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

Last change on this file since 236 was 233, checked in by Matthew Whiting, 17 years ago

Other than minor changes to code comments, major change is adding atrous_transform.cc back for development purposes. Also updating menus.cc and devel.hh.
File size: 15.5 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)
50{
51  extern Filter reconFilter;
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)
84{
85  extern Filter reconFilter;
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  extern Filter reconFilter;
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
262void atrousTransform2D(long &xdim, long &ydim, int &numScales, float *input, double *coeffs, double *wavelet)
263{
264  extern Filter reconFilter;
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  extern Filter reconFilter;
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
425void 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.