|
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 |
|