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);
    }
}