遥感影像显示相关的技术总结

前言

    从事遥感影像和图像处理有一段时间了，今天就把遥感影像显示相关的技术和大家分享一下。

    平常我们用的GIS软件或者说遥感软件都能讲遥感影像的数据显示在屏幕上，并且有些显示效果还不错，其中ENVI的显示效果是业界做得比较好的，尤其是ENVI5.0之后的大视图，能够根据真彩色的波段自动选择波段进行显示。遥感影像显示其实就是图像显示，和我们生活中常见的图像显示没有太大区别，基本上原理是一样的，但也有自己独特的地方，之后会讲到。首先，我们来看看图像显示的基本过程，即图像是怎么显示在计算机屏幕上的。

1、影像显示基本过程

    首先，CPU从存储介质中读取数据，存储介质可以是本地磁盘、移动存储设备以及远程服务器。然后将数据存储在内存缓冲区里面，然后读数器将这里面的数据根据彩色查找表转换为RGB的形式。然后模数转换器将RGB形式的数据值转换为适当的电信号，这样就得到模拟信号，该模拟信号调整RGB电子枪的强度，控制着每个像素在视频CRT屏幕上显示的亮度。

1.1全色显示（灰度显示）

    比如我们的遥感影像数据源，有很多数据时全色波段，例如高分一号的全色波段就是2米分辨率的，那么这个数据时如何显示在屏幕上？如果数据时8位的，即byte型，那么可以做拉伸也可以不用拉伸直接显示在屏幕上。这里就需要用到8位的图像处理器，它是一个8位连续的256个元素的查找表，并且RGB每个分量都是同一个值。比如黑色就是RGB（0,0,0）白色就是RGB（255,255,255）。最后通过数模转换显示在屏幕上。

1、2彩色显示

    说到彩色显示，可能更复合我们人眼的观察，分辨的能力更强。彩色显示主要用两种方案，一种是电子显示法，在数字图像显示中主要用这种方法。它是用彩色监视器显示，有时候将这种显示叫做软拷贝。另外一种是用彩色硬拷贝设备进行显示。这种显示方法主要用再印刷中，通常需要将RGB颜色空间转换到CYMK颜色空间，这种显示不在我们的讨论范围之内。

    对于遥感影像，如果选择的波段与RGB是一一对应的，那么显示出来的效果即为真彩色效果，接近于人类肉眼观察到的效果。如果输入波段与RGB不相对应，那么僵产生假彩色效果，就和我们人眼观察的效果有较大区别。

2、图像的拉伸

    图像拉伸主要是增强图像的显示效果，如果图像的对比度比较低，那么就无法看清楚地物，这时需要拉伸。另外一个方面，对于非8位的数据，我们一般也做拉伸处理后再显示。一般的拉伸有线性拉伸和非线性拉伸两种方法。

    线性拉伸主要有全域线性拉伸，分段线性拉伸等。具体有最大-最小值拉伸、标准差拉伸。非线性拉伸主要有指数拉伸，对数拉伸等。对于多波段的遥感影像，一般是将各个波段分别拉伸后，然后进行波段组合后显示。

3、编程方面的一些细节

3.1影像IO   

    对于遥感影像的读取，绝对推荐GDAL来做，你自己来解析的话很累。在某些需要提高IO速度的地方可以先将数据转储为一个裸数据，然后通过内存映射文件的方式来做，这样的效率很高。需要注意的一个地方是每次映射的时候要分块映射，在windows系统中，系统分配的粒度是64k，可以每次映射这么大的数据，分多次映射，每次映射的起始偏移量必须是64k的整数倍，偏移量是以字节为单位。

3、2关于拉伸的实现

首先是最大最小值拉伸，先在原始影像中统计最大最小值。最大和最小值拉伸的代码如下：

void MinMaxStreh(unsigned char* poData,int nLen,double dbMin,double dbMax)
{
	#pragma omp parallel for
	for (int i = 0;i < nLen; i ++)
	{
		poData[i] = double(poData[i]-dbMin)/(double)(dbMax-dbMin) * 255;

	}
}

另外一种比较常用的拉伸方法是标准差拉伸，这种拉伸方法实际上也是最大最小值拉伸，只不过在拉伸强需要先将最大最小值计算出来。代码如下：

//图像的标准差拉伸
template <typename T>
void StandDevStreh(const T* poData,unsigned char* puBytes,int nLen,double dbDev,int nDevCount,double dbMean)
{
	#pragma omp parallel for
	for (int i = 0;i < nLen; i ++)
	{
		double ucMax  = dbMean + nDevCount * dbDev;
		double ucMin  = dbMean - nDevCount * dbDev;

		if (poData[i] < ucMin)
		{
			puBytes[i] = 0;
		}

		else if (poData[i] > ucMax)
		{
			puBytes[i] = 255;
		}

		else if (poData[i] >= ucMin && poData[i] <= ucMax)
		{
			puBytes[i] = double(poData[i]-ucMin)/(double)(ucMax-ucMin) * 255;
		}

	}
}

    这里的拉伸都需要对原始影像做统计，如果是对整个波段的原始数据逐个统计，如果图像很大的话，速度绝对的慢，那么可以将原始影像的数据读到一个256*256的缓冲区里面或者512*512的缓冲区，也可以按照比例缩放，比如X方向上定死512个像素，那么Y方向的像素的个数就可以根据原始影像波段的高度和宽度之比来计算得到。通过一个小的缓冲区统计出来的值和原始数据统计出来的值差别很小，因为是对原始影像的抽样，概率密度分布基本上是一样的。

