fusionalgrithm.cpp

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

					*(	temp+y*WIDE+x	)=255;
			}

//3、与原全色图像合并???是加权平均呢，还是？？系数数如何确定？？
		
			for (y=0;y<height;y++)
						for (x=0;x<WIDE;x++)
						{
							panimage[y*WIDE+x]=panimage[y*WIDE+x]-temp[y*WIDE+x];
						}
					
		
//4、HIS融合??小波融合？？？



				
    delete []temp;
	delete []canny;
	delete []sobel;
	return HPAN;

}

HDIB CFusionAlgrithm::HPFFUSION(HDIB HCOLOR, HDIB HPAN)
{
	/*
	1、高分辨率全色影像与低分辨率多光谱影像各波段进行直方图匹配
	2、对进行直方图匹配后的全色影像进行高通滤波；
	3、将高通滤波后的影像分别加入多光谱各个波段；
	4、将多光谱各个波段的影像进行彩色合成为融合影像；
	*/
	
	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);
	float * temp;
	float *	R;
	float * G;
    float * B;
	float * Im;
    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];//存全色图像
    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] = 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;
		   }
		 //2、对进行直方图匹配后的全色影像进行高通滤波；
	//	 int K;
		 int r,c;
		 int SUM;
		 int HP[3][3]={-1,-1,-1,
			           -1,17,-1,
			           -1,-1,-1};
		 for(j=0;j<height-1;j++)
			 for(i=0;i<WIDE-1;i++)
			 {
				 SUM = 0;
				 /* image boundaries */
				 if(j==0 ||i==0 || j==height-1|| i==WIDE-1) 
					 SUM=panimage[j*WIDE+i];
				 
				 /* Convolution starts here */
				 else for( r=-1;r<=1;r++)//模板循环
					 for( c=-1;c<=1;c++)
					 {
						 SUM+=panimage[(j+r)*WIDE+i+c]*HP[r+1][c+1]/9;
					 }
				temp[j*WIDE+i]  =(float)SUM;//TEMP用于临时存放高斯滤波处理后的图像数据
			}
//3、将高通滤波后的影像分别加入多光谱各个波段；不加权，取两个图像和平均
			 for( j=0;j<height;j++)
				 for ( i=0;i<WIDE;i++)
				 {
					 B [j*WIDE+i]=(float)(B [j*WIDE+i]+temp[j*WIDE+i]) /2 ; 
					 G [j*WIDE+i]= (float)(G [j*WIDE+i]+temp[j*WIDE+i])/2 ;  
					 R [j*WIDE+i]= (float)(R [j*WIDE+i]+temp[j*WIDE+i])/2 ; 
					 if (*(B+j*WIDE+i)>255) *(B+j*WIDE+i)=255;
					 if (*(G+j*WIDE+i)>255) *(G+j*WIDE+i)=255;
					 if (*(R+j*WIDE+i)>255) *(R+j*WIDE+i)=255;
					 if (*(B+j*WIDE+i)<0) *(B+j*WIDE+i)=0;
					 if (*(G+j*WIDE+i)<0) *(G+j*WIDE+i)=0;
					 if (*(R+j*WIDE+i)<0) *(R+j*WIDE+i)=0;

				 } 
//4、将多光谱各个波段的影像进行彩色合成为融合影像；
				 for( j=0;j<height;j++)
					 for ( i=0;i<DibWidth;i+=3)
					 {
						 outimage[j*DibWidth+i  ]=(BYTE)B [j*WIDE+i/3];
						 outimage[j*DibWidth+i+1]=(BYTE)G [j*WIDE+i/3];
						 outimage[j*DibWidth+i+2]=(BYTE)R [j*WIDE+i/3];
					 }

    return HOUT;
}

