Newton, work in progress.
authorgburri <gregory.burri@master.hes-so.ch>
Sat, 6 Dec 2014 17:36:42 +0000 (18:36 +0100)
committergburri <gregory.burri@master.hes-so.ch>
Sat, 6 Dec 2014 17:36:42 +0000 (18:36 +0100)
WCudaMSE/Student_Cuda_Image/src/cpp/core/02_Mandelbrot_Julia/provider/FractalProvider.cpp
WCudaMSE/Student_Cuda_Image/src/cpp/core/02_Mandelbrot_Julia/provider/FractalProvider.h
WCudaMSE/Student_Cuda_Image/src/cpp/core/03_Newton/moo/device/NewtonDevice.cu
WCudaMSE/Student_Cuda_Image/src/cpp/core/03_Newton/moo/device/math/NewtonMath.h
WCudaMSE/Student_Cuda_Image/src/cpp/core/03_Newton/moo/host/Newton.cu
WCudaMSE/Student_Cuda_Image/src/cpp/core/mainGL.cpp

index 47df25e..5282d85 100755 (executable)
@@ -1,35 +1,6 @@
 #include "FractalProvider.h"\r
 \r
-\r
-/*----------------------------------------------------------------------*\\r
- |*                    Declaration                                     *|\r
- \*---------------------------------------------------------------------*/\r
-\r
-/*--------------------------------------*\\r
- |*            Imported                *|\r
- \*-------------------------------------*/\r
-\r
-/*--------------------------------------*\\r
- |*            Public                  *|\r
- \*-------------------------------------*/\r
-\r
-/*--------------------------------------*\\r
- |*            Private                 *|\r
- \*-------------------------------------*/\r
-\r
-/*----------------------------------------------------------------------*\\r
- |*                    Implementation                                  *|\r
- \*---------------------------------------------------------------------*/\r
-\r
-/*--------------------------------------*\\r
- |*            Public                  *|\r
- \*-------------------------------------*/\r
-\r
-/*-----------------*\\r
- |*    static     *|\r
- \*----------------*/\r
-\r
-Fractal* FractalProvider::createMandelbrot()\r
+Fractal* MandelbrotProvider::create()\r
     {\r
     int dw = 16 * 50;\r
     int dh = 16 * 30;\r
@@ -37,7 +8,7 @@ Fractal* FractalProvider::createMandelbrot()
     return new FractalMandelbrot(dw, dh, 0.2);\r
     }\r
 \r
-Fractal* FractalProvider::createJulia()\r
+Fractal* JuliaProvider::create()\r
     {\r
     int dw = 16 * 50;\r
     int dh = 16 * 30;\r
@@ -45,23 +16,15 @@ Fractal* FractalProvider::createJulia()
     return new FractalJulia(dw, dh, 0.01, -0.745, -0.32, -0.09, 0.1);\r
     }\r
 \r
-ImageFonctionel* FractalProvider::createMandelbrotGL()\r
+ImageFonctionel* MandelbrotProvider::createGL()\r
     {\r
     ColorRGB_01* ptrColorTitre = new ColorRGB_01(0, 0, 0);\r
-    return new ImageFonctionel(createMandelbrot(), ptrColorTitre); // both ptr destroy by destructor of ImageFonctionel\r
+    return new ImageFonctionel(create(), ptrColorTitre); // both ptr destroy by destructor of ImageFonctionel\r
     }\r
 \r
 \r
-ImageFonctionel* FractalProvider::createJuliaGL()\r
+ImageFonctionel* JuliaProvider::createGL()\r
     {\r
     ColorRGB_01* ptrColorTitre = new ColorRGB_01(0, 0, 0);\r
-    return new ImageFonctionel(createJulia(), ptrColorTitre); // both ptr destroy by destructor of ImageFonctionel\r
+    return new ImageFonctionel(create(), ptrColorTitre); // both ptr destroy by destructor of ImageFonctionel\r
     }\r
