7 #include "ColorTools.h"
13 B
, // (-1/2 sqrt(3)/2)
14 C
// (-1/2 -sqrt(3)/2)
19 * n est le nombre d'iteration.
22 NewtonMath(int n
, float epsilon
)
23 : n(n
), epsilon(epsilon
)
28 virtual ~NewtonMath() {}
31 void colorXY(uchar4
* ptrColor
, float x1
, float x2
) const
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;
40 //const float epsilon = 0.1;
42 float nearest_current_solution_distance
= FLT_MAX
;
43 Solution nearest_current_solution
= A
;
45 for (int i
= 0; i
< this->n
; i
++)
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
);
51 if (distance_to_A
< nearest_current_solution_distance
&& distance_to_A
< distance_to_B
&& distance_to_A
< distance_to_C
)
53 nearest_current_solution
= A
;
54 nearest_current_solution_distance
= distance_to_A
;
56 else if (distance_to_B
< nearest_current_solution_distance
&& distance_to_B
< distance_to_A
&& distance_to_B
< distance_to_C
)
58 nearest_current_solution
= B
;
59 nearest_current_solution_distance
= distance_to_B
;
61 else if (distance_to_C
< nearest_current_solution_distance
&& distance_to_C
< distance_to_A
&& distance_to_C
< distance_to_B
)
63 nearest_current_solution
= C
;
64 nearest_current_solution_distance
= distance_to_C
;
67 /*if (Indice2D::tid() == 0)
69 printf("nearest_current_solution_distance: %f\n", nearest_current_solution_distance);
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);
78 if (nearest_current_solution_distance
< this->epsilon
)
81 nextX(x1
, x2
, x1
, x2
);
84 switch (nearest_current_solution
)
110 * Renvoie la valeur (x1, x2) de l'itération suivante (i+1).
113 static void nextX(float x1
, float x2
, float& x1_next
, float& x2_next
)
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
;
118 // La matrice est représentée comme cela :
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
;
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
;
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
);
137 * Renvoie la distance entre deux vecteurs a et b.
140 static float distance_carre(float a1
, float a2
, float b1
, float b2
)
142 const float delta1
= a1
- b1
;
143 const float delta2
= a2
- b2
;
144 return delta1
* delta1
+ delta2
* delta2
;
148 static float distance(float a1
, float a2
, float b1
, float b2
)
150 const float delta1
= a1
- b1
;
151 const float delta2
= a2
- b2
;
153 return (delta1
* delta1
+ delta2
* delta2
) / (b1
* b1
+ b2
* b2
);