# 1. 原理

## 2) 日照方向

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

DEM渲染中一般将太阳高度角设置成45度，太阳方位角设置成315度，即西北照东南。

### (2) 计算过程

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

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

## 3) 晕渲强度

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

# 2. 实现

#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";

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 };
{
}

//申请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 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);
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;
}


# 3. 参考

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

posted @ 2019-07-14 17:37  charlee44  阅读(1761)  评论(0编辑  收藏  举报