当前位置：首页 > news >正文

osgEarth 图像融合正片叠底

news 2025/8/23 10:33:34

* 需求：
* 高程渲染图 RGB.tif、山体阴影图 AMP.tif
*
* 高程渲染图 rgb波段分别乘以山体阴影图r波段，然后除以255(AI说读取的纹理就已经归一化到了 0~1 范围，不用除以 255)。

本人遥感知识匮乏。

问了AI,以上需求在许多商业软件上已实现。在 ArcGIS 和 QGIS 中制作 Hillshade Overlay 的核心逻辑是通过数字高程模型（DEM）生成地形阴影，再将其与其他图层叠加以增强立体感。

环境：Win10 64bit, ARM R9 3900X ,Qt 5.15.2, C++，osgEarth3.7.0

1、C++ 结合 OpenCV实现。适合静态展示；

2、glsl实现。效率高。glsl 我也是抄、描、问，零零碎碎拼凑出来的，这里就想抛砖引玉，希望openGL大佬们来指导哈。

这完全是个笔记，借助AI拼凑出来的代码，希望遥感数据处理大牛们多多指教，

特别是osgEarth中的着色器编程。

无图无真相。先图后码。

原彩色浮雕图:

山体阴影

在Qt 里面使用 C++ OpenCV运行出来的效果，源码在后面。

osgearth着色器, minLumFactor = 0.125时的效果，明显偏暗。

minLumFactor = 0.618

好，把之前实现的代码和效果做个笔记。

代码如下：

