| MaximumFinder |

    package ij.plugin.filter;
    import ij.*;
    import ij.gui.*;
    import ij.measure.*;
    import ij.process.*;
    import java.awt.*;
    import java.util.*;

    /** This ImageJ plug-in filter creates a mask where the local maxima of the
     * current image are marked (255; unmarked pixels 0).
     * The plug-in can also create watershed-segmented particles (assume a
     * landscape of inverted heights: maxima of the images are now water sinks.
     * For each point in the image, the sink that the water goes to determines which
     * particle it belongs to.
     * Pixels with a level below the lower threshold can be left unprocessed.
     *
     * This plugin works with ROIs, including non-rectangular ROIs.
     * Since this plug-in creates a separate output image it processes
     *    only single images or slices, no stacks.
     *
     * version 09-Nov-2006 Michael Schmid
     * version 21-Nov-2006 Wayne Rasband. Adds "Display Point Selection" option and "Count" output type.
     * version 28-May-2007 Michael Schmid. Preview added, bugfix: minima of calibrated images, uses Arrays.sort (Java2 required)
     */

    public class MaximumFinder implements ExtendedPlugInFilter, DialogListener {
        //filter params
        /** maximum height difference between points that are not counted as separate maxima */
        private static double tolerance = 10;
        /** what type of output to create (see constants below)*/
        private static int outputType;
        /** what type of output to create was chosen in the dialog (see constants below)*/
        private static int dialogOutputType;
        /** Output type single points */
        public final static int SINGLE_POINTS=0;
        /** Output type all points around the maximum within the tolerance */
        public final static int IN_TOLERANCE=1;
        /** Output type watershed-segmented image */
        public final static int SEGMENTED=2;
        /** Do not create image, only mark points */
        public final static int POINT_SELECTION=3;
        /** Do not create an image, just count maxima and add count to Results table*/
        public final static int COUNT=4;
        /** output type names */
        final static String[] outputTypeNames = new String[]
            {"Single Points", "Maxima Within Tolerance", "Segmented Particles", "Point Selection", "Count"};
        /** whether to exclude maxima at the edge of the image*/
        private static boolean excludeOnEdges;
        /** whether to accept maxima only in the thresholded height range*/
        private static boolean useMinThreshold;
        /** whether to find darkest points on light background */
        private static boolean lightBackground;
        private ImagePlus imp;                          // the ImagePlus of the setup call
        private int flags = DOES_ALL|NO_CHANGES|NO_UNDO;// the flags (see interfaces PlugInFilter & ExtendedPlugInFilter)
        private boolean thresholded;                    // whether the input image has a threshold
        private boolean roiSaved;                       // whether the filter has changed the roi and saved the original roi
        private boolean preview;                        // true while dialog is displayed (processing for preview)
        private Vector checkboxes;                      // a reference to the Checkboxes of the dialog
        private boolean thresholdWarningShown = false;  // whether the warning "can't find minima with thresholding" has been shown
        private Label messageArea;                      // reference to the textmessage area for displaying the number of maxima
        int [] dirOffset, dirXoffset, dirYoffset;       // offsets of neighbor pixels for addressing
        final static int IS_LINE=1;                     // a point on a line (as a return type of isLineOrDot)
        final static int IS_DOT=2;                      // an isolated point (as a return type of isLineOrDot)
        /** the following constants are used to set bits corresponding to pixel types */
        final static byte MAXIMUM = (byte)1;            // marks local maxima (irrespective of noise tolerance)
        final static byte LISTED = (byte)2;             // marks points currently in the list
        final static byte PROCESSED = (byte)4;          // marks points processed previously
        final static byte MAX_AREA = (byte)8;           // marks areas near a maximum, within the tolerance
        final static byte EQUAL = (byte)16;             // marks contigous maximum points of equal level
        final static byte MAX_POINT = (byte)32;         // marks a single point standing for a maximum
        final static byte ELIMINATED = (byte)64;        // marks maxima that have been eliminated before watershed
        /** type masks corresponding to the output types */
        final static byte[] outputTypeMasks = new byte[] {MAX_POINT, MAX_AREA, MAX_AREA, MAX_POINT};

        /** Method to return types supported
         * @param arg   Not used by this plugin-filter
         * @param imp   The image to be filtered
         * @return      Code describing supported formats etc.
         * (see ij.plugin.filter.PlugInFilter & ExtendedPlugInFilter)
         */
        public int setup(String arg, ImagePlus imp) {
            this.imp = imp;
            return flags;
        }

        public int showDialog(ImagePlus imp, String command, PlugInFilterRunner pfr) {
            ImageProcessor ip = imp.getProcessor();
            ip.resetBinaryThreshold(); // remove invisible threshold set by MakeBinary and Convert to Mask
            thresholded = ip.getMinThreshold()!=ImageProcessor.NO_THRESHOLD;
            GenericDialog gd = new GenericDialog(command);
            int digits = (ip instanceof FloatProcessor)?2:0;
            String unit = (imp.getCalibration()!=null)?imp.getCalibration().getValueUnit():null;
            unit = (unit==null||unit.equals("Gray Value"))?":":" ("+unit+"):";
            gd.addNumericField("Noise Tolerance"+unit,tolerance, digits);
            gd.addChoice("Output type:", outputTypeNames, outputTypeNames[dialogOutputType]);
            gd.addCheckbox("Exclude Edge Maxima", excludeOnEdges);
            if (thresholded)
                gd.addCheckbox("Above Lower Threshold", useMinThreshold);
            gd.addCheckbox("Light Background", lightBackground);
            gd.addPreviewCheckbox(pfr, "Preview Point Selection");
            gd.addMessage("                        "); //space for number of maxima
            messageArea = (Label)gd.getMessage();
            gd.addDialogListener(this);
            checkboxes = gd.getCheckboxes();
            preview = true;
            gd.showDialog();          //input by the user (or macro) happens here
            if (gd.wasCanceled())
                return DONE;
            preview = false;
            if (!dialogItemChanged(gd, null))   //read parameters
                return DONE;
            IJ.register(this.getClass());       //protect static class variables (parameters) from garbage collection
            return flags;
        } // boolean showDialog

        /** Read the parameters (during preview or after showing the dialog) */
        public boolean dialogItemChanged(GenericDialog gd, AWTEvent e) {
            tolerance = gd.getNextNumber();
            if (tolerance<0) tolerance = 0;
            dialogOutputType = gd.getNextChoiceIndex();
            outputType = preview ? POINT_SELECTION : dialogOutputType;
            excludeOnEdges = gd.getNextBoolean();
            if (thresholded)
                useMinThreshold = gd.getNextBoolean();
            else
                useMinThreshold = false;
            lightBackground = gd.getNextBoolean();
            boolean invertedLut = imp.isInvertedLut();
            if (useMinThreshold && ((invertedLut&&!lightBackground) || (!invertedLut&&lightBackground))) {
                if (!thresholdWarningShown)
                    if (!IJ.showMessageWithCancel(
                        "Find Maxima",
                        "\"Above Lower Threshold\" option cannot be used\n"+
                        "when finding minima (image with light background\n"+
                        "or image with dark background and inverting LUT).")
                        && !preview)
                    return false;               // if faulty input is not detected during preview, "cancel" quits
                thresholdWarningShown = true;
                useMinThreshold = false;
                ((Checkbox)(checkboxes.elementAt(1))).setState(false); //set checkbox
            }
            if (!gd.getPreviewCheckbox().getState())
                messageArea.setText("");        // no "nnn Maxima" message when not previewing
            return (!gd.invalidNumber());
        } // public boolean DialogItemChanged

        /** unused method required by interface ExtendedPlugInFilter */
        public void setNPasses(int nPasses) {}

        /** The plugin is inferred from ImageJ by this method
         * @param ip The image where maxima (or minima) should be found
         */
        public void run(ImageProcessor ip) {
            Roi roi = imp.getRoi();
            if (outputType == POINT_SELECTION && !roiSaved) {
                imp.saveRoi(); // save previous selection so user can restore it
                roiSaved = true;
            }
            if (roi!=null && (!roi.isArea() || outputType==SEGMENTED)) {
                imp.killRoi();
                roi = null;
            }
            boolean invertedLut = imp.isInvertedLut();
            double threshold = useMinThreshold?ip.getMinThreshold():ImageProcessor.NO_THRESHOLD;
            if ((invertedLut&&!lightBackground) || (!invertedLut&&lightBackground)) {
                threshold = ImageProcessor.NO_THRESHOLD;    //don't care about threshold when finding minima
                float[] cTable = ip.getCalibrationTable();
                ip = ip.duplicate();
                if (cTable==null) {                 //invert calibration table for calibrated images
                    ip.invert();
                } else {                            //we are using getPixelValue, so the CalibrationTable must be inverted
                    float[] invertedCTable = new float[cTable.length];
                    for (int i=cTable.length-1; i>=0; i--)
                        invertedCTable[i] = -cTable[i];
                    ip.setCalibrationTable(invertedCTable);
                }
                ip.setRoi(roi);
            }
            ByteProcessor outIp = null;
            outIp = findMaxima(ip, tolerance, threshold, outputType, excludeOnEdges, false); //process the image
            if (outIp == null) return;              //cancelled by user or previewing or no output image
            if (!Prefs.blackBackground)             //normally, we use an inverted LUT, "active" pixels black (255) - like a mask
                outIp.invertLut();
            String resultName;
            if (outputType == SEGMENTED)            //Segmentation required
                resultName = " Segmented";
            else
                resultName = " Maxima";
            //IJ.write ("lowest/highest value"+maxPoints[0].value+","+maxPoints[nMax-1].value);
            String outname = imp.getTitle();
            if (imp.getNSlices()>1)
                outname += "("+imp.getCurrentSlice()+")";
            outname += resultName;
            if (WindowManager.getImage(outname)!=null)
                outname = WindowManager.getUniqueName(outname);
            ImagePlus maxImp = new ImagePlus(outname, outIp);
            Calibration cal = imp.getCalibration().copy();
            cal.disableDensityCalibration();
            maxImp.setCalibration(cal);             //keep the spatial calibration
            maxImp.show();
        } //public void run

        /** Here the processing is done: Find the maxima of an image (does not find minima)
         * @param ip             The input image
         * @param tolerance      Height tolerance: maxima are accepted only if protruding more than this value from the ridge to a higher maximum
         * @param threshold      minimum height of a maximum (uncalibrated); for no minimum height set it to ImageProcessor.NO_THRESHOLD
         * @param outputType     What to mark in output image: SINGLE_POINTS, MAXIMA_EXACT, IN_TOLERANCE or SEGMENTED
         * @param excludeOnEdges Whether to exclude edge maxima
         * @param isEDM          Whether the image is a 16-bit Euclidian Distance Map
         * @return               A new byteProcessor with a normal (uninverted) LUT where the marked points
         *                       are set to 255 (Background 0). Pixels outside of the roi of the input ip are not set.
         *                       Returns null if outputType does not require an output or if cancelled.
         */
        public ByteProcessor findMaxima(ImageProcessor ip, double tolerance, double threshold, int outputType, boolean excludeOnEdges, boolean isEDM) {
            //new ImagePlus("find maxima input", ip.duplicate()).show();
            int width = ip.getWidth();
            int height = ip.getHeight();
            Rectangle roi = ip.getRoi();
            byte[] mask = ip.getMaskArray();
            if (threshold!=ImageProcessor.NO_THRESHOLD && ip.getCalibrationTable()!=null && threshold>0 && threshold<ip.getCalibrationTable().length)
                threshold = ip.getCalibrationTable()[(int)threshold];   //convert threshold to calibrated
            ByteProcessor typeP = new ByteProcessor(width, height);     //will be a notepad for pixel types
            byte[] types = (byte[])typeP.getPixels();
            makeDirectionOffsets(ip);
            float globalMin = Float.MAX_VALUE;
            float globalMax = -Float.MAX_VALUE;
            for (int y=roi.y; y<roi.y+roi.height; y++) {         //find local minimum/maximum now
                for (int x=roi.x; x<roi.x+roi.width; x++) {      //ImageStatistics won't work if we have no ImagePlus
                    float v = ip.getPixelValue(x, y);
                    if (globalMin>v) globalMin = v;
                    if (globalMax<v) globalMax = v;
                }
            }
            if (threshold !=ImageProcessor.NO_THRESHOLD)
                threshold -= (globalMax-globalMin)*1e-6;//avoid rounding errors
            boolean excludeEdgesNow = excludeOnEdges && outputType!=SEGMENTED; //for segmentation, exclusion of edge maxima cannot be done now but has to be done after segmentation

            if (Thread.currentThread().isInterrupted()) return null;
            IJ.showStatus("Getting sorted maxima...");
            MaxPoint[] maxPoints = getSortedMaxPoints(ip, typeP, excludeEdgesNow, isEDM, globalMin, threshold);
            if (Thread.currentThread().isInterrupted()) return null;
            IJ.showStatus("Analyzing  maxima...");
            analyzeAndMarkMaxima(ip, typeP, maxPoints, excludeEdgesNow, isEDM, globalMin, tolerance);
            //new ImagePlus("Pixel types",typeP.duplicate()).show();
            if (outputType==COUNT || outputType==POINT_SELECTION)
                return null;

            ByteProcessor outIp;
            byte[] pixels;
            if (outputType == SEGMENTED) {                  //Segmentation required, convert to 8bit (also for 8-bit images, since the calibration may have a negative slope)
                //create 8-bit image from ip with background 0 and maximum areas 255
                outIp = make8bit(ip, typeP, isEDM, globalMin, globalMax, threshold);
                cleanupMaxima(outIp, typeP, maxPoints);     //eliminate all the small maxima (i.e. those outside MAX_AREA)
                //new ImagePlus("pixel types", typeP).show();
                //if (IJ.debugMode)
                //new ImagePlus("pre-watershed", outIp.duplicate()).show();
                if (!watershedSegment(outIp))               //do watershed segmentation
                    return null;                            //if user-cancelled, return
                if (!isEDM) cleanupExtraLines(outIp);       //eliminate lines due to local minima (none in EDM)
                watershedPostProcess(outIp);                //levels to binary image
                if (excludeOnEdges) deleteEdgeParticles(outIp, typeP);
            } else {                                        //outputType other than SEGMENTED
                for (int i=0; i<width*height; i++)
                    types[i] = (byte)(((types[i]&outputTypeMasks[outputType])!=0)?255:0);
                outIp = typeP;
            }
            byte[] outPixels = (byte[])outIp.getPixels();
            //IJ.write("roi: "+roi.toString());
            if (roi!=null) {
                for (int y=0, i=0; y<outIp.getHeight(); y++) { //delete everything outside roi
                    for (int x=0; x<outIp.getWidth(); x++, i++) {
                        if (x<roi.x || x>=roi.x+roi.width || y<roi.y || y>=roi.y+roi.height) outPixels[i] = (byte)0;
                        else if (mask !=null && (mask[x-roi.x + roi.width*(y-roi.y)]==0)) outPixels[i] = (byte)0;
                    }
                }
            }

            return outIp;
        } // public ByteProcessor MaximumFinder

        /** find all local maxima (irrespective whether they finally qualify as maxima or not)
         * @param ip the image to be analyzed
         * @param typeP A byte image, same size as ip, where the maximum points are mrked as MAXIMUM
         * @param excludeEdgesNow whether to exclude edge pixels
         * @param isEDM whether ip is a 16-bit Euclidian distance map
         * @param globalMin the minimum value of the image or roi
         * @param threshold the threshold (calibrated) below which no pixels are processed. Ignored if ImageProcessor.NO_THRESHOLD
         * @return an array of x, y coordinates and heights, sorted by height. Do not use the positions of these points in typeP (marked as MAXIMUM there).
         * returns null if the thread is interrupted.
         */
       MaxPoint[] getSortedMaxPoints(ImageProcessor ip, ByteProcessor typeP, boolean excludeEdgesNow, boolean isEDM, float globalMin, double threshold) {
            int width = ip.getWidth();
            int height = ip.getHeight();
            Rectangle roi = ip.getRoi();
            byte[] types =  (byte[])typeP.getPixels();
            int nMax = 0;  //counts local maxima
            boolean checkThreshold = threshold!=ImageProcessor.NO_THRESHOLD;
            Thread thread = Thread.currentThread();
            for (int y=roi.y, i=0; y<roi.y+roi.height; y++) {         // find local maxima now
                if (y%50==0 && thread.isInterrupted()) return null;
                for (int x=roi.x; x<roi.x+roi.width; x++, i++) {      //  (don't care for roi now; output will be restricted to the roi at the very end)
                    float v = ip.getPixelValue(x, y);
                    if (v==globalMin) continue;
                    if (excludeEdgesNow && (x==0 || x==width-1 || y==0 || y==height-1)) continue;
                    if (checkThreshold && v<threshold) continue;
                    boolean isMax = true;
                    for (int d=0; d<8; d++) {
                        if (isWithin(ip, x, y, d)) {
                            if (ip.getPixelValue(x+dirXoffset[d], y+dirYoffset[d]) > v)
                                isMax = false;
                        }
                    }
                    if (isMax) {
                        types[i] = MAXIMUM;
                        nMax++;
                    }
                } // for x
            } // for y
            if (thread.isInterrupted()) return null;

            MaxPoint [] maxPoints = new MaxPoint[nMax];
            int iMax = 0;
            for (int y=roi.y, i=0; y<roi.y+roi.height; y++)           //enter all maxima into an array
                for (int x=roi.x; x<roi.x+roi.width; x++, i++)
                    if (types[i]==MAXIMUM) {
                        maxPoints[iMax] = new MaxPoint((short)x, (short)y, isEDM?trueEdmHeight(x,y,ip):ip.getPixelValue(x,y));
                        //if(ip.getPixelValue(x,y)>=250) IJ.write("x,y,value="+maxPoints[iMax].x+","+maxPoints[iMax].y+","+maxPoints[iMax].value);
                        iMax++;
                    }
            if (thread.isInterrupted()) return null;
            Arrays.sort(maxPoints);                                 //sort the maxima by height
            return maxPoints;
        } //getSortedMaxPoints

       /** Check all maxima in list maxPoints, mark type of the points in typeP
        * @param ip             the image to be analyzed
        * @param typeP          here the point types are marked.
        * @param maxPoints      input: a list of all local maxima, sorted by height
        * @param excludeEdgesNow whether to avoid edge maxima
        * @param isEDM          whether ip is a 16-bit Euclidian distance map
        * @param globalMin      minimum pixel value in ip
        * @param tolerance      minimum pixel value difference for two separate maxima
        */
       void analyzeAndMarkMaxima(ImageProcessor ip, ByteProcessor typeP, MaxPoint[] maxPoints, boolean excludeEdgesNow, boolean isEDM, float globalMin, double tolerance) {
            int width = ip.getWidth();
            int height = ip.getHeight();
            byte[] types =  (byte[])typeP.getPixels();
            int nMax = maxPoints.length;
            short[] xList = new short[width*height];    //here we enter points starting from a maximum
            short[] yList = new short[width*height];
            Vector xyVector = null;
            Roi roi = null;
            boolean displayOrCount = imp!=null && (outputType==POINT_SELECTION||outputType==COUNT);
            for (int iMax=nMax-1; iMax>=0; iMax--) {    //process all maxima now, starting from the highest
                if (iMax%100 == 0 && Thread.currentThread().isInterrupted()) return;
                float v = maxPoints[iMax].value;
                if (v==globalMin) break;
                int offset = maxPoints[iMax].x + width*maxPoints[iMax].y;
                if ((types[offset]&PROCESSED)!=0)       //this maximum has been reached from another one, skip it
                    continue;
                xList[0] = maxPoints[iMax].x;           //we create a list of connected points and start the list
                yList[0] = maxPoints[iMax].y;           //  at the current maximum
                types[offset] |= (EQUAL|LISTED);        //mark first point as equal height (to itself) and listed
                int listLen = 1;                        //number of elements in the list
                int listI = 0;                          //index of current element in the list
                boolean isEdgeMaximum = (xList[0]==0 || xList[0]==width-1 || yList[0]==0 || yList[0]==height-1);
                boolean maxPossible = true;             //it may be a true maxmum
                double xEqual = xList[0];               //for creating a single point: determine average over the
                double yEqual = yList[0];               //  coordinates of contiguous equal-height points
                int nEqual = 1;                         //counts xEqual/yEqual points that we use for averaging
                do {
                    offset = xList[listI] + width*yList[listI];
                    for (int d=0; d<8; d++) {           //analyze all neighbors (in 8 directions) at the same level
                        int offset2 = offset+dirOffset[d];
                        if (isWithin(ip, xList[listI], yList[listI], d) && (types[offset2]&LISTED)==0) {
                            if ((types[offset2]&PROCESSED)!=0) {
                                maxPossible = false;    //we have reached a point processed previously, thus it is no maximum now
                                //if(xList[0]>510&&xList[0]<77)IJ.write("stop at processed neighbor from x,y="+xList[listI]+","+yList[listI]+", dir="+d);
                                break;
                            }
                            int x2 = xList[listI]+dirXoffset[d];
                            int y2 = yList[listI]+dirYoffset[d];
                            float v2 = ip.getPixelValue(x2, y2);
                            if (isEDM && (v2 <=v-(float)tolerance)) v2 = trueEdmHeight(x2, y2, ip); //correcting for EDM peaks may move the point up
                            if (v2 > v) {
                                maxPossible = false;    //we have reached a higher point, thus it is no maximum
                                //if(xList[0]>510&&xList[0]<77)IJ.write("stop at higher neighbor from x,y="+xList[listI]+","+yList[listI]+", dir,value,value2,v2-v="+d+","+v+","+v2+","+(v2-v));
                                break;
                            } else if (v2 >= v-(float)tolerance) {
                                xList[listLen] = (short)(x2);
                                yList[listLen] = (short)(y2);
                                listLen++;              //we have found a new point within the tolerance
                                types[offset2] |= LISTED;
                                if (x2==0 || x2==width-1 || y2==0 || y2==height-1) {
                                    isEdgeMaximum = true;
                                    if (excludeEdgesNow) {
                                        maxPossible = false;
                                        break;          //we have an edge maximum;
                                    }
                                }
                                if (v2==v) {            //prepare finding center of equal points (in case single point needed)
                                    types[offset2] |= EQUAL;
                                    xEqual += x2;
                                    yEqual += y2;
                                    nEqual ++;
                                }
                            }
                        } // if isWithin & not LISTED
                    } // for directions d
                    listI++;
                } while (listI < listLen);
                byte resetMask = (byte)~(maxPossible?LISTED:(LISTED|EQUAL));
                xEqual /= nEqual;
                yEqual /= nEqual;
                double minDist2 = 1e20;
                int nearestI = 0;
                for (listI=0; listI<listLen; listI++) {
                    offset = xList[listI] + width*yList[listI];
                    types[offset] &= resetMask;     //reset attributes no longer needed
                    types[offset] |= PROCESSED;     //mark as processed
                    if (maxPossible) {
                        types[offset] |= MAX_AREA;
                        if ((types[offset]&EQUAL)!=0) {
                            double dist2 = (xEqual-xList[listI])*(xEqual-xList[listI]) + (yEqual-yList[listI])*(yEqual-yList[listI]);
                            if (dist2 < minDist2) {
                                minDist2 = dist2;   //this could be the best "single maximum" point
                                nearestI = listI;
                            }
                        }
                    }
                } // for listI
                if (maxPossible) {
                    types[xList[nearestI] + width*yList[nearestI]] |= MAX_POINT;
                    if (displayOrCount && !(this.excludeOnEdges && isEdgeMaximum)) {
                        if (xyVector==null) {
                            xyVector = new Vector();
                            roi = imp.getRoi();
                        }
                        short mpx = xList[nearestI];
                        short mpy = yList[nearestI];
                        if (roi==null || roi.contains(mpx, mpy))
                            xyVector.addElement(new MaxPoint(mpx, mpy, 0f));
                    }
                }
            } // for all maxima iMax
            if (Thread.currentThread().isInterrupted()) return;
            if (displayOrCount && xyVector!=null) {
                int npoints = xyVector.size();
                if (outputType == POINT_SELECTION) {
                    int[] xpoints = new int[npoints];
                    int[] ypoints = new int[npoints];
                    MaxPoint mp;
                    for (int i=0; i<npoints; i++) {
                        mp = (MaxPoint)xyVector.elementAt(i);
                        xpoints[i] = mp.x;
                        ypoints[i] = mp.y;
                    }
                    imp.setRoi(new PointRoi(xpoints, ypoints, npoints));
                }
                if (outputType==COUNT) {
                    ResultsTable rt = ResultsTable.getResultsTable();
                    rt.incrementCounter();
                    rt.setValue("Count", rt.getCounter()-1, npoints);
                    rt.show("Results");
                }
            }
            if (preview)
                messageArea.setText((xyVector==null ? 0 : xyVector.size())+" Maxima");
        } //void analyzeAndMarkMaxima

       /** Create an 8-bit image by scaling the pixel values of ip (background 0) and mark maximum areas as 255.
        * For use as input for watershed segmentation
        * @param ip         The original image that should be segmented
        * @param typeP      Pixel types in ip
        * @param isEDM      Whether ip is an Euclidian distance map
        * @param globalMin  The minimum pixel value of ip
        * @param globalMax  The maximum pixel value of ip
        * @param threshold  Pixels of ip below this value (calibrated) are considered background. Ignored if ImageProcessor.NO_THRESHOLD
        * @return           The 8-bit output image.
        */
        ByteProcessor make8bit(ImageProcessor ip, ByteProcessor typeP, boolean isEDM, float globalMin, float globalMax, double threshold) {
            int width = ip.getWidth();
            int height = ip.getHeight();
            byte[] types = (byte[])typeP.getPixels();
            double minValue = (threshold == ImageProcessor.NO_THRESHOLD)?globalMin:threshold;
            double offset = minValue - (globalMax-minValue)*(1./253/2-1e-6); //everything above minValue should become >(byte)0
            double factor = 253/(globalMax-minValue);
            //IJ.write("min,max="+minValue+","+globalMax+"; offset, 1/factor="+offset+", "+(1./factor));
            if (isEDM && factor >1./EDM.ONE)
                factor = 1./EDM.ONE;   // in EDM, no better resolution is needed
            ByteProcessor outIp = new ByteProcessor(width, height);
            byte[] pixels = (byte[])outIp.getPixels();  //convert possibly calibrated image to byte without damaging threshold (setMinAndMax would kill threshold)
            long v;
            for (int y=0, i=0; y<height; y++) {
                for (int x=0; x<width; x++, i++) {
                    v = Math.round((ip.getPixelValue(x, y)-offset)*factor);
                    if (v<0) pixels[i] = (byte)0;
                    else if (v<=255) pixels[i] = (byte)(v&255);
                    else pixels[i] = (byte) 255;
                }
            }

            pixels = (byte[])outIp.getPixels();
            //new ImagePlus("pre-threshold", outIp.duplicate()).show();
            if (threshold != ImageProcessor.NO_THRESHOLD) {
                for (int y=0; y<height; y++) {          //Set out-of-threshold to background (inactive)
                    for (int x=0; x<width; x++) {
                        if (ip.getPixelValue(x,y) < threshold)
                            pixels[x+y*width] = (byte)0;
                    }
                }
            }
            for (int i=0; i<width*height; i++) {
                if ((pixels[i]&255) == 255)             //pixel value 255 is reserved
                    pixels[i]--;
                if((types[i]&MAX_AREA)!=0)
                    pixels[i] = (byte)255;              //prepare watershed by setting "true" maxima+surroundings to 255
            }
            return outIp;
        } // byteProcessor make8bit

       /** Get estimated "true" height of a maximum or saddle point of a Euclidian Distance Map.
         * This is needed since the point sampled is not necessarily at the highest position.
         * Note: for simplicity, in contrast to EDM we don't care about the Sqrt(5) distance here.
         * @param x     x-position of the point
         * @param y     y-position of the point
         * @param ip    the EDM (16-bit, scaling as defined in EDM class)
         * @return      estimated height
         */
        int trueEdmHeight(int x, int y, ImageProcessor ip) {
            int xmax = ip.getWidth()-1;
            int ymax = ip.getHeight()-1;
            int v =  ip.getPixel(x, y);
            if (x==0 || y==0 || x==xmax || y==ymax || v==0) {
                return v;                               //don't recalculate for edge pixels or background
            } else {
                int one = EDM.ONE;
                int sqrt2 = EDM.SQRT2;
                int trueH = v+sqrt2/2;                  //true height can never by higher than this
                boolean ridgeOrMax = false;
                for (int d=0; d<4; d++) {               //for all directions halfway around:
                    int d2 = (d+4)%8;                   //get the opposite direction and neighbors
                    int v1 = ip.getPixel(x+dirXoffset[d], y+dirYoffset[d]);
                    int v2 = ip.getPixel(x+dirXoffset[d2], y+dirYoffset[d2]);
                    int h;
                    if (v>=v1 && v>=v2) {
                        ridgeOrMax = true;
                        h = (v1 + v2)/2;
                    } else {
                        h = Math.min(v1, v2);
                    }
                    h += (d%2==0) ? one : sqrt2;        //in diagonal directions, distance is sqrt2
                    if (trueH > h) trueH = h;
                }
                if (!ridgeOrMax) trueH = v;
                return trueH;
            }
        }

        /** eliminate unmarked maxima for use by watershed. Starting from each previous maximum,
         * explore the surrounding down to successively lower levels until a marked maximum is
         * touched (or the plateau of a previously eliminated maximum leads to a marked maximum).
         * Then set all the points above this value to this value
         * @param outIp     the image containing the pixel values
         * @param typeP     the types of the pixels are mekred here
         * @param maxPoints array containing the coordinates of all maxima that might be relevant
         */
        void cleanupMaxima(ByteProcessor outIp, ByteProcessor typeP, MaxPoint[] maxPoints) {
            byte[] pixels = (byte[])outIp.getPixels();
            byte[] types = (byte[])typeP.getPixels();
            int width = outIp.getWidth();
            int height = outIp.getHeight();
            int nMax = maxPoints.length;
            short[] xList = new short[width*height];
            short[] yList = new short[width*height];
            for (int iMax = nMax-1; iMax>=0; iMax--) {
                int offset = maxPoints[iMax].x + width*maxPoints[iMax].y;
                //if (maxPoints[iMax].x==122) IJ.write("max#"+iMax+" at x,y="+maxPoints[iMax].x+","+maxPoints[iMax].y+": type="+types[offset]);
                if ((types[offset]&(MAX_AREA|ELIMINATED))!=0) continue;
                int level = pixels[offset]&255;
                int loLevel = level+1;
                xList[0] = maxPoints[iMax].x;           //we start the list at the current maximum
                yList[0] = maxPoints[iMax].y;
                //if (xList[0]==122) IJ.write("max#"+iMax+" at x,y="+xList[0]+","+yList[0]+"; level="+level);
                types[offset] |= LISTED;                //mark first point as listed
                int listLen = 1;                        //number of elements in the list
                int lastLen = 1;
                int listI = 0;                          //index of current element in the list
                boolean saddleFound = false;
                while (!saddleFound && loLevel >0) {
                    loLevel--;
                    lastLen = listLen;                  //remember end of list for previous level
                    listI = 0;                          //in each level, start analyzing the neighbors of all pixels
                    do {                                //for all pixels listed so far
                        offset = xList[listI] + width*yList[listI];
                        for (int d=0; d<8; d++) {       //analyze all neighbors (in 8 directions) at the same level
                            int offset2 = offset+dirOffset[d];
                            if (isWithin(typeP, xList[listI], yList[listI], d) && (types[offset2]&LISTED)==0) {
                                if ((types[offset2]&MAX_AREA)!=0 || (((types[offset2]&ELIMINATED)!=0) && (pixels[offset2]&255)>=loLevel)) {
                                    saddleFound = true; //we have reached a point touching a "true" maximum...
                                    //if (xList[0]==122) IJ.write("saddle found at level="+loLevel+"; x,y="+xList[listI]+","+yList[listI]+", dir="+d);
                                    break;              //...or a level not lower, but touching a "true" maximum
                                } else if ((pixels[offset2]&255)>=loLevel && (types[offset2]&ELIMINATED)==0) {
                                    xList[listLen] = (short)(xList[listI]+dirXoffset[d]);
                                    yList[listLen] = (short)(yList[listI]+dirYoffset[d]);
                                    listLen++;          //we have found a new point to be processed
                                    types[offset2] |= LISTED;
                                }
                            } // if isWithin & not LISTED
                        } // for directions d
                        if (saddleFound) break;         //no reason to search any further
                        listI++;
                    } while (listI < listLen);
                } // while !levelFound && loLevel>=0
                for (listI=0; listI<listLen; listI++)   //reset attribute since we may come to this place again
                    types[xList[listI] + width*yList[listI]] &= ~LISTED;
                for (listI=0; listI<lastLen; listI++) { //for all points higher than the level of the saddle point
                    offset = xList[listI] + width*yList[listI];
                    pixels[offset] = (byte)loLevel;     //set pixel value to the level of the saddle point
                    types[offset] |= ELIMINATED;        //mark as processed: there can't be a local maximum in this area
                }
            } // for all maxima iMax
        }

        /** Delete  single dots and lines ending somewhere within a segmented particle
         * Needed for post-processing watershed-segmented images that can have local minima
         * @param ip 8-bit image with background = 0, lines between 1 and 254 and segmented particles = 255
         */
        void cleanupExtraLines(ImageProcessor ip) {
            int width = ip.getWidth();
            int height = ip.getHeight();
            byte[] pixels =  (byte[])ip.getPixels();
            for (int y=0, i=0; y<height; y++) {
                for (int x=0; x<width; x++,i++) {
                    int v = pixels[i]&255;
                    if (v<255 && v>0) {
                        int type = isLineOrDot(ip, x, y);
                        if (type==IS_DOT) {
                            pixels[i] = (byte)255;                  //delete the point;
                        } else if (type==IS_LINE) {
                            int xEnd = x;
                            int yEnd = y;
                            boolean endFound = true;
                            while (endFound) {
                                pixels[xEnd + width*yEnd] = (byte)255;  //delete the point
                                //if(movie.getSize()<100)movie.addSlice("part-cleaned", ip.duplicate());
                                endFound = false;
                                for (int d=0; d<8; d++) {               //analyze neighbors of the point
                                    if (isWithin(ip, xEnd, yEnd, d)) {
                                        v = pixels[xEnd + width*yEnd + dirOffset[d]]&255;
                                        //if(x>210&&x<215&&y==13)IJ.write("x,y start="+x+","+y+": look at="+xEnd+","+yEnd+"+dir "+d+": v="+v);
                                        if (v<255 && v>0 && isLineOrDot(ip, xEnd+dirXoffset[d], yEnd+dirYoffset[d])==IS_LINE) {
                                            xEnd += dirXoffset[d];
                                            yEnd += dirYoffset[d];
                                            endFound = true;
                                            break;
                                        }
                                    }
                                } // for directions d
                            } // while (endFound)
                        } // if IS_LINE
                    } // if v<255 && v>0
                } // for x
            } // for y
        } //cleanupExtraLines

        /** analyze the neighbors of a pixel (x, y) in a byte image; pixels <255 are
         * considered part of lines.
         * @param ip
         * @param x
         * @param y
         * @return  IS_LINE if pixel is part of a line, IS_DOT if a single dot
         */
        int isLineOrDot(ImageProcessor ip, int x, int y) {
            int result = 0;
            int width = ip.getWidth();
            int height = ip.getHeight();
            byte[] pixels = (byte[])ip.getPixels();
            int offset = x + y*width;
            int whiteNeighbors = 0;             //counts neighbors that are not part of a line
            int countTransitions = 0;
            boolean pixelSet;
            boolean prevPixelSet = true;
            boolean firstPixelSet = true;       //initialize to make the compiler happy
            for (int d=0; d<8; d++) {           //walk around the point and note every no-line->line transition
                if (isWithin(ip, x, y, d)) {
                    pixelSet = (pixels[offset+dirOffset[d]]!=(byte)255);
                    if (!pixelSet) whiteNeighbors++;
                } else {
                    pixelSet = true;
                }
                if (pixelSet && !prevPixelSet)
                    countTransitions ++;
                prevPixelSet = pixelSet;
                if (d==0)
                    firstPixelSet = pixelSet;
            }
            if (firstPixelSet && !prevPixelSet)
                countTransitions ++;
                //if (x>=210&&x<=215 && y>=10 && y<=17)IJ.write("x,y="+x+","+y+": transitions="+countTransitions);
            if (countTransitions==1 && whiteNeighbors>=5)
                result = IS_LINE;
            else if (whiteNeighbors==8)
                result = IS_DOT;
            return result;
        } //isLineEnd

        /** after watershed, set all pixels in the background and segmentation lines to 0
         */
        private void watershedPostProcess(ImageProcessor ip) {
            //new ImagePlus("before postprocess",ip.duplicate()).show();
            byte[] pixels = (byte[])ip.getPixels();
            int size = ip.getWidth()*ip.getHeight();
            for (int i=0; i<size; i++) {
               if ((pixels[i]&255)<255)
                    pixels[i] = (byte)0;
            }
            //new ImagePlus("after postprocess",ip.duplicate()).show();
        }

        /** delete particles corresponding to edge maxima
         * @param typeP Here the pixel types of the original image are noted,
         * pixels with bit MAX_AREA at the edge are considered indicators of an edge maximum.
         * @param ip the image resulting from watershed segmentaiton
         * (foreground pixels, i.e. particles, are 255, background 0)
         */
        void deleteEdgeParticles(ByteProcessor ip, ByteProcessor typeP) {
            int width = ip.getWidth();
            int height = ip.getHeight();
            byte[] pixels = (byte[])ip.getPixels();
            byte[] types = (byte[])typeP.getPixels();
            width = ip.getWidth();
            height = ip.getHeight();
            ip.setValue(0);
            Wand wand = new Wand(ip);
            for (int x=0; x<width; x++) {
                int y = 0;
                if ((types[x+y*width]&MAX_AREA) != 0 && pixels[x+y*width] != 0)
                    deleteParticle(x,y,ip,wand);
                y = height - 1;
                if ((types[x+y*width]&MAX_AREA) != 0 && pixels[x+y*width] != 0)
                    deleteParticle(x,y,ip,wand);
            }
            for (int y=1; y<height-1; y++) {
                int x = 0;
                if ((types[x+y*width]&MAX_AREA) != 0 && pixels[x+y*width] != 0)
                    deleteParticle(x,y,ip,wand);
                x = width - 1;
                if ((types[x+y*width]&MAX_AREA) != 0 && pixels[x+y*width] != 0)
                    deleteParticle(x,y,ip,wand);
            }
        }
        /** delete a particle (set from value 255 to current fill value).
         * Position x,y must be within the particle
         */
        void deleteParticle(int x, int y, ByteProcessor ip, Wand wand) {
            wand.autoOutline(x, y, 255, 255);
            if (wand.npoints==0) {
                IJ.log("wand error selecting edge particle at x, y = "+x+", "+y);
                return;
            }
            Roi roi = new PolygonRoi(wand.xpoints, wand.ypoints, wand.npoints, Roi.TRACED_ROI);
            ip.snapshot();                  //prepare for reset outside of mask
            ip.setRoi(roi);
            ip.fill();
            ip.reset(ip.getMask());
        }

            /** Do watershed segmentation on a byte image, with the start points (maxima)
         * set to 255 and the background set to 0. The image should not have any local maxima
         * other than the marked ones. Local minima will lead to artifacts that can be removed
         * later. On output, all particles will be set to 255, segmentation lines remain at their
         * old value.
         * @param ip  The byteProcessor containing the image, with size given by the class variables width and height
         * @return    false if canceled by the user (note: can be cancelled only if called by "run" with a known ImagePlus)
         */
        private boolean watershedSegment(ByteProcessor ip) {
            int width = ip.getWidth();
            int height = ip.getHeight();
            boolean debug = IJ.debugMode;
            ImageStack movie=null;
            if (debug) {
                movie = new ImageStack(ip.getWidth(), ip.getHeight());
                movie.addSlice("pre-watershed EDM", ip.duplicate());
            }
            byte[] pixels = (byte[])ip.getPixels();
            // create arrays with the coordinates of all points between value 1 and 254
            //This method, suggested by Stein Roervik (stein_at_kjemi-dot-unit-dot-no),
            //greatly speeds up the watershed segmentation routine.
            int[] histogram = ip.getHistogram();
            int arraySize = width*height - histogram[0] -histogram[255];
            short[] xCoordinate = new short[arraySize];
            short[] yCoordinate = new short[arraySize];
            int highestValue = 0;
            int offset = 0;
            int[] levelStart = new int[256];
            for (int v=1; v<255; v++) {
                levelStart[v] = offset;
                offset += histogram[v];
                if (histogram[v]>0) highestValue = v;
            }
            int[] levelOffset = new int[highestValue + 1];
            for (int y=0, i=0; y<height; y++) {
                for (int x=0; x<width; x++, i++) {
                    int v = pixels[i]&255;
                    if (v>0 && v<255) {
                        offset = levelStart[v] + levelOffset[v];
                        xCoordinate[offset] = (short) x;
                        yCoordinate[offset] = (short) y;
                        levelOffset[v] ++;
                    }
               } //for x
            } //for y
            // now do the segmentation, starting at the highest level and working down.
            // At each level, dilate the particle, constrained to pixels whose values are
            // at that level and also constrained to prevent features from merging.
            int[] table = makeFateTable();
            IJ.showStatus("Segmenting (Esc to cancel)");
            ImageProcessor ip2 = ip.duplicate();
            for (int level=highestValue; level>=1; level--) {
                IJ.showProgress(highestValue-level, highestValue);
                int idle = 1;
                while(true) {                   // break the loop at any point after 8 idle processLevel calls
                    if (processLevel(level, 7, ip, ip2, table, histogram, levelStart, xCoordinate, yCoordinate))
                        idle = 0;
                    if (idle++ >=8) break;
                    if (processLevel(level, 3, ip, ip2, table, histogram, levelStart, xCoordinate, yCoordinate))
                        idle = 0;
                    if (idle++ >=8) break;
                    if (processLevel(level, 1, ip, ip2, table, histogram, levelStart, xCoordinate, yCoordinate))
                        idle = 0;
                    if (idle++ >=8) break;
                    if (processLevel(level, 5, ip, ip2, table, histogram, levelStart, xCoordinate, yCoordinate))
                        idle = 0;
                    if (idle++ >=8) break;
                    //IJ.write("diagonal only; level="+level+"; countW="+countW);
                    if (processLevel(level, 0, ip, ip2, table, histogram, levelStart, xCoordinate, yCoordinate))
                        idle = 0;
                    if (idle++ >=8) break;
                    if (processLevel(level, 4, ip, ip2, table, histogram, levelStart, xCoordinate, yCoordinate))
                        idle = 0;
                    if (idle++ >=8) break;
                    if (processLevel(level, 2, ip, ip2, table, histogram, levelStart, xCoordinate, yCoordinate))
                        idle = 0;
                    if (idle++ >=8) break;
                    if (processLevel(level, 6, ip, ip2, table, histogram, levelStart, xCoordinate, yCoordinate))
                        idle = 0;
                    if (idle++ >=8) break;
                    //IJ.write("All directions; level="+level+"; countW="+countW);
                }
                if (IJ.escapePressed()) {  // cancelled by the user
                    IJ.beep();
                    IJ.showProgress(1.0);
                    return false;
                }
                if (debug && (level>245 || level<30))
                    movie.addSlice("level "+level, ip.duplicate());
            }
            if (debug)
                new ImagePlus("Segmentation Movie", movie).show();
            IJ.showProgress(1.0);
            return true;
        }
        /** Creates the lookup table used by the watershed function for dilating the particles.
         * The algorithm allows dilation in both straight and diagonal directions.
         * There is an entry in the table for each possible 3x3 neighborhood:
         *          x-1          x          x+1
         *  y-1    128            1          2
         *  y       64     pxl_unset_yet     4
         *  y+1     32           16          8
         * (to find throws entry, sum up the numbers of the neighboring pixels set; e.g.
         * entry 6=2+4 if only the pixels (x,y-1) and (x+1, y-1) are set.
         * A pixel is added on the 1st pass if bit 0 (2^0 = 1) is set,
         * on the 2nd pass if bit 1 (2^1 = 2) is set, etc.
         * pass gives the direction of rotation, with 0 = to top left (x--,y--), 1 to top,
         * and clockwise up to 7 = to the left (x--).
         * E.g. 4 = add on 3rd pass, 3 = add on either 1st or 2nd pass.
         */
        private int[] makeFateTable() {
            int[] table = new int[256];
            boolean[] isSet = new boolean[8];
            for (int item=0; item<256; item++) {        //dissect into pixels
                for (int i=0, mask=1; i<8; i++) {
                    isSet[i] = (item&mask)==mask;
                    mask*=2;
                }
                for (int i=0, mask=1; i<8; i++) {       //we dilate in the direction opposite to the direction of the existing neighbors
                    if (isSet[(i+4)%8]) table[item] |= mask;
                    mask*=2;
                }
                for (int i=0; i<8; i+=2)                //if side pixels are set, for counting transitions it is as good as if the adjacent edges were also set
                    if (isSet[i]) {
                        isSet[(i+1)%8] = true;
                        isSet[(i+7)%8] = true;
                    }
                int transitions=0;
                for (int i=0, mask=1; i<8; i++) {
                    if (isSet[i] != isSet[(i+1)%8])
                        transitions++;
                }
                if (transitions>=4) {                   //if neighbors contain more than one region, dilation ito this pixel is forbidden
                    table[item] = 0;
                } else {
                }
            }
            return table;
        } // makeFateTable

        /** dilate the UEP on one level by one pixel in the direction specified by step, i.e., set pixels to 255
         * @param level the level of the EDM that should be processed (all other pixels are untouched)
         * @param pass gives direction of dilation, see makeFateTable
         * @param ip1 the EDM with the segmeted blobs successively getting set to 255
         * @param ip2 a processor used as scratch storage, must hold the same data as ip1 upon entry.
         *            This method ensures that ip2 equals ip1 upon return
         * @param table the fateTable
         * class variables used:
         *  xCoordinate[], yCoordinate[]    list of pixel coordinates sorted by level (in sequence of y, x within each level)
         *  levelStart[]                    offsets of the level in xCoordinate[], yCoordinate[]
         * @return                          true if pixels have been changed
         */
        private boolean processLevel(int level, int pass, ImageProcessor ip1, ImageProcessor ip2, int[] table,
        int[] histogram, int[] levelStart, short[] xCoordinate, short[] yCoordinate) {
            int rowSize = ip1.getWidth();
            int height = ip1.getHeight();
            int xmax = rowSize-1;
            int ymax = ip1.getHeight()-1;
            byte[] pixels1 = (byte[])ip1.getPixels();
            byte[] pixels2 = (byte[])ip2.getPixels();
            boolean changed = false;
            for (int i=0; i<histogram[level]; i++) {
                int coordOffset = levelStart[level] + i;
                int x = xCoordinate[coordOffset];
                int y = yCoordinate[coordOffset];
                int offset = x + y*rowSize;
                int index;
                if ((pixels2[offset]&255)!=255) {
                    index = 0;
                    if (y>0 && (pixels2[offset-rowSize]&255)==255)
                        index ^= 1;
                    if (x<xmax && y>0 && (pixels2[offset-rowSize+1]&255)==255)
                        index ^= 2;
                    if (x<xmax && (pixels2[offset+1]&255)==255)
                        index ^= 4;
                    if (x<xmax && y<ymax && (pixels2[offset+rowSize+1]&255)==255)
                        index ^= 8;
                    if (y<ymax && (pixels2[offset+rowSize]&255)==255)
                        index ^= 16;
                    if (x>0 && y<ymax && (pixels2[offset+rowSize-1]&255)==255)
                        index ^= 32;
                    if (x>0 && (pixels2[offset-1]&255)==255)
                        index ^= 64;
                    if (x>0 && y>0 && (pixels2[offset-rowSize-1]&255)==255)
                        index ^= 128;
                    switch (pass) {
                        case 0: if ((table[index]&1)==1) {
                            pixels1[offset] = (byte)255;
                            changed = true;
                        }
                        break;
                        case 1: if ((table[index]&2)==2) {
                            pixels1[offset] = (byte)255;
                            changed = true;
                        }
                        break;
                        case 2: if ((table[index]&4)==4) {
                            pixels1[offset] = (byte)255;
                            changed = true;
                        }
                        break;
                        case 3: if ((table[index]&8)==8) {
                            pixels1[offset] = (byte)255;
                            changed = true;
                        }
                        break;
                        case 4: if ((table[index]&16)==16) {
                            pixels1[offset] = (byte)255;
                            changed = true;
                        }
                        break;
                        case 5: if ((table[index]&32)==32) {
                            pixels1[offset] = (byte)255;
                            changed = true;
                        }
                        break;
                        case 6: if ((table[index]&64)==64) {
                            pixels1[offset] = (byte)255;
                            changed = true;
                        }
                        break;
                        case 7: if ((table[index]&128)==128) {
                            pixels1[offset] = (byte)255;
                            changed = true;
                        }
                        break;
                    } // switch
                } // if .. !=255
            } // for pixel i
            //IJ.write("level="+level+", pass="+pass+", changed="+changed);
            if (changed)                    //make sure that ip2 is updated for the next time
                System.arraycopy(pixels1, 0, pixels2, 0, rowSize*height);

            return changed;
        } //processLevel

        /** create an array of offsets within a pixel array for directions in clockwise order: 0=(x,y-1), 1=(x+1,y-1), ... 7=(x-1,y)
         * uses class variable width: width of the image where the pixels should be addressed
         * returns as class variables: the arrays of the offsets to the 8 neighboring pixels
         */
        private void makeDirectionOffsets(ImageProcessor ip) {
            int width = ip.getWidth();
            int height = ip.getHeight();
            dirOffset = new int[] { -width, -width+1, +1, +width+1, +width, +width-1,   -1, -width-1 };
            dirXoffset = new int[] {    0,       1,     1,     1,       0,      -1,      -1,    -1    };
            dirYoffset = new int[] {   -1,      -1,     0,     1,       1,       1,       0,    -1,   };
        }

        /** returns whether the neighbor in a given direction is within the image
         * NOTE: it is assumed that the pixel x,y itself is within the image!
         * Uses class variables width, height: dimensions of the image
         * @param x         x-coordinate of the pixel that has a neighbor in the given direction
         * @param y         y-coordinate of the pixel that has a neighbor in the given direction
         * @param direction the direction from the pixel towards the neighbor (see makeDirectionOffsets)
         * @return          true if the neighbor is within the image (provided that x, y is within)
         */
        boolean isWithin(ImageProcessor ip, int x, int y, int direction) {
            int width = ip.getWidth();
            int height = ip.getHeight();
            int xmax = width - 1;
            int ymax = height -1;
            switch(direction) {
                case 0:
                    return (y>0);
                case 1:
                    return (x<xmax && y>0);
                case 2:
                    return (x<xmax);
                case 3:
                    return (x<xmax && y<ymax);
                case 4:
                    return (y<ymax);
                case 5:
                    return (x>0 && y<ymax);
                case 6:
                    return (x>0);
                case 7:
                    return (x>0 && y>0);
            }
            return false;   //to make the compiler happy :-)
        } // isWithin

        /** Coordinates and the value of a local maximum.
         *  Implemented as a class for easy sorting of an array of maxima by pixel value
         **/
        class MaxPoint implements Comparable {
            float value;
            short x;
            short y;

            /** a constructor filling in the data */
            MaxPoint(short x, short y, float value) {
                this.x = x;
                this.y = y;
                this.value = value;
            }
            /** a comparator required for sorting (interface Comparable) */
            public int compareTo(Object o) {
                //return Float.compare(value,((MaxPoint)o).value); //not possible since this requires Java 1.4
                float diff = value-((MaxPoint)o).value;
                if (diff > 0f) return 1;
                else if (diff == 0f) return 0;
                else return -1;
            }
        }

    } // end of class MaximumFinder

| MaximumFinder |