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 37d8089..94e0714 100644 --- a/code/CMakeLists.txt +++ b/code/CMakeLists.txt @@ -22,14 +22,31 @@ endif() 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 nur im Debug (nicht für Tests, da gtest_discover_tests + # die Binary zur Buildzeit ausführt und ASan dabei hängt) + 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() +# ----- Testing (Google Test via FetchContent) ----- +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) +# -------------------------------------------------- + add_subdirectory(deps) add_subdirectory(deps/glfw-3.4) @@ -38,7 +55,7 @@ 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) # --------- VIEWER Lib (separat → nur bei src/viewer/* Änderung gebaut!) --------- add_library(viewer STATIC src/viewer/simple_viewer.cpp) -target_include_directories(VIEWER PUBLIC +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 @@ -52,7 +69,7 @@ target_include_directories(${PROJECT_NAME} SYSTEM PRIVATE ${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 # Core libigl + ${CMAKE_CURRENT_SOURCE_DIR}/deps/libigl-2.6.0/include ) target_include_directories(${PROJECT_NAME} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/include) #-------- Compiler-Definitionen --------- @@ -61,5 +78,5 @@ target_link_libraries(${PROJECT_NAME} PRIVATE viewer) #--------- Output-Ordner --------- set_target_properties(${PROJECT_NAME} PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/bin) - +add_subdirectory(tests) 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/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/src/apps/v0/conformallab_cli.cpp b/code/src/apps/v0/conformallab_cli.cpp index e404966..c002a28 100644 --- a/code/src/apps/v0/conformallab_cli.cpp +++ b/code/src/apps/v0/conformallab_cli.cpp @@ -2,6 +2,7 @@ #include #include + #include "viewer_utils.h" #include "mesh_utils.hpp" @@ -34,7 +35,6 @@ bool parse_arguments(int argc, char* argv[], std::string& input_file, std::strin } catch (const CLI::ParseError &e) { return false; } - return true; if (input_file.empty()) { std::cerr << "Input file is required.\n"; @@ -45,12 +45,18 @@ bool parse_arguments(int argc, char* argv[], std::string& input_file, std::strin std::cerr << "Unsupported file format. Please provide an OFF file.\n"; 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; @@ -63,17 +69,15 @@ 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"; - 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; diff --git a/code/tests/CMakeLists.txt b/code/tests/CMakeLists.txt new file mode 100644 index 0000000..f8e25d7 --- /dev/null +++ b/code/tests/CMakeLists.txt @@ -0,0 +1,17 @@ +add_executable(conformallab_tests + test_clausen.cpp + test_hyper_ideal_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); +}