////////////////////////////////////////////////////////////////////////

//-正變換，放在標頭檔案中

#pragma once

#include "cv.h"
#include "highgui.h"

using namespace cv;

Mat MWT(Mat img,int nlevs)
{
 int halfH = img.rows/2;
 int halfW = img.cols/2;

 Mat img64f;
 img.convertTo(img64f,CV_64F);

 Mat temp = Mat(img.rows,img.cols,CV_64F,Scalar::all(0));

 for(int n=nlevs;n>0;n--)
 {
  for(int i=0;i<halfH;i++)
  {
   for(int j=0;j<halfW;j++)
   {
    //結構元素B
    double b11 = img64f.at<double>(2*i,2*j);
    double b12 = img64f.at<double>(2*i,2*j+1);
    double b21 = img64f.at<double>(2*i+1,2*j);
    double b22 = img64f.at<double>(2*i+1,2*j+1);
    
    double M = min(min(min(b11,b12),b21),b22);
    double Yh = (b11-b21+b12-b22)/2;
    double Yv = (b11-b12+b21-b22)/2;    
    double Yd = (b11-b12-b21+b22)/2;
    
    temp.at<double>(i,j) = M;
    temp.at<double>(i,j+halfW) = Yh;
    temp.at<double>(i+halfH,j) = Yv;    
    temp.at<double>(i+halfH,j+halfW) = Yd;
   }
  }
  img64f = temp.clone();
  halfW = halfW/2;
  halfH = halfH/2;
 } 
 return temp;
}
////////////////////////////////////////////////////////////////////////////////////////////

//反變換，放在標頭檔案中

#ifndef _IMWT_H_
#define _IMWT_H_

#include "cv.h"
#include "highgui.h"

using namespace cv;

// 功能計演算法x的n次方
template <class T>
T power( T x, int m)
{
 int i;       //迴圈次數
 T p=1;   //存放結果
 for (i=1;i<=m;i++)
  p=p*x;
 return p;
}

Mat IMWT(Mat img,int nlevs)
{
 // 臨時矩陣,用於儲存矩陣資料
 Mat temp = Mat(img.rows,img.cols,CV_64F,Scalar::all(0));
 
 int p = power(2,nlevs);
 int iw = img.cols/p;
 int ih = img.rows/p;

 for(int n= nlevs;n>0;n--)
 {
  // 獲得尺度訊號的寬和高
  for(int i=0;i<ih;i++)
  {
   for(int j=0;j<iw;j++)
   {
    double m = img.at<double >(i,j);
    double h = img.at<double >(i,j+iw);
    double v = img.at<double >(i+ih,j);
    double d = img.at<double >(i+ih,j+iw);

    double  rm  = max(max(max(v+h,v+d),h+d),0.0);
    double  rh =  max(max(max(h-v,-v-d),h-d),0.0);
    double  rv =  max(max(max(v-h,v-d),-h-d),0.0);
    double  rd =  max(max(max(-v-h,d-v),d-h),0.0);

       temp.at<double>(2*i,2*j) = m+rm;
    temp.at<double>(2*i,2*j+1) = m+rh;
    temp.at<double>(2*i+1,2*j) = m+rv;
       temp.at<double>(2*i+1,2*j+1) = m+rd;
   }
  }
  iw = iw*2;
  ih = ih*2;
  for(int i=0;i<iw;i++)
  {
   for(int j=0;j<ih;j++)
   {
    img.at<double>(i,j) = temp.at<double>(i,j);
   }
  }
 }
 temp.convertTo(temp,CV_8U,1.0,0.0);
 return temp;
}
#endif

/////////////////////////////////////////

// -主程式

#include "cv.h"
#include "highgui.h"
#include "MWT.h"
#include "IMWT.h"

#include <iostream>
using namespace cv;
using namespace std;

int main()
{
 Mat img1 = imread("zoneplate.png",0);  //單通道
 int nlev;
 cout<<"Please input nlev:";
 cin>>nlev;

 Mat imgdec = Mat(img1.rows,img1.cols,CV_64F);
 imgdec = MWT(img1,nlev);

 Mat imgdec8u;
 imgdec.convertTo(imgdec8u,CV_8U);

 Mat imgrec = Mat(img1.rows,img1.cols,CV_8U,Scalar::all(0));
 imgrec = IMWT(imgdec,nlev);

 namedWindow("img",1);
 imshow("img",img1);

 namedWindow("imgdec",1);
 imshow("imgdec",imgdec8u);

 namedWindow("imgrec",1);
 imshow("imgrec",imgrec);

 subtract(img1,imgrec,imgrec);
 namedWindow("deta",1);
 imshow("deta",imgrec);

 Scalar m = mean(imgrec);

 cout<<"原圖-重建影象的均值："<<m.val[0]<<endl;
 cout<<"完美重建";

 waitKey();
 return 0;

}