showdibview.cpp

用vc实现的小波分析的mallat算法用以实现图像融合

// ShowDIBView.cpp : implementation of the CShowDIBView class
//

#include "stdafx.h"
#include "ShowDIB.h"

#include "ShowDIBDoc.h"
#include "ShowDIBView.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

/////////////////////////////////////////////////////////////////////////////
// CShowDIBView

IMPLEMENT_DYNCREATE(CShowDIBView, CView)

BEGIN_MESSAGE_MAP(CShowDIBView, CView)
	//{{AFX_MSG_MAP(CShowDIBView)
	ON_COMMAND(ID_TRANS_WAVE, OnTransWave)
	//}}AFX_MSG_MAP
	// Standard printing commands
	ON_COMMAND(ID_FILE_PRINT, CView::OnFilePrint)
	ON_COMMAND(ID_FILE_PRINT_DIRECT, CView::OnFilePrint)
	ON_COMMAND(ID_FILE_PRINT_PREVIEW, CView::OnFilePrintPreview)
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CShowDIBView construction/destruction

CShowDIBView::CShowDIBView()
{
	// TODO: add construction code here

}

CShowDIBView::~CShowDIBView()
{
}

BOOL CShowDIBView::PreCreateWindow(CREATESTRUCT& cs)
{
	// TODO: Modify the Window class or styles here by modifying
	//  the CREATESTRUCT cs

	return CView::PreCreateWindow(cs);
}

/////////////////////////////////////////////////////////////////////////////
// CShowDIBView drawing

void CShowDIBView::OnDraw(CDC* pDC)
{
	CShowDIBDoc* pDoc = GetDocument();
	ASSERT_VALID(pDoc);
	// TODO: add draw code for native data here
	if(!pDoc->m_pDib->IsEmpty())
	{
		CSize sizeDisp;
		sizeDisp = pDoc->m_pDib->GetDimensions();
		pDoc->m_pDib->Draw(pDC, CPoint(0, 0), sizeDisp);
	}
}

/////////////////////////////////////////////////////////////////////////////
// CShowDIBView printing

BOOL CShowDIBView::OnPreparePrinting(CPrintInfo* pInfo)
{
	// default preparation
	return DoPreparePrinting(pInfo);
}

void CShowDIBView::OnBeginPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
	// TODO: add extra initialization before printing
}

void CShowDIBView::OnEndPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
	// TODO: add cleanup after printing
}

/////////////////////////////////////////////////////////////////////////////
// CShowDIBView diagnostics

#ifdef _DEBUG
void CShowDIBView::AssertValid() const
{
	CView::AssertValid();
}

void CShowDIBView::Dump(CDumpContext& dc) const
{
	CView::Dump(dc);
}

CShowDIBDoc* CShowDIBView::GetDocument() // non-debug version is inline
{
	ASSERT(m_pDocument->IsKindOf(RUNTIME_CLASS(CShowDIBDoc)));
	return (CShowDIBDoc*)m_pDocument;
}
#endif //_DEBUG

/////////////////////////////////////////////////////////////////////////////
// CShowDIBView message handlers

