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