基于C++、GDAL、OpenCV的矢量数据骨架线提取算法

基于C++、GDAL、OpenCV的矢量数据骨架线提取算法

CGAL已经实现了该功能，但由于CGAL依赖于Boost库，编译后过大，因此本文所采用的这套方式实现骨架线提取功能。

效果：
骨架线提取效果

思路：
1、将导入shp按照要素逐一拆分成新的shp
2、将所有拆分后的shp分别转栅格，利用OpenCV提取骨架线
3、将所有骨架线转为shp，并合并输出

详细代码如下：

调用

    basePolygonAlgorithm::SkeletonExtractor extract2; extract2.polygon2Skelton("originFile.shp", "outputFile.shp");

.h

#include "opencv2/core/core.hpp"
#include "opencv2/opencv.hpp"
#include"opencv2/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp" 
#include "ogrsf_frmts.h"
#include "gdal_priv.h"
#include "ogr_geometry.h"
#include "ogr_attrind.h"
#include "ogr_srs_api.h" //提取骨架线class SkeletonExtractor{public:SkeletonExtractor();void polygon2Skelton(string src_path, string dst_path);//src_path待提取数据  dst_path：提取后数据;private:cv::Mat thinning(cv::Mat img);//Mat骨架线提取cv::Mat readRasterData(GDALDataset* pDstDataset);//pDstDataset转Matvoid polygon2subPolygon(const string& src_path, int* featureNum, double adfGeoTransform[6]);void subPolygon2skelton(const string& dst_path, const int& featureNum, double adfGeoTransform[6]);void mergeShp(vector<string>vecSrcFiles, string pszDstFile);void writeShp(string dst_path, cv::Mat skel, double adfGeoTransform[], OGRSpatialReference* oSRS);};