-\r
-/*--------------------------------------*\\r
- |*            Private                 *|\r
- \*-------------------------------------*/\r
-\r
-/*----------------------------------------------------------------------*\\r
- |*                    End                                             *|\r
- \*---------------------------------------------------------------------*/\r
index 27bd73e..eca7b3c 100755 (executable)
@@ -4,27 +4,19 @@
 #include "Fractal.h"\r
 #include "ImageFonctionel.h"\r
 \r
-/*----------------------------------------------------------------------*\\r
- |*                    Declaration                                     *|\r
- \*---------------------------------------------------------------------*/\r
+class MandelbrotProvider\r
+    {\r
+    public:\r
+       static Fractal* create();\r
+       static ImageFonctionel* createGL();\r
+    };\r
 \r
-/*--------------------------------------*\\r
- |*            Public                  *|\r
- \*-------------------------------------*/\r
 \r
-class FractalProvider\r
+class JuliaProvider\r
     {\r
     public:\r
-       static Fractal* createMandelbrot();\r
-       static Fractal* createJulia();\r
-\r
-       static ImageFonctionel* createMandelbrotGL();\r
-       static ImageFonctionel* createJuliaGL();\r
+        static Fractal* create();\r
+        static ImageFonctionel* createGL();\r
     };\r
 \r
 #endif\r
-\r
-/*----------------------------------------------------------------------*\\r
- |*                    End                                             *|\r
- \*---------------------------------------------------------------------*/\r
-\r
index 536b75b..3399d92 100755 (executable)
@@ -1,4 +1,6 @@
 #include <iostream>\r
+#include <stdio.h>\r
+using namespace std;\r
 \r
 #include "Indice2D.h"\r
 #include "IndiceTools.h"\r
@@ -8,9 +10,6 @@
 \r
 #include "NewtonMath.h"\r
 \r
