|
EDM |
|
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() }
|
EDM |
|