package ij.plugin.filter;
import ij.*;
import ij.gui.GenericDialog;
import ij.gui.DialogListener;
import ij.gui.Roi;
import ij.process.*;
import ij.plugin.ContrastEnhancer;
import java.awt.*;
import java.awt.event.*;
import java.util.Arrays;
/** This plugin implements the Mean, Minimum, Maximum, Variance, Median, Open Maxima, Close Maxima,
* Remove Outliers, Remove NaNs and Despeckle commands.
*/
public class RankFilters 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*(double)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;
if (gd.invalidNumber() || radius<0 || radius>maxRadius || (filterType==OUTLIERS && threshold <0))
return false;
return true;
}
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);
boolean isImagePart = (roi.width<ip.getWidth()) || (roi.height<ip.getHeight());
boolean[] aborted = new boolean[1]; // returns whether interrupted during preview or ESC pressed
for (int ch=0; ch<ip.getNChannels(); ch++) {
int filterType1 = filterType;
if (isMultiStepFilter(filterType)) {
filterType1 = (filterType==OPEN) ? MIN : MAX;
if (isImagePart) { //composite filters ('open maxima' etc.) need larger area in first step
int kRadius = kRadius(lineRadii);
int kHeight = kHeight(lineRadii);
Rectangle roiClone = (Rectangle)roi.clone();
roiClone.grow(kRadius, kHeight/2);
roi1 = roiClone.intersection(new Rectangle(ip.getWidth(), ip.getHeight()));
ip.setRoi(roi1);
}
}
doFiltering(ip, lineRadii, filterType1, whichOutliers, threshold, ch, aborted);
if (aborted[0]) break;
if (isMultiStepFilter(filterType)) {
ip.setRoi(roi);
ip.setMask(mask);
int filterType2 = (filterType==OPEN) ? MAX : MIN;
doFiltering(ip, lineRadii, filterType2, whichOutliers, threshold, ch, aborted);
if (aborted[0]) break;
if (isImagePart)
resetRoiBoundary(ip, roi, roi1);
}
}
}
// Filter a grayscale image or one channel of an RGB image with several threads
// Implementation: each thread uses the same input buffer (cache), always works on the next unfiltered line
// Usually, one thread reads reads several lines into the cache, while the others are processing the data.
// 'aborted[0]' is set if the main thread has been interrupted (during preview) or ESC pressed.
// 'aborted' must not be a class variable because it signals the other threads to stop; and this may be caused
// by an interrupted preview thread after the main calculation has been started.
private void doFiltering(final ImageProcessor ip, final int[] lineRadii, final int filterType,
final int whichOutliers, final float threshold, final int colorChannel, final boolean[] aborted) {
Rectangle roi = ip.getRoi();
int width = ip.getWidth();
Object pixels = ip.getPixels();
int numThreads = Math.min(roi.height, this.numThreads);
int kHeight = kHeight(lineRadii);
int kRadius = kRadius(lineRadii);
final int cacheWidth = roi.width+2*kRadius;
final int cacheHeight = kHeight + (numThreads>1 ? 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
final Thread[] threads = new Thread[numThreads];
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, whichOutliers, threshold, colorChannel,
yForThread, ti, aborted);
}
},
"RankFilters-"+t);
thread.setPriority(Thread.currentThread().getPriority());
thread.start();
threads[ti] = thread;
}
doFiltering(ip, lineRadii, cache, cacheWidth, cacheHeight,
filterType, whichOutliers, threshold, colorChannel,
yForThread, 0, aborted);
try {
for (final Thread thread : threads)
if (thread != null) thread.join();
}
catch (InterruptedException e) {
aborted[0] = true;
}
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, int whichOutliers, 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 = xmax<width ? xmax : width;
int widthInside = xmaxInside - xminInside;
boolean minOrMax = filterType == MIN || filterType == MAX;
boolean minOrMaxOrOutliers = minOrMax || filterType == OUTLIERS;
boolean sumFilter = filterType == MEAN || filterType == VARIANCE;
boolean medianFilter = filterType == MEDIAN || filterType == OUTLIERS;
double[] sums = sumFilter ? new double[2] : null;
float[] medianBuf1 = (medianFilter||filterType==REMOVE_NAN) ? new float[kNPoints] : null;
float[] medianBuf2 = (medianFilter||filterType==REMOVE_NAN) ? new float[kNPoints] : null;
float sign = filterType==MIN ? -1f : 1f;
if (filterType == OUTLIERS) //sign is -1 for high outliers: compare number with minimum
sign = (ip.isInvertedLut()==(whichOutliers==DARK_OUTLIERS)) ? -1f : 1f;
boolean smallKernel = kRadius < 2;
Object pixels = ip.getPixels();
boolean isFloat = pixels instanceof float[];
float maxValue = isFloat ? Float.NaN : (float)ip.maxValue();
float[] values = isFloat ? (float[])pixels : new float[roi.width];
int numThreads = yForThread.length;
long lastTime = System.currentTimeMillis();
int previousY = kHeight/2-cacheHeight;
boolean rgb = ip instanceof ColorProcessor;
while (!aborted[0]) {
int y = arrayMax(yForThread) + 1; // y of the next line that needs processing
yForThread[threadNumber] = y;
//IJ.log("thread "+threadNumber+" @y="+y+" needs"+(y-kHeight/2)+"-"+(y+kHeight/2)+" highestYinC="+highestYinCache);
boolean threadFinished = y >= roi.y+roi.height;
if (numThreads>1 && (threadWaiting || threadFinished)) // 'if' is not synchronized to avoid overhead
synchronized(this) {
notifyAll(); // we may have blocked another thread
//IJ.log("thread "+threadNumber+" @y="+y+" notifying");
}
if (threadFinished) return; // all done, break the loop
if (threadNumber==0) { // main thread checks for abort and ProgressBar
long time = System.currentTimeMillis();
if (time-lastTime>100) {
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; i<cachePointers.length; i++) //shift kernel pointers to new line
cachePointers[i] = (cachePointers[i] + cacheWidth*(y-previousY))%cache.length;
previousY = y;
if (numThreads>1) { // thread synchronization
int slowestThreadY = arrayMinNonNegative(yForThread);
if (y - slowestThreadY + kHeight > cacheHeight) { // we would overwrite data needed by another thread
synchronized(this) {
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();
return;
}
slowestThreadY = arrayMinNonNegative(yForThread);
} while (y - slowestThreadY + kHeight > cacheHeight);
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, sign, 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<array.length; i++)
if (array[i] > max) max = array[i];
return max;
}
//only checks non-negative numbers in the array. Returns Integer.MAX_VALUE if no such values.
private int arrayMinNonNegative(int[] array) {
int min = Integer.MAX_VALUE;
for (int i=0; i<array.length; i++)
if (array[i]>=0 && array[i]<min) min = array[i];
return min;
}
private void filterLine(float[] values, int width, float[] cache, int[] cachePointers, int kNPoints, int cacheLineP, Rectangle roi, int y,
double[] sums, float[] medianBuf1, float[] medianBuf2, float sign, float maxValue, boolean isFloat, int filterType,
boolean smallKernel, boolean sumFilter, boolean minOrMax, boolean minOrMaxOrOutliers) {
int valuesP = isFloat ? roi.x+y*width : 0;
float max = 0f;
float median = Float.isNaN(cache[cacheLineP]) ? 0 : cache[cacheLineP]; // a first guess
boolean fullCalculation = true;
for (int x=0; x<roi.width; x++, valuesP++) { // x is with respect to roi.x
if (fullCalculation) {
fullCalculation = smallKernel; //for small kernel, always use the full area, not incremental algorithm
if (minOrMaxOrOutliers)
max = getAreaMax(cache, x, cachePointers, 0, -Float.MAX_VALUE, sign);
if (minOrMax) {
values[valuesP] = max*sign;
continue;
}
else if (sumFilter)
getAreaSums(cache, x, cachePointers, sums);
} else {
if (minOrMaxOrOutliers) {
float newPointsMax = getSideMax(cache, x, cachePointers, true, sign);
if (newPointsMax >= max) { //compare with previous maximum 'max'
max = newPointsMax;
} else {
float removedPointsMax = getSideMax(cache, x, cachePointers, false, sign);
if (removedPointsMax >= max)
max = getAreaMax(cache, x, cachePointers, 1, newPointsMax, sign);
}
if (minOrMax) {
values[valuesP] = max*sign;
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*sign+threshold < max) { //for low outliers: median can't be higher than max (sign is +1)
median = getMedian(cache, x, cachePointers, medianBuf1, medianBuf2, kNPoints, median);
if (v*sign+threshold < median*sign)
v = median; //beyond threshold (below if sign=+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<pixelLineP+xminInside+widthInside; pp++,cp++)
cache[cp] = bPixels[pp]&0xff;
} else if (pixels instanceof short[]){
short[] sPixels = (short[])pixels;
for (int pp=pixelLineP+xminInside, cp=cacheLineP+padLeft; pp<pixelLineP+xminInside+widthInside; pp++,cp++)
cache[cp] = sPixels[pp]&0xffff;
} else if (pixels instanceof float[]) {
System.arraycopy(pixels, pixelLineP+xminInside, cache, cacheLineP+padLeft, widthInside);
} else { //RGB
int[] cPixels = (int[])pixels;
int shift = 16 - 8*colorChannel;
int byteMask = 255<<shift;
for (int pp=pixelLineP+xminInside, cp=cacheLineP+padLeft; pp<pixelLineP+xminInside+widthInside; pp++,cp++)
cache[cp] = (cPixels[pp]&byteMask)>>shift;
}
for (int cp=cacheLineP; cp<cacheLineP+padLeft; cp++)
cache[cp] = cache[cacheLineP+padLeft];
for (int cp=cacheLineP+padLeft+widthInside; cp<cacheLineP+padLeft+widthInside+padRight; cp++)
cache[cp] = cache[cacheLineP+padLeft+widthInside-1];
}
/** Write a line to pixels arrax, converting from float (not for float data!)
* No checking for overflow/underflow
*/
private static void writeLineToPixels(float[] values, Object pixels, int pixelP, int length, int colorChannel) {
if (pixels instanceof byte[]) {
byte[] bPixels = (byte[])pixels;
for (int i=0, p=pixelP; i<length; i++,p++)
bPixels[p] = (byte)(((int)(values[i] + 0.5f))&0xff);
} else if (pixels instanceof short[]) {
short[] sPixels = (short[])pixels;
for (int i=0, p=pixelP; i<length; i++,p++)
sPixels[p] = (short)(((int)(values[i] + 0.5f))&0xffff);
} else { //RGB
int[] cPixels = (int[])pixels;
int shift = 16 - 8*colorChannel;
int resetMask = 0xffffffff^(0xff<<shift);
for (int i=0, p=pixelP; i<length; i++,p++)
cPixels[p] = (cPixels[p]&resetMask) | (((int)(values[i] + 0.5f))<<shift);
}
}
/** Get max (or -min if sign=-1) within the kernel area.
* @param x between 0 and cacheWidth-1
* @param ignoreRight should be 0 for analyzing all data or 1 for leaving out the row at the right
* @param max should be -Float.MAX_VALUE or the smallest value the maximum can be */
private static float getAreaMax(float[] cache, int xCache0, int[] kernel, int ignoreRight, float max, float sign) {
for (int kk=0; kk<kernel.length; kk++) { // y within the cache stripe (we have 2 kernel pointers per cache line)
for (int p=kernel[kk++]+xCache0; p<=kernel[kk]+xCache0-ignoreRight; p++) {
float v = cache[p]*sign;
if (max < v) max = v;
}
}
return max;
}
/** Get max (or -min if sign=-1) at the right border inside or left border outside the kernel area.
* x between 0 and cacheWidth-1 */
private static float getSideMax(float[] cache, int xCache0, int[] kernel, boolean isRight, float sign) {
float max = -Float.MAX_VALUE;
if (!isRight) xCache0--;
for (int kk= isRight ? 1 : 0; kk<kernel.length; kk+=2) { // y within the cache stripe (we have 2 kernel pointers per cache line)
float v = cache[xCache0 + kernel[kk]]*sign;
if (max < v) max = v;
}
return max;
}
/** Get sum of values and values squared within the kernel area.
* x between 0 and cacheWidth-1
* Output is written to array sums[0] = sum; sums[1] = sum of squares */
private static void getAreaSums(float[] cache, int xCache0, int[] kernel, double[] sums) {
double sum=0, sum2=0;
for (int kk=0; kk<kernel.length; kk++) { // y within the cache stripe (we have 2 kernel pointers per cache line)
for (int p=kernel[kk++]+xCache0; p<=kernel[kk]+xCache0; p++) {
float v = cache[p];
sum += v;
sum2 += v*v;
}
}
sums[0] = sum;
sums[1] = sum2;
return;
}
/** Add all values and values squared at the right border inside minus at the left border outside the kernal area.
* Output is added or subtracted to/from array sums[0] += sum; sums[1] += sum of squares when at
* the right border, minus when at the left border */
private static void addSideSums(float[] cache, int xCache0, int[] kernel, double[] sums) {
double sum=0, sum2=0;
for (int kk=0; kk<kernel.length; /*k++;k++ below*/) {
float v = cache[kernel[kk++]+(xCache0-1)];
sum -= v;
sum2 -= v*v;
v = cache[kernel[kk++]+xCache0];
sum += v;
sum2 += v*v;
}
sums[0] += sum;
sums[1] += sum2;
return;
}
/** Get median of values within kernel-sized neighborhood. Kernel size kNPoints should be odd.
*/
private static float getMedian(float[] cache, int xCache0, int[] kernel,
float[] aboveBuf, float[]belowBuf, int kNPoints, float guess) {
int nAbove = 0, nBelow = 0;
for (int kk=0; kk<kernel.length; kk++) {
for (int p=kernel[kk++]+xCache0; p<=kernel[kk]+xCache0; p++) {
float v = cache[p];
if (v > 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<kernel.length; kk++) {
for (int p=kernel[kk++]+xCache0; p<=kernel[kk]+xCache0; p++) {
float v = cache[p];
if (Float.isNaN(v)) {
kNPoints--;
} else if (v > 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. <code>buf</code> will be modified.
* @param bufLength Number of values in <code>buf</code> 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<m) {
i=l ;
j=m ;
do {
while (buf[i]<med) i++ ;
while (med<buf[j]) j-- ;
dum=buf[j];
buf[j]=buf[i];
buf[i]=dum;
i++ ; j-- ;
} while ((j>=n) && (i<=n)) ;
if (j<n) l=i ;
if (n<i) m=j ;
med=buf[n] ;
}
return med ;
}
/** Reset region between inner rectangle 'roi' and outer rectangle 'roi1' to the snapshot */
private void resetRoiBoundary(ImageProcessor ip, Rectangle roi, Rectangle roi1) {
int width = ip.getWidth();
Object pixels = ip.getPixels();
Object snapshot = ip.getSnapshotPixels();
for (int y=roi1.y, p = roi1.x+y*width; y<roi.y; y++,p+=width)
System.arraycopy(snapshot, p, pixels, p, roi1.width);
int leftWidth = roi.x - roi1.x;
int rightWidth = roi1.x+roi1.width - (roi.x+roi.width);
for (int y=roi.y, pL=roi1.x+y*width, pR=roi.x+roi.width+y*width; y<roi.y+roi.height; y++,pL+=width,pR+=width) {
if (leftWidth > 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<roi1.y+roi1.height; y++,p+=width)
System.arraycopy(snapshot, p, pixels, p, roi1.width);
}
/** @deprecated
* Not needed any more, use the rank(ip, ...) method, which creates the kernel */
public void makeKernel(double radius) {
this.radius = radius;
}
/** Create a circular kernel (structuring element) of a given radius.
* @param radius:
* Radius = 0.5 includes the 4 neighbors of the pixel in the center,
* radius = 1 corresponds to a 3x3 kernel size.
* @param width: width of the roi (or image if roi width=image width) for filtering.
* @return:
* The output is an array that gives the length of each line of the structuring element
* (kernel) to the left (negative) and to the right (positive):
* [0] left in line 0, [1] right in line 0,
* [2] left in line 2, ...
* The maximum (absolute) value should be kernelRadius.
* Array elements at the end:
* length-2: nPoints, number of pixels in the kernel area
* length-1: kernelRadius in x direction (kernel width is 2*kernelRadius+1)
* Kernel height can be calculated as (array length - 1)/2 (odd number);
* Kernel radius in y direction is kernel height/2 (truncating integer division).
* Note that kernel width and height are the same for the circular kernels used here,
* but treated separately for the case of future extensions with non-circular kernels.
*/
protected int[] makeLineRadii(double radius) {
if (radius>=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<kHeight;i++)IJ.log(i+": "+kernel[2*i]+"-"+kernel[2*i+1]);
return kernel;
}
//kernel height
private int kHeight(int[] lineRadii) {
return (lineRadii.length-2)/2;
}
//kernel radius in x direction. width is 2+kRadius+1
private int kRadius(int[] lineRadii) {
return lineRadii[lineRadii.length-1];
}
//number of points in kernal area
private int kNPoints(int[] lineRadii) {
return lineRadii[lineRadii.length-2];
}
//cache pointers for a given kernel
private int[] makeCachePointers(int[] lineRadii, int cacheWidth) {
int kRadius = kRadius(lineRadii);
int kHeight = kHeight(lineRadii);
int[] cachePointers = new int[2*kHeight];
for (int i=0; i<kHeight; i++) {
cachePointers[2*i] = i*cacheWidth+kRadius + lineRadii[2*i];
cachePointers[2*i+1] = i*cacheWidth+kRadius + lineRadii[2*i+1];
}
return cachePointers;
}
void showMasks() {
int w=150, h=150;
ImageStack stack = new ImageStack(w, h);
//for (double r=0.1; r<3; r+=0.01) {
for (double r=0.5; r<50; r+=0.5) {
ImageProcessor ip = new FloatProcessor(w,h,new int[w*h]);
float[] pixels = (float[])ip.getPixels();
int[] lineRadii = makeLineRadii(r);
int kHeight = kHeight(lineRadii);
int kRadius = kRadius(lineRadii);
int y0 = h/2-kHeight/2;
for (int i = 0, y = y0; i<kHeight; i++, y++)
for (int x = w/2+lineRadii[2*i], p = x+y*w; x <= w/2+lineRadii[2*i+1]; x++, p++)
pixels[p] = 1f;
stack.addSlice("radius="+r+", size="+(2*kRadius+1), ip);
}
new ImagePlus("Masks", stack).show();
}
/** This method is called by ImageJ to set the number of calls to run(ip)
* corresponding to 100% of the progress bar */
public void setNPasses (int nPasses) {
this.nPasses = nPasses;
pass = 0;
}
private void showProgress(double percent, boolean rgb) {
int nPasses2 = rgb?nPasses*3:nPasses;
percent = (double)pass/nPasses2 + percent/nPasses2;
IJ.showProgress(percent);
}
}