#include <iostream>
#include <fstream>
#include <cstdlib>
#include <cstdio>
#include <iomanip>
#include <cmath>

#include "Daubechies.h"
#include "Coefficients.h"

char          Daubechies::TheP;
char          Daubechies::diese;
int           Daubechies::type;
int           Daubechies::dimX;
int           Daubechies::maxIntensity;
std::string   Daubechies::line_of_comments;

double      **Daubechies::v;
double      **Daubechies::workAAA;
double      **Daubechies::workBBB; 
double       *Daubechies::workCCC;
double       *Daubechies::h;
double       *Daubechies::g;
int           Daubechies::max_iter;

void Daubechies::Do_Discrete_Matrix_Transformation(char File_in[], 
           char File_out[], int Total_Iterations, int Dorder)
{

   ReadMatrix(File_in);

   AssignCoefficients(Dorder);

   char cumulMatrix[] = "cumulMatrice.dat"; 
   CumulativeEnergy(cumulMatrix);

   VerifyMaxIterations(Total_Iterations); 

   DWT2D(Dorder);

   int type;
   double threshold;

   std::cout << "Data Compression: press 1\n";
   std::cout << "Edge Detection: press 0\n";
   std::cin  >> type;
   
   if (type == 1) {
      std::cout << "Data Compression:  type in the threshold:";
      std::cin >> threshold;
   }
   AlterTransform(threshold, type);

   WriteMatrix(File_out);

   char cumulDaubechies[] = "cumulDaubechies.dat"; 
   CumulativeEnergy(cumulDaubechies);
  
   char image[] = "image.pgm"; Write_Matrix_as_PGM(image);

   for (int i=0; i < Daubechies::dimX+1; i++) { delete[] Daubechies::v[i]; } 
   delete[] Daubechies::v;

}

void Daubechies::Do_Inverse_Matrix_Transformation(char File_in[], 
           char File_out[], int Total_Iterations, int Dorder)
{

   ReadMatrix(File_in);

   AssignCoefficients(Dorder);

   VerifyMaxIterations(Total_Iterations); 

   IDWT2D(Dorder);
   WriteMatrix(File_out);

   char image[] = "InverseImage.pgm"; Write_Matrix_as_PGM(image);

   for (int i=0; i < Daubechies::dimX+1; i++) { delete[] Daubechies::v[i]; } 
   delete[] Daubechies::v;

}

// *****************************************************************************
// ******************************** DWT ****************************************
// *****************************************************************************
  
void Daubechies::DWT2D(int Dorder)
{

   int dim_matrix;

   for (int Jtrans=0; Jtrans < Daubechies::max_iter; Jtrans++) {

      dim_matrix = int(Daubechies::dimX/pow(2,Jtrans));

      AllocateWork(dim_matrix);

      GetCorner(Daubechies::workAAA,                      dim_matrix        );
      DWT2D1   (Daubechies::workAAA, Daubechies::workBBB, dim_matrix, Dorder);
      PutCorner(Daubechies::workBBB,                      dim_matrix        );

      DeallocateWork(dim_matrix);
   }

}

void Daubechies::DWT2D1(double **A, double **B, int dim_matrix, int Dorder)
{

   LeftDWT        (A, B, dim_matrix, Dorder);
   CopyMatrix     (B, A, dim_matrix        );
   TransposeMatrix(A,    dim_matrix        );

   LeftDWT        (A, B, dim_matrix, Dorder);
   TransposeMatrix(   B, dim_matrix        );

}

void Daubechies::LeftDWT(double **A, double **B, int dim_matrix, int Dorder)
{

   for (int j=1; j<= dim_matrix; j++) {
      DWT1D1(A[j], B[j], dim_matrix, Dorder);
   }

}

