plot routine added.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Tue, 12 Dec 2017 19:44:22 +0000 (19:44 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Tue, 12 Dec 2017 19:44:22 +0000 (19:44 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15308 ba61647d-9d00-f842-95cd-605cb4296b96

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

index cceff8a..2873be2 100644 (file)
@@ -1,9 +1,15 @@
 package gov.nih.mipav.model.algorithms;\r
 \r
+import java.awt.Color;\r
+import java.awt.Graphics;\r
+import java.awt.Rectangle;\r
 import java.io.IOException;\r
+import java.util.Vector;\r
 \r
 import gov.nih.mipav.model.structures.ModelImage;\r
 import gov.nih.mipav.view.MipavUtil;\r
+import gov.nih.mipav.view.ViewJComponentGraph;\r
+import gov.nih.mipav.view.ViewJFrameGraph;\r
 \r
 public class DoublyConnectedSC extends AlgorithmBase {\r
        // This is a port of the FORTRAN ACM TOMS Algorithm 785, collected\r
@@ -122,7 +128,11 @@ public class DoublyConnectedSC extends AlgorithmBase {
        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
+       private int ISHAPE = 0;\r
+       \r
+       private Vector<Double> linhx[][] = null;\r
+       private Vector<Double> linhy[][] = null;\r
+\r
        \r
        public DoublyConnectedSC() {\r
                \r
@@ -385,6 +395,8 @@ public class DoublyConnectedSC extends AlgorithmBase {
        // 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
+        plot(U[0], null, null, 200, 140, yDimSource-1, W0, W1, Z0, Z1,\r
+                       ALFA0, ALFA1, QWORK, C);\r
         double xcenterDest = (xDimDest-1.0)/2.0;\r
                double ycenterDest = (yDimDest-1.0)/2.0;\r
                double maxDistance = Math.min(xcenterDest,ycenterDest);\r
@@ -512,6 +524,423 @@ public class DoublyConnectedSC extends AlgorithmBase {
        return;\r
        }\r
        \r
+       private void plot(double U,\r
+                       double R[], double theta[], int num1draw, int num2draw, int yInvert,\r
+                       double W0[][], double W1[][], double Z0[][], double Z1[][], double ALFA0[], double ALFA1[],\r
+                       double QWORK[], double C[]) {\r
+               // Image of polar grid under Schwarz-Christoffel disk map.\r
+               // dplot will adaptively plot the images under the Schwarz-Christoffel \r
+               // disk map of circles and rays in the unit disk. \r
+               // If R.length == 1 and theta.length == 1, R[0] is the\r
+               // number of evenly spaced circles and theta[0] is the number of\r
+               // evenly spaced rays.  If R.length > 1 and theta.length > 1, then the\r
+               // circles are plotted at radii given by the entries of R and rays at\r
+               // the angles specified in theta.\r
+               \r
+               // nqpts Number of quadrature points per integration.\r
+               // Approximately equals -log10(error).  Increase if plot\r
+           // has false little zigzags in curves. \r
+               \r
+               // Original MATLAB dplot routine copyright 1998 by Toby Driscoll.\r
+               \r
+               // In original code num1draw = 20 and num2draw = 14\r
+               \r
+               int i;\r
+               int j;\r
+               int k;\r
+               int m;\r
+               int Rlength;\r
+               int thetalength;\r
+               double spacing;\r
+               double R2[];\r
+               double theta2[] = null;\r
+               int INEAR[] = new int[1];\r
+               int KNEAR[] = new int[1];\r
+               double zc[] = new double[2];\r
+               double wa[] = new double[2];\r
+               double wout[];\r
+               double cr[] = new double[1];\r
+               double ci[] = new double[1];\r
+               double posx1;\r
+               double posy1;\r
+               double posx2;\r
+               double posy2;\r
+               int x1;\r
+               int y1;\r
+               int x2;\r
+               int y2;\r
+               boolean drawTheta = true;\r
+               int n = Z0.length;\r
+               // Minimum line segment length, as a proportion of the axes box\r
+        double minlen = 0.005;\r
+        // Maximum line segment length, as a proportion of the axes box\r
+        double maxlen = 0.02;\r
+        // Max allowed number of adaptive refinements made to meet other requirements \r
+        int maxrefn = 2;\r
+               double axlim[] = new double[4];\r
+               Vector<Double>x1Vector = new Vector<Double>();\r
+               Vector<Double>y1Vector = new Vector<Double>();\r
+               Vector<Double>x2Vector = new Vector<Double>();\r
+               Vector<Double>y2Vector = new Vector<Double>();\r
+               Vector<Double>tpReal = new Vector<Double>();\r
+               Vector<Double>tpImag = new Vector<Double>();\r
+               Vector<Double>RpReal = new Vector<Double>();\r
+               Vector<Double>RpImag = new Vector<Double>();\r
+               Vector<Double>zpReal = new Vector<Double>();\r
+               Vector<Double>zpImag = new Vector<Double>();\r
+               Vector<Double>wpReal = new Vector<Double>();\r
+               Vector<Double>wpImag = new Vector<Double>();\r
+               Vector<Boolean>newlog = new Vector<Boolean>();\r
+               double zp[][];\r
+               double neww[][];\r
+               // Empty arguments default to 10.\r
+               if ((R == null) || (R.length == 0)) {\r
+                       R = new double[1];\r
+                       R[0] = 10;\r
+               }\r
+               if ((theta == null) || (theta.length == 0)) {\r
+                       theta = new double[1];\r
+                       theta[0] = 10;\r
+               }\r
+               \r
+               // Integer arguments must be converted to specific values\r
+               if ((R.length == 1) && (R[0] == Math.round(R[0]))) {\r
+                       Rlength = (int)R[0];\r
+                       R2 = new double[Rlength];\r
+                   spacing = (1.0-U)/(R[0] + 1.0);\r
+                   for (i = 1; i <= Rlength; i++) {\r
+                       R2[i-1] = U + i*spacing;        \r
+                   }\r
+               } // if ((R.length == 1) && (R[0] == Math.round(R[0])))\r
+               else {\r
+                       R2 = new double[R.length];\r
+                       for (i = 0; i < R.length; i++) {\r
+                               R2[i] = R[i];\r
+                       }\r
+               }\r
+               if ((theta.length == 1) && (theta[0] == 0.0)) {\r
+                       drawTheta = false;\r
+               }\r
+               else if ((theta.length == 1) && (theta[0] == Math.round(theta[0]))) {\r
+                       thetalength = (int)theta[0];\r
+                       theta2 = new double[thetalength];\r
+                   spacing = ((2.0*Math.PI)/theta[0]);\r
+                   for (i = 0; i < thetalength; i++) {\r
+                       theta2[i] = i*spacing;  \r
+                   }\r
+               } // if ((theta.length == 1) && (theta[0] == Math.round(theta[0])))\r
+               else {\r
+                       theta2 = new double[theta.length];\r
+                       for (i = 0; i < theta.length; i++) {\r
+                               theta2[i] = theta[i];\r
+                       }\r
+               }\r
+               int numinfinite = 0;\r
+               // Put 2 finite points at every infinity in plotpoly\r
+               for (i = 0; i < Z0.length; i++) {\r
+                   if (Double.isInfinite(Z0[i][0]) || Double.isInfinite(Z0[i][1]))     {\r
+                       numinfinite++;  \r
+                   }\r
+               }\r
+               float xPointArray[] = new float[n+1+numinfinite];\r
+               float yPointArray[] = new float[n+1+numinfinite];\r
+               double beta[] = new double[ALFA0.length];\r
+               for (i = 0; i < ALFA0.length; i++) {\r
+                       beta[i] = ALFA0[i] - 1.0;\r
+               }\r
+               ViewJFrameGraph pointGraph =scm. plotpoly(xPointArray, yPointArray, Z0, beta, \r
+                               false, axlim, yInvert, true, null);\r
+               ViewJComponentGraph graph = pointGraph.getGraph();\r
+               Rectangle graphBounds = graph.getGraphBounds();\r
+               Graphics g = graph.getGraphics();\r
+               double xScale = graphBounds.width / (axlim[1] - axlim[0]);\r
+        double yScale = graphBounds.height / (axlim[3] - axlim[2]);\r
+        double len = Math.max(axlim[1] - axlim[0], axlim[3] - axlim[2]);\r
+               minlen = len * minlen;\r
+               maxlen = len * maxlen;\r
+               \r
+               // Plot circles...\r
+           linhx = new Vector[R2.length][2];\r
+               linhy = new Vector[R2.length][2];\r
+               for (i = 0; i < R2.length; i++) {\r
+                       for (j = 0; j < 2; j++) {\r
+                               linhx[i][j] = new Vector<Double>();\r
+                               linhy[i][j] = new Vector<Double>();\r
+                       }\r
+               }\r
+               for (j = 0; j < R2.length; j++) {\r
+                   // Start with evenly spaced theta\r
+                       tpReal.clear();\r
+                       tpImag.clear();\r
+                       for (i = 0; i < num1draw; i++) {\r
+                               tpReal.add(i*(2.0*Math.PI)/(num1draw-1.0));\r
+                               tpImag.add(0.0);\r
+                       }\r
+                       newlog.clear();\r
+                       for (i = 0; i < num1draw; i++) {\r
+                               newlog.add(true);\r
+                       }\r
+                       wpReal.clear();\r
+                       wpImag.clear();\r
+                       for (i = 0; i < num1draw; i++) {\r
+                               wpReal.add(Double.NaN);\r
+                               wpImag.add(0.0);\r
+                       }\r
+                       \r
+                       // The individual points will be shown as they are found.\r
+                       \r
+                       // Adaptive refinement to make curve smooth\r
+                       int iter = 0;\r
+                       while (iter < maxrefn) {\r
+                               int numnewlog = 0;\r
+                               for (i = 0; i < newlog.size(); i++) {\r
+                                   if (newlog.get(i)) {\r
+                                       numnewlog++;\r
+                                   }\r
+                               } // for (i = 0; i < newlog.size(); i++)\r
+                               if (numnewlog == 0) {\r
+                                       break;\r
+                               }\r
+                               zp = new double[numnewlog][2];\r
+                               for (i = 0, m = 0; i < newlog.size(); i++) {\r
+                                   if (newlog.get(i)) {\r
+                                       zp[m][0] = R2[j]*Math.cos(tpReal.get(i));\r
+                                       zp[m++][1] = R2[j]*Math.sin(tpReal.get(i));\r
+                                   }\r
+                               } // for (i = 0, m = 0; i < newlog.length; i++)\r
+                               neww = new double[numnewlog][2];\r
+                               for (k = 0; k < numnewlog; k++) {\r
+                                       NEARW(M[0],N[0],W0,W1,ALFA0,zp[k],KNEAR,INEAR); \r
+                                       if (INEAR[0] == 0) {\r
+                                               zc[0] = Z0[KNEAR[0]-1][0];\r
+                                               zc[1] = Z0[KNEAR[0]-1][1];\r
+                                               wa[0] = W0[KNEAR[0]-1][0];\r
+                                               wa[1] = W0[KNEAR[0]-1][1];\r
+                                       }\r
+                                       else {\r
+                                               zc[0] = Z1[KNEAR[0]-1][0];\r
+                                               zc[1] = Z1[KNEAR[0]-1][1];\r
+                                               wa[0] = W1[KNEAR[0]-1][0];\r
+                                               wa[1] = W1[KNEAR[0]-1][1];      \r
+                                       }\r
+                                       wout = WQUAD1(wa, 0, KNEAR[0], INEAR[0], zp[k], 0, 0, M[0], N[0], U, W0, W1,ALFA0, ALFA1, NPTQ, QWORK, 0);\r
+                                       scm.zmlt(C[0], C[1], wout[0], wout[1], cr, ci);\r
+                                       neww[k][0] = zc[0] + cr[0];\r
+                                       neww[k][1] = zc[1] + ci[0];\r
+                               } // for (k = 0; k < numnewlog; k++)\r
+                               for (i = 0, m = 0; i < newlog.size(); i++) {\r
+                                       if (newlog.get(i)) {\r
+                                           wpReal.set(i, neww[m][0]);\r
+                                           wpImag.set(i, neww[m++][1]);\r
+                                       } \r
+                               } // for (i = 0, m = 0; i < newlog.size(); i++)\r
+                               iter = iter + 1;\r
+                               \r
+                               linhx[j][0].clear();\r
+                               linhy[j][0].clear();\r
+                               linhx[j][1].clear();\r
+                               linhy[j][1].clear();\r
+                               // Update the points to show progress\r
+                               for (i = 0; i < wpReal.size(); i++) {\r
+                                   linhx[j][0].add(wpReal.get(i));\r
+                                   linhy[j][0].add(wpImag.get(i));\r
+                               }\r
+                               for (i = 0; i < zp.length; i++) {\r
+                                       linhx[j][1].add(zp[i][0]);\r
+                                   linhy[j][1].add(zp[i][1]);  \r
+                               }\r
+                               \r
+                               // Add points to tp where necessary\r
+                               scm.scpadapt(tpReal, tpImag, wpReal, wpImag, newlog, minlen, maxlen, axlim);\r
+                       } // while (iter < maxrefn)\r
+               } // for (j = 0; j < R2.length; j++)\r
+               \r
+               for (i = 0; i < Z1.length-1; i++) {\r
+                   posx1 = Z1[i][0];\r
+                   posy1 = Z1[i][1];\r
+                   posx2 = Z1[i+1][0];\r
+                   posy2 = Z1[i+1][1];\r
+                   x1Vector.add(posx1);\r
+               y1Vector.add(posy1);\r
+               x2Vector.add(posx2);\r
+               y2Vector.add(posy2);\r
+           x1 =  (int)Math.round(graphBounds.x + xScale*(posx1 - axlim[0]));\r
+                   y1 =  (int)Math.round(graphBounds.y + yScale*(posy1 - axlim[2]));\r
+                   y1 = -y1 + 2*graphBounds.y + graphBounds.height;\r
+                   x2 =  (int)Math.round(graphBounds.x + xScale*(posx2 - axlim[0]));\r
+                   y2 =  (int)Math.round(graphBounds.y + yScale*(posy2 - axlim[2]));\r
+                   y2 = -y2 + 2*graphBounds.y + graphBounds.height;\r
+                   graph.drawLine(g, x1, y1, x2, y2);\r
+               }\r
+               posx1 = Z1[Z1.length-1][0];\r
+           posy1 = Z1[Z1.length-1][1];\r
+           posx2 = Z1[0][0];\r
+           posy2 = Z1[0][1];\r
+           x1Vector.add(posx1);\r
+               y1Vector.add(posy1);\r
+               x2Vector.add(posx2);\r
+               y2Vector.add(posy2);\r
+           x1 =  (int)Math.round(graphBounds.x + xScale*(posx1 - axlim[0]));\r
+           y1 =  (int)Math.round(graphBounds.y + yScale*(posy1 - axlim[2]));\r
+           y1 = -y1 + 2*graphBounds.y + graphBounds.height;\r
+           x2 =  (int)Math.round(graphBounds.x + xScale*(posx2 - axlim[0]));\r
+           y2 =  (int)Math.round(graphBounds.y + yScale*(posy2 - axlim[2]));\r
+           y2 = -y2 + 2*graphBounds.y + graphBounds.height;\r
+           graph.drawLine(g, x1, y1, x2, y2);\r
+               \r
+               for (i = 0; i < R2.length; i++) {\r
+                   if ((linhx[i][0] != null)&& (linhy[i][0] != null)) {\r
+                       for (j = 0; j < linhx[i][0].size()-1; j++) {\r
+                               posx1 = linhx[i][0].get(j);\r
+                               posy1 = linhy[i][0].get(j);\r
+                               posx2 = linhx[i][0].get(j+1);\r
+                               posy2 = linhy[i][0].get(j+1);\r
+                               x1Vector.add(posx1);\r
+                               y1Vector.add(posy1);\r
+                               x2Vector.add(posx2);\r
+                               y2Vector.add(posy2);\r
+                           x1 =  (int)Math.round(graphBounds.x + xScale*(posx1 - axlim[0]));\r
+                           y1 =  (int)Math.round(graphBounds.y + yScale*(posy1 - axlim[2]));\r
+                           y1 = -y1 + 2*graphBounds.y + graphBounds.height;\r
+                           x2 =  (int)Math.round(graphBounds.x + xScale*(posx2 - axlim[0]));\r
+                           y2 =  (int)Math.round(graphBounds.y + yScale*(posy2 - axlim[2]));\r
+                           y2 = -y2 + 2*graphBounds.y + graphBounds.height;\r
+                           graph.drawLine(g, x1, y1, x2, y2);\r
+                       }\r
+                   }\r
+               } // for (i = 0; i < R2.length; i++)\r
+               \r
+               // Plot radii\r
+               if (drawTheta) {\r
+                       linhx = new Vector[theta2.length][2];\r
+                       linhy = new Vector[theta2.length][2];\r
+                       for (i = 0; i < theta2.length; i++) {\r
+                               for (j = 0; j < 2; j++) {\r
+                                       linhx[i][j] = new Vector<Double>();\r
+                                       linhy[i][j] = new Vector<Double>();\r
+                               }\r
+                       }\r
+                       for (j = 0; j < theta2.length; j++) {\r
+                               RpReal.clear();\r
+                               RpImag.clear();\r
+                               zpReal.clear();\r
+                               zpImag.clear();\r
+                               for (i = 0; i < num2draw; i++) {\r
+                                       RpReal.add(U + (1-U)*i/(num2draw-1.0));\r
+                                       RpImag.add(0.0);\r
+                                       zpReal.add(RpReal.get(i)*Math.cos(theta2[j]));\r
+                                       zpImag.add(RpReal.get(i)*Math.sin(theta2[j]));\r
+                               }\r
+                               newlog.clear();\r
+                               for (i = 0; i < num2draw; i++) {\r
+                                       newlog.add(true);\r
+                               }\r
+                               wpReal.clear();\r
+                               wpImag.clear();\r
+                               for (i = 0; i < num2draw; i++) {\r
+                                       wpReal.add(Double.NaN);\r
+                                       wpImag.add(0.0);\r
+                               }\r
+                               \r
+                               // The individual points will be shown as they are found.\r
+                               \r
+                               // Adaptive refinement to make curve smooth\r
+                               int iter = 0;\r
+                               while (iter < maxrefn) {\r
+                                       int numnewlog = 0;\r
+                                       for (i = 0; i < newlog.size(); i++) {\r
+                                           if (newlog.get(i)) {\r
+                                               numnewlog++;\r
+                                           }\r
+                                       } // for (i = 0; i < newlog.size(); i++)\r
+                                       if (numnewlog == 0) {\r
+                                               break;\r
+                                       }\r
+                                       double zpnew[][] = new double[numnewlog][2];\r
+                                       for (i = 0, m = 0; i < newlog.size(); i++) {\r
+                                           if (newlog.get(i)) {\r
+                                               zpnew[m][0] = zpReal.get(i);\r
+                                               zpnew[m++][1] = zpImag.get(i);\r
+                                           }\r
+                                       } // for (i = 0, m = 0; i < newlog.length; i++)\r
+                                       neww = new double[numnewlog][2];\r
+                                       for (k = 0; k < numnewlog; k++) {\r
+                                               NEARW(M[0],N[0],W0,W1,ALFA0,zpnew[k],KNEAR,INEAR);      \r
+                                               if (INEAR[0] == 0) {\r
+                                                       zc[0] = Z0[KNEAR[0]-1][0];\r
+                                                       zc[1] = Z0[KNEAR[0]-1][1];\r
+                                                       wa[0] = W0[KNEAR[0]-1][0];\r
+                                                       wa[1] = W0[KNEAR[0]-1][1];\r
+                                               }\r
+                                               else {\r
+                                                       zc[0] = Z1[KNEAR[0]-1][0];\r
+                                                       zc[1] = Z1[KNEAR[0]-1][1];\r
+                                                       wa[0] = W1[KNEAR[0]-1][0];\r
+                                                       wa[1] = W1[KNEAR[0]-1][1];      \r
+                                               }\r
+                                               wout = WQUAD1(wa, 0, KNEAR[0], INEAR[0], zpnew[k], 0, 0, M[0], N[0], U, W0, W1,ALFA0, ALFA1, NPTQ, QWORK, 0);\r
+                                               scm.zmlt(C[0], C[1], wout[0], wout[1], cr, ci);\r
+                                               neww[k][0] = zc[0] + cr[0];\r
+                                               neww[k][1] = zc[1] + ci[0];\r
+                                       } // for (k = 0; k < numnewlog; k++)\r
+                                       for (i = 0, m = 0; i < newlog.size(); i++) {\r
+                                               if (newlog.get(i)) {\r
+                                                   wpReal.set(i, neww[m][0]);\r
+                                                   wpImag.set(i, neww[m++][1]);\r
+                                               } \r
+                                       } // for (i = 0, m = 0; i < newlog.size(); i++)\r
+                                       iter = iter + 1;\r
+                                       \r
+                                       linhx[j][0].clear();\r
+                                       linhy[j][0].clear();\r
+                                       linhx[j][1].clear();\r
+                                       linhy[j][1].clear();\r
+                                       // Update the points to show progress\r
+                                       for (i = 0; i < wpReal.size(); i++) {\r
+                                           linhx[j][0].add(wpReal.get(i));\r
+                                           linhy[j][0].add(wpImag.get(i));\r
+                                           linhx[j][1].add(zpReal.get(i));\r
+                                           linhy[j][1].add(zpImag.get(i));\r
+                                       }\r
+                                       \r
+                                       scm.scpadapt(zpReal, zpImag, wpReal, wpImag, newlog, minlen, maxlen, axlim);\r
+                               } // while (iter < maxrefn)\r
+                       } // for (j = 0; j < theta2.length; j++)\r
+                       \r
+                       for (i = 0; i < theta2.length; i++) {\r
+                           if ((linhx[i][0] != null)&& (linhy[i][0] != null)) {\r
+                               for (j = 0; j < linhx[i][0].size()-1; j++) {\r
+                                       posx1 = linhx[i][0].get(j);     \r
+                                       posy1 = linhy[i][0].get(j);\r
+                                       posx2 = linhx[i][0].get(j+1);\r
+                                       posy2 = linhy[i][0].get(j+1);\r
+                                       x1Vector.add(posx1);\r
+                                       y1Vector.add(posy1);\r
+                                       x2Vector.add(posx2);\r
+                                       y2Vector.add(posy2);\r
+                                   x1 =  (int)Math.round(graphBounds.x + xScale*(posx1 - axlim[0]));\r
+                                   y1 =  (int)Math.round(graphBounds.y + yScale*(posy1 - axlim[2]));\r
+                                   y1 = -y1 + 2*graphBounds.y + graphBounds.height;\r
+                                   x2 =  (int)Math.round(graphBounds.x + xScale*(posx2 - axlim[0]));\r
+                                   y2 =  (int)Math.round(graphBounds.y + yScale*(posy2 - axlim[2]));\r
+                                   y2 = -y2 + 2*graphBounds.y + graphBounds.height;\r
+                                   graph.drawLine(g, x1, y1, x2, y2);\r
+                               }\r
+                           }\r
+                       } // for (i = 0; i < theta2.length; i++)\r
+               } // if (drawTheta)\r
+               \r
+               graph.setX1Vector(x1Vector);\r
+               graph.setY1Vector(y1Vector);\r
+               graph.setX2Vector(x2Vector);\r
+               graph.setY2Vector(y2Vector);\r
+               graph.setAddSchwarzChristoffelLines(true);\r
+               graph.paintComponent(g);\r
+               return;\r
+       }\r
+       \r
+       \r
+       \r
        private void DSCDATA(int IPOLY, int M[], int N[], double Z0[][], double Z1[][], double ALFA0[], double ALFA1[]) {\r
                // DSCDATA generates data.\r
                double Q;\r