From e4632ce0ecbd609b6a624280bb1ac8378138e9a4 Mon Sep 17 00:00:00 2001 From: Markus Lanxinger Date: Sat, 21 Mar 2026 09:14:51 +0100 Subject: [PATCH 1/3] Vendor MeshSDFilter from Pete's fork, replacing AliceVision FetchContent Drop the AliceVision FetchContent + override mechanism in favor of directly vendored algorithm files from Pete's standalone MeshDenoiser fork. This unlocks CHOLMOD solver support, spatial hashing, deterministic mode, early-stop displacement, and full pipeline metrics (JSON/CSV). New features: - Flexible CLI: 2-arg mode with built-in defaults, --help, --write-default-options, --deterministic, --metrics-json/csv, --export-precision - Input validation rejects meshes with NaN/Inf coordinates - Pipeline timing printed to stdout on every run - RunStatistics with detailed solver/timing metrics Multi-format I/O preserved: OBJ/PLY/OFF/STL (OpenMesh), glTF (tinygltf v2.9.7), USD (tinyusdz v0.9.1 pinned with SHA256). Co-Authored-By: Claude Opus 4.6 (1M context) --- .github/workflows/build.yml | 14 +- CMakeLists.txt | 87 +- MeshDenoiserDefaults.txt | 25 +- README.md | 207 ++--- cmake/FetchDependencies.cmake | 103 +++ cmake/FetchMeshSDFilter.cmake | 104 --- src/AppMetrics.h | 131 +++ src/EigenTypes.h | 76 ++ .../MeshSDFilter => src}/MeshDenoiser.cpp | 817 +++++++---------- src/MeshNormalDenoising.h | 292 ++++++ src/MeshNormalFilter.h | 810 +++++++++++++++++ src/MeshTypes.h | 246 +++++ src/SDFilter.h | 848 ++++++++++++++++++ 13 files changed, 3018 insertions(+), 742 deletions(-) create mode 100644 cmake/FetchDependencies.cmake delete mode 100644 cmake/FetchMeshSDFilter.cmake create mode 100644 src/AppMetrics.h create mode 100644 src/EigenTypes.h rename {overrides/MeshSDFilter => src}/MeshDenoiser.cpp (52%) create mode 100644 src/MeshNormalDenoising.h create mode 100644 src/MeshNormalFilter.h create mode 100644 src/MeshTypes.h create mode 100644 src/SDFilter.h diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 2188119..b4baefc 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -68,11 +68,9 @@ jobs: if: runner.os != 'Windows' run: | mkdir -p artifacts - find build -name MeshSDFilter -type f -exec cp {} artifacts/ \; find build -name MeshDenoiser -type f -exec cp {} artifacts/ \; cp MeshDenoiserDefaults.txt artifacts/ cp README.md artifacts/ - chmod +x artifacts/MeshSDFilter chmod +x artifacts/MeshDenoiser # Bundle OpenMesh libraries and fix rpath @@ -80,22 +78,16 @@ jobs: mkdir -p artifacts/lib find build -name "libOpenMesh*.dylib" -exec cp {} artifacts/lib/ \; - # Fix rpath for both executables - for exe in artifacts/MeshSDFilter artifacts/MeshDenoiser; do + for exe in artifacts/MeshDenoiser; do if [ -f "$exe" ]; then - # Get list of OpenMesh libraries this executable depends on otool -L "$exe" | grep libOpenMesh | awk '{print $1}' | while read lib; do libname=$(basename "$lib") - # Change the rpath to look in the lib directory next to the executable install_name_tool -change "$lib" "@executable_path/lib/$libname" "$exe" done - - # Ad-hoc code sign to avoid some security warnings codesign --force --deep --sign - "$exe" fi done - # Also sign the libraries for lib in artifacts/lib/*.dylib; do if [ -f "$lib" ]; then codesign --force --sign - "$lib" @@ -105,8 +97,7 @@ jobs: mkdir -p artifacts/lib find build -name "libOpenMesh*.so*" -exec cp {} artifacts/lib/ \; - # Set rpath for both executables to look in lib directory - for exe in artifacts/MeshSDFilter artifacts/MeshDenoiser; do + for exe in artifacts/MeshDenoiser; do if [ -f "$exe" ]; then patchelf --set-rpath '$ORIGIN/lib' "$exe" || true fi @@ -117,7 +108,6 @@ jobs: if: runner.os == 'Windows' run: | mkdir artifacts - Get-ChildItem -Path build -Recurse -Filter MeshSDFilter.exe | ForEach-Object { Copy-Item $_.FullName -Destination artifacts\ } Get-ChildItem -Path build -Recurse -Filter MeshDenoiser.exe | ForEach-Object { Copy-Item $_.FullName -Destination artifacts\ } copy MeshDenoiserDefaults.txt artifacts\ copy README.md artifacts\ diff --git a/CMakeLists.txt b/CMakeLists.txt index 9605d17..7b3ebb0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,37 +1,94 @@ cmake_minimum_required(VERSION 3.18) -project(meshsdfilter-builder LANGUAGES CXX) +project(MeshDenoiser LANGUAGES CXX) -# Set C++17 standard for filesystem support set(CMAKE_CXX_STANDARD 17) -set(CMAKE_CXX_STANDARD_REQUIRED OFF) +set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) include(FetchContent) list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake") -# Option: enable/disable OpenMP (MeshSDFilter uses it if available) +# Compiler flags (applied per-target below, not globally, to avoid breaking fetched deps) + +# OpenMP (optional) option(ENABLE_OPENMP "Enable OpenMP if available" ON) +if(ENABLE_OPENMP) + find_package(OpenMP QUIET) +endif() + +# Eigen (required) +find_package(Eigen3 3.3 QUIET) +if(NOT TARGET Eigen3::Eigen) + message(STATUS "Eigen3 not found via find_package, trying EIGEN3_INCLUDE_DIR") + if(DEFINED EIGEN3_INCLUDE_DIR AND EXISTS "${EIGEN3_INCLUDE_DIR}") + add_library(Eigen3::Eigen INTERFACE IMPORTED) + set_target_properties(Eigen3::Eigen PROPERTIES + INTERFACE_INCLUDE_DIRECTORIES "${EIGEN3_INCLUDE_DIR}") + else() + message(FATAL_ERROR "Eigen3 not found. Install via system package manager or set -DEIGEN3_INCLUDE_DIR=...") + endif() +endif() -# --- Fetch the AliceVision source at the requested commit --- -set(ALICEVISION_COMMIT "14b0b8f8b3026765d165dfc3f219ed8c53635f52") -include(FetchMeshSDFilter) +# Fetch OpenMesh, tinygltf, tinyusdz, and optionally detect CHOLMOD +include(FetchDependencies) -# On MSVC, use static runtime if desired (optional) +# On MSVC, optionally link static runtime option(USE_STATIC_MSVC_RUNTIME "Link static MSVC runtime (/MT)" OFF) if(MSVC AND USE_STATIC_MSVC_RUNTIME) foreach(flag_var CMAKE_C_FLAGS CMAKE_C_FLAGS_DEBUG CMAKE_C_FLAGS_RELEASE - CMAKE_C_FLAGS_RELWITHDEBINFO CMAKE_C_FLAGS_MINSIZEREL - CMAKE_CXX_FLAGS CMAKE_CXX_FLAGS_DEBUG CMAKE_CXX_FLAGS_RELEASE - CMAKE_CXX_FLAGS_RELWITHDEBINFO CMAKE_CXX_FLAGS_MINSIZEREL) + CMAKE_C_FLAGS_RELWITHDEBINFO CMAKE_C_FLAGS_MINSIZEREL + CMAKE_CXX_FLAGS CMAKE_CXX_FLAGS_DEBUG CMAKE_CXX_FLAGS_RELEASE + CMAKE_CXX_FLAGS_RELWITHDEBINFO CMAKE_CXX_FLAGS_MINSIZEREL) if(${flag_var} MATCHES "/MD") string(REGEX REPLACE "/MD" "/MT" ${flag_var} "${${flag_var}}") endif() endforeach() endif() -if(ENABLE_OPENMP) - find_package(OpenMP) - if(OpenMP_CXX_FOUND) - message(STATUS "OpenMP found: enabling for targets built in MeshSDFilter") +# ---------- MeshDenoiser executable ---------- +add_executable(MeshDenoiser + src/EigenTypes.h + src/MeshTypes.h + src/SDFilter.h + src/MeshNormalFilter.h + src/MeshNormalDenoising.h + src/AppMetrics.h + src/MeshDenoiser.cpp +) + +# Per-target compiler flags +if(CMAKE_CXX_COMPILER_ID MATCHES "Clang|GNU|AppleClang") + target_compile_options(MeshDenoiser PRIVATE -Wall -Wextra) + target_compile_options(MeshDenoiser PRIVATE "$<$:-march=native;-mtune=native>") +elseif(MSVC) + target_compile_definitions(MeshDenoiser PRIVATE USE_MSVC _USE_MATH_DEFINES) +endif() + +target_include_directories(MeshDenoiser PRIVATE + "${CMAKE_CURRENT_SOURCE_DIR}/src" + "${tinygltf_SOURCE_DIR}" + "${tinyusdz_SOURCE_DIR}/src" +) + +target_link_libraries(MeshDenoiser PRIVATE + Eigen3::Eigen + OpenMeshCore + tinyusdz_static +) + +# OpenMP +if(OpenMP_CXX_FOUND) + target_compile_definitions(MeshDenoiser PRIVATE USE_OPENMP) + target_link_libraries(MeshDenoiser PRIVATE OpenMP::OpenMP_CXX) +endif() + +# CHOLMOD +if(SDFILTER_ENABLE_CHOLMOD) + target_compile_definitions(MeshDenoiser PRIVATE USE_CHOLMOD) + if(SDFILTER_CHOLMOD_INCLUDE_DIRS) + target_include_directories(MeshDenoiser PRIVATE ${SDFILTER_CHOLMOD_INCLUDE_DIRS}) + endif() + if(SDFILTER_CHOLMOD_LINK_LIBS) + target_link_libraries(MeshDenoiser PRIVATE ${SDFILTER_CHOLMOD_LINK_LIBS}) endif() endif() diff --git a/MeshDenoiserDefaults.txt b/MeshDenoiserDefaults.txt index 1f69ccb..b04da02 100644 --- a/MeshDenoiserDefaults.txt +++ b/MeshDenoiserDefaults.txt @@ -1,11 +1,34 @@ # Default options for MeshDenoiser tuned for detail-preserving cleanup. # -# Reduce OuterIterations or Lambda for even lighter touch; raise them (or Eta) if you need more smoothing. +# Reduce OuterIterations or Lambda for even lighter touch; +# raise them (or Eta) if you need more smoothing. +## Regularization weight. Higher = more smoothing per iteration. Lambda 0.15 + +## Spatial Gaussian sigma, scaled by average face-centroid distance. Eta 2.2 + +## Guidance normal difference weight. Mu 0.2 + +## Signal normal difference weight. Nu 0.25 + +## Closeness weight for mesh update (vertex position fidelity). MeshUpdateClosenessWeight 0.001 + +## Max iterations for mesh vertex update per outer iteration. MeshUpdateIterations 20 + +## Early-stop threshold for mesh update RMS displacement (<=0 disables). +MeshUpdateDisplacementEps 0.1 + +## Number of full filtering passes. OuterIterations 1 + +## Force single-threaded execution for reproducible output (0/1). +DeterministicMode 0 + +## Linear solver backend: 0=CG, 1=Eigen LDLT, 2=CHOLMOD (falls back to 1 if unavailable). +LinearSolverType 1 diff --git a/README.md b/README.md index 76190ca..92bf28f 100644 --- a/README.md +++ b/README.md @@ -1,163 +1,128 @@ -# meshdenoiser +# MeshDenoiser -A tiny, cross‑platform CMake repo that **uses the MeshSDFilter code from AliceVision** at commit -`14b0b8f8b3026765d165dfc3f219ed8c53635f52` and builds the two binaries: +A cross-platform mesh normal denoising tool based on the **Static/Dynamic Filtering** algorithm +([Zhang et al., arXiv:1712.03574](https://arxiv.org/abs/1712.03574)). -- `MeshSDFilter` -- `MeshDenoiser` +Smooths noisy 3D meshes while preserving geometric detail — useful for cleaning up +photogrammetry output, 3D scans, and other noisy mesh data. -Dependencies are handled per‑platform (OpenMesh sources are bundled by the upstream code; Eigen is installed by the workflow/package manager). +## Features -## Dependencies +- **Multi-format input:** OBJ, PLY, OFF, STL, glTF (.gltf/.glb), USD (.usd/.usda/.usdc/.usdz) +- **Multi-format output:** Any format supported by OpenMesh (OBJ, PLY, OFF, STL, etc.) +- **Flexible CLI:** Run with built-in defaults or a custom options file +- **Pipeline timing** printed to stdout; optional JSON/CSV metrics export +- **Deterministic mode** for reproducible results +- **Optional CHOLMOD** solver backend (auto-detected, falls back to Eigen LDLT) +- **OpenMP** parallelization (optional) +- **Input validation** — rejects meshes with NaN/Inf coordinates + +## Quick Start + +```bash +# Denoise with built-in defaults +MeshDenoiser input.obj output.obj + +# With a custom options file +MeshDenoiser options.txt input.obj output.obj -MeshSDFilter requires: -- **Eigen 3.3+** (header-only) -- **OpenMP** (optional, if compiler supports it) -- **OpenMesh 11.0** – fetched and built automatically from source -- **tinygltf** v2.9.6 (header-only, fetched automatically) to allow loading `.gltf` and `.glb` meshes when running `MeshDenoiser` -- **tinyusdz** (fetched automatically) to allow loading USD formats (`.usd`, `.usda`, `.usdc`, `.usdz`) when running `MeshDenoiser` +# Denoise a glTF or USD file +MeshDenoiser scan.glb denoised.obj +MeshDenoiser scene.usdz clean.ply -OpenMP is an open standard for shared-memory parallelism; compilers that support it (e.g. GCC, Clang with `libomp`, MSVC) let MeshSDFilter run heavy loops across multiple CPU cores. +# Generate a default options template +MeshDenoiser --write-default-options my_options.txt -> Note: MeshSDFilter expects Eigen to be discoverable via `find_package(Eigen3)`. Our CI installs Eigen per-platform so you don't have to. -> You can also point CMake at a local Eigen install with `-DEIGEN3_INCLUDE_DIR=/path/to/eigen` if needed. +# Show all options +MeshDenoiser --help +``` + +### Optional flags + +| Flag | Description | +|------|-------------| +| `--export-precision N` | Vertex coordinate precision for output (default: 16) | +| `--metrics-json PATH` | Write JSON timing and solver metrics | +| `--metrics-csv PATH` | Append CSV timing and solver metrics | +| `--deterministic` | Force single-threaded execution for reproducible output | +| `--write-default-options PATH` | Write the default options template to a file and exit | + +## Denoising Parameters + +The defaults are tuned for detail-preserving cleanup. Generate a commented template with +`MeshDenoiser --write-default-options options.txt`, or see `MeshDenoiserDefaults.txt`. + +| Parameter | Default | Description | +|-----------|---------|-------------| +| `Lambda` | 0.15 | Regularization weight. Higher = more smoothing per iteration | +| `Eta` | 2.2 | Spatial Gaussian sigma, scaled by average face-centroid distance | +| `Mu` | 0.2 | Guidance normal difference weight | +| `Nu` | 0.25 | Signal normal difference weight | +| `MeshUpdateClosenessWeight` | 0.001 | Vertex position fidelity during mesh update | +| `MeshUpdateIterations` | 20 | Max iterations for mesh vertex update per outer iteration | +| `MeshUpdateDisplacementEps` | 0.1 | Early-stop threshold for mesh update RMS displacement (<=0 disables) | +| `OuterIterations` | 1 | Number of full filtering passes. More = more smoothing | +| `DeterministicMode` | 0 | Force single-threaded execution (0/1) | +| `LinearSolverType` | 1 | 0=CG, 1=Eigen LDLT, 2=CHOLMOD (falls back to 1 if unavailable) | + +## Dependencies -## Build (local) +| Library | Version | Type | Notes | +|---------|---------|------|-------| +| **Eigen** | 3.3+ | Header-only | Must be findable via `find_package(Eigen3)` or `-DEIGEN3_INCLUDE_DIR=...` | +| **OpenMesh** | 11.0 | Fetched automatically | Mesh data structure and traditional format I/O | +| **tinygltf** | 2.9.7 | Fetched automatically | Header-only glTF 2.0 parser | +| **tinyusdz** | 0.9.1 | Fetched automatically | USD format support | +| **OpenMP** | — | Optional | Multi-threaded performance (auto-detected) | +| **SuiteSparse CHOLMOD** | — | Optional | Faster sparse solver (auto-detected, falls back to Eigen LDLT) | + +## Build ### Linux (Ubuntu/Debian) ```bash sudo apt-get update && sudo apt-get install -y build-essential cmake libeigen3-dev cmake -S . -B build -DCMAKE_BUILD_TYPE=Release cmake --build build -j -./build/MeshSDFilter --help ./build/MeshDenoiser --help ``` ### macOS -Using Homebrew: ```bash brew update && brew install cmake eigen cmake -S . -B build -DCMAKE_BUILD_TYPE=Release cmake --build build -j -./build/MeshSDFilter --help ./build/MeshDenoiser --help ``` -If you want OpenMP on macOS, install `gcc` and configure CMake with `-DCMAKE_C_COMPILER=gcc-14 -DCMAKE_CXX_COMPILER=g++-14` (or the version you installed). + +For OpenMP on macOS, install `gcc` and configure with `-DCMAKE_C_COMPILER=gcc-14 -DCMAKE_CXX_COMPILER=g++-14`. ### Windows (MSVC + vcpkg) ```powershell -# One-time: install vcpkg and integrate -# https://github.com/microsoft/vcpkg#quick-start-windows vcpkg install eigen3 - -# Configure with vcpkg toolchain so find_package(Eigen3) works cmake -S . -B build -DCMAKE_BUILD_TYPE=Release -DCMAKE_TOOLCHAIN_FILE="C:/path/to/vcpkg/scripts/buildsystems/vcpkg.cmake" cmake --build build --config Release -# Binaries will be under build/Release/ ``` +## Input Format Notes + +- For **glTF** and **USD** files with multiple meshes or transforms, all geometry is combined and transforms are applied automatically before filtering. +- Output format is determined by file extension (e.g., `.obj`, `.ply`, `.stl`). + ## CI / Pre-built Binaries -- GitHub Actions build and upload artifacts for Ubuntu, macOS, Windows. -- Artifacts include the two binaries per platform (`MeshSDFilter`, `MeshDenoiser`). -- Download the latest release binaries from the [Releases page](../../releases). -### macOS Security Note -When running downloaded binaries on macOS, you may see a security warning. To fix this: +GitHub Actions builds and uploads artifacts for Ubuntu, macOS, and Windows. +Download from the [Releases page](../../releases). -**Option 1 - Right-click method (Recommended):** -1. Right-click (or Control+click) on `MeshDenoiser` or `MeshSDFilter` -2. Select "Open" from the menu -3. Click "Open" in the security dialog -4. After doing this once, you can run the binary normally from Terminal +### macOS Security Note -**Option 2 - Command line:** +Downloaded binaries may be quarantined. To fix: ```bash -# Remove quarantine attribute from all files in the download xattr -cr meshdenoiser-macos/ ``` -**Option 3 - System Settings:** -1. Try to run the binary (it will be blocked) -2. Go to System Settings → Privacy & Security -3. Scroll down and click "Open Anyway" next to the blocked app message - ## Licensing -- **MeshSDFilter** code is BSD-3-Clause (see `LICENSES/MeshSDFilter-BSD-3-Clause.txt`). -- **OpenMesh** license is included as `LICENSES/OpenMesh-LICENSE.txt`. -- This wrapper repo is MIT by default (you can change it), and preserves upstream notices. - -## Using the tools - -### Single File Processing -```bash -# Filter -MeshSDFilter FilteringOptions.txt input_mesh.ply output_mesh.ply -# Denoise -MeshDenoiser DenoisingOptions.txt input_mesh.ply output_mesh.ply -``` - -### Batch Processing -Process multiple mesh files at once by providing directories instead of individual files: - -```bash -# Process all mesh files in input_dir/ and save to output_dir/ -MeshDenoiser DenoisingOptions.txt input_dir/ output_dir/ -``` - -Batch processing features: -- **Automatically discovers** all supported mesh files in the input directory -- **Preserves filenames** - each output file has the same name as its input -- **Creates output directory** if it doesn't exist -- **Progress indicators** - animated spinner showing loading, denoising, and saving progress -- **Real-time feedback** - displays mesh statistics (vertices, faces) and elapsed time -- **Error handling** - continues processing remaining files if one fails -- **Summary report** - displays success/failure counts and average processing time - -Example: -```bash -# Process a directory of scanned meshes -MeshDenoiser MeshDenoiserDefaults.txt scanned_meshes/ cleaned_meshes/ -``` - -Output example: -``` -Batch processing mode -Input directory: scanned_meshes/ -Output directory: cleaned_meshes/ - -Found 3 mesh file(s) to process - -[1/3] scan001.obj - Loaded mesh (12543 vertices, 25086 faces) (0.12s) - Complete (2.45s) - -[2/3] scan002.ply - Loaded mesh (8421 vertices, 16842 faces) (0.08s) - Complete (1.89s) - -[3/3] model.gltf - Loaded mesh (15632 vertices, 31264 faces) (0.15s) - Complete (3.12s) - -================================== -Batch processing complete - Successful: 3 - Failed: 0 - Total: 3 - Time: 7.56s (avg: 2.52s per file) -``` - -The progress indicator shows: -- **Loading phase**: Spinner animation with vertex/face count on completion -- **Denoising phase**: Spinner animation with elapsed time -- **Saving phase**: Quick spinner during file write -- Terminal animation works in interactive shells; falls back to simple dots in pipes/logs - -### Supported Formats -- A detail-preserving MeshDenoiser preset is in `MeshDenoiserDefaults.txt` (outer iterations 1, lambda 0.15, eta 2.2, mu 0.2, nu 0.25). Copy it to your working folder or pass it directly; raise lambda/eta or the iteration count only if you want stronger smoothing. -- `MeshDenoiser` accepts: - - Traditional formats: OBJ, PLY, OFF, STL (via OpenMesh) - - glTF formats: `.gltf`, `.glb` (via tinygltf) - - USD formats: `.usd`, `.usda`, `.usdc`, `.usdz` (via tinyusdz) - For glTF and USD files with multiple meshes or transforms, all geometry is combined and transforms are applied automatically before filtering. +- **MeshSDFilter algorithm** is BSD-3-Clause (see `LICENSES/MeshSDFilter-BSD-3-Clause.txt`). +- **OpenMesh** license: `LICENSES/OpenMesh-LICENSE.txt`. +- **tinygltf** license: `LICENSES/tinygltf-LICENSE.txt`. +- This project is licensed under the Mozilla Public License v2.0. diff --git a/cmake/FetchDependencies.cmake b/cmake/FetchDependencies.cmake new file mode 100644 index 0000000..4cabcd0 --- /dev/null +++ b/cmake/FetchDependencies.cmake @@ -0,0 +1,103 @@ +# Fetch external dependencies for MeshDenoiser. +include(FetchContent) + +if(POLICY CMP0169) + cmake_policy(SET CMP0169 OLD) +endif() + +if(NOT DEFINED CMAKE_POLICY_VERSION_MINIMUM OR CMAKE_POLICY_VERSION_MINIMUM VERSION_LESS 3.5) + set(CMAKE_POLICY_VERSION_MINIMUM 3.5) +endif() + +# ---------- OpenMesh ---------- +set(_OPENMESH_URL "https://www.graphics.rwth-aachen.de/media/openmesh_static/Releases/11.0/OpenMesh-11.0.0.tar.gz") +FetchContent_Declare( + openmesh + URL "${_OPENMESH_URL}" + DOWNLOAD_EXTRACT_TIMESTAMP TRUE +) +FetchContent_GetProperties(openmesh) +if(NOT openmesh_POPULATED) + FetchContent_Populate(openmesh) + set(BUILD_APPS OFF CACHE BOOL "" FORCE) + set(OPENMESH_DOCS OFF CACHE BOOL "" FORCE) + add_subdirectory("${openmesh_SOURCE_DIR}" "${openmesh_BINARY_DIR}") +endif() + +set(OpenMesh_DIR "${openmesh_BINARY_DIR}/src/cmake" CACHE PATH "" FORCE) +list(PREPEND CMAKE_PREFIX_PATH "${openmesh_BINARY_DIR}") +list(PREPEND CMAKE_MODULE_PATH "${openmesh_SOURCE_DIR}/cmake") + +# ---------- tinygltf (header-only) ---------- +set(_TINYGLTF_URL "https://github.com/syoyo/tinygltf/archive/refs/tags/v2.9.7.zip") +FetchContent_Declare( + tinygltf + URL "${_TINYGLTF_URL}" + DOWNLOAD_EXTRACT_TIMESTAMP TRUE + URL_HASH SHA256=1015f306721257fdcee602c2995542232042b4feda2ebb2e0c323c7d769ccd0e +) +FetchContent_GetProperties(tinygltf) +if(NOT tinygltf_POPULATED) + FetchContent_Populate(tinygltf) +endif() + +# ---------- tinyusdz ---------- +set(_TINYUSDZ_URL "https://github.com/lighttransport/tinyusdz/archive/refs/tags/v0.9.1.zip") + +set(TINYUSDZ_PRODUCTION_BUILD ON CACHE BOOL "" FORCE) +set(TINYUSDZ_WITH_OPENSUBDIV OFF CACHE BOOL "" FORCE) +set(TINYUSDZ_WITH_AUDIO OFF CACHE BOOL "" FORCE) +set(TINYUSDZ_WITH_EXR OFF CACHE BOOL "" FORCE) +set(TINYUSDZ_BUILD_TESTS OFF CACHE BOOL "" FORCE) +set(TINYUSDZ_BUILD_EXAMPLES OFF CACHE BOOL "" FORCE) +set(TINYUSDZ_NO_WERROR ON CACHE BOOL "" FORCE) + +FetchContent_Declare( + tinyusdz + URL "${_TINYUSDZ_URL}" + DOWNLOAD_EXTRACT_TIMESTAMP TRUE + URL_HASH SHA256=0cdbb15147783b916fb4a487a283495d271724c101df3c98303de96a3ef1af66 +) +FetchContent_MakeAvailable(tinyusdz) + +# ---------- Optional: SuiteSparse CHOLMOD ---------- +option(ENABLE_CHOLMOD "Enable SuiteSparse CHOLMOD backend when available" ON) +set(SDFILTER_ENABLE_CHOLMOD OFF) +set(SDFILTER_CHOLMOD_LINK_LIBS "") +set(SDFILTER_CHOLMOD_INCLUDE_DIRS "") + +if(ENABLE_CHOLMOD) + find_package(CHOLMOD CONFIG QUIET) + if(CHOLMOD_FOUND OR TARGET CHOLMOD::CHOLMOD) + set(SDFILTER_ENABLE_CHOLMOD ON) + message(STATUS "CHOLMOD found. Enabling CHOLMOD linear-solver backend.") + if(TARGET CHOLMOD::CHOLMOD) + set(SDFILTER_CHOLMOD_LINK_LIBS CHOLMOD::CHOLMOD) + elseif(DEFINED CHOLMOD_LIBRARIES) + set(SDFILTER_CHOLMOD_LINK_LIBS ${CHOLMOD_LIBRARIES}) + endif() + if(DEFINED CHOLMOD_INCLUDE_DIRS) + set(SDFILTER_CHOLMOD_INCLUDE_DIRS ${CHOLMOD_INCLUDE_DIRS}) + elseif(DEFINED CHOLMOD_INCLUDE_DIR) + set(SDFILTER_CHOLMOD_INCLUDE_DIRS ${CHOLMOD_INCLUDE_DIR}) + endif() + else() + find_package(SuiteSparse QUIET) + if(SuiteSparse_FOUND) + set(SDFILTER_ENABLE_CHOLMOD ON) + message(STATUS "SuiteSparse found. Enabling CHOLMOD linear-solver backend.") + if(TARGET SuiteSparse::CHOLMOD) + set(SDFILTER_CHOLMOD_LINK_LIBS SuiteSparse::CHOLMOD) + elseif(DEFINED CHOLMOD_LIBRARIES) + set(SDFILTER_CHOLMOD_LINK_LIBS ${CHOLMOD_LIBRARIES}) + endif() + if(DEFINED SuiteSparse_INCLUDE_DIRS) + set(SDFILTER_CHOLMOD_INCLUDE_DIRS ${SuiteSparse_INCLUDE_DIRS}) + elseif(DEFINED CHOLMOD_INCLUDE_DIRS) + set(SDFILTER_CHOLMOD_INCLUDE_DIRS ${CHOLMOD_INCLUDE_DIRS}) + endif() + else() + message(STATUS "SuiteSparse/CHOLMOD not found. CHOLMOD backend disabled; using Eigen LDLT.") + endif() + endif() +endif() diff --git a/cmake/FetchMeshSDFilter.cmake b/cmake/FetchMeshSDFilter.cmake deleted file mode 100644 index cc9a775..0000000 --- a/cmake/FetchMeshSDFilter.cmake +++ /dev/null @@ -1,104 +0,0 @@ -# Fetch the AliceVision tarball at a specific commit and add only MeshSDFilter. -include(FetchContent) - -if(POLICY CMP0169) - cmake_policy(SET CMP0169 OLD) -endif() - -if(NOT DEFINED CMAKE_POLICY_VERSION_MINIMUM OR CMAKE_POLICY_VERSION_MINIMUM VERSION_LESS 3.5) - set(CMAKE_POLICY_VERSION_MINIMUM 3.5) -endif() - -if(NOT ALICEVISION_COMMIT) - message(FATAL_ERROR "ALICEVISION_COMMIT must be set before including FetchMeshSDFilter.cmake") -endif() - -set(_AV_URL "https://github.com/alicevision/AliceVision/archive/${ALICEVISION_COMMIT}.zip") -FetchContent_Declare( - alicevision_src - URL "${_AV_URL}" - DOWNLOAD_EXTRACT_TIMESTAMP TRUE - URL_HASH SHA256=9f3d1718e5cff1fc81a21d2a498557c342913208d8479cf2565075cc0be71d9b -) -FetchContent_GetProperties(alicevision_src) -if(NOT alicevision_src_POPULATED) - FetchContent_Populate(alicevision_src) -endif() - -set(_OPENMESH_URL "https://www.graphics.rwth-aachen.de/media/openmesh_static/Releases/11.0/OpenMesh-11.0.0.tar.gz") -FetchContent_Declare( - openmesh - URL "${_OPENMESH_URL}" - DOWNLOAD_EXTRACT_TIMESTAMP TRUE -) -FetchContent_GetProperties(openmesh) -if(NOT openmesh_POPULATED) - FetchContent_Populate(openmesh) - set(BUILD_APPS OFF CACHE BOOL "" FORCE) - set(OPENMESH_DOCS OFF CACHE BOOL "" FORCE) - add_subdirectory("${openmesh_SOURCE_DIR}" "${openmesh_BINARY_DIR}") -endif() - -set(OpenMesh_DIR "${openmesh_BINARY_DIR}/src/cmake" CACHE PATH "" FORCE) -list(PREPEND CMAKE_PREFIX_PATH "${openmesh_BINARY_DIR}") -list(PREPEND CMAKE_MODULE_PATH "${openmesh_SOURCE_DIR}/cmake") - -set(_TINYGLTF_URL "https://github.com/syoyo/tinygltf/archive/refs/tags/v2.9.6.zip") -FetchContent_Declare( - tinygltf - URL "${_TINYGLTF_URL}" - DOWNLOAD_EXTRACT_TIMESTAMP TRUE - URL_HASH SHA256=a3eade1a5f00a756f81d43e1ef6b23c3bde9cda0c52a7a69cca5dd6ad8f1dd95 -) -FetchContent_GetProperties(tinygltf) -if(NOT tinygltf_POPULATED) - FetchContent_Populate(tinygltf) -endif() - -set(_TINYUSDZ_URL "https://github.com/lighttransport/tinyusdz/archive/refs/heads/release.zip") - -# TinyUSDZ build options must be set BEFORE FetchContent -set(TINYUSDZ_PRODUCTION_BUILD ON CACHE BOOL "" FORCE) -set(TINYUSDZ_WITH_OPENSUBDIV OFF CACHE BOOL "" FORCE) -set(TINYUSDZ_WITH_AUDIO OFF CACHE BOOL "" FORCE) -set(TINYUSDZ_WITH_EXR OFF CACHE BOOL "" FORCE) -set(TINYUSDZ_BUILD_TESTS OFF CACHE BOOL "" FORCE) -set(TINYUSDZ_BUILD_EXAMPLES OFF CACHE BOOL "" FORCE) - -FetchContent_Declare( - tinyusdz - URL "${_TINYUSDZ_URL}" - DOWNLOAD_EXTRACT_TIMESTAMP TRUE -) -FetchContent_MakeAvailable(tinyusdz) - -set(MESHSD_DIR "${alicevision_src_SOURCE_DIR}/src/dependencies/MeshSDFilter") - -if(NOT EXISTS "${MESHSD_DIR}/CMakeLists.txt") - message(FATAL_ERROR "MeshSDFilter CMakeLists.txt not found at ${MESHSD_DIR}. Commit may have changed structure.") -endif() - -set(_override_dir "${PROJECT_SOURCE_DIR}/overrides/MeshSDFilter") -message(STATUS "Checking override dir: ${_override_dir}") -if(EXISTS "${_override_dir}") - file(GLOB _override_sources "${_override_dir}/*") - foreach(_src IN LISTS _override_sources) - get_filename_component(_fname "${_src}" NAME) - message(STATUS "Overriding MeshSDFilter source: ${_fname}") - configure_file("${_src}" "${MESHSD_DIR}/${_fname}" COPYONLY) - endforeach() -endif() - -# Add include directories globally before building MeshSDFilter -# This ensures tinygltf and tinyusdz headers are available during compilation -include_directories("${tinygltf_SOURCE_DIR}") -include_directories("${tinyusdz_SOURCE_DIR}/src") - -add_subdirectory("${MESHSD_DIR}" "${CMAKE_CURRENT_BINARY_DIR}/MeshSDFilter-build") - -# Link tinyusdz to the targets after they're created -foreach(_mesh_target MeshSDFilter MeshDenoiser) - if(TARGET ${_mesh_target}) - target_link_libraries(${_mesh_target} tinyusdz_static) - endif() -endforeach() diff --git a/src/AppMetrics.h b/src/AppMetrics.h new file mode 100644 index 0000000..05dbf2a --- /dev/null +++ b/src/AppMetrics.h @@ -0,0 +1,131 @@ +#ifndef APPMETRICS_H +#define APPMETRICS_H + +#include "SDFilter.h" + +#include +#include + +namespace SDFilter +{ + +inline std::string json_escape(const std::string &s) +{ + std::string out; + out.reserve(s.size()); + for(char c : s) + { + switch(c) + { + case '\\': out += "\\\\"; break; + case '\"': out += "\\\""; break; + case '\n': out += "\\n"; break; + case '\r': out += "\\r"; break; + case '\t': out += "\\t"; break; + default: out.push_back(c); break; + } + } + return out; +} + +struct PipelineMetrics +{ + std::string mode; + std::string input_mesh; + std::string output_mesh; + double import_secs = 0.0; + double normalize_secs = 0.0; + double algorithm_secs = 0.0; + double restore_secs = 0.0; + double export_secs = 0.0; + double total_secs = 0.0; + int obj_export_precision = 16; +}; + +inline bool write_metrics_json(const std::string &path, const PipelineMetrics &m, const RunStatistics &s) +{ + std::ofstream out(path.c_str()); + if(!out.is_open()){ + return false; + } + + out << "{\n"; + out << " \"mode\": \"" << json_escape(m.mode) << "\",\n"; + out << " \"input_mesh\": \"" << json_escape(m.input_mesh) << "\",\n"; + out << " \"output_mesh\": \"" << json_escape(m.output_mesh) << "\",\n"; + out << " \"timing\": {\n"; + out << " \"import_secs\": " << m.import_secs << ",\n"; + out << " \"normalize_secs\": " << m.normalize_secs << ",\n"; + out << " \"algorithm_secs\": " << m.algorithm_secs << ",\n"; + out << " \"restore_secs\": " << m.restore_secs << ",\n"; + out << " \"export_secs\": " << m.export_secs << ",\n"; + out << " \"total_secs\": " << m.total_secs << ",\n"; + out << " \"preprocessing_secs\": " << s.preprocessing_secs << ",\n"; + out << " \"filtering_secs\": " << s.filtering_secs << ",\n"; + out << " \"mesh_update_secs\": " << s.mesh_update_secs << ",\n"; + out << " \"mesh_filter_total_secs\": " << s.mesh_filter_total_secs << ",\n"; + out << " \"denoise_total_secs\": " << s.denoise_total_secs << "\n"; + out << " },\n"; + out << " \"solver\": {\n"; + out << " \"iterations\": " << s.solver_iterations << ",\n"; + out << " \"converged\": " << (s.solver_converged ? "true" : "false") << ",\n"; + out << " \"outer_iterations\": " << s.outer_iterations << "\n"; + out << " },\n"; + out << " \"io\": {\n"; + out << " \"obj_export_precision\": " << m.obj_export_precision << "\n"; + out << " }\n"; + out << "}\n"; + + return out.good(); +} + +inline bool write_metrics_csv( + const std::string &path, + const PipelineMetrics &m, + const RunStatistics &s, + bool write_header) +{ + std::ofstream out; + if(write_header){ + out.open(path.c_str(), std::ios::out | std::ios::trunc); + } + else{ + out.open(path.c_str(), std::ios::out | std::ios::app); + } + + if(!out.is_open()){ + return false; + } + + if(write_header){ + out << "mode,input_mesh,output_mesh,import_secs,normalize_secs,algorithm_secs,restore_secs,export_secs,total_secs," + << "preprocessing_secs,filtering_secs,mesh_update_secs,mesh_filter_total_secs,denoise_total_secs," + << "solver_iterations,solver_converged,outer_iterations," + << "obj_export_precision\n"; + } + + out << m.mode << "," + << m.input_mesh << "," + << m.output_mesh << "," + << m.import_secs << "," + << m.normalize_secs << "," + << m.algorithm_secs << "," + << m.restore_secs << "," + << m.export_secs << "," + << m.total_secs << "," + << s.preprocessing_secs << "," + << s.filtering_secs << "," + << s.mesh_update_secs << "," + << s.mesh_filter_total_secs << "," + << s.denoise_total_secs << "," + << s.solver_iterations << "," + << (s.solver_converged ? 1 : 0) << "," + << s.outer_iterations << "," + << m.obj_export_precision << "\n"; + + return out.good(); +} + +} + +#endif // APPMETRICS_H diff --git a/src/EigenTypes.h b/src/EigenTypes.h new file mode 100644 index 0000000..c449120 --- /dev/null +++ b/src/EigenTypes.h @@ -0,0 +1,76 @@ +// BSD 3-Clause License +// +// Copyright (c) 2017, Bailin Deng +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// * Redistributions of source code must retain the above copyright notice, this +// list of conditions and the following disclaimer. +// +// * Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// * Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +// OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +#ifndef EIGENTYPES_H +#define EIGENTYPES_H + +#include +#include +#include + + +namespace SDFilter +{ + +// Define eigen matrix types +typedef Eigen::Matrix Matrix3X; +typedef Eigen::Matrix Matrix2X; +typedef Eigen::Matrix Matrix2Xi; +typedef Eigen::Matrix MatrixX3; +typedef Eigen::Matrix MatrixX2; +typedef Eigen::Matrix Matrix3Xi; +typedef Eigen::Matrix Matrix2XIdx; +typedef Eigen::Matrix VectorXIdx; +typedef Eigen::SparseMatrix SparseMatrixXd; +typedef Eigen::Triplet Triplet; + +// Conversion between a 3d vector type to Eigen::Vector3d +template +inline Eigen::Vector3d to_eigen_vec3d(const Vec_T &vec) +{ + return Eigen::Vector3d(vec[0], vec[1], vec[2]); +} + + +template +inline Vec_T from_eigen_vec3d(const Eigen::Vector3d &vec) +{ + Vec_T v; + v[0] = vec(0); + v[1] = vec(1); + v[2] = vec(2); + + return v; +} + +} + + +#endif // EIGENTYPES_H diff --git a/overrides/MeshSDFilter/MeshDenoiser.cpp b/src/MeshDenoiser.cpp similarity index 52% rename from overrides/MeshSDFilter/MeshDenoiser.cpp rename to src/MeshDenoiser.cpp index 01d21f7..ce569c6 100644 --- a/overrides/MeshSDFilter/MeshDenoiser.cpp +++ b/src/MeshDenoiser.cpp @@ -31,56 +31,26 @@ #include "MeshTypes.h" #include "MeshNormalDenoising.h" +#include "AppMetrics.h" #include -#include #include #include #include +#include +#include #include -#include +#include #include +#include #include #include #include #include -#include -#include -#include - -#if defined(_WIN32) || defined(_WIN64) - #include - #define ISATTY _isatty - #define FILENO _fileno -#else - #include - #define ISATTY isatty - #define FILENO fileno -#endif - -// Filesystem support -// Note: MSVC's std::experimental::filesystem has incomplete API, so we use native Windows API for MSVC -#if __cplusplus >= 201703L && !defined(_MSC_VER) - #include - namespace fs = std::filesystem; - #define HAS_FILESYSTEM -#elif defined(__has_include) && !defined(_MSC_VER) - #if __has_include() - #include - namespace fs = std::filesystem; - #define HAS_FILESYSTEM - #endif -#endif -#ifndef HAS_FILESYSTEM - #include - #include - #if defined(_WIN32) || defined(_WIN64) - #include - #else - #include - #endif +#ifdef USE_OPENMP +#include #endif #define TINYGLTF_IMPLEMENTATION @@ -93,113 +63,163 @@ #include "tydra/render-data.hh" #include "tydra/scene-access.hh" -// Simple progress indicator for terminal output -class ProgressIndicator +namespace { -public: - ProgressIndicator(const std::string &message, bool show_spinner = true) - : message_(message) - , show_spinner_(show_spinner && ISATTY(FILENO(stdout))) - , start_time_(std::chrono::steady_clock::now()) - , spinner_index_(0) - , active_(true) - { - if(show_spinner_){ - std::cout << message_ << " " << spinner_chars_[0] << std::flush; - } - else{ - std::cout << message_ << "..." << std::flush; - } + +using namespace tinyusdz; // For tydra namespace access + +// --------------------------------------------------------------------------- +// Utility helpers +// --------------------------------------------------------------------------- + +std::string to_lower_copy(std::string str) +{ + std::transform(str.begin(), str.end(), str.begin(), [](unsigned char c) { return static_cast(std::tolower(c)); }); + return str; +} + +bool has_extension(const std::string &path, const std::string &ext) +{ + if(path.size() < ext.size()){ + return false; } + return to_lower_copy(path.substr(path.size() - ext.size())) == ext; +} - ~ProgressIndicator() - { - finish(); +// --------------------------------------------------------------------------- +// CLI helpers +// --------------------------------------------------------------------------- + +void configure_deterministic_runtime(bool deterministic) +{ + if(!deterministic){ + return; } +#ifdef USE_OPENMP + omp_set_dynamic(0); + omp_set_num_threads(1); +#endif + Eigen::setNbThreads(1); +} - void update() - { - if(!active_ || !show_spinner_){ - return; - } +const char* default_options_text() +{ + return +R"(#### Option file for SD filter based mesh denoising +#### Lines starting with '#' are comments - spinner_index_ = (spinner_index_ + 1) % spinner_chars_.size(); - auto now = std::chrono::steady_clock::now(); - auto elapsed = std::chrono::duration_cast(now - start_time_).count(); +## Regularization weight, must be positive. +## Higher values = more smoothing per iteration. +Lambda 0.15 - // Clear current line and redraw - std::cout << "\r" << message_ << " " << spinner_chars_[spinner_index_]; - if(elapsed > 0){ - std::cout << " (" << elapsed << "s)"; - } - std::cout << std::flush; - } +## Gaussian standard deviation for spatial weight, scaled by the average +## distance between adjacent face centroids. Must be positive. +## Higher values = wider spatial neighborhood. +Eta 2.2 - void finish(const std::string &final_message = "") - { - if(!active_){ - return; - } +## Gaussian standard deviation for guidance weight, must be positive. +Mu 0.2 - active_ = false; +## Gaussian standard deviation for signal weight, must be positive. +Nu 0.25 - auto now = std::chrono::steady_clock::now(); - auto elapsed = std::chrono::duration_cast(now - start_time_).count(); +## Closeness weight for mesh update, must be positive. +MeshUpdateClosenessWeight 0.001 - if(show_spinner_){ - std::cout << "\r"; - // Clear the line - std::cout << std::string(message_.length() + 20, ' ') << "\r"; - } +## Iterations for mesh update, must be positive. +MeshUpdateIterations 20 - if(!final_message.empty()){ - std::cout << final_message; - } - else{ - std::cout << message_ << " - done"; - } +## Early-stop threshold for per-iteration mesh update RMS displacement (<=0 disables). +MeshUpdateDisplacementEps 1e-1 - // Show elapsed time if significant - if(elapsed >= 100){ - std::cout << " (" << std::fixed << std::setprecision(2) << (elapsed / 1000.0) << "s)"; - } +## Outer iterations for denoising, must be positive integer. +## More iterations = more smoothing overall. +OuterIterations 1 - std::cout << std::endl; - } +## Force deterministic mode (single-threaded runtime) for reproducible outputs (0/1). +DeterministicMode 0 - void set_message(const std::string &message) - { - message_ = message; +## Linear solver backend: 0=CG, 1=Eigen LDLT, 2=CHOLMOD (falls back to 1 if unavailable). +LinearSolverType 1 +)"; +} + +bool write_default_options_file(const std::string &path) +{ + std::ofstream out(path.c_str(), std::ios::out | std::ios::trunc); + if(!out.is_open()){ + return false; } -private: - std::string message_; - bool show_spinner_; - std::chrono::steady_clock::time_point start_time_; - size_t spinner_index_; - bool active_; - const std::vector spinner_chars_ = {'|', '/', '-', '\\'}; -}; + out << default_options_text(); + return out.good(); +} -namespace +void apply_default_options(SDFilter::MeshDenoisingParameters ¶m) { + param.lambda = 0.15; + param.eta = 2.2; + param.mu = 0.2; + param.nu = 0.25; + param.mesh_update_closeness_weight = 0.001; + param.mesh_update_iter = 20; + param.mesh_update_disp_eps = 1e-1; + param.outer_iterations = 1; + param.deterministic_mode = false; + param.linear_solver_type = static_cast(SDFilter::DEFAULT_LINEAR_SOLVER_TYPE); +} -using namespace tinyusdz; // For tydra namespace access - -std::string to_lower_copy(std::string str) +void print_help(const char *exe_name) { - std::transform(str.begin(), str.end(), str.begin(), [](unsigned char c) { return static_cast(std::tolower(c)); }); - return str; + std::cout << "MeshDenoiser - Mesh normal denoising filter" << std::endl; + std::cout << std::endl; + std::cout << "Usage:" << std::endl; + std::cout << " " << exe_name << " INPUT_MESH OUTPUT_MESH [optional flags]" << std::endl; + std::cout << " " << exe_name << " OPTION_FILE INPUT_MESH OUTPUT_MESH [optional flags]" << std::endl; + std::cout << " " << exe_name << " --write-default-options PATH" << std::endl; + std::cout << " " << exe_name << " --help" << std::endl; + std::cout << std::endl; + std::cout << "Positional arguments:" << std::endl; + std::cout << " OPTION_FILE Optional. Path to denoising options text file." << std::endl; + std::cout << " If omitted, built-in defaults are used." << std::endl; + std::cout << " INPUT_MESH Required. Input mesh to denoise." << std::endl; + std::cout << " OUTPUT_MESH Required. Output mesh path." << std::endl; + std::cout << std::endl; + std::cout << "Supported input formats:" << std::endl; + std::cout << " OBJ, PLY, OFF, STL (via OpenMesh)" << std::endl; + std::cout << " glTF (.gltf, .glb)" << std::endl; + std::cout << " USD (.usd, .usda, .usdc, .usdz)" << std::endl; + std::cout << std::endl; + std::cout << "Optional flags:" << std::endl; + std::cout << " --export-precision N Vertex coordinate precision for output. Default: 16." << std::endl; + std::cout << " --metrics-json PATH Write JSON timing and solver metrics." << std::endl; + std::cout << " --metrics-csv PATH Append CSV timing and solver metrics." << std::endl; + std::cout << " --deterministic Force deterministic single-threaded runtime." << std::endl; + std::cout << " --write-default-options PATH" << std::endl; + std::cout << " Write the default options template to PATH and exit." << std::endl; + std::cout << " --help Show this help text and exit." << std::endl; } -bool has_extension(const std::string &path, const std::string &ext) +// --------------------------------------------------------------------------- +// Input validation +// --------------------------------------------------------------------------- + +bool validate_mesh_vertices(const TriMesh &mesh) { - if(path.size() < ext.size()){ - return false; + for(auto v_it = mesh.vertices_begin(); v_it != mesh.vertices_end(); ++v_it){ + const auto &pt = mesh.point(*v_it); + if(!std::isfinite(pt[0]) || !std::isfinite(pt[1]) || !std::isfinite(pt[2])){ + std::cerr << "Error: non-finite vertex coordinate at index " << v_it->idx() << std::endl; + return false; + } } - return to_lower_copy(path.substr(path.size() - ext.size())) == ext; + return true; } +// --------------------------------------------------------------------------- +// glTF loader +// --------------------------------------------------------------------------- + struct PrimitiveInstance { const tinygltf::Primitive *primitive; @@ -266,11 +286,6 @@ void collect_primitives(const tinygltf::Model &model, } } -Eigen::Matrix4d default_transform() -{ - return Eigen::Matrix4d::Identity(); -} - template const unsigned char *accessor_element_ptr(const tinygltf::Model &model, const tinygltf::Accessor &accessor, @@ -372,7 +387,7 @@ bool load_gltf_mesh(const std::string &filename, TriMesh &mesh, std::string &war } std::vector instances; - Eigen::Matrix4d identity = default_transform(); + Eigen::Matrix4d identity = Eigen::Matrix4d::Identity(); if(model.scenes.empty()){ for(size_t i = 0; i < model.nodes.size(); ++i){ @@ -487,9 +502,7 @@ bool load_gltf_mesh(const std::string &filename, TriMesh &mesh, std::string &war Eigen::Vector4d p(position.x(), position.y(), position.z(), 1.0); Eigen::Vector4d transformed = instance.transform * p; - TriMesh::Point point(static_cast(transformed.x()), - static_cast(transformed.y()), - static_cast(transformed.z())); + TriMesh::Point point(transformed.x(), transformed.y(), transformed.z()); TriMesh::VertexHandle vh = mesh.add_vertex(point); cache_it = vertex_cache.emplace(key, vh).first; @@ -519,12 +532,15 @@ bool load_gltf_mesh(const std::string &filename, TriMesh &mesh, std::string &war return true; } +// --------------------------------------------------------------------------- +// USD loader +// --------------------------------------------------------------------------- + bool load_usd_mesh(const std::string &filename, TriMesh &mesh, std::string &warning, std::string &error) { tinyusdz::Stage stage; std::string warn, err; - // Load USD file based on extension bool success = false; if(has_extension(filename, ".usdz")){ success = tinyusdz::LoadUSDZFromFile(filename, &stage, &warn, &err); @@ -536,7 +552,6 @@ bool load_usd_mesh(const std::string &filename, TriMesh &mesh, std::string &warn success = tinyusdz::LoadUSDAFromFile(filename, &stage, &warn, &err); } else if(has_extension(filename, ".usd")){ - // Generic .usd extension - auto-detect format (ASCII or binary) success = tinyusdz::LoadUSDFromFile(filename, &stage, &warn, &err); } else{ @@ -557,7 +572,6 @@ bool load_usd_mesh(const std::string &filename, TriMesh &mesh, std::string &warn return false; } - // Convert USD stage to RenderScene using Tydra tydra::RenderScene render_scene; tydra::RenderSceneConverter converter; tydra::RenderSceneConverterEnv env(stage); @@ -572,18 +586,15 @@ bool load_usd_mesh(const std::string &filename, TriMesh &mesh, std::string &warn return false; } - // Clear and prepare the output mesh mesh.clear(); mesh.request_face_normals(); mesh.request_vertex_normals(); - // Process all meshes in the render scene for(const auto &render_mesh : render_scene.meshes){ if(render_mesh.points.empty()){ continue; } - // Add vertices std::vector vertex_handles; vertex_handles.reserve(render_mesh.points.size()); @@ -592,17 +603,14 @@ bool load_usd_mesh(const std::string &filename, TriMesh &mesh, std::string &warn vertex_handles.push_back(mesh.add_vertex(p)); } - // Determine which indices to use const std::vector *face_indices = nullptr; const std::vector *face_counts = nullptr; if(!render_mesh.triangulatedFaceVertexIndices.empty()){ - // Use triangulated indices if available face_indices = &render_mesh.triangulatedFaceVertexIndices; face_counts = &render_mesh.triangulatedFaceVertexCounts; } else{ - // Use original indices face_indices = &render_mesh.usdFaceVertexIndices; face_counts = &render_mesh.usdFaceVertexCounts; } @@ -611,11 +619,9 @@ bool load_usd_mesh(const std::string &filename, TriMesh &mesh, std::string &warn continue; } - // Add faces size_t index_offset = 0; for(uint32_t count : *face_counts){ if(count == 3){ - // Triangle - add directly std::vector face_verts; face_verts.reserve(3); @@ -634,7 +640,6 @@ bool load_usd_mesh(const std::string &filename, TriMesh &mesh, std::string &warn } } else if(count > 3){ - // Polygon - triangulate by fan method std::vector poly_indices; poly_indices.reserve(count); @@ -649,7 +654,6 @@ bool load_usd_mesh(const std::string &filename, TriMesh &mesh, std::string &warn } if(valid){ - // Simple fan triangulation from first vertex for(uint32_t i = 1; i + 1 < count; ++i){ std::vector face_verts; face_verts.push_back(vertex_handles[poly_indices[0]]); @@ -673,415 +677,250 @@ bool load_usd_mesh(const std::string &filename, TriMesh &mesh, std::string &warn return true; } -// Directory and file utilities -bool is_directory(const std::string &path) -{ -#ifdef HAS_FILESYSTEM - return fs::is_directory(path); -#else - #if defined(_WIN32) || defined(_WIN64) - DWORD attrib = GetFileAttributesA(path.c_str()); - return (attrib != INVALID_FILE_ATTRIBUTES && (attrib & FILE_ATTRIBUTE_DIRECTORY)); - #else - struct stat st; - if(stat(path.c_str(), &st) != 0){ - return false; - } - return S_ISDIR(st.st_mode); - #endif -#endif -} +// --------------------------------------------------------------------------- +// Multi-format mesh loading +// --------------------------------------------------------------------------- -bool is_regular_file(const std::string &path) +bool load_mesh(const std::string &filename, TriMesh &mesh) { -#ifdef HAS_FILESYSTEM - return fs::is_regular_file(path); -#else - #if defined(_WIN32) || defined(_WIN64) - DWORD attrib = GetFileAttributesA(path.c_str()); - if(attrib == INVALID_FILE_ATTRIBUTES){ - return false; - } - return !(attrib & FILE_ATTRIBUTE_DIRECTORY); - #else - struct stat st; - if(stat(path.c_str(), &st) != 0){ - return false; - } - return S_ISREG(st.st_mode); - #endif -#endif -} - -bool is_supported_mesh_extension(const std::string &path) -{ - // List of supported extensions - static const std::vector supported = { - ".obj", ".ply", ".off", ".stl", // OpenMesh formats - ".gltf", ".glb", // glTF formats - ".usd", ".usda", ".usdc", ".usdz" // USD formats - }; + std::string warning; + std::string error; + bool success = false; - for(const auto &ext : supported){ - if(has_extension(path, ext)){ - return true; + if(has_extension(filename, ".gltf") || has_extension(filename, ".glb")){ + success = load_gltf_mesh(filename, mesh, warning, error); + } + else if(has_extension(filename, ".usd") || has_extension(filename, ".usda") || + has_extension(filename, ".usdc") || has_extension(filename, ".usdz")){ + success = load_usd_mesh(filename, mesh, warning, error); + } + else{ + success = OpenMesh::IO::read_mesh(mesh, filename); + if(!success){ + error = "unable to read input mesh with OpenMesh"; } } - return false; -} -std::vector list_mesh_files(const std::string &directory) -{ - std::vector files; - -#ifdef HAS_FILESYSTEM - try{ - for(const auto &entry : fs::directory_iterator(directory)){ - if(entry.is_regular_file()){ - std::string filename = entry.path().filename().string(); - if(is_supported_mesh_extension(filename)){ - files.push_back(filename); - } - } - } + if(!warning.empty()){ + std::cout << "Warning while loading mesh: " << warning << std::endl; } - catch(const fs::filesystem_error &e){ - std::cerr << "Error reading directory: " << e.what() << std::endl; + + if(!success){ + std::cerr << "Error: " << error << std::endl; } -#else - #if defined(_WIN32) || defined(_WIN64) - WIN32_FIND_DATAA find_data; - std::string search_path = directory + "\\*"; - HANDLE handle = FindFirstFileA(search_path.c_str(), &find_data); - - if(handle != INVALID_HANDLE_VALUE){ - do{ - if(!(find_data.dwFileAttributes & FILE_ATTRIBUTE_DIRECTORY)){ - std::string filename = find_data.cFileName; - if(is_supported_mesh_extension(filename)){ - files.push_back(filename); - } - } - }while(FindNextFileA(handle, &find_data)); - FindClose(handle); - } - #else - DIR *dir = opendir(directory.c_str()); - if(dir){ - struct dirent *entry; - while((entry = readdir(dir)) != nullptr){ - if(entry->d_type == DT_REG || entry->d_type == DT_UNKNOWN){ - std::string filename = entry->d_name; - if(is_supported_mesh_extension(filename)){ - // Check if it's actually a regular file (for DT_UNKNOWN) - std::string full_path = directory + "/" + filename; - struct stat st; - if(stat(full_path.c_str(), &st) == 0 && S_ISREG(st.st_mode)){ - files.push_back(filename); - } - } - } - } - closedir(dir); - } - #endif -#endif - std::sort(files.begin(), files.end()); - return files; + return success; } -std::string join_path(const std::string &dir, const std::string &filename) -{ -#ifdef HAS_FILESYSTEM - return (fs::path(dir) / filename).string(); -#else - #if defined(_WIN32) || defined(_WIN64) - return dir + "\\" + filename; - #else - return dir + "/" + filename; - #endif -#endif -} +} // anonymous namespace -bool create_directory_if_not_exists(const std::string &path) -{ -#ifdef HAS_FILESYSTEM - try{ - if(!fs::exists(path)){ - return fs::create_directories(path); - } - return true; - } - catch(const fs::filesystem_error &e){ - std::cerr << "Error creating directory: " << e.what() << std::endl; - return false; - } -#else - #if defined(_WIN32) || defined(_WIN64) - return CreateDirectoryA(path.c_str(), NULL) != 0 || GetLastError() == ERROR_ALREADY_EXISTS; - #else - struct stat st; - if(stat(path.c_str(), &st) == 0){ - return S_ISDIR(st.st_mode); - } - return mkdir(path.c_str(), 0755) == 0; - #endif -#endif -} -// Process a single mesh file -bool process_single_mesh(const std::string &input_mesh_path, - const std::string &output_mesh_path, - const SDFilter::MeshDenoisingParameters ¶m, - bool show_progress = true, - const std::string &indent = "") +int main(int argc, char **argv) { - TriMesh mesh; - bool load_success = false; - std::string warning; - std::string error; - - // Load mesh based on file extension - if(show_progress){ - ProgressIndicator progress(indent + "Loading mesh"); - if(has_extension(input_mesh_path, ".gltf") || has_extension(input_mesh_path, ".glb")){ - load_success = load_gltf_mesh(input_mesh_path, mesh, warning, error); - } - else if(has_extension(input_mesh_path, ".usd") || has_extension(input_mesh_path, ".usda") || - has_extension(input_mesh_path, ".usdc") || has_extension(input_mesh_path, ".usdz")){ - load_success = load_usd_mesh(input_mesh_path, mesh, warning, error); - } - else{ - load_success = OpenMesh::IO::read_mesh(mesh, input_mesh_path); - if(!load_success){ - error = "unable to read input mesh with OpenMesh"; - } - } - if(load_success){ - std::ostringstream oss; - oss << indent << "Loaded mesh (" << mesh.n_vertices() << " vertices, " << mesh.n_faces() << " faces)"; - progress.finish(oss.str()); - } - else{ - progress.finish(indent + "Loading failed"); - } + if(argc == 1) + { + print_help(argv[0]); + return 0; } - else{ - if(has_extension(input_mesh_path, ".gltf") || has_extension(input_mesh_path, ".glb")){ - load_success = load_gltf_mesh(input_mesh_path, mesh, warning, error); - } - else if(has_extension(input_mesh_path, ".usd") || has_extension(input_mesh_path, ".usda") || - has_extension(input_mesh_path, ".usdc") || has_extension(input_mesh_path, ".usdz")){ - load_success = load_usd_mesh(input_mesh_path, mesh, warning, error); - } - else{ - load_success = OpenMesh::IO::read_mesh(mesh, input_mesh_path); - if(!load_success){ - error = "unable to read input mesh with OpenMesh"; - } + + if(argc == 2) + { + const std::string arg = argv[1]; + if(arg == "--help" || arg == "-h"){ + print_help(argv[0]); + return 0; } } - if(!warning.empty()){ - std::cout << indent << "Warning: " << warning << std::endl; - } + if(argc == 3) + { + const std::string arg = argv[1]; + if(arg == "--write-default-options") + { + if(!write_default_options_file(argv[2])){ + std::cerr << "Error: unable to write default options file " << argv[2] << std::endl; + return 1; + } - if(!load_success){ - std::cerr << indent << "Error: " << error << std::endl; - return false; + std::cout << "Wrote default options to " << argv[2] << std::endl; + return 0; + } } - // Normalize the input mesh - Eigen::Vector3d original_center; - double original_scale; - SDFilter::normalize_mesh(mesh, original_center, original_scale); - - // Filter the normals and construct the output mesh + // Count positional arguments (everything before the first --) + int positional_argc = 0; + for(int i = 1; i < argc; ++i) { - ProgressIndicator progress(indent + "Denoising mesh", show_progress); - SDFilter::MeshNormalDenoising denoiser(mesh); - TriMesh output_mesh; - if(!denoiser.denoise(param, output_mesh)){ - progress.finish(indent + "Denoising failed"); - std::cerr << indent << "Error: denoising failed" << std::endl; - return false; + const std::string arg = argv[i]; + if(arg.rfind("--", 0) == 0){ + break; } - SDFilter::restore_mesh(output_mesh, original_center, original_scale); - - // Save output mesh - progress.set_message(indent + "Saving mesh"); - if(!SDFilter::write_mesh_high_accuracy(output_mesh, output_mesh_path)){ - progress.finish(indent + "Saving failed"); - std::cerr << indent << "Error: unable to save result mesh to " << output_mesh_path << std::endl; - return false; - } - progress.finish(indent + "Complete"); + ++positional_argc; } - return true; -} - -// Batch process all mesh files in a directory -int process_batch(const std::string &option_file, - const std::string &input_dir, - const std::string &output_dir) -{ - // Load option file - SDFilter::MeshDenoisingParameters param; - if(!param.load(option_file.c_str())){ - std::cerr << "Error: unable to load option file " << option_file << std::endl; - return 1; - } - if(!param.valid_parameters()){ - std::cerr << "Invalid filter options. Aborting..." << std::endl; + if(positional_argc != 2 && positional_argc != 3) + { + print_help(argv[0]); return 1; } - std::cout << "Batch processing mode" << std::endl; - std::cout << "Input directory: " << input_dir << std::endl; - std::cout << "Output directory: " << output_dir << std::endl; - std::cout << std::endl; - param.output(); - std::cout << std::endl; + std::string options_path; + std::string input_mesh_path; + std::string output_mesh_path; + int flag_index = 1; + bool using_builtin_defaults = false; - // Create output directory if it doesn't exist - if(!create_directory_if_not_exists(output_dir)){ - std::cerr << "Error: unable to create output directory " << output_dir << std::endl; - return 1; + if(positional_argc == 2) + { + using_builtin_defaults = true; + input_mesh_path = argv[1]; + output_mesh_path = argv[2]; + flag_index = 3; } - - // List all mesh files in input directory - std::vector mesh_files = list_mesh_files(input_dir); - - if(mesh_files.empty()){ - std::cerr << "Error: no supported mesh files found in " << input_dir << std::endl; - return 1; + else + { + options_path = argv[1]; + input_mesh_path = argv[2]; + output_mesh_path = argv[3]; + flag_index = 4; } - std::cout << "Found " << mesh_files.size() << " mesh file(s) to process" << std::endl; - std::cout << std::endl; + int export_precision = 16; + std::string metrics_json_path; + std::string metrics_csv_path; + bool deterministic_cli = false; - int success_count = 0; - int failure_count = 0; - auto batch_start = std::chrono::steady_clock::now(); - - // Process each file - for(size_t i = 0; i < mesh_files.size(); ++i){ - const std::string &filename = mesh_files[i]; - std::string input_path = join_path(input_dir, filename); - std::string output_path = join_path(output_dir, filename); - - std::cout << "[" << (i + 1) << "/" << mesh_files.size() << "] " << filename << std::endl; - - if(process_single_mesh(input_path, output_path, param, true, " ")){ - ++success_count; + for(int i = flag_index; i < argc; ++i) + { + std::string arg = argv[i]; + if(arg == "--export-precision" && i + 1 < argc){ + export_precision = std::stoi(argv[++i]); + } + else if(arg == "--metrics-json" && i + 1 < argc){ + metrics_json_path = argv[++i]; + } + else if(arg == "--metrics-csv" && i + 1 < argc){ + metrics_csv_path = argv[++i]; + } + else if(arg == "--deterministic"){ + deterministic_cli = true; } else{ - ++failure_count; + std::cerr << "Unknown or incomplete argument: " << arg << std::endl; + return 1; } - std::cout << std::endl; - } - - // Summary with total elapsed time - auto batch_end = std::chrono::steady_clock::now(); - auto total_elapsed = std::chrono::duration_cast(batch_end - batch_start).count(); - - std::cout << "==================================" << std::endl; - std::cout << "Batch processing complete" << std::endl; - std::cout << " Successful: " << success_count << std::endl; - std::cout << " Failed: " << failure_count << std::endl; - std::cout << " Total: " << mesh_files.size() << std::endl; - std::cout << " Time: " << std::fixed << std::setprecision(2) << (total_elapsed / 1000.0) << "s"; - if(success_count > 0){ - std::cout << " (avg: " << std::fixed << std::setprecision(2) << (total_elapsed / 1000.0 / success_count) << "s per file)"; } - std::cout << std::endl; - - return failure_count > 0 ? 1 : 0; -} - -} // anonymous namespace + using Clock = std::chrono::steady_clock; + auto app_begin = Clock::now(); -int main(int argc, char **argv) -{ - if(argc != 4) + // Load input mesh + TriMesh mesh; + auto import_begin = Clock::now(); + if(!load_mesh(input_mesh_path, mesh)) { - std::cout << "MeshDenoiser - Mesh normal denoising filter" << std::endl; - std::cout << std::endl; - std::cout << "Usage: MeshDenoiser OPTION_FILE INPUT OUTPUT" << std::endl; - std::cout << std::endl; - std::cout << "Single file mode:" << std::endl; - std::cout << " MeshDenoiser OPTION_FILE INPUT_MESH OUTPUT_MESH" << std::endl; - std::cout << std::endl; - std::cout << "Batch processing mode:" << std::endl; - std::cout << " MeshDenoiser OPTION_FILE INPUT_DIR/ OUTPUT_DIR/" << std::endl; - std::cout << " (processes all mesh files in INPUT_DIR)" << std::endl; - std::cout << std::endl; - std::cout << "Supported input formats:" << std::endl; - std::cout << " - OBJ, PLY, OFF, STL (via OpenMesh)" << std::endl; - std::cout << " - glTF (.gltf, .glb)" << std::endl; - std::cout << " - USD (.usd, .usda, .usdc, .usdz)" << std::endl; - std::cout << std::endl; - std::cout << "Examples:" << std::endl; - std::cout << " MeshDenoiser options.txt input.obj output.obj" << std::endl; - std::cout << " MeshDenoiser options.txt model.gltf denoised.obj" << std::endl; - std::cout << " MeshDenoiser options.txt scene.usdz clean.ply" << std::endl; - std::cout << " MeshDenoiser options.txt input_dir/ output_dir/" << std::endl; return 1; } + auto import_end = Clock::now(); - const std::string option_file = argv[1]; - const std::string input_path = argv[2]; - const std::string output_path = argv[3]; + std::cout << "Loaded " << mesh.n_vertices() << " vertices, " << mesh.n_faces() << " faces." << std::endl; + + // Validate vertex coordinates + if(!validate_mesh_vertices(mesh)){ + return 1; + } #ifdef USE_OPENMP Eigen::initParallel(); #endif - // Determine processing mode based on input type - bool input_is_dir = is_directory(input_path); - - if(input_is_dir){ - // Batch processing mode: input is a directory - // Check if output exists and is a regular file (error case) - if(is_regular_file(output_path)){ - std::cerr << "Error: Input is a directory but output is an existing file." << std::endl; - std::cerr << "For batch processing, output must be a directory path (will be created if needed)." << std::endl; - return 1; - } - // Output is either a directory or doesn't exist yet (will be created in process_batch) - return process_batch(option_file, input_path, output_path); - } - else{ - // Single file processing mode: input is a file - // Check if output is an existing directory (error case - need filename) - if(is_directory(output_path)){ - std::cerr << "Error: Input is a file but output is a directory." << std::endl; - std::cerr << "For single file mode, output must be a file path (e.g., output_dir/file.obj)." << std::endl; - return 1; - } - } - - // Single file processing mode - // Load option file + // Load options SDFilter::MeshDenoisingParameters param; - if(!param.load(option_file.c_str())){ - std::cerr << "Error: unable to load option file " << option_file << std::endl; + if(using_builtin_defaults) + { + apply_default_options(param); + std::cout << "No option file specified; using built-in default denoising options." << std::endl; + } + else if(!param.load(options_path.c_str())){ + std::cerr << "Error: unable to load option file " << options_path << std::endl; return 1; } if(!param.valid_parameters()){ std::cerr << "Invalid filter options. Aborting..." << std::endl; return 1; } + if(deterministic_cli){ + param.deterministic_mode = true; + } + configure_deterministic_runtime(param.deterministic_mode); param.output(); - // Process the single mesh file - if(!process_single_mesh(input_path, output_path, param)){ + // Normalize the input mesh + Eigen::Vector3d original_center; + double original_scale; + auto normalize_begin = Clock::now(); + SDFilter::normalize_mesh(mesh, original_center, original_scale); + auto normalize_end = Clock::now(); + + // Filter the normals and construct the output mesh + SDFilter::MeshNormalDenoising denoiser(mesh); + TriMesh output_mesh; + auto denoise_begin = Clock::now(); + if(!denoiser.denoise(param, output_mesh)){ + std::cerr << "Error: denoising failed." << std::endl; return 1; } + auto denoise_end = Clock::now(); + + auto restore_begin = Clock::now(); + SDFilter::restore_mesh(output_mesh, original_center, original_scale); + auto restore_end = Clock::now(); + + // Save output mesh + auto export_begin = Clock::now(); + if(!SDFilter::write_mesh(output_mesh, output_mesh_path, export_precision)){ + std::cerr << "Error: unable to save the result mesh to file " << output_mesh_path << std::endl; + return 1; + } + auto export_end = Clock::now(); + auto app_end = Clock::now(); + + // Collect and report metrics + SDFilter::PipelineMetrics metrics; + metrics.mode = "denoise"; + metrics.input_mesh = input_mesh_path; + metrics.output_mesh = output_mesh_path; + metrics.import_secs = std::chrono::duration(import_end - import_begin).count(); + metrics.normalize_secs = std::chrono::duration(normalize_end - normalize_begin).count(); + metrics.algorithm_secs = std::chrono::duration(denoise_end - denoise_begin).count(); + metrics.restore_secs = std::chrono::duration(restore_end - restore_begin).count(); + metrics.export_secs = std::chrono::duration(export_end - export_begin).count(); + metrics.total_secs = std::chrono::duration(app_end - app_begin).count(); + metrics.obj_export_precision = export_precision; + + std::cout << std::endl; + std::cout << "====== Timing =========" << std::endl; + std::cout << "Import: " << metrics.import_secs << " s" << std::endl; + std::cout << "Normalize: " << metrics.normalize_secs << " s" << std::endl; + std::cout << "Denoise: " << metrics.algorithm_secs << " s" << std::endl; + std::cout << "Restore: " << metrics.restore_secs << " s" << std::endl; + std::cout << "Export: " << metrics.export_secs << " s" << std::endl; + std::cout << "Total: " << metrics.total_secs << " s" << std::endl; + std::cout << "=======================" << std::endl; + + const SDFilter::RunStatistics &stats = denoiser.run_stats(); + if(!metrics_json_path.empty()){ + if(!SDFilter::write_metrics_json(metrics_json_path, metrics, stats)){ + std::cerr << "Warning: unable to write metrics JSON file " << metrics_json_path << std::endl; + } + } + + if(!metrics_csv_path.empty()){ + bool write_header = !std::filesystem::exists(metrics_csv_path); + if(!SDFilter::write_metrics_csv(metrics_csv_path, metrics, stats, write_header)){ + std::cerr << "Warning: unable to write metrics CSV file " << metrics_csv_path << std::endl; + } + } - std::cout << "Denoising complete!" << std::endl; return 0; } diff --git a/src/MeshNormalDenoising.h b/src/MeshNormalDenoising.h new file mode 100644 index 0000000..e8d68bc --- /dev/null +++ b/src/MeshNormalDenoising.h @@ -0,0 +1,292 @@ +// BSD 3-Clause License +// +// Copyright (c) 2017, Bailin Deng +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// * Redistributions of source code must retain the above copyright notice, this +// list of conditions and the following disclaimer. +// +// * Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// * Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +// OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + +#ifndef MESHNORMALDENOISING_H_ +#define MESHNORMALDENOISING_H_ + +#include "MeshNormalFilter.h" +#include + + +namespace SDFilter +{ + +class MeshDenoisingParameters: public MeshFilterParameters +{ +public: + MeshDenoisingParameters() + :outer_iterations(5){} + + virtual ~MeshDenoisingParameters(){} + + int outer_iterations; // Number of outer iterations where the SD filter is applied + + virtual bool valid_parameters() const + { + if(!MeshFilterParameters::valid_parameters()){ + return false; + } + + if(outer_iterations <= 0){ + std::cerr << "Error: outer_iterations must be positive" << std::endl; + return false; + } + + return true; + } + +protected: + + virtual bool load_option(const OptionInterpreter &opt) + { + return MeshFilterParameters::load_option(opt) || + opt.load("OuterIterations", outer_iterations); + } + + virtual void output_options() + { + Parameters::output_options(); + std::cout << "Mesh denoising outer iterations: " << outer_iterations << std::endl; + } +}; + + +class MeshNormalDenoising : public MeshNormalFilter +{ +public: + MeshNormalDenoising(const TriMesh &mesh) + :MeshNormalFilter(mesh) + { + print_progress_ = false; + print_diagnostic_info_ = false; + print_timing_ = false; + print_error_evaluation_ = false; + } + + + bool denoise(const MeshDenoisingParameters ¶m, TriMesh &output_mesh) + { + assert(param.valid_parameters()); + + std::cout << "Denoising started" << std::endl; + + Timer timer; + Timer::EventID denoise_begin_time = timer.get_time(); + double preprocessing_sum = 0.0; + double filtering_sum = 0.0; + double mesh_update_sum = 0.0; + double mesh_filter_total_sum = 0.0; + + for(int i = 0; i < param.outer_iterations; ++ i) + { + std::cout << "Outer iteration " << (i+1) << "..." << std::endl; + + if(!filter(param, output_mesh)){ + std::cerr << "Unable to perform normal filter. Denoising aborted." << std::endl; + return false; + } + const RunStatistics &iter_stats = run_stats(); + preprocessing_sum += iter_stats.preprocessing_secs; + filtering_sum += iter_stats.filtering_secs; + mesh_update_sum += iter_stats.mesh_update_secs; + mesh_filter_total_sum += iter_stats.mesh_filter_total_secs; + + // Use the filtered mesh as the input mesh for the next outer iteration + set_mesh(output_mesh, false); + } + + Timer::EventID denoise_end_time = timer.get_time(); + stats_.preprocessing_secs = preprocessing_sum; + stats_.filtering_secs = filtering_sum; + stats_.mesh_update_secs = mesh_update_sum; + stats_.mesh_filter_total_secs = mesh_filter_total_sum; + stats_.denoise_total_secs = timer.elapsed_time(denoise_begin_time, denoise_end_time); + stats_.outer_iterations = param.outer_iterations; + std::cout << "Denoising completed, timing: " << + stats_.denoise_total_secs << std::endl; + + return true; + } + + +protected: + + virtual void get_initial_data(Eigen::MatrixXd &guidance_normals, Eigen::MatrixXd &init_normals, Eigen::VectorXd &area_weights) + { + // Call the base class method to fill in initial data + MeshNormalFilter::get_initial_data(guidance_normals, init_normals, area_weights); + + // Update the guidance to patch-based normals + + int face_count = mesh_.n_faces(), edge_count = mesh_.n_edges(), vtx_count = mesh_.n_vertices(); + + std::vector edge_saliency(edge_count, 0); // Pre-computed edge saliency, defined as difference between adjacent normals + + std::vector< std::vector > adj_faces_per_vtx(vtx_count); + std::vector< std::vector > adj_nonboundary_edges_per_vtx(vtx_count); + std::vector< std::vector > neighborhood_faces_per_face(face_count); + + std::vector patch_normal_consistency(face_count, 0); + Eigen::Matrix3Xd patch_avg_normal(3, face_count); + std::vector face_stamp(face_count, -1), edge_stamp(edge_count, -1); + + double epsilon = 1e-9; + + OMP_PARALLEL + { + OMP_FOR + for(int i = 0; i < edge_count; ++ i) + { + TriMesh::EdgeHandle eh(i); + + if(!mesh_.is_boundary(eh)) + { + int f1 = mesh_.face_handle(mesh_.halfedge_handle(eh, 0)).idx(); + int f2 = mesh_.face_handle(mesh_.halfedge_handle(eh, 1)).idx(); + edge_saliency[i] = (init_normals.col(f1) - init_normals.col(f2)).norm(); + } + } + + OMP_FOR + for(int i = 0; i < vtx_count; ++ i) + { + // Collect neighboring faces and non-boundary edges for each vertex + + TriMesh::VertexHandle vh(i); + + for(TriMesh::ConstVertexFaceIter cvf_it = mesh_.cvf_iter(vh); cvf_it.is_valid(); ++ cvf_it){ + adj_faces_per_vtx[i].push_back(cvf_it->idx()); + } + + for(TriMesh::ConstVertexEdgeIter cve_it = mesh_.cve_iter(vh); cve_it.is_valid(); ++ cve_it){ + if(!mesh_.is_boundary(*cve_it)){ + adj_nonboundary_edges_per_vtx[i].push_back(cve_it->idx()); + } + } + } + } + + for(int i = 0; i < face_count; ++ i) + { + // Each candidate patch is associated with a face at its center. + // We collect all the faces and non-boundary edges within such a patch. + TriMesh::FaceHandle fh(i); + std::vector faces_in_patch, edges_in_patch; + + for(TriMesh::ConstFaceVertexIter cfv_it = mesh_.cfv_iter(fh); cfv_it.is_valid(); ++ cfv_it){ + int vtx_idx = cfv_it->idx(); + const std::vector &adj_faces = adj_faces_per_vtx[vtx_idx]; + for(std::size_t k = 0; k < adj_faces.size(); ++k) + { + int f = adj_faces[k]; + if(face_stamp[f] != i){ + face_stamp[f] = i; + faces_in_patch.push_back(f); + } + } + + const std::vector &adj_edges = adj_nonboundary_edges_per_vtx[vtx_idx]; + for(std::size_t k = 0; k < adj_edges.size(); ++k) + { + int e = adj_edges[k]; + if(edge_stamp[e] != i){ + edge_stamp[e] = i; + edges_in_patch.push_back(e); + } + } + } + + neighborhood_faces_per_face[i] = faces_in_patch; + + // Collect face normals and edge saliency values from the patch, and compute patch normal consistency value and average normal + int n_faces_in_patch = faces_in_patch.size(); + int n_edges_in_patch = edges_in_patch.size(); + + Eigen::Matrix3Xd face_normals(3, n_faces_in_patch); + Eigen::VectorXd face_area(n_faces_in_patch); + Eigen::VectorXd edge_saliency_values(n_edges_in_patch); + + for(int k = 0; k < n_edges_in_patch; ++ k){ + edge_saliency_values(k) = edge_saliency[ edges_in_patch[k] ]; + } + + if(n_edges_in_patch > 0){ + patch_normal_consistency[i] = edge_saliency_values.maxCoeff() / (epsilon + edge_saliency_values.sum()); + } + + for(int k = 0; k < n_faces_in_patch; ++ k) + { + int f_k = faces_in_patch[k]; + face_normals.col(k) = init_normals.col(f_k); + face_area(k) = area_weights(f_k); + } + + // Find the max normal difference within a patch + double max_normal_diff = 0; + for(int k = 0; k < n_faces_in_patch - 1; ++ k) + { + Eigen::Matrix3Xd N = face_normals.block(0, k+1, 3, n_faces_in_patch-k-1); + N.colwise() -= face_normals.col(k); + double max_diff = N.colwise().norm().maxCoeff(); + max_normal_diff = std::max(max_normal_diff, max_diff); + } + + patch_normal_consistency[i] *= max_normal_diff; + patch_avg_normal.col(i) = (face_normals * face_area).normalized(); + } + + OMP_PARALLEL + { + OMP_FOR + for(int i = 0; i < face_count; ++ i) + { + // For each face, select the patch with the most consistent normals to construct its guidance + std::vector &neighborhood_faces = neighborhood_faces_per_face[i]; + assert(!neighborhood_faces.empty()); + + Eigen::VectorXd patch_scores(neighborhood_faces.size()); + for(int k= 0; k < static_cast(neighborhood_faces.size()); ++ k){ + patch_scores(k) = patch_normal_consistency[neighborhood_faces[k]]; + } + + int best_score_idx = -1; + patch_scores.minCoeff(&best_score_idx); + guidance_normals.col(i) = patch_avg_normal.col( neighborhood_faces[best_score_idx] ); + } + } + } +}; + +} + + + +#endif /* MESHNORMALDENOISING_H_ */ diff --git a/src/MeshNormalFilter.h b/src/MeshNormalFilter.h new file mode 100644 index 0000000..7b63de8 --- /dev/null +++ b/src/MeshNormalFilter.h @@ -0,0 +1,810 @@ +// BSD 3-Clause License +// +// Copyright (c) 2017, Bailin Deng +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// * Redistributions of source code must retain the above copyright notice, this +// list of conditions and the following disclaimer. +// +// * Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// * Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +// OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + +#ifndef MESHNORMALFILTER_H_ +#define MESHNORMALFILTER_H_ + +#include "MeshTypes.h" +#include "SDFilter.h" +#include +#include +#include +#include + +namespace SDFilter +{ + +class MeshFilterParameters : public Parameters +{ +public: + MeshFilterParameters() + :mesh_update_method(ITERATIVE_UPDATE), mesh_update_closeness_weight(0.001), mesh_update_iter(20), + mesh_update_disp_eps(0.0) + { + // Use an threshold value corresponding to eps_angle_degree degrees change of normal vectors between two iterations + double eps_angle_degree = 0.2; + avg_disp_eps = 2 * std::sin(eps_angle_degree * 0.5 * M_PI / 180); + } + + virtual ~MeshFilterParameters(){} + + // Methods for mesh vertex update according to filtered normals + enum MeshUpdateMethod + { + ITERATIVE_UPDATE, // ShapeUp styled iterative solver + POISSON_UPDATE, // Poisson-based update from [Wang et al. 2015] + MESH_UPDATE_METHOD_COUNT + }; + + MeshUpdateMethod mesh_update_method; + + double mesh_update_closeness_weight; // Weight for the closeness term in mesh update + + int mesh_update_iter; // Number of mesh update iterations + + double mesh_update_disp_eps; // Early-stop threshold on RMS vertex displacement per mesh-update iteration; <=0 disables early stop + + virtual bool valid_parameters() const + { + if(!Parameters::valid_parameters()){ + return false; + } + + if(mesh_update_iter <= 0){ + std::cerr << "Error: MeshUpdateIterations must be positive" << std::endl; + return false; + } + + if(mesh_update_closeness_weight < 0){ + std::cerr << "Error: MeshUpdateClosenessWeight must be positive" << std::endl; + return false; + } + + if(mesh_update_disp_eps < 0.0){ + std::cerr << "Error: MeshUpdateDisplacementEps must be non-negative" << std::endl; + return false; + } + + if(mesh_update_method < ITERATIVE_UPDATE || mesh_update_method >= MESH_UPDATE_METHOD_COUNT){ + std::cerr << "Error: invalid MeshUpdateMethod" << std::endl; + return false; + } + + return true; + } + +protected: + virtual bool load_option(const OptionInterpreter &opt) + { + return Parameters::load_option(opt) || + opt.load_enum("MeshUpdateMethod", MESH_UPDATE_METHOD_COUNT, mesh_update_method) || + opt.load("MeshUpdateClosenessWeight", mesh_update_closeness_weight) || + opt.load("MeshUpdateIterations", mesh_update_iter) || + opt.load("MeshUpdateDisplacementEps", mesh_update_disp_eps); + } + + virtual void output_options() + { + Parameters::output_options(); + std::cout << "Mesh update method: " << static_cast(mesh_update_method) << std::endl; + std::cout << "Mesh update closeness weight: " << mesh_update_closeness_weight << std::endl; + std::cout << "Mesh update iterations: " << mesh_update_iter << std::endl; + std::cout << "Mesh update displacement eps: " << mesh_update_disp_eps << std::endl; + } +}; + +class MeshNormalFilter : public SDFilter +{ +public: + + MeshNormalFilter(const TriMesh &mesh) + :mesh_(mesh), print_error_evaluation_(false), linear_solver_(Parameters::LDLT), system_matrix_factorized_(false), + mesh_update_system_closeness_weight_(std::numeric_limits::quiet_NaN()), + mesh_update_system_n_faces_(-1), mesh_update_system_n_vtx_(-1), + cached_face_vtx_idx_valid_(false), cached_face_vtx_n_faces_(-1), cached_face_vtx_n_vtx_(-1){} + + virtual ~MeshNormalFilter() {} + + // Pass param by value, to allow changing the value of eta + bool filter(MeshFilterParameters param, TriMesh &output_mesh) + { + assert(param.valid_parameters()); + linear_solver_.set_solver_type(param.linear_solver_type); + linear_solver_.set_solver_params(param.linear_solver_max_iterations, param.linear_solver_tolerance); + + Timer timer; + Timer::EventID mesh_flter_begin_time = timer.get_time(); + + // Rescale the eta parameter, according to average distance between neighboring face centroids + param.eta *= average_neighbor_face_centroid_dist(mesh_); + + if(!SDFilter::filter(param)){ + std::cerr << "Error in performing SD filter" << std::endl; + return false; + } + + + Timer::EventID update_begin_time = timer.get_time(); + + // Normalize the filtered face normals + Matrix3X target_normals = signals_.block(0, 0, 3, signals_.cols()); + target_normals.colwise().normalize(); + + if(param.mesh_update_method == MeshFilterParameters::ITERATIVE_UPDATE){ + if(!iterative_mesh_update(param, target_normals, output_mesh)){ + std::cerr << "Error in iteative mesh update" << std::endl; + return false; + } + } + else{ + if(!Poisson_mesh_update(target_normals, output_mesh)){ + std::cerr << "Error in Poisson mesh update" << std::endl; + return false; + } + } + + Timer::EventID update_end_time = timer.get_time(); + stats_.mesh_update_secs = timer.elapsed_time(update_begin_time, update_end_time); + stats_.mesh_filter_total_secs = timer.elapsed_time(mesh_flter_begin_time, update_end_time); + + if(print_timing_){ + std::cout << "Mesh udpate timing: " << stats_.mesh_update_secs << " secs" << std::endl; + std::cout << "Mesh filter total timing: " << stats_.mesh_filter_total_secs << " secs" << std::endl; + } + + if(print_error_evaluation_) + { + std::cout << std::endl; + show_normalized_mesh_displacement_norm(output_mesh); + show_normal_error_statistics(output_mesh, target_normals, 2, 10); + } + + return true; + } + +protected: + + TriMesh mesh_; + + bool print_error_evaluation_; // The printing of mesh update error + + bool get_neighborhood(const Parameters ¶m, Eigen::Matrix2Xi &neighbor_pairs, Eigen::VectorXd &neighbor_dist) + { + const double radius = 3.0 * param.eta; + const double radius_sqr = radius * radius; + const int n_faces = mesh_.n_faces(); + if(n_faces <= 0){ + return false; + } + + Matrix3X face_centroids(3, n_faces); + OMP_PARALLEL + { + OMP_FOR + for(int i = 0; i < n_faces; ++i){ + face_centroids.col(i) = to_eigen_vec3d(mesh_.calc_face_centroid(TriMesh::FaceHandle(i))); + } + } + + using CellKey = std::array; + struct CellKeyHash + { + std::size_t operator()(const CellKey &k) const + { + std::size_t h = static_cast(k[0]) * 73856093u; + h ^= static_cast(k[1]) * 19349663u; + h ^= static_cast(k[2]) * 83492791u; + return h; + } + }; + + Eigen::Vector3d min_c = face_centroids.rowwise().minCoeff(); + std::unordered_map, CellKeyHash> cells; + cells.reserve(static_cast(n_faces)); + + auto centroid_cell = [&](const Eigen::Vector3d &p) -> CellKey + { + Eigen::Vector3d q = (p - min_c) / radius; + return CellKey{ + static_cast(std::floor(q[0])), + static_cast(std::floor(q[1])), + static_cast(std::floor(q[2])) + }; + }; + + for(int i = 0; i < n_faces; ++i){ + cells[centroid_cell(face_centroids.col(i))].push_back(i); + } + + std::vector> pairs; + pairs.reserve(static_cast(n_faces) * 12u); + + for(int i = 0; i < n_faces; ++i) + { + const Eigen::Vector3d ci = face_centroids.col(i); + const CellKey c = centroid_cell(ci); + + for(int dx = -1; dx <= 1; ++dx) + { + for(int dy = -1; dy <= 1; ++dy) + { + for(int dz = -1; dz <= 1; ++dz) + { + CellKey nkey{c[0] + dx, c[1] + dy, c[2] + dz}; + auto it = cells.find(nkey); + if(it == cells.end()){ + continue; + } + + const std::vector &bucket = it->second; + for(std::size_t k = 0; k < bucket.size(); ++k) + { + const int j = bucket[k]; + if(j <= i){ + continue; + } + + if((ci - face_centroids.col(j)).squaredNorm() < radius_sqr){ + pairs.emplace_back(i, j); + } + } + } + } + } + } + + const Eigen::Index n_neighbor_pairs = static_cast(pairs.size()); + if(n_neighbor_pairs <= 0){ + return false; + } + + neighbor_pairs.resize(2, n_neighbor_pairs); + neighbor_dist.resize(n_neighbor_pairs); + for(Eigen::Index p = 0; p < n_neighbor_pairs; ++p) + { + const int i = pairs[static_cast(p)].first; + const int j = pairs[static_cast(p)].second; + neighbor_pairs(0, p) = i; + neighbor_pairs(1, p) = j; + neighbor_dist(p) = (face_centroids.col(i) - face_centroids.col(j)).norm(); + } + + return true; + } + + + void get_initial_data(Eigen::MatrixXd &guidance, Eigen::MatrixXd &init_signals, Eigen::VectorXd &area_weights) + { + init_signals.resize(3, mesh_.n_faces()); + + for(TriMesh::ConstFaceIter cf_it = mesh_.faces_begin(); cf_it != mesh_.faces_end(); ++ cf_it) + { + Eigen::Vector3d f_normal = to_eigen_vec3d(mesh_.calc_face_normal(*cf_it)).normalized(); + init_signals.col(cf_it->idx()) = f_normal; + } + + guidance = init_signals; + + get_face_area_weights(mesh_, area_weights); + } + + + void reset_mesh_update_system() + { + system_matrix_factorized_ = false; + mesh_update_system_closeness_weight_ = std::numeric_limits::quiet_NaN(); + mesh_update_system_n_faces_ = -1; + mesh_update_system_n_vtx_ = -1; + } + + + void set_mesh(const TriMesh &mesh, bool invalidate_update_system) + { + mesh_ = mesh; + cached_face_vtx_idx_valid_ = false; + cached_face_vtx_n_faces_ = -1; + cached_face_vtx_n_vtx_ = -1; + + if(invalidate_update_system){ + reset_mesh_update_system(); + } + } + + + +private: + + // Pre-computed information for mesh update problem + // \min_X ||AX - B||^2 + w ||X - X0||^2, + // where X are new mesh vertex positions, A is a mean-centering matrix, B are the target mean-centered positions, + // X0 are initial vertex positions, and w > 0 is a closeness weight. + // This amounts to solving a linear system + // (A^T A + w I) X = A^T B + w X0. + // We precompute matrix A^T, and pre-factorize A^T A + w I + LinearSolver linear_solver_; // Linear system solver for mesh update + SparseMatrixXd At_; // Transpose of part of the linear least squares matrix that corresponds to mean centering of face vertices + bool system_matrix_factorized_; // Whether the matrix + double mesh_update_system_closeness_weight_; + int mesh_update_system_n_faces_; + int mesh_update_system_n_vtx_; + Matrix3Xi cached_face_vtx_idx_; + bool cached_face_vtx_idx_valid_; + int cached_face_vtx_n_faces_; + int cached_face_vtx_n_vtx_; + + + + + // Set up and pre-factorize the linear system for iterative mesh update + bool setup_mesh_udpate_system(const Matrix3Xi &face_vtx_idx, double w_closeness) + { + const int n_faces = mesh_.n_faces(); + const int n_vtx = mesh_.n_vertices(); + const bool same_problem_size = (mesh_update_system_n_faces_ == n_faces && mesh_update_system_n_vtx_ == n_vtx); + const bool same_closeness = (std::isfinite(mesh_update_system_closeness_weight_) && + std::abs(mesh_update_system_closeness_weight_ - w_closeness) <= std::numeric_limits::epsilon()); + if(system_matrix_factorized_ && same_problem_size && same_closeness) + { + return true; + } + + std::vector A_triplets(9 * n_faces); + std::vector I_triplets(n_vtx); + + // Matrix for mean centering of three vertices + Eigen::Matrix3d mean_centering_mat; + get_mean_centering_matrix(mean_centering_mat); + + OMP_PARALLEL + { + OMP_FOR + for(int i = 0; i < n_faces; ++ i) + { + Eigen::Vector3i vtx_idx = face_vtx_idx.col(i); + + int triplet_addr = 9 * i; + int row_idx = 3 * i; + for(int j = 0; j < 3; ++ j) + { + for(int k = 0; k < 3; ++ k){ + A_triplets[triplet_addr++] = Triplet(row_idx, vtx_idx(k), mean_centering_mat(j, k)); + } + + row_idx++; + } + } + + OMP_FOR + for(int i = 0; i < n_vtx; ++ i) + { + I_triplets[i] = Triplet(i, i, w_closeness); + } + } + + + SparseMatrixXd A(3 * n_faces, n_vtx); + A.setFromTriplets(A_triplets.begin(), A_triplets.end()); + At_ = A.transpose(); + At_.makeCompressed(); + + SparseMatrixXd wI(n_vtx, n_vtx); + wI.setFromTriplets(I_triplets.begin(), I_triplets.end()); + SparseMatrixXd M = At_ * A + wI; + + linear_solver_.reset_pattern(); + if(!linear_solver_.compute(M)){ + std::cerr << "Error: failed to pre-factorize mesh update system" << std::endl; + return false; + } + + system_matrix_factorized_ = true; + mesh_update_system_closeness_weight_ = w_closeness; + mesh_update_system_n_faces_ = n_faces; + mesh_update_system_n_vtx_ = n_vtx; + return true; + } + + void get_cached_face_vertex_indices(const TriMesh &mesh, Matrix3Xi &face_vtx_idx) + { + const int n_faces = mesh.n_faces(); + const int n_vtx = mesh.n_vertices(); + if(!cached_face_vtx_idx_valid_ || cached_face_vtx_n_faces_ != n_faces || cached_face_vtx_n_vtx_ != n_vtx) + { + get_face_vertex_indices(mesh, cached_face_vtx_idx_); + cached_face_vtx_idx_valid_ = true; + cached_face_vtx_n_faces_ = n_faces; + cached_face_vtx_n_vtx_ = n_vtx; + } + face_vtx_idx = cached_face_vtx_idx_; + } + + + void get_face_area_weights(const TriMesh &mesh, Eigen::VectorXd &face_area_weights) const + { + face_area_weights.resize(mesh.n_faces()); + + for(TriMesh::ConstFaceIter cf_it = mesh.faces_begin(); cf_it != mesh.faces_end(); ++ cf_it) + { + face_area_weights(cf_it->idx()) = mesh.calc_sector_area(mesh.halfedge_handle(*cf_it)); + } + + face_area_weights /= face_area_weights.mean(); + } + + + bool iterative_mesh_update(const MeshFilterParameters ¶m, const Matrix3X &target_normals, TriMesh &output_mesh) + { + // Rescale closeness weight using the ratio between face number and vertex number, and take its square root + double w_closeness = param.mesh_update_closeness_weight * double(mesh_.n_faces()) / mesh_.n_vertices(); + + output_mesh = mesh_; + + Matrix3Xi face_vtx_idx; + get_cached_face_vertex_indices(output_mesh, face_vtx_idx); + + if(!setup_mesh_udpate_system(face_vtx_idx, w_closeness)){ + return false; + } + + std::cout << "Starting iterative mesh update......" << std::endl; + + Matrix3X vtx_pos; + get_vertex_points(output_mesh, vtx_pos); + + int n_faces = output_mesh.n_faces(); + Eigen::Matrix3Xd target_plane_local_frames(3, 2 * n_faces); // Local frame for the target plane of each face + std::vector local_frame_initialized(n_faces, false); + + using MatrixX3dRM = Eigen::Matrix; + MatrixX3dRM wX0 = vtx_pos.transpose() * w_closeness; // Part of the linear system right-hand-side that corresponds to initial vertex positions + MatrixX3dRM B(3 * n_faces, 3); // Per-face target position of the new vertices + + int n_vtx = output_mesh.n_vertices(); + MatrixX3dRM rhs(n_vtx, 3), sol(n_vtx, 3); + + for(int iter = 0; iter < param.mesh_update_iter; ++ iter) + { + OMP_PARALLEL + { + OMP_FOR + for(int i = 0; i < n_faces; ++ i) + { + Eigen::Vector3d current_normal = to_eigen_vec3d(output_mesh.calc_face_normal(TriMesh::FaceHandle(i))); + Eigen::Vector3d target_normal = target_normals.col(i); + + Eigen::Matrix3d face_vtx_pos; + get_mean_centered_face_vtx_pos(vtx_pos, face_vtx_idx.col(i), face_vtx_pos); + + Eigen::Matrix3d target_pos; + + // If the current normal is not pointing away from the target normal, simply project the points onto the target plane + if(current_normal.dot(target_normal) >= 0){ + target_pos = face_vtx_pos - target_normal * (target_normal.transpose() * face_vtx_pos); + } + else{ + // Otherwise, project the points onto a line in the target plane + typedef Eigen::Matrix Matrix32d; + Matrix32d current_local_frame; + if(local_frame_initialized[i]){ + current_local_frame = target_plane_local_frames.block(0, 2*i, 3, 2); + } + else{ + Eigen::JacobiSVD jSVD_normal(target_normal, Eigen::ComputeFullU); + current_local_frame = jSVD_normal.matrixU().block(0, 1, 3, 2); + target_plane_local_frames.block(0, 2*i, 3, 2) = current_local_frame; + local_frame_initialized[i] = true; + } + + Matrix32d local_coord = face_vtx_pos.transpose() * current_local_frame; + Eigen::JacobiSVD jSVD_coord(local_coord, Eigen::ComputeFullV); + Eigen::Vector2d fitting_line_direction = jSVD_coord.matrixV().col(0); + Eigen::Vector3d line_direction_3d = current_local_frame * fitting_line_direction; + target_pos = line_direction_3d * (line_direction_3d.transpose() * face_vtx_pos); + } + + B.block(3 * i, 0, 3, 3) = target_pos.transpose(); + } + } + + // Solver linear system + rhs = wX0; + rhs.noalias() += At_ * B; + if(!linear_solver_.solve(rhs, sol)){ + std::cerr << "Error: failed to solve mesh update system" << std::endl; + return false; + } + + const Matrix3X prev_vtx_pos = vtx_pos; + vtx_pos = sol.transpose(); + set_vertex_points(output_mesh, vtx_pos); + + if(param.mesh_update_disp_eps > 0.0) + { + const double rms_disp = std::sqrt((vtx_pos - prev_vtx_pos).colwise().squaredNorm().mean()); + if(rms_disp <= param.mesh_update_disp_eps) + { + if(print_progress_){ + std::cout << "Mesh update converged after " << (iter + 1) + << " iterations (rms displacement " << rms_disp << ")" << std::endl; + } + break; + } + } + } + + return true; + } + + + bool Poisson_mesh_update(const Matrix3X &target_normals, TriMesh &output_mesh) + { + output_mesh = mesh_; + + Matrix3Xi face_vtx_idx; + get_face_vertex_indices(output_mesh, face_vtx_idx); + + std::cout << "Starting Poisson mesh update......" << std::endl; + + Matrix3X vtx_pos; + get_vertex_points(output_mesh, vtx_pos); + + int n_faces = output_mesh.n_faces(); + int n_vtx = output_mesh.n_vertices(); + Eigen::VectorXd face_area_weights; + get_face_area_weights(output_mesh, face_area_weights); + + // Compute the initial centroid + Eigen::Vector3d initial_centroid = compute_centroid(face_vtx_idx, face_area_weights, vtx_pos); + + // Set up the linear least squares system + SparseMatrixXd A(3*n_faces + 1, n_vtx); + Eigen::MatrixX3d B(3*n_faces + 1, 3); + std::vector A_triplets(9 * n_faces + 1); + + // Set the target position of the first vertex at the origin + A_triplets.back() = Triplet(3*n_faces, 0, 1.0); + B.row(3*n_faces).setZero(); + + OMP_PARALLEL + { + OMP_FOR + for(int i = 0; i < n_faces; ++ i) + { + // Compute rotation from current face to the target plane + Eigen::Vector3d init_normal = to_eigen_vec3d(output_mesh.calc_face_normal(TriMesh::FaceHandle(i))); + Eigen::Vector3d target_normal = target_normals.col(i); + Eigen::Quaternion rot = Eigen::Quaternion::FromTwoVectors(init_normal, target_normal); + + Eigen::Vector3i vtx_idx = face_vtx_idx.col(i); + Eigen::Matrix3d face_vtx_pos; + for(int j = 0; j < 3; ++ j){ + face_vtx_pos.col(j) = vtx_pos.col(vtx_idx(j)); + } + + double area_weight = std::sqrt(face_area_weights(i)); + + for(int j = 0; j < 3; ++ j) + { + // Search for coefficients such that + // [v1 - (a * v2 + (1 - a) * v3)] * (v2 - v3) = 0 + // ==> a = [(v1 - v3) * (v2 - v3)] / ||v2 - v3||^2 + // Then the gradient coefficients become (1, -a, a - 1) + + int i1 = j, i2 = (j+1)%3, i3 = (j+2)%3; + + Eigen::Vector3d v1 = face_vtx_pos.col(i1), + v2 = face_vtx_pos.col(i2), + v3 = face_vtx_pos.col(i3); + + double a = 0.5; + if((v2 - v3).norm() > 1e-12){ + a = (v1 - v3).dot(v2 - v3) / (v2 - v3).squaredNorm(); + } + + // Compute gradient coefficient w.r.t vertex positions + Eigen::Vector3d current_grad_coef; + current_grad_coef(i1) = 1.0; + current_grad_coef(i2) = -a; + current_grad_coef(i3) = a - 1.0; + current_grad_coef *= area_weight; + + // Fill gradient coefficients into matrix + int triplet_addr = i*9 + j*3; + int row_idx = 3*i + j; + A_triplets[triplet_addr++] = Triplet(row_idx, vtx_idx(i1), current_grad_coef(i1)); + A_triplets[triplet_addr++] = Triplet(row_idx, vtx_idx(i2), current_grad_coef(i2)); + A_triplets[triplet_addr] = Triplet(row_idx, vtx_idx(i3), current_grad_coef(i3)); + + // Put the target gradient on the right-hand-side + Eigen::Vector3d target_grad = rot * (face_vtx_pos * current_grad_coef); + B.row(row_idx) = target_grad.transpose(); + } + } + } + + // Solver linear system + A.setFromTriplets(A_triplets.begin(), A_triplets.end()); + SparseMatrixXd At = A.transpose(); + SparseMatrixXd M = At * A; + Eigen::MatrixX3d rhs = At * B; + Eigen::MatrixX3d sol(n_vtx, 3); + + linear_solver_.reset_pattern(); + if(!(linear_solver_.compute(M) && linear_solver_.solve(rhs, sol))){ + std::cerr << "Error: failed to solve linear system for Poisson mesh update" << std::endl; + return false; + } + + // Align the new mesh with the intial mesh + vtx_pos = sol.transpose(); + Eigen::Vector3d new_centroid = compute_centroid(face_vtx_idx, area_weights_, vtx_pos); + vtx_pos.colwise() += initial_centroid - new_centroid; + set_vertex_points(output_mesh, vtx_pos); + + + return true; + } + + + + // Generate the matrix for mean-centering of the vertices of a triangle + void get_mean_centering_matrix(Eigen::Matrix3d &mat) + { + mat = Eigen::Matrix3d::Identity() - Eigen::Matrix3d::Constant(1.0/3); + } + + void get_mean_centered_face_vtx_pos(const Eigen::Matrix3Xd &vtx_pos, const Eigen::Vector3i &face_vtx, Eigen::Matrix3d &face_vtx_pos) + { + for(int i = 0; i < 3; ++ i){ + face_vtx_pos.col(i) = vtx_pos.col(face_vtx(i)); + } + + Eigen::Vector3d mean_pt = face_vtx_pos.rowwise().mean(); + face_vtx_pos.colwise() -= mean_pt; + } + + // Compute the centroid of a mesh given its vertex positions and face areas + Eigen::Vector3d compute_centroid(const Eigen::Matrix3Xi &face_vtx_idx, const Eigen::VectorXd &face_areas, const Eigen::Matrix3Xd &vtx_pos) + { + int n_faces = face_vtx_idx.cols(); + Eigen::Matrix3Xd face_centroids(3, n_faces); + + OMP_PARALLEL + { + OMP_FOR + for(int i = 0; i < n_faces; ++ i) + { + Eigen::Vector3d c = Eigen::Vector3d::Zero(); + Eigen::Vector3i face_vtx = face_vtx_idx.col(i); + + for(int j = 0; j < 3; ++ j){ + c += vtx_pos.col(face_vtx(j)); + } + + face_centroids.col(i) = c / 3.0; + } + } + + return (face_centroids * face_areas) / face_areas.sum(); + } + + + ////////// Methods for evaluating the quality of the updated mesh ///////////// + + // Compute the L2 norm between the initial mesh and filtered mesh + void show_normalized_mesh_displacement_norm(const TriMesh &filtered_mesh) + { + Eigen::Matrix3Xd init_vtx_pos, new_vtx_pos; + get_vertex_points(mesh_, init_vtx_pos); + get_vertex_points(filtered_mesh, new_vtx_pos); + Eigen::VectorXd vtx_disp_sqr_norm = (init_vtx_pos - new_vtx_pos).colwise().squaredNorm(); + + // Computer normalized vertex area weights from the original mesh + Eigen::VectorXd face_area_weights; + get_face_area_weights(mesh_, face_area_weights); + Eigen::Matrix3Xi face_vtx_indices; + get_face_vertex_indices(mesh_, face_vtx_indices); + int n_faces = mesh_.n_faces(); + + Eigen::VectorXd vtx_area(mesh_.n_vertices()); + vtx_area.setZero(); + for(int i = 0; i < n_faces; ++ i){ + for(int j = 0; j < 3; ++ j){ + vtx_area(face_vtx_indices(j, i)) += face_area_weights(i); + } + } + vtx_area /= vtx_area.sum(); + + std::cout << "Normalized mesh displacement norm: " << + std::sqrt(vtx_area.dot(vtx_disp_sqr_norm)) / average_edge_length(mesh_) << std::endl; + } + + + void show_error_statistics(const Eigen::VectorXd &err_values, double bin_size, int n_bins) + { + int n_elems = err_values.size(); + + Eigen::VectorXi error_bin_idx(n_elems); + OMP_PARALLEL + { + OMP_FOR + for(int i = 0; i < n_elems; ++ i){ + error_bin_idx(i) = std::min(n_bins, static_cast(std::floor(err_values(i) / bin_size))); + } + } + + Eigen::VectorXd bin_count(n_bins + 1); + bin_count.setZero(); + + for(int i = 0; i < n_elems; ++ i) + { + bin_count( error_bin_idx(i) ) += 1; + } + + bin_count /= bin_count.sum(); + + for(int i = 0; i < n_bins; ++ i) + { + double lower_val = bin_size * i; + double upper_val = bin_size * (i+1); + std::cout << lower_val << " to " << upper_val << ": " << bin_count(i) * 100 << "%" << std::endl; + } + + std::cout << "Over " << bin_size * n_bins << ": " << bin_count(n_bins) * 100 << "%" << std::endl; + } + + // Show statistics of the deviation between the new normals and target normals (in degrees) + void show_normal_error_statistics(const TriMesh &mesh, const Matrix3X &target_normals, int bin_size_in_degrees, int n_bins) + { + // Compute the normal deviation angle, and the number of flipped normals + int n_faces = mesh.n_faces(); + Eigen::VectorXd face_normal_error_angle(n_faces); + + for(int i = 0; i < n_faces; ++ i) + { + Eigen::Vector3d normal = to_eigen_vec3d(mesh.calc_face_normal(TriMesh::FaceHandle(i))); + double error_angle_cos = std::max(-1.0, std::min(1.0, normal.dot(target_normals.col(i)))); + face_normal_error_angle(i) = std::acos(error_angle_cos); + } + + face_normal_error_angle *= (180 / M_PI); + + std::cout << "Statistics of deviation between new normals and target normals:" << std::endl; + std::cout << "===============================================================" << std::endl; + show_error_statistics(face_normal_error_angle, bin_size_in_degrees, n_bins); + } +}; + +} + + + +#endif /* MESHNORMALFILTER_H_ */ diff --git a/src/MeshTypes.h b/src/MeshTypes.h new file mode 100644 index 0000000..da06ec3 --- /dev/null +++ b/src/MeshTypes.h @@ -0,0 +1,246 @@ +// BSD 3-Clause License +// +// Copyright (c) 2017, Bailin Deng +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// * Redistributions of source code must retain the above copyright notice, this +// list of conditions and the following disclaimer. +// +// * Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// * Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +// OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + +#ifndef MESHTYPES_H +#define MESHTYPES_H + +#include +#include +#include +#include "EigenTypes.h" + + +/// Default traits for triangle mesh +struct TriTraits : public OpenMesh::DefaultTraits +{ + /// Use double precision points + typedef OpenMesh::Vec3d Point; + /// Use double precision Normals + typedef OpenMesh::Vec3d Normal; + + /// Use RGBA Color + typedef OpenMesh::Vec4f Color; + +}; + +/// Simple name for triangle mesh +typedef OpenMesh::TriMesh_ArrayKernelT TriMesh; + + +namespace SDFilter +{ + +// Collect all point positions into a matrix. +// Argument points must be have been initialized with the correct number of columns. +template +void get_vertex_points(const MeshT &mesh, Matrix3X &points) +{ + points.resize(3, mesh.n_vertices()); + for(typename MeshT::ConstVertexIter cv_it = mesh.vertices_begin(); cv_it != mesh.vertices_end(); ++ cv_it) + { + points.col(cv_it->idx()) = to_eigen_vec3d(mesh.point(*cv_it)); + } +} + + +template +void set_vertex_points(MeshT &mesh, const Matrix3X &pos) +{ + assert(static_cast(mesh.n_vertices()) * 3 == static_cast(pos.size())); + + for(typename MeshT::ConstVertexIter cv_it = mesh.vertices_begin(); + cv_it != mesh.vertices_end(); ++ cv_it) + { + Eigen::Vector3d pt = pos.col(cv_it->idx()); + mesh.set_point(*cv_it, from_eigen_vec3d(pt)); + } +} + + +template +void get_vertex_points(const MeshT &mesh, std::vector &points) +{ + points.assign(3 * mesh.n_vertices(), 0.0); + for(typename MeshT::ConstVertexIter cv_it = mesh.vertices_begin(); cv_it != mesh.vertices_end(); ++ cv_it) + { + typename MeshT::Point pt = mesh.point(*cv_it); + int v_idx = cv_it->idx(); + + assert(v_idx >=0 && v_idx < static_cast(mesh.n_vertices())); + + for(int i = 0; i < 3; ++ i){ + points[3 * v_idx + i] = pt[i]; + } + } +} + +inline void get_face_vertex_indices(const TriMesh &mesh, Matrix3Xi &face_vertex_indices) +{ + face_vertex_indices.resize(3, mesh.n_faces()); + + for(TriMesh::ConstFaceIter cf_it = mesh.faces_begin(); cf_it != mesh.faces_end(); ++ cf_it) + { + int i = 0; + + for(TriMesh::ConstFaceVertexIter cfv_it = mesh.cfv_iter(*cf_it); cfv_it.is_valid(); ++ cfv_it) + { + face_vertex_indices(i++, cf_it->idx()) = cfv_it->idx(); + } + } +} + + +template +void set_vertex_points(MeshT &mesh, const std::vector &pos) +{ + assert(mesh.n_vertices() * 3 == pos.size()); + + for(typename MeshT::ConstVertexIter cv_it = mesh.vertices_begin(); + cv_it != mesh.vertices_end(); ++ cv_it) + { + int addr = cv_it->idx() * 3; + mesh.set_point(*cv_it, typename MeshT::Point(pos[addr], pos[addr+1], pos[addr+2])); + } +} + +template +inline Eigen::Vector3d bbox_dimension(const MeshT &mesh) +{ + Matrix3X vtx_pos; + get_vertex_points(mesh, vtx_pos); + + return vtx_pos.rowwise().maxCoeff() - vtx_pos.rowwise().minCoeff(); +} + + +template +inline double bbox_diag_length(const MeshT &mesh) +{ + return bbox_dimension(mesh).norm(); +} + +template +inline Eigen::Vector3d mesh_center(const MeshT &mesh) +{ + Matrix3X vtx_pos; + get_vertex_points(mesh, vtx_pos); + + return vtx_pos.rowwise().mean(); +} + + +template +inline double average_neighbor_face_centroid_dist(const MeshT &mesh) +{ + Matrix3X centroid_pos(3, mesh.n_faces()); + for(typename MeshT::ConstFaceIter cf_it = mesh.faces_begin(); cf_it != mesh.faces_end(); ++ cf_it) + { + centroid_pos.col(cf_it->idx()) = to_eigen_vec3d(mesh.calc_face_centroid(*cf_it)); + } + + double centroid_dist = 0; + int n = 0; + for(typename MeshT::ConstEdgeIter ce_it = mesh.edges_begin(); ce_it != mesh.edges_end(); ++ ce_it) + { + if(!mesh.is_boundary(*ce_it)) + { + int f1 = mesh.face_handle(mesh.halfedge_handle(*ce_it, 0)).idx(), + f2 = mesh.face_handle(mesh.halfedge_handle(*ce_it, 1)).idx(); + + centroid_dist += (centroid_pos.col(f1) - centroid_pos.col(f2)).norm(); + n++; + } + } + + return centroid_dist / n; + +} + +template +inline double average_edge_length(const MeshT &mesh) +{ + if(static_cast(mesh.n_edges()) == 0){ + return 0.0; + } + + double length = 0; + + for(typename MeshT::ConstEdgeIter ce_it = mesh.edges_begin(); ce_it != mesh.edges_end(); ++ ce_it) + { + length += mesh.calc_edge_length(*ce_it); + } + + return length / mesh.n_edges(); +} + + +// Center the mesh at the origin, and normalize its scale +// Return the orginal mesh center, and the orginal scale +template +inline void normalize_mesh(MeshT &mesh, Eigen::Vector3d &original_center, double &original_scale) +{ + original_center = mesh_center(mesh); + original_scale = bbox_diag_length(mesh); + + Matrix3X vtx_pos; + get_vertex_points(mesh, vtx_pos); + + vtx_pos.colwise() -= original_center; + vtx_pos *= (1.0 / original_scale); + + set_vertex_points(mesh, vtx_pos); +} + +// Scale and move a normalized mesh to according to its original scale and center +template +inline void restore_mesh(MeshT &mesh, const Eigen::Vector3d &original_center, double original_scale) +{ + Matrix3X vtx_pos; + get_vertex_points(mesh, vtx_pos); + + Eigen::Vector3d center = vtx_pos.rowwise().mean(); + vtx_pos.colwise() -= center; + + vtx_pos *= original_scale; + vtx_pos.colwise() += original_center; + + set_vertex_points(mesh, vtx_pos); +} + +template +inline bool write_mesh(const MeshT &mesh, const std::string &filename, int precision = 16) +{ + return OpenMesh::IO::write_mesh(mesh, filename, OpenMesh::IO::Options::Default, precision); +} + +} + +#endif // MESHTYPES_H diff --git a/src/SDFilter.h b/src/SDFilter.h new file mode 100644 index 0000000..3e24400 --- /dev/null +++ b/src/SDFilter.h @@ -0,0 +1,848 @@ +// BSD 3-Clause License +// +// Copyright (c) 2017, Bailin Deng +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// * Redistributions of source code must retain the above copyright notice, this +// list of conditions and the following disclaimer. +// +// * Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// * Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +// OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + +#ifndef ITERATIVESDFILTER_H_ +#define ITERATIVESDFILTER_H_ + +#include "EigenTypes.h" +#ifdef USE_CHOLMOD +#include +#endif +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef USE_OPENMP +#include +#ifdef USE_MSVC +#define OMP_PARALLEL __pragma(omp parallel) +#define OMP_FOR __pragma(omp for) +#define OMP_SINGLE __pragma(omp single) +#else +#define OMP_PARALLEL _Pragma("omp parallel") +#define OMP_FOR _Pragma("omp for") +#define OMP_SINGLE _Pragma("omp single") +#endif +#else +#include +#define OMP_PARALLEL +#define OMP_FOR +#define OMP_SINGLE +#endif + +namespace SDFilter +{ + +#ifdef USE_CHOLMOD +static constexpr int DEFAULT_LINEAR_SOLVER_TYPE = 2; +#else +static constexpr int DEFAULT_LINEAR_SOLVER_TYPE = 1; +#endif + +struct RunStatistics +{ + double preprocessing_secs = 0.0; + double filtering_secs = 0.0; + double mesh_update_secs = 0.0; + double mesh_filter_total_secs = 0.0; + double denoise_total_secs = 0.0; + int solver_iterations = 0; + bool solver_converged = false; + int outer_iterations = 0; +}; + +class Parameters +{ +public: + Parameters() + :lambda(10), eta(1.0), mu(1.0), nu(1.0), max_iter(100), avg_disp_eps(1e-6), + normalize_iterates(true), linear_solver_type(static_cast(DEFAULT_LINEAR_SOLVER_TYPE)), linear_solver_max_iterations(2000), + linear_solver_tolerance(1e-8) {} + + virtual ~Parameters() {} + + enum LinearSolverType + { + CG, + LDLT, + CHOLMOD, + LINEAR_SOLVER_TYPE_COUNT + }; + + double lambda; // Regularization weight + double eta; // Gaussian standard deviation for spatial weight, relative to the bounding box diagonal length + double mu; // Gaussian standard deviation for guidance weight + double nu; // Gaussian standard deviation for signal weight + + // Parameters related to termination criteria + int max_iter; // Max number of iterations + double avg_disp_eps; // Max average per-signal displacement threshold between two iterations for determining convergence + bool deterministic_mode = false; // Force deterministic execution mode (single-threaded OpenMP/Eigen) + bool normalize_iterates; // Normalization of the filtered normals in each iteration + LinearSolverType linear_solver_type; // Linear solver backend for sparse solves + int linear_solver_max_iterations; // Iteration cap for iterative sparse solvers + double linear_solver_tolerance; // Tolerance for iterative sparse solvers + + + // Load options from file + bool load(const char* filename) + { + std::ifstream ifile(filename); + if (!ifile.is_open()){ + std::cerr << "Error while opening file " << filename << std::endl; + return false; + } + + std::string line; + while(std::getline(ifile, line)) + { + std::string::size_type pos = line.find_first_not_of(' '); + if(pos == std::string::npos){ + continue; + } + + // Check for comment line + else if(line.at(pos) == '#'){ + continue; + } + + std::string::size_type end_pos = line.find_first_of(' '); + std::string option_str = line.substr(pos, end_pos - pos); + std::string value_str = line.substr(end_pos + 1, std::string::npos); + OptionInterpreter opt(option_str, value_str); + + load_option(opt); + } + + std::cout << "Successfully loaded options from file " << filename << std::endl; + + return true; + } + + + void output() + { + std::cout << std::endl; + std::cout << "====== Filter parameters =========" << std::endl; + output_options(); + std::cout << "==================================" << std::endl; + std::cout << std::endl; + } + + // Check whether the parameter values are valid + virtual bool valid_parameters() const + { + if(lambda <= 0.0){ + std::cerr << "Error: Lambda must be positive" << std::endl; + return false; + } + + if(eta <= 0.0){ + std::cerr << "Error: Eta must be positive" << std::endl; + return false; + } + + if(mu <= 0.0){ + std::cerr << "Error: Mu must be positive" << std::endl; + return false; + } + + if(nu <= 0.0){ + std::cerr << "Error: Nu must be positive" << std::endl; + return false; + } + + if(max_iter < 1){ + std::cerr << "Error: MaxIterations must be at least 1" << std::endl; + return false; + } + + if(avg_disp_eps <= 0.0){ + std::cerr << "Error: average displacement threshold must be positive" << std::endl; + return false; + } + + if(linear_solver_type < CG || linear_solver_type >= LINEAR_SOLVER_TYPE_COUNT){ + std::cerr << "Error: invalid LinearSolverType" << std::endl; + return false; + } + + if(linear_solver_max_iterations <= 0){ + std::cerr << "Error: LinearSolverMaxIterations must be positive" << std::endl; + return false; + } + + if(linear_solver_tolerance <= 0.0 || !std::isfinite(linear_solver_tolerance)){ + std::cerr << "Error: LinearSolverTolerance must be finite and positive" << std::endl; + return false; + } + + return true; + } + +protected: + + class OptionInterpreter + { + public: + OptionInterpreter(const std::string &option_str, const std::string &value_str) + :option_str_(option_str), value_str_(value_str){} + + template + bool load(const std::string &target_option_name, T &target_option_value) const + { + if(option_str_ == target_option_name){ + if(!load_value(value_str_, target_option_value)){ + std::cerr << "Error loading option: " << target_option_name << std::endl; + return false; + } + + return true; + } + else{ + return false; + } + } + + template + bool load_enum(const std::string &target_option_name, EnumT enum_value_count, EnumT &value) const + { + if(option_str_ == target_option_name) + { + int enum_int = 0; + if(load_value(value_str_, enum_int)) + { + if(enum_int >= 0 && enum_int < static_cast(enum_value_count)){ + value = static_cast(enum_int); + return true; + } + } + + std::cerr << "Error loading option: " << target_option_name << std::endl; + return false; + } + else{ + return false; + } + } + + private: + std::string option_str_, value_str_; + + bool load_value(const std::string &str, double &value) const + { + try{ + value = std::stod(str); + } + catch (const std::invalid_argument& ia){ + std::cerr << "Invalid argument: " << ia.what() << std::endl; + return false; + } + catch (const std::out_of_range &oor){ + std::cerr << "Out of Range error: " << oor.what() << std::endl; + return false; + } + + return true; + } + + + bool load_value(const std::string &str, int &value) const + { + try{ + value = std::stoi(str); + } + catch (const std::invalid_argument& ia){ + std::cerr << "Invalid argument: " << ia.what() << std::endl; + return false; + } + catch (const std::out_of_range &oor){ + std::cerr << "Out of Range error: " << oor.what() << std::endl; + return false; + } + + return true; + } + + bool load_value(const std::string &str, bool &value) const + { + int bool_value = 0; + if(load_value(str, bool_value)){ + value = (bool_value != 0); + return true; + } + else{ + return false; + } + } + }; + + virtual bool load_option(const OptionInterpreter &opt) + { + return opt.load("Lambda", lambda) || + opt.load("Eta", eta) || + opt.load("Mu", mu) || + opt.load("Nu", nu) || + opt.load("MaxFilterIterations", max_iter) || + opt.load("DeterministicMode", deterministic_mode) || + opt.load_enum("LinearSolverType", LINEAR_SOLVER_TYPE_COUNT, linear_solver_type) || + opt.load("LinearSolverMaxIterations", linear_solver_max_iterations) || + opt.load("LinearSolverTolerance", linear_solver_tolerance); + } + + virtual void output_options() + { + std::cout << "Lambda: " << lambda << std::endl; + std::cout << "Eta:" << eta << std::endl; + std::cout << "Mu:" << mu << std::endl; + std::cout << "Nu:" << nu << std::endl; + std::cout << "Deterministic mode: " << (deterministic_mode ? 1 : 0) << std::endl; + std::cout << "Linear solver type: " << static_cast(linear_solver_type) << std::endl; + std::cout << "Linear solver max iterations: " << linear_solver_max_iterations << std::endl; + std::cout << "Linear solver tolerance: " << linear_solver_tolerance << std::endl; + } + + +}; + + +class Timer +{ +public: + + typedef int EventID; + + EventID get_time() + { + EventID id = time_values_.size(); + + #ifdef USE_OPENMP + time_values_.push_back(omp_get_wtime()); + #else + time_values_.push_back(clock()); + #endif + + return id; + } + + double elapsed_time(EventID event1, EventID event2) + { + assert(event1 >= 0 && event1 < static_cast(time_values_.size())); + assert(event2 >= 0 && event2 < static_cast(time_values_.size())); + + #ifdef USE_OPENMP + return time_values_[event2] - time_values_[event1]; + #else + return double(time_values_[event2] - time_values_[event1]) / CLOCKS_PER_SEC; + #endif + } + +private: + #ifdef USE_OPENMP + std::vector time_values_; + #else + std::vector time_values_; + #endif +}; + + +class SDFilter +{ +public: + + SDFilter() + :signal_dim_(-1), signal_count_(-1), print_progress_(true), + print_timing_(true), print_diagnostic_info_(false){} + + virtual ~SDFilter(){} + + const RunStatistics& run_stats() const + { + return stats_; + } + +protected: + + bool filter(Parameters param) + { + stats_.preprocessing_secs = 0.0; + stats_.filtering_secs = 0.0; + stats_.solver_iterations = 0; + stats_.solver_converged = false; + + std::cout << "Preprocessing......" << std::endl; + + Timer timer; + Timer::EventID begin_time = timer.get_time(); + + if(!initialize_filter(param)) + { + std::cerr << "Error: unable to initialize filter" << std::endl; + return false; + } + + Timer::EventID preprocess_end_time = timer.get_time(); + + std::cout << "Filtering......" << std::endl; + + fixedpoint_solver(param); + + Timer::EventID filter_end_time = timer.get_time(); + stats_.preprocessing_secs = timer.elapsed_time(begin_time, preprocess_end_time); + stats_.filtering_secs = timer.elapsed_time(preprocess_end_time, filter_end_time); + + if(print_timing_){ + std::cout << "Preprocessing timing: " << stats_.preprocessing_secs << " secs" << std::endl; + std::cout << "Filtering timing: " << stats_.filtering_secs << " secs" << std::endl; + } + + return true; + } + + + void fixedpoint_solver(const Parameters ¶m) + { + // Store signals in the previous iteration + Eigen::MatrixXd init_signals = signals_; + Eigen::MatrixXd prev_signals; + + // Weighted initial signals, as used in the fixed-point solver + Eigen::MatrixXd weighted_init_signals = init_signals * (area_weights_ * (2 * param.nu * param.nu / param.lambda)).asDiagonal(); + + Eigen::MatrixXd filtered_signals; + double h = -0.5 / (param.nu * param.nu); + + Eigen::Index n_neighbor_pairs = neighboring_pairs_.cols(); + + // The weights for neighboring pairs that are used for convex combination of neighboring signals in the fixed-point solver + Eigen::VectorXd neighbor_pair_weights(n_neighbor_pairs); + + // Compute the termination threshold for area weighted squread norm of signal change between two iterations + double disp_sqr_norm_threshold = area_weights_.sum() * param.avg_disp_eps * param.avg_disp_eps; + + int output_frequency = 10; + + for(int num_iter = 1; num_iter <= param.max_iter; ++ num_iter) + { + prev_signals = signals_; + filtered_signals = weighted_init_signals; + + OMP_PARALLEL + { + OMP_FOR + for(Eigen::Index i = 0; i < n_neighbor_pairs; ++ i) + { + int idx1 = neighboring_pairs_(0, i), idx2 = neighboring_pairs_(1, i); + + neighbor_pair_weights(i) = + precomputed_area_spatial_guidance_weights_(i) * + std::exp( h * (signals_.col(idx1) - signals_.col(idx2)).squaredNorm()); + } + + OMP_FOR + for(int i = 0; i < signal_count_; ++ i) + { + Eigen::Index neighbor_info_start_idx = neighborhood_info_boundaries_(i); + Eigen::Index neighbor_info_end_idx = neighborhood_info_boundaries_(i+1); + + for(Eigen::Index j = neighbor_info_start_idx; j < neighbor_info_end_idx; ++ j) + { + Eigen::Index neighbor_idx = neighborhood_info_(0, j); + Eigen::Index coef_idx = neighborhood_info_(1, j); + + filtered_signals.col(i) += signals_.col(neighbor_idx) * neighbor_pair_weights(coef_idx); + } + + if(param.normalize_iterates){ + filtered_signals.col(i).normalize(); + } + else{ + filtered_signals.col(i) /= filtered_signals(signal_dim_, i); + } + } + } + + signals_ = filtered_signals; + + double var_disp_sqrnorm = area_weights_.dot((signals_ - prev_signals).colwise().squaredNorm()); + + if(print_diagnostic_info_){ + std::cout << "Iteration " << num_iter << ", Target function value " << target_function(param, init_signals) << std::endl; + } + else if(print_progress_ && num_iter % output_frequency == 0){ + std::cout << "Iteration "<< num_iter << "..." << std::endl; + } + + + if(var_disp_sqrnorm <= disp_sqr_norm_threshold){ + stats_.solver_iterations = num_iter; + stats_.solver_converged = true; + std::cout << "Solver converged after " << num_iter << " iterations" << std::endl; + break; + } + else if(num_iter == param.max_iter){ + stats_.solver_iterations = param.max_iter; + stats_.solver_converged = false; + std::cout << "Solver terminated after " << param.max_iter << " iterations" << std::endl; + break; + } + } + } + + + // Linear solver for symmetric positive definite matrix, + class LinearSolver + { + public: + LinearSolver(Parameters::LinearSolverType solver_type) + :solver_type_(solver_type), solver_max_iterations_(2000), solver_tolerance_(1e-8), pattern_analyzed(false){} + + // Initialize the solver with matrix + bool compute(const SparseMatrixXd &A) + { + cached_matrix_ = A; + if(solver_type_ == Parameters::LDLT) + { + if(!pattern_analyzed) + { + LDLT_solver_.analyzePattern(A); + if(!check_error(LDLT_solver_, "Cholesky analyzePattern failed")){ + return false; + } + + pattern_analyzed = true; + } + + LDLT_solver_.factorize(A); + return check_error(LDLT_solver_, "Cholesky factorization failed"); + } + else if(solver_type_ == Parameters::CHOLMOD) + { +#ifdef USE_CHOLMOD + if(!pattern_analyzed) + { + cholmod_solver_.analyzePattern(A); + if(!check_error(cholmod_solver_, "CHOLMOD analyzePattern failed")){ + return false; + } + pattern_analyzed = true; + } + cholmod_solver_.factorize(A); + return check_error(cholmod_solver_, "CHOLMOD factorization failed"); +#else + std::cerr << "CHOLMOD solver requested but not available in this build. Falling back to LDLT." << std::endl; + solver_type_ = Parameters::LDLT; + reset_pattern(); + return compute(A); +#endif + } + else if(solver_type_ == Parameters::CG){ + CG_solver_.setMaxIterations(solver_max_iterations_); + CG_solver_.setTolerance(solver_tolerance_); + CG_solver_.compute(A); + return check_error(CG_solver_, "CG solver compute failed"); + } + else{ + return false; + } + } + + template + bool solve(const MatrixT &rhs, MatrixT &sol) + { + if(solver_type_ == Parameters::LDLT || solver_type_ == Parameters::CHOLMOD) + { + return solve_cpu_direct(rhs, sol, nullptr); + } + else if(solver_type_ == Parameters::CG) + { + sol = CG_solver_.solve(rhs); + return check_error(CG_solver_, "CG solve failed"); + } + else{ + return false; + } + } + + void reset_pattern() + { + pattern_analyzed = false; + } + + void set_solver_type(Parameters::LinearSolverType type) + { + solver_type_ = type; +#ifndef USE_CHOLMOD + if(solver_type_ == Parameters::CHOLMOD){ + solver_type_ = Parameters::LDLT; + } +#endif + if(solver_type_ != Parameters::CG){ + reset_pattern(); + } + } + + void set_solver_params(int max_iterations, double tolerance) + { + solver_max_iterations_ = std::max(1, max_iterations); + solver_tolerance_ = std::max(tolerance, std::numeric_limits::epsilon()); + } + + private: + Parameters::LinearSolverType solver_type_; + Eigen::SimplicialLDLT LDLT_solver_; +#ifdef USE_CHOLMOD + Eigen::CholmodSupernodalLLT cholmod_solver_; +#endif + Eigen::ConjugateGradient > CG_solver_; + int solver_max_iterations_; + double solver_tolerance_; + SparseMatrixXd cached_matrix_; + + bool pattern_analyzed; // Flag for symbolic factorization + + template + bool solve_cpu_direct(const MatrixT &rhs, MatrixT &sol, double *out_secs) + { + const auto begin = std::chrono::steady_clock::now(); + bool ok = false; + if(solver_type_ == Parameters::LDLT) + { +#ifdef USE_OPENMP + const int n_cols = static_cast(rhs.cols()); + sol.resize(rhs.rows(), rhs.cols()); + OMP_PARALLEL + { + OMP_FOR + for(int i = 0; i < n_cols; ++ i){ + sol.col(i) = LDLT_solver_.solve(rhs.col(i)); + } + } +#else + sol = LDLT_solver_.solve(rhs); +#endif + ok = check_error(LDLT_solver_, "LDLT solve failed"); + } + else if(solver_type_ == Parameters::CHOLMOD) + { +#ifdef USE_CHOLMOD + sol = cholmod_solver_.solve(rhs); + ok = check_error(cholmod_solver_, "CHOLMOD solve failed"); +#else + ok = false; +#endif + } + + if(out_secs){ + const auto end = std::chrono::steady_clock::now(); + *out_secs = std::chrono::duration(end - begin).count(); + } + return ok; + } + + template + bool check_error(const SolverT &solver, const std::string &error_message){ + if(solver.info() != Eigen::Success){ + std::cerr << error_message << std::endl; + } + + return solver.info() == Eigen::Success; + } + }; + + + +protected: + + int signal_dim_; // Dimension of the signals + int signal_count_; // Number of signals + + Eigen::MatrixXd signals_; // Signals to be filtered. Represented in homogeneous form when there is no normalization constraint + Eigen::VectorXd area_weights_; // Area weights for each element + + Eigen::Matrix2Xi neighboring_pairs_; // Each column stores the indices for a pair of neighboring elements + Eigen::VectorXd precomputed_area_spatial_guidance_weights_; // Precomputed weights (area, spatial Gaussian and guidance Gaussian) for neighboring pairs + + // The neighborhood information for each signal element is stored as contiguous columns within the neighborhood_info_ matrix + // For each column, the first element is the index of a neighboring element, the second one is the corresponding address within array neighboring_pairs_ + Matrix2XIdx neighborhood_info_; + VectorXIdx neighborhood_info_boundaries_; // Boundary positions for the neighborhood information segments + + bool print_progress_; + bool print_timing_; + bool print_diagnostic_info_; + RunStatistics stats_; + + // Overwrite this in a subclass to provide the initial spatial positions, guidance, and signals. + virtual void get_initial_data(Eigen::MatrixXd &guidance, Eigen::MatrixXd &init_signals, Eigen::VectorXd &area_weights) = 0; + + bool initialize_filter(Parameters ¶m) + { + // Retrive input signals and their area weights + Eigen::MatrixXd guidance, init_signals; + get_initial_data(guidance, init_signals, area_weights_); + + signal_dim_ = init_signals.rows(); + signal_count_ = init_signals.cols(); + if(signal_count_ <= 0){ + return false; + } + + if(param.normalize_iterates){ + signals_ = init_signals; + } + else{ + signals_.resize(signal_dim_ + 1, signal_count_); + signals_.block(0, 0, signal_dim_, signal_count_) = init_signals; + signals_.row(signal_dim_).setOnes(); + } + + Eigen::VectorXd neighbor_dists; + if(!get_neighborhood(param, neighboring_pairs_, neighbor_dists)){ + std::cerr << "Unable to get neighborhood information, no filtering done..." << std::endl; + return false; + } + + // Pre-compute filtering weights, and rescale the lambda parameter + + Eigen::Index n_neighbor_pairs = neighboring_pairs_.cols(); + if(n_neighbor_pairs <= 0){ + return false; + } + + precomputed_area_spatial_guidance_weights_.resize(n_neighbor_pairs); + double h_spatial = - 0.5 / (param.eta * param.eta); + double h_guidance = - 0.5 / (param.mu * param.mu); + Eigen::VectorXd area_spatial_weights(n_neighbor_pairs); // Area-integrated spatial weights, used for rescaling lambda + + + OMP_PARALLEL + { + OMP_FOR + for(Eigen::Index i = 0; i < n_neighbor_pairs; ++ i) + { + // Read the indices of a neighboring pair, and their distance + int idx1 = neighboring_pairs_(0, i), idx2 = neighboring_pairs_(1, i); + double d = neighbor_dists(i); + + // Compute the weights associated with the pair + area_spatial_weights(i) = (area_weights_(idx1) + area_weights_(idx2)) * std::exp( h_spatial * d * d ); + precomputed_area_spatial_guidance_weights_(i) = (area_weights_(idx1) + area_weights_(idx2)) * + std::exp( h_guidance * (guidance.col(idx1) - guidance.col(idx2)).squaredNorm() + h_spatial * d * d ); + } + + } + + + assert(neighbor_dists.size() > 0); + param.lambda *= ( area_weights_.sum() / area_spatial_weights.sum() ); // Rescale lambda to make regularization and fidelity terms comparable + + // Pre-compute neighborhood_info_ in CSR-like contiguous arrays. + std::vector degree(static_cast(signal_count_), Eigen::Index(0)); + for(Eigen::Index i = 0; i < n_neighbor_pairs; ++i) + { + const int idx1 = neighboring_pairs_(0, i), idx2 = neighboring_pairs_(1, i); + ++degree[static_cast(idx1)]; + ++degree[static_cast(idx2)]; + } + + neighborhood_info_boundaries_.resize(signal_count_ + 1); + neighborhood_info_boundaries_(0) = 0; + for(int i = 0; i < signal_count_; ++i){ + neighborhood_info_boundaries_(i + 1) = neighborhood_info_boundaries_(i) + degree[static_cast(i)]; + } + + neighborhood_info_.resize(2, 2 * n_neighbor_pairs); + std::vector cursor(static_cast(signal_count_)); + for(int i = 0; i < signal_count_; ++i){ + cursor[static_cast(i)] = neighborhood_info_boundaries_(i); + } + + for(Eigen::Index i = 0; i < n_neighbor_pairs; ++i) + { + const int idx1 = neighboring_pairs_(0, i), idx2 = neighboring_pairs_(1, i); + + Eigen::Index pos1 = cursor[static_cast(idx1)]++; + neighborhood_info_(0, pos1) = idx2; + neighborhood_info_(1, pos1) = i; + + Eigen::Index pos2 = cursor[static_cast(idx2)]++; + neighborhood_info_(0, pos2) = idx1; + neighborhood_info_(1, pos2) = i; + } + + return true; + } + + // Find out all neighboring paris, as well as their distance + virtual bool get_neighborhood(const Parameters ¶m, Eigen::Matrix2Xi &neighbor_pairs, Eigen::VectorXd &neighbor_dist) = 0; + + double target_function(const Parameters ¶m, const Eigen::MatrixXd &init_signals) + { + // Compute regularizer term, using the contribution from each neighbor pair + Eigen::Index n_neighbor_pairs = neighboring_pairs_.cols(); + Eigen::VectorXd pair_values(n_neighbor_pairs); + pair_values.setZero(); + double h = - 0.5 / (param.nu * param.nu); + + OMP_PARALLEL + { + OMP_FOR + for(Eigen::Index i = 0; i < n_neighbor_pairs; ++ i) + { + int idx1 = neighboring_pairs_(0, i), idx2 = neighboring_pairs_(1, i); + pair_values[i] = precomputed_area_spatial_guidance_weights_(i) + * std::max(0.0, 1.0 - std::exp( h * (signals_.col(idx1) - signals_.col(idx2)).squaredNorm())); + } + } + + double reg = pair_values.sum(); + + // Compute the fidelity term, which is the squared difference between current and initial signals, weighted by the areas + double fid = area_weights_.dot((signals_ - init_signals).colwise().squaredNorm()); + + return fid + reg * param.lambda; + } +}; + +} + + + + +#endif /* ITERATIVESDFILTER_H_ */ From af7ba1e28895738d5f21a1d2747c9b42180c3aba Mon Sep 17 00:00:00 2001 From: Markus Lanxinger Date: Sat, 21 Mar 2026 09:28:12 +0100 Subject: [PATCH 2/3] Fix Eigen3 detection with fallback for CI environments find_package(Eigen3 3.3) can fail on some CI runners even when Eigen is installed. Add cascading fallbacks: retry without version, probe common system include paths (/usr/include/eigen3, /usr/local/include/eigen3, /opt/homebrew/include/eigen3), then try default search paths. Co-Authored-By: Claude Opus 4.6 (1M context) --- CMakeLists.txt | 22 +++++++++++++++++++--- 1 file changed, 19 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 7b3ebb0..ceda257 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,14 +16,30 @@ if(ENABLE_OPENMP) find_package(OpenMP QUIET) endif() -# Eigen (required) +# Eigen (required, header-only) find_package(Eigen3 3.3 QUIET) if(NOT TARGET Eigen3::Eigen) - message(STATUS "Eigen3 not found via find_package, trying EIGEN3_INCLUDE_DIR") + # Retry without version constraint + find_package(Eigen3 QUIET) +endif() +if(NOT TARGET Eigen3::Eigen) + # Try user-provided include dir if(DEFINED EIGEN3_INCLUDE_DIR AND EXISTS "${EIGEN3_INCLUDE_DIR}") + set(_eigen_dir "${EIGEN3_INCLUDE_DIR}") + else() + # Probe common system locations + find_path(_eigen_dir NAMES Eigen/Dense + PATHS /usr/include/eigen3 /usr/local/include/eigen3 /opt/homebrew/include/eigen3 + NO_DEFAULT_PATH) + if(NOT _eigen_dir) + find_path(_eigen_dir NAMES Eigen/Dense) + endif() + endif() + if(_eigen_dir) + message(STATUS "Using Eigen from ${_eigen_dir}") add_library(Eigen3::Eigen INTERFACE IMPORTED) set_target_properties(Eigen3::Eigen PROPERTIES - INTERFACE_INCLUDE_DIRECTORIES "${EIGEN3_INCLUDE_DIR}") + INTERFACE_INCLUDE_DIRECTORIES "${_eigen_dir}") else() message(FATAL_ERROR "Eigen3 not found. Install via system package manager or set -DEIGEN3_INCLUDE_DIR=...") endif() From 84d533db344d82bf5f893e967f98aeff3370fd05 Mon Sep 17 00:00:00 2001 From: Markus Lanxinger Date: Sat, 21 Mar 2026 09:57:00 +0100 Subject: [PATCH 3/3] Bump GitHub Actions to Node.js 24 compatible versions MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - actions/checkout v4 → v5 - actions/upload-artifact v4 → v5 - softprops/action-gh-release v1 → v2 Fixes Node.js 20 deprecation warnings. Co-Authored-By: Claude Opus 4.6 (1M context) --- .github/workflows/build.yml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index b4baefc..47c6472 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -31,7 +31,7 @@ jobs: steps: - name: Checkout code - uses: actions/checkout@v4 + uses: actions/checkout@v5 - name: Install dependencies (Ubuntu) if: runner.os == 'Linux' @@ -113,7 +113,7 @@ jobs: copy README.md artifacts\ - name: Upload artifacts - uses: actions/upload-artifact@v4 + uses: actions/upload-artifact@v5 with: name: ${{ matrix.artifact_name }} path: artifacts/ @@ -133,14 +133,14 @@ jobs: - name: Upload to release (Linux/macOS) if: github.event_name == 'release' && runner.os != 'Windows' - uses: softprops/action-gh-release@v1 + uses: softprops/action-gh-release@v2 with: files: ${{ matrix.artifact_name }}.tar.gz token: ${{ secrets.GITHUB_TOKEN }} - name: Upload to release (Windows) if: github.event_name == 'release' && runner.os == 'Windows' - uses: softprops/action-gh-release@v1 + uses: softprops/action-gh-release@v2 with: files: ${{ matrix.artifact_name }}.zip token: ${{ secrets.GITHUB_TOKEN }}