rsipview.cpp

基于小波的SAR斑点处理

		int MaxClassNum = characterWinClsDlg.m_nClassNum;
		float DeltaRange = (float)characterWinClsDlg.m_nClusterTime/100;
		int FileNum = characterWinClsDlg.m_nImageCount;
		int HasClassedNum = 0;
		
		DWORD  *ClassPoint = new DWORD[MaxClassNum+1];
		double *ClassDim = new double[(MaxClassNum+1)*FileNum];
		double *TotalClassDim = new double[(MaxClassNum+1)*FileNum];
		double *CoMatrix = new double[(MaxClassNum+1)*FileNum*FileNum];
		double *UpLimit = new double[FileNum*MaxClassNum];
		double *DownLimit = new double[FileNum*MaxClassNum];
		//Alloc the multi-Memory
		int BmpCol=(Width+3)/4*4;
		DWORD size = (DWORD)BmpCol*Height;

		BYTE *pNewData = (BYTE *)malloc(size);

		BOOL bInRange;
		//Initialize
		for(int i=0; i<MaxClassNum+1; i++)
		{
			ClassPoint[i] = 0L;
			for(int j=0; j<FileNum; j++)
				TotalClassDim[i*FileNum+j] = 0.0;
			for(j=0; j<FileNum*FileNum; j++)
				CoMatrix[i*FileNum*FileNum+j] = 0.0;
		}

		BYTE *Block = new BYTE[Width*FileNum];
		
		WORD RPro;
		WORD Process = 0;
		CString *pszFileName;
		CFile File;
		
		//Begin Classification
		BeginWaitCursor();
		for(i=0; i<Height; i++)
		{
			RPro=(int)(100.0 * (i+1) / Height);
			if(RPro>Process)
			{
				for(int j=0; j<RPro-Process; j++)
					UpdateStatusBar();
				Process=RPro;
			}
			//Read File
			for(int j=0; j<FileNum; j++)
			{
				pszFileName = (CString *)(characterWinClsDlg.m_aImageArray[j]);
				File.Open(*pszFileName,CFile::modeRead);
				
				File.Seek((DWORD)i*Width,CFile::begin);
				File.Read(Block+(DWORD)j*Width,(DWORD)Width);
				File.Close();
			}

			for(j=0; j<Width; j++)
			{
				for(int k=0; k<HasClassedNum; k++)
				{
					bInRange = TRUE;
					for(int l=0; l<FileNum; l++)
						if(Block[l*Width+j]<DownLimit[k*FileNum+l]||
							Block[l*Width+j]>UpLimit[k*FileNum+l])
							{
								bInRange = FALSE;
								break;
							}
					if(bInRange)  //Has searched the perfect class
						break;
				}
				if(k<HasClassedNum)
				{
					ClassPoint[k]++;
					for(int l=0; l<FileNum; l++)
						{
							TotalClassDim[k*FileNum+l] += Block[l*Width+j];
							ClassDim[k*FileNum+l] = TotalClassDim[k*FileNum+l]/ClassPoint[k];
							DownLimit[k*FileNum+l] = ClassDim[k*FileNum+l]*(1-DeltaRange/2);
							UpLimit[k*FileNum+l] = ClassDim[k*FileNum+l]*(1+DeltaRange/2);
							pNewData[(DWORD)(Height-i-1)*BmpCol+j] = k;
							for(int n=0;n<FileNum;n++)
								CoMatrix[k*FileNum*FileNum+l*FileNum+n]+=
									(double)Block[l*Width+j]*Block[n*Width+j];
						}
				}
				else if(k==HasClassedNum && HasClassedNum<MaxClassNum)
				{
					HasClassedNum++;
					ClassPoint[k]=1;
					for(int l=0; l<FileNum; l++)
						{
							TotalClassDim[k*FileNum+l] = Block[l*Width+j];
							ClassDim[k*FileNum+l] = Block[l*Width+j];
							DownLimit[k*FileNum+l] = ClassDim[k*FileNum+l]*(1-DeltaRange/2);
							UpLimit[k*FileNum+l] = ClassDim[k*FileNum+l]*(1+DeltaRange/2);
							pNewData[(DWORD)(Height-i-1)*BmpCol+j] = k;
							for(int n=0; n<FileNum; n++)
								CoMatrix[k*FileNum*FileNum+l*FileNum+n]+=
									(double)Block[l*Width+j]*Block[n*Width+j];
						}
				}
				else if(HasClassedNum == MaxClassNum)
				{
					//assert(HasClassedNum==MaxClassNum+1);
					ClassPoint[MaxClassNum]++;
					for(int l=0; l<FileNum; l++)
						{
							TotalClassDim[MaxClassNum*FileNum+l] += Block[l*Width+j];
							ClassDim[MaxClassNum*FileNum+l] = TotalClassDim[k*FileNum+l]/ClassPoint[k];
							pNewData[(DWORD)(Height-i-1)*BmpCol+j] = MaxClassNum;
							for(int n=0;n<FileNum;n++)
								CoMatrix[MaxClassNum*FileNum*FileNum+l*FileNum+n]+=
									(double)Block[l*Width+j]*Block[n*Width+j];
						}
				}
			}//j from 0 to width ended
		}//Class ended;
		delete []Block;
		delete []UpLimit;
		delete []DownLimit;

		MaxClassNum++;

		CRSIPApp *theApp = (CRSIPApp *)AfxGetApp();
		CRSImage *pRSImage = new CRSImage();
		pRSImage->m_pPal = new COLORREF[256];
		for(i=0; i<256; i++)
		{
			pRSImage->m_pPal[i] = theApp->m_pSystemPal[i];
		}
		pRSImage->m_nImgType = IMAGE_THEMATIC;
		pRSImage->m_nImgHeight = Height;
		pRSImage->m_nImgWidth = Width;
		pRSImage->m_pImgData = pNewData;

		pRSImage->m_nClassNumber = MaxClassNum;
		CRsThematicInfo *pThemeInfo;
		CString szNoNameClass;
		for(i=0; i<MaxClassNum; i++)
		{
			pThemeInfo = new CRsThematicInfo();
			pThemeInfo->m_nValue = i;		//1--255
			szNoNameClass.LoadString(IDS_NONAMECLASS);
			pThemeInfo->m_szClassName.Format("%s %d",szNoNameClass,i);
			pThemeInfo->m_clrColor = theApp->m_pSystemPal[i];
			pThemeInfo->m_pdCenter = new double[FileNum];
			for(int j=0; j<FileNum; j++)
			{
				pThemeInfo->m_pdCenter[j] = ClassDim[i*FileNum+j];
			}
			pThemeInfo->m_pdCovariance = new double[FileNum*FileNum];
			for(j=0; j<FileNum; j++)
				for(int k=0; k<FileNum; k++)
				{
					if(ClassPoint[i]==0L)
						pThemeInfo->m_pdCovariance[j*FileNum+k]=0.0;
					else
						pThemeInfo->m_pdCovariance[j*FileNum+k]=
								CoMatrix[i*FileNum*FileNum+j*FileNum+k]/
								ClassPoint[i]-
								ClassDim[i*FileNum+j]*ClassDim[i*FileNum+k];
				}
			pRSImage->m_aThemeInfoArray.Add(pThemeInfo);
		}
		pRSImage->m_nLinkedImage = FileNum;
		CString *pszTemp;
		for(i=0; i<FileNum; i++)
		{
			pszTemp = new CString();
			*pszTemp = *(CString *)(characterWinClsDlg.m_aImageArray[i]);
			pRSImage->m_aLinkedImageNameArray.Add(pszTemp);
		}

		pRSImage->m_pdwHistogram = new DWORD[256];
		for(i=0; i<256; i++)
			pRSImage->m_pdwHistogram[i] = 0;
		for(i=0; i<MaxClassNum; i++)
			pRSImage->m_pdwHistogram[i] = ClassPoint[i];

		delete []ClassDim;
		delete []CoMatrix;
		delete []ClassPoint;
		delete []TotalClassDim;
		
		theApp->CreateNewImageDoc(pRSImage,"特征窗口曲线法分类结果");
		EndWaitCursor();
	}	
}

