package ij.plugin.filter;
import ij.*;
import ij.plugin.*;
import ij.process.*;
import ij.gui.*;

/**
 * This plugin implements the Euclidean Distance Map (EDM), Ultimate Eroded Points and
 * Watershed commands in the Process/Binary submenu.
 * Note: These functions do not take ROIs into account (any ROI gets deselected).
 * setup is called with no argument for EDM, "points" for ultimate eroded points and "watershed"
 * for watershed segmentation.
 * Ultimate Eroded Points and Watershed are handled by the MaximumFinder 
 * plugin applied to the EDM
 * 
 * version 09-Nov-2006 Michael Schmid
 */
	public class EDM implements PlugInFilter {
    /** unit in 16-bit EDM image: this value corresponds to a distance of one pixel */
    public static final int ONE = 41;
    /** in 16-bit EDM image this value corresponds to a pixel distance of sqrt(2) */
    public static final int SQRT2 = 58; // ~ 41 * sqrt(2)
    /** in 16-bit EDM image this value corresponds to a pixel distance of sqrt(2) */
    public static final int SQRT5 = 92; // ~ 41 * sqrt(5)

    ImagePlus imp;
    String arg;
    int slice;                                  // when processing a stack
    boolean invertImage;
    
    public int setup(String arg, ImagePlus imp) {
        this.imp = imp;
        this.arg = arg;
        if (imp != null) {
            boolean invertedLut = imp.isInvertedLut();
            invertImage = (invertedLut && Prefs.blackBackground) || (!invertedLut && !Prefs.blackBackground);
        }
        return IJ.setupDialog(imp, DOES_8G);
    }
    
    public void run(ImageProcessor ip) {
        imp.killRoi();                          //otherwise invert won't work as expected
        int[] histogram;
        slice++;
        ImageStatistics stats = imp.getStatistics();
        if (slice==1 && stats.histogram[0]+stats.histogram[255]!=stats.pixelCount) {
            IJ.error("8-bit binary image (0 and 255) required.");
            return;
        }

        if (invertImage)
            ip.invert();                            // independent of settings, now 255=foreground, 0=background
        if (arg.equals("points")||arg.equals("watershed")) {
            int outputType;
            if (arg.equals("watershed")) {
                outputType = MaximumFinder.SEGMENTED;
            } else {                                // output ultimate points: one per maximum
                outputType = MaximumFinder.SINGLE_POINTS;
            }
            ImageProcessor ip16 = make16bitEDM(ip);
            //new ImagePlus("16-bit EDM",ip16).show();
            MaximumFinder fm = new MaximumFinder();

            //setting the tolerance to lower values such as 0.3*ONE creates more segments,
            //larger values like 0.7*ONE give fewer segments
            //reasonable values are 0.3 ... 0.8 * ONE.
            ByteProcessor maxIp = fm.findMaxima(ip16, 0.5*ONE, ImageProcessor.NO_THRESHOLD, outputType, false, true);
            if (maxIp != null) {                    //no further action if cancelled by user
                if (arg.equals("watershed")) {
                    ip.copyBits(maxIp, 0, 0, Blitter.COPY);
                    
                } else {                            // output ultimate points
                    byte[] image8 = (byte[])ip.getPixels();
                    short[] image16 = (short[])ip16.getPixels();
                    convertToBytes(image16, image8);
                    ip.copyBits(maxIp, 0, 0, Blitter.AND); // maxima are 255 and keep the EDM value
                }
            }
        }
        else {                                      // if output should be EDM
            toEDM(ip);
        }
        if (invertImage)
            ip.invert();                            // back to black foreground (unless settings say reverse)
    } //public void run
    
    
    
    /**	Converts a binary image into a 8-bit grayscale Euclidean Distance Map
     * (EDM). Each foreground (zero) pixel in the binary image is
     * assigned a value equal to its distance from the nearest
     * background (255) pixel.
     */
    public void toEDM(ImageProcessor ip) {
        ImageProcessor ip16 = make16bitEDM(ip);
        byte[] image8 = (byte[])ip.getPixels();
        short[] image16 = (short[])ip16.getPixels();
        convertToBytes(image16, image8);
    }
        
    /**	Calculates a 16-bit grayscale Euclidean Distance Map for a binary 8-bit image.
     * Each foreground (black) pixel in the binary image is
     * assigned a value equal to its distance from the nearest
     * background (white) pixel.  Uses the two-pass EDM algorithm
     * from the "Image Processing Handbook" by John Russ.
     */
    public ShortProcessor make16bitEDM(ImageProcessor ip) {
        int xmax, ymax;
        int offset, rowsize;
        
        IJ.showStatus("Generating EDM");
        int width = imp.getWidth();
        int height = imp.getHeight();
        rowsize = width;
        xmax    = width - 2;
        ymax    = height - 2;
        ShortProcessor ip16 = (ShortProcessor)ip.convertToShort(false);
        ip16.multiply(128); //foreground pixels set (almost) as high as possible for signed short
        short[] image16 = (short[])ip16.getPixels();
        int inc = Math.max(height/50,1);
        
        for (int y=0; y<height; y++) {
            for (int x=0; x<width; x++) {
                offset = x + y * rowsize;
                if (image16[offset] > 0) {
                    if ((x<=1) || (x>=xmax) || (y<=1) || (y>=ymax))
                        setEdgeValue(offset, rowsize, image16, x, y, xmax, ymax);
                    else
                        setValue(offset, rowsize, image16);
                }
            } // for x
            if ((inc&0)==0) IJ.showProgress(y/2, height);
        } // for y
        
        for (int y=height-1; y>=0; y--) {
            for (int x=width-1; x>=0; x--) {
                offset = x + y * rowsize;
                if (image16[offset] > 0) {
                    if ((x<=1) || (x>=xmax) || (y<=1) || (y>=ymax))
                        setEdgeValue(offset, rowsize, image16, x, y, xmax, ymax);
                    else
                        setValue(offset, rowsize, image16);
                }
            } // for x
            if ((inc&0)==0) IJ.showProgress(height/2+(height-y)/2, height);
        } // for y
        IJ.showProgress(1,1);
        //(new ImagePlus("EDM16", ip16.duplicate())).show();
        
        return ip16;
    } // make16bitEDM(ip)
    
    void setValue(int offset, int rowsize, short[] image16) {
        int  v;
        int r1  = offset - rowsize - rowsize - 2;
        int r2  = r1 + rowsize;
        int r3  = r2 + rowsize;
        int r4  = r3 + rowsize;
        int r5  = r4 + rowsize;
        int min = 32767;
        
        v = image16[r2 + 2] + ONE;
        if (v < min)
            min = v;
        v = image16[r3 + 1] + ONE;
        if (v < min)
            min = v;
        v = image16[r3 + 3] + ONE;
        if (v < min)
            min = v;
        v = image16[r4 + 2] + ONE;
        if (v < min)
            min = v;
        
        v = image16[r2 + 1] + SQRT2;
        if (v < min)
            min = v;
        v = image16[r2 + 3] + SQRT2;
        if (v < min)
            min = v;
        v = image16[r4 + 1] + SQRT2;
        if (v < min)
            min = v;
        v = image16[r4 + 3] + SQRT2;
        if (v < min)
            min = v;
        
        v = image16[r1 + 1] + SQRT5;
        if (v < min)
            min = v;
        v = image16[r1 + 3] + SQRT5;
        if (v < min)
            min = v;
        v = image16[r2 + 4] + SQRT5;
        if (v < min)
            min = v;
        v = image16[r4 + 4] + SQRT5;
        if (v < min)
            min = v;
        v = image16[r5 + 3] + SQRT5;
        if (v < min)
            min = v;
        v = image16[r5 + 1] + SQRT5;
        if (v < min)
            min = v;
        v = image16[r4] + SQRT5;
        if (v < min)
            min = v;
        v = image16[r2] + SQRT5;
        if (v < min)
            min = v;
        
        image16[offset] = (short)min;
        
    } // setValue()
    
    void setEdgeValue(int offset, int rowsize, short[] image16, int x, int y, int xmax, int ymax) {
        int  v;
        int r1 = offset - rowsize - rowsize - 2;
        int r2 = r1 + rowsize;
        int r3 = r2 + rowsize;
        int r4 = r3 + rowsize;
        int r5 = r4 + rowsize;
        int min = 32767;
        int offimage = image16[r3 + 2];
        
        if (y<1)
            v = offimage + ONE;
        else
            v = image16[r2 + 2] + ONE;
        if (v < min)
            min = v;
        
        if (x<1)
            v = offimage + ONE;
        else
            v = image16[r3 + 1] + ONE;
        if (v < min)
            min = v;
        
        if (x>xmax)
            v = offimage + ONE;
        else
            v = image16[r3 + 3] + ONE;
        if (v < min)
            min = v;
        
        if (y>ymax)
            v = offimage + ONE;
        else
            v = image16[r4 + 2] + ONE;
        if (v < min)
            min = v;
        
        if ((x<1) || (y<1))
            v = offimage + SQRT2;
        else
            v = image16[r2 + 1] + SQRT2;
        if (v < min)
            min = v;
        
        if ((x>xmax) || (y<1))
            v = offimage + SQRT2;
        else
            v = image16[r2 + 3] + SQRT2;
        if (v < min)
            min = v;
        
        if ((x<1) || (y>ymax))
            v = offimage + SQRT2;
        else
            v = image16[r4 + 1] + SQRT2;
        if (v < min)
            min = v;
        
        if ((x>xmax) || (y>ymax))
            v = offimage + SQRT2;
        else
            v = image16[r4 + 3] + SQRT2;
        if (v < min)
            min = v;
        
        if ((x<1) || (y<=1))
            v = offimage + SQRT5;
        else
            v = image16[r1 + 1] + SQRT5;
        if (v < min)
            min = v;
        
        if ((x>xmax) || (y<=1))
            v = offimage + SQRT5;
        else
            v = image16[r1 + 3] + SQRT5;
        if (v < min)
            min = v;
        
        if ((x>=xmax) || (y<1))
            v = offimage + SQRT5;
        else
            v = image16[r2 + 4] + SQRT5;
        if (v < min)
            min = v;
        
        if ((x>=xmax) || (y>ymax))
            v = offimage + SQRT5;
        else
            v = image16[r4 + 4] + SQRT5;
        if (v < min)
            min = v;
        
        if ((x>xmax) || (y>=ymax))
            v = offimage + SQRT5;
        else
            v = image16[r5 + 3] + SQRT5;
        if (v < min)
            min = v;
        
        if ((x<1) || (y>=ymax))
            v = offimage + SQRT5;
        else
            v = image16[r5 + 1] + SQRT5;
        if (v < min)
            min = v;
        
        if ((x<=1) || (y>ymax))
            v = offimage + SQRT5;
        else
            v = image16[r4] + SQRT5;
        if (v < min)
            min = v;
        
        if ((x<=1) || (y<1))
            v = offimage + SQRT5;
        else
            v = image16[r2] + SQRT5;
        if (v < min)
            min = v;
        
        image16[offset] = (short)min;
        
    } // setEdgeValue()
    
    /**convert 16-bit EDM to 8 bits*/
    void convertToBytes(short[] image16, byte[] image8) {
        int v;
        int round = ONE / 2;

        for (int i=0; i<image16.length; i++) {
                v = (image16[i] + round) / ONE;
                if (v > 255)
                    v = 255;
                image8[i] = (byte) v;
        } // for i
    } // end ConvertToBytes()
}
