模式识别第三次作业-KL变换.doc

遥感信息工程学院PRS032班 黄寅捷 200332590083模式识别作业三一、编程要求：编程实现KL变换，并对TM六波段图像进行演算。KL变换的思想是：从n维特征选取m维特征，删去的n-m维特征不一定就是无用的信息，如何在信息损失最小的情况下选取特征，在理论上就显得更严密些。通常采用正交变换，得到新的经变换的模式， 以保证信息损失最小情况下获得有利于分类的特征。二、编程思想：将6副图象依次输入获得灰度值存在一个6*size(size为一副图象的像素数)的二维数组中，计算每个波段的灰度均值，然后计算协方差矩阵，得出特征值和特征向量，再重新计算要输出的处理的图象的灰度。三、核心代码：void CKLchangeView:Oninput() /TM图象数据读入AfxMessageBox(注意:请一次将6副图象都读入!);/读图象的方法int flag=0;while(flag6)CFileDialog dlg(TRUE,bmp,.bmp);if(dlg.DoModal()=IDOK)m_NewFileNameflag=dlg.GetFileName();for(int j=0;j5;j+)if(flag=j|m_NewFileNameflag!=m_NewFileNamej) continue;else AfxMessageBox(已打开过的图象!请按方法重新打开!);flag=0; for(int i=0;i6;i+) m_NewFileNamei=;return;NewFileName=m_NewFileNameflag;flag+;/判断是否是已经打开的图象，若已经打开，则要重新按方法打开else /不读图象，或是只是要看图象for(int i=0;ipDib=&m_pDib;OnDraw(pDC);/调用OnDraw函数绘图sizeImage = m_pDib.GetDimensions(); sizeSaveImage = m_pDib.GetDibSaveDim(); a=sizeSaveImage.cx * sizeSaveImage.cy;/象素数if(flag=1) tm=CMatrix(6,a);/开第一副图象时为矩阵分配内存 /将每幅图象象素灰度值放入矩阵中if(flag=1)for(int i=0;ia;i+)tm.m_pDatai=*(m_pDib.m_lpImage+sizeSaveImage.cx*(i/sizeImage.cx)+i%sizeImage.cx);if(flag=2)for(int i=0;ia;i+)tm.m_pDataa+i=*(m_pDib.m_lpImage+sizeSaveImage.cx*(i/sizeImage.cx)+i%sizeImage.cx);if(flag=3)for(int i=0;ia;i+)tm.m_pData2*a+i=*(m_pDib.m_lpImage+sizeSaveImage.cx*(i/sizeImage.cx)+i%sizeImage.cx);if(flag=4)for(int i=0;ia;i+)tm.m_pData3*a+i=*(m_pDib.m_lpImage+sizeSaveImage.cx*(i/sizeImage.cx)+i%sizeImage.cx);if(flag=5)for(int i=0;ia;i+)tm.m_pData4*a+i=*(m_pDib.m_lpImage+sizeSaveImage.cx*(i/sizeImage.cx)+i%sizeImage.cx);if(flag=6)for(int i=0;ia;i+)tm.m_pData5*a+i=*(m_pDib.m_lpImage+sizeSaveImage.cx*(i/sizeImage.cx)+i%sizeImage.cx);void CKLchangeView:OnKl() /KL变换计算int i,j;xfea.Init(6,1);xfeat.Init(1,6);temp=CMatrix(6,6);feature=new double6;fea.Init(6,6);double tmave6;/求每幅图象的灰度均值for (i=0;i6;i+)tmavei=0;for(j=0;ja;j+)tmavei+=tm.m_pDatai*a+j;tmavei/=(double)a;/求方差矩阵Rfor(j=0;j36;j+) temp.m_pDataj=0.0;for(i=0;ia;i+) temp.m_pDataj+=(tm.m_pData(j/6)*a+i%6-tmavej/6)*(tm.m_pData(j%6)*a+i%6-tmavej%6);temp.m_pDataj/=(double)a; temp.JacobiEigenv2(feature, fea, 0.000001);/求特征值和特征向量 double t=0.0;double featuretemp6;/临时变量 for(i=0;i6;i+)featuretempi=featurei;/给临时数组赋值 for(i=0;i6;i+)/冒泡法将特征值排序 for(j=0;jfeaturej+1) t=featurej; featurej=featurej+1; featurej+1=t; for(i=0;i6;i+)/求出特征值排序后对应在特征值数组中的位置 for(j=0;j6;j+) if(featurei=featuretempj) indexi=j; void CKLchangeView:OnFirst() /第一主分量 if(NewFileName=) AfxMessageBox(请先读入图象数据!);return;if(fea.m_pData0=0) OnKl();result.Init(1,a);for(int i=0;ia;i+)for(int j=0;j6;j+) result.m_pDatai+=fea.m_pDataindex5+j*6*tm.m_pDataj*a+i; Draw();void CKLchangeView:OnFive() /第五主分量if(NewFileName=) AfxMessageBox(请先读入图象数据!);return;if(fea.m_pData0=0) OnKl();result.Init(1,a);for(int i=0;ia;i+)for(int j=0;j6;j+) result.m_pDatai+=fea.m_pDataindex1+j*6*tm.m_pDataj*a+i; Draw();void CKLchangeView:OnFour() /第四主分量if(NewFileName=) AfxMessageBox(请先读入图象数据!);return;if(fea.m_pData0=0) OnKl();result.Init(1,a);for(int i=0;ia;i+)for(int j=0;j6;j+) result.m_pDatai+=fea.m_pDataindex2+j*6*tm.m_pDataj*a+i; Draw();void CKLchangeView:OnSecond() /第二主分量if(NewFileName=) AfxMessageBox(请先读入图象数据!);return;if(fea.m_pData0=0) OnKl();result.Init(1,a);for(int i=0;ia;i+)for(int j=0;j6;j+) result.m_pDatai+=fea.m_pDataindex4+j*6*tm.m_pDataj*a+i;Draw();void CKLchangeView:OnSix() /第六主分量if(NewFileName=) AfxMessageBox(请先读入图象数据!);return;if(fea.m_pData0=0) OnKl();result.Init(1,a);for(int i=0;ia;i+)for(int j=0;j6;j+) result.m_pDatai+=fea.m_pDataindex0+j*6*tm.m_pDataj*a+i; Draw();void CKLchangeView:OnThird() /第三主分量if(NewFileName=) AfxMessageBox(请先读入图象数据!);return;if(fea.m_pData0=0) OnKl(); result.Init(1,a);for(int i=0;ia;i+)for(int j=0;j6;j+) result.m_pDatai+=fea.m_