void CRSIPView::OnClassFisherrule() 
{
	// 读出分类信息
	CFileDialog fileDlg(TRUE,"*.CLS",NULL,
		OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT,
		"Classes Files(*.cls)|*.cls||",this);
	
	if(fileDlg.DoModal()==IDOK)
	{
		CFile file;
		if(0 == file.Open(fileDlg.GetPathName(),CFile::modeRead))
		{
			AfxMessageBox("打开影象分类信息文件失败!");
			return;
		}
		else
		{
			// 读影象分类信息文件
			CRSClassifyInfoFile classifyInfoFile;
			classifyInfoFile.ReadFromFile(file);
			file.Close();

			int ClassNum = classifyInfoFile.m_nClassCount;
			int FileNum = classifyInfoFile.m_nImageCount;
			int FileNum2 = FileNum*FileNum;

			double *ClassDim = new double[ClassNum*FileNum];
			double *CoMatrix = new double[ClassNum*FileNum2];

			CRSClassInfo *pClassInfo;
			for(int i=0; i<ClassNum; i++)
			{
				pClassInfo = (CRSClassInfo *)(classifyInfoFile.m_aClassInfoArray[i]);
				for(int j=0; j<FileNum; j++)
					ClassDim[i*FileNum+j] = pClassInfo->m_pdCenter[j];
				for(j=0; j<FileNum2; j++)
					CoMatrix[i*FileNum2+j] = pClassInfo->m_pdCovariance[j];
			}
	
			double *MatrixE = new double[FileNum2];
			double *InvMatrixE = new double[FileNum2];
			double *MatrixB = new double[FileNum2];
			double *MatrixQ = new double[ClassNum*ClassNum];
			double *TClassDim = new double[ClassNum*FileNum];
			double *ResultMat = new double[FileNum2];
			double *FResultMat = new double[FileNum2];
			double *TempMatrix = new double[FileNum*ClassNum];
	
			double *NewCoMatrix = new double[ClassNum*FileNum2];
			double *NewCenter = new double[ClassNum*FileNum];
			DWORD *ClassPoint= new DWORD[ClassNum];
	
			for(i=0; i<ClassNum; i++)
			{
				ClassPoint[i] = 0L;
				for(int j=0; j<FileNum; j++)
					NewCenter[i*FileNum+j] = 0.0;
				for(j=0; j<FileNum2; j++)
					NewCoMatrix[i*FileNum2+j] = 0.0;
			}	
			//Get Matrix E
			for(i=0; i<FileNum2; i++)
				MatrixE[i]=0.0;
			for(i=0; i<ClassNum; i++)
				for(int j=0; j<FileNum2; j++)
					MatrixE[j] += CoMatrix[i*FileNum2+j];
	
			//Get Matrix Q
			for(i=0; i<ClassNum; i++)
				for(int j=0; j<ClassNum; j++)
				{
					if(i==j)
						MatrixQ[i*ClassNum+j]=(double)(ClassNum-1)/ClassNum;
					else
						MatrixQ[i*ClassNum+j]=(double)1/(-ClassNum);
				}
		
			//Trans ClassDim
			trn(ClassDim,TClassDim,ClassNum,FileNum);
			//Mul Matrix,Get B
			mul(TClassDim,MatrixQ,TempMatrix,FileNum,ClassNum,ClassNum);
			mul(TempMatrix,ClassDim,MatrixB,FileNum,ClassNum,FileNum);
	
			//Inv E
			inv(MatrixE,InvMatrixE,FileNum);
			//mul,Get Result Matrix
			mul(InvMatrixE,MatrixB,ResultMat,FileNum,FileNum,FileNum);
	
			delete []MatrixE;
			delete []InvMatrixE;
			delete []MatrixB;
			delete []MatrixQ;
			delete []TClassDim;
			delete []TempMatrix;
	
			double *Lamda = new double[FileNum];
			Yakbi(ResultMat,FResultMat,Lamda,FileNum);
	
			int *m=new int[FileNum];
			double *ResultVec=new double[FileNum];
	
			sortd(Lamda,FileNum,m);
			for(i=0; i<FileNum; i++)
				ResultVec[i]=FResultMat[m[0]*FileNum+i];

			delete []ResultMat;
			delete []FResultMat;
			delete []Lamda;
			delete []m;

			double *SamCheck = new double[ClassNum];
			for(i=0; i<ClassNum; i++)
				mul(ResultVec,ClassDim+i*FileNum,SamCheck+i,1,FileNum,1);
	
			int Height = classifyInfoFile.m_nImageHeight;
			int Width = classifyInfoFile.m_nImageWidth;
	
			BeginWaitCursor();
			//Alloc the multi-Memory
			int BmpCol=(Width+3)/4*4;
			DWORD size = (DWORD)BmpCol*Height;
	
			BYTE *pNewData = (BYTE *)malloc(size);
	
			if(pNewData==NULL)
			{
				MessageBox("没有足够内存进行分类!",
					NULL,MB_OK|MB_ICONQUESTION);
				return;
			}

			BYTE *Block = new BYTE[Width*FileNum];
	
			CFile File;
			double *X_M=new double[FileNum];
			float Min;
			BYTE ThisClass;
			double Result;
			double TResult;
	
			int Process=0;
			int RPro;
			CString *pszFileName;
		
			for(i=0; i<Height; i++)
			{
				RPro=(int)(100.0 * (i+1) / Height);
				if(RPro>Process)
				{
					for(int j=0;j<RPro-Process;j++)
						UpdateStatusBar();
					Process=RPro;
				}
				//Read File
				for(int j=0; j<FileNum; j++)
				{
					pszFileName = (CString *)(classifyInfoFile.m_aImageNameArray[j]);
					File.Open(*pszFileName,CFile::modeRead);	

					File.Seek((DWORD)i*Width,CFile::begin);
					File.Read((Block+(DWORD)j*Width),(DWORD)Width);
					File.Close();
				}	
				for(j=0; j<Width; j++)
				{
					Min = MAXFLOAT;
					for(int k=0; k<FileNum; k++)
						X_M[k] = Block[k*Width+j];
					mul(ResultVec,X_M,&Result,1,FileNum,1);
					//Search The Min Function Value
					for(int l=0; l<ClassNum; l++)
					{
						TResult=fabs(Result-SamCheck[l]);
						if(Min>TResult)
						{
							Min=float(TResult);
							ThisClass=l;
						}
					}
					*(pNewData+(DWORD)(Height-i-1)*BmpCol+j)=ThisClass;
					ClassPoint[ThisClass]++;
					for(l=0; l<FileNum; l++)
					{
						NewCenter[ThisClass*FileNum+l]+=
							(double)Block[l*Width+j];
						for(int k=0;k<FileNum;k++)
							NewCoMatrix[ThisClass*FileNum2+l*FileNum+k]+=
									(double)Block[l*Width+j]*
									(double)Block[k*Width+j];
					}
				}
			}
			delete []X_M;
			delete []SamCheck;
			delete []ResultVec;
			delete []ClassDim;
			delete []CoMatrix;
			delete []Block;
	
			//Create new document
			CRSIPApp *theApp = (CRSIPApp *)AfxGetApp();
			CRSImage *pRSImage = new CRSImage();
	
			pRSImage->m_nImgType = IMAGE_THEMATIC;
			pRSImage->m_nImgHeight = Height;
			pRSImage->m_nImgWidth = Width;
			pRSImage->m_pImgData = pNewData;
	
			pRSImage->m_nClassNumber = ClassNum;
			CRsThematicInfo *pThemeInfo;
			CString szNoNameClass;
			for(i=0; i<ClassNum; i++)
			{
				pThemeInfo = new CRsThematicInfo();
	
				pClassInfo = (CRSClassInfo *)(classifyInfoFile.m_aClassInfoArray[i]);
	
				pThemeInfo->m_nValue = pClassInfo->m_nValue;
				pThemeInfo->m_szClassName = pClassInfo->m_szClassName;
				pThemeInfo->m_clrColor = pClassInfo->m_clrColor;
				pThemeInfo->m_pdCenter = new double[FileNum];
				for(int j=0; j<FileNum; j++)
				{
					if(ClassPoint[i]==0)
						pThemeInfo->m_pdCenter[j] = 0.0;
					else
						pThemeInfo->m_pdCenter[j] = NewCenter[i*FileNum+j]/ClassPoint[i];
				}
				pThemeInfo->m_pdCovariance = new double[FileNum2];
				for(j=0; j<FileNum; j++)
					for(int k=0; k<FileNum; k++)
					{
						if(ClassPoint[i]==0L)
							pThemeInfo->m_pdCovariance[j*FileNum+k]=0.0;
						else
							pThemeInfo->m_pdCovariance[j*FileNum+k]=
								NewCoMatrix[i*FileNum2+j*FileNum+k]/ClassPoint[i]-
								ClassDim[i*FileNum+j]/ClassPoint[i]*ClassDim[i*FileNum+k]/ClassPoint[i];
					}
				pRSImage->m_aThemeInfoArray.Add(pThemeInfo);
			}
			pRSImage->m_nLinkedImage = FileNum;
			CString *pszTemp;
			for(i=0;i<FileNum;i++)
			{
				pszTemp = new CString();
				*pszTemp = *(CString *)(classifyInfoFile.m_aImageNameArray[i]);
				pRSImage->m_aLinkedImageNameArray.Add(pszTemp);
			}
	
			pRSImage->m_pdwHistogram = new DWORD[256];
			for(i=0; i<256; i++)
				pRSImage->m_pdwHistogram[i] = 0;
			for(i=0; i<ClassNum; i++)
				pRSImage->m_pdwHistogram[i] = ClassPoint[i];
	
			delete []NewCenter;
			delete []NewCoMatrix;
			delete []ClassPoint;
			
			theApp->CreateNewImageDoc(pRSImage,"费歇准则分类结果");

			EndWaitCursor();
		}
	}
}

void CRSIPView::OnClassIsodistance() 
{
	CEqualMixedDistDlg equalMixedDistDlg(this);
	if(equalMixedDistDlg.DoModal() == IDOK)
	{
		struct Class
		{
			int ClassOrder;
			double *ClassCenter;
			double