HDIB CFusionAlgrithm::D4Fusion(HDIB grayDib, HDIB colorDib)
{
	//指向原图像的指针
	LPSTR lpgrayDIB;
	LPSTR lpgrayDIBBits;
	
	
	LPSTR lpcolorDIB;
	LPSTR lpcolorDIBBits;
	
	lpgrayDIB=(LPSTR) ::GlobalLock((HGLOBAL) grayDib);
	lpgrayDIBBits= ::FindDIBBits(lpgrayDIB);
	
	lpcolorDIB=(LPSTR) ::GlobalLock((HGLOBAL) colorDib);
	lpcolorDIBBits= ::FindDIBBits(lpcolorDIB);
	
	if(colorDib==NULL || grayDib==NULL)
	{
		AfxMessageBox("必须选择两幅影像!",MB_ICONSTOP);
		return NULL;
	}
	if(colorDib==grayDib)
	{
		AfxMessageBox("不可以出现两幅相同的影像!",MB_ICONSTOP);
		return NULL;
	}
	
	int WidthgrayDIB = (int) ::DIBWidth(lpgrayDIB);
	int HeightgrayDIB = (int) ::DIBHeight(lpgrayDIB);
	int WidthcolorDIB = (int) ::DIBWidth(lpcolorDIB);
	int HeightcolorDIB = (int) ::DIBHeight(lpcolorDIB);
	
	if(WidthgrayDIB!=WidthcolorDIB || HeightgrayDIB!=HeightcolorDIB)
	{
		AfxMessageBox("两幅影像的宽度高度必须相等!",MB_ICONSTOP);
		return NULL;
	}
	

	
	HDIB tmphDIB=(HDIB)::GlobalAlloc(GHND,WidthcolorDIB*HeightcolorDIB*3+
		sizeof(BITMAPINFOHEADER));
	if (tmphDIB == NULL)
	{
		// 分配内存失败
		return NULL;
	}
	LPSTR lpTmpDIB=(LPSTR)::GlobalLock((HGLOBAL) tmphDIB);
	memcpy(lpTmpDIB,lpcolorDIB,sizeof(BITMAPINFOHEADER)+
		WidthcolorDIB*HeightcolorDIB*3);
	LPSTR lpTmpDIBBits = ::FindDIBBits(lpTmpDIB);//save the color image
	
	float *tempR,*tempG,*tempB,*gray,*img1,*img2,*img3;
	tempR=new float[WidthcolorDIB*HeightcolorDIB];
	tempG=new float[WidthcolorDIB*HeightcolorDIB];
	tempB=new float[WidthcolorDIB*HeightcolorDIB];
	gray= new float[WidthcolorDIB*HeightcolorDIB];
	img1= new float[WidthcolorDIB*HeightcolorDIB];
	img2= new float[WidthcolorDIB*HeightcolorDIB];
	img3= new float[WidthcolorDIB*HeightcolorDIB];
	
	int x,y;
	unsigned char pixel;
	
	for(x=0; x<HeightcolorDIB; x++)
	{
		for(y=0; y<WidthcolorDIB; y++)
		{
			lpcolorDIB = (char*)lpcolorDIBBits+3*(WidthcolorDIB)*x+3*y;
			*(tempB+x*WidthcolorDIB+y) = float((unsigned char)*lpcolorDIB);
			*(tempG+x*WidthcolorDIB+y) = float((unsigned char)*(lpcolorDIB+1));
			*(tempR+x*WidthcolorDIB+y) =float((unsigned char) *(lpcolorDIB+2));
			
			lpgrayDIB = (char*)lpgrayDIBBits+WidthcolorDIB*x+y;
			pixel=(unsigned char)*lpgrayDIB;
			*(gray+x*WidthcolorDIB+y) = float(pixel);
		}
	}

	CFusionAlgrithm wave;
	wave.imageWidth=WidthcolorDIB;
	
	wave.CWaveDecompose(tempR,WidthcolorDIB,HeightcolorDIB);
	wave.CWaveDecompose(tempR,WidthcolorDIB/2,HeightcolorDIB/2);

	wave.CWaveDecompose(tempG,WidthcolorDIB,HeightcolorDIB);
	wave.CWaveDecompose(tempG,WidthcolorDIB/2,HeightcolorDIB/2);

	wave.CWaveDecompose(tempB,WidthcolorDIB,HeightcolorDIB);
	wave.CWaveDecompose(tempB,WidthcolorDIB/2,HeightcolorDIB/2);

	wave.CWaveDecompose(gray,WidthcolorDIB,HeightcolorDIB);
	wave.CWaveDecompose(gray,WidthcolorDIB/2,HeightcolorDIB/2);

	for(x=0;x<HeightcolorDIB/4;x++)
		for(y=0;y<WidthcolorDIB/4;y++)
			img1[x*WidthcolorDIB+y]=tempR[x*WidthcolorDIB+y];
	for(x=0;x<HeightcolorDIB;x++)
		for(y=0;y<WidthcolorDIB;y++)
			tempR[x*WidthcolorDIB+y]=gray[x*WidthcolorDIB+y];
	for(x=0;x<HeightcolorDIB/4;x++)
		for(y=0;y<WidthcolorDIB/4;y++)
			tempR[x*WidthcolorDIB+y]=img1[x*WidthcolorDIB+y];
				
	for(x=0;x<HeightcolorDIB/4;x++)
		for(y=0;y<WidthcolorDIB/4;y++)
			img1[x*WidthcolorDIB+y]=tempG[x*WidthcolorDIB+y];
	for(x=0;x<HeightcolorDIB;x++)
		for(y=0;y<WidthcolorDIB;y++)
			tempG[x*WidthcolorDIB+y]=gray[x*WidthcolorDIB+y];
	for(x=0;x<HeightcolorDIB/4;x++)
		for(y=0;y<WidthcolorDIB/4;y++)
			tempG[x*WidthcolorDIB+y]=img1[x*WidthcolorDIB+y];
							
	for(x=0;x<HeightcolorDIB/4;x++)
		for(y=0;y<WidthcolorDIB/4;y++)
			img1[x*WidthcolorDIB+y]=tempB[x*WidthcolorDIB+y];
	for(x=0;x<HeightcolorDIB;x++)
		for(y=0;y<WidthcolorDIB;y++)
			tempB[x*WidthcolorDIB+y]=gray[x*WidthcolorDIB+y];
	for(x=0;x<HeightcolorDIB/4;x++)
		for(y=0;y<WidthcolorDIB/4;y++)
			tempB[x*WidthcolorDIB+y]=img1[x*WidthcolorDIB+y];
		
		wave.CWaveCompose(tempR,WidthcolorDIB/2,HeightcolorDIB/2);
		wave.CWaveCompose(tempR,WidthcolorDIB,HeightcolorDIB);
		wave.CWaveCompose(tempG,WidthcolorDIB/2,HeightcolorDIB/2);
		wave.CWaveCompose(tempG,WidthcolorDIB,HeightcolorDIB);
		wave.CWaveCompose(tempB,WidthcolorDIB/2,HeightcolorDIB/2);
		wave.CWaveCompose(tempB,WidthcolorDIB,HeightcolorDIB);
		
	for(x=0; x<HeightcolorDIB; x++)
		{
			for(y=0; y<WidthcolorDIB; y++)
				{
					lpTmpDIB = (char*)lpTmpDIBBits+3*(WidthcolorDIB)*x+3*y;
					if(*(tempB+x*WidthcolorDIB+y)>255)			
						*lpTmpDIB =char(255);
					else if(*(tempB+x*WidthcolorDIB+y)<0)
						*lpTmpDIB =char(0);
					else
						*lpTmpDIB = char(*(tempB+x*WidthcolorDIB+y));
												
					if(*(tempG+x*WidthcolorDIB+y)>255)			
						*(lpTmpDIB+1)=char(255);
					else if(*(tempG+x*WidthcolorDIB+y)<0)
						*(lpTmpDIB+1) =char(0);
					else
						*(lpTmpDIB+1) = char(*(tempG+x*WidthcolorDIB+y));
												
					if(*(tempR+x*WidthcolorDIB+y)>255)
						*(lpTmpDIB+2) =char(255);
					else if(*(tempR+x*WidthcolorDIB+y)<0)
						*(lpTmpDIB+2) =char(0);
					else
						*(lpTmpDIB+2) = char(*(tempR+x*WidthcolorDIB+y));
			}
	}
	delete []img1;
	delete []img2;
	delete []img3;
	delete []tempR;
	delete []tempG;
	delete []tempB;
	delete []gray;
	::GlobalUnlock((HGLOBAL)grayDib);
	::GlobalUnlock((HGLOBAL)colorDib);
	::GlobalUnlock((HGLOBAL)tmphDIB);
	
	return tmphDIB;

}