void Daubechies::DWT1D1(double *A, double *B, int N, int Dorder)
{
// Reference: Van Fleet, Discrete Wavelet Transformations, page 270-271

   int p, L;

   L = 2*Dorder - 1;
   p = (L - 1)/2; // number of wrapping rows.

// Compute lowpass and highpass values for the non-wrapping rows.
// To index highpass values, add N/2 to the subscript

   for (int k=1; k <= N/2 -p; k++) {
      B[      k] = 0.0;
      B[N/2 + k] = 0.0;
      for (int j=0; j <= L; j++) {
         B[      k] += Daubechies::h[L - j]*A[2*k - 1 + j];
         B[N/2 + k] += Daubechies::g[L - j]*A[2*k - 1 + j];
      }
   }

// compute lowpass and highpass values for the wrapping rows.
   for (int k=1; k <= p; k++) {
      B[N/2 - p + k] = 0.0;
      B[N   - p + k] = 0.0;
      // This loop is the end sum (7.94)
      for (int j=0; j <= L-2*k; j++) {
         B[N/2 - p + k] += Daubechies::h[L - j]*A[N - L + 2*k + j];
         B[N   - p + k] += Daubechies::g[L - j]*A[N - L + 2*k + j];
      }
      // This loop adds the beginning sum (7.93)
      for (int j=1; j <= 2*k; j++) {
         B[N/2 - p + k] += Daubechies::h[2*k - j]*A[j];
         B[N   - p + k] += Daubechies::g[2*k - j]*A[j];
      }
   }

}

// *****************************************************************************
// ******************************* IDWT ****************************************
// *****************************************************************************
  
void Daubechies::IDWT2D(int Dorder)
{
// Reference: Van Fleet, "Discrete Wavelet Transformations", 2007, 
// page 169 *********************************
// Careful, there are several errors in the algorithm given in the book by 
// Van Fleet
// Follow the argument in the book and some index become clear

   int dim_matrix;

   for (int Jtrans=1; Jtrans <= Daubechies::max_iter; Jtrans++) {

      dim_matrix = int(Daubechies::dimX/pow(2,Daubechies::max_iter - Jtrans));

      AllocateWork(dim_matrix);

      GetCorner(Daubechies::workAAA,                      dim_matrix        );
      IDWT2D1  (Daubechies::workAAA, Daubechies::workBBB, dim_matrix, Dorder);
      PutCorner(Daubechies::workBBB,                      dim_matrix        );

      DeallocateWork(dim_matrix);
   }
}

void Daubechies::IDWT2D1(double **A, double **B, int dim_matrix, int Dorder)
{
// Reference: Van Fleet, "Discrete Wavelet Transformations", 2007, 
//  page 196 *********************************

   LeftIDWT       (A, B, dim_matrix, Dorder);
   CopyMatrix     (B, A, dim_matrix        );
   TransposeMatrix(   A, dim_matrix        );

   LeftIDWT       (A, B, dim_matrix, Dorder);
   TransposeMatrix(   B, dim_matrix);

}

void Daubechies::LeftIDWT(double **A, double **B, int dim_matrix, int Dorder)
{
// Reference: Van Fleet, "Discrete Wavelet Transformations", 2007, 
// page 192 *********************************

   for (int j=1; j<= dim_matrix; j++) {
      IDWT1D1(A[j], B[j], dim_matrix, Dorder);
   }
 
}

void Daubechies::IDWT1D1(double *A, double *B, int N, int Dorder)
{
// algorithm page 277-278 Van Fleet

   double  lowPass[N/2 + 1];
   double highPass[N/2 + 1];
   double  InverseTransform[N + 1];

// Separate the matrix vector into low and high pass
   for (int i = 1; i <= N/2; i++) {
       lowPass[i] = A[i      ];
      highPass[i] = A[i + N/2];
   }

// nullify the vector, just in case.
   for (int i = 1; i <= N ; i++) { 
      InverseTransform[i] = 0.0;
   }

// Perform the inverse transform on the low pass section
   IDWTht( lowPass, InverseTransform, Daubechies::h, N, Dorder);

// copy result in work space B
   for (int i = 1; i <= N ; i++) { 
      B[i] = InverseTransform[i];
   }

// Perform the inverse transform on the high pass section
   IDWTht(highPass, InverseTransform, Daubechies::g, N, Dorder);

// add value to work space B
   for (int i = 1; i <= N ; i++) { 
      B[i] += InverseTransform[i];
   }

// B is return

}

