source: tags/release-1.2.2/src/ATrous/atrous_transform.cc

Last change on this file was 1069, checked in by MatthewWhiting, 12 years ago

Adding the licensing text to files that didn't have it.
File size: 18.0 KB
Line 
1// -----------------------------------------------------------------------
2// atrous_transform.cc: Simplified a trous transform functions - only used
3//                      for testing purposes.
4// -----------------------------------------------------------------------
5// Copyright (C) 2006, Matthew Whiting, ATNF
6//
7// This program is free software; you can redistribute it and/or modify it
8// under the terms of the GNU General Public License as published by the
9// Free Software Foundation; either version 2 of the License, or (at your
10// option) any later version.
11//
12// Duchamp is distributed in the hope that it will be useful, but WITHOUT
13// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14// FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
15// for more details.
16//
17// You should have received a copy of the GNU General Public License
18// along with Duchamp; if not, write to the Free Software Foundation,
19// Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
20//
21// Correspondence concerning Duchamp may be directed to:
22//    Internet email: Matthew.Whiting [at] atnf.csiro.au
23//    Postal address: Dr. Matthew Whiting
24//                    Australia Telescope National Facility, CSIRO
25//                    PO Box 76
26//                    Epping NSW 1710
27//                    AUSTRALIA
28// -----------------------------------------------------------------------
29#include <iostream>
30#include <math.h>
31#include <duchamp/param.hh>
32#include <duchamp/ATrous/atrous.hh>
33#include <duchamp/Utils/utils.hh>
34#include <duchamp/Devel/devel.hh>
35
36namespace duchamp{
37
38  /***********************************************************************/
39  /////  1-DIMENSIONAL TRANSFORM
40  /***********************************************************************/
41
42  // template <class T>
43  // void atrousTransform(long &length, T *&spectrum, T *&coeffs, T *&wavelet)
44  // {
45  //   int filterHW = filterwidth/2;
46  //   int numScales = getNumScales(length);
47
48  //   delete [] coeffs;
49  //   coeffs = new T[(numScales+1)*length];
50  //   delete [] wavelet;
51  //   wavelet = new T[(numScales+1)*length];
52  //   for(int i=0;i<length;i++)  coeffs[i] = wavelet[i] = spectrum[i];
53
54  //   int spacing = 1;
55  //   for(int scale = 0; scale<numScales; scale++){
56
57  //     for(int pos = 0; pos<length; pos++){
58  //       coeffs[(scale+1)*length+pos] = 0;
59  //       for(int offset=-filterHW; offset<=filterHW; offset++){
60  //    int x = pos + spacing*offset;
61  //    if(x<0) x = -x;                    // boundary conditions are
62  //    if(x>=length) x = 2*(length-1) - x;    //    reflection.
63  // //         if(x<0) x = x+length;                // boundary conditions are
64  // //         if(x>=length) x = x-length;            //    continuous.
65
66  //    coeffs[(scale+1)*length+pos] += filter1D[offset+filterHW] *
67  //      coeffs[scale*length+x];
68  //       }
69  //       wavelet[(scale+1)*length+pos] = coeffs[scale*length+pos] -
70  //    coeffs[(scale+1)*length+pos];
71  //     }
72
73  //     spacing *= 2;
74
75  //   }
76
77  // }
78
79 
80  void atrousTransform(size_t &length, int &numScales, float *spectrum, double *coeffs, double *wavelet, duchamp::Param &par)
81  {
82    duchamp::Filter reconFilter = par.filter();
83    int filterHW = reconFilter.width()/2;
84
85    for(int i=0;i<length;i++)  coeffs[i] = wavelet[i] = spectrum[i];
86
87    int spacing = 1;
88    for(int scale = 0; scale<numScales; scale++){
89
90      for(int pos = 0; pos<length; pos++){
91        coeffs[(scale+1)*length+pos] = 0;
92        for(int offset=-filterHW; offset<=filterHW; offset++){
93          int x = pos + spacing*offset;
94          if(x<0) x = -x;                    // boundary conditions are
95          if(x>=length) x = 2*(length-1) - x;    //    reflection.
96          //    if(x<0) x = x+length;                // boundary conditions are
97          //    if(x>=length) x = x-length;            //    continuous.
98
99          //    coeffs[(scale+1)*length+pos] += filter1D[offset+filterHW] *
100          //      coeffs[scale*length+x];
101          coeffs[(scale+1)*length+pos] += reconFilter.coeff(offset+filterHW) *
102            coeffs[scale*length+x];
103        }
104        wavelet[(scale+1)*length+pos] = coeffs[scale*length+pos] -
105          coeffs[(scale+1)*length+pos];
106      }
107
108      spacing *= 2;
109
110    }
111
112  }
113
114  void atrousTransform(size_t &length, float *spectrum, float *coeffs, float *wavelet, duchamp::Param &par)
115  {
116    duchamp::Filter reconFilter = par.filter();
117    int filterHW = reconFilter.width()/2;
118    int numScales = reconFilter.getNumScales(length);
119
120    delete [] coeffs;
121    coeffs = new float[(numScales+1)*length];
122    delete [] wavelet;
123    wavelet = new float[(numScales+1)*length];
124    for(int i=0;i<length;i++) coeffs[i] = wavelet[i] = spectrum[i];
125
126    int spacing = 1;
127    for(int scale = 0; scale<numScales; scale++){
128
129      for(int pos = 0; pos<length; pos++){
130        coeffs[(scale+1)*length+pos] = 0;
131        for(int offset=-filterHW; offset<=filterHW; offset++){
132          int x = pos + spacing*offset;
133          if(x<0) x = -x;                    // boundary conditions are
134          if(x>=length) x = 2*(length-1) - x;    //    reflection.
135          //    if(x<0) x = x+length;                // boundary conditions are
136          //    if(x>=length) x = x-length;            //    continuous.
137
138          //    coeffs[(scale+1)*length+pos] += filter1D[offset+filterHW] *
139          //      coeffs[scale*length+x];
140          coeffs[(scale+1)*length+pos] += reconFilter.coeff(offset+filterHW) *
141            coeffs[scale*length+x];
142        }
143        wavelet[(scale+1)*length+pos] = coeffs[scale*length+pos] -
144          coeffs[(scale+1)*length+pos];
145      }
146
147      spacing *= 2;
148
149    }
150
151  }
152
153
154
155
156  /***********************************************************************/
157  /////  2-DIMENSIONAL TRANSFORM
158  /***********************************************************************/
159
160  void atrousTransform2D(size_t &xdim, size_t &ydim, int &numScales, float *input, double *coeffs, double *wavelet, duchamp::Param &par)
161  {
162    duchamp::Filter reconFilter = par.filter();
163    float blankPixValue = par.getBlankPixVal();
164    int filterHW = reconFilter.width()/2;
165
166    double *filter = new double[reconFilter.width()*reconFilter.width()];
167    for(int i=0;i<reconFilter.width();i++){
168      for(int j=0;j<reconFilter.width();j++){
169        filter[i*reconFilter.width()+j] = reconFilter.coeff(i) * reconFilter.coeff(j);
170      }
171    }
172
173    size_t size = xdim * ydim;
174    float *oldcoeffs = new float[size];
175
176    // locating the borders of the image -- ignoring BLANK pixels
177    //   int xLim1=0, yLim1=0, xLim2=xdim-1, yLim2=ydim-1;
178    //   for(int row=0;row<ydim;row++){
179    //     while((xLim1<xLim2)&&(input[row*xdim+xLim1]==blankPixValue)) xLim1++;
180    //     while((xLim2>xLim1)&&(input[row*xdim+xLim1]==blankPixValue)) xLim2--;
181    //   }
182    //   for(int col=0;col<xdim;col++){
183    //     while((yLim1<yLim2)&&(input[col+xdim*yLim1]==blankPixValue)) yLim1++;
184    //     while((yLim2>yLim1)&&(input[col+xdim*yLim1]==blankPixValue)) yLim2--;
185    //   }
186    //   std::cerr << "X Limits: "<<xLim1<<" "<<xLim2<<std::endl;
187    //   std::cerr << "Y Limits: "<<yLim1<<" "<<yLim2<<std::endl;
188 
189    int *xLim1 = new int[ydim];
190    int *yLim1 = new int[xdim];
191    int *xLim2 = new int[ydim];
192    int *yLim2 = new int[xdim];
193    for(int row=0;row<ydim;row++){
194      int ct1 = 0;
195      int ct2 = xdim - 1;
196      while((ct1<ct2)&&(input[row*xdim+ct1]==blankPixValue) ) ct1++;
197      while((ct2>ct1)&&(input[row*xdim+ct2]==blankPixValue) ) ct2--;
198      xLim1[row] = ct1;
199      xLim2[row] = ct2;
200      std::cerr<<row<<":"<<xLim1[row]<<","<<xLim2[row]<<" ";
201    }
202    std::cerr<<std::endl;
203
204    for(int col=0;col<xdim;col++){
205      int ct1=0;
206      int ct2=ydim-1;
207      while((ct1<ct2)&&(input[col+xdim*ct1]==blankPixValue) ) ct1++;
208      while((ct2>ct1)&&(input[col+xdim*ct2]==blankPixValue) ) ct2--;
209      yLim1[col] = ct1;
210      yLim2[col] = ct2;
211    }
212
213
214    bool *isGood = new bool[size];
215    for(int pos=0;pos<size;pos++) //isGood[pos] = (!flagBlank) || (input[pos]!=blankPixValue);
216      isGood[pos] = !par.isBlank(input[pos]);
217
218    for(int i=0;i<size;i++)  coeffs[i] = wavelet[i] = input[i];
219
220    int spacing = 1;
221    for(int scale = 0; scale<numScales; scale++){
222
223      for(int i=0;i<size;i++) oldcoeffs[i] = coeffs[i];
224
225      //std::cerr << numScales<<" "<<scale<<" "<<spacing<<" "<<reconFilter.width()*spacing<<std::endl;
226      for(int ypos = 0; ypos<ydim; ypos++){
227        for(int xpos = 0; xpos<xdim; xpos++){
228          // loops over each pixel in the image
229          int pos = ypos*xdim + xpos;
230          coeffs[pos] = 0;
231
232          //    if((par.getFlagBlankPix())&&(oldcoeffs[pos] == blankPixValue) )
233          if(par.isBlank(oldcoeffs[pos]) )
234            coeffs[pos] = oldcoeffs[pos];
235          else{
236            for(int yoffset=-filterHW; yoffset<=filterHW; yoffset++){
237              int y = ypos + spacing*yoffset;
238              int newy;
239              //          if(y<0) y = -y;                 // boundary conditions are
240              //          if(y>=ydim) y = 2*(ydim-1) - y; //    reflection.
241              //          while((y<yLim1)||(y>yLim2)){
242              //            if(y<yLim1) y = 2*yLim1 - y;      // boundary conditions are
243              //            if(y>yLim2) y = 2*yLim2 - y;      //    reflection.
244              //          }
245              // boundary conditions are reflection.
246         
247              for(int xoffset=-filterHW; xoffset<=filterHW; xoffset++){
248                int x = xpos + spacing*xoffset;
249                int newx;
250                //if(x<0) x = -x;                 // boundary conditions are
251                // if(x>=xdim) x = 2*(xdim-1) - x; //    reflection.
252                //while((x<xLim1)||(x>xLim2)){
253                //            if(x<xLim1) x = 2*xLim1 - x;      // boundary conditions are
254                //            if(x>xLim2) x = 2*xLim2 - x;      //    reflection.
255                //          }
256                // boundary conditions are reflection.
257                if(y<yLim1[xpos]) newy = 2*yLim1[xpos] - y;     
258                else if(y>yLim2[xpos]) newy = 2*yLim2[xpos] - y;     
259                else newy = y;
260                if(x<xLim1[ypos]) newx = 2*xLim1[ypos] - x;
261                else if(x>xLim2[ypos]) newx = 2*xLim2[ypos] - x;     
262                else newx=x;     
263             
264                x = newx;
265                y = newy;
266
267                int filterpos = (yoffset+filterHW)*reconFilter.width() + (xoffset+filterHW);
268                int oldpos = y*xdim + x;
269
270                if(// (x>=0)&&(x<xdim)&&(y>=0)&&(y<ydim)&&
271                   (isGood[pos]))
272                  coeffs[pos] += filter[filterpos] * oldcoeffs[oldpos];
273              }
274            }
275     
276          }
277          wavelet[(scale+1)*size+pos] = oldcoeffs[pos] - coeffs[pos];
278
279
280        }
281      }
282 
283      spacing *= 2;
284
285    }
286
287    delete [] filter;
288    delete [] oldcoeffs;
289
290  }
291
292  // void atrousTransform2D(long &xdim, long &ydim, int &numScales, float *input, double *coeffs, double *wavelet)
293  // {
294  //   Filter reconFilter = par.filter();
295  //   int filterHW = reconFilter.width()/2;
296
297  //   double *filter = new double[reconFilter.width()*reconFilter.width()];
298  //   for(int i=0;i<reconFilter.width();i++){
299  //     for(int j=0;j<reconFilter.width();j++){
300  //       filter[i*reconFilter.width()+j] = reconFilter.coeff(i) * reconFilter.coeff(j);
301  //     }
302  //   }
303
304  //   long size = xdim * ydim;
305  //   float *oldcoeffs = new float[size];
306
307  //   for(int i=0;i<size;i++)  coeffs[i] = wavelet[i] = input[i];
308
309 
310  //   int spacing = 1;
311  //   for(int scale = 0; scale<numScales; scale++){
312
313  //     for(int i=0;i<size;i++) oldcoeffs[i] = coeffs[i];
314
315  //     std::cerr << numScales<<" "<<scale<<" "<<spacing<<std::endl;
316  //     for(int ypos = 0; ypos<ydim; ypos++){
317  //       for(int xpos = 0; xpos<xdim; xpos++){
318  //    // loops over each pixel in the image
319  //    int pos = ypos*xdim + xpos;
320  //    coeffs[pos] = 0;
321
322  //    for(int yoffset=-filterHW; yoffset<=filterHW; yoffset++){
323  //      int y = ypos + spacing*yoffset;
324  //      if(y<0) y = -y;                 // boundary conditions are
325  //      if(y>=ydim) y = 2*(ydim-1) - y; //    reflection.
326         
327  //      for(int xoffset=-filterHW; xoffset<=filterHW; xoffset++){
328  //        int x = xpos + spacing*xoffset;
329  //        if(x<0) x = -x;                 // boundary conditions are
330  //        if(x>=xdim) x = 2*(xdim-1) - x; //    reflection.
331
332  //        int filterpos = (yoffset+filterHW)*reconFilter.width() + (xoffset+filterHW);
333  //        int oldpos = y*xdim + x;
334  //        coeffs[pos] += filter[filterpos] * oldcoeffs[oldpos];
335  //      }
336  //    }
337       
338  //    wavelet[(scale+1)*size+pos] = oldcoeffs[pos] - coeffs[pos];
339
340  //       }
341  //     }
342 
343  //     spacing *= 2;
344
345  //   }
346
347  //   delete [] filter;
348  //   delete [] oldcoeffs;
349
350  // }
351
352  /***********************************************************************/
353  /////  3-DIMENSIONAL TRANSFORM
354  /***********************************************************************/
355
356  void atrousTransform3D(size_t &xdim, size_t &ydim, size_t &zdim, int &numScales, float *&input, float *&coeffs, float *&wavelet, duchamp::Param &par)
357  {
358    duchamp::Filter reconFilter = par.filter();
359    float blankPixValue = par.getBlankPixVal();
360    int filterHW = reconFilter.width()/2;
361
362    size_t size = xdim * ydim * zdim;
363    float *oldcoeffs = new float[size];
364
365    std::cerr << "%";
366    int fsize = reconFilter.width()*reconFilter.width()*reconFilter.width();
367    std::cerr << "%";
368    double *filter = new double[fsize];
369    for(int i=0;i<reconFilter.width();i++){
370      for(int j=0;j<reconFilter.width();j++){
371        for(int k=0;k<reconFilter.width();k++){
372          filter[i +j*reconFilter.width() + k*reconFilter.width()*reconFilter.width()] =
373            reconFilter.coeff(i) * reconFilter.coeff(j) * reconFilter.coeff(k);
374        }
375      }
376    }
377    // locating the borders of the image -- ignoring BLANK pixels
378    // HAVE NOT DONE THIS FOR Z --> ASSUMING NO TRIMMING IN SPECTRAL DIRECTION
379    int xLim1 = 0, yLim1 = 0, xLim2 = xdim-1, yLim2 = ydim-1;
380    for(int col=0;col<xdim;col++){
381      while((yLim1<yLim2)&&(input[col+xdim*yLim1]==blankPixValue) ) yLim1++;
382      while((yLim2>yLim1)&&(input[col+xdim*yLim1]==blankPixValue) ) yLim2--;
383    }
384    for(int row=0;row<ydim;row++){
385      while((xLim1<xLim2)&&(input[row*xdim+xLim1]==blankPixValue) ) xLim1++;
386      while((xLim2>xLim1)&&(input[row*xdim+xLim1]==blankPixValue) ) xLim2--;
387    }
388
389    for(int i=0;i<size;i++)  coeffs[i] = wavelet[i] = input[i];
390
391    int spacing = 1;
392    std::cerr<<xdim<<"x"<<ydim<<"x"<<zdim<<"x"<<numScales;
393    for(int scale = 0; scale<numScales; scale++){
394      std::cerr << ".";
395      for(int i=0;i<size;i++) oldcoeffs[i] = coeffs[i];
396
397      for(int zpos = 0; zpos<zdim; zpos++){
398        for(int ypos = 0; ypos<ydim; ypos++){
399          for(int xpos = 0; xpos<xdim; xpos++){
400            // loops over each pixel in the image
401            int pos = zpos*xdim*ydim + ypos*xdim + xpos;
402            coeffs[pos] = 0;
403
404            if(par.isBlank(oldcoeffs[pos]) )
405              coeffs[pos] = oldcoeffs[pos];
406            else{
407              for(int zoffset=-filterHW; zoffset<=filterHW; zoffset++){
408                int z = zpos + spacing*zoffset;
409                if(z<0) z = -z;                 // boundary conditions are
410                if(z>=zdim) z = 2*(zdim-1) - z; //    reflection.
411         
412                for(int yoffset=-filterHW; yoffset<=filterHW; yoffset++){
413                  int y = ypos + spacing*yoffset;
414                  //if(y<0) y = -y;                 // boundary conditions are
415                  // if(y>=ydim) y = 2*(ydim-1) - y; //    reflection.
416                  if(y<yLim1) y = 2*yLim1 - y;      // boundary conditions are
417                  if(y>yLim2) y = 2*yLim2 - y;      //    reflection.
418         
419                  for(int xoffset=-filterHW; xoffset<=filterHW; xoffset++){
420                    int x = xpos + spacing*xoffset;
421                    //if(x<0) x = -x;                 // boundary conditions are
422                    //if(x>=xdim) x = 2*(xdim-1) - x; //    reflection.
423                    if(x<xLim1) x = 2*xLim1 - x;      // boundary conditions are
424                    if(x>xLim2) x = 2*xLim2 - x;      //    reflection.
425
426                    int filterpos = (zoffset+filterHW)*reconFilter.width()*reconFilter.width() +
427                      (yoffset+filterHW)*reconFilter.width() + (xoffset+filterHW);
428                    int oldpos = z*xdim*ydim + y*xdim + x;
429
430                    if(!par.isBlank(oldcoeffs[oldpos]))
431                      coeffs[pos] += filter[filterpos] * oldcoeffs[oldpos];
432
433                  }
434                }
435              }
436            }           
437 
438            wavelet[(scale+1)*size+pos] = oldcoeffs[pos] - coeffs[pos];
439
440          }
441        }
442      }
443 
444      spacing *= 2;
445
446    }
447    std::cerr << "|";
448
449    delete [] filter;
450    delete [] oldcoeffs;
451
452  }
453
454
455  // void atrousTransform3D(long &xdim, long &ydim, long &zdim, int &numScales, float *input, float *coeffs, float *wavelet)
456  // {
457  //   extern Filter reconFilter;
458  //   int filterHW = reconFilter.width()/2;
459
460  //   long size = xdim * ydim * zdim;
461  //   float *oldcoeffs = new float[size];
462
463  //   std::cerr << "%";
464  //   int fsize = reconFilter.width()*reconFilter.width()*reconFilter.width();
465  //   std::cerr << "%";
466  //   double *filter = new double[fsize];
467  //   for(int i=0;i<reconFilter.width();i++){
468  //     for(int j=0;j<reconFilter.width();j++){
469  //       for(int k=0;k<reconFilter.width();k++){
470  //       filter[i +j*reconFilter.width() + k*reconFilter.width()*reconFilter.width()] =
471  //    reconFilter.coeff(i) * reconFilter.coeff(j) * reconFilter.coeff(k);
472  //       }
473  //     }
474  //   }
475
476  //   for(int i=0;i<size;i++)  coeffs[i] = wavelet[i] = input[i];
477
478  //   int spacing = 1;
479  //   std::cerr<<xdim<<"x"<<ydim<<"x"<<zdim<<"x"<<numScales;
480  //   for(int scale = 0; scale<numScales; scale++){
481  //     std::cerr << ".";
482  //     for(int i=0;i<size;i++) oldcoeffs[i] = coeffs[i];
483
484  //     for(int zpos = 0; zpos<zdim; zpos++){
485  //       for(int ypos = 0; ypos<ydim; ypos++){
486  //    for(int xpos = 0; xpos<xdim; xpos++){
487  //      // loops over each pixel in the image
488  //      int pos = zpos*xdim*ydim + ypos*xdim + xpos;
489  //      coeffs[pos] = 0;
490
491  //      for(int zoffset=-filterHW; zoffset<=filterHW; zoffset++){
492  //        int z = zpos + spacing*zoffset;
493  //        if(z<0) z = -z;                 // boundary conditions are
494  //        if(z>=zdim) z = 2*(zdim-1) - z; //    reflection.
495         
496  //        for(int yoffset=-filterHW; yoffset<=filterHW; yoffset++){
497  //          int y = ypos + spacing*yoffset;
498  //          if(y<0) y = -y;                 // boundary conditions are
499  //          if(y>=ydim) y = 2*(ydim-1) - y; //    reflection.
500         
501  //          for(int xoffset=-filterHW; xoffset<=filterHW; xoffset++){
502  //            int x = xpos + spacing*xoffset;
503  //            if(x<0) x = -x;                 // boundary conditions are
504  //            if(x>=xdim) x = 2*(xdim-1) - x; //    reflection.
505
506  //            int filterpos = (zoffset+filterHW)*reconFilter.width()*reconFilter.width() +
507  //              (yoffset+filterHW)*reconFilter.width() + (xoffset+filterHW);
508  //            int oldpos = z*xdim*ydim + y*xdim + x;
509
510  //            coeffs[pos] += filter[filterpos] * oldcoeffs[oldpos];
511
512  //          }
513  //        }
514  //      }
515               
516  //      wavelet[(scale+1)*size+pos] = oldcoeffs[pos] - coeffs[pos];
517
518  //    }
519  //       }
520  //     }
521 
522  //     spacing *= 2;
523
524  //   }
525  //   std::cerr << "|";
526
527  //   delete [] filter;
528  //   delete [] oldcoeffs;
529
530  // }
531
532
533}
Note: See TracBrowser for help on using the repository browser.