fusionalgrithm.cpp

影像融合与融合精度评价源码影像融合与融合精度评价源码

				//将灰度值赋给新数组的该像素* /
				gray = (BYTE)img1[i*width+j];
				imgResult[i*width+j] = (double)cpd[gray];
			}
			
			delete []histo1;
			delete []histo2;
			delete []cpd;
			delete []p1;
			delete[]p2;
			delete []f;
			delete []g;
			
			return imgResult;
}

HDIB CFusionAlgrithm::AVERAGEFUSION(HDIB HCOLOR, HDIB HPAN)
{

	LPSTR  muldib,pandib,outdib; //指向图像的指针
	BYTE * colorimage;//指向彩色图像的数据区位置
	BYTE * panimage;//指向全色图像的数据区位置
	muldib= (LPSTR) ::GlobalLock((HGLOBAL) HCOLOR);
	colorimage=(BYTE *)FindDIBBits(muldib);
	pandib= (LPSTR) ::GlobalLock((HGLOBAL) HPAN);
	panimage =(BYTE *)FindDIBBits(pandib);	

	int wide,height,DibWidth,WIDE;
	wide=::DIBWidth(muldib);//取宽的像素个数
	height=::DIBHeight(muldib);
	DibWidth=3*wide+( sizeof(DWORD)- 3*wide %sizeof(DWORD) )%sizeof(DWORD);//取得原图的每行总位数
	WIDE=((wide*8)+31)/32*4;//取宽为4的整数倍
	LONG dwPixelNum=WIDE*height;  //数据像素数
	BYTE * outimage;  //输出图像的数据区位置
	HDIB HOUT;
	HOUT=(HDIB)::GlobalAlloc(GHND,sizeof(BITMAPINFOHEADER)+DibWidth*height);//生成一个新的图像句柄
	outdib	= (LPSTR)::GlobalLock((HGLOBAL)HOUT);
	memcpy(outdib,muldib,sizeof(BITMAPINFOHEADER));
	outimage= (BYTE *) FindDIBBits (outdib);
		
//参数检测	
	ASSERT(panimage!=NULL);
	ASSERT(colorimage!=NULL);
	ASSERT(outimage!=NULL);

//统计平均灰度
	int i,j;
	int dPanchroMeanGray = 0; //全色平均值
	LONG red=0,green=0,blue=0;
	
	for( i=0;i<dwPixelNum;i++)
	{
		dPanchroMeanGray += panimage[i];// 全色求和
	}
	for( j=0;j<dwPixelNum;j+=3)
	{
		blue  += colorimage[3*j];   //蓝求和
		green += colorimage[3*j+1];  //绿求和
		red   += colorimage[3*j+2];    //红求和
	}
	dPanchroMeanGray/= dwPixelNum;   //全色中值
	blue /= dwPixelNum;   //蓝色中值
	green /= dwPixelNum;  //绿色中值
	red /= dwPixelNum;    //红色中值
	
	
//计算相关系数			
	double dRelationRed=0 ,dRelationGreen=0, dRelationBlue=0;//相关系数初始化
	double dTemp1=0,dTemp2=0,dTemp3=0;
	double temp1=0,temp2=0;

/////////蓝色相关系数
	for(i=0;i<dwPixelNum;i++)
	{
		dTemp1 += (panimage[i] - blue) * 
			(colorimage[i] - blue);
		dTemp2 += (panimage[i] - blue) * 
			(panimage[i] - blue);
		dTemp3 += (colorimage[i] - blue) * 
			(colorimage[i] - blue);
	}		 

	temp1=sqrt(dTemp2);
	temp2=sqrt(dTemp3);
	dRelationBlue = (fabs)(dTemp1/(temp1*temp2));
	temp1=temp2=0;
	dTemp1 = dTemp2 = dTemp3 = 0;//临时系数归零
/////////绿色相关系数
	for(i=0;i<dwPixelNum;i++)
	{
		dTemp1 += (panimage[i] - green) * 
			(colorimage[i] - green);
		dTemp2 += (panimage[i] - green) * 
			(panimage[i] - green);
		dTemp3 += (colorimage[i] - green) * 
			(colorimage[i] - green);
	}	
	temp1=sqrt(dTemp2);
	temp2=sqrt(dTemp3);
	dRelationGreen = (fabs)(dTemp1/(temp1*temp2));
	temp1=temp2=0;
	dTemp1 = dTemp2 = dTemp3 = 0;//临时系数归零
////////红色相关系数
	for(i=0;i<dwPixelNum;i++)
	{
		dTemp1 += (panimage[i] - red) * 
			(colorimage[i] - red);
		dTemp2 += (panimage[i] - red) * 
			(panimage[i] - red);
		dTemp3 += (colorimage[i] - red) * 
			(colorimage[i] - red);
	}
	temp1=sqrt(dTemp2);
	temp2=sqrt(dTemp3);
	dRelationRed = (fabs)(dTemp1/(temp1*temp2));
	temp1=temp2=0;
	dTemp1 = dTemp2 = dTemp3 = 0;//临时系数归零

/////////融合/////////////
	
	for(i=0;i<dwPixelNum;i++)
	{
	    outimage[3*i]  = (BYTE)((1-dRelationBlue)*panimage[i] + 
			(0.5+dRelationBlue/2)*colorimage[3*i]);
		outimage[3*i+1]=(BYTE)((1-dRelationGreen)*panimage[i] + 
			(0.5+dRelationGreen/2)*colorimage[3*i+1]);
		outimage[3*i+2]=(BYTE)((1-dRelationRed)*panimage[i] + 
			(0.5+dRelationRed/2)*colorimage[3*i+2]);
	}
	return (HOUT);

}