void Daubechies::IDWTht(double *signal, double *InverseTransform, 
                        double *filter, int N, int Dorder)
{

   int L, p, r, n;

   L = 2*Dorder - 1;
   p = (L - 1)/2; // number of wrapping rows.
   r = p + 1;

   double even[r], odd[r];

   n = N/2;

// Initialize the vectors e and o.
   for (int k=0; k <= p; k++) {
      even[k+1] = filter[2*k  ];
       odd[k+1] = filter[2*k+1];
   }

// Compute the nonwrapping row values
   for (int k=1; k <= n-p; k++) {
      InverseTransform[L + 2*k - 2] = 0.0;
      InverseTransform[L + 2*k - 1] = 0.0;
//------------------------------------------------------------------
//    incorrect index in algorithm 7.2, it should be p+k           |
//    and equation (7.99) is wrong, it is (L+1)/2+k-1 or (L-1)/2+k |
//------------------------------------------------------------------
      for (int j=0; j<=p; j++) {
//       odd-even index start at 1, not 0
         InverseTransform[L + 2*k - 2] +=  odd[j + 1]*signal[j + k]; 
         InverseTransform[L + 2*k - 1] += even[j + 1]*signal[j + k]; 
      }
   }

// Compute the wrapping row values
   for (int k=1; k <= p; k++) {
      InverseTransform[2*k - 1] = 0.0;
      InverseTransform[2*k  ] = 0.0;

      // This loop is the beginning sum (7.101), page 275
      for (int j=1; j<= k ; j++) {
         InverseTransform[2*k - 1] +=  odd[r - k + j]*signal[j];
         InverseTransform[2*k    ] += even[r - k + j]*signal[j];
      }

      // This loop adds the end sm (7.102), page 275
//------------------------------------------------------------------
// End sum is incorrect in the textbook, it should be O[j]         |
//------------------------------------------------------------------
      for (int j=1; j<= r-k ; j++) {
         InverseTransform[2*k - 1] +=  odd[j]*signal[n - r + k + j];
         InverseTransform[2*k    ] += even[j]*signal[n - r + k + j];
      }
   }

}

// *****************************************************************************
// *************************** UTIL ********************************************
// *****************************************************************************
  
void Daubechies::AssignCoefficients(int Dorder)
{

   Daubechies::h = new double[2*Dorder];
   Daubechies::g = new double[2*Dorder];

   if (2*Dorder > Daubechies::dimX) {
      std::cout << "\n\nReduce order to a value less than or equal to " 
                << Daubechies::dimX/2 << "; Program is STOPPING\n";
      exit(8);
   }
      
   for (int i=0; i< 2*Dorder; i++) {
      switch (Dorder) {
         case  1: Daubechies::h[i] = Daub01[i]; break;
         case  2: Daubechies::h[i] = Daub02[i]; break;
         case  3: Daubechies::h[i] = Daub03[i]; break;
         case  4: Daubechies::h[i] = Daub04[i]; break;
         case  5: Daubechies::h[i] = Daub05[i]; break;
         case  6: Daubechies::h[i] = Daub06[i]; break;
         case  7: Daubechies::h[i] = Daub07[i]; break;
         case  8: Daubechies::h[i] = Daub08[i]; break;
         case  9: Daubechies::h[i] = Daub09[i]; break;
         case 10: Daubechies::h[i] = Daub10[i]; break;
         case 11: Daubechies::h[i] = Daub11[i]; break;
         case 12: Daubechies::h[i] = Daub12[i]; break;
         case 13: Daubechies::h[i] = Daub13[i]; break;
         case 14: Daubechies::h[i] = Daub14[i]; break;
         case 15: Daubechies::h[i] = Daub15[i]; break;
         case 16: Daubechies::h[i] = Daub16[i]; break;
         case 17: Daubechies::h[i] = Daub17[i]; break;
         case 18: Daubechies::h[i] = Daub18[i]; break;
         case 19: Daubechies::h[i] = Daub19[i]; break;
         case 20: Daubechies::h[i] = Daub20[i]; break;
         case 21: Daubechies::h[i] = Daub21[i]; break;
         case 22: Daubechies::h[i] = Daub22[i]; break;
         case 23: Daubechies::h[i] = Daub23[i]; break;
         case 24: Daubechies::h[i] = Daub24[i]; break;
         case 25: Daubechies::h[i] = Daub25[i]; break;
         case 26: Daubechies::h[i] = Daub26[i]; break;
         case 27: Daubechies::h[i] = Daub27[i]; break;
         case 28: Daubechies::h[i] = Daub28[i]; break;
         case 29: Daubechies::h[i] = Daub29[i]; break;
         case 30: Daubechies::h[i] = Daub30[i]; break;
         case 31: Daubechies::h[i] = Daub31[i]; break;
         case 32: Daubechies::h[i] = Daub32[i]; break;
         case 33: Daubechies::h[i] = Daub33[i]; break;
         case 34: Daubechies::h[i] = Daub34[i]; break;
         case 35: Daubechies::h[i] = Daub35[i]; break;
         case 36: Daubechies::h[i] = Daub36[i]; break;
         case 37: Daubechies::h[i] = Daub37[i]; break;
         case 38: Daubechies::h[i] = Daub38[i]; break;

         default: std::cout << "order must be between 1 and 38, inclusive\n"; 
                  exit(8); break;
         
         }
   }
   for (int k = 0; k <  2*Dorder; k++) {
      g[k] = pow(-1, k)*Daubechies::h[(2*Dorder-1) - k];
//    std::cout << "h[" << k << "]=" << Daubechies::h[k] <<"\n";
//    std::cout << "g[" << k << "]=" << Daubechies::g[k] <<"\n";
   }


}

