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

Last change on this file since 393 was 393, checked in by MatthewWhiting, 17 years ago

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