Merge branch 'vegetation-mix' into 'master'

Vegetation false positives reduction

See merge request gislab/cloudtools!10

CloudTools.Vegetation.Verify/main.cpp

@@ -306,9 +306,12 @@ int main(int argc, char* argv[]) try
 	          << "Verification completed!" << std::endl
 	          << std::endl << std::fixed << std::setprecision(2)
 	          << "[Basic statistic]" << std::endl
+			  << "Detected trees: " << inputTrees.size() << std::endl
 	          << "Reference trees matched: " << matched.size() << std::endl
 	          << "Reference trees failed: " << missed.size() << std::endl
 	          << "Match ratio: " << ((100.f * matched.size()) / (referenceTrees.size())) << "%" << std::endl
+			  << "False Positives: " << inputTrees.size() - matched.size()
+			  << "  " << ((100.f * (inputTrees.size()-matched.size())) / (inputTrees.size())) << "%" <<std::endl
 	          << std::endl
 	          << "[Miss statistic]" << std::endl
 	          << "Average radius: " << ((std::accumulate(missed.begin(), missed.end(),

CloudTools.Vegetation/BuildingFacadeSeedRemoval.hpp

@@ -0,0 +1,78 @@
+#pragma once
+
+#include <string>
+#include <vector>
+#include <cmath>
+#include <iostream>
+#include <CloudTools.DEM/SweepLineCalculation.hpp>
+#include <CloudTools.DEM/DatasetCalculation.hpp>
+
+
+namespace CloudTools::DEM
+{
+/// <summary>
+/// Represents the removal of false positive seed points close to buildings
+/// </summary>
+template <typename DataType = float>
+class BuildingFacadeSeedRemoval : public DatasetCalculation<DataType>
+{
+	std::vector<OGRPoint> &seedPoints;
+public:
+	int threshold;
+	/// <summary>
+	/// Initializes a new instance of the class. Loads input metadata and defines calculation.
+	/// </summary>
+	/// <param name="sourcePaths">The source files of the difference comparison.</param>
+	/// <param name="progress">The callback method to report progress.</param>
+	BuildingFacadeSeedRemoval(std::vector<OGRPoint>& seedPoints,
+	                          const std::vector<std::string>& sourcePaths,
+	                          Operation::ProgressType progress = nullptr,
+							  int threshold = 20)
+		: DatasetCalculation<DataType>(sourcePaths, nullptr, progress),
+		    seedPoints(seedPoints),
+			threshold(threshold)
+	{
+		initialize();
+	}
+
+	BuildingFacadeSeedRemoval(const BuildingFacadeSeedRemoval&) = delete;
+	BuildingFacadeSeedRemoval& operator=(const BuildingFacadeSeedRemoval&) = delete;
+
+private:
+	/// <summary>
+	/// Initializes the new instance of the class.
+	/// </summary>
+	void initialize();
+};
+
+template <typename DataType>
+void BuildingFacadeSeedRemoval<DataType>::initialize()
+{
+	this->computation = [this](int x, int y)
+	{
+		std::vector<int> idxs;
+		int c = 0;
+		for (auto& point: seedPoints)
+		{
+			int counter = 0;
+			int ws = 3; // window size = ws*2 + 1
+			int px = point.getX();
+			int py = point.getY();
+			for(int i = px - ws; i <= px + ws; i++){
+				for(int j = py - ws; j <= py + ws; j++){
+					if(!this->hasSourceData(0,i,j) && this->sourceData(1,i,j) > 10){
+						counter++;
+					}
+				}
+			}
+			if(counter > threshold){
+				idxs.push_back(c);
+			}
+			c++;
+		}
+		for (int i = idxs.size() - 1; i >= 0; --i) {
+			seedPoints.erase(seedPoints.begin() + idxs[i]);
+		}
+	};
+}
+} // CloudTools

CloudTools.Vegetation/CMakeLists.txt

@@ -13,13 +13,20 @@ add_executable(vegetation
 	DistanceCalculation.h
 	HeightDifference.h HeightDifference.cpp
 	PreProcess.cpp PreProcess.h
-	PostProcess.cpp PostProcess.h)
+	PostProcess.cpp PostProcess.h
+		BuildingFacadeSeedRemoval.hpp
+	RiverMask.hpp)
 
 target_link_libraries(vegetation
+	PRIVATE ${PDAL_LIBRARIES}
 	dem
 	common
 	Threads::Threads)
 
+target_include_directories(vegetation PRIVATE
+		${PDAL_INCLUDE_DIRS}
+		${PDAL_INCLUDE_DIRS}/pdal)
+
 install(TARGETS vegetation
 	DESTINATION ${CMAKE_INSTALL_PREFIX})
 

CloudTools.Vegetation/PreProcess.cpp

@@ -13,6 +13,8 @@
 #include "InterpolateNoData.h"
 #include "TreeCrownSegmentation.h"
 #include "MorphologyClusterFilter.h"
+#include "RiverMask.hpp"
+#include "BuildingFacadeSeedRemoval.hpp"
 
 using namespace CloudTools::DEM;
 using namespace CloudTools::IO;
@@ -38,13 +40,36 @@ void PreProcess::onPrepare()
 
 void PreProcess::onExecute()
 {
-	_progressMessage = "Creating CHM (" + _prefix + ")";
-	newResult("CHM");
+	if(_processingMethod == PreProcess::SeedRemoval)
 	{
-		Difference<float> comparison({_dtmInputPath, _dsmInputPath}, result("CHM").path(), _progress);
-		comparison.execute();
-		result("CHM").dataset = comparison.target();
-		_targetMetadata = comparison.targetMetadata();
+		_progressMessage = "Creating River Map (" + _prefix + ")";
+		newResult("RM");
+		{
+			RiverMask<float> riverMap({_dtmInputPath, _dsmInputPath}, result("RM").path(), _progress);
+			riverMap.execute();
+			result("RM").dataset = riverMap.target();
+			_targetMetadata = riverMap.targetMetadata();
+		}
+
+		_progressMessage = "Creating CHM (" + _prefix + ")";
+		newResult("CHM");
+		{
+			auto dsm = static_cast<GDALDataset*>(GDALOpen(_dsmInputPath.c_str(), GA_ReadOnly));
+			Difference<float> comparison({{result("RM").dataset},
+			                              {dsm}}, result("CHM").path(), _progress);
+			comparison.execute();
+			result("CHM").dataset = comparison.target();
+			_targetMetadata = comparison.targetMetadata();
+		}
+	} else {
+		_progressMessage = "Creating CHM (" + _prefix + ")";
+		newResult("CHM");
+		{
+			Difference<float> comparison({_dtmInputPath, _dsmInputPath}, result("CHM").path(), _progress);
+			comparison.execute();
+			result("CHM").dataset = comparison.target();
+			_targetMetadata = comparison.targetMetadata();
+		}
 	}
 
 	_progressMessage = "Matrix transformation (" + _prefix + ")";
@@ -76,7 +101,16 @@ void PreProcess::onExecute()
 	_progressMessage = "Seed points collection (" + _prefix + ")";
 	std::vector<OGRPoint> seedPoints = collectSeedPoints(result("interpol").dataset);
 	if (debug)
+	{
 		writePointsToFile(seedPoints, (fs::path(_outputDir) / (_prefix + "_seedpoints.json")).string());
+	}
+
+	if(_processingMethod == PreProcess::SeedRemoval)
+	{
+		_progressMessage = "Seed Removal(" + _prefix + ")";
+		::BuildingFacadeSeedRemoval<float> seedRemoval(seedPoints, {_dtmInputPath, _dsmInputPath}, _progress);
+		seedRemoval.execute();
+	}
 
 	_progressMessage = "Tree crown segmentation (" + _prefix + ")";
 	{

CloudTools.Vegetation/PreProcess.h

@@ -12,6 +12,11 @@ namespace Vegetation
 class PreProcess : public CloudTools::Operation, protected CloudTools::IO::ResultCollection
 {
 public:
+	enum ProcessingMethod
+	{
+		Standard,
+		SeedRemoval
+	};
 	/// <summary>
 	/// Callback function for reporting progress.
 	/// </summary>
@@ -51,11 +56,13 @@ class PreProcess : public CloudTools::Operation, protected CloudTools::IO::Resul
 	PreProcess(const std::string& prefix,
 	           const std::string& dtmInputPath,
 	           const std::string& dsmInputPath,
-	           const std::string& outputDir)
+	           const std::string& outputDir,
+			   const ProcessingMethod processingMethod = ProcessingMethod::Standard)
 		: _prefix(prefix),
 		  _dtmInputPath(dtmInputPath),
 		  _dsmInputPath(dsmInputPath),
-		  _outputDir(outputDir)
+		  _outputDir(outputDir),
+		  _processingMethod(processingMethod)
 	{
 	}
 
@@ -93,6 +100,8 @@ class PreProcess : public CloudTools::Operation, protected CloudTools::IO::Resul
 	CloudTools::IO::Result* createResult(const std::string& name, bool isFinal = false) override;
 
 private:
+	ProcessingMethod _processingMethod;
+
 	std::string _prefix, _dtmInputPath, _dsmInputPath, _outputDir;
 	CloudTools::DEM::RasterMetadata _targetMetadata;
 	CloudTools::DEM::ClusterMap _targetCluster;

CloudTools.Vegetation/RiverMask.hpp

@@ -0,0 +1,60 @@
+#pragma once
+
+#include <string>
+#include <vector>
+#include <cmath>
+#include <iostream>
+
+#include <CloudTools.DEM/SweepLineTransformation.hpp>
+
+namespace CloudTools
+{
+namespace DEM
+{
+/// <summary>
+/// Represents a DTM in which rivers are not nodata points, but instead have a constant value of -0.4
+/// </summary>
+template <typename DataType = float>
+class RiverMask : public SweepLineTransformation<DataType>
+{
+public:
+	RiverMask(const std::vector<std::string>& sourcePaths,
+	          const std::string& targetPath,
+	          Operation::ProgressType progress = nullptr)
+		: SweepLineTransformation<DataType>(sourcePaths, targetPath, nullptr, progress)
+	{
+		initialize();
+	}
+
+	RiverMask(const RiverMask&) = delete;
+	RiverMask& operator=(const RiverMask&) = delete;
+
+private:
+	/// <summary>
+	/// Initializes the new instance of the class.
+	/// </summary>
+	void initialize();
+        /// <summary>
+	/// Extremal low value for rivers.
+	/// </summary>
+	const int riverHeight = -1000;
+};
+
+template <typename DataType>
+void RiverMask<DataType>::initialize()
+{
+	/// returns the DTM value, if both the DTM and the DSM value is given;
+	/// returns -0.4 if only the DTM value or neither the DTM and the DSM value is given;
+	/// return the DSM value, if only the DSM value is given.
+	this->computation = [this](int x, int y, const std::vector<Window<DataType>>& sources)
+	{
+		if (!sources[0].hasData() && sources[1].hasData())
+			return static_cast<DataType>(this->nodataValue);
+		if (!sources[1].hasData())
+			return static_cast<DataType>(riverHeight);
+		if (sources[0].hasData())
+			return static_cast<DataType>(sources[0].data());
+	};
+}
+} // DEM
+} // CloudTools

CloudTools.Vegetation/main.cpp

@@ -38,6 +38,7 @@ int main(int argc, char* argv[])
 		("dtm-input-path-B,t", po::value<std::string>(&dtmInputPathB), "Epoch-B DTM input path")
 		("output-dir,o", po::value<std::string>(&outputDir)->default_value(outputDir), "result directory path")
 		("hausdorff-distance", "use Hausdorff-distance")
+		("srm", "removes trees possibly to close to buildings")
 		("parallel,p", "parallel execution for A & B epochs")
 		("debug,d", "keep intermediate results on disk after progress")
 		("verbose,v", "verbose output")
@@ -136,8 +137,14 @@ int main(int argc, char* argv[])
 	};
 
 	// Create preprocessors
-	PreProcess preProcessA("a", dtmInputPathA, dsmInputPathA, outputDir);
-	PreProcess preProcessB("b", dtmInputPathB, dsmInputPathB, outputDir);
+	PreProcess preProcessA("a", dtmInputPathA, dsmInputPathA, outputDir,
+						   vm.count("srm")
+						   ? PreProcess::ProcessingMethod::SeedRemoval
+						   : PreProcess::ProcessingMethod::Standard);
+	PreProcess preProcessB("b", dtmInputPathB, dsmInputPathB, outputDir,
+						   vm.count("srm")
+						   ? PreProcess::ProcessingMethod::SeedRemoval
+						   : PreProcess::ProcessingMethod::Standard);
 
 	preProcessA.debug = vm.count("debug");
 	preProcessB.debug = vm.count("debug");