void Daubechies::VerifyMaxIterations(int Total_Iterations)
{
   int ratio_of_matrix;
   int test_modulo;
   int i;

   Daubechies::max_iter = Total_Iterations;

   i = 1;
   while (i <= Total_Iterations) {
      ratio_of_matrix = int(pow(2,i));
      test_modulo = Daubechies::dimX % ratio_of_matrix;
      if (test_modulo != 0) {
         Daubechies::max_iter = i-1;
         break;
      }
      i++;
   }
   std::cout << "# iterations =" << Daubechies::max_iter << "\n";
}

void Daubechies::AllocateWork(int dim_matrix)
{
   Daubechies::workAAA = new double*[dim_matrix+1];
   Daubechies::workBBB = new double*[dim_matrix+1];
   for (int i=0; i < dim_matrix+1; i++) {
      Daubechies::workAAA[i] = new double[dim_matrix+1];
      Daubechies::workBBB[i] = new double[dim_matrix+1];
   }
}

void Daubechies::DeallocateWork(int dim_matrix)
{
   for (int i=0; i < dim_matrix+1; i++) {
      delete[] Daubechies::workAAA[i];
      delete[] Daubechies::workBBB[i];
   }
   delete[] Daubechies::workAAA;
   delete[] Daubechies::workBBB;
}

void Daubechies::GetCorner(double **workAAA, int dim_matrix)
{

   for (int i=1; i<= dim_matrix; i++) {
      for (int j=1; j<= dim_matrix; j++) {
         workAAA[i][j] = v[i][j];
      }
   }

}

void Daubechies::PutCorner(double **workBBB, int dim_matrix)
{

   for (int i=1; i<= dim_matrix; i++) {
      for (int j=1; j<= dim_matrix; j++) {
         v[i][j] = workBBB[i][j];
      }
   }
}

void Daubechies::TransposeMatrix(double **A, int dim_matrix)
{
   double tempo;

   for (int i=1; i<= dim_matrix; i++) {
      for (int j=i+1; j<= dim_matrix; j++) {
         tempo   = A[i][j];
         A[i][j] = A[j][i];
         A[j][i] = tempo;
      }
   }
}

void Daubechies::CopyMatrix(double **A, double **B, int dim_matrix)
{
   for (int i=1; i<= dim_matrix; i++) {
      for (int j=1; j<= dim_matrix; j++) {
      B[i][j] = A[i][j];
      }
   }
}

void Daubechies::CopyVector(double *A, double *B, int dim_vector)
{
   for (int i=1; i<= dim_vector; i++) {
      B[i] = A[i];
   }
}

void Daubechies::ReadMatrix(char File_in[])
{

// Notice that this version accepts only square matrices
   int dimY, dummy;
   std::ifstream file_in_A;

   file_in_A.open(File_in);
   if (file_in_A.fail() ) {
      std::cerr << "Unable to open input files\n";
      exit(8);
   }

   file_in_A >> Daubechies::TheP >> Daubechies::type;
   if (Daubechies::type == 999) {
      std::cout << "Real matrix\n";
      file_in_A >> Daubechies::dimX >> dimY;
   }
   if (Daubechies::type == 2) {
      std::cout << "PGM B&W ascii image\n";
      file_in_A >> Daubechies::diese;
      getline(file_in_A, Daubechies::line_of_comments);
      file_in_A >> Daubechies::dimX >> dimY;
   }
   if (Daubechies::dimX != dimY) {
      std::cout << "Something's wrong, PARSEC has square matrices; STOP\n";
      exit(8);
   }
   Daubechies::v = new double*[dimX+1];
   for (int i=0; i < dimX+1; i++) {
      Daubechies::v[i] = new double[dimX+1];
   }
// Notice that the index ij are reversed so that lines are consecutive for a 
// given column
   for (int i=1; i<= Daubechies::dimX; i++) {
      for (int j=1; j<= Daubechies::dimX; j++) {
         file_in_A  >> Daubechies::v[j][i];
      }
   }

   file_in_A.close();
}

