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