图像处理——优化PCA图像融合文中问题解答及相关优化

Posted 2021-12-18 nanke_yh

        在很早之前就发过一篇文章：《图像处理系列——图像融合之主成分分析（PCA）》，现在这篇文就是在它的内容上进行的。上文没看过的可以去看看，这样再看这一篇文章时才更容易明白和理解。

        上一文章中基本实现了基于主成分分析的图像融合主要的功能，也获得了结果图，但在文章中提出了很多的问题。一直以来偶尔思考但也没花时间和精力给他重新梳理完善。这篇文章的诞生还是前些天写完图像处理系列——直方图之直方图规定化(Histogram Specification)一文后，才想着重新梳理一下当初写过的PCA融合内容的。

        在PCA融合文章中曾提到过几个问题：

1、根据 Jacory_Gao作者提供的代码，实现影像融合时，出现结果影像色彩紊乱的情况（具体描述为：水彩画，涂鸦）；

2、对于大影像而言，直方图匹配部分报错；

3、逆变换后得到的多波段影像，按照不同顺序输出产生的结果影像均不相同，产生的具体原因是什么呢？

当时文中也有对应的解释，这里再重新的梳理和解答一下。

第一问：

        在融合过程中，数据类型的存储将会是造成输出结果不正确的主要原因。由于PCA融合算法中涉及到大量数值计算，里面像素灰度值计算取值将会发生变化。一般影像读入数据类型不是unsigned char 就是unsigned short等，如果中间计算部分也用unsigned char等数据类型将会对计算结果进行取整截断（unsigned char是存储整型数据），那么将导致结果的光谱色彩不正确。因此在中间的计算过程均采用double类型或float数据类型来保存传递内存值。（这一点前文中也说明白了，在此多提一句而已）。

第二问：

        前文中对于直方图匹配的处理方式：采用小影像试验，发现加直方图匹配的结果并不比不加直方图匹配的结果好，相反有所不如。因此，在这里直接没有加上这一步骤。

这样虽然能够得到结果，但不是特别的严谨。如果对于影像质量相差特别大的影像来说，如此处理就会出现问题。为了重新验证和排除问题，继续在融合工程中进行了完善：

先放原图：

全色影像	多光谱影像

 强烈吐槽一下：为什么写文章内上传图像不能上传tif格式的图片了，那搞的我只能截图了（不想转格式，这样还有直方图信息）。

        我们知道，多光谱影像经过PCA正变换后，得到double类型的各个成分图层，那么将全色影像与第一主成分进行直方图匹配就需要分别统计两者的直方图累积概率分布。在这个过程对于double图层则需要保留精度，factor选择100，因为选值太大会消耗大量的时间，同时色彩并无太大影响。

首先，对于大规模的数组问题：不可定义成临时数组，因为临时变量存在栈区，是有着最大容量限制的，一般是2M。无法满足大规模数组的存储，为此需要使用动态内存来实现，否则会报错；

再次，对于正变换后的第一主成分，其值域范围已经发生了变化，我们可以看到输出最值情况如下，不可还是按照0~255来，应该先找到其最值，最大值max*factor才是值域右侧；

然后，创建数组内存空间时，不可按照max*factor的大小来，否则如下报错，应该length=max*factor+1，因为需要考虑0值索引值就要大一个;

 这个问题是真的不容易发现，困扰我很久，多方调试后最终才定位排查出来的。其现象是，开辟内存后，向其中写入过多数值后，内存释放时报错堆被损坏。记住：内存方面报错，往往就是内存创建的空间小于实际内存，增大创建空间即可。

最后，解决这些问题后，实现一下直方图匹配：

        这里面设计是全色源图像色阶不变，仅仅扩充第一主成分色阶，再分别统计直方图，累积分布，最后映射即可。而区别于C++实现遥感图像PCA融合_Jacory_Gao的专栏-CSDN博客中将源图像和目标图像求最值，拉伸至一致，再拷贝扩大色阶统计直方图，匹配。如此可减少大量的计算过程和内存消耗。

double* HistImgMatch(unsigned char*srcImg,double *destImg,int Width,int Height)

	if (srcImg == NULL || destImg == NULL)
	
		return NULL;
	
	double min,max;
	max_min_value(destImg,Width*Height,max,min);
	cout<<max<<","<<min<<endl;
	int accuracy = 100;//设置精度
	double* ImageMatch = new double[Width*Height];
	memset(ImageMatch,0,sizeof(double)*Width*Height);
	int Len = (int)(max*accuracy+0.5);//目标图像能够取到到的最大值
	int Length = Len+1;//创建内存，需要考虑到总数（包括0），不然报错

	double *srchist = CalImgHist(srcImg,Width,Height,ColorLen,1);
	double*desthist = CalImgHist1(destImg,Width,Height,Length,accuracy);
	
	//原图累计分布直方图
	//double inc_srchist[ColorLen*1000] =  0 ;
	double *inc_srchist = new double[ColorLen];
	memset(inc_srchist,0,sizeof(double)*ColorLen);
	double temp = 0;
	for (int i = 0; i < ColorLen; i++)
	
		temp += srchist[i];
		inc_srchist[i] = temp;
	
	delete []srchist;srchist=NULL;
	//目标图累计分布直方图
	double *inc_desthist = new double[Length];
	memset(inc_desthist,0,sizeof(double)*Length);
	temp = 0;
	for (int i = 0; i < Length; i++)
	
		temp += desthist[i];
		inc_desthist[i] = temp;
	
	delete []desthist;desthist=NULL;
	
	//单映射关系
	int *SML = new int[ColorLen];//映射关系内存大小与源图像空间大小一致
	memset(SML,0,sizeof(int)*ColorLen);
	for (int i = 0; i < ColorLen; i++)
	
		double minvalue = 1;
		double orig_value = inc_srchist[i];
		for (int j = 0; j < Length; j++)
		
			double stan_value = inc_desthist[j];
			double srcmin = fabs(stan_value - orig_value);//与目标的累积直方图最近的点
			if (srcmin < minvalue)
			
				minvalue = srcmin;
				SML[i] = j;
			
		
	
	delete []inc_srchist;inc_srchist=NULL;
	delete []inc_desthist;inc_desthist=NULL;
	for (int i = 0; i < Height; i++)
	
		for (int j = 0; j < Width; j++)
		
			ImageMatch[i*Width+j] = SML[srcImg[i*Width+j]]/accuracy;
		
	

	delete []SML;SML=NULL;
	return ImageMatch;