void Daubechies::Write_Matrix_as_PGM(char File_out[])
{
   std::ofstream file_out_A;
   int white = 255;

   file_out_A.open(File_out);
   if (file_out_A.fail() ) {
      std::cerr << "Unable to open output files\n";
      exit(8);
   }

   file_out_A  << "P2\n"; // I don't put colors
   file_out_A << "#file.pgm\n";
   file_out_A << Daubechies::dimX << " " << Daubechies::dimX << "\n";
   file_out_A << 255 << "\n";

// find the min and the max
   double min =  1.0e+10;
   double max = -1.0e+10;
   double tempo;

   for (int i=1; i<= dimX; i++) {
     for (int j=1; j<= dimX; j++) {
       if (v[j][i] < min) min = v[j][i];
       if (v[j][i] > max) max = v[j][i];
     }
   }
   std::cout << "All Matrix: min=" << min << " max=" << max << "\n";

   for (int i=1; i<= dimX; i++) {
      for (int j=1; j<= dimX; j++) {
         tempo = (v[j][i]-min)/double(max-min);
         if (tempo < 0.0) { std::cout << "negative value\n"; exit(8); }
  
//       gamma transform
//       if (tempo > 1.0e-15) {
//          tempo = pow(tempo,0.1); 
//       } else {
//          tempo = 0;
//       }
         tempo *= white;
//       reverse video if wanted:
//       tempo = white - tempo;
         if (int(tempo) > 255) std::cout << tempo << i << j;
         file_out_A << int(tempo) << " ";
      }
      file_out_A << "\n";
   }

   file_out_A.close();
} 

void Daubechies::Write_Transform_as_PGM(char File_out[])
{
   std::ofstream file_out_A;
   int white = 255;

   file_out_A.open(File_out);
   if (file_out_A.fail() ) {
      std::cerr << "Unable to open output files\n";
      exit(8);
   }

   file_out_A  << "P2\n"; // I don't put colors
   file_out_A << "#file.pgm\n";
   file_out_A << Daubechies::dimX << " " << Daubechies::dimX << "\n";
   file_out_A << 255 << "\n";

// find the min and the max
   double min =  1.0e+10;
   double max = -1.0e+10;
   double tempo;

   for (int i=1; i<= dimX; i++) {
     for (int j=1; j<= dimX; j++) {
       if (v[j][i] < min) min = v[j][i];
       if (v[j][i] > max) max = v[j][i];
     }
   }
   std::cout << "All Matrix: min=" << min << " max=" << max << "\n";

   for (int i=1; i<= dimX; i++) {
      for (int j=1; j<= dimX; j++) {
         tempo = (v[j][i]-min)/double(max-min);
         if (tempo < 0.0) { std::cout << "negative value\n"; exit(8); }
  
//       gamma transform
//       if (tempo > 1.0e-15) {
//          tempo = pow(tempo,0.1); 
//       } else {
//          tempo = 0;
//       }
         tempo *= white;
//       reverse video if wanted:
//       tempo = white - tempo;
         if (int(tempo) > 255) std::cout << tempo << i << j;
         file_out_A << int(tempo) << " ";
      }
      file_out_A << "\n";
   }

   file_out_A.close();
} 

void Daubechies::WriteMatrix(char File_out[])
{
   std::ofstream file_out_A;

   file_out_A.open(File_out);
   if (file_out_A.fail() ) {
      std::cerr << "Unable to open output files\n";
      exit(8);
   }

   file_out_A  << Daubechies::TheP << 999 << "\n";
   file_out_A << Daubechies::dimX << " " << Daubechies::dimX << "\n";

   for (int i=1; i<= Daubechies::dimX; i++) {
      for (int j=1; j<= Daubechies::dimX; j++) {
         file_out_A  << std::setw(15) << std::fixed << std::setprecision(8) 
                     << Daubechies::v[j][i];
         if (dimX > 64) file_out_A << "\n";
    }
    if (dimX <= 64) file_out_A << "\n";
  }

  file_out_A.close();
} 

void Daubechies::Entropy()
{
// Calculated entropy
// Not used

   std::cout << "Entropy: not done yet, need practical definition for reals\n";
}

