From 822a3a3e1a8ebb781cc3149dfbb337e370eff5e3 Mon Sep 17 00:00:00 2001 From: Michele Morisco <43170306+michisco@users.noreply.github.com> Date: Thu, 25 Aug 2022 10:30:17 +0200 Subject: [PATCH] Cubic Stylization filter added --- .../filter_cubization/CMakeLists.txt | 28 ++ .../CubicStylizationFiles/README.md | 18 + .../CubicStylizationFiles/conversionMeshes.h | 49 +++ .../CubicStylizationFiles/cubic_stylizing.h | 53 +++ .../CubicStylizationFiles/src/LICENSE.MPL2 | 373 ++++++++++++++++++ .../src/cube_style_data.h | 78 ++++ .../src/cube_style_precomputation.cpp | 70 ++++ .../src/cube_style_precomputation.h | 24 ++ .../src/cube_style_single_iteration.cpp | 41 ++ .../src/cube_style_single_iteration.h | 23 ++ .../src/fit_rotations_l1.cpp | 100 +++++ .../src/fit_rotations_l1.h | 20 + .../src/normalize_unitbox.cpp | 8 + .../src/normalize_unitbox.h | 15 + .../src/orthogonal_procrustes.cpp | 34 ++ .../src/orthogonal_procrustes.h | 12 + .../CubicStylizationFiles/src/shrinkage.cpp | 32 ++ .../CubicStylizationFiles/src/shrinkage.h | 11 + .../filter_cubization/curvdata.h | 93 +++++ .../filter_cubization/curvedgeflip.h | 339 ++++++++++++++++ .../filter_cubization/filter_cubization.cpp | 310 +++++++++++++++ .../filter_cubization/filter_cubization.h | 73 ++++ 22 files changed, 1804 insertions(+) create mode 100644 src/meshlabplugins/filter_cubization/CMakeLists.txt create mode 100644 src/meshlabplugins/filter_cubization/CubicStylizationFiles/README.md create mode 100644 src/meshlabplugins/filter_cubization/CubicStylizationFiles/conversionMeshes.h create mode 100644 src/meshlabplugins/filter_cubization/CubicStylizationFiles/cubic_stylizing.h create mode 100644 src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/LICENSE.MPL2 create mode 100644 src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/cube_style_data.h create mode 100644 src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/cube_style_precomputation.cpp create mode 100644 src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/cube_style_precomputation.h create mode 100644 src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/cube_style_single_iteration.cpp create mode 100644 src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/cube_style_single_iteration.h create mode 100644 src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/fit_rotations_l1.cpp create mode 100644 src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/fit_rotations_l1.h create mode 100644 src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/normalize_unitbox.cpp create mode 100644 src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/normalize_unitbox.h create mode 100644 src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/orthogonal_procrustes.cpp create mode 100644 src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/orthogonal_procrustes.h create mode 100644 src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/shrinkage.cpp create mode 100644 src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/shrinkage.h create mode 100644 src/meshlabplugins/filter_cubization/curvdata.h create mode 100644 src/meshlabplugins/filter_cubization/curvedgeflip.h create mode 100644 src/meshlabplugins/filter_cubization/filter_cubization.cpp create mode 100644 src/meshlabplugins/filter_cubization/filter_cubization.h diff --git a/src/meshlabplugins/filter_cubization/CMakeLists.txt b/src/meshlabplugins/filter_cubization/CMakeLists.txt new file mode 100644 index 000000000..64473bedd --- /dev/null +++ b/src/meshlabplugins/filter_cubization/CMakeLists.txt @@ -0,0 +1,28 @@ +# Copyright 2019-2020, Collabora, Ltd. +# SPDX-License-Identifier: BSL-1.0 + + +set(SOURCES + filter_cubization.cpp + CubicStylizationFiles/src/cube_style_precomputation.cpp + CubicStylizationFiles/src/cube_style_single_iteration.cpp + CubicStylizationFiles/src/fit_rotations_l1.cpp + CubicStylizationFiles/src/normalize_unitbox.cpp + CubicStylizationFiles/src/orthogonal_procrustes.cpp + CubicStylizationFiles/src/shrinkage.cpp) + +set(HEADERS + filter_cubization.h + curvdata.h + curvedgeflip.h + CubicStylizationFiles/conversionMeshes.h + CubicStylizationFiles/cubic_stylizing.h + CubicStylizationFiles/src/cube_style_data.h + CubicStylizationFiles/src/cube_style_precomputation.h + CubicStylizationFiles/src/cube_style_single_iteration.h + CubicStylizationFiles/src/fit_rotations_l1.h + CubicStylizationFiles/src/normalize_unitbox.h + CubicStylizationFiles/src/orthogonal_procrustes.h + CubicStylizationFiles/src/shrinkage.h) + +add_meshlab_plugin(filter_cubization ${SOURCES} ${HEADERS}) diff --git a/src/meshlabplugins/filter_cubization/CubicStylizationFiles/README.md b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/README.md new file mode 100644 index 000000000..8a8694e26 --- /dev/null +++ b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/README.md @@ -0,0 +1,18 @@ +# Cubic Stylization filter + +Cubic stylization is a 3D stylization tool. Unlike image stylization (2D to 2D) and non-photorealistic rendering (3D to 2D), cubic stylization is a 3D to 3D stylization algorithm which takes a manifold triangle mesh as the input and outputs a cubified triangle mesh. + +This is a MeshLab filter based on "[Cubic Stylization](https://www.dgp.toronto.edu/projects/cubic-stylization/)" by [Hsueh-Ti Derek Liu](https://www.dgp.toronto.edu/~hsuehtil/) and [Alec Jacobson](https://www.cs.toronto.edu/~jacobson/). + +The folder `/src` contains the original C++ implementation of Cubic stylization performed by [Hsueh-Ti Derek Liu](https://www.dgp.toronto.edu/~hsuehtil/) and [Alec Jacobson](https://www.cs.toronto.edu/~jacobson/). +Their code is licensed under [MPL2](https://www.mozilla.org/en-US/MPL/2.0/). + +## Bibtex +``` +@article{Liu:CubicStyle:2019, + title = {Cubic Stylization}, + author = {Hsueh-Ti Derek Liu and Alec Jacobson}, + year = {2019}, + journal = {ACM Transactions on Graphics}, +} +``` diff --git a/src/meshlabplugins/filter_cubization/CubicStylizationFiles/conversionMeshes.h b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/conversionMeshes.h new file mode 100644 index 000000000..75fc0e3be --- /dev/null +++ b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/conversionMeshes.h @@ -0,0 +1,49 @@ +#ifndef CONVERSIONMESHES_H +#define CONVERSIONMESHES_H + +#include +#include + +namespace vcg{ +namespace tri{ + template + void Mesh2Matrix(MeshType& convertingMesh, Eigen::MatrixXd &output_vertexes, Eigen::MatrixXi &output_faces){ + + typedef typename MeshType::ScalarType ScalarType; + typedef typename Eigen::Matrix MatrixXm; + + MatrixXm temp_verts; + //convert mesh in matrix + vcg::tri::MeshToMatrix< MeshType >::GetTriMeshData( convertingMesh, temp_verts, output_faces); + output_vertexes = temp_verts.template cast(); + } + + template + void Matrix2Mesh(MeshType& convertedMesh, Eigen::MatrixXd vertexes, Eigen::MatrixXi faces){ + typedef typename MeshType::VertexPointer VertexPointer; + typedef typename MeshType::FaceIterator FaceIterator; + + //reconvert V and F matrixes in a mesh + convertedMesh.Clear(); + //vcg::tri::Allocator::AddVertices(convertedMesh,vertexes); + //vcg::tri::Allocator::AddFaces(convertedMesh,faces); + + Allocator::AddVertices(convertedMesh,vertexes.rows()); + Allocator::AddFaces(convertedMesh,faces.rows()); + VertexPointer ivp[vertexes.rows()]; + + int i; + for (i=0; i < vertexes.rows(); i++){ + for (int j = 0; j < 3; j++) + convertedMesh.vert[i].P()[j] = vertexes(i,j); + ivp[i]=&convertedMesh.vert[i]; + } + + FaceIterator fi; + for (i=0,fi=convertedMesh.face.begin();fi!=convertedMesh.face.end();i++,fi++) + for (int j = 0; j < 3; j++) + (*fi).V(j) = ivp[faces(i, j)]; + } +}} + +#endif // CONVERSIONMESHES_H diff --git a/src/meshlabplugins/filter_cubization/CubicStylizationFiles/cubic_stylizing.h b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/cubic_stylizing.h new file mode 100644 index 000000000..5dbb5b529 --- /dev/null +++ b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/cubic_stylizing.h @@ -0,0 +1,53 @@ +#ifndef CUBIC_STYLIZING_H +#define CUBIC_STYLIZING_H + +#include +#include +#include +#include + +#include +#include + +#include + +namespace vcg{ +namespace tri{ + +template +class Cubization{ + public: + static void Init(MeshType& m, Eigen::MatrixXd& V, Eigen::MatrixXd& U, Eigen::MatrixXi& F, cube_style_data& data){ + + // check requirements + vcg::tri::VertexVectorHasPerVertexTexCoord( m.vert ); + vcg::tri::VertexVectorHasPerVertexFlags( m.vert ); + + Mesh2Matrix(m, V, F); + + normalize_unitbox(V); + Eigen::RowVector3d meanV = V.colwise().mean(); + V = V.rowwise() - meanV; + U = V; + + data.bc.resize(1,3); + data.bc << V.row(F(0,0)); + + data.b.resize(1); + data.b << F(0,0); + + // precomputation ARAP and initialize ADMM parameters + cube_style_precomputation(V,F,data); + } + + static double Stylize(Eigen::MatrixXd &V, Eigen::MatrixXd &U, Eigen::MatrixXi &F, cube_style_data& data, Eigen::VectorXd& energy_verts, double& energyCubic) + { + // apply cubic stylization + cube_style_single_iteration(V,U,data, energy_verts); + energyCubic = data.objVal; + + return data.reldV; + } +}; +}} +#endif // CUBIC_STYLIZING_H diff --git a/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/LICENSE.MPL2 b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/LICENSE.MPL2 new file mode 100644 index 000000000..f4bbcd200 --- /dev/null +++ b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/LICENSE.MPL2 @@ -0,0 +1,373 @@ +Mozilla Public License Version 2.0 +================================== + +1. Definitions +-------------- + +1.1. "Contributor" + means each individual or legal entity that creates, contributes to + the creation of, or owns Covered Software. + +1.2. "Contributor Version" + means the combination of the Contributions of others (if any) used + by a Contributor and that particular Contributor's Contribution. + +1.3. "Contribution" + means Covered Software of a particular Contributor. + +1.4. "Covered Software" + means Source Code Form to which the initial Contributor has attached + the notice in Exhibit A, the Executable Form of such Source Code + Form, and Modifications of such Source Code Form, in each case + including portions thereof. + +1.5. "Incompatible With Secondary Licenses" + means + + (a) that the initial Contributor has attached the notice described + in Exhibit B to the Covered Software; or + + (b) that the Covered Software was made available under the terms of + version 1.1 or earlier of the License, but not also under the + terms of a Secondary License. + +1.6. "Executable Form" + means any form of the work other than Source Code Form. + +1.7. "Larger Work" + means a work that combines Covered Software with other material, in + a separate file or files, that is not Covered Software. + +1.8. "License" + means this document. + +1.9. "Licensable" + means having the right to grant, to the maximum extent possible, + whether at the time of the initial grant or subsequently, any and + all of the rights conveyed by this License. + +1.10. "Modifications" + means any of the following: + + (a) any file in Source Code Form that results from an addition to, + deletion from, or modification of the contents of Covered + Software; or + + (b) any new file in Source Code Form that contains any Covered + Software. + +1.11. "Patent Claims" of a Contributor + means any patent claim(s), including without limitation, method, + process, and apparatus claims, in any patent Licensable by such + Contributor that would be infringed, but for the grant of the + License, by the making, using, selling, offering for sale, having + made, import, or transfer of either its Contributions or its + Contributor Version. + +1.12. "Secondary License" + means either the GNU General Public License, Version 2.0, the GNU + Lesser General Public License, Version 2.1, the GNU Affero General + Public License, Version 3.0, or any later versions of those + licenses. + +1.13. "Source Code Form" + means the form of the work preferred for making modifications. + +1.14. "You" (or "Your") + means an individual or a legal entity exercising rights under this + License. For legal entities, "You" includes any entity that + controls, is controlled by, or is under common control with You. For + purposes of this definition, "control" means (a) the power, direct + or indirect, to cause the direction or management of such entity, + whether by contract or otherwise, or (b) ownership of more than + fifty percent (50%) of the outstanding shares or beneficial + ownership of such entity. + +2. License Grants and Conditions +-------------------------------- + +2.1. Grants + +Each Contributor hereby grants You a world-wide, royalty-free, +non-exclusive license: + +(a) under intellectual property rights (other than patent or trademark) + Licensable by such Contributor to use, reproduce, make available, + modify, display, perform, distribute, and otherwise exploit its + Contributions, either on an unmodified basis, with Modifications, or + as part of a Larger Work; and + +(b) under Patent Claims of such Contributor to make, use, sell, offer + for sale, have made, import, and otherwise transfer either its + Contributions or its Contributor Version. + +2.2. Effective Date + +The licenses granted in Section 2.1 with respect to any Contribution +become effective for each Contribution on the date the Contributor first +distributes such Contribution. + +2.3. Limitations on Grant Scope + +The licenses granted in this Section 2 are the only rights granted under +this License. No additional rights or licenses will be implied from the +distribution or licensing of Covered Software under this License. +Notwithstanding Section 2.1(b) above, no patent license is granted by a +Contributor: + +(a) for any code that a Contributor has removed from Covered Software; + or + +(b) for infringements caused by: (i) Your and any other third party's + modifications of Covered Software, or (ii) the combination of its + Contributions with other software (except as part of its Contributor + Version); or + +(c) under Patent Claims infringed by Covered Software in the absence of + its Contributions. + +This License does not grant any rights in the trademarks, service marks, +or logos of any Contributor (except as may be necessary to comply with +the notice requirements in Section 3.4). + +2.4. Subsequent Licenses + +No Contributor makes additional grants as a result of Your choice to +distribute the Covered Software under a subsequent version of this +License (see Section 10.2) or under the terms of a Secondary License (if +permitted under the terms of Section 3.3). + +2.5. Representation + +Each Contributor represents that the Contributor believes its +Contributions are its original creation(s) or it has sufficient rights +to grant the rights to its Contributions conveyed by this License. + +2.6. Fair Use + +This License is not intended to limit any rights You have under +applicable copyright doctrines of fair use, fair dealing, or other +equivalents. + +2.7. Conditions + +Sections 3.1, 3.2, 3.3, and 3.4 are conditions of the licenses granted +in Section 2.1. + +3. Responsibilities +------------------- + +3.1. Distribution of Source Form + +All distribution of Covered Software in Source Code Form, including any +Modifications that You create or to which You contribute, must be under +the terms of this License. You must inform recipients that the Source +Code Form of the Covered Software is governed by the terms of this +License, and how they can obtain a copy of this License. You may not +attempt to alter or restrict the recipients' rights in the Source Code +Form. + +3.2. Distribution of Executable Form + +If You distribute Covered Software in Executable Form then: + +(a) such Covered Software must also be made available in Source Code + Form, as described in Section 3.1, and You must inform recipients of + the Executable Form how they can obtain a copy of such Source Code + Form by reasonable means in a timely manner, at a charge no more + than the cost of distribution to the recipient; and + +(b) You may distribute such Executable Form under the terms of this + License, or sublicense it under different terms, provided that the + license for the Executable Form does not attempt to limit or alter + the recipients' rights in the Source Code Form under this License. + +3.3. Distribution of a Larger Work + +You may create and distribute a Larger Work under terms of Your choice, +provided that You also comply with the requirements of this License for +the Covered Software. If the Larger Work is a combination of Covered +Software with a work governed by one or more Secondary Licenses, and the +Covered Software is not Incompatible With Secondary Licenses, this +License permits You to additionally distribute such Covered Software +under the terms of such Secondary License(s), so that the recipient of +the Larger Work may, at their option, further distribute the Covered +Software under the terms of either this License or such Secondary +License(s). + +3.4. Notices + +You may not remove or alter the substance of any license notices +(including copyright notices, patent notices, disclaimers of warranty, +or limitations of liability) contained within the Source Code Form of +the Covered Software, except that You may alter any license notices to +the extent required to remedy known factual inaccuracies. + +3.5. Application of Additional Terms + +You may choose to offer, and to charge a fee for, warranty, support, +indemnity or liability obligations to one or more recipients of Covered +Software. However, You may do so only on Your own behalf, and not on +behalf of any Contributor. You must make it absolutely clear that any +such warranty, support, indemnity, or liability obligation is offered by +You alone, and You hereby agree to indemnify every Contributor for any +liability incurred by such Contributor as a result of warranty, support, +indemnity or liability terms You offer. You may include additional +disclaimers of warranty and limitations of liability specific to any +jurisdiction. + +4. Inability to Comply Due to Statute or Regulation +--------------------------------------------------- + +If it is impossible for You to comply with any of the terms of this +License with respect to some or all of the Covered Software due to +statute, judicial order, or regulation then You must: (a) comply with +the terms of this License to the maximum extent possible; and (b) +describe the limitations and the code they affect. Such description must +be placed in a text file included with all distributions of the Covered +Software under this License. Except to the extent prohibited by statute +or regulation, such description must be sufficiently detailed for a +recipient of ordinary skill to be able to understand it. + +5. Termination +-------------- + +5.1. The rights granted under this License will terminate automatically +if You fail to comply with any of its terms. However, if You become +compliant, then the rights granted under this License from a particular +Contributor are reinstated (a) provisionally, unless and until such +Contributor explicitly and finally terminates Your grants, and (b) on an +ongoing basis, if such Contributor fails to notify You of the +non-compliance by some reasonable means prior to 60 days after You have +come back into compliance. Moreover, Your grants from a particular +Contributor are reinstated on an ongoing basis if such Contributor +notifies You of the non-compliance by some reasonable means, this is the +first time You have received notice of non-compliance with this License +from such Contributor, and You become compliant prior to 30 days after +Your receipt of the notice. + +5.2. If You initiate litigation against any entity by asserting a patent +infringement claim (excluding declaratory judgment actions, +counter-claims, and cross-claims) alleging that a Contributor Version +directly or indirectly infringes any patent, then the rights granted to +You by any and all Contributors for the Covered Software under Section +2.1 of this License shall terminate. + +5.3. In the event of termination under Sections 5.1 or 5.2 above, all +end user license agreements (excluding distributors and resellers) which +have been validly granted by You or Your distributors under this License +prior to termination shall survive termination. + +************************************************************************ +* * +* 6. Disclaimer of Warranty * +* ------------------------- * +* * +* Covered Software is provided under this License on an "as is" * +* basis, without warranty of any kind, either expressed, implied, or * +* statutory, including, without limitation, warranties that the * +* Covered Software is free of defects, merchantable, fit for a * +* particular purpose or non-infringing. The entire risk as to the * +* quality and performance of the Covered Software is with You. * +* Should any Covered Software prove defective in any respect, You * +* (not any Contributor) assume the cost of any necessary servicing, * +* repair, or correction. This disclaimer of warranty constitutes an * +* essential part of this License. No use of any Covered Software is * +* authorized under this License except under this disclaimer. * +* * +************************************************************************ + +************************************************************************ +* * +* 7. Limitation of Liability * +* -------------------------- * +* * +* Under no circumstances and under no legal theory, whether tort * +* (including negligence), contract, or otherwise, shall any * +* Contributor, or anyone who distributes Covered Software as * +* permitted above, be liable to You for any direct, indirect, * +* special, incidental, or consequential damages of any character * +* including, without limitation, damages for lost profits, loss of * +* goodwill, work stoppage, computer failure or malfunction, or any * +* and all other commercial damages or losses, even if such party * +* shall have been informed of the possibility of such damages. This * +* limitation of liability shall not apply to liability for death or * +* personal injury resulting from such party's negligence to the * +* extent applicable law prohibits such limitation. Some * +* jurisdictions do not allow the exclusion or limitation of * +* incidental or consequential damages, so this exclusion and * +* limitation may not apply to You. * +* * +************************************************************************ + +8. Litigation +------------- + +Any litigation relating to this License may be brought only in the +courts of a jurisdiction where the defendant maintains its principal +place of business and such litigation shall be governed by laws of that +jurisdiction, without reference to its conflict-of-law provisions. +Nothing in this Section shall prevent a party's ability to bring +cross-claims or counter-claims. + +9. Miscellaneous +---------------- + +This License represents the complete agreement concerning the subject +matter hereof. If any provision of this License is held to be +unenforceable, such provision shall be reformed only to the extent +necessary to make it enforceable. Any law or regulation which provides +that the language of a contract shall be construed against the drafter +shall not be used to construe this License against a Contributor. + +10. Versions of the License +--------------------------- + +10.1. New Versions + +Mozilla Foundation is the license steward. Except as provided in Section +10.3, no one other than the license steward has the right to modify or +publish new versions of this License. Each version will be given a +distinguishing version number. + +10.2. Effect of New Versions + +You may distribute the Covered Software under the terms of the version +of the License under which You originally received the Covered Software, +or under the terms of any subsequent version published by the license +steward. + +10.3. Modified Versions + +If you create software not governed by this License, and you want to +create a new license for such software, you may create and use a +modified version of this License if you rename the license and remove +any references to the name of the license steward (except to note that +such modified license differs from this License). + +10.4. Distributing Source Code Form that is Incompatible With Secondary +Licenses + +If You choose to distribute Source Code Form that is Incompatible With +Secondary Licenses under the terms of this version of the License, the +notice described in Exhibit B of this License must be attached. + +Exhibit A - Source Code Form License Notice +------------------------------------------- + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + +If it is not possible or desirable to put the notice in a particular +file, then You may include the notice in a location (such as a LICENSE +file in a relevant directory) where a recipient would be likely to look +for such a notice. + +You may add additional accurate notices of copyright ownership. + +Exhibit B - "Incompatible With Secondary Licenses" Notice +--------------------------------------------------------- + + This Source Code Form is "Incompatible With Secondary Licenses", as + defined by the Mozilla Public License, v. 2.0. \ No newline at end of file diff --git a/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/cube_style_data.h b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/cube_style_data.h new file mode 100644 index 000000000..f60eacccf --- /dev/null +++ b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/cube_style_data.h @@ -0,0 +1,78 @@ +#ifndef CUBE_STYLE_DATA_H +#define CUBE_STYLE_DATA_H + +#include +#include +#include +#include + +struct cube_style_data +{ + // user should tune these + double lambda = 0.0; + double rhoInit = 1e-3; + double ABSTOL = 1e-6; + double RELTOL = 1e-3; + + // usually these don't need to tune + double mu = 10; + double tao = 2; + double maxIter_ADMM = 100; + double objVal = 0; + double reldV = std::numeric_limits::max(); + + std::vector objHis; + std::vector hEList; + std::vector dVList, UHis; + std::vector WVecList; + + Eigen::SparseMatrix K, L; + Eigen::MatrixXd N, VA, zAll, uAll; + Eigen::VectorXd rhoAll, objValVec; + + // bool useBc = false; + Eigen::MatrixXd bc; + Eigen::VectorXi b; + + igl::min_quad_with_fixed_data solver_data; + + // for plane constraints + Eigen::VectorXi bx, by, bz; + igl::min_quad_with_fixed_data solver_data_x, solver_data_y, solver_data_z; + double xPlane = 0.0; + double yPlane = 0.0; + double zPlane = 0.0; + // for plane constraints + + void reset() + { + // user should tune these + ABSTOL = 1e-5; + rhoInit = 1e-3; + RELTOL = 1e-3; + + // usually these don't need to tune + double mu = 10; + double tao = 2; + double maxIter_ADMM = 100; + double objVal = 0; + double reldV = std::numeric_limits::max(); + + objHis.clear(); + hEList.clear(); + dVList.clear(); + UHis.clear(); + WVecList.clear(); + K = Eigen::SparseMatrix(); + L = Eigen::SparseMatrix(); + N = Eigen::MatrixXd(); + VA = Eigen::MatrixXd(); + zAll =Eigen::MatrixXd(); + uAll = Eigen::MatrixXd(); + rhoAll = Eigen::VectorXd(); + objValVec = Eigen::VectorXd(); + + igl::min_quad_with_fixed_data solver_data; + } +}; +#endif // CUBE_STYLE_DATA_H diff --git a/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/cube_style_precomputation.cpp b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/cube_style_precomputation.cpp new file mode 100644 index 000000000..1e648ad84 --- /dev/null +++ b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/cube_style_precomputation.cpp @@ -0,0 +1,70 @@ +#include "cube_style_precomputation.h" + +void cube_style_precomputation( + const Eigen::MatrixXd & V, + const Eigen::MatrixXi & F, + cube_style_data & data) +{ + using namespace Eigen; + using namespace std; + + data.reset(); + + igl::per_vertex_normals(V,F, data.N); + + igl::cotmatrix(V,F,data.L); + + SparseMatrix M; + igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_BARYCENTRIC,M); + data.VA = M.diagonal(); + + vector> adjFList, VI; + igl::vertex_triangle_adjacency(V.rows(),F,adjFList,VI); + + igl::arap_rhs(V,F,V.cols(),igl::ARAP_ENERGY_TYPE_SPOKES_AND_RIMS, data.K); + + data.hEList.resize(V.rows()); + data.WVecList.resize(V.rows()); + data.dVList.resize(V.rows()); + vector adjF; + for (int ii=0; ii(),false,data.solver_data); + + data.zAll.resize(3, V.rows()); data.zAll.setRandom(); + data.uAll.resize(3, V.rows()); data.uAll.setRandom(); + data.rhoAll.resize(V.rows()); data.rhoAll.setConstant(data.rhoInit); +} diff --git a/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/cube_style_precomputation.h b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/cube_style_precomputation.h new file mode 100644 index 000000000..f299f8d8d --- /dev/null +++ b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/cube_style_precomputation.h @@ -0,0 +1,24 @@ +#ifndef CUBE_STYLE_PRECOMPUTATION_H +#define CUBE_STYLE_PRECOMPUTATION_H + +#include +#include +#include +#include + +// include libigl functions +#include +#include +#include +#include +#include +#include +#include +#include + +#include +void cube_style_precomputation( + const Eigen::MatrixXd & V, + const Eigen::MatrixXi & F, + cube_style_data & data); +#endif diff --git a/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/cube_style_single_iteration.cpp b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/cube_style_single_iteration.cpp new file mode 100644 index 000000000..06bc6428a --- /dev/null +++ b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/cube_style_single_iteration.cpp @@ -0,0 +1,41 @@ +#include "cube_style_single_iteration.h" + +void cube_style_single_iteration( + const Eigen::MatrixXd & V, + Eigen::MatrixXd & U, + cube_style_data & data, + Eigen::VectorXd & energyVects) +{ + using namespace Eigen; + using namespace std; + + // local step + MatrixXd RAll(3,V.rows()*3); + { + fit_rotations_l1(V, U, RAll, data, energyVects); + } + + //R.setZero(); + //R = RAll; + + // global step + MatrixXd Upre = U; + { + VectorXd Rcol; + igl::columnize(RAll, V.rows(), 2, Rcol); + VectorXd Bcol = data.K * Rcol; + for(int dim=0; dim +#include +#include +#include + +// include libigl functions +#include +#include +#include + +// include cube flow functions +#include +#include + +void cube_style_single_iteration( + const Eigen::MatrixXd & V, + Eigen::MatrixXd & U, + cube_style_data & data, + Eigen::VectorXd & energyVects); +#endif // CUBE_STYLE_SINGLE_ITERATION_H diff --git a/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/fit_rotations_l1.cpp b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/fit_rotations_l1.cpp new file mode 100644 index 000000000..af71ad15f --- /dev/null +++ b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/fit_rotations_l1.cpp @@ -0,0 +1,100 @@ +#include "fit_rotations_l1.h" + +void fit_rotations_l1( + const Eigen::MatrixXd & V, + Eigen::MatrixXd & U, + Eigen::MatrixXd & RAll, + cube_style_data & data, + Eigen::VectorXd & energy_vects) +{ + using namespace Eigen; + using namespace std; + + data.objValVec.setZero(V.rows()); + igl::parallel_for( + V.rows(), + [&data, &RAll, &U](const int ii) + { + // warm start parameters + VectorXd z = data.zAll.col(ii); + VectorXd u = data.uAll.col(ii); + VectorXd n = data.N.row(ii).transpose(); + double rho = data.rhoAll(ii); + Matrix3d R; + + // get energy parameters + // Note: dVn = [dV n], dUn = [dU z-u] + MatrixXi hE = data.hEList[ii]; + MatrixXd dU(3,hE.rows()); + { + MatrixXd U_hE0, U_hE1; + igl::slice(U,hE.col(0),1,U_hE0); + igl::slice(U,hE.col(1),1,U_hE1); + dU = (U_hE1 - U_hE0).transpose(); + } + + // Note: + // S = [dV n] * [W 0; 0 rho] * [dU (z-u)]' + // = dV * W * dU' + n * rho * (z-u)' + // = Spre + n * rho * (z-u)' + MatrixXd dV = data.dVList[ii]; + VectorXd WVec = data.WVecList[ii]; + Matrix3d Spre = dV * WVec.asDiagonal() * dU.transpose(); + + // ADMM + for (int k=0; k data.mu * s_norm) + { + rho = data.tao * rho; + u = u / data.tao; + } + else if (s_norm > data.mu * r_norm) + { + rho = rho / data.tao; + u = u * data.tao; + } + + // stopping criteria + double nz = double(z.size()); + double eps_pri = sqrt(2.0*nz)*data.ABSTOL + data.RELTOL*max( (R*n).norm(),z.norm() ); + double eps_dual = sqrt(1.0*nz)*data.ABSTOL + data.RELTOL* ((rho*u).norm()); + if ( (r_norm +#include +#include +#include +#include +#include +#include +#include +#include + +void fit_rotations_l1( + const Eigen::MatrixXd & V, + Eigen::MatrixXd & U, + Eigen::MatrixXd & RAll, + cube_style_data & data, + Eigen::VectorXd & energy_vects); +#endif // FIT_ROTATIONS_L1_H diff --git a/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/normalize_unitbox.cpp b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/normalize_unitbox.cpp new file mode 100644 index 000000000..91c4d03c1 --- /dev/null +++ b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/normalize_unitbox.cpp @@ -0,0 +1,8 @@ +#include "normalize_unitbox.h" + +void normalize_unitbox( + Eigen::MatrixXd & V) +{ + V = V.rowwise() - V.colwise().minCoeff(); + V /= V.maxCoeff(); +} diff --git a/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/normalize_unitbox.h b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/normalize_unitbox.h new file mode 100644 index 000000000..97a619353 --- /dev/null +++ b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/normalize_unitbox.h @@ -0,0 +1,15 @@ +#ifndef NORMALIZE_UNITBOX_H +#define NORMALIZE_UNITBOX_H + +#include +#include + +// Inputs: +// V a matrix of vertex positions +// Outputs: +// V a matrix of vertex positions (in a unit box) +void normalize_unitbox( + Eigen::MatrixXd & V); + +#endif // NORMALIZE_UNITBOX_H + diff --git a/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/orthogonal_procrustes.cpp b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/orthogonal_procrustes.cpp new file mode 100644 index 000000000..bbe3fc846 --- /dev/null +++ b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/orthogonal_procrustes.cpp @@ -0,0 +1,34 @@ +#include "orthogonal_procrustes.h" + +void orthogonal_procrustes( + const Eigen::Matrix3d & S, + Eigen::Matrix3d & R) +{ + using namespace Eigen; + using namespace std; + + // using svd3x3 which has problems for double precision + // Matrix3d SU, SV; + // Matrix SS; + // igl::svd3x3(S, SU, SS, SV); + // R = SV * SU.transpose(); + // if (R.determinant() < 0) + // { + // SU.col(2) = -SU.col(2); + // R = SV * SU.transpose(); + // } + + // using Eigen svd + JacobiSVD svd; + svd.compute(S, Eigen::ComputeFullU | Eigen::ComputeFullV ); + Matrix3d SU = svd.matrixU(); + Matrix3d SV = svd.matrixV(); + R = SV * SU.transpose(); + if (R.determinant() < 0) + { + SU.col(2) = -SU.col(2); + R = SV * SU.transpose(); + } + + assert(R.determinant() > 0); +} diff --git a/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/orthogonal_procrustes.h b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/orthogonal_procrustes.h new file mode 100644 index 000000000..e81315abb --- /dev/null +++ b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/orthogonal_procrustes.h @@ -0,0 +1,12 @@ +#ifndef ORTHOGONAL_PROCRUSTES_H +#define ORTHOGONAL_PROCRUSTES_H + +#include +#include +#include +#include + +void orthogonal_procrustes( + const Eigen::Matrix3d & S, + Eigen::Matrix3d & R); +#endif // ORTHOGONAL_PROCRUSTES_H diff --git a/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/shrinkage.cpp b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/shrinkage.cpp new file mode 100644 index 000000000..5ab44bbeb --- /dev/null +++ b/src/meshlabplugins/filter_cubization/CubicStylizationFiles/src/shrinkage.cpp @@ -0,0 +1,32 @@ +#include "shrinkage.h" + +void shrinkage( + const Eigen::VectorXd & x, + const double & k, + Eigen::VectorXd & z) +{ + using namespace Eigen; + using namespace std; + + VectorXd tmp1 = x.array() - k; + VectorXd posMax = tmp1.array().max(0.0); + + VectorXd tmp2 = -x.array() - k; + VectorXd negMax = tmp2.array().max(0.0); + + // VectorXd posMax = x.array() - k; + // for (int ii=0; ii +#include + +void shrinkage( + const Eigen::VectorXd & x, + const double & k, + Eigen::VectorXd & z); +#endif // SHRINKAGE_H diff --git a/src/meshlabplugins/filter_cubization/curvdata.h b/src/meshlabplugins/filter_cubization/curvdata.h new file mode 100644 index 000000000..cb1d4a42f --- /dev/null +++ b/src/meshlabplugins/filter_cubization/curvdata.h @@ -0,0 +1,93 @@ +/**************************************************************************** +* VCGLib o o * +* Visual and Computer Graphics Library o o * +* _ O _ * +* Copyright(C) 2004 \/)\/ * +* Visual Computing Lab /\/| * +* ISTI - Italian National Research Council | * +* \ * +* All rights reserved. * +* * +* This program is free software; you can redistribute it and/or modify * +* it under the terms of the GNU General Public License as published by * +* the Free Software Foundation; either version 2 of the License, or * +* (at your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, * +* but WITHOUT ANY WARRANTY; without even the implied warranty of * +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) * +* for more details. * +* * +****************************************************************************/ +#ifndef __CURVDATA +#define __CURVDATA + +namespace vcg +{ +class CurvData +{ +public: + CurvData() + { + A = 0; // area + H = 0; // absolute mean curvature + K = 0; // gaussian curvature + } + + virtual ~CurvData() {} + + friend const CurvData operator+(const CurvData &lhs, const CurvData &rhs) + { + CurvData result; + result.A = lhs.A + rhs.A; + result.H = lhs.H + rhs.H; + result.K = lhs.K + rhs.K; + return result; + } + + friend CurvData &operator+=(CurvData &lhs, const CurvData &rhs) + { + lhs.A += rhs.A; + lhs.H += rhs.H; + lhs.K += rhs.K; + return lhs; + } + + float A; + float H; + float K; +}; + + +class NSMCEval +{ +public: + float operator() (const CurvData& c) + { + return (powf(c.H / 4.0, 2.0f) / c.A); + } +}; + +class MeanCEval +{ +public: + float operator() (const CurvData& c) + { + return (c.H / 4.0); + } +}; + +class AbsCEval +{ +public: + float operator() (const CurvData& c) + { + float k = 2 * M_PI - c.K; + if(k > 0) return 2.0 * (c.H / 4.0); + else return 2.0 * math::Sqrt(powf(c.H / 4.0, 2.0) - (c.A * k)); + } +}; + +} //end namespace vcg +#endif // __CURVDATA diff --git a/src/meshlabplugins/filter_cubization/curvedgeflip.h b/src/meshlabplugins/filter_cubization/curvedgeflip.h new file mode 100644 index 000000000..ebc988fb9 --- /dev/null +++ b/src/meshlabplugins/filter_cubization/curvedgeflip.h @@ -0,0 +1,339 @@ +/**************************************************************************** + * MeshLab o o * + * A versatile mesh processing toolbox o o * + * _ O _ * + * Copyright(C) 2005 \/)\/ * + * Visual Computing Lab /\/| * + * ISTI - Italian National Research Council | * + * \ * + * All rights reserved. * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License (http://www.gnu.org/licenses/gpl.txt) * + * for more details. * + * * + ****************************************************************************/ + +#ifndef __CURVEDGEFLIP +#define __CURVEDGEFLIP + +#include +#include +#include +#include + +#include + +#include "curvdata.h" + +namespace vcg +{ +namespace tri +{ + +/* + * This flip happens only if decreases the curvature of the surface + * Edge flip optimization based on the paper + * "optimizing 3d triangulations using discrete curvature analysis" + * http://www2.in.tu-clausthal.de/~hormann/papers/Dyn.2001.OTU.pdf + */ +template +class CurvEdgeFlip : public TopoEdgeFlip +{ +protected: + typedef typename TRIMESH_TYPE::FaceType FaceType; + typedef typename TRIMESH_TYPE::FacePointer FacePointer; + typedef typename TRIMESH_TYPE::FaceIterator FaceIterator; + typedef typename TRIMESH_TYPE::VertexIterator VertexIterator; + typedef typename TRIMESH_TYPE::VertexType VertexType; + typedef typename TRIMESH_TYPE::ScalarType ScalarType; + typedef typename TRIMESH_TYPE::VertexPointer VertexPointer; + typedef typename TRIMESH_TYPE::CoordType CoordType; + typedef vcg::face::Pos PosType; + typedef vcg::face::VFIterator VFIteratorType; + typedef typename LocalOptimization::HeapElem HeapElem; + typedef typename LocalOptimization::HeapType HeapType; + typedef TriEdgeFlip Parent; + typedef typename vcg::Triangle3 TriangleType; + + typedef typename TRIMESH_TYPE::VertContainer VertContainer; + + // New curvature precomputed for the vertices of the faces + // adjacent to edge to be flipped + ScalarType _cv0, _cv1, _cv2, _cv3; + + + static CurvData FaceCurv(VertexPointer v0, + VertexPointer v1, + VertexPointer v2, + CoordType fNormal) + { + CurvData res; + + float ang0 = math::Abs(Angle(v1->P() - v0->P(), v2->P() - v0->P())); + float ang1 = math::Abs(Angle(v0->P() - v1->P(), v2->P() - v1->P())); + float ang2 = M_PI - ang0 - ang1; + + float s01 = SquaredDistance(v1->P(), v0->P()); + float s02 = SquaredDistance(v2->P(), v0->P()); + + // voronoi cell of vertex + if (ang0 >= M_PI/2) { + TriangleType triangle(v0->P(), v1->P(), v2->P()); + res.A += (0.5f * DoubleArea(triangle) - + (s01 * tan(ang1) + s02 * tan(ang2)) / 8.0 ); + } else if (ang1 >= M_PI/2) + res.A += (s01 * tan(ang0)) / 8.0; + else if (ang2 >= M_PI/2) + res.A += (s02 * tan(ang0) / 8.0); + else // non obctuse triangle + res.A += ((s02 / tan(ang1)) + (s01 / tan(ang2))) / 8.0; + + res.K += ang0; + + ang1 = math::Abs(Angle(fNormal, v1->N())); + ang2 = math::Abs(Angle(fNormal, v2->N())); + res.H += ( (math::Sqrt(s01) / 2.0) * ang1 + + (math::Sqrt(s02) / 2.0) * ang2 ); + + return res; + } + + /* Compute vertex curvature walking around a verte with VF adjacency + * f1, f2 --> this faces are to be ignored + * */ + static CurvData Curvature(VertexPointer v, FacePointer f1 = NULL, FacePointer f2 = NULL) + { + CurvData curv; + VFIteratorType vfi(v); + + while (!vfi.End()) { + if (vfi.F() != f1 && vfi.F() != f2 && !vfi.F()->IsD()) { + int i = vfi.I(); + curv += FaceCurv(vfi.F()->V0(i), + vfi.F()->V1(i), + vfi.F()->V2(i), + vfi.F()->N()); + } + ++vfi; + } + + return curv; + } + +public: + CurvEdgeFlip() {} + + CurvEdgeFlip(PosType pos, int mark,BaseParameterClass *pp) + { + this->_pos = pos; + this->_localMark = mark; + this->_priority = ComputePriority(pp); + } + + + CurvEdgeFlip(CurvEdgeFlip &par) + { + this->_pos = par.GetPos(); + this->_localMark = par.GetMark(); + this->_priority = par.Priority(); + } + + // temporary empty (flip is already done in constructor) + void Execute(TRIMESH_TYPE& m, BaseParameterClass *) + { + int i = this->_pos.E(); + FacePointer f1 = this->_pos.F(); + int j = f1->FFi(i); + FacePointer f2 = f1->FFp(i); + VertexPointer v0 = f1->V0(i), v1 = f1->V1(i), + v2 = f1->V2(i), v3 = f2->V2(j); + + // save precomputed curvature in vertex quality + v0->Q() = _cv0; + v1->Q() = _cv1; + v2->Q() = _cv2; + v3->Q() = _cv3; + + CoordType n1 = Normal(v0->P(), v3->P(), v2->P()); + CoordType n2 = Normal(v1->P(), v2->P(), v3->P()); + + // compute and update normals on vertices after flip + v0->N() = v0->N() - f1->N() - f2->N() + n1; + v1->N() = v1->N() - f1->N() - f2->N() + n2; + v2->N() = v2->N() - f1->N() + n1 + n2; + v3->N() = v3->N() - f2->N() + n1 + n2; + + // fix VF adjacency - detach + assert(f1->V1(i) == v1); + vcg::face::VFDetach(*f1, (i + 1) % 3); // detach v1 from f1 + assert(f2->V1(j) == v0); + vcg::face::VFDetach(*f2, (j + 1) % 3); // detach v0 from f2 + + // do the flip + vcg::face::FlipEdge(*this->_pos.F(), this->_pos.E()); + + // fix VF adjacency - append + assert(f2->V1(j) == v2); + vcg::face::VFAppend(f2, (j + 1) % 3); // attach v2 with f2 + assert(f1->V1(i) == v3); + vcg::face::VFAppend(f1, (i + 1) % 3); // attach v3 with f1 + + // update face normals after the flip + f1->N() = n1; + f2->N() = n2; + + // avoid texture coordinates swap after flip + if (tri::HasPerWedgeTexCoord(m)) { + f2->WT((j + 1) % 3) = f1->WT((i + 2) % 3); + f1->WT((i + 1) % 3) = f2->WT((j + 2) % 3); + } + } + + + virtual bool IsFeasible(BaseParameterClass *_pp) + { + PlanarEdgeFlipParameter *pp=(PlanarEdgeFlipParameter *)_pp; + // First the flip must be topologically correct. + if(!vcg::face::CheckFlipEdge(*this->_pos.F(), this->_pos.E())) + return false; + + // then the angle between the involved normals must be greater??? + if (math::ToDeg(Angle(this->_pos.FFlip()->cN(), this->_pos.F()->cN()) ) <= pp->CoplanarAngleThresholdDeg ) + return false; + + CoordType v0, v1, v2, v3; + int i = this->_pos.E(); + FacePointer f = this->_pos.F(); + v0 = f->P0(i); + v1 = f->P1(i); + v2 = f->P2(i); + v3 = f->FFp(i)->P2(f->FFi(i)); + + // Take the parallelogram formed by the adjacent faces of edge + // If a corner of the parallelogram on extreme of edge to flip is >= 180 + // the flip will produce two identical faces - avoid this + if( (Angle(v2 - v0, v1 - v0) + Angle(v3 - v0, v1 - v0) >= M_PI) || + (Angle(v2 - v1, v0 - v1) + Angle(v3 - v1, v0 - v1) >= M_PI)) + return false; + + // if any of two faces adj to edge in non writable, the flip is unfeasible + if (!this->_pos.F()->IsW() || !this->_pos.F()->FFp(i)->IsW()) + return false; + + return true; + } + + + ScalarType ComputePriority(BaseParameterClass *pp) + { + /* + 1 + /|\ + / | \ + 2 | 3 + \ | / + \|/ + 0 + */ + + if(!this->IsFeasible(pp)) + return std::numeric_limits::infinity(); + + VertexPointer v0, v1, v2, v3; + int i = this->_pos.E(); + + FacePointer f1 = this->_pos.F(); + v0 = f1->V0(i); + v1 = f1->V1(i); + v2 = f1->V2(i); + + FacePointer f2 = this->_pos.F()->FFp(i); + v3 = f2->V2(f1->FFi(i)); + + // save sum of curvatures of vertices + float cbefore = v0->Q() + v1->Q() + v2->Q() + v3->Q(); + + // saving current vertex normals + CoordType nv0orig = v0->N(); + CoordType nv1orig = v1->N(); + CoordType nv2orig = v2->N(); + CoordType nv3orig = v3->N(); + + CurvData cd0, cd1, cd2, cd3; + CoordType n1 = Normal(v0->P(), v3->P(), v2->P()); + CoordType n2 = Normal(v1->P(), v2->P(), v3->P()); + + // new vertex normal after the flip + // set this now to evaluate new curvatures + v0->N() = nv0orig - f1->N() - f2->N() + n1; + v1->N() = nv1orig - f1->N() - f2->N() + n2; + v2->N() = nv2orig - f1->N() + n1 + n2; + v3->N() = nv3orig - f2->N() + n1 + n2; + + cd0 = FaceCurv(v0, v3, v2, n1) + Curvature(v0, f1, f2); + cd1 = FaceCurv(v1, v2, v3, n2) + Curvature(v1, f1, f2); + cd2 = FaceCurv(v2, v0, v3, n1) + FaceCurv(v2, v3, v1, n2) + Curvature(v2, f1, f2); + cd3 = FaceCurv(v3, v2, v0, n1) + FaceCurv(v3, v1, v2, n2) + Curvature(v3, f1, f2); + + // restore original pervertex normals after computation + v0->N() = nv0orig; + v1->N() = nv1orig; + v2->N() = nv2orig; + v3->N() = nv3orig; + + CURVEVAL curveval; + + _cv0 = curveval(cd0); + _cv1 = curveval(cd1); + _cv2 = curveval(cd2); + _cv3 = curveval(cd3); + float cafter = _cv0 + _cv1 + _cv2 + _cv3; + + // The priority of an edge flip is **how much we lower the overall curvature**. + // If after the flip the sum of the curvature is decreased it is a good move; + // good flips have: cafter < cbefore + // Since the local optimization is designed to make the minimum cost move we put inside + // negative values (the more negative the better). + this->_priority = (cafter - cbefore); + //qDebug("computed curvature change, %f->%f (priority = %f)", cbefore,cafter,this->_priority); + return this->_priority; + } + + + static void Init(TRIMESH_TYPE &m, HeapType &heap, BaseParameterClass *pp) + { + CURVEVAL curveval; + heap.clear(); + + // comuputing edge flip priority require non normalized vertex normals AND non normalized face normals. + vcg::tri::UpdateNormal::PerVertexPerFace(m); + + VertexIterator vi; + for (vi = m.vert.begin(); vi != m.vert.end(); ++vi) + if (!(*vi).IsD() && (*vi).IsW()) + (*vi).Q() = curveval(Curvature(&(*vi))); + + FaceIterator fi; + for (fi = m.face.begin(); fi != m.face.end(); ++fi) + if (!(*fi).IsD()) + for (unsigned int i = 0; i < 3; i++) + if ((*fi).V1(i) - (*fi).V0(i) > 0) { + PosType newpos(&*fi, i); + TopoEdgeFlip::Insert(heap, newpos, tri::IMark(m),pp); + } + } +}; // end CurvEdgeFlip class + + +}; // namespace tri +}; // namespace vcg + +#endif // __CURVEDGEFLIP diff --git a/src/meshlabplugins/filter_cubization/filter_cubization.cpp b/src/meshlabplugins/filter_cubization/filter_cubization.cpp new file mode 100644 index 000000000..decacb751 --- /dev/null +++ b/src/meshlabplugins/filter_cubization/filter_cubization.cpp @@ -0,0 +1,310 @@ +/**************************************************************************** +* MeshLab o o * +* A versatile mesh processing toolbox o o * +* _ O _ * +* Copyright(C) 2005 \/)\/ * +* Visual Computing Lab /\/| * +* ISTI - Italian National Research Council | * +* \ * +* All rights reserved. * +* * +* This program is free software; you can redistribute it and/or modify * +* it under the terms of the GNU General Public License as published by * +* the Free Software Foundation; either version 2 of the License, or * +* (at your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, * +* but WITHOUT ANY WARRANTY; without even the implied warranty of * +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) * +* for more details. * +* * +****************************************************************************/ + +#include +#include +#include + +#include "filter_cubization.h" + +#include +#include + +#include +#include + +#include +#include + +#include "curvedgeflip.h" +#include "curvdata.h" + +using namespace vcg; + +// forward declarations +class MeanCEFlip; +typedef Histogram Histogramm; + +class MeanCEFlip : public vcg::tri::CurvEdgeFlip +{ +public: + MeanCEFlip(PosType pos, int mark,BaseParameterClass *_pp) : + vcg::tri::CurvEdgeFlip(pos, mark,_pp) {} +}; + + +// Constructor usually performs only two simple tasks of filling the two lists +// - typeList: with all the possible id of the filtering actions +// - actionList with the corresponding actions. If you want to add icons to +// your filtering actions you can do here by construction the QActions accordingly +CubizationPlugin::CubizationPlugin() +{ + typeList = { + FP_CUBIZATION, + FP_CUBIZATION_COLORIZE + }; + + for(ActionIDType tt: types()) + actionList.push_back(new QAction(filterName(tt), this)); + + cubic_ApplyEdgeFlip = false; +} + +QString CubizationPlugin::pluginName() const +{ + return "FilterCubization"; +} + +QString CubizationPlugin::pythonFilterName(ActionIDType f) const +{ + switch (f) { + case FP_CUBIZATION: return tr("applying_cubic_stylization"); + case FP_CUBIZATION_COLORIZE: return tr("applying_cubization_filter_colorizing_vertices"); + default: assert(0); return QString(); + } +} + +QString CubizationPlugin::filterName(ActionIDType filterId) const +{ + switch (filterId) { + case FP_CUBIZATION: return tr("Cubic stylization"); + case FP_CUBIZATION_COLORIZE: return tr("Cubic stylization and colorize"); + default: assert(0); return QString(); + } +} + +FilterPlugin::FilterArity CubizationPlugin::filterArity(const QAction*) const +{ + return SINGLE_MESH; +} + + +int CubizationPlugin::getPreConditions(const QAction*) const +{ + return MeshModel::MM_VERTCOORD | MeshModel::MM_FACEVERT; +} + +int CubizationPlugin::getRequirements(const QAction *action) +{ + switch (ID(action)) { + case FP_CUBIZATION: + case FP_CUBIZATION_COLORIZE: + return MeshModel::MM_VERTCOORD | + MeshModel::MM_FACEFACETOPO | + MeshModel::MM_VERTFACETOPO | + MeshModel::MM_VERTMARK; + } + + return 0; +} + +QString CubizationPlugin::filterInfo(ActionIDType filterId) const +{ + switch(filterId) { + case FP_CUBIZATION: + return tr("Cubic stylization of the mesh. For all detailed about cubic stylization see:
Hsueh-Ti Derek Liu and Alec Jacobson, 'Cubic Stylization', ACM Transactions on Graphics 2019
"); + case FP_CUBIZATION_COLORIZE: + return tr("Colorize the vertices of a mesh with cubic stylization. For all detailed about cubic stylization see:
Hsueh-Ti Derek Liu and Alec Jacobson, 'Cubic Stylization', ACM Transactions on Graphics 2019
"); + default : assert(0); + } + return {}; +} + + CubizationPlugin::FilterClass CubizationPlugin::getClass(const QAction *action) const +{ + switch(ID(action)) { + case FP_CUBIZATION: return FilterPlugin::Remeshing; + case FP_CUBIZATION_COLORIZE: return FilterPlugin::Remeshing; + } + return FilterPlugin::Generic; +} + +int CubizationPlugin::postCondition(const QAction *a) const +{ + switch(ID(a)) + { + case FP_CUBIZATION : return MeshModel::MM_ALL; + case FP_CUBIZATION_COLORIZE : return MeshModel::MM_ALL | + MeshModel::MM_VERTQUALITY; + default : assert(0); + } + return {}; +} + +// This function define the needed parameters for each filter. +// Return true if the filter has some parameters +// it is called every time, so you can set the default value of parameters according to the mesh +// For each parameter you need to define, +// - the name of the parameter, +// - the string shown in the dialog +// - the default value +// - a possibly long string describing the meaning of that parameter (shown as a popup help in the dialog) +RichParameterList CubizationPlugin::initParameterList(const QAction *action, const MeshModel &m) +{ + RichParameterList parlst; + if (ID(action) == FP_CUBIZATION) { + parlst.addParam(RichFloat("lcubeness", 0.2f, + tr("Cubeness parameter (λ)"), + tr("Control the cubeness of the mesh. Generally, the higher the cubeness parameter, the more cubic the mesh is. λ ∈ [0, 1] "))); + + + parlst.addParam(RichBool("applyef", cubic_ApplyEdgeFlip, tr("Apply edge flipping"), tr("Apply edge flip optimization on cubic stylization."))); + + } + if (ID(action) == FP_CUBIZATION_COLORIZE) { + parlst.addParam(RichFloat("lcubeness", 0.2f, + tr("Cubeness parameter (λ)"), + tr("Control the cubeness of the mesh. Generally, the higher the cubeness parameter, the more cubic the mesh is. λ ∈ [0, 1] "))); + + + parlst.addParam(RichBool("applyef", cubic_ApplyEdgeFlip, tr("Apply edge flipping"), tr("Apply edge flip optimization on cubic stylization."))); + } + return parlst; +} + +// The Real Core Function doing the actual mesh processing. +// Run mesh optimization +std::map CubizationPlugin::applyFilter( + const QAction *filter, + const RichParameterList & par, + MeshDocument &md, + unsigned int& /*postConditionMask*/, + vcg::CallBackPos *cb) +{ + + MeshModel &m=*(md.mm()); + double energyTotal = 0.f; + time_t start = clock(); + + if (ID(filter) == FP_CUBIZATION) { + m = ComputeCubicStylization(md, par, energyTotal); + } + + if (ID(filter) == FP_CUBIZATION_COLORIZE) { + m = ComputeCubicStylization(md, par, energyTotal, true); + + m.updateDataMask(MeshModel::MM_FACEFACETOPO); + m.updateDataMask(MeshModel::MM_VERTCOLOR); + m.updateDataMask(MeshModel::MM_VERTQUALITY); + + Histogramm H; + vcg::tri::Stat::ComputePerVertexQualityHistogram(m.cm, H); + vcg::tri::UpdateColor::PerVertexQualityRamp( + m.cm, H.Percentile(0.01f), H.Percentile(0.99f)); + } + + m.updateBoxAndNormals(); + + log( "cubic stylization performed in %.2f sec. with cubic energy equal to %.5f", (clock() - start) / (float) CLOCKS_PER_SEC, energyTotal); + return std::map(); +} + +MeshModel CubizationPlugin::ComputeCubicStylization( + MeshDocument& md, + const RichParameterList& par, + double& totalEnergy, + bool isColorizing){ + + MeshModel &m=*(md.mm()); + float limit = -std::numeric_limits::epsilon(); + + int delvert = tri::Clean::RemoveUnreferencedVertex(m.cm); + if (delvert) + log( "Pre-Curvature Cleaning: Removed %d unreferenced vertices", + delvert); + + + tri::Allocator::CompactVertexVector(m.cm); + tri::Allocator::CompactFaceVector(m.cm); + m.updateDataMask(MeshModel::MM_FACEFACETOPO); + vcg::tri::UpdateFlags::FaceBorderFromFF(m.cm); + + bool isApplyEdgeFlip = par.getBool("applyef"); + + if ( tri::Clean::CountNonManifoldEdgeFF(m.cm) >0 && isApplyEdgeFlip) { + throw MLException("Mesh has some not 2-manifold faces, edge flips requires manifoldness"); + } + + Eigen::MatrixXd verts; + Eigen::MatrixXi faces; + Eigen::MatrixXd u_verts; + Eigen::VectorXd energy_verts; + + float lambda = par.getFloat("lcubeness"); + cube_style_data data; + + data.lambda = lambda; + + vcg::tri::Cubization::Init(m.cm, verts, u_verts, faces, data); + + // apply cubic stylization + int maxIter = 1000; + double stopReldV = 1e-3; // stopping criteria for relative displacement + double reldV = 0; + + for (int iter=0; iter::Stylize(verts, u_verts, faces, data, energy_verts, totalEnergy); + + //apply Edge Flips + if((iter%30 == 0 || reldV < stopReldV) && isApplyEdgeFlip){ + Matrix2Mesh(m.cm, u_verts, faces); + + vcg::tri::PlanarEdgeFlipParameter pp; + float limit = -std::numeric_limits::epsilon(); + + vcg::LocalOptimization optimiz(m.cm,&pp); + float pthr = 10; + + // VF adjacency needed for edge flips based on vertex curvature + vcg::tri::UpdateTopology::VertexFace(m.cm); + vcg::tri::UpdateTopology::TestVertexFace(m.cm); + + pp.CoplanarAngleThresholdDeg = pthr; + optimiz.Init(); + + // stop when flips become harmful + optimiz.SetTargetMetric(limit); + optimiz.DoOptimization(); + optimiz.h.clear(); + + Mesh2Matrix(m.cm, u_verts, faces); + } + + if (reldV < stopReldV) break; + } + + Matrix2Mesh(m.cm, u_verts, faces); + + if(isColorizing){ + assert(energy_verts.size() == m.cm.vert.size()); + for(int i = 0; i < energy_verts.size(); i++){ + m.cm.vert[i].Q() = energy_verts[i]; + } + + } + + return m; +} + +MESHLAB_PLUGIN_NAME_EXPORTER(CubizationPlugin) diff --git a/src/meshlabplugins/filter_cubization/filter_cubization.h b/src/meshlabplugins/filter_cubization/filter_cubization.h new file mode 100644 index 000000000..a3dcb4e2f --- /dev/null +++ b/src/meshlabplugins/filter_cubization/filter_cubization.h @@ -0,0 +1,73 @@ +/**************************************************************************** +* MeshLab o o * +* A versatile mesh processing toolbox o o * +* _ O _ * +* Copyright(C) 2005 \/)\/ * +* Visual Computing Lab /\/| * +* ISTI - Italian National Research Council | * +* \ * +* All rights reserved. * +* * +* This program is free software; you can redistribute it and/or modify * +* it under the terms of the GNU General Public License as published by * +* the Free Software Foundation; either version 2 of the License, or * +* (at your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, * +* but WITHOUT ANY WARRANTY; without even the implied warranty of * +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) * +* for more details. * +* * +****************************************************************************/ + +#ifndef FILTER_CUBIZATION_H +#define FILTER_CUBIZATION_H + +#include +#include + +class CubizationPlugin : public QObject, public FilterPlugin +{ + Q_OBJECT + MESHLAB_PLUGIN_IID_EXPORTER(FILTER_PLUGIN_IID) + Q_INTERFACES(FilterPlugin) + +public: + enum { + // mesh improvement by edge flipping + FP_CUBIZATION, + FP_CUBIZATION_COLORIZE + }; + + CubizationPlugin(); + + QString pythonFilterName(ActionIDType f) const; + QString pluginName() const; + QString filterName(ActionIDType filter) const; + QString filterInfo(ActionIDType filter) const; + RichParameterList initParameterList(const QAction*, const MeshModel &/*m*/); + std::map applyFilter( + const QAction* action, + const RichParameterList & parameters, + MeshDocument &md, + unsigned int& postConditionMask, + vcg::CallBackPos * cb); + int getPreConditions(const QAction *) const; + int getRequirements(const QAction*); + FilterClass getClass(const QAction *) const; + int postCondition(const QAction* ) const; + FilterArity filterArity(const QAction*) const; + +private: + MeshModel ComputeCubicStylization( + MeshDocument& md, + const RichParameterList& par, + double& totalEnergy, + bool isColorizing = false); + +protected: + bool cubic_ApplyEdgeFlip; +}; + +#endif // FILTER_CUBIZATION_H