研究方向需要运用C/C++语言实现航空影像的处理工作，主要包括有配准和融合处理。在此主要讲一下自己学习到实现PCA融合算法的过程。起初限于自己水平，只能不断的摸索和学习。主要是借助于CSDN的资源，认真地学习了一些博主的文章，了解了PCA的原理，当自己实现时，还是会出现各种问题，下面就开始讲核心东西吧（原理、实现代码、注意事项和结果图）。

PCA融合算法原理：

       PCA变换，也称为主成分分析，是基于K-L(Kathunen-Loeve)变换来实现的。PCA变换融合处理流程如下图所示，就是将N个波段的多光谱影像进行K-L变换，并根据向量特征值排序依次得到N个主分量，再将高分辨率全色影像与第一主分量进行直方图匹配，使全色影像灰度均值和方差与第一主分量影像的保持一致；然后将匹配后的全色影像直接替换掉第一主分量影像，最后通过K-L逆变换处理将可获得高分辨率的多光谱融合影像。

以上为自己简单总结概括的原理，网上或书上对PCA融合原理介绍的会更加详细，在这里主要也可以参考以下三篇博文：

1、https://blog.csdn.net/u013165921/article/details/78339942

2、https://blog.csdn.net/deirjie/article/details/41410163

3、https://blog.csdn.net/liminlu0314/article/details/8957009

在我实现该算法过程中，这几篇博文对我有着较大的帮助。

特别是第2篇博文，里面也是用C++实现遥感影像PCA融合的，其中有完整的实现代码，所以这里我将不贴出全部代码，主要讲一下我在实现过程遇到的问题，以及需要注意的东西。

问题：

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

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

开始附代码解释：在这里很多函数基本内部有注释，也是承接上面参考博客2的。

第一步：特征值求解排序。

////////////第一步：多光谱影像之间求特征值和特征向量//////////////////////////
int width =0;
int height = 0;
int nRastercount=0;
//读入多光谱影像
unsigned char **pMSImage = GDALreadImage(width,height,nRastercount,filepath1);
//1.1计算波段均值并计算协方差矩阵  
double* covMatrix = CovMatrix(filepath1);
//1.2计算协方差所形成的矩阵的特征值与特征向量  
double eps = 0.000001;   //控制精度要求  
double *eigenVector = new double[nRastercount * nRastercount];  
eejcb( covMatrix, nRastercount, eigenVector, eps, 100000);//雅可比法求解  
//1.3按特征值由大到小的顺序排列特征向量  
sortEigenVector(eigenVector, covMatrix, nRastercount);
delete []covMatrix;covMatrix = NULL;

第二步：主成分变换并线性拉伸

///////////////////////////第二步：主成分正变换//////////////////////////
//2.1将多光谱图像进行主分量变换  
transMatrix(eigenVector,nRastercount,nRastercount);//转置
double *dPixelX = new double[nRastercount];
memset(dPixelX,0,sizeof(double)*nRastercount);
double *dPixelY = new double[nRastercount];
memset(dPixelY,0,sizeof(double)*nRastercount);
double ** result = new double*[nRastercount]; //☆☆☆
for ( int i = 0; i < nRastercount; i++ )  
{  result[i] = new double[Num];  for ( int j = 0; j < Num ; j++)  {  result[i][j] = 0;  //初始化}  
}   
unsigned char* pHRImage = OpenTif(filepath2);
//2.2将高光谱分辨率的图像进行主分量变换
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];}}
}
int length = ColorLen - 1;
double *dMax = new double[nRastercount];
double *dMin = new double[nRastercount];
for (int i = 0;i < nRastercount;i++)
{dMax[i] = -length;dMin[i] = length;
}
//统计极值
for(int i=0;i<height;i++)
{for(int j=0;j<width;j++){for(int k=0;k<nRastercount;k++){double dImageBits = (double)result[k][i*width+j];if(dMax[k]<(dImageBits))dMax[k] = dImageBits;if(dMin[k]>(dImageBits))dMin[k] = dImageBits;}}
}
//线性拉伸
for(int i=0;i<height;i++)
{for(int j=0;j<width;j++){for(int k=0;k<nRastercount;k++){double dImageBits = (double)result[k][i*width+j];if((dMax[k]-dMin[k])<length)result[k][i*width+j] = (dImageBits-dMin[k]);elseresult[k][i*width+j] = (dImageBits-dMin[k])*length/(dMax[k]-dMin[k]);}			}
}

