diff --git a/CloudTools.Vegetation.Verify/main.cpp b/CloudTools.Vegetation.Verify/main.cpp index 72a7b31..61237d0 100644 --- a/CloudTools.Vegetation.Verify/main.cpp +++ b/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())) << "%" < +#include +#include +#include +#include +#include + + +namespace CloudTools::DEM +{ +/// +/// Represents the removal of false positive seed points close to buildings +/// +template +class BuildingFacadeSeedRemoval : public DatasetCalculation +{ + std::vector &seedPoints; +public: + int threshold; + /// + /// Initializes a new instance of the class. Loads input metadata and defines calculation. + /// + /// The source files of the difference comparison. + /// The callback method to report progress. + BuildingFacadeSeedRemoval(std::vector& seedPoints, + const std::vector& sourcePaths, + Operation::ProgressType progress = nullptr, + int threshold = 20) + : DatasetCalculation(sourcePaths, nullptr, progress), + seedPoints(seedPoints), + threshold(threshold) + { + initialize(); + } + + BuildingFacadeSeedRemoval(const BuildingFacadeSeedRemoval&) = delete; + BuildingFacadeSeedRemoval& operator=(const BuildingFacadeSeedRemoval&) = delete; + +private: + /// + /// Initializes the new instance of the class. + /// + void initialize(); +}; + +template +void BuildingFacadeSeedRemoval::initialize() +{ + this->computation = [this](int x, int y) + { + std::vector 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 diff --git a/CloudTools.Vegetation/CMakeLists.txt b/CloudTools.Vegetation/CMakeLists.txt index 4170baf..fa5ea5c 100644 --- a/CloudTools.Vegetation/CMakeLists.txt +++ b/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}) diff --git a/CloudTools.Vegetation/PreProcess.cpp b/CloudTools.Vegetation/PreProcess.cpp index 08aab90..c1639ca 100644 --- a/CloudTools.Vegetation/PreProcess.cpp +++ b/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 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 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(GDALOpen(_dsmInputPath.c_str(), GA_ReadOnly)); + Difference 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 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 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 seedRemoval(seedPoints, {_dtmInputPath, _dsmInputPath}, _progress); + seedRemoval.execute(); + } _progressMessage = "Tree crown segmentation (" + _prefix + ")"; { diff --git a/CloudTools.Vegetation/PreProcess.h b/CloudTools.Vegetation/PreProcess.h index 88fefc9..e132947 100644 --- a/CloudTools.Vegetation/PreProcess.h +++ b/CloudTools.Vegetation/PreProcess.h @@ -12,6 +12,11 @@ namespace Vegetation class PreProcess : public CloudTools::Operation, protected CloudTools::IO::ResultCollection { public: + enum ProcessingMethod + { + Standard, + SeedRemoval + }; /// /// Callback function for reporting progress. /// @@ -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; diff --git a/CloudTools.Vegetation/RiverMask.hpp b/CloudTools.Vegetation/RiverMask.hpp new file mode 100644 index 0000000..d50345c --- /dev/null +++ b/CloudTools.Vegetation/RiverMask.hpp @@ -0,0 +1,60 @@ +#pragma once + +#include +#include +#include +#include + +#include + +namespace CloudTools +{ +namespace DEM +{ +/// +/// Represents a DTM in which rivers are not nodata points, but instead have a constant value of -0.4 +/// +template +class RiverMask : public SweepLineTransformation +{ +public: + RiverMask(const std::vector& sourcePaths, + const std::string& targetPath, + Operation::ProgressType progress = nullptr) + : SweepLineTransformation(sourcePaths, targetPath, nullptr, progress) + { + initialize(); + } + + RiverMask(const RiverMask&) = delete; + RiverMask& operator=(const RiverMask&) = delete; + +private: + /// + /// Initializes the new instance of the class. + /// + void initialize(); + /// + /// Extremal low value for rivers. + /// + const int riverHeight = -1000; +}; + +template +void RiverMask::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>& sources) + { + if (!sources[0].hasData() && sources[1].hasData()) + return static_cast(this->nodataValue); + if (!sources[1].hasData()) + return static_cast(riverHeight); + if (sources[0].hasData()) + return static_cast(sources[0].data()); + }; +} +} // DEM +} // CloudTools diff --git a/CloudTools.Vegetation/main.cpp b/CloudTools.Vegetation/main.cpp index 4c3b7ee..27cffbe 100644 --- a/CloudTools.Vegetation/main.cpp +++ b/CloudTools.Vegetation/main.cpp @@ -38,6 +38,7 @@ int main(int argc, char* argv[]) ("dtm-input-path-B,t", po::value(&dtmInputPathB), "Epoch-B DTM input path") ("output-dir,o", po::value(&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");