#include <iostream>
#include <iomanip>
#include <ATrous/atrous.hh>
#define DEBUG true

void getSigmaFactors1DNew(int &numScales)
{
//   extern Filter reconFilter;
//   std::cerr<<"-->"<<reconFilter.getName()<<"<--"<<std::endl;  
  std::cerr << "-->"<<Filter::getName()<<"<--"<<std::endl;
  std::cerr << "-->"<<getNumScales(125)<<"<--"<<std::endl;
  
//   int filterHW = reconFilter.width()/2;
  int filterwidth = Filter::width();
  int filterHW = Filter::width()/2;

  // starting values:
  long xdim = 1;
  long size = xdim;
  double *c = new double[size];
  double *w = new double[size];
  w[0] = 1.;

  double *filter = new double[filterwidth];
//   for(int i=0;i<filterwidth;i++) filter[i] = reconFilter.coeff(i);
  for(int i=0;i<filterwidth;i++) filter[i] = Filter::coeff(i);
  
  int spacing = 1;
//   const float FACTOR = 16.;
//   const float FACTOR = 1./reconFilter.coeff(0);
  const float FACTOR = 1.;
  const float FACFAC = FACTOR*FACTOR;

  if(DEBUG){
    std::cout << "divide by " << pow(FACTOR,0) << std::endl;
    for(int i=0;i<size;i++){
      std::cout<<i<<" "<<w[i]<<std::endl;
    }
    std::cout<<"===============\n";
  }

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

    // re-size and incorporate previous w into new c
    long oldsize = size;
    long oldxdim = xdim;
//     xdim = 2 * xdim + 3;
//     xdim = reconFilter.getMaxSize(scale);
    xdim += filterHW*long(pow(2,scale));
    size = xdim;
    double *newc = new double[size];
    for(int i=0;i<size;i++){
      int oldpos = i - xdim/2;
      if(// (oldpos<oldsize)&&(oldpos>=0)
	 abs(oldpos)<=oldsize/2) newc[i] = w[oldpos+oldsize/2] * FACTOR;
      else newc[i] = 0.;
      
      if(DEBUG) std::cout<<i<<" "<<oldpos<<" "<<size<<" "<<oldsize<<" "<<oldsize/2<<"   "
			 <<newc[i]// *pow(16,scale-1)
			 <<std::endl;
    }
    if(DEBUG) std::cout<<"---------------\nspacing = " << spacing << std::endl;
    delete [] c;
    c = newc;
    double *neww = new double[size];
    delete [] w;
    w = neww;
    for(int i=0;i<size;i++) w[i] = 0.;

    long double sum = 0.;

    for(int xpos = 0; xpos<xdim; xpos++){
      // loops over each pixel in the image
      w[xpos] = 0;

      for(int xoffset=-filterHW; xoffset<=filterHW; xoffset++){
	int x = xpos + spacing*xoffset;
	int filterpos = (xoffset+filterHW);
	if((x>=0)&&(x<xdim))  // if outside boundary --> set to 0.
	  w[xpos] += filter[filterpos] * c[x];
      } //-> end of xoffset loop
       
      sum += (c[xpos]-w[xpos])*(c[xpos]-w[xpos]) / pow(FACFAC,scale);

    } //-> end of xpos loop

    if(DEBUG){
      std::cout << "divide by " << pow(FACTOR,scale) << std::endl;
      for(int i=0;i<size;i++){
	std::cout<<i<<" "<<w[i]// *pow(16,scale)
		 <<std::endl;
      }
      std::cout<<"===============\n";
    }

    std::cerr << "Scale "<<std::setw(2)<<scale<<":";
    std::cerr.setf(std::ios::scientific);
    std::cerr.setf(std::ios::left);
    std::cerr <<"\t1D = "<<std::setw(15)<<std::setprecision(9)<<sqrt(sum)<<std::endl;
    std::cerr.unsetf(std::ios::left);
      

    spacing *= 2;

  } //-> end of scale loop 

}

void getSigmaFactors2DNew(int &numScales)
{
/*

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

  
  // starting values:
  long xdim = 1;
  long size = xdim*xdim*xdim;
  double *c = new double[size];
  double *w = new double[size];
  w[0] = 1.;

  int spacing = 1;
//   const float FACTOR = 1./(reconFilter.coeff(0)*reconFilter.coeff(0));
//   const float FACTOR = 16.*16.;
  const float FACTOR = 1.;
  const float FACFAC = FACTOR * FACTOR;

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

    // re-size and incorporate previous w into new c
    long oldsize = size;
    long oldxdim = xdim;
//     xdim = 2 * oldxdim + 3;
    xdim += filterHW*long(pow(2,scale));
    size = xdim * xdim;
    double *newc = new double[size];
    for(int i=0;i<size;i++){
      int x = i%xdim;
      int y = i/xdim;
      int centrepos=0;
      if((abs(x-xdim/2)<=oldxdim/2)&&(abs(y-xdim/2)<=oldxdim/2)){
	centrepos = x-xdim/2+oldxdim/2 + oldxdim*(y-xdim/2+oldxdim/2);
	newc[i] = w[centrepos] * FACTOR;
      }
      else newc[i] = 0.;
//       std::cout<<i<<" "<<x<<" "<<y<<" "<<centrepos<<" "<<xdim<<" "<<size
// 	       <<" "<<oldxdim<<" "<<oldsize<<"   "<<newc[i]*pow(256,scale-1)<<std::endl;
    }
//     std::cout<<"---------------\n";
    delete [] c;
    c = newc;
//     std::cerr<<".";
    double *neww = new double[size];
    delete [] w;
    w = neww;
    for(int i=0;i<size;i++) w[i] = 0.;
//     std::cerr<<".";

    long double sum = 0.;

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

	for(int yoffset=-filterHW; yoffset<=filterHW; yoffset++){
	  int y = ypos + spacing*yoffset;

	  for(int xoffset=-filterHW; xoffset<=filterHW; xoffset++){
	    int x = xpos + spacing*xoffset;

	    int filterpos = (yoffset+filterHW)*reconFilter.width() + (xoffset+filterHW);
	    int oldpos = y*xdim + x;
		    
	    if((y>=0)&&(y<xdim)&&(x>=0)&&(x<xdim))
	      w[pos] += filter[filterpos] * c[oldpos];

	  } //-> end of xoffset loop
	} //-> end of yoffset loop
 
	sum += (c[pos]-w[pos])*(c[pos]-w[pos]) / pow(FACFAC,scale);
	    
      } //-> end of xpos loop
    } //-> end of ypos loop

//     for(int i=0;i<size;i++){
//       std::cout<<i<<" "<<i%xdim<<" "<<i/xdim<<" "<<w[i]*pow(256,scale)<<std::endl;
//     }
//     std::cout<<"===============\n";

    std::cerr << "Scale "<<std::setw(2)<<scale<<":";
    std::cerr.setf(std::ios::scientific);
    std::cerr.setf(std::ios::left);
    std::cerr <<"\t2D = "<<std::setw(11)<<std::setprecision(9)<<sqrt(sum)<<std::endl;
    std::cerr.unsetf(std::ios::left);
      
    spacing *= 2;

  } //-> end of scale loop 

*/

}