//邢帅的四系数小波分解子程序
void CFusionAlgrithm::CWaveDecompose(float *TranImage, long xs, long ys)
{
	
	int i,j,k,m,n;
    float  *bufferx,*buffery,*temp_bufx,*temp_bufx1;
	long xsize=imageWidth;
	double k1=0,k2=0,temp1=0,temp2=0;
    static	float h[4]={0.341506351f,0.591506350f,0.158493649f,(-0.091506350f)};
    static  float g[4]={(-0.091506350f),(-0.1584936490f),0.591506350f,(-0.341506351f)};
    bufferx=new float [xs*1];
    temp_bufx=new float [(xs+4)*1];
    temp_bufx1=new float [(xs+2)*1];
	
	if(!(bufferx)||!(temp_bufx)||!(temp_bufx1))
	{
		AfxMessageBox("WT:1Memory allocation error!");
		return;
	}
	/*----- To filter in the horizonal direction-----*/
	k1=(h[2]/h[0]);
	k2=(h[3]/h[0]);
	for(m=0;m<1;m++)
	{
		temp_bufx[1+m]=0.0;
		temp_bufx[(xs+2)*1+m]=temp_bufx[(xs+3)*1+m]=0.0;
		temp_bufx1[xs*1+m]=temp_bufx1[(xs+1)*1+m]=0.0;
	}
	for(m=0;m<1;m++)
	{
		for(n=0;n<ys;n++)
		{
			temp_bufx[m]=(float)( (-k1*TranImage[(1L*n*xsize)*1+m]-k2*TranImage[(1L*n*xsize+1)*1+m]));
			for( i=0;i<xs;i++)
			{
				temp_bufx[(i+2)*1+m]=TranImage[(1L*n*xsize+i)*1+m];
				temp_bufx1[i*1+m]=temp_bufx[(i+2)*1+m];
			}
			for(i=0;i<xs;i+=2)
			{
				temp1=0.0;
				temp2=0.0;
				for(j=0;j<4;j++)
				{
					temp1=temp1+h[j]*temp_bufx1[(i+j)*1+m];
					temp2=temp2+g[j]*temp_bufx[(i+j)*1+m];
				}
				bufferx[(i/2)*1+m]=(float)( temp1);
				bufferx[(i/2+xs/2)*1+m]= (float)(temp2);
			}
			for( k=0;k<xs;k++)
				TranImage[(1L*n*xsize+k)*1+m]=bufferx[k*1+m];
		}
	}
	
	delete [xs*1] bufferx;
	delete [(xs+4)*1] temp_bufx;
	delete [(xs+2)*1] temp_bufx1;
	
	buffery=new float [ys*1];
	temp_bufx=new float [(ys+4)*1];
	temp_bufx1=new float [(ys+2)*1];
	if(!(buffery)||!(temp_bufx)||!(temp_bufx1))
	{
		AfxMessageBox("WT:2Memory allocation error!");
		return;
	}
	/*----- To filter in the vertical direction-----*/
	for(m=0;m<1;m++)
	{ 
		temp_bufx[(ys+2)*1+m]=temp_bufx[(ys+3)*1+m]=0.0;
        temp_bufx1[ys*1+m]=temp_bufx1[(ys+1)*1+m]=0.0;
	}
	
	for(m=0;m<1;m++)
	{
		for(n=0;n<xs;n++)
		{
			temp_bufx[m]= (float)((-k1*TranImage[(1L*n)*1+m]-k2*TranImage[(1L*xsize+n)*1+m])); 
			temp_bufx[1+m]=0;
			for(i=0;i<ys;i++)
			{
				temp_bufx[(i+2)*1+m]=TranImage[(1L*i*xsize+n)*1+m];
				temp_bufx1[i*1+m]=temp_bufx[(i+2)*1+m];
			}
			for(i=0;i<ys;i+=2)
			{
				temp1=0.0;
				temp2=0.0;
				for(j=0;j<4;j++)
				{
					temp1+=h[j]*temp_bufx1[(i+j)*1+m];
					temp2+=g[j]*temp_bufx[(i+j)*1+m];
				}
				buffery[(i/2)*1+m]=(float)( temp1);
				buffery[(i/2+ys/2)*1+m]= (float)(temp2);
			}
			for(k=0;k<ys;k++)
				TranImage[(1L*k*xsize+n)*1+m]=buffery[k*1+m];
		}
	}
	delete [ys*1] buffery;
	delete [(ys+4)*1] temp_bufx;
	delete [(ys+2)*1] temp_bufx1;
}


