package ini.trakem2.imaging;
import ij.IJ;
import ij.ImageJ;
import ij.ImagePlus;
import ij.process.FloatProcessor;
import ij.process.ShortProcessor;
/** An integral histogram for 2d images.
*
* @author Albert Cardona
*
*/
public class IntegralHistogram2d {
/** Returns a byte[] of dimensions (width + 1) * (height + 1) * nBins,
* representing the integral histogram of the 2d image.
*
* @param pixels The image.
* @param width The width of the image.
* @param height The height of the image.
* @param nBins The desired number of histogram bins.
* @param min The lowest value, under which all values are stashed into the first bin.
* @param max The highest value, above which all values are stashed into the last bin.
* @return The integral histogram, with the histogram as the fastest coordinate.
*/
static public final long[] integralHistogram2d(
final short[] pixels,
final int width,
final int height,
final int nBins,
final double min,
final double max)
{
final long[] hist = new long[(width + 1) * (height + 1) * nBins];
final double K = nBins / (max - min + 1.0d);
final int histWidth = (width + 1) * nBins;
for (int y=0, i=0; y<height; ++y) {
final int offset = histWidth * (y + 1);
// Index of first bin in hist for histogram on pixel i
int histIndex = offset + nBins;
for (int x=0; x<width; ++x, ++i, histIndex += nBins) {
// Compute bin index for pixels[i]
int binIndex = (int)(((pixels[i]&0xffff) - min) * K);
// Check for values over the max
if (binIndex >= nBins) binIndex = nBins -1;
// Copy histogram on pixel i-1
for (int b=0; b<nBins; b++) {
hist[histIndex + b] = hist[histIndex + b - nBins];
}
// Increase corresponding bin
hist[histIndex + binIndex] += 1;
}
}
// Sum columns
final int w = (width + 1) * nBins;
final int h = (height + 1);
// Add the previous row to the current row
// (Skip the first row, which is all zeros,
// and skip the first nBins columns, which are all zeros)
for (int y=1; y<h; ++y) {
for (int x=nBins; x<w; ++x) {
hist[y * w + x] += hist[(y - 1) * w + x];
}
}
/*
// print last 16 values
for (int i=hist.length-16; i < hist.length; ++i) {
System.out.print(hist[i] + ", ");
}
System.out.print('\n');
*/
return hist;
}
/**
*
* @param width The width of the image.
* @param height The height of the image.
* @param hist The array containing the integral histogram as generated by {@link IntegralHistogram2d#integralHistogram2d(short[], int, int, int, double, double)}.
* @param nBins The number of bins of the integral histogram.
* @param min The lowest value, under which all values are stashed into the first bin.
* @param max The highest value, above which all values are stashed into the last bin.
* @param radius The radius to approximate the median.
* @return The pixels of the image in 16-bit.
*/
static public final short[] median(
final int width,
final int height,
final long[] hist,
final int nBins,
final double min,
final double max,
final int radius)
{
final short[] median = new short[width * height];
final int w = (width + 1) * nBins;
final double binWorth = (max - min + 1.0) / nBins;
for (int y = 0; y < height; ++y) {
final int offset = y * width;
for (int x = 0; x < width; ++x) {
// Naive
// To convert to hist array coordinates:
// x,y must be shifted by 1 to account for the leading empty row and column
// and x must be multiplied by nBins
int x0 = Math.max(0, x + 1 - radius) * nBins;
int y0 = Math.max(0, y + 1 - radius);
int x1 = Math.min(width, x + 1 + radius) * nBins;
int y1 = Math.max(0, y + 1 - radius);
int x2 = Math.min(width , x + 1 + radius) * nBins;
int y2 = Math.min(height, y + 1 + radius);
int x3 = Math.max(0, x + 1 - radius) * nBins;
int y3 = Math.min(height, y + 1 + radius);
// Half the area
final int halfHistCount = (int)((x2 - x0 + 1) * (y2 - y0 + 1) / nBins / 2.0);
long sum = 0;
float val = 0;
for (int b = 0; b < nBins; ++b, val += binWorth) {
long binCount = hist[y0 * w + x0 + b]
- hist[y1 * w + x1 + b]
+ hist[y2 * w + x2 + b]
- hist[y3 * w + x3 + b];
sum += binCount;
if (sum > halfHistCount) {
// median value has been reached
// Could set the median to val + binWorth, but why not attempt an approximation,
// which requires a fraction of binWorth according to the exact position of the half count:
median[offset + x] = (short)(val + binWorth * ((halfHistCount - (sum - binCount)) / (float)binCount));
break;
}
}
}
}
return median;
}
@SuppressWarnings("unused")
static private final void view(
final long[] hist,
final int w1,
final int h1,
final int nBins) {
final float[] pixels = new float[w1 * h1 * nBins];
for (int i=0; i<hist.length; ++i) {
pixels[i] = hist[i];
}
new ImagePlus("Integral Histogram", new FloatProcessor(w1 * nBins, h1, pixels)).show();
}
static public final ShortProcessor median(
final ShortProcessor sp,
final int nBins,
final double min,
final double max,
final int radius) {
final long[] hist = integralHistogram2d((short[])sp.getPixels(), sp.getWidth(), sp.getHeight(), nBins, min, max);
final short[] median = median(sp.getWidth(), sp.getHeight(), hist, nBins, min, max, radius);
return new ShortProcessor(sp.getWidth(), sp.getHeight(), median, null);
}
@SuppressWarnings("unused")
static public final void main(String[] args) {
final ShortProcessor sp = (ShortProcessor) IJ.openImage("/home/albert/Desktop/t2/bridge-16bit.tif").getProcessor();
final ShortProcessor filtered = median(sp, 64, 0, 65535, 100);
new ImageJ();
new ImagePlus("median", filtered).show();
}
}