7 #include "ColorTools.h"
13 B
, // (-1/2 sqrt(3)/2)
14 C
// (-1/2 -sqrt(3)/2)
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
27 * n est le nombre d'iteration.
30 NewtonMath(int n
, float epsilon
)
31 : n(n
), epsilon(epsilon
)
36 virtual ~NewtonMath() {}
39 void colorXY(uchar4
* ptrColor
, float x1
, float x2
) const
41 float nearest_current_solution_distance
= FLT_MAX
;
42 Solution nearest_current_solution
= A
;
44 for (int i
= 0; i
< this->n
; i
++)
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
);
50 if (distance_to_A
< nearest_current_solution_distance
&& distance_to_A
< distance_to_B
&& distance_to_A
< distance_to_C
)
52 nearest_current_solution
= A
;
53 nearest_current_solution_distance
= distance_to_A
;
55 else if (distance_to_B
< nearest_current_solution_distance
&& distance_to_B
< distance_to_A
&& distance_to_B
< distance_to_C
)
57 nearest_current_solution
= B
;
58 nearest_current_solution_distance
= distance_to_B
;
60 else if (distance_to_C
< nearest_current_solution_distance
&& distance_to_C
< distance_to_A
&& distance_to_C
< distance_to_B
)
62 nearest_current_solution
= C
;
63 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);
79 if (nearest_current_solution_distance
< this->epsilon
)
82 nextX(x1
, x2
, x1
, x2
);
85 switch (nearest_current_solution
)
111 * Renvoie la valeur (x1, x2) de l'itération suivante (i+1).
114 static void nextX(float x1
, float x2
, float& x1_next
, float& x2_next
)
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
;
119 // La matrice est représentée comme cela :
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
;
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
;
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
);
138 * Renvoie la distance entre deux vecteurs a et b.
141 static float distance_carre(float a1
, float a2
, float b1
, float b2
)
143 const float delta1
= a1
- b1
;
144 const float delta2
= a2
- b2
;
145 return delta1
* delta1
+ delta2
* delta2
;
149 static float distance(float a1
, float a2
, float b1
, float b2
)
151 const float delta1
= a1
- b1
;
152 const float delta2
= a2
- b2
;
154 return (delta1
* delta1
+ delta2
* delta2
) / (b1
* b1
+ b2
* b2
);