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

1、模式识别作业三一、编程要求：编程实现KL变换，并对TM六波段图像进行演算。KL变换的思想是：从n维特征选取m维特征，删去的n-m维特征不一定就是无用的信息，如何在信息损失最小的情况下选取特征，在理论上就显得更严密些。通常采用正交变换，得到新的经变换的模式， 以保证信息损失最小情况下获得有利于分类的特征。二、编程思想：将6副图象依次输入获得灰度值存在一个6*size(size为一副图象的像素数)的二维数组中，计算每个波段的灰度均值，然后计算协方差矩阵，得出特征值和特征向量，再重新计算要输出的处理的图象的灰度。三、核心代码：void CKLchangeView:Oninput() /TM图象数据读

2、入AfxMessageBox("注意:请一次将6副图象都读入!");/读图象的方法int flag=0;while(flag<6)CFileDialog dlg(TRUE,"bmp",".bmp");if(dlg.DoModal()=IDOK)m_NewFileNameflag=dlg.GetFileName();for(int j=0;j<5;j+)if(flag=j|m_NewFileNameflag!=m_NewFileNamej) continue;else AfxMessageBox("已打开过的图象!

3、请按方法重新打开!");flag=0; for(int i=0;i<6;i+) m_NewFileNamei=""return;NewFileName=m_NewFileNameflag;flag+;/判断是否是已经打开的图象，若已经打开，则要重新按方法打开else /不读图象，或是只是要看图象for(int i=0;i<6;i+) m_NewFileNamei=""return;CDC *pDC = GetWindowDC();CFile f(NewFileName,0);m_pDib.Read(&f);/读图象if(m_

4、pDib.m_nColorTableEntries!=256)/判断是否是灰度图象AfxMessageBox("请打开灰度图象!");return;CKLchangeDoc* pDoc = GetDocument(); pDoc->pDib=&m_pDib;OnDraw(pDC);/调用OnDraw函数绘图sizeImage = m_pDib.GetDimensions(); sizeSaveImage = m_pDib.GetDibSaveDim(); a=sizeSaveImage.cx * sizeSaveImage.cy;/象素数if(flag=1)

5、tm=CMatrix(6,a);/开第一副图象时为矩阵分配内存 /将每幅图象象素灰度值放入矩阵中if(flag=1)for(int i=0;i<a;i+)tm.m_pDatai=*(m_pDib.m_lpImage+sizeSaveImage.cx*(i/sizeImage.cx)+i%sizeImage.cx);if(flag=2)for(int i=0;i<a;i+)tm.m_pDataa+i=*(m_pDib.m_lpImage+sizeSaveImage.cx*(i/sizeImage.cx)+i%sizeImage.cx);if(flag=3)for(int i=0;i&

6、lt;a;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;i<a;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;i<a;i+)tm.m_pData4*a+i=*(m_pDib.m_lpImage+sizeSaveImage.cx*(i/siz

7、eImage.cx)+i%sizeImage.cx);if(flag=6)for(int i=0;i<a;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=

8、0;i<6;i+)tmavei=0;for(j=0;j<a;j+)tmavei+=tm.m_pDatai*a+j;tmavei/=(double)a;/求方差矩阵Rfor(j=0;j<36;j+) temp.m_pDataj=0.0;for(i=0;i<a;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);/求特征值

9、和特征向量 double t=0.0;double featuretemp6;/临时变量 for(i=0;i<6;i+)featuretempi=featurei;/给临时数组赋值 for(i=0;i<6;i+)/"冒泡法"将特征值排序 for(j=0;j<5;j+) if(featurej>featurej+1) t=featurej; featurej=featurej+1; featurej+1=t; for(i=0;i<6;i+)/求出特征值排序后对应在特征值数组中的位置 for(j=0;j<6;j+) if(featurei=f

10、eaturetempj) indexi=j; void CKLchangeView:OnFirst() /第一主分量 if(NewFileName="") AfxMessageBox("请先读入图象数据!");return;if(fea.m_pData0=0) OnKl();result.Init(1,a);for(int i=0;i<a;i+)for(int j=0;j<6;j+) result.m_pDatai+=fea.m_pDataindex5+j*6*tm.m_pDataj*a+i; Draw();void CKLchangeVie

11、w:OnFive() /第五主分量if(NewFileName="") AfxMessageBox("请先读入图象数据!");return;if(fea.m_pData0=0) OnKl();result.Init(1,a);for(int i=0;i<a;i+)for(int j=0;j<6;j+) result.m_pDatai+=fea.m_pDataindex1+j*6*tm.m_pDataj*a+i; Draw();void CKLchangeView:OnFour() /第四主分量if(NewFileName="&quo

12、t;) AfxMessageBox("请先读入图象数据!");return;if(fea.m_pData0=0) OnKl();result.Init(1,a);for(int i=0;i<a;i+)for(int j=0;j<6;j+) result.m_pDatai+=fea.m_pDataindex2+j*6*tm.m_pDataj*a+i; Draw();void CKLchangeView:OnSecond() /第二主分量if(NewFileName="") AfxMessageBox("请先读入图象数据!")

13、;return;if(fea.m_pData0=0) OnKl();result.Init(1,a);for(int i=0;i<a;i+)for(int j=0;j<6;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.Ini

14、t(1,a);for(int i=0;i<a;i+)for(int j=0;j<6;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;i<a;i+)for(int j=0;j<6;j+) result.m_pDatai+=fea.m_pDataindex3+j*6*tm.m_pDataj*a+i; Draw();void CKLchangeView:Draw() /图象输出int i;for(i=0;i<a;i+)/处理越界的灰度值if(result.m_pDatai<0) result.m_pDatai=0; if(result.m_pDatai>255) result.m_pDatai=255; for (i = 0; i <a; i +)/重新给图象每个象素赋灰度值 un