| BackgroundSubtracter.java |
package ij.plugin.filter;
import ij.*;
import ij.gui.*;
import ij.process.*;
import ij.measure.*;
import ij.util.Tools;
import java.awt.*;
/** Implements ImageJ's Subtract Background command. Based on the concept of the
rolling ball algorithm described in Stanley Sternberg's article, "Biomedical Image
Processing", IEEE Computer, January 1983.
Imagine that the 2D grayscale image has a third (height) dimension by the image
value at every point in the image, creating a surface. A ball of given radius is
rolled over the bottom side of this surface; the hull of the volume reachable by
the ball is the background.
With "Sliding Parabvoloid", the rolling ball is replaced by a sliding paraboloid
of rotation with the same curvature at its apex as a ball of a given radius.
A paraboloid has the advantage that suitable paraboloids can be found for any image
values, even if the pixel values are much larger than a typical object size (in pixels).
The paraboloid of rotation is approximated as parabolae in 4 directions: x, y and
the two 45-degree directions. Lines of the image in these directions are processed
by sliding a parabola against them. Obtaining the hull needs the parabola for a
given direction to be applied multiple times (after doing the other directions);
in this respect the current code is a compromise between accuracy and speed.
For noise rejection, with the sliding paraboloid algorithm, a 3x3 maximum of the
background is applied. With both, rolling ball and sliding paraboloid,
the image used for calculating the background is slightly smoothened (3x3 average).
This can result in negative values after background subtraction. This preprocessing
can be disabled.
In the sliding paraboloid algorithm, additional code has been added to avoid
subtracting corner objects as a background (note that a paraboloid or ball would
always touch the 4 corner pixels and thus make them background pixels).
This code assumes that corner particles reach less than 1/4 of the image size
into the image.
Rolling ball code based on the NIH Image Pascal version by Michael Castle and Janice
Keller of the University of Michigan Mental Health Research Institute.
Sliding Paraboloid by Michael Schmid, 2007.
Version 10-Jan-2008
*/
public class BackgroundSubtracter implements ExtendedPlugInFilter, DialogListener {
/* parameters from the dialog: */
private static double staticRadius = 50; // default rolling ball radius
private static boolean staticLightBackground = Prefs.get("bs.background", true);
private static boolean staticSeparateColors; // whether to create a separate background for each color channel
private static boolean staticCreateBackground; // don't subtract background (e.g., for processing the background before subtracting)
private static boolean staticUseParaboloid; // use "Sliding Paraboloid" instead of rolling ball algorithm
private static boolean staticDoPresmooth = true; // smoothen the image before creating the background
private double radius = staticRadius;
private boolean lightBackground = staticLightBackground;
private boolean separateColors = staticSeparateColors;
private boolean createBackground = staticCreateBackground;
private boolean useParaboloid = staticUseParaboloid;
private boolean doPresmooth = staticDoPresmooth;
/* more class variables */
private boolean isRGB; // whether we have an RGB image
private boolean previewing;
private final static int MAXIMUM = 0, MEAN = 1; //filter types of filter3x3
private final static int X_DIRECTION = 0, Y_DIRECTION = 1,
DIAGONAL_1A = 2, DIAGONAL_1B = 3, DIAGONAL_2A = 4, DIAGONAL_2B = 5; //filter directions
private final static int DIRECTION_PASSES = 9; //number of passes for different directions
private int nPasses = DIRECTION_PASSES;
private int pass;
private int flags = DOES_ALL|FINAL_PROCESSING|KEEP_PREVIEW|PARALLELIZE_STACKS;
private boolean calledAsPlugin;
public int setup(String arg, ImagePlus imp) {
if (arg.equals("final")) {
imp.getProcessor().resetMinAndMax();
return DONE;
} else
return flags;
}
public int showDialog(ImagePlus imp, String command, PlugInFilterRunner pfr) {
isRGB = imp.getProcessor() instanceof ColorProcessor;
calledAsPlugin = true;
String options = Macro.getOptions();
if (options!=null) { //macro
Macro.setOptions(options.replaceAll("white", "light"));
radius = 50;
lightBackground = false;
separateColors = false;
createBackground = false;
useParaboloid = false;
doPresmooth = true;
}
GenericDialog gd = new GenericDialog(command);
gd.addNumericField("Rolling ball radius:", radius, 1, 6, "pixels");
gd.addCheckbox("Light background", lightBackground);
if (isRGB) gd.addCheckbox("Separate colors", separateColors);
gd.addCheckbox("Create background (don't subtract)", createBackground);
gd.addCheckbox("Sliding paraboloid", useParaboloid);
gd.addCheckbox("Disable smoothing", !doPresmooth);
gd.addPreviewCheckbox(pfr);
gd.addDialogListener(this);
previewing = true;
gd.addHelp(IJ.URL+"/docs/menus/process.html#background");
gd.showDialog();
previewing = false;
if (gd.wasCanceled()) return DONE;
IJ.register(this.getClass()); //protect static class variables (filter parameters) from garbage collection
if ((imp.getProcessor() instanceof FloatProcessor) && !createBackground)
flags |= SNAPSHOT; //FloatProcessors need the original to subtract it from the background
if (options==null) { // not a macro
staticRadius = radius;
staticLightBackground = lightBackground;
staticSeparateColors = separateColors;
staticCreateBackground = createBackground;
staticUseParaboloid = useParaboloid;
staticDoPresmooth = doPresmooth;
Prefs.set("bs.background", lightBackground);
}
return IJ.setupDialog(imp, flags); //ask whether to process all slices of stack (if a stack)
}
public boolean dialogItemChanged(GenericDialog gd, AWTEvent e) {
radius = gd.getNextNumber();
if (radius <= 0.0001 || gd.invalidNumber())
return false;
lightBackground = gd.getNextBoolean();
if (isRGB) separateColors = gd.getNextBoolean();
createBackground = gd.getNextBoolean();
useParaboloid = gd.getNextBoolean();
doPresmooth = !gd.getNextBoolean();
return true;
}
/** Background for any image type */
public void run(ImageProcessor ip) {
if (isRGB && !separateColors)
rollingBallBrightnessBackground((ColorProcessor)ip, radius, createBackground, lightBackground, useParaboloid, doPresmooth, true);
else
rollingBallBackground(ip, radius, createBackground, lightBackground, useParaboloid, doPresmooth, true);
if (previewing && (ip instanceof FloatProcessor || ip instanceof ShortProcessor)) {
ip.resetMinAndMax();
}
}
/** Depracated. For compatibility with previous ImageJ versions */
public void subtractRGBBackround(ColorProcessor ip, int ballRadius) {
rollingBallBrightnessBackground(ip, (double)ballRadius, false, lightBackground, false, true, true);
}
/** Depracated. For compatibility with previous ImageJ versions */
public void subtractBackround(ImageProcessor ip, int ballRadius) {
rollingBallBackground(ip, (double)ballRadius, false, lightBackground, false, true, true);
}
/** Create or subtract a background, based on the brightness of an RGB image (keeping
* the hue of each pixel unchanged)
* @param ip The RGB image. On output, it will become the background-subtracted image or
* the background (depending on <code>createBackground</code>).
* @param radius Radius of the rolling ball creating the background (actually a
* paraboloid of rotation with the same curvature)
* @param createBackground Whether to create a background, not to subtract it.
* @param lightBackground Whether the image has a light background.
* @param doPresmooth Whether the image should be smoothened (3x3 mean) before creating
* the background. With smoothing, the background will not necessarily
* be below the image data.
* @param correctCorners Whether the algorithm should try to detect corner particles to avoid
* subtracting them as a background.
*/
public void rollingBallBrightnessBackground(ColorProcessor ip, double radius, boolean createBackground,
boolean lightBackground, boolean useParaboloid, boolean doPresmooth, boolean correctCorners) {
int width = ip.getWidth();
int height = ip.getHeight();
byte[] H = new byte[width*height];
byte[] S = new byte[width*height];
byte[] B = new byte[width*height];
ip.getHSB(H, S, B);
ByteProcessor bp = new ByteProcessor(width, height, B, null);
rollingBallBackground(bp, radius, createBackground, lightBackground, useParaboloid, doPresmooth, correctCorners);
ip.setHSB(H, S, (byte[])bp.getPixels());
}
/** Create or subtract a background, works for all image types. For RGB images, the
* background is subtracted from each channel separately
* @param ip The image. On output, it will become the background-subtracted image or
* the background (depending on <code>createBackground</code>).
* @param radius Radius of the rolling ball creating the background (actually a
* paraboloid of rotation with the same curvature)
* @param createBackground Whether to create a background, not to subtract it.
* @param lightBackground Whether the image has a light background.
* @param useParaboloid Whether to use the "sliding paraboloid" algorithm.
* @param doPresmooth Whether the image should be smoothened (3x3 mean) before creating
* the background. With smoothing, the background will not necessarily
* be below the image data.
* @param correctCorners Whether the algorithm should try to detect corner particles to avoid
* subtracting them as a background.
*/
public void rollingBallBackground(ImageProcessor ip, double radius, boolean createBackground,
boolean lightBackground, boolean useParaboloid, boolean doPresmooth, boolean correctCorners) {
boolean invertedLut = ip.isInvertedLut();
boolean invert = (invertedLut && !lightBackground) || (!invertedLut && lightBackground);
RollingBall ball = null;
if (!useParaboloid) ball = new RollingBall(radius);
FloatProcessor fp = null;
for (int channelNumber=0; channelNumber<ip.getNChannels(); channelNumber++) {
fp = ip.toFloat(channelNumber, fp);
if ((ip instanceof FloatProcessor) && !calledAsPlugin && !createBackground)
fp.snapshot(); //float images need a snapshot to subtract
if (useParaboloid)
slidingParaboloidFloatBackground(fp, (float)radius, invert, doPresmooth, correctCorners);
else
rollingBallFloatBackground(fp, (float)radius, invert, doPresmooth, ball);
if (createBackground)
ip.setPixels(channelNumber, fp);
else { //subtract the background now
float[] bgPixels = (float[])fp.getPixels(); //currently holds the background
if (ip instanceof FloatProcessor) { //here ip and fp are the same (bgPixels will be output)
float[] snapshotPixels = (float[])fp.getSnapshotPixels(); //original data in the snapshot
for (int p=0; p<bgPixels.length; p++)
bgPixels[p] = snapshotPixels[p]-bgPixels[p];
//for all others, the image data are in ip, the background is in fp
} else if (ip instanceof ShortProcessor) {
float offset = invert ? 65535.5f : 0.5f;//includes 0.5 for rounding when converting float to short
short[] pixels = (short[])ip.getPixels();
for (int p=0; p<bgPixels.length; p++) {
float value = (pixels[p]&0xffff) - bgPixels[p] + offset;
if (value<0f) value = 0f;
if (value>65535f) value = 65535f;
pixels[p] = (short)(value);
}
} else if (ip instanceof ByteProcessor) {
float offset = invert ? 255.5f : 0.5f; //includes 0.5 for rounding when converting float to byte
byte[] pixels = (byte[])ip.getPixels();
for (int p=0; p<bgPixels.length; p++) {
float value = (pixels[p]&0xff) - bgPixels[p] + offset;
if (value<0f) value = 0f;
if (value>255f) value = 255f;
pixels[p] = (byte)(value);
}
} else if (ip instanceof ColorProcessor) {
float offset = invert ? 255.5f : 0.5f;
int[] pixels = (int[])ip.getPixels();
int shift = 16 - 8*channelNumber;
int byteMask = 255<<shift;
int resetMask = 0xffffffff^(255<<shift);
for (int p=0; p<bgPixels.length; p++) {
int pxl = pixels[p];
float value = ((pxl&byteMask)>>shift) - bgPixels[p] + offset;
if (value<0f) value = 0f;
if (value>255f) value = 255f;
pixels[p] = (pxl&resetMask) | ((int)value<<shift);
}
}
}
}
}
// S L I D E P A R A B O L O I D S E C T I O N
/** Create background for a float image by sliding a paraboloid over
* the image. */
void slidingParaboloidFloatBackground(FloatProcessor fp, float radius, boolean invert,
boolean doPresmooth, boolean correctCorners) {
float[] pixels = (float[])fp.getPixels(); //this will become the background
int width = fp.getWidth();
int height = fp.getHeight();
float[] cache = new float[Math.max(width, height)]; //work array for lineSlideParabola
int[] nextPoint = new int[Math.max(width, height)]; //work array for lineSlideParabola
float coeff2 = 0.5f/radius; //2nd-order coefficient of the polynomial approximating the ball
float coeff2diag = 1.f/radius; //same for diagonal directions where step is sqrt2
showProgress(0.000001); //start the progress bar (only filter1D will increment it)
if (invert)
for (int i=0; i<pixels.length; i++)
pixels[i] = -pixels[i];
float shiftBy = 0;
if (doPresmooth) {
shiftBy = (float)filter3x3(fp, MAXIMUM);//3x3 maximum to remove dust etc.
showProgress(0.5);
filter3x3(fp, MEAN); //smoothing to remove noise
pass++;
}
if (correctCorners)
correctCorners(fp, coeff2, cache, nextPoint); //modify corner data, avoids subtracting corner particles
/* Slide the parabola over the image in different directions */
/* Doing the diagonal directions at the end is faster (diagonal lines are denser,
* so there are more such lines, and the algorithm gets faster with each iteration) */
filter1D(fp, X_DIRECTION, coeff2, cache, nextPoint);
filter1D(fp, Y_DIRECTION, coeff2, cache, nextPoint);
filter1D(fp, X_DIRECTION, coeff2, cache, nextPoint); //redo for better accuracy
filter1D(fp, DIAGONAL_1A, coeff2diag, cache, nextPoint);
filter1D(fp, DIAGONAL_1B, coeff2diag, cache, nextPoint);
filter1D(fp, DIAGONAL_2A, coeff2diag, cache, nextPoint);
filter1D(fp, DIAGONAL_2B, coeff2diag, cache, nextPoint);
filter1D(fp, DIAGONAL_1A, coeff2diag, cache, nextPoint);//redo for better accuracy
filter1D(fp, DIAGONAL_1B, coeff2diag, cache, nextPoint);
if (invert)
for (int i=0; i<pixels.length; i++)
pixels[i] = -(pixels[i] - shiftBy);
else if (doPresmooth)
for (int i=0; i<pixels.length; i++)
pixels[i] -= shiftBy; //correct for shift by 3x3 maximum
}
/** Filter by subtracting a sliding parabola for all lines in one direction, x, y or one of
* the two diagonal directions (diagonals are processed only for half the image per call). */
void filter1D(FloatProcessor fp, int direction, float coeff2, float[] cache, int[] nextPoint) {
float[] pixels = (float[])fp.getPixels(); //this will become the background
int width = fp.getWidth();
int height = fp.getHeight();
int startLine = 0; //index of the first line to handle
int nLines = 0; //index+1 of the last line to handle (initialized to avoid compile-time error)
int lineInc = 0; //increment from one line to the next in pixels array
int pointInc = 0; //increment from one point to the next along the line
int length = 0; //length of the line
switch (direction) {
case X_DIRECTION: //lines parallel to x direction
nLines = height;
lineInc = width;
pointInc = 1;
length = width;
break;
case Y_DIRECTION: //lines parallel to y direction
nLines = width;
lineInc = 1;
pointInc = width;
length = height;
break;
case DIAGONAL_1A: //lines parallel to x=y, starting at x axis
nLines = width-2; //the algorithm makes no sense for lines shorter than 3 pixels
lineInc = 1;
pointInc = width + 1;
break;
case DIAGONAL_1B: //lines parallel to x=y, starting at y axis
startLine = 1;
nLines = height-2;
lineInc = width;
pointInc = width + 1;
break;
case DIAGONAL_2A: //lines parallel to x=-y, starting at x axis
startLine = 2;
nLines = width;
lineInc = 1;
pointInc = width - 1;
break;
case DIAGONAL_2B: //lines parallel to x=-y, starting at x=width-1, y=variable
startLine = 0;
nLines = height-2;
lineInc = width;
pointInc = width - 1;
break;
}
for (int i=startLine; i<nLines; i++) {
if (i%50==0) {
if (Thread.currentThread().isInterrupted()) return;
showProgress(i/(double)nLines);
}
int startPixel = i*lineInc;
if (direction == DIAGONAL_2B) startPixel += width-1;
switch (direction) {
case DIAGONAL_1A: length = Math.min(height, width-i); break;
case DIAGONAL_1B: length = Math.min(width, height-i); break;
case DIAGONAL_2A: length = Math.min(height, i+1); break;
case DIAGONAL_2B: length = Math.min(width, height-i); break;
}
lineSlideParabola(pixels, startPixel, pointInc, length, coeff2, cache, nextPoint, null);
}
pass++;
} //void filter1D
/** Process one straight line in the image by sliding a parabola along the line
* (from the bottom) and setting the values to make all points reachable by
* the parabola
* @param pixels Image data, will be modified by parabolic interpolation
* where the parabola does not touch.
* @param start Index of first pixel of the line in pixels array
* @param inc Increment of index in pixels array
* @param length Number of points the line consists of
* @param coeff2 2nd order coefficient of the polynomial describing the parabola,
* must be positive (although a parabola with negative curvature is
* actually used)
* @param cache Work array, length at least <code>length</code>. Will usually remain
* in the CPU cache and may therefore speed up the code.
* @param nextPoint Work array. Will hold the index of the next point with sufficient local
* curvature to get touched by the parabola.
* @param correctedEdges Should be a 2-element array used for output or null.
* @return The correctedEdges array (if non-null on input) with the two estimated
* edge pixel values corrected for edge particles.
*/
static float[] lineSlideParabola(float[] pixels, int start, int inc, int length, float coeff2, float[] cache, int[] nextPoint, float[] correctedEdges) {
float minValue = Float.MAX_VALUE;
int lastpoint = 0;
int firstCorner = length-1; // the first point except the edge that is touched
int lastCorner = 0; // the last point except the edge that is touched
float vPrevious1 = 0f;
float vPrevious2 = 0f;
float curvatureTest = 1.999f*coeff2; //not 2: numeric scatter of 2nd derivative
/* copy data to cache, determine the minimum, and find points with local curvature such
* that the parabola can touch them - only these need to be examined futher on */
for (int i=0, p=start; i<length; i++, p+=inc) {
float v = pixels[p];
cache[i] = v;
if (v < minValue) minValue = v;
if (i >= 2 && vPrevious1+vPrevious1-vPrevious2-v < curvatureTest) {
nextPoint[lastpoint] = i-1; // point i-1 may be touched
lastpoint = i-1;
}
vPrevious2 = vPrevious1;
vPrevious1 = v;
}
nextPoint[lastpoint] = length-1;
nextPoint[length-1] = Integer.MAX_VALUE;// breaks the search loop
int i1 = 0; // i1 and i2 will be the two points where the parabola touches
while (i1<length-1) {
float v1 = cache[i1];
float minSlope = Float.MAX_VALUE;
int i2 = 0; //(initialized to avoid compile-time error)
int searchTo = length;
int recalculateLimitNow = 0; // when 0, limits for searching will be recalculated
/* find the second point where the parabola through point i1,v1 touches: */
for (int j=nextPoint[i1]; j<searchTo; j=nextPoint[j], recalculateLimitNow++) {
float v2 = cache[j];
float slope = (v2-v1)/(j-i1)+coeff2*(j-i1);
if (slope < minSlope) {
minSlope = slope;
i2 = j;
recalculateLimitNow = -3;
}
if (recalculateLimitNow==0) { //time-consuming recalculation of search limit: wait a bit after slope is updated
double b = 0.5f*minSlope/coeff2;
int maxSearch = i1+(int)(b+Math.sqrt(b*b+(v1-minValue)/coeff2)+1); //(numeric overflow may make this negative)
if (maxSearch < searchTo && maxSearch > 0) searchTo = maxSearch;
}
}
if (i1 == 0) firstCorner = i2;
if (i2 == length-1) lastCorner = i1;
/* interpolate between the two points where the parabola touches: */
for (int j=i1+1, p=start+j*inc; j<i2; j++, p+=inc)
pixels[p] = v1 + (j-i1)*(minSlope - (j-i1)*coeff2);
i1 = i2; // continue from this new point
} //while (i1<length-1)
/* Now calculate estimated edge values without an edge particle, allowing for vignetting
* described as a 6th-order polynomial: */
if (correctedEdges != null) {
if (4*firstCorner >= length) firstCorner = 0; // edge particles must be < 1/4 image size
if (4*(length - 1 - lastCorner) >= length) lastCorner = length - 1;
float v1 = cache[firstCorner];
float v2 = cache[lastCorner];
float slope = (v2-v1)/(lastCorner-firstCorner); // of the line through the two outermost non-edge touching points
float value0 = v1 - slope * firstCorner; // offset of this line
float coeff6 = 0; // coefficient of 6th order polynomial
float mid = 0.5f * (lastCorner + firstCorner);
for (int i=(length+2)/3; i<=(2*length)/3; i++) {// compare with mid-image pixels to detect vignetting
float dx = (i-mid)*2f/(lastCorner-firstCorner);
float poly6 = dx*dx*dx*dx*dx*dx - 1f; // the 6th order polynomial, zero at firstCorner and lastCorner
if (cache[i] < value0 + slope*i + coeff6*poly6) {
coeff6 = -(value0 + slope*i - cache[i])/poly6;
}
}
float dx = (firstCorner-mid)*2f/(lastCorner-firstCorner);
correctedEdges[0] = value0 + coeff6*(dx*dx*dx*dx*dx*dx - 1f) + coeff2*firstCorner*firstCorner;
dx = (lastCorner-mid)*2f/(lastCorner-firstCorner);
correctedEdges[1] = value0 + (length-1)*slope + coeff6*(dx*dx*dx*dx*dx*dx - 1f) + coeff2*(length-1-lastCorner)*(length-1-lastCorner);
}
return correctedEdges;
} //void lineSlideParabola
/** Detect corner particles and adjust corner pixels if a particle is there.
* Analyzing the directions parallel to the edges and the diagonals, we
* average over the 3 correction values (found for the 3 directions)
*/
void correctCorners(FloatProcessor fp, float coeff2, float[] cache, int[] nextPoint) {
int width = fp.getWidth();
int height = fp.getHeight();
float[] pixels = (float[])fp.getPixels();
float[] corners = new float[4]; //(0,0); (xmax,0); (ymax,0); (xmax,ymax)
float[] correctedEdges = new float[2];
correctedEdges = lineSlideParabola(pixels, 0, 1, width, coeff2, cache, nextPoint, correctedEdges);
corners[0] = correctedEdges[0];
corners[1] = correctedEdges[1];
correctedEdges = lineSlideParabola(pixels, (height-1)*width, 1, width, coeff2, cache, nextPoint, correctedEdges);
corners[2] = correctedEdges[0];
corners[3] = correctedEdges[1];
correctedEdges = lineSlideParabola(pixels, 0, width, height, coeff2, cache, nextPoint, correctedEdges);
corners[0] += correctedEdges[0];
corners[2] += correctedEdges[1];
correctedEdges = lineSlideParabola(pixels, width-1, width, height, coeff2, cache, nextPoint, correctedEdges);
corners[1] += correctedEdges[0];
corners[3] += correctedEdges[1];
int diagLength = Math.min(width,height); //length of a 45-degree line from a corner
float coeff2diag = 2 * coeff2;
correctedEdges = lineSlideParabola(pixels, 0, 1+width, diagLength, coeff2diag, cache, nextPoint, correctedEdges);
corners[0] += correctedEdges[0];
correctedEdges = lineSlideParabola(pixels, width-1, -1+width, diagLength, coeff2diag, cache, nextPoint, correctedEdges);
corners[1] += correctedEdges[0];
correctedEdges = lineSlideParabola(pixels, (height-1)*width, 1-width, diagLength, coeff2diag, cache, nextPoint, correctedEdges);
corners[2] += correctedEdges[0];
correctedEdges = lineSlideParabola(pixels, width*height-1, -1-width, diagLength, coeff2diag, cache, nextPoint, correctedEdges);
corners[3] += correctedEdges[0];
if (pixels[0] > corners[0]/3) pixels[0] = corners[0]/3;
if (pixels[width-1] > corners[1]/3) pixels[width-1] = corners[1]/3;
if (pixels[(height-1)*width] > corners[2]/3) pixels[(height-1)*width] = corners[2]/3;
if (pixels[width*height-1] > corners[3]/3) pixels[width*height-1] = corners[3]/3;
//new ImagePlus("corner corrected",fp.duplicate()).show();
} //void correctCorners
// R O L L B A L L S E C T I O N
/** Create background for a float image by rolling a ball over
* the image. */
void rollingBallFloatBackground(FloatProcessor fp, float radius, boolean invert,
boolean doPresmooth, RollingBall ball) {
float[] pixels = (float[])fp.getPixels(); //this will become the background
boolean shrink = ball.shrinkFactor >1;
showProgress(0.0);
if (invert)
for (int i=0; i<pixels.length; i++)
pixels[i] = -pixels[i];
if (doPresmooth)
filter3x3(fp, MEAN);
double[] minmax = Tools.getMinMax(pixels);
if (Thread.currentThread().isInterrupted()) return;
FloatProcessor smallImage = shrink ? shrinkImage(fp, ball.shrinkFactor) : fp;
if (Thread.currentThread().isInterrupted()) return;
rollBall(ball, smallImage);
if (Thread.currentThread().isInterrupted()) return;
showProgress(0.9);
if (shrink)
enlargeImage(smallImage, fp, ball.shrinkFactor);
if (Thread.currentThread().isInterrupted()) return;
if (invert)
for (int i=0; i<pixels.length; i++)
pixels[i] = -pixels[i];
pass++;
}
/** Creates a lower resolution image for ball-rolling. */
FloatProcessor shrinkImage(FloatProcessor ip, int shrinkFactor) {
int width = ip.getWidth();
int height = ip.getHeight();
float[] pixels = (float[])ip.getPixels();
int sWidth = (width+shrinkFactor-1)/shrinkFactor;
int sHeight = (height+shrinkFactor-1)/shrinkFactor;
showProgress(0.1);
FloatProcessor smallImage = new FloatProcessor(sWidth, sHeight);
float[] sPixels = (float[])smallImage.getPixels();
float min, thispixel;
for (int ySmall=0; ySmall<sHeight; ySmall++) {
for (int xSmall=0; xSmall<sWidth; xSmall++) {
min = Float.MAX_VALUE;
for (int j=0, y=shrinkFactor*ySmall; j<shrinkFactor&&y<height; j++, y++) {
for (int k=0, x=shrinkFactor*xSmall; k<shrinkFactor&&x<width; k++, x++) {
thispixel = pixels[x+y*width];
if (thispixel<min)
min = thispixel;
}
}
sPixels[xSmall+ySmall*sWidth] = min; // each point in small image is minimum of its neighborhood
}
}
//new ImagePlus("smallImage", smallImage).show();
return smallImage;
}
/** 'Rolls' a filtering object over a (shrunken) image in order to find the
image's smooth continuous background. For the purpose of explaining this
algorithm, imagine that the 2D grayscale image has a third (height) dimension
defined by the intensity value at every point in the image. The center of
the filtering object, a patch from the top of a sphere having radius BallRadius,
is moved along each scan line of the image so that the patch is tangent to the
image at one or more points with every other point on the patch below the
corresponding (x,y) point of the image. Any point either on or below the patch
during this process is considered part of the background. Shrinking the image
before running this procedure is advised for large ball radii because the
processing time increases with ball radius^2.
*/
void rollBall(RollingBall ball, FloatProcessor fp) {
float[] pixels = (float[])fp.getPixels(); //the input pixels
int width = fp.getWidth();
int height = fp.getHeight();
float[] zBall = ball.data;
int ballWidth = ball.width;
int radius = ballWidth/2;
float[] cache = new float[width*ballWidth]; //temporarily stores the pixels we work on
Thread thread = Thread.currentThread();
long lastTime = System.currentTimeMillis();
for (int y=-radius; y<height+radius; y++) { //for all positions of the ball center:
long time = System.currentTimeMillis();
if (time-lastTime > 100) {
lastTime = time;
if (thread.isInterrupted()) return;
showProgress(0.1+0.8*y/(height+ballWidth));
}
int nextLineToWriteInCache = (y+radius)%ballWidth;
int nextLineToRead = y + radius; //line of the input not touched yet
if (nextLineToRead<height) {
System.arraycopy(pixels, nextLineToRead*width, cache, nextLineToWriteInCache*width, width);
for (int x=0, p=nextLineToRead*width; x<width; x++,p++)
pixels[p] = -Float.MAX_VALUE; //unprocessed pixels start at minus infinity
}
int y0 = y-radius; //the first line to see whether the ball touches
if (y0 < 0) y0 = 0;
int yBall0 = y0-y+radius; //y coordinate in the ball corresponding to y0
int yend = y+radius; //the last line to see whether the ball touches
if (yend>=height) yend = height-1;
for (int x=-radius; x<width+radius; x++) {
float z = Float.MAX_VALUE; //the height of the ball (ball is in position x,y)
int x0 = x-radius;
if (x0 < 0) x0 = 0;
int xBall0 = x0-x+radius;
int xend = x+radius;
if (xend>=width) xend = width-1;
for (int yp=y0, yBall=yBall0; yp<=yend; yp++,yBall++) { //for all points inside the ball
int cachePointer = (yp%ballWidth)*width+x0;
for (int xp=x0, bp=xBall0+yBall*ballWidth; xp<=xend; xp++, cachePointer++, bp++) {
float zReduced = cache[cachePointer] - zBall[bp];
if (z > zReduced) //does this point imply a greater height?
z = zReduced;
}
}
for (int yp=y0, yBall=yBall0; yp<=yend; yp++,yBall++) //raise pixels to ball surface
for (int xp=x0, p=xp+yp*width, bp=xBall0+yBall*ballWidth; xp<=xend; xp++, p++, bp++) {
float zMin = z + zBall[bp];
if (pixels[p] < zMin)
pixels[p] = zMin;
}
// if (x>=0&&y>=0&&x<width&&y<height) bgPixels[x+y*width] = z; //debug, ball height output
}
}
//new ImagePlus("bg rolled", fp.duplicate()).show();
}
/** Uses bilinear interpolation to find the points in the full-scale background
given the points from the shrunken image background. (At the edges, it is
actually extrapolation.)
*/
void enlargeImage(FloatProcessor smallImage, FloatProcessor fp, int shrinkFactor) {
int width = fp.getWidth();
int height = fp.getHeight();
int smallWidth = smallImage.getWidth();
int smallHeight = smallImage.getHeight();
float[] pixels = (float[])fp.getPixels();
float[] sPixels = (float[])smallImage.getPixels();
int[] xSmallIndices = new int[width]; //index of first point in smallImage
float[] xWeights = new float[width]; //weight of this point
makeInterpolationArrays(xSmallIndices, xWeights, width, smallWidth, shrinkFactor);
int[] ySmallIndices = new int[height];
float[] yWeights = new float[height];
makeInterpolationArrays(ySmallIndices, yWeights, height, smallHeight, shrinkFactor);
float[] line0 = new float[width];
float[] line1 = new float[width];
for (int x=0; x<width; x++) //x-interpolation of the first smallImage line
line1[x] = sPixels[xSmallIndices[x]] * xWeights[x] +
sPixels[xSmallIndices[x]+1] * (1f - xWeights[x]);
int ySmallLine0 = -1; //line0 corresponds to this y of smallImage
for (int y=0; y<height; y++) {
if (ySmallLine0 < ySmallIndices[y]) {
float[] swap = line0; //previous line1 -> line0
line0 = line1;
line1 = swap; //keep the other array for filling with new data
ySmallLine0++;
int sYPointer = (ySmallIndices[y]+1)*smallWidth; //points to line0 + 1 in smallImage
for (int x=0; x<width; x++) //x-interpolation of the new smallImage line -> line1
line1[x] = sPixels[sYPointer+xSmallIndices[x]] * xWeights[x] +
sPixels[sYPointer+xSmallIndices[x]+1] * (1f - xWeights[x]);
}
float weight = yWeights[y];
for (int x=0, p=y*width; x<width; x++,p++)
pixels[p] = line0[x]*weight + line1[x]*(1f - weight);
}
}
/** Create arrays of indices and weigths for interpolation.
<pre>
Example for shrinkFactor = 4:
small image pixel number | 0 | 1 | 2 | ...
full image pixel number | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 |10 |11 | ...
smallIndex for interpolation(0) | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 2 | 2 | ...
(0) Note: This is smallIndex for the left pixel; for the right pixel used for interpolation
it is higher by one
</pre>
*/
void makeInterpolationArrays(int[] smallIndices, float[] weights, int length, int smallLength, int shrinkFactor) {
for (int i=0; i<length; i++) {
int smallIndex = (i - shrinkFactor/2)/shrinkFactor;
if (smallIndex >= smallLength-1) smallIndex = smallLength - 2;
smallIndices[i] = smallIndex;
float distance = (i + 0.5f)/shrinkFactor - (smallIndex + 0.5f); //distance of pixel centers (in smallImage pixels)
weights[i] = 1f - distance;
}
}
// C O M M O N S E C T I O N F O R B O T H A L G O R I T H M S
/** Replace the pixels by the mean or maximum in a 3x3 neighborhood.
* No snapshot is required (less memory needed than e.g., fp.smooth()).
* When used as maximum filter, it returns the average change of the
* pixel value by this operation
*/
double filter3x3(FloatProcessor fp, int type) {
int width = fp.getWidth();
int height = fp.getHeight();
double shiftBy = 0;
float[] pixels = (float[])fp.getPixels();
for (int y=0; y<height; y++)
shiftBy += filter3(pixels, width, y*width, 1, type);
for (int x=0; x<width; x++)
shiftBy += filter3(pixels, height, x, width, type);
return shiftBy/width/height;
}
/** Filter a line: maximum or average of 3-pixel neighborhood */
double filter3(float[] pixels, int length, int pixel0, int inc, int type) {
double shiftBy = 0;
float v3 = pixels[pixel0]; //will be pixel[i+1]
float v2 = v3; //will be pixel[i]
float v1; //will be pixel[i-1]
for (int i=0, p=pixel0; i<length; i++,p+=inc) {
v1 = v2;
v2 = v3;
if (i<length-1) v3 = pixels[p+inc];
if (type == MAXIMUM) {
float max = v1 > v3 ? v1 : v3;
if (v2 > max) max = v2;
shiftBy += max - v2;
pixels[p] = max;
} else
pixels[p] = (v1+v2+v3)*0.33333333f;
}
return shiftBy;
}
public void setNPasses(int nPasses) {
if (isRGB && separateColors)
nPasses *= 3;
if (useParaboloid)
nPasses*= (doPresmooth) ? DIRECTION_PASSES+2 : DIRECTION_PASSES;
this.nPasses = nPasses;
pass = 0;
}
private void showProgress(double percent) {
if (nPasses <= 0) return;
percent = (double)pass/nPasses + percent/nPasses;
IJ.showProgress(percent);
}
}
// C L A S S R O L L I N G B A L L
/** A rolling ball (or actually a square part thereof)
* Here it is also determined whether to shrink the image
*/
class RollingBall {
float[] data;
int width;
int shrinkFactor;
RollingBall(double radius) {
int arcTrimPer;
if (radius<=10) {
shrinkFactor = 1;
arcTrimPer = 24; // trim 24% in x and y
} else if (radius<=30) {
shrinkFactor = 2;
arcTrimPer = 24; // trim 24% in x and y
} else if (radius<=100) {
shrinkFactor = 4;
arcTrimPer = 32; // trim 32% in x and y
} else {
shrinkFactor = 8;
arcTrimPer = 40; // trim 40% in x and y
}
buildRollingBall(radius, arcTrimPer);
}
/** Computes the location of each point on the rolling ball patch relative to the
center of the sphere containing it. The patch is located in the top half
of this sphere. The vertical axis of the sphere passes through the center of
the patch. The projection of the patch in the xy-plane below is a square.
*/
void buildRollingBall(double ballradius, int arcTrimPer) {
double rsquare; // rolling ball radius squared
int xtrim; // # of pixels trimmed off each end of ball to make patch
int xval, yval; // x,y-values on patch relative to center of rolling ball
double smallballradius; // radius of rolling ball (downscaled in x,y and z when image is shrunk)
int halfWidth; // distance in x or y from center of patch to any edge (patch "radius")
this.shrinkFactor = shrinkFactor;
smallballradius = ballradius/shrinkFactor;
if (smallballradius<1)
smallballradius = 1;
rsquare = smallballradius*smallballradius;
xtrim = (int)(arcTrimPer*smallballradius)/100; // only use a patch of the rolling ball
halfWidth = (int)Math.round(smallballradius - xtrim);
width = 2*halfWidth+1;
data = new float[width*width];
for (int y=0, p=0; y<width; y++)
for (int x=0; x<width; x++, p++) {
xval = x - halfWidth;
yval = y - halfWidth;
double temp = rsquare - xval*xval - yval*yval;
data[p] = temp>0. ? (float)(Math.sqrt(temp)) : 0f;
//-Float.MAX_VALUE might be better than 0f, but gives different results than earlier versions
}
}
}