8 #include "ColorTools.h"
14 B
, // (-1/2 sqrt(3)/2)
15 C
// (-1/2 -sqrt(3)/2)
19 * Renvoie la valeur (x1, x2) de l'itération suivante (i+1).
22 static void nextX(float x1
, float x2
, float& x1_next
, float& x2_next
)
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
;
27 // La matrice est représentée comme cela :
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);
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
;
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
);
46 * Renvoie la distance entre deux vecteurs a et b.
49 static float distance_carre(float a1
, float a2
, float b1
, float b2
)
51 return powf(a1
- b1
, 2.0) + powf(a2
- b2
, 2.0);
55 static float distance(float a1
, float a2
, float b1
, float b2
)
57 return (powf(a1
- b1
, 2.0) + powf(a2
- b2
, 2.0)) / (powf(b1
, 2.0) + powf(b2
, 2.0));
62 * n est le nombre d'iteration.
65 NewtonMath(int n
= 1000)
71 virtual ~NewtonMath() {}
74 void colorXY(uchar4
* ptrColor
, float x1
, float x2
) const
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;
83 const float epsilon
= 0.001;
85 float nearest_current_solution_distance
= FLT_MAX
;
86 Solution nearest_current_solution
= A
;
88 for (int i
= 0; i
< this->n
; i
++)
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
);
94 if (distance_to_A
< nearest_current_solution_distance
&& distance_to_A
< distance_to_B
&& distance_to_A
< distance_to_C
)
96 nearest_current_solution
= A
;
97 nearest_current_solution_distance
= distance_to_A
;
99 else if (distance_to_B
< nearest_current_solution_distance
&& distance_to_B
< distance_to_A
&& distance_to_B
< distance_to_C
)
101 nearest_current_solution
= B
;
102 nearest_current_solution_distance
= distance_to_B
;
104 else if (distance_to_C
< nearest_current_solution_distance
&& distance_to_C
< distance_to_A
&& distance_to_C
< distance_to_B
)
106 nearest_current_solution
= C
;
107 nearest_current_solution_distance
= distance_to_C
;
110 if (nearest_current_solution_distance
< epsilon
)
113 nextX(x1
, x2
, x1
, x2
);
116 switch (nearest_current_solution
)