3.3、关于影像黑边的显示处理

影像黑边的话，可以默认显示，也可以将它过滤掉，比如将黑边显示为白色，其效果如下：

有黑边的显示

无黑边的显示

3.4、关于和矢量数据的叠加

    这个需要注意的是要根据地理坐标进行计算显示的像素范围。首先将屏幕当前显示的地理范围和影像的地理范围求交集，然后将这个交集映射到影像的行列号，即可以取到数据进行显示了。待显示的数据缓冲区大小是地理范围求交集所占的屏幕像素的大小。具体可以参考我的另外一篇博客：http://blog.csdn.net/zhouxuguang236/article/details/21882801。

3.5、关于优化显示效果

    这里讲一讲如何处理影像黑边对波段统计结果的影像。如果将黑边统计进去，然后进行拉伸的话，势必会对影像的显示效果产生影像，一般影像的黑边是原始数据经过正射校正或者匹配校正后产生的，总之，对图像进行纠正后都会产生黑边。其实，校正后的黑边面积比较大，可以按照3.2中提到的方法，进行统计时过滤掉黑边，具体操作可以在进行直方图统计时将比例最大的那一部分看成黑边，然后找到这个值，最后计算均值标准差的时候过滤掉就可以了。这只是我想到的一种方法，有没有更好的方法，敬请期待！下面两幅图是没有过滤黑边和过滤掉黑边显示效果的对比，明显过滤掉黑边后进行标准差拉伸更接近真彩色效果。

没有对黑边处理的效果

对黑边处理的效果

对比这两幅图，效果明显的是统计时过滤掉黑边的效果较好。

对于过滤掉黑边的统计函数如下：

#define TILE_SIZE 256

bool GeoGdalImageLayer::GetBandStatisInfoEx(
	int nChannelIndex, 
	float* pfHist, 
	float& fMax, 
	float& fMin, 
	float& fMean, 
	float& fVar, 
	float fExclude) const
{
	GDALDataset *poDataset = (GDALDataset*)m_poDataset;
	if (NULL == poDataset)
	{
		return false;
	}

	assert(NULL != poDataset);
	int nXsize = poDataset->GetRasterXSize();
	int nYsize = poDataset->GetRasterYSize();
	int pBandList[1] = {nChannelIndex+1};

	int nTileHeight = ( nYsize/float(nXsize) ) * TILE_SIZE;

	float *pBytes = new float[TILE_SIZE*nTileHeight];

	poDataset->RasterIO(GF_Read,0,0,nXsize,nYsize,pBytes,TILE_SIZE,nTileHeight,GDT_Float32,1,pBandList,0,0,0);

	std::vector<float> vecPixels;
	vecPixels.resize(TILE_SIZE*nTileHeight);

	int nCount = 0;
	double dbSum = 0;

	//开始统计
	fMax = DBL_MIN;
	fMin = DBL_MAX;
	for (int i = 0; i < TILE_SIZE*nTileHeight; i ++)
	{
		if (fabs( pBytes[i] - fExclude) <= 0.000001)
		{
			continue;
		}

		else
		{
			if (pBytes[i] <= fMin)
			{
				fMin = pBytes[i];
			}

			if (pBytes[i] >= fMax)
			{
				fMax = pBytes[i];
			}

			dbSum += pBytes[i];
			vecPixels[nCount] = pBytes[i];
			nCount ++;
		}
	}

	//计算均值
	fMean = dbSum/nCount;

	//计算标准差
	dbSum = 0;
	for (int i = 0; i < nCount; i ++)
	{
		dbSum += (vecPixels[i]-fMean)*(vecPixels[i]-fMean);
	}
	dbSum /= nCount;
	fVar = sqrt(dbSum);

	if (pBytes != NULL)
	{
		delete []pBytes;
		pBytes = NULL;
	}
	return 1;
}

3.6、关于影像的自动校正显示

    这个以前有碰到过这个需求，但是目前还没有实现。例如高分一号数据的自动校正显示，SPOT原始数据的自动校正显示等。这个可以再显示前做重采样实现。

4、总结

    以上都是个人总结，如有不妥之处，也请指出来共同探讨。影像显示的界面框架是QT界面库实现，具体的绘制通过OpenGL来绘制。其实说实在的，这个显示效果和ENVI的显示效果差太多了，还得研究ENVI的显示机制。

遥感影像显示相关的技术总结