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

Last change on this file since 1441 was 1393, checked in by MatthewWhiting, 10 years ago

A large changeset that moves all use of bool arrays (ie. bool *) to std::vector<bool>. This will be a bit safer, plus potentially cheaper in memory usage. Interfaces to all stats functions with masks have changed accordingly, as well as the Gaussian smoothing functions and some spectral utility functions.
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    std::vector<bool> isGood(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.