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