主分量变换.txt

主分量变换

bool Zhufenliangchange(CPBDib * m_Dib[2])
{
	CSize size;
    int m_palsize=m_Dib[0]->m_nColorTableEntries*sizeof(RGBQUAD);//条色板大小
    size=m_Dib[0]->GetDimensions();

	int nImageScanWidth=m_Dib[0]->m_ScanWidth;
	int nImageHeight=size.cy;
	int nImageWidth=size.cx;

    int i,ii,jj; 

     int nInputNum=2;
     int nOutputNum=2;
	 double* average;//存放各波段的均值
     average=new double[nInputNum];
	 
	 double tempall;

	 for(i=0;i<nInputNum;i++)
	 {
	     tempall=0.0;

	     for(ii=0;ii<nImageHeight;ii++)
	        for(jj=0;jj<nImageWidth;jj++)
			  tempall+=*(m_Dib[i]->m_lpImage+ii*nImageScanWidth+jj);
		 average[i]=tempall/(nImageHeight*nImageWidth);
	}

// 计算图象的协方差矩阵
    double  **jz; //图象的协方差矩阵
    int i1,j1;

	CMathBase math;
    jz=math.TwoArrayAlloc(nInputNum,nInputNum);

/*以下可能有误
	for(int ii=0;ii<nImageHeight;ii++)
	    for(int jj=0;jj<nImageWidth;jj++)
		{
		    for(j=0;j<nInputNum;j++)  //every band
		        vector[j]=*(m_Dib[j]->m_lpImage+ii*nImageScanWidth+jj)-average[j];//一个像素在各波段的离差
	        for(i1=0;i1<nInputNum;i1++)
		        for(j1=0;j1<nInputNum;j1++)
		           jz[i1][j1]+=vector[i1]*vector[j1];
		 }
*/
//计算协方差矩阵
    for(i1=0;i1<nInputNum;i1++)  //every band
        for(j1=0;j1<nInputNum;j1++)
		{
			jz[i1][j1]=0.0;
			for(int ii=0;ii<nImageHeight;ii++)
        	    for(int jj=0;jj<nImageWidth;jj++)
		           jz[i1][j1]+=(*(m_Dib[i1]->m_lpImage+ii*nImageScanWidth+jj)-average[i1])*
					           (*(m_Dib[j1]->m_lpImage+ii*nImageScanWidth+jj)-average[j1]);
		
			jz[i1][j1]/=nImageHeight*nImageWidth;
		}

		if(average) delete average;

//计算协方差矩阵的特征值和特征向量
	double **resultjz;
	double  *character;

	character=(double*)calloc(nInputNum,sizeof(double));
	
	resultjz=math.TwoArrayAlloc(nInputNum,nInputNum);//申请内存，大小nInputNum×nInputNum
	
/*	if(!math.JACOBI(nInputNum,jz,character,resultjz))//计算出的特征值存于character, 特征向量存于resultjz
	{
		 MessageBox("Failed!","数据错误",MB_OK);
		 return;
	}
	*/
	math.TwoArrayFree(jz);

/////////////////////////////////////////////////////////////////

// 特征值和特征向量从大到小排序
    #define MAX  20
	
	struct
	{
	  double value;
	  int  number;
	}index[MAX],exc;
	
	for (i=0;i<nInputNum;i++)
	{
	  index[i].value=character[i];
	  index[i].number=i;
	}
	
	int a,b;
	
	for(a=0;a<nInputNum;a++)    //特征值从大到小排序
	  for (b=0;b<a;b++)
	  {
		 if (index[b].value<index[b+1].value)
		 {
			exc=index[b];index[b]=index[b+1];index[b+1]=exc;
		 }
	  }


	double **result;

	result=math.TwoArrayAlloc(nInputNum,nInputNum);

    int j;
	
	for(i=0;i<nInputNum;i++)   //按照特征值重排特征向量;结果存于result中
	{
		for(j=0;j<nInputNum;j++)
		{
			 result[j][i]=resultjz[j][index[i].number];
		}
	}
    
	math.TwoArrayFree(resultjz);
    
	free(character);

//将特征向量矩阵result转置构成变换矩阵(K)
	for(i=0;i<nInputNum;i++)
	   for (j=i+1;j<nInputNum;j++)
	   {
		   tempall=result[i][j];result[i][j]=result[j][i];result[j][i]=tempall;
	   }


//进行K__L变换，结果存入文件名以kl开头,扩展名为tmp的m个文件中
//////////////////////////////////////////////////////////////////////////
	struct
	{
	  double max;
	  double min;
	}scale[100];

	int aa;
	unsigned int bb;
	
	int **hist;//直方图，用于进行灰度调整

	hist=new int*[nOutputNum];
	for(aa=0;aa<nOutputNum;aa++)
	{
	   hist[aa]=new int[65535];
	}

	for (aa=0;aa<nOutputNum;aa++)
		for(bb=0;bb<65535;bb++)
				 hist[aa][bb]=0;

//////////////////////////////////////////////////////////////////////////////

	double  last;
	int ci;
	FILE ** fpout;
	
	fpout=new FILE*[nOutputNum];
	char  na[15];

	for(i=0;i<nOutputNum;i++)  //创建中间结果文件
	{
		sprintf(na,"kl%d.tmp",i);
	  /*  if((fpout[i]=fopen(na,"wb"))==NULL)
		{
			MessageBox("文件打开错误 ","文件错误",MB_OK );
			return;
		 }*/
	}


	for(ii=0;ii<nImageHeight;ii++)
	    for(jj=0;jj<nImageWidth;jj++)
		{
	       for(ci=0;ci<nOutputNum;ci++)
		   {
		      last=0.0;
		      for(j=0;j<nInputNum;j++)
	    	     last+=result[ci][j] * *(m_Dib[j]->m_lpImage+ii*nImageScanWidth+jj)/*vect[j]*/;

			  if(ii==0&&jj==0)
			     scale[ci].max=scale[ci].min=last;

		      if(last>scale[ci].max)
			     scale[ci].max=last;
		      if(last<scale[ci].min)
	  	         scale[ci].min=last;
		      
			  hist[ci][(int)last]++;     // 统计直方图，存于数组
		      
			  fwrite(&last,sizeof(double),1,fpout[ci]);
		  }
		}

     math.TwoArrayFree(result);

	for(i=0;i<nOutputNum;i++)
		fclose(fpout[i]);
    
///////////////////////////////////////////////////////////////////////////
//拉伸
//1.处理直方图
/*
	int csr;
	for(aa=0;aa<nOutputNum;aa++) // n幅图象
	{
		csr=scale[aa].min;
		while(hist[aa][csr]+hist[aa][csr+1]+hist[aa][csr+2]+hist[aa][csr+3]+hist[aa][csr+4]<=10)
		   csr++;
		scale[aa].min=csr;
		csr=scale[aa].max;
		while(hist[aa][csr]+hist[aa][csr-1]+hist[aa][csr-2]+hist[aa][csr-3]+hist[aa][csr-4]<=10)
  	 	   csr--;
		scale[aa].max=csr;
	}
*/
//2。拉伸图象
	double sc;
	
	for(i=0;i<nOutputNum;i++)  //打开中间结果文件
	{
		sprintf(na,"kl%d.tmp",i);

	}

    for(aa=0;aa<nOutputNum;aa++)
	{
	    for(int ii=0;ii<nImageHeight;ii++)
	        for(int jj=0;jj<nImageWidth;jj++)
			{
			    fread(&last,sizeof(double),1,fpout[aa]);
			    if(last>scale[aa].max)
					last=scale[aa].max;
			    if(last<scale[aa].min)
					last=scale[aa].min;
			    sc=scale[aa].max-scale[aa].min;
			 
				last=(last-scale[aa].min)*245/sc+10.0;  //拉伸
			 
                *(m_Dib[aa]->m_lpImage+ii*nImageScanWidth+jj)=(unsigned char)last;
				//				fwrite(&bb,sizeof(char),1,fpoutend[aa]);
			}
	}

//	for(aa=0;aa<nOutputNum;aa++)
//		delete hist[aa];
	delete hist;

	for(a=0;a<nOutputNum;a++)  //关闭结果文件
		fclose(fpout[a]);

	for(i=0;i<nOutputNum;i++)  //删除中间结果文件
	{
	   sprintf(na,"kl%d.tmp",i);
	   remove(na);
	}

      char docname[10];
	  for(i=0;i<nOutputNum;i++)
	  {
		  sprintf(docname,"KL%i",i+1);
          //ShowResultImage(m_DibOutput[i],(LPCTSTR)docname);
	  }
      


	return true;
}