HDIB CFusionAlgrithm::DAVERAGEFUSION(HDIB HCOLOR, HDIB HPAN)
{
	LPSTR  muldib,pandib,outdib; //指向图像的指针
	BYTE * colorimage;//指向彩色图像的数据区位置
	BYTE * panimage;//指向全色图像的数据区位置
	muldib= (LPSTR) ::GlobalLock((HGLOBAL) HCOLOR);
	colorimage=(BYTE *)FindDIBBits(muldib);
	pandib= (LPSTR) ::GlobalLock((HGLOBAL) HPAN);
	panimage =(BYTE *)FindDIBBits(pandib);	
	
	int wide,height,DibWidth,WIDE;
	wide=::DIBWidth(muldib);//取宽的像素个数
	height=::DIBHeight(muldib);
	DibWidth=3*wide+( sizeof(DWORD)- 3*wide %sizeof(DWORD) )%sizeof(DWORD);//取得原图的每行总位数
	WIDE=((wide*8)+31)/32*4;//取宽为4的整数倍
	LONG dwPixelNum=WIDE*height;  //数据像素数
	BYTE * outimage;  //输出图像的数据区位置
	HDIB HOUT;
	HOUT=(HDIB)::GlobalAlloc(GHND,sizeof(BITMAPINFOHEADER)+DibWidth*height);//生成一个新的图像句柄
	outdib	= (LPSTR)::GlobalLock((HGLOBAL)HOUT);
	memcpy(outdib,muldib,sizeof(BITMAPINFOHEADER));
	outimage= (BYTE *) FindDIBBits (outdib);
	for(int j=0; j<height; j++)
		for(int i=0; i<WIDE; i++)
		{
			outimage[j*DibWidth+i*3]  =(colorimage[j*DibWidth+i*3]  +panimage[j*WIDE+i])/2;
			outimage[j*DibWidth+i*3+1]=(colorimage[j*DibWidth+i*3+1]+panimage[j*WIDE+i])/2;
			outimage[j*DibWidth+i*3+2]=(colorimage[j*DibWidth+i*3+2]+panimage[j*WIDE+i])/2;

		}
return HOUT;


}

