imageJ.RankFiltersOrbit Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of orbit-image-analysis Show documentation
Show all versions of orbit-image-analysis Show documentation
Orbit, a versatile image analysis software for biological image-based quantification
package imageJ;
import ij.*;
import ij.gui.DialogListener;
import ij.gui.GenericDialog;
import ij.gui.Roi;
import ij.plugin.ContrastEnhancer;
import ij.plugin.filter.ExtendedPlugInFilter;
import ij.plugin.filter.PlugInFilterRunner;
import ij.process.ColorProcessor;
import ij.process.FloatProcessor;
import ij.process.ImageProcessor;
import java.awt.*;
import java.util.Arrays;
/** This plugin implements the Mean, Minimum, Maximum, Variance, Median, Open Maxima, Close Maxima,
* Remove Outliers, Remove NaNs and Despeckle commands.
*/
// Version 2012-07-15 M. Schmid: Fixes a bug that could cause preview not to work correctly
// Version 2012-12-23 M. Schmid: Test for inverted LUT only once (not in each slice)
public class RankFiltersOrbit implements ExtendedPlugInFilter, DialogListener {
public static final int MEAN=0, MIN=1, MAX=2, VARIANCE=3, MEDIAN=4, OUTLIERS=5, DESPECKLE=6, REMOVE_NAN=7,
OPEN=8, CLOSE=9;
private static int HIGHEST_FILTER = CLOSE;
private static final int BRIGHT_OUTLIERS = 0, DARK_OUTLIERS = 1;
private static final String[] outlierStrings = {"Bright","Dark"};
// Filter parameters
private double radius;
private double threshold;
private int whichOutliers;
private int filterType;
// Remember filter parameters for the next time
private static double[] lastRadius = new double[HIGHEST_FILTER+1]; //separate for each filter type
private static double lastThreshold = 50.;
private static int lastWhichOutliers = BRIGHT_OUTLIERS;
//
// F u r t h e r c l a s s v a r i a b l e s
int flags = DOES_ALL|SUPPORTS_MASKING|KEEP_PREVIEW;
private ImagePlus imp;
private int nPasses = 1; // The number of passes (color channels * stack slices)
private PlugInFilterRunner pfr;
private int pass;
// M u l t i t h r e a d i n g - r e l a t e d
private int numThreads = Prefs.getThreads();
// Current state of processing is in class variables. Thus, stack parallelization must be done
// ONLY with one thread for the image (not using these class variables):
private int highestYinCache; // the highest line read into the cache so far
private boolean threadWaiting; // a thread waits until it may read data
private boolean copyingToCache; // whether a thread is currently copying data to the cache
private boolean isMultiStepFilter(int filterType) {
return filterType>=OPEN;
}
/** Setup of the PlugInFilter. Returns the flags specifying the capabilities and needs
* of the filter.
*
* @param arg Defines type of filter operation
* @param imp The ImagePlus to be processed
* @return Flags specifying further action of the PlugInFilterRunner
*/
public int setup(String arg, ImagePlus imp) {
this.imp = imp;
if (arg.equals("mean"))
filterType = MEAN;
else if (arg.equals("min"))
filterType = MIN;
else if (arg.equals("max"))
filterType = MAX;
else if (arg.equals("variance")) {
filterType = VARIANCE;
flags |= FINAL_PROCESSING;
} else if (arg.equals("median"))
filterType = MEDIAN;
else if (arg.equals("outliers"))
filterType = OUTLIERS;
else if (arg.equals("despeckle"))
filterType = DESPECKLE;
else if (arg.equals("close"))
filterType = CLOSE;
else if (arg.equals("open"))
filterType = OPEN;
else if (arg.equals("nan")) {
filterType = REMOVE_NAN;
if (imp!=null && imp.getBitDepth()!=32) {
IJ.error("RankFilters","\"Remove NaNs\" requires a 32-bit image");
return DONE;
}
} else if (arg.equals("final")) { //after variance filter, adjust brightness&contrast
if (imp!=null && imp.getBitDepth()!=8 && imp.getBitDepth()!=24 && imp.getRoi()==null)
new ContrastEnhancer().stretchHistogram(imp.getProcessor(), 0.5);
} else if (arg.equals("masks")) {
showMasks();
return DONE;
} else {
IJ.error("RankFilters","Argument missing or undefined: "+arg);
return DONE;
}
if (isMultiStepFilter(filterType) && imp!=null) { //composite filter: 'open maxima' etc:
Roi roi = imp.getRoi();
if (roi!=null && !roi.getBounds().contains(new Rectangle(imp.getWidth(), imp.getHeight())))
//Roi < image? (actually tested: NOT (Roi>=image))
flags |= SNAPSHOT; //snapshot for resetRoiBoundary
}
return flags;
}
public int showDialog(ImagePlus imp, String command, PlugInFilterRunner pfr) {
if (filterType == DESPECKLE) {
filterType = MEDIAN;
radius = 1.0;
} else {
GenericDialog gd = new GenericDialog(command+"...");
radius = lastRadius[filterType]<=0 ? 2 : lastRadius[filterType];
gd.addNumericField("Radius", radius, 1, 6, "pixels");
int digits = imp.getType() == ImagePlus.GRAY32 ? 2 : 0;
if (filterType==OUTLIERS) {
gd.addNumericField("Threshold", lastThreshold, digits);
gd.addChoice("Which outliers", outlierStrings, outlierStrings[lastWhichOutliers]);
gd.addHelp(IJ.URL+"/docs/menus/process.html#outliers");
} else if (filterType==REMOVE_NAN)
gd.addHelp(IJ.URL+"/docs/menus/process.html#nans");
gd.addPreviewCheckbox(pfr); //passing pfr makes the filter ready for preview
gd.addDialogListener(this); //the DialogItemChanged method will be called on user input
gd.showDialog(); //display the dialog; preview runs in the now
if (gd.wasCanceled()) return DONE;
IJ.register(this.getClass()); //protect static class variables (filter parameters) from garbage collection
if (Macro.getOptions() == null) { //interactive only: remember parameters entered
lastRadius[filterType] = radius;
if (filterType == OUTLIERS) {
lastThreshold = threshold;
lastWhichOutliers = whichOutliers;
}
}
}
this.pfr = pfr;
flags = IJ.setupDialog(imp, flags); //ask whether to process all slices of stack (if a stack)
if ((flags&DOES_STACKS)!=0) {
int size = imp.getWidth() * imp.getHeight();
Roi roi = imp.getRoi();
if (roi != null) {
Rectangle roiRect = roi.getBounds();
size = roiRect.width * roiRect.height;
}
double workToDo = size * radius; //estimate computing time (arb. units)
if (filterType==MEAN || filterType==VARIANCE) workToDo *= 0.5;
else if (filterType==MEDIAN) workToDo *= radius*0.5;
if (workToDo < 1e6 && imp.getImageStackSize()>=numThreads) {
numThreads = 1; //for fast operations, avoid overhead of multi-threading in each image
flags |= PARALLELIZE_STACKS;
}
}
return flags;
}
public boolean dialogItemChanged(GenericDialog gd, AWTEvent e) {
radius = gd.getNextNumber();
if (filterType == OUTLIERS) {
threshold = gd.getNextNumber();
whichOutliers = gd.getNextChoiceIndex();
}
int maxRadius = (filterType==MEDIAN || filterType==OUTLIERS || filterType==REMOVE_NAN) ? 100 : 1000;
return !(gd.invalidNumber() || radius < 0 || radius > maxRadius || (filterType == OUTLIERS && threshold < 0));
}
/**
* Manuel: added to call this filter from external
* @param ip
* @param radius
* @return
*/
public ImagePlus doRemoveOutliers(ImagePlus ip, double radius) {
this.threshold = 0.5f;
setup("outliers",ip);
this.whichOutliers=BRIGHT_OUTLIERS;
makeKernel(radius);
//ip.setProcessor(ip.getProcessor().convertToFloat());
run(ip.getProcessor());
//ip.setProcessor(ip.getProcessor().convertToRGB());
return ip;
}
public void run(ImageProcessor ip) {
rank(ip, radius, filterType, whichOutliers, (float)threshold);
if (IJ.escapePressed()) // interrupted by user?
ip.reset();
}
/** Filters an image by any method except 'despecle' or 'remove outliers'.
* @param ip The ImageProcessor that should be filtered (all 4 types supported)
* @param radius Determines the kernel size, see Process>Filters>Show Circular Masks.
* Must not be negative. No checking is done for large values that would
* lead to excessive computing times.
* @param filterType May be MEAN, MIN, MAX, VARIANCE, or MEDIAN.
*/
public void rank(ImageProcessor ip, double radius, int filterType) {
rank(ip, radius, filterType, 0, 50f);
}
/** Filters an image by any method except 'despecle' (for 'despeckle', use 'median' and radius=1)
* @param ip The image subject to filtering
* @param radius The kernel radius
* @param filterType as defined above; DESPECKLE is not a valid type here; use median and
* a radius of 1.0 instead
* @param whichOutliers BRIGHT_OUTLIERS or DARK_OUTLIERS for 'outliers' filter
* @param threshold Threshold for 'outliers' filter
*/
public void rank(ImageProcessor ip, double radius, int filterType, int whichOutliers, float threshold) {
Rectangle roi = ip.getRoi();
ImageProcessor mask = ip.getMask();
Rectangle roi1 = null;
int[] lineRadii = makeLineRadii(radius);
float minMaxOutliersSign = filterType==MIN ? -1f : 1f;
if (filterType == OUTLIERS) //sign is -1 for high outliers: compare number with minimum
minMaxOutliersSign = (ip.isInvertedLut()==(whichOutliers==DARK_OUTLIERS)) ? -1f : 1f;
boolean isImagePart = (roi.width1 ? 2*numThreads : 0);
// 'cache' is the input buffer. Each line y in the image is mapped onto cache line y%cacheHeight
final float[] cache = new float[cacheWidth*cacheHeight];
highestYinCache = Math.max(roi.y-kHeight/2, 0) - 1; //this line+1 will be read into the cache first
final int[] yForThread = new int[numThreads]; //threads announce here which line they currently process
Arrays.fill(yForThread, -1);
yForThread[numThreads-1] = roi.y-1; //first thread started should begin at roi.y
//IJ.log("going to filter lines "+roi.y+"-"+(roi.y+roi.height-1)+"; cacheHeight="+cacheHeight);
final Thread[] threads = new Thread[numThreads-1]; //thread number 0 is this one, not in the array
for (int t=numThreads-1; t>0; t--) {
final int ti=t;
final Thread thread = new Thread(
new Runnable() {
final public void run() {
doFiltering(ip, lineRadii, cache, cacheWidth, cacheHeight,
filterType, minMaxOutliersSign, threshold, colorChannel,
yForThread, ti, aborted);
}
},
"RankFilters-"+t);
thread.setPriority(Thread.currentThread().getPriority());
thread.start();
threads[ti-1] = thread;
}
doFiltering(ip, lineRadii, cache, cacheWidth, cacheHeight,
filterType, minMaxOutliersSign, threshold, colorChannel,
yForThread, 0, aborted);
for (final Thread thread : threads)
try {
if (thread != null) thread.join();
} catch (InterruptedException e) {
aborted[0] = true;
Thread.currentThread().interrupt(); //keep interrupted status (PlugInFilterRunner needs it)
}
showProgress(1.0, ip instanceof ColorProcessor);
pass++;
}
// Filter a grayscale image or one channel of an RGB image using one thread
//
// Synchronization: unless a thread is waiting, we avoid the overhead of 'synchronized'
// statements. That's because a thread waiting for another one should be rare.
//
// Data handling: The area needed for processing a line is written into the array 'cache'.
// This is a stripe of sufficient width for all threads to have each thread processing one
// line, and some extra space if one thread is finished to start the next line.
// This array is padded at the edges of the image so that a surrounding with radius kRadius
// for each pixel processed is within 'cache'. Out-of-image
// pixels are set to the value of the nearest edge pixel. When adding a new line, the lines in
// 'cache' are not shifted but rather the smaller array with the start and end pointers of the
// kernel area is modified to point at the addresses for the next line.
//
// Algorithm: For mean and variance, except for very small radius, usually do not calculate the
// sum over all pixels. This sum is calculated for the first pixel of every line only. For the
// following pixels, add the new values and subtract those that are not in the sum any more.
// For min/max, also first look at the new values, use their maximum if larger than the old
// one. The look at the values not in the area any more; if it does not contain the old
// maximum, leave the maximum unchanged. Otherwise, determine the maximum inside the area.
// For outliers, calculate the median only if the pixel deviates by more than the threshold
// from any pixel in the area. Therfore min or max is calculated; this is a much faster
// operation than the median.
private void doFiltering(ImageProcessor ip, int[] lineRadii, float[] cache, int cacheWidth, int cacheHeight,
int filterType, float minMaxOutliersSign, float threshold, int colorChannel,
int [] yForThread, int threadNumber, boolean[] aborted) {
if (aborted[0] || Thread.currentThread().isInterrupted()) return;
int width = ip.getWidth();
int height = ip.getHeight();
Rectangle roi = ip.getRoi();
int kHeight = kHeight(lineRadii);
int kRadius = kRadius(lineRadii);
int kNPoints = kNPoints(lineRadii);
int xmin = roi.x - kRadius;
int xmax = roi.x + roi.width + kRadius;
int[]cachePointers = makeCachePointers(lineRadii, cacheWidth);
int padLeft = xmin<0 ? -xmin : 0;
int padRight = xmax>width? xmax-width : 0;
int xminInside = xmin>0 ? xmin : 0;
int xmaxInside = xmax100) {
lastTime = time;
showProgress((y-roi.y)/(double)(roi.height), rgb);
if (Thread.currentThread().isInterrupted() || (imp!= null && IJ.escapePressed())) {
aborted[0] = true;
synchronized(this) {notifyAll();}
return;
}
}
}
for (int i=0; i1) { // thread synchronization
int slowestThreadY = arrayMinNonNegative(yForThread); // non-synchronized check to avoid overhead
if (y - slowestThreadY + kHeight > cacheHeight) { // we would overwrite data needed by another thread
synchronized(this) {
slowestThreadY = arrayMinNonNegative(yForThread); //recheck whether we have to wait
if (y - slowestThreadY + kHeight > cacheHeight) {
do {
notifyAll(); // avoid deadlock: wake up others waiting
threadWaiting = true;
//IJ.log("Thread "+threadNumber+" waiting @y="+y+" slowest@y="+slowestThreadY);
try {
wait();
if (aborted[0]) return;
} catch (InterruptedException e) {
aborted[0] = true;
notifyAll();
Thread.currentThread().interrupt(); //keep interrupted status (PlugInFilterRunner needs it)
return;
}
slowestThreadY = arrayMinNonNegative(yForThread);
} while (y - slowestThreadY + kHeight > cacheHeight);
} //if
threadWaiting = false;
}
}
}
if (numThreads==1) { // R E A D
int yStartReading = y==roi.y ? Math.max(roi.y-kHeight/2, 0) : y+kHeight/2;
for (int yNew = yStartReading; yNew<=y+kHeight/2; yNew++) { //only 1 line except at start
readLineToCacheOrPad(pixels, width, height, roi.y, xminInside, widthInside,
cache, cacheWidth, cacheHeight, padLeft, padRight, colorChannel, kHeight, yNew);
}
} else {
if (!copyingToCache || highestYinCache < y+kHeight/2) synchronized(cache) {
copyingToCache = true; // copy new line(s) into cache
while (highestYinCache < arrayMinNonNegative(yForThread) - kHeight/2 + cacheHeight - 1) {
int yNew = highestYinCache + 1;
readLineToCacheOrPad(pixels, width, height, roi.y, xminInside, widthInside,
cache, cacheWidth, cacheHeight, padLeft, padRight, colorChannel, kHeight, yNew);
highestYinCache = yNew;
}
copyingToCache = false;
}
}
int cacheLineP = cacheWidth * (y % cacheHeight) + kRadius; //points to pixel (roi.x, y)
filterLine(values, width, cache, cachePointers, kNPoints, cacheLineP, roi, y, // F I L T E R
sums, medianBuf1, medianBuf2, minMaxOutliersSign, maxValue, isFloat, filterType,
smallKernel, sumFilter, minOrMax, minOrMaxOrOutliers);
if (!isFloat) //Float images: data are written already during 'filterLine'
writeLineToPixels(values, pixels, roi.x+y*width, roi.width, colorChannel); // W R I T E
//IJ.log("thread "+threadNumber+" @y="+y+" line done");
} // while (!aborted[0]); loop over y (lines)
}
private int arrayMax(int[] array) {
int max = Integer.MIN_VALUE;
for (int i=0; i max) max = array[i];
return max;
}
//returns the minimum of the array, but not less than 0
private int arrayMinNonNegative(int[] array) {
int min = Integer.MAX_VALUE;
for (int i=0; i= max) { //compare with previous maximum 'max'
max = newPointsMax;
} else {
float removedPointsMax = getSideMax(cache, x, cachePointers, false, minMaxOutliersSign);
if (removedPointsMax >= max)
max = getAreaMax(cache, x, cachePointers, 1, newPointsMax, minMaxOutliersSign);
}
if (minOrMax) {
values[valuesP] = max*minMaxOutliersSign;
continue;
}
} else if (sumFilter) {
addSideSums(cache, x, cachePointers, sums);
if (Double.isNaN(sums[0])) //avoid perpetuating NaNs into remaining line
fullCalculation = true;
}
}
if (sumFilter) {
if (filterType == MEAN)
values[valuesP] = (float)(sums[0]/kNPoints);
else {// Variance: sum of squares - square of sums
float value = (float)((sums[1] - sums[0]*sums[0]/kNPoints)/kNPoints);
if (value>maxValue) value = maxValue;
values[valuesP] = value;
}
} else if (filterType == MEDIAN) {
median = getMedian(cache, x, cachePointers, medianBuf1, medianBuf2, kNPoints, median);
values[valuesP] = median;
} else if (filterType == OUTLIERS) {
float v = cache[cacheLineP+x];
if (v*minMaxOutliersSign+threshold < max) { //for low outliers: median can't be higher than max (minMaxOutliersSign is +1)
median = getMedian(cache, x, cachePointers, medianBuf1, medianBuf2, kNPoints, median);
if (v*minMaxOutliersSign+threshold < median*minMaxOutliersSign)
v = median; //beyond threshold (below if minMaxOutliersSign=+1), replace outlier by median
}
values[valuesP] = v;
} else if (filterType == REMOVE_NAN) { //float only; then 'values' is pixels array
if (Float.isNaN(values[valuesP]))
values[valuesP] = getNaNAwareMedian(cache, x, cachePointers, medianBuf1, medianBuf2, kNPoints, median);
else
median = values[valuesP]; //initial guess for the next point
}
} // for x
}
/** Read a line into the cache (including padding in x).
* If y>=height, instead of reading new data, it duplicates the line y=height-1.
* If y==0, it also creates the data for y<0, as far as necessary, thus filling the cache with
* more than one line (padding by duplicating the y=0 row).
*/
private static void readLineToCacheOrPad(Object pixels, int width, int height, int roiY, int xminInside, int widthInside,
float[]cache, int cacheWidth, int cacheHeight, int padLeft, int padRight, int colorChannel,
int kHeight, int y) {
int lineInCache = y%cacheHeight;
if (y < height) {
readLineToCache(pixels, y*width, xminInside, widthInside,
cache, lineInCache*cacheWidth, padLeft, padRight, colorChannel);
if (y==0) for (int prevY = roiY-kHeight/2; prevY<0; prevY++) { //for y<0, pad with y=0 border pixels
int prevLineInCache = cacheHeight+prevY;
System.arraycopy(cache, 0, cache, prevLineInCache*cacheWidth, cacheWidth);
}
} else
System.arraycopy(cache, cacheWidth*((height-1)%cacheHeight), cache, lineInCache*cacheWidth, cacheWidth);
}
/** Read a line into the cache (includes conversion to flaot). Pad with edge pixels in x if necessary */
private static void readLineToCache(Object pixels, int pixelLineP, int xminInside, int widthInside,
float[] cache, int cacheLineP, int padLeft, int padRight, int colorChannel) {
if (pixels instanceof byte[]) {
byte[] bPixels = (byte[])pixels;
for (int pp=pixelLineP+xminInside, cp=cacheLineP+padLeft; pp>shift;
}
for (int cp=cacheLineP; cp guess) {
aboveBuf[nAbove] = v;
nAbove++;
}
else if (v < guess) {
belowBuf[nBelow] = v;
nBelow++;
}
}
}
int half = kNPoints/2;
if (nAbove>half)
return findNthLowestNumber(aboveBuf, nAbove, nAbove-half-1);
else if (nBelow>half)
return findNthLowestNumber(belowBuf, nBelow, half);
else
return guess;
}
/** Get median of values within kernel-sized neighborhood.
* NaN data values are ignored; the output is NaN only if there are only NaN values in the
* kernel-sized neighborhood */
private static float getNaNAwareMedian(float[] cache, int xCache0, int[] kernel,
float[] aboveBuf, float[]belowBuf, int kNPoints, float guess) {
int nAbove = 0, nBelow = 0;
for (int kk=0; kk guess) {
aboveBuf[nAbove] = v;
nAbove++;
}
else if (v < guess) {
belowBuf[nBelow] = v;
nBelow++;
}
}
}
if (kNPoints == 0) return Float.NaN; //only NaN data in the neighborhood?
int half = kNPoints/2;
if (nAbove>half)
return findNthLowestNumber(aboveBuf, nAbove, nAbove-half-1);
else if (nBelow>half)
return findNthLowestNumber(belowBuf, nBelow, half);
else
return guess;
}
/** Find the n-th lowest number in part of an array
* @param buf The input array. Only values 0 ... bufLength are read. buf will be modified.
* @param bufLength Number of values in buf that should be read
* @param n which value should be found; n=0 for the lowest, n=bufLength-1 for the highest
* @return the value */
public final static float findNthLowestNumber(float[] buf, int bufLength, int n) {
// Hoare's find, algorithm, based on http://www.geocities.com/zabrodskyvlada/3alg.html
// Contributed by Heinz Klar
int i,j;
int l=0;
int m=bufLength-1;
float med=buf[n];
float dum ;
while (l=n) && (i<=n)) ;
if (j 0)
System.arraycopy(snapshot, pL, pixels, pL, leftWidth);
if (rightWidth > 0)
System.arraycopy(snapshot, pR, pixels, pR, rightWidth);
}
for (int y=roi.y+roi.height, p = roi1.x+y*width; y=1.5 && radius<1.75) //this code creates the same sizes as the previous RankFilters
radius = 1.75;
else if (radius>=2.5 && radius<2.85)
radius = 2.85;
int r2 = (int) (radius*radius) + 1;
int kRadius = (int)(Math.sqrt(r2+1e-10));
int kHeight = 2*kRadius + 1;
int[] kernel = new int[2*kHeight + 2];
kernel[2*kRadius] = -kRadius;
kernel[2*kRadius+1] = kRadius;
int nPoints = 2*kRadius+1;
for (int y=1; y<=kRadius; y++) { //lines above and below center together
int dx = (int)(Math.sqrt(r2-y*y+1e-10));
kernel[2*(kRadius-y)] = -dx;
kernel[2*(kRadius-y)+1] = dx;
kernel[2*(kRadius+y)] = -dx;
kernel[2*(kRadius+y)+1] = dx;
nPoints += 4*dx+2; //2*dx+1 for each line, above&below
}
kernel[kernel.length-2] = nPoints;
kernel[kernel.length-1] = kRadius;
//for (int i=0; i