source: tags/release-1.1.12/src/ATrous/atrous_transform.cc

Last change on this file was 641, checked in by MatthewWhiting, 15 years ago

Some changes to the development code, with some image plotting functions, and fixing namespaces for the atrous transform code.
File size: 16.5 KB
Line 
1#include <iostream>
2#include <math.h>
3#include <duchamp/param.hh>
4#include <duchamp/ATrous/atrous.hh>
5#include <duchamp/Utils/utils.hh>
6
7namespace duchamp{
8
9  /***********************************************************************/
10  /////  1-DIMENSIONAL TRANSFORM
11  /***********************************************************************/
12
13  // template <class T>
14  // void atrousTransform(long &length, T *&spectrum, T *&coeffs, T *&wavelet)
15  // {
16  //   int filterHW = filterwidth/2;
17  //   int numScales = getNumScales(length);
18
19  //   delete [] coeffs;
20  //   coeffs = new T[(numScales+1)*length];
21  //   delete [] wavelet;
22  //   wavelet = new T[(numScales+1)*length];
23  //   for(int i=0;i<length;i++)  coeffs[i] = wavelet[i] = spectrum[i];
24
25  //   int spacing = 1;
26  //   for(int scale = 0; scale<numScales; scale++){
27
28  //     for(int pos = 0; pos<length; pos++){
29  //       coeffs[(scale+1)*length+pos] = 0;
30  //       for(int offset=-filterHW; offset<=filterHW; offset++){
31  //    int x = pos + spacing*offset;
32  //    if(x<0) x = -x;                    // boundary conditions are
33  //    if(x>=length) x = 2*(length-1) - x;    //    reflection.
34  // //         if(x<0) x = x+length;                // boundary conditions are
35  // //         if(x>=length) x = x-length;            //    continuous.
36
37  //    coeffs[(scale+1)*length+pos] += filter1D[offset+filterHW] *
38  //      coeffs[scale*length+x];
39  //       }
40  //       wavelet[(scale+1)*length+pos] = coeffs[scale*length+pos] -
41  //    coeffs[(scale+1)*length+pos];
42  //     }
43
44  //     spacing *= 2;
45
46  //   }
47
48  // }
49
50 
51  void atrousTransform(long &length, int &numScales, float *spectrum, double *coeffs, double *wavelet, duchamp::Param &par)
52  {
53    duchamp::Filter reconFilter = par.filter();
54    int filterHW = reconFilter.width()/2;
55
56    for(int i=0;i<length;i++)  coeffs[i] = wavelet[i] = spectrum[i];
57
58    int spacing = 1;
59    for(int scale = 0; scale<numScales; scale++){
60
61      for(int pos = 0; pos<length; pos++){
62        coeffs[(scale+1)*length+pos] = 0;
63        for(int offset=-filterHW; offset<=filterHW; offset++){
64          int x = pos + spacing*offset;
65          if(x<0) x = -x;                    // boundary conditions are
66          if(x>=length) x = 2*(length-1) - x;    //    reflection.
67          //    if(x<0) x = x+length;                // boundary conditions are
68          //    if(x>=length) x = x-length;            //    continuous.
69
70          //    coeffs[(scale+1)*length+pos] += filter1D[offset+filterHW] *
71          //      coeffs[scale*length+x];
72          coeffs[(scale+1)*length+pos] += reconFilter.coeff(offset+filterHW) *
73            coeffs[scale*length+x];
74        }
75        wavelet[(scale+1)*length+pos] = coeffs[scale*length+pos] -
76          coeffs[(scale+1)*length+pos];
77      }
78
79      spacing *= 2;
80
81    }
82
83  }
84
85  void atrousTransform(long &length, float *spectrum, float *coeffs, float *wavelet, duchamp::Param &par)
86  {
87    duchamp::Filter reconFilter = par.filter();
88    int filterHW = reconFilter.width()/2;
89    int numScales = reconFilter.getNumScales(length);
90
91    delete [] coeffs;
92    coeffs = new float[(numScales+1)*length];
93    delete [] wavelet;
94    wavelet = new float[(numScales+1)*length];
95    for(int i=0;i<length;i++) coeffs[i] = wavelet[i] = spectrum[i];
96
97    int spacing = 1;
98    for(int scale = 0; scale<numScales; scale++){
99
100      for(int pos = 0; pos<length; pos++){
101        coeffs[(scale+1)*length+pos] = 0;
102        for(int offset=-filterHW; offset<=filterHW; offset++){
103          int x = pos + spacing*offset;
104          if(x<0) x = -x;                    // boundary conditions are
105          if(x>=length) x = 2*(length-1) - x;    //    reflection.
106          //    if(x<0) x = x+length;                // boundary conditions are
107          //    if(x>=length) x = x-length;            //    continuous.
108
109          //    coeffs[(scale+1)*length+pos] += filter1D[offset+filterHW] *
110          //      coeffs[scale*length+x];
111          coeffs[(scale+1)*length+pos] += reconFilter.coeff(offset+filterHW) *
112            coeffs[scale*length+x];
113        }
114        wavelet[(scale+1)*length+pos] = coeffs[scale*length+pos] -
115          coeffs[(scale+1)*length+pos];
116      }
117
118      spacing *= 2;
119
120    }
121
122  }
123
124
125
126
127  /***********************************************************************/
128  /////  2-DIMENSIONAL TRANSFORM
129  /***********************************************************************/
130
131  void atrousTransform2D(long &xdim, long &ydim, int &numScales, float *input, double *coeffs, double *wavelet, duchamp::Param &par)
132  {
133    duchamp::Filter reconFilter = par.filter();
134    float blankPixValue = par.getBlankPixVal();
135    int filterHW = reconFilter.width()/2;
136
137    double *filter = new double[reconFilter.width()*reconFilter.width()];
138    for(int i=0;i<reconFilter.width();i++){
139      for(int j=0;j<reconFilter.width();j++){
140        filter[i*reconFilter.width()+j] = reconFilter.coeff(i) * reconFilter.coeff(j);
141      }
142    }
143
144    long size = xdim * ydim;
145    float *oldcoeffs = new float[size];
146
147    // locating the borders of the image -- ignoring BLANK pixels
148    //   int xLim1=0, yLim1=0, xLim2=xdim-1, yLim2=ydim-1;
149    //   for(int row=0;row<ydim;row++){
150    //     while((xLim1<xLim2)&&(input[row*xdim+xLim1]==blankPixValue)) xLim1++;
151    //     while((xLim2>xLim1)&&(input[row*xdim+xLim1]==blankPixValue)) xLim2--;
152    //   }
153    //   for(int col=0;col<xdim;col++){
154    //     while((yLim1<yLim2)&&(input[col+xdim*yLim1]==blankPixValue)) yLim1++;
155    //     while((yLim2>yLim1)&&(input[col+xdim*yLim1]==blankPixValue)) yLim2--;
156    //   }
157    //   std::cerr << "X Limits: "<<xLim1<<" "<<xLim2<<std::endl;
158    //   std::cerr << "Y Limits: "<<yLim1<<" "<<yLim2<<std::endl;
159 
160    int *xLim1 = new int[ydim];
161    int *yLim1 = new int[xdim];
162    int *xLim2 = new int[ydim];
163    int *yLim2 = new int[xdim];
164    for(int row=0;row<ydim;row++){
165      int ct1 = 0;
166      int ct2 = xdim - 1;
167      while((ct1<ct2)&&(input[row*xdim+ct1]==blankPixValue) ) ct1++;
168      while((ct2>ct1)&&(input[row*xdim+ct2]==blankPixValue) ) ct2--;
169      xLim1[row] = ct1;
170      xLim2[row] = ct2;
171      std::cerr<<row<<":"<<xLim1[row]<<","<<xLim2[row]<<" ";
172    }
173    std::cerr<<std::endl;
174
175    for(int col=0;col<xdim;col++){
176      int ct1=0;
177      int ct2=ydim-1;
178      while((ct1<ct2)&&(input[col+xdim*ct1]==blankPixValue) ) ct1++;
179      while((ct2>ct1)&&(input[col+xdim*ct2]==blankPixValue) ) ct2--;
180      yLim1[col] = ct1;
181      yLim2[col] = ct2;
182    }
183
184
185    bool *isGood = new bool[size];
186    for(int pos=0;pos<size;pos++) //isGood[pos] = (!flagBlank) || (input[pos]!=blankPixValue);
187      isGood[pos] = !par.isBlank(input[pos]);
188
189    for(int i=0;i<size;i++)  coeffs[i] = wavelet[i] = input[i];
190
191    int spacing = 1;
192    for(int scale = 0; scale<numScales; scale++){
193
194      for(int i=0;i<size;i++) oldcoeffs[i] = coeffs[i];
195
196      //std::cerr << numScales<<" "<<scale<<" "<<spacing<<" "<<reconFilter.width()*spacing<<std::endl;
197      for(int ypos = 0; ypos<ydim; ypos++){
198        for(int xpos = 0; xpos<xdim; xpos++){
199          // loops over each pixel in the image
200          int pos = ypos*xdim + xpos;
201          coeffs[pos] = 0;
202
203          //    if((par.getFlagBlankPix())&&(oldcoeffs[pos] == blankPixValue) )
204          if(par.isBlank(oldcoeffs[pos]) )
205            coeffs[pos] = oldcoeffs[pos];
206          else{
207            for(int yoffset=-filterHW; yoffset<=filterHW; yoffset++){
208              int y = ypos + spacing*yoffset;
209              int newy;
210              //          if(y<0) y = -y;                 // boundary conditions are
211              //          if(y>=ydim) y = 2*(ydim-1) - y; //    reflection.
212              //          while((y<yLim1)||(y>yLim2)){
213              //            if(y<yLim1) y = 2*yLim1 - y;      // boundary conditions are
214              //            if(y>yLim2) y = 2*yLim2 - y;      //    reflection.
215              //          }
216              // boundary conditions are reflection.
217         
218              for(int xoffset=-filterHW; xoffset<=filterHW; xoffset++){
219                int x = xpos + spacing*xoffset;
220                int newx;
221                //if(x<0) x = -x;                 // boundary conditions are
222                // if(x>=xdim) x = 2*(xdim-1) - x; //    reflection.
223                //while((x<xLim1)||(x>xLim2)){
224                //            if(x<xLim1) x = 2*xLim1 - x;      // boundary conditions are
225                //            if(x>xLim2) x = 2*xLim2 - x;      //    reflection.
226                //          }
227                // boundary conditions are reflection.
228                if(y<yLim1[xpos]) newy = 2*yLim1[xpos] - y;     
229                else if(y>yLim2[xpos]) newy = 2*yLim2[xpos] - y;     
230                else newy = y;
231                if(x<xLim1[ypos]) newx = 2*xLim1[ypos] - x;
232                else if(x>xLim2[ypos]) newx = 2*xLim2[ypos] - x;     
233                else newx=x;     
234             
235                x = newx;
236                y = newy;
237
238                int filterpos = (yoffset+filterHW)*reconFilter.width() + (xoffset+filterHW);
239                int oldpos = y*xdim + x;
240
241                if(// (x>=0)&&(x<xdim)&&(y>=0)&&(y<ydim)&&
242                   (isGood[pos]))
243                  coeffs[pos] += filter[filterpos] * oldcoeffs[oldpos];
244              }
245            }
246     
247          }
248          wavelet[(scale+1)*size+pos] = oldcoeffs[pos] - coeffs[pos];
249
250
251        }
252      }
253 
254      spacing *= 2;
255
256    }
257
258    delete [] filter;
259    delete [] oldcoeffs;
260
261  }
262
263  // void atrousTransform2D(long &xdim, long &ydim, int &numScales, float *input, double *coeffs, double *wavelet)
264  // {
265  //   Filter reconFilter = par.filter();
266  //   int filterHW = reconFilter.width()/2;
267
268  //   double *filter = new double[reconFilter.width()*reconFilter.width()];
269  //   for(int i=0;i<reconFilter.width();i++){
270  //     for(int j=0;j<reconFilter.width();j++){
271  //       filter[i*reconFilter.width()+j] = reconFilter.coeff(i) * reconFilter.coeff(j);
272  //     }
273  //   }
274
275  //   long size = xdim * ydim;
276  //   float *oldcoeffs = new float[size];
277
278  //   for(int i=0;i<size;i++)  coeffs[i] = wavelet[i] = input[i];
279
280 
281  //   int spacing = 1;
282  //   for(int scale = 0; scale<numScales; scale++){
283
284  //     for(int i=0;i<size;i++) oldcoeffs[i] = coeffs[i];
285
286  //     std::cerr << numScales<<" "<<scale<<" "<<spacing<<std::endl;
287  //     for(int ypos = 0; ypos<ydim; ypos++){
288  //       for(int xpos = 0; xpos<xdim; xpos++){
289  //    // loops over each pixel in the image
290  //    int pos = ypos*xdim + xpos;
291  //    coeffs[pos] = 0;
292
293  //    for(int yoffset=-filterHW; yoffset<=filterHW; yoffset++){
294  //      int y = ypos + spacing*yoffset;
295  //      if(y<0) y = -y;                 // boundary conditions are
296  //      if(y>=ydim) y = 2*(ydim-1) - y; //    reflection.
297         
298  //      for(int xoffset=-filterHW; xoffset<=filterHW; xoffset++){
299  //        int x = xpos + spacing*xoffset;
300  //        if(x<0) x = -x;                 // boundary conditions are
301  //        if(x>=xdim) x = 2*(xdim-1) - x; //    reflection.
302
303  //        int filterpos = (yoffset+filterHW)*reconFilter.width() + (xoffset+filterHW);
304  //        int oldpos = y*xdim + x;
305  //        coeffs[pos] += filter[filterpos] * oldcoeffs[oldpos];
306  //      }
307  //    }
308       
309  //    wavelet[(scale+1)*size+pos] = oldcoeffs[pos] - coeffs[pos];
310
311  //       }
312  //     }
313 
314  //     spacing *= 2;
315
316  //   }
317
318  //   delete [] filter;
319  //   delete [] oldcoeffs;
320
321  // }
322
323  /***********************************************************************/
324  /////  3-DIMENSIONAL TRANSFORM
325  /***********************************************************************/
326
327  void atrousTransform3D(long &xdim, long &ydim, long &zdim, int &numScales, float *&input, float *&coeffs, float *&wavelet, duchamp::Param &par)
328  {
329    duchamp::Filter reconFilter = par.filter();
330    float blankPixValue = par.getBlankPixVal();
331    int filterHW = reconFilter.width()/2;
332
333    long size = xdim * ydim * zdim;
334    float *oldcoeffs = new float[size];
335
336    std::cerr << "%";
337    int fsize = reconFilter.width()*reconFilter.width()*reconFilter.width();
338    std::cerr << "%";
339    double *filter = new double[fsize];
340    for(int i=0;i<reconFilter.width();i++){
341      for(int j=0;j<reconFilter.width();j++){
342        for(int k=0;k<reconFilter.width();k++){
343          filter[i +j*reconFilter.width() + k*reconFilter.width()*reconFilter.width()] =
344            reconFilter.coeff(i) * reconFilter.coeff(j) * reconFilter.coeff(k);
345        }
346      }
347    }
348    // locating the borders of the image -- ignoring BLANK pixels
349    // HAVE NOT DONE THIS FOR Z --> ASSUMING NO TRIMMING IN SPECTRAL DIRECTION
350    int xLim1 = 0, yLim1 = 0, xLim2 = xdim-1, yLim2 = ydim-1;
351    for(int col=0;col<xdim;col++){
352      while((yLim1<yLim2)&&(input[col+xdim*yLim1]==blankPixValue) ) yLim1++;
353      while((yLim2>yLim1)&&(input[col+xdim*yLim1]==blankPixValue) ) yLim2--;
354    }
355    for(int row=0;row<ydim;row++){
356      while((xLim1<xLim2)&&(input[row*xdim+xLim1]==blankPixValue) ) xLim1++;
357      while((xLim2>xLim1)&&(input[row*xdim+xLim1]==blankPixValue) ) xLim2--;
358    }
359
360    for(int i=0;i<size;i++)  coeffs[i] = wavelet[i] = input[i];
361
362    int spacing = 1;
363    std::cerr<<xdim<<"x"<<ydim<<"x"<<zdim<<"x"<<numScales;
364    for(int scale = 0; scale<numScales; scale++){
365      std::cerr << ".";
366      for(int i=0;i<size;i++) oldcoeffs[i] = coeffs[i];
367
368      for(int zpos = 0; zpos<zdim; zpos++){
369        for(int ypos = 0; ypos<ydim; ypos++){
370          for(int xpos = 0; xpos<xdim; xpos++){
371            // loops over each pixel in the image
372            int pos = zpos*xdim*ydim + ypos*xdim + xpos;
373            coeffs[pos] = 0;
374
375            if(par.isBlank(oldcoeffs[pos]) )
376              coeffs[pos] = oldcoeffs[pos];
377            else{
378              for(int zoffset=-filterHW; zoffset<=filterHW; zoffset++){
379                int z = zpos + spacing*zoffset;
380                if(z<0) z = -z;                 // boundary conditions are
381                if(z>=zdim) z = 2*(zdim-1) - z; //    reflection.
382         
383                for(int yoffset=-filterHW; yoffset<=filterHW; yoffset++){
384                  int y = ypos + spacing*yoffset;
385                  //if(y<0) y = -y;                 // boundary conditions are
386                  // if(y>=ydim) y = 2*(ydim-1) - y; //    reflection.
387                  if(y<yLim1) y = 2*yLim1 - y;      // boundary conditions are
388                  if(y>yLim2) y = 2*yLim2 - y;      //    reflection.
389         
390                  for(int xoffset=-filterHW; xoffset<=filterHW; xoffset++){
391                    int x = xpos + spacing*xoffset;
392                    //if(x<0) x = -x;                 // boundary conditions are
393                    //if(x>=xdim) x = 2*(xdim-1) - x; //    reflection.
394                    if(x<xLim1) x = 2*xLim1 - x;      // boundary conditions are
395                    if(x>xLim2) x = 2*xLim2 - x;      //    reflection.
396
397                    int filterpos = (zoffset+filterHW)*reconFilter.width()*reconFilter.width() +
398                      (yoffset+filterHW)*reconFilter.width() + (xoffset+filterHW);
399                    int oldpos = z*xdim*ydim + y*xdim + x;
400
401                    if(!par.isBlank(oldcoeffs[oldpos]))
402                      coeffs[pos] += filter[filterpos] * oldcoeffs[oldpos];
403
404                  }
405                }
406              }
407            }           
408 
409            wavelet[(scale+1)*size+pos] = oldcoeffs[pos] - coeffs[pos];
410
411          }
412        }
413      }
414 
415      spacing *= 2;
416
417    }
418    std::cerr << "|";
419
420    delete [] filter;
421    delete [] oldcoeffs;
422
423  }
424
425
426  // void atrousTransform3D(long &xdim, long &ydim, long &zdim, int &numScales, float *input, float *coeffs, float *wavelet)
427  // {
428  //   extern Filter reconFilter;
429  //   int filterHW = reconFilter.width()/2;
430
431  //   long size = xdim * ydim * zdim;
432  //   float *oldcoeffs = new float[size];
433
434  //   std::cerr << "%";
435  //   int fsize = reconFilter.width()*reconFilter.width()*reconFilter.width();
436  //   std::cerr << "%";
437  //   double *filter = new double[fsize];
438  //   for(int i=0;i<reconFilter.width();i++){
439  //     for(int j=0;j<reconFilter.width();j++){
440  //       for(int k=0;k<reconFilter.width();k++){
441  //       filter[i +j*reconFilter.width() + k*reconFilter.width()*reconFilter.width()] =
442  //    reconFilter.coeff(i) * reconFilter.coeff(j) * reconFilter.coeff(k);
443  //       }
444  //     }
445  //   }
446
447  //   for(int i=0;i<size;i++)  coeffs[i] = wavelet[i] = input[i];
448
449  //   int spacing = 1;
450  //   std::cerr<<xdim<<"x"<<ydim<<"x"<<zdim<<"x"<<numScales;
451  //   for(int scale = 0; scale<numScales; scale++){
452  //     std::cerr << ".";
453  //     for(int i=0;i<size;i++) oldcoeffs[i] = coeffs[i];
454
455  //     for(int zpos = 0; zpos<zdim; zpos++){
456  //       for(int ypos = 0; ypos<ydim; ypos++){
457  //    for(int xpos = 0; xpos<xdim; xpos++){
458  //      // loops over each pixel in the image
459  //      int pos = zpos*xdim*ydim + ypos*xdim + xpos;
460  //      coeffs[pos] = 0;
461
462  //      for(int zoffset=-filterHW; zoffset<=filterHW; zoffset++){
463  //        int z = zpos + spacing*zoffset;
464  //        if(z<0) z = -z;                 // boundary conditions are
465  //        if(z>=zdim) z = 2*(zdim-1) - z; //    reflection.
466         
467  //        for(int yoffset=-filterHW; yoffset<=filterHW; yoffset++){
468  //          int y = ypos + spacing*yoffset;
469  //          if(y<0) y = -y;                 // boundary conditions are
470  //          if(y>=ydim) y = 2*(ydim-1) - y; //    reflection.
471         
472  //          for(int xoffset=-filterHW; xoffset<=filterHW; xoffset++){
473  //            int x = xpos + spacing*xoffset;
474  //            if(x<0) x = -x;                 // boundary conditions are
475  //            if(x>=xdim) x = 2*(xdim-1) - x; //    reflection.
476
477  //            int filterpos = (zoffset+filterHW)*reconFilter.width()*reconFilter.width() +
478  //              (yoffset+filterHW)*reconFilter.width() + (xoffset+filterHW);
479  //            int oldpos = z*xdim*ydim + y*xdim + x;
480
481  //            coeffs[pos] += filter[filterpos] * oldcoeffs[oldpos];
482
483  //          }
484  //        }
485  //      }
486               
487  //      wavelet[(scale+1)*size+pos] = oldcoeffs[pos] - coeffs[pos];
488
489  //    }
490  //       }
491  //     }
492 
493  //     spacing *= 2;
494
495  //   }
496  //   std::cerr << "|";
497
498  //   delete [] filter;
499  //   delete [] oldcoeffs;
500
501  // }
502
503
504}
Note: See TracBrowser for help on using the repository browser.