Newton, work in progress.
[GPU.git] / 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