-using std::cout;\r
-using std::endl;\r
-\r
 __global__ void newton(uchar4* ptrDevPixels, int w, int h, DomaineMath domaineMath)\r
     {\r
     const int TID = Indice2D::tid();\r
@@ -20,7 +19,7 @@ __global__ void newton(uchar4* ptrDevPixels, int w, int h, DomaineMath domaineMa
     NewtonMath newtonMath;\r
 \r
     uchar4 color;\r
-    color.z = 255; // Par défaut, l'image est opaque.\r
+    color.w = 255; // Par défaut, l'image est opaque.\r
 \r
     double x, y;\r
     int pixelI, pixelJ;\r
@@ -30,8 +29,8 @@ __global__ void newton(uchar4* ptrDevPixels, int w, int h, DomaineMath domaineMa
         {\r
         IndiceTools::toIJ(s, w, &pixelI, &pixelJ); // update (pixelI, pixelJ)\r
 \r
-        // (i,j) domaine ecran\r
-        // (x,y) domaine math\r
+        // (i,j) domaine écran.\r
+        // (x,y) domaine math.\r
         domaineMath.toXY(pixelI, pixelJ, &x, &y); //  (i,j) -> (x,y)\r
 \r
         newtonMath.colorXY(&color, x, y);\r
index 5edf6ee..3c76170 100755 (executable)
 #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
+        /*\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()\r
-           : calibreur(IntervalF(1, 100), IntervalF(0, 1))\r
+       NewtonMath(int n = 1000)\r
+           : n(n)\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
-           }\r
+            const float A1 = 1.0;\r
+            const float A2 = 0.0;\r
+            const float B1 = -1.0 / 2.0;\r
+            const float B2 = sqrt(3.0) / 2.0;\r
+            const float C1 = -1.0 / 2.0;\r
+            const float C2 = -sqrt(3.0) / 2.0;\r
 \r
-    private:\r
+            const float epsilon = 0.001;\r
+\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
+                if (nearest_current_solution_distance < 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
-       CalibreurF calibreur;\r
+       int n;\r
     };\r
 \r
 #endif\r
index 028f924..685fb5e 100755 (executable)
@@ -1,12 +1,11 @@
 #include <iostream>\r
 #include <assert.h>\r
+#include <stdio.h>\r
+using namespace std;\r
 \r
 #include "Newton.h"\r
 #include "Device.h"\r
 \r
-using std::cout;\r
-using std::endl;\r
-\r
 extern __global__ void newton(uchar4* ptrDevPixels, int w, int h, DomaineMath domaineMath);\r
 \r
 Newton::Newton(int w, int h)\r
@@ -14,9 +13,10 @@ Newton::Newton(int w, int h)
       w(w), h(h),\r
       dg(8, 8, 1),\r
       db(16, 16, 1),\r
+      ptrDomaineMathInit(new DomaineMath(-2, -2, 2, 2)),\r
       title("Fractal Newton")\r
     {\r
-    //print(dg, db);\r
+    // print(dg, db);\r
     Device::assertDim(dg, db);\r
     }\r
 \r
index 3d4d02f..d2b6adf 100755 (executable)
@@ -1,6 +1,7 @@
 #include <iostream>\r
 #include <stdlib.h>\r
 #include <string.h>\r
+using namespace std;\r
 \r
 #include "GLUTImageViewers.h"\r
 \r
 \r
 #include "Rippling0Provider.h"\r
 #include "RipplingProvider.h"\r
-\r
 #include "FractalProvider.h"\r
+#include "NewtonProvider.h"\r
 \r
-using std::cout;\r
-using std::endl;\r
-using std::string;\r
 \r
-class RipplingViewer\r
+template <class TOutput, class TProvider>\r
+class Viewer\r
     {\r
+    private:\r
+        TOutput* ptrProvider;\r
+        GLUTImageViewers viewer;\r
+\r
     public:\r
-        RipplingViewer()\r
-            : ptrRippling0(Rippling0Provider::createGL()), ptrRippling(RipplingProvider::createGL()),\r
-              rippling0Viewer(this->ptrRippling0, true, true, 0, 0),\r
-              ripplingViewer(this->ptrRippling, true, true, 20, 20)\r
-            {}\r
-        ~RipplingViewer()\r
+        Viewer(bool isAnimation, bool isSelection, int pxFrame, int pyFrame):\r
+            ptrProvider(TProvider::createGL()),\r
+            viewer(ptrProvider, isAnimation, isSelection, pxFrame, pyFrame)\r
             {\r
-                delete this->ptrRippling0;\r
-                delete this->ptrRippling;\r
             }\r
-    private:\r
-        Rippling0Image* ptrRippling0;\r
-        Image* ptrRippling;\r
-        GLUTImageViewers rippling0Viewer, ripplingViewer;\r
-    };\r
 \r
-class FractalViewer\r
-    {\r
-    public:\r
-        FractalViewer()\r
-            : ptrMandelbrot(FractalProvider::createMandelbrotGL()), ptrJulia(FractalProvider::createJuliaGL()),\r
-              mandelbrotViewer(this->ptrMandelbrot, true, true, 0, 0),\r
-              juliaViewer(this->ptrJulia, true, true, 20, 20)\r
-            {}\r
-        ~FractalViewer()\r
+        ~Viewer()\r
             {\r
-                delete this->ptrMandelbrot;\r
-                delete this->ptrJulia;\r
+            delete this->ptrProvider;\r
             }\r
-    private:\r
-        ImageFonctionel* ptrMandelbrot;\r
-        ImageFonctionel* ptrJulia;\r
-        GLUTImageViewers mandelbrotViewer, juliaViewer;\r
     };\r
 \r
 int mainGL(void)\r
     {\r
-    // RipplingViewer rippling;\r
-    // FractalViewer fractals;\r
+    //Viewer<Rippling0Image, Rippling0Provider> rippling0(true, true, 10, 10);\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
+\r
+\r
+    Viewer<ImageFonctionel, NewtonProvider> newtown(true, true, 20, 20);\r
+\r
 \r
     GLUTImageViewers::runALL(); // Bloquant, Tant qu'une fenetre est ouverte\r
 \r