From 8927b706ca9d99830f70068ace9ef6308e7f3599 Mon Sep 17 00:00:00 2001 From: Alocias Date: Thu, 20 Oct 2022 16:06:11 +0200 Subject: [PATCH 1/2] fix test pipeline, add build-push image running when Dockerfile changes, move tests CmakeLists to root, add Catch2 using lib/Catch2, remove date lib, update deps in Dockerfile, use cmake binary instead of source compile, compile clang+llvm 16 from source, cleanup cpp_tests, add land-cover-code option(defaulting to clc18_kode) --- .github/workflows/main.yaml | 78 +- .gitignore | 16 +- CMakeLists.txt | 46 +- Dockerfile | 41 +- Pipfile | 18 - README.md | 24 +- cpp_test/CMakeLists.txt | 10 - cpp_test/test_sun_position.cpp | 43 +- setup.cfg | 56 ++ setup.py | 80 +- src/rasputin/application.py | 18 +- src/rasputin/bindings.cpp | 412 --------- src/rasputin/cpp/bindings.cpp | 21 + src/rasputin/cpp/mesh.h | 269 ++++++ src/rasputin/cpp/mesh_bindings.cpp | 81 ++ src/rasputin/cpp/polygon.h | 54 ++ src/rasputin/cpp/polygon_bindings.cpp | 159 ++++ src/rasputin/cpp/raster_bindings.cpp | 69 ++ src/rasputin/cpp/raster_data.h | 110 +++ src/rasputin/cpp/shade.h | 156 ++++ src/rasputin/cpp/slope.h | 225 +++++ src/rasputin/cpp/slope_bindings.cpp | 22 + src/rasputin/{ => cpp}/solar_position.h | 541 +++++------ src/rasputin/cpp/solar_position_bindings.cpp | 41 + src/rasputin/cpp/triangulate_dem.h | 83 ++ src/rasputin/cpp/types.h | 63 ++ src/rasputin/cpp/vec_bindings.cpp | 84 ++ src/rasputin/geometry.py | 4 +- src/rasputin/gml_repository.py | 23 +- src/rasputin/triangulate_dem.h | 916 ------------------- tests/test_mesh.py | 6 +- tests/test_polygons.py | 39 +- 32 files changed, 2021 insertions(+), 1787 deletions(-) delete mode 100644 Pipfile delete mode 100644 cpp_test/CMakeLists.txt create mode 100644 setup.cfg delete mode 100644 src/rasputin/bindings.cpp create mode 100644 src/rasputin/cpp/bindings.cpp create mode 100644 src/rasputin/cpp/mesh.h create mode 100644 src/rasputin/cpp/mesh_bindings.cpp create mode 100644 src/rasputin/cpp/polygon.h create mode 100644 src/rasputin/cpp/polygon_bindings.cpp create mode 100644 src/rasputin/cpp/raster_bindings.cpp create mode 100644 src/rasputin/cpp/raster_data.h create mode 100644 src/rasputin/cpp/shade.h create mode 100644 src/rasputin/cpp/slope.h create mode 100644 src/rasputin/cpp/slope_bindings.cpp rename src/rasputin/{ => cpp}/solar_position.h (56%) create mode 100644 src/rasputin/cpp/solar_position_bindings.cpp create mode 100644 src/rasputin/cpp/triangulate_dem.h create mode 100644 src/rasputin/cpp/types.h create mode 100644 src/rasputin/cpp/vec_bindings.cpp delete mode 100644 src/rasputin/triangulate_dem.h diff --git a/.github/workflows/main.yaml b/.github/workflows/main.yaml index a9debb1..436125a 100644 --- a/.github/workflows/main.yaml +++ b/.github/workflows/main.yaml @@ -1,4 +1,3 @@ - name: CI on: @@ -7,14 +6,73 @@ on: pull_request: branches: [ master ] +env: + REGISTRY: ghcr.io + RASPUTIN_DATA_DIR: examples/example_data + jobs: + check_changes: + runs-on: ubuntu-latest + outputs: + run_job: ${{ steps.check_dockerfile.outputs.run_job }} + steps: + - name: Checkout + uses: actions/checkout@v3 + with: + fetch-depth: 2 + + - name: Check Dockerfile + id: check_dockerfile + run: | + if git diff --name-only HEAD^ HEAD | grep -q "^Dockerfile$"; then + echo "run_job=true" >> "$GITHUB_OUTPUT" + else + echo "run_job=false" >> "$GITHUB_OUTPUT" + fi + + build-and-push-image: + runs-on: ubuntu-latest + needs: check_changes + if: ${{ needs.check_changes.outputs.run_job == 'true' }} + permissions: + contents: read + packages: write + + steps: + - uses: actions/checkout@v3 + + - name: Login + uses: docker/login-action@v2 + with: + registry: ${{ env.REGISTRY }} + username: ${{ github.actor }} + password: ${{ secrets.GITHUB_TOKEN }} + + - name: Extract metadata + id: meta + uses: docker/metadata-action@v4 + with: + images: ${{ env.REGISTRY }}/${{ github.repository }}_${{ github.ref_name }} + tags: | + type=raw,value=latest,enable=true + + - name: Build and Push + uses: docker/build-push-action@v4 + with: + context: . + push: true + tags: ${{ steps.meta.outputs.tags }} + labels: ${{ steps.meta.outputs.labels }} + build: runs-on: ubuntu-latest + needs: [check_changes, build-and-push-image] + if: ${{ !failure() && !cancelled() }} container: - image: dnumgis/rasputin:latest + image: ghcr.io/${{ github.repository }}_${{ github.ref_name }} steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v3 - name: Copy dependencies run: | @@ -22,11 +80,7 @@ jobs: - name: Install run: | - python3 setup.py install - - - name: Set Data Dir - run: | - echo "::set-env name=RASPUTIN_DATA_DIR::/__w/rasputin/rasputin/examples/example_data" + pip -v install -e . - name: Rasputin Store run: | @@ -34,11 +88,9 @@ jobs: - name: Rasputin Web run: | - rasputin_web -output examples/example_data/web/ingoya_test -material examples/example_data/materials/material.yaml ingoya_test + rasputin_web -output examples/example_data/web/ingoya_test -material examples/example_data/materials/material.yaml ingoya_test - name: Run tests run: | - pytest tests/test_mesh.py - pytest tests/test_polygons.py - pytest tests/test_read_raster_file.py - pytest tests/test_tin_repository.py + pytest tests + ./build/test_sun_position 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..4ee809e 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") @@ -40,15 +47,15 @@ add_definitions(-DHAS_BOOST) # date # ---- # TODO: Check CXX version if this is really needed? -if (EXISTS "${RASPUTIN_LIB_DIR}/date") - add_library(date INTERFACE) - target_include_directories(date INTERFACE - "${RASPUTIN_LIB_DIR}/date/include") -else() - find_package(date REQUIRED) - add_library(date INTERFACE IMPORTED) -endif() -list(APPEND RASPUTIN_DEPENDENCIES date) +# if (EXISTS "${RASPUTIN_LIB_DIR}/date") +# add_library(date INTERFACE) +# target_include_directories(date INTERFACE +# "${RASPUTIN_LIB_DIR}/date/include") +# else() +# find_package(date REQUIRED) +# add_library(date INTERFACE IMPORTED) +# endif() +# list(APPEND RASPUTIN_DEPENDENCIES date) # CGAL @@ -83,21 +90,24 @@ find_package(LAPACK) # Header-only Rasputin library # ---------------------------- -set(RASPUTIN_SOURCE "${CMAKE_CURRENT_SOURCE_DIR}/src/rasputin") +set(RASPUTIN_SOURCE "${CMAKE_CURRENT_SOURCE_DIR}/src/rasputin/cpp") add_library(rasputin INTERFACE) target_include_directories(rasputin PRIVATE INTERFACE ${RASPUTIN_SOURCE}) -if (ENABLE_TESTS) - add_subdirectory(cpp_test) -endif() - # Python bindnigs # --------------- # Use pybind11 as a system library to avoid problems with include order -pybind11_add_module(triangulate_dem MODULE SYSTEM - "${RASPUTIN_SOURCE}/bindings.cpp") +file(GLOB RASPUTIN_SOURCES ${RASPUTIN_SOURCE}/*.cpp) +pybind11_add_module(triangulate_dem MODULE SYSTEM ${RASPUTIN_SOURCES}) # Need to link BLAS for header-only Aramdillo target_link_libraries(triangulate_dem PRIVATE ${RASPUTIN_DEPENDENCIES} - ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES}) + ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES}) + +if (ENABLE_TESTS) + enable_testing() + add_subdirectory(lib/Catch2) + add_executable(test_sun_position cpp_test/test_sun_position.cpp) + target_link_libraries(test_sun_position PRIVATE rasputin Catch2::Catch2WithMain ${RASPUTIN_DEPENDENCIES}) +endif() diff --git a/Dockerfile b/Dockerfile index 3693b5b..f9af7ae 100644 --- a/Dockerfile +++ b/Dockerfile @@ -1,8 +1,10 @@ -FROM ubuntu:20.04 +FROM ubuntu:22.04 WORKDIR /root -# Fetch packages Dependancies +ENV DEBIAN_FRONTEND=noninteractive + +# Fetch packages Dependencies RUN apt-get update &&\ apt-get upgrade -y &&\ apt-get install -y\ @@ -17,7 +19,7 @@ RUN apt-get update &&\ liblapack-dev \ libatlas3-base \ libatlas-base-dev \ - python3.6 \ + python3-pip \ git \ pybind11-dev \ libarmadillo-dev \ @@ -25,24 +27,31 @@ RUN apt-get update &&\ python3-pip \ libcurl4 \ libssl-dev \ - libcurl4-openssl-dev &&\ + libcurl4-openssl-dev \ + libtinfo5 && \ apt-get clean && \ rm -rf /var/lib/apt/lists/* /tmp/* /var/tmp/* -# Python Depndancies -RUN pip3 install numpy pyproj Pillow h5py lxml shapely descartes meshio pytest setuptools pipenv pyyaml requests +# Python Dependencies +RUN python3 -m pip install numpy pyproj Pillow h5py lxml shapely descartes meshio pytest setuptools pipenv pyyaml requests # Cmake -RUN wget https://github.com/Kitware/CMake/releases/download/v3.16.5/cmake-3.16.5.tar.gz \ - && tar -zxvf cmake-3.16.5.tar.gz \ - && rm cmake-3.16.5.tar.gz \ - && cd cmake-3.16.5 \ - && ./bootstrap \ - && make \ - && make install +RUN wget https://github.com/Kitware/CMake/releases/download/v3.25.2/cmake-3.25.2-linux-x86_64.tar.gz \ + && tar -zxf cmake-3.25.2-linux-x86_64.tar.gz \ + && rm cmake-3.25.2-linux-x86_64.tar.gz \ + && cp -R cmake-3.25.2-linux-x86_64/* /usr/ \ + && rm -r cmake-3.25.2-linux-x86_64 + +# clang+llvm 16 +RUN wget https://github.com/llvm/llvm-project/releases/download/llvmorg-16.0.0/clang+llvm-16.0.0-x86_64-linux-gnu-ubuntu-18.04.tar.xz \ + && tar -xf clang+llvm-16.0.0-x86_64-linux-gnu-ubuntu-18.04.tar.xz \ + && rm clang+llvm-16.0.0-x86_64-linux-gnu-ubuntu-18.04.tar.xz \ + && cp -R clang+llvm-16.0.0-x86_64-linux-gnu-ubuntu-18.04/* /usr/ \ + && rm -r clang+llvm-16.0.0-x86_64-linux-gnu-ubuntu-18.04 + +ENV CXX=/usr/bin/clang++ # Github resources RUN mkdir lib && cd lib && \ - git clone https://github.com/boostorg/geometry.git --branch=boost-1.71.0 &&\ - git clone https://github.com/HowardHinnant/date.git &&\ - git clone https://github.com/catchorg/Catch2.git --branch=v2.10.2 catch2 + git clone https://github.com/boostorg/geometry.git --branch=boost-1.81.0 --depth=1 &&\ + git clone https://github.com/catchorg/Catch2.git --branch=v3.3.1 --depth=1 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 deleted file mode 100644 index 94cc2bf..0000000 --- a/cpp_test/CMakeLists.txt +++ /dev/null @@ -1,10 +0,0 @@ -enable_testing() - -# Find Catch2 -# ----------- -find_package(Catch2 3 REQUIRED) - -# Build tests -# ----------- -add_executable(test_sun_position test_sun_position.cpp) -target_link_libraries(test_sun_position PRIVATE rasputin Catch2::Catch2WithMain ${RASPUTIN_DEPENDENCIES}) diff --git a/cpp_test/test_sun_position.cpp b/cpp_test/test_sun_position.cpp index c7966fb..c0a8dfc 100644 --- a/cpp_test/test_sun_position.cpp +++ b/cpp_test/test_sun_position.cpp @@ -1,11 +1,8 @@ #define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one cpp file -//#include +// #include #include #include #include -#ifndef __clang__ -#include -#endif #include #include @@ -35,12 +32,12 @@ TEST_CASE("Reference example test", "[reference]") { const double T = 11; using namespace rasputin::test_utils; using namespace rasputin::solar_position; - 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); + const auto [Phi, e0] = calendar_solar_position( + year, month, day, lat, lon, masl, collectors::azimuth_and_elevation(), fixed_cal_delta_t_calc() + ); + 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]") { @@ -67,16 +64,13 @@ TEST_CASE("UTC cal test 1", "[utccalc]") { using namespace std::chrono; -#ifndef __clang__ - using namespace date; -#endif - const auto [Phi0, elevation0] = calendar_solar_position(year, mon, dom + h/24.0, lat, lon, masl, - collectors::azimuth_and_elevation(), - fixed_cal_delta_t_calc()); - const auto tp = sys_days{month(mon)/dom/year} + hours(h); - const auto [Phi1, elevation1] = time_point_solar_position(tp, lat, lon, masl, - collectors::azimuth_and_elevation(), - fixed_time_point_delta_t_calc()); + const auto [Phi0, elevation0] = calendar_solar_position( + year, mon, dom + h/24.0, lat, lon, masl, collectors::azimuth_and_elevation(), fixed_cal_delta_t_calc() + ); + const auto tp = sys_days{std::chrono::year(year)/month(mon)/dom} + hours(h); + const auto [Phi1, elevation1] = time_point_solar_position( + tp, lat, lon, masl, collectors::azimuth_and_elevation(), fixed_time_point_delta_t_calc() + ); REQUIRE(elevation0 == elevation1); REQUIRE(Phi0 == Phi1); } @@ -95,13 +89,10 @@ TEST_CASE("DeltaT test", "[DT]") { const unsigned int h = 19; using namespace std::chrono; -#ifndef __clang__ - using namespace date; -#endif using namespace rasputin::solar_position; - auto tp = sys_days{month(mon)/dom/year}; - const auto result = time_point_solar_position(tp, lat, lon, masl, - collectors::azimuth_and_elevation(), - delta_t_calculator::coarse_timestamp_calc()); + auto tp = sys_days{std::chrono::year(year)/month(mon)/dom}; + const auto result = time_point_solar_position( + tp, lat, lon, masl, collectors::azimuth_and_elevation(), delta_t_calculator::coarse_timestamp_calc() + ); } diff --git a/setup.cfg b/setup.cfg new file mode 100644 index 0000000..452e689 --- /dev/null +++ b/setup.cfg @@ -0,0 +1,56 @@ +[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 +package_dir= + =src +packages=find: +install_requires = + cython + numpy + pyproj + Pillow + h5py + lxml + shapely + descartes + meshio + pybind11 +setup_requires = + cython + h5py + pybind11 +python_requires=>=3.10 + +[options.packages.find] +where=src + +[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..eb7f827 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,8 +40,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,55 +57,35 @@ def build_extension(self, ext): cmake_args += ['-A', 'x64'] build_args += ['--', '/m'] else: - cmake_args += ['-DCMAKE_BUILD_TYPE=' + cfg] - build_args += ['--', '-j2'] + cmake_args += ['-DCMAKE_BUILD_TYPE=' + cfg, f"-DCMAKE_PREFIX_PATH={Path.home() / '.local'};{Path().cwd() / 'lib'}"] + build_args += ['--', '-j4'] 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"] + cmake_args + _compile_and_check(cmake_build_command, env=env) + _compile_and_check(["cmake", "--build", "build"] + build_args) + 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')], # 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/application.py b/src/rasputin/application.py index 6037460..f315368 100644 --- a/src/rasputin/application.py +++ b/src/rasputin/application.py @@ -18,7 +18,7 @@ from rasputin import globcov_repository -def store_tin(): +def store_tin(args=None): """ Avalance Forecast Visualization Example. @@ -29,10 +29,12 @@ def store_tin(): * Use more of the information from varsom.no (and perhaps alpha blending) to better display avalanche dangers """ + if args is None: + args = sys.argv[:1] + logging.basicConfig(level=logging.CRITICAL, format='Rasputin[%(levelname)s]: %(message)s') logger = logging.getLogger() - arg_parser = argparse.ArgumentParser() arg_parser.add_argument("-x", nargs="+", type=float, help="x-coordinates of polygon", default=None) arg_parser.add_argument("-y", nargs="+", type=float, help="y-coordinates of polygon", default=None) @@ -46,9 +48,13 @@ def store_tin(): default="", choices=["corine", "globcov"], help="Partition mesh by land type") + arg_parser.add_argument( + "--land-cover-code", type=str, default="clc18_kode", + help="Corine LandCover type code (in ogr iteration). Default to clc18_kode" + ) arg_parser.add_argument("uid", type=str, help="Unique ID for the result TIN") arg_parser.add_argument("-silent", action="store_true", help="Run in silent mode") - res = arg_parser.parse_args(sys.argv[1:]) + res = arg_parser.parse_args(args) if not res.silent: logger.setLevel(logging.INFO) @@ -124,14 +130,12 @@ def store_tin(): if res.land_type_partition: if res.land_type_partition == "corine": - lt_repo = gml_repository.GMLRepository(path=corine_archive) + lt_repo = gml_repository.GMLRepository(path=corine_archive, land_cover_code=res.land_cover_code) else: lt_repo = globcov_repository.GlobCovRepository(path=gc_archive) geo_cell_centers = GeoPoints(xy=mesh.cell_centers[:, :2], crs=target_crs) - terrain_cover = lt_repo.land_cover(land_types=None, - geo_points=geo_cell_centers, - domain=target_domain) + terrain_cover = lt_repo.land_cover(land_types=None, geo_points=geo_cell_centers, domain=target_domain) terrain_colors = np.empty((terrain_cover.shape[0], 3), dtype='d') extracted_terrain_types = set() for i, cell in enumerate(terrain_cover): diff --git a/src/rasputin/bindings.cpp b/src/rasputin/bindings.cpp deleted file mode 100644 index d1c2c77..0000000 --- a/src/rasputin/bindings.cpp +++ /dev/null @@ -1,412 +0,0 @@ -#include "triangulate_dem.h" -#include "solar_position.h" -#include -#include -#include -#include - -#include - -// Surface mesh simplication policies -#include -#include -#include -#include -#include -#include -#include -#include - -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>); - -template -py::buffer_info vecarray_buffer(std::vector> &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 */ - std::vector { v.size(), n }, /* Buffer dimensions */ - { sizeof(T) * n, sizeof(T) } /* Strides (in bytes) for each index */ - ); -} - -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 -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); -} - - -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; - }) - ; -} diff --git a/src/rasputin/cpp/bindings.cpp b/src/rasputin/cpp/bindings.cpp new file mode 100644 index 0000000..4aa29ad --- /dev/null +++ b/src/rasputin/cpp/bindings.cpp @@ -0,0 +1,21 @@ +#include + +namespace py = pybind11; + +void bind_vectors(py::module&); +void bind_raster(py::module&); +void bind_mesh(py::module&); +void bind_polygons(py::module&); +void bind_solars(py::module&); +void bind_shades(py::module&); +void bind_slopes(py::module&); + +PYBIND11_MODULE(triangulate_dem, m) { + bind_vectors(m); + bind_raster(m); + bind_mesh(m); + bind_polygons(m); + bind_slopes(m); + bind_shades(m); + bind_solars(m); +} diff --git a/src/rasputin/cpp/mesh.h b/src/rasputin/cpp/mesh.h new file mode 100644 index 0000000..3fe348a --- /dev/null +++ b/src/rasputin/cpp/mesh.h @@ -0,0 +1,269 @@ +#pragma once +#include + +#include "types.h" +#include "raster_data.h" +#include "polygon.h" + +namespace rasputin { +inline 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; +}; + +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]}); + } + } +}; + +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; + } + + // 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)); + } + } + 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; + + 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; + }; + + 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)); + } + } + pvm.clear(); + return Mesh(mesh, proj4_str); +}; + + +inline 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(); + + 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())); + } + 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 +Mesh mesh_from_raster(const RasterData& raster, const std::string proj4_str) { + return make_mesh(raster.raster_points(), proj4_str); +} + +} diff --git a/src/rasputin/cpp/mesh_bindings.cpp b/src/rasputin/cpp/mesh_bindings.cpp new file mode 100644 index 0000000..55252c0 --- /dev/null +++ b/src/rasputin/cpp/mesh_bindings.cpp @@ -0,0 +1,81 @@ +#include + +#include +#include +#include + +// Surface mesh simplication policies +#include +#include +#include +#include + +#include "raster_data.h" +#include "mesh.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(CGAL::MultiPolygon); +PYBIND11_MAKE_OPAQUE(std::vector>); +PYBIND11_MAKE_OPAQUE(std::vector>); + +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); +} + +void bind_mesh(py::module& m) { + 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); + + 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); + + m.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); +} diff --git a/src/rasputin/cpp/polygon.h b/src/rasputin/cpp/polygon.h new file mode 100644 index 0000000..3ca52f8 --- /dev/null +++ b/src/rasputin/cpp/polygon.h @@ -0,0 +1,54 @@ +#pragma once + +#include "types.h" + +namespace CGAL { +inline bool point_inside_polygon(const Point2 &x, const SimplePolygon &polygon) { + return polygon.has_on_bounded_side(x); +} + +inline 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; +} + +inline 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; +} + +inline std::vector extract_boundaries(const SimplePolygon &polygon) { + std::vector ret; + ret.emplace_back(polygon); + + return ret; +} + +inline 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; +} + +inline 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; +} +} diff --git a/src/rasputin/cpp/polygon_bindings.cpp b/src/rasputin/cpp/polygon_bindings.cpp new file mode 100644 index 0000000..c7c6c3d --- /dev/null +++ b/src/rasputin/cpp/polygon_bindings.cpp @@ -0,0 +1,159 @@ +#include +#include + +#include +#include +#include + +#include "polygon.h" + +namespace py = pybind11; + +PYBIND11_MAKE_OPAQUE(CGAL::MultiPolygon); + +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; +} + + +void bind_polygons(py::module& 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); +} + diff --git a/src/rasputin/cpp/raster_bindings.cpp b/src/rasputin/cpp/raster_bindings.cpp new file mode 100644 index 0000000..c861ee1 --- /dev/null +++ b/src/rasputin/cpp/raster_bindings.cpp @@ -0,0 +1,69 @@ +#include +#include +#include + +#include "raster_data.h" + +namespace py = pybind11; + +PYBIND11_MAKE_OPAQUE(std::vector>); +PYBIND11_MAKE_OPAQUE(std::vector>); + +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); +} + +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); +} + +void bind_raster(py::module& m) { + 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"); +} diff --git a/src/rasputin/cpp/raster_data.h b/src/rasputin/cpp/raster_data.h new file mode 100644 index 0000000..872c098 --- /dev/null +++ b/src/rasputin/cpp/raster_data.h @@ -0,0 +1,110 @@ +#pragma once + +#include + +#include "types.h" + +namespace rasputin { +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; + } + + // 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); + + 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 { + // 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; + } + + // 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)); + + return rectangle; + } + + // Compute intersection of raster rectangle with polygon + template + CGAL::MultiPolygon compute_intersection(const P& polygon) const { + CGAL::SimplePolygon rectangle = exterior(); + + CGAL::MultiPolygon intersection_polygon; + CGAL::intersection(rectangle, polygon, std::back_inserter(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)); + } +}; +} diff --git a/src/rasputin/cpp/shade.h b/src/rasputin/cpp/shade.h new file mode 100644 index 0000000..c4011ab --- /dev/null +++ b/src/rasputin/cpp/shade.h @@ -0,0 +1,156 @@ +#include +#include +#include +#include + +#include +#include + +#include + +#include "types.h" +#include "mesh.h" +#include "solar_position.h" + +namespace rasputin { +inline 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}; +} + +inline 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); + } + ++i; + } + return shade; +}; + +inline 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; +} + +inline std::vector 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; +} + +inline 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]}); +}; + +inline 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]); + } + } + result.emplace_back(std::move(shade)); + } + return result; + */ +}; +} diff --git a/src/rasputin/cpp/slope.h b/src/rasputin/cpp/slope.h new file mode 100644 index 0000000..a2d7160 --- /dev/null +++ b/src/rasputin/cpp/slope.h @@ -0,0 +1,225 @@ +#include + +#include + +#include "types.h" + +namespace rasputin { +inline 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 result; +}; + +inline double compute_slope(const point3 &normal) { + return std::atan2(pow(pow(normal[0], 2) + pow(normal[1], 2), 0.5), normal[2]); +} + +inline 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; +}; + +inline double compute_aspect(const point3 &normal) { + return std::atan2(normal[0], normal[1]); +} + +inline 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; +}; + +inline 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]}; +} + +inline 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 void iadd(T& v, const T& o) { + v[0] += o[0]; + v[1] += o[1]; + v[2] += o[2]; +} + +inline 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)); +} + +inline 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; + }); +} + +inline 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; + }); +} + +inline 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; +} + +inline 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; +} + +inline 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)); +} +} diff --git a/src/rasputin/cpp/slope_bindings.cpp b/src/rasputin/cpp/slope_bindings.cpp new file mode 100644 index 0000000..5321b2c --- /dev/null +++ b/src/rasputin/cpp/slope_bindings.cpp @@ -0,0 +1,22 @@ +#include + +#include "slope.h" + +namespace py = pybind11; + +void bind_slopes(py::module& m) { + m.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."); + +} diff --git a/src/rasputin/solar_position.h b/src/rasputin/cpp/solar_position.h similarity index 56% rename from src/rasputin/solar_position.h rename to src/rasputin/cpp/solar_position.h index 109952d..63e3c59 100644 --- a/src/rasputin/solar_position.h +++ b/src/rasputin/cpp/solar_position.h @@ -1,20 +1,20 @@ // // Created by Ola Skavhaug on 26/08/2019. // - #pragma once #include #include -#ifndef __clang__ -#include -#endif #include +#include +#include namespace rasputin::solar_position::collectors { -auto azimuth_and_elevation() { - return [] (double e0, double Gamma, double Phi, const double alpha_mark, const double delta_mark, const double H_mark) { +inline auto azimuth_and_elevation() { + return [] ( + double e0, double Gamma, double Phi, const double alpha_mark, const double delta_mark, const double H_mark + ) { return std::make_tuple(Phi, e0); }; } @@ -24,16 +24,18 @@ auto azimuth_and_elevation() { namespace rasputin::solar_position::delta_t_calculator { const std::vector _dts{31.1, 33.2, 35.7, 40.2, 45.5, 50.5, 54.3, 56.9, 60.8, 63.8, 64.7}; -auto parabolic_dt_calc() { +inline auto parabolic_dt_calc() { return [] (const double year) { return -20 + 32*pow((year - 1820)/20.0, 2); }; } -auto coarse_date_calc() { - return [] (const unsigned int year, - const unsigned int month, - const unsigned int day) { +inline auto coarse_date_calc() { + return [] ( + const unsigned int year, + const unsigned int month, + const unsigned int day + ) { if (year < 1955 || year > 2009) return parabolic_dt_calc()(year); const auto idx = (year - 1955) / 5; @@ -41,13 +43,10 @@ auto coarse_date_calc() { }; } -auto coarse_timestamp_calc() { +inline auto coarse_timestamp_calc() { using namespace std::chrono; -#ifndef __clang__ - using namespace date; -#endif - auto ep = sys_days{January/1/1970}; + auto ep = sys_days{1970y / January / 1}; return [ep] (system_clock::time_point tp) { const auto year = round(duration_cast(tp - ep).count()*1.0/duration_cast(years(1)).count()) + 1970; @@ -62,19 +61,19 @@ auto coarse_timestamp_calc() { namespace rasputin::solar_position { -auto limit_degrees(const double degrees) { +inline double limit_degrees(const double degrees) { return fmod(fmod(degrees, 360) + 360, 360); } -auto d2r(const double d) { +inline double d2r(const double d) { return d*M_PI/180.0; } -auto r2d(const double r) { +inline double r2d(const double r) { return r*180.0/M_PI; } -auto jd_from_cal(unsigned int year, unsigned int month, const double day){ +inline double jd_from_cal(unsigned int year, unsigned int month, const double day){ if (month < 3) { month += 12; @@ -82,40 +81,39 @@ auto jd_from_cal(unsigned int year, unsigned int month, const double day){ } // Equation (4) - const auto jd_temp = trunc(365.25*(year + 4716)) + trunc(30.6001*(month + 1)) + day - 1524.5; + const double jd_temp = trunc(365.25*(year + 4716)) + trunc(30.6001*(month + 1)) + day - 1524.5; if (jd_temp < 2299160) return jd_temp; - const auto A = trunc(year/100); + const double A = trunc(year/100.); return jd_temp + 2 - A + trunc(A/4); } -double jd_from_clock(const std::chrono::system_clock::time_point tp) { +inline double jd_from_clock(const std::chrono::system_clock::time_point tp) { using namespace std::chrono; -#ifndef __clang__ - using namespace date; -#endif - const auto dt = tp - sys_days{January/1/1970}; - const auto t = duration_cast(dt).count(); + + const auto dt = tp - sys_days{1970y / January / 1}; + const long t = duration_cast(dt).count(); + return t/86400.0 + 2440587.5; } -double jde(const double julian_day, const double dt){ +inline double jde(const double julian_day, const double dt){ return julian_day + dt/86400.0; } -double jc(const double julian_day) { +inline double jc(const double julian_day) { return (julian_day - 2451545.0)/36525.0; } -double jce(const double julian_ephemeris_day) { +inline double jce(const double julian_ephemeris_day) { return (julian_ephemeris_day - 2451545.0)/36525.0; } -double jme(const double julian_ephemeris_centry) { +inline double jme(const double julian_ephemeris_centry) { return julian_ephemeris_centry/10.0; } -double heliocentric_longitude(const double julian_ephemeris_millennium) { +inline double heliocentric_longitude(const double julian_ephemeris_millennium) { const double jme = julian_ephemeris_millennium; const double L00 = 175347046.0; const double L01 = 3341656*cos(4.6692568 + 6283.07585 *jme); @@ -280,7 +278,7 @@ double heliocentric_longitude(const double julian_ephemeris_millennium) { return limit_degrees(L_deg); } -double heliocentric_latitude(double julian_ephemeris_millennium) { +inline double heliocentric_latitude(double julian_ephemeris_millennium) { const double jme = julian_ephemeris_millennium; const double B00 = 280*cos(3.199 + 84334.662*jme); @@ -289,7 +287,7 @@ double heliocentric_latitude(double julian_ephemeris_millennium) { const double B03 = 44*cos(3.7 + 2352.87 *jme); const double B04 = 32*cos(4 + 1577.34 *jme); const double B0 = B00 + B01 + B02 + B03 + B04; - + const double B10 = 9*cos(3.9 + 5507.55*jme); const double B11 = 6*cos(1.73 + 5223.69*jme); const double B1 = B10 + B11; @@ -299,7 +297,7 @@ double heliocentric_latitude(double julian_ephemeris_millennium) { return B_deg; } -double heliocentric_radius_vector(double julian_ephemeris_millennium) { +inline double heliocentric_radius_vector(double julian_ephemeris_millennium) { const double jme = julian_ephemeris_millennium; const double R00 = 100013989; @@ -375,130 +373,130 @@ double heliocentric_radius_vector(double julian_ephemeris_millennium) { return (R0 + R1*jme + R2*pow(jme, 2) + R3*pow(jme, 3) + R4*pow(jme, 4))*1e-8; } -double geocentric_longitude(const double heliocentric_longitude) { +inline double geocentric_longitude(const double heliocentric_longitude) { return heliocentric_longitude + 180; } -double geocentric_latitude(const double heliocentric_latitude) { +inline double geocentric_latitude(const double heliocentric_latitude) { return -heliocentric_latitude; } -auto nutation(const double julian_ephemeris_centry) { - const auto jce = julian_ephemeris_centry; - const auto jce2 = jce*jce; - const auto jce3 = jce2*jce; - - const auto X0 = 297.85036 + 445267.111480*jce - 0.0019142*jce2 + jce3/189474; - const auto X1 = 357.52772 + 35999.050340*jce - 0.0001603*jce2 - jce3/300000; - const auto X2 = 134.96298 + 477198.867398*jce + 0.0086972*jce2 + jce3/56250; - const auto X3 = 93.27191 + 483202.017538*jce - 0.0036825*jce2 + jce3/327270; - const auto X4 = 125.04452 - 1934.136261*jce + 0.0020708*jce2 + jce3/450000; - - const auto dksi_00 = (-171996 - 174.2*jce)*sin(d2r( X4 )); - const auto dksi_01 = ( -13187 - 1.6*jce)*sin(d2r(-X0*2 +X3*2 +X4*2)); - const auto dksi_02 = ( -2274 - 0.2*jce)*sin(d2r( X3*2 +X4*2)); - const auto dksi_03 = ( 2062 + 0.2*jce)*sin(d2r( X4*2)); - const auto dksi_04 = ( 1426 - 3.4*jce)*sin(d2r( X1 )); - const auto dksi_05 = ( 712 + 0.1*jce)*sin(d2r( X2 )); - const auto dksi_06 = ( -517 + 1.2*jce)*sin(d2r(-X0*2 +X1 +X3*2 +X4*2)); - const auto dksi_07 = ( -386 - 0.4*jce)*sin(d2r( X3*2 +X4 )); - const auto dksi_08 = -301 *sin(d2r( X2 +X3*2 +X4*2)); - const auto dksi_09 = ( 217 - 0.5*jce)*sin(d2r(-X0*2 -X1 +X3*2 +X4*2)); - const auto dksi_10 = -158 *sin(d2r(-X0*2 +X2 )); - const auto dksi_11 = ( 129 + 0.1*jce)*sin(d2r(-X0*2 +X3*2 +X4 )); - const auto dksi_12 = 123 *sin(d2r( -X2 +X3*2 +X4*2)); - const auto dksi_13 = 63 *sin(d2r( X0*2 )); - const auto dksi_14 = ( 63 + 0.1*jce)*sin(d2r( X2 +X4 )); - const auto dksi_15 = -59 *sin(d2r( X0*2 -X2 +X3*2 +X4*2)); - const auto dksi_16 = ( -58 - 0.1*jce)*sin(d2r( -X2 +X4 )); - const auto dksi_17 = -51 *sin(d2r( X2 +X3*2 +X4 )); - const auto dksi_18 = 48 *sin(d2r(-X0*2 +X2*2 )); - const auto dksi_19 = 46 *sin(d2r( -X2*2 +X3*2 +X4 )); - const auto dksi_20 = -38 *sin(d2r( X0*2 +X3*2 +X4*2)); - const auto dksi_21 = -31 *sin(d2r( X2*2 +X3*2 +X4*2)); - const auto dksi_22 = 29 *sin(d2r( X2*2 )); - const auto dksi_23 = 29 *sin(d2r(-X0*2 +X2 +X3*2 +X4*2)); - const auto dksi_24 = 26 *sin(d2r( X3*2 )); - const auto dksi_25 = -22 *sin(d2r(-X0*2 +X3*2 )); - const auto dksi_26 = 21 *sin(d2r( -X2 +X3*2 +X4 )); - const auto dksi_27 = ( 17 - 0.1*jce)*sin(d2r( X1*2 )); - const auto dksi_28 = 16 *sin(d2r( X0*2 -X2 +X4 )); - const auto dksi_29 = ( -16 + 0.1*jce)*sin(d2r(-X0*2 +X1*2 +X3*2 +X4*2)); - const auto dksi_30 = -15 *sin(d2r( X1 +X4 )); - const auto dksi_31 = -13 *sin(d2r(-X0*2 +X2 +X4 )); - const auto dksi_32 = -12 *sin(d2r( -X1 +X4 )); - const auto dksi_33 = 11 *sin(d2r( X2*2 -X3*2 )); - const auto dksi_34 = -10 *sin(d2r( X0*2 -X2 +X3*2 +X4 )); - const auto dksi_35 = -8 *sin(d2r( X0*2 +X2 +X3*2 +X4*2)); - const auto dksi_36 = 7 *sin(d2r( X1 +X3*2 +X4*2)); - const auto dksi_37 = -7 *sin(d2r(-X0*2 +X1 +X2 )); - const auto dksi_38 = -7 *sin(d2r( -X1 +X3*2 +X4*2)); - const auto dksi_39 = -7 *sin(d2r( X0*2 +X3*2 +X4 )); - const auto dksi_40 = 6 *sin(d2r( X0*2 +X2 )); - const auto dksi_41 = 6 *sin(d2r(-X0*2 +X2*2 +X3*2 +X4*2)); - const auto dksi_42 = 6 *sin(d2r(-X0*2 +X2 +X3*2 +X4 )); - const auto dksi_43 = -6 *sin(d2r( X0*2 -X2*2 +X4 )); - const auto dksi_44 = -6 *sin(d2r( X0*2 +X4 )); - const auto dksi_45 = 5 *sin(d2r( -X1 +X2 )); - const auto dksi_46 = -5 *sin(d2r(-X0*2 -X1 +X3*2 +X4 )); - const auto dksi_47 = -5 *sin(d2r(-X0*2 +X4 )); - const auto dksi_48 = -5 *sin(d2r( X2*2 +X3*2 +X4 )); - const auto dksi_49 = 4 *sin(d2r(-X0*2 +X2*2 +X4 )); - const auto dksi_50 = 4 *sin(d2r(-X0*2 +X1 +X3*2 +X4 )); - const auto dksi_51 = 4 *sin(d2r( X2 -X3*2 )); - const auto dksi_52 = -4 *sin(d2r(-X0 +X2 )); - const auto dksi_53 = -4 *sin(d2r(-X0*2 +X1 )); - const auto dksi_54 = -4 *sin(d2r( X0 )); - const auto dksi_55 = 3 *sin(d2r( X2 +X3*2 )); - const auto dksi_56 = -3 *sin(d2r( -X2*2 +X3*2 +X4*2)); - const auto dksi_57 = -3 *sin(d2r(-X0 -X1 +X2 )); - const auto dksi_58 = -3 *sin(d2r( X1 +X2 )); - const auto dksi_59 = -3 *sin(d2r( -X1 +X2 +X3*2 +X4*2)); - const auto dksi_60 = -3 *sin(d2r( X0*2 -X1 -X2 +X3*2 +X4*2)); - const auto dksi_61 = -3 *sin(d2r( X2*3 +X3*2 +X4*2)); - const auto dksi_62 = -3 *sin(d2r( X0*2 -X1 +X3*2 +X4*2)); - - const auto deps_00 = (92025 + 8.9*jce)*cos(d2r( X4 )); - const auto deps_01 = ( 5736 - 3.1*jce)*cos(d2r(-X0*2 +X3*2 +X4*2)); - const auto deps_02 = ( 977 - 0.5*jce)*cos(d2r( X3*2 +X4*2)); - const auto deps_03 = ( -895 + 0.5*jce)*cos(d2r( X4*2)); - const auto deps_04 = ( 54 - 0.1*jce)*cos(d2r( X1 )); - const auto deps_05 = -7 *cos(d2r( X2 )); - const auto deps_06 = ( 224 - 0.6*jce)*cos(d2r(-X0*2 +X1 +X3*2 +X4*2)); - const auto deps_07 = 200 *cos(d2r( X3*2 +X4 )); - const auto deps_08 = ( 129 - 0.1*jce)*cos(d2r( X2 +X3*2 +X4*2)); - const auto deps_09 = ( -95 + 0.3*jce)*cos(d2r(-X0*2 -X1 +X3*2 +X4*2)); - const auto deps_10 = -70 *cos(d2r(-X0*2 +X3*2 +X4 )); - const auto deps_11 = -53 *cos(d2r( -X2 +X3*2 +X4*2)); - const auto deps_12 = -33 *cos(d2r( X2 +X4 )); - const auto deps_13 = 26 *cos(d2r( X0*2 -X2 +X3*2 +X4*2)); - const auto deps_14 = 32 *cos(d2r( -X2 +X4 )); - const auto deps_15 = 27 *cos(d2r( X2 +X3*2 +X4 )); - const auto deps_16 = -24 *cos(d2r( -X2*2 +X3*2 +X4 )); - const auto deps_17 = 16 *cos(d2r( X0*2 +X3*2 +X4*2)); - const auto deps_18 = 13 *cos(d2r( X2*2 +X3*2 +X4*2)); - const auto deps_19 = -12 *cos(d2r(-X0*2 +X2 +X3*2 +X4*2)); - const auto deps_20 = -10 *cos(d2r( -X2 +X3*2 +X4 )); - const auto deps_21 = -8 *cos(d2r( X0*2 -X2 +X4 )); - const auto deps_22 = 7 *cos(d2r(-X0*2 +X1*2 +X3*2 +X4*2)); - const auto deps_23 = 9 *cos(d2r( X1 +X4 )); - const auto deps_24 = 7 *cos(d2r(-X0*2 +X2 +X4 )); - const auto deps_25 = 6 *cos(d2r( -X1 +X4 )); - const auto deps_26 = 5 *cos(d2r( X0*2 -X2 +X3*2 +X4 )); - const auto deps_27 = 3 *cos(d2r( X0*2 +X2 +X3*2 +X4*2)); - const auto deps_28 = -3 *cos(d2r( X1 +X3*2 +X4*2)); - const auto deps_29 = 3 *cos(d2r( -X1 +X3*2 +X4*2)); - const auto deps_30 = 3 *cos(d2r( X0*2 +X3*2 +X4 )); - const auto deps_31 = -3 *cos(d2r(-X0*2 +X2*2 +X3*2 +X4*2)); - const auto deps_32 = -3 *cos(d2r(-X0*2 +X2 +X3*2 +X4 )); - const auto deps_33 = 3 *cos(d2r( X0*2 -X2*2 +X4 )); - const auto deps_34 = 3 *cos(d2r( X0*2 +X4 )); - const auto deps_35 = 3 *cos(d2r(-X0*2 -X1 +X3*2 +X4 )); - const auto deps_36 = 3 *cos(d2r(-X0*2 +X4 )); - const auto deps_37 = 3 *cos(d2r( X2*2 +X3*2 +X4 )); - - const auto dksi = (dksi_00 + dksi_01 + dksi_02 + dksi_03 + dksi_04 + dksi_05 + dksi_06 + dksi_07 + dksi_08 + dksi_09 +inline std::pair nutation(const double julian_ephemeris_centry) { + const double jce = julian_ephemeris_centry; + const double jce2 = jce*jce; + const double jce3 = jce2*jce; + + const double X0 = 297.85036 + 445267.111480*jce - 0.0019142*jce2 + jce3/189474; + const double X1 = 357.52772 + 35999.050340*jce - 0.0001603*jce2 - jce3/300000; + const double X2 = 134.96298 + 477198.867398*jce + 0.0086972*jce2 + jce3/56250; + const double X3 = 93.27191 + 483202.017538*jce - 0.0036825*jce2 + jce3/327270; + const double X4 = 125.04452 - 1934.136261*jce + 0.0020708*jce2 + jce3/450000; + + const double dksi_00 = (-171996 - 174.2*jce)*sin(d2r( X4 )); + const double dksi_01 = ( -13187 - 1.6*jce)*sin(d2r(-X0*2 +X3*2 +X4*2)); + const double dksi_02 = ( -2274 - 0.2*jce)*sin(d2r( X3*2 +X4*2)); + const double dksi_03 = ( 2062 + 0.2*jce)*sin(d2r( X4*2)); + const double dksi_04 = ( 1426 - 3.4*jce)*sin(d2r( X1 )); + const double dksi_05 = ( 712 + 0.1*jce)*sin(d2r( X2 )); + const double dksi_06 = ( -517 + 1.2*jce)*sin(d2r(-X0*2 +X1 +X3*2 +X4*2)); + const double dksi_07 = ( -386 - 0.4*jce)*sin(d2r( X3*2 +X4 )); + const double dksi_08 = -301 *sin(d2r( X2 +X3*2 +X4*2)); + const double dksi_09 = ( 217 - 0.5*jce)*sin(d2r(-X0*2 -X1 +X3*2 +X4*2)); + const double dksi_10 = -158 *sin(d2r(-X0*2 +X2 )); + const double dksi_11 = ( 129 + 0.1*jce)*sin(d2r(-X0*2 +X3*2 +X4 )); + const double dksi_12 = 123 *sin(d2r( -X2 +X3*2 +X4*2)); + const double dksi_13 = 63 *sin(d2r( X0*2 )); + const double dksi_14 = ( 63 + 0.1*jce)*sin(d2r( X2 +X4 )); + const double dksi_15 = -59 *sin(d2r( X0*2 -X2 +X3*2 +X4*2)); + const double dksi_16 = ( -58 - 0.1*jce)*sin(d2r( -X2 +X4 )); + const double dksi_17 = -51 *sin(d2r( X2 +X3*2 +X4 )); + const double dksi_18 = 48 *sin(d2r(-X0*2 +X2*2 )); + const double dksi_19 = 46 *sin(d2r( -X2*2 +X3*2 +X4 )); + const double dksi_20 = -38 *sin(d2r( X0*2 +X3*2 +X4*2)); + const double dksi_21 = -31 *sin(d2r( X2*2 +X3*2 +X4*2)); + const double dksi_22 = 29 *sin(d2r( X2*2 )); + const double dksi_23 = 29 *sin(d2r(-X0*2 +X2 +X3*2 +X4*2)); + const double dksi_24 = 26 *sin(d2r( X3*2 )); + const double dksi_25 = -22 *sin(d2r(-X0*2 +X3*2 )); + const double dksi_26 = 21 *sin(d2r( -X2 +X3*2 +X4 )); + const double dksi_27 = ( 17 - 0.1*jce)*sin(d2r( X1*2 )); + const double dksi_28 = 16 *sin(d2r( X0*2 -X2 +X4 )); + const double dksi_29 = ( -16 + 0.1*jce)*sin(d2r(-X0*2 +X1*2 +X3*2 +X4*2)); + const double dksi_30 = -15 *sin(d2r( X1 +X4 )); + const double dksi_31 = -13 *sin(d2r(-X0*2 +X2 +X4 )); + const double dksi_32 = -12 *sin(d2r( -X1 +X4 )); + const double dksi_33 = 11 *sin(d2r( X2*2 -X3*2 )); + const double dksi_34 = -10 *sin(d2r( X0*2 -X2 +X3*2 +X4 )); + const double dksi_35 = -8 *sin(d2r( X0*2 +X2 +X3*2 +X4*2)); + const double dksi_36 = 7 *sin(d2r( X1 +X3*2 +X4*2)); + const double dksi_37 = -7 *sin(d2r(-X0*2 +X1 +X2 )); + const double dksi_38 = -7 *sin(d2r( -X1 +X3*2 +X4*2)); + const double dksi_39 = -7 *sin(d2r( X0*2 +X3*2 +X4 )); + const double dksi_40 = 6 *sin(d2r( X0*2 +X2 )); + const double dksi_41 = 6 *sin(d2r(-X0*2 +X2*2 +X3*2 +X4*2)); + const double dksi_42 = 6 *sin(d2r(-X0*2 +X2 +X3*2 +X4 )); + const double dksi_43 = -6 *sin(d2r( X0*2 -X2*2 +X4 )); + const double dksi_44 = -6 *sin(d2r( X0*2 +X4 )); + const double dksi_45 = 5 *sin(d2r( -X1 +X2 )); + const double dksi_46 = -5 *sin(d2r(-X0*2 -X1 +X3*2 +X4 )); + const double dksi_47 = -5 *sin(d2r(-X0*2 +X4 )); + const double dksi_48 = -5 *sin(d2r( X2*2 +X3*2 +X4 )); + const double dksi_49 = 4 *sin(d2r(-X0*2 +X2*2 +X4 )); + const double dksi_50 = 4 *sin(d2r(-X0*2 +X1 +X3*2 +X4 )); + const double dksi_51 = 4 *sin(d2r( X2 -X3*2 )); + const double dksi_52 = -4 *sin(d2r(-X0 +X2 )); + const double dksi_53 = -4 *sin(d2r(-X0*2 +X1 )); + const double dksi_54 = -4 *sin(d2r( X0 )); + const double dksi_55 = 3 *sin(d2r( X2 +X3*2 )); + const double dksi_56 = -3 *sin(d2r( -X2*2 +X3*2 +X4*2)); + const double dksi_57 = -3 *sin(d2r(-X0 -X1 +X2 )); + const double dksi_58 = -3 *sin(d2r( X1 +X2 )); + const double dksi_59 = -3 *sin(d2r( -X1 +X2 +X3*2 +X4*2)); + const double dksi_60 = -3 *sin(d2r( X0*2 -X1 -X2 +X3*2 +X4*2)); + const double dksi_61 = -3 *sin(d2r( X2*3 +X3*2 +X4*2)); + const double dksi_62 = -3 *sin(d2r( X0*2 -X1 +X3*2 +X4*2)); + + const double deps_00 = (92025 + 8.9*jce)*cos(d2r( X4 )); + const double deps_01 = ( 5736 - 3.1*jce)*cos(d2r(-X0*2 +X3*2 +X4*2)); + const double deps_02 = ( 977 - 0.5*jce)*cos(d2r( X3*2 +X4*2)); + const double deps_03 = ( -895 + 0.5*jce)*cos(d2r( X4*2)); + const double deps_04 = ( 54 - 0.1*jce)*cos(d2r( X1 )); + const double deps_05 = -7 *cos(d2r( X2 )); + const double deps_06 = ( 224 - 0.6*jce)*cos(d2r(-X0*2 +X1 +X3*2 +X4*2)); + const double deps_07 = 200 *cos(d2r( X3*2 +X4 )); + const double deps_08 = ( 129 - 0.1*jce)*cos(d2r( X2 +X3*2 +X4*2)); + const double deps_09 = ( -95 + 0.3*jce)*cos(d2r(-X0*2 -X1 +X3*2 +X4*2)); + const double deps_10 = -70 *cos(d2r(-X0*2 +X3*2 +X4 )); + const double deps_11 = -53 *cos(d2r( -X2 +X3*2 +X4*2)); + const double deps_12 = -33 *cos(d2r( X2 +X4 )); + const double deps_13 = 26 *cos(d2r( X0*2 -X2 +X3*2 +X4*2)); + const double deps_14 = 32 *cos(d2r( -X2 +X4 )); + const double deps_15 = 27 *cos(d2r( X2 +X3*2 +X4 )); + const double deps_16 = -24 *cos(d2r( -X2*2 +X3*2 +X4 )); + const double deps_17 = 16 *cos(d2r( X0*2 +X3*2 +X4*2)); + const double deps_18 = 13 *cos(d2r( X2*2 +X3*2 +X4*2)); + const double deps_19 = -12 *cos(d2r(-X0*2 +X2 +X3*2 +X4*2)); + const double deps_20 = -10 *cos(d2r( -X2 +X3*2 +X4 )); + const double deps_21 = -8 *cos(d2r( X0*2 -X2 +X4 )); + const double deps_22 = 7 *cos(d2r(-X0*2 +X1*2 +X3*2 +X4*2)); + const double deps_23 = 9 *cos(d2r( X1 +X4 )); + const double deps_24 = 7 *cos(d2r(-X0*2 +X2 +X4 )); + const double deps_25 = 6 *cos(d2r( -X1 +X4 )); + const double deps_26 = 5 *cos(d2r( X0*2 -X2 +X3*2 +X4 )); + const double deps_27 = 3 *cos(d2r( X0*2 +X2 +X3*2 +X4*2)); + const double deps_28 = -3 *cos(d2r( X1 +X3*2 +X4*2)); + const double deps_29 = 3 *cos(d2r( -X1 +X3*2 +X4*2)); + const double deps_30 = 3 *cos(d2r( X0*2 +X3*2 +X4 )); + const double deps_31 = -3 *cos(d2r(-X0*2 +X2*2 +X3*2 +X4*2)); + const double deps_32 = -3 *cos(d2r(-X0*2 +X2 +X3*2 +X4 )); + const double deps_33 = 3 *cos(d2r( X0*2 -X2*2 +X4 )); + const double deps_34 = 3 *cos(d2r( X0*2 +X4 )); + const double deps_35 = 3 *cos(d2r(-X0*2 -X1 +X3*2 +X4 )); + const double deps_36 = 3 *cos(d2r(-X0*2 +X4 )); + const double deps_37 = 3 *cos(d2r( X2*2 +X3*2 +X4 )); + + const double dksi = (dksi_00 + dksi_01 + dksi_02 + dksi_03 + dksi_04 + dksi_05 + dksi_06 + dksi_07 + dksi_08 + dksi_09 + dksi_10 + dksi_11 + dksi_12 + dksi_13 + dksi_14 + dksi_15 + dksi_16 + dksi_17 + dksi_18 + dksi_19 + dksi_20 + dksi_21 + dksi_22 + dksi_23 + dksi_24 + dksi_25 + dksi_26 + dksi_27 + dksi_28 + dksi_29 + dksi_30 + dksi_31 + dksi_32 + dksi_33 + dksi_34 + dksi_35 + dksi_36 + dksi_37 + dksi_38 + dksi_39 @@ -506,7 +504,7 @@ auto nutation(const double julian_ephemeris_centry) { + dksi_50 + dksi_51 + dksi_52 + dksi_53 + dksi_54 + dksi_55 + dksi_56 + dksi_57 + dksi_58 + dksi_59 + dksi_60 + dksi_61 + dksi_62)/36000000.0; - const auto deps = (deps_00 + deps_01 + deps_02 + deps_03 + deps_04 + deps_05 + deps_06 + deps_07 + deps_08 + deps_09 + const double deps = (deps_00 + deps_01 + deps_02 + deps_03 + deps_04 + deps_05 + deps_06 + deps_07 + deps_08 + deps_09 + deps_10 + deps_11 + deps_12 + deps_13 + deps_14 + deps_15 + deps_16 + deps_17 + deps_18 + deps_19 + deps_20 + deps_21 + deps_22 + deps_23 + deps_24 + deps_25 + deps_26 + deps_27 + deps_28 + deps_29 + deps_30 + deps_31 + deps_32 + deps_33 + deps_34 + deps_35 + deps_36 + deps_37)/36000000.0; @@ -514,8 +512,7 @@ auto nutation(const double julian_ephemeris_centry) { return std::make_pair(dksi, deps); } -double true_ecliptic_obliquity(const double julian_ephemeris_millennium, - const double obliquity_nutation) { +inline double true_ecliptic_obliquity(const double julian_ephemeris_millennium, const double obliquity_nutation) { const double U = julian_ephemeris_millennium/10; const double U2 = U*U; const double U3 = U2*U; @@ -527,70 +524,81 @@ double true_ecliptic_obliquity(const double julian_ephemeris_millennium, const double U9 = U8*U; const double U10 = U9*U; // Equation (24) - const auto eps_0 = 84381.448 - 4680.93*U - 1.55*U2 + 1999.25*U3 + const double eps_0 = 84381.448 - 4680.93*U - 1.55*U2 + 1999.25*U3 - 51.38*U4 -249.67*U5 - 39.05*U6 + 7.12*U7 + 27.87*U8 + 5.79*U9 + 2.45*U10; // Equation (25) return eps_0/3600.0 + obliquity_nutation; } -double aberration_correction(const double R) { +inline double aberration_correction(const double R) { // Equation (26) return -20.4898/(3600.0*R); } -double apparent_sun_longitude(const double geoc_long, const double nut_long, const double ab_corr) { +inline double apparent_sun_longitude(const double geoc_long, const double nut_long, const double ab_corr) { // Equation (27) return geoc_long + nut_long + ab_corr; } -double apparent_Greenwich_sidereal_time(const double julian_day, - const double julian_ephemeris_centry, - const double nut_long, - const double true_ecliptic_obliquity) { + +inline double apparent_Greenwich_sidereal_time( + const double julian_day, + const double julian_ephemeris_centry, + const double nut_long, + const double true_ecliptic_obliquity +) { // Equation (28) - const double nu0 = 280.46061837 + 360.98564736629*(julian_day - 2451545) + 0.000387933*pow(julian_ephemeris_centry, 2) - - pow(julian_ephemeris_centry, 3)/38710000; + const double nu0 = + 280.46061837 + + 360.98564736629*(julian_day - 2451545) + + 0.000387933*pow(julian_ephemeris_centry, 2) + - pow(julian_ephemeris_centry, 3)/38710000; // Equation (29) with [0, 360) range limit return limit_degrees(nu0) + nut_long*cos(d2r(true_ecliptic_obliquity)); } -double geocentric_sun_right_ascension(const double app_sun_longitude, - const double true_ecliptic_obliq, - const double geocentric_lat) { +inline double geocentric_sun_right_ascension( + const double app_sun_longitude, const double true_ecliptic_obliq, const double geocentric_lat +) { // Equation (30) - const auto lambda_rad = d2r(app_sun_longitude); - const auto eps_rad = d2r(true_ecliptic_obliq); - const auto beta_rad = d2r(geocentric_lat); - const auto alpha = atan2(sin(lambda_rad)*cos(eps_rad) - tan(beta_rad)*sin(eps_rad), cos(lambda_rad)); + const double lambda_rad = d2r(app_sun_longitude); + const double eps_rad = d2r(true_ecliptic_obliq); + const double beta_rad = d2r(geocentric_lat); + const double alpha = atan2(sin(lambda_rad)*cos(eps_rad) - tan(beta_rad)*sin(eps_rad), cos(lambda_rad)); + return limit_degrees(r2d(alpha)); } -double geocentric_sun_declination(const double app_sun_longitude, - const double true_ecliptic_obliq, - const double geocentric_lat) { +inline double geocentric_sun_declination( + const double app_sun_longitude, const double true_ecliptic_obliq, const double geocentric_lat +) { // Equation (31) - const auto lambda_rad = d2r(app_sun_longitude); - const auto eps_rad = d2r(true_ecliptic_obliq); - const auto beta_rad = d2r(geocentric_lat); - const auto delta = asin(sin(beta_rad)*cos(eps_rad) + cos(beta_rad)*sin(eps_rad)*sin(lambda_rad)); + const double lambda_rad = d2r(app_sun_longitude); + const double eps_rad = d2r(true_ecliptic_obliq); + const double beta_rad = d2r(geocentric_lat); + const double delta = asin(sin(beta_rad)*cos(eps_rad) + cos(beta_rad)*sin(eps_rad)*sin(lambda_rad)); + return r2d(delta); } -double observer_local_hour_angle(const double app_Greenwich_sidereal_time, - const double geocentric_sun_right_asc, - const double geographical_lon) { +inline double observer_local_hour_angle( + const double app_Greenwich_sidereal_time, const double geocentric_sun_right_asc, const double geographical_lon +) { // Equation (32) const double H = app_Greenwich_sidereal_time + geographical_lon - geocentric_sun_right_asc; + return limit_degrees(H); } -auto topocentric_values(const double R, - const double geographical_latitude, - const double masl, - const double geocentric_sun_right_asc, - const double geocentric_sun_decl, - const double obs_local_hour_angle) { +inline std::tuple topocentric_values( + const double R, + const double geographical_latitude, + const double masl, + const double geocentric_sun_right_asc, + const double geocentric_sun_decl, + const double obs_local_hour_angle +) { // Convert to rad const double lat_rad = d2r(geographical_latitude); @@ -622,138 +630,133 @@ auto topocentric_values(const double R, return std::make_tuple(topo_sun_right_asc , topo_sun_declination, topo_local_hour_angle); } -auto uncorrected_topocentric_elevation_angle(const double geographical_latitude, - const double topo_sun_declination, - const double topo_local_hour_angle) { +inline double uncorrected_topocentric_elevation_angle( + const double geographical_latitude, const double topo_sun_declination, const double topo_local_hour_angle +) { // Equation (41) const double lat_rad = d2r(geographical_latitude); const double sun_decl_rad = d2r(topo_sun_declination); const double hour_angle_rad = d2r(topo_local_hour_angle); - const double eps_0 = r2d(asin(sin(lat_rad)*sin(sun_decl_rad) + - cos(lat_rad)*cos(sun_decl_rad)*cos(hour_angle_rad))); + const double eps_0 = r2d(asin(sin(lat_rad)*sin(sun_decl_rad) + cos(lat_rad)*cos(sun_decl_rad)*cos(hour_angle_rad))); + return eps_0; } -auto topocentric_zenith_angle(const double eps_0, - const double P, - const double T) { +inline std::pair topocentric_zenith_angle(const double eps_0, const double P, const double T) { // Equation (42) const double deps = P/1010*283/(273 + T)*1.02/(60.0*tan(d2r(eps_0 + 10.3/(eps_0 + 5.11)))); // Equation (43) const double topocentric_elevation_angle = eps_0 + deps; // Equation (44) const double topocentric_zenith_angle = 90 - topocentric_elevation_angle; + return std::make_pair(topocentric_elevation_angle, topocentric_zenith_angle); } -auto topocentric_astronomers_and_azimuth_angle(const double topo_local_hour_angle, - const double geographical_latitude, - const double topo_sun_declination) { +inline std::pair topocentric_astronomers_and_azimuth_angle( + const double topo_local_hour_angle, const double geographical_latitude, const double topo_sun_declination +) { // Equation (45) const double h_rad = d2r(topo_local_hour_angle); const double lat_rad = d2r(geographical_latitude); const double sun_decl_rad = d2r(topo_sun_declination); - const double Gamma = limit_degrees(r2d(atan2(sin(h_rad), - cos(h_rad)*sin(lat_rad) - tan(sun_decl_rad)*cos(lat_rad)))); + const double Gamma = limit_degrees( + r2d(atan2(sin(h_rad), cos(h_rad)*sin(lat_rad) - tan(sun_decl_rad)*cos(lat_rad))) + ); // Equation (46) const double Phi = limit_degrees(Gamma + 180); - return std::make_tuple(Gamma, Phi); + + return std::make_pair(Gamma, Phi); } -double indicence_angle(const double topo_zenith_angle, const double surface_slope, const double surface_azimuth, const double astronomers_azimuth_angle) { +inline double indicence_angle( + const double topo_zenith_angle, + const double surface_slope, + const double surface_azimuth, + const double astronomers_azimuth_angle +) { // Equation (47) const double zenit_rad = d2r(topo_zenith_angle); const double slope_rad = d2r(surface_slope); const double astro_azimuth_rad = d2r(astronomers_azimuth_angle); const double surf_azimuth_rad = d2r(surface_azimuth); - const double I = r2d(acos(cos(zenit_rad)*cos(slope_rad) - + sin(slope_rad)*sin(zenit_rad)*cos(astro_azimuth_rad - surf_azimuth_rad))); + const double I = r2d( + acos(cos(zenit_rad)*cos(slope_rad) + sin(slope_rad)*sin(zenit_rad)*cos(astro_azimuth_rad - surf_azimuth_rad)) + ); + return I; } template -auto solar_position(const double julian_day, - const double DT, - const double geographic_latitude, - const double geographic_longitude, - const double masl, - collector_t collector) { +auto solar_position( + const double julian_day, + const auto DT, + const double geographic_latitude, + const double geographic_longitude, + const double masl, + collector_t collector +) { // Note that the args should be in UT, and that |UT - UTC| < 1.0 //const auto DT = Delta_T(year); - const auto julian_ephemeris_day = jde(julian_day, DT); - const auto julian_ephemeris_centry = jce(julian_ephemeris_day); - const auto julian_ephemeris_millennium = jme(julian_ephemeris_centry); - const auto L = heliocentric_longitude(julian_ephemeris_millennium); - const auto B = heliocentric_latitude(julian_ephemeris_millennium); - const auto R = heliocentric_radius_vector(julian_ephemeris_millennium); - const auto beta = geocentric_latitude(B); - const auto Theta = geocentric_longitude(L); + const double julian_ephemeris_day = jde(julian_day, DT); + const double julian_ephemeris_centry = jce(julian_ephemeris_day); + const double julian_ephemeris_millennium = jme(julian_ephemeris_centry); + const double L = heliocentric_longitude(julian_ephemeris_millennium); + const double B = heliocentric_latitude(julian_ephemeris_millennium); + const double R = heliocentric_radius_vector(julian_ephemeris_millennium); + const double beta = geocentric_latitude(B); + const double Theta = geocentric_longitude(L); const auto [longitude_nutation, obliquity_nutation] = nutation(julian_ephemeris_centry); - const auto epsilon = true_ecliptic_obliquity(julian_ephemeris_millennium, obliquity_nutation); - const auto delta_tau = aberration_correction(R); - const auto lambda = apparent_sun_longitude(Theta, longitude_nutation, delta_tau); - const auto nu = apparent_Greenwich_sidereal_time(julian_day, - julian_ephemeris_centry, - longitude_nutation, - epsilon); - const auto alpha = geocentric_sun_right_ascension(lambda, epsilon, beta); - const auto delta = geocentric_sun_declination(lambda, epsilon, beta); - const auto H = observer_local_hour_angle(nu, alpha, geographic_longitude); - const auto [alpha_mark, delta_mark, H_mark] = topocentric_values(R, - geographic_latitude, - masl, - alpha, - delta, - H); - const double e0 = uncorrected_topocentric_elevation_angle(geographic_latitude, - delta_mark, - H_mark); - const auto [Gamma, Phi] = topocentric_astronomers_and_azimuth_angle(H_mark, - geographic_latitude, - delta_mark); + const double epsilon = true_ecliptic_obliquity(julian_ephemeris_millennium, obliquity_nutation); + const double delta_tau = aberration_correction(R); + const double lambda = apparent_sun_longitude(Theta, longitude_nutation, delta_tau); + const double nu = apparent_Greenwich_sidereal_time(julian_day, julian_ephemeris_centry, longitude_nutation, epsilon); + const double alpha = geocentric_sun_right_ascension(lambda, epsilon, beta); + const double delta = geocentric_sun_declination(lambda, epsilon, beta); + const double H = observer_local_hour_angle(nu, alpha, geographic_longitude); + const auto [alpha_mark, delta_mark, H_mark] = topocentric_values(R, geographic_latitude, masl, alpha, delta, H); + const double e0 = uncorrected_topocentric_elevation_angle(geographic_latitude, delta_mark, H_mark); + const auto [Gamma, Phi] = topocentric_astronomers_and_azimuth_angle(H_mark, geographic_latitude, delta_mark); + return collector(e0, Gamma, Phi, alpha_mark, delta_mark, H_mark); } -auto corrected_solar_elevation(const double e0, const double P, const double T) { +inline std::pair corrected_solar_elevation(const double e0, const double P, const double T) { return topocentric_zenith_angle(e0, P, T); } template -auto calendar_solar_position(unsigned int year, - unsigned int month, - double day, - const double geographic_latitude, - const double geographic_longitude, - const double masl, - collector_t collector, - dt_calc_t cal_calc) { - const double DT = cal_calc(year, month, day); +auto calendar_solar_position( + unsigned int year, + unsigned int month, + double day, + const double geographic_latitude, + const double geographic_longitude, + const double masl, + collector_t collector, + dt_calc_t cal_calc +) { + const auto DT = cal_calc(year, month, day); const double julian_day = jd_from_cal(year, month, day); - return solar_position(julian_day, - DT, - geographic_latitude, - geographic_longitude, - masl, collector); + return solar_position(julian_day, DT, geographic_latitude, geographic_longitude, masl, collector); } template -auto time_point_solar_position(const std::chrono::system_clock::time_point time_point, - const double geographic_latitude, - const double geographic_longitude, - const double masl, - collector_t collector, - dt_calc_t time_point_calc) { - +auto time_point_solar_position( + const std::chrono::system_clock::time_point time_point, + const double geographic_latitude, + const double geographic_longitude, + const double masl, + collector_t collector, + dt_calc_t time_point_calc +) { const auto DT = time_point_calc(time_point); - const auto julian_day = jd_from_clock(time_point); - return solar_position(julian_day, - DT, - geographic_latitude, - geographic_longitude, - masl, collector); + const double julian_day = jd_from_clock(time_point); + + return solar_position(julian_day, DT, geographic_latitude, geographic_longitude, masl, collector); } } diff --git a/src/rasputin/cpp/solar_position_bindings.cpp b/src/rasputin/cpp/solar_position_bindings.cpp new file mode 100644 index 0000000..4f7b555 --- /dev/null +++ b/src/rasputin/cpp/solar_position_bindings.cpp @@ -0,0 +1,41 @@ +#include +#include + +#include "solar_position.h" +#include "shade.h" + +namespace py = pybind11; + +using namespace std::chrono; + +void bind_solars(py::module& m) { + m.def("timestamp_solar_position", [] (const double timestamp, const double geographic_latitude, const double geographic_longitude, const double masl) { + const auto tp = sys_days{1970y /January / 1} + 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."); +} + +void bind_shades(py::module& m) { + 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("shade", [] (const rasputin::Mesh& mesh, const double timestamp) { + const auto secs = seconds(int(std::round(timestamp))); + const auto millisecs = milliseconds(int(round(1000*fmod(timestamp, 1)))); + const auto tp = sys_days{1970y / January / 1} + 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; + }); +} diff --git a/src/rasputin/cpp/triangulate_dem.h b/src/rasputin/cpp/triangulate_dem.h new file mode 100644 index 0000000..b9d4e2b --- /dev/null +++ b/src/rasputin/cpp/triangulate_dem.h @@ -0,0 +1,83 @@ +// +// Created by Ola Skavhaug on 08/10/2018. +// + +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include +#include + +#include "types.h" + +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}; } + }; +} diff --git a/src/rasputin/cpp/types.h b/src/rasputin/cpp/types.h new file mode 100644 index 0000000..e68c25f --- /dev/null +++ b/src/rasputin/cpp/types.h @@ -0,0 +1,63 @@ +#pragma once + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +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; +} + +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; +} diff --git a/src/rasputin/cpp/vec_bindings.cpp b/src/rasputin/cpp/vec_bindings.cpp new file mode 100644 index 0000000..4b5a3a1 --- /dev/null +++ b/src/rasputin/cpp/vec_bindings.cpp @@ -0,0 +1,84 @@ +#include +#include +#include +#include + +#include "types.h" + +namespace py = pybind11; + +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); + +template +py::buffer_info vecarray_buffer(std::vector> &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 */ + std::vector { v.size(), n }, /* Buffer dimensions */ + { sizeof(T) * n, sizeof(T) } /* Strides (in bytes) for each index */ + ); +} + +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 +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) }); +} + +void bind_vectors(py::module& 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"); +} diff --git a/src/rasputin/geometry.py b/src/rasputin/geometry.py index 21d53d2..131bd60 100644 --- a/src/rasputin/geometry.py +++ b/src/rasputin/geometry.py @@ -264,8 +264,8 @@ def to_cpp(self) -> Union[td.simple_polygon, td.polygon]: from shapely.geometry.polygon import orient # Get point sequences as arrays - exterior = np.asarray(orient(self.polygon).exterior)[:-1] - interiors = [np.asarray(hole)[:-1] + exterior = np.asarray(orient(self.polygon).exterior.coords)[:-1] + interiors = [np.asarray(hole.coords)[:-1] for hole in orient(self.polygon, -1).interiors] cgal_polygon = td.simple_polygon(exterior) diff --git a/src/rasputin/gml_repository.py b/src/rasputin/gml_repository.py index 9326ada..a38d73a 100644 --- a/src/rasputin/gml_repository.py +++ b/src/rasputin/gml_repository.py @@ -132,13 +132,15 @@ class GMLRepository(LandCoverRepository): land_cover_type = LandCoverType land_cover_meta_info_type = LandCoverMetaInfo - def __init__(self, path: Path) -> None: + + def __init__(self, path: Path, land_cover_code: str="clc18_kode") -> None: self.path = path files = list(path.glob("*.gml")) assert len(files) == 1 self.fn = files[0] self.parser = etree.XMLParser(encoding="utf-8", recover=True, huge_tree=True) self.data_crs: Optional[CRS] = None + self.land_cover_code = land_cover_code def _parse_polygon(self, polygon: etree._Element, nsmap: Dict[str, str]) -> Polygon: gml = f"{{{nsmap['gml']}}}" @@ -170,7 +172,15 @@ def read(self, domain: GeoPolygon) -> Dict[LandCoverType, List[Polygon]]: ogr = f"{{{root.nsmap['ogr']}}}" result = {} for elm in root.iter(f"{gml}featureMember"): - code = LandCoverType(int(next(elm.iter(f"{ogr}clc18_kode")).text)) + try: + code = LandCoverType(int(next(elm.iter(f"{ogr}{self.land_cover_code}")).text)) + except ValueError as err: + if len(err.args) == 1 and err.args[0].endswith("is not a valid LandCoverType"): + continue + else: + raise ValueError from err + except StopIteration as err: + raise StopIteration(f"LandCoverCode {self.land_cover_code} invalid") from err polygon = self._parse_polygon(next(elm.iter(f"{gml}Polygon")), root.nsmap) if polygon.intersects(domain.polygon): if code not in result: @@ -182,10 +192,11 @@ def constraints(self, *, domain: GeoPolygon) -> List[GeoPolygon]: return [] def land_cover(self, - *, - land_types: Optional[List[LandCoverType]], - geo_points: GeoPoints, - domain: GeoPolygon) -> np.ndarray: + *, + land_types: Optional[List[LandCoverType]], + geo_points: GeoPoints, + domain: GeoPolygon, + ) -> np.ndarray: tree = etree.parse(str(self.fn), self.parser) root = tree.getroot() self._assign_data_crs(root) diff --git a/src/rasputin/triangulate_dem.h b/src/rasputin/triangulate_dem.h deleted file mode 100644 index 3dc346d..0000000 --- a/src/rasputin/triangulate_dem.h +++ /dev/null @@ -1,916 +0,0 @@ -// -// 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 "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 -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; - } - - // 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); - - 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 { - // 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; - } - - // 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)); - - return rectangle; - } - - // Compute intersection of raster rectangle with polygon - template - CGAL::MultiPolygon compute_intersection(const P& polygon) const { - CGAL::SimplePolygon rectangle = exterior(); - - CGAL::MultiPolygon intersection_polygon; - CGAL::intersection(rectangle, polygon, std::back_inserter(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)); - } -}; - -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; - } - - // 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)); - } - } - 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; - - 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; - }; - - 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)); - } - } - 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(); - - 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())); - } - 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 -Mesh mesh_from_raster(const RasterData& raster, const std::string proj4_str) { - return make_mesh(raster.raster_points(), proj4_str); -} - -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); - } - ++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]}); -}; - -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]); - } - } - 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 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; -}; - -double compute_aspect(const point3 &normal) { - return std::atan2(normal[0], normal[1]); -} - -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; -}; - -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 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)); -} -} diff --git a/tests/test_mesh.py b/tests/test_mesh.py index 8fdec5b..caa501d 100644 --- a/tests/test_mesh.py +++ b/tests/test_mesh.py @@ -1,9 +1,9 @@ import pytest -from numpy import array, cos, sin, linspace, pi, float32, zeros, ndindex, sqrt, meshgrid +from numpy import array, cos, sin, linspace, pi, float32, zeros, ndindex, sqrt, meshgrid, int64 from numpy.linalg import norm from datetime import datetime, timedelta -import pyproj +import pyproj from shapely.geometry import Polygon, Point from rasputin.geometry import GeoPolygon @@ -152,7 +152,7 @@ def dist(mesh, pt): def test_mesh_sub_mesh(raster, polygon): mesh = Mesh.from_raster(data=raster, domain=polygon) - sub_mesh = mesh.extract_sub_mesh(array([0, 1])) + sub_mesh = mesh.extract_sub_mesh(array([0, 1], dtype=int64)) assert len(sub_mesh.faces) == 2 def test_mesh_w_hole(raster, polygon_w_hole): diff --git a/tests/test_polygons.py b/tests/test_polygons.py index edbeaf5..38bb60d 100644 --- a/tests/test_polygons.py +++ b/tests/test_polygons.py @@ -1,4 +1,5 @@ from numpy import array, cos, sin, linspace, pi + import pytest import rasputin.triangulate_dem as td @@ -9,7 +10,12 @@ def rectangle_polygon(): [1.0, 0.0], [1.0, 1.0], [0.48, 1.0]]) - +@pytest.fixture +def square_polygon(): + return td.simple_polygon([[0.48, 0.0], + [0.96, 0.0], + [0.96, 0.96], + [0.48, 0.96]]) @pytest.fixture def circle_polygon(): @@ -20,8 +26,35 @@ def circle_polygon(): return td.simple_polygon([(x + r*cos(t), y + r*sin(t)) for t in linspace(0, 2*pi, N+1)[:-1]]) - @pytest.mark.skip # Segfaults -def test_intersection(rectangle_polygon, circle_polygon): +def test_intersection_rect_circ(rectangle_polygon, circle_polygon): polygon = rectangle_polygon.intersection(circle_polygon) assert polygon.num_parts() == 1 + +@pytest.mark.skip # Segfaults +def test_intersection_square_circ(square_polygon, circle_polygon): + polygon = square_polygon.intersection(circle_polygon) + assert len(polygon) == 1 + +def test_intersection_rect_square(rectangle_polygon, square_polygon): + polygon = rectangle_polygon.intersection(square_polygon) + assert isinstance(polygon, td.multi_polygon) and polygon.num_parts() == 1 + assert isinstance(polygon[0], td.polygon) + assert len(polygon[0].holes()) == 0 + assert isinstance(polygon[0].exterior(), td.simple_polygon) + + polygon = square_polygon.intersection(rectangle_polygon) + assert isinstance(polygon, td.multi_polygon) and polygon.num_parts() == 1 + assert isinstance(polygon[0], td.polygon) + assert len(polygon[0].holes()) == 0 + assert isinstance(polygon[0].exterior(), td.simple_polygon) + +def test_difference_rect_square(rectangle_polygon, square_polygon): + polygon = rectangle_polygon.difference(square_polygon) + assert isinstance(polygon, td.multi_polygon) and polygon.num_parts() == 1 + assert isinstance(polygon[0], td.polygon) + assert len(polygon[0].holes()) == 0 + assert isinstance(polygon[0].exterior(), td.simple_polygon) + + polygon = square_polygon.difference(rectangle_polygon) + assert isinstance(polygon, td.multi_polygon) and polygon.num_parts() == 0 From aa224630bf25e97e2b5057b4b96fbf68fd8b0e44 Mon Sep 17 00:00:00 2001 From: Alocias Date: Thu, 20 Oct 2022 16:06:11 +0200 Subject: [PATCH 2/2] Add general RasterBase type with derived implementations using CGAL, add General MeshBase and implement with CGAL, add trait Triangulation, add MeshInit trait for initialization specific functions used in MeshBase, add Shade trait, add trait PolygonUtils and use in Triangulation trait special for CGAL, make binding code general using traits with no need for specialization --- src/rasputin/cpp/bindings.cpp | 8 +- src/rasputin/cpp/cgal_mesh.h | 296 +++++++++++++ src/rasputin/cpp/cgal_polygon.h | 149 +++++++ src/rasputin/cpp/cgal_raster_data.h | 65 +++ src/rasputin/cpp/cgal_shade.h | 109 +++++ src/rasputin/cpp/mesh.h | 442 ++++++++++--------- src/rasputin/cpp/mesh_bindings.cpp | 71 +-- src/rasputin/cpp/polygon.h | 228 ++++++++-- src/rasputin/cpp/polygon_bindings.cpp | 152 +------ src/rasputin/cpp/raster_bindings.cpp | 59 +-- src/rasputin/cpp/raster_data.h | 184 ++++++-- src/rasputin/cpp/shade.h | 254 +++++------ src/rasputin/cpp/shade_bindings.cpp | 7 + src/rasputin/cpp/solar_position_bindings.cpp | 19 - 14 files changed, 1321 insertions(+), 722 deletions(-) create mode 100644 src/rasputin/cpp/cgal_mesh.h create mode 100644 src/rasputin/cpp/cgal_polygon.h create mode 100644 src/rasputin/cpp/cgal_raster_data.h create mode 100644 src/rasputin/cpp/cgal_shade.h create mode 100644 src/rasputin/cpp/shade_bindings.cpp diff --git a/src/rasputin/cpp/bindings.cpp b/src/rasputin/cpp/bindings.cpp index 4aa29ad..46b7093 100644 --- a/src/rasputin/cpp/bindings.cpp +++ b/src/rasputin/cpp/bindings.cpp @@ -3,17 +3,17 @@ namespace py = pybind11; void bind_vectors(py::module&); -void bind_raster(py::module&); -void bind_mesh(py::module&); +void bind_raster(py::module& m); +void bind_meshes(py::module& m); void bind_polygons(py::module&); -void bind_solars(py::module&); void bind_shades(py::module&); +void bind_solars(py::module&); void bind_slopes(py::module&); PYBIND11_MODULE(triangulate_dem, m) { bind_vectors(m); bind_raster(m); - bind_mesh(m); + bind_meshes(m); bind_polygons(m); bind_slopes(m); bind_shades(m); diff --git a/src/rasputin/cpp/cgal_mesh.h b/src/rasputin/cpp/cgal_mesh.h new file mode 100644 index 0000000..7265322 --- /dev/null +++ b/src/rasputin/cpp/cgal_mesh.h @@ -0,0 +1,296 @@ +#pragma once + +#include +#include +#include +#include +#include + +#include "cgal_polygon.h" +#include "cgal_raster_data.h" +#include "mesh.h" + +namespace SMS = CGAL::Surface_mesh_simplification; + +namespace rasputin { +using VertexIndex = CGAL::VertexIndex; +using face_descriptor = CGAL::face_descriptor; +using Constraints = CGAL::DelaunayConstraints; + +using Base = traits::MeshBase< + CGAL::Mesh, VertexIndex, face_descriptor, Constraints, PointList, Polygon, SimplePolygon, MultiPolygon +>; + +struct Mesh: Base { + Mesh(CGAL::Mesh mesh, const std::string proj4_str): Base(mesh, proj4_str) { } + + static constexpr char coarsen_ratio_name[] = "lindstrom_turk_by_ratio"; + static constexpr char coarsen_size_name[] = "lindstrom_turk_by_size"; + static constexpr char coarsen_ratio_description[] = + "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."; + static constexpr char coarsen_size_description[] = + "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."; + + Mesh copy() const { + // Call CGAL::Mesh copy constructor to do a deep copy + return Mesh(CGAL::Mesh(this->mesh), this->proj4_str); + } + + template + requires std::same_as> + Mesh _coarsen( + const Stop& stop, + const SMS::LindstromTurk_placement& placement, + const SMS::LindstromTurk_cost& cost + ) const { + CGAL::Mesh new_cgal_mesh = CGAL::Mesh(this->mesh); + CGAL::Surface_mesh_simplification::edge_collapse( + new_cgal_mesh, stop, CGAL::parameters::get_cost(cost).get_placement(placement) + ); + return Mesh(new_cgal_mesh, this->proj4_str); + } + + template + Mesh coarsen(PredRatio ratio) const { + return this->_coarsen>( + SMS::Count_ratio_stop_predicate(ratio), + SMS::LindstromTurk_placement(), + SMS::LindstromTurk_cost() + ); + } + + Mesh extract_sub_mesh(const std::vector &face_indices) const { + return Mesh(this->construct_sub_mesh(face_indices), this->proj4_str); + } + + size_t num_edges() const { + return this->mesh.number_of_edges(); + } + size_t num_vertices() const { + return this->mesh.number_of_vertices(); + } + size_t num_faces() const { + return this->mesh.number_of_faces(); + } +}; + +namespace traits { +template<> +struct MeshInit { + using Mesh_t = CGAL::Mesh; + using VertexIndex_t = VertexIndex; + using face_descriptor_t = face_descriptor; + using VertexIndexMap = std::map; + using FaceDescrMap = std::map; + + static void set_points_faces(Mesh_t& mesh, point3_vector& points, face_vector& faces) { + points.clear(); + faces.clear(); + + points.reserve(mesh.num_vertices()); + faces.reserve(mesh.num_faces()); + + int n = 0; + std::map reindex; + for (auto f: mesh.faces()) { + std::array fl; + size_t idx = 0; + for (auto v: mesh.vertices_around_face(mesh.halfedge(f))) { + if (reindex.count(v) == 0) { + reindex.emplace(v, n++); + const auto pt = 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]}); + } + } + + static Mesh_t construct_mesh( + const point3_vector &pts, const face_vector &faces, VertexIndexMap& index_map, FaceDescrMap& face_map + ) { + Mesh_t 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 +struct Triangulation { + using PolygonUtils_t = PolygonUtils; + using T = traits::Triangulation; + using RasterData_t = RasterData; + using Mesh_t = Mesh; + using VertexIndex_t = typename Mesh_t::VertexIndex_t; + using face_descriptor_t = typename Mesh_t::face_descriptor_t; + using Constraints_t = typename Mesh_t::Constraints_t; + + using VertexIndexMap = std::map; + using FaceDescrMap = std::map; + + template P> + static Constraints 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 : PolygonUtils_t::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; + } + + // 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)); + } + } + return interpolated_points; + } + + template Pgn> + static Mesh_t make_mesh( + const typename RasterData::PointList_t &pts, + const Pgn& inclusion_polygon, + const Constraints &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; + + 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; + }; + + 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 (PolygonUtils_t::point_inside_polygon(face_midpoint, inclusion_polygon)) { + mesh.add_face(index(u), index(v), index(w)); + } + } + pvm.clear(); + return Mesh_t(mesh, proj4_str); + } + + static Mesh_t make_mesh(const typename RasterData::PointList_t &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(); + + return Mesh_t(mesh, proj4_str); + } + + template Pgn> + static Mesh_t 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); + } + + static Mesh_t 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())); + } + return make_mesh(raster_points, proj4_str); + } + + template Pgn> + static Mesh_t 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); + } + + static Mesh_t mesh_from_raster(const RasterData& raster, const std::string proj4_str) { + return make_mesh(raster.raster_points(), proj4_str); + } +}; +} // namespace traits +} // namespace rasputin diff --git a/src/rasputin/cpp/cgal_polygon.h b/src/rasputin/cpp/cgal_polygon.h new file mode 100644 index 0000000..90b406b --- /dev/null +++ b/src/rasputin/cpp/cgal_polygon.h @@ -0,0 +1,149 @@ +#pragma once + +#include +#include +#include + +#include +#include + +#include "types.h" +#include "polygon.h" + +namespace rasputin::traits { +using Point2 = CGAL::Point2; +using Polygon = CGAL::Polygon; +using SimplePolygon = CGAL::SimplePolygon; + +template<> +struct PolygonUtils { + using Polygon_t = Polygon; + using SimplePolygon_t = SimplePolygon; + using MultiPolygon_t = std::vector; + + static bool point_inside_polygon(const Point2 &x, const SimplePolygon &polygon) { + return polygon.has_on_bounded_side(x); + } + + static 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; + } + + static 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; + } + + static std::vector extract_boundaries(const MultiPolygon_t &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; + } + + static SimplePolygon exterior(const Polygon& polygon) { + return polygon.outer_boundary(); + } + + static py::list holes(const Polygon& polygon) { + py::list result; + for (auto h = polygon.holes_begin(); h != polygon.holes_end(); ++h) + result.append(*h); + return result; + } + + static py::array_t polygon_array(const SimplePolygon& polygon) { + py::array_t result(polygon.size() * 2); + result.resize(py::array::ShapeContainer({static_cast(polygon.size()), 2}), true); + auto info = result.request(); + double* data = static_cast(info.ptr); + for (auto v = polygon.vertices_begin(); v != polygon.vertices_end(); ++v) { + data[0] = v->x(); + data[1] = v->y(); + data += 2; + } + return result; + } + + template + static MultiPolygon_t difference_polygons(const P0& polygon0, const P1& polygon1) { + MultiPolygon_t difference_result; + CGAL::difference(polygon0, polygon1, std::back_inserter(difference_result)); + + return difference_result; + } + + template + static MultiPolygon_t intersect_polygons(const P0& polygon0, const P1& polygon1) { + MultiPolygon_t intersection_result; + CGAL::intersection(polygon0, polygon1, std::back_inserter(intersection_result)); + + return intersection_result; + } + + template + static MultiPolygon_t join_polygons(const P0& polygon0, const P1& polygon1) { + Polygon joined; + MultiPolygon_t join_result; + if (CGAL::join(polygon0, polygon1, joined)) + join_result.push_back(std::move(joined)); + else { + join_result.push_back(Polygon(polygon0)); + join_result.push_back(Polygon(polygon1)); + } + return join_result; + } + + static MultiPolygon_t join_multipolygons(const MultiPolygon_t& polygon0, const MultiPolygon_t& polygon1) { + MultiPolygon_t join_result; + CGAL::join(polygon0.begin(), polygon0.end(), polygon1.begin(), polygon1.end(), std::back_inserter(join_result)); + + return join_result; + } + + static 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; + + 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; + } + + static std::size_t num_vertices(const SimplePolygon& polygon) { + return polygon.size(); + } +}; +} // namespace rasputin::traits diff --git a/src/rasputin/cpp/cgal_raster_data.h b/src/rasputin/cpp/cgal_raster_data.h new file mode 100644 index 0000000..c93fd89 --- /dev/null +++ b/src/rasputin/cpp/cgal_raster_data.h @@ -0,0 +1,65 @@ +#pragma once + +#include + +#include "raster_data.h" + +namespace rasputin { +using PointList = CGAL::PointList; +using Polygon = CGAL::Polygon; +using SimplePolygon = CGAL::SimplePolygon; +using MultiPolygon = CGAL::MultiPolygon; + +template +requires std::floating_point +struct RasterData: public traits::RasterBase> { + 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 + ): + traits::RasterBase>( + x_min, y_max, delta_x, delta_y, num_points_x, num_points_y, data + ) + { + } + + PointList raster_points() const { + PointList points; + points.reserve(this->num_points_x * this->num_points_y); + + for (std::size_t i = 0; i < this->num_points_y; ++i) + for (std::size_t j = 0; j < this->num_points_x; ++j) + points.emplace_back( + this->x_min + j*this->delta_x, this->y_max - i*this->delta_y, this->data[i*this->num_points_x + j] + ); + return points; + } + + // Boundary of the raster domain as a CGAL polygon + SimplePolygon exterior() const { + SimplePolygon rectangle; + rectangle.push_back(CGAL::Point2(this->x_min, this->get_y_min())); + rectangle.push_back(CGAL::Point2(this->get_x_max(), this->get_y_min())); + rectangle.push_back(CGAL::Point2(this->get_x_max(), this->y_max)); + rectangle.push_back(CGAL::Point2(this->x_min, this->y_max)); + + return rectangle; + } + + // Compute intersection of raster rectangle with polygon + template P> + MultiPolygon compute_intersection(const P& polygon) const { + SimplePolygon rectangle = exterior(); + + MultiPolygon intersection_polygon; + intersection(rectangle, polygon, std::back_inserter(intersection_polygon)); + + return intersection_polygon; + } +}; +} diff --git a/src/rasputin/cpp/cgal_shade.h b/src/rasputin/cpp/cgal_shade.h new file mode 100644 index 0000000..05879a0 --- /dev/null +++ b/src/rasputin/cpp/cgal_shade.h @@ -0,0 +1,109 @@ +#pragma once + +#include +#include +#include +#include + + +#include +#include + +#include "shade.h" +#include "cgal_mesh.h" + +namespace rasputin::traits { +template<> +struct Shade { + using Mesh_t = Mesh; + using S = traits::Shade; + using face_descriptor_t = typename Mesh_t::face_descriptor_t; + + static CGAL::Point3 centroid(const Mesh_t& mesh, const face_descriptor_t &face) { + CGAL::Point3 c{0, 0, 0}; + for (auto v: mesh.mesh.vertices_around_face(mesh.mesh.halfedge(face))) { + const auto pt = mesh.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}; + } + + static bool is_shaded( + const CGAL::Tree& tree, + const face_descriptor_t& 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; + } + + static std::vector compute_shadow(const Mesh_t& mesh, const point3& sun_direction) { + std::vector shade; + auto cgal_mesh = mesh.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 CGAL::Point3 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); + } + ++i; + } + return shade; + }; + + static std::vector shade(const Mesh_t& 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.mesh).first, CGAL::faces(mesh.mesh).second, mesh.mesh); + for (auto fd: mesh.mesh.faces()) { + const CGAL::Point3 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.mesh); + shade_vec.emplace_back(is_shaded(tree, fd, face_normal, c, azimuth, elevation)); + } + return shade_vec; + } +}; +} // namespace rasputin::traits diff --git a/src/rasputin/cpp/mesh.h b/src/rasputin/cpp/mesh.h index 3fe348a..8f2f33d 100644 --- a/src/rasputin/cpp/mesh.h +++ b/src/rasputin/cpp/mesh.h @@ -1,60 +1,70 @@ #pragma once -#include + +#include + +#include +#include +#include #include "types.h" #include "raster_data.h" -#include "polygon.h" + +namespace py = pybind11; namespace rasputin { -inline 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; +namespace traits { +template +struct MeshInit { + static_assert("Trait MeshInit not implemented!"); }; -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); - } +template +struct Triangulation { + static_assert("Trait Trianguation not implemented!"); +}; - Mesh copy() const { - // Call CGAL::Mesh copy constructor to do a deep copy - return Mesh(CGAL::Mesh(cgal_mesh), proj4_str); - } +template +concept PredRatioType = + std::same_as || (std::unsigned_integral && !std::same_as); + +template< + typename Mesh, + typename VertexIndex, + typename face_descriptor, + typename Constraints, + typename PointList, + typename Polygon, + typename SimplePolygon, + typename MultiPolygon +> +class MeshBase { + public: + using Mesh_t = Mesh; + using VertexIndex_t = VertexIndex; + using face_descriptor_t = face_descriptor; + using PointList_t = PointList; + using Polygon_t = Polygon; + using SimplePolygon_t = SimplePolygon; + using MultiPolygon_t = MultiPolygon; + using Constraints_t = Constraints; + + const Mesh mesh; + const std::string proj4_str; + point3_vector points; + face_vector faces; - const point3_vector& get_points() const { - return points; + MeshBase(Mesh mesh, const std::string proj4_str): mesh(mesh), proj4_str(proj4_str) { + MeshInit::set_points_faces(mesh, points, faces); } - const face_vector& get_faces() const { - return faces; - } + template + Mesh coarsen(PredRatio ratio) const; - 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();} + size_t num_edges() const; + size_t num_vertices() const; + size_t num_faces() const; - Mesh extract_sub_mesh(const std::vector &face_indices) const { + Mesh_t construct_sub_mesh(const std::vector& face_indices) const { std::map remap; point3_vector new_points; face_vector new_faces; @@ -71,199 +81,217 @@ struct Mesh { } 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); + + return MeshInit::construct_mesh( + new_points, new_faces, index_map, face_map + ); } - private: - point3_vector points; - face_vector faces; + const point3_vector& get_points() const { + return points; + } - 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 face_vector& get_faces() const { + return faces; } }; -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; - } - - // 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)); - } +template +concept MeshInitType = + requires { + typename Mesh::Mesh_t; + typename Mesh::VertexIndex_t; + typename Mesh::face_descriptor_t; + typename Mesh::VertexIndexMap; + typename Mesh::FaceDescrMap; } - 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; - - 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; + && requires { + { Mesh::VertexIndexMap } -> std::same_as>; + { Mesh::FaceDescrMap } -> std::same_as>; }; - 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(); +template +concept MeshType = + requires(Mesh m) { + { m.copy() } -> std::same_as; + { m.coarsen_ratio_name[0] } -> std::same_as; + { m.coarsen_ratio_description[0] } -> std::same_as; + { m.coarsen_size_name[0] } -> std::same_as; + { m.coarsen_size_description[0] } -> std::same_as; + } + && requires { + typename Mesh::Mesh_t; + typename Mesh::VertexIndex_t; + typename Mesh::face_descriptor_t; + typename Mesh::Constraints_t; + typename Mesh::PointList_t; + typename Mesh::Polygon_t; + typename Mesh::SimplePolygon_t; + typename Mesh::MultiPolygon_t; + } + && requires(Mesh m, const std::vector& face_indices) { + { m.extract_sub_mesh(face_indices) } -> std::same_as; + } + && std::derived_from< + Mesh, + MeshBase< + typename Mesh::Mesh_t, + typename Mesh::VertexIndex_t, + typename Mesh::face_descriptor_t, + typename Mesh::Constraints_t, + typename Mesh::PointList_t, + typename Mesh::Polygon_t, + typename Mesh::SimplePolygon_t, + typename Mesh::MultiPolygon_t + > + >; + +} // namespace traits + +template +struct MeshInit { + using M = traits::MeshInit; + using Mesh_t = Mesh; + using VertexIndex_t = VertexIndex; + using face_descriptor_t = face_descriptor; + using VertexIndexMap = std::map; + using FaceDescrMap = std::map; + + static void set_points_faces(Mesh_t const& mesh, point3_vector& points, face_vector& faces) { + M::set_points_faces(mesh); + } - // 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)); - } + static Mesh_t construct_mesh ( + const point3_vector &pts, const face_vector &faces, VertexIndexMap& index_map, FaceDescrMap& face_map + ) { + return M::construct_mesh(pts, faces, index_map, face_map); } - pvm.clear(); - return Mesh(mesh, proj4_str); }; +template +struct Triangulation { + using RasterData_t = RasterData; + using Mesh_t = Mesh; + using VertexIndex_t = typename Mesh::VertexIndex_t; + using face_descriptor_t = typename Mesh::face_descriptor_t; + using Constraints_t = typename Mesh::Constraints_t; + using T = traits::Triangulation; + using VertexIndexMap = std::map; + using FaceDescrMap = std::map; + + template P> + static Constraints_t interpolate_boundary_points(const RasterData_t& raster, const P& boundary_polygon) { + return T::template interpolate_boundary_points(raster, boundary_polygon); + } -inline 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; + template Pgn> + static Mesh make_mesh( + const typename RasterData_t::PointList_t &pts, + const Pgn& inclusion_polygon, + const Constraints_t &constraints, + const std::string proj4_str + ) { + return T::template make_mesh(pts, inclusion_polygon, constraints, proj4_str); + } - 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()))); + static Mesh make_mesh(const typename RasterData_t::PointList_t &pts, const std::string proj4_str) { + return T::template make_mesh(pts, proj4_str); } - 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]); + template Pgn> + static Mesh mesh_from_raster( + const std::vector& raster_list, const Pgn& boundary_polygon, const std::string proj4_str + ) { + return T::template mesh_from_raster(raster_list, boundary_polygon, proj4_str); } - pvm.clear(); - return Mesh(mesh, proj4_str); -}; + static Mesh mesh_from_raster(const std::vector& raster_list, const std::string proj4_str) { + return T::mesh_from_raster(raster_list, 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())); + template Pgn> + static Mesh mesh_from_raster(const RasterData_t& raster, const Pgn& boundary_polygon, const std::string proj4_str) { + return T::template mesh_from_raster(raster, proj4_str); } - return make_mesh(raster_points, boundary_polygon, boundary_points, proj4_str); -} + static Mesh mesh_from_raster(const RasterData_t& raster, const std::string proj4_str) { + return T::mesh_from_raster(raster, 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())); +concept TriangulationType = + requires { + typename T::RasterData_t; + typename T::Mesh_t; + typename T::VertexIndex_t; + typename T::face_descriptor_t; + typename T::Constraints_t; } - return make_mesh(raster_points, proj4_str); + && traits::MeshType; + +template P> +void bind_make_mesh(py::module &m) { + using R = typename T::RasterData_t; + using Mesh_t = typename T::Mesh_t::Mesh_t; + using M = MeshInit; + + m.def("make_mesh", + [] (const R& raster_data, const P polygon, const std::string proj4_str) { + return T::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 T::mesh_from_raster(raster_data, proj4_str); + }, py::return_value_policy::take_ownership) + .def("make_mesh", + [] (const std::vector& raster_data, const P polygon, const std::string proj4_str) { + return T::mesh_from_raster(raster_data, polygon, proj4_str); + }, py::return_value_policy::take_ownership) + .def("make_mesh", + [] (const std::vector& raster_data, const std::string proj4_str) { + return T::mesh_from_raster(raster_data, proj4_str); + }, py::return_value_policy::take_ownership) + .def("construct_mesh", + [] (const rasputin::point3_vector& points, const rasputin::face_vector & faces, const std::string proj4_str) { + typename T::VertexIndexMap index_map; + typename T::FaceDescrMap face_map; + return typename T::Mesh_t(M::construct_mesh(points, faces, index_map, face_map), proj4_str); + }, py::return_value_policy::take_ownership); } - -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 -Mesh mesh_from_raster(const RasterData& raster, const std::string proj4_str) { - return make_mesh(raster.raster_points(), proj4_str); +template +void bind_make_mesh(py::module& m) { + using R = typename T::RasterData_t; + rasputin::bind_make_mesh(m); + rasputin::bind_make_mesh(m); } +template +void bind_mesh(py::module& m) { + py::class_>(m, "Mesh") + .def(Mesh::coarsen_ratio_name, + [] (const Mesh& self, double ratio) { return self.template coarsen(ratio); }, + py::return_value_policy::take_ownership, + Mesh::coarsen_ratio_description + ) + .def(Mesh::coarsen_size_name, + [] (const Mesh& self, std::size_t max_size) { return self.template coarsen(max_size); }, + py::return_value_policy::take_ownership, + Mesh::coarsen_size_description + ) + .def("copy", &Mesh::copy, py::return_value_policy::take_ownership) + .def("extract_sub_mesh", &Mesh::extract_sub_mesh, py::return_value_policy::take_ownership) + .def_property_readonly("num_vertices", &Mesh::num_vertices) + .def_property_readonly("num_edges", &Mesh::num_edges) + .def_property_readonly("num_faces", &Mesh::num_faces) + .def_property_readonly("points", &Mesh::get_points, py::return_value_policy::reference_internal) + .def_property_readonly("faces", &Mesh::get_faces, py::return_value_policy::reference_internal); + + ( + bind_make_mesh>(m) ,...); } +} // namespace rasputin diff --git a/src/rasputin/cpp/mesh_bindings.cpp b/src/rasputin/cpp/mesh_bindings.cpp index 55252c0..80d2a6e 100644 --- a/src/rasputin/cpp/mesh_bindings.cpp +++ b/src/rasputin/cpp/mesh_bindings.cpp @@ -1,81 +1,16 @@ -#include - #include -#include #include -// Surface mesh simplication policies -#include -#include -#include -#include - -#include "raster_data.h" -#include "mesh.h" +#include "cgal_mesh.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(CGAL::MultiPolygon); PYBIND11_MAKE_OPAQUE(std::vector>); PYBIND11_MAKE_OPAQUE(std::vector>); -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); -} - -void bind_mesh(py::module& m) { - 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); - - 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); - - m.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); +void bind_meshes(py::module& m) { + rasputin::bind_mesh, rasputin::RasterData>(m); } diff --git a/src/rasputin/cpp/polygon.h b/src/rasputin/cpp/polygon.h index 3ca52f8..2e8a0ac 100644 --- a/src/rasputin/cpp/polygon.h +++ b/src/rasputin/cpp/polygon.h @@ -1,54 +1,206 @@ #pragma once +#include +#include -#include "types.h" +#include +#include -namespace CGAL { -inline bool point_inside_polygon(const Point2 &x, const SimplePolygon &polygon) { - return polygon.has_on_bounded_side(x); -} +namespace py = pybind11; + +namespace rasputin { +namespace traits { +template +struct PolygonUtils { + static_assert("PolygonUtils trait not implemented"); +}; + +template +concept PolygonType = std::same_as || std::same_as; +} // namespace traits + +template +struct PolygonUtils { + using P = traits::PolygonUtils; + using Polygon_t = Polygon; + using SimplePolygon_t = SimplePolygon; + using MultiPolygon_t = std::vector; + + template + static bool point_inside_polygon(const Point2 &x, const SimplePolygon &polygon) { + return P::point_inside_polygon(x, polygon); + } + + template + static bool point_inside_polygon(const Point2 &x, const Polygon &polygon) { + return P::point_inside_polygon(x, polygon); + } -inline bool point_inside_polygon(const Point2 &x, const Polygon &polygon) { - if (not polygon.outer_boundary().has_on_bounded_side(x)) + template + static bool point_inside_polygon(const Point2 &x, const MultiPolygon_t &polygon) { + for (const auto& part: polygon) + if (P::point_inside_polygon(x, part)) + return true; 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; -} + } -inline 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; -} + static std::vector extract_boundaries(const SimplePolygon &polygon) { + std::vector ret; + ret.emplace_back(polygon); + + return ret; + } + + static py::array_t polygon_array(const SimplePolygon& polygon) { + return P::polygon_array(polygon); + } + + static std::vector extract_boundaries(const Polygon &polygon) { + return P::extract_boundaries(polygon); + } + + static std::vector extract_boundaries(const MultiPolygon_t &polygon) { + return P::extract_boundaries(polygon); + } + + static SimplePolygon polygon_from_numpy(py::array_t& buf) { + return P::polygon_from_numpy(buf); + } + + static MultiPolygon_t join_multipolygons(const MultiPolygon_t& polygon0, const MultiPolygon_t& polygon1) { + return P::join_multipolygons(polygon0, polygon1); + } + + template + requires traits::PolygonType && traits::PolygonType + static MultiPolygon_t join_polygons(const P0& polygon0, const P1& polygon1) { + return P::join_polygons(polygon0, polygon1); + } + + template + requires traits::PolygonType && traits::PolygonType + static MultiPolygon_t intersect_polygons(const P0& polygon0, const P1& polygon1) { + return P::intersect_polygons(polygon0, polygon1); + } + + template + requires traits::PolygonType && traits::PolygonType + static MultiPolygon_t difference_polygons(const P0& polygon0, const P1& polygon1) { + return P::difference_polygons(polygon0, polygon1); + } + + static SimplePolygon exterior(const Polygon& polygon) { + return P::exterior(polygon); + } + + static py::list holes(const Polygon& polygon) { + return P::holes(polygon); + } + + static std::size_t num_vertices(const SimplePolygon& polygon) { + return P::num_vertices(polygon); + } +}; -inline std::vector extract_boundaries(const SimplePolygon &polygon) { - std::vector ret; - ret.emplace_back(polygon); +template +concept PolygonUtilsType = + requires { + typename P::Polygon_t; + typename P::SimplePolygon_t; + typename P::MultiPolygon_t; + } + && std::same_as>; + +template +void bind_polygon_class(py::module& m) { + using Polygon = typename P::Polygon_t; + using SimplePolygon = typename P::SimplePolygon_t; - return ret; + py::class_>(m, "polygon") + .def(py::init([] (py::array_t& buf) { + const SimplePolygon exterior = P::polygon_from_numpy(buf); + return Polygon(exterior);})) + .def("holes", &P::holes) + .def("extract_boundaries", [] (const Polygon& self) { return P::extract_boundaries(self); }, + py::return_value_policy::take_ownership) + .def("exterior", &P::exterior, py::return_value_policy::take_ownership) + .def("join", &P::template join_polygons) + .def("join", &P::template join_polygons) + .def("difference", &P::template difference_polygons) + .def("difference", &P::template difference_polygons) + .def("intersection", &P::template intersect_polygons) + .def("intersection", &P::template intersect_polygons); } -inline 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); +template +void bind_simple_polygon_class(py::module& m) { + using Polygon = typename P::Polygon_t; + using SimplePolygon = typename P::SimplePolygon_t; - return ret; + py::class_>(m, "simple_polygon") + .def(py::init(&P::polygon_from_numpy)) + .def("num_vertices", &P::num_vertices) + .def("array", + [] (const SimplePolygon& self) { return P::polygon_array(self); }, + py::return_value_policy::move) + .def("join", &P::template join_polygons) + .def("join", &P::template join_polygons) + .def("difference", &P::template difference_polygons) + .def("difference", &P::template difference_polygons) + .def("intersection", &P::template intersect_polygons) + .def("intersection", &P::template intersect_polygons); } -inline 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); - } +template +void bind_multi_polygon_class(py::module& m) { + using Polygon = typename P::Polygon_t; + using SimplePolygon = typename P::SimplePolygon_t; + using MultiPolygon = typename P::MultiPolygon_t; - return ret; + py::class_>(m, "multi_polygon") + .def(py::init( + [] (const Polygon& polygon) { + MultiPolygon self; + self.push_back(Polygon(polygon)); + return self; + })) + .def(py::init( + [] (const SimplePolygon& polygon) { + MultiPolygon self; + self.push_back(Polygon(polygon)); + return self; + })) + .def("num_parts", &MultiPolygon::size) + .def("parts", + [] (const MultiPolygon& self) { + py::list result; + for (auto p = self.begin(); p != self.end(); ++p) + result.append(*p); + return result; + }) + .def("extract_boundaries", [] (const MultiPolygon& self) { return P::extract_boundaries(self); }, + py::return_value_policy::take_ownership) + .def("join", + [] (const MultiPolygon a, SimplePolygon b) { + return P::join_multipolygons(a, MultiPolygon({static_cast(b)})); + }) + .def("join", + [] (const MultiPolygon a, Polygon b) { + return P::join_multipolygons(a, MultiPolygon({b})); + }) + .def("join", + [] (const MultiPolygon a, MultiPolygon b) { + return P::join_multipolygons(a, b); + }) + .def("__getitem__", [] (const MultiPolygon& self, int idx) { + return self.at(idx); + }, py::return_value_policy::reference_internal); } + +template +void bind_polygons(py::module& m) { + bind_polygon_class

(m); + bind_simple_polygon_class

(m); + bind_multi_polygon_class

(m); } +} // namespace rasputin + diff --git a/src/rasputin/cpp/polygon_bindings.cpp b/src/rasputin/cpp/polygon_bindings.cpp index c7c6c3d..5741a20 100644 --- a/src/rasputin/cpp/polygon_bindings.cpp +++ b/src/rasputin/cpp/polygon_bindings.cpp @@ -1,159 +1,15 @@ #include #include -#include -#include -#include - -#include "polygon.h" +#include "cgal_polygon.h" namespace py = pybind11; +PYBIND11_MAKE_OPAQUE(CGAL::Polygon); +PYBIND11_MAKE_OPAQUE(CGAL::SimplePolygon); PYBIND11_MAKE_OPAQUE(CGAL::MultiPolygon); -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; -} - - void bind_polygons(py::module& 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); + rasputin::bind_polygons>(m); } diff --git a/src/rasputin/cpp/raster_bindings.cpp b/src/rasputin/cpp/raster_bindings.cpp index c861ee1..14dd05c 100644 --- a/src/rasputin/cpp/raster_bindings.cpp +++ b/src/rasputin/cpp/raster_bindings.cpp @@ -2,68 +2,13 @@ #include #include -#include "raster_data.h" +#include "cgal_raster_data.h" namespace py = pybind11; PYBIND11_MAKE_OPAQUE(std::vector>); PYBIND11_MAKE_OPAQUE(std::vector>); -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); -} - -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); -} - void bind_raster(py::module& m) { - 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"); + rasputin::bind_raster(m); } diff --git a/src/rasputin/cpp/raster_data.h b/src/rasputin/cpp/raster_data.h index 872c098..45ebbc0 100644 --- a/src/rasputin/cpp/raster_data.h +++ b/src/rasputin/cpp/raster_data.h @@ -1,13 +1,45 @@ #pragma once -#include +#include + +#include +#include +#include #include "types.h" +#include "polygon.h" + +namespace py = pybind11; namespace rasputin { -template -struct RasterData { - RasterData( +namespace traits { + +template< + typename FT, + typename PointList, + typename PolygonUtils +> +requires std::floating_point +class RasterBase { + public: + using FT_t = FT; + using PointList_t = PointList; + using PolygonUtils_t = PolygonUtils; + using Polygon_t = typename PolygonUtils::Polygon_t; + using SimplePolygon_t = typename PolygonUtils::SimplePolygon_t; + using MultiPolygon_t = typename PolygonUtils::MultiPolygon_t; + + 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; + + RasterBase( double x_min, double y_max, double delta_x, @@ -25,27 +57,11 @@ struct RasterData { 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; + 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; } // For every point inside the raster rectangle we identify indices (i, j) of the upper-left vertex of the cell containing the point @@ -78,33 +94,105 @@ struct RasterData { 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)); - - return rectangle; - } - - // Compute intersection of raster rectangle with polygon - template - CGAL::MultiPolygon compute_intersection(const P& polygon) const { - CGAL::SimplePolygon rectangle = exterior(); + PointList raster_points() const; + SimplePolygon_t exterior() const; - CGAL::MultiPolygon intersection_polygon; - CGAL::intersection(rectangle, polygon, std::back_inserter(intersection_polygon)); - - return intersection_polygon; - } + template P> + MultiPolygon_t compute_intersection(const P& polygon) const; // 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)); + double eps = pow(pow(this->delta_x, 2) + pow(this->delta_y, 2), 0.5) * 1e-10; + return ((x > this->x_min + eps) and (x < this->get_x_max() - eps) + and (y > this->get_y_min() + eps) and (y < this->y_max - eps)); } }; + +template +concept RasterType = + requires { + typename RasterData::FT_t; + typename RasterData::PointList_t; + typename RasterData::Polygon_t; + typename RasterData::SimplePolygon_t; + typename RasterData::MultiPolygon_t; + typename RasterData::PolygonUtils_t; + } + && std::floating_point + && std::convertible_to< + RasterData, + RasterBase< + typename RasterData::FT_t, + typename RasterData::PointList_t, + typename RasterData::PolygonUtils_t + > + >; +} // namespace traits + +template +concept PolygonType = + traits::RasterType + && traits::PolygonType; + +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, 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); +} + +template +void bind_rasterdata(py::module &m, const std::string& pyname) { + using FT = typename RasterData::FT_t; + 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( + 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", &RasterData::x_min) + .def_readwrite("y_max", &RasterData::y_max) + .def_readwrite("delta_x", &RasterData::delta_x) + .def_readwrite("delta_y", &RasterData::delta_y) + .def_readonly("num_points_x", &RasterData::num_points_x) + .def_readonly("num_points_y", &RasterData::num_points_y) + .def_property_readonly("x_max", &RasterData::get_x_max) + .def_property_readonly("y_min", &RasterData::get_y_min) + .def("__getitem__", [](RasterData& self, std::pair idx) { + auto [i, j] = idx; + return self.data[self.num_points_x * i + j]; }) + .def("get_indices", &RasterData::get_indices) + .def("exterior", &RasterData::exterior, py::return_value_policy::take_ownership) + .def("contains", &RasterData::contains) + .def("get_interpolated_value_at_point", &RasterData::get_interpolated_value_at_point); +} + +template class RasterData> +requires traits::RasterType> && traits::RasterType> +void bind_raster(py::module& m) { + 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"); } +} // namespace rasputin diff --git a/src/rasputin/cpp/shade.h b/src/rasputin/cpp/shade.h index c4011ab..d8e7e13 100644 --- a/src/rasputin/cpp/shade.h +++ b/src/rasputin/cpp/shade.h @@ -1,156 +1,144 @@ -#include -#include -#include -#include - -#include -#include +#include #include +#include +#include +#include + #include "types.h" #include "mesh.h" #include "solar_position.h" +namespace py = pybind11; + namespace rasputin { -inline 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}; -} +namespace traits { +template +struct Shade { + static_assert("Shade trait not implemented!"); +}; +} // namespace traits -inline 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]); +template +struct Shade { + using S = traits::Shade; + using Mesh_t = Mesh; + using face_descriptor_t = typename Mesh::face_descriptor_t; - 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); - } - ++i; + template + static Point3 centroid(const Mesh_t& mesh, const face_descriptor_t& face) { + return S::centroid(mesh, face); } - return shade; -}; -inline 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; -} + static std::vector compute_shadow(const Mesh_t & mesh, const point3 &sun_direction) { + return S::compute_shadow(mesh, sun_direction); + } -inline std::vector 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)); + template + static bool is_shaded( + const Tree& tree, + const face_descriptor_t& fd, + const Vector& face_normal, + const Point3& face_center, + const double azimuth, + const double elevation + ) { + // return S::template is_shaded(tree, fd, face_normal, face_center, azimuth, elevation); + return S::is_shaded(tree, fd, face_normal, face_center, azimuth, elevation); } - return shade_vec; -} -inline 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]}); -}; + static std::vector compute_shadow(const Mesh_t &mesh, const double azimuth, const double 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]}); + } -inline 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]); + static std::vector> compute_shadows( + const Mesh_t& mesh, const std::vector>& sun_rays + ) { + // TODO: Rewrite totally! + std::vector> result; + return result; + /* + result.reserve(sun_rays.size()); + CGAL::Mesh cgal_mesh = mesh.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) + 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]); + } } + result.emplace_back(std::move(shade)); } - result.emplace_back(std::move(shade)); + return result; + */ } - return result; - */ }; + +template +concept ShadeType = + requires { + typename S::Mesh_t; + typename S::face_descriptor_t; + } + && + traits::MeshType; + +template +void bind_shades(py::module& m) { + using namespace std::chrono; + using Mesh = typename Shade::Mesh_t; + m.def("compute_shadow", + (std::vector (*)(const Mesh&, const rasputin::point3 &))&Shade::compute_shadow, + "Compute shadows for given topocentric sun position.") + .def("compute_shadow", + (std::vector (*)(const Mesh&, const double, const double))&Shade::compute_shadow, + "Compute shadows for given azimuth and elevation.") + // .def("compute_shadows", + // (std::vector> (*)(const Mesh&, const std::vector>&))&Shade::compute_shadows, + // "Compute shadows for a series of times and ray directions.") + .def("shade", [] (const Mesh& mesh, const double timestamp) { + const auto secs = seconds(int(std::round(timestamp))); + const auto millisecs = milliseconds(int(round(1000*fmod(timestamp, 1)))); + const auto tp = sys_days{1970y / January / 1} + secs + millisecs; + const auto shade_vec = Shade::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; + }); +} + +template +void bind_shades(py::module& m) { + bind_shades>(m); } +} // namespace rasputin diff --git a/src/rasputin/cpp/shade_bindings.cpp b/src/rasputin/cpp/shade_bindings.cpp new file mode 100644 index 0000000..277cfdc --- /dev/null +++ b/src/rasputin/cpp/shade_bindings.cpp @@ -0,0 +1,7 @@ +#include + +#include "cgal_shade.h" + +void bind_shades(py::module& m) { + rasputin::bind_shades(m); +} diff --git a/src/rasputin/cpp/solar_position_bindings.cpp b/src/rasputin/cpp/solar_position_bindings.cpp index 4f7b555..33f916f 100644 --- a/src/rasputin/cpp/solar_position_bindings.cpp +++ b/src/rasputin/cpp/solar_position_bindings.cpp @@ -1,8 +1,6 @@ #include -#include #include "solar_position.h" -#include "shade.h" namespace py = pybind11; @@ -22,20 +20,3 @@ void bind_solars(py::module& m) { }, "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."); } - -void bind_shades(py::module& m) { - 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("shade", [] (const rasputin::Mesh& mesh, const double timestamp) { - const auto secs = seconds(int(std::round(timestamp))); - const auto millisecs = milliseconds(int(round(1000*fmod(timestamp, 1)))); - const auto tp = sys_days{1970y / January / 1} + 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; - }); -}