source: branches/fitshead-branch/ATrous/atrous_2d_reconstruct.cc @ 95

#include <iostream>
#include <iomanip>
#include <math.h>
#include <ATrous/atrous.hh>
#include <Utils/utils.hh>

using std::endl;
using std::setw;

void atrous2DReconstruct(long &xdim, long &ydim, float *&input,float *&output, Param &par)
{
  /**
   *  atrous2DReconstruct(xdim, ydim, input, output, par)
   *
   *  A routine that uses the a trous wavelet method to reconstruct a 
   *   2-dimensional image.
   *  The Param object "par" contains all necessary info about the filter and 
   *   reconstruction parameters, although a Filter object has to be declared
   *   elsewhere previously.
   *  The input array is in "input", of dimensions "xdim"x"ydim", and the reconstructed
   *   array is in "output".
   */

  extern Filter reconFilter;
  bool flagBlank=par.getFlagBlankPix();
  float blankPixValue = par.getBlankPixVal();
  long size = xdim * ydim;
  long mindim = xdim;
  if (ydim<mindim) mindim = ydim;
  int numScales = reconFilter.getNumScales(mindim);
  double *sigmaFactors = new double[numScales+1];
  for(int i=0;i<=numScales;i++){
    if(i<=reconFilter.maxFactor(2)) sigmaFactors[i] = reconFilter.sigmaFactor(2,i);
    else sigmaFactors[i] = sigmaFactors[i-1] / 2.;
  }

  float mean,sigma,originalSigma,originalMean,oldsigma,newsigma;
  bool *isGood = new bool[size];
  for(int pos=0;pos<size;pos++) 
    isGood[pos] = !par.isBlank(input[pos]);
 
  float *array = new float[size];
  int goodSize=0;
  for(int i=0;i<size;i++) if(isGood[i]) array[goodSize++] = input[i];
  findMedianStats(array,goodSize,originalMean,originalSigma);
  delete [] array;
  
  float *coeffs    = new float[size];
  float *wavelet   = new float[size];

  for(int pos=0;pos<size;pos++) output[pos]=0.;

  int filterHW = reconFilter.width()/2;
  double *filter = new double[reconFilter.width()*reconFilter.width()];
  for(int i=0;i<reconFilter.width();i++){
    for(int j=0;j<reconFilter.width();j++){
      filter[i*reconFilter.width()+j] = reconFilter.coeff(i) * reconFilter.coeff(j);
    }
  }

  int *xLim1 = new int[ydim];
  int *yLim1 = new int[xdim];
  int *xLim2 = new int[ydim];
  int *yLim2 = new int[xdim];
  float avGapX = 0, avGapY = 0;
  for(int row=0;row<ydim;row++){
    int ct1 = 0;
    int ct2 = xdim - 1;
    while((ct1<ct2)&&(input[row*xdim+ct1]==blankPixValue) ) ct1++;
    while((ct2>ct1)&&(input[row*xdim+ct2]==blankPixValue) ) ct2--;
    xLim1[row] = ct1;
    xLim2[row] = ct2;
    avGapX += ct2 - ct1;
  }
  avGapX /= float(ydim);

  for(int col=0;col<xdim;col++){
    int ct1=0;
    int ct2=ydim-1;
    while((ct1<ct2)&&(input[col+xdim*ct1]==blankPixValue) ) ct1++;
    while((ct2>ct1)&&(input[col+xdim*ct2]==blankPixValue) ) ct2--;
    yLim1[col] = ct1;
    yLim2[col] = ct2;
    avGapY += ct2 - ct1;
  }
  avGapY /= float(ydim);
 
  mindim = int(avGapX);
  if(avGapY < avGapX) mindim = int(avGapY);
  numScales = reconFilter.getNumScales(mindim);
  

  /***********************************************************************/
  /***********************************************************************/

  float threshold;
  int iteration=0;
  newsigma = 1.e9;
  for(int i=0;i<size;i++) output[i] = 0;
  do{
    if(par.isVerbose()) std::cout << "Iteration #"<<++iteration<<":             ";
    // first, get the value of oldsigma and set it to the previous newsigma value
    oldsigma = newsigma;
    // we are transforming the residual array
    for(int i=0;i<size;i++)  coeffs[i] = input[i] - output[i];  

    int spacing = 1;
    for(int scale = 1; scale<numScales; scale++){

      if(par.isVerbose()) {
        std::cout << "\b\b\b\b\b\b\b\b\b\b\b\bScale ";
        std::cout << setw(2)<<scale<<" /"<<setw(2)<<numScales<<std::flush;
      }

      for(int ypos = 0; ypos<ydim; ypos++){
        for(int xpos = 0; xpos<xdim; xpos++){
          // loops over each pixel in the image
          int pos = ypos*xdim + xpos;
          
          wavelet[pos] = coeffs[pos];

          if(!isGood[pos]) wavelet[pos] = 0.;
          else{

            for(int yoffset=-filterHW; yoffset<=filterHW; yoffset++){
              int y = ypos + spacing*yoffset;
              int newy;
              //          if(y<0) y = -y;                 // boundary conditions are 
              //          if(y>=ydim) y = 2*(ydim-1) - y; //    reflection.
              //          while((y<yLim1)||(y>yLim2)){
              //            if(y<yLim1) y = 2*yLim1 - y;      // boundary conditions are 
              //            if(y>yLim2) y = 2*yLim2 - y;      //    reflection.
              //          }
              // boundary conditions are reflection.
                while((y<yLim1[xpos])||(y>yLim2[xpos])){
                  if(y<yLim1[xpos]) y = 2*yLim1[xpos] - y;      
                  else if(y>yLim2[xpos]) y = 2*yLim2[xpos] - y;      
                }
          
              for(int xoffset=-filterHW; xoffset<=filterHW; xoffset++){
                int x = xpos + spacing*xoffset;
                int newx;
                //if(x<0) x = -x;                 // boundary conditions are 
                // if(x>=xdim) x = 2*(xdim-1) - x; //    reflection.
                //while((x<xLim1)||(x>xLim2)){
                //            if(x<xLim1) x = 2*xLim1 - x;      // boundary conditions are 
                //            if(x>xLim2) x = 2*xLim2 - x;      //    reflection.
                //          }
                // boundary conditions are reflection.
                while((x<xLim1[ypos])||(x>xLim2[ypos])){
                  if(x<xLim1[ypos]) x = 2*xLim1[ypos] - x;      
                  else if(x>xLim2[ypos]) x = 2*xLim2[ypos] - x;      
                }

//              if(y<yLim1[xpos]) newy = 2*yLim1[xpos] - y;      
//              else if(y>yLim2[xpos]) newy = 2*yLim2[xpos] - y;      
//              else newy = y;
//              if(x<xLim1[ypos]) newx = 2*xLim1[ypos] - x;
//              else if(x>xLim2[ypos]) newx = 2*xLim2[ypos] - x;      
//              else newx=x;      
              
//              x = newx;
//              y = newy;

                int filterpos = (yoffset+filterHW)*reconFilter.width() + (xoffset+filterHW);
                int oldpos = y*xdim + x;

                if(isGood[pos]) 
                  wavelet[pos] -= filter[filterpos] * coeffs[oldpos];

              } //-> end of xoffset loop
            } //-> end of yoffset loop
          } //-> end of else{ ( from if(!isGood[pos])  )
        
        } //-> end of xpos loop
      } //-> end of ypos loop

      // Need to do this after we've done *all* the convolving
      for(int pos=0;pos<size;pos++) coeffs[pos] = coeffs[pos] - wavelet[pos];

      // Have found wavelet coeffs for this scale -- now threshold    
      if(scale>=par.getMinScale()){
        array = new float[size];
        goodSize=0;
        for(int pos=0;pos<size;pos++) if(isGood[pos]) array[goodSize++] = wavelet[pos];
        findMedianStats(array,goodSize,mean,sigma);
        delete [] array;

        threshold = mean + par.getAtrousCut() * originalSigma * sigmaFactors[scale];
        for(int pos=0;pos<size;pos++){
          if(!isGood[pos]) output[pos] = blankPixValue; // preserve the Blank pixel values in the output.
          else if( fabs(wavelet[pos]) > threshold ) output[pos] += wavelet[pos];
        }
      }
      spacing *= 2;

    } // END OF LOOP OVER SCALES

    for(int pos=0;pos<size;pos++) if(isGood[pos]) output[pos] += coeffs[pos];

    array = new float[size];
    goodSize=0;
    for(int i=0;i<size;i++) if(isGood[i]) array[goodSize++] = input[i] - output[i];
    findMedianStats(array,goodSize,mean,newsigma);
    delete [] array;
    
    if(par.isVerbose()) std::cout << "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b";

  } while( (iteration==1) || 
           (fabsf(oldsigma-newsigma)/newsigma > reconTolerance) );

  if(par.isVerbose()) std::cout << "Completed "<<iteration<<" iterations. ";

  delete [] coeffs;
  delete [] wavelet;
  delete [] isGood;
  delete [] filter;
  delete [] sigmaFactors;
}