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

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

#pragma once

#include "cv.h"
#include "highgui.h"
#include <iostream>
#include <algorithm>

using namespace cv;
using namespace std;

// fuction：morphological transformation by the mearns of a paper
// named:Ishita De,Bhabatosh.A simple and efficient algorithm for
// multifocus image fusion uisng morphological wavelets
// Input:
//      img: an image
//      nlevs: levels of decomposition
// Output:
//      return a Mat(int)

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

 // 將img由uchar轉換為int
 Mat imgint;
 img.convertTo(imgint,CV_32S,1.0,0.0);

 Mat imgf = Mat(img.rows,img.cols,CV_32S,Scalar::all(0));

 for(int n=nlevs;n>0;n--)
 {  
  for(int i=0;i<nHeight;i=i+2)
  {
   for(int j=0;j<nWidth;j=j+2)
   {
    //結構元素B
    int b11 = imgint.at<int>(Point(i,j));
    int b12 = imgint.at<int>(Point(i,j+1));
    int b21 = imgint.at<int>(Point(i+1,j));
    int b22 = imgint.at<int>(Point(i+1,j+1));
    
    int m = max(max(max(b11,b12),b21),b22);
    int t1 = m-b12;
    int t2 = m-b21;
    int t3 = m-b22;
    
    imgf.at<int>(Point(i/2,j/2)) = m;
    
    if(t1>0)
     imgf.at<int>(Point(i/2,j/2+halfW)) = t1;
    else   
     t1 = b11-b12;
        imgf.at<int>(Point(i/2,j/2+halfW)) = t1;
     
    if(t2>0)
     imgf.at<int>(Point(i/2+halfH,j/2)) = t2;
    else
     t2 = b11-b21;
        imgf.at<int>(Point(i/2+halfH,j/2)) = t2;
     
    if(t3>0)
     imgf.at<int>(Point(i/2+halfH,j/2+halfW)) = t3;
    else
     t3 = b11-b22;
        imgf.at<int>(Point(i/2+halfH,j/2+halfW)) = t3;
   }
  }
  imgint = imgf.clone();
  nHeight = nHeight/2;
  nWidth  = nWidth/2;
  halfW = halfW/2;
  halfH = halfH/2;
 } 
 return imgf;
}
//------------------------------------------------------

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

#pragma once
#include "cv.h"
#include "highgui.h"
#include <iostream>
#include <algorithm>

using namespace cv;

// fuction：inverse morphological transformation by the mearns of a paper
// named:Ishita De,Bhabatosh.A simple and efficient algorithm for
// multifocus image fusion uisng morphological wavelets
// Input:
//      img: a Mat
//      nlevs: levels of decomposition
// Output:
//      return a Mat(int)

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

Mat IMWT(Mat img,int nlevs)
{
 // 臨時矩陣,用於儲存矩陣資料
 Mat temp = Mat(img.rows,img.cols,CV_32S,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++)
   {
    int m = img.at<int>(Point(i,j));
    int h = img.at<int>(Point(i,j+iw));
    int v = img.at<int>(Point(i+ih,j));
    int d = img.at<int>(Point(i+ih,j+iw));

    int rm = min(min(min(h,v),d),0);
    int rh = min(-h,0);
    int rv = min(-v,0);
    int rd = min(-d,0);

       temp.at<int>(Point(2*i,2*j)) = m+rm;
    temp.at<int>(Point(2*i,2*j+1)) = m+rh;
    temp.at<int>(Point(2*i+1,2*j)) = m+rv;
       temp.at<int>(Point(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<int>(Point(i,j)) = temp.at<int>(Point(i,j));
   }
  }
 }
 temp.convertTo(temp,CV_8U);
 return temp;
}

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

//主程式
#include "MWT.h"
#include "IMWT.h"

int main()
{
 Mat img1 = imread("medA-t.jpg",0);  //單通道

 cout<<"請輸入分解層數："<<endl;
 int nlevs;
 cin>>nlevs;
 Mat dec1(img1.rows,img1.cols,CV_32S);
 dec1 = MWT(img1,nlevs);
 Mat imgrec1(img1.rows,img1.cols,CV_8U);
 imgrec1 = IMWT(dec1,nlevs);

 dec1.convertTo(dec1,CV_8U);
 
 namedWindow("imgdec",1);
 imshow("imgdec",dec1);
 namedWindow("imgrec",1);
 imshow("imgrec",imgrec1);
 waitKey(0);

}