source: trunk/src/ATrous/atrous_transform.cc @ 861

Last change on this file since 861 was 861, checked in by MatthewWhiting, 13 years ago

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