下面计算影像直方图函数：

double* CalImgHist1(double* Image,int Width,int Height,int length,int offset)

	if (Image == NULL)
	
		return NULL;
	
	double* ImgHist = new double[length];
	memset(ImgHist,0,sizeof(double)*length);//涉及值的传递问题，数据存储四区，必须创建内存存储

	for (int i = 0; i < Height; i++)
	
		for (int j = 0; j < Width; j++)
		
			int temp = Image[i*Width+j]*offset+0.5;
			ImgHist[temp]++;//统计各个灰度级的数量
			//cout<<temp << ":"<<ImgHist[temp]<<endl;
		
	
	for (int i = 0; i < length; i++)
	
		ImgHist[i] = ImgHist[i] / (Width*Height);//归一化求各个灰度级出现频率
	
	return ImgHist;

double* CalImgHist(unsigned char* Image,int Width,int Height,int length,int offset)

	if (Image == NULL)
	
		return NULL;
	
	//double ImgHist[ColorLen] =  0 ;//0~255每个像素值出现的概率
	double* ImgHist = new double[length];
	memset(ImgHist,0,sizeof(double)*length);//涉及值的传递问题，数据存储四区，必须创建内存存储
	for (int i = 0; i < Height; i++)
	
		for (int j = 0; j < Width; j++)
		
			ImgHist[Image[i*Width+j]*offset]++;//统计各个灰度级的数量
		
	
	for (int i = 0; i < length; i++)
	
		ImgHist[i] = ImgHist[i] / (Width*Height);//归一化求各个灰度级出现频率
	
	return ImgHist;

可以看出，上面两个函数很相似，因为参数数据类型不同而需要定义多个函数，这个其实可通过函数模板来解决优化，后续将继续研究。

最后替换第一主分量，再逆变换就完成了PCA工作了。

上一篇文章中替换：是边正变换边替换的
for(int i=0;i<height;i++)

    for(int j=0;j<width;j++)
    
        for(int k=0;k<nRastercount;k++)
        
            dPixelX[k] = pMSImage[k][i * width + j];
        
        //主分量变换
        multimatrix(eigenVector,dPixelX,nRastercount,nRastercount,1,dPixelY);//矩阵相乘
        dPixelY[0] = pHRImage[i*width+j];//将高分辨率图像替代第一主分量
        for(int k=0;k<nRastercount;k++)
        
            result[k][i*width+j] = dPixelY[k];
        
    

————————————————
版权声明：本文为CSDN博主「nanke_yh」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/nanke_yh/article/details/89481095

下面则改为：

for(int i=0;i<height;i++)

	for(int j=0;j<width;j++)
	
		for(int k=0;k<nRastercount;k++)
		
			dPixelX[k] = pMSImage[k][i * width + j];
		
		//主分量变换
		multimatrix(eigenVector,dPixelX,nRastercount,nRastercount,1,dPixelY);//矩阵相乘
		//dPixelY[0] = pHRImage[i*width+j];//将高分辨率图像替代第一主分量
		for(int k=0;k<nRastercount;k++)
		
			result[k][i*width+j] = dPixelY[k];
		
	

double*matchimg = HistImgMatch(pHRImage,result[0],width,height);//直接替换第一分量
for(int i=0;i<height;i++)

	for(int j=0;j<width;j++)
	
		result[0][i*width+j] = matchimg[i*width+j];
	

delete []matchimg;matchimg=NULL;
delete []pHRImage;pHRImage=NULL;

得到的结果：

全色直接替换	直方图匹配后替换

         可以看出，两者的色彩还是有些区别的。从直方图上看：直方图匹配后替换得到的影像波峰偏右整体色彩更亮。从色彩上看：右图色彩的对比度也更好一些，比如再右上角的山区部分，右图山脊山沟更明显；左上角的水域与林地色差区别更好。最后从与原图对比来看：右图色彩分布更贴合多光谱影像，色彩保留较好。而左图除了提高了分辨率外，色彩也受到了全色影像的影响。 

 第三问：

        不管按怎么样的波段顺序读入的影像，经过PCA分析得到的是各个主成分，然后再排序，得到依次递减的“波段”，如此不同的输出顺序对应的不同的RGB通道上就是会出现不同的结果影像。（这一点在上文的评论中也给出了解释，在此提上一句），一般是不用纠结其如何输出的。一般是按照顺序输出到RGB通道上即可。

result1[i][j] = (unsigned char)result[i][j];