使用GDAL实现DEM的地貌晕渲图(一)

目录

• 1. 原理
  • 1) 点法向量
  • 2) 日照方向
    • (1) 太阳高度角和太阳方位角
    • (2) 计算过程
  • 3) 晕渲强度
• 2. 实现
• 3. 参考

1. 原理

以前一直以为对DEM的渲染就是简单的根据DEM的高度不同赋予不同的颜色就可以。后来实际这么做的时候获取的效果跟别的软件相比，根本体现不出地形起伏的变化。如果要体现出地形的起伏变化，需要得到地貌晕渲图才行。晕渲法假设地形接受固定于某一位置光源的平行光线，随坡面与光源方向的夹角不同，产生不同色调明暗效果。
根据文献[1][2],可以通过计算DEM格网点的法向量与日照方向的的夹角，来确定该格网点的像素值。

1) 点法向量

我们知道三点成面，面的法向量就是其三个顶点的法向量（三点成面计算其法向量可参看《已知三点求平面法向量》）。但是一个顶点可能会构成多个不同的面，那么这种存在多个面的顶点的法向量怎么求呢？其实很简单，只需要把该点对应面的法向量相加就可以了。可以不用求平均，因为反正最后是要正规化的。
具体到DEM上来说，可以将一个DEM的矩形网格分成两个同样顺序排列的三角形，每个点涉及1到6个不等的面法向量。将这些面法向量相加并正则化，就得到了每个点的法向量。如下图所示。

2) 日照方向

关于日照方向，我在《通过OSG实现对模型的日照模拟》这篇文章里面有过详细的表述，那么这里就直接搬运过来。

(1) 太阳高度角和太阳方位角

对于太阳光照来说，其方向并不是随便设置的。这里需要引入太阳高度角和太阳方位角两个概念，通过这两个角度，可以确定日照的方向。
太阳高度角指的就是太阳光的入射方向和地平面之间的夹角；而太阳方位角略微复杂点，指的是太阳光线在地平面上的投影与当地子午线的夹角，可近似地看作是竖立在地面上的直线在阳光下的阴影与正南方向的夹角。其中方位角以正南方向为0，由南向东向北为负，有南向西向北为正。例如太阳在正东方，则其方位角为-90度；在正东北方时，方位角为-135度；在正西方时，方位角是90度，在正西北方为135度；当然在正北方时方位角可以表示为正负180度。
DEM渲染中一般将太阳高度角设置成45度，太阳方位角设置成315度，即西北照东南。

(2) 计算过程

根据上述定义，对于空间某一点的日照光线，可以有如下示意图。

令太阳光线长度L1=1，有如下推算过程：

α是太阳高度角，则日照方向Z长度L3=sin(α)；
L1在地平面（XY）平面的长度L2 = cos(α)；
β是太阳方位角，则日照方向X长度L5 = L2cos(β);
同时日照方向Y长度L4 = L2sin(β)。

因此，对于太阳高度角α和太阳方位角β，日照光线的单位向量n(x,y,z)为：

X = cos(α)cos(β);
Y = cos(α)sin(β);
Z = sin(α);

3) 晕渲强度

在文献[1][2]中提出由格网点法向量与光源方向的夹角，确定当前格网点的晕渲强度值。其晕渲图像素值i_cellvalue_hillshade计算公式如下所示（其中d_vectorvalue是法向量，a_rayvector是日照方向向量）：

这里的夹角d_raytovector_angle的计算公式略微奇怪。其实夹角计算远那么复杂，如果法向量和日照方向向量都已经正则化，那么其夹角可以直接用向量点积公式：

d_vectorvalue * a_rayvector = |d_vectorvalue| * |a_rayvector|* cos(d_raytovector_angle) = cos(d_raytovector_angle)

即其夹角的余弦值为两个正则化向量的点积。经过验证，两者计算出来的夹角值是一致的。

2. 实现

根据上述原理，其具体实现如下。我这里用到了GDAL来读写DEM和图像，此外还有向量计算用到了osg库里面的内容，如果没有osg，可以自己简单实现下，都是很简单的数学知识。

#include <iostream>
#include <algorithm>
#include <gdal_priv.h>
#include <osg/Vec3d>
#include <fstream>

