source: branches/pixel-map-branch/src/ATrous/atrous_transform.cc

Last change on this file was 253, checked in by Matthew Whiting, 17 years ago
  • Added a mergeList(vector<Scan>) function to the PixelMap? namespace. This is usefu

l for combining lists of Scans into a set of independent ones.

  • Wrote a new spectralSelection function to deal with selection of a spectrum for development applications.
File size: 15.9 KB
Line 
1#include <iostream>
2#include <math.h>
3#include <param.hh>
4#include <ATrous/atrous.hh>
5#include <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, Param &par)
50{
51  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, Param &par)
84{
85  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, Param &par)
130{
131  Filter reconFilter = par.filter();
132  bool flagBlank=par.getFlagBlankPix();
133  float blankPixValue = par.getBlankPixVal();
134  int filterHW = reconFilter.width()/2;
135
136  double *filter = new double[reconFilter.width()*reconFilter.width()];
137  for(int i=0;i<reconFilter.width();i++){
138    for(int j=0;j<reconFilter.width();j++){
139      filter[i*reconFilter.width()+j] = reconFilter.coeff(i) * reconFilter.coeff(j);
140    }
141  }
142
143  long size = xdim * ydim;
144  float *oldcoeffs = new float[size];
145
146  // locating the borders of the image -- ignoring BLANK pixels
147//   int xLim1=0, yLim1=0, xLim2=xdim-1, yLim2=ydim-1;
148//   for(int row=0;row<ydim;row++){
149//     while((xLim1<xLim2)&&(input[row*xdim+xLim1]==blankPixValue)) xLim1++;
150//     while((xLim2>xLim1)&&(input[row*xdim+xLim1]==blankPixValue)) xLim2--;
151//   }
152//   for(int col=0;col<xdim;col++){
153//     while((yLim1<yLim2)&&(input[col+xdim*yLim1]==blankPixValue)) yLim1++;
154//     while((yLim2>yLim1)&&(input[col+xdim*yLim1]==blankPixValue)) yLim2--;
155//   }
156//   std::cerr << "X Limits: "<<xLim1<<" "<<xLim2<<std::endl;
157//   std::cerr << "Y Limits: "<<yLim1<<" "<<yLim2<<std::endl;
158 
159  int *xLim1 = new int[ydim];
160  int *yLim1 = new int[xdim];
161  int *xLim2 = new int[ydim];
162  int *yLim2 = new int[xdim];
163  for(int row=0;row<ydim;row++){
164    int ct1 = 0;
165    int ct2 = xdim - 1;
166    while((ct1<ct2)&&(input[row*xdim+ct1]==blankPixValue) ) ct1++;
167    while((ct2>ct1)&&(input[row*xdim+ct2]==blankPixValue) ) ct2--;
168    xLim1[row] = ct1;
169    xLim2[row] = ct2;
170    std::cerr<<row<<":"<<xLim1[row]<<","<<xLim2[row]<<" ";
171  }
172  std::cerr<<std::endl;
173
174  for(int col=0;col<xdim;col++){
175    int ct1=0;
176    int ct2=ydim-1;
177    while((ct1<ct2)&&(input[col+xdim*ct1]==blankPixValue) ) ct1++;
178    while((ct2>ct1)&&(input[col+xdim*ct2]==blankPixValue) ) ct2--;
179    yLim1[col] = ct1;
180    yLim2[col] = ct2;
181  }
182
183
184  bool *isGood = new bool[size];
185  for(int pos=0;pos<size;pos++) //isGood[pos] = (!flagBlank) || (input[pos]!=blankPixValue);
186    isGood[pos] = !par.isBlank(input[pos]);
187
188  for(int i=0;i<size;i++)  coeffs[i] = wavelet[i] = input[i];
189
190  int spacing = 1;
191  for(int scale = 0; scale<numScales; scale++){
192
193    for(int i=0;i<size;i++) oldcoeffs[i] = coeffs[i];
194
195    //std::cerr << numScales<<" "<<scale<<" "<<spacing<<" "<<reconFilter.width()*spacing<<std::endl;
196    for(int ypos = 0; ypos<ydim; ypos++){
197      for(int xpos = 0; xpos<xdim; xpos++){
198        // loops over each pixel in the image
199        int pos = ypos*xdim + xpos;
200        coeffs[pos] = 0;
201
202//      if((par.getFlagBlankPix())&&(oldcoeffs[pos] == blankPixValue) )
203        if(par.isBlank(oldcoeffs[pos]) )
204          coeffs[pos] = oldcoeffs[pos];
205        else{
206          for(int yoffset=-filterHW; yoffset<=filterHW; yoffset++){
207            int y = ypos + spacing*yoffset;
208            int newy;
209            //    if(y<0) y = -y;                 // boundary conditions are
210            //    if(y>=ydim) y = 2*(ydim-1) - y; //    reflection.
211            //    while((y<yLim1)||(y>yLim2)){
212//          if(y<yLim1) y = 2*yLim1 - y;      // boundary conditions are
213//          if(y>yLim2) y = 2*yLim2 - y;      //    reflection.
214            //    }
215              // boundary conditions are reflection.
216         
217            for(int xoffset=-filterHW; xoffset<=filterHW; xoffset++){
218              int x = xpos + spacing*xoffset;
219              int newx;
220              //if(x<0) x = -x;                 // boundary conditions are
221              // if(x>=xdim) x = 2*(xdim-1) - x; //    reflection.
222              //while((x<xLim1)||(x>xLim2)){
223//            if(x<xLim1) x = 2*xLim1 - x;      // boundary conditions are
224//            if(x>xLim2) x = 2*xLim2 - x;      //    reflection.
225              //            }
226              // boundary conditions are reflection.
227              if(y<yLim1[xpos]) newy = 2*yLim1[xpos] - y;     
228              else if(y>yLim2[xpos]) newy = 2*yLim2[xpos] - y;     
229              else newy = y;
230              if(x<xLim1[ypos]) newx = 2*xLim1[ypos] - x;
231              else if(x>xLim2[ypos]) newx = 2*xLim2[ypos] - x;     
232              else newx=x;     
233             
234              x = newx;
235              y = newy;
236
237              int filterpos = (yoffset+filterHW)*reconFilter.width() + (xoffset+filterHW);
238              int oldpos = y*xdim + x;
239
240              if(// (x>=0)&&(x<xdim)&&(y>=0)&&(y<ydim)&&
241                 (isGood[pos]))
242                coeffs[pos] += filter[filterpos] * oldcoeffs[oldpos];
243            }
244          }
245     
246        }
247        wavelet[(scale+1)*size+pos] = oldcoeffs[pos] - coeffs[pos];
248
249
250      }
251    }
252 
253    spacing *= 2;
254
255  }
256
257  delete [] filter;
258  delete [] oldcoeffs;
259
260}
261
262// void atrousTransform2D(long &xdim, long &ydim, int &numScales, float *input, double *coeffs, double *wavelet)
263// {
264//   Filter reconFilter = par.filter();
265//   int filterHW = reconFilter.width()/2;
266
267//   double *filter = new double[reconFilter.width()*reconFilter.width()];
268//   for(int i=0;i<reconFilter.width();i++){
269//     for(int j=0;j<reconFilter.width();j++){
270//       filter[i*reconFilter.width()+j] = reconFilter.coeff(i) * reconFilter.coeff(j);
271//     }
272//   }
273
274//   long size = xdim * ydim;
275//   float *oldcoeffs = new float[size];
276
277//   for(int i=0;i<size;i++)  coeffs[i] = wavelet[i] = input[i];
278
279 
280//   int spacing = 1;
281//   for(int scale = 0; scale<numScales; scale++){
282
283//     for(int i=0;i<size;i++) oldcoeffs[i] = coeffs[i];
284
285//     std::cerr << numScales<<" "<<scale<<" "<<spacing<<std::endl;
286//     for(int ypos = 0; ypos<ydim; ypos++){
287//       for(int xpos = 0; xpos<xdim; xpos++){
288//      // loops over each pixel in the image
289//      int pos = ypos*xdim + xpos;
290//      coeffs[pos] = 0;
291
292//      for(int yoffset=-filterHW; yoffset<=filterHW; yoffset++){
293//        int y = ypos + spacing*yoffset;
294//        if(y<0) y = -y;                 // boundary conditions are
295//        if(y>=ydim) y = 2*(ydim-1) - y; //    reflection.
296         
297//        for(int xoffset=-filterHW; xoffset<=filterHW; xoffset++){
298//          int x = xpos + spacing*xoffset;
299//          if(x<0) x = -x;                 // boundary conditions are
300//          if(x>=xdim) x = 2*(xdim-1) - x; //    reflection.
301
302//          int filterpos = (yoffset+filterHW)*reconFilter.width() + (xoffset+filterHW);
303//          int oldpos = y*xdim + x;
304//          coeffs[pos] += filter[filterpos] * oldcoeffs[oldpos];
305//        }
306//      }
307       
308//      wavelet[(scale+1)*size+pos] = oldcoeffs[pos] - coeffs[pos];
309
310//       }
311//     }
312 
313//     spacing *= 2;
314
315//   }
316
317//   delete [] filter;
318//   delete [] oldcoeffs;
319
320// }
321
322/***********************************************************************/
323/////  3-DIMENSIONAL TRANSFORM
324/***********************************************************************/
325
326void atrousTransform3D(long &xdim, long &ydim, long &zdim, int &numScales, float *&input, float *&coeffs, float *&wavelet, Param &par)
327{
328  Filter reconFilter = par.filter();
329  float blankPixValue = par.getBlankPixVal();
330  int filterHW = reconFilter.width()/2;
331
332  long size = xdim * ydim * zdim;
333  float *oldcoeffs = new float[size];
334
335  std::cerr << "%";
336  int fsize = reconFilter.width()*reconFilter.width()*reconFilter.width();
337  std::cerr << "%";
338  double *filter = new double[fsize];
339  for(int i=0;i<reconFilter.width();i++){
340    for(int j=0;j<reconFilter.width();j++){
341      for(int k=0;k<reconFilter.width();k++){
342      filter[i +j*reconFilter.width() + k*reconFilter.width()*reconFilter.width()] =
343        reconFilter.coeff(i) * reconFilter.coeff(j) * reconFilter.coeff(k);
344      }
345    }
346  }
347  // locating the borders of the image -- ignoring BLANK pixels
348  // HAVE NOT DONE THIS FOR Z --> ASSUMING NO TRIMMING IN SPECTRAL DIRECTION
349  int xLim1 = 0, yLim1 = 0, xLim2 = xdim-1, yLim2 = ydim-1;
350  for(int col=0;col<xdim;col++){
351    while((yLim1<yLim2)&&(input[col+xdim*yLim1]==blankPixValue) ) yLim1++;
352    while((yLim2>yLim1)&&(input[col+xdim*yLim1]==blankPixValue) ) yLim2--;
353  }
354  for(int row=0;row<ydim;row++){
355    while((xLim1<xLim2)&&(input[row*xdim+xLim1]==blankPixValue) ) xLim1++;
356    while((xLim2>xLim1)&&(input[row*xdim+xLim1]==blankPixValue) ) xLim2--;
357  }
358
359  for(int i=0;i<size;i++)  coeffs[i] = wavelet[i] = input[i];
360
361  int spacing = 1;
362  std::cerr<<xdim<<"x"<<ydim<<"x"<<zdim<<"x"<<numScales;
363  for(int scale = 0; scale<numScales; scale++){
364    std::cerr << ".";
365    for(int i=0;i<size;i++) oldcoeffs[i] = coeffs[i];
366
367    for(int zpos = 0; zpos<zdim; zpos++){
368      for(int ypos = 0; ypos<ydim; ypos++){
369        for(int xpos = 0; xpos<xdim; xpos++){
370          // loops over each pixel in the image
371          int pos = zpos*xdim*ydim + ypos*xdim + xpos;
372          coeffs[pos] = 0;
373
374          if(par.isBlank(oldcoeffs[pos]) )
375            coeffs[pos] = oldcoeffs[pos];
376          else{
377            for(int zoffset=-filterHW; zoffset<=filterHW; zoffset++){
378              int z = zpos + spacing*zoffset;
379              if(z<0) z = -z;                 // boundary conditions are
380              if(z>=zdim) z = 2*(zdim-1) - z; //    reflection.
381         
382              for(int yoffset=-filterHW; yoffset<=filterHW; yoffset++){
383                int y = ypos + spacing*yoffset;
384                //if(y<0) y = -y;                 // boundary conditions are
385                // if(y>=ydim) y = 2*(ydim-1) - y; //    reflection.
386                if(y<yLim1) y = 2*yLim1 - y;      // boundary conditions are
387                if(y>yLim2) y = 2*yLim2 - y;      //    reflection.
388         
389                for(int xoffset=-filterHW; xoffset<=filterHW; xoffset++){
390                  int x = xpos + spacing*xoffset;
391                  //if(x<0) x = -x;                 // boundary conditions are
392                  //if(x>=xdim) x = 2*(xdim-1) - x; //    reflection.
393                  if(x<xLim1) x = 2*xLim1 - x;      // boundary conditions are
394                  if(x>xLim2) x = 2*xLim2 - x;      //    reflection.
395
396                  int filterpos = (zoffset+filterHW)*reconFilter.width()*reconFilter.width() +
397                    (yoffset+filterHW)*reconFilter.width() + (xoffset+filterHW);
398                  int oldpos = z*xdim*ydim + y*xdim + x;
399
400                 if(!par.isBlank(oldcoeffs[oldpos]))
401                    coeffs[pos] += filter[filterpos] * oldcoeffs[oldpos];
402
403                }
404              }
405            }
406          }             
407 
408          wavelet[(scale+1)*size+pos] = oldcoeffs[pos] - coeffs[pos];
409
410        }
411      }
412    }
413 
414    spacing *= 2;
415
416  }
417  std::cerr << "|";
418
419  delete [] filter;
420  delete [] oldcoeffs;
421
422}
423
424
425// void atrousTransform3D(long &xdim, long &ydim, long &zdim, int &numScales, float *input, float *coeffs, float *wavelet)
426// {
427//   extern Filter reconFilter;
428//   int filterHW = reconFilter.width()/2;
429
430//   long size = xdim * ydim * zdim;
431//   float *oldcoeffs = new float[size];
432
433//   std::cerr << "%";
434//   int fsize = reconFilter.width()*reconFilter.width()*reconFilter.width();
435//   std::cerr << "%";
436//   double *filter = new double[fsize];
437//   for(int i=0;i<reconFilter.width();i++){
438//     for(int j=0;j<reconFilter.width();j++){
439//       for(int k=0;k<reconFilter.width();k++){
440//       filter[i +j*reconFilter.width() + k*reconFilter.width()*reconFilter.width()] =
441//      reconFilter.coeff(i) * reconFilter.coeff(j) * reconFilter.coeff(k);
442//       }
443//     }
444//   }
445
446//   for(int i=0;i<size;i++)  coeffs[i] = wavelet[i] = input[i];
447
448//   int spacing = 1;
449//   std::cerr<<xdim<<"x"<<ydim<<"x"<<zdim<<"x"<<numScales;
450//   for(int scale = 0; scale<numScales; scale++){
451//     std::cerr << ".";
452//     for(int i=0;i<size;i++) oldcoeffs[i] = coeffs[i];
453
454//     for(int zpos = 0; zpos<zdim; zpos++){
455//       for(int ypos = 0; ypos<ydim; ypos++){
456//      for(int xpos = 0; xpos<xdim; xpos++){
457//        // loops over each pixel in the image
458//        int pos = zpos*xdim*ydim + ypos*xdim + xpos;
459//        coeffs[pos] = 0;
460
461//        for(int zoffset=-filterHW; zoffset<=filterHW; zoffset++){
462//          int z = zpos + spacing*zoffset;
463//          if(z<0) z = -z;                 // boundary conditions are
464//          if(z>=zdim) z = 2*(zdim-1) - z; //    reflection.
465         
466//          for(int yoffset=-filterHW; yoffset<=filterHW; yoffset++){
467//            int y = ypos + spacing*yoffset;
468//            if(y<0) y = -y;                 // boundary conditions are
469//            if(y>=ydim) y = 2*(ydim-1) - y; //    reflection.
470         
471//            for(int xoffset=-filterHW; xoffset<=filterHW; xoffset++){
472//              int x = xpos + spacing*xoffset;
473//              if(x<0) x = -x;                 // boundary conditions are
474//              if(x>=xdim) x = 2*(xdim-1) - x; //    reflection.
475
476//              int filterpos = (zoffset+filterHW)*reconFilter.width()*reconFilter.width() +
477//                (yoffset+filterHW)*reconFilter.width() + (xoffset+filterHW);
478//              int oldpos = z*xdim*ydim + y*xdim + x;
479
480//              coeffs[pos] += filter[filterpos] * oldcoeffs[oldpos];
481
482//            }
483//          }
484//        }
485               
486//        wavelet[(scale+1)*size+pos] = oldcoeffs[pos] - coeffs[pos];
487
488//      }
489//       }
490//     }
491 
492//     spacing *= 2;
493
494//   }
495//   std::cerr << "|";
496
497//   delete [] filter;
498//   delete [] oldcoeffs;
499
500// }
501
502
Note: See TracBrowser for help on using the repository browser.