.cpp

        SkeletonExtractor::SkeletonExtractor(){}void SkeletonExtractor::polygon2Skelton(string src_path, string dst_path){int featureNum = 0;double adfGeoTransform[6];polygon2subPolygon(src_path, &featureNum, adfGeoTransform);subPolygon2skelton(dst_path, featureNum, adfGeoTransform);return;}cv::Mat SkeletonExtractor::thinning(cv::Mat img){cv::Mat skel(img.size(), CV_8UC1, cv::Scalar(0));cv::Mat temp(img.size(), CV_8UC1);cv::Mat element = cv::getStructuringElement(cv::MORPH_CROSS, cv::Size(3, 3));bool done;do{cv::morphologyEx(img, temp, cv::MORPH_OPEN, element);cv::bitwise_not(temp, temp);cv::bitwise_and(img, temp, temp);cv::bitwise_or(skel, temp, skel);cv::erode(img, img, element);double maxVal = 0;cv::minMaxLoc(img, nullptr, &maxVal);done = (maxVal == 0);} while (!done);img.release();temp.release();return skel;  }cv::Mat SkeletonExtractor::readRasterData(GDALDataset* pDstDataset){// Read the first band of the raster datasetGDALRasterBand* pBand = pDstDataset->GetRasterBand(1);// Allocate memory for the raster dataint nXSize = pBand->GetXSize();int nYSize = pBand->GetYSize();GByte* pData = new GByte[nXSize * nYSize];// Read the raster datapBand->RasterIO(GF_Read, 0, 0, nXSize, nYSize, pData, nXSize, nYSize, GDT_Byte, 0, 0);//转化为Matcv::Mat _Mymat(nYSize, nXSize, CV_8UC1);for (int i = 0; i < nXSize; i++){for (int j = 0; j < nYSize; j++){_Mymat.at<uchar>(j, i) = (uchar)pData[j * nXSize + i];}}// 执行骨架化cv::Mat skel;skel = thinning(_Mymat);// 定义结构元素（如3x3的矩形）int morph_size = 1;cv::Mat element = cv::getStructuringElement(cv::MORPH_RECT,cv::Size(2 * morph_size + 1, 2 * morph_size + 1),cv::Point(morph_size, morph_size));// 膨胀操作cv::dilate(skel, skel, element);// 腐蚀操作cv::erode(skel, skel, element);delete[] pData;return skel;}void SkeletonExtractor::polygon2subPolygon(const string& src_path, int* featureNum, double adfGeoTransform[6]){// Register GDAL driversGDALAllRegister();OGRRegisterAll();//	NULL, NULL, NULL)); GDALDataset* pSrcDataset = (GDALDataset*)GDALOpenEx(src_path.c_str(), GDAL_OF_VECTOR, NULL, NULL, NULL);// The one layer.OGRLayer* pSrcLayer = pSrcDataset->GetLayer(0);OGRSpatialReference* oSRS = NULL;oSRS = new OGRSpatialReference;oSRS = pSrcLayer->GetSpatialRef();*featureNum = pSrcLayer->GetFeatureCount();pSrcDataset->GetGeoTransform(adfGeoTransform);for (int i = 0; i < *featureNum; i++) {// 获取要素OGRFeature* pSrcFeature = pSrcLayer->GetFeature(i);// 创建新的GDALDataset和图层char outputFilename[256];sprintf(outputFilename, "temp_%d.shp", i);GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("ESRI Shapefile");GDALDataset* temp_pDstDataset = pDriver->Create(outputFilename, 0, 0, 0, GDT_Unknown, NULL);OGRLayer* pDstLayer = temp_pDstDataset->CreateLayer(pSrcLayer->GetName(), oSRS, pSrcLayer->GetGeomType(), NULL);// 复制要素到新的图层中OGRFeature* pDstFeature = pSrcFeature->Clone();pDstLayer->CreateFeature(pDstFeature);OGRFeature::DestroyFeature(pDstFeature);GDALFlushCache(temp_pDstDataset);// 释放资源GDALClose(temp_pDstDataset);OGRFeature::DestroyFeature(pSrcFeature);}GDALClose(pSrcDataset); }void SkeletonExtractor::mergeShp(vector<string> vecSrcFiles, string pszDstFile){GDALAllRegister();// 获取Shapefile驱动GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ESRI Shapefile");if (poDriver == nullptr) {cerr << "ESRI Shapefile driver not available." << endl;return;}// 创建目标shapefileGDALDataset* poDstDS = poDriver->Create(pszDstFile.c_str(), 0, 0, 0, GDT_Unknown, nullptr);if (poDstDS == nullptr) {cerr << "Creation of output file failed." << endl;return;}OGRLayer* poDstLayer = nullptr;// 遍历所有源shapefilefor (const auto& srcFile : vecSrcFiles) {GDALDataset* poSrcDS = static_cast<GDALDataset*>(GDALOpenEx(srcFile.c_str(), GDAL_OF_VECTOR, nullptr, nullptr, nullptr));if (poSrcDS == nullptr) {cerr << "Open failed: " << srcFile << endl;continue;}OGRLayer* poSrcLayer = poSrcDS->GetLayer(0);if (poSrcLayer == nullptr) {cerr << "Failed to get layer from source file: " << srcFile << endl;GDALClose(poSrcDS);continue;}// 如果目标图层不存在，则根据源图层创建一个if (poDstLayer == nullptr) {poDstLayer = poDstDS->CreateLayer(poSrcLayer->GetName(), poSrcLayer->GetSpatialRef(), poSrcLayer->GetGeomType(), nullptr);if (poDstLayer == nullptr) {cerr << "Failed to create destination layer." << endl;GDALClose(poSrcDS);GDALClose(poDstDS);return;}// 复制属性表结构OGRFeatureDefn* poSrcFDefn = poSrcLayer->GetLayerDefn();for (int i = 0; i < poSrcFDefn->GetFieldCount(); ++i) {poDstLayer->CreateField(poSrcFDefn->GetFieldDefn(i));}}// 复制要素OGRFeature* poFeature;poSrcLayer->ResetReading();while ((poFeature = poSrcLayer->GetNextFeature()) != nullptr) {OGRFeature* poDstFeature = OGRFeature::CreateFeature(poDstLayer->GetLayerDefn());poDstFeature->SetFrom(poFeature);poDstFeature->SetGeometry(poFeature->GetGeometryRef());if (poDstLayer->CreateFeature(poDstFeature) != OGRERR_NONE) {cerr << "Failed to create feature in destination shapefile." << endl;OGRFeature::DestroyFeature(poFeature);OGRFeature::DestroyFeature(poDstFeature);GDALClose(poSrcDS);GDALClose(poDstDS);return;}OGRFeature::DestroyFeature(poFeature);OGRFeature::DestroyFeature(poDstFeature);}GDALClose(poSrcDS);}GDALClose(poDstDS);}void SkeletonExtractor::writeShp(string dst_path, cv::Mat skel, double adfGeoTransform[], OGRSpatialReference* oSRS){GDALAllRegister();OGRRegisterAll();GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ESRI Shapefile");GDALDataset* poDS = poDriver->Create(dst_path.c_str(), 0, 0, 0, GDT_Unknown, NULL);OGRLayer* poLayer = poDS->CreateLayer("layer", oSRS, wkbUnknown, NULL);std::vector<cv::Vec4i> hierarchy;std::vector<std::vector<cv::Point>> contours;cv::findContours(skel, contours, hierarchy, cv::RETR_EXTERNAL, cv::CHAIN_APPROX_NONE);OGRMultiLineString multiLineString;OGRFeature* poFeature = OGRFeature::CreateFeature(poLayer->GetLayerDefn());for (int i = 0; i < contours.size(); i++){OGRLineString lineString;for (int j = 0; j < contours[i].size(); j++){double x = adfGeoTransform[0] + contours[i][j].x * adfGeoTransform[1] + contours[i][j].y * adfGeoTransform[2];double y = adfGeoTransform[3] + contours[i][j].x * adfGeoTransform[4] + contours[i][j].y * adfGeoTransform[5];lineString.addPoint(x, y);}multiLineString.addGeometry(&lineString);}poFeature->SetGeometry(&multiLineString);skelton_result = poFeature->GetGeometryRef()->clone();//clone将skelton_result复制到一个新的 OGRGeometry 对象中，防止因poFeature释放而变成空poLayer->CreateFeature(poFeature);OGRFeature::DestroyFeature(poFeature);GDALClose(poDS);hierarchy.clear();contours.clear(); }void SkeletonExtractor::subPolygon2skelton(const string& dst_path, const int& featureNum, double adfGeoTransform[6]){GDALAllRegister();OGRRegisterAll();//栅格参数设置char** argv = NULL;//nodataargv = CSLAddString(argv, "-a_nodata");argv = CSLAddString(argv, "-9999");//初始值argv = CSLAddString(argv, "-init");argv = CSLAddString(argv, "0");//分辨率argv = CSLAddString(argv, "-ts");argv = CSLAddString(argv, "1024");argv = CSLAddString(argv, "1024");//存储类型argv = CSLAddString(argv, "-ot");argv = CSLAddString(argv, "Byte");GDALRasterizeOptions* pOptions = GDALRasterizeOptionsNew(argv, NULL);OGRSpatialReference* oSRS = NULL;vector<string>filePaths;for (int i = 0; i < featureNum; i++) { //重新打开新的SHP文件char outputFilename[256];sprintf(outputFilename, "temp_%d.shp", i);GDALDataset* pNewDS = (GDALDataset*)GDALOpenEx(outputFilename, GDAL_OF_VECTOR, NULL, NULL, NULL);OGRLayer* pNewLayer = pNewDS->GetLayer(0);oSRS = new OGRSpatialReference;oSRS = pNewLayer->GetSpatialRef();//获取转换6参数int pbUsageError = FALSE;GDALDataset* pImageDataset = static_cast<GDALDataset*>(GDALRasterize("output.tif", NULL, pNewDS, pOptions, &pbUsageError));double adfGeoTransform[6];pImageDataset->GetGeoTransform(adfGeoTransform);//提取骨架线cv::Mat skel = readRasterData(pImageDataset);//写入结果outputFilename[256];sprintf(outputFilename, "res_%d.shp", i);writeShp(outputFilename, skel, adfGeoTransform, oSRS);filePaths.push_back(outputFilename);GDALClose(pNewDS);GDALClose(pImageDataset);skel.release();}// Close the datasetsGDALRasterizeOptionsFree(pOptions);//合并shpmergeShp(filePaths, dst_path);//删除临时文件for (int i = 0; i < featureNum; i++) {char tempFile[256];char resFile[256];sprintf(tempFile, "temp_%d.shp", i);std::remove(tempFile);sprintf(tempFile, "temp_%d.shx", i);std::remove(tempFile);sprintf(tempFile, "temp_%d.dbf", i);std::remove(tempFile);sprintf(tempFile, "temp_%d.prj", i);std::remove(tempFile);sprintf(resFile, "res_%d.shp", i);std::remove(resFile);sprintf(resFile, "res_%d.shx", i);std::remove(resFile);sprintf(resFile, "res_%d.dbf", i);std::remove(resFile);sprintf(resFile, "res_%d.prj", i);std::remove(resFile);}std::remove("output.tif");return;}

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.hqwc.cn/news/14891.html

如若内容造成侵权/违法违规/事实不符，请联系编程知识网进行投诉反馈email:809451989@qq.com，一经查实，立即删除！

相关文章

类Twitter风格的RSS阅读器

类Twitter风格的RSS阅读器

本文完成于 2 月中旬&#xff0c;其中的反代还是 frp npm 方案&#xff1b; 什么是 RSS ? RSS 是用 PHP、Laravel、Inertia.js、Tailwind 和 Vue.js 编写的简单的类Twitter 风格的 RSS阅读器&#xff0c;支持 RSS和ATOM 格式。 命令行安装 在群晖上以 Docker 方式安装。 官…
阅读更多...
Nacos2.3.0源码启动报错找不到符号com.alibaba.nacos.consistency.entity

Nacos2.3.0源码启动报错找不到符号com.alibaba.nacos.consistency.entity

一. 源码下载编译&#xff1a;找不到符号com.alibaba.nacos.consistency.entity 如果报错找不到符号com.alibaba.nacos.consistency.entity Nacos\consistency\src\main\java\com\alibaba\nacos\consistency\entity 这个包下没有相关的java文件&#xff0c;其实是我们没有编译…
阅读更多...
hadoop --- MapReduce

hadoop --- MapReduce

MapReduce定义&#xff1a; MapReduce可以分解为Map (映射) Reduce (规约) &#xff0c; 具体过程&#xff1a; Map : 输入数据集被切分成多个小块&#xff0c;并分配给不同的计算节点进行处理Shuffle and Sort&#xff1a;洗牌和排序&#xff0c;在 Map 阶段结束后&#xf…
阅读更多...
15.1 BP神经网络实现图像压缩——了解神经网络在图像处理方面的应用（matlab程序）

15.1 BP神经网络实现图像压缩——了解神经网络在图像处理方面的应用（matlab程序）

1.简述 BP神经网络现在来说是一种比较成熟的网络模型了,因为神经网络对于数字图像处理的先天优势,特别是在图像压缩方面更具有先天的优势,因此,我这一段时间在研究神经网络的时候同时研究了一下关于BP网络实现图像压缩的原理和过程,并且是在MATLAB上进行了仿真的实验,结果发现设…
阅读更多...
TinyStories: How Small Can Language Models Be and Still Speak Coherent English?

TinyStories: How Small Can Language Models Be and Still Speak Coherent English?

本文是LLM系列的文章之一&#xff0c;针对《TinyStories: How Small Can Language Models Be and Still Speak Coherent English?》的翻译。 TinyStories&#xff1a;语言模型能有多小&#xff0c;还能说连贯的英语&#xff1f; 摘要1 引言2 TinyStories数据集的描述2.1 Tiny…
阅读更多...
3D模型轻量化开发工具HOOPS与WebGL的对比分析

3D模型轻量化开发工具HOOPS与WebGL的对比分析

HOOPS是一种商业级的3D开发平台&#xff0c;由Tech Soft 3D公司开发。它提供了一套全面的工具和API&#xff0c;用于构建和展示高度复杂的3D场景和模型&#xff0c;可以在多个平台和环境中使用&#xff0c;包括Web、移动设备和桌面&#xff0c;这使得开发者能够在不同的设备上展…
阅读更多...
UE5接入在线直播视频源，如hls(m3u8)格式

UE5接入在线直播视频源，如hls(m3u8)格式

文章目录 1.实现目标2.实现过程2.1 VlcMedia插件重编译2.2 UE5接入在线直播2.3 创建材质3.参考资料1.实现目标 通过重编译VlcMedia插件，以支持在线直播视频在UE5中的播放，GIF动图如下： 2.实现过程 本文主要包括插件的重编译、在线直播视频的接入，以及材质的创建三个部分。…
阅读更多...
ELK部署安装

ELK部署安装

目录 一、环境准备 1.准备三台服务器&#xff08;带图形化的linuxCentOS7&#xff0c;最小化缺少很多环境&#xff09; 2.修改主机名 3.关闭防火墙 4.elk-node1、elk-node2 用系统自带的java 5.上传软件包到node1和node2 二、部署elasticsearch 1、node1、node2操作 2.no…
阅读更多...
Coggle 30 Days of ML（23年7月）任务四：线性模型训练与预测

Coggle 30 Days of ML（23年7月）任务四：线性模型训练与预测

Coggle 30 Days of ML&#xff08;23年7月&#xff09;任务四&#xff1a;线性模型训练与预测 任务四&#xff1a;使用TFIDF特征和线性模型完成训练和预测 说明&#xff1a;在这个任务中&#xff0c;你需要使用TFIDF特征和线性模型&#xff08;如逻辑回归&#xff09;完成训练…
阅读更多...
Jmeter做单接口测试-超详细步骤讲解

Jmeter做单接口测试-超详细步骤讲解

测试项目&#xff1a;本章节将以此测试项目为大家讲解怎么使用jmeter做一个接口测试 CSDN - 专业开发者社区CSDN是全球知名中文IT技术交流平台,创建于1999年,包含原创博客、精品问答、职业培训、技术论坛、资源下载等产品服务,提供原创、优质、完整内容的专业IT技术开发社区.h…
阅读更多...
基于梯度下降算法的无约束函数极值问题求解

基于梯度下降算法的无约束函数极值问题求解

基于梯度下降算法的无约束函数极值问题求解 1 知识预警1.1导数1.2偏导数1.3方向导数1.4梯度 2 梯度下降算法3 无约束函数极值问题求解3.1 算例13.1.1 Python编程求解3.1.2 求解结果与可视化 3.2 算例2 Rosenbrock函数3.2.1 Python编程求解3.2.2 求解结果与可视化 1 知识预警 1…
阅读更多...
多元分类预测 | Matlab 麻雀算法(SSA)优化xgboost的分类预测，多特征输入模型，SSA-xgboost分类预测模型

多元分类预测 | Matlab 麻雀算法(SSA)优化xgboost的分类预测，多特征输入模型，SSA-xgboost分类预测模型

文章目录 效果一览文章概述部分源码参考资料效果一览 文章概述 多元分类预测 | Matlab 麻雀算法(SSA)优化xgboost的分类预测，多特征输入模型，SSA-xgboost分类预测模型 多特征输入单输出的二分类及多分类模型。程序
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023