using namespace std;
using namespace osg;

//计算三点成面的法向量
void Cal_Normal_3D(const Vec3d& v1, const Vec3d& v2, const Vec3d& v3, Vec3d &vn)
{
	//v1(n1,n2,n3);
	//平面方程: na * (x – n1) + nb * (y – n2) + nc * (z – n3) = 0 ;
	double na = (v2.y() - v1.y())*(v3.z() - v1.z()) - (v2.z() - v1.z())*(v3.y() - v1.y());
	double nb = (v2.z() - v1.z())*(v3.x() - v1.x()) - (v2.x() - v1.x())*(v3.z() - v1.z());
	double nc = (v2.x() - v1.x())*(v3.y() - v1.y()) - (v2.y() - v1.y())*(v3.x() - v1.x());

	//平面法向量
	vn.set(na, nb, nc);
}

int main()
{
	GDALAllRegister();          //GDAL所有操作都需要先注册格式
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");  //支持中文路径

	const char* demPath = "D:/CloudSpace/我的技术文章/素材/DEM的渲染/dst.tif";
	//const char* demPath = "D:/Data/imgDemo/K51E001022/k51e001022dem/w001001.adf";
	
	GDALDataset* img = (GDALDataset *)GDALOpen(demPath, GA_ReadOnly);
	if (!img)
	{
		cout << "Can't Open Image!" << endl;
		return 1;
	}

	int imgWidth = img->GetRasterXSize();   //图像宽度
	int imgHeight = img->GetRasterYSize();  //图像高度
	int bandNum = img->GetRasterCount();    //波段数
	int depth = GDALGetDataTypeSize(img->GetRasterBand(1)->GetRasterDataType()) / 8;    //图像深度

	GDALDriver *pDriver = GetGDALDriverManager()->GetDriverByName("GTIFF"); //图像驱动
	char** ppszOptions = NULL;
	ppszOptions = CSLSetNameValue(ppszOptions, "BIGTIFF", "IF_NEEDED"); //配置图像信息
	const char* dstPath = "D:\dst.tif";
	int bufWidth = imgWidth;
	int bufHeight = imgHeight;
	int dstBand = 1;
	int dstDepth = 1;
	GDALDataset* dst = pDriver->Create(dstPath, bufWidth, bufHeight, dstBand, GDT_Byte, ppszOptions);
	if (!dst)
	{
		printf("Can't Write Image!");
		return false;
	}

	dst->SetProjection(img->GetProjectionRef());
	double padfTransform[6] = { 0 };
	if (CE_None == img->GetGeoTransform(padfTransform))
	{
		dst->SetGeoTransform(padfTransform);
	}

	//申请buf
        depth = 4;
	size_t imgBufNum = (size_t)bufWidth * bufHeight * bandNum;
	float *imgBuf = new float[imgBufNum];
	//读取
	img->RasterIO(GF_Read, 0, 0, bufWidth, bufHeight, imgBuf, bufWidth, bufHeight,
		GDT_Float32, bandNum, nullptr, bandNum*depth, bufWidth*bandNum*depth, depth);

	if (bandNum != 1)
	{
		return 1;
	}

	double startX = padfTransform[0];			//左上角点坐标X
	double dx = padfTransform[1];			//X方向的分辨率
	double startY = padfTransform[3]; 			//左上角点坐标Y
	double dy = padfTransform[5];			//Y方向的分辨率

	//
	double minZ = DBL_MAX;
	double maxZ = -DBL_MAX;
	double noValue = img->GetRasterBand(1)->GetNoDataValue();

	vector<Vec3d> dotList;			//所有的顶点	
	for (int yi = 0; yi < bufHeight; yi++)
	{
		for (int xi = 0; xi < bufWidth; xi++)
		{
			size_t m = (size_t)bufWidth * yi + xi;
			double x = startX + xi * dx;
			double y = startY + yi * dy;
			double z = imgBuf[m];
			dotList.push_back(Vec3d(x, y, z));

			if (abs(z - noValue) < 0.01 || z<-11034 || z>8844.43)
			{
				continue;
			}

			minZ = (std::min)(minZ, z);
			maxZ = (std::max)(maxZ, z);
		}
	}

	//计算每个面的法向量
	multimap<size_t, size_t> dot_face;
	vector<Vec3d> faceNomalList;
	for (int yi = 0; yi < bufHeight - 1; yi++)
	{
		for (int xi = 0; xi < bufWidth - 1; xi++)
		{
			size_t y0x0 = (size_t)bufWidth * yi + xi;
			size_t y1x0 = (size_t)bufWidth *(yi + 1) + xi;
			size_t y0x1 = (size_t)bufWidth *yi + xi + 1;
			size_t y1x1 = (size_t)bufWidth *(yi + 1) + xi + 1;

			Vec3d vn;
			Cal_Normal_3D(dotList[y0x0], dotList[y1x0], dotList[y0x1], vn);
			dot_face.insert(make_pair(y0x0, faceNomalList.size()));
			dot_face.insert(make_pair(y1x0, faceNomalList.size()));
			dot_face.insert(make_pair(y0x1, faceNomalList.size()));
			faceNomalList.push_back(vn);

			Cal_Normal_3D(dotList[y1x0], dotList[y1x1], dotList[y0x1], vn);
			dot_face.insert(make_pair(y1x0, (int)faceNomalList.size()));
			dot_face.insert(make_pair(y1x1, (int)faceNomalList.size()));
			dot_face.insert(make_pair(y0x1, (int)faceNomalList.size()));
			faceNomalList.push_back(vn);
		}
	}


	//申请buf
	size_t dstBufNum = (size_t)bufWidth * bufHeight * dstBand * dstDepth;
	GByte *dstBuf = new GByte[dstBufNum];
	memset(dstBuf, 255, dstBufNum*sizeof(GByte));

	//设置方向：平行光
	double solarAltitude = 45.0;
	double solarAzimuth = 315.0;
	osg::Vec3d arrayvector(0.0f, 0.0f, -1.0f);
	double fAltitude = osg::DegreesToRadians(solarAltitude);				//光源高度角
	double fAzimuth = osg::DegreesToRadians(solarAzimuth);			//光源方位角
	arrayvector[0] = cos(fAltitude)*cos(fAzimuth);
	arrayvector[1] = cos(fAltitude)*sin(fAzimuth);
	arrayvector[2] = sin(fAltitude);
	
	vector<Vec3d> normalList;
	double alpha = 0.5;		//A不透明度 α*A+(1-α)*B
	for (int yi = 0; yi < bufHeight; yi++)
	{
		for (int xi = 0; xi < bufWidth; xi++)
		{
			size_t m = (size_t)bufWidth * yi + xi;

			auto beg = dot_face.lower_bound(m);
			auto end = dot_face.upper_bound(m);

			Vec3d n(0, 0, 0);
			int ci = 0;
			for (auto it = beg; it != end; ++it)
			{
				n = n + faceNomalList[it->second];
				ci++;
			}

			n.normalize();
			normalList.push_back(n);

			double angle = osg::RadiansToDegrees(acos(n * arrayvector));
			//double d_tmp = (n - arrayvector).length();
			//double angle = osg::RadiansToDegrees(asin(d_tmp / 2.0)) * 2;

			double value = (std::min)((std::max)(angle / 90 * 255, 0.0), 255.0);
			dstBuf[m] = (GByte)value;
		}
	}

	//写入
	dst->RasterIO(GF_Write, 0, 0, bufWidth, bufHeight, dstBuf, bufWidth, bufHeight,
		GDT_Byte, dstBand, nullptr, dstBand*dstDepth, bufWidth*dstBand*dstDepth, dstDepth);
	
	//释放
	delete[] imgBuf;
	imgBuf = nullptr;

	//释放
	delete[] dstBuf;
	dstBuf = nullptr;

	//
	GDALClose(dst);
	dst = nullptr;

	GDALClose(img);
	img = nullptr;

	return 0;
}

最后得到的效果与ArcMap里面生成的晕渲效果比较如下，应该还是比较接近的：

这里只是得到了晕渲的灰白强度图，后续会继续实现彩色晕渲图的实现。

3. 参考

[1].地貌晕渲图的生成原理与实现.丁宇萍,蒋球伟
[2].DEM-地貌晕渲图的生成原理