第三步：逆变换、拉伸和输出

//////////////第三步：主分量逆变换/////////////////////////////
inverseMatrix(eigenVector,nRastercount);//求特征向量矩阵的逆矩阵。	
for(int i=0;i<height;i++)
{for(int j=0;j<width;j++){for(int k=0;k<nRastercount;k++){dPixelX[k] = result[k][i * width + j];}		multimatrix(eigenVector,dPixelX,nRastercount,nRastercount,1,dPixelY);for(int k=0;k<nRastercount;k++){result[k][i*width+j] = dPixelY[k];}}
}
delete []eigenVector;
delete []dPixelX;
delete []dPixelY;
//重新赋值
for (int i = 0;i < nRastercount;i++)
{dMax[i] = -length;dMin[i] = length;
}
//统计极值
for(int i=0;i<height;i++)
{for(int j=0;j<width;j++){for(int k=0;k<nRastercount;k++){double dImageBits = (double)result[k][i*width+j];if(dMax[k]<(dImageBits))dMax[k] = dImageBits;if(dMin[k]>(dImageBits))dMin[k] = dImageBits;}}
}
//线性拉伸
for(int i=0;i<height;i++)
{for(int j=0;j<width;j++){for(int k=0;k<nRastercount;k++){double dImageBits = (double)result[k][i*width+j];if((dMax[k]-dMin[k])<length)result[k][i*width+j] = (dImageBits-dMin[k]);elseresult[k][i*width+j] = (dImageBits-dMin[k])*length/(dMax[k]-dMin[k]);}			}
}
unsigned char** result1 = new unsigned char*[nRastercount];//采用unsigned char类型写出  
for ( int i = 0; i < nRastercount; i++ )  
{  result1[i] = new unsigned char[Num];  for ( int j = 0; j < Num ; j++ )  {  result1[i][j] = (unsigned char)result[nRastercount-i-1][j]; //????? }  
} 
delete []dMax;
delete []dMin;
for (int i = 0;i < nRastercount;i++)
{delete []result[i];result[i] = NULL;
}
delete []result;
result = NULL;

注意：

1、在融合过程中，数据类型的存储将会是造成输出结果不正确的主要原因。由于PCA融合算法中涉及到大量数值计算，里面像素灰度值值得取值将会发生变化。参考博文2内之所以不出现问题1是因为作者全部采用的float型数据进行计算，而本文读入unsigned char ，如果中间计算部分也用unsigned char数据类型将会对计算结果进行截断（unsigned char是存储整型数据），那么将导致结果的光谱色彩不正确。因此在计算中均采用double类型或float数据类型参与计算（如第二步代码中☆☆☆）。

2、该算法中理论上应当采用直方图匹配的，但是大规模影像的直方图匹配一直报错，然后采用小影像试验，发现加直方图匹配的结果并不比不加直方图匹配的结果好，相反有所不如。因此，在这里直接没有加上这一步骤。

3、上述代码是基于8bit源影像数据实验的，如果不同位深影像，则要相应的更改数据类型（比如我自己处理的航空影像为16bit影像，那么就要采用unsigned short），即影像的输入输出则根据影像自身bit数采用float、unsigned char、unsigned short等。

采用小规模遥感影像进行融合试验，得到PCA融合影像。

问题发现：

文中将中间过程输出到结果影像的过程中（数据类型转换）需要将波段颠倒输出，该问题最终发现是因为借助GDAL库读影像时，是逆波段读取的。应该根据影像相应波段进行读写操作。如果顺序读入那么写出时，代码//？？？？？处应该改为：result1[i][j] = (unsigned char)result[i][j];

但又出现了新问题：1顺序读取波段顺序输出、2顺序读取波段逆波段输出、3逆波段读取再颠倒波段输出、4逆波段读取顺序输出，产生的结果影像均不相同，理论上情况1与情况3产生的结果影像应当一样，但是结果不一致（分别为下图和上结果图）。产生的具体原因是什么呢？