void CShowDIBView::OnTransWave() 
{
	// TODO: Add your command handler code here
	CShowDIBDoc * pDoc = GetDocument();

	CDib * pDib = pDoc->m_pDib;	
	LPBITMAPINFOHEADER lpBMIH = pDib->m_lpBMIH;
	if (lpBMIH->biBitCount != 8)
	{
		MessageBox("目前只支持256色位图!", "系统提示" ,
			MB_ICONINFORMATION | MB_OK);
		return;
	}
	
	CSize sizeImage	= pDib->GetDimensions();
	int nWidth = sizeImage.cx		;
	int nHeight = sizeImage.cy		;
	LPBYTE lpImage = pDib->m_lpImage;
	//CSize sizeImageSave	= pDib->GetDibSaveDim();
	//int nSaveWidth = sizeImageSave.cx;
	
	//W97-1 wavelet
	double Lo[20] = {-0.125,0.25,0.75,0.25,-0.125,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
	double Hi[20]={0.25,-0.5,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};

	int i, j;
	double* LLLH;
	double* HHHL;
	double* tempImage;
	int tempHeight = nHeight;
	int tempWidth = nWidth;
	double* resImage = new double[nHeight * nWidth];

	for(int n = 0; n < 1 ; n++)//此处的n为分解的层数
	{
		if(n == 0)
		{
			tempImage = new double[nHeight * nWidth];

			for(i = 0; i < nHeight; i++)
				for(j = 0; j < nWidth; j++)
					tempImage[i * nWidth + j] = lpImage[i * nWidth + j];
		}

		LLLH = new double[tempHeight * tempWidth / 2];
    	HHHL = new double[tempHeight * tempWidth / 2];
		LLLH = RealMallat(Lo, 5, 9, Hi, 4, 7, tempImage, tempHeight, tempWidth, TRUE);
		HHHL = RealMallat(Hi, 4, 7, Lo, 5, 9, tempImage, tempHeight, tempWidth, FALSE);

		for(i = (nHeight - tempHeight); i < nHeight; i++)
			for(j = 0; j < tempWidth / 2; j++)
			{
				resImage[i * nWidth + j] = LLLH[(i - (nHeight - tempHeight)) * (tempWidth / 2) + j];
				resImage[i * nWidth + j + tempWidth / 2] = HHHL[(i - (nHeight - tempHeight)) * (tempWidth / 2) + j];
			}

		delete[] tempImage;
		tempHeight /= 2;
		tempWidth /= 2;
		tempImage = new double[tempHeight * tempWidth];

		for(i = 0; i < tempHeight; i++)
			for(j = 0; j < tempWidth; j++)
				tempImage[i * tempWidth + j] = LLLH[(i + tempHeight) * tempWidth + j];
	}
	
	//Gray strain
	double max, min;
	max = -2000;
	min = 2000;
	for(i = 0; i < nHeight; i++)
		for(j = 0; j < nWidth; j++)
		{
			if(resImage[i * nWidth + j] < min)
				min = resImage[i * nWidth + j];
			if(resImage[i * nWidth + j] > max)
				max = resImage[i * nWidth + j];
		}

	for(i = 0; i < nHeight; i++)
		for(j = 0; j < nWidth; j++)
			resImage[i * nWidth + j] = (resImage[i * nWidth + j] - min) / (max-min) * 255;

	for(i = 0; i < nHeight; i++)
		for(j = 0; j < nWidth; j++)
			pDoc->m_pDib->m_lpImage[i * nWidth + j] = FloatToByte(resImage[i * nWidth + j]);
	pDoc->UpdateAllViews(NULL);
	pDoc->SetModifiedFlag(TRUE);
	//DisplayImage(nHeight, nWidth, resImage);
	delete[] resImage;
}

double* CShowDIBView::RealMallat(double h[], int h0, int ln_h, double g[], int g0, int ln_g, double *Image, int hh, int ww, bool Low)
{
	double *tmp=new double[ww*hh];
	double *resultImage=new double[ww*hh/2];//transformed image
	int i,j,k;
	double var;
	for(i=0;i<hh;i++)
	{
		for(j=0;j<h0-1;j++)
		{
			var=0;
			for(k=0;k<ln_h;k++)
			{	
				if((j+k-(h0-1)<0))
					var=var+h[k]*Image[i*ww+(abs(j+k-(h0-1))-(ln_h-1)%2)];
				else
					var=var+h[k]*Image[i*ww+(j+k-(h0-1))];
			}
			tmp[i*ww+j]=var;
		}

		for(j=h0-1;j<ww-(ln_h-h0);j++)
		{
			var=0;
			for(k=0;k<ln_h;k++)
				var=var+h[k]*Image[i*ww+j+k-(h0-1)];
			tmp[i*ww+j]=var;
		}
		for(j=ww-(ln_h-h0);j<ww;j++)
		{
			var=0;
			for(k=0;k<ln_h;k++)
			{
				/*if((j+k-h0+1)<ww)
					var=var+h[k]*Image[(j+k-h0+1)*ww+j];
				else
					var=var+h[k]*Image[(2*ww-(j+k-h0+1)-2+(ln_h-1)%2)*ww+j];*/
				if((j + k - h0 + 1) >= ww)
					var = var + h[k] * Image[i * ww + ww - (j + k - h0 + 1 - ww) - (ln_h - 1) % 2 - 2];
				else
					var = var + h[k] * Image[i * ww + (j + k - h0 + 1)];
			}
			tmp[i*ww+j]=var;
		}
	}

	for(i=0;i<hh;i++)
	{
		for(j=0;j<ww/2;j++)
			tmp[i*ww+j]=tmp[i*ww+2*j];
	}
//行低通计算完成

//列高通。LH HL
	
	for(j=0;j<ww/2;j++)
	{
		for(i=0;i<g0-1;i++)
		{
			var=0;
			for(k=0;k<ln_g;k++)
			{
				if((i+k-(g0-1)<0))
					var=var+g[k]*tmp[(abs(i+k-(g0-1))-(ln_g-1)%2)*ww+j];
				else
					var=var+g[k]*tmp[(i+k-(g0-1))*ww+j];
			}
			tmp[i*ww+j+ww/2]=var;
		}

		for(i=g0-1;i<hh-(ln_g-g0);i++)
		{
			var=0;
			for(k=0;k<ln_g;k++)
				var=var+g[k]*tmp[(i+k-(g0-1))*ww+j];
			tmp[i*ww+j+ww/2]=var;
		}

		for(i=hh-(ln_g-g0);i<hh;i++)
		{
			var=0;
			for(k=0;k<ln_g;k++)
			{
				/*if((i+k-g0+1)<hh)
					var=var+g[k]*tmp[(i+k-g0+1)*ww+j];
				else
					var=var+g[k]*tmp[(2*hh-(i+k-g0+1)-2+(ln_g-1)%2)*ww+j];*/
				if((i + k - g0 + 1) >= hh)
					var = var + g[k] * tmp[(hh - ((i + k - g0 + 1) - hh) + (ln_g - 1) % 2 - 2) * ww + j];
				else
					var = var + g[k] * tmp[(i + k - g0 + 1) * ww + j];
			}
			tmp[i*ww+j+ww/2]=var;
		}
	}//return tmp;

	for(i=0;i<hh/2;i++)
		for(j=0;j<ww/2;j++)
		{
			if(Low)
				resultImage[i*(ww/2)+j]=tmp[2*i*ww+j+ww/2];
			else
				resultImage[(i+hh/2)*(ww/2)+j]=tmp[2*i*ww+j+ww/2];
		}
//计算LL HH
	for(j=0;j<ww/2;j++)
	{
		for(i=0;i<h0-1;i++)
		{
			var=0;
			for(k=0;k<ln_h;k++)
			{
				if((i+k-(h0-1)<0))
					var=var+h[k]*tmp[(abs(i+k-(h0-1))-(ln_h-1)%2)*ww+j];
				else
					var=var+h[k]*tmp[(i+k-(h0-1))*ww+j];
			}
			tmp[i*ww+j+ww/2]=var;
		}

		for(i=h0-1;i<hh-(ln_h-h0);i++)
		{
			var=0;
			for(k=0;k<ln_h;k++)
				var=var+h[k]*tmp[(i+k-(h0-1))*ww+j];
			tmp[i*ww+j+ww/2]=var;
		}

		for(i=hh-(ln_h-h0);i<hh;i++)
		{
			var=0;
			for(k=0;k<ln_h;k++)
			{
				if((i+k-h0+1)<hh)
					var=var+h[k]*tmp[(i+k-h0+1)*ww+j];
				else
					var=var+h[k]*tmp[(2*hh-(i+k-h0+1)-2+(ln_h-1)%2)*ww+j];//*/
				//if((i + k - h0 + 1))
			}
			tmp[i*ww+j+ww/2]=var;
		}
	}//return tmp;
   
	for(i=hh/2;i<hh;i++)
		for(j=0;j<ww/2;j++)
		{
			if(Low)
				resultImage[i*(ww/2)+j]=tmp[2*(i - hh/2)*ww+j+ww/2];
			else
				resultImage[(i-hh/2)*(ww/2)+j]=tmp[2*(i - hh/2)*ww+j+ww/2];
		}
//	DisplayImage(m_orgWidth/2,m_orgHeight,0,resultImage);
	return resultImage;
}