Merge remote-tracking branch 'origin/master'

README.md

@@ -7,16 +7,18 @@ After the manuscript is accepted the software will be released free and open
 source.
 
 USAGE
-run src/MATLAB/demoPR.m to see the generic pixel replication usage
+=====
 
-src/ImageJ - the ImageJ plugin
++src/MATLAB/demoPR.m -- generic 2-D pixel replication example
 	
-src/Python - Python implementation
++src/MATLAB/demoPR_3D.m -- includes 3-D ellipse generator. Uses same PixelReplicate.m as 2-D example, only the simulated ellipse generation is different.
 
-src/MATLAB/3D - a 3D version
++src/ImageJ -- ImageJ plugin.
 
-src/MATLAB/DistributionSampling - a more efficient implementation (useful
-    for larger connected components
++src/Python -- Python implementation.
 
-src/PerformanceCharacterization - quantitive comparison of different algorithms on simulated 
-    elllipse data
\ No newline at end of file
++src/MATLAB/DistributionSampling -- a more efficient implementation (useful
+    for larger connected components.
+    
++src/PerformanceCharacterization -- quantitive comparison of different algorithms on simulated 
+    elllipse data.
\ No newline at end of file

src/ImageJ/Demo_Pixel_Replication.java

+import ij.*;
+import ij.plugin.filter.PlugInFilter;
+import ij.process.*;
+import ij.gui.*;
+import java.awt.*;
+import java.awt.image.*;
+import java.nio.ByteBuffer;
+import java.nio.FloatBuffer;
+import java.util.ArrayList;
+
+import static org.bytedeco.javacpp.opencv_core.*;
+import static org.bytedeco.javacpp.opencv_imgproc.*;
+import org.bytedeco.javacpp.indexer.*;
+import org.bytedeco.javacpp.indexer.DoubleIndexer;
+import org.bytedeco.javacpp.opencv_ml.*;
+
+public class Demo_Pixel_Replication implements PlugInFilter {
+
+    int k;
+
+    public int setup(String arg, ImagePlus imp) {
+        if (IJ.versionLessThan("1.37j"))
+            return DONE;
+        else    
+            return DOES_ALL+DOES_STACKS+SUPPORTS_MASKING;
+    }
+
+    public boolean showDialog() {
+        String[] kChoices = {"2", "3", "4", "5", "6"};
+
+        GenericDialog gd = new GenericDialog("Number of ellipses");
+        gd.addChoice("Number of ellipses:", kChoices, kChoices[0]);
+        gd.showDialog();
+        if (gd.wasCanceled()) {
+            return false;
+        } else {
+            k = gd.getNextChoiceIndex() + 2;
+            return true;
+        }
+    }
+
+    public void run(ImageProcessor ip) {
+        if (!showDialog())
+            return;
+
+        IplImage ipl = ip2ipl(ip);
+        Rectangle r = ip.getRoi();
+
+        // foreground points in ip
+        ArrayList<Point2f> UniqPts = new ArrayList<>();
+        // pixel-replicated points in ip
+        ArrayList<Point2f> PrPts = new ArrayList<>();
+
+        // use distance transform to populate UniqPts and PrPts
+        CvMat distMat = DistTransform(ipl);
+        cvReleaseImage(ipl);
+        FloatIndexer distMatIdx = distMat.createIndexer();
+
+        for (int y=r.y; y<(r.y+r.height); y++) {
+            for (int x=r.x; x<(r.x+r.width); x++) {
+                int n = (int) Math.round(distMatIdx.get(y,x));
+                if (n > 0) {
+                    Point2f p = new Point2f(x, y);
+                    UniqPts.add(p);
+                    for (int i = 0; i < n; i++) {
+                        PrPts.add(p);
+                    }
+                }
+            }
+        }
+
+        // compute the gmm
+        EM em = computeGMM(PrPts, k);
+
+        // create a new image showing the partition
+        ImageProcessor ipOutPr = ip.duplicate();
+        Mat PrMat = new Mat(1, 2, CV_32FC1);
+        FloatIndexer PrMatIndexer = PrMat.createIndexer();
+        int clrOffset = 256 / k;
+        for (Point2f p : UniqPts) {
+            PrMatIndexer.put(0, 0, p.x());
+            PrMatIndexer.put(0, 1, p.y());
+            Point2d pr = em.predict(PrMat);
+            int val = ((int) pr.y() + 1) * clrOffset - 1;
+            ipOutPr.set((int) p.x(), (int) p.y(), val);
+        }
+
+        ImagePlus segIP = new ImagePlus("Segmentation", ipOutPr);
+        segIP.show();
+
+        // draw ellipses around the regions found by the gmms
+        Mat means = em.getMat("means");
+        MatVector covs = em.getMatVector("covs");
+        Overlay ellipses = new Overlay();
+        for (int i = 0; i < k; i++) {
+            // compute eigenvalues and eigenvectors of the covariance matrix
+            Mat eigVal = new Mat();
+            Mat eigVec = new Mat();
+            int lo, hi;
+            eigen(covs.get(i), eigVal, eigVec);
+            DoubleIndexer valIdx = eigVal.createIndexer();
+            if (valIdx.get(0,0) > valIdx.get(1,0)) {
+                hi = 0;
+                lo = 1;
+            } else {
+                hi = 1;
+                lo = 0;
+            }
+
+            double A = Math.sqrt(valIdx.get(hi,0) * 20 / 3);
+            double B = Math.sqrt(valIdx.get(lo,0) * 20 / 3);
+            // double aspectRatio = B / A;
+
+            EllipseRoi elRoi = makeEllipseRoi(means.row(i), eigVec.row(hi), A, B);
+            // Mat end1 = new Mat();
+            // Mat end2 = new Mat();
+            // scaleAdd(eigVec.row(hi), A, means.row(i), end1);
+            // scaleAdd(eigVec.row(hi), -A, means.row(i), end2);
+            // DoubleIndexer end1Idx = end1.createIndexer();
+            // DoubleIndexer end2Idx = end2.createIndexer();
+            // EllipseRoi elRoi = new EllipseRoi(end1Idx.get(0,0), end1Idx.get(0,1), end2Idx.get(0,0), end2Idx.get(0,1), aspectRatio);
+            elRoi.setStrokeWidth(2);
+            elRoi.setStrokeColor(Color.red);
+            ellipses.add(elRoi);
+        }
+        segIP.setOverlay(ellipses);
+    }
+
+    private IplImage ip2ipl(ImageProcessor src) {
+
+        BufferedImage bi = src.getBufferedImage();
+        return IplImage.createFrom(bi);
+    }
+
+    private CvMat DistTransform(IplImage iplIn) {
+        IplImage iplOut = cvCreateImage(cvGetSize(iplIn), IPL_DEPTH_32F, 1);
+        cvDistTransform(iplIn, iplOut);
+        return iplOut.asCvMat();
+    }
+
+    private Mat PrPts2Mat(ArrayList<Point2f> PrPts) {
+        Mat emMat = new Mat(PrPts.size(), 2, CV_32FC1);
+        FloatIndexer emMatIndexer = emMat.createIndexer();
+        for (int i = 0; i < PrPts.size(); i++) {
+            Point2f p = PrPts.get(i);
+            emMatIndexer.put(i, 0, p.x());
+            emMatIndexer.put(i, 1, p.y());
+        }
+
+        return emMat;
+    }
+
+    private EM computeGMM(ArrayList<Point2f> PrPts, int k) {
+        EM em = new EM();
+        em.set("nclusters", k);
+        em.set("covMatType", EM.COV_MAT_GENERIC);
+
+        Mat emMat = PrPts2Mat(PrPts);
+        em.train(emMat);
+
+        return em;
+    }
+
+    private EllipseRoi makeEllipseRoi(Mat center, Mat unitVec, double A, double B) {
+        Mat end1 = new Mat();
+        Mat end2 = new Mat();
+
+        scaleAdd(unitVec, A, center, end1);
+        scaleAdd(unitVec, -A, center, end2);
+        double aspectRatio = B / A;
+
+        DoubleIndexer end1Idx = end1.createIndexer();
+        DoubleIndexer end2Idx = end2.createIndexer();
+
+        return new EllipseRoi(end1Idx.get(0,0), end1Idx.get(0,1), end2Idx.get(0,0), end2Idx.get(0,1), aspectRatio);
+    }
+}