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 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