GDAL计算栅格图像统计值相关说明

一、        简介

在图像处理的过程中,会经常获取图像的统计值来进行一些计算。这里说的图像的统计值指:直方图,最大值,最小值,均值和方差这几个。下面分别对怎么使用GDAL来计算这些值做一个简单的说明,希望对大家有用。

下面将分为三个部分进行说明,分别是:直方图统计、极值统计、均值标准方差统计。

二、        统计直方图

关于直方图,这里不多介绍。这里主要对如何使用GDAL获取直方图进行说明。使用GDAL获取直方图的函数叫做GDALRasterBand::GetHistogram(),下面对这个函数的参数进行一个大致的说明。函数GDALRasterBand::GetHistogram的定义如下:

CPLErrGDALRasterBand::GetHistogram   (
double dfMin,
double dfMax,
int nBuckets,
int *     panHistogram,
int bIncludeOutOfRange,
int bApproxOK,
GDALProgressFunc   pfnProgress,
void *   pProgressData  
)
为了更形象的说明直方图,将结合图1进行说明。图1是一个8bit的图像的全部直方图像素值的一个列表,每个横向的线段代表16个灰度值间隔。图中的数字表示右下角交叉处的灰度值(图中红色圈的位置表示的是0和256的位置,可能有点不清楚,瞪大眼睛仔细看~~)。


图1 直方图示例

第一个参数和第二个参数,即dfMin和dfMax,这两个参数的意思是指统计直方图的最小值和最大值,在图1中的话,就是这个横向数轴的一个区间。比如一个8bit的图像,你想统计所有的直方图,那么这两个值可以设置为-0.5和255.5。至于为什么不是0和255,是因为GDAL的这个函数是个开区间,就是如果设置为0和255的话,那么0和255这两个值是不会被统计的,所以需要把最小值往小减小一点,最大值往大增加一点。-0.5和255.5只是我的习惯,你也可以用-0.1和255.1,只要把0和255包括进去即可。如果你只想统计128到255之间的直方图,那么把最小值设置为127.5即可;如果想统计0到128之间的直方图,那么把最大值设置为128.5即可。

第三个参数nBuckets,表示直方图统计的份数,什么意思呢?举例说明一下,我要统计一个8bit的图像的某个波段的直方图,我希望每个灰度值都单独统计一下,就是0有多少个,1有多少个,2有多少个,等等,那么这样的话,我就要把0~255(8bit图像的像素存储范围)一共256个数分成256份,那么这个参数就是256。如果我希望没两个灰度值统计在一起,就是0和1统计一个,2和3统计一个,一直到254和255统计一个,那么这样的话,份数就是128个,那么这个值就是128。其他情况自己类比就出来了。这个参数还有一个作用,就是用来指定第四个参数的数组个数。

第四个参数panHistogram就是用来存储直方图的数组,这个数组的个数由第三个参数确定,如果这个数组的个数没有第三个参数大,程序会崩溃滴,切忌,最好设置为一样的。

第五个参数bIncludeOutOfRange,从字面意思就能看出来,就是是否统计区间外的值。从第一个和第二个参数可以得知,统计直方图是有一个区间的。如果这个参数设置为TRUE,那么图像中的像元值小于最小值的像元值将被统计到直方图数组中的第一个里面去,图像中的像元值大于最大值的像元会被统计到直方图数组中的最后一个里面去。如果设置为FALSE,那么图像中的像元值小于最小值的像元不进行统计,同样,像元值超过最大值的像元值也不进行统计。

第六个参数bApproxOK,表示是否进行粗略统计。由于图像大的话,统计直方图会很慢,为了解决这个问题,GDAL提供了一个粗略的直方图统计方式,就是用这个参数来进行指定。如果这个参数设置为TRUE,那么统计直方图的时候为了速度快,但是只会统计一个大致的直方图,而不是遍历图像中的所有像素来进行统计。如果该参数设置为FALSE,那么统计的直方图就是图像的准确的直方图,速度比较慢。

最后两个参数pfnProgress和pProgressData表示进度条的回调函数和进度条数据。关于这两个参数,可以参考我之前的博客《GDAL算法进度条使用说明》。

