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

这篇具有很好参考价值的文章主要介绍了基于C++、GDAL、OpenCV的矢量数据骨架线提取算法。希望对大家有所帮助。如果存在错误或未考虑完全的地方,请大家不吝赐教,您也可以点击"举报违法"按钮提交疑问。

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

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

效果:
基于C++、GDAL、OpenCV的矢量数据骨架线提取算法,计算机图形学,图像处理,数值运算,opencv,c++,算法,gdal

思路:
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转Mat
			void 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文章来源地址https://www.toymoban.com/news/detail-523244.html

        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 dataset
			GDALRasterBand* pBand = pDstDataset->GetRasterBand(1);

			// Allocate memory for the raster data
			int nXSize = pBand->GetXSize();
			int nYSize = pBand->GetYSize();
			GByte* pData = new GByte[nXSize * nYSize];

			// Read the raster data
			pBand->RasterIO(GF_Read, 0, 0, nXSize, nYSize, pData, nXSize, nYSize, GDT_Byte, 0, 0);
			//转化为Mat
			cv::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 drivers
			GDALAllRegister();
			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;
			}

			// 创建目标shapefile
			GDALDataset* 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;

			// 遍历所有源shapefile
			for (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;
			//nodata
			argv = 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 datasets
			GDALRasterizeOptionsFree(pOptions);

			//合并shp
			mergeShp(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;
		}

到了这里,关于基于C++、GDAL、OpenCV的矢量数据骨架线提取算法的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处: 如若内容造成侵权/违法违规/事实不符,请点击违法举报进行投诉反馈,一经查实,立即删除!

领支付宝红包 赞助服务器费用

相关文章

  • 创建大量栅格文件并分别写入像元数据:C++ GDAL代码实现

      本文介绍基于 C++ 语言 GDAL 库, 批量创建 大量栅格遥感影像文件,并将数据批量写入其中的方法。   首先,我们来明确一下本文所需实现的需求。已知我们对 大量遥感影像 进行了批量读取与数据处理操作——具体过程可以参考文章C++ GDAL提取多时相遥感影像中像素随

    2024年02月19日
    浏览(43)
  • c++ OpenCV——提取图像的局部区域

    有时候整幅图像需要采取局部,如何进行采取呢。方案如下: 确定裁剪区域的大下 将原图copy一份,原图备用,防止损坏 将copy的图进行裁剪 原图: copy图及需要裁剪的部分: 裁剪图:

    2024年02月09日
    浏览(43)
  • 基于图像的裂缝分割与裂缝宽度计算(正交骨架线法)

    From https://subce.gitee.io/htmls/essays/crack_width_calculation.html 可以使用 opencv 的阈值分割,或者使用 CNN 等深度学习的方法进行裂缝分割,一般得到的分割结果如下,这里不再赘述。 寻找边缘有两种方法,如下 scikit-image skimage.measure.find_contours, 官方的示例如下: 另一种可以使用 Dela

    2024年02月05日
    浏览(102)
  • 学习笔记:Opencv实现图像特征提取算法SIFT

    2023.8.19 为了在暑假内实现深度学习的进阶学习,特意学习一下传统算法,分享学习心得,记录学习日常 SIFT的百科: SIFT = Scale Invariant Feature Transform, 尺度不变特征转换 全网最详细SIFT算法原理实现_ssift算法_Tc.小浩的博客-CSDN博客 在环境配置中要配置opencv: pip install opencv-c

    2024年02月12日
    浏览(48)
  • opencv提取中心线坐标点算法(python)

    一、通过opencv方法提取骨线, 二、通过骨线,提取每个线的像素坐标 三、循环每个点根据距离连接 循环思路: 1、寻找一个点的八个方向,确定点存在数组里面,就放入线中 2、交点处,在交点位置就新建一个线,通过读取一种条线在接着往下寻找点,没有读取的线存在数组

    2024年04月09日
    浏览(42)
  • OpenCV图像特征提取学习五,HOG特征检测算法

    一、HOG向梯度直方图概述   向梯度直方图(Histogram of Oriented Gradient, HOG)特征是基于对稠密网格中归一化的局部方向梯度直方图的计算。此方法的基本观点是:局部目标的外表和形状可以被局部梯度或边缘方向的分布很好的描述,即使我们不知道对应的梯度和边缘的位置。在

    2024年02月04日
    浏览(48)
  • C++如何用OpenCV中实现图像的边缘检测和轮廓提取?

    主要实现代码:

    2024年02月07日
    浏览(40)
  • 基于 NCC/灰度信息 的模板匹配算法(QT + Opencv + C++),10ms内获取匹配结果,部分源码

    文后代码,优化效果图结尾处,最快3ms得到匹配结果 NCC,全称为Normalized Cross Correlation,即归一化互相关系数, 在模板匹配中使用的非常非常广泛,也是众多模板匹配方法中非常耀眼的存在, 这个匹配的理论核心基础公式如下: 其实Opencv的matchTemplate函数使用的就是这个公式

    2024年02月08日
    浏览(46)
  • C++ 图像线特征提取【HoughLinesP算法】

       HoughLinesP :是一种基于Hough变换的直线检测算法。它可以识别图像中的直线,并返回它们的端点坐标。其函数接口如下: cv::HoughLinesP(   InputArray src,   // 输入图像,必须 8-bit 的灰度图像   OutputArray lines,  // 输出的极坐标来表示直线   double rho,    // 生成极

    2024年02月07日
    浏览(40)
  • Lesson4-3:OpenCV图像特征提取与描述---SIFT/SURF算法

    学习目标 理解 S I F T / S U R F SIFT/SURF S I FT / S U RF 算法的原理, 能够使用 S I F T / S U R F SIFT/SURF S I FT / S U RF 进行关键点的检测 1.1 SIFT原理 前面两节我们介绍了 H a r r i s Harris H a rr i s 和 S h i − T o m a s i Shi-Tomasi S hi − T o ma s i 角点检测算法,这两种算法具有旋转不变性,但不具

    2024年02月09日
    浏览(52)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包,优惠每天领

二维码1

领取红包

二维码2

领红包