diff --git a/.gitignore b/.gitignore index 7251aa5..00c17ba 100644 --- a/.gitignore +++ b/.gitignore @@ -2,15 +2,15 @@ cmake-* .idea .env *.swp -lib/CGAL* -lib/pybind11 -lib/armadillo* -src/rasputin/*so -src/rasputin.egg-info +lib/**/* +**/*.egg-info +**/*.so __pycache__ build dist .DS_Store -lib/boo* -lib/geometry -lib/build_* +.envrc +compile_commands.json +.cache/ +.python-version +*.eggs diff --git a/CMakeLists.txt b/CMakeLists.txt index 57508a9..1155419 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,9 +4,17 @@ 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(CMAKE_CXX_FLAGS_RELEASE_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") @@ -54,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. @@ -77,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 @@ -100,4 +108,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/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/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/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_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/cpp_test/test_triangulate.cpp b/cpp_test/test_triangulate.cpp new file mode 100644 index 0000000..1facc4d --- /dev/null +++ b/cpp_test/test_triangulate.cpp @@ -0,0 +1,246 @@ +#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]") { + 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, ""); + + 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 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]") { + // 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.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 567b1d6..310c2bd 100644 --- a/setup.py +++ b/setup.py @@ -1,16 +1,14 @@ 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 +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) @@ -39,10 +40,19 @@ 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' + ] + + if self.debug: + cfg = 'Debug' + build_path = "build/debug" + else: + cfg = 'Release' + build_path = "build/release" - cfg = 'Debug' if self.debug else 'Release' build_args = ['--config', cfg] if platform.system() == "Windows": @@ -53,55 +63,35 @@ 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()) + 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_path] + cmake_args + _compile_and_check(cmake_build_command, env=env) + _compile_and_check(["cmake", "--build", build_path]) + 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', - 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'), - install_requires=[ - 'numpy', - 'pyproj', - 'Pillow', - 'h5py', - 'lxml', - 'shapely', - 'descartes', - 'meshio' - ], - # 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')], + 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']} ) - diff --git a/src/rasputin/bindings.cpp b/src/rasputin/bindings.cpp index d1c2c77..03e6abe 100644 --- a/src/rasputin/bindings.cpp +++ b/src/rasputin/bindings.cpp @@ -1,47 +1,55 @@ -#include "triangulate_dem.h" -#include "solar_position.h" +#include +#include + +#include +#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 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>(); } + }; +} -template -py::buffer_info vecarray_buffer(std::vector> &v) { +namespace rasputin { +template +requires std::floating_point +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) @@ -74,339 +82,66 @@ std::vector>* vecarray_from_numpy(py::array_t buf) { return vec.release(); } -template +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()) +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 rasputin::RasterData(x_min, y_max, delta_x, delta_y, n, m, static_cast(buffer.ptr) ); + 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([] (rasputin::RasterData& self) { + .def_buffer([] (RasterData& self) { return py::buffer_info( - self.data, + (void*)(&(self.get_points()[0])), 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) } + bg::dimension(), + std::vector { self.get_num_y(), self.get_num_y() }, + { sizeof(FT) * self.get_num_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) { + .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.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); + 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); } -template +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;})) + 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, rasputin::RasterData raster_data) { + [] (std::vector& self, RasterData raster_data) { self.push_back(raster_data); }, py::keep_alive<1,2>()) .def("__getitem__", - [] (std::vector>& self, int index) { + [] (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(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::MultiPoint>(m, "points", py::buffer_protocol()) + .def_buffer(&rasputin::vecarray_buffer::value>); 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; - }) - ; -} + rasputin::bind_raster_data(m, "raster_data"); + rasputin::bind_raster_list(m, "raster_list"); +}; diff --git a/src/rasputin/solar_position.h b/src/rasputin/solar_position.h index 109952d..38aa2d3 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 { diff --git a/src/rasputin/triangulate_dem.h b/src/rasputin/triangulate_dem.h index 3dc346d..0f2acc1 100644 --- a/src/rasputin/triangulate_dem.h +++ b/src/rasputin/triangulate_dem.h @@ -1,916 +1,376 @@ -// -// 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; - } +namespace bg = boost::geometry; + +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 { + public: + using MultiPoint = bg::model::multi_point; + using LineString = bg::model::linestring; + using Box = bg::model::box; + using Polygon = bg::model::polygon; - const face_vector& get_faces() const { - return faces; - } + private: - 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]); + 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]) + ); } - 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]}); + const MultiPoint& get_points() const { + return points; } - } -}; - -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; -}; + const std::size_t& get_num_x() const { + return this->num_points_x; + } -template -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; - std::size_t num_points_x; - - double y_max; - double 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); - - 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]); - return points; - } + const std::size_t& get_num_y() const { + return this->num_points_y; + } - // 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); + const FT& get_xmin() const { + return this->x_min; + } - return std::make_pair(i,j); - } + const FT& get_ymin() const { + return this->y_min; + } - // Interpolate data using using a bilinear interpolation rule on each cell - FT get_interpolated_value_at_point(double x, double 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; - - // 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) - - return h; - } + const FT& get_xmax() const { + return this->x_max; + } - // 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)); + const FT& get_ymax() const { + return this->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 - CGAL::MultiPolygon compute_intersection(const P& polygon) const { - CGAL::SimplePolygon rectangle = exterior(); + return std::make_pair(i, j); + } - CGAL::MultiPolygon intersection_polygon; - CGAL::intersection(rectangle, polygon, std::back_inserter(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 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)); - } -}; + // 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) -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; - } + return h; + } - // 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); - } - interpolated_points.push_back(std::move(interpolated_points_on_edge)); + // 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; } - } - return interpolated_points; -} + // Compute intersection of raster rectangle with polygon + template + requires std::same_as + Polygon compute_intersection(const P &polygon) const { + Box rectangle = exterior(); -template -Mesh make_mesh(const CGAL::PointList &pts, - const Pgn& inclusion_polygon, - const CGAL::DelaunayConstraints &constraints, - const std::string proj4_str) { + Polygon intersection_polygon; + bg::intersection(rectangle, polygon, intersection_polygon); - CGAL::ConstrainedDelaunay dtin; - for (const auto p: pts) - dtin.insert(p); + return intersection_polygon; + } - for (auto point_sequence: constraints) - dtin.insert_constraint(point_sequence.begin(), point_sequence.end(), false); + // 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)); + } +}; - CGAL::Mesh mesh; - CGAL::PointVertexMap pvm; +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) + { + } - 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*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); + inline static Face _sort_points(const PointPtr& a, const PointPtr& b, const PointPtr& c) { + using tri_type = std::pair; - CGAL::PointVertexMap pvm; - CGAL::Mesh mesh; + 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); - 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(); + FT ccw = (bx - ax)*(cy - ay) - (cx - ax)*(by - ay); - return Mesh(mesh, proj4_str); -}; - - -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.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); + 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