Squiz Matrix  4.12.2
 All Data Structures Namespaces Functions Variables Pages
BackgroundSubtracter.java
1 package ij.plugin.filter;
2 import ij.*;
3 import ij.gui.*;
4 import ij.process.*;
5 import ij.measure.*;
6 import java.awt.*;
7 
15 public class BackgroundSubtracter implements PlugInFilter {
16 
17  private static int radius = 50; // default rolling ball radius
18  private static boolean whiteBackground = true;
19  private ImagePlus imp;
20  private boolean canceled;
21  private int slice;
22  private boolean invert;
23 
24  public int setup(String arg, ImagePlus imp) {
25  IJ.register(BackgroundSubtracter.class);
26  this.imp = imp;
27  if (imp!=null) {
28  showDialog();
29  if (canceled)
30  return DONE;
31  }
32  IJ.register(BackgroundSubtracter.class);
33  return IJ.setupDialog(imp, DOES_8G+DOES_RGB);
34  }
35 
36  public void run(ImageProcessor ip) {
37  if (canceled)
38  return;
39  slice++;
40  if (slice>1)
41  IJ.showStatus("Subtract Background: "+slice+"/"+imp.getStackSize());
42  if (ip instanceof ColorProcessor)
43  subtractRGBBackround((ColorProcessor)ip, radius);
44  else
45  subtractBackround(ip, radius);
46  }
47 
48  public void showDialog() {
49  GenericDialog gd = new GenericDialog("Subtract Background");
50  gd.addNumericField("Rolling Ball Radius:", radius, 0);
51  gd.addCheckbox("White Background", whiteBackground);
52  gd.showDialog();
53  if (gd.wasCanceled())
54  canceled = true;
55  else {
56  radius = (int)gd.getNextNumber();
57  whiteBackground = gd.getNextBoolean();
58  }
59  boolean invertedLut = imp.isInvertedLut();
60  invert = (invertedLut && !whiteBackground) || (!invertedLut && whiteBackground);
61  }
62 
63  public void subtractRGBBackround(ColorProcessor ip, int ballRadius) {
64  int width = ip.getWidth();
65  int height = ip.getHeight();
66  byte[] H = new byte[width*height];
67  byte[] S = new byte[width*height];
68  byte[] B = new byte[width*height];
69  ip.getHSB(H, S, B);
70  ByteProcessor brightness = new ByteProcessor(width, height, B, null);
71  subtractBackround(brightness, radius);
72  ip.setHSB(H, S, (byte[])brightness.getPixels());
73  }
74 
81  public void subtractBackround(ImageProcessor ip, int ballRadius) {
82  if (imp!=null)
83  imp.killRoi();
84  else
85  ip.resetRoi();
86  ip.setProgressBar(null);
87  IJ.showProgress(0.0);
88  if (invert)
89  ip.invert();
90  RollingBall ball = new RollingBall(ballRadius);
91  //new ImagePlus("ball", new ByteProcessor(ball.patchwidth+1, ball.patchwidth+1, ball.data, null)).show();
92  //ImageProcessor smallImage = ip.resize(ip.getWidth()/ball.shrinkfactor, ip.getHeight()/ball.shrinkfactor);
93  ImageProcessor smallImage = shrinkImage(ip, ball.shrinkfactor);
94  //new ImagePlus("small image", smallImage).show();
95  if (slice==1)
96  IJ.showStatus("Rolling ball ("+ball.shrinkfactor+")...");
97  ImageProcessor background = rollBall(ball, ip, smallImage);
98  interpolateBackground(background, ball);
99  extrapolateBackground(background, ball);
100  IJ.showProgress(0.9);
101  if (IJ.altKeyDown())
102  new ImagePlus("background", background).show();
103  ip.copyBits(background, 0, 0, Blitter.SUBTRACT);
104  if (invert)
105  ip.invert();
106  IJ.showProgress(1.0);
107  }
108 
121  ImageProcessor rollBall(RollingBall ball, ImageProcessor image, ImageProcessor smallImage) {
122 
123  int halfpatchwidth; //distance in x or y from patch center to any edge
124  int ptsbelowlastpatch; //number of points we may ignore because they were below last patch
125  int xpt2, ypt2; // current (x,y) point in the patch relative to upper left corner
126  int xval, yval; // location in ball in shrunken image coordinates
127  int zdif; // difference in z (height) between point on ball and point on image
128  int zmin; // smallest zdif for ball patch with center at current point
129  int zctr; // current height of the center of the sphere of which the patch is a part
130  int zadd; // height of a point on patch relative to the xy-plane of the shrunken image
131  int ballpt; // index to array storing the precomputed ball patch
132  int imgpt; // index to array storing the shrunken image
133  int backgrpt; // index to array storing the calculated background
134  int ybackgrpt; // displacement to current background scan line
135  int p1, p2; // temporary indexes to background, ball, or small image
136  int ybackgrinc; // distance in memory between two shrunken y-points in background
137  int smallimagewidth; // length of a scan line in shrunken image
138  int left, right, top, bottom;
139  byte[] pixels = (byte[])smallImage.getPixels();
140  byte[] patch = ball.data;
141  int width = image.getWidth();
142  int height = image.getHeight();
143  int swidth = smallImage.getWidth();
144  int sheight = smallImage.getHeight();
145  ImageProcessor background = new ByteProcessor(width, height);
146  byte[] backgroundpixels = (byte[])background.getPixels();
147  int shrinkfactor = ball.shrinkfactor;
148  int leftroll = 0;
149  int rightroll = width/shrinkfactor-1;
150  int toproll = 0;
151  int bottomroll = height/shrinkfactor-1;
152 
153  left = 1;
154  right = rightroll - leftroll - 1;
155  top = 1;
156  bottom = bottomroll - toproll - 1;
157  smallimagewidth = swidth;
158  int patchwidth = ball.patchwidth;
159  halfpatchwidth = patchwidth / 2;
160  ybackgrinc = shrinkfactor*width; // real dist btwn 2 adjacent (dy=1) shrunk pts
161  zctr = 0; // start z-center in the xy-plane
162  for (int ypt=top; ypt<=(bottom+patchwidth); ypt++) {
163  for (int xpt=left; xpt<=(right+patchwidth); xpt++) {// while patch is tangent to edges or within image...
164  // xpt is far right edge of ball patch
165  // do we have to move the patch up or down to make it tangent to but not above image?...
166  zmin = 255; // highest could ever be 255
167  ballpt = 0;
168  ypt2 = ypt - patchwidth; // ypt2 is top edge of ball patch
169  imgpt = ypt2*smallimagewidth + xpt - patchwidth;
170  while (ypt2<=ypt) {
171  xpt2 = xpt-patchwidth; // xpt2 is far left edge of ball patch
172  while (xpt2<=xpt) { // check every point on ball patch
173  // only examine points on
174  if ((xpt2>=left) && (xpt2<=right) && (ypt2>=top) && (ypt2<=bottom)) {
175  p1 = ballpt;
176  p2 = imgpt;
177  zdif = (pixels[p2]&255) - (zctr + (patch[p1]&255)); //curve - circle points
178  //if (xpt==50 && ypt==50) IJ.write(zdif
179  //+" "+(pixels[p2]&255)
180  //+" "+zctr
181  //+" "+(patch[p1]&255)
182  //+" "+p2
183  //+" "+p1
184  //);
185  if (zdif<zmin) // keep most negative, since ball should always be below curve
186  zmin = zdif;
187  } // if xpt2,ypt2
188  ballpt++;
189  xpt2++;
190  imgpt++;
191  } // while xpt2
192  ypt2++;
193  imgpt = imgpt - patchwidth - 1 + smallimagewidth;
194  } // while ypt2
195  if (zmin!=0)
196  zctr += zmin; // move ball up or down if we find a new minimum
197  if (zmin<0)
198  ptsbelowlastpatch = halfpatchwidth; // ignore left half of ball patch when dz < 0
199  else
200  ptsbelowlastpatch = 0;
201  // now compare every point on ball with background, and keep highest number
202  yval = ypt - patchwidth;
203  ypt2 = 0;
204  ballpt = 0;
205  ybackgrpt = (yval - top + 1) * ybackgrinc;
206  while (ypt2<=patchwidth) {
207  xval = xpt - patchwidth + ptsbelowlastpatch;
208  xpt2 = ptsbelowlastpatch;
209  ballpt += ptsbelowlastpatch;
210  backgrpt = ybackgrpt + (xval - left + 1) * shrinkfactor;
211  while (xpt2<=patchwidth) { // for all the points in the ball patch
212  if ((xval >= left) && (xval <= right) && (yval >= top) && (yval <= bottom)) {
213  p1 = ballpt;
214  zadd = zctr + (patch[p1]&255);
215  p1 = backgrpt;
216  //if (backgrpt>=backgroundpixels.length) backgrpt = 0; //(debug)
217  if (zadd>(backgroundpixels[p1]&255)) //keep largest adjustment}
218  backgroundpixels[p1] = (byte)zadd;
219  }
220  ballpt++;
221  xval++;
222  xpt2++;
223  backgrpt += shrinkfactor; // move to next point in x
224  } // while xpt2
225  yval++;
226  ypt2++;
227  ybackgrpt += ybackgrinc; // move to next point in y
228  } // while ypt2
229  } // for xpt
230  if (ypt%20==0)
231  IJ.showProgress(0.2+0.6*ypt/(bottom+patchwidth));
232  } // for ypt
233  return background;
234  }
235 
237  ImageProcessor shrinkImage(ImageProcessor ip, int shrinkfactor) {
238  int width = ip.getWidth();
239  int height = ip.getHeight();
240  int swidth = width/shrinkfactor;
241  int sheight = height/shrinkfactor;
242  ImageProcessor ip2 = ip.duplicate();
243  ip2.smooth();
244  IJ.showProgress(0.1);
245  ImageProcessor smallImage = new ByteProcessor(swidth, sheight);
246  byte[] pixels = (byte[])ip2.getPixels();
247  byte[] spixels = (byte[])smallImage.getPixels();
248  int xmaskmin, ymaskmin, min, nextrowoffset, paddr, thispixel;
249  for (int i=0; i<(swidth*sheight); i++) {
250  xmaskmin = shrinkfactor*(i%swidth);
251  ymaskmin = shrinkfactor*(i/swidth);
252  min = 255;
253  nextrowoffset = width-shrinkfactor;
254  paddr = ymaskmin*width+xmaskmin;
255  for (int j=1; j<=shrinkfactor; j++) {
256  for (int k=1; k<=shrinkfactor; k++) {
257  thispixel = pixels[paddr++]&255;
258  if (thispixel<min)
259  min = thispixel;
260  }
261  paddr += nextrowoffset;
262  }
263  spixels[i] = (byte)min; // each point in small image is minimum of its neighborhood
264  }
265  return smallImage;
266  }
267 
274  void interpolateBackground(ImageProcessor background, RollingBall ball) {
275  int hloc, vloc; // position of current pixel in calculated background
276  int vinc; // memory offset from current calculated pos to current interpolated pos
277  int lastvalue, nextvalue; // calculated pixel values between which we are interpolating
278  int p; // pointer to current interpolated pixel value
279  int bglastptr, bgnextptr; // pointers to calculated pixel values between which we are interpolating
280 
281  int width = background.getWidth();
282  int height = background.getHeight();
283  int shrinkfactor = ball.shrinkfactor;
284  int leftroll = 0;
285  int rightroll = width/shrinkfactor-1;
286  int toproll = 0;
287  int bottomroll = height/shrinkfactor-1;
288  byte[] pixels = (byte[])background.getPixels();
289 
290  vloc = 0;
291  for (int j=1; j<=(bottomroll-toproll-1); j++) { //interpolate to find background interior
292  hloc = 0;
293  vloc += shrinkfactor;
294  for (int i=1; i<=(rightroll-leftroll); i++) {
295  hloc += shrinkfactor;
296  bgnextptr = vloc*width + hloc;
297  bglastptr = bgnextptr - shrinkfactor;
298  nextvalue = pixels[bgnextptr]&255;
299  lastvalue = pixels[bglastptr]&255;
300  for (int ii=1; ii<=(shrinkfactor-1); ii++) { //interpolate horizontally
301  p = bgnextptr - ii;
302  pixels[p] = (byte)(lastvalue+(shrinkfactor-ii)*(nextvalue-lastvalue)/shrinkfactor);
303  }
304  for (int ii=0; ii<=(shrinkfactor-1); ii++) { //interpolate vertically
305  bglastptr = (vloc-shrinkfactor)*width+hloc-ii;
306  bgnextptr = vloc*width+hloc-ii;
307  lastvalue = pixels[bglastptr]&255;
308  nextvalue = pixels[bgnextptr]&255;
309  vinc = 0;
310  for (int jj=1; jj<=(shrinkfactor-1); jj++) {
311  vinc = vinc-width;
312  p = bgnextptr+vinc;
313  pixels[p] = (byte)(lastvalue+(shrinkfactor-jj)*(nextvalue-lastvalue)/shrinkfactor);
314  } // for jj
315  } // for ii
316  } // for i
317  } // for j
318  }
319 
328  void extrapolateBackground(ImageProcessor background, RollingBall ball) {
329 
330  int edgeslope; // difference of last two consecutive pixel values on an edge
331  int pvalue; // current extrapolated pixel value
332  int lastvalue, nextvalue; //calculated pixel values from which we are extrapolating
333  int p; // pointer to current extrapolated pixel value
334  int bglastptr, bgnextptr; // pointers to calculated pixel values from which we are extrapolating
335 
336  int width = background.getWidth();
337  int height = background.getHeight();
338  int shrinkfactor = ball.shrinkfactor;
339  int leftroll = 0;
340  int rightroll = width/shrinkfactor-1;
341  int toproll = 0;
342  int bottomroll = height/shrinkfactor-1;
343  byte[] pixels = (byte[])background.getPixels();
344 
345  for (int hloc=shrinkfactor; hloc<=(shrinkfactor*(rightroll-leftroll)-1); hloc++) {
346  // extrapolate on top and bottom
347  bglastptr = shrinkfactor*width+hloc;
348  bgnextptr = (shrinkfactor+1)*width+hloc;
349  lastvalue = pixels[bglastptr]&255;
350  nextvalue = pixels[bgnextptr]&255;
351  edgeslope = nextvalue-lastvalue;
352  p = bglastptr;
353  pvalue = lastvalue;
354  for (int jj=1; jj<=shrinkfactor; jj++) {
355  p = p-width;
356  pvalue = pvalue-edgeslope;
357  if (pvalue<0)
358  pixels[p] = 0;
359  else if (pvalue>255)
360  pixels[p] = (byte)255;
361  else
362  pixels[p] = (byte)pvalue;
363  } // for jj
364  bglastptr = (shrinkfactor*(bottomroll-toproll-1)-1)*width+hloc;
365  bgnextptr = shrinkfactor*(bottomroll-toproll-1)*width+hloc;
366  lastvalue = pixels[bglastptr]&255;
367  nextvalue = pixels[bgnextptr]&255;
368  edgeslope = nextvalue-lastvalue;
369  p = bgnextptr;
370  pvalue = nextvalue;
371  for (int jj=1; jj<=((height-1)-shrinkfactor*(bottomroll-toproll-1)); jj++) {
372  p += width;
373  pvalue += edgeslope;
374  if (pvalue<0)
375  pixels[p] = 0;
376  else if (pvalue>255)
377  pixels[p] = (byte)255;
378  else
379  pixels[p] = (byte)pvalue;
380  } // for jj
381  } // for hloc
382  for (int vloc=0; vloc<height; vloc++) {
383  // extrapolate on left and right
384  bglastptr = vloc*width+shrinkfactor;
385  bgnextptr = bglastptr+1;
386  lastvalue = pixels[bglastptr]&255;
387  nextvalue = pixels[bgnextptr]&255;
388  edgeslope = nextvalue-lastvalue;
389  p = bglastptr;
390  pvalue = lastvalue;
391  for (int ii=1; ii<=shrinkfactor; ii++) {
392  p--;
393  pvalue = pvalue - edgeslope;
394  if (pvalue<0)
395  pixels[p] = 0;
396  else if (pvalue>255)
397  pixels[p] = (byte)255;
398  else
399  pixels[p] = (byte)pvalue;
400  } // for ii
401  bgnextptr = vloc*width+shrinkfactor*(rightroll-leftroll-1)-1;
402  bglastptr = bgnextptr-1;
403  lastvalue = pixels[bglastptr]&255;
404  nextvalue = pixels[bgnextptr]&255;
405  edgeslope = nextvalue-lastvalue;
406  p = bgnextptr;
407  pvalue = nextvalue;
408  for (int ii=1; ii<=((width-1)-shrinkfactor*(rightroll-leftroll-1)+1); ii++) {
409  p++;
410  pvalue = pvalue+edgeslope;
411  if (pvalue<0)
412  pixels[p] = 0;
413  else if (pvalue>255)
414  pixels[p] = (byte)255;
415  else
416  pixels[p] = (byte)pvalue;
417  } // for ii
418  } // for vloc
419  }
420 
421 }
422 
423 
424 class RollingBall {
425 
426  byte[] data;
427  int patchwidth;
428  int shrinkfactor;
429 
430  RollingBall(int radius) {
431  int arcTrimPer;
432  if (radius<=10) {
433  shrinkfactor = 1;
434  arcTrimPer = 12; // trim 24% in x and y
435  } else if (radius<=30) {
436  shrinkfactor = 2;
437  arcTrimPer = 12; // trim 24% in x and y
438  } else if (radius<=100) {
439  shrinkfactor = 4;
440  arcTrimPer = 16; // trim 32% in x and y
441  } else {
442  shrinkfactor = 8;
443  arcTrimPer = 20; // trim 40% in x and y
444  }
445  buildRollingBall(radius, arcTrimPer);
446  }
447 
453  void buildRollingBall(int ballradius, int arcTrimPer) {
454  int rsquare; // rolling ball radius squared
455  int xtrim; // # of pixels trimmed off each end of ball to make patch
456  int xval, yval; // x,y-values on patch relative to center of rolling ball
457  int smallballradius, diam; // radius and diameter of rolling ball
458  int temp; // value must be >=0 to take square root
459  int halfpatchwidth; // distance in x or y from center of patch to any edge
460  int ballsize; // size of rolling ball array
461 
462  this.shrinkfactor = shrinkfactor;
463  smallballradius = ballradius/shrinkfactor;
464  if (smallballradius<1)
465  smallballradius = 1;
466  rsquare = smallballradius*smallballradius;
467  diam = smallballradius*2;
468  xtrim = (arcTrimPer*diam)/100; // only use a patch of the rolling ball
469  patchwidth = diam - xtrim - xtrim;
470  halfpatchwidth = smallballradius - xtrim;
471  ballsize = (patchwidth+1)*(patchwidth+1);
472  data = new byte[ballsize];
473 
474  for (int i=0; i<ballsize; i++) {
475  xval = i % (patchwidth+1) - halfpatchwidth;
476  yval = i / (patchwidth+1) - halfpatchwidth;
477  temp = rsquare - (xval*xval) - (yval*yval);
478  if (temp >= 0)
479  data[i] = (byte)Math.round(Math.sqrt(temp));
480  else
481  data[i] = 0;
482  }
483 
484  }
485 }
486