diff --git a/.gitmodules b/.gitmodules deleted file mode 100644 index b01fc8bfc0..0000000000 --- a/.gitmodules +++ /dev/null @@ -1,3 +0,0 @@ -[submodule "src/dependencies/MeshSDFilter"] - path = src/dependencies/MeshSDFilter - url = https://github.com/alicevision/MeshSDFilter.git diff --git a/src/dependencies/MeshSDFilter b/src/dependencies/MeshSDFilter deleted file mode 160000 index f63adc897d..0000000000 --- a/src/dependencies/MeshSDFilter +++ /dev/null @@ -1 +0,0 @@ -Subproject commit f63adc897d7e2ea01b49554846cc1ca7ac8784b1 diff --git a/src/dependencies/MeshSDFilter/CMakeLists.txt b/src/dependencies/MeshSDFilter/CMakeLists.txt new file mode 100644 index 0000000000..b43a972041 --- /dev/null +++ b/src/dependencies/MeshSDFilter/CMakeLists.txt @@ -0,0 +1,94 @@ +cmake_minimum_required(VERSION 3.1) + +project(SDFilter) + +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +set(CMAKE_CXX_EXTENSIONS OFF) +set(CMAKE_VERBOSE_MAKEFILE ON) + + +# Additional compiler flags +if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang") + message("Clang compiler found.") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}") +elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "AppleClang") + message("AppleClang compiler found.") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}") +elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU") + message("GNU compiler found.") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}") +elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "MSVC") + message("MSVC compiler found.") + add_definitions(/DUSE_MSVC) + add_definitions(/D_USE_MATH_DEFINES) +endif() + +# Detect OpenMP environment +set(OPENMP ON CACHE BOOL "OpenMP") +if(OPENMP) + find_package(OpenMP QUIET) + if(OPENMP_FOUND) + message("OpenMP found. OpenMP activated in release.") + add_definitions(-DUSE_OPENMP) + + else() + message("OpenMP not found.") + endif() +endif() + + +# Detect Eigen3 +if(NOT EIGEN3_FOUND) + find_package(Eigen3 REQUIRED) + if(EIGEN3_FOUND) + message("Found external Eigen. include: ${EIGEN3_INCLUDE_DIR}, version: ${EIGEN3_VERSION_STRING}.") + endif() +endif() + +find_package(OpenMesh REQUIRED) + +add_library(MeshSDLibrary INTERFACE) +target_include_directories(MeshSDLibrary INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}) +target_link_libraries(MeshSDLibrary INTERFACE Eigen3::Eigen OpenMeshCore) + + +# Executable for filtering +add_executable(MeshSDFilter + EigenTypes.h + MeshTypes.h + SDFilter.h + MeshNormalFilter.h + MeshSDFilter.cpp +) +target_link_libraries(MeshSDFilter MeshSDLibrary) + + +# Executable for denoising +add_executable(MeshDenoiser + EigenTypes.h + MeshTypes.h + SDFilter.h + MeshNormalFilter.h + MeshNormalDenoising.h + MeshDenoiser.cpp +) +target_link_libraries(MeshDenoiser MeshSDLibrary) + + +if(OPENMP_FOUND) + #target_compile_options(MeshSDLibrary PUBLIC "$<$:${OpenMP_CXX_FLAGS}>") + #target_compile_definitions(MeshSDLibrary PUBLIC "$<$:USE_OPENMP>") + #target_link_libraries(MeshSDLibrary "$<$:${OpenMP_CXX_FLAGS}>") + + target_compile_options(MeshSDFilter PUBLIC "$<$:${OpenMP_CXX_FLAGS}>") + target_compile_definitions(MeshSDFilter PUBLIC "$<$:USE_OPENMP>") + target_link_libraries(MeshSDFilter "$<$:${OpenMP_CXX_FLAGS}>") + + target_compile_options(MeshDenoiser PUBLIC "$<$:${OpenMP_CXX_FLAGS}>") + target_compile_definitions(MeshDenoiser PUBLIC "$<$:USE_OPENMP>") + target_link_libraries(MeshDenoiser "$<$:${OpenMP_CXX_FLAGS}>") +endif() diff --git a/src/dependencies/MeshSDFilter/DenoisingOptions.txt b/src/dependencies/MeshSDFilter/DenoisingOptions.txt new file mode 100644 index 0000000000..df845a5ea6 --- /dev/null +++ b/src/dependencies/MeshSDFilter/DenoisingOptions.txt @@ -0,0 +1,23 @@ +#### Option file for SD filter based mesh denoising +#### Lines starting with '#' are comments + +## Regularization weight, must be positive. +Lambda 2 + +## Gaussian standard deviation for spatial weight, scaled by the average distance between adjacent face cetroids. Must be positive. +Eta 1.5 + +## Gaussian standard deviation for guidance weight, must be positive. +Mu 1.5 + +## Gaussian standard deviation for signal weight, must be positive +Nu 0.3 + +## Closeness weight for mesh update, must be positive +MeshUpdateClosenessWeight 0.001 + +## Iterations for mesh update, must be positive +MeshUpdateIterations 20 + +## Outer iteration for denoising, must be positive integers +OuterIterations 5 \ No newline at end of file diff --git a/src/dependencies/MeshSDFilter/EigenTypes.h b/src/dependencies/MeshSDFilter/EigenTypes.h new file mode 100644 index 0000000000..b9392913f4 --- /dev/null +++ b/src/dependencies/MeshSDFilter/EigenTypes.h @@ -0,0 +1,75 @@ +// 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 + + +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/src/dependencies/MeshSDFilter/FilteringOptions.txt b/src/dependencies/MeshSDFilter/FilteringOptions.txt new file mode 100644 index 0000000000..4ae65f98c9 --- /dev/null +++ b/src/dependencies/MeshSDFilter/FilteringOptions.txt @@ -0,0 +1,20 @@ +#### Option file for SD filtering of mesh face normals +#### Lines starting with '#' are comments + +## Regularization weight, must be positive. +Lambda 10 + +## Gaussian standard deviation for spatial weight, scaled by the average distance between adjacent face cetroids. Must be positive. +Eta 2.5 + +## Gaussian standard deviation for guidance weight, must be positive. +Mu 20 + +## Gaussian standard deviation for signal weight, must be positive +Nu 0.26 + +## Closeness weight for mesh update, must be positive +MeshUpdateClosenessWeight 0.001 + +## Iterations for mesh update, must be positive +MeshUpdateIterations 20 \ No newline at end of file diff --git a/src/dependencies/MeshSDFilter/LICENSE b/src/dependencies/MeshSDFilter/LICENSE new file mode 100644 index 0000000000..bf12463607 --- /dev/null +++ b/src/dependencies/MeshSDFilter/LICENSE @@ -0,0 +1,29 @@ +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. diff --git a/src/dependencies/MeshSDFilter/MeshDenoiser.cpp b/src/dependencies/MeshSDFilter/MeshDenoiser.cpp new file mode 100644 index 0000000000..248f126800 --- /dev/null +++ b/src/dependencies/MeshSDFilter/MeshDenoiser.cpp @@ -0,0 +1,96 @@ +// 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. + + +#include "MeshTypes.h" +#include "MeshNormalDenoising.h" +#include +#include +#include + + + +int main(int argc, char **argv) +{ + if(argc != 4) + { + std::cout << "Usage: MeshDenoiser OPTION_FILE INPUT_MESH OUTPUT_MESH" << std::endl; + return 1; + } + + TriMesh mesh; + if(!OpenMesh::IO::read_mesh(mesh, argv[2])) + { + std::cerr << "Error: unable to read input mesh from the file " << argv[2] << std::endl; + return 1; + } + + #ifdef USE_OPENMP + Eigen::initParallel(); + #endif + + // Load option file + SDFilter::MeshDenoisingParameters param; + if(!param.load(argv[1])){ + std::cerr << "Error: unable to load option file " << argv[1] << std::endl; + return 1; + } + if(!param.valid_parameters()){ + std::cerr << "Invalid filter options. Aborting..." << std::endl; + return 1; + } + param.output(); + + + // 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 + SDFilter::MeshNormalDenoising denoiser(mesh); + TriMesh output_mesh; + denoiser.denoise(param, output_mesh); + + SDFilter::restore_mesh(output_mesh, original_center, original_scale); + + // Save output mesh + if(!SDFilter::write_mesh_high_accuracy(output_mesh, argv[3])){ + std::cerr << "Error: unable to save the result mesh to file " << argv[3] << std::endl; + return 1; + } + + return 0; +} + + + + + diff --git a/src/dependencies/MeshSDFilter/MeshNormalDenoising.h b/src/dependencies/MeshSDFilter/MeshNormalDenoising.h new file mode 100644 index 0000000000..7343c4d40a --- /dev/null +++ b/src/dependencies/MeshSDFilter/MeshNormalDenoising.h @@ -0,0 +1,264 @@ +// 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(); + + 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; + } + + // 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(); + std::cout << "Denoising completed, timing: " << + timer.elapsed_time(denoise_begin_time, denoise_end_time) << 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); + + 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()); + } + } + } + + OMP_FOR + 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(); + faces_in_patch.insert(faces_in_patch.end(), adj_faces_per_vtx[vtx_idx].begin(), adj_faces_per_vtx[vtx_idx].end()); + edges_in_patch.insert(edges_in_patch.end(), adj_nonboundary_edges_per_vtx[vtx_idx].begin(), adj_nonboundary_edges_per_vtx[vtx_idx].end()); + } + + // Sort and remove duplicates + std::sort(faces_in_patch.begin(), faces_in_patch.end()); + faces_in_patch.erase(std::unique(faces_in_patch.begin(), faces_in_patch.end()), faces_in_patch.end()); + neighborhood_faces_per_face[i] = faces_in_patch; + + std::sort(edges_in_patch.begin(), edges_in_patch.end()); + edges_in_patch.erase(std::unique(edges_in_patch.begin(), edges_in_patch.end()), edges_in_patch.end()); + + // 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_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/dependencies/MeshSDFilter/MeshNormalFilter.h b/src/dependencies/MeshSDFilter/MeshNormalFilter.h new file mode 100644 index 0000000000..fd99766e78 --- /dev/null +++ b/src/dependencies/MeshSDFilter/MeshNormalFilter.h @@ -0,0 +1,786 @@ +// 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 + +namespace SDFilter +{ + +class MeshFilterParameters : public Parameters +{ +public: + MeshFilterParameters() + :mesh_update_method(ITERATIVE_UPDATE), mesh_update_closeness_weight(0.001), mesh_update_iter(20) + { + // 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 [Want et al. 2015] + }; + + 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 + + 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; + } + + return true; + } + +protected: + virtual bool load_option(const OptionInterpreter &opt) + { + return Parameters::load_option(opt) || + opt.load("MeshUpdateClosenessWeight", mesh_update_closeness_weight) || + opt.load("MeshUpdateIterations", mesh_update_iter); + } + + virtual void output_options() + { + Parameters::output_options(); + std::cout << "Mesh update closeness weight: " << mesh_update_closeness_weight << std::endl; + std::cout << "Mesh update iterations: " << mesh_update_iter << 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){} + + virtual ~MeshNormalFilter() {} + + // Pass param by value, to allow changing the value of eta + bool filter(MeshFilterParameters param, TriMesh &output_mesh) + { + assert(param.valid_parameters()); + + 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(); + + if(print_timing_){ + std::cout << "Mesh udpate timing: " << timer.elapsed_time(update_begin_time, update_end_time) << " secs" << std::endl; + std::cout << "Mesh filter total timing: " << timer.elapsed_time(mesh_flter_begin_time, update_end_time) << " 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) + { + // Neighborhood radius set to three times of eta + double radius = 3.0 * param.eta; + int n_faces = mesh_.n_faces(); + int n_vtx = mesh_.n_vertices(); + + Matrix3X face_centroids(3, n_faces); + std::vector< std::vector > upper_neighbor_lists(n_faces), lower_neighbor_lists(n_faces); // Lists of neighbor faces with indices larger/smaller than the current face + + OMP_PARALLEL + { + // Pre-compute face centroids + OMP_FOR + for(int i = 0; i < n_faces; ++ i) + { + face_centroids.col(i) = to_eigen_vec3d(mesh_.calc_face_centroid(TriMesh::FaceHandle(i))); + } + + // Store neighbor lists for each face + OMP_FOR + for(int i = 0; i < n_faces; ++ i) + { + std::queue face_queue; + std::vector face_visited(n_faces, false), vtx_visited(n_vtx, false); + face_visited[i] = true; + face_queue.push(i); + + Eigen::Vector3d center_face_centroid = face_centroids.col(i); + + while(!face_queue.empty()) + { + TriMesh::FaceHandle fh(face_queue.front()); + face_queue.pop(); + + // Enumerate the three vertices of the current face, and check their neighboring faces + for(TriMesh::ConstFaceVertexIter cfv_it = mesh_.cfv_iter(fh); cfv_it.is_valid(); ++ cfv_it) + { + int current_vtx = cfv_it->idx(); + + if(!vtx_visited[current_vtx]) + { + vtx_visited[current_vtx] = true; + + for(TriMesh::ConstVertexFaceIter cvf_it = mesh_.cvf_iter(*cfv_it); cvf_it.is_valid(); ++ cvf_it) + { + int current_face = cvf_it->idx(); + + if(!face_visited[current_face]) + { + face_visited[current_face] = true; + + double dist = (face_centroids.col(current_face) - center_face_centroid).norm(); + + if(dist < radius) + { + face_queue.push(current_face); + + if(current_face < i){ + // Make sure each neigbor list only stores indices smaller than its face index + lower_neighbor_lists[i].push_back(current_face); + } + else if(current_face > i){ + upper_neighbor_lists[i].push_back(current_face); + } + } + } + } + } + } + } + } + + // Sort the neighbor lists + OMP_FOR + for(int i = 0; i < n_faces; ++ i) + { + std::sort(upper_neighbor_lists[i].begin(), upper_neighbor_lists[i].end()); + std::sort(lower_neighbor_lists[i].begin(), lower_neighbor_lists[i].end()); + } + } + + // For each upper neighbor list, check the lower neighbor list of the candidate neighbor to make sure the neighboring relationship is symmetric + + Eigen::VectorXi lower_neighbor_front_idx(n_faces); // Indices for keeping track of elements that have been checked in each lower neighbor list + lower_neighbor_front_idx.setZero(); + VectorXIdx segment_start_addr(n_faces); // Starting address for the segment of each face within the neighbor pair array + Eigen::Index n_neighbor_pairs = 0; + + for(int i = 0; i < n_faces; ++ i) + { + std::vector ¤t_upper_list = upper_neighbor_lists[i]; + int j = 0; + + while(j < static_cast(current_upper_list.size())) + { + // Try to locate the current face in the lower list of the current upper neighbor. + // If it can be found, then the neighboring relationship is symmetric + int neighbor_idx = current_upper_list[j]; + std::vector &lower_list = lower_neighbor_lists[neighbor_idx]; + + int lower_list_size = lower_list.size(); + int current_lower_list_idx = lower_neighbor_front_idx(neighbor_idx); + + while(current_lower_list_idx < lower_list_size && lower_list[current_lower_list_idx] < i){ + current_lower_list_idx ++; + } + + lower_neighbor_front_idx(neighbor_idx) = current_lower_list_idx; + + // If the current upper neighbor is not symmetric, remove it from the upper neighbor list. + // Otherwise, move on to check the next upper neighbor. + if(current_lower_list_idx < lower_list_size && lower_list[current_lower_list_idx] == i){ + j++; + } + else{ + current_upper_list.erase(current_upper_list.begin() + j); + } + } + + segment_start_addr(i) = n_neighbor_pairs; + n_neighbor_pairs += current_upper_list.size(); + } + + neighbor_pairs.resize(2, n_neighbor_pairs); + neighbor_dist.resize(n_neighbor_pairs); + + OMP_PARALLEL + { + OMP_FOR + for(int i = 0; i < n_faces; ++ i) + { + std::vector &upper_neighbors = upper_neighbor_lists[i]; + + if(!upper_neighbors.empty()){ + Eigen::Index start_col = segment_start_addr(i); + Eigen::Index n_cols = upper_neighbors.size(); + neighbor_pairs.block(0, start_col, 1, n_cols).setConstant(i); + neighbor_pairs.block(1, start_col, 1, n_cols) = Eigen::Map(upper_neighbors.data(), upper_neighbors.size()).transpose(); + } + } + + OMP_FOR + for(Eigen::Index i = 0; i < n_neighbor_pairs; ++ i) + { + int idx1 = neighbor_pairs(0, i), idx2 = neighbor_pairs(1, i); + neighbor_dist(i) = (face_centroids.col(idx1) - face_centroids.col(idx2)).norm(); + } + } + + //std::cout << "n_neighbor_pairs:" << n_neighbor_pairs << std::endl; + //std::cout << "Average neighborhood size: " << 2.0 * n_neighbor_pairs / n_faces << std::endl; + + return n_neighbor_pairs > Eigen::Index(0); + } + + + 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; + } + + + void set_mesh(const TriMesh &mesh, bool invalidate_update_system) + { + mesh_ = mesh; + + 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 + + + + + // 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) + { + if(system_matrix_factorized_) + { + return true; + } + + int n_faces = mesh_.n_faces(); + int n_vtx = mesh_.n_vertices(); + 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; + return true; + } + + + 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_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); + + Eigen::MatrixX3d wX0 = vtx_pos.transpose() * w_closeness; // Part of the linear system right-hand-side that corresponds to initial vertex positions + Eigen::MatrixX3d B(3 * n_faces, 3); // Per-face target position of the new vertices + + int n_vtx = output_mesh.n_vertices(); + Eigen::MatrixX3d 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::Matrix3Xd 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 = At_ * B + wX0; + if(!linear_solver_.solve(rhs, sol)){ + std::cerr << "Error: failed to solve mesh update system" << std::endl; + return false; + } + + vtx_pos = sol.transpose(); + set_vertex_points(output_mesh, vtx_pos); + } + + 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/dependencies/MeshSDFilter/MeshSDFilter.cpp b/src/dependencies/MeshSDFilter/MeshSDFilter.cpp new file mode 100644 index 0000000000..ccfcfec405 --- /dev/null +++ b/src/dependencies/MeshSDFilter/MeshSDFilter.cpp @@ -0,0 +1,93 @@ +// 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. + + +#include "MeshTypes.h" +#include "MeshNormalFilter.h" +#include +#include +#include + + + +int main(int argc, char **argv) +{ + if(argc != 4) + { + std::cout << "Usage: MeshSDFilter OPTION_FILE INPUT_MESH OUTPUT_MESH" << std::endl; + return 1; + } + + + + TriMesh mesh; + if(!OpenMesh::IO::read_mesh(mesh, argv[2])) + { + std::cerr << "Error: unable to read input mesh from the file " << argv[2] << std::endl; + return 1; + } + + #ifdef USE_OPENMP + Eigen::initParallel(); + #endif + + // Load option file + SDFilter::MeshFilterParameters param; + if(!param.load(argv[1])){ + std::cerr << "Error: unable to load option file " << argv[1] << std::endl; + return 1; + } + if(!param.valid_parameters()){ + std::cerr << "Invalid filter options. Aborting..." << std::endl; + return 1; + } + param.output(); + + + // 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 + TriMesh output_mesh; + SDFilter::MeshNormalFilter mesh_filter(mesh); + mesh_filter.filter(param, output_mesh); + + SDFilter::restore_mesh(output_mesh, original_center, original_scale); + + // Save output mesh + if(!SDFilter::write_mesh_high_accuracy(output_mesh, argv[3])){ + std::cerr << "Error: unable to save the result mesh to file " << argv[3] << std::endl; + return 1; + } + + return 0; +} diff --git a/src/dependencies/MeshSDFilter/MeshTypes.h b/src/dependencies/MeshSDFilter/MeshTypes.h new file mode 100644 index 0000000000..2928ad1170 --- /dev/null +++ b/src/dependencies/MeshSDFilter/MeshTypes.h @@ -0,0 +1,248 @@ +// 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 "EigenTypes.h" + +#include + + +/// 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]; + } + } +} + +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])); + } +} + +// Write a mesh to an ASCII file with high accuracy +template +inline bool write_mesh_high_accuracy(const MeshT &mesh, const std::string &filename) +{ + return OpenMesh::IO::write_mesh(mesh, filename, OpenMesh::IO::Options::Default, 16); +} + +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); +} + +} + +#endif // MESHTYPES_H diff --git a/src/dependencies/MeshSDFilter/README b/src/dependencies/MeshSDFilter/README new file mode 100644 index 0000000000..a4aa637962 --- /dev/null +++ b/src/dependencies/MeshSDFilter/README @@ -0,0 +1,64 @@ +MeshSDFilter -- Static/Dynamic Filtering for Mesh Geometry +=========================================================== + + +This code implements the mesh normal filtering algorithm from the following paper: + +- Juyong Zhang, Bailin Deng, Yang Hong, Yue Peng, Wenjie Qin, Ligang Liu. Static/Dynamic Filtering for Mesh Geometry. arXiv:1712.03574. + + + +1. Compiling + +The code requires c++ compiplers with c++11 support. It has been tested on the following systems and compilers: +- Windows 10 using Visual Studio 2015; +- macOS sierra, using homebrew GCC 7, or clang compiler from Xcode 8.3; +- Debian 9.1 using GCC 6.3. + + +Follow the following steps to compile the code: + +1) Make sure Eigen is installed. We recommend version 3.3+. + - Download Eigen from eigen.tuxfamily.org and extract it into a folder 'eigen' within the 'external' folder. Make sure the file 'external/eigen/Eigen/Dense' can be found + - Alternatively: 1) On Ubuntu and Debian, use the command "apt-get install libeigen3-dev" to install Eigen; 2) On macOS, install homebrew (https://brew.sh/) and run "brew install cmake open-mesh eigen". + +2) Create a folder “build” within the root directory of the code + +3) Run cmake to generate the build files inside the build folder, and compile the source code: + +- On linux or mac, run the following commands within the build folder: + $ cmake -DCMAKE_BUILD_TYPE=Release .. + $ make + +- On mac, the default Apple clang compiler does not support OpenMP. To enable OpenMP, first install the homebrew GCC: + $ brew install gcc + This should install gcc 7, with compiler commands gcc-7 and g++7. Then run the following command within the build folder: + $ cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_C_COMPILER=gcc-7 -DCMAKE_CXX_COMPILER=g++-7 .. + $ make + +- On windows, use the cmake GUI to generate a visual studio solution file, and build the solution. + +4) Afterwards, there should be two exectuable files 'MeshSDFilter' and 'MeshDenoiser' generated. + + + +2. Usage + +Use one of the following commands to filter or denoise a mesh: + $ MeshSDFilter OPTION_FILE INPUT_MESH OUTPUT_MESH + $ MeshDenoiser OPTION_FILE INPUT_MESH OUTPUT_MESH + +- INPUT_MESH and OUTPUT_MESH are mesh files for input and result. +- OPTION_FILE is a file that specifies the parameters. Example files for filtering and denoising are provided in FilteringOptions.txt and DenoisingOptions.txt . + + + +3. License + +The code is released under BSD 3-Clause License. + + + +4. Contact + +Feel free to contact Bailin Deng if you have any comments or questions. \ No newline at end of file diff --git a/src/dependencies/MeshSDFilter/SDFilter.h b/src/dependencies/MeshSDFilter/SDFilter.h new file mode 100644 index 0000000000..387e1df3c8 --- /dev/null +++ b/src/dependencies/MeshSDFilter/SDFilter.h @@ -0,0 +1,714 @@ +// 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" +#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 +{ + +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) {} + + virtual ~Parameters() {} + + enum LinearSolverType + { + CG, + LDLT + }; + + 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 normalize_iterates; // Normalization of the filtered normals in each iteration + + + // 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; + } + + 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); + } + + 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; + } + + +}; + + +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 +{ +protected: + + SDFilter() + :signal_dim_(-1), signal_count_(-1), print_progress_(true), + print_timing_(true), print_diagnostic_info_(false){} + + virtual ~SDFilter(){} + + + bool filter(Parameters param) + { + 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(); + + if(print_timing_){ + std::cout << "Preprocessing timing: " << timer.elapsed_time(begin_time, preprocess_end_time) << " secs" << std::endl; + std::cout << "Filtering timing: " << timer.elapsed_time(preprocess_end_time, filter_end_time) << " 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){ + std::cout << "Solver converged after " << num_iter << " iterations" << std::endl; + break; + } + else if(num_iter == param.max_iter){ + 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), pattern_analyzed(false){} + + // Initialize the solver with matrix + bool compute(const SparseMatrixXd &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::CG){ + 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) + { +#ifdef USE_OPENMP + int n_cols = 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 + + return check_error(LDLT_solver_, "LDLT solve failed"); + } + 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; + if(solver_type_ == Parameters::LDLT){ + reset_pattern(); + } + } + + private: + Parameters::LinearSolverType solver_type_; + Eigen::SimplicialLDLT LDLT_solver_; + Eigen::ConjugateGradient CG_solver_; + + bool pattern_analyzed; // Flag for symbolic factorization + + 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_; + + // 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_ + std::vector< std::vector< Eigen::Index > > neighbors(signal_count_); + for(Eigen::Index i = 0; i < n_neighbor_pairs; ++ i) + { + int idx1 = neighboring_pairs_(0, i), idx2 = neighboring_pairs_(1, i); + + neighbors[idx1].push_back(idx2); + neighbors[idx1].push_back(i); + + neighbors[idx2].push_back(idx1); + neighbors[idx2].push_back(i); + } + + 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) + neighbors[i].size() / 2; + } + + neighborhood_info_.resize(2, 2 * n_neighbor_pairs); + + OMP_PARALLEL + { + OMP_FOR + for(int i = 0; i < signal_count_; ++ i) + { + std::vector ¤t_neighbor_info = neighbors[i]; + + if(!current_neighbor_info.empty()) + { + Eigen::Index n_cols = current_neighbor_info.size() / 2; + neighborhood_info_.block(0, neighborhood_info_boundaries_(i), 2, n_cols) = + Eigen::Map< Matrix2XIdx >(current_neighbor_info.data(), 2, n_cols); + } + } + } + + 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_ */