有人会问,如果图像是浮点数的话,这个直方图怎么处理。比如一个32F的数据,其数据范围大致是-3.40E+38~ +3.40E+38,如果是一个64F的数据,其数据范围大致是-1.79E+308~ +1.79E+308。由于浮点型图像的像素值是个小数,不像8bit或者16bit之类的整数类型,其像元值是个整数,所以在统计直方图的时候都是按照一个区间来进行统计的。比如,一个浮点数图像(NDVI之类的图像就是浮点数),他的大致范围是-100.5到200.9,我要把这个范围的值统计为100份,那么调用这个函数的时候,应该这么写:

int panHistogram[100] = {0};
pBand->GetHistogram    (-101, 201, 100, panHistogram, TRUE, FALSE,NULL, NULL);
最后在说明一点,就是如果调用了这个GetHistogram()函数,GDAL会把直方图信息写入一个xml文件,这个xml的位置与图像位置在一起,文件名是在图像的名字后面加了一个.aux.xml。下一次再调用这个函数是,GDAL会先从这个xml文件中读取,如果读取失败再进行统计。所以看到这个.aux.xml的文件不要删,可以大幅度的减少图像的统计时间,尤其是对大图像更为明显。

三、        统计极值

极值的统计就是计算图像某个波段的最大值和最小值。在GDAL中,里面有两个函数可以来进行计算,分别是GDALRasterBand::GetMinimum()、GDALRasterBand::GetMaximum() 和GDALRasterBand::ComputeRasterMinMax ()。前两个函数的原型如下:

doubleGDALRasterBand::GetMinimum(int * pbSuccess = NULL )
double GDALRasterBand::GetMaximum(int* pbSuccess = NULL )

这两个函数很像,返回值是一个double值,表示获取的最大值或者最小值。同时这两个函数还有个参数叫pbSuccess,从字面意思看就是是否成功的意思。通过传递一个变量的地址来进行返回,如果这个pbSuccess的值是TRUE,那么就说明获取的最大值或者最小值是正确的。如果pbSuccess的值是FALSE,那么说明获取的最大值或者最小值是错的。这个时候要统计图像的极值就需要函数ComputeRasterMinMax ()出马了。函数ComputeRasterMinMax ()的原型是:

CPLErr GDALRasterBand::ComputeRasterMinMax(int bApproxOK, double *pdfMinMax)

这个函数的第一个参数bApproxOK和直方图函数里面的意义一样,FALSE表示精确统计,速度慢,图像的所有像元都遍历一边,TRUE表示粗略统计,速度快,但是不一定准确。函数还有个返回值是CPLErr类型,其实就是个int值,0就是CE_None,就是表示计算成功,其他的可以参考GDAL的相关说明。

第二个参数pdfMinMax,就是用来存储统计出来的最小值和最大值。这个数组一般是一个double [2],第0个表示最小值,第1个表示最大值,代码示例如下:

double dMinMax[2] = {0.0,255.0};
m_pDS->GetRasterBand(1)->ComputeRasterMinMax(FALSE,dMinMax);
m_dMin = dMinMax[0];
m_dMax = dMinMax[1];

同样的,调用了函数GDALRasterBand::ComputeRasterMinMax ()之后,GDAL也会把最大值和最小值写入.aux.xml中,下次就不需要再统计一遍了。如果有这个.aux.xml文件,直接调用函数GDALRasterBand::GetMinimum()和GDALRasterBand::GetMaximum()就能获取到准确的结果。

四、        统计均值标准方差

均值和标准方差的统计的函数叫GDALRasterBand::ComputeStatistics()。这个函数的原型是:

CPLErrGDALRasterBand::ComputeStatistics(
int bApproxOK,
double * pdfMin,
double * pdfMax,
double * pdfMean,
double * pdfStdDev,
GDALProgressFunc pfnProgress,
void * pProgressData )

第一个参数就不说了,和上面两个一样。

接下来的四个参数,分别是用来存储计算的最小值,最大值,均值和标准方差的值。

这个函数和前面的一样,也会将统计的结果存储在.aux.xml文件中,下次获取的时候就不用再统计了。

最后两个参数就是进度条函数指针和进度条参数,具体参考《GDAL算法进度条使用说明》。

五、        参考资料

[1] http://www.gdal.org/classGDALRasterBand.html

[2] http://blog.csdn.net/liminlu0314/article/details/7276954

原文地址:https://www.cnblogs.com/xiaowangba/p/6313991.html