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.URL2+"/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 createBackground
).
* @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 createBackground
).
* @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; channelNumber65535f) 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; p255f) 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) - bgPixels[p] + offset;
if (value<0f) value = 0f;
if (value>255f) value = 255f;
pixels[p] = (pxl&resetMask) | ((int)value<length. 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= 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 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= 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 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) yend = height-1;
for (int x=-radius; x=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 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 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
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
*/
void makeInterpolationArrays(int[] smallIndices, float[] weights, int length, int smallLength, int shrinkFactor) {
for (int i=0; i= 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 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; y0. ? (float)(Math.sqrt(temp)) : 0f;
//-Float.MAX_VALUE might be better than 0f, but gives different results than earlier versions
}
}
}