#include "GenerateColorReliefThread.h"
#include "Global.h"
#include "MyFile/MyRaster.h"
#include "gdal_utils.h"
#include "osg/BlendFunc"
#include "qcoreapplication.h"#include "SettingsManager.h"#include "ProjectManage/MyProject.h"GenerateColorReliefThread::GenerateColorReliefThread( const QString &filePath, GDALDataset *dataset, const std::string &colorFilename, QObject *parent): QThread(parent), m_filePath(filePath), m_colorFilename(colorFilename), m_dataset(dataset)
{}GenerateColorReliefThread::~GenerateColorReliefThread()
{wait();
}bool GenerateColorReliefThread::createColorRelief( std::string &filePathColorRelief )
{CPLSetConfigOption("GDAL_NUM_THREADS", "ALL_CPUS"); // 启用所有CPU核心并行处理// 3) 调 GDALDEMProcessing 做 color-relief（或 hillshade）const char* dArgvStr = "-alpha -of GTiff";char**papszArgv = CSLTokenizeString(dArgvStr);// 处理选项GDALDEMProcessingOptions *psOptions = GDALDEMProcessingOptionsNew(papszArgv, nullptr);if (psOptions == nullptr){qDebugV5() <<"Failed to create processing options.";return false;}static std::atomic<int> tempFileCounter(0); // 避免内存文件名冲突filePathColorRelief = "/vsimem/color_relief_" + std::to_string(tempFileCounter++) + ".tif";// std::cout<<"memFileName: "<<filePathColorRelief<<std::endl;int pbUsageError = 0;GDALDatasetH hColorRelief = GDALDEMProcessing(filePathColorRelief.c_str(), m_dataset, "color-relief", m_colorFilename.c_str(), psOptions, &pbUsageError);if (!hColorRelief || pbUsageError != 0){qDebugV5() << "Processing failed with error code: " << pbUsageError;GDALDEMProcessingOptionsFree(psOptions);return false;}cv::Mat matHillShade;FileCategory fileCategory = Global::getFileCategory(m_filePath);switch (fileCategory.rasterSubCategory) {case RASTER_OPT:case RASTER_AMP:break;case RASTER_DEM:{switch (fileCategory.demSubSubCategory) {case DEM_COPERNICUS:matHillShade = MyProject::instance()->matHillShadeCopernicus().clone();break;case DEM_RAW:matHillShade = MyProject::instance()->matHillShadeRaw().clone();break;case DEM_CACHE:{double zFactor  = SettingsManager::instance().value("HILL_SHADE_ZFACTOR", 0.00001).toDouble();double scale    = SettingsManager::instance().value("HILL_SHADE_SCALE", 1.0).toDouble();double azimuth  = SettingsManager::instance().value("HILL_SHADE_AZIMUTH", 315.0).toDouble();double altitude = SettingsManager::instance().value("HILL_SHADE_ALTITUDE", 45.0).toDouble();bool   combined = SettingsManager::instance().value("HILL_SHADE_COMBINED", false).toBool();Global::generateHillShade(matHillShade, m_filePath, zFactor, scale, azimuth, altitude, combined);if(fileCategory.demSubSubCategory == DEM_CACHE){MyProject::instance()->setMatHillShadeCache( matHillShade );}}break;            }}break;case RASTER_MASK:matHillShade = MyProject::instance()->matHillShadeCache().clone();break;}// 获取图像尺寸int width  = m_dataset->GetRasterXSize();int height = m_dataset->GetRasterYSize();int panBandMap [3]= {1, 2, 3};// 转换为C++对象GDALDataset* poColorRelief = (GDALDataset*)hColorRelief;if(!poColorRelief){return false;}// 创建三通道矩阵一次性读取RGB数据cv::Mat matRGB(height, width, CV_8UC3);// 一次性读取三个波段CPLErr readErr = poColorRelief->RasterIO(GF_Read,0, 0,width, height,(void*)matRGB.data,width, height,GDT_Byte, 3, panBandMap,3, width*3, 1);if (readErr != CE_None) {qDebugV5() << "Failed to read RGB bands!";GDALClose(hColorRelief);GDALDEMProcessingOptionsFree(psOptions);return false;}// 将单通道hillshade转换为三通道以匹配RGBcv::Mat matHillShade3C;cv::cvtColor(matHillShade, matHillShade3C, cv::COLOR_GRAY2BGR);// 2. 分别对每个通道做 multiply（与你原来的分通道逻辑完全一致）cv::multiply(matRGB, matHillShade3C, matRGB, 1.0 / 255.0, CV_8UC3); // R通道// 一次性写入三个波段CPLErr writeErr = poColorRelief->RasterIO(GF_Write,0, 0,width, height,(void*)matRGB.data,width, height,GDT_Byte, 3, panBandMap,3, width*3, 1);if (writeErr != CE_None) {std::cerr << "GDAL 写入失败: " << CPLGetLastErrorMsg() << std::endl;GDALClose(hColorRelief);GDALDEMProcessingOptionsFree(psOptions);return false;}GDALClose(hColorRelief);GDALDEMProcessingOptionsFree(psOptions);return true;
}void GenerateColorReliefThread::run()
{// 创建并启动计时器QElapsedTimer timer;timer.start();std::string filePathColorRelief;bool ret = createColorRelief(filePathColorRelief);// qDebugV0() <<m_filePath<< "  createColorRelief() 代码执行时间: " << timer.elapsed() << " 毫秒";  // 计算执行时间（毫秒）if(!ret){return;}osg::ref_ptr<GDALImageLayer> layer = nullptr;if (!filePathColorRelief.empty()) {timer.restart();Global::buildOverviews(filePathColorRelief, -1); // 减少金字塔级数// qDebugV0() <<QString::fromStdString(filePathColorRelief)<< "  buildOverviewsTime() 代码执行时间: " << timer.elapsed() << " 毫秒";  // 计算执行时间（毫秒）osgEarth::GDALImageLayer::Options options;options.set_url(osgEarth::URI(filePathColorRelief));// 基础优化options.set_async( true );  // 异步加载layer = new GDALImageLayer(options);layer->options().set_name(m_filePath.toStdString());layer->setAsyncLoading(true);}emit sigProcessingFinished(layer);
}

