source: trunk/src/ATrous/atrous_2d_reconstruct.cc @ 848

Last change on this file since 848 was 846, checked in by MatthewWhiting, 13 years ago

Lots of small changes to the reconstruction code, getting the type definitions appropriate (signed vs unsigned ints, with size_t parameters in there as well.)

File size: 9.3 KB
Line 
1// -----------------------------------------------------------------------
2// atrous_2d_reconstruct.cc: 2-dimensional wavelet reconstruction.
3// -----------------------------------------------------------------------
4// Copyright (C) 2006, Matthew Whiting, ATNF
5//
6// This program is free software; you can redistribute it and/or modify it
7// under the terms of the GNU General Public License as published by the
8// Free Software Foundation; either version 2 of the License, or (at your
9// option) any later version.
10//
11// Duchamp is distributed in the hope that it will be useful, but WITHOUT
12// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13// FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
14// for more details.
15//
16// You should have received a copy of the GNU General Public License
17// along with Duchamp; if not, write to the Free Software Foundation,
18// Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
19//
20// Correspondence concerning Duchamp may be directed to:
21//    Internet email: Matthew.Whiting [at] atnf.csiro.au
22//    Postal address: Dr. Matthew Whiting
23//                    Australia Telescope National Facility, CSIRO
24//                    PO Box 76
25//                    Epping NSW 1710
26//                    AUSTRALIA
27// -----------------------------------------------------------------------
28#include <iostream>
29#include <iomanip>
30#include <math.h>
31#include <duchamp/duchamp.hh>
32#include <duchamp/param.hh>
33#include <duchamp/ATrous/atrous.hh>
34#include <duchamp/ATrous/filter.hh>
35#include <duchamp/Utils/utils.hh>
36#include <duchamp/Utils/feedback.hh>
37#include <duchamp/Utils/Statistics.hh>
38using Statistics::madfmToSigma;
39
40namespace duchamp
41{
42
43  void atrous2DReconstruct(unsigned long &xdim, unsigned long &ydim, float *&input, float *&output, Param &par)
44  {
45    ///  A routine that uses the a trous wavelet method to reconstruct a
46    ///   2-dimensional image.
47    ///
48    ///  If there are no non-BLANK pixels (and we are testing for
49    ///  BLANKs), the reconstruction cannot be done, so we return the
50    ///  input array as the output array and give a warning message.
51    ///
52    ///  \param xdim The length of the x-axis of the image.
53    ///  \param ydim The length of the y-axis of the image.
54    ///  \param input The input spectrum.
55    ///  \param output The returned reconstructed spectrum. This array
56    ///  needs to be declared beforehand.
57    ///  \param par The Param set:contains all necessary info about the
58    ///  filter and reconstruction parameters.
59
60    size_t size = xdim * ydim;
61    unsigned long mindim = xdim;
62    if (ydim<mindim) mindim = ydim;
63    size_t numScales = par.filter().getNumScales(mindim);
64    double *sigmaFactors = new double[numScales+1];
65    for(size_t i=0;i<=numScales;i++){
66      if(i<=par.filter().maxFactor(2))
67        sigmaFactors[i] = par.filter().sigmaFactor(2,i);
68      else sigmaFactors[i] = sigmaFactors[i-1] / 2.;
69    }
70
71    float mean,sigma,originalSigma,originalMean,oldsigma,newsigma;
72    size_t goodSize=0;
73    bool *isGood = new bool[size];
74    for(size_t pos=0;pos<size;pos++){
75      isGood[pos] = !par.isBlank(input[pos]);
76      if(isGood[pos]) goodSize++;
77    }
78
79    if(goodSize == 0){
80      // There are no good pixels -- everything is BLANK for some reason.
81      // Return the input array as the output, and give a warning message.
82
83      for(size_t pos=0;pos<size; pos++) output[pos] = input[pos];
84
85      duchampWarning("2D Reconstruction","\
86There are no good pixels to be reconstructed -- all are BLANK.\n\
87Returning input array.\n");
88    }
89    else{
90      // Otherwise, all is good, and we continue.
91
92      findMedianStats(input,goodSize,isGood,originalMean,originalSigma);
93      originalSigma = madfmToSigma(originalSigma);
94 
95      float *coeffs    = new float[size];
96      float *wavelet   = new float[size];
97      float *residual  = new float[size];
98
99      for(size_t pos=0;pos<size;pos++) output[pos]=0.;
100
101      unsigned int filterwidth = par.filter().width();
102      int filterHW = filterwidth/2;
103      double *filter = new double[filterwidth*filterwidth];
104      for(size_t i=0;i<filterwidth;i++){
105        for(size_t j=0;j<filterwidth;j++){
106          filter[i*filterwidth+j] = par.filter().coeff(i) * par.filter().coeff(j);
107        }
108      }
109
110      long *xLim1 = new long[ydim];
111      for(size_t i=0;i<ydim;i++) xLim1[i] = 0;
112      long *yLim1 = new long[xdim];
113      for(size_t i=0;i<xdim;i++) yLim1[i] = 0;
114      long *xLim2 = new long[ydim];
115      for(size_t i=0;i<ydim;i++) xLim2[i] = xdim-1;
116      long *yLim2 = new long[xdim];
117      for(size_t i=0;i<xdim;i++) yLim2[i] = ydim-1;
118
119      if(par.getFlagBlankPix()){
120        float avGapX = 0, avGapY = 0;
121        for(size_t row=0;row<ydim;row++){
122          size_t ct1 = 0;
123          size_t ct2 = xdim - 1;
124          while((ct1<ct2)&&(par.isBlank(input[row*xdim+ct1]))) ct1++;
125          while((ct2>ct1)&&(par.isBlank(input[row*xdim+ct2]))) ct2--;
126          xLim1[row] = ct1;
127          xLim2[row] = ct2;
128          avGapX += ct2 - ct1;
129        }
130        avGapX /= float(ydim);
131   
132        for(size_t col=0;col<xdim;col++){
133          size_t ct1=0;
134          size_t ct2=ydim-1;
135          while((ct1<ct2)&&(par.isBlank(input[col+xdim*ct1]))) ct1++;
136          while((ct2>ct1)&&(par.isBlank(input[col+xdim*ct2]))) ct2--;
137          yLim1[col] = ct1;
138          yLim2[col] = ct2;
139          avGapY += ct2 - ct1;
140        }
141        avGapY /= float(xdim);
142   
143        if(avGapX < mindim) mindim = int(avGapX);
144        if(avGapY < mindim) mindim = int(avGapY);
145        numScales = par.filter().getNumScales(mindim);
146      }
147
148      float threshold;
149      int iteration=0;
150      newsigma = 1.e9;
151      for(size_t i=0;i<size;i++) output[i] = 0;
152      do{
153        if(par.isVerbose()) {
154          std::cout << "Iteration #"<<std::setw(2)<<++iteration<<":";
155          printBackSpace(13);
156        }
157
158        // first, get the value of oldsigma and set it to the previous
159        //   newsigma value
160        oldsigma = newsigma;
161        // we are transforming the residual array
162        for(size_t i=0;i<size;i++)  coeffs[i] = input[i] - output[i]; 
163
164        int spacing = 1;
165        for(unsigned int scale = 1; scale<numScales; scale++){
166
167          if(par.isVerbose()){
168            std::cout << "Scale ";
169            std::cout << std::setw(2)<<scale<<" / "<<std::setw(2)<<numScales;
170            printBackSpace(13);
171            std::cout <<std::flush;
172          }
173
174          for(unsigned long ypos = 0; ypos<ydim; ypos++){
175            for(unsigned long xpos = 0; xpos<xdim; xpos++){
176              // loops over each pixel in the image
177              size_t pos = ypos*xdim + xpos;
178         
179              wavelet[pos] = coeffs[pos];
180
181              if(!isGood[pos]) wavelet[pos] = 0.;
182              else{
183
184                size_t filterpos = -1;
185                for(int yoffset=-filterHW; yoffset<=filterHW; yoffset++){
186                  long y = long(ypos) + spacing*yoffset;
187                  // Boundary conditions -- assume reflection at boundaries.
188                  // Use limits as calculated above
189                  //          if(yLim1[xpos]!=yLim2[xpos]){
190                  //            // if these are equal we will get into an infinite loop here
191                  //            while((y<yLim1[xpos])||(y>yLim2[xpos])){
192                  //              if(y<yLim1[xpos]) y = 2*yLim1[xpos] - y;     
193                  //              else if(y>yLim2[xpos]) y = 2*yLim2[xpos] - y;     
194                  //            }
195                  //          }
196                  size_t oldrow = y * xdim;
197         
198                  for(int xoffset=-filterHW; xoffset<=filterHW; xoffset++){
199                    long x = long(xpos) + spacing*xoffset;
200                    // Boundary conditions -- assume reflection at boundaries.
201                    // Use limits as calculated above
202                    //          if(xLim1[ypos]!=xLim2[ypos]){
203                    //            // if these are equal we will get into an infinite loop here
204                    //            while((x<xLim1[ypos])||(x>xLim2[ypos])){
205                    //              if(x<xLim1[ypos]) x = 2*xLim1[ypos] - x;     
206                    //              else if(x>xLim2[ypos]) x = 2*xLim2[ypos] - x;     
207                    //            }
208                    //          }
209
210                    size_t oldpos = oldrow + x;
211
212                    float oldCoeff;
213                    if((y>=yLim1[xpos])&&(y<=yLim2[xpos])&&
214                       (x>=xLim1[ypos])&&(x<=xLim2[ypos])  )
215                      oldCoeff = coeffs[oldpos];
216                    else oldCoeff = 0.;
217
218                    filterpos++;
219
220                    if(isGood[pos]) wavelet[pos] -= filter[filterpos] * oldCoeff;
221                    //            wavelet[pos] -= filter[filterpos] * coeffs[oldpos];
222
223                  } //-> end of xoffset loop
224                } //-> end of yoffset loop
225              } //-> end of else{ ( from if(!isGood[pos])  )
226       
227            } //-> end of xpos loop
228          } //-> end of ypos loop
229
230          // Need to do this after we've done *all* the convolving
231          for(size_t pos=0;pos<size;pos++) coeffs[pos] = coeffs[pos] - wavelet[pos];
232
233          // Have found wavelet coeffs for this scale -- now threshold   
234          if(scale>=par.getMinScale()){
235            findMedianStats(wavelet,goodSize,isGood,mean,sigma);
236
237            threshold = mean +
238              par.getAtrousCut() * originalSigma * sigmaFactors[scale];
239            for(size_t pos=0;pos<size;pos++){
240              if(!isGood[pos]) output[pos] = input[pos];
241              // preserve the Blank pixel values in the output.
242              else if( fabs(wavelet[pos]) > threshold )
243                output[pos] += wavelet[pos];
244            }
245          }
246          spacing *= 2;
247
248        } // END OF LOOP OVER SCALES
249
250        for(size_t pos=0;pos<size;pos++) if(isGood[pos]) output[pos] += coeffs[pos];
251
252        for(size_t i=0;i<size;i++) residual[i] = input[i] - output[i];
253        findMedianStats(residual,goodSize,isGood,mean,newsigma);
254        newsigma = madfmToSigma(newsigma);
255   
256        if(par.isVerbose()) printBackSpace(15);
257
258      } while( (iteration==1) ||
259               (fabs(oldsigma-newsigma)/newsigma > reconTolerance) );
260
261      if(par.isVerbose()) std::cout << "Completed "<<iteration<<" iterations. ";
262
263      delete [] xLim1;
264      delete [] xLim2;
265      delete [] yLim1;
266      delete [] yLim2;
267      delete [] filter;
268      delete [] coeffs;
269      delete [] wavelet;
270      delete [] residual;
271
272    }
273
274    delete [] isGood;
275    delete [] sigmaFactors;
276  }
277   
278}
Note: See TracBrowser for help on using the repository browser.