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

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

using namespace std;
using namespace osg;

// a b c
// d e f
// g h i
double CalHillshade(float *tmpBuf, double Zenith_rad, double Azimuth_rad, double dx, double dy, double z_factor)
{
    double dzdx = ((tmpBuf[2] + 2 * tmpBuf[5] + tmpBuf[8]) - (tmpBuf[0] + 2 * tmpBuf[3] + tmpBuf[6])) / (8 * dx);
    double dzdy = ((tmpBuf[6] + 2 * tmpBuf[7] + tmpBuf[8]) - (tmpBuf[0] + 2 * tmpBuf[1] + tmpBuf[2])) / (8 * dy);

    double Slope_rad = atan(z_factor * sqrt(dzdx*dzdx + dzdy*dzdy));
    double Aspect_rad = 0;
    if (abs(dzdx) > 1e-9)
    {
        Aspect_rad = atan2(dzdy, -dzdx);
        if (Aspect_rad < 0)
        {
            Aspect_rad = 2 * PI + Aspect_rad;
        }
    }
    else
    {
        if (dzdy > 0)
        {
            Aspect_rad = PI / 2;
        }
        else if (dzdy < 0)
        {
            Aspect_rad = 2 * PI - PI / 2;
        }
        else
        {
            Aspect_rad = Aspect_rad;
        }
    }

    double Hillshade = 255.0 * ((cos(Zenith_rad) * cos(Slope_rad)) + (sin(Zenith_rad) * sin(Slope_rad) * cos(Azimuth_rad - Aspect_rad)));
    return Hillshade;
}


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
    size_t imgBufNum = (size_t)bufWidth * bufHeight * bandNum * depth;
    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方向的分辨率
        
    //申请buf
    size_t dstBufNum = (size_t)bufWidth * bufHeight * dstBand * dstDepth;
    GByte *dstBuf = new GByte[dstBufNum];
    memset(dstBuf, 0, dstBufNum*sizeof(GByte));

    //设置方向：平行光
    double solarAltitude = 45.0;
    double solarAzimuth = 315.0;
    
    //
    double Zenith_rad = osg::DegreesToRadians(90 - solarAltitude);
    double Azimuth_math = 360.0 - solarAzimuth + 90;
    if (Azimuth_math >= 360.0)
    {
        Azimuth_math = Azimuth_math - 360.0;
    }   
    double Azimuth_rad = osg::DegreesToRadians(Azimuth_math);

    //a b c
    //d e f
    //g h i
    double z_factor = 1;
    for (int yi = 1; yi < bufHeight-1; yi++)
    {
        for (int xi = 1; xi < bufWidth-1; xi++)
        {
            size_t e = (size_t)bufWidth * yi + xi;
            size_t f = e + 1;
            size_t d = e - 1;

            size_t b = e - bufWidth;
            size_t c = b + 1;
            size_t a = b - 1;

            size_t h = e + bufWidth;
            size_t i = h + 1;
            size_t g = h - 1;
            
            float tmpBuf[9] = { imgBuf[a], imgBuf[b], imgBuf[c], imgBuf[d], imgBuf[e], imgBuf[f], imgBuf[g],imgBuf[h], imgBuf[i] };
            double Hillshade = CalHillshade(tmpBuf, Zenith_rad, Azimuth_rad, dx, -dy, z_factor);
    
            dstBuf[e] = (GByte)(std::min(std::max(Hillshade, 0.0),255.0));
        }
    }

    //写入
    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;
}