基于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转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

        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