source: branches/fitshead-branch/ATrous/atrous_transform.cc @ 1441

Last change on this file since 1441 was 3, checked in by Matthew Whiting, 18 years ago

This is the first full import of all working code to
the Duchamp repository.
Made three directories at top level:

branches/ tags/ trunk/

and trunk/ has the full set of code:
ATrous/ Cubes/ Detection/ InputComplete? InputExample? README Utils/ docs/ mainDuchamp.cc param.cc param.hh
File size: 15.6 KB
Line 
1#include <iostream>
2#include <math.h>
3#include <ATrous/atrous.hh>
4//#include <Cubes/cubes.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.