HDIB CFusionAlgrithm::BROVEYFUSION(HDIB HCOLOR, HDIB HPAN)
{
	LPSTR  muldib,pandib,outdib; //指向图像的指针
	BYTE * colorimage;//指向彩色图像的数据区位置
	BYTE * panimage;//指向全色图像的数据区位置
	float * Temp;
	float *	R;
	float * G;
    float * B;
	float *	RP;//存R'
	float * GP;//存G'
    float * BP;//存B'
	float * Im;
	muldib= (LPSTR) ::GlobalLock((HGLOBAL) HCOLOR);
	colorimage=(BYTE *)FindDIBBits(muldib);
	pandib= (LPSTR) ::GlobalLock((HGLOBAL) HPAN);
	panimage =(BYTE *)FindDIBBits(pandib);	
	
	int wide,height,DibWidth,WIDE;
	wide=::DIBWidth(muldib);//取宽的像素个数
	height=::DIBHeight(muldib);
	DibWidth=3*wide+( sizeof(DWORD)- 3*wide %sizeof(DWORD) )%sizeof(DWORD);//取得原图的每行总位数
	WIDE=((wide*8)+31)/32*4;//取宽为4的整数倍
	LONG dwPixelNum=WIDE*height;  //数据像素数
	BYTE * outimage;  //输出图像的数据区位置
	HDIB HOUT;
	HOUT=(HDIB)::GlobalAlloc(GHND,sizeof(BITMAPINFOHEADER)+DibWidth*height);//生成一个新的图像句柄
	outdib	= (LPSTR)::GlobalLock((HGLOBAL)HOUT);
	memcpy(outdib,muldib,sizeof(BITMAPINFOHEADER));
	outimage= (BYTE *) FindDIBBits (outdib);
	// compute here
     int i,j;
	Temp = new float[WIDE*height];
  	R   = new float[WIDE*height];
	G   = new float[WIDE*height];
	B   = new float[WIDE*height];
	Im  = new float[WIDE*height];//存全色图像
	RP   = new float[WIDE*height];
	GP   = new float[WIDE*height];
	BP   = new float[WIDE*height];
	for( j=0;j<height;j++)
		for ( i=0;i<DibWidth;i+=3)
		{
			B [j*WIDE+i/3]=colorimage[j*DibWidth+i]  ;  //得到蓝值   
			G [j*WIDE+i/3]=colorimage[j*DibWidth+i+1];  //得到绿值  
			R [j*WIDE+i/3]=colorimage[j*DibWidth+i+2];  //得到红值
			Temp[j*WIDE+i/3]=B [j*WIDE+i/3]+G [j*WIDE+i/3]+R [j*WIDE+i/3];
			Im[(DWORD)j*WIDE+i/3] = (float)panimage[(DWORD)j*WIDE+i/3];
			
		}
		float * imB=new float[WIDE*height];
			HistogramMatch(R,B,WIDE,height,imB);//用B分量图像对全色图像进行直方图均衡
		float * imG= new float[WIDE*height];
			HistogramMatch(G,imB,WIDE,height,imG);
		float * imh= new float[WIDE*height];
			HistogramMatch(Im,imG,WIDE,height,imh);
		if(imh == NULL)
		{
			AfxMessageBox("直方图均衡失败!",MB_ICONSTOP);
			return NULL;
		}
		
		for( j=0;j<height;j++)
			for ( i=0;i<DibWidth;i+=3)
			{ 
				BP[j*WIDE+i/3]=3*B [j*WIDE+i/3]*imh[j*WIDE+i/3]/Temp[j*WIDE+i/3];
				GP[j*WIDE+i/3]=3*G [j*WIDE+i/3]*imh[j*WIDE+i/3]/Temp[j*WIDE+i/3];
				RP[j*WIDE+i/3]=3*R [j*WIDE+i/3]*imh[j*WIDE+i/3]/Temp[j*WIDE+i/3];
				if(BP[j*WIDE+i/3]>255) BP[j*WIDE+i/3]=255;
                if(GP[j*WIDE+i/3]>255) GP[j*WIDE+i/3]=255;
				if(RP[j*WIDE+i/3]>255) RP[j*WIDE+i/3]=255;
                if(BP[j*WIDE+i/3]<0) BP[j*WIDE+i/3]=(-1)*BP[j*WIDE+i/3];
				if(GP[j*WIDE+i/3]<0) GP[j*WIDE+i/3]=(-1)*GP[j*WIDE+i/3];
				if(RP[j*WIDE+i/3]<0) RP[j*WIDE+i/3]=(-1)*RP[j*WIDE+i/3];
                outimage[j*DibWidth+i]  =(BYTE)BP[j*WIDE+i/3];
				outimage[j*DibWidth+i+1]=(BYTE)GP[j*WIDE+i/3];
				outimage[j*DibWidth+i+2]=(BYTE)RP[j*WIDE+i/3];

			}
			free(BP);
			free(RP);
			free(GP);
			free(Temp);
			free(R);
			free(G);
			free(B);
			free(Im);

	return HOUT;

}

