From 0c5eb9add5fd3c84d5c43b183b13687433d2f0ad Mon Sep 17 00:00:00 2001 From: Pedro Maciel Date: Wed, 12 Jun 2024 09:19:13 +0100 Subject: [PATCH] interpolation=voronoi-average --- CMakeLists.txt | 16 ++ src/tools/CMakeLists.txt | 8 + src/tools/mir-voronoi-average.cc | 262 +++++++++++++++++++++++++++++++ 3 files changed, 286 insertions(+) create mode 100644 src/tools/mir-voronoi-average.cc diff --git a/CMakeLists.txt b/CMakeLists.txt index 4498386df..5486aef3a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -34,6 +34,22 @@ ecbuild_find_package(NAME eckit VERSION 1.22 REQUIRED) ecbuild_add_option(FEATURE ATLAS DEFAULT ON DESCRIPTION "support for Atlas" REQUIRED_PACKAGES "atlas VERSION 0.35.1") +ecbuild_add_option( FEATURE VORONOI + DEFAULT OFF + DESCRIPTION "eckit::maths library convex hull/Delaunay triangulation" ) + # REQUIRED_PACKAGES "Qhull COMPONENTS C CXX" + +if( mir_HAVE_VORONOI ) + find_package( Qhull REQUIRED CONFIG ) + + if( NOT TARGET Qhull::qhullcpp OR NOT TARGET Qhull::qhullstatic_r ) + message( FATAL_ERROR "eckit::maths ENABLE_VORONOI requires Qhull C/C++ libraries" ) + endif() + + add_library(Qhull::Qhull INTERFACE IMPORTED) + target_link_libraries(Qhull::Qhull INTERFACE Qhull::qhullcpp Qhull::qhullstatic_r ) +endif() + set(MIR_LIBRARIES mir) diff --git a/src/tools/CMakeLists.txt b/src/tools/CMakeLists.txt index 3a1419f0f..95761bc2a 100644 --- a/src/tools/CMakeLists.txt +++ b/src/tools/CMakeLists.txt @@ -32,6 +32,14 @@ if(eckit_HAVE_ECKIT_CODEC) endforeach() endif() +if(mir_HAVE_VORONOI) + ecbuild_add_executable( + TARGET mir-voronoi-average + SOURCES mir-voronoi-average.cc + LIBS mir Qhull::Qhull + ${mir_INSTALL_TOOLS}) +endif() + if(mir_HAVE_ATLAS) foreach(tool IN ITEMS mir-load-legendre diff --git a/src/tools/mir-voronoi-average.cc b/src/tools/mir-voronoi-average.cc new file mode 100644 index 000000000..078beabec --- /dev/null +++ b/src/tools/mir-voronoi-average.cc @@ -0,0 +1,262 @@ +/* + * (C) Copyright 1996- ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation nor + * does it submit to any jurisdiction. + */ + + +#include +#include +#include +#include + +#include "eckit/option/CmdArgs.h" +#include "libqhullcpp/Qhull.h" +#include "libqhullcpp/QhullFacetList.h" +#include "libqhullcpp/QhullVertexSet.h" + +#include "mir/data/MIRField.h" +#include "mir/input/GribFileInput.h" +#include "mir/repres/Iterator.h" +#include "mir/repres/Representation.h" +#include "mir/tools/MIRTool.h" +#include "mir/util/Exceptions.h" +#include "mir/util/Log.h" +#include "mir/util/Trace.h" + + +namespace mir::tools { + + +class Voronoi { +public: + using coord_t = std::vector; + + Voronoi(const coord_t& coord) { + constexpr size_t N = 3; + + trace::Timer trace("Voronoi"); + + ASSERT(0 < N && coord.size() % N == 0); + + auto pointDimension = static_cast(N); + auto pointCount = static_cast(coord.size() / N); + + qh_ = new orgQhull::Qhull; + ASSERT(qh_ != nullptr); + + std::ostringstream err; + qh_->setErrorStream(&err); + qh_->setOutputStream(&Log::info()); + qh_->enableOutputStream(); + + try { + qh_->runQhull("", pointDimension, pointCount, coord.data(), "v Qz"); + ASSERT(qh_->qhullStatus() == 0); + } + catch (...) { + ASSERT(false); + } + } + + Voronoi(const Voronoi&) = delete; + Voronoi(Voronoi&&) = delete; + + ~Voronoi() { delete qh_; } + + Voronoi& operator=(const Voronoi&) = delete; + Voronoi& operator=(Voronoi&&) = delete; + + coord_t list_voronoi_vertices() const { + coord_t vertices; + vertices.reserve(3 * qh_->facetCount()); + + for (auto f : qh_->facetList()) { + /*if (!f.isUpperDelaunay())*/ { + if (auto centre = f.voronoiVertex(); centre.size() == 3) { + vertices.insert(vertices.end(), centre.begin(), centre.end()); + } + } + } + + return vertices; + } + + std::vector> list_voronoi_facets() const { + std::vector> facets; + + auto facet = [](const auto& f) { + decltype(facets)::value_type facet; + for (const auto& v : f.vertices()) { + facet.emplace_back(v.id()); + } + return facet; + }; + + for (const auto& f : qh_->facetList()) { + if (!f.isUpperDelaunay()) { + facets.emplace_back(facet(f)); + } + } + + return facets; + } + +private: + orgQhull::Qhull* qh_; +}; + + +#if 0 +struct Voronoi { + using coord_t = eckit::maths::Qhull::coord_t; + + explicit Voronoi(const coord_t& coord) { + trace::Timer timer("Voronoi"); + qh_ = std::make_unique(3, coord, "v Qz"); + ASSERT(qh_); + } + + std::vector list_vertices() const { return qh_->list_vertices(); } + + std::vector> list_facets() const { + + std::vector> facets; + facets.reserve(qh_->facetCount()); + + for (const auto& facet : qh_->facetList()) { + const auto vertices = facet.vertices(); + + std::vector f; + f.reserve(vertices.size()); + + for (const auto& vertex : vertices) { + f.emplace_back(vertex.point().id()); + } + + facets.emplace_back(std::move(f)); + } + + return facets; + + + return qh_->list_facets(); } + +#if 0 + void writeGmshFile(const std::vector& points, const std::vector& voronoiVertices, const std::vector>& voronoiRegions, std::ostream&out) { + out << "$MeshFormat" + "\n2.2 0 8" + "\n$EndMeshFormat" + "\n"; + + // Write points + out << "$Nodes\n"; + out << points.size() / 3 + voronoiVertices.size() << "\n"; + int nodeId = 1; + for (size_t i = 0; i < points.size(); i += 3) { + out << nodeId++ << " " << points[i] << " " << points[i + 1] << " " << points[i + 2] << "\n"; + } + for (const auto& v : voronoiVertices) { + out << nodeId++ << " " << v[0] << " " << v[1] << " " << v[2] << "\n"; + } + out << "$EndNodes\n"; + + // Write elements + out << "$Elements\n"; + int numElements = 0; + for (const auto& region : voronoiRegions) { + numElements += region.size(); + } + out << numElements << "\n"; + + int elementId = 1; + int offset = points.size() / 3; + for (size_t i = 0; i < voronoiRegions.size(); ++i) { + for (int index : voronoiRegions[i]) { + out << elementId++ << " 2 2 " << i + 1 << " " << i + 1 << " " << i + 1 << " " << index + offset + 1 << "\n"; + } + } + out << "$EndElements\n"; + } +#endif + +private: + std::unique_ptr qh_; +}; +#endif + + +struct MIRVoronoiAverage : MIRTool { + using MIRTool::MIRTool; + + int numberOfPositionalArguments() const override { return 1; } + + void usage(const std::string& tool) const override { + Log::info() << "\nVoronoi regions." + "\n" + "\nUsage:" + "\n\t" + << tool << " file.grib" << std::endl; + } + + void execute(const eckit::option::CmdArgs& args) override { + ASSERT(args.count() == 1); + + + // Input + std::unique_ptr input(new input::GribFileInput(args(0))); + ASSERT(input->next()); + + + // Get field representation + auto field = input->field(); + repres::RepresentationHandle rep{field.representation()}; + + + // Get field points + Voronoi::coord_t coord(3 * rep->numberOfPoints()); + + size_t i = 0; + for (std::unique_ptr it{rep->iterator()}; it->next();) { + const auto p = it->point3D(); + coord.at(i++) = p[0]; + coord.at(i++) = p[1]; + coord.at(i++) = p[2]; + } + + + // Calculate Voronoi regions + Voronoi v(coord); + + for (const auto& v : v.list_voronoi_vertices()) { + Log::info() << "Vertex: " << v << std::endl; + } + + for (const auto& f : v.list_voronoi_facets()) { + Log::info() << "Facet: "; + for (const auto& v : f) { + Log::info() << v << " "; + } + Log::info() << std::endl; + } + + + // Write Gmsh file + std::ofstream out("filename.gmsh"); + ASSERT(!out); + } +}; + + +} // namespace mir::tools + + +int main(int argc, char** argv) { + mir::tools::MIRVoronoiAverage tool(argc, argv); + return tool.start(); +}