#include <stdio.h>\r
using namespace std;\r
\r
+#include "NewtonDevice.h"\r
+\r
#include "Indice2D.h"\r
#include "IndiceTools.h"\r
-#include "DomaineMath.h"\r
#include "cudaTools.h"\r
#include "Device.h"\r
\r
#include "NewtonMath.h"\r
\r
-__global__ void newton(uchar4* ptrDevPixels, int w, int h, DomaineMath domaineMath)\r
+__global__\r
+void newton(uchar4* ptrDevPixels, int w, int h, DomaineMath domaineMath, int n, float epsilon)\r
{\r
const int TID = Indice2D::tid();\r
const int NB_THREAD = Indice2D::nbThread();\r
const int WH = w * h;\r
\r
- NewtonMath newtonMath;\r
+ NewtonMath newtonMath(n, epsilon);\r
\r
uchar4 color;\r
color.w = 255; // Par défaut, l'image est opaque.\r
\r
// (i,j) domaine écran.\r
// (x,y) domaine math.\r
- domaineMath.toXY(pixelI, pixelJ, &x, &y); // (i,j) -> (x,y)\r
+ domaineMath.toXY(pixelI, pixelJ, &x, &y); // (i,j) -> (x,y).\r
\r
newtonMath.colorXY(&color, x, y);\r
\r
#ifndef NEWTON_MATH_H_\r
#define NEWTON_MATH_H_\r
\r
-#include <cmath>\r
#include <float.h>\r
#include <stdio.h>\r
\r
C // (-1/2 -sqrt(3)/2)\r
};\r
\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 = powf(x1, 3.0) - 3.0 * x1 * powf(x2, 2.0) - 1.0;\r
- float f_x2 = powf(x2, 3.0) - 3.0 * powf(x1, 2.0) * 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 * powf(x1, 2.0) - 3.0 * powf(x2, 2.0);\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 * powf(x1, 2.0) + 3.0 * powf(x2, 2.0);\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
- return powf(a1 - b1, 2.0) + powf(a2 - b2, 2.0);\r
- }\r
-\r
- __device__\r
- static float distance(float a1, float a2, float b1, float b2)\r
- {\r
- return (powf(a1 - b1, 2.0) + powf(a2 - b2, 2.0)) / (powf(b1, 2.0) + powf(b2, 2.0));\r
- }\r
-\r
public:\r
/*\r
* n est le nombre d'iteration.\r
*/\r
__device__\r
- NewtonMath(int n = 1000)\r
- : n(n)\r
+ NewtonMath(int n, float epsilon)\r
+ : n(n), epsilon(epsilon)\r
{\r
}\r
\r
const float C1 = -1.0 / 2.0;\r
const float C2 = -sqrt(3.0) / 2.0;\r
\r
- const float epsilon = 0.001;\r
+ //const float epsilon = 0.1;\r
\r
float nearest_current_solution_distance = FLT_MAX;\r
Solution nearest_current_solution = A;\r
nearest_current_solution_distance = distance_to_C;\r
}\r
\r
- if (nearest_current_solution_distance < epsilon)\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
+ if (nearest_current_solution_distance < this->epsilon)\r
break;\r
\r
nextX(x1, x2, x1, x2);\r
}\r
\r
private:\r
+\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 * x2 - 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
\r
#include "Newton.h"\r
#include "Device.h"\r
-\r
-extern __global__ void newton(uchar4* ptrDevPixels, int w, int h, DomaineMath domaineMath);\r
+#include "NewtonDevice.h"\r
\r
Newton::Newton(int w, int h)\r
- : variateurAnimation(IntervalF(30, 100), 1),\r
+ : variateurN(IntervalI(5, 1000), 1),\r
+ variateurEpsilon(IntervalF(0.01, 10), 0.01),\r
w(w), h(h),\r
dg(8, 8, 1),\r
- db(16, 16, 1),\r
+ db(32, 32, 1),\r
ptrDomaineMathInit(new DomaineMath(-2, -2, 2, 2)),\r
title("Fractal Newton")\r
{\r
- // print(dg, db);\r
Device::assertDim(dg, db);\r
}\r
\r
\r
void Newton::runGPU(uchar4* ptrDevPixels, const DomaineMath& domaineMath)\r
{\r
- newton<<<dg,db>>>(ptrDevPixels, this->w, this->h, domaineMath);\r
+ newton<<<dg,db>>>(ptrDevPixels, this->w, this->h, domaineMath, this->t, this->epsilon);\r
+\r
+ HANDLE_ERROR(cudaDeviceSynchronize()); // Pour flusher les 'printf' (pour le DEBUG).\r
}\r
\r
void Newton::animationStep()\r
{\r
- this->t = this->variateurAnimation.varierAndGet();\r
+ this->t = this->variateurN.varierAndGet();\r
+ this->epsilon = this->variateurEpsilon.varierAndGet();\r
}\r
\r
int Newton::getW()\r
\r
float Newton::getT()\r
{\r
- return 0;\r
+ return this->t;\r
}\r
\r
string Newton::getTitle()\r
//Viewer<Image, RipplingProvider> rippling0(true, true, 10, 10);\r
//Viewer<ImageFonctionel, MandelbrotProvider> fractalMandelbrot(true, true, 20, 20);\r
//Viewer<ImageFonctionel, JuliaProvider> fractalJulia(true, true, 30, 30);\r
- //Viewer<ImageFonctionel, NewtonProvider> newtown(true, true, 20, 20);\r
- Viewer<Image, HeatTransfertProvider> heatTransfert(true, false, 20, 20);\r
+ Viewer<ImageFonctionel, NewtonProvider> newtown(true, true, 20, 20);\r
+ //Viewer<Image, HeatTransfertProvider> heatTransfert(true, false, 20, 20);\r
\r
GLUTImageViewers::runALL(); // Bloquant, Tant qu'une fenetre est ouverte\r
\r