pixelfusiondoc.cpp

本程序采用Visual C++6.0对全色图像和多光谱遥感图像进行了基于相关系数的图像融合

//+++++++++++++++++++++++++++++++++++++++++++++
void CPixelFusionDoc::OnCorrelationFusion() 
{
	// 存储文件对话框
	CFileDialog dlg(FALSE, "Bmp", NULL, OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT, 
		"BMP文件(*.bmp)|*.bmp|所有文件(*.*)|*.*||",NULL);
	if (dlg.DoModal() == IDCANCEL)
		return ;

	// 获取存储文件路径
	lpszResultPath = dlg.GetPathName();

	// TODO: Add your command handler code here
	BITMAPFILEHEADER bmfHeaderSpot,bmfHeaderMulspec;
	BITMAPINFOHEADER bmiHeaderSpot,bmiHeaderMulspec;

	////////////////////////////////////////////////
	//获取高空间分辨率影像的信息	
	CFile fileSpot;
	if(!fileSpot.Open(lpszSpotPath,CFile::modeRead))
	{
		AfxMessageBox("当前文件不能打开!",MB_OK);
		return ;
	}
	//读文件头
	if(fileSpot.Read(&bmfHeaderSpot,sizeof(BITMAPFILEHEADER)) != sizeof(BITMAPFILEHEADER))
		return ;
	//判断读入的文件是否是BMP文件
	if(bmfHeaderSpot.bfType != DIB_HEADER_MARKER)
		return ;
	//读信息头
	if(fileSpot.Read(&bmiHeaderSpot,sizeof(BITMAPINFOHEADER)) != sizeof(BITMAPINFOHEADER))
	{
		return ;
	}	
	WORD wBitCountSpot = bmiHeaderSpot.biBitCount; 
	LONG  lSpotWidth = bmiHeaderSpot.biWidth;
	LONG  lSpotHeight = bmiHeaderSpot.biHeight;
	DWORD dwSpotSizeOff = bmfHeaderSpot.bfOffBits; //文件头到象素的偏移量
	//计算高分辨率图象的行宽度
	LONG lLineBytesSpot = WIDTHBYTES(bmiHeaderSpot.biWidth*bmiHeaderSpot.biBitCount*bmiHeaderSpot.biPlanes);

	//判断图像是否满足SPOT全色影像的条件 
	if(wBitCountSpot != 8)
	{
		AfxMessageBox("SPOT全色图像必须为8-bit,请重新输入图像!",MB_OK);
		return ;
	}
	fileSpot.Close();

	////////////////////////////////////////////////
	//获取高光谱分辨率影像的信息
	CFile fileMulspec;
	if(!fileMulspec.Open(lpszMulspecPath,CFile::modeRead))
	{
		AfxMessageBox("当前文件不能打开!",MB_OK);
		return ;
	}	
	//读文件头
	if(fileMulspec.Read(&bmfHeaderMulspec,sizeof(BITMAPFILEHEADER)) != sizeof(BITMAPFILEHEADER))
		return ;
	//判断读入的文件是否是BMP文件
	if(bmfHeaderMulspec.bfType != DIB_HEADER_MARKER)
		return ;
	//读信息头
	if(fileMulspec.Read(&bmiHeaderMulspec,sizeof(BITMAPINFOHEADER)) != sizeof(BITMAPINFOHEADER))
	{
		return ;
	}	
	WORD wBitCountMulspec = bmiHeaderMulspec.biBitCount; 
	LONG  lMulspecWidth = bmiHeaderMulspec.biWidth;
	LONG  lMulspecHeight = bmiHeaderMulspec.biHeight;
	DWORD dwMulspecSizeOff = bmfHeaderMulspec.bfOffBits; //文件头到象素的偏移量
	//计算高分辨率图象的行宽度
	LONG lLineBytesMulspec = WIDTHBYTES(bmiHeaderMulspec.biWidth*bmiHeaderMulspec.biBitCount*bmiHeaderMulspec.biPlanes);
	int m = bmiHeaderMulspec.biWidth*bmiHeaderMulspec.biBitCount*bmiHeaderMulspec.biPlanes;
	//判断图像是否满足SPOT全色影像的条件 
	if(wBitCountMulspec != 24)
	{
		AfxMessageBox("多光谱图像必须为24-bit,请重新输入图像!",MB_OK);
		return ;
	}
	fileMulspec.Close();

	LONG i,j;
	////////////////////////////////////////////////
	//判断两幅图像尺寸是否一致。
	if((lSpotWidth != lMulspecWidth) || (lSpotHeight != lMulspecHeight))
	{
		AfxMessageBox("图像大小不匹配，请重新输入!");
		return ;
	}
	/////////////////////////////////////////////////
	//建立文件映射
	//高空间分辨率图象
	HANDLE hFileSpot = ::CreateFile(lpszSpotPath,GENERIC_READ,FILE_SHARE_READ|FILE_SHARE_WRITE,NULL,
		                            OPEN_EXISTING,FILE_ATTRIBUTE_NORMAL,NULL);
	HANDLE hMapSpot = ::CreateFileMapping(hFileSpot,NULL,PAGE_READONLY,0,0,NULL);

	if(hMapSpot != NULL && GetLastError() == ERROR_ALREADY_EXISTS)
	{
		CloseHandle(hMapSpot);
		hMapSpot = NULL;
		return ;
	}
	LPSTR lpSpotHead = (LPSTR)::MapViewOfFile(hMapSpot,FILE_MAP_READ,0,0,0);
	LPSTR lpSpotDIBBits = lpSpotHead + dwSpotSizeOff;
	unsigned char* lpSpotBits;       //指向SPOT源图象象素的指针
		
	//高光谱分辨率图象
	HANDLE hFileMulspec = ::CreateFile(lpszMulspecPath,GENERIC_READ,FILE_SHARE_READ|FILE_SHARE_WRITE,NULL,
		OPEN_EXISTING,FILE_ATTRIBUTE_NORMAL,NULL);
	HANDLE hMapMulspec = ::CreateFileMapping(hFileMulspec,NULL,PAGE_READONLY,0,0,NULL);
	
	if(hMapMulspec != NULL && GetLastError() == ERROR_ALREADY_EXISTS)
	{
		CloseHandle(hMapMulspec);
		hMapMulspec = NULL;
		return ;
	}
	LPSTR lpMulspecHead = (LPSTR)::MapViewOfFile(hMapMulspec,FILE_MAP_READ,0,0,0);
	LPSTR lpMulspecDIBBits = lpMulspecHead + dwMulspecSizeOff;
	unsigned char* lpMulspecBits;       //指向多光谱源图象象素的指针
	
	////////////////////////////////////////////////////////////////////////////
	
	//计算SPOT影像的平均灰度
	double SpotMeanGray=0;
	double temp = 0;
	for(i=0;i<lSpotHeight;i++)
	{
		for(j=0;j<lSpotWidth;j++)
		{
			lpSpotBits = (unsigned char*)lpSpotDIBBits 
				             + lLineBytesSpot*(lSpotHeight-1-i) + j;
			temp += *lpSpotBits;
		}
	}
	SpotMeanGray = temp/((double)lSpotHeight*lSpotWidth);
	temp = 0;

    double dMulspecMeanGrayR=0,dMulspecMeanGrayG=0,dMulspecMeanGrayB=0;
	//R
	for(i=0;i<lMulspecHeight;i++)
	{
		for(j=0;j<lMulspecWidth;j++)
		{
			lpMulspecBits = (unsigned char*)lpMulspecDIBBits 
				+ lLineBytesMulspec*(lSpotHeight-1-i) + j*3;
			temp += *lpMulspecBits;
		}
	}
	
	dMulspecMeanGrayR = temp/((double)lMulspecHeight*lMulspecWidth);
	temp = 0;
	//G
	for(i=0;i<lMulspecHeight;i++)
	{
		for(j=0;j<lMulspecWidth;j++)
		{
			lpMulspecBits = (unsigned char*)lpMulspecDIBBits 
				+ lLineBytesMulspec*(lSpotHeight-1-i) + j*3+1;
			temp += *lpMulspecBits;
		}
	}
	dMulspecMeanGrayG = temp/((double)lMulspecHeight*lMulspecWidth);
	temp = 0;
	//B
	for(i=0;i<lMulspecHeight;i++)
	{
		for(j=0;j<lMulspecWidth;j++)
		{
			lpMulspecBits = (unsigned char*)lpMulspecDIBBits 
				+ lLineBytesMulspec*(lSpotHeight-1-i) + j*3+2;
			temp += *lpMulspecBits;
		}
	}
	dMulspecMeanGrayB = temp/((double)lMulspecHeight*lMulspecWidth);
	temp = 0;

////////////////////////////////////////////////////////////////////////////
	//计算SPOT全色影像和多光谱影像每个波段的相关系数。
	double dRelation[3];
	for(i=0;i<3;i++)
	{
		dRelation[i]=0;
	}
	double temp1=0,temp2=0,temp3=0;	
	//计算SPOT全色影像和多光谱影像各个波段的相关系数,
	//dRelation[1]，dRelation[2],dRelation[3]分别代表与红波段，绿波段，兰波段的相关系数。
	for(int k=0;k<3;k++)
	{
		for(i=0;i<lSpotHeight;i++)
		{
			for(j=0;j<lSpotWidth;j++)
			{
				
				lpSpotBits = (unsigned char*)lpSpotDIBBits 
						+ lLineBytesSpot*(lSpotHeight-1-i) + j;
				lpMulspecBits = (unsigned char*)lpMulspecDIBBits 
						+ lLineBytesMulspec*(lSpotHeight-1-i) + j*3+k;
				temp1 += (*lpSpotBits - SpotMeanGray)*(*lpMulspecBits - dMulspecMeanGrayR);
				temp2 += (*lpSpotBits - SpotMeanGray)*(*lpSpotBits - SpotMeanGray);
				temp3 += (*lpMulspecBits - dMulspecMeanGrayR)*(*lpMulspecBits - dMulspecMeanGrayR);		
			}
		}
		dRelation[k] = temp1/sqrt(temp2*temp3);
		temp1=temp2=temp3=0;
	}
	/////////////////////////////////////////////////////////
	//创建融合后影像的文件头
	//根据高光谱分辨率图象的头信息写融合生成图象的文件
	CFile fileResult(lpszResultPath,CFile::modeCreate|CFile::modeWrite);
	fileResult.Write(&bmfHeaderMulspec,sizeof(BITMAPFILEHEADER));
	fileResult.Write(&bmiHeaderMulspec,sizeof(BITMAPINFO));
	BYTE* data;
	if((data = new BYTE[lLineBytesMulspec]) == NULL)
	{
		AfxMessageBox("Can't alloc enough memory!",MB_OK);
		return ;
	}
	for(i = 0;i<lLineBytesMulspec;i++)
	{
		data[i] = 0;
	}
	fileResult.Seek(dwMulspecSizeOff,CFile::begin);
	for(i=0;i<lMulspecHeight;i++)
	{
		fileResult.Write(data,lLineBytesMulspec);
	}

	fileResult.Close();
	//融合结果影像的文件映射
	HANDLE hFileResult = ::CreateFile(lpszResultPath,GENERIC_READ|GENERIC_WRITE,FILE_SHARE_WRITE,NULL,OPEN_ALWAYS,
		                      FILE_ATTRIBUTE_NORMAL,NULL);
	HANDLE hMapResult = ::CreateFileMapping(hFileResult,NULL,PAGE_READWRITE,0,0,NULL);
	
	if(hMapResult != NULL && GetLastError() == ERROR_ALREADY_EXISTS)
	{
		CloseHandle(hMapResult);
		hMapResult = NULL;
		return ;
	}
	LPSTR lpResultHead = (LPSTR)::MapViewOfFile(hMapResult,FILE_MAP_WRITE,0,0,0);
	LPSTR lpResultDIBBits = lpResultHead + dwMulspecSizeOff;
	unsigned char* lpResultBits;       //指向融合源图象象素的指针

	LONG  lLineBytesResult = lLineBytesMulspec;

	//把SPOT全色图像和多光谱影像加权融合生成的像素值填入相应位置的像素。
	for(i=0;i<lSpotHeight;i++)
	{
		for(j=0;j<lSpotWidth;j++)
		{
			for(int k=0;k<3;k++)
			{
				lpSpotBits = (unsigned char*)lpSpotDIBBits 
					+ lLineBytesSpot*(lSpotHeight-1-i) + j;
				lpMulspecBits = (unsigned char*)lpMulspecDIBBits 
					+ lLineBytesMulspec*(lMulspecHeight-1-i) + j*3+k;
				lpResultBits = (unsigned char*)lpResultDIBBits 
					+ lLineBytesResult*(lSpotHeight-1-i) + j*3+k;
				*lpResultBits =(unsigned char) ((((1+abs(dRelation[k]))*(*lpSpotBits) 
					          + (1-abs(dRelation[k]))*(*lpMulspecBits)))/2);		
			
			}
		}
	}
	::UnmapViewOfFile(lpSpotBits);
	::CloseHandle(hMapSpot);
	::CloseHandle(hFileSpot);

	::UnmapViewOfFile(lpMulspecBits);
	::CloseHandle(hMapMulspec);
	::CloseHandle(hFileMulspec);

	::UnmapViewOfFile(lpResultBits);
	::CloseHandle(hMapResult);
	::CloseHandle(hFileResult);

	//+++++++++++++++++++++++++++++
	// 打开存储文件，用于显示
	//++++++++++++++++++++++++++++
	CFile fileM;
	CFileException feM;

	if (!fileM.Open(lpszResultPath, CFile::modeRead | CFile::shareDenyWrite, &feM))
	{
		// 失败
		ReportSaveLoadException(lpszResultPath, &feM,
			FALSE, AFX_IDP_FAILED_TO_OPEN_DOC);
		
		// 返回FALSE
		return ;
	}
	
	// 尝试调用ReadDIBFile()读取图像
	TRY
	{
		m_hFusionDIB = ::ReadDIBFile(fileM);
	}
	CATCH (CFileException, eLoad)
	{
		// 读取失败
		fileM.Abort();
		
		// 恢复光标形状
		EndWaitCursor();
		
		// 报告失败
		ReportSaveLoadException(lpszResultPath, eLoad,
			FALSE, AFX_IDP_FAILED_TO_OPEN_DOC);
		
		// 设置DIB为空
		m_hFusionDIB = NULL;
		
		fileM.Close();

		// 返回FALSE
		return ;
	}
	END_CATCH
	
	// 判断读取文件是否成功
	if (m_hFusionDIB == NULL)
	{
		// 失败，可能非BMP格式
		CString strMsg;
		strMsg = "读取图像时出错！可能是不支持该类型的图像文件！";
		
		// 提示出错
		//MessageBox(NULL, strMsg, "系统提示", MB_ICONINFORMATION | MB_OK);
		
		fileM.Close();

		// 返回FALSE
		return ;
	}

	fileM.Close();

	UpdateAllViews(NULL);
}