部分源码选自GDAL库的官方网址：www.gdal.org，其余的代码为笔者自己编写。
// readfile.cpp : 定义控制台应用程序的入口点。
//
/*
part of the codes were cite from: http://www.gdal.org/gdal_tutorial.html
and remaining of code were created :by www.cnblogs.com/AmatVictorialCuram/
and please mark the site if you cite the program fo commercial of publication.
*/
#include "stdafx.h"
#include <iostream>
#include "gdal.h"
#include "gdal_priv.h"
/**/
#include <iomanip>
#include <fstream>
using namespace std;
/**/
int minLabel(int a,int b); int maxLabel(int a,int b); double csd(int *Parea,int length); int main()
{
/*
part of the codes were cite from: http://www.gdal.org/gdal_tutorial.html
and remaining of code were created :by www.cnblogs.com/AmatVictorialCuram/
and please mark the site if you cite the program fo commercial of publication.
*/
const char *pszFilename="E:\\tif/fragstats/sample.tif";
GDALDataset *poDataset; GDALAllRegister(); poDataset = (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );
if( poDataset == NULL )
{
cout<<"file read error!"<<endl;;
}
double adfGeoTransform[]; printf( "Driver: %s/%s\n",
poDataset->GetDriver()->GetDescription(),
poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) ); printf( "Size is %dx%dx%d\n",
poDataset->GetRasterXSize(), poDataset->GetRasterYSize(),
poDataset->GetRasterCount() ); if( poDataset->GetProjectionRef() != NULL )
printf( "Projection is `%s'\n", poDataset->GetProjectionRef() ); if( poDataset->GetGeoTransform( adfGeoTransform ) == CE_None )
{
printf( "Origin = (%.6f,%.6f)\n",
adfGeoTransform[], adfGeoTransform[] ); printf( "Pixel Size = (%.6f,%.6f)\n",
adfGeoTransform[], adfGeoTransform[] );
} GDALRasterBand *poBand;
int nBlockXSize, nBlockYSize;
int bGotMin, bGotMax;
double adfMinMax[]; poBand = poDataset->GetRasterBand( );
poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
printf( "Block=%dx%d Type=%s, ColorInterp=%s\n",
nBlockXSize, nBlockYSize,
GDALGetDataTypeName(poBand->GetRasterDataType()),
GDALGetColorInterpretationName(
poBand->GetColorInterpretation()) ); adfMinMax[] = poBand->GetMinimum( &bGotMin );
adfMinMax[] = poBand->GetMaximum( &bGotMax );
if( ! (bGotMin && bGotMax) )
GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax); printf( "Min=%.3fd, Max=%.3f\n", adfMinMax[], adfMinMax[] ); if( poBand->GetOverviewCount() > )
printf( "Band has %d overviews.\n", poBand->GetOverviewCount() ); if( poBand->GetColorTable() != NULL )
printf( "Band has a color table with %d entries.\n",
poBand->GetColorTable()->GetColorEntryCount() );
/*****
Reading Raster Data
*****/
float *pafScanline;
int nXSize = poBand->GetXSize();
int nYSize=poBand->GetYSize(); //读取图像的nXSize*nYSize数据
pafScanline = (float*) CPLMalloc(sizeof(float)*nXSize*nYSize);//创建指针
poBand->RasterIO(
GF_Read,//第一个参数表示要读入数据还是写入数据
, ,//nXOff, nYOff表示读取或者写入图像数据的起始坐标图像的左上角坐标为（0,0）
nXSize, nYSize,/*nXSize, nYSize表示读取或者写入图像数据的窗口大小,nXSize表示宽度,
nYSize表示高度,均使用像素为单位,该宽度和高度是从第二个和第三个参数处开始计算。
这两个参数和第二第三个参数一起表示就是,读取和写入图像的窗口位置和大小。*/
pafScanline, //指向存储数据的一个指针
nXSize, nYSize,//指定缓冲区的大小
GDT_Float32,
, );
// poBand->RasterIO( GF_Read, 0, 0, nXSize, nYSize, pafScanline, nXSize, nYSize, GDT_Float32, 0, 0 );
cout<<"列数：nXSize="<<nXSize<<endl;
cout<<"行数：nYSize="<<nYSize<<endl;
//system("pause"); //创建nXSize*nYSize的float数组
float **data=new float *[nXSize];//每列有nXSize列,data数组的大小与dataLabel数组的大小相同,类型不同
int **dataLabel=new int *[nXSize];//创建标签数组,有nYSize行,nXSize列
int t=,i=,j=;
for(int i=;i<nYSize;i++)//nYSize行
{
dataLabel[i]=new int [nXSize];
data[i] = new float [nXSize];
} cout<<"输出元数据数组data："<<endl;
for(int i=;i<nYSize;i++)
{
for(int j=;j<nXSize;j++)
{
data[i][j]=pafScanline[t++];
cout<<setw()<<data[i][j];
}
cout<<endl;
}
cout<<endl;
//system("pause");
cout<<"输出标签数组dataLabel数组："<<endl;
for(int i=;i<nYSize;i++)
{
cout<<endl;
for(int j=;j<nXSize;j++)
{
dataLabel[i][j]=nYSize*nXSize;
cout<<setw()<<dataLabel[i][j];
}
}
//system("pause");
cout<<endl;
cout<<"联通区域标记算法开始："<<endl;
dataLabel[][]=;//把左上角的第一个标签赋值为1
int indexcolor=;
for(i=;i<nYSize;i++)
{
for(j=;j<nXSize;j++)
{
if(i==)
{
if(j==)
continue;
if(data[i][j]==data[i][j-])
dataLabel[i][j]=dataLabel[i][j-];
else
{
dataLabel[i][j]=++indexcolor;
}
continue;
}
if(j==)
{
if(data[i][j]==data[i-][j])
{
if(dataLabel[i-][j]<dataLabel[i][j])
dataLabel[i][j]=dataLabel[i-][j];
}
if(data[i][j]==data[i-][j+])
{
if(dataLabel[i-][j+]<dataLabel[i][j])
dataLabel[i][j]=dataLabel[i-][j+];
}
if(dataLabel[i][j]==nYSize*nXSize)
{
dataLabel[i][j]=++indexcolor;
}
continue;
}
//if(edgeCheckL(i,j) && data[i][j]==data[i][j-1])//左边的,
if(j> && data[i][j]==data[i][j-])
{
if(dataLabel[i][j-]<dataLabel[i][j])
dataLabel[i][j]=dataLabel[i][j-]; }
//if(edgeCheckLU(i,j) && data[i][j]==data[i-1][j-1])//左上角的
if((i> && j>)&& data[i][j]==data[i-][j-])
{
if(dataLabel[i-][j-]<dataLabel[i][j])
dataLabel[i][j]=dataLabel[i-][j-];
}
//if(edgeCheckU(i,j) && data[i][j]==data[i-1][j])//上面的
if(i> && data[i][j]==data[i-][j])
{
if(dataLabel[i-][j]<dataLabel[i][j])
dataLabel[i][j]=dataLabel[i-][j];
}
//if(edgeCheckUR(i,j) && data[i][j]==data[i-1][j+1])//右上角的
if((i> && j<nYSize-) && data[i][j]==data[i-][j+])
{
if(dataLabel[i-][j+]<dataLabel[i][j])
dataLabel[i][j]=dataLabel[i-][j+];
}
if(dataLabel[i][j]==nYSize*nXSize)
dataLabel[i][j]=++indexcolor;
}
}
//system("pause");
cout<<endl;
cout<<"输出第一次标记数组"<<endl<<endl;
for(i=;i<nYSize;i++)
{
for(j=;j<nXSize;j++)
{
cout<<setw()<<dataLabel[i][j];
}
cout<<endl;
}
//system("pause");
cout<<"合并首次生成的标签数组。。。。"<<endl<<endl; //合并：像素值相同，但是标签不同，就把标签值大的变为小的。
for(i=;i<nYSize;i++)//行
{
for(j=;j<nXSize;j++)//列
{
if(i==)//第0行,只判左边的！
{
if(j==)//第一个元素
{
j=;//跳过第一个元素，直接从第二个元素：data[0][1]判断
//判断并执行合并：
if(data[i][j-]==data[i][j] && dataLabel[i][j-]!=dataLabel[i][j])//像素值相同 && 标签不同：合并！
{
//把所有的大标签替换为当前两个标签中的一个较小的值
//执行一次遍历
//cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
for(int m=;m<nYSize;m++)
{
for(int n=;n<nXSize;n++)
{//把大的标签替换为小的标签
if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i][j-]))
{
dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i][j-]);
}
}
}
}
}
}
else//非0行
{
if(j==)//第0列，但不是第一行的：只判断上、右上两个方向
{
if(data[i-][j]==data[i][j] && dataLabel[i-][j]!=dataLabel[i][j])//像素值相同 && 标签不同：合并！
{
//把所有的大标签替换为当前两个标签中的一个较小的值
//执行一次遍历
//cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
for(int m=;m<nYSize;m++)
{
for(int n=;n<nXSize;n++)
{//把大的标签替换为小的标签
if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i-][j]))
{
dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i-][j]);
}
}
}
}
if(data[i-][j+]==data[i][j] && dataLabel[i-][j=]!=dataLabel[i][j])//像素值相同 && 标签不同：合并！
{
//把所有的大标签替换为当前两个标签中的一个较小的值
//执行一次遍历
//cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
for(int m=;m<nYSize;m++)
{
for(int n=;n<nXSize;n++)
{//把大的标签替换为小的标签
if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i-][j+]))
{
dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i-][j+]);
}
}
}
}
}
else if(j==nYSize-)//非0行且最后一列的：判断左、左上、上三个方向
{
if(data[i][j-]==data[i][j] && dataLabel[i][j-]!=dataLabel[i][j])//像素值相同 && 标签不同：合并！
{
//把所有的大标签替换为当前两个标签中的一个较小的值
//执行一次遍历
//cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
for(int m=;m<nYSize;m++)
{
for(int n=;n<nXSize;n++)
{//把大的标签替换为小的标签
if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i][j-]))
{
dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i][j-]);
}
}
}
}
if(data[i-][j-]==data[i][j] && dataLabel[i-][j-]!=dataLabel[i][j])//像素值相同 && 标签不同：合并！
{
//把所有的大标签替换为当前两个标签中的一个较小的值
//执行一次遍历
//cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
for(int m=;m<nYSize;m++)
{
for(int n=;n<nXSize;n++)
{//把大的标签替换为小的标签
if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i-][-]))
{
dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i-][j-]);
}
}
}
}
if(data[i-][j]==data[i][j] && dataLabel[i-][j]!=dataLabel[i][j])//像素值相同 && 标签不同：合并！
{
//把所有的大标签替换为当前两个标签中的一个较小的值
//执行一次遍历
//cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
for(int m=;m<nYSize;m++)
{
for(int n=;n<nXSize;n++)
{//把大的标签替换为小的标签
if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i-][j]))
{
dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i-][j]);
}
}
}
} }
else//非0行且（既不是第一列，也不是最后一列的）：判断左、左上、上、右上四个方向
{
if(data[i][j-]==data[i][j] && dataLabel[i][j-]!=dataLabel[i][j])//像素值相同 && 标签不同：合并！
{
//把所有的大标签替换为当前两个标签中的一个较小的值
//执行一次遍历
//cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
for(int m=;m<nYSize;m++)
{
for(int n=;n<nXSize;n++)
{//把大的标签替换为小的标签
if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i][j-]))
{
dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i][j-]);
}
}
}
}
if(data[i-][j-]==data[i][j] && dataLabel[i-][j-]!=dataLabel[i][j])//像素值相同 && 标签不同：合并！
{
//把所有的大标签替换为当前两个标签中的一个较小的值
//执行一次遍历
//cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
for(int m=;m<nYSize;m++)
{
for(int n=;n<nXSize;n++)
{//把大的标签替换为小的标签
if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i-][-]))
{
dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i-][j-]);
}
}
}
}
if(data[i-][j]==data[i][j] && dataLabel[i-][j]!=dataLabel[i][j])//像素值相同 && 标签不同：合并！
{
//把所有的大标签替换为当前两个标签中的一个较小的值
//执行一次遍历
//cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
for(int m=;m<nYSize;m++)
{
for(int n=;n<nXSize;n++)
{//把大的标签替换为小的标签
if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i-][j]))
{
dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i-][j]);
}
}
}
}
if(data[i-][j+]==data[i][j] && dataLabel[i-][j+]!=dataLabel[i][j])//像素值相同 && 标签不同：合并！
{
//把所有的大标签替换为当前两个标签中的一个较小的值
//执行一次遍历
//cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
for(int m=;m<nYSize;m++)
{
for(int n=;n<nXSize;n++)
{//把大的标签替换为小的标签
if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i-][j+]))
{
dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i-][j+]);
}
}
}
}
}
} }
}
t=;
//system("pause");
cout<<"输出合并后的标签数组dataLabel:"<<endl;
for(i=;i<nYSize;i++)
{
for(j=;j<nXSize;j++)
{
if(dataLabel[i][j]>t)
{
t++;
}
cout<<setw()<<dataLabel[i][j]; }
cout<<endl;
}
cout<<endl<<"t="<<t<<endl;
cout<<"联通区域标记算法结束end！！！"<<endl;
//system("pause");
int max=;
for(int i=;i<nYSize;i++)
{
for(int j=;j<nXSize;j++)
{
if(dataLabel[i][j]>max)
{
max=dataLabel[i][j];
}
}
}
cout<<"maxLabel="<<max<<endl;
int *Pid=new int [max+];
int *Parea=new int [max+];
for(i=;i<max+;i++)
{
Pid[i]=i;
Parea[i]=;
}
/////////////////////////////
for(int i=;i<nYSize;i++)
{
for(int j=;j<nXSize;j++)
{
Parea[dataLabel[i][j]]++;
}
}
t=max+;
for(i=;i<t;i++)
{ cout<<"dataLabel为"<<i<<"的面积为:"<<Parea[i]<<endl;
} double Xi=,Si=;
int NPatch=;
for(i=;i<t;i++)
{
if(Parea[i]!=)
{
NPatch++;
Xi=Xi+Parea[i];//求出板块总面积
}
/*cout<<"NPatch="<<NPatch<<endl;
cout<<"Xi="<<Xi<<endl;*/
}
cout<<"NPatch="<<NPatch<<endl;
//cout<<"Xi="<<Xi<<endl;
Xi=Xi/NPatch;//面积平均数
cout<<"Initial Value:Xi="<<Xi<<endl;
cout<<"Initial Value:Si="<<Si<<endl;
//计算出所有
for(i=;i<t;i++)
{
if(Parea[i]!=)
{
Si=Si+((Parea[i]-Xi)*(Parea[i]-Xi))/NPatch;
}
}
cout<<"方差="<<Si<<endl;
Si=sqrt(Si);
cout<<"Standard Deviation(标准差)="<<Si<<endl;
cout<<endl; /****输出内容写出到文件当中****/
ofstream ocout;
//ocout.open("result.csv");
//ocout.open("result.txt");
ocout.open("result.csv");
/****将计算出的结果输出到屏幕****/
ocout<<"PatchID"<<","<<"Area"<<","<<"CSD"<<endl;
for(i=;i<t;i++)
{
if(Parea[i]!=)
{ //ocout<<"Patch: Label="<<i<<setw(4)<<" Area="<<Parea[i]<<setw(10)<<"CSD="<<(Parea[i]-Xi)/Si<<endl;
ocout<<i<<",";
ocout<<Parea[i]<<",";
ocout<<(Parea[i]-Xi)/Si<<","<<endl;
}
}
/*******后面的这一部分没有太大的实质性用处，可有可无**********/
const char *pszFormat = "GTiff";
GDALDriver *poDriver;
char **papszMetadata; poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat); if( poDriver == NULL )
exit( ); papszMetadata = poDriver->GetMetadata();
if( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATE, FALSE ) )
printf( "Driver %s supports Create() method.\n", pszFormat );
if( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATECOPY, FALSE ) )
printf( "Driver %s supports CreateCopy() method.\n", pszFormat ); /*释放读取文件用的指针*/
CPLFree(pafScanline);
//关闭文件
GDALClose((GDALDatasetH)poDataset);
free(dataLabel);
free(data);
free(Parea);
free(Pid);
//
ocout.close();
}
int minLabel(int a,int b)
{
return (a>b)?b:a;
}
int maxLabel(int a,int b)
{
return (a<b)?b:a;
}
double csd(int *Parea,int length)
{
return ;
}