void CFusionAlgrithm::CWaveCompose(float *TranImage, long xs, long ys)
{
//邢帅的四系数小波合成子程序
	 int i,j,k,m,n;
	float  *buffer1,*buffer2,*buffer12;
	long xsize,ysize;
	xsize=imageWidth;ysize=imageWidth;
	double k1=0,temp1=0,temp2=0;
    static float hi[4]={(-0.091506350f),0.158493649f,0.591506350f,0.341506351f};
    static float gi[4]={(-0.341506351f),0.591506350f,(-0.1584936490f),(-0.091506350f)};
    buffer1=new float [(ysize+3)*1];
    buffer12=new float [(ysize+3)*1];
    buffer2=new float [ysize*1];
    if(!(buffer1)||!(buffer2)||!(buffer12))
	{
	    AfxMessageBox("iwt:1Memory allocation error!");
	    return;
	}
	k1=hi[0]/hi[3];
	for(m=0;m<1;m++)
	{
		buffer1[m]=buffer1[1+m]=buffer1[2*1+m]=0.0;
		buffer12[m]=0.0;
		buffer12[(ys+2)*1+m]=0.0;
	}
	for(m=0;m<1;m++)
	{
		for(n=0;n<xs;n++)
		{
			for( k=0;k<(ys/2);k++)
			{
			   buffer1[(k*2+3)*1+m]=TranImage[(1L*k*xsize+n)*1+m];
			   buffer1[(k*2+3+1)*1+m]=0;
			}
			for( i=0;i<ys;i++)
			{
				 temp1=0;
				 for( j=0;j<4;j++)
					temp1+=hi[j]*buffer1[(i+j)*1+m];
				 buffer2[i*1+m]= (float)(temp1);
			}
			buffer12[(ys+1)*1+m]=(float)( (k1*TranImage[(1L*(ys/2-1)*xsize+n)*1+m]));
			for(k=0;k<(ys/2);k++)
			{
				buffer12[(k*2+1)*1+m]=TranImage[(1L*(k+ys/2)*xsize+n)*1+m];
				buffer12[(k*2+1+1)*1+m]=0.0;
			}
			for(i=0;i<ys;i++)
			{
				temp1=0;
				for(j=0;j<4;j++)
				temp1+=gi[j]*buffer12[(i+j)*1+m];
				temp2=temp1+buffer2[i*1+m];
				buffer2[i*1+m]= (float)(temp2);
			}
			for(k=0;k<ys;k++)
			TranImage[(1L*k*xsize+n)*1+m]=buffer2[k*1+m];
		}
	}
	delete [(ysize+3)*1] buffer1;
	delete [(ysize+3)*1] buffer2;
	delete [ysize*1] buffer12;
	buffer1=(float *)malloc((xsize+3)*sizeof(float)*1);
	buffer12=(float *)malloc((xsize+3)*sizeof(float)*1 );
	buffer2=(float *)malloc(xsize*sizeof(float) *1);
	if(!(buffer1)||!(buffer2)||!(buffer12))
	{
		 AfxMessageBox("iwt:2Memory allocation error!");
		 return;
	}
	for(m=0;m<1;m++)
	{
		buffer1[m]=buffer1[1+m]=buffer1[2*1+m]=0.0;
		buffer12[m]=0.0;
		buffer12[(xs+2)*1+m]=0.0;
	}
	for(m=0;m<1;m++)
	{
		for(n=0;n<ys;n++)
		{
			for(k=0;k<(xs/2);k++)
			{
				buffer1[(k*2+3)*1+m]=TranImage[(1L*n*xsize+k)*1+m];
				buffer1[(k*2+3+1)*1+m]=0;
			}
			for(i=0;i<xs;i++)
			{
				temp1=0;
				for(j=0;j<4;j++)
				 temp1+=hi[j]*buffer1[(i+j)*1+m];
				buffer2[i*1+m]=(float)( temp1);
			}
			buffer12[(xs+1)*1+m]= (float)((k1*TranImage[(1L*n*xsize+xs/2-1)*1+m]));
			for(k=0;k<(xs/2);k++)
			{
				buffer12[(k*2+1)*1+m]=TranImage[(1L*n*xsize+xs/2+k)*1+m];
				buffer12[(k*2+1+1)*1+m]=0;
			}
			for(i=0;i<xs;i++)
			{
				temp1=0;temp2=0.0;
				for(j=0;j<4;j++)
				temp1+=gi[j]*buffer12[(i+j)*1+m];
				temp2=temp1+buffer2[i*1+m];
				buffer2[i*1+m]=(float)( temp2);
			}

			for(k=0;k<xs;k++)
			TranImage[(1L*n*xsize+k)*1+m]=buffer2[k*1+m]*4;
		}
	}
	free(buffer1);
	free(buffer2);
	free(buffer12);
}

HDIB CFusionAlgrithm::PCAFUSION(HDIB HCOLOR, HDIB HPAN)
{
return NULL;
}