高光谱图像ISODATA分类算法

//#include <iostream> #include <string.h> #include "conio.h" #include "math.h" #include "isodata.h" isodata::isodata( ) { n_K=4; //预期聚类中心数目 n_TypeNum=7; //初始分类数目 n_Iteration=3; n_ThetaD=22; n_Threshold=5; //判断收敛阈值 n_ThetaS=30; //同一聚类域中样本标准差的最大值 n_ThetaN=10000; //每个聚类域中最少的样本数 n_L=2; //输入一次可以合并的聚类中心的最多对数 } isodata::~isodata() { for (int j=0;j<n_TypeNum;j++) { delete []CenterDis[j]; } delete []CenterDis; delete []EucliDistance; delete []EachTypePixelSum; delete []CopyEachTypePixelSum; CPLFree(Max); CPLFree(ImageData); CPLFree(ResultImageData); CPLFree(ClusterCenter); delete(EachTypePixelSum); delete(LastEachTypePixelSum); } //计算初始聚类中心 void isodata::gdComputeInitCenter() { int x1=1*nXSize/4,x2=2*nXSize/4,x3=3*nXSize/4; int y1=1*nYSize/4,y2=2*nYSize/4,y3=3*nYSize/4; int A[9][2]={x1,y1,x2,y2,x3,y3,x1,y3,x3,y1,x2,y1,x2,y3,x1,y2,x3,y2}; ClusterCenter=(IDTcentStT *)malloc(sizeof(IDTcentStT)*bandCounts*n_TypeNum); for(int i=0;i<bandCounts;i++) { for(int j=0;j<n_TypeNum;j++) { if(!ROIFlag) { ClusterCenter[i*n_TypeNum+j].d=ImageData[i*nXSize*nYSize+A[j][1]*nXSize+A[j][0]]; ClusterCenter[i*n_TypeNum+j].b=1; } else { ClusterCenter[i*n_TypeNum+j].d=ImageData[i*nXSize*nYSize+A[j][1]*nXSize+A[j][0]]; ClusterCenter[i*n_TypeNum+j].b=1; } } } } void isodata::gdSetColorEntry(GDALRasterBand *poDataBand,int TypeNum) { GDALColorTable *colorTable=new GDALColorTable; GDALColorEntry *colorEntry=new GDALColorEntry; poDataBand->SetColorInterpretation(GCI_PaletteIndex); for(int i=0;i<TypeNum;i++) { colorEntry=&ColorIndexPalatte[i]; colorTable->SetColorEntry(i+1,colorEntry); } poDataBand->SetColorTable(colorTable); } /*--------------------------------------------------------------- 函数说明：根据图像数据，计算每个像元与每个类别中心的欧式距离 输入参数：这个数据块的X、Y大小尺寸，以及图像的像元数据 返回为标记了类别的数据及像元值为1.2.3...N ----------------------------------------------------------------*/ void isodata::gdComputeResult(int nXSize,int nYSize) { double temp; double result; int index; for(int i=0;i<nYSize;i++) { for(int j=0;j<nXSize;j++) { temp=1000000000000000; for(int k=0;k<n_TypeNum;k++) { result=gdComputeEucliDistance(j,i,ImageData,k); if(result<temp) { temp=result; index=k+1;//将像元分类到最小欧氏距离类别 每个循环得出最小距离result } } ResultImageData[i*nXSize+j]=index; } } } //计算每个像元与各聚类中心的距离 //输入有抽象数据集，图像数据 //像元的坐标位置，分类的类别数，以及聚类中心的数据 ///计算每个像元与聚类中心的距离，并返回类别号 double isodata::gdComputeEucliDistance(int x,int y,float* ImageData,int nTypeIndex)//掩膜处理不确定 { double result=0.0; if(!BandMaskFlag) { result=pow((ImageData[y*nXSize+x]-ClusterCenter[nTypeIndex].d),2); for(int i=1;i<bandCounts;i++) { result=result+pow((ImageData[y*nXSize+x+i*nXSize*nYSize]-ClusterCenter[i*n_TypeNum+nTypeIndex].d),2); } } else { int bandNum=nbandMask.bandNum; result=pow((ImageData[y*nXSize+x]-ClusterCenter[nTypeIndex].d),2); for(int i=1;i<bandNum;i++) { result=result+pow((ImageData[y*nXSize+x+i*nXSize*nYSize]-ClusterCenter[i*n_TypeNum+nTypeIndex].d),2); } } return result; } /*根据分类的数据将样本个数少于ThetaN的类划入其他类*/ bool isodata::gdInsert() { int b=0; int CopyTypeNum=n_TypeNum; for (int j=0;j<n_TypeNum;j++) { int temp; if(EachTypePixelSum[j]<n_ThetaN) { if (j==n_TypeNum-1) //若判断是第5类小于样本个数则把它归并到第1类 { temp=EachTypePixelSum[0]; IDTclassStT *copy=new IDTclassStT[temp]; memcpy(copy,Type[0],sizeof(IDTclassStT)*temp); //内存临时存放变量 EachTypePixelSum[0]+=EachTypePixelSum[j]; delete []Type[0]; Type[0]=new IDTclassStT[EachTypePixelSum[0]]; //重新分配大内存 memcpy(Type[0],copy,sizeof(IDTclassStT)*temp); //将临时内存放回 for (int m1=temp,m2=0;m1<EachTypePixelSum[0],m2<EachTypePixelSum[j];m1++,m2++) { Type[0][m1].data.x=Type[j][m2].data.x; Type[0][m1].data.y=Type[j][m2].data.y; ResultImageData[Type[j][m2].data.y*nXSize+Type[j][m2].data.x]=1; } delete []copy; delete []Type[j]; } else { temp=EachTypePixelSum[j+1]; IDTclassStT *copy=new IDTclassStT[temp]; memcpy(copy,Type[j+1],sizeof(IDTclassStT)*temp);//内存临时存放变量 EachTypePixelSum[j+1]+=EachTypePixelSum[j]; delete []Type[j+1]; Type[j+1]=new IDTclassStT[EachTypePixelSum[j+1]]; //重新分配大内存 memcpy(Type[j+1],copy,sizeof(IDTclassStT)*temp);//将临时内存放回 for (int m1=temp,m2=0;m1<EachTypePixelSum[j+1],m2<EachTypePixelSum[j];m1++,m2++) { Type[j+1][m1].data.x=Type[j][m2].data.x; Type[j+1][m1].data.y=Type[j][m2].data.y; ResultImageData[Type[j][m2].data.y*nXSize+Type[j][m2].data.x]=j+1; } delete []copy; /* delete []Type[j];*/ for (int z=j+1;z<n_TypeNum;z++) { Type[z]=Type[z-1]; EachTypePixelSum[z]=EachTypePixelSum[z-1]; for (int t=0;t<bandCounts;t++) { ClusterCenter[t*n_TypeNum+j]=ClusterCenter[t*n_TypeNum+j-1]; } } delete []Type[n_TypeNum-1]; } n_TypeNum--; } } if (n_TypeNum<CopyTypeNum) { b=1; } return b; } //计算类各分量的标准差及最大分量标准差 void isodata::gdComputeQji() { double **Q=new double*[n_TypeNum]; Max=new double[n_TypeNum]; for (int i=0;i<n_TypeNum;i++) { Q[i]=new double[bandCounts]; Max[i]=0.0; } for (int j=0;j<n_TypeNum;j++) { for (i=0;i<bandCounts;i++) { Q[j][i]=0.0; } } for (i=0;i<bandCounts;i++) { for (int j=0;j<n_TypeNum;j++) { for (int m=0;m<EachTypePixelSum[j];m++) { Q[j][i]+=pow(ImageData[i*nXSize*nYSize+Type[j][m].data.y*nXSize+Type[j][m].data.x]-ClusterCenter[i*n_TypeNum+j].d,2); } Q[j][i]/=EachTypePixelSum[j]; Q[j][i]=sqrt(Q[j][i]); if (Max[j]<Q[j][i])// { Max[j]=Q[j][i]; } } } for(i=0;i<n_TypeNum;i++) { delete []Q[i]; } delete []Q; } //计算新的聚类中心 //输入有抽象数据集，图像数据 //上一次分类的结果数据，类别数 //还有上一次的聚类中心数据 void isodata::gdComputeNewClusterCenter() { //每个类别的像素值总和 double* Sum=new double[n_TypeNum]; for(int i=0;i<bandCounts;i++) { for(int j=0;j<n_TypeNum;j++) { Sum[j]=0; } for(int y=0;y<nYSize;y++) { for(int x=0;x<nXSize;x++) { Sum[ResultImageData[y*nXSize+x]-1]=Sum[ResultImageData[y*nXSize+x]-1]+ImageData[i*nXSize*nYSize+y*nXSize+x]; } } for(int z=0;z<n_TypeNum;z++) { ClusterCenter[i*n_TypeNum+z].d=Sum[z]/EachTypePixelSum[z]; } } delete []Sum; } //计算所有样本偏离均值的平均距离 void isodata::gdComputeALLDAV() { ALLDAV=0.0; EucliDistance=new double[n_TypeNum]; for (int i=0;i<n_TypeNum;i++) { EucliDistance[i]=0; } //类样本偏移均值的平均距离EucliDistance[j] for (i=0;i<n_TypeNum;i++) { for (int j=0;j<EachTypePixelSum[i];j++) { for (int t=0;t<bandCounts;t++) { EucliDistance[i]+=pow(ImageData[Type[i][j].data.y*nXSize+Type[i][j].data.x+t*nXSize*nYSize]-ClusterCenter[t*n_TypeNum+i].d,2); } } } for(int j=0;j<n_TypeNum;j++) { ALLDAV+=EachTypePixelSum[j]*EucliDistance[j]; } ALLDAV/=nXSize*nYSize; } //初始化每个类别的像元个数 void isodata::gdInitLastEachTypePixelSum() { LastEachTypePixelSum=new int[n_TypeNum]; for(int i=0;i<n_TypeNum;i++) { LastEachTypePixelSum[i]=1; } } //根据分类的数据计算各个类别像元个数 void isodata::gdComputeEachTypePixelSum(int nXSize,int nYSize) { int m=0; int t=0; EachT