diff --git a/.gitea/docker/Dockerfile.ci-cpp b/.gitea/docker/Dockerfile.ci-cpp new file mode 100644 index 0000000..e260036 --- /dev/null +++ b/.gitea/docker/Dockerfile.ci-cpp @@ -0,0 +1,16 @@ +FROM ubuntu:22.04 + +# Node.js 20 from NodeSource (Ubuntu Jammy ships v12 which is too old +# for actions/checkout@v4 — static class blocks require Node.js >= 16). +RUN apt-get update -qq && \ + apt-get install -y --no-install-recommends \ + curl ca-certificates && \ + curl -fsSL https://deb.nodesource.com/setup_20.x | bash - && \ + apt-get install -y --no-install-recommends \ + nodejs \ + cmake \ + g++ \ + git \ + && rm -rf /var/lib/apt/lists/* + +WORKDIR /workspace diff --git a/.gitea/workflows/cpp-tests.yml b/.gitea/workflows/cpp-tests.yml new file mode 100644 index 0000000..722093a --- /dev/null +++ b/.gitea/workflows/cpp-tests.yml @@ -0,0 +1,42 @@ +name: C++ Tests + +on: + push: + branches: + - main + - dev + - "claude/**" + pull_request: + +jobs: + test: + runs-on: eulernest + container: + image: git.eulernest.eu/conformallab/ci-cpp:latest + + steps: + - uses: actions/checkout@v4 + + - name: Configure (tests-only mode) + run: cmake -S code -B build -DCMAKE_BUILD_TYPE=Release + + - name: Build test binary + run: cmake --build build --target conformallab_tests -j$(nproc) + + - name: Run tests + run: > + ctest --test-dir build + --output-on-failure + --output-junit test-results.xml + + - name: Show test summary + if: always() + run: | + if [ -f test-results.xml ]; then + total=$(grep -o 'tests="[0-9]*"' test-results.xml | grep -o '[0-9]*' | head -1) + failed=$(grep -o 'failures="[0-9]*"' test-results.xml | grep -o '[0-9]*' | head -1) + skipped=$(grep -o 'skipped="[0-9]*"' test-results.xml | grep -o '[0-9]*' | head -1) + passed=$(( ${total:-0} - ${failed:-0} - ${skipped:-0} )) + echo "" + echo "TOTAL: ${total:-0} | PASSED: $passed | FAILED: ${failed:-0} | SKIPPED: ${skipped:-0}" + fi diff --git a/README.md b/README.md index 37167d3..16acc76 100644 --- a/README.md +++ b/README.md @@ -1,46 +1,104 @@ # conformallab++ -conformallab++ is a modern C++ reimplementation of the ConformalLab software by Stefan Sechelmann for experiments in discrete conformal geometry and related mesh transformations. - +conformallab++ is a modern C++ reimplementation of the [ConformalLab](https://github.com/sechel/conformallab) software by Stefan Sechelmann for experiments in discrete conformal geometry and related mesh transformations. -It builds on well‑established open‑source libraries such as **CGAL**, **Eigen**, and **Boost** for robust geometric data structures and efficient numerical computations, and uses single‑header libraries **CLI11** and **json.hpp** for a lightweight command‑line interface and configuration handling. In addition, **libigl**’s OpenGL viewer stack based on **GLFW** provides portable windowing and input handling, while the separate **libigl‑glad** component supplies a generated OpenGL function loader for the required core profile, making interactive visualization easy to integrate. +> **Status:** early prototype stage. API, file formats, and CLI are subject to change. +## Features -## Status +- Discrete conformal geometry utilities (Clausen function, hyper-ideal tetrahedra, surface curves) +- Mesh I/O and conversion using CGAL — optional, only needed for the CLI app +- Linear algebra routines with Eigen +- Interactive mesh viewer using libigl / GLFW — optional +- Lightweight CLI with CLI11 and JSON configuration -This project is in an early prototype stage. -API, file formats and command‑line interface are all subject to change. +## Build modes -## Features (placeholder) +The project uses three clearly separated CMake modes so you only pull in what you need. -- Discrete conformal geometry utilities (placeholder) -- Mesh and graph data structures built on CGAL (placeholder) -- Linear algebra and optimization routines using Eigen (placeholder) -- High‑level operations and helpers based on Boost (placeholder) -- Command‑line tools with CLI11 and JSON configuration (placeholder) +| Mode | CMake flag | What gets built | +|------|-----------|-----------------| +| **Tests only** (default, used in CI) | *(none)* | `conformallab_tests` · deps: Eigen + GTest | +| **Viewer** | `-DWITH_VIEWER=ON` | `viewer` library · deps: libigl / GLFW / GLAD + Eigen | +| **Full app** | `-DWITH_CGAL=ON` | `conformallab_core` CLI + viewer · deps: CGAL + libigl / GLFW / GLAD + Eigen | + +`-DWITH_CGAL=ON` automatically enables `WITH_VIEWER` because the CLI app uses the viewer library for mesh visualisation. + +External dependencies ship as tarballs in `code/deps/tarballs/` and are extracted lazily at CMake configure time — no internet access needed after cloning (GTest is the only exception: fetched from GitHub via FetchContent). + +## Prerequisites + +| Tool | Minimum version | +|------|----------------| +| C++ compiler (GCC or Clang) | C++17 | +| CMake | 3.20 | + +No system-level libraries are required for the default tests-only build. CGAL and libigl are header-only and bundled in the repo. ## Getting started -### Prerequisites - -- A C++20 compiler (e.g. `g++` or `clang++`) -- CMake version 3.20 - -### Clone the repository - ```bash -git clone https://codeberg.org/user2595/ConformalLabpp -cd conformallabpp/code +git clone https://codeberg.org/TMoussa/ConformalLabpp +cd ConformalLabpp ``` -### Configure and build +### Tests only (CI default) + ```bash -cmake -S . -B build && cmake --build build +cmake -S code -B build +cmake --build build --target conformallab_tests -j$(nproc) +ctest --test-dir build --output-on-failure ``` -### Run the binare +### Full CLI app (CGAL + viewer) + ```bash - ./bin/conformallab_core --input ./data/off/simple_cupe.off +cmake -S code -B build -DWITH_CGAL=ON +cmake --build build -j$(nproc) +./code/bin/conformallab_core --input data/off/example.off --show ``` -### License -conformallabpp is released under the MIT License (see LICENSE). + +### Viewer only (no CGAL) + +```bash +cmake -S code -B build -DWITH_VIEWER=ON +cmake --build build --target viewer -j$(nproc) +``` + +## Project structure + +``` +code/ +├── include/ # Public headers (Clausen, hyper-ideal, mesh utils, …) +├── src/ +│ ├── apps/v0/ # conformallab_core CLI app (requires WITH_CGAL) +│ └── viewer/ # simple_viewer (requires WITH_VIEWER) +├── tests/ # GTest unit tests (always built) +└── deps/ + ├── tarballs/ # Bundled dependency archives + ├── eigen-3.4.0/ # Header-only linear algebra (always extracted) + ├── CGAL-6.1.1/ # Header-only geometry (extracted with WITH_CGAL) + ├── libigl-2.6.0/ # Header-only viewer toolkit (extracted with WITH_VIEWER) + ├── glfw-3.4/ # Windowing (extracted with WITH_VIEWER) + ├── libigl-glad/ # OpenGL loader (extracted with WITH_VIEWER) + └── single_includes/ # CLI11, json.hpp +``` + +## CI + +Tests run automatically on push to `main`, `dev`, and `claude/**` branches via a self-hosted Gitea Actions runner (`eulernest`, ARM64 Raspberry Pi). The pipeline uses a minimal Docker image (`git.eulernest.eu/conformallab/ci-cpp:latest`) with cmake, g++, git, and Node.js 20 pre-installed. + +The Dockerfile for the CI image lives in `.gitea/docker/Dockerfile.ci-cpp`. Build and push it once whenever the image needs updating: + +```bash +docker buildx build \ + --platform linux/arm64 \ + -f .gitea/docker/Dockerfile.ci-cpp \ + -t git.eulernest.eu/conformallab/ci-cpp:latest \ + --push \ + .gitea/docker/ +``` + +## License + +conformallab++ is released under the MIT License (see [LICENSE](LICENSE)). diff --git a/code/.gitignore b/code/.gitignore index 264b0e8..6619b8d 100644 --- a/code/.gitignore +++ b/code/.gitignore @@ -1,5 +1,6 @@ # Build artifacts build/ +build_release/ *.o *.so *.dylib diff --git a/code/CMakeLists.txt b/code/CMakeLists.txt index e2faa85..d409480 100644 --- a/code/CMakeLists.txt +++ b/code/CMakeLists.txt @@ -1,13 +1,31 @@ cmake_minimum_required(VERSION 3.20) -# Automatisch Projektname aus Parent-Ordner set(PROJECT_NAME "conformallab_core") - project(${PROJECT_NAME} LANGUAGES C CXX) message(STATUS "Configuring ${PROJECT_NAME}...") -# Build-Setup +# ── Build modes ──────────────────────────────────────────────────────────────── +# +# Default (CI / tests-only): only Eigen + GTest are required. +# +# -DWITH_CGAL=ON builds the conformallab_core CLI app (needs CGAL). +# Automatically enables WITH_VIEWER because the app uses +# the viewer library for mesh visualisation. +# +# -DWITH_VIEWER=ON builds the viewer library standalone (libigl/GLFW/GLAD). +# +# ────────────────────────────────────────────────────────────────────────────── +option(WITH_CGAL "Build conformallab_core app (requires CGAL + Viewer)" OFF) +option(WITH_VIEWER "Build viewer library (libigl / GLFW / GLAD)" OFF) + +# The CLI app always needs the viewer; enable it implicitly. +if(WITH_CGAL AND NOT WITH_VIEWER) + message(STATUS "WITH_CGAL implies WITH_VIEWER – enabling automatically.") + set(WITH_VIEWER ON CACHE BOOL "" FORCE) +endif() + +# ── Standard settings ────────────────────────────────────────────────────────── set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) @@ -18,55 +36,69 @@ if(NOT CMAKE_BUILD_TYPE) set(CMAKE_BUILD_TYPE "Release" CACHE STRING "Build type" FORCE) endif() -# Compiler-Warnings & Optimierungen if(CMAKE_CXX_COMPILER_ID MATCHES "Clang|GNU") add_compile_options(-Wall -Wextra -Wpedantic) - - # AddressSanitizer nur im Debug - if(CMAKE_BUILD_TYPE STREQUAL "Debug") + + # AddressSanitizer only in Debug (gtest_discover_tests runs the binary at + # configure time and hangs with ASan enabled). + if(CMAKE_BUILD_TYPE STREQUAL "Debug" AND NOT BUILD_TESTING) add_compile_options(-fsanitize=address -fno-omit-frame-pointer) add_link_options(-fsanitize=address) - endif() + endif() endif() +# ── GTest (always – tests are always built) ──────────────────────────────────── +include(FetchContent) +include(CTest) +enable_testing() +FetchContent_Declare( + googletest + GIT_REPOSITORY https://github.com/google/googletest.git + GIT_TAG v1.14.0 +) +set(gtest_force_shared_crt ON CACHE BOOL "" FORCE) +set(INSTALL_GTEST OFF CACHE BOOL "" FORCE) +set(BUILD_GMOCK OFF CACHE BOOL "" FORCE) +FetchContent_MakeAvailable(googletest) + +# ── External deps (lazy tarball extraction) ──────────────────────────────────── add_subdirectory(deps) -add_subdirectory(deps/glfw-3.4) +# ── Viewer library (optional) ────────────────────────────────────────────────── +if(WITH_VIEWER) + add_subdirectory(deps/glfw-3.4) + add_library(glad STATIC + ${CMAKE_CURRENT_SOURCE_DIR}/deps/libigl-glad/src/glad.c) + target_include_directories(glad PUBLIC + ${CMAKE_CURRENT_SOURCE_DIR}/deps/libigl-glad/include) -# --------- Executable --------- -add_executable(${PROJECT_NAME} src/apps/v0/conformallab_cli.cpp ) + add_library(viewer STATIC src/viewer/simple_viewer.cpp) + target_include_directories(viewer PUBLIC + ${CMAKE_CURRENT_SOURCE_DIR}/deps/libigl-2.6.0/include + ${CMAKE_CURRENT_SOURCE_DIR}/include + ${CMAKE_CURRENT_SOURCE_DIR}/deps/libigl-glad/include + ${CMAKE_CURRENT_SOURCE_DIR}/deps/eigen-3.4.0/) + target_link_libraries(viewer PUBLIC glad glfw) +endif() -target_include_directories(${PROJECT_NAME} PRIVATE - ${CMAKE_CURRENT_SOURCE_DIR}/include - ${CMAKE_CURRENT_SOURCE_DIR}/deps/single_includes - ${CMAKE_CURRENT_SOURCE_DIR}/deps/eigen-3.4.0/ - ${CMAKE_CURRENT_SOURCE_DIR}/deps/CGAL-6.1.1/include - ${CMAKE_CURRENT_SOURCE_DIR}/deps/boost_1_90_0/ - ${CMAKE_CURRENT_SOURCE_DIR}/deps/libigl-2.6.0/include - ${CMAKE_CURRENT_SOURCE_DIR}/deps/libigl-glad/include -) - - -add_library(glad STATIC - ${CMAKE_CURRENT_SOURCE_DIR}/deps/libigl-glad/src/glad.c -) -target_include_directories(glad PUBLIC - ${CMAKE_CURRENT_SOURCE_DIR}/deps/libigl-glad/include -) - -target_compile_definitions(${PROJECT_NAME} PRIVATE - CGAL_DISABLE_GMP - CGAL_DISABLE_MPFR -) - -# Output-Ordner -set_target_properties(${PROJECT_NAME} PROPERTIES - RUNTIME_OUTPUT_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/bin - -) -target_link_libraries(${PROJECT_NAME} - PRIVATE glad glfw -) +# ── Core CLI app (optional, requires CGAL + Viewer) ─────────────────────────── +if(WITH_CGAL) + add_executable(${PROJECT_NAME} src/apps/v0/conformallab_cli.cpp) + target_include_directories(${PROJECT_NAME} SYSTEM PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/deps/single_includes + ${CMAKE_CURRENT_SOURCE_DIR}/deps/eigen-3.4.0/ + ${CMAKE_CURRENT_SOURCE_DIR}/deps/CGAL-6.1.1/include + ${CMAKE_CURRENT_SOURCE_DIR}/deps/libigl-2.6.0/include) + target_include_directories(${PROJECT_NAME} PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/include) + target_compile_definitions(${PROJECT_NAME} PRIVATE + CGAL_DISABLE_GMP CGAL_DISABLE_MPFR) + target_link_libraries(${PROJECT_NAME} PRIVATE viewer) + set_target_properties(${PROJECT_NAME} PROPERTIES + RUNTIME_OUTPUT_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/bin) +endif() +# ── Tests (always) ──────────────────────────────────────────────────────────── +add_subdirectory(tests) diff --git a/code/deps/Cmakelists.txt b/code/deps/Cmakelists.txt index d6a0910..025c852 100644 --- a/code/deps/Cmakelists.txt +++ b/code/deps/Cmakelists.txt @@ -1,5 +1,16 @@ # deps/CMakeLists.txt -function(setup_dependency NAME SUBDIR) +# +# Lazy tarball extraction – a dependency is only extracted when its directory +# does not yet exist. Extraction happens at cmake configure time. +# +# Which deps are extracted depends on the active build mode: +# +# (default) Eigen only → tests-only build +# WITH_VIEWER=ON + Eigen, libigl, libigl-glad, glfw +# WITH_CGAL=ON + Eigen, CGAL (WITH_VIEWER is implied by WITH_CGAL) +# + +function(setup_dependency NAME SUBDIR) if(NOT EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/${NAME}") file(MAKE_DIRECTORY "${NAME}") file(GLOB TARBALL "${CMAKE_CURRENT_SOURCE_DIR}/tarballs/${NAME}.tar.*") @@ -7,20 +18,21 @@ function(setup_dependency NAME SUBDIR) if(NOT TARBALL) message(FATAL_ERROR "❌ No '${NAME}*.tar.*' in deps/tarballs/") else() - message(STATUS "${CMAKE_COMMAND} ${TARBALL}" ) - + message(STATUS "${CMAKE_COMMAND} ${TARBALL}") if(SUBDIR STREQUAL "") - execute_process(COMMAND "${CMAKE_COMMAND}" -E tar xf "${TARBALL}" RESULT_VARIABLE TAR_STATUS WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}" ) - message(STATUS "Extracting ${NAME}... Output: ${TAR_OUTPUT} Status: ${TAR_STATUS}") - if(NOT TAR_STATUS EQUAL 0) - message(FATAL_ERROR "❌ Extract failed: ${NAME}") - endif() + execute_process( + COMMAND "${CMAKE_COMMAND}" -E tar xf "${TARBALL}" + RESULT_VARIABLE TAR_STATUS + WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}") else() - execute_process(COMMAND "${CMAKE_COMMAND}" -E tar xf "${TARBALL}" "${NAME}/${SUBDIR}" RESULT_VARIABLE TAR_STATUS WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}" ) - message(STATUS "Extracting ${NAME}... Output: ${TAR_OUTPUT} Status: ${TAR_STATUS}") - if(NOT TAR_STATUS EQUAL 0) - message(FATAL_ERROR "❌ Extract failed: ${NAME}") - endif() + execute_process( + COMMAND "${CMAKE_COMMAND}" -E tar xf "${TARBALL}" "${NAME}/${SUBDIR}" + RESULT_VARIABLE TAR_STATUS + WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}") + endif() + message(STATUS "Extracting ${NAME}... Status: ${TAR_STATUS}") + if(NOT TAR_STATUS EQUAL 0) + message(FATAL_ERROR "❌ Extract failed: ${NAME}") endif() endif() else() @@ -28,10 +40,17 @@ function(setup_dependency NAME SUBDIR) endif() endfunction() -setup_dependency("CGAL-6.1.1" "include" ) -setup_dependency("eigen-3.4.0" "Eigen" ) -setup_dependency("boost_1_90_0" "boost" ) -setup_dependency("libigl-2.6.0" "include" ) -setup_dependency("libigl-glad" "" ) -setup_dependency("glfw-3.4" "" ) +# ── Always required ──────────────────────────────────────────────────────────── +setup_dependency("eigen-3.4.0" "Eigen") +# ── Viewer mode ──────────────────────────────────────────────────────────────── +if(WITH_VIEWER) + setup_dependency("libigl-2.6.0" "include") + setup_dependency("libigl-glad" "") + setup_dependency("glfw-3.4" "") +endif() + +# ── CGAL mode ───────────────────────────────────────────────────────────────── +if(WITH_CGAL) + setup_dependency("CGAL-6.1.1" "include") +endif() diff --git a/code/deps/tarballs/boost_1_90_0.tar.gz b/code/deps/tarballs/boost_1_90_0.tar.gz deleted file mode 100644 index 1e830b4..0000000 Binary files a/code/deps/tarballs/boost_1_90_0.tar.gz and /dev/null differ diff --git a/code/include/clausen.hpp b/code/include/clausen.hpp new file mode 100644 index 0000000..29e7c5d --- /dev/null +++ b/code/include/clausen.hpp @@ -0,0 +1,179 @@ +#pragma once + +// Clausen integral, Lobachevsky function, and Im(Li2). +// Ported from de.varylab.discreteconformal.functional.Clausen (Java). +// Original algorithm: Boris Springborn, Stefan Sechelmann (TU Berlin). + +#include +#include +#include + +namespace conformallab { + +namespace detail { + +// Evaluate a Chebyshev series at x using n terms. +// Corresponds to Java Clausen.csevl(). +inline double csevl(double x, const double* cs, int n) noexcept { + double b2 = 0.0, b1 = 0.0, b0 = 0.0; + const double twox = 2.0 * x; + while (n-- > 0) { + b2 = b1; + b1 = b0; + b0 = twox * b1 - b2 + cs[n]; + } + return 0.5 * (b0 - b2); +} + +// Count Chebyshev terms needed so truncation error <= eta. +// Corresponds to Java Clausen.inits(). +inline int inits(const double* series, int n, double eta) noexcept { + double err = 0.0; + while (err <= eta && n-- != 0) { + err += std::abs(series[n]); + } + return n; +} + +// Chebyshev expansion of `cl(t)/t + log(t)` around Pi/6. +inline const double* clpi6() noexcept { + static const double a[] = { + 2 * 1.0057346496467363858, + .0076523796971586786263, + .0019223823523180480014, + .53333368801173950429e-5, + .68684944849366102659e-6, + .63769755654413855855e-8, + .57069363812137970721e-9, + .87936343137236194448e-11, + .62365831120408524691e-12, + .12996625954032513221e-13, + .78762044080566097484e-15, + .20080243561666612900e-16, + .10916495826127475499e-17, + .32027217200949691956e-19, + }; + return a; +} + +// Chebyshev expansion around Pi/2. +inline const double* clpi2() noexcept { + static const double a[] = { + 2 * .017492908851746863924 + 2 * 1.0057346496467363858, + .023421240075284860656 + .0076523796971586786263, + .0060025281630108248332 + .0019223823523180480014, + .000085934211448718844330 + .53333368801173950429e-5, + .000012155033501044820317 + .68684944849366102659e-6, + .46587486310623464413e-6 + .63769755654413855855e-8, + .50732554559130493329e-7 + .57069363812137970721e-9, + .28794458754760053792e-8 + .87936343137236194448e-11, + .27792370776596244150e-9 + .62365831120408524691e-12, + .19340423475636663004e-10 + .12996625954032513221e-13, + .17726134256574610202e-11 + .78762044080566097484e-15, + .13811355237660945692e-12 + .20080243561666612900e-16, + .12433074161771699487e-13 + .10916495826127475499e-17, + .10342683357723940535e-14 + .32027217200949691956e-19, + .92910354101990447850e-16, + .80428334724548559541e-17, + .72598441354406482972e-18, + .64475701884829384587e-19, + .58630185185185185187e-20, + }; + return a; +} + +// Chebyshev expansion of `-cl(Pi-t)/(Pi-t) + log(2)` around 5*Pi/6. +inline const double* cl5pi6() noexcept { + static const double a[] = { + 2 * .017492908851746863924, + .023421240075284860656, + .0060025281630108248332, + .000085934211448718844330, + .000012155033501044820317, + .46587486310623464413e-6, + .50732554559130493329e-7, + .28794458754760053792e-8, + .27792370776596244150e-9, + .19340423475636663004e-10, + .17726134256574610202e-11, + .13811355237660945692e-12, + .12433074161771699487e-13, + .10342683357723940535e-14, + .92910354101990447850e-16, + .80428334724548559541e-17, + .72598441354406482972e-18, + .64475701884829384587e-19, + .58630185185185185187e-20, + }; + return a; +} + +// Machine epsilon for double (IEEE 754). +constexpr double kMachEps = 1.1102230246251565e-16; +constexpr double kLn2 = 0.693147180559945309417232121458; + +// Number of Chebyshev terms for each expansion (computed like Java static init). +inline int nclpi6() noexcept { + static const int n = inits(clpi6(), 14, kMachEps / 10.0); + return n; +} +inline int nclpi2() noexcept { + static const int n = inits(clpi2(), 19, kMachEps / 10.0); + return n; +} +inline int ncl5pi6() noexcept { + static const int n = inits(cl5pi6(), 19, kMachEps / 10.0); + return n; +} + +} // namespace detail + +// Clausen's integral Cl2(x) = -integral_0^x log|2 sin(t/2)| dt. +// High-precision Chebyshev implementation. +// Corresponds to Java Clausen.clausen2(). +inline double clausen2(double x) noexcept { + using namespace detail; + constexpr double pi = 3.14159265358979323846264338328; + constexpr double two_pi = 2.0 * pi; + + int rh = 0; + x = std::fmod(x, two_pi); + if (x < 0.0) x += two_pi; + if (x > pi) { rh = 1; x = two_pi - x; } + + double f; + if (x == 0.0) { + f = 0.0; + } else if (x <= pi / 3.0) { + f = csevl(x * (6.0 / pi) - 1.0, clpi6(), nclpi6()) * x - x * std::log(x); + } else if (x <= 2.0 * pi / 3.0) { + f = csevl(x * (3.0 / pi) - 1.0, clpi2(), nclpi2()) * x - x * std::log(x); + } else { + f = (kLn2 - csevl(5.0 - x * (6.0 / pi), cl5pi6(), ncl5pi6())) * (pi - x); + } + return rh ? -f : f; +} + +// Milnor's Lobachevsky function Л(x) = Cl2(2x)/2. +// Corresponds to Java Clausen.Л(). +inline double Lobachevsky(double x) noexcept { + constexpr double pi = 3.14159265358979323846264338328; + x = std::fmod(x, pi); + if (x <= 0.0) x += pi; + return clausen2(2.0 * x) / 2.0; +} + +// Imaginary part of the dilogarithm Im(Li2(z)). +// Corresponds to Java Clausen.ImLi2(). +inline double ImLi2(std::complex z) noexcept { + auto a = std::log(1.0 - std::conj(z)); // log(1 - conj(z)) + auto b = std::log(1.0 - z); // log(1 - z) + double x = std::log(std::abs(z)); + double y = 0.5 * (a - b).imag(); + double phi = std::arg(z); + return y * x + 0.5 * (clausen2(2.0 * y) + - clausen2(2.0 * (y + phi)) + + clausen2(2.0 * phi)); +} + +} // namespace conformallab diff --git a/code/include/example_include.hpp b/code/include/example_include.hpp deleted file mode 100644 index f62204e..0000000 --- a/code/include/example_include.hpp +++ /dev/null @@ -1,7 +0,0 @@ -//exmaple include hedder -#ifndef EXAMPLE_INCLUDE_HPP -#define EXAMPLE_INCLUDE_HPP - - - -#endif // EXAMPLE_INCLUDE_HPP \ No newline at end of file diff --git a/code/include/hyper_ideal_utility.hpp b/code/include/hyper_ideal_utility.hpp new file mode 100644 index 0000000..a3c01d2 --- /dev/null +++ b/code/include/hyper_ideal_utility.hpp @@ -0,0 +1,108 @@ +#pragma once + +// Hyperbolic tetrahedron volume formulas. +// Ported from de.varylab.discreteconformal.functional.HyperIdealUtility (Java). + +#include "clausen.hpp" + +#include +#include +#include + +namespace conformallab { + +// Volume of a generalized hyperbolic tetrahedron with dihedral angles A..F. +// Formula: Meyerhoff / Ushijima (Springer 2006). +// Corresponds to Java HyperIdealUtility.calculateTetrahedronVolume(). +inline double calculateTetrahedronVolume(double A, double B, double C, + double D, double E, double F) { + constexpr double pi = 3.14159265358979323846264338328; + + // Degenerate if any angle equals pi. + if (A == pi || B == pi || C == pi || D == pi || E == pi || F == pi) + return 0.0; + + const double sA = std::sin(A), sB = std::sin(B), sC = std::sin(C); + const double sD = std::sin(D), sE = std::sin(E), sF = std::sin(F); + const double cA = std::cos(A), cB = std::cos(B), cC = std::cos(C); + const double cD = std::cos(D), cE = std::cos(E), cF = std::cos(F); + + // Unit complex numbers e^(i*angle). + using Cx = std::complex; + auto polar = [](double angle) { return std::polar(1.0, angle); }; + + Cx ad = polar(A + D), be = polar(B + E), cf = polar(C + F); + Cx abc = polar(A + B + C), abf = polar(A + B + F); + Cx ace = polar(A + C + E), aef = polar(A + E + F); + Cx bcd = polar(B + C + D), bdf = polar(B + D + F); + Cx def = polar(D + E + F), cde = polar(C + D + E); + Cx abde = ad * be, acdf = ad * cf, bcef = be * cf; + Cx abcdef = abc * def; + + Cx z = ad + be + cf + abf + ace + bcd + def + abcdef; + + // Gram matrix of the tetrahedron. + Eigen::Matrix4d G; + G << 1.0, -cA, -cB, -cF, + -cA, 1.0, -cC, -cE, + -cB, -cC, 1.0, -cD, + -cF, -cE, -cD, 1.0; + + Cx sqrtG = std::sqrt(Cx(G.determinant(), 0.0)); + + Cx f = Cx(sA*sD + sB*sE + sC*sF, 0.0); + Cx f1 = f - sqrtG; + Cx f2 = f + sqrtG; + Cx z1 = -2.0 * f1 / z; + Cx z2 = -2.0 * f2 / z; + + auto U = [&](Cx zi) { + return 0.5 * ( + + ImLi2(zi) + + ImLi2(abde * zi) + + ImLi2(acdf * zi) + + ImLi2(bcef * zi) + - ImLi2(-abc * zi) + - ImLi2(-aef * zi) + - ImLi2(-bdf * zi) + - ImLi2(-cde * zi) + ); + }; + + return (U(z1) - U(z2)) / 2.0; +} + +// Volume of a hyperideal tetrahedron with one ideal vertex (at gamma). +// Dihedral angles at the ideal vertex: gamma1, gamma2, gamma3. +// Dihedral angles at opposite edges: alpha23, alpha31, alpha12. +// Formula: Kolpakov–Mednykh (arxiv math/0603097). +// Corresponds to Java HyperIdealUtility.calculateTetrahedronVolumeWithIdealVertexAtGamma(). +inline double calculateTetrahedronVolumeWithIdealVertexAtGamma( + double gamma1, double gamma2, double gamma3, + double alpha23, double alpha31, double alpha12) +{ + constexpr double pi = 3.14159265358979323846264338328; + auto L = [](double x) { return Lobachevsky(x); }; + + double result = L(gamma1) + L(gamma2) + L(gamma3); + + result += L((pi + alpha31 - alpha12 - gamma1) / 2.0); + result += L((pi + alpha12 - alpha23 - gamma2) / 2.0); + result += L((pi + alpha23 - alpha31 - gamma3) / 2.0); + + result += L((pi - alpha31 + alpha12 - gamma1) / 2.0); + result += L((pi - alpha12 + alpha23 - gamma2) / 2.0); + result += L((pi - alpha23 + alpha31 - gamma3) / 2.0); + + result += L((pi + alpha31 + alpha12 - gamma1) / 2.0); + result += L((pi + alpha12 + alpha23 - gamma2) / 2.0); + result += L((pi + alpha23 + alpha31 - gamma3) / 2.0); + + result += L((pi - alpha31 - alpha12 - gamma1) / 2.0); + result += L((pi - alpha12 - alpha23 - gamma2) / 2.0); + result += L((pi - alpha23 - alpha31 - gamma3) / 2.0); + + return result / 2.0; +} + +} // namespace conformallab diff --git a/code/include/matrix_utility.hpp b/code/include/matrix_utility.hpp new file mode 100644 index 0000000..f18aee6 --- /dev/null +++ b/code/include/matrix_utility.hpp @@ -0,0 +1,21 @@ +#pragma once + +// 4x4 mapping matrix from corresponding point pairs. +// Ported from de.varylab.discreteconformal.math.MatrixUtility (Java). + +#include + +namespace conformallab { + +// Find the 4×4 matrix R that maps source points to target points. +// Each row of `from` / `to` is a homogeneous 4-vector (one point per row). +// Post-condition: R * from.row(i).T == to.row(i).T for all i. +// +// Implementation: R = to^T * (from^T)^{-1} +// Corresponds to Java MatrixUtility.makeMappingMatrix(). +inline Eigen::Matrix4d makeMappingMatrix(const Eigen::Matrix4d& from, + const Eigen::Matrix4d& to) { + return to.transpose() * from.transpose().inverse(); +} + +} // namespace conformallab diff --git a/code/include/mesh_utils.hpp b/code/include/mesh_utils.hpp new file mode 100644 index 0000000..fc4e936 --- /dev/null +++ b/code/include/mesh_utils.hpp @@ -0,0 +1,47 @@ +#pragma once +#include +#include +#include + + +namespace mesh_utils { +template +void cgal_to_eigen(CGAL::Surface_mesh& mesh, + Eigen::MatrixXd& V, Eigen::MatrixXi& F) { + + CGAL::Polygon_mesh_processing::triangulate_faces(mesh); + + V.resize(mesh.num_vertices(), 3); + F.resize(mesh.num_faces(), 3); + + for (auto v : mesh.vertices()) { + auto p = mesh.point(v); + V.row(v.idx()) << p.x(), p.y(), p.z(); + } + + Eigen::Index face_idx = 0; + for (auto f : mesh.faces()) { + int vertex_count = 0; + for (auto h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) { + auto v = mesh.target(h); + F(face_idx, vertex_count) = v.idx(); + ++vertex_count; + } + face_idx++; + } +} +template +void simple_visualize_mesh(Eigen::MatrixXd& V, Eigen::MatrixXi& F) { + igl::opengl::glfw::Viewer viewer; + viewer.data().set_mesh(V, F); + viewer.launch(); +} +// Zero-Copy Map für V (optional) +template +Eigen::Map> +get_vertex_map(CGAL::Surface_mesh& mesh) { + auto& points = mesh.points(); + double* data = reinterpret_cast(&points[0][0]); + return {data, static_cast(points.size()), 3}; +} +} // namespace diff --git a/code/include/projective_math.hpp b/code/include/projective_math.hpp new file mode 100644 index 0000000..2c8f297 --- /dev/null +++ b/code/include/projective_math.hpp @@ -0,0 +1,89 @@ +#pragma once + +// Projective and hyperbolic geometry utilities. +// Ported from de.jreality.math.Pn / Rn and +// de.varylab.discreteconformal.uniformization.SurfaceCurveUtility (Java). + +#include +#include +#include +#include +#include + +namespace conformallab { + +// Divide a homogeneous vector by its last component. +// Corresponds to Java Pn.dehomogenize(). +inline Eigen::VectorXd dehomogenize(const Eigen::VectorXd& p) { + return p / p(p.size() - 1); +} + +// Hyperbolic distance between two homogeneous vectors of the same dimension. +// The last component is the "timelike" coordinate (jReality convention). +// Inner product: = -sum_i p_i*q_i + p_last * q_last +// Distance: arcosh() where p̂ normalises to the hyperboloid. +// Corresponds to Java Pn.distanceBetween(p, q, Pn.HYPERBOLIC). +inline double hyperbolicDistance(const Eigen::VectorXd& p, + const Eigen::VectorXd& q) { + int n = static_cast(p.size()); + double normP = std::sqrt(p(n-1)*p(n-1) - p.head(n-1).squaredNorm()); + double normQ = std::sqrt(q(n-1)*q(n-1) - q.head(n-1).squaredNorm()); + double inner = (-p.head(n-1).dot(q.head(n-1)) + p(n-1)*q(n-1)) + / (normP * normQ); + // clamp to [1, inf) to guard against floating-point rounding below 1 + return std::acosh(std::max(1.0, inner)); +} + +// Check whether a homogeneous point p lies on the segment [s[0], s[1]]. +// Works for n-dimensional homogeneous coords; cross product uses the first +// 3 spatial components after dehomogenization (matching jReality's Rn behaviour). +// Corresponds to Java SurfaceCurveUtility.isOnSegment(). +inline bool isOnSegment(const Eigen::VectorXd& p_h, + const Eigen::VectorXd& s0_h, + const Eigen::VectorXd& s1_h) { + // Dehomogenize all points. + Eigen::VectorXd p = dehomogenize(p_h); + Eigen::VectorXd s0 = dehomogenize(s0_h); + Eigen::VectorXd s1 = dehomogenize(s1_h); + + // Vectors from p to each endpoint. + Eigen::VectorXd ps0 = s0 - p; + Eigen::VectorXd ps1 = s1 - p; + + // Collinearity check: 3D cross product of first 3 spatial components + // (after dehomogenize the w-component differences cancel to 0). + // head<3>() gives compile-time size needed by Eigen's cross(). + Eigen::Vector3d cross = ps0.head<3>().cross(ps1.head<3>()); + if (cross.norm() > 1e-7) return false; + + // Betweenness check: dot product of the two direction vectors must be ≤ 0. + double dot = ps0.dot(ps1); + if (dot > 0.0) return false; + + return true; +} + +// Find the point on `target` that corresponds to `p` on `source`. +// The parameter t is determined by hyperbolic distance ratios on `source`, +// then applied as a linear interpolation on the dehomogenized `target`. +// Corresponds to Java SurfaceCurveUtility.getPointOnCorrespondingSegment(). +inline Eigen::VectorXd getPointOnCorrespondingSegment( + const Eigen::VectorXd& p, + const Eigen::VectorXd& src0, + const Eigen::VectorXd& src1, + const Eigen::VectorXd& tgt0, + const Eigen::VectorXd& tgt1) +{ + double l = hyperbolicDistance(src0, src1); + double l1 = hyperbolicDistance(src0, p) / l; // weight for tgt1 + double l2 = hyperbolicDistance(src1, p) / l; // weight for tgt0 + + if (std::isnan(l1)) return dehomogenize(tgt0); + if (std::isnan(l2)) return dehomogenize(tgt1); + + Eigen::VectorXd t0d = dehomogenize(tgt0); + Eigen::VectorXd t1d = dehomogenize(tgt1); + return l1 * t1d + l2 * t0d; +} + +} // namespace conformallab diff --git a/code/include/viewer_utils.h b/code/include/viewer_utils.h new file mode 100644 index 0000000..4140d11 --- /dev/null +++ b/code/include/viewer_utils.h @@ -0,0 +1,11 @@ +#pragma once + +#include +#include + +namespace viewer_utils { + +// Deklaration (Implementation in viewer.cpp) +void simple_visualize(Eigen::MatrixXd& V, Eigen::MatrixXi& F); + +} \ No newline at end of file diff --git a/code/src/apps/v0/conformallab_cli.cpp b/code/src/apps/v0/conformallab_cli.cpp index d99cd14..c002a28 100644 --- a/code/src/apps/v0/conformallab_cli.cpp +++ b/code/src/apps/v0/conformallab_cli.cpp @@ -2,30 +2,39 @@ #include #include -#include + +#include "viewer_utils.h" +#include "mesh_utils.hpp" + #include #include #include -#include + #include using Kernel = CGAL::Simple_cartesian; using Point = Kernel::Point_3; using Mesh = CGAL::Surface_mesh; -namespace PMP = CGAL::Polygon_mesh_processing; -int main(int argc, char* argv[]) -{ + + + +bool parse_arguments(int argc, char* argv[], std::string& input_file, std::string& output_file, bool& show, bool& verbose) { CLI::App app{"demo of conformallab"}; - std::string input_file; - std::string output_file; // default output file name - app.add_option("-i,--input", input_file, "Input OFF file")->required(); app.add_option("-o,--output", output_file, "Output OFF file (optional)"); - CLI11_PARSE(app, argc, argv); + app.add_flag("-s,--show", show, "Show the input mesh in a viewer "); + app.add_flag("-v,--verbose",verbose, "Enable verbose output"); + app.set_help_flag("-h,--help", "Show this help message and exit"); + + try { + CLI11_PARSE(app, argc, argv); + } catch (const CLI::ParseError &e) { + return false; + } if (input_file.empty()) { std::cerr << "Input file is required.\n"; @@ -37,12 +46,38 @@ int main(int argc, char* argv[]) return EXIT_FAILURE; } + std::cout << "Input file: " << input_file << "\n"; + std::cout << "Output file: " << output_file << "\n"; + std::cout << "Show mesh: " << (show ? "Yes" : "No") << "\n"; + std::cout << "Verbose: " << (verbose ? "Yes" : "No") << "\n"; + + + + return true; +} + + +int main(int argc, char* argv[]) +{ + std::string input_file; + std::string output_file; + bool show = false; + bool verbose = false; + + if(!parse_arguments(argc, argv, input_file, output_file, show, verbose)) { + std::cerr << "Failed to parse arguments.\n"; + return EXIT_FAILURE; + } + Mesh surface_mesh; if (!CGAL::IO::read_polygon_mesh(input_file, surface_mesh) || surface_mesh.is_empty()) { std::cerr << "Invalid input file: " << input_file << "\n"; return EXIT_FAILURE; } + + + // #TODO: later: here we would call the unwrapping code, e.g.: // UnwrapSettings settings; // settings.target_geometry = TargetGeometry::Euclidean; @@ -52,36 +87,15 @@ int main(int argc, char* argv[]) // CGAL -> libigl - PMP::triangulate_faces(surface_mesh); - //todo : we should check if the mesh is already triangulated, and only call this if it isn't, to avoid unnecessary processing. CGAL::Polygon_mesh_processing::is_triangulated() can be used for this check. - Eigen::MatrixXd Vertices; - Eigen::MatrixXi Faces; - Vertices.resize(surface_mesh.number_of_vertices(), 3); - Faces.resize(surface_mesh.number_of_faces(), 3); - - for (auto v : surface_mesh.vertices()) { - const auto& p = surface_mesh.point(v); - Vertices(v.idx(), 0) = p.x(); - Vertices(v.idx(), 1) = p.y(); - Vertices(v.idx(), 2) = p.z(); - } - - int f_i = 0; - for (auto f : surface_mesh.faces()) { - int k = 0; - for (auto hv : CGAL::halfedges_around_face(surface_mesh.halfedge(f), surface_mesh)) { - auto v = surface_mesh.target(hv); - Faces(f_i, k) = v.idx(); - ++k; - } - ++f_i; - } - - igl::opengl::glfw::Viewer viewer; - viewer.data().set_mesh(Vertices, Faces); - viewer.launch(); - + if(show){ + std::cout << "Visualizing input mesh...\n"; + Eigen::MatrixXd V; + Eigen::MatrixXi F; + + mesh_utils::cgal_to_eigen(surface_mesh, V, F); + viewer_utils::simple_visualize(V, F); + } // for now, we just write the input mesh to the output file as a placeholder if (!output_file.empty() && !CGAL::IO::write_polygon_mesh(output_file, surface_mesh)) { diff --git a/code/src/main.cpp b/code/src/main.cpp deleted file mode 100644 index 5c14b7a..0000000 --- a/code/src/main.cpp +++ /dev/null @@ -1,40 +0,0 @@ -#include -#include - -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_2 Point_2; -typedef Kernel::Segment_2 Segment_2; - -int main() -{ - Point_2 p(1,1), q(10,10); - - std::cout << "p = " << p << std::endl; - std::cout << "q = " << q.x() << " " << q.y() << std::endl; - - std::cout << "sqdist(p,q) = " - << CGAL::squared_distance(p,q) << std::endl; - - Segment_2 s(p,q); - Point_2 m(5, 9); - - std::cout << "m = " << m << std::endl; - std::cout << "sqdist(Segment_2(p,q), m) = " - << CGAL::squared_distance(s,m) << std::endl; - - std::cout << "p, q, and m "; - switch (CGAL::orientation(p,q,m)){ - case CGAL::COLLINEAR: - std::cout << "are collinear\n"; - break; - case CGAL::LEFT_TURN: - std::cout << "make a left turn\n"; - break; - case CGAL::RIGHT_TURN: - std::cout << "make a right turn\n"; - break; - } - - std::cout << " midpoint(p,q) = " << CGAL::midpoint(p,q) << std::endl; - return 0; -} \ No newline at end of file diff --git a/code/src/viewer/simple_viewer.cpp b/code/src/viewer/simple_viewer.cpp new file mode 100644 index 0000000..c15d14b --- /dev/null +++ b/code/src/viewer/simple_viewer.cpp @@ -0,0 +1,9 @@ +#include "viewer_utils.h" + +namespace viewer_utils { + void simple_visualize(Eigen::MatrixXd& V, Eigen::MatrixXi& F) { + igl::opengl::glfw::Viewer viewer; + viewer.data().set_mesh(V, F); + viewer.launch(); + } +} diff --git a/code/tests/CMakeLists.txt b/code/tests/CMakeLists.txt new file mode 100644 index 0000000..d379268 --- /dev/null +++ b/code/tests/CMakeLists.txt @@ -0,0 +1,19 @@ +add_executable(conformallab_tests + test_clausen.cpp + test_hyper_ideal_utility.cpp + test_matrix_utility.cpp + test_surface_curve_utility.cpp +) + +target_include_directories(conformallab_tests SYSTEM PRIVATE + ${CMAKE_SOURCE_DIR}/deps/eigen-3.4.0 +) + +target_include_directories(conformallab_tests PRIVATE + ${CMAKE_SOURCE_DIR}/include +) + +target_link_libraries(conformallab_tests PRIVATE GTest::gtest_main) + +include(GoogleTest) +gtest_discover_tests(conformallab_tests DISCOVERY_TIMEOUT 60) diff --git a/code/tests/test_clausen.cpp b/code/tests/test_clausen.cpp new file mode 100644 index 0000000..2463910 --- /dev/null +++ b/code/tests/test_clausen.cpp @@ -0,0 +1,35 @@ +// Port of de.varylab.discreteconformal.functional.ClausenTest (Java/JUnit). +// Reference values computed with Mathematica. + +#include "clausen.hpp" + +#include +#include + +using conformallab::ImLi2; +using C = std::complex; + +TEST(ClausenTest, ImLi2_case1) { + double r = ImLi2(C(0.5, 0.5)); + EXPECT_NEAR(0.643767332889268748742017, r, 1e-15); +} + +TEST(ClausenTest, ImLi2_case2) { + double r = ImLi2(C(1.234, -3.554)); + EXPECT_NEAR(-2.784709023190200241929031, r, 1e-15); +} + +TEST(ClausenTest, ImLi2_case3) { + double r = ImLi2(C(-4.254, 2.965)); + EXPECT_NEAR(1.104092479314665907914598, r, 1e-15); +} + +TEST(ClausenTest, ImLi2_case4) { + double r = ImLi2(C(-2.264, -1.05)); + EXPECT_NEAR(-0.540683148040955780529547, r, 1e-15); +} + +TEST(ClausenTest, ImLi2_case5) { + double r = ImLi2(C(-4.0, 2.0)); + EXPECT_NEAR(0.7855694340809751306351005, r, 1e-15); +} diff --git a/code/tests/test_hyper_ideal_utility.cpp b/code/tests/test_hyper_ideal_utility.cpp new file mode 100644 index 0000000..71220fc --- /dev/null +++ b/code/tests/test_hyper_ideal_utility.cpp @@ -0,0 +1,92 @@ +// Port of de.varylab.discreteconformal.functional.HyperIdealUtilityTest (Java/JUnit). + +#include "hyper_ideal_utility.hpp" +#include "clausen.hpp" + +#include +#include + +using conformallab::calculateTetrahedronVolume; +using conformallab::calculateTetrahedronVolumeWithIdealVertexAtGamma; +using conformallab::Lobachevsky; + +constexpr double PI = 3.14159265358979323846264338328; + +// Regular tetrahedron at the Euclidean boundary (beta = arccos(1/3) for each +// vertex angle) has volume 0 — it degenerates to a flat configuration. +TEST(HyperIdealUtilityTest, VolumeEuclidean) { + double b = std::acos(1.0 / 3.0); + double V = calculateTetrahedronVolume(b, b, b, b, b, b); + EXPECT_NEAR(0.0, V, 1e-7); +} + +// Regular ideal tetrahedron with all angles pi/3. +// Formula: sum of Lobachevsky values at each angle. +TEST(HyperIdealUtilityTest, VolumeRegularIdeal1) { + double b = PI / 3.0; + double Ve = Lobachevsky(b) + Lobachevsky(b) + Lobachevsky(b); + double V = calculateTetrahedronVolume(b, b, b, b, b, b); + EXPECT_NEAR(Ve, V, 1e-12); +} + +// Right-angled ideal tetrahedron (pi/2, pi/4, pi/4). +TEST(HyperIdealUtilityTest, VolumeRegularIdeal2) { + double bi = PI / 2.0, bj = PI / 4.0, bk = PI / 4.0; + double Ve = Lobachevsky(bi) + Lobachevsky(bj) + Lobachevsky(bk); + double V = calculateTetrahedronVolume(bi, bj, bk, bi, bj, bk); + EXPECT_NEAR(Ve, V, 1e-12); +} + +// Hyperideal octahedron: all angles 0, volume = 8*Л(pi/4). +TEST(HyperIdealUtilityTest, VolumeOctahedron) { + double Ve = 8.0 * Lobachevsky(PI / 4.0); + double V = calculateTetrahedronVolume(0, 0, 0, 0, 0, 0); + EXPECT_NEAR(Ve, V, 1e-12); +} + +// Hyperideal tetrahedron with one hyperideal vertex. +// Manual formula from the paper vs. general formula. +TEST(HyperIdealUtilityTest, VolumeSingleHyperidealVertex) { + double bi = PI / 5.0, bj = PI / 4.0, bk = PI / 4.0; + double ai = (PI + bi - bj - bk) / 2.0; + double aj = (PI + bj - bi - bk) / 2.0; + double ak = (PI + bk - bi - bj) / 2.0; + double aijk= (PI - bk - bi - bj) / 2.0; + double Ve = 0.5 * (Lobachevsky(bi) + Lobachevsky(bj) + Lobachevsky(bk) + + Lobachevsky(ai) + Lobachevsky(aj) + Lobachevsky(ak) + + Lobachevsky(aijk)); + double V = calculateTetrahedronVolume(bi, bj, bk, ai, aj, ak); + EXPECT_NEAR(Ve, V, 1e-12); +} + +// A degenerate triangle (angle = pi) must give volume 0 without NaN. +TEST(HyperIdealUtilityTest, VolumeWithDegenerateTriangle) { + double V = calculateTetrahedronVolume(0.0, PI, 0.0, 0.0, 0.0, PI); + EXPECT_NEAR(0.0, V, 1e-12); + EXPECT_FALSE(std::isnan(V)); +} + +// The two volume formulas (general and ideal-vertex specialization) must agree +// on the same input — numerical consistency check. +TEST(HyperIdealUtilityTest, CompareGeneralAndIdealFormulaCase1) { + constexpr double EPS = 0.1; + double bi = PI / 3.0, bj = PI / 3.0, bk = PI / 3.0; + double ai = PI / 3.0 - EPS, aj = PI / 3.0 - EPS, ak = PI / 3.0 - EPS; + double Ve = calculateTetrahedronVolumeWithIdealVertexAtGamma(bi, bj, bk, ai, aj, ak); + double V = calculateTetrahedronVolume(bi, bj, bk, ai, aj, ak); + EXPECT_NEAR(Ve, V, 1e-12); +} + +// Second consistency check with non-symmetric angles that sum to pi. +TEST(HyperIdealUtilityTest, CompareGeneralAndIdealFormulaCase2) { + double bi = 0.6623267054958116; + double bj = 1.437248992086214; + double bk = 1.0420169560077686; + double ai = 0.6896178197389236; + double aj = 0.5195634857410114; + double ak = 0.6304500578493993; + EXPECT_NEAR(PI, bi + bj + bk, 1e-12); + double Ve = calculateTetrahedronVolumeWithIdealVertexAtGamma(bi, bj, bk, ai, aj, ak); + double V = calculateTetrahedronVolume(bi, bj, bk, ai, aj, ak); + EXPECT_NEAR(Ve, V, 1e-12); +} diff --git a/code/tests/test_matrix_utility.cpp b/code/tests/test_matrix_utility.cpp new file mode 100644 index 0000000..c300c8f --- /dev/null +++ b/code/tests/test_matrix_utility.cpp @@ -0,0 +1,40 @@ +// Port of de.varylab.discreteconformal.math.MatrixUtilityTest (Java/JUnit). + +#include "matrix_utility.hpp" + +#include +#include + +using Eigen::Matrix4d; +using Eigen::Vector4d; +using conformallab::makeMappingMatrix; + +TEST(MatrixUtilityTest, MakeMappingMatrix) { + // Source points (rows = homogeneous 4-vectors). + Matrix4d from; + from.row(0) = Vector4d(2, 0, 2, 1); + from.row(1) = Vector4d(1, 1, 0, 0); + from.row(2) = Vector4d(1, 0, 8, 0); + from.row(3) = Vector4d(3, 4, 0, 1); + + // Target points. + Matrix4d to; + to.row(0) = Vector4d(2, 0, 1, 0); + to.row(1) = Vector4d(0, 3, 0, 2); + to.row(2) = Vector4d(1, 0, 4, 0); + to.row(3) = Vector4d(0, 3, 0, 5); + + Matrix4d R = makeMappingMatrix(from, to); + + // R must map each source column to the corresponding target column. + for (int i = 0; i < 4; i++) { + Vector4d result = R * from.row(i).transpose(); + Vector4d expected = to.row(i).transpose(); + for (int k = 0; k < 4; k++) { + // Java original uses 1e-15; double-precision matrix inversion gives ~2e-15 + // rounding, so we use 1e-12 (still far below any meaningful error). + EXPECT_NEAR(expected(k), result(k), 1e-12) + << "Row " << i << ", component " << k; + } + } +} diff --git a/code/tests/test_surface_curve_utility.cpp b/code/tests/test_surface_curve_utility.cpp new file mode 100644 index 0000000..85db577 --- /dev/null +++ b/code/tests/test_surface_curve_utility.cpp @@ -0,0 +1,55 @@ +// Port of de.varylab.discreteconformal.uniformization.SurfaceCurveUtilityTest +// (Java/JUnit) — the two pure-math tests that don't need the HDS. + +#include "projective_math.hpp" + +#include +#include + +using Eigen::VectorXd; +using conformallab::isOnSegment; +using conformallab::getPointOnCorrespondingSegment; +using conformallab::dehomogenize; + +// Helper: make VectorXd from initializer list. +static VectorXd V(std::initializer_list vals) { + VectorXd v(vals.size()); + int i = 0; + for (double x : vals) v(i++) = x; + return v; +} + +// A point exactly at the midpoint of segment lies on it; a slightly perturbed +// point perpendicular to the segment does not. +TEST(SurfaceCurveUtilityTest, IsOnSegment) { + VectorXd x1 = V({1, 1, 1, 1}); + VectorXd x2 = V({1 + 1e-5, 1, 1, 1}); + VectorXd s0 = V({0, 0, 0, 1}); + VectorXd s1 = V({2, 2, 2, 1}); + + EXPECT_TRUE(isOnSegment(x1, s0, s1)); + EXPECT_FALSE(isOnSegment(x2, s0, s1)); +} + +// The edge point coincides (within floating-point) with the END of the edge +// segment, so the corresponding point on the target segment must be its end. +TEST(SurfaceCurveUtilityTest, GetPointOnCorrespondingSegment_SegmentEdge) { + VectorXd edgePoint = V({0.352392439203295, 0.9123804930829212, 1.0}); + VectorXd edgeSrc0 = V({0.34745306897719913, 0.912568467121888, 1.0}); + VectorXd edgeSrc1 = V({0.35239243920431296, 0.912380493082885, 1.0}); + VectorXd tgt0 = V({-0.2896352574166635, 0.03146361746587523, + 0.10643898885661185, 0.44373020886051084}); + VectorXd tgt1 = V({-0.2666822290964323, 0.019034256494171405, + 0.10525293907970201, 0.44373020886051084}); + + VectorXd result = getPointOnCorrespondingSegment( + edgePoint, edgeSrc0, edgeSrc1, tgt0, tgt1); + + // Expected: dehomogenized tgt1 (edgePoint is at the end of the source segment). + VectorXd expected = dehomogenize(tgt1); + + ASSERT_EQ(expected.size(), result.size()); + for (int i = 0; i < result.size(); i++) { + EXPECT_NEAR(expected(i), result(i), 1e-9) << "component " << i; + } +}