使用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-地貌暈渲圖的生成原理