From ec1e8e79b24ea5973ef842d793deecb508d70ef0 Mon Sep 17 00:00:00 2001 From: Alocias Date: Thu, 20 Oct 2022 16:06:11 +0200 Subject: [PATCH 1/4] fix install for clang 14 install using local libs, add tests enable env var --- .gitignore | 11 +++++------ CMakeLists.txt | 11 +++++++++-- README.md | 24 +++++++++++------------ cpp_test/test_sun_position.cpp | 4 ++-- setup.py | 36 +++++++++++++++++++++++----------- src/rasputin/bindings.cpp | 8 ++++---- src/rasputin/solar_position.h | 4 +++- 7 files changed, 60 insertions(+), 38 deletions(-) diff --git a/.gitignore b/.gitignore index 7251aa5..bd73c85 100644 --- a/.gitignore +++ b/.gitignore @@ -2,15 +2,14 @@ cmake-* .idea .env *.swp -lib/CGAL* -lib/pybind11 -lib/armadillo* +lib/**/* src/rasputin/*so src/rasputin.egg-info __pycache__ build dist .DS_Store -lib/boo* -lib/geometry -lib/build_* +.envrc +compile_commands.json +.cache/ +.python-version diff --git a/CMakeLists.txt b/CMakeLists.txt index 57508a9..56455d9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,9 +4,16 @@ project (Rasputin) set(CMAKE_CXX_STANDARD 20) set(CMAKE_CXX_STANDARD_REQUIRED ON) -# set(CMAKE_VERBOSE_MAKEFILE ON) +set(CMAKE_CXX_FLAGS_DEBUG_INIT "-Wall") set(ENABLE_TESTS ON) +if (DEFINED ENV{ENABLE_TESTS}) + if ($ENV{ENABLE_TESTS} STREQUAL "ON") + set(ENABLE_TESTS ON) + else() + set(ENABLE_TESTS OFF) + endif() +endif() # Default location for dependencies set(RASPUTIN_LIB_DIR "${CMAKE_CURRENT_SOURCE_DIR}/lib") @@ -100,4 +107,4 @@ pybind11_add_module(triangulate_dem MODULE SYSTEM # Need to link BLAS for header-only Aramdillo target_link_libraries(triangulate_dem PRIVATE ${RASPUTIN_DEPENDENCIES} - ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES}) + ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES}) diff --git a/README.md b/README.md index b958460..79dcad7 100644 --- a/README.md +++ b/README.md @@ -31,16 +31,16 @@ Installing Rasputin is easy, as the C++ dependencies are header only. Simply, ei For examaple, in some location where you have your source code, type: ``` -git clone https://github.com/pybind/pybind11.git -git clone https://gitlab.com/conradsnicta/armadillo-code.git -git clone https://github.com/boostorg/geometry.git -git clone https://github.com/catchorg/Catch2.git -git clone https://github.com/CGAL/cgal.git +git clone https://github.com/pybind/pybind11.git +git clone https://gitlab.com/conradsnicta/armadillo-code.git +git clone https://github.com/boostorg/geometry.git +git clone https://github.com/catchorg/Catch2.git +git clone https://github.com/CGAL/cgal.git ``` For Howard Hinnant's `date` library to work, enter the rasputin source root directory and checkout the source under the `lib` folder: ``` cd lib -git clone https://github.com/HowardHinnant/date.git +git clone https://github.com/HowardHinnant/date.git ``` Rasputin does not aim at being backwards compatible with older compilers. @@ -108,11 +108,11 @@ from rasputin.reader import Rasterdata from rasputin.mesh import Mesh def construct_rasterdata(): - raster = np.array([0, 0, 0, - 0, 1, 0, + raster = np.array([0, 0, 0, + 0, 1, 0, 0, 0, 0], dtype=np.float32).reshape(3,3) cs = pyproj.CRS.from_epsg(32633) - return Rasterdata(shape=(raster.shape[1], raster.shape[0]), x_min=0, + return Rasterdata(shape=(raster.shape[1], raster.shape[0]), x_min=0, y_max=20, delta_x=10, delta_y=10, array=raster, coordinate_system=cs.to_proj4(), info={}) @@ -179,9 +179,9 @@ Choose "Nedlasting" from the left hand side of the map, and choose "Landsdekkend and finally click DTM10. Download and unpack in, for instance, `$HOME/rasputin_data/dem_archive`, and `export RASPUTIN_DATA_DIR=$HOME/rasputin_data`. -It is possible to include land cover types in your triangulation, through the -[GlobCover dataset](http://due.esrin.esa.int/page_globcover.php) from ESA. It is a raster based -300m (approx) resolution data set that contains 23 different land cover types. +It is possible to include land cover types in your triangulation, through the +[GlobCover dataset](http://due.esrin.esa.int/page_globcover.php) from ESA. It is a raster based +300m (approx) resolution data set that contains 23 different land cover types. Download the data set and unpack it in `$RASPUTIN_DATA_DIR/globcov` to access the land types using the `rasputin.globcov_repository.GlobCovRepository` class. diff --git a/cpp_test/test_sun_position.cpp b/cpp_test/test_sun_position.cpp index c7966fb..82002c7 100644 --- a/cpp_test/test_sun_position.cpp +++ b/cpp_test/test_sun_position.cpp @@ -38,9 +38,9 @@ TEST_CASE("Reference example test", "[reference]") { const auto [Phi, e0] = calendar_solar_position(year, month, day, lat, lon, masl, collectors::azimuth_and_elevation(), fixed_cal_delta_t_calc()); - REQUIRE(abs(Phi - 194.34024) < 1.0e-4); + REQUIRE(std::abs(Phi - 194.34024) < 1.0e-4); const auto [e, Theta] = corrected_solar_elevation(e0, P, T); - REQUIRE(abs(Theta - 50.11162) < 1.0e-4); + REQUIRE(std::abs(Theta - 50.11162) < 1.0e-4); } TEST_CASE("JD test 1", "[jd1]") { diff --git a/setup.py b/setup.py index 567b1d6..188e2d5 100644 --- a/setup.py +++ b/setup.py @@ -1,9 +1,9 @@ import os import re import sys -import sysconfig import platform import subprocess +from pathlib import Path from distutils.version import LooseVersion from setuptools import setup, find_packages, Extension @@ -39,8 +39,11 @@ def run(self): def build_extension(self, ext): extdir = os.path.abspath( os.path.dirname(self.get_ext_fullpath(ext.name))) - cmake_args = ['-DCMAKE_LIBRARY_OUTPUT_DIRECTORY=' + extdir, - '-DPYTHON_EXECUTABLE=' + sys.executable] + cmake_args = [ + '-DCMAKE_LIBRARY_OUTPUT_DIRECTORY=' + extdir, + '-DPYTHON_EXECUTABLE=' + sys.executable, + '-DCMAKE_EXPORT_COMPILE_COMMANDS=ON' + ] cfg = 'Debug' if self.debug else 'Release' build_args = ['--config', cfg] @@ -53,21 +56,30 @@ def build_extension(self, ext): cmake_args += ['-A', 'x64'] build_args += ['--', '/m'] else: - cmake_args += ['-DCMAKE_BUILD_TYPE=' + cfg] + cmake_args += ['-DCMAKE_BUILD_TYPE=' + cfg, f"-DCMAKE_PREFIX_PATH={Path.home() / '.local'};{Path().cwd() / 'lib'}"] build_args += ['--', '-j2'] env = os.environ.copy() env['CXXFLAGS'] = '{} -DVERSION_INFO=\\"{}\\"'.format( env.get('CXXFLAGS', ''), self.distribution.get_version()) + if not os.path.exists(self.build_temp): os.makedirs(self.build_temp) - subprocess.check_call(['cmake', ext.sourcedir] + cmake_args, - cwd=self.build_temp, env=env) - subprocess.check_call(['cmake', '--build', '.'] + build_args, - cwd=self.build_temp) + + cmake_build_command = ['cmake', "-S", ".", "-B", "build"] + cmake_args + _compile_and_check(cmake_build_command, env=env) + _compile_and_check(["cmake", "--build", "build"]) + print() # Add an empty line for cleaner output +def _compile_and_check(cmd, *args, **kwargs): + proc = subprocess.run(cmd, *args, **kwargs, check=True) + try: + proc.check_returncode() + except subprocess.CalledProcessError as err: + raise RuntimeError(f"{cmd} exited with {proc.returncode}\n{proc.stderr}") from err + setup( name='rasputin', version='0.1', @@ -77,7 +89,8 @@ def build_extension(self, ext): long_description='', # tell setuptools to look for any packages under 'src' packages=find_packages('src'), - install_requires=[ + requires=[ + 'Cython', 'numpy', 'pyproj', 'Pillow', @@ -85,15 +98,16 @@ def build_extension(self, ext): 'lxml', 'shapely', 'descartes', - 'meshio' + 'meshio', ], + extras_require={'all': ['pytest']}, # tell setuptools that all packages will be under the 'src' directory # and nowhere else package_dir={'':'src'}, data_files=[("rasputin/web", ["web/index.js", "web/index.html", "web/data.js", "web/favicon.ico"]), ("rasputin/web/js", ["web/js/three.min.js", "web/js/dat.min.js"]), ("rasputin/web/js/controls", ["web/js/controls/OrbitControls.js", "web/js/controls/PointerLockControls.js"])], - # add an extension module named 'python_cpp_example' to the package + # add an extension module named 'python_cpp_example' to the package # 'python_cpp_example' ext_modules=[CMakeExtension('rasputin/rasputin')], # add custom build_ext command diff --git a/src/rasputin/bindings.cpp b/src/rasputin/bindings.cpp index d1c2c77..6af9891 100644 --- a/src/rasputin/bindings.cpp +++ b/src/rasputin/bindings.cpp @@ -383,13 +383,13 @@ PYBIND11_MODULE(triangulate_dem, m) { using namespace date; #endif const auto tp = sys_days{January/1/1970} + seconds(long(std::round(timestamp))); - return rasputin::solar_position::time_point_solar_position(tp, geographic_latitude, geographic_longitude, masl, - rasputin::solar_position::collectors::azimuth_and_elevation(), + return rasputin::solar_position::time_point_solar_position(tp, geographic_latitude, geographic_longitude, masl, + rasputin::solar_position::collectors::azimuth_and_elevation(), rasputin::solar_position::delta_t_calculator::coarse_timestamp_calc()); }, "Compute azimuth and elevation of sun for given UTC timestamp.") .def("calendar_solar_position", [] (unsigned int year,unsigned int month, double day, const double geographic_latitude, const double geographic_longitude, const double masl) { - return rasputin::solar_position::calendar_solar_position(year, month, day, geographic_latitude, geographic_longitude, masl, - rasputin::solar_position::collectors::azimuth_and_elevation(), + return rasputin::solar_position::calendar_solar_position(year, month, day, geographic_latitude, geographic_longitude, masl, + rasputin::solar_position::collectors::azimuth_and_elevation(), rasputin::solar_position::delta_t_calculator::coarse_date_calc()); }, "Compute azimuth and elevation of sun for given UT calendar coordinate.") .def("solar_elevation_correction", &rasputin::solar_position::corrected_solar_elevation, "Correct elevation based on pressure and temperature.") diff --git a/src/rasputin/solar_position.h b/src/rasputin/solar_position.h index 109952d..e6d8caf 100644 --- a/src/rasputin/solar_position.h +++ b/src/rasputin/solar_position.h @@ -10,6 +10,8 @@ #include #endif #include +#include +#include namespace rasputin::solar_position::collectors { @@ -47,7 +49,7 @@ auto coarse_timestamp_calc() { using namespace date; #endif - auto ep = sys_days{January/1/1970}; + auto ep = sys_days(January/1/1970); return [ep] (system_clock::time_point tp) { const auto year = round(duration_cast(tp - ep).count()*1.0/duration_cast(years(1)).count()) + 1970; From d2bd3c64d958fbd6089d50eb212383c97bf4ac76 Mon Sep 17 00:00:00 2001 From: Alocias Date: Thu, 1 Dec 2022 14:19:37 +0000 Subject: [PATCH 2/4] add setup.cfg --- .gitignore | 5 +++-- Pipfile | 18 ------------------ setup.cfg | 52 ++++++++++++++++++++++++++++++++++++++++++++++++++++ setup.py | 52 +++++++++++----------------------------------------- 4 files changed, 66 insertions(+), 61 deletions(-) delete mode 100644 Pipfile create mode 100644 setup.cfg diff --git a/.gitignore b/.gitignore index bd73c85..00c17ba 100644 --- a/.gitignore +++ b/.gitignore @@ -3,8 +3,8 @@ cmake-* .env *.swp lib/**/* -src/rasputin/*so -src/rasputin.egg-info +**/*.egg-info +**/*.so __pycache__ build dist @@ -13,3 +13,4 @@ dist compile_commands.json .cache/ .python-version +*.eggs diff --git a/Pipfile b/Pipfile deleted file mode 100644 index f800235..0000000 --- a/Pipfile +++ /dev/null @@ -1,18 +0,0 @@ -[[source]] -url = 'https://pypi.python.org/simple' -verify_ssl = true -name = 'pypi' - -[packages] -numpy = '*' -pyproj = '*' -Pillow = '*' -h5py = '*' -lxml = '*' -shapely = '*' -descartes = '*' -meshio = '*' -rasputin = {path = "."} - -[dev-packages] -pytest = '*' diff --git a/setup.cfg b/setup.cfg new file mode 100644 index 0000000..09867c5 --- /dev/null +++ b/setup.cfg @@ -0,0 +1,52 @@ +[metadata] +name = rasputin +version = 0.1 +author = Ola Skavhaug +author_email = ola@xal.no +description = A simple tool for constructing TINs from DEMs +readme = README.md + +[options] +zip_safe = False +include_package_data = True +install_requires = + cython + numpy + pyproj + Pillow + h5py + lxml + shapely + descartes + meshio +setup_requires = + cython + h5py +python_requires = >=3.10 + +[options.packages.find] +where = src +include = rasputin/* + +[options.extras_require] +all = pytest +ipy = ipython + +[options.package_data] +web = + index.js + index.html + data.js + favicon.ico +web.js = + three.min.js + dat.min.js +web.js.controls = + OrbitControls.js + PointerLockControls.js + +[options.entry_points] +console_scripts = + rasputin_web = rasputin.web_visualize:visualize_tin + rasputin_store = rasputin.application:store_tin + rasputin_shade = rasputin.application:compute_shades diff --git a/setup.py b/setup.py index 188e2d5..a540a5f 100644 --- a/setup.py +++ b/setup.py @@ -6,11 +6,9 @@ from pathlib import Path from distutils.version import LooseVersion -from setuptools import setup, find_packages, Extension +from setuptools import setup, Extension from setuptools.command.build_ext import build_ext -web_templates = "src/rasputin/web/*" - class CMakeExtension(Extension): def __init__(self, name, sourcedir=''): @@ -28,10 +26,13 @@ def run(self): ", ".join(e.name for e in self.extensions)) if platform.system() == "Windows": - cmake_version = LooseVersion(re.search(r'version\s*([\d.]+)', - out.decode()).group(1)) - if cmake_version < '3.1.0': - raise RuntimeError("CMake >= 3.1.0 is required on Windows") + cmake_version = re.search(r'version\s*([\d.]+)', out.decode()) + if cmake_version is not None: + cmake_version = LooseVersion(cmake_version.group(1)) + if cmake_version < '3.1.0': + raise RuntimeError("CMake >= 3.1.0 is required on Windows") + else: + raise ValueError("cmake version not found") for ext in self.extensions: self.build_extension(ext) @@ -62,7 +63,8 @@ def build_extension(self, ext): env = os.environ.copy() env['CXXFLAGS'] = '{} -DVERSION_INFO=\\"{}\\"'.format( env.get('CXXFLAGS', ''), - self.distribution.get_version()) + self.distribution.get_version() + ) if not os.path.exists(self.build_temp): os.makedirs(self.build_temp) @@ -81,41 +83,9 @@ def _compile_and_check(cmd, *args, **kwargs): raise RuntimeError(f"{cmd} exited with {proc.returncode}\n{proc.stderr}") from err setup( - name='rasputin', - version='0.1', - author='Ola Skavhaug', - author_email='ola@xal.no', - description='A simple tool for constructing TINs from DEMs', - long_description='', - # tell setuptools to look for any packages under 'src' - packages=find_packages('src'), - requires=[ - 'Cython', - 'numpy', - 'pyproj', - 'Pillow', - 'h5py', - 'lxml', - 'shapely', - 'descartes', - 'meshio', - ], - extras_require={'all': ['pytest']}, - # tell setuptools that all packages will be under the 'src' directory - # and nowhere else - package_dir={'':'src'}, - data_files=[("rasputin/web", ["web/index.js", "web/index.html", "web/data.js", "web/favicon.ico"]), - ("rasputin/web/js", ["web/js/three.min.js", "web/js/dat.min.js"]), - ("rasputin/web/js/controls", ["web/js/controls/OrbitControls.js", "web/js/controls/PointerLockControls.js"])], # add an extension module named 'python_cpp_example' to the package # 'python_cpp_example' - ext_modules=[CMakeExtension('rasputin/rasputin')], + ext_modules=[CMakeExtension('src/rasputin')], # add custom build_ext command cmdclass=dict(build_ext=CMakeBuild), - zip_safe=False, - entry_points={ - 'console_scripts': ['rasputin_web = rasputin.web_visualize:visualize_tin', - 'rasputin_store = rasputin.application:store_tin', - 'rasputin_shade = rasputin.application:compute_shades']} ) - From 11b951e88eb6b767bd7b0f637866633012a60c55 Mon Sep 17 00:00:00 2001 From: Alocias Date: Thu, 20 Oct 2022 16:06:11 +0200 Subject: [PATCH 3/4] fix install for clang 14 install using local libs, add tests enable env lar, remove CGAL, add namespace Triangulation, split debug and release builds to seperate folders, add tests for point within inscribed circle, add tests for ccw sorting. --- CMakeLists.txt | 13 +- cpp_test/CMakeLists.txt | 2 + cpp_test/test_triangulate.cpp | 205 ++++++ setup.py | 12 +- src/rasputin/bindings.cpp | 467 +++---------- src/rasputin/solar_position.h | 2 +- src/rasputin/triangulate_dem.h | 1127 ++++++++------------------------ 7 files changed, 589 insertions(+), 1239 deletions(-) create mode 100644 cpp_test/test_triangulate.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 56455d9..1155419 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -5,6 +5,7 @@ set(CMAKE_CXX_STANDARD 20) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_FLAGS_DEBUG_INIT "-Wall") +set(CMAKE_CXX_FLAGS_RELEASE_INIT "-Wall") set(ENABLE_TESTS ON) if (DEFINED ENV{ENABLE_TESTS}) @@ -61,11 +62,11 @@ list(APPEND RASPUTIN_DEPENDENCIES date) # CGAL # ---- # Required for now but should be optional in the future -find_package(CGAL REQUIRED) -if (CGAL_FOUND) - add_definitions(-DHAS_CGAL) # For future use when CGAL is optional - list(APPEND RASPUTIN_DEPENDENCIES CGAL::CGAL) -endif() +# find_package(CGAL REQUIRED) +# if (CGAL_FOUND) +# add_definitions(-DHAS_CGAL) # For future use when CGAL is optional +# list(APPEND RASPUTIN_DEPENDENCIES CGAL::CGAL) +# endif() # NOTE: Calling`include(${CGAL_USE_FILE})` will modify INCLUDE_DIRECTORIES, # with the consquence that directories added by CGAL will be searched before # anything else, potentiall messing up the include order. @@ -84,8 +85,8 @@ list(APPEND RASPUTIN_DEPENDENCIES Armadillo) # BLAS # ---- # Need to link to BLAS libraries when not using armadillo wrappers -find_package(BLAS REQUIRED) find_package(LAPACK) +find_package(BLAS REQUIRED) # Header-only Rasputin library diff --git a/cpp_test/CMakeLists.txt b/cpp_test/CMakeLists.txt index 94cc2bf..e7bcc83 100644 --- a/cpp_test/CMakeLists.txt +++ b/cpp_test/CMakeLists.txt @@ -8,3 +8,5 @@ find_package(Catch2 3 REQUIRED) # ----------- add_executable(test_sun_position test_sun_position.cpp) target_link_libraries(test_sun_position PRIVATE rasputin Catch2::Catch2WithMain ${RASPUTIN_DEPENDENCIES}) +add_executable(test_triangulate test_triangulate.cpp) +target_link_libraries(test_triangulate PRIVATE rasputin Catch2::Catch2WithMain ${RASPUTIN_DEPENDENCIES}) diff --git a/cpp_test/test_triangulate.cpp b/cpp_test/test_triangulate.cpp new file mode 100644 index 0000000..aeaf018 --- /dev/null +++ b/cpp_test/test_triangulate.cpp @@ -0,0 +1,205 @@ +#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one cpp file +#include + +#include +#include + +#include +#include +#include + +#include + +namespace bg = boost::geometry; +using Point = bg::model::point; +namespace boost::geometry::traits { + using shared_type = std::shared_ptr; + BOOST_GEOMETRY_DETAIL_SPECIALIZE_POINT_TRAITS(shared_type, 3, double, cs::cartesian) + template<> + struct access { + static inline double get(shared_type const& p) { return p->get<0>(); } + }; + template<> + struct access { + static inline double get(shared_type const& p) { return p->get<1>(); } + }; + template<> + struct access { + static inline double get(shared_type const& p) { return p->get<2>(); } + }; +} + +namespace rasputin::test_utils { + namespace bg = boost::geometry; + using Point = bg::model::point; + using Triangulation = rasputin::Triangulation; + + void check_points_exists(const std::vector& triangles, const std::vector& points) { + for (const Point& point : points) { + bool found = false; + for (const Triangulation::Face& triangle : triangles) { + const auto& [a, b, c] = triangle; + if (&point == a || &point == b || &point == c) { + found = true; + break; + } + } + + INFO("Point (" << bg::get<0>(point) << " " << bg::get<1>(point) << " " << bg::get<2>(point) << ") not found"); + CHECK(found); + } + } + + void check_points_inside_ccirc( + const std::vector& triangles, const std::vector& points + ) { + for ( + Triangulation::Faces::const_iterator triangle1=triangles.begin(); + triangle1 != triangles.end(); + ++triangle1 + ) { + auto [t1a, t1b, t1c] = *triangle1; + for ( + Triangulation::Faces::const_iterator triangle2=triangle1 + 1; + triangle2 != triangles.end(); + ++triangle2 + ) { + bool ta_inside = Triangulation::_is_inside_ccirc(t1a, *triangle2); + bool tb_inside = Triangulation::_is_inside_ccirc(t1b, *triangle2); + bool tc_inside = Triangulation::_is_inside_ccirc(t1c, *triangle2); + INFO("Point (" << bg::get<0>(t1a) << " " << bg::get<1>(t1a) << " " << bg::get<2>(t1a) << ") inside"); + CHECK(!ta_inside); + INFO("Point (" << bg::get<0>(t1a) << " " << bg::get<1>(t1a) << " " << bg::get<2>(t1a) << ") inside"); + CHECK(!tb_inside); + INFO("Point (" << bg::get<0>(t1c) << " " << bg::get<1>(t1c) << " " << bg::get<2>(t1c) << ") inside"); + CHECK(!tc_inside); + } + } + } + + void check_triangles_unique(const std::vector& triangles) { + for ( + Triangulation::Faces::const_iterator triangle1=triangles.begin(); + triangle1 != triangles.end(); + ++triangle1 + ) { + const auto& [a1, b1, c1] = *triangle1; + for ( + Triangulation::Faces::const_iterator triangle2=triangle1 + 1; + triangle2 != triangles.end(); + ++triangle2 + ) { + const auto& [a2, b2, c2] = *triangle2; + bool same = + (&a1 == &a2 && &b1 == &b2 && &c1 == &c2) + || (&a1 == &c2 && &b1 == &a2 && &c1 == &b2) + || (&a1 == &b2 && &b1 == &c2 && &c1 == &a2); + INFO( + "Duplicate Triangle [(" + << bg::get<0>(a1) << " " << bg::get<1>(a1) << " " << bg::get<2>(a1) << ") " + << bg::get<0>(b1) << " " << bg::get<1>(b1) << " " << bg::get<2>(b1) << ") " + << bg::get<0>(c1) << " " << bg::get<1>(c1) << " " << bg::get<2>(c1) << ")]" + ); + CHECK(!same); + } + } + } + + TEST_CASE("[Test sort polar]", "[Triangulation]") { + std::random_device rd; + std::default_random_engine engine(rd()); + std::uniform_real_distribution distr(0.0, 1.0); + + const Point a = {0.0, 0.0, 0.0}; + const Point b = {distr(engine), distr(engine), 0.0}; + const Point c = {distr(engine), 0.0, 0.0}; + + Triangulation::Face triangle = Triangulation::_sort_points(&a, &b, &c); + REQUIRE(std::get<0>(triangle) == &a); + REQUIRE(std::get<1>(triangle) == &c); + REQUIRE(std::get<2>(triangle) == &b); + + triangle = Triangulation::_sort_points(&a, &c, &b); + REQUIRE(std::get<0>(triangle) == &a); + REQUIRE(std::get<1>(triangle) == &c); + REQUIRE(std::get<2>(triangle) == &b); + + triangle = Triangulation::_sort_points(&c, &a, &b); + REQUIRE(std::get<0>(triangle) == &c); + REQUIRE(std::get<1>(triangle) == &b); + REQUIRE(std::get<2>(triangle) == &a); + + triangle = Triangulation::_sort_points(&c, &b, &a); + REQUIRE(std::get<0>(triangle) == &c); + REQUIRE(std::get<1>(triangle) == &b); + REQUIRE(std::get<2>(triangle) == &a); + } + + TEST_CASE("[Test inside circumcircle]", "[Triangulation]") { + const Point a = {0.0, 0.0, 0.0}; + const Point b = {0.5, 1.0, 0.0}; + const Point c = {1.0, 0.0, 0.0}; + Triangulation::Face triangle = {&a, &b, &c}; + + Point inside_point = {0.25, 0.25, 0.0}; + bool is_inside = Triangulation::_is_inside_ccirc(&inside_point, triangle); + REQUIRE(is_inside == true); + + inside_point = {100.25, 100.25, 0.0}; + is_inside = Triangulation::_is_inside_ccirc(&inside_point, triangle); + REQUIRE(is_inside == false); + + inside_point = {0.85, 0.25, 0.0}; + is_inside = Triangulation::_is_inside_ccirc(&inside_point, triangle); + REQUIRE(is_inside == true); + + inside_point = {0.15, 0.25, 0.0}; + is_inside = Triangulation::_is_inside_ccirc(&inside_point, triangle); + REQUIRE(is_inside == true); + + inside_point = {0.95, 0.25, 0.0}; + is_inside = Triangulation::_is_inside_ccirc(&inside_point, triangle); + REQUIRE(is_inside == true); + + inside_point = {0.05, 0.25, 0.0}; + is_inside = Triangulation::_is_inside_ccirc(&inside_point, triangle); + REQUIRE(is_inside == true); + + inside_point = {1.15, 0.25, 0.0}; + is_inside = Triangulation::_is_inside_ccirc(&inside_point, triangle); + REQUIRE(is_inside == false); + + inside_point = {-0.85, 0.25, 0.0}; + is_inside = Triangulation::_is_inside_ccirc(&inside_point, triangle); + REQUIRE(is_inside == false); + } + + TEST_CASE("[Test unconstrained Delaunay triangulation]", "[Triangulation]") { + std::vector raw_data = {0, 0, 0, 0, 0, 0, 0, 0}; + RasterData data = {0.0, 1.0, 1.0, 1.0, 2, 2, raw_data.data()}; + + auto mesh = Triangulation::delaunay(data, ""); + + check_points_exists(mesh.triangles, data.get_points()); + check_points_inside_ccirc(mesh.triangles, data.get_points()); + check_triangles_unique(mesh.triangles); + + REQUIRE(mesh.num_faces() == 2); + } + + // TEST_CASE("[Test random unconstrained Delaunay triangulation]", "[Triangulation]") { + // std::random_device rd; + // std::default_random_engine engine(rd()); + // std::uniform_real_distribution distr(0.0, 1.0); + // + // std::vector raw_data = {0, 0, 0, 0, 0, 0, 0, 0}; + // RasterData data = {0.0, 1.0, 1.0, 1.0, 2, 2, raw_data.data()}; + // + // auto mesh = Triangulation::delaunay(data, ""); + // + // REQUIRE(mesh.num_faces() == 2); + // + // check_points_exists(mesh.triangles, data.points); + // check_points_inside_ccirc(mesh.triangles, data.points); + // } +} diff --git a/setup.py b/setup.py index a540a5f..310c2bd 100644 --- a/setup.py +++ b/setup.py @@ -46,7 +46,13 @@ def build_extension(self, ext): '-DCMAKE_EXPORT_COMPILE_COMMANDS=ON' ] - cfg = 'Debug' if self.debug else 'Release' + if self.debug: + cfg = 'Debug' + build_path = "build/debug" + else: + cfg = 'Release' + build_path = "build/release" + build_args = ['--config', cfg] if platform.system() == "Windows": @@ -69,9 +75,9 @@ def build_extension(self, ext): if not os.path.exists(self.build_temp): os.makedirs(self.build_temp) - cmake_build_command = ['cmake', "-S", ".", "-B", "build"] + cmake_args + cmake_build_command = ['cmake', "-S", ".", "-B", build_path] + cmake_args _compile_and_check(cmake_build_command, env=env) - _compile_and_check(["cmake", "--build", "build"]) + _compile_and_check(["cmake", "--build", build_path]) print() # Add an empty line for cleaner output diff --git a/src/rasputin/bindings.cpp b/src/rasputin/bindings.cpp index 6af9891..f7bc40e 100644 --- a/src/rasputin/bindings.cpp +++ b/src/rasputin/bindings.cpp @@ -1,35 +1,18 @@ -#include "triangulate_dem.h" -#include "solar_position.h" +#include + +#include + #include -#include #include -#include - -#include +#include -// Surface mesh simplication policies -#include -#include -#include -#include -#include -#include -#include -#include +#include "triangulate_dem.h" namespace py = pybind11; -namespace SMS = CGAL::Surface_mesh_simplification; - -PYBIND11_MAKE_OPAQUE(rasputin::point3_vector); -PYBIND11_MAKE_OPAQUE(rasputin::point2_vector); -PYBIND11_MAKE_OPAQUE(rasputin::face_vector); -PYBIND11_MAKE_OPAQUE(rasputin::double_vector); -PYBIND11_MAKE_OPAQUE(rasputin::index_vector); -PYBIND11_MAKE_OPAQUE(CGAL::MultiPolygon); -PYBIND11_MAKE_OPAQUE(std::vector>); -PYBIND11_MAKE_OPAQUE(std::vector>); +namespace rasputin { template +requires std::floating_point py::buffer_info vecarray_buffer(std::vector> &v) { return py::buffer_info( &v[0], /* Pointer to buffer */ @@ -42,371 +25,83 @@ py::buffer_info vecarray_buffer(std::vector> &v) { } template -std::vector>* vecarray_from_numpy(py::array_t buf) { - auto info = buf.request(); - if (info.ndim < 1 or info.ndim > 2) - throw py::type_error("Can only convert one- and two-dimensional arrays."); - - // Make sure total size can be preserved - auto size = info.shape[0]; - if (info.ndim == 1 and size % n) - throw py::type_error("Size of one-dimensional array must be divisible by!" + std::to_string(n) + "."); - - if (info.ndim == 2 and info.shape[1] != n) - throw py::type_error("Second dimension does not have size equal to " + std::to_string(n) + "."); - - if (info.ndim == 2) - size *= n; - - auto vec = std::unique_ptr>> (new std::vector> ()); - vec->reserve(size / n); - - // Copy values - double* p = static_cast(info.ptr); - double* end = p + size; - for (; p < end; p += n) { - auto a = std::array {}; - for (auto k = 0; k < n; k++) - a[k] = *(p + k); - vec->push_back(a); - } - - return vec.release(); -} - -template +requires std::floating_point py::buffer_info vector_buffer(std::vector &v) { return py::buffer_info(&v[0], sizeof(T), py::format_descriptor::format(), 1, { v.size() }, { sizeof(T) }); } - -template -CGAL::MultiPolygon difference_polygons(const P0& polygon0, const P1& polygon1) { - CGAL::MultiPolygon difference_result; - CGAL::difference(polygon0, polygon1, std::back_inserter(difference_result)); - - return difference_result; -} - -template -CGAL::MultiPolygon intersect_polygons(const P0& polygon0, const P1& polygon1) { - CGAL::MultiPolygon intersection_result; - CGAL::intersection(polygon0, polygon1, std::back_inserter(intersection_result)); - - return intersection_result; -} - -template -CGAL::MultiPolygon join_polygons(const P0& polygon0, const P1& polygon1) { - CGAL::Polygon joined; - CGAL::MultiPolygon join_result; - if (CGAL::join(polygon0, polygon1, joined)) - join_result.push_back(std::move(joined)); - else { - join_result.push_back(CGAL::Polygon(polygon0)); - join_result.push_back(CGAL::Polygon(polygon1)); - } - return join_result; -} - - -CGAL::MultiPolygon join_multipolygons(const CGAL::MultiPolygon& polygon0, const CGAL::MultiPolygon& polygon1) { - CGAL::MultiPolygon join_result; - CGAL::join(polygon0.begin(), polygon0.end(), polygon1.begin(), polygon1.end(), std::back_inserter(join_result)); - - return join_result; -} - - - -CGAL::SimplePolygon polygon_from_numpy(py::array_t& buf) { - auto info = buf.request(); - if (info.ndim < 1 or info.ndim > 2) - throw py::type_error("Can only convert one- and two-dimensional arrays."); - - // Make sure total size can be preserved - auto size = info.shape[0]; - if (info.ndim == 1 and size % 2) - throw py::type_error("Size of one-dimensional array must be divisible by 2."); - - if (info.ndim == 2 and info.shape[1] != 2) - throw py::type_error("Second dimension does not have size equal to 2."); - - if (info.ndim == 2) - size *= 2; - - CGAL::SimplePolygon exterior; - - // Copy values - double* p = static_cast(info.ptr); - double* end = p + size; - for (; p < end; p += 2) - exterior.push_back(CGAL::Point2(*p, *(p+1))); - return exterior; -} - -template -void bind_rasterdata(py::module &m, const std::string& pyname) { - py::class_, std::unique_ptr>>(m, pyname.c_str(), py::buffer_protocol()) - .def(py::init([] (py::array_t& data_array, double x_min, double y_max, double delta_x, double delta_y) { - auto buffer = data_array.request(); - int m = buffer.shape[0], n = buffer.shape[1]; - - return rasputin::RasterData(x_min, y_max, delta_x, delta_y, n, m, static_cast(buffer.ptr) ); - }), py::return_value_policy::take_ownership, py::keep_alive<1, 2>(), - py::arg("data_array").noconvert(), py::arg("x_min"), py::arg("y_max"), py::arg("delta_x"), py::arg("delta_y")) - .def_buffer([] (rasputin::RasterData& self) { - return py::buffer_info( - self.data, - sizeof(FT), - py::format_descriptor::format(), - 2, - std::vector { self.num_points_y, self.num_points_x }, - { sizeof(FT) * self.num_points_x, sizeof(FT) } - ); - }) - .def_readwrite("x_min", &rasputin::RasterData::x_min) - .def_readwrite("y_max", &rasputin::RasterData::y_max) - .def_readwrite("delta_x", &rasputin::RasterData::delta_x) - .def_readwrite("delta_y", &rasputin::RasterData::delta_y) - .def_readonly("num_points_x", &rasputin::RasterData::num_points_x) - .def_readonly("num_points_y", &rasputin::RasterData::num_points_y) - .def_property_readonly("x_max", &rasputin::RasterData::get_x_max) - .def_property_readonly("y_min", &rasputin::RasterData::get_y_min) - .def("__getitem__", [](rasputin::RasterData& self, std::pair idx) { - auto [i, j] = idx; - return self.data[self.num_points_x * i + j]; }) - .def("get_indices", &rasputin::RasterData::get_indices) - .def("exterior", &rasputin::RasterData::exterior, py::return_value_policy::take_ownership) - .def("contains", &rasputin::RasterData::contains) - .def("get_interpolated_value_at_point", &rasputin::RasterData::get_interpolated_value_at_point); -} - -template -void bind_make_mesh(py::module &m) { - m.def("make_mesh", - [] (const R& raster_data, const P polygon, const std::string proj4_str) { - return rasputin::mesh_from_raster(raster_data, polygon, proj4_str); - }, py::return_value_policy::take_ownership) - .def("make_mesh", - [] (const R& raster_data, const std::string proj4_str) { - return rasputin::mesh_from_raster(raster_data, proj4_str); - }, py::return_value_policy::take_ownership); -} - -template -void bind_raster_list(py::module &m, const std::string& pyname) { - py::class_>, std::unique_ptr>>> (m, pyname.c_str()) - .def(py::init( [] () {std::vector> self; return self;})) - .def("add_raster", - [] (std::vector>& self, rasputin::RasterData raster_data) { - self.push_back(raster_data); - }, py::keep_alive<1,2>()) - .def("__getitem__", - [] (std::vector>& self, int index) { - return self.at(index); - }, py::return_value_policy::reference_internal); +#define MAKE_OPAQUE(T, Point) (PYBIND11_MAKE_OPAQUE(rasputin::Points)) +#define BIND_MODULE(T, Point) {\ + PYBIND10_MODULE(triangulate_dem, m) {\ + py::bind_vector>(m, "points", py::buffer_protocol())\ + .def_buffer(&vecarray_buffer::value>);\ + py::bind_vector>(m, "int_vector");\ + py::bind_vector>>(m, "shadow_vector");\ + };\ } - -PYBIND11_MODULE(triangulate_dem, m) { - py::bind_vector(m, "point3_vector", py::buffer_protocol()) - .def_buffer(&vecarray_buffer) - .def("from_numpy", &vecarray_from_numpy); - py::bind_vector(m, "point2_vector", py::buffer_protocol()) - .def_buffer(&vecarray_buffer) - .def("from_numpy", &vecarray_from_numpy); - py::bind_vector(m, "face_vector", py::buffer_protocol()) - .def_buffer(&vecarray_buffer) - .def("from_numpy", &vecarray_from_numpy); - py::bind_vector(m, "index_vector", py::buffer_protocol()) - .def_buffer(&vecarray_buffer) - .def("from_numpy", &vecarray_from_numpy); - py::bind_vector(m, "double_vector", py::buffer_protocol()) - .def_buffer(&vector_buffer); - - py::bind_vector>(m, "int_vector"); - py::bind_vector>>(m, "shadow_vector"); - - - bind_rasterdata(m, "raster_data_float"); - bind_rasterdata(m, "raster_data_double"); - - bind_raster_list(m, "raster_list_float"); - bind_raster_list(m, "raster_list_double"); - - bind_make_mesh>, CGAL::SimplePolygon>(m); - bind_make_mesh>, CGAL::SimplePolygon>(m); - bind_make_mesh>, CGAL::Polygon>(m); - bind_make_mesh>, CGAL::Polygon>(m); - // bind_make_mesh>, CGAL::MultiPolygon>(m); - // bind_make_mesh>, CGAL::MultiPolygon>(m); - bind_make_mesh, CGAL::SimplePolygon>(m); - bind_make_mesh, CGAL::SimplePolygon>(m); - bind_make_mesh, CGAL::Polygon>(m); - bind_make_mesh, CGAL::Polygon>(m); - // bind_make_mesh, CGAL::MultiPolygon>(m); - // bind_make_mesh, CGAL::MultiPolygon>(m); - - py::class_>(m, "simple_polygon") - .def(py::init(&polygon_from_numpy)) - .def("num_vertices", &CGAL::SimplePolygon::size) - .def("array", [] (const CGAL::SimplePolygon& self) { - py::array_t result(self.size() * 2); - result.resize(py::array::ShapeContainer({static_cast(self.size()), 2}), true); - auto info = result.request(); - double* data = static_cast(info.ptr); - for (auto v = self.vertices_begin(); v != self.vertices_end(); ++v) { - data[0] = v->x(); - data[1] = v->y(); - data += 2; - } - return result; - }, py::return_value_policy::move) - .def("join", &join_polygons) - .def("join", &join_polygons) - .def("difference", &difference_polygons) - .def("difference", &difference_polygons) - .def("intersection", &intersect_polygons) - .def("intersection", &intersect_polygons); - - py::class_>(m, "polygon") - .def(py::init([] (py::array_t& buf) { - const CGAL::SimplePolygon exterior = polygon_from_numpy(buf); - return CGAL::Polygon(exterior);})) - .def("holes", [] (const CGAL::Polygon& self) { - py::list result; - for (auto h = self.holes_begin(); h != self.holes_end(); ++h) - result.append(*h); - return result; - }) - .def("extract_boundaries", [] (const CGAL::Polygon& self) {return CGAL::extract_boundaries(self);}, - py::return_value_policy::take_ownership) - .def("exterior", [] (const CGAL::Polygon& self) {return self.outer_boundary();}) - .def("join", &join_polygons) - .def("join", &join_polygons) - .def("difference", &difference_polygons) - .def("difference", &difference_polygons) - .def("intersection", &intersect_polygons) - .def("intersection", &intersect_polygons); - - py::class_>(m, "multi_polygon") - .def(py::init( - [] (const CGAL::Polygon& polygon) { - CGAL::MultiPolygon self; - self.push_back(CGAL::Polygon(polygon)); - return self; - })) - .def(py::init( - [] (const CGAL::SimplePolygon& polygon) { - CGAL::MultiPolygon self; - self.push_back(CGAL::Polygon(polygon)); - return self; - })) - .def("num_parts", &CGAL::MultiPolygon::size) - .def("parts", - [] (const CGAL::MultiPolygon& self) { - py::list result; - for (auto p = self.begin(); p != self.end(); ++p) - result.append(*p); - return result; - }) - .def("extract_boundaries", [] (const CGAL::MultiPolygon& self) {return CGAL::extract_boundaries(self);}, - py::return_value_policy::take_ownership) - .def("join", - [] (const CGAL::MultiPolygon a, CGAL::SimplePolygon b) { - return join_multipolygons(a, CGAL::MultiPolygon({static_cast(b)})); - }) - .def("join", - [] (const CGAL::MultiPolygon a, CGAL::Polygon b) { - return join_multipolygons(a, CGAL::MultiPolygon({b})); - }) - .def("join", - [] (const CGAL::MultiPolygon a, CGAL::MultiPolygon b) { - return join_multipolygons(a, b); - }) - .def("__getitem__", [] (const CGAL::MultiPolygon& self, int idx) { - return self.at(idx); - }, py::return_value_policy::reference_internal); - - py::class_>(m, "Mesh") - .def("lindstrom_turk_by_ratio", - [] (const rasputin::Mesh& self, double ratio) { - return self.coarsen(SMS::Count_ratio_stop_predicate(ratio), - SMS::LindstromTurk_placement(), - SMS::LindstromTurk_cost()); - }, py::return_value_policy::take_ownership, - "Simplify the mesh.\n\nThe LindstromTurk cost and placement strategy is used, and simplification process stops when the number of undirected edges drops below the size threshold.") - .def("lindstrom_turk_by_size", - [] (const rasputin::Mesh& self, int max_size) { - return self.coarsen(SMS::Count_stop_predicate(max_size), - SMS::LindstromTurk_placement(), - SMS::LindstromTurk_cost()); - }, py::return_value_policy::take_ownership, - "Simplify the mesh.\n\nThe LindstromTurk cost and placement strategy is used, and simplification process stops when the number of undirected edges drops below the ratio threshold.") - .def("copy", &rasputin::Mesh::copy, py::return_value_policy::take_ownership) - .def("extract_sub_mesh", &rasputin::Mesh::extract_sub_mesh, py::return_value_policy::take_ownership) - - .def_property_readonly("num_vertices", &rasputin::Mesh::num_vertices) - .def_property_readonly("num_edges", &rasputin::Mesh::num_edges) - .def_property_readonly("num_faces", &rasputin::Mesh::num_faces) - - .def_property_readonly("points", &rasputin::Mesh::get_points, py::return_value_policy::reference_internal) - .def_property_readonly("faces", &rasputin::Mesh::get_faces, py::return_value_policy::reference_internal); - - m.def("compute_shadow", (std::vector (*)(const rasputin::Mesh &, const rasputin::point3 &))&rasputin::compute_shadow, "Compute shadows for given topocentric sun position.") - .def("compute_shadow", (std::vector (*)(const rasputin::Mesh &, const double, const double))&rasputin::compute_shadow, "Compute shadows for given azimuth and elevation.") - .def("compute_shadows", &rasputin::compute_shadows, "Compute shadows for a series of times and ray directions.") - .def("construct_mesh", - [] (const rasputin::point3_vector& points, const rasputin::face_vector & faces, const std::string proj4_str) { - rasputin::VertexIndexMap index_map; - rasputin::FaceDescrMap face_map; - return rasputin::Mesh(rasputin::construct_mesh(points, faces, index_map, face_map), proj4_str); - }, py::return_value_policy::take_ownership) - .def("surface_normals", &rasputin::surface_normals, - "Compute surface normals for all faces in the mesh.", - py::return_value_policy::take_ownership) - .def("point_normals", &rasputin::point_normals, "Compute surface normals for all vertices in the mesh.") - .def("orient_tin", &rasputin::orient_tin, "Orient all triangles in the TIN and returns their surface normals.") - .def("extract_lakes", &rasputin::extract_lakes, "Extract lakes as separate face list.") - .def("compute_slopes", &rasputin::compute_slopes,"Compute slopes (i.e. angles relative to xy plane) for the all the vectors in list.") - .def("compute_aspects", &rasputin::compute_aspects, "Compute aspects for the all the vectors in list.") - .def("extract_avalanche_expositions", &rasputin::extract_avalanche_expositions, "Extract avalanche exposed cells.") - .def("consolidate", &rasputin::consolidate, "Make a stand alone consolidated tin.") - .def("cell_centers", &rasputin::cell_centers, "Compute cell centers for triangulation.") - .def("coordinates_to_indices", &rasputin::coordinates_to_indices, "Transform from coordinate space to index space.") - .def("extract_uint8_buffer_values", &rasputin::extract_buffer_values, "Extract raster data for given indices.") - // From solar_position.h - .def("timestamp_solar_position", [] (const double timestamp, const double geographic_latitude, const double geographic_longitude, const double masl) { - using namespace std::chrono; -#ifndef __clang__ - using namespace date; -#endif - const auto tp = sys_days{January/1/1970} + seconds(long(std::round(timestamp))); - return rasputin::solar_position::time_point_solar_position(tp, geographic_latitude, geographic_longitude, masl, - rasputin::solar_position::collectors::azimuth_and_elevation(), - rasputin::solar_position::delta_t_calculator::coarse_timestamp_calc()); - }, "Compute azimuth and elevation of sun for given UTC timestamp.") - .def("calendar_solar_position", [] (unsigned int year,unsigned int month, double day, const double geographic_latitude, const double geographic_longitude, const double masl) { - return rasputin::solar_position::calendar_solar_position(year, month, day, geographic_latitude, geographic_longitude, masl, - rasputin::solar_position::collectors::azimuth_and_elevation(), - rasputin::solar_position::delta_t_calculator::coarse_date_calc()); - }, "Compute azimuth and elevation of sun for given UT calendar coordinate.") - .def("solar_elevation_correction", &rasputin::solar_position::corrected_solar_elevation, "Correct elevation based on pressure and temperature.") - .def("shade", [] (const rasputin::Mesh& mesh, const double timestamp) { - using namespace std::chrono; -#ifndef __clang__ - using namespace date; -#endif - const auto secs = seconds(int(std::round(timestamp))); - const auto millisecs = milliseconds(int(round(1000*fmod(timestamp, 1)))); - const auto tp = sys_days{January / 1 / 1970} + secs + millisecs; - const auto shade_vec = rasputin::shade(mesh, tp); - py::array_t result(shade_vec.size()); - auto info = result.request(); - auto *data = static_cast(info.ptr); - std::copy(std::begin(shade_vec), std::end(shade_vec), data); - return result; - }) - ; +// static constexpr void make_bindings() { + // PYBIND10_MODULE(triangulate_dem, m) { + // py::bind_vector>(m, "points", py::buffer_protocol()) + // .def_buffer(&vecarray_buffer::value>); + // py::bind_vector(m, "point2_vector", py::buffer_protocol()) + // .def_buffer(&vecarray_buffer); + // py::bind_vector(m, "face_vector", py::buffer_protocol()) + // .def_buffer(&vecarray_buffer); + // py::bind_vector(m, "index_vector", py::buffer_protocol()) + // .def_buffer(&vecarray_buffer); + // py::bind_vector(m, "double_vector", py::buffer_protocol()) + // .def_buffer(&vector_buffer); + + // py::bind_vector>(m, "int_vector"); + // py::bind_vector>>(m, "shadow_vector"); + // m + // .def( + // "lindstrom_turk_by_size", + // [] (const rasputin::point3_vector& raster_coordinates, size_t result_mesh_size) { + // return rasputin::make_tin(raster_coordinates) // use own implementation + // }, + // "Construct a TIN based on the points provided.\n\nThe LindstromTurk cost and placement strategy is used, and simplification process stops when the number of undirected edges drops below the size threshold." + // ) + // .def( + // "lindstrom_turk_by_ratio", + // [] (const rasputin::point3_vector& raster_coordinates, double ratio) { + // return rasputin::make_tin(raster_coordinates) // use own implementation + // }, + // "Construct a TIN based on the points provided.\n\nThe LindstromTurk cost and placement strategy is used, and simplification process stops when the number of undirected edges drops below the ratio threshold." + // ) + // .def("compute_shadow", &rasputin::compute_shadow, "Compute shadows for given sun ray direction.") + // .def("compute_shadows", &rasputin::compute_shadows, "Compute shadows for a series of times and ray directions.") + // .def("surface_normals", &rasputin::surface_normals, "Compute surface normals for all faces in the mesh.") + // .def("point_normals", &rasputin::point_normals, "Compute surface normals for all vertices in the mesh.") + // .def("orient_tin", &rasputin::orient_tin, "Orient all triangles in the TIN and returns their surface normals.") + // .def("extract_lakes", &rasputin::extract_lakes, "Extract lakes as separate face list.") + // .def("compute_slopes", &rasputin::compute_slopes,"Compute slopes (i.e. angles relative to xy plane) for the all the vectors in list.") + // .def("compute_aspects", &rasputin::compute_aspects, "Compute aspects for the all the vectors in list.") + // .def("cell_centers", &rasputin::cell_centers, "Compute cell centers for triangulation.") + // .def("extract_avalanche_expositions", &rasputin::extract_avalanche_expositions, "Extract avalanche exposed cells.") + // .def("coordinates_to_indices", &rasputin::coordinates_to_indices, "Transform from coordinate space to index space.") + // .def("extract_uint8_buffer_values", &rasputin::extract_buffer_values, "Extract raster data for given indices.") + // .def( + // "rasterdata_to_pointvector", + // [] (py::array_t array, double x0, double y1, double dx, double dy) { + // rasputin::point3_vector raster_coordinates; + // auto buffer = array.request(); + // unsigned long M = (unsigned long)buffer.shape[0]; + // unsigned long N = (unsigned long)buffer.shape[1]; + // double* ptr = (double*) buffer.ptr; + // + // //double dx = (x1 - x0)/(n - 1), dy = (y1 - y0)/(m - 1); + // + // raster_coordinates.reserve(M*N); + // for (std::size_t i = 0; i < M; ++i) + // for (std::size_t j = 0; j < N; ++j) + // raster_coordinates.emplace_back(std::array{x0 + j*dx, y1 - i*dy, ptr[i*N + j]}); + // + // return raster_coordinates; + // }, + // "point3_vector from raster data" + // ); + // } +// } } diff --git a/src/rasputin/solar_position.h b/src/rasputin/solar_position.h index e6d8caf..38aa2d3 100644 --- a/src/rasputin/solar_position.h +++ b/src/rasputin/solar_position.h @@ -49,7 +49,7 @@ auto coarse_timestamp_calc() { using namespace date; #endif - auto ep = sys_days(January/1/1970); + auto ep = sys_days{January/1/1970}; return [ep] (system_clock::time_point tp) { const auto year = round(duration_cast(tp - ep).count()*1.0/duration_cast(years(1)).count()) + 1970; diff --git a/src/rasputin/triangulate_dem.h b/src/rasputin/triangulate_dem.h index 3dc346d..9056567 100644 --- a/src/rasputin/triangulate_dem.h +++ b/src/rasputin/triangulate_dem.h @@ -1,916 +1,357 @@ -// -// Created by Ola Skavhaug on 08/10/2018. -// - #pragma once -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include +#include +#include +#include +#include +#include #include +#include #include #include -#include #include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include "solar_position.h" - - - -namespace CGAL { -using K = Exact_predicates_inexact_constructions_kernel; -using Gt = Projection_traits_xy_3; -using Delaunay = Delaunay_triangulation_2; -using ConstrainedDelaunay = Constrained_Delaunay_triangulation_2; - -using Point = Gt::Point_2; -using Point3 = K::Point_3; -using Point2 = K::Point_2; - -using Vector = K::Vector_3; -using PointList = std::vector; -using Mesh = Surface_mesh; -using VertexIndex = Mesh::Vertex_index; -using FaceIndex = Mesh::Face_index; -using PointVertexMap = std::map; -using Ray = K::Ray_3; -using Primitive = CGAL::AABB_face_graph_triangle_primitive; -using Traits = CGAL::AABB_traits; -using Tree = CGAL::AABB_tree; -using Ray_intersection = boost::optional::Type>; -using face_descriptor = boost::graph_traits::face_descriptor; - -using PointSequence = std::vector; -using DelaunayConstraints = std::vector; - -using SimplePolygon = Polygon_2; -using Polygon = Polygon_with_holes_2; -using MultiPolygon = std::vector; - -bool point_inside_polygon(const Point2 &x, const SimplePolygon &polygon) { - return polygon.has_on_bounded_side(x); -} - -bool point_inside_polygon(const Point2 &x, const Polygon &polygon) { - if (not polygon.outer_boundary().has_on_bounded_side(x)) - return false; - if (not polygon.has_holes()) - return true; - for (auto it = polygon.holes_begin(); it != polygon.holes_end(); ++it) - if (it->has_on_bounded_side(x)) - return false; - return true; -} - -bool point_inside_polygon(const Point2 &x, const MultiPolygon &polygon) { - for (const auto& part: polygon) - if (point_inside_polygon(x, part)) - return true; - return false; -} - -std::vector extract_boundaries(const SimplePolygon &polygon) { - std::vector ret; - ret.emplace_back(polygon); - - return ret; -} -std::vector extract_boundaries(const Polygon &polygon) { - std::vector ret; - ret.emplace_back(polygon.outer_boundary()); - for (auto it = polygon.holes_begin(); it != polygon.holes_end(); ++it) - ret.emplace_back(*it); - - return ret; -} - -std::vector extract_boundaries(const MultiPolygon &polygon) { - std::vector ret; - for (const auto& part: polygon) { - ret.emplace_back(part.outer_boundary()); - for (auto it = part.holes_begin(); it != part.holes_end(); ++it) - ret.emplace_back(*it); - } - - return ret; -} -} - -namespace boost::geometry::traits { - - template<> struct tag - { typedef point_tag type; }; - - template<> struct coordinate_type - { typedef double type; }; - - template<> struct coordinate_system - {typedef cs::cartesian type; }; - - template<> struct dimension : boost::mpl::int_<2> {}; - - template <> - struct access - { - static double get(CGAL::Point2 const& p) { return p.x(); } - static void set(CGAL::Point2& p, double const& value) { p = CGAL::Point2{value, p.y()}; } - }; - - template <> - struct access - { - static double get(CGAL::Point2 const& p) { return p.y(); } - static void set(CGAL::Point2& p, double const& value) { p = CGAL::Point2{p.x(), value}; } - }; - - template<> struct tag - { typedef point_tag type; }; - - template<> struct coordinate_type - { typedef double type; }; - - template<> struct coordinate_system - {typedef cs::cartesian type; }; - - template<> struct dimension : boost::mpl::int_<3> {}; - - template <> - struct access - { - static double get(CGAL::Point3 const& p) { return p.x(); } - static void set(CGAL::Point3& p, double const& value) { p = CGAL::Point3{value, p.y(), p.z()}; } - }; - - template <> - struct access - { - static double get(CGAL::Point3 const& p) { return p.y(); } - static void set(CGAL::Point3& p, double const& value) { p = CGAL::Point3{p.x(), value, p.z()}; } - }; - - template <> - struct access - { - static double get(CGAL::Point3 const& p) { return p.z(); } - static void set(CGAL::Point3& p, double const& value) { p = CGAL::Point3{p.x(), p.y(), value}; } - }; -} - - namespace rasputin { -using point3 = std::array; -using point3_vector = std::vector>; -using point2 = std::array; -using point2_vector= std::vector; -using face = std::array; -using index = std::array; -using face_vector = std::vector; -using index_vector = std::vector; -using double_vector = std::vector; -using uint8_vector = std::vector; - -// Clean up below -using VertexIndexMap = std::map; -using FaceDescrMap = std::map; - -CGAL::Mesh construct_mesh(const point3_vector &pts, - const face_vector &faces, - VertexIndexMap& index_map, - FaceDescrMap& face_map); -struct Mesh { - const CGAL::Mesh cgal_mesh; - const std::string proj4_str; - Mesh(CGAL::Mesh cgal_mesh, const std::string proj4_str) - : cgal_mesh(cgal_mesh), proj4_str(proj4_str) {set_points_faces();} - - template - Mesh coarsen(const S& stop, const P& placement, const C& cost) const { - CGAL::Mesh new_cgal_mesh = CGAL::Mesh(this->cgal_mesh); - CGAL::Surface_mesh_simplification::edge_collapse(new_cgal_mesh, - stop, - CGAL::parameters::get_cost(cost) - .get_placement(placement)); - return Mesh(new_cgal_mesh, proj4_str); - } - - Mesh copy() const { - // Call CGAL::Mesh copy constructor to do a deep copy - return Mesh(CGAL::Mesh(cgal_mesh), proj4_str); - } - - const point3_vector& get_points() const { - return points; - } - - const face_vector& get_faces() const { - return faces; - } - - size_t num_edges() const {return cgal_mesh.number_of_edges();} - size_t num_vertices() const {return cgal_mesh.number_of_vertices();} - size_t num_faces() const {return cgal_mesh.number_of_faces();} - - Mesh extract_sub_mesh(const std::vector &face_indices) const { - std::map remap; - point3_vector new_points; - face_vector new_faces; - int counter = 0; - for (auto face_idx: face_indices) { - std::array new_face; - int i = 0; - for (auto idx: faces[face_idx]) { - if (remap.count(idx) == 0) { - remap[idx] = counter++; - new_points.emplace_back(points[idx]); - } - new_face[i++] = remap[idx]; - } - new_faces.emplace_back(new_face); - } - VertexIndexMap index_map; - FaceDescrMap face_map; - return Mesh(construct_mesh(new_points, new_faces, index_map, face_map), proj4_str); - } - - private: - point3_vector points; - face_vector faces; - - void set_points_faces() { - points.clear(); - faces.clear(); - - points.reserve(cgal_mesh.num_vertices()); - faces.reserve(cgal_mesh.num_faces()); - - int n = 0; - std::map reindex; - for (auto f: cgal_mesh.faces()) { - std::array fl; - size_t idx = 0; - for (auto v: cgal_mesh.vertices_around_face(cgal_mesh.halfedge(f))) { - if (reindex.count(v) == 0) { - reindex.emplace(v, n++); - const auto pt = cgal_mesh.point(v); - points.emplace_back(point3{pt.x(), pt.y(), pt.z()}); - } - fl[idx++] = reindex[v]; - } - faces.emplace_back(face{fl[0], fl[1], fl[2]}); - } - } -}; - -CGAL::Mesh construct_mesh(const point3_vector &pts, - const face_vector &faces, - VertexIndexMap& index_map, - FaceDescrMap& face_map){ - CGAL::Mesh mesh; - index_map.clear(); - face_map.clear(); - size_t i = 0; - size_t j = 0; - for (auto p: pts) - index_map[i++] = mesh.add_vertex(CGAL::Point(p[0], p[1], p[2])); - for (auto f: faces) - face_map[mesh.add_face(index_map[f[0]], index_map[f[1]], index_map[f[2]])] = j++; - return mesh; -}; - - -template +namespace bg = boost::geometry; +// namespace bgi = boost::geometry::index; + +template +constexpr bool CLOCKWISE = std::same_as< + std::integral_constant::value>, + std::integral_constant +>; + +template +using DimConst = std::integral_constant; + +template +requires + std::same_as::value>, DimConst<3>> + && std::floating_point + && std::same_as::type, FT> struct RasterData { - RasterData(double x_min, double y_max, double delta_x, double delta_y, - std::size_t num_points_x, std::size_t num_points_y, - FT* data) - : x_min(x_min), delta_x(delta_x), num_points_x(num_points_x), - y_max(y_max), delta_y(delta_y), num_points_y(num_points_y), - data(data) {} - - double x_min; - double delta_x; + using MultiPoint = bg::model::multi_point; + using LineString = bg::model::linestring; + using Box = bg::model::box; + using Polygon = bg::model::polygon; + + FT x_min; + FT x_max; + FT delta_x; std::size_t num_points_x; - double y_max; - double delta_y; + FT y_min; + FT y_max; + FT delta_y; std::size_t num_points_y; - FT* data; - - double get_x_max() const {return x_min + (num_points_x - 1)*delta_x; } - double get_y_min() const {return y_max - (num_points_y - 1)*delta_y; } - - CGAL::PointList raster_points() const { - CGAL::PointList points; - points.reserve(num_points_x * num_points_y); + MultiPoint points; + + RasterData( + FT x_min, FT y_max, FT delta_x, FT delta_y, std::size_t num_points_x, std::size_t num_points_y, FT* data + ) : + x_min(x_min), + x_max(x_min + (num_points_x - 1) * delta_x), + y_min(y_max - (num_points_y - 1) * delta_y), + y_max(y_max), + delta_x(delta_x), + delta_y(delta_y), + num_points_x(num_points_x), + num_points_y(num_points_y) + { + for (std::size_t i = 0; i < num_points_x; ++i) { + for (std::size_t j = 0; j < num_points_y; ++j) { + bg::append( + this->points, + Point(x_min + i * delta_x, y_max - j * delta_y, data[j * num_points_x + i]) + ); + } + } + } - for (std::size_t i = 0; i < num_points_y; ++i) - for (std::size_t j = 0; j < num_points_x; ++j) - points.emplace_back(x_min + j*delta_x, y_max - i*delta_y, data[i*num_points_x + j]); + const MultiPoint& get_points() const { return points; } - // For every point inside the raster rectangle we identify indices (i, j) of the upper-left vertex of the cell containing the point - std::pair get_indices(double x, double y) const { - int i = std::min(std::max(static_cast((y_max-y) /delta_y), 0), num_points_y - 1); - int j = std::min(std::max(static_cast((x-x_min) /delta_x), 0), num_points_x - 1); + // For every point inside the raster rectangle we identify indices (i, j) of + // the upper-left vertex of the cell containing the point + std::pair get_indices(FT x, FT y) const { + int i = std::min(std::max(static_cast((y_max - y) / delta_y), 0), num_points_y - 1); + int j = std::min(std::max(static_cast((x - x_min) / delta_x), 0), num_points_x - 1); - return std::make_pair(i,j); + return std::make_pair(i, j); } // Interpolate data using using a bilinear interpolation rule on each cell - FT get_interpolated_value_at_point(double x, double y) const { + FT get_interpolated_value_at_point(FT x, FT y) const { // Determine indices of the cell containing (x, y) auto [i, j] = get_indices(x, y); // Determine the cell corners // (x0, y0) -- upper left // (x1, y1) -- lower right - double x_0 = x_min + j * delta_x, - y_0 = y_max - i * delta_y, - x_1 = x_min + (j+1) * delta_x, - y_1 = y_max - (i+1) * delta_y; + FT x_0 = x_min + j * delta_x, y_0 = y_max - i * delta_y; + FT x_1 = x_min + (j + 1) * delta_x, y_1 = y_max - (i + 1) * delta_y; - // Using bilinear interpolation on the celll - double h = data[(i + 0)*num_points_x + j + 0] * (x_1 - x)/delta_x * (y - y_1)/delta_y // (x0, y0) - + data[(i + 0)*num_points_x + j + 1] * (x - x_0)/delta_x * (y - y_1)/delta_y // (x1, y0) - + data[(i + 1)*num_points_x + j + 0] * (x_1 - x)/delta_x * (y_0 - y)/delta_y // (x0, y1) - + data[(i + 1)*num_points_x + j + 1] * (x - x_0)/delta_x * (y_0 - y)/delta_y; // (x1, y1) + // Using bilinear interpolation on the cell + FT h = points[(i + 0) * num_points_x + j + 0] * (x_1 - x) / delta_x * (y - y_1) / delta_y // (x0, y0) + + points[(i + 0) * num_points_x + j + 1] * (x - x_0) / delta_x * (y - y_1) / delta_y // (x1, y0) + + points[(i + 1) * num_points_x + j + 0] * (x_1 - x) / delta_x * (y_0 - y) / delta_y // (x0, y1) + + points[(i + 1) * num_points_x + j + 1] * (x - x_0) / delta_x * (y_0 - y) / delta_y; // (x1, y1) return h; } - // Boundary of the raster domain as a CGAL polygon - CGAL::SimplePolygon exterior() const { - CGAL::SimplePolygon rectangle; - rectangle.push_back(CGAL::Point2(x_min, get_y_min())); - rectangle.push_back(CGAL::Point2(get_x_max(), get_y_min())); - rectangle.push_back(CGAL::Point2(get_x_max(), y_max)); - rectangle.push_back(CGAL::Point2(x_min, y_max)); - + // Boundary of the raster domain as a boost geometry polygon + Polygon exterior() const { + Polygon rectangle{{x_min, y_min}, {x_max, y_min}, {x_max, y_max}, {x_min, y_max}}; return rectangle; } // Compute intersection of raster rectangle with polygon - template - CGAL::MultiPolygon compute_intersection(const P& polygon) const { - CGAL::SimplePolygon rectangle = exterior(); + template + requires std::same_as + Polygon compute_intersection(const P &polygon) const { + Box rectangle = exterior(); - CGAL::MultiPolygon intersection_polygon; - CGAL::intersection(rectangle, polygon, std::back_inserter(intersection_polygon)); + Polygon intersection_polygon; + bg::intersection(rectangle, polygon, intersection_polygon); return intersection_polygon; } - // Determine if a point (x, y) is is strictly inside the raster domain - bool contains(double x, double y) const { - double eps = pow(pow(delta_x, 2) + pow(delta_y, 2), 0.5) * 1e-10; - return ((x > x_min + eps) and (x < get_x_max() - eps) - and (y > get_y_min() + eps) and (y < y_max - eps)); + // Determine if a point (x, y) is strictly inside the raster domain + bool contains(FT x, FT y) const { + FT eps = pow(pow(delta_x, 2) + pow(delta_y, 2), 0.5) * 1e-10; + return ((x > x_min + eps) and (x < x_max - eps) and (y > y_min + eps) and (y < y_max - eps)); } }; -template -CGAL::DelaunayConstraints interpolate_boundary_points(const RasterData& raster, - const P& boundary_polygon) { - // First we need to determine intersection points between the raster domain and the polygon - CGAL::MultiPolygon intersection_polygon = raster.compute_intersection(boundary_polygon); - - // Iterate over edges of the intersection polygon and interpolate points - // TODO: Handle holes. Only needed if the intersecting polygon has holes. - CGAL::DelaunayConstraints interpolated_points; - for (const auto& part : CGAL::extract_boundaries(intersection_polygon)) { - for (auto e = part.edges_begin(); e != part.edges_end(); ++e) { - CGAL::Point2 first_vertex = e->vertex(0); - CGAL::Point2 second_vertex = e->vertex(1); - - // We can skip edges that are aligend with the raster boundary - // TODO: Consider checking using CGAL and exact arithmentic - bool edge_is_aligned = not raster.contains((first_vertex.x() + second_vertex.x())/2, - (first_vertex.y() + second_vertex.y())/2); - if (edge_is_aligned) { - continue; +template < + typename Point, + typename FT, + typename MultiPoint=bg::model::multi_point, + typename Segment=bg::model::segment, + typename LineString=bg::model::linestring, + typename MultiLineString=bg::model::multi_linestring, + typename Box=bg::model::box, + typename Polygon=bg::model::polygon, + typename MultiPolygon=bg::model::multi_polygon +> +requires std::same_as::value>, DimConst<3>> && std::floating_point +class Triangulation { +public: + using PointPtr = const Point*; + using Face = std::tuple; + using Faces = std::vector; + using Edge = std::pair; + + struct Mesh { + private: + const std::string proj4_str; + + public: + const Faces triangles; + + explicit Mesh(Faces& triangles, const std::string proj4_str) : + proj4_str(proj4_str), + triangles(triangles) + { } - // Sample with the approximately the same resolution as the raster data along the boundary edges - double edge_len_x = second_vertex.x() - first_vertex.x(); - double edge_len_y = second_vertex.y() - first_vertex.y(); - std::size_t num_subedges = static_cast(std::max(std::fabs(edge_len_x/raster.delta_x), - std::fabs(edge_len_y/raster.delta_y))); - num_subedges = std::max(1, num_subedges); - - double edge_dx = edge_len_x / num_subedges; // signed distance - double edge_dy = edge_len_y / num_subedges; - - // Iterate over edge samples - std::vector interpolated_points_on_edge; - for (std::size_t k=0; k < num_subedges + 1; ++k) { - double x = first_vertex.x() + k * edge_dx; - double y = first_vertex.y() + k * edge_dy; - double z = raster.get_interpolated_value_at_point(x, y); - interpolated_points_on_edge.emplace_back(x, y, z); + // TODO: fix iterators? + size_t num_vertices() const { + return this->triangles.size(); } - interpolated_points.push_back(std::move(interpolated_points_on_edge)); - } - } - return interpolated_points; -} - - -template -Mesh make_mesh(const CGAL::PointList &pts, - const Pgn& inclusion_polygon, - const CGAL::DelaunayConstraints &constraints, - const std::string proj4_str) { - - CGAL::ConstrainedDelaunay dtin; - for (const auto p: pts) - dtin.insert(p); - - for (auto point_sequence: constraints) - dtin.insert_constraint(point_sequence.begin(), point_sequence.end(), false); - CGAL::Mesh mesh; - CGAL::PointVertexMap pvm; + size_t num_edges() const { + return this->triangles.size(); + } - auto index = [&pvm, &mesh] (CGAL::Point v) -> CGAL::VertexIndex { - // Find index of point in mesh, adding it if needed - auto iter = pvm.find(v); - if (iter == pvm.end()) - iter = pvm.emplace(std::make_pair(v, mesh.add_vertex(v))).first; - return iter->second; + size_t num_faces() const { + return this->triangles.size(); + } }; - for (auto f = dtin.finite_faces_begin(); f != dtin.finite_faces_end(); ++f) { - CGAL::Point u = f->vertex(0)->point(); - CGAL::Point v = f->vertex(1)->point(); - CGAL::Point w = f->vertex(2)->point(); - - // Add face if midpoint is contained - CGAL::Point2 face_midpoint(u.x()/3 + v.x()/3 + w.x()/3, - u.y()/3 + v.y()/3 + w.y()/3); - if (CGAL::point_inside_polygon(face_midpoint, inclusion_polygon)) { - mesh.add_face(index(u), index(v), index(w)); + inline static bool _is_inside_ccirc(const PointPtr& point, const Face& triangle) { + const auto [a, b, c] = triangle; + const FT px = bg::get<0>(point); + const FT py = bg::get<1>(point); + const FT pxa = bg::get<0>(a) - px; + const FT pya = bg::get<1>(a) - py; + const FT pxb = bg::get<0>(b) - px; + const FT pyb = bg::get<1>(b) - py; + const FT pxc = bg::get<0>(c) - px; + const FT pyc = bg::get<1>(c) - py; + + // |pxa pya a2 pxa*pxa + pya*pya| + // |pxb pyb b2 pxb*pxb + pyb*pyb| + // |pxc pyc c2 pxc*pxc + pyc*pyc| + // const FT det = pxa*pyb*c2 - pxa*b2*pyc - pya*pxb*c2 + pya*b2*pxc + a2*pxb*pyc - a2*pyb*pxc; + const FT det = + (pxa*pxa + pya*pya)*(pxb*pyc - pxc*pyb) + - (pxb*pxb + pyb*pyb)*(pxa*pyc - pxc*pya) + + (pxc*pxc + pyc*pyc)*(pxa*pyb - pxb*pya); + + if constexpr(CLOCKWISE) { + return det < 0 ? true : false; + } else { + return det > 0 ? true : false; } } - pvm.clear(); - return Mesh(mesh, proj4_str); -}; - - -Mesh make_mesh(const CGAL::PointList &pts, const std::string proj4_str) { - - CGAL::ConstrainedDelaunay dtin; - for (const auto p: pts) - dtin.insert(p); - - CGAL::PointVertexMap pvm; - CGAL::Mesh mesh; - - for (auto v = dtin.finite_vertices_begin(); v != dtin.finite_vertices_end(); ++v) { - CGAL::Point2 pt(v->point().x(), v->point().y()); - pvm.emplace(std::make_pair(v->point(), mesh.add_vertex(v->point()))); - } - for (auto f = dtin.finite_faces_begin(); f != dtin.finite_faces_end(); ++f) { - CGAL::Point u = f->vertex(0)->point(); - CGAL::Point v = f->vertex(1)->point(); - CGAL::Point w = f->vertex(2)->point(); - mesh.add_face(pvm[u], pvm[v], pvm[w]); - } - pvm.clear(); + inline static Face _sort_points(const PointPtr& a, const PointPtr& b, const PointPtr& c) { + using tri_type = std::pair; - return Mesh(mesh, proj4_str); -}; + const FT ax = bg::get<0>(*a); + const FT ay = bg::get<1>(*a); + const FT bx = bg::get<0>(*b); + const FT by = bg::get<1>(*b); + const FT cx = bg::get<0>(*c); + const FT cy = bg::get<1>(*c); + FT ccw = (bx - ax)*(cy - ay) - (cx - ax)*(by - ay); -template -Mesh mesh_from_raster(const std::vector>& raster_list, - const Pgn& boundary_polygon, - const std::string proj4_str) { - CGAL::PointList raster_points; - CGAL::DelaunayConstraints boundary_points; - for (auto raster : raster_list) { - CGAL::PointList new_points = raster.raster_points(); - raster_points.insert(raster_points.end(), - std::make_move_iterator(new_points.begin()), - std::make_move_iterator(new_points.end())); - CGAL::DelaunayConstraints new_constraints = interpolate_boundary_points(raster, boundary_polygon); - boundary_points.insert(boundary_points.end(), - std::make_move_iterator(new_constraints.begin()), - std::make_move_iterator(new_constraints.end())); - } - return make_mesh(raster_points, boundary_polygon, boundary_points, proj4_str); -} - - -template -Mesh mesh_from_raster(const std::vector>& raster_list, const std::string proj4_str) { - CGAL::PointList raster_points; - for (auto raster : raster_list) { - CGAL::PointList new_points = raster.raster_points(); - raster_points.insert(raster_points.end(), - std::make_move_iterator(new_points.begin()), - std::make_move_iterator(new_points.end())); + if constexpr(CLOCKWISE) { + return ccw < 0 ? Face{a, c, b} : Face{a, b, c}; + } else { + return ccw > 0 ? Face{a, c, b} : Face{a, b, c}; + } } - return make_mesh(raster_points, proj4_str); -} - - -template -Mesh mesh_from_raster(const RasterData& raster, - const Pgn& boundary_polygon, - const std::string proj4_str) { - CGAL::PointList raster_points = raster.raster_points(); - CGAL::DelaunayConstraints boundary_points = interpolate_boundary_points(raster, boundary_polygon); - - return make_mesh(raster_points, boundary_polygon, boundary_points, proj4_str); -} + template + static Mesh delaunay(const RasterData& raster, const std::string proj4_str) { + auto add_new_triangle = [](Faces& triangles, const Edge& edge, const PointPtr& point) { + /* add new triangles formed by 2 points in edge and point. Sort points to be clockwise with respect to + centroid */ + Face triangle = _sort_points(edge.first, edge.second, point); + + triangles.push_back(triangle); + }; + + // add super triangle + const FT lower_x = raster.x_min; + const FT lower_y = raster.y_min; + const FT upper_x = raster.x_max; + const FT upper_y = raster.y_max; + + const FT diff_max = std::max(upper_x - lower_x, upper_y - lower_y); + const FT mid_x = 0.5 * (upper_x + lower_x); + const FT mid_y = 0.5 * (upper_y + lower_y); + + Faces triangles; + const Point sta = {mid_x - DN*diff_max, mid_y - diff_max}; + const Point stb = {mid_x, mid_y + DN*diff_max}; + const Point stc = {mid_x + DN*diff_max, mid_y - diff_max}; + const Face super_triangle = std::make_tuple(&sta, &stb, &stc); + triangles.push_back(super_triangle); + + const MultiPoint& points = raster.get_points(); + + add_new_triangle(triangles, {std::get<0>(super_triangle), std::get<1>(super_triangle)}, &(points[0])); + add_new_triangle(triangles, {std::get<1>(super_triangle), std::get<2>(super_triangle)}, &(points[0])); + add_new_triangle(triangles, {std::get<2>(super_triangle), std::get<0>(super_triangle)}, &(points[0])); + + for (typename MultiPoint::const_iterator point=points.begin() + 1; point != points.end(); ++point) { + auto start_bad = std::remove_if( + triangles.begin(), + triangles.end(), + [&point](const Face& triangle) -> bool { + return _is_inside_ccirc(&*point, triangle); + } + ); + + std::vector new_edges{}; + for (typename Faces::iterator triangle1=start_bad; triangle1 != triangles.end(); ++triangle1) { + const auto [a1, b1, c1] = *triangle1; + for (const Edge e1 : {std::make_pair(a1, b1), std::make_pair(b1, c1), std::make_pair(c1, a1)}) { + bool edge_good = true; + for ( + typename Faces::iterator triangle2=triangle1 + 1; + triangle2 != triangles.end() && edge_good; + ++triangle2 + ) { + const auto [a2, b2, c2] = *triangle2; + for (const Edge e2 : { + std::make_pair(a2, b2), std::make_pair(b2, c2), std::make_pair(c2, a2) + }) { + if ( + ((e1.first == e2.first) && (e1.second == e2.second)) + || ((e1.second == e2.first) && (e1.first == e2.second)) + ) { + edge_good = false; + break; + } + } + } + + if (edge_good) { + new_edges.push_back(e1); + } + } + } -template -Mesh mesh_from_raster(const RasterData& raster, const std::string proj4_str) { - return make_mesh(raster.raster_points(), proj4_str); -} + triangles.erase(start_bad, triangles.end()); -CGAL::Point3 centroid(const Mesh& mesh, const CGAL::face_descriptor &face) { - CGAL::Point3 c{0, 0, 0}; - for (auto v: mesh.cgal_mesh.vertices_around_face(mesh.cgal_mesh.halfedge(face))) { - const auto pt = mesh.cgal_mesh.point(v); - c = CGAL::Point3{c.x() + pt.x(), c.y() + pt.y(), c.z() + pt.z()}; - } - return CGAL::Point3{c.x()/3.0, c.y()/3.0, c.z()/3.0}; -} - -std::vector compute_shadow(const Mesh & mesh, - const point3 &sun_direction) { - std::vector shade; - auto cgal_mesh = mesh.cgal_mesh; - const CGAL::Vector sun_vec(-sun_direction[0], -sun_direction[1], -sun_direction[2]); - - int i = 0; - for (auto fd: cgal_mesh.faces()) { - auto v = CGAL::Polygon_mesh_processing::compute_face_normal(fd, cgal_mesh); - if ( v[0]*sun_vec[0] + v[1]*sun_vec[1] + v[2]*sun_vec[2] > 0.0 ) - shade.emplace_back(i); - else { - CGAL::Tree tree(CGAL::faces(cgal_mesh).first, CGAL::faces(cgal_mesh).second, cgal_mesh); - const auto c = centroid(mesh, fd); - CGAL::Ray sun_ray(c, -sun_vec); - auto intersection = tree.first_intersection(sun_ray, - [fd] (const CGAL::face_descriptor &t) { return (t == fd); }); - if (intersection) - shade.emplace_back(i); + for (const Edge& edge : new_edges) { + add_new_triangle(triangles, edge, &*point); + } } - ++i; - } - return shade; -}; -bool is_shaded(const CGAL::Tree &tree, - const CGAL::face_descriptor &fd, - const CGAL::Vector &face_normal, - const CGAL::Point3 &face_center, - const double azimuth, - const double elevation) { - if (elevation < 0.0) - return true; - const arma::vec::fixed<3> sd = arma::normalise(arma::vec::fixed<3>{sin(azimuth*M_PI/180.0), - cos(azimuth*M_PI/180.0), - tan(elevation*M_PI/180.0)}); - const CGAL::Vector sun_vec(-sd[0], -sd[1], -sd[2]); - if ( face_normal[0]*sun_vec[0] + face_normal[1]*sun_vec[1] + face_normal[2]*sun_vec[2] > 0.0 ) - return true; - CGAL::Ray sun_ray(face_center, -sun_vec); - if (tree.first_intersection(sun_ray, [fd] (const CGAL::face_descriptor &t) { return (t == fd); })) - return true; - return false; -} - -auto shade(const Mesh &mesh, - const std::chrono::system_clock::time_point tp) { - std::vector shade_vec; - shade_vec.reserve(mesh.num_faces()); - namespace bg = boost::geometry; - using point_car = bg::model::point; - using point_geo = bg::model::point>; - bg::srs::transformation<> tr{ - bg::srs::proj4(mesh.proj4_str), - bg::srs::epsg(4326) - }; - const CGAL::Tree tree(CGAL::faces(mesh.cgal_mesh).first, CGAL::faces(mesh.cgal_mesh).second, mesh.cgal_mesh); - for (auto fd: mesh.cgal_mesh.faces()) { - const auto c = centroid(mesh, fd); - const point_car x_car{c.x(), c.y()}; - point_geo x_geo; - tr.forward(x_car, x_geo); - const auto lat = bg::get<1>(x_geo); - const auto lon = bg::get<0>(x_geo); - const auto [azimuth, elevation] = solar_position::time_point_solar_position( - tp, - lat, - lon, - c.z(), - rasputin::solar_position::collectors::azimuth_and_elevation(), - rasputin::solar_position::delta_t_calculator::coarse_timestamp_calc() - ); - auto face_normal = CGAL::Polygon_mesh_processing::compute_face_normal(fd, mesh.cgal_mesh); - shade_vec.emplace_back(is_shaded(tree, fd, face_normal, c, azimuth, elevation)); - } - return shade_vec; -} - -std::vector compute_shadow(const Mesh &mesh, - const double azimuth, - const double elevation) { - // Topocentric azimuth and elevation - const arma::vec::fixed<3> sd = arma::normalise(arma::vec::fixed<3>{sin(azimuth*M_PI/180.0), - cos(azimuth*M_PI/180.0), - tan(elevation*M_PI/180.0)}); - return compute_shadow(mesh, point3{sd[0], sd[1], sd[2]}); -}; + // remove all triangles connected to super triangle + std::erase_if( + triangles, + [&super_triangle](const Face& triangle) -> bool { + const auto [a, b, c] = triangle; + const auto [s1, s2, s3] = super_triangle; + + if ( + (a == s1) + || (a == s2) + || (a == s3) + || (b == s1) + || (b == s2) + || (b == s3) + || (c == s1) + || (c == s2) + || (c == s3) + ) { + return true; + } -std::vector> compute_shadows(const Mesh &mesh, - const std::vector> & sun_rays) { - // TODO: Rewrite totally! - std::vector> result; - return result; - /* - result.reserve(sun_rays.size()); - CGAL::Mesh cgal_mesh = mesh.cgal_mesh; - std::map index_map; - std::map face_map; - size_t i = 0; - size_t j = 0; - CGAL::Tree tree(CGAL::faces(cgal_mesh).first, CGAL::faces(cgal_mesh).second, cgal_mesh); - for (auto item: sun_rays) { - std::vector shade; - shade.reserve(cgal_mesh.num_faces()); - auto utc_time = item.first; - const CGAL::Vector sun_ray(item.second[0], - item.second[1], - item.second[2]); - - for (auto fd: CGAL::faces(cgal_mesh)) { - auto hd = halfedge(fd, cgal_mesh); - auto p = CGAL::centroid(cgal_mesh.point(source(hd, cgal_mesh)), - cgal_mesh.point(target(hd, cgal_mesh)), - cgal_mesh.point(target(next(hd, cgal_mesh), cgal_mesh))); - auto v = CGAL::Polygon_mesh_processing::compute_face_normal(fd, cgal_mesh); - if ( v[0]*sun_ray[0] + v[1]*sun_ray[1] + v[2]*sun_ray[2] >= 0.0 ) - shade.emplace_back(face_map[fd]); - else { - CGAL::Ray ray_towards_sun(p, -sun_ray); - auto intersection = tree.first_intersection(ray_towards_sun, - [fd] (const CGAL::face_descriptor &t) { - return (t == fd); - }); - if (intersection) - shade.emplace_back(face_map[fd]); + return false; } - } - result.emplace_back(std::move(shade)); - } - return result; - */ -}; - -point3_vector orient_tin(const point3_vector &pts, face_vector &faces) { - point3_vector result; - result.reserve(faces.size()); - for (auto& face: faces) { - // Compute ccw normal - const auto p0 = pts[face[0]]; - const auto p1 = pts[face[1]]; - const auto p2 = pts[face[2]]; - const arma::vec::fixed<3> v0{p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2]}; - const arma::vec::fixed<3> v1{p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2]}; - const arma::vec::fixed<3> n = arma::cross(v0, v1); - double c = arma::norm(n); - - // Reverse triangle orientation if it is negatively oriented relative to xy plane - if (n[2] < 0.0) { - c *= -1.0; - std::reverse(face.begin(), face.end()); - } + ); - // Store normalised and correctly oriented normal vector - result.push_back(point3{n[0]/c, n[1]/c, n[2]/c}); + return Mesh(triangles, proj4_str); } - return result; -}; -double compute_slope(const point3 &normal) { - return std::atan2(pow(pow(normal[0], 2) + pow(normal[1], 2), 0.5), normal[2]); -} - -double_vector compute_slopes(const point3_vector &normals) { - double_vector result; - result.reserve(normals.size()); - - for (const auto &n : normals) - result.push_back(compute_slope(n)); - - return result; -}; + template + requires std::same_as + Mesh DelaunayConstrained(const RasterData &raster, P &bound); -double compute_aspect(const point3 &normal) { - return std::atan2(normal[0], normal[1]); -} + template + requires std::same_as + Mesh DelaunayConstrained(const RasterData &raster, P &bound); -double_vector compute_aspects(const point3_vector &normals) { - double_vector result; - result.reserve(normals.size()); - for (const auto &n: normals) - result.emplace_back(compute_aspect(n)); - return result; -}; + template + requires std::same_as + Mesh DelaunayConstrained(const RasterData &raster, P &bound); -point3 normal(const point3 &p0, const point3 &p1, const point3 &p2) { - const arma::vec::fixed<3> v0{p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2]}; - const arma::vec::fixed<3> v1{p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2]}; - arma::vec::fixed<3> n = arma::cross(v0, v1); - n /= arma::norm(n); - return n[2] >= 0.0 ? point3{n[0], n[1], n[2]} : point3{-n[0], -n[1], -n[2]}; -} - -point3_vector surface_normals(const point3_vector &pts, const face_vector &faces) { - point3_vector result; - result.reserve(faces.size()); - for (const auto face: faces) - result.emplace_back(normal(pts[face[0]], pts[face[1]], pts[face[2]])); - return result; + template + requires std::same_as + Mesh DelaunayConstrained(const RasterData &raster, P &bound); + + // template + // requires std::same_as && std::floating_point + // std::vector interpolate_boundary_points(const RasterData& raster, const P& boundary_polygon) { MultiPolygon intersection_polygon + // = raster.compute_intersection(boundary_polygon); + // + // std::vector interpolated_points; + // + // return interpolated_points; + // } }; - -template void iadd(T& v, const T& o) { - v[0] += o[0]; - v[1] += o[1]; - v[2] += o[2]; -} - -point3_vector point_normals(const point3_vector &pts, const face_vector &faces) { - point3_vector result(pts.size(), {0.0, 0.0, 0.0}); - for (auto face: faces) { - const auto p0 = pts[face[0]]; - const auto p1 = pts[face[1]]; - const auto p2 = pts[face[2]]; - const arma::vec::fixed<3> v0{p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2]}; - const arma::vec::fixed<3> v1{p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2]}; - arma::vec::fixed<3> n = arma::cross(v0, v1); - n /= arma::norm(n); - const auto v = (n[2] >= 0.0) ? point3{n[0], n[1], n[2]} : point3{-n[0], -n[1], -n[2]}; - iadd(result[face[0]], v); - iadd(result[face[1]], v); - iadd(result[face[2]], v); - } - for (int i = 0; i < result.size(); ++i) { - point3 &p = result[i]; - const double norm = std::sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]); - if (norm > 1.0e-16) { - p[0] /= norm; - p[1] /= norm; - p[2] /= norm; - } - } - return result; -} - -template -std::tuple partition(const point3_vector &pts, - const face_vector &faces, - CB criterion) { - face_vector part1, part2; - for (auto face: faces) { - if (criterion(pts[face[0]], pts[face[1]], pts[face[2]])) - part1.emplace_back(face); - else - part2.emplace_back(face); - } - return std::make_pair(std::move(part1), std::move(part2)); -} - -std::tuple extract_lakes(const point3_vector&pts, const face_vector &faces) { - return partition(pts, faces, [] (const point3 &p0, const point3 &p1, const point3 &p2){ - return compute_slope(normal(p0, p1, p2)) < 1.0e-2; - }); -} - -std::tuple extract_avalanche_expositions(const point3_vector &pts, - const face_vector &faces, - const point2_vector &exposed_intervals, - const point2_vector &height_intervals){ - return partition(pts, faces, [exposed_intervals, height_intervals](const point3 &p0, const point3 &p1, const point3 &p2){ - const auto max_height = std::max(p0[2], std::max(p1[2], p2[2])); - const auto min_height = std::min(p0[2], std::min(p1[2], p2[2])); - bool inside = false; - for (auto height_interval: height_intervals) { - if ((max_height <= height_interval[1] && max_height >= height_interval[0]) || - (min_height <= height_interval[1] && min_height >= height_interval[0])) - inside = true; - } - if (not inside) - return false; - const auto cell_normal = normal(p0, p1, p2); - const auto cell_slope = compute_slope(cell_normal); - if (cell_slope < 30./180.*M_PI) - return false; - const auto aspect = compute_aspect(cell_normal); - for (auto exposition: exposed_intervals) { - if ((exposition[0] < aspect) && (aspect < exposition[1])) - return true; - else if ((exposition[0] > exposition[1]) && ((exposition[0] < aspect) || (aspect < exposition[1]))) - return true; - } - return false; - }); -} - -point3_vector cell_centers(const point3_vector& points, const face_vector & faces) { - point3_vector result; - result.reserve(faces.size()); - for (auto f: faces) { - auto p0 = points[f[0]]; - auto p1 = points[f[1]]; - auto p2 = points[f[2]]; - auto x = (p0[0] + p1[0] + p2[0])/3.0; - auto y = (p0[1] + p1[1] + p2[1])/3.0; - auto z = (p0[2] + p1[2] + p2[2])/3.0; - result.emplace_back(point3{x, y, z}); - } - return result; -} - -index_vector coordinates_to_indices(double x0, - double y1, - double dx, - double dy, - unsigned int M, - unsigned int N, - point2_vector pts) { - - index_vector indices; - indices.reserve(pts.size()); - for (auto pt: pts) - indices.emplace_back(std::array{(unsigned int)((pt[0]-x0)/dx), (unsigned int)((y1 - pt[1])/dy)}); - return indices; -} - -template -std::vector extract_buffer_values(const index_vector& indices, pybind11::array_t& array) { - std::vector result; - result.reserve(indices.size()); - auto buffer = array.request(); - unsigned long M = (unsigned long)buffer.shape[0]; - unsigned long N = (unsigned long)buffer.shape[1]; - T* ptr = (T *)buffer.ptr; - for (auto idx: indices) - result.emplace_back(ptr[idx[0]*N + idx[1]]); // TODO: Implement range check? - return result; -} - -std::tuple consolidate(const point3_vector &points, const face_vector &faces){ - face_vector new_faces; - point3_vector new_points; - new_faces.reserve(faces.size()); - std::map point_map; - int n = 0; - for (const auto _face: faces) { - face new_face; - int i = 0; - for (const auto f: _face) { - if (not point_map.count(f)) { - point_map.insert(std::make_pair(f, n++)); - new_points.emplace_back(points[f]); - } - new_face[i++] = point_map[f]; - } - new_faces.emplace_back(new_face); - } - return std::make_pair(std::move(new_points), std::move(new_faces)); -} -} +}; // namespace rasputin From b33abc8ddb60ac7ca9fc6f106035182535fd312f Mon Sep 17 00:00:00 2001 From: Alocias Date: Tue, 29 Nov 2022 16:32:47 +0000 Subject: [PATCH 4/4] add bindings for Point and points vec types for 3D double --- cpp_test/test_triangulate.cpp | 55 +++++++-- src/rasputin/bindings.cpp | 192 ++++++++++++++++++------------ src/rasputin/triangulate_dem.h | 211 ++++++++++++++++++--------------- 3 files changed, 279 insertions(+), 179 deletions(-) diff --git a/cpp_test/test_triangulate.cpp b/cpp_test/test_triangulate.cpp index aeaf018..1facc4d 100644 --- a/cpp_test/test_triangulate.cpp +++ b/cpp_test/test_triangulate.cpp @@ -175,16 +175,57 @@ namespace rasputin::test_utils { } TEST_CASE("[Test unconstrained Delaunay triangulation]", "[Triangulation]") { - std::vector raw_data = {0, 0, 0, 0, 0, 0, 0, 0}; - RasterData data = {0.0, 1.0, 1.0, 1.0, 2, 2, raw_data.data()}; + SECTION("square") { + std::vector raw_data = {0, 0, 0, 0, 0, 0, 0, 0}; + RasterData data = {0.0, 1.0, 1.0, 1.0, 2, 2, raw_data.data()}; - auto mesh = Triangulation::delaunay(data, ""); + auto mesh = Triangulation::delaunay(data, ""); - check_points_exists(mesh.triangles, data.get_points()); - check_points_inside_ccirc(mesh.triangles, data.get_points()); - check_triangles_unique(mesh.triangles); + check_points_exists(mesh.triangles, data.get_points()); + check_points_inside_ccirc(mesh.triangles, data.get_points()); + check_triangles_unique(mesh.triangles); - REQUIRE(mesh.num_faces() == 2); + REQUIRE(mesh.num_faces() == 2); + } + + SECTION("rectangle 1") { + std::vector raw_data = {0, 0, 0, 0, 0, 0, 0, 0}; + RasterData data = {0.0, 1.0, 2.0, 1.0, 2, 2, raw_data.data()}; + + auto mesh = Triangulation::delaunay(data, ""); + + check_points_exists(mesh.triangles, data.get_points()); + check_points_inside_ccirc(mesh.triangles, data.get_points()); + check_triangles_unique(mesh.triangles); + + REQUIRE(mesh.num_faces() == 2); + } + + SECTION("rectangle 2") { + std::vector raw_data = {0, 0, 0, 0, 0, 0, 0, 0}; + RasterData data = {0.0, 2.0, 1.0, 2.0, 2, 2, raw_data.data()}; + + auto mesh = Triangulation::delaunay(data, ""); + + check_points_exists(mesh.triangles, data.get_points()); + check_points_inside_ccirc(mesh.triangles, data.get_points()); + check_triangles_unique(mesh.triangles); + + REQUIRE(mesh.num_faces() == 2); + } + + SECTION("trapezoid 1") { + std::vector raw_data = {0, 0, 0, 0, 0, 0, 0, 0}; + RasterData data = {0.0, 3.0, 1.0, 3.0, 2, 2, raw_data.data()}; + + auto mesh = Triangulation::delaunay(data, ""); + + check_points_exists(mesh.triangles, data.get_points()); + check_points_inside_ccirc(mesh.triangles, data.get_points()); + check_triangles_unique(mesh.triangles); + + REQUIRE(mesh.num_faces() == 2); + } } // TEST_CASE("[Test random unconstrained Delaunay triangulation]", "[Triangulation]") { diff --git a/src/rasputin/bindings.cpp b/src/rasputin/bindings.cpp index f7bc40e..03e6abe 100644 --- a/src/rasputin/bindings.cpp +++ b/src/rasputin/bindings.cpp @@ -1,107 +1,147 @@ #include +#include #include +#include +#include +#include #include -#include #include +#include +#include #include "triangulate_dem.h" namespace py = pybind11; +namespace bg = boost::geometry; +using Point = bg::model::point; +namespace boost::geometry::traits { + using shared_type = std::shared_ptr; + BOOST_GEOMETRY_DETAIL_SPECIALIZE_POINT_TRAITS(shared_type, 3, double, cs::cartesian) + template<> + struct access { + static inline double get(shared_type const& p) { return p->get<0>(); } + }; + template<> + struct access { + static inline double get(shared_type const& p) { return p->get<1>(); } + }; + template<> + struct access { + static inline double get(shared_type const& p) { return p->get<2>(); } + }; +} + namespace rasputin { -template +template requires std::floating_point -py::buffer_info vecarray_buffer(std::vector> &v) { +py::buffer_info vecarray_buffer(typename RasterData::MultiPoint &v) { return py::buffer_info( &v[0], /* Pointer to buffer */ sizeof(T), /* Size of one scalar */ py::format_descriptor::format(), /* Python struct-style format descriptor */ - 2, /* Number of dimensions */ + bg::dimension(), /* Number of dimensions */ std::vector { v.size(), n }, /* Buffer dimensions */ { sizeof(T) * n, sizeof(T) } /* Strides (in bytes) for each index */ ); } +template +requires std::floating_point +std::vector>* vecarray_from_numpy(py::array_t buf) { + auto info = buf.request(); + if (info.ndim < 1 or info.ndim > 2) + throw py::type_error("Can only convert one- and two-dimensional arrays."); + + // Make sure total size can be preserved + auto size = info.shape[0]; + if (info.ndim == 1 and size % n) + throw py::type_error("Size of one-dimensional array must be divisible by!" + std::to_string(n) + "."); + + if (info.ndim == 2 and info.shape[1] != n) + throw py::type_error("Second dimension does not have size equal to " + std::to_string(n) + "."); + + if (info.ndim == 2) + size *= n; + + auto vec = std::unique_ptr>> (new std::vector> ()); + vec->reserve(size / n); + + // Copy values + double* p = static_cast(info.ptr); + double* end = p + size; + for (; p < end; p += n) { + auto a = std::array {}; + for (auto k = 0; k < n; k++) + a[k] = *(p + k); + vec->push_back(a); + } + + return vec.release(); +} + template requires std::floating_point py::buffer_info vector_buffer(std::vector &v) { return py::buffer_info(&v[0], sizeof(T), py::format_descriptor::format(), 1, { v.size() }, { sizeof(T) }); } -#define MAKE_OPAQUE(T, Point) (PYBIND11_MAKE_OPAQUE(rasputin::Points)) -#define BIND_MODULE(T, Point) {\ - PYBIND10_MODULE(triangulate_dem, m) {\ - py::bind_vector>(m, "points", py::buffer_protocol())\ - .def_buffer(&vecarray_buffer::value>);\ - py::bind_vector>(m, "int_vector");\ - py::bind_vector>>(m, "shadow_vector");\ - };\ +template +void bind_raster_data(py::module& m, const std::string& pyname) { + using RasterData = rasputin::RasterData; + py::class_>(m, pyname.c_str(), py::buffer_protocol()) + .def(py::init([] (py::array_t& data_array, double x_min, double y_max, double delta_x, double delta_y) { + auto buffer = data_array.request(); + int m = buffer.shape[0], n = buffer.shape[1]; + + return RasterData(x_min, y_max, delta_x, delta_y, n, m, static_cast(buffer.ptr) ); + }), py::return_value_policy::take_ownership, py::keep_alive<1, 2>(), + py::arg("data_array").noconvert(), py::arg("x_min"), py::arg("y_max"), py::arg("delta_x"), py::arg("delta_y")) + .def_buffer([] (RasterData& self) { + return py::buffer_info( + (void*)(&(self.get_points()[0])), + sizeof(FT), + py::format_descriptor::format(), + bg::dimension(), + std::vector { self.get_num_y(), self.get_num_y() }, + { sizeof(FT) * self.get_num_x(), sizeof(FT) } + ); + }) + .def_property_readonly("x_max", &RasterData::get_xmax) + .def_property_readonly("y_min", &RasterData::get_ymin) + .def("__getitem__", [](RasterData& self, std::pair idx) { + auto [i, j] = idx; + return self.get_points()[self.get_num_x() * i + j]; }) + .def("get_indices", &RasterData::get_indices) + .def("exterior", &RasterData::exterior) + .def("contains", &RasterData::contains) + .def("get_interpolated_value_at_point", &RasterData::get_interpolated_value_at_point); } -// static constexpr void make_bindings() { - // PYBIND10_MODULE(triangulate_dem, m) { - // py::bind_vector>(m, "points", py::buffer_protocol()) - // .def_buffer(&vecarray_buffer::value>); - // py::bind_vector(m, "point2_vector", py::buffer_protocol()) - // .def_buffer(&vecarray_buffer); - // py::bind_vector(m, "face_vector", py::buffer_protocol()) - // .def_buffer(&vecarray_buffer); - // py::bind_vector(m, "index_vector", py::buffer_protocol()) - // .def_buffer(&vecarray_buffer); - // py::bind_vector(m, "double_vector", py::buffer_protocol()) - // .def_buffer(&vector_buffer); - - // py::bind_vector>(m, "int_vector"); - // py::bind_vector>>(m, "shadow_vector"); - // m - // .def( - // "lindstrom_turk_by_size", - // [] (const rasputin::point3_vector& raster_coordinates, size_t result_mesh_size) { - // return rasputin::make_tin(raster_coordinates) // use own implementation - // }, - // "Construct a TIN based on the points provided.\n\nThe LindstromTurk cost and placement strategy is used, and simplification process stops when the number of undirected edges drops below the size threshold." - // ) - // .def( - // "lindstrom_turk_by_ratio", - // [] (const rasputin::point3_vector& raster_coordinates, double ratio) { - // return rasputin::make_tin(raster_coordinates) // use own implementation - // }, - // "Construct a TIN based on the points provided.\n\nThe LindstromTurk cost and placement strategy is used, and simplification process stops when the number of undirected edges drops below the ratio threshold." - // ) - // .def("compute_shadow", &rasputin::compute_shadow, "Compute shadows for given sun ray direction.") - // .def("compute_shadows", &rasputin::compute_shadows, "Compute shadows for a series of times and ray directions.") - // .def("surface_normals", &rasputin::surface_normals, "Compute surface normals for all faces in the mesh.") - // .def("point_normals", &rasputin::point_normals, "Compute surface normals for all vertices in the mesh.") - // .def("orient_tin", &rasputin::orient_tin, "Orient all triangles in the TIN and returns their surface normals.") - // .def("extract_lakes", &rasputin::extract_lakes, "Extract lakes as separate face list.") - // .def("compute_slopes", &rasputin::compute_slopes,"Compute slopes (i.e. angles relative to xy plane) for the all the vectors in list.") - // .def("compute_aspects", &rasputin::compute_aspects, "Compute aspects for the all the vectors in list.") - // .def("cell_centers", &rasputin::cell_centers, "Compute cell centers for triangulation.") - // .def("extract_avalanche_expositions", &rasputin::extract_avalanche_expositions, "Extract avalanche exposed cells.") - // .def("coordinates_to_indices", &rasputin::coordinates_to_indices, "Transform from coordinate space to index space.") - // .def("extract_uint8_buffer_values", &rasputin::extract_buffer_values, "Extract raster data for given indices.") - // .def( - // "rasterdata_to_pointvector", - // [] (py::array_t array, double x0, double y1, double dx, double dy) { - // rasputin::point3_vector raster_coordinates; - // auto buffer = array.request(); - // unsigned long M = (unsigned long)buffer.shape[0]; - // unsigned long N = (unsigned long)buffer.shape[1]; - // double* ptr = (double*) buffer.ptr; - // - // //double dx = (x1 - x0)/(n - 1), dy = (y1 - y0)/(m - 1); - // - // raster_coordinates.reserve(M*N); - // for (std::size_t i = 0; i < M; ++i) - // for (std::size_t j = 0; j < N; ++j) - // raster_coordinates.emplace_back(std::array{x0 + j*dx, y1 - i*dy, ptr[i*N + j]}); - // - // return raster_coordinates; - // }, - // "point3_vector from raster data" - // ); - // } -// } +template +void bind_raster_list(py::module &m, const std::string& pyname) { + using RasterData = rasputin::RasterData; + py::class_, std::unique_ptr>> (m, pyname.c_str()) + .def(py::init( [] () {std::vector self; return self;})) + .def("add_raster", + [] (std::vector& self, RasterData raster_data) { + self.push_back(raster_data); + }, py::keep_alive<1,2>()) + .def("__getitem__", + [] (std::vector& self, int index) { + return self.at(index); + }, py::return_value_policy::reference_internal); } +} + +PYBIND11_MAKE_OPAQUE(rasputin::RasterData::MultiPoint) +PYBIND11_MODULE(triangulate_dem, m) { + py::bind_vector::MultiPoint>(m, "points", py::buffer_protocol()) + .def_buffer(&rasputin::vecarray_buffer::value>); + py::bind_vector>(m, "int_vector"); + py::bind_vector>>(m, "shadow_vector"); + rasputin::bind_raster_data(m, "raster_data"); + rasputin::bind_raster_list(m, "raster_list"); +}; diff --git a/src/rasputin/triangulate_dem.h b/src/rasputin/triangulate_dem.h index 9056567..0f2acc1 100644 --- a/src/rasputin/triangulate_dem.h +++ b/src/rasputin/triangulate_dem.h @@ -14,7 +14,6 @@ namespace rasputin { namespace bg = boost::geometry; -// namespace bgi = boost::geometry::index; template constexpr bool CLOCKWISE = std::same_as< @@ -31,101 +30,131 @@ requires && std::floating_point && std::same_as::type, FT> struct RasterData { - using MultiPoint = bg::model::multi_point; - using LineString = bg::model::linestring; - using Box = bg::model::box; - using Polygon = bg::model::polygon; - - FT x_min; - FT x_max; - FT delta_x; - std::size_t num_points_x; - - FT y_min; - FT y_max; - FT delta_y; - std::size_t num_points_y; - - MultiPoint points; - - RasterData( - FT x_min, FT y_max, FT delta_x, FT delta_y, std::size_t num_points_x, std::size_t num_points_y, FT* data - ) : - x_min(x_min), - x_max(x_min + (num_points_x - 1) * delta_x), - y_min(y_max - (num_points_y - 1) * delta_y), - y_max(y_max), - delta_x(delta_x), - delta_y(delta_y), - num_points_x(num_points_x), - num_points_y(num_points_y) - { - for (std::size_t i = 0; i < num_points_x; ++i) { - for (std::size_t j = 0; j < num_points_y; ++j) { - bg::append( - this->points, - Point(x_min + i * delta_x, y_max - j * delta_y, data[j * num_points_x + i]) - ); + public: + using MultiPoint = bg::model::multi_point; + using LineString = bg::model::linestring; + using Box = bg::model::box; + using Polygon = bg::model::polygon; + + private: + + FT x_min; + FT x_max; + FT delta_x; + std::size_t num_points_x; + + FT y_min; + FT y_max; + FT delta_y; + std::size_t num_points_y; + + FT* data; + MultiPoint points; + + public: + RasterData( + FT x_min, FT y_max, FT delta_x, FT delta_y, std::size_t num_points_x, std::size_t num_points_y, FT* data + ) : + x_min(x_min), + x_max(x_min + (num_points_x - 1) * delta_x), + y_min(y_max - (num_points_y - 1) * delta_y), + y_max(y_max), + delta_x(delta_x), + delta_y(delta_y), + num_points_x(num_points_x), + num_points_y(num_points_y), + data(data) + { + for (std::size_t i = 0; i < num_points_x; ++i) { + for (std::size_t j = 0; j < num_points_y; ++j) { + bg::append( + this->points, + Point(x_min + i * delta_x, y_max - j * delta_y, data[j * num_points_x + i]) + ); + } } } - } - const MultiPoint& get_points() const { - return points; - } + const MultiPoint& get_points() const { + return points; + } - // For every point inside the raster rectangle we identify indices (i, j) of - // the upper-left vertex of the cell containing the point - std::pair get_indices(FT x, FT y) const { - int i = std::min(std::max(static_cast((y_max - y) / delta_y), 0), num_points_y - 1); - int j = std::min(std::max(static_cast((x - x_min) / delta_x), 0), num_points_x - 1); + const std::size_t& get_num_x() const { + return this->num_points_x; + } - return std::make_pair(i, j); - } + const std::size_t& get_num_y() const { + return this->num_points_y; + } - // Interpolate data using using a bilinear interpolation rule on each cell - FT get_interpolated_value_at_point(FT x, FT y) const { - // Determine indices of the cell containing (x, y) - auto [i, j] = get_indices(x, y); + const FT& get_xmin() const { + return this->x_min; + } - // Determine the cell corners - // (x0, y0) -- upper left - // (x1, y1) -- lower right - FT x_0 = x_min + j * delta_x, y_0 = y_max - i * delta_y; - FT x_1 = x_min + (j + 1) * delta_x, y_1 = y_max - (i + 1) * delta_y; + const FT& get_ymin() const { + return this->y_min; + } - // Using bilinear interpolation on the cell - FT h = points[(i + 0) * num_points_x + j + 0] * (x_1 - x) / delta_x * (y - y_1) / delta_y // (x0, y0) - + points[(i + 0) * num_points_x + j + 1] * (x - x_0) / delta_x * (y - y_1) / delta_y // (x1, y0) - + points[(i + 1) * num_points_x + j + 0] * (x_1 - x) / delta_x * (y_0 - y) / delta_y // (x0, y1) - + points[(i + 1) * num_points_x + j + 1] * (x - x_0) / delta_x * (y_0 - y) / delta_y; // (x1, y1) + const FT& get_xmax() const { + return this->x_max; + } - return h; - } + const FT& get_ymax() const { + return this->y_max; + } - // Boundary of the raster domain as a boost geometry polygon - Polygon exterior() const { - Polygon rectangle{{x_min, y_min}, {x_max, y_min}, {x_max, y_max}, {x_min, y_max}}; - return rectangle; - } + // For every point inside the raster rectangle we identify indices (i, j) of + // the upper-left vertex of the cell containing the point + std::pair get_indices(FT x, FT y) const { + int i = std::min(std::max(static_cast((y_max - y) / delta_y), 0), num_points_y - 1); + int j = std::min(std::max(static_cast((x - x_min) / delta_x), 0), num_points_x - 1); - // Compute intersection of raster rectangle with polygon - template - requires std::same_as - Polygon compute_intersection(const P &polygon) const { - Box rectangle = exterior(); + return std::make_pair(i, j); + } - Polygon intersection_polygon; - bg::intersection(rectangle, polygon, intersection_polygon); + // Interpolate data using using a bilinear interpolation rule on each cell + FT get_interpolated_value_at_point(FT x, FT y) const { + // Determine indices of the cell containing (x, y) + auto [i, j] = get_indices(x, y); - return intersection_polygon; - } + // Determine the cell corners + // (x0, y0) -- upper left + // (x1, y1) -- lower right + const FT x_0 = x_min + j * delta_x, y_0 = y_max - i * delta_y; + const FT x_1 = x_min + (j + 1) * delta_x, y_1 = y_max - (i + 1) * delta_y; - // Determine if a point (x, y) is strictly inside the raster domain - bool contains(FT x, FT y) const { - FT eps = pow(pow(delta_x, 2) + pow(delta_y, 2), 0.5) * 1e-10; - return ((x > x_min + eps) and (x < x_max - eps) and (y > y_min + eps) and (y < y_max - eps)); - } + // Using bilinear interpolation on the cell + const FT h = data[(i + 0) * num_points_x + j + 0] * (x_1 - x) / delta_x * (y - y_1) / delta_y // (x0, y0) + + data[(i + 0) * num_points_x + j + 1] * (x - x_0) / delta_x * (y - y_1) / delta_y // (x1, y0) + + data[(i + 1) * num_points_x + j + 0] * (x_1 - x) / delta_x * (y_0 - y) / delta_y // (x0, y1) + + data[(i + 1) * num_points_x + j + 1] * (x - x_0) / delta_x * (y_0 - y) / delta_y; // (x1, y1) + + return h; + } + + // Boundary of the raster domain as a boost geometry polygon + Polygon exterior() const { + Polygon rectangle{{{x_min, y_min}, {x_max, y_min}, {x_max, y_max}, {x_min, y_max}}}; + return rectangle; + } + + // Compute intersection of raster rectangle with polygon + template + requires std::same_as + Polygon compute_intersection(const P &polygon) const { + Box rectangle = exterior(); + + Polygon intersection_polygon; + bg::intersection(rectangle, polygon, intersection_polygon); + + return intersection_polygon; + } + + // Determine if a point (x, y) is strictly inside the raster domain + bool contains(FT x, FT y) const { + FT eps = pow(pow(delta_x, 2) + pow(delta_y, 2), 0.5) * 1e-10; + return ((x > x_min + eps) and (x < x_max - eps) and (y > y_min + eps) and (y < y_max - eps)); + } }; template < @@ -160,15 +189,6 @@ class Triangulation { { } - // TODO: fix iterators? - size_t num_vertices() const { - return this->triangles.size(); - } - - size_t num_edges() const { - return this->triangles.size(); - } - size_t num_faces() const { return this->triangles.size(); } @@ -188,7 +208,6 @@ class Triangulation { // |pxa pya a2 pxa*pxa + pya*pya| // |pxb pyb b2 pxb*pxb + pyb*pyb| // |pxc pyc c2 pxc*pxc + pyc*pyc| - // const FT det = pxa*pyb*c2 - pxa*b2*pyc - pya*pxb*c2 + pya*b2*pxc + a2*pxb*pyc - a2*pyb*pxc; const FT det = (pxa*pxa + pya*pya)*(pxb*pyc - pxc*pyb) - (pxb*pxb + pyb*pyb)*(pxa*pyc - pxc*pya) @@ -231,10 +250,10 @@ class Triangulation { }; // add super triangle - const FT lower_x = raster.x_min; - const FT lower_y = raster.y_min; - const FT upper_x = raster.x_max; - const FT upper_y = raster.y_max; + const FT& lower_x = raster.get_xmin(); + const FT& lower_y = raster.get_ymin(); + const FT& upper_x = raster.get_xmax(); + const FT& upper_y = raster.get_ymax(); const FT diff_max = std::max(upper_x - lower_x, upper_y - lower_y); const FT mid_x = 0.5 * (upper_x + lower_x);