HDIB CFusionAlgrithm::CANNYSOBELFUSION(HDIB HCOLOR, HDIB HPAN)
{
/*1、对全色图像分别进行CANNY和SOBEL特征提取，
  2、对提取的结果进行逻辑与运算
  3、与原全色图像合并
  4、合并图像与多光谱图像进行融合（小波融合？HIS融合？？？）；
  */
	LPSTR  muldib,pandib,outdib; //指向图像的指针
	BYTE * colorimage;//指向彩色图像的数据区位置
	BYTE * panimage;//指向全色图像的数据区位置
	//锁定DIB
	muldib= (LPSTR) ::GlobalLock((HGLOBAL) HCOLOR);
	colorimage=(BYTE *)FindDIBBits(muldib);
	pandib= (LPSTR) ::GlobalLock((HGLOBAL) HPAN);
	panimage =(BYTE *)FindDIBBits(pandib);
	
	int wide,height,DibWidth,WIDE;
	wide=::DIBWidth(muldib);//取宽的像素个数
	height=::DIBHeight(muldib);
    DibWidth=3*wide+( sizeof(DWORD)- 3*wide %sizeof(DWORD) )%sizeof(DWORD);//取得原图的每行总位数
	WIDE=((wide*8)+31)/32*4;//取宽为4的整数倍
	
	BYTE * outimage;  //输出图像的数据区位置
	HDIB HOUT;
	
	HOUT=(HDIB)::GlobalAlloc(GHND,sizeof(BITMAPINFOHEADER)+DibWidth*height);//生成一个新的图像句柄
	if(!HOUT){	return NULL;	}
	
	outdib	= (LPSTR)::GlobalLock((HGLOBAL)HOUT);
	memcpy(outdib,muldib,sizeof(BITMAPINFOHEADER));//sizeof(BITMAPINFOHEADER)=40
	outimage= (BYTE *) FindDIBBits (outdib);
	
	if(muldib==NULL || pandib==NULL)
	{
		AfxMessageBox("必须选择两幅影像!",MB_ICONSTOP);	
		return NULL;
	}
	if(muldib==pandib)
	{
		AfxMessageBox("不可以出现两幅相同的影像!",MB_ICONSTOP);
		return NULL;
	}
	if(::DIBWidth(muldib)!=::DIBWidth(pandib) || ::DIBHeight(muldib)!= ::DIBHeight(pandib))
	{
		AfxMessageBox("两幅影像的宽度高度必须相等!",MB_ICONSTOP);
		return NULL;
	}
	
   	BYTE *	canny;//用于CANNY 算子处理后图像
	BYTE * sobel;//用于SOBEL 算子处理后图像
    BYTE * temp;//用于逻辑与处理后图像
	temp  = new BYTE[WIDE*height]     ;
	canny = new BYTE[WIDE*height]     ;
	sobel = new BYTE[WIDE*height]     ;

	// 指向梯度数据的指针
	double * pdGrad;
	
	// 按照图像的大小开辟内存空间，存储梯度计算的结果
	pdGrad=new double[height*WIDE];

	//1、对全色图像分别进行CANNY和SOBEL特征提取
    CFeatureAlgrithm feature;
    feature.CCanny(panimage,WIDE,height,0.4,0.4,0.79,canny);
	feature.SobelOperator(pandib,pdGrad);
    int x,y;
	for( y=0; y<height ; y++ )
		for( x=0 ; x<WIDE ; x++ )
		{
			canny[y*WIDE+x]=(unsigned char)(255-canny[y*WIDE+x]);
			if(*(pdGrad+y*WIDE+x)>50)
				*(	sobel+y*WIDE+x	)=0;
			else
				*(	sobel+y*WIDE+x	)=255;			
		}
//2、对提取的结果进行逻辑与运算
		for(y=0;y<height;y++)
			for(x=0;x<WIDE;x++)
			{
				if(canny[y*WIDE+x]==0&&sobel[y*WIDE+x]==0)
					*(	temp+y*WIDE+x	)=0;
				else