WCudaMSE/Student_Cuda_Image/src/cpp/core/03_Newton/moo/device/math/NewtonMath.h

   1 #ifndef NEWTON_MATH_H_
   2 #define NEWTON_MATH_H_
   3
   4 #include <float.h>
   5 #include <stdio.h>
   6
   7 #include "ColorTools.h"
   8
   9 class NewtonMath
  10     {
  11         enum Solution {
  12             A, // (1 0)
  13             B, // (-1/2 sqrt(3)/2)
  14             C // (-1/2 -sqrt(3)/2)
  15         };
  16
  17         // Les trois solutions du système d'équation : A, B et C (vecteur 2D).
  18         static const float A1 = 1.0;
  19         static const float A2 = 0.0;
  20         static const float B1 = -0.5; // -1.0 / 2.0
  21         static const float B2 = 0.866025403785; // sqrt(3.0) / 2.0
  22         static const float C1 = -0.5; // -1.0 / 2.0
  23         static const float C2 = -0.866025403785; // -sqrt(3.0) / 2.0
  24
  25     public:
  26         /*
  27          * n est le nombre d'iteration.
  28          */
  29         __device__
  30         NewtonMath(int n, float epsilon)
  31             : n(n), epsilon(epsilon)
  32             {
  33             }
  34
  35         __device__
  36         virtual ~NewtonMath() {}
  37
  38         __device__
  39         void colorXY(uchar4* ptrColor, float x1, float x2) const
  40             {
  41             float nearest_current_solution_distance = FLT_MAX;
  42             Solution nearest_current_solution = A;
  43
  44             for (int i = 0; i < this->n; i++)
  45                 {
  46                 float distance_to_A = distance(x1, x2, A1, A2);
  47                 float distance_to_B = distance(x1, x2, B1, B2);
  48                 float distance_to_C = distance(x1, x2, C1, C2);
  49
  50                 if (distance_to_A < nearest_current_solution_distance && distance_to_A < distance_to_B && distance_to_A < distance_to_C)
  51                     {
  52                     nearest_current_solution = A;
  53                     nearest_current_solution_distance = distance_to_A;
  54                     }
  55                 else if (distance_to_B < nearest_current_solution_distance && distance_to_B < distance_to_A && distance_to_B < distance_to_C)
  56                     {
  57                     nearest_current_solution = B;
  58                     nearest_current_solution_distance = distance_to_B;
  59                     }
  60                 else if (distance_to_C < nearest_current_solution_distance && distance_to_C < distance_to_A && distance_to_C < distance_to_B)
  61                     {
  62                     nearest_current_solution = C;
  63                     nearest_current_solution_distance = distance_to_C;
  64                     }
  65
  66                 /* For DEBUG.
  67                  * if (Indice2D::tid() == 0)
  68                     {
  69                     printf("nearest_current_solution_distance: %f\n", nearest_current_solution_distance);
  70                     }*/
  71
  72                     /*printf("x1: %f, x2: %f\n", x1, x2);
  73                     printf("d to a: %f\n", distance_to_A);
  74                     printf("d to a: %f\n", distance_to_B);
  75                     printf("d to a: %f\n", distance_to_C);
  76                     }
  77                  */
  78
  79                 if (nearest_current_solution_distance < this->epsilon)
  80                     break;
  81
  82                 nextX(x1, x2, x1, x2);
  83                 }
  84
  85             switch (nearest_current_solution)
  86                 {
  87                 // Noir.
  88                 case A :
  89                     ptrColor->x = 0;
  90                     ptrColor->y = 0;
  91                     ptrColor->z = 0;
  92                     break;
  93                 // Gris.
  94                 case B :
  95                     ptrColor->x = 128;
  96                     ptrColor->y = 128;
  97                     ptrColor->z = 128;
  98                     break;
  99                 // Blanc.
 100                 case C :
 101                     ptrColor->x = 255;
 102                     ptrColor->y = 255;
 103                     ptrColor->z = 255;
 104                     break;
 105                 }
 106             }
 107
 108     private:
 109
 110         /*
 111          * Renvoie la valeur (x1, x2) de l'itération suivante (i+1).
 112          */
 113         __device__
 114         static void nextX(float x1, float x2, float& x1_next, float& x2_next)
 115             {
 116             float f_x1 = x1 * x1 * x1 - 3.0 * x1 * x2 * x2 - 1.0;
 117             float f_x2 = x2 * x2 * x2 - 3.0 * x1 * x1 * x2;
 118
 119             // La matrice est représentée comme cela :
 120             // a b
 121             // c d
 122             float jacobienne_f_x_a = 3.0 * x1 * x1 - 3.0 * x2 * x2;
 123             float jacobienne_f_x_b = -6.0 * x1 * x2;
 124             float jacobienne_f_x_c = -6.0 * x1 * x2;
 125             float jacobienne_f_x_d = -3.0 * x1 * x1 + 3.0 * x2 * x2;
 126
 127             float det_inverse_jacobienne = 1.0 / (jacobienne_f_x_a * jacobienne_f_x_d - jacobienne_f_x_b * jacobienne_f_x_c);
 128             float jacobienne_f_x_a_inverse = jacobienne_f_x_d * det_inverse_jacobienne;
 129             float jacobienne_f_x_b_inverse = -jacobienne_f_x_b * det_inverse_jacobienne;
 130             float jacobienne_f_x_c_inverse = -jacobienne_f_x_c * det_inverse_jacobienne;
 131             float jacobienne_f_x_d_inverse = jacobienne_f_x_a * det_inverse_jacobienne;
 132
 133             x1_next = x1 - (jacobienne_f_x_a_inverse * f_x1 + jacobienne_f_x_b_inverse * f_x2);
 134             x2_next = x2 - (jacobienne_f_x_c_inverse * f_x1 + jacobienne_f_x_d_inverse * f_x2);
 135             }
 136
 137         /*
 138          * Renvoie la distance entre deux vecteurs a et b.
 139          */
 140         __device__
 141         static float distance_carre(float a1, float a2, float b1, float b2)
 142             {
 143             const float delta1 = a1 - b1;
 144             const float delta2 = a2 - b2;
 145             return delta1 * delta1 + delta2 * delta2;
 146             }
 147
 148         __device__
 149         static float distance(float a1, float a2, float b1, float b2)
 150             {
 151             const float delta1 = a1 - b1;
 152             const float delta2 = a2 - b2;
 153
 154             return (delta1 * delta1 + delta2 * delta2) / (b1 * b1 + b2 * b2);
 155             }
 156
 157         int n;
 158         float epsilon;
 159     };
 160
 161 #endif