C++ 里面嵌套着色器代码实现：


/******************************************************************************* 需求：* 高程渲染图 RGB.tif、 山体阴影图 AMP.tif** 高程渲染图 rgb波段分别  乘以 山体阴影图r波段， 然后除以255(AI说  读取的纹理就已经归一化到了 0~1 范围，不用除以 255)。*/void MyWidget::testGlsl()
{osgEarth::GDALImageLayer::Options ampOpt;// 2. AMP 阴影图层，开启 sharedampOpt.url() = "D:/Demo/glsl/AMP.tif";ampOpt.shared() = true;                   // 关键点：让它暴露共享采样器和矩阵ampOpt.shareTexUniformName() = "ampTex";        // Uniform 采样器名ampOpt.shareTexMatUniformName() = "ampTexMatrix";  // Uniform 矩阵名osg::ref_ptr<GDALImageLayer> ampLayer = new GDALImageLayer(ampOpt);ampLayer->setOpacity(0.00f);_map->addLayer(ampLayer.get());Global::buildOverviews("D:/Demo/glsl/RGB.tif", 20, "AVERAGE"); // 减少金字塔级数osgEarth::GDALImageLayer::Options rgbOpt;rgbOpt.tileSize() = 256;rgbOpt.maxLevel() = 20;rgbOpt.url() = "D:/Demo/glsl/RGB.tif";osg::ref_ptr<GDALImageLayer> rgbLayer = new GDALImageLayer(rgbOpt);_map->addLayer(rgbLayer.get());osg::ref_ptr<osg::StateSet> ss = rgbLayer->getOrCreateStateSet();ss->addUniform(new osg::Uniform("minLumFactor", 0.125f));//虚拟程序设置auto vp = osgEarth::VirtualProgram::getOrCreate(ss);// 获取可执行文件目录QString appDirQt = QCoreApplication::applicationDirPath();// GLSL 文件路径（同级 glsl 文件夹）QString glslPathQt = appDirQt + "/glsl/colorRelief_hillShade.glsl";std::string glslFile = glslPathQt.toStdString();// 打印检查std::cout << "glslFile: " << glslFile << std::endl;// 从文件里加载一个 GLSL 源码到 osg::Shader 对象osg::ref_ptr<osg::Shader> fragShader = osgDB::readRefShaderFile(osg::Shader::FRAGMENT, glslFile);if (fragShader.valid()) {std::string src = fragShader->getShaderSource();vp->setFunction("customFragment", src , osgEarth::ShaderComp::LOCATION_FRAGMENT_COLORING);}else{qDebugV5()<<"fragShader ×";}
}

以下是着色器glsl代码：

uniform sampler2D ampTex;
uniform mat4 ampTexMatrix;uniform float minLumFactor;// 新增 uniform，用滑条控制in vec2 oe_layer_texc;void customFragment(inout vec4 color)
{// RGB 原始颜色vec3 baseColor=color.rgb;// 转到 AMP 图层坐标vec2 uv_amp=(ampTexMatrix*vec4(oe_layer_texc,0.,1.)).st;// Hillshade 灰度float amp=texture(ampTex,uv_amp).r;// --- 关键：用 hillshade 调制亮度 ---// 1. 先算彩色图的亮度（感知加权公式）float lum=dot(baseColor,vec3(.299,.587,.114));// 2. 用 hillshade 去调制亮度// float lumMod = lum * amp;float lumMod=lum*(minLumFactor+amp*(1.-minLumFactor));// 把 hillshade 映射到 minLumFactor~1.0// 3. 保留原色相和饱和度：缩放 baseColor，使亮度变成 lumModfloat lumBase=max(lum,1e-4);// 防止除零vec3 resultColor=baseColor*(lumMod/lumBase);color.rgb=clamp(resultColor,0.,1.);
}

查看全文

http://www.dtcms.com/a/344925.html