Entered interface of mapping to MIPAV image in runAlgorithm.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Fri, 8 Dec 2017 22:58:17 +0000 (22:58 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Fri, 8 Dec 2017 22:58:17 +0000 (22:58 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15303 ba61647d-9d00-f842-95cd-605cb4296b96

mipav/src/gov/nih/mipav/model/algorithms/DoublyConnectedSC.java

index 53c5dc3..cceff8a 100644 (file)
@@ -1,5 +1,9 @@
 package gov.nih.mipav.model.algorithms;\r
 \r
+import java.io.IOException;\r
+\r
+import gov.nih.mipav.model.structures.ModelImage;\r
+import gov.nih.mipav.view.MipavUtil;\r
 \r
 public class DoublyConnectedSC extends AlgorithmBase {\r
        // This is a port of the FORTRAN ACM TOMS Algorithm 785, collected\r
@@ -77,48 +81,80 @@ public class DoublyConnectedSC extends AlgorithmBase {
     // COMMON /PARAM5/ISPRT,ICOUNT\r
        \r
        private SchwarzChristoffelMapping scm = new SchwarzChristoffelMapping();\r
-       \r
+       // NOTE2:    USERS ARE RESPONSIBLE FOR CHECKING THE FOLLOWING:\r
+    //   1. EACH COMPONENT OF THE OUTER POLYGON CONTAINS AT LEAST ONE\r
+       //      FINITE VERTEX.\r
+       //   2. Z1[N-1] IS THE CLOSEST (OR ALMOST CLOSEST)\r
+       //      INNER VERTEX TO Z0[M-1].\r
+    //   3. COUNTERCLOCKWISE ORIENTATION OF THE VERTICES.\r
+\r
+       // M[0] points of outer polygon\r
+       private double Z0[][] = null;\r
+       // N[0] points of inner polygon\r
+       private double Z1[][] = null;\r
        // Specifies the geometry of the polygon region for 7 test examples\r
-       private int IPOLY;\r
+       private int IPOLY = -1;\r
        // The number of Gauss-Jacobi points\r
        // Recommended values for NPTQ are 2-8.\r
-       private int NPTQ;\r
+       private int NPTQ = 8;\r
        // 1 for solving the nonlinear system, 2 for not solving the nonlinear system\r
-       private int ISOLV;\r
+       private int ISOLV = 1;\r
+       //IGUESS    (=0 )A NON-EQUALLY SPACED INITIAL GUESS OR(=1)THE\r
+       //C            OTHER  EQUALLY-SPACED  INITIAL GUESS PROVIDED, OR\r
+       //C            (=2)USER-SUPPLIED INITIAL GUESS WHICH IS REQUIRED\r
+       //C            TO BE THE ARGUMENTS OF THE INITIAL PREVERTICES.\r
+       //C            ROUTINE ARGUM MAY BE USED FOR COMPUTING ARGUMENTS\r
+       //C            IF NEEDED. NOTE: C WILL BE COMPUTED IN THE ROUTINE\r
+       //C            (NOT SUPPLIED!)\r
+    private int IGUESS = 1;\r
+    // LINEARC   INTEGER VARIABLE TO CONTROL INTEGRATION PATH.IN PATICULAR:\r
+    //C            LINEARC=0:INTEGRATING ALONG LINE SEGMENT WHENEVER POSSIBLE\r
+    //C            LINEARC=1:INTEGRATING ALONG CIRCULAR ARC WHENEVER POSSIBLE\r
+    private int LINEARC = 1;\r
+    // TOL       A TOLERANCE TO CONTROL THE CONVERGENCE IN HYBRD\r
+    private double TOL = 1.0E-10;\r
        // Inverse points\r
        private double INVERSE_POINTS[][] = null;\r
        private double FORWARD_POINTS[][] = null;\r
        private boolean testRoutine = false;\r
        private double MACHEP = 2.2204460E-16;\r
+       private int M[] = new int[1];\r
+       private int N[] = new int[1];\r
+       // ISHAPE = 0 outer polygon has no infinite vertices\r
+    // ISHAPE = 1 outer polygon has some infinite vertices\r
+       int ISHAPE = 0;\r
        \r
        public DoublyConnectedSC() {\r
                \r
        }\r
        \r
-       public DoublyConnectedSC(int IPOLY, int NPTQ, int ISOLV, int ISPRT, double INVERSE_POINTS[][]) {\r
+       public DoublyConnectedSC(int IPOLY, int NPTQ, int ISOLV, int ISPRT) {\r
                this.IPOLY = IPOLY;\r
                this.NPTQ = NPTQ;\r
                this.ISOLV = ISOLV;\r
                this.ISPRT = ISPRT;\r
-               this.INVERSE_POINTS = INVERSE_POINTS;\r
                testRoutine = true;\r
        }\r
        \r
+       public DoublyConnectedSC(ModelImage destImg, ModelImage srcImg, double Z0[][], double Z1[][], int NPTQ, int ISPRT,\r
+                       int IGUESS, int LINEARC, double TOL) {\r
+               super(destImg, srcImg); \r
+               this.Z0 = Z0;\r
+               this.Z1 = Z1;\r
+               this.NPTQ = NPTQ;\r
+               this.ISPRT = ISPRT;\r
+               this.IGUESS = IGUESS;\r
+               this.LINEARC = LINEARC;\r
+               this.TOL = TOL;\r
+               M[0] = Z0.length;\r
+               N[0] = Z1.length;\r
+       }\r
+       \r
        public void runAlgorithm() {\r
                int I;\r
-               int M[] = new int[1];\r
-               int N[] = new int[1];\r
-               double Z0[][] = new double[30][2];\r
-               double Z1[][] = new double[30][2];\r
                double ALFA0[] = new double[30];\r
                double ALFA1[] = new double[30];\r
                double QWORK[] = new double[1660];\r
-               // ISHAPE = 0 outer polygon has no infinite vertices\r
-               // ISHAPE = 1 outer polygon has some infinite vertices\r
-               int ISHAPE;\r
-               int IGUESS;\r
-               int LINEARC;\r
-               double TOL;\r
                double U[] = new double[1];\r
                double C[] = new double[2];\r
                double W0[][] = new double[30][2];\r
@@ -130,7 +166,37 @@ public class DoublyConnectedSC extends AlgorithmBase {
                double WW[];\r
                double ZZ0[];\r
                \r
-double neweps;\r
+               int i;\r
+               int cf;\r
+               int xDimSource;\r
+               int yDimSource;\r
+               int sourceSlice;\r
+               int xDimDest;\r
+               int yDimDest;\r
+               int destSlice;\r
+               double srcBuffer[];\r
+               double destBuffer[];\r
+               double imageMin;\r
+               int x;\r
+               int y; \r
+               int c;\r
+               double zp[];\r
+               int j;\r
+               int index;\r
+               double wp[][];\r
+               double xp;\r
+               double yp;\r
+               int x0;\r
+               int y0;\r
+               int deltaX;\r
+               int deltaY;\r
+               double dx;\r
+               double dy;\r
+               double dx1;\r
+               double dy1;\r
+               int position;\r
+               \r
+        double neweps;\r
                \r
                // MACHEP is a machine dependent parameter specifying the relative precision of\r
                // floating point arithmetic.\r
@@ -154,6 +220,8 @@ double neweps;
         } // while(true)\r
                \r
                if (testRoutine) {\r
+                       Z0 = new double[30][2];\r
+                       Z1 = new double[30][2];\r
                        DSCDATA(IPOLY, M, N, Z0, Z1, ALFA0, ALFA1);\r
                        ANGLES(N[0], Z1, ALFA1, 1);\r
                        if ((IPOLY == 1) || (IPOLY == 3) || (IPOLY == 4) || (IPOLY == 6) || (IPOLY == 8) || (IPOLY == 9)) {\r
@@ -239,9 +307,209 @@ double neweps;
                            System.out.println("INVERSE MAPPING YIELDS = (" + WW[0] + ", " + WW[1] + ")");\r
                        } // for (I = 0; I < FORWARD_POINTS.length; I++)\r
                    } // if ((FORWARD_POINTS != null) && (FORWARD_POINTS.length != 0))\r
+                   return;\r
                } // if (testRoutine)\r
                \r
+               if (srcImage == null) {\r
+            displayError("Source Image is null");\r
+            finalize();\r
+\r
+            return;\r
+        }\r
+\r
+        fireProgressStateChanged(srcImage.getImageName(), "Between 2 polygons to annulus ...");\r
+\r
+        if (srcImage.isColorImage()) {\r
+            cf = 4;\r
+        } else {\r
+            cf = 1;\r
+        }\r
+        xDimSource = srcImage.getExtents()[0];\r
+        yDimSource = srcImage.getExtents()[1];\r
+        sourceSlice = xDimSource * yDimSource;\r
+\r
+        xDimDest = destImage.getExtents()[0];\r
+        yDimDest = destImage.getExtents()[1];\r
+        destSlice = xDimDest * yDimDest;\r
+        srcBuffer = new double[cf * sourceSlice];\r
+\r
+        try {\r
+            srcImage.exportData(0, cf * sourceSlice, srcBuffer);\r
+        } catch (IOException e) {\r
+            MipavUtil.displayError("IOException " + e + " on srcImage.exportData");\r
+\r
+            setCompleted(false);\r
+\r
+            return;\r
+        }\r
+        \r
+     // Invert so y axis increases going up so can use ccw ordering\r
+        double srcBuffer2[] = new double[cf *sourceSlice];\r
+        for (y = 0; y < yDimSource; y++) {\r
+               for (x = 0; x < xDimSource; x++) {\r
+                       if (cf == 1) {\r
+                       srcBuffer2[x + (yDimSource - 1 - y)*xDimSource] = srcBuffer[ x + y*xDimSource];\r
+                       }\r
+                       else {\r
+                           for (c = 0; c < 4; c++) {\r
+                               srcBuffer2[4*(x + (yDimSource - 1 - y)*xDimSource) + c] = srcBuffer[4*(x + y*xDimSource) + c];  \r
+                           }\r
+                       }\r
+               }\r
+        }\r
+\r
+        destBuffer = new double[cf * destSlice];\r
+\r
+        if (!srcImage.isColorImage()) {\r
+            imageMin = srcImage.getMin();\r
+\r
+            for (i = 0; i < destSlice; i++) {\r
+                destBuffer[i] = imageMin;\r
+            }\r
+        } // if (!srcImage.isColorImage())\r
+        //double w[][] = new double[xSource.length][2];\r
+        //for (i = 0; i < xSource.length; i++) {\r
+        //w[i][0] = xSource[i];\r
+        // Invert so y axis increases going up so can use ccw ordering\r
+        //w[i][1] = yDimSource-1-ySource[i];\r
+        //}\r
+        ANGLES(N[0], Z1, ALFA1, 1);\r
+        ANGLES(M[0], Z0, ALFA0, 0);\r
+        // Generate the Gauss-Jacobi weights & nodes and check the input\r
+       QINIT(M[0],N[0],ALFA0, ALFA1, NPTQ, QWORK);\r
+       CHECK(ALFA0, ALFA1, M[0], N[0], ISHAPE);\r
+       DSCSOLV(TOL, IGUESS, M[0], N[0], U, C, W0, W1, PHI0, PHI1,\r
+                               Z0, Z1, ALFA0, ALFA1, NPTQ, QWORK, ISHAPE, LINEARC);\r
+        // Output will be arranged in DSCSOLV\r
+               \r
+       // Compute data for theta-function and test the accuracy\r
+       THDATA(U);\r
+        DSCTEST(M[0],N[0],U[0],C,W0,W1,Z0,Z1,ALFA0,ALFA1,NPTQ,QWORK);\r
+        double xcenterDest = (xDimDest-1.0)/2.0;\r
+               double ycenterDest = (yDimDest-1.0)/2.0;\r
+               double maxDistance = Math.min(xcenterDest,ycenterDest);\r
+               boolean idx[] = new boolean[destSlice];\r
+               j = 0;\r
+        for (y = 0; y < yDimDest; y++) {\r
+               for (x = 0; x < xDimDest; x++) {\r
+                       index = x + xDimDest*y;\r
+                       double distX = (x-xcenterDest)/maxDistance;\r
+                       double distY = (y-ycenterDest)/maxDistance;\r
+                       double distSquared = distX*distX + distY*distY;\r
+                       double dist = Math.sqrt(distSquared);\r
+                       if ((dist >= U[0]) && (dist < 1.0)) {\r
+                               idx[index] = true;\r
+                               j++;\r
+                       }\r
+               }\r
+        } // for (y = 0; y < yDimDest; y++)\r
+        zp = new double[2];\r
+        wp = new double[j][2];\r
+        j = 0;\r
+        for (y = 0; y < yDimDest; y++) {\r
+               for (x = 0; x < xDimDest; x++) {\r
+                       index = x + xDimDest*y;\r
+                       double distX = (x-xcenterDest)/maxDistance;\r
+                       double distY = (y-ycenterDest)/maxDistance;\r
+                       if (idx[index]) {\r
+                               zp[0] = distX;\r
+                               zp[1] = distY;\r
+                               ZZ0 = ZDSC(zp,0,2,M[0],N[0],U[0],C,W0,W1,Z0,Z1,ALFA0,ALFA1,PHI0,\r
+                               PHI1,NPTQ,QWORK,1);\r
+                               wp[j][0] = ZZ0[0];\r
+                               wp[j][1] = ZZ0[1];\r
+                               j++;\r
+                       }\r
+               }\r
+        } // for (y = 0; y < yDimDest; y++)\r
+        \r
+        j = 0;\r
+        for (y = 0; y < yDimDest; y++) {\r
+               for (x = 0; x < xDimDest; x++) {\r
+                       index = x + xDimDest*y;\r
+                       if (idx[index]) {\r
+                               xp = wp[j][0];\r
+                               yp = wp[j][1];\r
+                               j++;\r
+                                       if ((xp > -0.5) && (xp < xDimSource - 0.5) && (yp > -0.5) && (yp < yDimSource - 0.5)) {\r
+                                  if (xp <= 0) {\r
+                                   x0 = 0;\r
+                                   dx = 0;\r
+                                   deltaX = 0;\r
+                               } else if (xp >= (xDimSource - 1)) {\r
+                                   x0 = xDimSource - 1;\r
+                                   dx = 0;\r
+                                   deltaX = 0;\r
+                               } else {\r
+                                   x0 = (int) xp;\r
+                                   dx = xp - x0;\r
+                                   deltaX = 1;\r
+                               }\r
+               \r
+                               if (yp <= 0) {\r
+                                   y0 = 0;\r
+                                   dy = 0;\r
+                                   deltaY = 0;\r
+                               } else if (yp >= (yDimSource - 1)) {\r
+                                   y0 = yDimSource - 1;\r
+                                   dy = 0;\r
+                                   deltaY = 0;\r
+                               } else {\r
+                                   y0 = (int) yp;\r
+                                   dy = yp - y0;\r
+                                   deltaY = xDimSource;\r
+                               }\r
+               \r
+                               dx1 = 1 - dx;\r
+                               dy1 = 1 - dy;\r
                \r
+                               position = (y0 * xDimSource) + x0;\r
+               \r
+                               if (cf == 1) {\r
+                                   destBuffer[index] = (dy1 * ( (dx1 * srcBuffer2[position]) + (dx * srcBuffer2[position + deltaX])))\r
+                                       + (dy * ( (dx1 * srcBuffer2[position + deltaY]) + (dx * srcBuffer2[position + deltaY + deltaX])));   \r
+                               }\r
+                               else {\r
+                                        for (c = 0; c < 4; c++) {\r
+                                        destBuffer[4*index+c] = (dy1 * ( (dx1 * srcBuffer2[4*position]+c) + (dx * srcBuffer2[4*(position + deltaX) + c])))\r
+                                                + (dy * ( (dx1 * srcBuffer2[4*(position + deltaY) + c]) + (dx * srcBuffer2[4*(position + deltaY + deltaX) + c])));\r
+                                    } // for (c = 0; c < 4; c++)       \r
+                               }\r
+                              } // if ((xp > -0.5) && (xp < xDimSource - 0.5) && (yp > -0.5) && (yp < yDimSource - 0.5))\r
+                       }\r
+               }\r
+               } // for (y = 0; y < yDimDest; y++)\r
+               \r
+       // Undo y axis inversion\r
+    double destBuffer2[] = new double[cf * destSlice];\r
+       for (y = 0; y < yDimDest; y++) {\r
+       for (x = 0; x < xDimDest; x++) {\r
+               if (cf == 1) {\r
+               destBuffer2[x + (yDimDest - 1 - y)*xDimDest]    = destBuffer[ x + y*xDimDest];\r
+               }\r
+               else {\r
+                   for (c = 0; c < 4; c++) {\r
+                       destBuffer2[4*(x + (yDimDest - 1 - y)*xDimDest) + c]    = destBuffer[4*(x + y*xDimDest) + c];\r
+                   }\r
+               }\r
+       }\r
+    }\r
+    \r
+\r
+       \r
+       try {\r
+           destImage.importData(0, destBuffer2, true);\r
+       } catch (IOException e) {\r
+           MipavUtil.displayError("IOException " + e + " on destImage.importData");\r
+\r
+           setCompleted(false);\r
+\r
+           return;\r
+       }\r
+\r
+       setCompleted(true);\r
+\r
+       return;\r
        }\r
        \r
        private void DSCDATA(int IPOLY, int M[], int N[], double Z0[][], double Z1[][], double ALFA0[], double ALFA1[]) {\r