From 4edf5dede5a4430b86d4dc65705c59e643d955df Mon Sep 17 00:00:00 2001 From: Zhongshi Date: Sun, 14 Oct 2018 17:55:36 -0400 Subject: [PATCH] Adding 2D version of SCAF for Bijective Maps and a tutorial (#752) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Add scaf code Add the tutorial entry Tutorial entry for SCAF. still need polishing Put -4 in the energy computation of scaf:SymmD * :school: Insert SCAF to 711 Fixes in tutorial.md updating tutorial.html Fix for gcc 4.8 Remove 4.1 Refactoring scaf WIP Further refactoring WIP Code Cleaning * 🔨 Improve SLIM interface * :recycle: refactor SCAF with less lines * :recycle: refactor harmonic with boundary * :sparkles: Add direct signature for scaf * :memo: Change tutorial order wip * :school: Timer in SLIM tutorial makes more sense * :sparkles: Add direct interface for SLIM * :art: format SCAF code. * :snake: Python binding for slim and scaf * :snake: expose SLIMData (maybe elevate its status more in the future) * :alien: :snake: deprecate PYBIND11_PLUGIN * :sparkles: harmonic with holes * :recycle: SCAF, Decouple Initialization to tutorial * :school: SCAF tutorial improvement * :school: Tutorial order; insert scaf * :art: Removing Adjusted_grad * :school: Inserted SCAF, not messing around with tutorials order until the new page is online * :green_heart: Fix for static build * :recycle: Uniform SLIM and SCAF and python binding canonical * :bug: Change in ScafData * :apple: Fix namespace for build * :sparkles: Introduced Mapping Energy Type to replace adhoc SLIMEnergy * :recycle: Refactor subroutines from SLIM * Revert python changes * leave fill hole outside harmonic * Revert Python changes * hole fill * Minor Fixes for SCAF SCAF->IGL_SCAF Topological_hole_fill "" in slim tutorial --- include/igl/MappingEnergyType.h | 26 + include/igl/mapping_energy_with_jacobians.cpp | 141 + include/igl/mapping_energy_with_jacobians.h | 36 + include/igl/scaf.cpp | 707 ++ include/igl/scaf.h | 101 + include/igl/slice.cpp | 1 + include/igl/slim.cpp | 1004 +-- include/igl/slim.h | 30 +- include/igl/sparse_cached.h | 2 +- include/igl/topological_hole_fill.cpp | 55 + include/igl/topological_hole_fill.h | 44 + .../{710_SLIM => 709_SLIM}/CMakeLists.txt | 0 tutorial/{710_SLIM => 709_SLIM}/main.cpp | 33 +- tutorial/710_SCAF/CMakeLists.txt | 5 + tutorial/710_SCAF/main.cpp | 127 + tutorial/CMakeLists.txt | 3 +- tutorial/images/710_SCAF.png | Bin 0 -> 167471 bytes tutorial/shared/camel_b.obj | 7640 +++++++++++++++++ 18 files changed, 9347 insertions(+), 608 deletions(-) create mode 100644 include/igl/MappingEnergyType.h create mode 100644 include/igl/mapping_energy_with_jacobians.cpp create mode 100644 include/igl/mapping_energy_with_jacobians.h create mode 100644 include/igl/scaf.cpp create mode 100644 include/igl/scaf.h create mode 100644 include/igl/topological_hole_fill.cpp create mode 100644 include/igl/topological_hole_fill.h rename tutorial/{710_SLIM => 709_SLIM}/CMakeLists.txt (100%) rename tutorial/{710_SLIM => 709_SLIM}/main.cpp (92%) create mode 100644 tutorial/710_SCAF/CMakeLists.txt create mode 100755 tutorial/710_SCAF/main.cpp create mode 100644 tutorial/images/710_SCAF.png create mode 100644 tutorial/shared/camel_b.obj diff --git a/include/igl/MappingEnergyType.h b/include/igl/MappingEnergyType.h new file mode 100644 index 000000000..acf2dc0e5 --- /dev/null +++ b/include/igl/MappingEnergyType.h @@ -0,0 +1,26 @@ +// This file is part of libigl, a simple c++ geometry processing library. +// +// Copyright (C) 2018 Zhongshi Jiang +// +// 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/. +#ifndef IGL_MAPPINGENERGYTYPE_H +#define IGL_MAPPINGENERGYTYPE_H +namespace igl +{ + // Energy Types used for Parameterization/Mapping. + // Refer to SLIM [Rabinovich et al. 2017] for more details + // Todo: Integrate with ARAPEnergyType + + enum MappingEnergyType + { + ARAP, + LOG_ARAP, + SYMMETRIC_DIRICHLET, + CONFORMAL, + EXP_CONFORMAL, + EXP_SYMMETRIC_DIRICHLET + }; +} +#endif diff --git a/include/igl/mapping_energy_with_jacobians.cpp b/include/igl/mapping_energy_with_jacobians.cpp new file mode 100644 index 000000000..791e01f89 --- /dev/null +++ b/include/igl/mapping_energy_with_jacobians.cpp @@ -0,0 +1,141 @@ +// This file is part of libigl, a simple c++ geometry processing library. +// +// Copyright (C) 2018 Zhongshi Jiang +// +// 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/. + +#include "mapping_energy_with_jacobians.h" +#include "polar_svd.h" + +IGL_INLINE double igl::mapping_energy_with_jacobians( + const Eigen::MatrixXd &Ji, + const Eigen::VectorXd &areas, + igl::MappingEnergyType slim_energy, + double exp_factor){ + + double energy = 0; + if (Ji.cols() == 4) + { + Eigen::Matrix ji; + for (int i = 0; i < Ji.rows(); i++) + { + ji(0, 0) = Ji(i, 0); + ji(0, 1) = Ji(i, 1); + ji(1, 0) = Ji(i, 2); + ji(1, 1) = Ji(i, 3); + + typedef Eigen::Matrix Mat2; + typedef Eigen::Matrix Vec2; + Mat2 ri, ti, ui, vi; + Vec2 sing; + igl::polar_svd(ji, ri, ti, ui, sing, vi); + double s1 = sing(0); + double s2 = sing(1); + + switch (slim_energy) + { + case igl::MappingEnergyType::ARAP: + { + energy += areas(i) * (pow(s1 - 1, 2) + pow(s2 - 1, 2)); + break; + } + case igl::MappingEnergyType::SYMMETRIC_DIRICHLET: + { + energy += areas(i) * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2)); + break; + } + case igl::MappingEnergyType::EXP_SYMMETRIC_DIRICHLET: + { + energy += areas(i) * exp(exp_factor * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2))); + break; + } + case igl::MappingEnergyType::LOG_ARAP: + { + energy += areas(i) * (pow(log(s1), 2) + pow(log(s2), 2)); + break; + } + case igl::MappingEnergyType::CONFORMAL: + { + energy += areas(i) * ((pow(s1, 2) + pow(s2, 2)) / (2 * s1 * s2)); + break; + } + case igl::MappingEnergyType::EXP_CONFORMAL: + { + energy += areas(i) * exp(exp_factor * ((pow(s1, 2) + pow(s2, 2)) / (2 * s1 * s2))); + break; + } + + } + + } + } + else + { + Eigen::Matrix ji; + for (int i = 0; i < Ji.rows(); i++) + { + ji(0, 0) = Ji(i, 0); + ji(0, 1) = Ji(i, 1); + ji(0, 2) = Ji(i, 2); + ji(1, 0) = Ji(i, 3); + ji(1, 1) = Ji(i, 4); + ji(1, 2) = Ji(i, 5); + ji(2, 0) = Ji(i, 6); + ji(2, 1) = Ji(i, 7); + ji(2, 2) = Ji(i, 8); + + typedef Eigen::Matrix Mat3; + typedef Eigen::Matrix Vec3; + Mat3 ri, ti, ui, vi; + Vec3 sing; + igl::polar_svd(ji, ri, ti, ui, sing, vi); + double s1 = sing(0); + double s2 = sing(1); + double s3 = sing(2); + + switch (slim_energy) + { + case igl::MappingEnergyType::ARAP: + { + energy += areas(i) * (pow(s1 - 1, 2) + pow(s2 - 1, 2) + pow(s3 - 1, 2)); + break; + } + case igl::MappingEnergyType::SYMMETRIC_DIRICHLET: + { + energy += areas(i) * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2) + pow(s3, 2) + pow(s3, -2)); + break; + } + case igl::MappingEnergyType::EXP_SYMMETRIC_DIRICHLET: + { + energy += areas(i) * exp(exp_factor * + (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2) + pow(s3, 2) + pow(s3, -2))); + break; + } + case igl::MappingEnergyType::LOG_ARAP: + { + energy += areas(i) * (pow(log(s1), 2) + pow(log(std::abs(s2)), 2) + pow(log(std::abs(s3)), 2)); + break; + } + case igl::MappingEnergyType::CONFORMAL: + { + energy += areas(i) * ((pow(s1, 2) + pow(s2, 2) + pow(s3, 2)) / (3 * pow(s1 * s2 * s3, 2. / 3.))); + break; + } + case igl::MappingEnergyType::EXP_CONFORMAL: + { + energy += areas(i) * exp((pow(s1, 2) + pow(s2, 2) + pow(s3, 2)) / (3 * pow(s1 * s2 * s3, 2. / 3.))); + break; + } + } + } + } + + return energy; +} + + +#ifdef IGL_STATIC_LIBRARY +// Explicit template instantiation +#endif diff --git a/include/igl/mapping_energy_with_jacobians.h b/include/igl/mapping_energy_with_jacobians.h new file mode 100644 index 000000000..2952db60d --- /dev/null +++ b/include/igl/mapping_energy_with_jacobians.h @@ -0,0 +1,36 @@ +// This file is part of libigl, a simple c++ geometry processing library. +// +// Copyright (C) 2018 Zhongshi Jiang +// +// 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/. +#ifndef IGL_MAPPING_ENERGY_WITH_JACOBIANS_H +#define IGL_MAPPING_ENERGY_WITH_JACOBIANS_H + +#include "igl_inline.h" +#include +#include "MappingEnergyType.h" + +namespace igl +{ + // compute the rotation-invariant energy of a mapping (represented in Jacobians and areas) + // Input: + // Ji: #F by 4 (9 if 3D) entries of jacobians + // areas: #F by 1 face areas + // slim_energy: energy type as in igl::MappingEnergyType + // exp_factor: see igl::MappingEnergyType + // + // Output: + // energy value + IGL_INLINE double mapping_energy_with_jacobians(const Eigen::MatrixXd &Ji, + const Eigen::VectorXd &areas, + igl::MappingEnergyType slim_energy, + double exp_factor); + +} +#ifndef IGL_STATIC_LIBRARY +# include "mapping_energy_with_jacobians.cpp" +#endif + +#endif \ No newline at end of file diff --git a/include/igl/scaf.cpp b/include/igl/scaf.cpp new file mode 100644 index 000000000..a3f969a9c --- /dev/null +++ b/include/igl/scaf.cpp @@ -0,0 +1,707 @@ +// This file is part of libigl, a simple c++ geometry processing library. +// +// Copyright (C) 2018 Zhongshi Jiang +// +// 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/. + +#include "scaf.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "mapping_energy_with_jacobians.h" + +#include +#include +#include +#include +#include +namespace igl +{ + +namespace scaf +{ +void update_scaffold(igl::SCAFData &s) +{ + s.mv_num = s.m_V.rows(); + s.mf_num = s.m_T.rows(); + + s.v_num = s.w_uv.rows(); + s.sf_num = s.s_T.rows(); + + s.sv_num = s.v_num - s.mv_num; + s.f_num = s.sf_num + s.mf_num; + + s.s_M = Eigen::VectorXd::Constant(s.sf_num, s.scaffold_factor); +} + +void adjusted_grad(Eigen::MatrixXd &V, + Eigen::MatrixXi &F, + double area_threshold, + Eigen::SparseMatrix &Dx, + Eigen::SparseMatrix &Dy, + Eigen::SparseMatrix &Dz) +{ + Eigen::VectorXd M; + igl::doublearea(V, F, M); + std::vector degen; + for (int i = 0; i < M.size(); i++) + if (M(i) < area_threshold) + degen.push_back(i); + + Eigen::SparseMatrix G; + igl::grad(V, F, G); + + Dx = G.topRows(F.rows()); + Dy = G.block(F.rows(), 0, F.rows(), V.rows()); + Dz = G.bottomRows(F.rows()); + + // handcraft uniform gradient for faces area falling below threshold. + double sin60 = std::sin(M_PI / 3); + double cos60 = std::cos(M_PI / 3); + double deno = std::sqrt(sin60 * area_threshold); + Eigen::MatrixXd standard_grad(3, 3); + standard_grad << -sin60 / deno, sin60 / deno, 0, + -cos60 / deno, -cos60 / deno, 1 / deno, + 0, 0, 0; + + for (auto k : degen) + for (int j = 0; j < 3; j++) + { + Dx.coeffRef(k, F(k, j)) = standard_grad(0, j); + Dy.coeffRef(k, F(k, j)) = standard_grad(1, j); + Dz.coeffRef(k, F(k, j)) = standard_grad(2, j); + } +} + +void compute_scaffold_gradient_matrix(SCAFData &s, + Eigen::SparseMatrix &D1, + Eigen::SparseMatrix &D2) +{ + using namespace Eigen; + Eigen::SparseMatrix G; + MatrixXi F_s = s.s_T; + int vn = s.v_num; + MatrixXd V = MatrixXd::Zero(vn, 3); + V.leftCols(2) = s.w_uv; + + double min_bnd_edge_len = INFINITY; + int acc_bnd = 0; + for (int i = 0; i < s.bnd_sizes.size(); i++) + { + int current_size = s.bnd_sizes[i]; + + for (int e = acc_bnd; e < acc_bnd + current_size - 1; e++) + { + min_bnd_edge_len = (std::min)(min_bnd_edge_len, + (s.w_uv.row(s.internal_bnd(e)) - + s.w_uv.row(s.internal_bnd(e + 1))) + .squaredNorm()); + } + min_bnd_edge_len = (std::min)(min_bnd_edge_len, + (s.w_uv.row(s.internal_bnd(acc_bnd)) - + s.w_uv.row(s.internal_bnd(acc_bnd + current_size - 1))) + .squaredNorm()); + acc_bnd += current_size; + } + + double area_threshold = min_bnd_edge_len / 4.0; + Eigen::SparseMatrix Dx, Dy, Dz; + adjusted_grad(V, F_s, area_threshold, Dx, Dy, Dz); + + MatrixXd F1, F2, F3; + igl::local_basis(V, F_s, F1, F2, F3); + D1 = F1.col(0).asDiagonal() * Dx + F1.col(1).asDiagonal() * Dy + + F1.col(2).asDiagonal() * Dz; + D2 = F2.col(0).asDiagonal() * Dx + F2.col(1).asDiagonal() * Dy + + F2.col(2).asDiagonal() * Dz; +} + +void mesh_improve(igl::SCAFData &s) +{ + using namespace Eigen; + MatrixXd m_uv = s.w_uv.topRows(s.mv_num); + MatrixXd V_bnd; + V_bnd.resize(s.internal_bnd.size(), 2); + for (int i = 0; i < s.internal_bnd.size(); i++) // redoing step 1. + { + V_bnd.row(i) = m_uv.row(s.internal_bnd(i)); + } + + if (s.rect_frame_V.size() == 0) + { + Matrix2d ob; // = rect_corners; + { + VectorXd uv_max = m_uv.colwise().maxCoeff(); + VectorXd uv_min = m_uv.colwise().minCoeff(); + VectorXd uv_mid = (uv_max + uv_min) / 2.; + + Eigen::Array2d scaf_range(3, 3); + ob.row(0) = uv_mid.array() + scaf_range * ((uv_min - uv_mid).array()); + ob.row(1) = uv_mid.array() + scaf_range * ((uv_max - uv_mid).array()); + } + Vector2d rect_len; + rect_len << ob(1, 0) - ob(0, 0), ob(1, 1) - ob(0, 1); + int frame_points = 5; + + s.rect_frame_V.resize(4 * frame_points, 2); + for (int i = 0; i < frame_points; i++) + { + // 0,0;0,1 + s.rect_frame_V.row(i) << ob(0, 0), ob(0, 1) + i * rect_len(1) / frame_points; + // 0,0;1,1 + s.rect_frame_V.row(i + frame_points) + << ob(0, 0) + i * rect_len(0) / frame_points, + ob(1, 1); + // 1,0;1,1 + s.rect_frame_V.row(i + 2 * frame_points) << ob(1, 0), ob(1, 1) - i * rect_len(1) / frame_points; + // 1,0;0,1 + s.rect_frame_V.row(i + 3 * frame_points) + << ob(1, 0) - i * rect_len(0) / frame_points, + ob(0, 1); + // 0,0;0,1 + } + s.frame_ids = Eigen::VectorXi::LinSpaced(s.rect_frame_V.rows(), s.mv_num, s.mv_num + s.rect_frame_V.rows()); + } + + // Concatenate Vert and Edge + MatrixXd V; + MatrixXi E; + igl::cat(1, V_bnd, s.rect_frame_V, V); + E.resize(V.rows(), 2); + for (int i = 0; i < E.rows(); i++) + E.row(i) << i, i + 1; + int acc_bs = 0; + for (auto bs : s.bnd_sizes) + { + E(acc_bs + bs - 1, 1) = acc_bs; + acc_bs += bs; + } + E(V.rows() - 1, 1) = acc_bs; + assert(acc_bs == s.internal_bnd.size()); + + MatrixXd H = MatrixXd::Zero(s.component_sizes.size(), 2); + { + int hole_f = 0; + int hole_i = 0; + for (auto cs : s.component_sizes) + { + for (int i = 0; i < 3; i++) + H.row(hole_i) += m_uv.row(s.m_T(hole_f, i)); // redoing step 2 + hole_f += cs; + hole_i++; + } + } + H /= 3.; + + MatrixXd uv2; + igl::triangle::triangulate(V, E, H, std::basic_string("qYYQ"), uv2, s.s_T); + auto bnd_n = s.internal_bnd.size(); + + for (auto i = 0; i < s.s_T.rows(); i++) + for (auto j = 0; j < s.s_T.cols(); j++) + { + auto &x = s.s_T(i, j); + if (x < bnd_n) + x = s.internal_bnd(x); + else + x += m_uv.rows() - bnd_n; + } + + igl::cat(1, s.m_T, s.s_T, s.w_T); + s.w_uv.conservativeResize(m_uv.rows() - bnd_n + uv2.rows(), 2); + s.w_uv.bottomRows(uv2.rows() - bnd_n) = uv2.bottomRows(-bnd_n + uv2.rows()); + + update_scaffold(s); + + // after_mesh_improve + compute_scaffold_gradient_matrix(s, s.Dx_s, s.Dy_s); + + s.Dx_s.makeCompressed(); + s.Dy_s.makeCompressed(); + s.Dz_s.makeCompressed(); + s.Ri_s = MatrixXd::Zero(s.Dx_s.rows(), s.dim * s.dim); + s.Ji_s.resize(s.Dx_s.rows(), s.dim * s.dim); + s.W_s.resize(s.Dx_s.rows(), s.dim * s.dim); +} + +void add_new_patch(igl::SCAFData &s, const Eigen::MatrixXd &V_ref, + const Eigen::MatrixXi &F_ref, + const Eigen::RowVectorXd ¢er, + const Eigen::MatrixXd &uv_init) +{ + using namespace std; + using namespace Eigen; + + assert(uv_init.rows() != 0); + Eigen::VectorXd M; + igl::doublearea(V_ref, F_ref, M); + s.mesh_measure += M.sum() / 2; + + Eigen::VectorXi bnd; + Eigen::MatrixXd bnd_uv; + + std::vector> all_bnds; + igl::boundary_loop(F_ref, all_bnds); + int num_holes = all_bnds.size() - 1; + + s.component_sizes.push_back(F_ref.rows()); + + MatrixXd m_uv = s.w_uv.topRows(s.mv_num); + igl::cat(1, m_uv, uv_init, s.w_uv); + + s.m_M.conservativeResize(s.mf_num + M.size()); + s.m_M.bottomRows(M.size()) = M / 2; + + for (auto cur_bnd : all_bnds) + { + s.internal_bnd.conservativeResize(s.internal_bnd.size() + cur_bnd.size()); + s.internal_bnd.bottomRows(cur_bnd.size()) = Map(cur_bnd.data(), cur_bnd.size()) + s.mv_num; + s.bnd_sizes.push_back(cur_bnd.size()); + } + + s.m_T.conservativeResize(s.mf_num + F_ref.rows(), 3); + s.m_T.bottomRows(F_ref.rows()) = F_ref.array() + s.mv_num; + s.mf_num += F_ref.rows(); + + s.m_V.conservativeResize(s.mv_num + V_ref.rows(), 3); + s.m_V.bottomRows(V_ref.rows()) = V_ref; + s.mv_num += V_ref.rows(); + + s.rect_frame_V = MatrixXd(); + + mesh_improve(s); +} + +void compute_jacobians(SCAFData &s, const Eigen::MatrixXd &V_new, bool whole) +{ + auto comp_J2 = [](const Eigen::MatrixXd &uv, + const Eigen::SparseMatrix &Dx, + const Eigen::SparseMatrix &Dy, + Eigen::MatrixXd &Ji) { + // Ji=[D1*u,D2*u,D1*v,D2*v]; + Ji.resize(Dx.rows(), 4); + Ji.col(0) = Dx * uv.col(0); + Ji.col(1) = Dy * uv.col(0); + Ji.col(2) = Dx * uv.col(1); + Ji.col(3) = Dy * uv.col(1); + }; + + Eigen::MatrixXd m_V_new = V_new.topRows(s.mv_num); + comp_J2(m_V_new, s.Dx_m, s.Dy_m, s.Ji_m); + if (whole) + comp_J2(V_new, s.Dx_s, s.Dy_s, s.Ji_s); +} + +double compute_energy_from_jacobians(const Eigen::MatrixXd &Ji, + const Eigen::VectorXd &areas, + igl::MappingEnergyType energy_type) +{ + double energy = 0; + if (energy_type == igl::MappingEnergyType::SYMMETRIC_DIRICHLET) + energy = -4; // comply with paper description + return energy + igl::mapping_energy_with_jacobians(Ji, areas, energy_type, 0); +} + +double compute_soft_constraint_energy(const SCAFData &s) +{ + double e = 0; + for (auto const &x : s.soft_cons) + e += s.soft_const_p * (x.second - s.w_uv.row(x.first)).squaredNorm(); + + return e; +} + +double compute_energy(SCAFData &s, Eigen::MatrixXd &w_uv, bool whole) +{ + if (w_uv.rows() != s.v_num) + assert(!whole); + compute_jacobians(s, w_uv, whole); + double energy = compute_energy_from_jacobians(s.Ji_m, s.m_M, s.slim_energy); + + if (whole) + energy += compute_energy_from_jacobians(s.Ji_s, s.s_M, s.scaf_energy); + energy += compute_soft_constraint_energy(s); + return energy; +} + +void buildAm(const Eigen::VectorXd &sqrt_M, + const Eigen::SparseMatrix &Dx, + const Eigen::SparseMatrix &Dy, + const Eigen::MatrixXd &W, + Eigen::SparseMatrix &Am) +{ + std::vector> IJV; + Eigen::SparseMatrix Dz; + + Eigen::SparseMatrix MDx = sqrt_M.asDiagonal() * Dx; + Eigen::SparseMatrix MDy = sqrt_M.asDiagonal() * Dy; + igl::slim_buildA(MDx, MDy, Dz, W, IJV); + + Am.setFromTriplets(IJV.begin(), IJV.end()); + Am.makeCompressed(); +} + +void buildRhs(const Eigen::VectorXd &sqrt_M, + const Eigen::MatrixXd &W, + const Eigen::MatrixXd &Ri, + Eigen::VectorXd &f_rhs) +{ + const int dim = (W.cols() == 4) ? 2 : 3; + const int f_n = W.rows(); + f_rhs.resize(dim * dim * f_n); + + for (int i = 0; i < f_n; i++) + { + auto sqrt_area = sqrt_M(i); + f_rhs(i + 0 * f_n) = sqrt_area * (W(i, 0) * Ri(i, 0) + W(i, 1) * Ri(i, 1)); + f_rhs(i + 1 * f_n) = sqrt_area * (W(i, 0) * Ri(i, 2) + W(i, 1) * Ri(i, 3)); + f_rhs(i + 2 * f_n) = sqrt_area * (W(i, 2) * Ri(i, 0) + W(i, 3) * Ri(i, 1)); + f_rhs(i + 3 * f_n) = sqrt_area * (W(i, 2) * Ri(i, 2) + W(i, 3) * Ri(i, 3)); + } +} + +void get_complement(const Eigen::VectorXi &bnd_ids, int v_n, Eigen::ArrayXi &unknown_ids) +{ // get the complement of bnd_ids. + int assign = 0, i = 0; + for (int get = 0; i < v_n && get < bnd_ids.size(); i++) + { + if (bnd_ids(get) == i) + get++; + else + unknown_ids(assign++) = i; + } + while (i < v_n) + unknown_ids(assign++) = i++; + assert(assign + bnd_ids.size() == v_n); +} + +void build_surface_linear_system(const SCAFData &s, Eigen::SparseMatrix &L, Eigen::VectorXd &rhs) +{ + using namespace Eigen; + using namespace std; + + const int v_n = s.v_num - (s.frame_ids.size()); + const int dim = s.dim; + const int f_n = s.mf_num; + + // to get the complete A + Eigen::VectorXd sqrtM = s.m_M.array().sqrt(); + Eigen::SparseMatrix A(dim * dim * f_n, dim * v_n); + auto decoy_Dx_m = s.Dx_m; + decoy_Dx_m.conservativeResize(s.W_m.rows(), v_n); + auto decoy_Dy_m = s.Dy_m; + decoy_Dy_m.conservativeResize(s.W_m.rows(), v_n); + buildAm(sqrtM, decoy_Dx_m, decoy_Dy_m, s.W_m, A); + + const VectorXi &bnd_ids = s.fixed_ids; + auto bnd_n = bnd_ids.size(); + if (bnd_n == 0) + { + + Eigen::SparseMatrix At = A.transpose(); + At.makeCompressed(); + + Eigen::SparseMatrix id_m(At.rows(), At.rows()); + id_m.setIdentity(); + + L = At * A; + + Eigen::VectorXd frhs; + buildRhs(sqrtM, s.W_m, s.Ri_m, frhs); + rhs = At * frhs; + } + else + { + MatrixXd bnd_pos; + igl::slice(s.w_uv, bnd_ids, 1, bnd_pos); + ArrayXi known_ids(bnd_ids.size() * dim); + ArrayXi unknown_ids((v_n - bnd_ids.rows()) * dim); + get_complement(bnd_ids, v_n, unknown_ids); + VectorXd known_pos(bnd_ids.size() * dim); + for (int d = 0; d < dim; d++) + { + auto n_b = bnd_ids.rows(); + known_ids.segment(d * n_b, n_b) = bnd_ids.array() + d * v_n; + known_pos.segment(d * n_b, n_b) = bnd_pos.col(d); + unknown_ids.block(d * (v_n - n_b), 0, v_n - n_b, unknown_ids.cols()) = + unknown_ids.topRows(v_n - n_b) + d * v_n; + } + + Eigen::SparseMatrix Au, Ae; + igl::slice(A, unknown_ids, 2, Au); + igl::slice(A, known_ids, 2, Ae); + + Eigen::SparseMatrix Aut = Au.transpose(); + Aut.makeCompressed(); + + L = Aut * Au; + + Eigen::VectorXd frhs; + buildRhs(sqrtM, s.W_m, s.Ri_m, frhs); + + rhs = Aut * (frhs - Ae * known_pos); + } + + // add soft constraints. + for (auto const &x : s.soft_cons) + { + int v_idx = x.first; + + for (int d = 0; d < dim; d++) + { + rhs(d * (v_n) + v_idx) += s.soft_const_p * x.second(d); // rhs + L.coeffRef(d * v_n + v_idx, + d * v_n + v_idx) += s.soft_const_p; // diagonal + } + } +} + +void build_scaffold_linear_system(const SCAFData &s, Eigen::SparseMatrix &L, Eigen::VectorXd &rhs) +{ + using namespace Eigen; + + const int f_n = s.W_s.rows(); + const int v_n = s.Dx_s.cols(); + const int dim = s.dim; + + Eigen::VectorXd sqrtM = s.s_M.array().sqrt(); + Eigen::SparseMatrix A(dim * dim * f_n, dim * v_n); + buildAm(sqrtM, s.Dx_s, s.Dy_s, s.W_s, A); + + VectorXi bnd_ids; + igl::cat(1, s.fixed_ids, s.frame_ids, bnd_ids); + + auto bnd_n = bnd_ids.size(); + assert(bnd_n > 0); + MatrixXd bnd_pos; + igl::slice(s.w_uv, bnd_ids, 1, bnd_pos); + + ArrayXi known_ids(bnd_ids.size() * dim); + ArrayXi unknown_ids((v_n - bnd_ids.rows()) * dim); + + get_complement(bnd_ids, v_n, unknown_ids); + + VectorXd known_pos(bnd_ids.size() * dim); + for (int d = 0; d < dim; d++) + { + auto n_b = bnd_ids.rows(); + known_ids.segment(d * n_b, n_b) = bnd_ids.array() + d * v_n; + known_pos.segment(d * n_b, n_b) = bnd_pos.col(d); + unknown_ids.block(d * (v_n - n_b), 0, v_n - n_b, unknown_ids.cols()) = + unknown_ids.topRows(v_n - n_b) + d * v_n; + } + Eigen::VectorXd sqrt_M = s.s_M.array().sqrt(); + + // manual slicing for A(:, unknown/known)' + Eigen::SparseMatrix Au, Ae; + igl::slice(A, unknown_ids, 2, Au); + igl::slice(A, known_ids, 2, Ae); + + Eigen::SparseMatrix Aut = Au.transpose(); + Aut.makeCompressed(); + + L = Aut * Au; + + Eigen::VectorXd frhs; + buildRhs(sqrtM, s.W_s, s.Ri_s, frhs); + + rhs = Aut * (frhs - Ae * known_pos); +} + +void solve_weighted_arap(SCAFData &s, Eigen::MatrixXd &uv) +{ + using namespace Eigen; + using namespace std; + int dim = s.dim; + igl::Timer timer; + timer.start(); + + VectorXi bnd_ids; + igl::cat(1, s.fixed_ids, s.frame_ids, bnd_ids); + const auto v_n = s.v_num; + const auto bnd_n = bnd_ids.size(); + assert(bnd_n > 0); + MatrixXd bnd_pos; + igl::slice(s.w_uv, bnd_ids, 1, bnd_pos); + + ArrayXi known_ids(bnd_n * dim); + ArrayXi unknown_ids((v_n - bnd_n) * dim); + + get_complement(bnd_ids, v_n, unknown_ids); + + VectorXd known_pos(bnd_ids.size() * dim); + for (int d = 0; d < dim; d++) + { + auto n_b = bnd_ids.rows(); + known_ids.segment(d * n_b, n_b) = bnd_ids.array() + d * v_n; + known_pos.segment(d * n_b, n_b) = bnd_pos.col(d); + unknown_ids.block(d * (v_n - n_b), 0, v_n - n_b, unknown_ids.cols()) = + unknown_ids.topRows(v_n - n_b) + d * v_n; + } + + Eigen::SparseMatrix L; + Eigen::VectorXd rhs; + + // fixed frame solving: + // x_e as the fixed frame, x_u for unknowns (mesh + unknown scaffold) + // min ||(A_u*x_u + A_e*x_e) - b||^2 + // => A_u'*A_u*x_u = Au'* (b - A_e*x_e) := Au'* b_u + // + // separate matrix build: + // min ||A_m x_m - b_m||^2 + ||A_s x_all - b_s||^2 + soft + proximal + // First change dimension of A_m to fit for x_all + // (Not just at the end, since x_all is flattened along dimensions) + // L = A_m'*A_m + A_s'*A_s + soft + proximal + // rhs = A_m'* b_m + A_s' * b_s + soft + proximal + // + Eigen::SparseMatrix L_m, L_s; + Eigen::VectorXd rhs_m, rhs_s; + build_surface_linear_system(s, L_m, rhs_m); // complete Am, with soft + build_scaffold_linear_system(s, L_s, rhs_s); // complete As, without proximal + + L = L_m + L_s; + rhs = rhs_m + rhs_s; + L.makeCompressed(); + + Eigen::VectorXd unknown_Uc((v_n - s.frame_ids.size() - s.fixed_ids.size()) * dim), Uc(dim * v_n); + + SimplicialLDLT> solver; + unknown_Uc = solver.compute(L).solve(rhs); + igl::slice_into(unknown_Uc, unknown_ids.matrix(), 1, Uc); + igl::slice_into(known_pos, known_ids.matrix(), 1, Uc); + + uv = Map>(Uc.data(), v_n, dim); +} + +double perform_iteration(SCAFData &s) +{ + Eigen::MatrixXd V_out = s.w_uv; + compute_jacobians(s, V_out, true); + igl::slim_update_weights_and_closest_rotations_with_jacobians(s.Ji_m, s.slim_energy, 0, s.W_m, s.Ri_m); + igl::slim_update_weights_and_closest_rotations_with_jacobians(s.Ji_s, s.scaf_energy, 0, s.W_s, s.Ri_s); + solve_weighted_arap(s, V_out); + auto whole_E = [&s](Eigen::MatrixXd &uv) { return compute_energy(s, uv, true); }; + + Eigen::MatrixXi w_T; + if (s.m_T.cols() == s.s_T.cols()) + igl::cat(1, s.m_T, s.s_T, w_T); + else + w_T = s.s_T; + return igl::flip_avoiding_line_search(w_T, s.w_uv, V_out, + whole_E, -1) / + s.mesh_measure; +} +} +} + +IGL_INLINE void igl::scaf_precompute( + const Eigen::MatrixXd &V, + const Eigen::MatrixXi &F, + const Eigen::MatrixXd &V_init, + igl::SCAFData &data, + igl::MappingEnergyType slim_energy, + Eigen::VectorXi &b, + Eigen::MatrixXd &bc, + double soft_p) +{ + Eigen::MatrixXd CN; + Eigen::MatrixXi FN; + igl::scaf::add_new_patch(data, V, F, Eigen::RowVector2d(0, 0), V_init); + data.soft_const_p = soft_p; + for (int i = 0; i < b.rows(); i++) + data.soft_cons[b(i)] = bc.row(i); + data.slim_energy = slim_energy; + + auto &s = data; + + if (!data.has_pre_calc) + { + int v_n = s.mv_num + s.sv_num; + int f_n = s.mf_num + s.sf_num; + int dim = s.dim; + Eigen::MatrixXd F1, F2, F3; + igl::local_basis(s.m_V, s.m_T, F1, F2, F3); + auto face_proj = [](Eigen::MatrixXd& F){ + std::vector >IJV; + int f_num = F.rows(); + for(int i=0; i(i, i, F(i,0))); + IJV.push_back(Eigen::Triplet(i, i+f_num, F(i,1))); + IJV.push_back(Eigen::Triplet(i, i+2*f_num, F(i,2))); + } + Eigen::SparseMatrix P(f_num, 3*f_num); + P.setFromTriplets(IJV.begin(), IJV.end()); + return P; + }; + Eigen::SparseMatrix G; + igl::grad(s.m_V, s.m_T, G); + s.Dx_m = face_proj(F1) * G; + s.Dy_m = face_proj(F2) * G; + + igl::scaf::compute_scaffold_gradient_matrix(s, s.Dx_s, s.Dy_s); + + s.Dx_m.makeCompressed(); + s.Dy_m.makeCompressed(); + s.Ri_m = Eigen::MatrixXd::Zero(s.Dx_m.rows(), dim * dim); + s.Ji_m.resize(s.Dx_m.rows(), dim * dim); + s.W_m.resize(s.Dx_m.rows(), dim * dim); + + s.Dx_s.makeCompressed(); + s.Dy_s.makeCompressed(); + s.Ri_s = Eigen::MatrixXd::Zero(s.Dx_s.rows(), dim * dim); + s.Ji_s.resize(s.Dx_s.rows(), dim * dim); + s.W_s.resize(s.Dx_s.rows(), dim * dim); + + data.has_pre_calc = true; + } +} + +IGL_INLINE Eigen::MatrixXd igl::scaf_solve(SCAFData &s, int iter_num) +{ + using namespace std; + using namespace Eigen; + s.energy = igl::scaf::compute_energy(s, s.w_uv, false) / s.mesh_measure; + + for (int it = 0; it < iter_num; it++) + { + s.total_energy = igl::scaf::compute_energy(s, s.w_uv, true) / s.mesh_measure; + s.rect_frame_V = Eigen::MatrixXd(); + igl::scaf::mesh_improve(s); + + double new_weight = s.mesh_measure * s.energy / (s.sf_num * 100); + s.scaffold_factor = new_weight; + igl::scaf::update_scaffold(s); + + s.total_energy = igl::scaf::perform_iteration(s); + + s.energy = + igl::scaf::compute_energy(s, s.w_uv, false) / s.mesh_measure; + } + + return s.w_uv.topRows(s.mv_num); +} + +#ifdef IGL_STATIC_LIBRARY +#endif \ No newline at end of file diff --git a/include/igl/scaf.h b/include/igl/scaf.h new file mode 100644 index 000000000..425cc2bb3 --- /dev/null +++ b/include/igl/scaf.h @@ -0,0 +1,101 @@ +// This file is part of libigl, a simple c++ geometry processing library. +// +// Copyright (C) 2018 Zhongshi Jiang +// +// 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/. +#ifndef IGL_SCAF_H +#define IGL_SCAF_H + +#include "slim.h" +#include "igl_inline.h" +#include "MappingEnergyType.h" + +namespace igl +{ + // Use a similar interface to igl::slim + // Implement ready-to-use 2D version of the algorithm described in + // SCAF: Simplicial Complex Augmentation Framework for Bijective Maps + // Zhongshi Jiang, Scott Schaefer, Daniele Panozzo, ACM Trancaction on Graphics (Proc. SIGGRAPH Asia 2017) + // For a complete implementation and customized UI, please refer to https://github.com/jiangzhongshi/scaffold-map + + struct SCAFData + { + double scaffold_factor = 10; + igl::MappingEnergyType scaf_energy = igl::MappingEnergyType::SYMMETRIC_DIRICHLET; + igl::MappingEnergyType slim_energy = igl::MappingEnergyType::SYMMETRIC_DIRICHLET; + + // Output + int dim = 2; + double total_energy; // scaffold + isometric + double energy; // objective value + + long mv_num = 0, mf_num = 0; + long sv_num = 0, sf_num = 0; + long v_num{}, f_num = 0; + Eigen::MatrixXd m_V; // input initial mesh V + Eigen::MatrixXi m_T; // input initial mesh F/T + // INTERNAL + Eigen::MatrixXd w_uv; // whole domain uv: mesh + free vertices + Eigen::MatrixXi s_T; // scaffold domain tets: scaffold tets + Eigen::MatrixXi w_T; + + Eigen::VectorXd m_M; // mesh area or volume + Eigen::VectorXd s_M; // scaffold area or volume + Eigen::VectorXd w_M; // area/volume weights for whole + double mesh_measure; // area or volume + double proximal_p = 0; + + Eigen::VectorXi frame_ids; + Eigen::VectorXi fixed_ids; + + std::map soft_cons; + double soft_const_p = 1e4; + + Eigen::VectorXi internal_bnd; + Eigen::MatrixXd rect_frame_V; + // multi-chart support + std::vector component_sizes; + std::vector bnd_sizes; + + // reweightedARAP interior variables. + bool has_pre_calc = false; + Eigen::SparseMatrix Dx_s, Dy_s, Dz_s; + Eigen::SparseMatrix Dx_m, Dy_m, Dz_m; + Eigen::MatrixXd Ri_m, Ji_m, Ri_s, Ji_s; + Eigen::MatrixXd W_m, W_s; + }; + + +// Compute necessary information to start using SCAF +// Inputs: +// V #V by 3 list of mesh vertex positions +// F #F by 3/3 list of mesh faces (triangles/tets) +// data igl::SCAFData +// slim_energy Energy type to minimize +// b list of boundary indices into V (soft constraint) +// bc #b by dim list of boundary conditions (soft constraint) +// soft_p Soft penalty factor (can be zero) + IGL_INLINE void scaf_precompute( + const Eigen::MatrixXd &V, + const Eigen::MatrixXi &F, + const Eigen::MatrixXd &V_init, + SCAFData &data, + MappingEnergyType slim_energy, + Eigen::VectorXi& b, + Eigen::MatrixXd& bc, + double soft_p); + + +// Run iter_num iterations of SCAF, with precomputed data +// Outputs: +// V_o (in SLIMData): #V by dim list of mesh vertex positions + IGL_INLINE Eigen::MatrixXd scaf_solve(SCAFData &data, int iter_num); + } + +#ifndef IGL_STATIC_LIBRARY +# include "scaf.cpp" +#endif + +#endif //IGL_SCAF_H diff --git a/include/igl/slice.cpp b/include/igl/slice.cpp index df3a8e9ce..0d20df96e 100644 --- a/include/igl/slice.cpp +++ b/include/igl/slice.cpp @@ -355,6 +355,7 @@ template Eigen::Matrix igl::slice >, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase > const&, Eigen::DenseBase > const&, int, Eigen::Matrix&); template void igl::slice, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::DenseBase > const&, Eigen::DenseBase > const&, Eigen::DenseBase > const&, Eigen::PlainObjectBase >&); template void igl::slice(Eigen::SparseMatrix const&, Eigen::Matrix const&, Eigen::Matrix const&, Eigen::SparseMatrix&); +template void igl::slice, Eigen::Array, Eigen::SparseMatrix >(Eigen::SparseMatrix const&, Eigen::DenseBase > const&, int, Eigen::SparseMatrix&); #ifdef WIN32 template void igl::slice, class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, class Eigen::PlainObjectBase>>(class Eigen::Matrix<__int64, -1, 1, 0, -1, 1> const &, class Eigen::DenseBase> const &, int, class Eigen::PlainObjectBase> &); template void igl::slice>, class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, class Eigen::PlainObjectBase>>(class Eigen::PlainObjectBase> const &, class Eigen::DenseBase> const &, int, class Eigen::PlainObjectBase> &); diff --git a/include/igl/slim.cpp b/include/igl/slim.cpp index a88674274..b1d0f4dba 100644 --- a/include/igl/slim.cpp +++ b/include/igl/slim.cpp @@ -24,6 +24,7 @@ #include "volume.h" #include "polar_svd.h" #include "flip_avoiding_line_search.h" +#include "mapping_energy_with_jacobians.h" #include #include @@ -47,10 +48,6 @@ namespace igl namespace slim { // Definitions of internal functions - IGL_INLINE void compute_surface_gradient_matrix(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, - const Eigen::MatrixXd &F1, const Eigen::MatrixXd &F2, - Eigen::SparseMatrix &D1, Eigen::SparseMatrix &D2); - IGL_INLINE void buildA(igl::SLIMData& s, std::vector > & IJV); IGL_INLINE void buildRhs(igl::SLIMData& s, const Eigen::SparseMatrix &A); IGL_INLINE void add_soft_constraints(igl::SLIMData& s, Eigen::SparseMatrix &L); IGL_INLINE double compute_energy(igl::SLIMData& s, Eigen::MatrixXd &V_new); @@ -58,10 +55,7 @@ namespace igl const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, Eigen::MatrixXd &V_o); - IGL_INLINE double compute_energy_with_jacobians(igl::SLIMData& s, - const Eigen::MatrixXd &V, - const Eigen::MatrixXi &F, const Eigen::MatrixXd &Ji, - Eigen::MatrixXd &uv, Eigen::VectorXd &areas); + IGL_INLINE void solve_weighted_arap(igl::SLIMData& s, const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, @@ -69,28 +63,12 @@ namespace igl Eigen::VectorXi &soft_b_p, Eigen::MatrixXd &soft_bc_p); IGL_INLINE void update_weights_and_closest_rotations( igl::SLIMData& s, - const Eigen::MatrixXd &V, - const Eigen::MatrixXi &F, Eigen::MatrixXd &uv); IGL_INLINE void compute_jacobians(igl::SLIMData& s, const Eigen::MatrixXd &uv); IGL_INLINE void build_linear_system(igl::SLIMData& s, Eigen::SparseMatrix &L); IGL_INLINE void pre_calc(igl::SLIMData& s); // Implementation - IGL_INLINE void compute_surface_gradient_matrix(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, - const Eigen::MatrixXd &F1, const Eigen::MatrixXd &F2, - Eigen::SparseMatrix &D1, Eigen::SparseMatrix &D2) - { - - Eigen::SparseMatrix G; - igl::grad(V, F, G); - Eigen::SparseMatrix Dx = G.block(0, 0, F.rows(), V.rows()); - Eigen::SparseMatrix Dy = G.block(F.rows(), 0, F.rows(), V.rows()); - Eigen::SparseMatrix Dz = G.block(2 * F.rows(), 0, F.rows(), V.rows()); - - D1 = F1.col(0).asDiagonal() * Dx + F1.col(1).asDiagonal() * Dy + F1.col(2).asDiagonal() * Dz; - D2 = F2.col(0).asDiagonal() * Dx + F2.col(1).asDiagonal() * Dy + F2.col(2).asDiagonal() * Dz; - } IGL_INLINE void compute_jacobians(igl::SLIMData& s, const Eigen::MatrixXd &uv) { @@ -116,280 +94,14 @@ namespace igl } } - IGL_INLINE void update_weights_and_closest_rotations(igl::SLIMData& s, - const Eigen::MatrixXd &V, - const Eigen::MatrixXi &F, - Eigen::MatrixXd &uv) + IGL_INLINE void update_weights_and_closest_rotations(igl::SLIMData& s, Eigen::MatrixXd &uv) { compute_jacobians(s, uv); - - const double eps = 1e-8; - double exp_f = s.exp_factor; - - if (s.dim == 2) - { - for (int i = 0; i < s.Ji.rows(); ++i) - { - typedef Eigen::Matrix Mat2; - typedef Eigen::Matrix Vec2; - Mat2 ji, ri, ti, ui, vi; - Vec2 sing; - Vec2 closest_sing_vec; - Mat2 mat_W; - Vec2 m_sing_new; - double s1, s2; - - ji(0, 0) = s.Ji(i, 0); - ji(0, 1) = s.Ji(i, 1); - ji(1, 0) = s.Ji(i, 2); - ji(1, 1) = s.Ji(i, 3); - - igl::polar_svd(ji, ri, ti, ui, sing, vi); - - s1 = sing(0); - s2 = sing(1); - - // Update Weights according to energy - switch (s.slim_energy) - { - case igl::SLIMData::ARAP: - { - m_sing_new << 1, 1; - break; - } - case igl::SLIMData::SYMMETRIC_DIRICHLET: - { - double s1_g = 2 * (s1 - pow(s1, -3)); - double s2_g = 2 * (s2 - pow(s2, -3)); - m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))); - break; - } - case igl::SLIMData::LOG_ARAP: - { - double s1_g = 2 * (log(s1) / s1); - double s2_g = 2 * (log(s2) / s2); - m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))); - break; - } - case igl::SLIMData::CONFORMAL: - { - double s1_g = 1 / (2 * s2) - s2 / (2 * pow(s1, 2)); - double s2_g = 1 / (2 * s1) - s1 / (2 * pow(s2, 2)); - - double geo_avg = sqrt(s1 * s2); - double s1_min = geo_avg; - double s2_min = geo_avg; - - m_sing_new << sqrt(s1_g / (2 * (s1 - s1_min))), sqrt(s2_g / (2 * (s2 - s2_min))); - - // change local step - closest_sing_vec << s1_min, s2_min; - ri = ui * closest_sing_vec.asDiagonal() * vi.transpose(); - break; - } - case igl::SLIMData::EXP_CONFORMAL: - { - double s1_g = 2 * (s1 - pow(s1, -3)); - double s2_g = 2 * (s2 - pow(s2, -3)); - - double geo_avg = sqrt(s1 * s2); - double s1_min = geo_avg; - double s2_min = geo_avg; - - double in_exp = exp_f * ((pow(s1, 2) + pow(s2, 2)) / (2 * s1 * s2)); - double exp_thing = exp(in_exp); - - s1_g *= exp_thing * exp_f; - s2_g *= exp_thing * exp_f; - - m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))); - break; - } - case igl::SLIMData::EXP_SYMMETRIC_DIRICHLET: - { - double s1_g = 2 * (s1 - pow(s1, -3)); - double s2_g = 2 * (s2 - pow(s2, -3)); - - double in_exp = exp_f * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2)); - double exp_thing = exp(in_exp); - - s1_g *= exp_thing * exp_f; - s2_g *= exp_thing * exp_f; - - m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))); - break; - } - } - - if (std::abs(s1 - 1) < eps) m_sing_new(0) = 1; - if (std::abs(s2 - 1) < eps) m_sing_new(1) = 1; - mat_W = ui * m_sing_new.asDiagonal() * ui.transpose(); - - s.W_11(i) = mat_W(0, 0); - s.W_12(i) = mat_W(0, 1); - s.W_21(i) = mat_W(1, 0); - s.W_22(i) = mat_W(1, 1); - - // 2) Update local step (doesn't have to be a rotation, for instance in case of conformal energy) - s.Ri(i, 0) = ri(0, 0); - s.Ri(i, 1) = ri(1, 0); - s.Ri(i, 2) = ri(0, 1); - s.Ri(i, 3) = ri(1, 1); - } - } - else - { - typedef Eigen::Matrix Vec3; - typedef Eigen::Matrix Mat3; - Mat3 ji; - Vec3 m_sing_new; - Vec3 closest_sing_vec; - const double sqrt_2 = sqrt(2); - for (int i = 0; i < s.Ji.rows(); ++i) - { - ji(0, 0) = s.Ji(i, 0); - ji(0, 1) = s.Ji(i, 1); - ji(0, 2) = s.Ji(i, 2); - ji(1, 0) = s.Ji(i, 3); - ji(1, 1) = s.Ji(i, 4); - ji(1, 2) = s.Ji(i, 5); - ji(2, 0) = s.Ji(i, 6); - ji(2, 1) = s.Ji(i, 7); - ji(2, 2) = s.Ji(i, 8); - - Mat3 ri, ti, ui, vi; - Vec3 sing; - igl::polar_svd(ji, ri, ti, ui, sing, vi); - - double s1 = sing(0); - double s2 = sing(1); - double s3 = sing(2); - - // 1) Update Weights - switch (s.slim_energy) - { - case igl::SLIMData::ARAP: - { - m_sing_new << 1, 1, 1; - break; - } - case igl::SLIMData::LOG_ARAP: - { - double s1_g = 2 * (log(s1) / s1); - double s2_g = 2 * (log(s2) / s2); - double s3_g = 2 * (log(s3) / s3); - m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))), sqrt(s3_g / (2 * (s3 - 1))); - break; - } - case igl::SLIMData::SYMMETRIC_DIRICHLET: - { - double s1_g = 2 * (s1 - pow(s1, -3)); - double s2_g = 2 * (s2 - pow(s2, -3)); - double s3_g = 2 * (s3 - pow(s3, -3)); - m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))), sqrt(s3_g / (2 * (s3 - 1))); - break; - } - case igl::SLIMData::EXP_SYMMETRIC_DIRICHLET: - { - double s1_g = 2 * (s1 - pow(s1, -3)); - double s2_g = 2 * (s2 - pow(s2, -3)); - double s3_g = 2 * (s3 - pow(s3, -3)); - m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))), sqrt(s3_g / (2 * (s3 - 1))); - - double in_exp = exp_f * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2) + pow(s3, 2) + pow(s3, -2)); - double exp_thing = exp(in_exp); - - s1_g *= exp_thing * exp_f; - s2_g *= exp_thing * exp_f; - s3_g *= exp_thing * exp_f; - - m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))), sqrt(s3_g / (2 * (s3 - 1))); - - break; - } - case igl::SLIMData::CONFORMAL: - { - double common_div = 9 * (pow(s1 * s2 * s3, 5. / 3.)); - - double s1_g = (-2 * s2 * s3 * (pow(s2, 2) + pow(s3, 2) - 2 * pow(s1, 2))) / common_div; - double s2_g = (-2 * s1 * s3 * (pow(s1, 2) + pow(s3, 2) - 2 * pow(s2, 2))) / common_div; - double s3_g = (-2 * s1 * s2 * (pow(s1, 2) + pow(s2, 2) - 2 * pow(s3, 2))) / common_div; - - double closest_s = sqrt(pow(s1, 2) + pow(s3, 2)) / sqrt_2; - double s1_min = closest_s; - double s2_min = closest_s; - double s3_min = closest_s; - - m_sing_new << sqrt(s1_g / (2 * (s1 - s1_min))), sqrt(s2_g / (2 * (s2 - s2_min))), sqrt( - s3_g / (2 * (s3 - s3_min))); - - // change local step - closest_sing_vec << s1_min, s2_min, s3_min; - ri = ui * closest_sing_vec.asDiagonal() * vi.transpose(); - break; - } - case igl::SLIMData::EXP_CONFORMAL: - { - // E_conf = (s1^2 + s2^2 + s3^2)/(3*(s1*s2*s3)^(2/3) ) - // dE_conf/ds1 = (-2*(s2*s3)*(s2^2+s3^2 -2*s1^2) ) / (9*(s1*s2*s3)^(5/3)) - // Argmin E_conf(s1): s1 = sqrt(s1^2+s2^2)/sqrt(2) - double common_div = 9 * (pow(s1 * s2 * s3, 5. / 3.)); - - double s1_g = (-2 * s2 * s3 * (pow(s2, 2) + pow(s3, 2) - 2 * pow(s1, 2))) / common_div; - double s2_g = (-2 * s1 * s3 * (pow(s1, 2) + pow(s3, 2) - 2 * pow(s2, 2))) / common_div; - double s3_g = (-2 * s1 * s2 * (pow(s1, 2) + pow(s2, 2) - 2 * pow(s3, 2))) / common_div; - - double in_exp = exp_f * ((pow(s1, 2) + pow(s2, 2) + pow(s3, 2)) / (3 * pow((s1 * s2 * s3), 2. / 3)));; - double exp_thing = exp(in_exp); - - double closest_s = sqrt(pow(s1, 2) + pow(s3, 2)) / sqrt_2; - double s1_min = closest_s; - double s2_min = closest_s; - double s3_min = closest_s; - - s1_g *= exp_thing * exp_f; - s2_g *= exp_thing * exp_f; - s3_g *= exp_thing * exp_f; - - m_sing_new << sqrt(s1_g / (2 * (s1 - s1_min))), sqrt(s2_g / (2 * (s2 - s2_min))), sqrt( - s3_g / (2 * (s3 - s3_min))); - - // change local step - closest_sing_vec << s1_min, s2_min, s3_min; - ri = ui * closest_sing_vec.asDiagonal() * vi.transpose(); - } - } - if (std::abs(s1 - 1) < eps) m_sing_new(0) = 1; - if (std::abs(s2 - 1) < eps) m_sing_new(1) = 1; - if (std::abs(s3 - 1) < eps) m_sing_new(2) = 1; - Mat3 mat_W; - mat_W = ui * m_sing_new.asDiagonal() * ui.transpose(); - - s.W_11(i) = mat_W(0, 0); - s.W_12(i) = mat_W(0, 1); - s.W_13(i) = mat_W(0, 2); - s.W_21(i) = mat_W(1, 0); - s.W_22(i) = mat_W(1, 1); - s.W_23(i) = mat_W(1, 2); - s.W_31(i) = mat_W(2, 0); - s.W_32(i) = mat_W(2, 1); - s.W_33(i) = mat_W(2, 2); - - // 2) Update closest rotations (not rotations in case of conformal energy) - s.Ri(i, 0) = ri(0, 0); - s.Ri(i, 1) = ri(1, 0); - s.Ri(i, 2) = ri(2, 0); - s.Ri(i, 3) = ri(0, 1); - s.Ri(i, 4) = ri(1, 1); - s.Ri(i, 5) = ri(2, 1); - s.Ri(i, 6) = ri(0, 2); - s.Ri(i, 7) = ri(1, 2); - s.Ri(i, 8) = ri(2, 2); - } // for loop end - - } // if dim end - + slim_update_weights_and_closest_rotations_with_jacobians(s.Ji, s.slim_energy, s.exp_factor, s.W, s.Ri); } + + + IGL_INLINE void solve_weighted_arap(igl::SLIMData& s, const Eigen::MatrixXd &V, @@ -448,12 +160,25 @@ namespace igl s.dim = 2; Eigen::MatrixXd F1, F2, F3; igl::local_basis(s.V, s.F, F1, F2, F3); - compute_surface_gradient_matrix(s.V, s.F, F1, F2, s.Dx, s.Dy); + Eigen::SparseMatrix G; + igl::grad(s.V, s.F, G); + Eigen::SparseMatrix Face_Proj; - s.W_11.resize(s.f_n); - s.W_12.resize(s.f_n); - s.W_21.resize(s.f_n); - s.W_22.resize(s.f_n); + auto face_proj = [](Eigen::MatrixXd& F){ + std::vector >IJV; + int f_num = F.rows(); + for(int i=0; i(i, i, F(i,0))); + IJV.push_back(Eigen::Triplet(i, i+f_num, F(i,1))); + IJV.push_back(Eigen::Triplet(i, i+2*f_num, F(i,2))); + } + Eigen::SparseMatrix P(f_num, 3*f_num); + P.setFromTriplets(IJV.begin(), IJV.end()); + return P; + }; + + s.Dx = face_proj(F1) * G; + s.Dy = face_proj(F2) * G; } else { @@ -464,19 +189,9 @@ namespace igl s.Dx = G.block(0, 0, s.F.rows(), s.V.rows()); s.Dy = G.block(s.F.rows(), 0, s.F.rows(), s.V.rows()); s.Dz = G.block(2 * s.F.rows(), 0, s.F.rows(), s.V.rows()); - - - s.W_11.resize(s.f_n); - s.W_12.resize(s.f_n); - s.W_13.resize(s.f_n); - s.W_21.resize(s.f_n); - s.W_22.resize(s.f_n); - s.W_23.resize(s.f_n); - s.W_31.resize(s.f_n); - s.W_32.resize(s.f_n); - s.W_33.resize(s.f_n); } + s.W.resize(s.f_n, s.dim * s.dim); s.Dx.makeCompressed(); s.Dy.makeCompressed(); s.Dz.makeCompressed(); @@ -501,7 +216,7 @@ namespace igl std::vector > IJV; #ifdef SLIM_CACHED - buildA(s,IJV); + slim_buildA(s.Dx, s.Dy, s.Dz, s.W, IJV); if (s.A.rows() == 0) { s.A = Eigen::SparseMatrix(s.dim * s.dim * s.f_n, s.dim * s.v_n); @@ -511,7 +226,7 @@ namespace igl igl::sparse_cached(IJV,s.A_data,s.A); #else Eigen::SparseMatrix A(s.dim * s.dim * s.f_n, s.dim * s.v_n); - buildA(s,IJV); + slim_buildA(s.Dx, s.Dy, s.Dz, s.W, IJV); A.setFromTriplets(IJV.begin(),IJV.end()); A.makeCompressed(); #endif @@ -574,7 +289,7 @@ namespace igl IGL_INLINE double compute_energy(igl::SLIMData& s, Eigen::MatrixXd &V_new) { compute_jacobians(s,V_new); - return compute_energy_with_jacobians(s, s.V, s.F, s.Ji, V_new, s.M) + + return mapping_energy_with_jacobians(s.Ji, s.M, s.slim_energy, s.exp_factor) + compute_soft_const_energy(s, s.V, s.F, V_new); } @@ -591,255 +306,7 @@ namespace igl return e; } - IGL_INLINE double compute_energy_with_jacobians(igl::SLIMData& s, - const Eigen::MatrixXd &V, - const Eigen::MatrixXi &F, const Eigen::MatrixXd &Ji, - Eigen::MatrixXd &uv, Eigen::VectorXd &areas) - { - double energy = 0; - if (s.dim == 2) - { - Eigen::Matrix ji; - for (int i = 0; i < s.f_n; i++) - { - ji(0, 0) = Ji(i, 0); - ji(0, 1) = Ji(i, 1); - ji(1, 0) = Ji(i, 2); - ji(1, 1) = Ji(i, 3); - - typedef Eigen::Matrix Mat2; - typedef Eigen::Matrix Vec2; - Mat2 ri, ti, ui, vi; - Vec2 sing; - igl::polar_svd(ji, ri, ti, ui, sing, vi); - double s1 = sing(0); - double s2 = sing(1); - - switch (s.slim_energy) - { - case igl::SLIMData::ARAP: - { - energy += areas(i) * (pow(s1 - 1, 2) + pow(s2 - 1, 2)); - break; - } - case igl::SLIMData::SYMMETRIC_DIRICHLET: - { - energy += areas(i) * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2)); - break; - } - case igl::SLIMData::EXP_SYMMETRIC_DIRICHLET: - { - energy += areas(i) * exp(s.exp_factor * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2))); - break; - } - case igl::SLIMData::LOG_ARAP: - { - energy += areas(i) * (pow(log(s1), 2) + pow(log(s2), 2)); - break; - } - case igl::SLIMData::CONFORMAL: - { - energy += areas(i) * ((pow(s1, 2) + pow(s2, 2)) / (2 * s1 * s2)); - break; - } - case igl::SLIMData::EXP_CONFORMAL: - { - energy += areas(i) * exp(s.exp_factor * ((pow(s1, 2) + pow(s2, 2)) / (2 * s1 * s2))); - break; - } - - } - - } - } - else - { - Eigen::Matrix ji; - for (int i = 0; i < s.f_n; i++) - { - ji(0, 0) = Ji(i, 0); - ji(0, 1) = Ji(i, 1); - ji(0, 2) = Ji(i, 2); - ji(1, 0) = Ji(i, 3); - ji(1, 1) = Ji(i, 4); - ji(1, 2) = Ji(i, 5); - ji(2, 0) = Ji(i, 6); - ji(2, 1) = Ji(i, 7); - ji(2, 2) = Ji(i, 8); - - typedef Eigen::Matrix Mat3; - typedef Eigen::Matrix Vec3; - Mat3 ri, ti, ui, vi; - Vec3 sing; - igl::polar_svd(ji, ri, ti, ui, sing, vi); - double s1 = sing(0); - double s2 = sing(1); - double s3 = sing(2); - - switch (s.slim_energy) - { - case igl::SLIMData::ARAP: - { - energy += areas(i) * (pow(s1 - 1, 2) + pow(s2 - 1, 2) + pow(s3 - 1, 2)); - break; - } - case igl::SLIMData::SYMMETRIC_DIRICHLET: - { - energy += areas(i) * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2) + pow(s3, 2) + pow(s3, -2)); - break; - } - case igl::SLIMData::EXP_SYMMETRIC_DIRICHLET: - { - energy += areas(i) * exp(s.exp_factor * - (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2) + pow(s3, 2) + pow(s3, -2))); - break; - } - case igl::SLIMData::LOG_ARAP: - { - energy += areas(i) * (pow(log(s1), 2) + pow(log(std::abs(s2)), 2) + pow(log(std::abs(s3)), 2)); - break; - } - case igl::SLIMData::CONFORMAL: - { - energy += areas(i) * ((pow(s1, 2) + pow(s2, 2) + pow(s3, 2)) / (3 * pow(s1 * s2 * s3, 2. / 3.))); - break; - } - case igl::SLIMData::EXP_CONFORMAL: - { - energy += areas(i) * exp((pow(s1, 2) + pow(s2, 2) + pow(s3, 2)) / (3 * pow(s1 * s2 * s3, 2. / 3.))); - break; - } - } - } - } - - return energy; - } - - IGL_INLINE void buildA(igl::SLIMData& s, std::vector > & IJV) - { - // formula (35) in paper - if (s.dim == 2) - { - IJV.reserve(4 * (s.Dx.outerSize() + s.Dy.outerSize())); - - /*A = [W11*Dx, W12*Dx; - W11*Dy, W12*Dy; - W21*Dx, W22*Dx; - W21*Dy, W22*Dy];*/ - for (int k = 0; k < s.Dx.outerSize(); ++k) - { - for (Eigen::SparseMatrix::InnerIterator it(s.Dx, k); it; ++it) - { - int dx_r = it.row(); - int dx_c = it.col(); - double val = it.value(); - - IJV.push_back(Eigen::Triplet(dx_r, dx_c, val * s.W_11(dx_r))); - IJV.push_back(Eigen::Triplet(dx_r, s.v_n + dx_c, val * s.W_12(dx_r))); - - IJV.push_back(Eigen::Triplet(2 * s.f_n + dx_r, dx_c, val * s.W_21(dx_r))); - IJV.push_back(Eigen::Triplet(2 * s.f_n + dx_r, s.v_n + dx_c, val * s.W_22(dx_r))); - } - } - - for (int k = 0; k < s.Dy.outerSize(); ++k) - { - for (Eigen::SparseMatrix::InnerIterator it(s.Dy, k); it; ++it) - { - int dy_r = it.row(); - int dy_c = it.col(); - double val = it.value(); - - IJV.push_back(Eigen::Triplet(s.f_n + dy_r, dy_c, val * s.W_11(dy_r))); - IJV.push_back(Eigen::Triplet(s.f_n + dy_r, s.v_n + dy_c, val * s.W_12(dy_r))); - - IJV.push_back(Eigen::Triplet(3 * s.f_n + dy_r, dy_c, val * s.W_21(dy_r))); - IJV.push_back(Eigen::Triplet(3 * s.f_n + dy_r, s.v_n + dy_c, val * s.W_22(dy_r))); - } - } - } - else - { - - /*A = [W11*Dx, W12*Dx, W13*Dx; - W11*Dy, W12*Dy, W13*Dy; - W11*Dz, W12*Dz, W13*Dz; - W21*Dx, W22*Dx, W23*Dx; - W21*Dy, W22*Dy, W23*Dy; - W21*Dz, W22*Dz, W23*Dz; - W31*Dx, W32*Dx, W33*Dx; - W31*Dy, W32*Dy, W33*Dy; - W31*Dz, W32*Dz, W33*Dz;];*/ - IJV.reserve(9 * (s.Dx.outerSize() + s.Dy.outerSize() + s.Dz.outerSize())); - for (int k = 0; k < s.Dx.outerSize(); k++) - { - for (Eigen::SparseMatrix::InnerIterator it(s.Dx, k); it; ++it) - { - int dx_r = it.row(); - int dx_c = it.col(); - double val = it.value(); - - IJV.push_back(Eigen::Triplet(dx_r, dx_c, val * s.W_11(dx_r))); - IJV.push_back(Eigen::Triplet(dx_r, s.v_n + dx_c, val * s.W_12(dx_r))); - IJV.push_back(Eigen::Triplet(dx_r, 2 * s.v_n + dx_c, val * s.W_13(dx_r))); - - IJV.push_back(Eigen::Triplet(3 * s.f_n + dx_r, dx_c, val * s.W_21(dx_r))); - IJV.push_back(Eigen::Triplet(3 * s.f_n + dx_r, s.v_n + dx_c, val * s.W_22(dx_r))); - IJV.push_back(Eigen::Triplet(3 * s.f_n + dx_r, 2 * s.v_n + dx_c, val * s.W_23(dx_r))); - - IJV.push_back(Eigen::Triplet(6 * s.f_n + dx_r, dx_c, val * s.W_31(dx_r))); - IJV.push_back(Eigen::Triplet(6 * s.f_n + dx_r, s.v_n + dx_c, val * s.W_32(dx_r))); - IJV.push_back(Eigen::Triplet(6 * s.f_n + dx_r, 2 * s.v_n + dx_c, val * s.W_33(dx_r))); - } - } - - for (int k = 0; k < s.Dy.outerSize(); k++) - { - for (Eigen::SparseMatrix::InnerIterator it(s.Dy, k); it; ++it) - { - int dy_r = it.row(); - int dy_c = it.col(); - double val = it.value(); - - IJV.push_back(Eigen::Triplet(s.f_n + dy_r, dy_c, val * s.W_11(dy_r))); - IJV.push_back(Eigen::Triplet(s.f_n + dy_r, s.v_n + dy_c, val * s.W_12(dy_r))); - IJV.push_back(Eigen::Triplet(s.f_n + dy_r, 2 * s.v_n + dy_c, val * s.W_13(dy_r))); - - IJV.push_back(Eigen::Triplet(4 * s.f_n + dy_r, dy_c, val * s.W_21(dy_r))); - IJV.push_back(Eigen::Triplet(4 * s.f_n + dy_r, s.v_n + dy_c, val * s.W_22(dy_r))); - IJV.push_back(Eigen::Triplet(4 * s.f_n + dy_r, 2 * s.v_n + dy_c, val * s.W_23(dy_r))); - - IJV.push_back(Eigen::Triplet(7 * s.f_n + dy_r, dy_c, val * s.W_31(dy_r))); - IJV.push_back(Eigen::Triplet(7 * s.f_n + dy_r, s.v_n + dy_c, val * s.W_32(dy_r))); - IJV.push_back(Eigen::Triplet(7 * s.f_n + dy_r, 2 * s.v_n + dy_c, val * s.W_33(dy_r))); - } - } - - for (int k = 0; k < s.Dz.outerSize(); k++) - { - for (Eigen::SparseMatrix::InnerIterator it(s.Dz, k); it; ++it) - { - int dz_r = it.row(); - int dz_c = it.col(); - double val = it.value(); - - IJV.push_back(Eigen::Triplet(2 * s.f_n + dz_r, dz_c, val * s.W_11(dz_r))); - IJV.push_back(Eigen::Triplet(2 * s.f_n + dz_r, s.v_n + dz_c, val * s.W_12(dz_r))); - IJV.push_back(Eigen::Triplet(2 * s.f_n + dz_r, 2 * s.v_n + dz_c, val * s.W_13(dz_r))); - - IJV.push_back(Eigen::Triplet(5 * s.f_n + dz_r, dz_c, val * s.W_21(dz_r))); - IJV.push_back(Eigen::Triplet(5 * s.f_n + dz_r, s.v_n + dz_c, val * s.W_22(dz_r))); - IJV.push_back(Eigen::Triplet(5 * s.f_n + dz_r, 2 * s.v_n + dz_c, val * s.W_23(dz_r))); - - IJV.push_back(Eigen::Triplet(8 * s.f_n + dz_r, dz_c, val * s.W_31(dz_r))); - IJV.push_back(Eigen::Triplet(8 * s.f_n + dz_r, s.v_n + dz_c, val * s.W_32(dz_r))); - IJV.push_back(Eigen::Triplet(8 * s.f_n + dz_r, 2 * s.v_n + dz_c, val * s.W_33(dz_r))); - } - } - } - } IGL_INLINE void buildRhs(igl::SLIMData& s, const Eigen::SparseMatrix &A) { @@ -853,10 +320,10 @@ namespace igl W21*R12 + W22*R22];*/ for (int i = 0; i < s.f_n; i++) { - f_rhs(i + 0 * s.f_n) = s.W_11(i) * s.Ri(i, 0) + s.W_12(i) * s.Ri(i, 1); - f_rhs(i + 1 * s.f_n) = s.W_11(i) * s.Ri(i, 2) + s.W_12(i) * s.Ri(i, 3); - f_rhs(i + 2 * s.f_n) = s.W_21(i) * s.Ri(i, 0) + s.W_22(i) * s.Ri(i, 1); - f_rhs(i + 3 * s.f_n) = s.W_21(i) * s.Ri(i, 2) + s.W_22(i) * s.Ri(i, 3); + f_rhs(i + 0 * s.f_n) = s.W(i, 0) * s.Ri(i, 0) + s.W(i, 1) * s.Ri(i, 1); + f_rhs(i + 1 * s.f_n) = s.W(i, 0) * s.Ri(i, 2) + s.W(i, 1) * s.Ri(i, 3); + f_rhs(i + 2 * s.f_n) = s.W(i, 2) * s.Ri(i, 0) + s.W(i, 3) * s.Ri(i, 1); + f_rhs(i + 3 * s.f_n) = s.W(i, 2) * s.Ri(i, 2) + s.W(i, 3) * s.Ri(i, 3); } } else @@ -872,15 +339,15 @@ namespace igl W31*R13 + W32*R23 + W33*R33;];*/ for (int i = 0; i < s.f_n; i++) { - f_rhs(i + 0 * s.f_n) = s.W_11(i) * s.Ri(i, 0) + s.W_12(i) * s.Ri(i, 1) + s.W_13(i) * s.Ri(i, 2); - f_rhs(i + 1 * s.f_n) = s.W_11(i) * s.Ri(i, 3) + s.W_12(i) * s.Ri(i, 4) + s.W_13(i) * s.Ri(i, 5); - f_rhs(i + 2 * s.f_n) = s.W_11(i) * s.Ri(i, 6) + s.W_12(i) * s.Ri(i, 7) + s.W_13(i) * s.Ri(i, 8); - f_rhs(i + 3 * s.f_n) = s.W_21(i) * s.Ri(i, 0) + s.W_22(i) * s.Ri(i, 1) + s.W_23(i) * s.Ri(i, 2); - f_rhs(i + 4 * s.f_n) = s.W_21(i) * s.Ri(i, 3) + s.W_22(i) * s.Ri(i, 4) + s.W_23(i) * s.Ri(i, 5); - f_rhs(i + 5 * s.f_n) = s.W_21(i) * s.Ri(i, 6) + s.W_22(i) * s.Ri(i, 7) + s.W_23(i) * s.Ri(i, 8); - f_rhs(i + 6 * s.f_n) = s.W_31(i) * s.Ri(i, 0) + s.W_32(i) * s.Ri(i, 1) + s.W_33(i) * s.Ri(i, 2); - f_rhs(i + 7 * s.f_n) = s.W_31(i) * s.Ri(i, 3) + s.W_32(i) * s.Ri(i, 4) + s.W_33(i) * s.Ri(i, 5); - f_rhs(i + 8 * s.f_n) = s.W_31(i) * s.Ri(i, 6) + s.W_32(i) * s.Ri(i, 7) + s.W_33(i) * s.Ri(i, 8); + f_rhs(i + 0 * s.f_n) = s.W(i, 0) * s.Ri(i, 0) + s.W(i, 1) * s.Ri(i, 1) + s.W(i, 2) * s.Ri(i, 2); + f_rhs(i + 1 * s.f_n) = s.W(i, 0) * s.Ri(i, 3) + s.W(i, 1) * s.Ri(i, 4) + s.W(i, 2) * s.Ri(i, 5); + f_rhs(i + 2 * s.f_n) = s.W(i, 0) * s.Ri(i, 6) + s.W(i, 1) * s.Ri(i, 7) + s.W(i, 2) * s.Ri(i, 8); + f_rhs(i + 3 * s.f_n) = s.W(i, 3) * s.Ri(i, 0) + s.W(i, 4) * s.Ri(i, 1) + s.W(i, 5) * s.Ri(i, 2); + f_rhs(i + 4 * s.f_n) = s.W(i, 3) * s.Ri(i, 3) + s.W(i, 4) * s.Ri(i, 4) + s.W(i, 5) * s.Ri(i, 5); + f_rhs(i + 5 * s.f_n) = s.W(i, 3) * s.Ri(i, 6) + s.W(i, 4) * s.Ri(i, 7) + s.W(i, 5) * s.Ri(i, 8); + f_rhs(i + 6 * s.f_n) = s.W(i, 6) * s.Ri(i, 0) + s.W(i, 7) * s.Ri(i, 1) + s.W(i, 8) * s.Ri(i, 2); + f_rhs(i + 7 * s.f_n) = s.W(i, 6) * s.Ri(i, 3) + s.W(i, 7) * s.Ri(i, 4) + s.W(i, 8) * s.Ri(i, 5); + f_rhs(i + 8 * s.f_n) = s.W(i, 6) * s.Ri(i, 6) + s.W(i, 7) * s.Ri(i, 7) + s.W(i, 8) * s.Ri(i, 8); } } Eigen::VectorXd uv_flat(s.dim *s.v_n); @@ -894,14 +361,392 @@ namespace igl } } +IGL_INLINE void igl::slim_update_weights_and_closest_rotations_with_jacobians(const Eigen::MatrixXd &Ji, + igl::MappingEnergyType slim_energy, + double exp_factor, + Eigen::MatrixXd &W, + Eigen::MatrixXd &Ri) +{ + const double eps = 1e-8; + double exp_f = exp_factor; + const int dim = (Ji.cols()==4? 2:3); + + if (dim == 2) + { + for (int i = 0; i < Ji.rows(); ++i) + { + typedef Eigen::Matrix2d Mat2; + typedef Eigen::Matrix RMat2; + typedef Eigen::Vector2d Vec2; + Mat2 ji, ri, ti, ui, vi; + Vec2 sing; + Vec2 closest_sing_vec; + RMat2 mat_W; + Vec2 m_sing_new; + double s1, s2; + + ji(0, 0) = Ji(i, 0); + ji(0, 1) = Ji(i, 1); + ji(1, 0) = Ji(i, 2); + ji(1, 1) = Ji(i, 3); + + igl::polar_svd(ji, ri, ti, ui, sing, vi); + + s1 = sing(0); + s2 = sing(1); + + // Update Weights according to energy + switch (slim_energy) + { + case igl::MappingEnergyType::ARAP: + { + m_sing_new << 1, 1; + break; + } + case igl::MappingEnergyType::SYMMETRIC_DIRICHLET: + { + double s1_g = 2 * (s1 - pow(s1, -3)); + double s2_g = 2 * (s2 - pow(s2, -3)); + m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))); + break; + } + case igl::MappingEnergyType::LOG_ARAP: + { + double s1_g = 2 * (log(s1) / s1); + double s2_g = 2 * (log(s2) / s2); + m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))); + break; + } + case igl::MappingEnergyType::CONFORMAL: + { + double s1_g = 1 / (2 * s2) - s2 / (2 * pow(s1, 2)); + double s2_g = 1 / (2 * s1) - s1 / (2 * pow(s2, 2)); + + double geo_avg = sqrt(s1 * s2); + double s1_min = geo_avg; + double s2_min = geo_avg; + + m_sing_new << sqrt(s1_g / (2 * (s1 - s1_min))), sqrt(s2_g / (2 * (s2 - s2_min))); + + // change local step + closest_sing_vec << s1_min, s2_min; + ri = ui * closest_sing_vec.asDiagonal() * vi.transpose(); + break; + } + case igl::MappingEnergyType::EXP_CONFORMAL: + { + double s1_g = 2 * (s1 - pow(s1, -3)); + double s2_g = 2 * (s2 - pow(s2, -3)); + + double geo_avg = sqrt(s1 * s2); + double s1_min = geo_avg; + double s2_min = geo_avg; + + double in_exp = exp_f * ((pow(s1, 2) + pow(s2, 2)) / (2 * s1 * s2)); + double exp_thing = exp(in_exp); + + s1_g *= exp_thing * exp_f; + s2_g *= exp_thing * exp_f; + + m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))); + break; + } + case igl::MappingEnergyType::EXP_SYMMETRIC_DIRICHLET: + { + double s1_g = 2 * (s1 - pow(s1, -3)); + double s2_g = 2 * (s2 - pow(s2, -3)); + + double in_exp = exp_f * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2)); + double exp_thing = exp(in_exp); + + s1_g *= exp_thing * exp_f; + s2_g *= exp_thing * exp_f; + + m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))); + break; + } + } + + if (std::abs(s1 - 1) < eps) m_sing_new(0) = 1; + if (std::abs(s2 - 1) < eps) m_sing_new(1) = 1; + mat_W = ui * m_sing_new.asDiagonal() * ui.transpose(); + + W.row(i) = Eigen::Map>(mat_W.data()); + // 2) Update local step (doesn't have to be a rotation, for instance in case of conformal energy) + Ri.row(i) = Eigen::Map>(ri.data()); + } + } + else + { + typedef Eigen::Matrix Vec3; + typedef Eigen::Matrix Mat3; + typedef Eigen::Matrix RMat3; + Mat3 ji; + Vec3 m_sing_new; + Vec3 closest_sing_vec; + const double sqrt_2 = sqrt(2); + for (int i = 0; i < Ji.rows(); ++i) + { + ji << Ji(i,0), Ji(i,1), Ji(i,2), + Ji(i,3), Ji(i,4), Ji(i,5), + Ji(i,6), Ji(i,7), Ji(i,8); + + Mat3 ri, ti, ui, vi; + Vec3 sing; + igl::polar_svd(ji, ri, ti, ui, sing, vi); + + double s1 = sing(0); + double s2 = sing(1); + double s3 = sing(2); + + // 1) Update Weights + switch (slim_energy) + { + case igl::MappingEnergyType::ARAP: + { + m_sing_new << 1, 1, 1; + break; + } + case igl::MappingEnergyType::LOG_ARAP: + { + double s1_g = 2 * (log(s1) / s1); + double s2_g = 2 * (log(s2) / s2); + double s3_g = 2 * (log(s3) / s3); + m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))), sqrt(s3_g / (2 * (s3 - 1))); + break; + } + case igl::MappingEnergyType::SYMMETRIC_DIRICHLET: + { + double s1_g = 2 * (s1 - pow(s1, -3)); + double s2_g = 2 * (s2 - pow(s2, -3)); + double s3_g = 2 * (s3 - pow(s3, -3)); + m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))), sqrt(s3_g / (2 * (s3 - 1))); + break; + } + case igl::MappingEnergyType::EXP_SYMMETRIC_DIRICHLET: + { + double s1_g = 2 * (s1 - pow(s1, -3)); + double s2_g = 2 * (s2 - pow(s2, -3)); + double s3_g = 2 * (s3 - pow(s3, -3)); + m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))), sqrt(s3_g / (2 * (s3 - 1))); + + double in_exp = exp_f * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2) + pow(s3, 2) + pow(s3, -2)); + double exp_thing = exp(in_exp); + + s1_g *= exp_thing * exp_f; + s2_g *= exp_thing * exp_f; + s3_g *= exp_thing * exp_f; + + m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))), sqrt(s3_g / (2 * (s3 - 1))); + + break; + } + case igl::MappingEnergyType::CONFORMAL: + { + double common_div = 9 * (pow(s1 * s2 * s3, 5. / 3.)); + + double s1_g = (-2 * s2 * s3 * (pow(s2, 2) + pow(s3, 2) - 2 * pow(s1, 2))) / common_div; + double s2_g = (-2 * s1 * s3 * (pow(s1, 2) + pow(s3, 2) - 2 * pow(s2, 2))) / common_div; + double s3_g = (-2 * s1 * s2 * (pow(s1, 2) + pow(s2, 2) - 2 * pow(s3, 2))) / common_div; + + double closest_s = sqrt(pow(s1, 2) + pow(s3, 2)) / sqrt_2; + double s1_min = closest_s; + double s2_min = closest_s; + double s3_min = closest_s; + + m_sing_new << sqrt(s1_g / (2 * (s1 - s1_min))), sqrt(s2_g / (2 * (s2 - s2_min))), sqrt( + s3_g / (2 * (s3 - s3_min))); + + // change local step + closest_sing_vec << s1_min, s2_min, s3_min; + ri = ui * closest_sing_vec.asDiagonal() * vi.transpose(); + break; + } + case igl::MappingEnergyType::EXP_CONFORMAL: + { + // E_conf = (s1^2 + s2^2 + s3^2)/(3*(s1*s2*s3)^(2/3) ) + // dE_conf/ds1 = (-2*(s2*s3)*(s2^2+s3^2 -2*s1^2) ) / (9*(s1*s2*s3)^(5/3)) + // Argmin E_conf(s1): s1 = sqrt(s1^2+s2^2)/sqrt(2) + double common_div = 9 * (pow(s1 * s2 * s3, 5. / 3.)); + + double s1_g = (-2 * s2 * s3 * (pow(s2, 2) + pow(s3, 2) - 2 * pow(s1, 2))) / common_div; + double s2_g = (-2 * s1 * s3 * (pow(s1, 2) + pow(s3, 2) - 2 * pow(s2, 2))) / common_div; + double s3_g = (-2 * s1 * s2 * (pow(s1, 2) + pow(s2, 2) - 2 * pow(s3, 2))) / common_div; + + double in_exp = exp_f * ((pow(s1, 2) + pow(s2, 2) + pow(s3, 2)) / (3 * pow((s1 * s2 * s3), 2. / 3)));; + double exp_thing = exp(in_exp); + + double closest_s = sqrt(pow(s1, 2) + pow(s3, 2)) / sqrt_2; + double s1_min = closest_s; + double s2_min = closest_s; + double s3_min = closest_s; + + s1_g *= exp_thing * exp_f; + s2_g *= exp_thing * exp_f; + s3_g *= exp_thing * exp_f; + + m_sing_new << sqrt(s1_g / (2 * (s1 - s1_min))), sqrt(s2_g / (2 * (s2 - s2_min))), sqrt( + s3_g / (2 * (s3 - s3_min))); + + // change local step + closest_sing_vec << s1_min, s2_min, s3_min; + ri = ui * closest_sing_vec.asDiagonal() * vi.transpose(); + } + } + if (std::abs(s1 - 1) < eps) m_sing_new(0) = 1; + if (std::abs(s2 - 1) < eps) m_sing_new(1) = 1; + if (std::abs(s3 - 1) < eps) m_sing_new(2) = 1; + RMat3 mat_W; + mat_W = ui * m_sing_new.asDiagonal() * ui.transpose(); + + W.row(i) = Eigen::Map>(mat_W.data()); + // 2) Update closest rotations (not rotations in case of conformal energy) + Ri.row(i) = Eigen::Map>(ri.data()); + } // for loop end + + } // if dim end + +} + +IGL_INLINE void igl::slim_buildA(const Eigen::SparseMatrix &Dx, + const Eigen::SparseMatrix &Dy, + const Eigen::SparseMatrix &Dz, + const Eigen::MatrixXd &W, +std::vector > & IJV) +{ + const int dim = (W.cols() == 4) ? 2 : 3; + const int f_n = W.rows(); + const int v_n = Dx.cols(); + + // formula (35) in paper + if (dim == 2) + { + IJV.reserve(4 * (Dx.outerSize() + Dy.outerSize())); + + /*A = [W11*Dx, W12*Dx; + W11*Dy, W12*Dy; + W21*Dx, W22*Dx; + W21*Dy, W22*Dy];*/ + for (int k = 0; k < Dx.outerSize(); ++k) + { + for (Eigen::SparseMatrix::InnerIterator it(Dx, k); it; ++it) + { + int dx_r = it.row(); + int dx_c = it.col(); + double val = it.value(); + + IJV.push_back(Eigen::Triplet(dx_r, dx_c, val * W(dx_r, 0))); + IJV.push_back(Eigen::Triplet(dx_r, v_n + dx_c, val * W(dx_r, 1))); + + IJV.push_back(Eigen::Triplet(2 * f_n + dx_r, dx_c, val * W(dx_r, 2))); + IJV.push_back(Eigen::Triplet(2 * f_n + dx_r, v_n + dx_c, val * W(dx_r, 3))); + } + } + + for (int k = 0; k < Dy.outerSize(); ++k) + { + for (Eigen::SparseMatrix::InnerIterator it(Dy, k); it; ++it) + { + int dy_r = it.row(); + int dy_c = it.col(); + double val = it.value(); + + IJV.push_back(Eigen::Triplet(f_n + dy_r, dy_c, val * W(dy_r, 0))); + IJV.push_back(Eigen::Triplet(f_n + dy_r, v_n + dy_c, val * W(dy_r, 1))); + + IJV.push_back(Eigen::Triplet(3 * f_n + dy_r, dy_c, val * W(dy_r, 2))); + IJV.push_back(Eigen::Triplet(3 * f_n + dy_r, v_n + dy_c, val * W(dy_r, 3))); + } + } + } + else + { + + /*A = [W11*Dx, W12*Dx, W13*Dx; + W11*Dy, W12*Dy, W13*Dy; + W11*Dz, W12*Dz, W13*Dz; + W21*Dx, W22*Dx, W23*Dx; + W21*Dy, W22*Dy, W23*Dy; + W21*Dz, W22*Dz, W23*Dz; + W31*Dx, W32*Dx, W33*Dx; + W31*Dy, W32*Dy, W33*Dy; + W31*Dz, W32*Dz, W33*Dz;];*/ + IJV.reserve(9 * (Dx.outerSize() + Dy.outerSize() + Dz.outerSize())); + for (int k = 0; k < Dx.outerSize(); k++) + { + for (Eigen::SparseMatrix::InnerIterator it(Dx, k); it; ++it) + { + int dx_r = it.row(); + int dx_c = it.col(); + double val = it.value(); + + IJV.push_back(Eigen::Triplet(dx_r, dx_c, val * W(dx_r, 0))); + IJV.push_back(Eigen::Triplet(dx_r, v_n + dx_c, val * W(dx_r, 1))); + IJV.push_back(Eigen::Triplet(dx_r, 2 * v_n + dx_c, val * W(dx_r, 2))); + + IJV.push_back(Eigen::Triplet(3 * f_n + dx_r, dx_c, val * W(dx_r, 3))); + IJV.push_back(Eigen::Triplet(3 * f_n + dx_r, v_n + dx_c, val * W(dx_r, 4))); + IJV.push_back(Eigen::Triplet(3 * f_n + dx_r, 2 * v_n + dx_c, val * W(dx_r, 5))); + + IJV.push_back(Eigen::Triplet(6 * f_n + dx_r, dx_c, val * W(dx_r, 6))); + IJV.push_back(Eigen::Triplet(6 * f_n + dx_r, v_n + dx_c, val * W(dx_r, 7))); + IJV.push_back(Eigen::Triplet(6 * f_n + dx_r, 2 * v_n + dx_c, val * W(dx_r, 8))); + } + } + + for (int k = 0; k < Dy.outerSize(); k++) + { + for (Eigen::SparseMatrix::InnerIterator it(Dy, k); it; ++it) + { + int dy_r = it.row(); + int dy_c = it.col(); + double val = it.value(); + + IJV.push_back(Eigen::Triplet(f_n + dy_r, dy_c, val * W(dy_r, 0))); + IJV.push_back(Eigen::Triplet(f_n + dy_r, v_n + dy_c, val * W(dy_r, 1))); + IJV.push_back(Eigen::Triplet(f_n + dy_r, 2 * v_n + dy_c, val * W(dy_r, 2))); + + IJV.push_back(Eigen::Triplet(4 * f_n + dy_r, dy_c, val * W(dy_r, 3))); + IJV.push_back(Eigen::Triplet(4 * f_n + dy_r, v_n + dy_c, val * W(dy_r, 4))); + IJV.push_back(Eigen::Triplet(4 * f_n + dy_r, 2 * v_n + dy_c, val * W(dy_r, 5))); + + IJV.push_back(Eigen::Triplet(7 * f_n + dy_r, dy_c, val * W(dy_r, 6))); + IJV.push_back(Eigen::Triplet(7 * f_n + dy_r, v_n + dy_c, val * W(dy_r, 7))); + IJV.push_back(Eigen::Triplet(7 * f_n + dy_r, 2 * v_n + dy_c, val * W(dy_r, 8))); + } + } + + for (int k = 0; k < Dz.outerSize(); k++) + { + for (Eigen::SparseMatrix::InnerIterator it(Dz, k); it; ++it) + { + int dz_r = it.row(); + int dz_c = it.col(); + double val = it.value(); + + IJV.push_back(Eigen::Triplet(2 * f_n + dz_r, dz_c, val * W(dz_r, 0))); + IJV.push_back(Eigen::Triplet(2 * f_n + dz_r, v_n + dz_c, val * W(dz_r, 1))); + IJV.push_back(Eigen::Triplet(2 * f_n + dz_r, 2 * v_n + dz_c, val * W(dz_r, 2))); + + IJV.push_back(Eigen::Triplet(5 * f_n + dz_r, dz_c, val * W(dz_r, 3))); + IJV.push_back(Eigen::Triplet(5 * f_n + dz_r, v_n + dz_c, val * W(dz_r, 4))); + IJV.push_back(Eigen::Triplet(5 * f_n + dz_r, 2 * v_n + dz_c, val * W(dz_r, 5))); + + IJV.push_back(Eigen::Triplet(8 * f_n + dz_r, dz_c, val * W(dz_r, 6))); + IJV.push_back(Eigen::Triplet(8 * f_n + dz_r, v_n + dz_c, val * W(dz_r, 7))); + IJV.push_back(Eigen::Triplet(8 * f_n + dz_r, 2 * v_n + dz_c, val * W(dz_r, 8))); + } + } + } +} /// Slim Implementation IGL_INLINE void igl::slim_precompute( const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const Eigen::MatrixXd &V_init, - SLIMData &data, - SLIMData::SLIM_ENERGY slim_energy, + igl::SLIMData &data, + igl::MappingEnergyType slim_energy, Eigen::VectorXi &b, Eigen::MatrixXd &bc, double soft_p) @@ -934,7 +779,7 @@ IGL_INLINE void igl::slim_precompute( data.energy = igl::slim::compute_energy(data,data.V_o) / data.mesh_area; } -IGL_INLINE Eigen::MatrixXd igl::slim_solve(SLIMData &data, int iter_num) +IGL_INLINE Eigen::MatrixXd igl::slim_solve(igl::SLIMData &data, int iter_num) { for (int i = 0; i < iter_num; i++) { @@ -942,7 +787,7 @@ IGL_INLINE Eigen::MatrixXd igl::slim_solve(SLIMData &data, int iter_num) dest_res = data.V_o; // Solve Weighted Proxy - igl::slim::update_weights_and_closest_rotations(data,data.V, data.F, dest_res); + igl::slim::update_weights_and_closest_rotations(data, dest_res); igl::slim::solve_weighted_arap(data,data.V, data.F, dest_res, data.b, data.bc); double old_energy = data.energy; @@ -956,6 +801,7 @@ IGL_INLINE Eigen::MatrixXd igl::slim_solve(SLIMData &data, int iter_num) return data.V_o; } + #ifdef IGL_STATIC_LIBRARY // Explicit template instantiation #endif diff --git a/include/igl/slim.h b/include/igl/slim.h index 1965e2bba..aebfb2688 100644 --- a/include/igl/slim.h +++ b/include/igl/slim.h @@ -9,6 +9,7 @@ #define SLIM_H #include "igl_inline.h" +#include "MappingEnergyType.h" #include #include @@ -29,16 +30,7 @@ struct SLIMData // Input Eigen::MatrixXd V; // #V by 3 list of mesh vertex positions Eigen::MatrixXi F; // #F by 3/3 list of mesh faces (triangles/tets) - enum SLIM_ENERGY - { - ARAP, - LOG_ARAP, - SYMMETRIC_DIRICHLET, - CONFORMAL, - EXP_CONFORMAL, - EXP_SYMMETRIC_DIRICHLET - }; - SLIM_ENERGY slim_energy; + MappingEnergyType slim_energy; // Optional Input // soft constraints @@ -64,9 +56,7 @@ struct SLIMData Eigen::VectorXd WGL_M; Eigen::VectorXd rhs; Eigen::MatrixXd Ri,Ji; - Eigen::VectorXd W_11; Eigen::VectorXd W_12; Eigen::VectorXd W_13; - Eigen::VectorXd W_21; Eigen::VectorXd W_22; Eigen::VectorXd W_23; - Eigen::VectorXd W_31; Eigen::VectorXd W_32; Eigen::VectorXd W_33; + Eigen::MatrixXd W; Eigen::SparseMatrix Dx,Dy,Dz; int f_n,v_n; bool first_solve; @@ -94,7 +84,7 @@ IGL_INLINE void slim_precompute( const Eigen::MatrixXi& F, const Eigen::MatrixXd& V_init, SLIMData& data, - SLIMData::SLIM_ENERGY slim_energy, + MappingEnergyType slim_energy, Eigen::VectorXi& b, Eigen::MatrixXd& bc, double soft_p); @@ -104,6 +94,18 @@ IGL_INLINE void slim_precompute( // V_o (in SLIMData): #V by dim list of mesh vertex positions IGL_INLINE Eigen::MatrixXd slim_solve(SLIMData& data, int iter_num); +// Internal Routine. Exposed for Integration with SCAF +IGL_INLINE void slim_update_weights_and_closest_rotations_with_jacobians(const Eigen::MatrixXd &Ji, + igl::MappingEnergyType slim_energy, + double exp_factor, + Eigen::MatrixXd &W, + Eigen::MatrixXd &Ri); + +IGL_INLINE void slim_buildA(const Eigen::SparseMatrix &Dx, + const Eigen::SparseMatrix &Dy, + const Eigen::SparseMatrix &Dz, + const Eigen::MatrixXd &W, + std::vector > & IJV); } // END NAMESPACE #ifndef IGL_STATIC_LIBRARY diff --git a/include/igl/sparse_cached.h b/include/igl/sparse_cached.h index 40d288962..c5a66c99f 100644 --- a/include/igl/sparse_cached.h +++ b/include/igl/sparse_cached.h @@ -39,7 +39,7 @@ namespace igl // Example: // Eigen::SparseMatrix A; // std::vector > IJV; - // buildA(IJV); + // slim_buildA(IJV); // if (A.rows() == 0) // { // A = Eigen::SparseMatrix(rows,cols); diff --git a/include/igl/topological_hole_fill.cpp b/include/igl/topological_hole_fill.cpp new file mode 100644 index 000000000..c09963a7c --- /dev/null +++ b/include/igl/topological_hole_fill.cpp @@ -0,0 +1,55 @@ +// This file is part of libigl, a simple c++ geometry processing library. +// +// Copyright (C) 2018 Zhongshi Jiang +// +// 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/. + +#include "topological_hole_fill.h" + template < + typename DerivedF, + typename Derivedb, + typename VectorIndex, + typename DerivedF_filled> +IGL_INLINE void igl::topological_hole_fill( + const Eigen::MatrixBase & F, + const Eigen::MatrixBase & b, + const std::vector & holes, + Eigen::PlainObjectBase &F_filled) +{ + int n_filled_faces = 0; + int num_holes = holes.size(); + int real_F_num = F.rows(); + const int V_rows = F.maxCoeff()+1; + + for (int i = 0; i < num_holes; i++) + n_filled_faces += holes[i].size(); + F_filled.resize(n_filled_faces + real_F_num, 3); + F_filled.topRows(real_F_num) = F; + + int new_vert_id = V_rows; + int new_face_id = real_F_num; + + for (int i = 0; i < num_holes; i++, new_vert_id++) + { + int cur_bnd_size = holes[i].size(); + int it = 0; + int back = holes[i].size() - 1; + F_filled.row(new_face_id++) << holes[i][it], holes[i][back], new_vert_id; + while (it != back) + { + F_filled.row(new_face_id++) + << holes[i][(it + 1)], + holes[i][(it)], new_vert_id; + it++; + } + } + assert(new_face_id == F_filled.rows()); + assert(new_vert_id == V_rows + num_holes); + +} + +#ifdef IGL_STATIC_LIBRARY +template void igl::topological_hole_fill, Eigen::Matrix, std::vector > >(Eigen::MatrixBase > const&, Eigen::MatrixBase > const&, std::vector >, std::allocator > > > const&, Eigen::PlainObjectBase >&); +#endif \ No newline at end of file diff --git a/include/igl/topological_hole_fill.h b/include/igl/topological_hole_fill.h new file mode 100644 index 000000000..945e127e8 --- /dev/null +++ b/include/igl/topological_hole_fill.h @@ -0,0 +1,44 @@ +// This file is part of libigl, a simple c++ geometry processing library. +// +// Copyright (C) 2018 Zhongshi Jiang +// +// 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/. +#ifndef IGL_TOPOLOGICAL_HOLE_FILL_H +#define IGL_TOPOLOGICAL_HOLE_FILL_H +#include "igl_inline.h" +#include +#include +namespace igl +{ + + // Topological fill hole on a mesh, with one additional vertex each hole + // Index of new abstract vertices will be F.maxCoeff() + (index of hole) + // + // Inputs: + // F #F by simplex-size list of element indices + // b #b boundary indices to preserve + // holes vector of hole loops to fill + // Outputs: + // F_filled input F stacked with filled triangles. + // + template < + typename DerivedF, + typename Derivedb, + typename VectorIndex, + typename DerivedF_filled> +IGL_INLINE void topological_hole_fill( + const Eigen::MatrixBase & F, + const Eigen::MatrixBase & b, + const std::vector & holes, + Eigen::PlainObjectBase &F_filled); + +} + + +#ifndef IGL_STATIC_LIBRARY +# include "topological_hole_fill.cpp" +#endif + +#endif \ No newline at end of file diff --git a/tutorial/710_SLIM/CMakeLists.txt b/tutorial/709_SLIM/CMakeLists.txt similarity index 100% rename from tutorial/710_SLIM/CMakeLists.txt rename to tutorial/709_SLIM/CMakeLists.txt diff --git a/tutorial/710_SLIM/main.cpp b/tutorial/709_SLIM/main.cpp similarity index 92% rename from tutorial/710_SLIM/main.cpp rename to tutorial/709_SLIM/main.cpp index 0521dff88..eaa2ca3b5 100644 --- a/tutorial/710_SLIM/main.cpp +++ b/tutorial/709_SLIM/main.cpp @@ -1,15 +1,16 @@ #include -#include "igl/slim.h" +#include -#include "igl/vertex_components.h" -#include "igl/readOBJ.h" -#include "igl/writeOBJ.h" -#include "igl/Timer.h" +#include +#include +#include +#include -#include "igl/boundary_loop.h" -#include "igl/map_vertices_to_circle.h" -#include "igl/harmonic.h" +#include +#include +#include +#include #include #include #include @@ -96,9 +97,9 @@ void param_2d_demo_iter(igl::opengl::glfw::Viewer& viewer) { cout << "initialized parametrization" << endl; - sData.slim_energy = igl::SLIMData::SYMMETRIC_DIRICHLET; + sData.slim_energy = igl::MappingEnergyType::SYMMETRIC_DIRICHLET; Eigen::VectorXi b; Eigen::MatrixXd bc; - slim_precompute(V,F,uv_init,sData, igl::SLIMData::SYMMETRIC_DIRICHLET, b,bc,0); + slim_precompute(V,F,uv_init,sData, igl::MappingEnergyType::SYMMETRIC_DIRICHLET, b,bc,0); uv_scale_param = 15 * (1./sqrt(sData.mesh_area)); viewer.data().set_mesh(V, F); @@ -109,10 +110,12 @@ void param_2d_demo_iter(igl::opengl::glfw::Viewer& viewer) { first_iter = false; } else { + timer.start(); slim_solve(sData,1); // 1 iter viewer.data().set_uv(sData.V_o*uv_scale_param); - cout << "time = " << timer.getElapsedTime() << endl; } + cout << "time = " << timer.getElapsedTime() << endl; + cout << "energy = " << sData.energy << endl; } void soft_const_demo_iter(igl::opengl::glfw::Viewer& viewer) { @@ -127,7 +130,7 @@ void soft_const_demo_iter(igl::opengl::glfw::Viewer& viewer) { Eigen::VectorXi b; Eigen::MatrixXd bc; get_soft_constraint_for_circle(V_0,F,b,bc); double soft_const_p = 1e5; - slim_precompute(V,F,V_0,sData,igl::SLIMData::SYMMETRIC_DIRICHLET,b,bc,soft_const_p); + slim_precompute(V,F,V_0,sData,igl::MappingEnergyType::SYMMETRIC_DIRICHLET,b,bc,soft_const_p); viewer.data().set_mesh(V, F); viewer.core.align_camera_center(V,F); @@ -144,6 +147,7 @@ void soft_const_demo_iter(igl::opengl::glfw::Viewer& viewer) { void deform_3d_demo_iter(igl::opengl::glfw::Viewer& viewer) { if (first_iter) { + timer.start(); igl::readOBJ(TUTORIAL_SHARED_PATH "/cube_40k.obj", V, F); Eigen::MatrixXd V_0 = V; @@ -152,16 +156,19 @@ void deform_3d_demo_iter(igl::opengl::glfw::Viewer& viewer) { double soft_const_p = 1e5; sData.exp_factor = 5.0; - slim_precompute(V,F,V_0,sData,igl::SLIMData::EXP_CONFORMAL,b,bc,soft_const_p); + slim_precompute(V,F,V_0,sData,igl::MappingEnergyType::EXP_CONFORMAL,b,bc,soft_const_p); //cout << "precomputed" << endl; first_iter = false; display_3d_mesh(viewer); } else { + timer.start(); slim_solve(sData,1); // 1 iter display_3d_mesh(viewer); } + cout << "time = " << timer.getElapsedTime() << endl; + cout << "energy = " << sData.energy << endl; } void display_3d_mesh(igl::opengl::glfw::Viewer& viewer) { diff --git a/tutorial/710_SCAF/CMakeLists.txt b/tutorial/710_SCAF/CMakeLists.txt new file mode 100644 index 000000000..7261d33ee --- /dev/null +++ b/tutorial/710_SCAF/CMakeLists.txt @@ -0,0 +1,5 @@ +get_filename_component(PROJECT_NAME ${CMAKE_CURRENT_SOURCE_DIR} NAME) +project(${PROJECT_NAME}) + +add_executable(${PROJECT_NAME}_bin main.cpp) +target_link_libraries(${PROJECT_NAME}_bin igl::core igl::opengl igl::opengl_glfw igl::triangle tutorials) diff --git a/tutorial/710_SCAF/main.cpp b/tutorial/710_SCAF/main.cpp new file mode 100755 index 000000000..5a583c869 --- /dev/null +++ b/tutorial/710_SCAF/main.cpp @@ -0,0 +1,127 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "tutorial_shared_path.h" + +Eigen::MatrixXd V; +Eigen::MatrixXi F; +Eigen::MatrixXd V_uv; +igl::Timer timer; +igl::SCAFData scaf_data; + +bool show_uv = false; +float uv_scale = 0.2; + +bool key_down(igl::opengl::glfw::Viewer& viewer, unsigned char key, int modifier) +{ + if (key == '1') + show_uv = false; + else if (key == '2') + show_uv = true; + + if (key == ' ') + { + timer.start(); + igl::scaf_solve(scaf_data, 1); + std::cout << "time = " << timer.getElapsedTime() << std::endl; + } + + const auto& V_uv = uv_scale * scaf_data.w_uv.topRows(V.rows()); + if (show_uv) + { + viewer.data().clear(); + viewer.data().set_mesh(V_uv,F); + viewer.data().set_uv(V_uv); + viewer.core.align_camera_center(V_uv,F); + } + else + { + viewer.data().set_mesh(V,F); + viewer.data().set_uv(V_uv); + viewer.core.align_camera_center(V,F); + } + + viewer.data().compute_normals(); + + return false; +} + +int main(int argc, char *argv[]) +{ + using namespace std; + // Load a mesh in OFF format + igl::readOBJ(TUTORIAL_SHARED_PATH "/camel_b.obj", V, F); + + Eigen::MatrixXd bnd_uv, uv_init; + + Eigen::VectorXd M; + igl::doublearea(V, F, M); + std::vector> all_bnds; + igl::boundary_loop(F, all_bnds); + + // Heuristic primary boundary choice: longest + auto primary_bnd = std::max_element(all_bnds.begin(), all_bnds.end(), [](const std::vector &a, const std::vector &b) { return a.size()(primary_bnd->data(), primary_bnd->size()); + + igl::map_vertices_to_circle(V, bnd, bnd_uv); + bnd_uv *= sqrt(M.sum() / (2 * igl::PI)); + if (all_bnds.size() == 1) + { + if (bnd.rows() == V.rows()) // case: all vertex on boundary + { + uv_init.resize(V.rows(), 2); + for (int i = 0; i < bnd.rows(); i++) + uv_init.row(bnd(i)) = bnd_uv.row(i); + } + else + { + igl::harmonic(V, F, bnd, bnd_uv, 1, uv_init); + if (igl::flipped_triangles(uv_init, F).size() != 0) + igl::harmonic(F, bnd, bnd_uv, 1, uv_init); // fallback uniform laplacian + } + } + else + { + // if there is a hole, fill it and erase additional vertices. + all_bnds.erase(primary_bnd); + Eigen::MatrixXi F_filled; + igl::topological_hole_fill(F, bnd, all_bnds, F_filled); + igl::harmonic(F_filled, bnd, bnd_uv ,1, uv_init); + uv_init = uv_init.topRows(V.rows()); + } + + Eigen::VectorXi b; Eigen::MatrixXd bc; + igl::scaf_precompute(V, F, uv_init, scaf_data, igl::MappingEnergyType::SYMMETRIC_DIRICHLET, b, bc, 0); + + // Plot the mesh + igl::opengl::glfw::Viewer viewer; + viewer.data().set_mesh(V, F); + const auto& V_uv = uv_scale * scaf_data.w_uv.topRows(V.rows()); + viewer.data().set_uv(V_uv); + viewer.callback_key_down = &key_down; + + // Enable wireframe + viewer.data().show_lines = true; + + // Draw checkerboard texture + viewer.data().show_texture = true; + + + std::cerr << "Press space for running an iteration." << std::endl; + std::cerr << "Press 1 for Mesh 2 for UV" << std::endl; + + // Launch the viewer + viewer.launch(); +} diff --git a/tutorial/CMakeLists.txt b/tutorial/CMakeLists.txt index 587291606..947fb86cc 100644 --- a/tutorial/CMakeLists.txt +++ b/tutorial/CMakeLists.txt @@ -144,7 +144,8 @@ if(TUTORIALS_CHAPTER7) endif() add_subdirectory("707_SweptVolume") add_subdirectory("708_Picking") - add_subdirectory("710_SLIM") + add_subdirectory("709_SLIM") + add_subdirectory("710_SCAF") add_subdirectory("711_Subdivision") add_subdirectory("712_DataSmoothing") add_subdirectory("713_ShapeUp") diff --git a/tutorial/images/710_SCAF.png b/tutorial/images/710_SCAF.png new file mode 100644 index 0000000000000000000000000000000000000000..0d5f21bca3dff436750d4e1500fad9344c6f36c1 GIT binary patch literal 167471 zcmX_nbwE?!`|uPM73nZ28G_Q%4GN=6Iz~u0NR1w-bc%E<-J`p^ksLjGFiLVX@9^{e zy??-kd(VB&)90xh{7F#?_vx#r00012MjE6706ZiG0PYt*zK^-X^MNvBu=saq6_NZ?hMmUcm8zxVF6H6Bf+8{-0J+2 zhfGu)(zfg#Q?4DJvYSqD?~H*a~JRs+JnerSE+H;lO}^O2z>FRyoK zMB?fh0Pus^i=p|m`X3iPIj~?o8|Haw?5co{b*QwiF~r#%^S^HxJ83Lwe*V3bB+O<% zL&oYe%G-g!-yWQuVjq(ZOI%OAW-I-J-o0fK$`cKLDJ9}lv0JWcs?nUR!yle$9W1>Q z5t~n|26r&Z<^tFB<3h{#R`W>_&;5Cbm8)S3(63GSEH9EXl6~}jv*;wX4o*od+MI&l z*C?n0O2zp`^DG72y7^!I+Fhc+1w0|5Rpjdb{i*hgvGOjq+qefJ8sAf#Rhz6!+?u}U z$))pVV<=s&xN7kucrNROr-lx@bzx3H?Qq=TX{@*1Yadk?oDT9meGh|}0inT3bq1i3 z+y%EVtoLZcex&f>S!$xixiBvRy5_Fwazv8Gd$18bTz}!0Db;-3d7UEscY+wNSg+h` zC?)B>0xx%-l;XJ3#Yysr$l+K_1GRe**iI+8i3o5;nVogUOICuPWy+%gKw9|j20C#8 zDitN%gC#(3s!E&!#TMUEpGpQZD*)9vxoj=%W7hft#p4$~h!}M2*8QIF_py)H7Bh*m zwi+HPDuTSEtNgnki?vU)@MYl_dA)g`ch-Q7QH@5qa<rP8uD0VHACP6+ZD$TFhvk zG?YwJO>zgK-GQdi|B)B2;#(Fx18%q)MI6vpTPv|?B`;3XGg9>)6uVNA~c($ zV>RSWqS%d$CW8|apE1P7YUtJ-Qf@&Uzq`DJ-STDh+ZfRTJInS-B}_mhPM}x?{o4cf zhdq^B1K4M{FA2p>Kth)1&vm`Z)=z>7X)~mOzdtK0m+i_8*~QC%(kdI1?t|7q;^e%f z?$=*WcR9q>KMrhXqH?T!)H;@5GtC+Ia6a?@_Hp;*WD(hl1EF3=gUfzKNq-s7F&bDM z8L}M_>*iYy+D5wv_qwD?SALyyDIkbi1?kc7jApE`XwSGq@f}~?y9#)d6E!%;Cxun{ zBqj6SFBecJkVbWzBf)3*emOHyaC$=F2<)ICZauZ)*I{V@Te7Lzb?+8go#eCNiJ}5B zjzN>@^82Quas9!#Q>i)!-veigMc8`^UI8U+KthA6Oi(_Z2)hAI{J&r# zuY87I@`_g8F5v;oy$*_SM%=BP)MJ7IB1LWvLlu!0mI_}dTNIm)+H!ef)d-X$o4@li zCwUJm;P1|{aZjlgM;2TiT?{PDFY@XMM+J$2HPXEoJ2_deE`uT@o-h7AR%H`TwJWq2 z@w}ZePiMms1WjaL`hx8nOvfYl4KDwZD+%4|#woR2OtkH zt5X7SyIHZVuV-C0gduuzs-VFG&W#+Fxfqy#rnOAtLFh_VPAuV*C(#5_a2|~szVhr; zO=9JpP8(2=G&?y^;SDI6He7`50}&B>_Y2zYpm*#^UA5gWDm+dKL~qy2#qnpak6aC) zEB?G9RHc`OmyjURFR8sgRJUxG98@R8&U49%6Wr%)Z@o6^;H-Q{pIZ{SjQ5g8 z#RXaQ>LwEYPQZ|qKGRG7KUjcNa{Fp-wcI}j_;$oEo>H#gBgzN`C)r$mCUp*s)c76#5TMX8! zHs2l&%y6~{IH#)`J#07KMwYN$Y)_(1b8`tsvo-acPrKlgmcMhZmyIJmHnQ=X&+QJT zuRU)0toMG=h+iiaf;wLiYF~9A>5M2ndV@}^HKHb)HCyxt;3cJAyN`~yIkx9ur(52C zztLqIpHgu7oWl1KsN|3=B0G(ocUz24`HGd8y-wq! zvCEu0M`Axiby|ot#u5rpJhSevy!CzGsCnx(Y}|Lff+VpOUf6Twnwa{h^cO$Wg}>OE z%JweueM>8-o~4Mxy$fzW{T^_ddV0M!isjAuyrk}MgtwdRD>-dRu!j0e-m3|5d@V-( z#ptx<7h7}114ubdmYdVi50N6P5h6@!TbPf<>=(G2SM*)8XfJVWcZ;iUta5qLy8I_EqrpJrva|H%)nDYs9fb3{XWAW8>L&f4_ANKrb5AY_NBa{I@rjRg@*E~Mmw815 z+)>z!+%QCI+HQE~6exsdRYAt|)WGLbukIc`VK93;$mo>PsWmH->bey?TA_NolP1p5 zj^BJaC6#)K-+Qc+{3>7NXh(2K;lcFD3nJrJ*PY80xPR!sL_G4HVwYqF7Nov^&w(mk z+Tn-3jhm?0AmZ+^d(=en(nZVI{bT02acDd=Jyl&+mi;ABggRYzP*!b`UL{h?Dr@iK z47PX8h1A%3$V|{CYBeNL+9n$E>GWf`FveUtCf};7Ztqxy4PLI#JvZ3*UTez~fv)G* z`O^`P2-oKLo|^9qXQ4e|qJ5>JW*`p#Ty=EK$n@ZB0e^AmHc5Q--2A%ZvgEeMjFU_t4z(6KuW&&4U%Zg9P$`GQ zdR&fOp}K_cF9Txlv42O#J~-b5Z%xl>M0bada9k=1U8i^lI&&$SHu|z7^&Fb6H)#<* zM=q~fjy9u5JdpW4($zfbXxd+19f8kuI9-O9dniaer7b10BtY(*tSy_3*y-L|Fr>ok zjCdSx~84WU&d)N^~~&ip9qPAR{ZC z5%rkK6nx6l(L`CyZVakJj@WTdst9Li#YO0?C zDu2gJ#H?f}Kz*!K{V}?EHgp%*DLj6}RF{$Acz#GKSpOFA0d%|(jc&Z^8x)q>Ni^zX z6TOXzn!8L#FPbqej$YsQMaWoLSvic-r;D0Oxjq1(+A;1{tc|gA1QVKf4kldmRT?~x z{scu*Nm07g>P1@-10eSp7#JiZl2ufaD109IQr^S3w!DIZwUEOQQ@Qxv%`{A?3&?(O zUqXUh)$4`NBLF~JK>atPE)N#wFS0Tg7l`v+bo}cc(8$02_X#=A&VR3%B27Xc4-a4- z$-H9`{#>l3bM=Q-^nVX2D-Bu*FG-df?;>r#l7@}fgNQrGRNv0$FmLq9fCxk}@ZaFj za0I2_Mb|$+{tL=MeqsOp5%#|;7s&7Z%^i@G|1Lsr$iy%&emHEqz{KdjZ`$98{Rbe5 zVCC*c7L8bU0_6e@W)llv_D2Vtn6*S7Z{A}bt@&LU-eC)p2K@&s?0;7-&<|0}7O!P! zd|qIds}N$iiNW#7zn3}M&hB88i~WaR>CH>bN=w7*yPN;_to`2-BEtXHWc^1PsoxzG z8Ck#oQMzOB5JN`ge_QKzzW%ot2hsQcka}?Wj~@(S|1QsI{zF0f@jo2K2xy-DyZ?3# z|4u552`#u7SXrRX*->xICcb=_1gBIu%poL{jFgm3zrswvfbF{V8rylB+hf`2PI0E) zUqkN%Ko<679`J{e5+{YN<-%mAeCqYZjw#f1hi1Ens{G{h``VpL>RgvC zuFtK$t57)W{i}l~-T&kpNN~rp%J%P(7xZD&%E-@2@I<%u@(X`sI!?~1A+YXWPiMT( z=3h_lm<9m&xG`zuinyVdzD1~4;V(K|u3)F7Y2>c?9nju*+v5X)PleACE{D(;XcTSh zFyOrAPKN?MU%dqAh%;<|@6&d{3jhvXO^BZ6Wgv>A0D?4xJGXag1S<6*1yqh{BDB}n z=bq7H(+-&|?v2%Pu))t+1Rs{F{!a}boOb|zRRIa1tHpR#k+mh3Sw$h&?#tQV+QYSm zpZs{&c4d2~Z0#6(@twLrd6L@L4MW275gSdHQTx2R_^pZo?ElmjL?DRO3dt){vu|1W zeeMVfB{LT~VE$(AOV;6sk)iLKZ;xLAh?Wne?IIc+q9be?dfodq564Fa!X1oVx9zPY zK*TUzqHyP65Nd-Q%QwgH_Mesm{V;a)C#4N*&|qZYM=GU*Ur4CXh)8-1BvrfRYfO26 zzwkJe5nN}5*Hu%1QqSG%WQrF$YOh=rUPu$Z+zD#Nx{CS-3WI`ZNPmBLp@#P@DvPG9 zflO6jacy=$lTg{?=9Y^IwvI1Zm}K)ksMjS5LjH^UKV>Ei6T$l6-OVSFyg_EKU*>p1 z0R+Ap`^>H>uU^R!^+@mg*~?#D?|sNIrZV&T+aq7wBCVUrl`cfXnw@ZG)s#?;)u*Dep8&h zC8O4kvX-X1HCDJenfKqIDpfsKurSv21M~pv561C=B8a`N7Z^cp1cg*DPU#{Q$Z5#ZxI}?K6(GJ*FzXnYFw`A=rMhB|$KnHMV{daz*-%ivwsOEb7efr3F$mzg<%bT&JT^InO{m=JdR6ol0 zQfyg36iyzI(X=?_*Ou;`w0Xt@==^-=r|uxk8vhNjh8bwwo?%RAo5J$>W2HtT86aJ7 zUXICh-$y>cAK+Q*-BQ4tRE$>P!nckasZpY88g7!Qw4~&`&!1iV8g6^CDs^+&tnY%S z*l5(p$#(Stv~?5aeM&wK!El7`KMw<(zI_P@%9ArkkA?EQ(H9}Ia-iS6hZSjIN$RVY z67+8o#t}WGiliuwkQWt7Ya?;7FW&Zw#5f#JvuDpzo!YsyKIbR(VERSj88y z>dngy7xwnGS%Bi#ju1mLky3RG5#BWqN~Dw?Pa zX5vmv=R_=To7dll$+@rU*)z>zPoMvv76{PM75WjS3KoM!Zgjz>SsCSd8Et_fJTn8D zGjrp7t-}xgC$zDxnFZcoBfZuoU?c0%kV>*S-~iEFhLh5QI0%ynvID!NkrAkC_;WDX zT-+ELm~4j3wci0xG%ucgrQRrvuY|qZc@Pq+VK%={db&Ww8|Mmo#eFWw|G)h0L-{Z)>%093B7Sv#saH1 z9gZ#jGX{sqWoq8J&+m#NUOXQhJgllk00g=J#~vS>Q4260Aq>3Qg8aHO$oxzYYr_8w z6GmX1-WG@wD9iN?SJ0MhdpKccudPRw-JJZ_$N#fyMG1eI#xHMXc^bHsZn)4}G6@y9 zf8PBMu#EXk2)(4QHI-FzX=0G>a&FpSGRNe&2JZcx7y{A#N{Y8y^``2b`6eN?-NBOL zQPNu01N%y-iGTBrE`#8s@>zl zQXU`9cOKX=$nha@$Q+8>SrA5VVmUoMypn9}KJiPPU@_g*DB|ogy_uMpNBl0lBkGgE zn|3A9I%s=e2dI-9C>a*S1e9z(`;u>3VOnMzWm+~h!4KzO^QAf{F3?@Al>`ZxHJEOP zV#HEbT{)Fi*?%l0u#Qd91Iko}vqvzMu!GFnRG8f!Sx~JYQ^Uz;v})~8uu!r$DsDTQ z`nh346G_@-#hKMTH2IDx;hwvaQ9JPxb{4+Arm2OAfu#Bj&_3e=Pji-@RGrKvVf{a~pz&W!4;2dyHrKs1X&mnB} z=Qb*jHu2>t7OzVSidEQn?-dOLPFxI25}9>prcm9v@NCug5_;X-3}WPOLwv^u+rI5S zyUo84g-X0EE&?bLc`}>Ea}Poplg3hi!3J zb&PG9ZADKmw)_CXN?cMM!oh>@=VS^SOQ20;l{_W1rJ=;rWfFM)nh>wj3O#JtZa#NnQ z`N);&`g)r*u|Jw(!hP&oSSZOB+;!_F4K5NYssx7i(l~_k&DO36bi44Of>smqXnFM} z9H#AN;umi8zO&%;t&kqY7bUbHe|^;E&y1Pet8uUUxBHKY53HjBo&ApAl64-cxVZGW zl_hd=HZE4u7XUusc8>JDdxa|rTpzWNR^?-xlRc0LL6~G0?l2D+tLwxi$GtlAXwPph zMdSttze^8JJQ{AX9L+HG$E`kJD!xo{G$}BRT33$ivf-R#j+2mm0*GVmxx6DImL(O9hvY&Hv zW2zjYQ1vr1C!K=x zoA~umO|n2nIUy0^aKCSa?x^8~@%tK(xSCV=8%{#P=*R4Y?3OFjEE>!qTQdzF2HuhA zxpd7*_ih<9wOf0Xz9~a_unplZ`wF3o7zp?#L@WA7DNSaue@whOa|#c02~W2prQ_yt z^h7^*1LENCvA`~RC&{}%Dm7A+i@;$mI9}Fl__E4)kpoua8vC*}%w_bB^To%*$)>4C zQ(AH5!Zz0BHq~6z-qQBQosYYCi#W#U5DT?l-=~VI#&j5l;=0_m?QZoGOSbDvK$MOY z4ivlu1i-gE`jk(%J*j12MV-UyfMBt!t4j{j{>0UBCV}T}^{%I6b{pA+8HgL*+Rs}# z)2r%}q46!LR3b-bB*0+oUT&5zL!v}DX|X_~%mm}n6JfZs_x{nTt0+$<{a@(usq$Nj zt@ls+rG=?d2YK?!aM?_s*ph1xkGE1Ha_*1KB@ZhCj*@2l6nVE+CV5cmm zgh>RX=d_bY$pr4B$BRPjbCV68lcl;F6{yJ92iFOe#VrKLz7L*YG~}E!@NVq+(NtJZ z_JCC+bDRb{3zSLpggl#zDPF6tq*jt+{CE`G{B%?85cmO;iL8C!_rWUoRW=$&nEA;m@At-GS8lSx#iJcxzuC8 zYE?gVp+LItt?WLS0N^B)8FG9c@o;cb?L;p-C>ZokwjY^2`0O$(v*Om7BHh(?^ShFG z4@II_hy)=+@W_2c!LV48S+kG-X5FpbjP<9*aV5xV-txOLj=`)L!|I@U3$i zcKpTW?TOoA?8@^+Xq76ohN=lqQKrv&v;m%A!S&*)-4lP=3`LpuTKf9REA@&xuQ>^N z$7c7-IE5AgkUUJ@>?hd9xf| zP6Z#ns>5;E^TbXS(Z=h(0AjrG@EOMw0dtju)x2q%2VNtV^}wAv2pb>eEr010`01U8DoXUv12vH8mRsMwVyLByP! zP(5HFK7C)DF+-a@)ToPjP=JVAlBB|ro~N7Oc70-Jq(u%~Vc0c`tAySZsRYBlZIvWL z7o}EX_?Az)^MgSB8KVX53fr`!1cbKpy851SfeI3H6GlehnO!G=jg8~sK1vJi8!}-A zjAUE!-cBr?IbpSo6OlZVsE1HZl#K`NZjrXK;v5{=cPbn)4;+sS^uJGh2?`j&(OP`0!HgDZ z5}KL-M?b&vS_%NeUO;AF=7Xw#%wAMAC_GlEYL2jel#wVs#mv@lbzw^FQ|2;#wa|Ce za?R+2z~TtXl(hsWY7J6ESiD8?9efes8PV+yL!DkFSNNN39lm6W8$r#vx_#{NPI_X| z<&x2*_{CJ>93@@r7gT8m)68|7^j)-k`E~HvcewWtluh_{F@Ur%G`+!GNsh@3NKZx< zH@*EQM7BR<&wC@bWI?IS*6$C7Eu~J<7BvL|V0FT%k@05Zc!w6{bO6L7D19!e)^br@ zo8f5eepEX0k@kTmVubup?h=9*a^xPHy*^gKq3^~t+D+x{unX<*?E2~FAY@}^mXDX1 zm>3hn;uC@Yk%`{OaHa6J47xGT&8<+Zujx2UYH+nJ`cR5;c6lSYUd4v7wzi!)Yo)Gs zmF=UI!35XnCW6&KWUK1sNrUUcu3U@2XSOlm&U>a%(z4Z~eCY%RMpn&8F;2a+++!pn z4DQ%B2tMoRUw5iUzluv5V<^;%oV9(mTKhZp_Y$h$uu4Y|nr;{Wn`2CsbG$grIW^HN z`=N74q_@b)5L@;4lwM5e`fnB>76z?^OxFL@>MC-!1&0Lr`3EP;gt3bO=UZHp+^x93 z-6bPu9i#nox=|o9(VhWjgP&8Z_-&jYwrQYtf-CjECHy5S8X=n}2Jk7rOsq?Y*V@!x zwL%N@`FV&vR!Q?O_l>}470zpXbf(CQh(#&It@59PS z)efmT&{)-PNst_IwyFeNfq&W$>&?`@Y(s9^4s# z4Lvujo66#zA+sEc%dX#Ozz#>*>8`$FmkeVsVM=vcu+Yw*%Hm?`YWwA~Py0}pmlO|h zcS!Jj4pUX;fPGmV%q~O)6q%?0xouRP11B*L=`hJfn~qi?bDPc-Gr~ubYBr7sCTHG4 zs2f^HaX^VOR;Q`J;J`oJFY8FDAsQWRPPS|(hi>x^^XxlLci2HLV3(rN?IwuGR+NHXJPE56$#sH{X23 zVyQLaA!FT|ZZ3C4PE4q!#wGnsDo|Bb(<#8q%9|iDohu?`$YC)ZmZOgGIX(y0{I--W zv&QcJG{LSo5iXEX*)Ojd8W*_1QXr3GLVsv?SA2OByVq21-vosTV2>ZAl{lHSpEV(% z#?GnViD6oQn$c74k;3)3yW_cV`#(~?dXHC-yH#ahrY@>P?cdFQ63uDYPKv^gD6KRq zo++OmGtn-Jl~pL~p&UamC_Se?BN;alNx!O!*gFZnw(yuU0q&F z%C=U8ztMcKJl~QkTx*@%DjrQ%T^ZFiz!M?t-sT#O-L)6Oqo022-{*kK$r);BXcmvL zC;{FVqQp^zNmakTOnqygUT6w*Y2Ga#md&Rd_|DdBa6DNVtAJi^9KQ&pRW_Cm(v6?8 zdEA0_U$1=0B-c+1nrPX7*CJXd!FD&p^!nw4@4BMcTC_0B(bsFXPZ>tpGJss4V7y#l{1f#&FS<2%?FHZ4@77=0(PhP55BIvO6I||pKQA; zV56%(qDy(-h3K-|d_=Slu97*4=ZOp1jq!IlOKW?e?>IXsp>-GY;CIB_RWi+XZ{}ymDlmZJLjr!Hx~@|e=mOV zyPpWrhsr|H6Ph=gT<6T?| zLF7h1qN=pN%-s|CU@L{!MT?KRzPR<+@+BK!Pgw~W34-pP(?A=#!lCCKA+d6dGApVx zRZ^IUYo@iK#&RraM3h^VcjolRO{$2SM<<=RnRc8`73+i+J9s=Q5dRFVcUi2< zx7Uv2h$A2yqNZ$VYGN$WH!JO+y!MyghM?U)Q>yt2R(1ah8SV7xcleDn3Y4xa%?%&cAF@Cf( zyf0!t_u~31yt66zRt~p2or})bgkga%Fy(m+=~8dnBSp>%r!r!vF_pK=M%2h#Tv8fZ zqp;4Kg{Ic4vlA7$3ezYJlNlRp)NI>b>8Cr0JwvFs5%&WArIhnRbW zot;~TDOGFC=Bl-<$Kq7TqRnW43R$2SpqUCQQ#FYwRB6QfbntQ;SC!nR!YObtv&XsV zu<!a!s7A?e-z4q6uCAx= z={j@O8u76FgP{lhQv?h-rL+= z5b*|1>|OvR&LVWaqLYp`OQHE6S!hB=X`ROlo%D0=+@Pktawp=S7S%JuTEkG z=M~APdJ>)lUj!C!&E_&i7Twi;(07=TSP49Q%&tvn z*dZGuZl3wXB@P zoPs6QwoV%dP2pK&FEn6fe+Lp+6uLyEaH z!dt1%`Jho^+qjt0M=LU!ZHr`1j9=t%uwAwiXl5ochzT?3Rk&`>V%jr=?atwBLy*H| zpB>j7hp(t!{!>yLR8cm+hc2p?_o7A69U8ir%ht?PAM0rSQ-o1&hKkeEZoQB?K$gj^ zuRntts=xkte=2wBI#+;+-YlD>t^|~6=bo)y&!MZl#c1bN*fGB|NN24U*!by<>}Vc& z#HV*g-)NrrE4;T5BP4#B{6$O5u%?K^)_CSMZSeJ7_Z{$~NLEfmVvvWwvr$9|XKu+q z^dR%`FWc8p#h@|h!0?U+flis7(>UF&msm(orkYoKbr8^qunSaXh}8O9L`PLdUiAvoVUzXk$Br-2-mE|PzEXpXaTFHT;T8?9~x7|lpt@mZ^%$2D6f zh6IVl4&%DoB^`IeNZV#GY8<3VJ0SKd&*nRijDs|!ZL+f z&N8XA`fS9U%0f;Gz$sgfNpX5WsotIcNKYbv%u6gpten1hZy>YRX!K1;FawWZqSh04 z5J_1<)tZ)Hf(&Qb7&8go&nrcSeWAtMLm!jI!%z5dcF{`&4|6f`gZ@KiLLhVD6Ub-`C1Qb*TT4p z5)9GzvW}ZOtGU2M%L_W}SjcT#Ooa21JnK(7F_bGCDYOM|>1cPMDdyR8y-`vIokxX3 z;7Ys23>WWuBm_ke1RBc87+Ytm0X-vKL!%}f(2(6rO#v1Z8%8BQ9O_hgK~ia{8axud z7W%Mpzm$%T-|7neKLaV>4I?q&w2$vHa<6ntsFe0fg+51smWwrT;y?@zqZGa2jT&sie(X&0)FT3Pi@6%~dQN~1Bf0el|U)v3J9 zZ!jz6o4#`&df~qcOj35{B&5V=m+-uvmu~-SIv;kH9&HWQgLqi&sB6XWQ?I0OTGkOio z3>P~?fueFwkSS{I>m_Bi!s7)R-oG1JA1Bd9BX&V*72wArmtPQ2+T_Ue%aeXV^u*rt z>uky=*U0R72{~RKu<7ie`o7Po<*~Hj3&*Ft5-~tIrsBTp)#n=ci@8Tj&0k}!f#gpI zwl(2J^)BRN`QwXFAB~C;%=8D~Zy@oL;-vS`bOCI^^@PHbG+2yQQYDdGu&Ht@8*z6e z#^4NzfQI6fj`MK+1XZ0xNv;lLM!W%2)=IOX1Ci7ib`kc(>Pjw;cO}n)#-Fq3X@fvcWF>;f%-koFH*sCSq6GKt zQqHcp5AauzlGyC9goGR@$Cu;ftJT5SGIO`Hy=zCzNK;?m8eMtrKOZ@$x(pT~z-UZ&>qvz;Ftu+BWc2(xBq-+M^$?JhWELkQ_ zjgp%1gAF8dZDVc_I)G18Lq+b^C{bq= zJ63vaOLR}+icGJc^k&I6N4f{Y?>MG!d#Oq`A}D|;r=?l3Hj{_3V_)S~7u3K`G55Ps zm!sgDX;lo#zNsKSUHwK=@nWn41G~}ExzFd_wMO6CSQ22^f9TbYBz`50-8clAcnI0d z_3zh22y?9lw3szre7bcJrSKX$7jI71Z#FW`{-L>2m>T+oe?or;B_RRENZtn7`0%O+IVBG)&+O9+;l$bZJ}maIk!+-r9;8Y>b7Cez|ZnNQ~WM z=em_}(yN!*a^^fjLpx@?D^pxVDi1Q!hvs<#*%E_1oK4K8LSi=(aWY+nc`^-P^U{7{ zZv=$!Fl3Y-X*bw?JW(hSCL#{!sq&Hio}s8ccY^adKH!XBxj(UCKaGUkYbvgsrM#-# z7NdX0kaO4Z)qyYrch?&zHdM2=fSHOPhnwS27ycc;L=Tf7mr31`vr`vJ{sN+bDUNPT zd)evL3N)f=ICF37e+KP&xpntm=Ri*#`c50MB8A3UeD)ieMyPpwE;#V7P+82l|Au~FL60)P;)`vxJIXiOmtpg|e88xjDr^k1qSc+t#U}c-t_1U5z zc6c+Ubv3_VD6TP{Ef#UAtFX zL7(oRk;?KwVVwNYltyK>*Vzxen)xZITCqQ`9lA*gF4s z1kk|@W)T6S>O67H=IAjwluC8n^ilM)$rBj4fIc04|h7l z#hze1OdE0K(st|l*#nkltjIS?j##(fPC`!m7;)B33fg0j)U6d3e8Kv9AV@tz&RK2M z-@`K8cGe;z5w9moM?zs0b+4HoD zBRtfZ$QDQC7-@#7)!iJ&onk(n!gTOT_oZ*C#L;vaa%{I%Q-(Ad&3gGaIT|3YgTwmo z>^e?GZE=JXEO?dQ~b}0c?-L5PcrFs&MY@V zu<$e4=p*yuzcu{+YOx`0x|YRh=6$eRKK)e&rqBX0QiLylMZI&kX30T4oSalXK9JHf zM(ms)z~z!RmENG{(QAWR%p%So0`oc|tP9!|FN+`wEts#y006=4LHYVi>eOHXVr8AZ zOcyF))A`oiv{w>kqajwrZ0Xm7IEYCG_bvC!`5W|C(Aq13ZQWqd58X$fJq>RlM;zXA z7tX3>ZRP_S)e=--cj4J>Q@D$PmpNp7a-?Nrz;nMu@Ppo5I4=uA<&{jM=qlUQadAz? zWXEuU2h5VLlTBfE-uW-Z(LfO?HibZVP>_7nPKihl&|a_a4TxSN zt|mve(C()+2;}ESzkn#-8YUrA@A_t*E+H(C%@^@PoKzh%q?=-9a928O9Y(e)E3m*S zoQ2$NY^->Bo2CEkhmMdseI>(j*NO8S~wyoRfn>8uIo<8j5$h5c8{j;D$i#ehH?KNXb| zO3j@~End7u7&6grK<3+~3#V_pLbu)R9X;1tw`W5& zC3u3x=K|zD%enSVgv&u+)c6>m>sD^YUFO#0`D5G%wZ|3P=;zw0GPeg;2Z8IJd)Mm? z(MVnmbdI2A_gi~SHfV&>Vg#Z%t(!#*6tCF(sv3;E2tVFmSe_@YQ1#N`Bvk$w0Vjp% zJq<|tAZZt>e>cS2I&6sk9h0VC(Sd0|5mjkS;!i>=zy7pj$jQ+{m)OoL`RWm`AamT* zv|H>7g?Y5eJcVbCLwhnCsdkUbdJ1p$o7;mdKW<*`Mh-|)Uo<_DvW|G4h}5^#{sW%d zdhT<2(DN}Op+hm@7Qvg81YArw55>XfH~ z5f(C`+e7h#JggS)6@N1g=bm}UzGFd^lltz-2#_N5w?vR%*q4??OU)F{tQ@<*l-BU) z`-?R29(f%Z#Td8=A%jwc; z`>m%J`kr_NebdP7&v~w?_k~`r0TMP-Y~LmdLs7W&xdda%W0^fBLf5Yd zN3XI?6Z5I0dbd9wVy1Kf)9WP9N9v9bGiHM=>BA5;iis(Zjj;i9%v4ItT6^|hPy9%} zWpQH>eqn^|tbwl{-U?FkWUKZKcUQUa&cVC%K%weGk68GCesZj)==EWUMopu$#OQLa zm-|xl0@SLjFWg@!9wjpjo`agR00T`K7dF4eAL!J36^rkCPQknskZ&F7G; zo3&jx0E;|`zPcqgh+R$1N)!Ti1sZk95)xY~WF;g#IvoR%aFc}BiwxwFl7xlOc=k^37W^LR-!i-DuC`ps+8@j%hV{S6Aekx(Dd};2=_ljp=Zd z>OI#z0)zBu6gtavqoxuw$Bao+oeMQbWPg!hix0&ION>V4sa{$bE~^-)NMgo@#0W%r zUyByLgI$}YE1e7V8(Yf}bA|?OhY)j@(B|h)LH2&RqRRL?JDFV(GnZqc_%;$a&(lYF zDWX12<4VnAeg(|;*l$m@xZ(guFujn(3zRQJTg1^%CP?dCNK^OErv4D_h3RPl9Dzru z9ym;d=wZCEXWyQ}eydwywc|&fdsF!{h_O)D2vdX#vZf0&p$GQc$0wmqNY<*(#P0BK zrJnV@QNyGjKsxd82vPXLC?mxTy;AXynD}CATkLzjj}M%_GoMdZ#~6il4VK77FkvRU zjCJF~n25N?yg0*Wy+0!zdQ4~R{7aZ4ZHgKA$?O)aXE1G&LY0 z4AO#oOV7+yLP#o%3Z6Q;iG$d!6=cgavbe1k(7KZ@%VT)dUQO!*)EX0($mx!Y%cj%C zjz8KQ+-(WB$42s>O>$IpqIl0n*}}dNP?V(f!*AJ-HN7n9;-X7(_=t zLmUFe+(7S`gKEkN3L4=PMCKy}mC`E&4xYQ5Qz9bCD|`aKkq~NE)uU3ZHGD*kO^y9A zlMs~}Mt&eQ%=E_uKU;Ru{5$ewP0ajMYsor%#`0U>V zGxeVq@Gb}4`a;3Cxk*f5K(T^+R5^fBI5yp-CH=_D^W=O>b4#fQ?^vri9eq>7=?E{0 zRC}}#BCD#ruA!E6D4Sa#%UCc^jMta{v>IH^Bjjml8ECWL`w%#SL_H-cd^8FB8 zQe{Hs$M=Cibpj&U$1?A>On>&7^m`9q_M|4t(u6&Zy?s?T(8n2RLc3UK;Y;mK4L7!*&T$y#A#*a|&7=LrWdeQKSfbip|!l%R}&i67u z_#Ls8rq@_3jg8=}_*hmc!cWi?UfY?e#W64EW*A>AO|T>{eGIq;0X|NDNJ>zc3UoPG9Qd!4oRUBx*=8vG6WVCQ;Z zfTF&^Ub(hat>*kefZNB>afADLhuR_TeyQEsg|BX@w1U{%acb)#8-fT5eGup9-8YF-LlruPxtjQPp3Va35c+#pmu1w|+r;16{>?4>^n z+p%B0w|~gNHRxv|95JN_fbp{k)SljRv7rWCb|s-PZuOL8atfT}sb?Irz`xGiBX6me zlP0D%&=AB#0Fgw@Lqnl?6$bhzJal+nq&s>4jo z`OHVrV5)T$)q{GM`1I# z`h8^YaM)(2AFK`@>#nq9(`R?jAH1q3j>-ycN{g++T>f-1c*eyCr@00r2bu3ZLO}Qm z9c}RApLz?_&j(r;HfQ~8j%MjL;-GC|^c)YK6?6=!&5 zzfZhuiY_Q3^iLdFFn{7-@3wO0KJtt^dAK(}BXSINTvQ5Hw4kU#*z?xE^TO(O!p5r| zH_@0s+1;j1)gWk6m0R}TA!u1O+Up$T!aCStPW~JPk&@IA3j%($T^-Ns$rxaPr&zXU zeCRu{zI=i8@|XaMA|xz$2dy|bwqh4lnsX^!&ae=>cb5M|RoTpwT|cpLap%_pY3b6? zQHA|I?esT=p=Zy*p2u!n4hD8I3_m_@?$5-&e!Gu016$Pr|gCg zw6zrahcN9;m1pz)-S_p@C%D_FFOti+@1>(_<Ojl65?^}&V(|dzHhlR{}ShR!MJpt&YnTGtM~u+Xz)q_S~^3N@lc|&)(T!x@QWF*Lkx(*{j{y1mo9X z5w(;x&c2?s2gU`Hn&mDqR6!{1yV&>MkTk8l!Q3o5m?=-2on17}umJ-$NdD)w+ ztF5PYav5T&_rT!Zx_c#45kf#USeJz5akcVTYaHHdJTUlSAcbq#Q0Ca4E}UV(kSsj` z75)&gq;pWpS>#G>ZC+YZ#5x0=l8+zYyBMs%DbuXD9FWZCcQinhWEE=ZvygP>;f zRZ~wb)(h#5Ki)t(AT%yryVrBb_kF-0AEvczqY)9ywUM7QFoeFPe#^Vc&Z*H^Ihq;&_rI`YL zhSE0Xllv49I%ZMnqNLG8YMQ&pzh}}TJ2Y?feFr9m1&!t>WDTky+Ux^AG`o?YAKUw< zpjrW#a+wXYJI?MYC015~+UY4XImQeWhut-}l(x(zPM82-3>EgOAoqIg9lfdCYqH2n zr+k-zirM~}ycc8d_O4*plBdYGBJp6B@20s@w#uK5`k_RM6~`D1A?Z~i-iN*Wj16{o z#jW?s#4l>fSy+F$<2r>j2TF6+pZ3iw>yf02|trIUV zHQr6wd#f_x_A5u|+FT6HSoqmoPLh=XKkL^n+yo;*b!EzVdpi zk3afObCsq%;ot)!QZx-=R~(fOLg6)?K_DMt_m=7fdNxaI>8U65;tX1RaIrc@A(aG*n;PW>$%;)BBDc`kEhj zrz0d^n%w~QhFV}t3Rd~jWsw%QXm0gj5PpXE6xUCxRPmS+xdLW zcU8W&ThRqhpC5i%zNds>0s;Ulh};i~1Yy8L&65N&)8}uY?T`B{IC37lZ=t%x6_Omx zLIGJh7Rr4^?*!w*acB1vA8=<~KGh*tKi|QqA|De9AwaFDn~(nje;XE%2l=p(lA zEX?O>5O;Io+dmhV2TbrW=ccySpHuh-w+T*OQ1vBH=Nv~)FB@niGF<{W~nF7S{zIJrn};RjUOV*gT;h z&pS?S_y5BWfw#tYQm^cc#J1o7yzEG-J=iy{41Ti>|DXM+nWzyvYmH9{L00ta6s#v` zVA;tn_2)mIOp^(TX2*2UdzZQMJ(a{})B#fW*V1+eK`y-RNB|@xe1*A5U$hkaj}>)m ziB4~ygmpP6U}~0yVWb#PIxqyo*B9tbeKaantm|RL4ZzwJus!beo%ac_q@zPH{(sG{ zFg_?9Q4U{q`>blnV&l|scWRk5P%06(>q_HvkQcHsAw2iNhtR-}$v7?BB-sj*0dVSl-&USaAg@!6vhGa*o0ho9MpilsV zU&jc=K_*;&B`4Y?IHFxco%=>f;3E)#67hMLHV#E84)V50Z?U|UWA&U7mmxXoNJrnt z{6hwsnYW^ij)kiGL~^W!E@1?DGNN0zECgojZYcwiXz}ro@bNp`;|Um}C{O$P6FkFN zEh(({>>ijG+M-l>{@wq;!lD&IFJ}Lbzf^&al;S091)gul_1_i{zriB(zj91MCidJH z&wbb0WN4&S-L0xdxVU*u0Ngv?$$HiVhE7S6>V?qDESN<$JZWw- z?shU{*NCMwUfq=fUZVr>SnQ>~`vPcVm?M=4kYr>6(4Ah)jfHh(4rHP#-ItVHzoq!d z#`bZaoKC-tr)$B-y$!wI_}JS046IzU#3Fc?Xi0H86_*Su{id;$CA33~^a z=9{&CIL<{`YFb|XnQw+g?Ek!h2*crT;r&yU751OzgB5;`+tYfpg8>jmLwAXgL)}lbZp`gOdbVa~qf;!WRz22AZS6*fp!cQ4eS>Vu}3^q?2 z#k1?j!n#*a=Uy|d{Q31Oau68;d>NUz@%j7U!in8bE0*EtzF7Lya9nGQx#6SpLzv65 zkn8Wj2|fTkn%rLC;G8TJ5Cd{KcJOYP!$oOzjEhM*Zp&^j1_5tv8*aLZDXxi1_m*H^ zOlJp5bI1@)24Vc-;>t8#+|qqa^y}UCS>4KUDA2?ERjlqFRJ|G(O(JDUkDDTe+c7n5 zMZ3GmY`s7bpw&_7hD6=*4^NBvo1x*Q&hO~XOmuY8m^C`wI% zh90clq@uG!Hadzz>E883SM5s^z1)M@^F?c3ry(QVXP^``w^?7a)F z`@V1p4CW~>!*y(YF%aMZYs_|l&vG2hdI=iu_T>HTWPbaQj=&)As@}TUo6Xmp8gI`J zHo^(G71=wj4dBCE!dRu5&2LOUkpq=V&gVStxFApX>Enoe3S&j{u{SLqQXtgteWle=uTT4>Rf-;c zk>9j=_iTZ{iH1n%G{GTuqu}r}{J~fjIqvV(r=nQB+&kBO;>N0?Mt&c;DZfkxw*&wM&26`WzWkSZ)_Eo)T}qWI7C<);>J72#q$nkR-00OE2WQW_T=u zhylS^D)9uik);DQyF`;25D5F~Sb{U+-AxX4VDV}|sQ*?Q_%KZt!2ewR zAQ;^v@UX2@e(+j}fDj)~n(S0GK93fe+Zk?M_x1e3p<{=og5$x^9JYz79EYTT-_w>^ zzVP#8P`c_1zh}pxBO#F=Pz2?q{kOg>VXNVM?G$l$z|XaZpB0Zp>M4m3x7{M__C}Mx zjb<7c_E0_LuQ_pl6iQ#R&!6bFuGbu1%EHk`-)O&q-Kr1$gG|B{1@bos-e3W`r(eg{ zXZI=58_dB5i(QD7(g(__cf|n*O%b4cNWG5G=3#$`gz#I8TQJ{-W(Z*W%Ll`C&!vqY z@G{(X*n0GeHWl^h2V|relo&8WkOz3-=n(c)LmwyuSVEb9l!gSko6lZb|2F ztBF>kD*I94&}GJWP_wn_SCh+nkIeIwB9k(^s(iM&kn05=jD|p>~xGD9nAaabb ztHJv~+t&3c`f?GGBT2V!zplIAw|g9z;=`*cmug+3TqH<7hD7n4{@RZ+$j9og{TG2^ z<}*td1kUzT=e=_)cw!$y1mE7|{Dp*MAAydlJ#tHCJb(av8H^MNOQ%8i1=9@xpt^p0 zlo>TxUD2s-{07lRPhNPC1%YXI7v}c+bgGY2%knc}N5nv4)Ig{NtppM+6burl4J>93 zV^xtT3{z0eW-k1+t6g!tR>WI+SO05tPt1>5LG`z}3L#;SZ(CwbAET~=k6!LVgXc}Z z(mgEC)%RcL7DXZ>WB5mAK7-E|*^sNRnSRI@$FBd^!jw}fn<*(80D&>wAN83q^NRXL zD<-^}K7MIVum2#gZ#;#ve|z*_}TF@Wj-Slxgg_@*n-g?KU$3? zdXN2a(XrC*He1O4!CXhmdBw`n*m_1oBe9ti?IasifPJ1IfdQZ0XkPGONS6&j|0&%H z0e-_ME6<6*q=nV!t-X+zezIFc6>K61r5`8;`iB-FzU!vD&7~IkbC%Q z%9hoQ6u0)JL+4@5Q8)S$KCrACH3zCdB#tB+LkIGI{p?|kgrCN{`o`JoWUSqijITm! zG8un?F&7Echx3!Wj}g<)kJLvr&VDWvS8bLLxkH`e6W~wl763pfKU9~R@LH^NNA9Cu zLtS~Ly+o#vvd{8O(7%xNAqd7*WMlizEU?K|QKmPIT0U(!O-%?;64t*q$xziRa&|r= zVl^UVWJMJIwjv~TuqfA(BS z{ZxN%T#}FHD^DL_$wkV3hV`5Bg^=sFnpeUT2Fw)}>>LMKDZg{aXxyLgeEd8)Lzx5d ze1WKObWm^(6w}+C@H?fY$2xXHw{6pA(|S_6)zM719129{V;qWj9As|--GkvB5(HKE z4aeq1T@b!<#oaal!2Uo$P-^)s^h5V1V=ye@GQvQ(OS|j8q|d^H%j2~W^(MJRcpvkUhr+bYm_QMtT>ATVw4I{#2*qHBkr-r zCK|)z6Hz5^tF`V|>@xSq`rK~O!F=)dOYV9MjmCd$%kvV1;aknvOwO)ZL`5udpo+sd zKGhMkcSlFivipu_yG9u`OMJR6qDc_>^5!ZiBJ&GoO|%)u--TN+v9RF0y6tJnRX#jeOx69TjB+J%XBXvd!PgtDLj z2@!nIQm~nB@qKSV$fTr?PcPMhuD99HhP>7{Gl1LkRJUjrR8&MX2CtAN$`5Fa;55e3 z#c6lbzOCdA)BI6r95QeqFMQn9efzJ2dHhW36Q_tj;om@8u4TxP*lD>`_3B&{857mw znw_JGp__{eL3tenbd_eRRpmOK4LUz8h3@?B?BdsbYh>QMX`H`Q(Qf_s{zh)C+4Rsq zWC{(tGE*=+pli`5GSD!>rkj#pXpVZs^|yig3mm7MqOY1JoAfURgM8z}WZS)~?++XA zei*5Jcp9@g;TVf{m|l$Rl^R?Y&7!;PS_uRO(Gf|IjEREj(8xc&@M{>|%W_&CZB0iS%YThA_;FKz)7Z zMa0MbCQ~ZS@PFD_H(TxJscp4Q#*JLOqG{c3#dvW2eBU3h$rSxre7c0M`|p$u8orfF z>cTh&uV65eW5lnOy2z_kLaJbnJ7Qwq$n^KK2iSUrl z#&u?a82CV2GZdmQJU<_eP%z4{ke&S@k}C;6S;AMaZF}dxov(Ov;AOSuL#69uWaUmd zlC^?f6)7TO1DjYJeU)mEjF2TWaETLj`+JE9UOh(KkO_V>oY5WI2c z8td@HKJv=j$`U*Il4gf!k`Px*hJShd!m(qIy>YL8&bA2~pud_nfUO#PR6uFbZsVyU zHbc{ufe#pNH%g#Ca>Y00|KcqvHbd{1x9iLCZt@w=fX~8uv`t4H{vB|w!^*rA&EIs zcyKM^ZOHJ6pXcri_3W7NO9qC6K{nh&j-sK#2(|&Q$)ZNiCx?y22B{}oA98R`mes4T zScj<#MKZ1AZNDnzDrO@=*&4Rl{a;_G5Qz)oYJo8iQ;nzuOU)vXK1ENbL@*cP{az|c z_?O3q9r{mPVb3ceWevSQ_9@jB8U~pC(K?47uFD_UO_34EuA8BUbeq`#jyOAgF(gqC zi@?$u=Npoy(jNDn)?O05rq{FwN2$k}7lQ{BFECrn%DTs<$l@qEt>AI&sG@E(vDrwPx4w6^s0`;q*R zKybD`__M%=L~U?zSQv(mXg8{|K+LYTYxCht;aeOhm;cVo(Tj!s%bwF*h>dwkn9&ZG zYgR2=5C_7m^1rzN90jM6YKq)B9HztokFD`Se9E9f#Hl=IPF;!EvtbrumM+ZW-H34f z)VEVcboZKyDu16&KsswCNuT}&+giQ z>whm)Q534G0^wl|bqkbnJ=2Oy#W9-G2P`4HUYrtBON_b2X6q|nxhLwbIZ5>nhnlOKf0ULgwXX5(lqf0We z1T&-a^1$NsKv7gF1kUc+zk;r-`HS{&7iV-^ahkWDGK};QRVJ!Xr2TrEF+?l%g zZgB0s4QpGln!1*|Y}Gi|?7TR(4wpbefq)S_S7;y<5<`#&1``27vt(MH8T+lmjj+x6 zFvB&AoR3WqOb8-@%Inu0O@ga6Na`nlRdF^{+I-Z}rsU(T_lxQzwG@mb(fOly!P`VG zkh<2-#l%yp(Ili=x5YL~)LKXRDw)rG!tp1BLI0~iErUcIW07%$@Y@ve4hqmzW@Nzsn9z4?V9>JNJ@8x5h`xxcPhCu@Eo{*c*)& zgfZu%Tw~2)LC(wT*}kfBe|nT_nrQ`Fze*{7zg&MH-PtyHI?=YUPQKW0T{X?)&PGgs z^TFsJaGXtp)jJQp%x@JWz2wAm^1^;ObWv*5An|3U0cOG!M824i7C6O7vHBcp^oy*3 z@LznJS2Mg}iiwG^`RI*WmArvLW7)ejR9j_M{~=$$rQXMk`06_~5SrpywOR3O3&X1! zSY2@cF@LylF;HkA^KM0UM_x9@Ps#S{h~(H*MWXGZ;@7cAdhydGN~EAHm+##BFj34< za}i%I?`n&#-L$&YuuZ*zjV*x;pWO(?hA(Z^+a(-+Y`)0oR_&{99aFl~^ADLxI*w<) z+BhL-FrBGW6U35CMd10@(}({IIz7FF{SzSK&L`Sc&D_l zm?>ZhnBww8HFRWgBQejf2G6hi%YbMjR%i@nn{HIiJMz^Z-&7Z=ip}J%^8vo@39JS3 z*NJva^^epoHM?%8z8?Ok5GjwY$<~7TfYIGbalcO2<>wse7e z#QvLldajuDk2;h&Cf=s?Oi3(v==Dv`<(Ad{u&KctjPDW=5{zG$h09ZWC-lZq_J6_ZuJjQ?I(17Ci$*2@ z(dMQ587nl?0k^JqFQp8Ju9bXSW|y=Wu4;MApt}nkMVHb#J)pgT`na#IG_Po`k=I@3 zQr4RzAPn3iPK)LxC7PbQd*>!4S_A6p+;AU}neK}zJeg3?7FOr#il5;=5Fi2lfZUN# z!w!SH-eh0bkevv*d;SVFJpN8-IH9R^$=>fk`w;NBcD_SX?#)Q zlmd_{v1B}tM=#k=$m79LKsP!h57PBSpE_CHN7b$U-|B^5#$uOjgq;%JZz+!_=15sL z$c9DXudL+cSnbn|q#fV)^mzF?dxqRHc0quZ-m&OTOpCn!*s z+2Er$JF12?bRYs3M)5r}@{l0izLHModH+cfJ3T$=o22ZPI@tg-hN+ya)>)m)BI4nt zr)HGR&hFwN;N~q}U{IoR^w8rID`z`ZnGI$OQ}o?v8IW z!VNb3;bS_PXzE~}C2If?k^~_ds?InHV5NB>$`CVp$SRB=H8K>@8lV9RS<7q}(wse3 z=M566P^klGEiG&dlz9XicATv)KPB@p?hQPvywYxDTBNdzrtyI%b;}DSQ~jF250eFX z>Cn^d)vWgIv4*lBLh7d(w~}&u>0%JeS=hcdt#e3iT|zF`a?!AA9ji8(>l(| zIO#Q5Ql%acvrcdIRDb8Ret|~K|7%IFL(}Z;nwW)6t@q_RBD$aNP@uUTceQ%*mAIr}pKc&|AnCRaP{GEe*eY(I! z8YC~Nf=&UXJ>Ab@P1tUF>S6GJM|m~RE6cA(|jof+b)19k|EQ zCp){@UT(mLQ*jQ|uGXp}T(m)tkbB>$75;ouai)(gN$SX$-VB8xB_buFY*u}~Z&OmJ zhVYw&S&uL2pcnsZV0>4aljMh)6JC?orh(a7G{fAx6qIvTVF z_)?V5pUXBB^sL^utmd|v)0Q_&#(NHiQN>&Lc+11BOBzh_E zIAJ+5?O%|nNVe|)$TW91A^pqO!j}LXD%M)2WBHHAR6CXbft1`y{{DpY*cUMr%TP4r z$7_BXHXyNxsBf_^-k-{;mUr=5Lq;A)`D0kccnEF~Pf>B)+TY!=i_IGkfVFQn-*T-@ zlm*TLNB}U??6s-~TrhzE>F)q>NKs2A0F@b&&L0Z(%MBr2A;C$A3za6o_xt{JS_aKr z$Pn*xywz<7d5@U)h;P}D6wMd$IM=^p#X_B!W)2G7fqv3MAYihXw|;I*ayJMIg`T4) zF;k3x(U@;2XPfTx#9(#g8r9=dybN`;^)()t&c8Tkn@J)$>YVtE$YXg~)(kUdqAG7M zO9qm@yg{!rr}~>9s183CYa znt2+`anl>!$L<9x<6RY4`N_v-blCBQi+E2u!M_}wQ?o3M-v2<{i@SH@^KQ3i{UITU zM&I}DBB||me+NHtZtTbTp?|plm=m$=^>j1cuo7WpZx9Tp0V(Uw)yw0gJDi?fC*h)^ zYZwHpqd#iEH@(>Y{51C|DL4`h@b!tR#O~;>+BU{op3Hghi>hFXuXyk~{6AR>ucGFr z@H{9>Tgpnvy(y6J_RYqi$*xRv)?!|Dv|+QNcgwK-?O{$C^83A$jMdt+<^VvYLux#* zw@?AB@wsbmO{50>Uo#tp*9sA4Lp5{sfh{oPmz%*EW4IlL@krQn*6K#b&l``L|t9nUC<4XC3 z(y6?sM$HikdTP}MzGNM`aqAM=`4?r5Ut%WPt#&47F=QPGSSen9Td+1q_++1Eno$_c zU+zR?ua%*d+LR8AJ%igsNVgFmZNleZVppo)#cWNh2Hj;M`2rNe-wyT5)mbFvOs1Uv zl$Bc}lutThvDP}5@zLHKo3rYc1ENe`wM>Ni2(xULyQIeFxf-|oXK0%EwChd%K9Q1* zTjR_RY-LytAZ}CDBiT;`C9ilxO(ol>YvX$qSERsuD~*nbsgaDd`gRVjvjgpG`8&AY zf9zCGUa{>SwxB)?s+H~un&zcYLlMj(iVTYV4C;!orVi-|3K1v?3MC5k^Mpc( zDBo2|2=lvp^CoCFT_sO1T>TRR4dVt4saBcYF!5gD6B|S^mU2?o$kbvH?phNcS?OZD zB(|!G`X;2Z>-r)jib$)B?f&w|S9UhDv@Y(^Gb-9%E!MQ>hpeMsC4@1sFXgk@q={oQKs_c%Fp zlG-0iXpHQ06)h%)=bPR+LFjIqRdtS&#p@Q-Z;zZEe^zJ?ZT%u+BtS?jD9yy2cVCMZ zk}FrD@b5i|U}xnkh_AKO-d@$r1d?P58;(PlDvlR#u0z)P9i^YzFnu8?MHS+z3Lq}h zrkLl2=9A0#7Y~0y5+bUo5?RIm#jI=A7j`3drw?09Y!}|uTSJ`uOCFnt?6Sx#GV!V@ zYI{za9|x<4EaMM!6l=o8ex*HJ6WOBrqBYoAiL}Hn?WFVm8M`q?Ou)d$$Ds4ron2mN zaMluroob52mxPQ^yE9JM^`T7dVc<&x@d4}qYsNS;I$cTbkQBfet@{{2@V)_27!dJ} zE_&bCQeB&Qne%5g0TDZ&nar&}e7OF;D&WHFYE6#%vnCHL+;v1*ZfC?cA7i4z9Y29* zD)Sn?K%39uA6b#3DaBJ~i4JVqGA-K()vDzD7^xBbP7y(X^q#qCwGjZbEB!M6xW|h| zmIAMOn#>?Q)Kg&d?S%kBWlGo8w6rNs8DZ2lhiJ}G2AdVl;Bh?|I|?Rk)zGdPqCUCz zgrr?uu06A& zPC#Q`>0I`pcv6EmZP}epYJ-n701fFbNNPiKvyfKsq?Ok3Rz^(%j8J$;)A_~A)ubVs z^JoD}_nAcWj6(10n>~CAF{&Yd(Iyjh7gKZ^E`YQQ*X-L?K$++-)EJQI$qS@_dyT{j z%=zXIK4tFhmQ?c#z5nhV6T7ps1Vy{|gB}dFQH~xZDNGUh%W;P(NFvoxE9fVgtNCnH zm>RYvg&?IYX~PICDtlERaU64_)z#`o7{XeTQ{Q5kOo+E!Yj8sk^LlXO47jVh-9CWL z@S+`m^LbjEB0My2VEu++X&6_`$ehQIqk;!z;F4AX{fTPEZK<{w>v?MyH8gaTKGB@C zKAXx-6km1qEd%gpUFupGzR~wqy>2V_n{LJIAu)CIwu{wRi;}5a#zjEO>SpLPd6oMS zs&bWQK{0>40?O0FN-)EaYu`oPMTr$!roh1ZMf@195X~Jg6d?2GYGB(`Q`bm9X{D_C zRqH|~kI;M9SD!vs{~lE0p280DoU=$l2jC={v zPII4q3&oHQE&ScF)@`m5fo4#}V^(ln;(Gk5qCQd8jJo_yG6$KQnWBM%5D;rW{q>Ij z@U-#Nh&s0C^|MPYRf}0eo_o#1XtJ>q3IPCV4?7bh#tA~;zf!&m8!8ljq+fkG@;-t0 z7{=+e`j(0I{E1gYG4>!R=v4zcMq%gHz?ycOqiX_6;fhcsC=>*d0Rm#;^4cH^8`adW z=M`17)K#?50d_!uFO*cE$;=~6=i_cX2l~xPy`y}`fv#d3Cug+SuYb}UjmrV!b+ z#E2YU3xK|SdghM|&+`HG-go}rpSrrnFvIJCz9siGHrUoh?P2p@;b}!WN>Y%5l6$$$NsdS1P!_4rkPyVYwGCAGSjQEIfv0wuN<{v zdRd4c6?A0<(iu2dn`sDy-f=rQ>H4sKd>gkNQF9;BUAz=`E zlkO^*q&rz`Y2m)**R%+i8E6&w&{FPlQhp-r+;kuy?D(mL$Axkf+rZCf!t1vI*#rt; z*YdoN@@Mx~L)w7f`sX=uRE)}ow$wL``lJ9vDWym%$;WYNUo8TPfP?OLfxinR??(7} zX}tVQOh$2q*QmH!Ek7S+*s3>h!7Mt*{;HbKdQCf|CR<>(p-JF#I-n;kxfFHNLB9xFWKwB9Ng!n1I zSidQJ!~-(;0e(Gi`^~KqR`|%tzDf7ohUOwZMi3ZE{p)*Mm0B*k6>I z-Ua)N4!xdLJJGIx?B*M6R6bieGiH=3TqfNNan2M4_-|iVbas9E-byFxO3m$mopC+|mgZ+E0gP%p;Y^XF$Jw_d7wF?N3{1($UKWDsI5jwibdfH7UyS02=zrgCR% zB;RSPkxg&0xqZWwmSI_E(!&qxcG5&BOH5QQWhvu%o%D4*iq``Y%3Xz{O4?O4e(a9M zrp-F&n51t6)P?$On1t%rJ)RJ+t+qdh^E_m?wU$!UC$4%oFTqq2 z^d{+u-Tu-jWVP?c?T^K?kMlOWGb?8-i~RBLg{34tDR$TZ#I~LBp=1OD;ips=ZmFJy z-=&Dc&&;fbgxEhM@ez*jw8#N&;)` zAea~U!j7}|;SF1H>awg)jV2NRR0PhbE|%F6EU`qAZ?8R`CJ z^l*jEVa>0}8(oX;dd`leb%*n*9sXWQvCFCY7=}hi(4>MU3d9Ih0wN^f6A9c$^6hDF z4DC%n`fpuZmBAmD6rIw_oo9l?^zuBu(!@3$7wZDj+vwkYGy{321HTn^u>PvIJy zJWW9pp=k8iqkN;SZu;x7mWZ7A4|SuLKGg3`E3+H*)qZI@Txl5Cbf94$X?GcAiD$&d zu#H;qWqm3s@twp;>N%n*fvqSFIje_g8qhIk=SIvCM00GAa@1D3h<_E0d*c=Aa5q@^ zIs#z{>429sy;%+B@b)-7SzAg=H2iWdv5B79OZ){JFsgRdE!uM%2L=#K+MTq$9c6hX z0NXuMXX5khoH|@mHP^0kw@yQ+)=L{Du&+4zfl9WwpS)h*=Kj(Tf`Tbd1jY-3N+7BF zh1;z(JC3c(VQ0h|lEEx9dX2$lXCu*=0?w|k$Gyl*QC0j+wYcrpZ|%n$E11D~^80f~ z0}qa}rVd#Rq^y@xwZx zdNS$PFK;H!`w$XHzEM=Y*N(euGy@Pi7C(f203YqWI({64CI36*8`@`0DxSwo$tSBb!PrPs!J8ip!5IA)$csA1x*>fJ=o{o=FbSg5=Cy%OEMc05+4fkijp6ej&l z;@{h_*EnaXZd15&?)_=_TASwB?Rdi09E$%t0xge=#L4rfU`%)^Bg49=H1YuX^g2W! z!}~2E&~*V>feSq)22{il#_@j5mLtQ_Crm^l7YQ4Ur^mj!iR*2O z8YGkNmb){b=dJ*FwtWAgyD#_#3mR~ha80gXLm4jefeQ#`kmLh*-IXuTa}?(yU!7&4 z2gLj(vugF9ftevSYf9d@D{e>8DFa>fi*~HLNJ(*&iK6ou0B1Xu!%yyDH{0HK${QTR z$?X-CenJ|Gum>8(C7q>4hNs$FxNU`_L{n~-XF;}Ixyd1O2_6`aSfUVliu?^hwM_nE zznnmV!uMZ%&0lwuUKpkxfyQlVCMw3m*1684>Vv(BabRu z>(B!&5h0cgK*VX`cO%~fkcctcXi7^|ia}Jgp(OM5p(-BAE6v|WzNEJiNzNltlxUfyd@d1c5fDWgunC0r%QegB|e{} zr3C{YNk)*LBKZ2?509z_YY<*A^pP(GOO7RV6H&?jZKC6DsNEhWIq!LKOIrNWVb82^ z)wi`bPT~TIL@yIo3b14%ij06ZxSP&-le*g)IJJ9I?1gPV87CEd#kcav}b3=(>`9TQnoU+Vbj@r-uQrdZpE z&w}rT7JH6%u)`&Y0!6}30Zu2HyrT)C%{qln2}UQGMe8^`6Rsh?OTSaKQ&0I}IS(O{ z-VIt-02qSHsTDa_SN6SedR~AZulzoj6ggJiybxx7L#~pmsH)0EuAm4&_>~YAAwcCH z1S)Ubql(I{v6<_%jMieN{l5P1Ecm!JXg`+=IKjYtZ2C!QCMc+%32} z0fIYAfMAQeySqEf*}UKR=XUR&+3BvX>aKQMT>04T00`RHSc~L+#WCXf5r#-D60)wr zIQg0}E>N}8y&aLifX;#cD5D98dJJL^cT@+JY~&>SCYG|8C}omPtAMSR)x%b~EW~@AK2=|_{OgJB0NJq!@eq8Omm)@7Pdfm7RIB7q;yM7x>a=AJ}Pkeg7KPwBduM zjYyFIOC4@EYH21Ng^GYh;&ObnzZfz*-+_Ux@d1?~6UCaHY=kf<*ST91lJ~=GqA`w( z<+R_QrW}qk7yzptFCEP(M3mq3xmcQnCevIb`qR8!OFL{SW>y?JAePG`WA<<;RclJH zfbe^V>1Y>Z1X(Nc&GxD?ZY>cQQ79lACg66YLmNQ@8uogKunH^+2+qs~0Inkk=VS98 zd4%j09eWTXqQU;p-^*w}L=#rQUELo@bV?oVY04$%bx0Avr7ne~1+gCKzzpY^T0qn5N~Gk6ow6($%z6 ziA+nDPl_Y=ood(^%3dxssZ^yA_~xu@;;e55X|92k3%S}8^^8fffn#5Kwdl%UCkXQ} zYw|wy%treAyo&f9kdR8vy}i4h=Tskc0SZFbB@$`6D3Q3oks%$2$k!o%?_#_ULt>y{ zpg_Wf{RoR$5{^q&Koxo`@6$Xg1QfCKXjuBjowa%MCKTok#9f3)QWrT(7FNO!jp<|- zz8dBD^G3%@OtQ5SWcUHD(^nrI*IzEGzjHa>oh**)YnZ;BO42APpQU5#WAB_@>^fTXIIi_fm4;YQidUphA4e#3F zp%ed_Y8xUR)UYp}1@#~z1PSb9N6Xvi@?Z5W<5bsQZ(ce-aJjjP0*n(tw+H}_a}So$ zsb+g(ldsWfhjSZx5=9wBM+f^c{nS|`p@58hU*E4Y>HhYb^FGbIEaY4wIkNZKX_lXO zzqVdIr~_1US4euas*a&}lP&4zbq7=}5gKaH4bRiHq@*lGjlY@gKW-G3UrM@M=Uw%! zWW-vt*VT}Wo~;~`95nzUFcQi6lqR!0haZHVe&%!_jgJPbMSTx`8nYbkHgMmoJoC4|N@Z^{t&9Rx zCU6@cr5s#UY$$-dTtW{hefF{tp#83Q;3x!Bp1QHU6JPI;>r6*xxz;6Ov)XfJ&S# z0k}a&JdTgNwyB)nGUr0&cglCWe&MHZk%1pqssZ&fGzoMBBvq(940uxpzTDK+BXt-; z=ScFSC57sX=k!xWLIsg!bFKmnPFP8OL5!7!4CyNV8`oPxqT*e~gWZs(^m`d7ZG&Pn8WtcTp~V#aJF0Nx9z8b12}s?0PHH z)0^gCTJ6}>$z;An(Rm?IHDqyD-in$9s{ zEZ@M&T3tSN53W;~$<_i}>6{zZ$=leixF#1XMEJllV@LTX`Ga zEk0|5L*5Z_Dp&8)dYGve7ukf__Y`w~P9aAz27M