void getSigmaFactors3DNew(int &numScales)
{
  /*  
  extern Filter reconFilter;
  int filterHW = reconFilter.width()/2;
  int fsize = reconFilter.width()*reconFilter.width()*reconFilter.width();
  double *filter = new double[fsize];
  for(int i=0;i<reconFilter.width();i++){
    for(int j=0;j<reconFilter.width();j++){
      for(int k=0;k<reconFilter.width();k++){
	filter[i +j*reconFilter.width() + k*reconFilter.width()*reconFilter.width()] = 
	  reconFilter.coeff(i) * reconFilter.coeff(j) * reconFilter.coeff(k);
      }
    }
  }

  // starting values:
  long xdim = 1;
  long size = xdim*xdim*xdim;
  double *c = new double[size];
  double *w = new double[size];
  int x,y,z,centrepos;
  w[0] = 1.;
  
  int spacing = 1;
//   const float FACTOR = 16.*16.*16.;
  const float FACTOR = 1.;
  const float FACFAC = FACTOR * FACTOR;

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

    // re-size and incorporate previous w into new c
    long oldsize = size;
    long oldxdim = xdim;
//     xdim = 2 * xdim + 3;
    xdim += filterHW*long(pow(2,scale));
    size = xdim * xdim * xdim;
    double *newc = new double[size];
    for(int i=0;i<size;i++){
      x = i%xdim;
      y = (i/xdim) % xdim;
      z = i / (xdim*xdim);
      centrepos = 0;
      if((abs(x-xdim/2)<=oldxdim/2)&&
	 (abs(y-xdim/2)<=oldxdim/2)&&
	 (abs(z-xdim/2)<=oldxdim/2)){
	centrepos = x-xdim/2+oldxdim/2 + oldxdim*(y-xdim/2+oldxdim/2)
	  + oldxdim*oldxdim*(z-xdim/2+oldxdim/2);
	newc[i] = w[centrepos] * FACTOR;
      }	
      else newc[i] = 0.;
//       std::cout<<i<<" "<<x<<" "<<y<<" "<<z<<" "<<oldpos<<" "<<size<<" "<<oldsize
// 	       <<"   "<<newc[i]*pow(16,3*(scale-1))<<std::endl;
    }
//     std::cout<<"---------------\n";
    delete [] c;
    c = newc;
    double *neww = new double[size];
    delete [] w;
    w = neww;
    for(int i=0;i<size;i++) w[i] = 0.;
    long double sum = 0.;

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

	  for(int zoffset=-filterHW; zoffset<=filterHW; zoffset++){
	    z = zpos + spacing*zoffset;
	    for(int yoffset=-filterHW; yoffset<=filterHW; yoffset++){
	      y = ypos + spacing*yoffset;
	      for(int xoffset=-filterHW; xoffset<=filterHW; xoffset++){
		x = xpos + spacing*xoffset;

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

		// only consider points inside boundary -- ie. =0 outside
		if((z>=0)&&(z<xdim)&&(y>=0)&&(y<xdim)&&(x>=0)&&(x<xdim))
		  w[pos] += filter[filterpos] * c[oldpos];

	      } //-> end of xoffset loop
	    } //-> end of yoffset loop
	  } //-> end of zoffset loop
 
	  sum += (c[pos]-w[pos])*(c[pos]-w[pos]) / pow(FACFAC,scale);
	    
	} //-> end of xpos loop
      } //-> end of ypos loop
    } //-> end of zpos loop

//     for(int i=0;i<size;i++){
//       std::cout<<i<<" "<<w[i]*pow(16,3*scale)<<std::endl;
//     }
//     std::cout<<"===============\n";

    std::cerr << "Scale "<<std::setw(2)<<scale<<":";
    std::cerr.setf(std::ios::scientific);
    std::cerr.setf(std::ios::left);
    std::cerr <<"\t3D = "<<std::setw(11)<<std::setprecision(9)<<sqrt(sum)<<std::endl;
    std::cerr.unsetf(std::ios::left);

    spacing *= 2;

  } //-> end of scale loop 

*/


}