| BackgroundSubtracter |

    package ij.plugin.filter;
    import ij.*;
    import ij.gui.*;
    import ij.process.*;
    import ij.measure.*;
    import java.awt.*;

    /** Implements ImageJ's Subtract Background command. Based on
        the NIH Image Pascal version by Michael Castle and Janice
        Keller of the University of Michigan Mental Health Research
        Institute. Rolling ball algorithm inspired by Stanley
        Sternberg's article, "Biomedical Image Processing",
        IEEE Computer, January 1983.
    */
    public class BackgroundSubtracter implements PlugInFilter {

        private static int radius = 50; // default rolling ball radius
        private static boolean lightBackground = Prefs.get("bs.background", true);
        private ImagePlus imp;
        private boolean canceled;
        private int slice;
        private boolean invert;

        public int setup(String arg, ImagePlus imp) {
            IJ.register(BackgroundSubtracter.class);
            this.imp = imp;
            if (imp!=null) {
                showDialog();
                if (canceled)
                    return DONE;
            }
            IJ.register(BackgroundSubtracter.class);
            return IJ.setupDialog(imp, DOES_8G+DOES_16+DOES_RGB);
        }

        public void run(ImageProcessor ip) {
            if (canceled)
                return;
            slice++;
            if (slice>1)
                IJ.showStatus("Subtract Background: "+slice+"/"+imp.getStackSize());
            if (ip instanceof ColorProcessor)
                subtractRGBBackround((ColorProcessor)ip, radius);
            else
                subtractBackround(ip, radius);
            if (slice==1 && ip instanceof ShortProcessor)
                imp.getProcessor().resetMinAndMax();
        }

        public void showDialog() {
            String options = Macro.getOptions();
            if  (options!=null)
                Macro.setOptions(options.replaceAll("white", "light"));
            GenericDialog gd = new GenericDialog("Subtract Background");
            gd.addNumericField("Rolling Ball Radius:", radius, 0);
            gd.addCheckbox("Light Background", lightBackground);
            gd.showDialog();
            if (gd.wasCanceled())
                canceled = true;
            else {
                radius = (int)gd.getNextNumber();
                lightBackground = gd.getNextBoolean();
                Prefs.set("bs.background", lightBackground);
            }
            boolean invertedLut = imp.isInvertedLut();
            invert = (invertedLut && !lightBackground) || (!invertedLut && lightBackground);
        }

        public void subtractRGBBackround(ColorProcessor ip, int ballRadius) {
            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 brightness = new ByteProcessor(width, height, B, null);
            subtractBackround(brightness, radius);
            ip.setHSB(H, S, (byte[])brightness.getPixels());
        }

        /** Implements a rolling-ball algorithm for the removal of smooth continuous background
            from a two-dimensional gel image.  It rolls the ball (actually a square patch on the
            top of a sphere) on a low-resolution (by a factor of 'shrinkfactor' times) copy of
            the original image in order to increase speed with little loss in accuracy.  It uses
            interpolation and extrapolation to blow the shrunk image to full size.
        */
        public void subtractBackround(ImageProcessor ip, int ballRadius) {
            if (imp!=null)
                imp.killRoi();
            else
                ip.resetRoi();
            ip.setProgressBar(null);
            IJ.showProgress(0.0);
            if (invert)
                ip.invert();
            RollingBall ball = new RollingBall(ballRadius);
            //new ImagePlus("ball", new ByteProcessor(ball.patchwidth+1, ball.patchwidth+1, ball.data, null)).show();
            //ImageProcessor smallImage = ip.resize(ip.getWidth()/ball.shrinkfactor, ip.getHeight()/ball.shrinkfactor);
            ImageProcessor smallImage = shrinkImage(ip, ball.shrinkfactor);
            //new ImagePlus("small image", smallImage).show();
            if (slice==1)
                IJ.showStatus("Rolling ball ("+ball.shrinkfactor+")...");
            ImageProcessor background;
            if (ip instanceof ShortProcessor) {
                background = rollBall16(ball, ip, smallImage);
                interpolateBackground16(background, ball);
                extrapolateBackground16(background, ball);
            } else {
                background = rollBall(ball, ip, smallImage);
                interpolateBackground(background, ball);
                extrapolateBackground(background, ball);
            }
            IJ.showProgress(0.9);
            if (IJ.altKeyDown())
                new ImagePlus("background", background).show();
            ip.copyBits(background, 0, 0, Blitter.SUBTRACT);
            if (invert)
                ip.invert();
            IJ.showProgress(1.0);
        }

        /** 'Rolls' a filtering object over a (shrunken) image in order to find the
            image's smooth continuous background.  For the purpose of explaining this
            algorithm, imagine that the 2D grayscale image has a third (height) dimension
            defined by the intensity value at every point in the image.  The center of
            the filtering object, a patch from the top of a sphere having radius BallRadius,
            is moved along each scan line of the image so that the patch is tangent to the
            image at one or more points with every other point on the patch below the
            corresponding (x,y) point of the image.  Any point either on or below the patch
            during this process is considered part of the background.  Shrinking the image
            before running this procedure is advised due to the fourth-degree complexity
            of the algorithm.
        */
        ImageProcessor rollBall(RollingBall ball, ImageProcessor image, ImageProcessor smallImage) {

            int halfpatchwidth;     //distance in x or y from patch center to any edge
            int ptsbelowlastpatch;  //number of points we may ignore because they were below last patch
            int xpt2, ypt2;         // current (x,y) point in the patch relative to upper left corner
            int xval, yval;         // location in ball in shrunken image coordinates
            int zdif;               // difference in z (height) between point on ball and point on image
            int zmin;               // smallest zdif for ball patch with center at current point
            int zctr;               // current height of the center of the sphere of which the patch is a part
            int zadd;               // height of a point on patch relative to the xy-plane of the shrunken image
            int ballpt;             // index to array storing the precomputed ball patch
            int imgpt;              // index to array storing the shrunken image
            int backgrpt;           // index to array storing the calculated background
            int ybackgrpt;          // displacement to current background scan line
            int p1, p2;             // temporary indexes to background, ball, or small image
            int ybackgrinc;             // distance in memory between two shrunken y-points in background
            int smallimagewidth;    // length of a scan line in shrunken image
            int left, right, top, bottom;
            byte[] pixels = (byte[])smallImage.getPixels();
            byte[] patch = ball.data;
            int width = image.getWidth();
            int height = image.getHeight();
            int swidth = smallImage.getWidth();
            int sheight = smallImage.getHeight();
            ImageProcessor background = new ByteProcessor(width, height);
            byte[] backgroundpixels = (byte[])background.getPixels();
            int shrinkfactor = ball.shrinkfactor;
            int leftroll = 0;
            int rightroll = width/shrinkfactor-1;
            int toproll = 0;
            int bottomroll = height/shrinkfactor-1;

            left = 1;
            right = rightroll - leftroll - 1;
            top = 1;
            bottom = bottomroll - toproll - 1;
            smallimagewidth = swidth;
            int patchwidth = ball.patchwidth;
            halfpatchwidth = patchwidth / 2;
            ybackgrinc = shrinkfactor*width; // real dist btwn 2 adjacent (dy=1) shrunk pts
            zctr = 0; // start z-center in the xy-plane
            for (int ypt=top; ypt<=(bottom+patchwidth); ypt++) {
                for (int xpt=left; xpt<=(right+patchwidth); xpt++) {// while patch is tangent to edges or within image...
                    // xpt is far right edge of ball patch
                    // do we have to move the patch up or down to make it tangent to but not above image?...
                    zmin = 255; // highest could ever be 255
                    ballpt = 0;
                    ypt2 = ypt - patchwidth; // ypt2 is top edge of ball patch
                    imgpt = ypt2*smallimagewidth + xpt - patchwidth;
                    while (ypt2<=ypt) {
                        xpt2 = xpt-patchwidth; // xpt2 is far left edge of ball patch
                        while (xpt2<=xpt) { // check every point on ball patch
                            // only examine points on
                            if ((xpt2>=left) && (xpt2<=right) && (ypt2>=top) && (ypt2<=bottom)) {
                                p1 = ballpt;
                                p2 = imgpt;
                                zdif = (pixels[p2]&255) - (zctr + (patch[p1]&255));  //curve - circle points
                                //if (xpt==50 && ypt==50) IJ.write(zdif
                                //+" "+(pixels[p2]&255)
                                //+" "+zctr
                                //+" "+(patch[p1]&255)
                                //+" "+p2
                                //+" "+p1
                                //);
                                if (zdif<zmin) // keep most negative, since ball should always be below curve
                                    zmin = zdif;
                            } // if xpt2,ypt2
                            ballpt++;
                            xpt2++;
                            imgpt++;
                        } // while xpt2
                        ypt2++;
                        imgpt = imgpt - patchwidth - 1 + smallimagewidth;
                    }  // while ypt2
                    if (zmin!=0)
                        zctr += zmin; // move ball up or down if we find a new minimum
                    if (zmin<0)
                        ptsbelowlastpatch = halfpatchwidth; // ignore left half of ball patch when dz < 0
                    else
                        ptsbelowlastpatch = 0;
                    // now compare every point on ball with background,  and keep highest number
                    yval = ypt - patchwidth;
                    ypt2 = 0;
                    ballpt = 0;
                    ybackgrpt = (yval - top + 1) * ybackgrinc;
                    while (ypt2<=patchwidth) {
                        xval = xpt - patchwidth + ptsbelowlastpatch;
                        xpt2 = ptsbelowlastpatch;
                        ballpt += ptsbelowlastpatch;
                        backgrpt = ybackgrpt + (xval - left + 1) * shrinkfactor;
                        while (xpt2<=patchwidth) { // for all the points in the ball patch
                            if ((xval >= left) && (xval <= right) && (yval >= top) && (yval <= bottom)) {
                                p1 = ballpt;
                                zadd = zctr + (patch[p1]&255);
                                p1 = backgrpt;
                                //if (backgrpt>=backgroundpixels.length) backgrpt = 0; //(debug)
                                if (zadd>(backgroundpixels[p1]&255)) //keep largest adjustment}
                                    backgroundpixels[p1] = (byte)zadd;
                            }
                            ballpt++;
                            xval++;
                            xpt2++;
                            backgrpt += shrinkfactor; // move to next point in x
                        } // while xpt2
                        yval++;
                        ypt2++;
                        ybackgrpt += ybackgrinc; // move to next point in y
                    } // while ypt2
                } // for xpt
                if (ypt%20==0)
                    IJ.showProgress(0.2+0.6*ypt/(bottom+patchwidth));
            } // for ypt
            return background;
        }

        /** Creates a lower resolution image for ball-rolling. */
        ImageProcessor shrinkImage(ImageProcessor ip, int shrinkfactor) {
            int width = ip.getWidth();
            int height = ip.getHeight();
            int swidth = width/shrinkfactor;
            int sheight = height/shrinkfactor;
            ImageProcessor ip2 = ip.duplicate();
            ip2.smooth();
            IJ.showProgress(0.1);
            ImageProcessor smallImage = ip.createProcessor(swidth, sheight);
            int xmaskmin, ymaskmin, min, thispixel;
            for (int y=0; y<sheight; y++) {
                for (int x=0; x<swidth; x++) {
                    xmaskmin = shrinkfactor*x;
                    ymaskmin = shrinkfactor*y;
                    min = 65535;
                    for (int j=0; j<shrinkfactor; j++) {
                        for (int k=0; k<shrinkfactor; k++) {
                            thispixel = ip2.getPixel(xmaskmin+j, ymaskmin+k);
                            if (thispixel<min)
                                min = thispixel;
                        }
                    }
                    smallImage.putPixel(x,y,min); // each point in small image is minimum of its neighborhood
                }
            }
            //new ImagePlus("smallImage", smallImage).show();
            return smallImage;
        }

    /** Uses bilinear interpolation to find the points in the full-scale background
            given the points from the shrunken image background.  Since the shrunken background
            is found from an image composed of minima (over a sufficiently large mask), it
            is certain that no point in the full-scale interpolated background has a higher
            pixel value than the corresponding point in the original image
        */
        void interpolateBackground(ImageProcessor background, RollingBall ball) {
            int hloc, vloc; // position of current pixel in calculated background
            int vinc;       // memory offset from current calculated pos to current interpolated pos
            int lastvalue, nextvalue;   // calculated pixel values between which we are interpolating
            int p;          // pointer to current interpolated pixel value
            int bglastptr, bgnextptr;   // pointers to calculated pixel values between which we are interpolating

            int width = background.getWidth();
            int height = background.getHeight();
            int shrinkfactor = ball.shrinkfactor;
            int leftroll = 0;
            int rightroll = width/shrinkfactor-1;
            int toproll = 0;
            int bottomroll = height/shrinkfactor-1;
            byte[] pixels = (byte[])background.getPixels();

            vloc = 0;
            for (int j=1; j<=(bottomroll-toproll-1); j++) { //interpolate to find background interior
                hloc = 0;
                vloc += shrinkfactor;
                for (int i=1; i<=(rightroll-leftroll); i++) {
                    hloc += shrinkfactor;
                    bgnextptr = vloc*width + hloc;
                    bglastptr = bgnextptr - shrinkfactor;
                    nextvalue = pixels[bgnextptr]&255;
                    lastvalue = pixels[bglastptr]&255;
                    for (int ii=1; ii<=(shrinkfactor-1); ii++) { //interpolate horizontally
                        p = bgnextptr - ii;
                        pixels[p] = (byte)(lastvalue+(shrinkfactor-ii)*(nextvalue-lastvalue)/shrinkfactor);
                    }
                    for (int ii=0; ii<=(shrinkfactor-1); ii++) { //interpolate vertically
                        bglastptr = (vloc-shrinkfactor)*width+hloc-ii;
                        bgnextptr = vloc*width+hloc-ii;
                        lastvalue = pixels[bglastptr]&255;
                        nextvalue = pixels[bgnextptr]&255;
                        vinc = 0;
                        for (int jj=1; jj<=(shrinkfactor-1); jj++) {
                            vinc = vinc-width;
                            p = bgnextptr+vinc;
                            pixels[p] = (byte)(lastvalue+(shrinkfactor-jj)*(nextvalue-lastvalue)/shrinkfactor);
                        } // for jj
                    } // for ii
                } // for i
            } // for j
        }

        /** Uses linear extrapolation to find pixel values on the top, left, right,
            and bottom edges of the background.  First it finds the top and bottom
            edge points by extrapolating from the edges of the calculated and
            interpolated background interior.  Then it uses the edge points on the
            new calculated, interpolated, and extrapolated background to find all
            of the left and right edge points.  If extrapolation yields values
            below zero or above 255, then they are set to zero and 255 respectively.
        */
        void extrapolateBackground(ImageProcessor background, RollingBall ball) {

            int edgeslope;          // difference of last two consecutive pixel values on an edge
            int pvalue;             // current extrapolated pixel value
            int lastvalue, nextvalue;   //calculated pixel values from which we are extrapolating
            int p;                  // pointer to current extrapolated pixel value
            int bglastptr, bgnextptr;   // pointers to calculated pixel values from which we are extrapolating

            int width = background.getWidth();
            int height = background.getHeight();
            int shrinkfactor = ball.shrinkfactor;
            int leftroll = 0;
            int rightroll = width/shrinkfactor-1;
            int toproll = 0;
            int bottomroll = height/shrinkfactor-1;
            byte[] pixels = (byte[])background.getPixels();

            for (int hloc=shrinkfactor; hloc<=(shrinkfactor*(rightroll-leftroll)-1); hloc++) {
                // extrapolate on top and bottom
                bglastptr = shrinkfactor*width+hloc;
                bgnextptr = (shrinkfactor+1)*width+hloc;
                lastvalue = pixels[bglastptr]&255;
                nextvalue = pixels[bgnextptr]&255;
                edgeslope = nextvalue-lastvalue;
                p = bglastptr;
                pvalue = lastvalue;
                for (int jj=1; jj<=shrinkfactor; jj++) {
                    p = p-width;
                    pvalue = pvalue-edgeslope;
                    if (pvalue<0)
                        pixels[p] = 0;
                    else if (pvalue>255)
                        pixels[p] = (byte)255;
                    else
                        pixels[p] = (byte)pvalue;
                } // for jj
                bglastptr = (shrinkfactor*(bottomroll-toproll-1)-1)*width+hloc;
                bgnextptr = shrinkfactor*(bottomroll-toproll-1)*width+hloc;
                lastvalue = pixels[bglastptr]&255;
                nextvalue = pixels[bgnextptr]&255;
                edgeslope = nextvalue-lastvalue;
                p = bgnextptr;
                pvalue = nextvalue;
                for (int jj=1; jj<=((height-1)-shrinkfactor*(bottomroll-toproll-1)); jj++) {
                    p += width;
                    pvalue += edgeslope;
                    if (pvalue<0)
                        pixels[p] = 0;
                    else if (pvalue>255)
                        pixels[p] = (byte)255;
                    else
                        pixels[p] = (byte)pvalue;
                } // for jj
            } // for hloc
            for (int vloc=0; vloc<height; vloc++) {
                // extrapolate on left and right
                bglastptr = vloc*width+shrinkfactor;
                bgnextptr = bglastptr+1;
                lastvalue = pixels[bglastptr]&255;
                nextvalue = pixels[bgnextptr]&255;
                edgeslope = nextvalue-lastvalue;
                p = bglastptr;
                pvalue = lastvalue;
                for (int ii=1; ii<=shrinkfactor; ii++) {
                    p--;
                    pvalue = pvalue - edgeslope;
                    if (pvalue<0)
                        pixels[p] = 0;
                    else if (pvalue>255)
                        pixels[p] = (byte)255;
                    else
                        pixels[p] = (byte)pvalue;
                } // for ii
                bgnextptr = vloc*width+shrinkfactor*(rightroll-leftroll-1)-1;
                bglastptr = bgnextptr-1;
                lastvalue = pixels[bglastptr]&255;
                nextvalue = pixels[bgnextptr]&255;
                edgeslope = nextvalue-lastvalue;
                p = bgnextptr;
                pvalue = nextvalue;
                for (int ii=1; ii<=((width-1)-shrinkfactor*(rightroll-leftroll-1)+1); ii++) {
                    p++;
                    pvalue = pvalue+edgeslope;
                    if (pvalue<0)
                        pixels[p] = 0;
                    else if (pvalue>255)
                        pixels[p] = (byte)255;
                    else
                        pixels[p] = (byte)pvalue;
                } // for ii
            } // for vloc
        }

        /** This is a 16-bit version of the rollBall() method. */
        ImageProcessor rollBall16(RollingBall ball, ImageProcessor image, ImageProcessor smallImage) {

            int halfpatchwidth;     //distance in x or y from patch center to any edge
            int ptsbelowlastpatch;  //number of points we may ignore because they were below last patch
            int xpt2, ypt2;         // current (x,y) point in the patch relative to upper left corner
            int xval, yval;         // location in ball in shrunken image coordinates
            int zdif;               // difference in z (height) between point on ball and point on image
            int zmin;               // smallest zdif for ball patch with center at current point
            int zctr;               // current height of the center of the sphere of which the patch is a part
            int zadd;               // height of a point on patch relative to the xy-plane of the shrunken image
            int ballpt;             // index to array storing the precomputed ball patch
            int imgpt;              // index to array storing the shrunken image
            int backgrpt;           // index to array storing the calculated background
            int ybackgrpt;          // displacement to current background scan line
            int p1, p2;             // temporary indexes to background, ball, or small image
            int ybackgrinc;             // distance in memory between two shrunken y-points in background
            int smallimagewidth;    // length of a scan line in shrunken image
            int left, right, top, bottom;
            short[] pixels = (short[])smallImage.getPixels();
            byte[] patch = ball.data;
            int width = image.getWidth();
            int height = image.getHeight();
            int swidth = smallImage.getWidth();
            int sheight = smallImage.getHeight();
            ImageProcessor background = new ShortProcessor(width, height);
            short[] backgroundpixels = (short[])background.getPixels();
            int shrinkfactor = ball.shrinkfactor;
            int leftroll = 0;
            int rightroll = width/shrinkfactor-1;
            int toproll = 0;
            int bottomroll = height/shrinkfactor-1;

            left = 1;
            right = rightroll - leftroll - 1;
            top = 1;
            bottom = bottomroll - toproll - 1;
            smallimagewidth = swidth;
            int patchwidth = ball.patchwidth;
            halfpatchwidth = patchwidth / 2;
            ybackgrinc = shrinkfactor*width; // real dist btwn 2 adjacent (dy=1) shrunk pts
            zctr = 0; // start z-center in the xy-plane
            for (int ypt=top; ypt<=(bottom+patchwidth); ypt++) {
                for (int xpt=left; xpt<=(right+patchwidth); xpt++) {// while patch is tangent to edges or within image...
                    // xpt is far right edge of ball patch
                    // do we have to move the patch up or down to make it tangent to but not above image?...
                    zmin = 65535; // highest possible value
                    ballpt = 0;
                    ypt2 = ypt - patchwidth; // ypt2 is top edge of ball patch
                    imgpt = ypt2*smallimagewidth + xpt - patchwidth;
                    while (ypt2<=ypt) {
                        xpt2 = xpt-patchwidth; // xpt2 is far left edge of ball patch
                        while (xpt2<=xpt) { // check every point on ball patch
                            // only examine points on
                            if ((xpt2>=left) && (xpt2<=right) && (ypt2>=top) && (ypt2<=bottom)) {
                                p1 = ballpt;
                                p2 = imgpt;
                                zdif = (pixels[p2]&0xffff) - (zctr + (patch[p1]&255));  //curve - circle points
                                if (zdif<zmin) // keep most negative, since ball should always be below curve
                                    zmin = zdif;
                            } // if xpt2,ypt2
                            ballpt++;
                            xpt2++;
                            imgpt++;
                        } // while xpt2
                        ypt2++;
                        imgpt = imgpt - patchwidth - 1 + smallimagewidth;
                    }  // while ypt2
                    if (zmin!=0)
                        zctr += zmin; // move ball up or down if we find a new minimum
                    if (zmin<0)
                        ptsbelowlastpatch = halfpatchwidth; // ignore left half of ball patch when dz < 0
                    else
                        ptsbelowlastpatch = 0;
                    // now compare every point on ball with background,  and keep highest number
                    yval = ypt - patchwidth;
                    ypt2 = 0;
                    ballpt = 0;
                    ybackgrpt = (yval - top + 1) * ybackgrinc;
                    while (ypt2<=patchwidth) {
                        xval = xpt - patchwidth + ptsbelowlastpatch;
                        xpt2 = ptsbelowlastpatch;
                        ballpt += ptsbelowlastpatch;
                        backgrpt = ybackgrpt + (xval - left + 1) * shrinkfactor;
                        while (xpt2<=patchwidth) { // for all the points in the ball patch
                            if ((xval >= left) && (xval <= right) && (yval >= top) && (yval <= bottom)) {
                                p1 = ballpt;
                                zadd = zctr + (patch[p1]&255);
                                p1 = backgrpt;
                                //if (backgrpt>=backgroundpixels.length) backgrpt = 0; //(debug)
                                if (zadd>(backgroundpixels[p1]&0xffff)) //keep largest adjustment}
                                    backgroundpixels[p1] = (short)zadd;
                            }
                            ballpt++;
                            xval++;
                            xpt2++;
                            backgrpt += shrinkfactor; // move to next point in x
                        } // while xpt2
                        yval++;
                        ypt2++;
                        ybackgrpt += ybackgrinc; // move to next point in y
                    } // while ypt2
                } // for xpt
                if (ypt%20==0)
                    IJ.showProgress(0.2+0.6*ypt/(bottom+patchwidth));
            } // for ypt
            return background;
        }

        /** This is a 16-bit version of the interpolateBackground(0 method. */
        void interpolateBackground16(ImageProcessor background, RollingBall ball) {
            int hloc, vloc; // position of current pixel in calculated background
            int vinc;       // memory offset from current calculated pos to current interpolated pos
            int lastvalue, nextvalue;   // calculated pixel values between which we are interpolating
            int p;          // pointer to current interpolated pixel value
            int bglastptr, bgnextptr;   // pointers to calculated pixel values between which we are interpolating

            int width = background.getWidth();
            int height = background.getHeight();
            int shrinkfactor = ball.shrinkfactor;
            int leftroll = 0;
            int rightroll = width/shrinkfactor-1;
            int toproll = 0;
            int bottomroll = height/shrinkfactor-1;
            short[] pixels = (short[])background.getPixels();

            vloc = 0;
            for (int j=1; j<=(bottomroll-toproll-1); j++) { //interpolate to find background interior
                hloc = 0;
                vloc += shrinkfactor;
                for (int i=1; i<=(rightroll-leftroll); i++) {
                    hloc += shrinkfactor;
                    bgnextptr = vloc*width + hloc;
                    bglastptr = bgnextptr - shrinkfactor;
                    nextvalue = pixels[bgnextptr]&0xffff;
                    lastvalue = pixels[bglastptr]&0xffff;
                    for (int ii=1; ii<=(shrinkfactor-1); ii++) { //interpolate horizontally
                        p = bgnextptr - ii;
                        pixels[p] = (short)(lastvalue+(shrinkfactor-ii)*(nextvalue-lastvalue)/shrinkfactor);
                    }
                    for (int ii=0; ii<=(shrinkfactor-1); ii++) { //interpolate vertically
                        bglastptr = (vloc-shrinkfactor)*width+hloc-ii;
                        bgnextptr = vloc*width+hloc-ii;
                        lastvalue = pixels[bglastptr]&0xffff;
                        nextvalue = pixels[bgnextptr]&0xffff;
                        vinc = 0;
                        for (int jj=1; jj<=(shrinkfactor-1); jj++) {
                            vinc = vinc-width;
                            p = bgnextptr+vinc;
                            pixels[p] = (short)(lastvalue+(shrinkfactor-jj)*(nextvalue-lastvalue)/shrinkfactor);
                        } // for jj
                    } // for ii
                } // for i
            } // for j
        }

        /** This is a 16-bit version of the extrapolateBackground() method. */
        void extrapolateBackground16(ImageProcessor background, RollingBall ball) {

            int edgeslope;          // difference of last two consecutive pixel values on an edge
            int pvalue;             // current extrapolated pixel value
            int lastvalue, nextvalue;   //calculated pixel values from which we are extrapolating
            int p;                  // pointer to current extrapolated pixel value
            int bglastptr, bgnextptr;   // pointers to calculated pixel values from which we are extrapolating

            int width = background.getWidth();
            int height = background.getHeight();
            int shrinkfactor = ball.shrinkfactor;
            int leftroll = 0;
            int rightroll = width/shrinkfactor-1;
            int toproll = 0;
            int bottomroll = height/shrinkfactor-1;
            short[] pixels = (short[])background.getPixels();

            for (int hloc=shrinkfactor; hloc<=(shrinkfactor*(rightroll-leftroll)-1); hloc++) {
                // extrapolate on top and bottom
                bglastptr = shrinkfactor*width+hloc;
                bgnextptr = (shrinkfactor+1)*width+hloc;
                lastvalue = pixels[bglastptr]&0xffff;
                nextvalue = pixels[bgnextptr]&0xffff;
                edgeslope = nextvalue-lastvalue;
                p = bglastptr;
                pvalue = lastvalue;
                for (int jj=1; jj<=shrinkfactor; jj++) {
                    p = p-width;
                    pvalue = pvalue-edgeslope;
                    if (pvalue<0)
                        pixels[p] = 0;
                    else if (pvalue>65535)
                        pixels[p] = (short)65535;
                    else
                        pixels[p] = (short)pvalue;
                } // for jj
                bglastptr = (shrinkfactor*(bottomroll-toproll-1)-1)*width+hloc;
                bgnextptr = shrinkfactor*(bottomroll-toproll-1)*width+hloc;
                lastvalue = pixels[bglastptr]&0xffff;
                nextvalue = pixels[bgnextptr]&0xffff;
                edgeslope = nextvalue-lastvalue;
                p = bgnextptr;
                pvalue = nextvalue;
                for (int jj=1; jj<=((height-1)-shrinkfactor*(bottomroll-toproll-1)); jj++) {
                    p += width;
                    pvalue += edgeslope;
                    if (pvalue<0)
                        pixels[p] = 0;
                    else if (pvalue>65535)
                        pixels[p] = (short)65535;
                    else
                        pixels[p] = (short)pvalue;
                } // for jj
            } // for hloc
            for (int vloc=0; vloc<height; vloc++) {
                // extrapolate on left and right
                bglastptr = vloc*width+shrinkfactor;
                bgnextptr = bglastptr+1;
                lastvalue = pixels[bglastptr]&0xffff;
                nextvalue = pixels[bgnextptr]&0xffff;
                edgeslope = nextvalue-lastvalue;
                p = bglastptr;
                pvalue = lastvalue;
                for (int ii=1; ii<=shrinkfactor; ii++) {
                    p--;
                    pvalue = pvalue - edgeslope;
                    if (pvalue<0)
                        pixels[p] = 0;
                    else if (pvalue>65535)
                        pixels[p] = (short)65535;
                    else
                        pixels[p] = (short)pvalue;
                } // for ii
                bgnextptr = vloc*width+shrinkfactor*(rightroll-leftroll-1)-1;
                bglastptr = bgnextptr-1;
                lastvalue = pixels[bglastptr]&0xffff;
                nextvalue = pixels[bgnextptr]&0xffff;
                edgeslope = nextvalue-lastvalue;
                p = bgnextptr;
                pvalue = nextvalue;
                for (int ii=1; ii<=((width-1)-shrinkfactor*(rightroll-leftroll-1)+1); ii++) {
                    p++;
                    pvalue = pvalue+edgeslope;
                    if (pvalue<0)
                        pixels[p] = 0;
                    else if (pvalue>65535)
                        pixels[p] = (short)65535;
                    else
                        pixels[p] = (short)pvalue;
                } // for ii
            } // for vloc
        }

    }

    class RollingBall {

        byte[] data;
        int patchwidth;
        int shrinkfactor;

        RollingBall(int radius) {
            int arcTrimPer;
            if (radius<=10) {
                shrinkfactor = 1;
                arcTrimPer = 12; // trim 24% in x and y
            } else if (radius<=30) {
                shrinkfactor = 2;
                arcTrimPer = 12; // trim 24% in x and y
            } else if (radius<=100) {
                shrinkfactor = 4;
                arcTrimPer = 16; // trim 32% in x and y
            } else {
                shrinkfactor = 8;
                arcTrimPer = 20; // 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(int ballradius, int arcTrimPer) {
            int 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
            int smallballradius, diam; // radius and diameter of rolling ball
            int temp;           // value must be >=0 to take square root
            int halfpatchwidth; // distance in x or y from center of patch to any edge
            int ballsize;       // size of rolling ball array

            this.shrinkfactor = shrinkfactor;
            smallballradius = ballradius/shrinkfactor;
            if (smallballradius<1)
                smallballradius = 1;
            rsquare = smallballradius*smallballradius;
            diam = smallballradius*2;
            xtrim = (arcTrimPer*diam)/100; // only use a patch of the rolling ball
            patchwidth = diam - xtrim - xtrim;
            halfpatchwidth = smallballradius - xtrim;
            ballsize = (patchwidth+1)*(patchwidth+1);
            data = new byte[ballsize];

            for (int i=0; i<ballsize; i++) {
                xval = i % (patchwidth+1) - halfpatchwidth;
                yval = i / (patchwidth+1) - halfpatchwidth;
                temp = rsquare - (xval*xval) - (yval*yval);
                if (temp >= 0)
                    data[i] = (byte)Math.round(Math.sqrt(temp));
                else
                    data[i] = 0;
            }
        }

    }

| BackgroundSubtracter |