void Daubechies::CumulativeEnergy(char cumul[])
{
// Calculated cumulative energy page 81 Van Fleet

   std::ofstream file_out_cumul;

   file_out_cumul.open(cumul);
   if (file_out_cumul.fail() ) {
      std::cerr << "Unable to open output files\n";
      exit(8);
   }

   double norm2;
   double cumul_E;

   int total_mat;
   total_mat = (Daubechies::dimX)*(Daubechies::dimX);
   Daubechies::workCCC = new double[total_mat];

   for (int i=1; i<= dimX; i++) {
      for (int j=1; j<= dimX; j++) {
         if (v[i][j] <  0.0) workCCC[(i-1) + (dimX)*(j-1)] = -v[i][j];
         if (v[i][j] >= 0.0) workCCC[(i-1) + (dimX)*(j-1)] =  v[i][j];
      }
   }

   quickSort(workCCC,total_mat);

// At this point the elements are sorted in increasing order, 
// which means that element workCCC[dim] is the largest

// calculate square of the norm
   norm2 = 0.0;
   for (int i = total_mat - 1; i >= 0; i--) {
      norm2 += workCCC[i]*workCCC[i];
   }

   cumul_E = 0.0;
   for (int i = total_mat - 1; i >= 0; i--) {
      cumul_E += workCCC[i]*workCCC[i]/norm2;
//    The conditions here reduces the size of the vector
//    Those numbers are arbitrarily set
      if ( (total_mat - i) <= 671) 
        file_out_cumul << total_mat -i << " " << cumul_E << "\n";
      if ( (total_mat-i)%1001 == 1) 
        file_out_cumul << total_mat -i << " " << cumul_E << "\n";
   }

   delete[] Daubechies::workCCC;
  file_out_cumul.close();
}

void Daubechies::PrintMatrixDEBUG(std::string message, double **matrix, 
          int dim_matrix)
{
// Not used, only for debugging purposes

   for (int i=1; i<=dim_matrix; i++) {
      for (int j=1; j<=dim_matrix; j++) {
         std::cout << matrix[i][j] << " ";
      }
      std::cout << "MATRIX " << message << "\n";
   }

}

void Daubechies::PrintVectorDEBUG(std::string message, double *vector, 
         int dim_vector)
{
// Not used, only for debugging purposes

   for (int i=1; i<=dim_vector; i++) {
      std::cout << vector[i] << " ";
   }
   std::cout << "VECTOR " << message << "\n";

}

void Daubechies::quickSort(double *numbers, int array_size)
{
// Reference: http://linux.wku.edu/~lamonml/algor/sort/quick.html
  q_sort(numbers, 0, array_size - 1);
}


void Daubechies::q_sort(double *numbers, int left, int right)
{
// Reference: http://linux.wku.edu/~lamonml/algor/sort/quick.html
  double pivot; 
  int l_hold, r_hold, pivot_index;

  l_hold = left;
  r_hold = right;
  pivot = numbers[left];
  while (left < right)
  {
    while ((numbers[right] >= pivot) && (left < right))
      right--;
    if (left != right)
    {
      numbers[left] = numbers[right];
      left++;
    }
    while ((numbers[left] <= pivot) && (left < right))
      left++;
    if (left != right)
    {
      numbers[right] = numbers[left];
      right--;
    }
  }
  numbers[left] = pivot;
  pivot_index   = left;
  left          = l_hold;
  right         = r_hold;
  if (left < pivot_index)
    q_sort(numbers, left, pivot_index-1);
  if (right > pivot_index)
    q_sort(numbers, pivot_index+1, right);
}

void Daubechies::AlterTransform(double threshold, int type)
{
// example: set to zero the blur matrix

   int dim_blur;

   dim_blur = int(Daubechies::dimX/pow(2,Daubechies::max_iter));

   if (type == 1) {
      std::cout << "WARNING: Compression: zeroing some of the complement\n"
                << " of the blur sub-matrix of dim=" 
                << dim_blur << "\n";
      for (int i=1; i<= dimX; i++) {
         for (int j=1; j<= dimX; j++) {
            if (!(i <= dim_blur && j <= dim_blur)) {  
               if (abs(v[j][i]) < threshold) Daubechies::v[j][i] = 0.0;
            }
         }
      } 
   }
   if (type == 0) {
      std::cout << "WARNING: Edge-detection: Zeroing the blur sub-matrix of dim=" 
                << dim_blur << "\n";
      for (int i=1; i<= dimX; i++) {
         for (int j=1; j<= dimX; j++) {
            if ( (i <= dim_blur && j <= dim_blur)) {  
               Daubechies::v[j][i] = 0.0;
            }
         }
      } 
   }
      
}