Squiz Matrix  4.12.2
 All Data Structures Namespaces Functions Variables Pages
FFTFilter.java
1 package ij.plugin.filter;
2 import ij.*;
3 import ij.process.*;
4 import ij.gui.*;
5 import ij.measure.*;
6 import ij.plugin.ContrastEnhancer;
7 import java.awt.*;
8 import java.util.*;
9 
10 
15 public class FFTFilter implements PlugInFilter, Measurements {
16 
17  private ImagePlus imp;
18  private String arg;
19  private static int filterIndex = 1;
20  private FHT fht;
21  private int slice;
22  private int stackSize = 1;
23 
24  private static double filterLargeDia = 40.0;
25  private static double filterSmallDia = 3.0;
26  private static int choiceIndex = 0;
27  private static String[] choices = {"None","Horizontal","Vertical"};
28  private static String choiceDia = choices[0];
29  private static double toleranceDia = 5.0;
30  private static boolean doScalingDia = true;
31  private static boolean saturateDia = true;
32  private static boolean displayFilter;
33  private static boolean processStack;
34 
35  public int setup(String arg, ImagePlus imp) {
36  this.arg = arg;
37  this.imp = imp;
38  if (imp==null)
39  {IJ.noImage(); return DONE;}
40  stackSize = imp.getStackSize();
41  fht = (FHT)imp.getProperty("FHT");
42  if (fht!=null) {
43  IJ.showMessage("FFT Filter", "Spatial domain image required");
44  return DONE;
45  }
46  if (!showBandpassDialog(imp))
47  return DONE;
48  else
49  return processStack?DOES_ALL+DOES_STACKS:DOES_ALL;
50  }
51 
52  public void run(ImageProcessor ip) {
53  slice++;
54  filter(ip);
55  }
56 
57  void filter(ImageProcessor ip) {
58  ImageProcessor ip2 = ip;
59  if (ip2 instanceof ColorProcessor) {
60  showStatus("Extracting brightness");
61  ip2 = ((ColorProcessor)ip2).getBrightness();
62  }
63  Rectangle roiRect = ip2.getRoi();
64  int maxN = Math.max(roiRect.width, roiRect.height);
65  // scale filterLarge and filterSmall to fractions of image size
66  double filterLarge = filterLargeDia / maxN;
67  double filterSmall = filterSmallDia / maxN;
68  double sharpness = (100.0 - toleranceDia) / 100.0;
69  boolean doScaling = doScalingDia;
70  boolean saturate = saturateDia;
71 
72  IJ.showProgress(1,20);
73 
74  /* tile mirrored image to power of 2 size
75  first determine smallest power 2 >= 1.5 * image width/height
76  factor of 1.5 to avoid wrap-around effects of Fourier Trafo */
77 
78  int i=2;
79  while(i<1.5 * maxN) i *= 2;
80 
81  // fit image into power of 2 size
82  Rectangle fitRect = new Rectangle();
83  fitRect.x = (int) Math.round( (i - roiRect.width) / 2.0 );
84  fitRect.y = (int) Math.round( (i - roiRect.height) / 2.0 );
85  fitRect.width = roiRect.width;
86  fitRect.height = roiRect.height;
87 
88  // put image (ROI) into power 2 size image
89  // mirroring to avoid wrap around effects
90  showStatus("Pad to "+i+"x"+i);
91  ip2 = tileMirror(ip2, i, i, fitRect.x, fitRect.y);
92  IJ.showProgress(2,20);
93 
94  // transform forward
95  showStatus(i+"x"+i+" forward transform");
96  FHT fht = new FHT(ip2);
97  fht.setShowProgress(false);
98  fht.transform();
99  System.gc();
100  IJ.showProgress(9,20);
101  //new ImagePlus("after fht",ip2.crop()).show();
102 
103  // filter out large and small structures
104  showStatus("Filter in frequency domain");
105  filterLargeSmall(fht, filterLarge, filterSmall, choiceIndex, sharpness);
106  //new ImagePlus("filter",ip2.crop()).show();
107  IJ.showProgress(11,20);
108 
109  // transform backward
110  showStatus("Inverse transform");
111  fht.inverseTransform();
112  IJ.showProgress(19,20);
113  //new ImagePlus("after inverse",ip2).show();
114 
115  // crop to original size and do scaling if selected
116  showStatus("Crop and convert to original type");
117  fht.setRoi(fitRect);
118  ip2 = fht.crop();
119  if (doScaling) {
120  ImagePlus imp2 = new ImagePlus(imp.getTitle()+"-filtered", ip2);
121  new ContrastEnhancer().stretchHistogram(imp2, saturate?1.0:0.0);
122  ip2 = imp2.getProcessor();
123  }
124 
125  // convert back to original data type
126  int bitDepth = imp.getBitDepth();
127  switch (bitDepth) {
128  case 8: ip2 = ip2.convertToByte(doScaling); break;
129  case 16: ip2 = ip2.convertToShort(doScaling); break;
130  case 24:
131  ip.snapshot();
132  showStatus("Setting brightness");
133  ((ColorProcessor)ip).setBrightness((FloatProcessor)ip2);
134  break;
135  case 32: break;
136  }
137 
138  // copy filtered image back into original image
139  if (bitDepth!=24) {
140  ip.snapshot();
141  ip.copyBits(ip2, roiRect.x, roiRect.y, Blitter.COPY);
142  }
143  ip.resetMinAndMax();
144  System.gc();
145  IJ.showProgress(20,20);
146  }
147 
148  void showStatus(String msg) {
149  if (stackSize>1 && processStack)
150  IJ.showStatus("FFT Filter: "+slice+"/"+stackSize);
151  else
152  IJ.showStatus(msg);
153  }
154 
157  ImageProcessor tileMirror(ImageProcessor ip, int width, int height, int x, int y) {
158 
159  if (x < 0 || x > (width -1) || y < 0 || y > (height -1)) {
160  IJ.error("Image to be tiled is out of bounds.");
161  return null;
162  }
163 
164  ImageProcessor ipout = ip.createProcessor(width, height);
165 
166  ImageProcessor ip2 = ip.crop();
167  int w2 = ip2.getWidth();
168  int h2 = ip2.getHeight();
169 
170  //how many times does ip2 fit into ipout?
171  int i1 = (int) Math.ceil(x / (double) w2);
172  int i2 = (int) Math.ceil( (width - x) / (double) w2);
173  int j1 = (int) Math.ceil(y / (double) h2);
174  int j2 = (int) Math.ceil( (height - y) / (double) h2);
175 
176  //tile
177  if ( (i1%2) > 0.5)
178  ip2.flipHorizontal();
179  if ( (j1%2) > 0.5)
180  ip2.flipVertical();
181 
182  for (int i=-i1; i<i2; i += 2) {
183  for (int j=-j1; j<j2; j += 2) {
184  ipout.insert(ip2, x-i*w2, y-j*h2);
185  }
186  }
187 
188  ip2.flipHorizontal();
189  for (int i=-i1+1; i<i2; i += 2) {
190  for (int j=-j1; j<j2; j += 2) {
191  ipout.insert(ip2, x-i*w2, y-j*h2);
192  }
193  }
194 
195  ip2.flipVertical();
196  for (int i=-i1+1; i<i2; i += 2) {
197  for (int j=-j1+1; j<j2; j += 2) {
198  ipout.insert(ip2, x-i*w2, y-j*h2);
199  }
200  }
201 
202  ip2.flipHorizontal();
203  for (int i=-i1; i<i2; i += 2) {
204  for (int j=-j1+1; j<j2; j += 2) {
205  ipout.insert(ip2, x-i*w2, y-j*h2);
206  }
207  }
208 
209  return ipout;
210  }
211 
212 
213  /*
214  filterLarge: down to which size are large structures suppressed?
215  filterSmall: up to which size are small structures suppressed?
216  filterLarge and filterSmall are given as fraction of the image size
217  in the original (untransformed) image.
218  stripesHorVert: filter out: 0) nothing more 1) horizontal 2) vertical stripes
219  (i.e. frequencies with x=0 / y=0)
220  scaleStripes: width of the stripe filter, same unit as filterLarge
221  */
222  void filterLargeSmall(ImageProcessor ip, double filterLarge, double filterSmall, int stripesHorVert, double scaleStripes) {
223 
224  int maxN = ip.getWidth();
225 
226  float[] fht = (float[])ip.getPixels();
227  float[] filter = new float[maxN*maxN];
228  for (int i=0; i<maxN*maxN; i++)
229  filter[i]=1f;
230 
231  int row;
232  int backrow;
233  float rowFactLarge;
234  float rowFactSmall;
235 
236  int col;
237  int backcol;
238  float factor;
239  float colFactLarge;
240  float colFactSmall;
241 
242  float factStripes;
243 
244  // calculate factor in exponent of Gaussian from filterLarge / filterSmall
245 
246  double scaleLarge = filterLarge*filterLarge;
247  double scaleSmall = filterSmall*filterSmall;
248  scaleStripes = scaleStripes*scaleStripes;
249  //float FactStripes;
250 
251  // loop over rows
252  for (int j=1; j<maxN/2; j++) {
253  row = j * maxN;
254  backrow = (maxN-j)*maxN;
255  rowFactLarge = (float) Math.exp(-(j*j) * scaleLarge);
256  rowFactSmall = (float) Math.exp(-(j*j) * scaleSmall);
257 
258 
259  // loop over columns
260  for (col=1; col<maxN/2; col++){
261  backcol = maxN-col;
262  colFactLarge = (float) Math.exp(- (col*col) * scaleLarge);
263  colFactSmall = (float) Math.exp(- (col*col) * scaleSmall);
264  factor = (1 - rowFactLarge*colFactLarge) * rowFactSmall*colFactSmall;
265  switch (stripesHorVert) {
266  case 1: factor *= (1 - (float) Math.exp(- (col*col) * scaleStripes)); break;// hor stripes
267  case 2: factor *= (1 - (float) Math.exp(- (j*j) * scaleStripes)); // vert stripes
268  }
269 
270  fht[col+row] *= factor;
271  fht[col+backrow] *= factor;
272  fht[backcol+row] *= factor;
273  fht[backcol+backrow] *= factor;
274  filter[col+row] *= factor;
275  filter[col+backrow] *= factor;
276  filter[backcol+row] *= factor;
277  filter[backcol+backrow] *= factor;
278  }
279  }
280 
281  //process meeting points (maxN/2,0) , (0,maxN/2), and (maxN/2,maxN/2)
282  int rowmid = maxN * (maxN/2);
283  rowFactLarge = (float) Math.exp(- (maxN/2)*(maxN/2) * scaleLarge);
284  rowFactSmall = (float) Math.exp(- (maxN/2)*(maxN/2) * scaleSmall);
285  factStripes = (float) Math.exp(- (maxN/2)*(maxN/2) * scaleStripes);
286 
287  fht[maxN/2] *= (1 - rowFactLarge) * rowFactSmall; // (maxN/2,0)
288  fht[rowmid] *= (1 - rowFactLarge) * rowFactSmall; // (0,maxN/2)
289  fht[maxN/2 + rowmid] *= (1 - rowFactLarge*rowFactLarge) * rowFactSmall*rowFactSmall; // (maxN/2,maxN/2)
290  filter[maxN/2] *= (1 - rowFactLarge) * rowFactSmall; // (maxN/2,0)
291  filter[rowmid] *= (1 - rowFactLarge) * rowFactSmall; // (0,maxN/2)
292  filter[maxN/2 + rowmid] *= (1 - rowFactLarge*rowFactLarge) * rowFactSmall*rowFactSmall; // (maxN/2,maxN/2)
293 
294  switch (stripesHorVert) {
295  case 1: fht[maxN/2] *= (1 - factStripes);
296  fht[rowmid] = 0;
297  fht[maxN/2 + rowmid] *= (1 - factStripes);
298  filter[maxN/2] *= (1 - factStripes);
299  filter[rowmid] = 0;
300  filter[maxN/2 + rowmid] *= (1 - factStripes);
301  break; // hor stripes
302  case 2: fht[maxN/2] = 0;
303  fht[rowmid] *= (1 - factStripes);
304  fht[maxN/2 + rowmid] *= (1 - factStripes);
305  filter[maxN/2] = 0;
306  filter[rowmid] *= (1 - factStripes);
307  filter[maxN/2 + rowmid] *= (1 - factStripes);
308  break; // vert stripes
309  }
310 
311  //loop along row 0 and maxN/2
312  rowFactLarge = (float) Math.exp(- (maxN/2)*(maxN/2) * scaleLarge);
313  rowFactSmall = (float) Math.exp(- (maxN/2)*(maxN/2) * scaleSmall);
314  for (col=1; col<maxN/2; col++){
315  backcol = maxN-col;
316  colFactLarge = (float) Math.exp(- (col*col) * scaleLarge);
317  colFactSmall = (float) Math.exp(- (col*col) * scaleSmall);
318 
319  switch (stripesHorVert) {
320  case 0:
321  fht[col] *= (1 - colFactLarge) * colFactSmall;
322  fht[backcol] *= (1 - colFactLarge) * colFactSmall;
323  fht[col+rowmid] *= (1 - colFactLarge*rowFactLarge) * colFactSmall*rowFactSmall;
324  fht[backcol+rowmid] *= (1 - colFactLarge*rowFactLarge) * colFactSmall*rowFactSmall;
325  filter[col] *= (1 - colFactLarge) * colFactSmall;
326  filter[backcol] *= (1 - colFactLarge) * colFactSmall;
327  filter[col+rowmid] *= (1 - colFactLarge*rowFactLarge) * colFactSmall*rowFactSmall;
328  filter[backcol+rowmid] *= (1 - colFactLarge*rowFactLarge) * colFactSmall*rowFactSmall;
329  break;
330  case 1:
331  factStripes = (float) Math.exp(- (col*col) * scaleStripes);
332  fht[col] *= (1 - colFactLarge) * colFactSmall * (1 - factStripes);
333  fht[backcol] *= (1 - colFactLarge) * colFactSmall * (1 - factStripes);
334  fht[col+rowmid] *= (1 - colFactLarge*rowFactLarge) * colFactSmall*rowFactSmall * (1 - factStripes);
335  fht[backcol+rowmid] *= (1 - colFactLarge*rowFactLarge) * colFactSmall*rowFactSmall * (1 - factStripes);
336  filter[col] *= (1 - colFactLarge) * colFactSmall * (1 - factStripes);
337  filter[backcol] *= (1 - colFactLarge) * colFactSmall * (1 - factStripes);
338  filter[col+rowmid] *= (1 - colFactLarge*rowFactLarge) * colFactSmall*rowFactSmall * (1 - factStripes);
339  filter[backcol+rowmid] *= (1 - colFactLarge*rowFactLarge) * colFactSmall*rowFactSmall * (1 - factStripes);
340  break;
341  case 2:
342  factStripes = (float) Math.exp(- (maxN/2)*(maxN/2) * scaleStripes);
343  fht[col] = 0;
344  fht[backcol] = 0;
345  fht[col+rowmid] *= (1 - colFactLarge*rowFactLarge) * colFactSmall*rowFactSmall * (1 - factStripes);
346  fht[backcol+rowmid] *= (1 - colFactLarge*rowFactLarge) * colFactSmall*rowFactSmall * (1 - factStripes);
347  filter[col] = 0;
348  filter[backcol] = 0;
349  filter[col+rowmid] *= (1 - colFactLarge*rowFactLarge) * colFactSmall*rowFactSmall * (1 - factStripes);
350  filter[backcol+rowmid] *= (1 - colFactLarge*rowFactLarge) * colFactSmall*rowFactSmall * (1 - factStripes);
351  }
352  }
353 
354  // loop along column 0 and maxN/2
355  colFactLarge = (float) Math.exp(- (maxN/2)*(maxN/2) * scaleLarge);
356  colFactSmall = (float) Math.exp(- (maxN/2)*(maxN/2) * scaleSmall);
357  for (int j=1; j<maxN/2; j++) {
358  row = j * maxN;
359  backrow = (maxN-j)*maxN;
360  rowFactLarge = (float) Math.exp(- (j*j) * scaleLarge);
361  rowFactSmall = (float) Math.exp(- (j*j) * scaleSmall);
362 
363  switch (stripesHorVert) {
364  case 0:
365  fht[row] *= (1 - rowFactLarge) * rowFactSmall;
366  fht[backrow] *= (1 - rowFactLarge) * rowFactSmall;
367  fht[row+maxN/2] *= (1 - rowFactLarge*colFactLarge) * rowFactSmall*colFactSmall;
368  fht[backrow+maxN/2] *= (1 - rowFactLarge*colFactLarge) * rowFactSmall*colFactSmall;
369  filter[row] *= (1 - rowFactLarge) * rowFactSmall;
370  filter[backrow] *= (1 - rowFactLarge) * rowFactSmall;
371  filter[row+maxN/2] *= (1 - rowFactLarge*colFactLarge) * rowFactSmall*colFactSmall;
372  filter[backrow+maxN/2] *= (1 - rowFactLarge*colFactLarge) * rowFactSmall*colFactSmall;
373  break;
374  case 1:
375  factStripes = (float) Math.exp(- (maxN/2)*(maxN/2) * scaleStripes);
376  fht[row] = 0;
377  fht[backrow] = 0;
378  fht[row+maxN/2] *= (1 - rowFactLarge*colFactLarge) * rowFactSmall*colFactSmall * (1 - factStripes);
379  fht[backrow+maxN/2] *= (1 - rowFactLarge*colFactLarge) * rowFactSmall*colFactSmall * (1 - factStripes);
380  filter[row] = 0;
381  filter[backrow] = 0;
382  filter[row+maxN/2] *= (1 - rowFactLarge*colFactLarge) * rowFactSmall*colFactSmall * (1 - factStripes);
383  filter[backrow+maxN/2] *= (1 - rowFactLarge*colFactLarge) * rowFactSmall*colFactSmall * (1 - factStripes);
384  break;
385  case 2:
386  factStripes = (float) Math.exp(- (j*j) * scaleStripes);
387  fht[row] *= (1 - rowFactLarge) * rowFactSmall * (1 - factStripes);
388  fht[backrow] *= (1 - rowFactLarge) * rowFactSmall * (1 - factStripes);
389  fht[row+maxN/2] *= (1 - rowFactLarge*colFactLarge) * rowFactSmall*colFactSmall * (1 - factStripes);
390  fht[backrow+maxN/2] *= (1 - rowFactLarge*colFactLarge) * rowFactSmall*colFactSmall * (1 - factStripes);
391  filter[row] *= (1 - rowFactLarge) * rowFactSmall * (1 - factStripes);
392  filter[backrow] *= (1 - rowFactLarge) * rowFactSmall * (1 - factStripes);
393  filter[row+maxN/2] *= (1 - rowFactLarge*colFactLarge) * rowFactSmall*colFactSmall * (1 - factStripes);
394  filter[backrow+maxN/2] *= (1 - rowFactLarge*colFactLarge) * rowFactSmall*colFactSmall * (1 - factStripes);
395  }
396  }
397  if (displayFilter && slice==1) {
398  FHT f = new FHT(new FloatProcessor(maxN, maxN, filter, null));
399  f.swapQuadrants();
400  new ImagePlus("Filter", f).show();
401  }
402  }
403 
404  boolean showBandpassDialog(ImagePlus imp) {
405  GenericDialog gd = new GenericDialog("FFT Bandpass Filter");
406  gd.addNumericField("Filter_Large Structures Down to", filterLargeDia, 0, 4, "pixels");
407  gd.addNumericField("Filter_Small Structures Up to", filterSmallDia, 0, 4, "pixels");
408  gd.addChoice("Suppress Stripes:", choices, choiceDia);
409  gd.addNumericField("Tolerance of Direction:", toleranceDia, 0, 2, "%");
410  gd.addCheckbox("Autoscale After Filtering", doScalingDia);
411  gd.addCheckbox("Saturate Image when Autoscaling", saturateDia);
412  gd.addCheckbox("Display Filter", displayFilter);
413  if (stackSize>1)
414  gd.addCheckbox("Process Entire Stack", processStack);
415  gd.showDialog();
416  if(gd.wasCanceled())
417  return false;
418  if(gd.invalidNumber()) {
419  IJ.showMessage("Error", "Invalid input number");
420  return false;
421  }
422  filterLargeDia = gd.getNextNumber();
423  filterSmallDia = gd.getNextNumber();
424  choiceIndex = gd.getNextChoiceIndex();
425  choiceDia = choices[choiceIndex];
426  toleranceDia = gd.getNextNumber();
427  doScalingDia = gd.getNextBoolean();
428  saturateDia = gd.getNextBoolean();
429  displayFilter = gd.getNextBoolean();
430  if (stackSize>1)
431  processStack = gd.getNextBoolean();
432  return true;
433  }
434 
435 }
436