#ifndef NEWTON_MATH_H_\r
#define NEWTON_MATH_H_\r
\r
-#include <cmath>\r
+#include <float.h>\r
+#include <stdio.h>\r
\r
-#include "CalibreurF.h"\r
#include "ColorTools.h"\r
\r
class NewtonMath\r
{\r
+ enum Solution {\r
+ A, // (1 0)\r
+ B, // (-1/2 sqrt(3)/2)\r
+ C // (-1/2 -sqrt(3)/2)\r
+ };\r
+\r
+ // Les trois solutions du système d'équation : A, B et C (vecteur 2D).\r
+ static const float A1 = 1.0;\r
+ static const float A2 = 0.0;\r
+ static const float B1 = -0.5; // -1.0 / 2.0\r
+ static const float B2 = 0.866025403785; // sqrt(3.0) / 2.0\r
+ static const float C1 = -0.5; // -1.0 / 2.0\r
+ static const float C2 = -0.866025403785; // -sqrt(3.0) / 2.0\r
+\r
public:\r
+ /*\r
+ * n est le nombre d'iteration.\r
+ */\r
__device__\r
- NewtonMath()\r
- : calibreur(IntervalF(1, 100), IntervalF(0, 1))\r
+ NewtonMath(int n, float epsilon)\r
+ : n(n), epsilon(epsilon)\r
{\r
}\r
\r
__device__\r
virtual ~NewtonMath() {}\r
\r
- public:\r
- /**\r
- * x=pixelI\r
- * y=pixelJ\r
- */\r
__device__\r
- void colorXY(uchar4* ptrColor, float x, float y) const\r
+ void colorXY(uchar4* ptrColor, float x1, float x2) const\r
{\r
- ptrColor->x = 0;\r
- ptrColor->y = 0;\r
- ptrColor->z = 0;\r
-\r
- int i = 0;\r
- float s = static_cast<float>(i);\r
- this->calibreur.calibrer(s);\r
- ColorTools::HSB_TO_RVB(s, ptrColor);\r
+ float nearest_current_solution_distance = FLT_MAX;\r
+ Solution nearest_current_solution = A;\r
+\r
+ for (int i = 0; i < this->n; i++)\r
+ {\r
+ float distance_to_A = distance(x1, x2, A1, A2);\r
+ float distance_to_B = distance(x1, x2, B1, B2);\r
+ float distance_to_C = distance(x1, x2, C1, C2);\r
+\r
+ if (distance_to_A < nearest_current_solution_distance && distance_to_A < distance_to_B && distance_to_A < distance_to_C)\r
+ {\r
+ nearest_current_solution = A;\r
+ nearest_current_solution_distance = distance_to_A;\r
+ }\r
+ else if (distance_to_B < nearest_current_solution_distance && distance_to_B < distance_to_A && distance_to_B < distance_to_C)\r
+ {\r
+ nearest_current_solution = B;\r
+ nearest_current_solution_distance = distance_to_B;\r
+ }\r
+ else if (distance_to_C < nearest_current_solution_distance && distance_to_C < distance_to_A && distance_to_C < distance_to_B)\r
+ {\r
+ nearest_current_solution = C;\r
+ nearest_current_solution_distance = distance_to_C;\r
+ }\r
+\r
+ /* For DEBUG.\r
+ * if (Indice2D::tid() == 0)\r
+ {\r
+ printf("nearest_current_solution_distance: %f\n", nearest_current_solution_distance);\r
+ }*/\r
+\r
+ /*printf("x1: %f, x2: %f\n", x1, x2);\r
+ printf("d to a: %f\n", distance_to_A);\r
+ printf("d to a: %f\n", distance_to_B);\r
+ printf("d to a: %f\n", distance_to_C);\r
+ }\r
+ */\r
+\r
+ if (nearest_current_solution_distance < this->epsilon)\r
+ break;\r
+\r
+ nextX(x1, x2, x1, x2);\r
+ }\r
+\r
+ switch (nearest_current_solution)\r
+ {\r
+ // Noir.\r
+ case A :\r
+ ptrColor->x = 0;\r
+ ptrColor->y = 0;\r
+ ptrColor->z = 0;\r
+ break;\r
+ // Gris.\r
+ case B :\r
+ ptrColor->x = 128;\r
+ ptrColor->y = 128;\r
+ ptrColor->z = 128;\r
+ break;\r
+ // Blanc.\r
+ case C :\r
+ ptrColor->x = 255;\r
+ ptrColor->y = 255;\r
+ ptrColor->z = 255;\r
+ break;\r
+ }\r
}\r
\r
private:\r
\r
- private:\r
- CalibreurF calibreur;\r
+ /*\r
+ * Renvoie la valeur (x1, x2) de l'itération suivante (i+1).\r
+ */\r
+ __device__\r
+ static void nextX(float x1, float x2, float& x1_next, float& x2_next)\r
+ {\r
+ float f_x1 = x1 * x1 * x1 - 3.0 * x1 * x2 * x2 - 1.0;\r
+ float f_x2 = x2 * x2 * x2 - 3.0 * x1 * x1 * x2;\r
+\r
+ // La matrice est représentée comme cela :\r
+ // a b\r
+ // c d\r
+ float jacobienne_f_x_a = 3.0 * x1 * x1 - 3.0 * x2 * x2;\r
+ float jacobienne_f_x_b = -6.0 * x1 * x2;\r
+ float jacobienne_f_x_c = -6.0 * x1 * x2;\r
+ float jacobienne_f_x_d = -3.0 * x1 * x1 + 3.0 * x2 * x2;\r
+\r
+ float det_inverse_jacobienne = 1.0 / (jacobienne_f_x_a * jacobienne_f_x_d - jacobienne_f_x_b * jacobienne_f_x_c);\r
+ float jacobienne_f_x_a_inverse = jacobienne_f_x_d * det_inverse_jacobienne;\r
+ float jacobienne_f_x_b_inverse = -jacobienne_f_x_b * det_inverse_jacobienne;\r
+ float jacobienne_f_x_c_inverse = -jacobienne_f_x_c * det_inverse_jacobienne;\r
+ float jacobienne_f_x_d_inverse = jacobienne_f_x_a * det_inverse_jacobienne;\r
+\r
+ x1_next = x1 - (jacobienne_f_x_a_inverse * f_x1 + jacobienne_f_x_b_inverse * f_x2);\r
+ x2_next = x2 - (jacobienne_f_x_c_inverse * f_x1 + jacobienne_f_x_d_inverse * f_x2);\r
+ }\r
+\r
+ /*\r
+ * Renvoie la distance entre deux vecteurs a et b.\r
+ */\r
+ __device__\r
+ static float distance_carre(float a1, float a2, float b1, float b2)\r
+ {\r
+ const float delta1 = a1 - b1;\r
+ const float delta2 = a2 - b2;\r
+ return delta1 * delta1 + delta2 * delta2;\r
+ }\r
+\r
+ __device__\r
+ static float distance(float a1, float a2, float b1, float b2)\r
+ {\r
+ const float delta1 = a1 - b1;\r
+ const float delta2 = a2 - b2;\r
+\r
+ return (delta1 * delta1 + delta2 * delta2) / (b1 * b1 + b2 * b2);\r
+ }\r
+\r
+ int n;\r
+ float epsilon;\r
};\r
\r
#endif\r
projects / GPU.git / blobdiff