Skip to content

Commit

Permalink
clean function signatures
Browse files Browse the repository at this point in the history
  • Loading branch information
David Jourdan committed Jul 20, 2024
1 parent 0b351fc commit b543f08
Show file tree
Hide file tree
Showing 6 changed files with 84 additions and 87 deletions.
20 changes: 8 additions & 12 deletions src/cli.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,10 +110,7 @@ int main(int argc, char* argv[])
std::vector<int> fixedIdx = findCenterFaceIndices(P, F);

// Save the input mesh DOFs
Eigen::VectorXd xTarget(V.size());
for(int i = 0; i < V.rows(); ++i)
for(int j = 0; j < 3; ++j)
xTarget(3 * i + j) = V(i, j);
Eigen::MatrixXd targetV = V;

// Simulation
std::cout << "**********\nSIMULATION\n**********\n";
Expand All @@ -127,13 +124,12 @@ int main(int argc, char* argv[])
optimTimer.stop();

// Display distance with target mesh
Eigen::MatrixXd VTarget = xTarget.reshaped<Eigen::RowMajor>(V.rows(), 3);
Eigen::VectorXd d = (V - VTarget).cwiseProduct((V - VTarget)).rowwise().sum();
Eigen::VectorXd d = (V - targetV).cwiseProduct((V - targetV)).rowwise().sum();
d = d.array().sqrt();
std::cout << "Avg distance = "
<< 100 * d.sum() / d.size() / (VTarget.colwise().maxCoeff() - VTarget.colwise().minCoeff()).norm() << "\n";
<< 100 * d.sum() / d.size() / (targetV.colwise().maxCoeff() - targetV.colwise().minCoeff()).norm() << "\n";
std::cout << "Max distance = "
<< 100 * d.lpNorm<Eigen::Infinity>() / (VTarget.colwise().maxCoeff() - VTarget.colwise().minCoeff()).norm()
<< 100 * d.lpNorm<Eigen::Infinity>() / (targetV.colwise().maxCoeff() - targetV.colwise().minCoeff()).norm()
<< "\n";

// Directions optimizatiom
Expand All @@ -143,17 +139,17 @@ int main(int argc, char* argv[])
// Define optimization function
auto adjointFunc = adjointFunction(geometry, F, MrInv, theta1, E1, lambda1, lambda2, deltaLambda, thickness);
// Optimize this energy function using SGN [Zehnder et al. 2021]
sparse_gauss_newton(geometry, V, xTarget, MrInv, theta1, theta2, adjointFunc, fixedIdx, n_iter, lim, wM, wL, E1,
sparse_gauss_newton(geometry, targetV, MrInv, theta1, theta2, adjointFunc, fixedIdx, n_iter, lim, wM, wL, E1,
lambda1, lambda2, deltaLambda, thickness);
optimTimer.stop();

// Display distance with target mesh
d = (V - VTarget).cwiseProduct((V - VTarget)).rowwise().sum();
d = (V - targetV).cwiseProduct((V - targetV)).rowwise().sum();
d = d.array().sqrt();
std::cout << "Avg distance = "
<< 100 * d.sum() / d.size() / (VTarget.colwise().maxCoeff() - VTarget.colwise().minCoeff()).norm() << "\n";
<< 100 * d.sum() / d.size() / (targetV.colwise().maxCoeff() - targetV.colwise().minCoeff()).norm() << "\n";
std::cout << "Max distance = "
<< 100 * d.lpNorm<Eigen::Infinity>() / (VTarget.colwise().maxCoeff() - VTarget.colwise().minCoeff()).norm()
<< 100 * d.lpNorm<Eigen::Infinity>() / (targetV.colwise().maxCoeff() - targetV.colwise().minCoeff()).norm()
<< "\n";

// Generate trajectories
Expand Down
4 changes: 2 additions & 2 deletions src/functions.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,8 @@ using namespace geometrycentral::surface;

TinyAD::ScalarFunction<3, double, Vertex> simulationFunction(IntrinsicGeometryInterface& geometry,
const FaceData<Eigen::Matrix2d>& MrInv,
FaceData<double>& theta1,
VertexData<double>& theta2,
const FaceData<double>& theta1,
const VertexData<double>& theta2,
double E1,
double lambda1,
double lambda2,
Expand Down
4 changes: 2 additions & 2 deletions src/functions.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,8 @@
TinyAD::ScalarFunction<3, double, geometrycentral::surface::VertexRangeF::Etype>
simulationFunction(geometrycentral::surface::IntrinsicGeometryInterface& geometry,
const geometrycentral::surface::FaceData<Eigen::Matrix2d>& MrInv,
geometrycentral::surface::FaceData<double>& theta1,
geometrycentral::surface::VertexData<double>& theta2,
const geometrycentral::surface::FaceData<double>& theta1,
const geometrycentral::surface::VertexData<double>& theta2,
double E1,
double lambda1,
double lambda2,
Expand Down
22 changes: 8 additions & 14 deletions src/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -240,10 +240,7 @@ int main(int argc, char* argv[])
std::vector<int> fixedIdx = findCenterFaceIndices(P, F);

// Save the input mesh DOFs
Eigen::VectorXd xTarget(V.size());
for(int i = 0; i < V.rows(); ++i)
for(int j = 0; j < 3; ++j)
xTarget(3 * i + j) = V(i, j);
Eigen::MatrixXd targetV = V;

// Init timers for GUI
Timer optimTimer("Directions optimization");
Expand Down Expand Up @@ -285,7 +282,7 @@ int main(int argc, char* argv[])
V *= scaleFactor;
P *= scaleFactor;
geometry.inputVertexPositions *= scaleFactor;
xTarget *= scaleFactor;
targetV *= scaleFactor;
geometry.refreshQuantities();
MrInv = precomputeSimData(mesh, P, F);
}
Expand All @@ -295,8 +292,7 @@ int main(int argc, char* argv[])
{
// Load input mesh to Polyscope (not displayed by default)
polyscope::removeAllStructures();
Eigen::MatrixXd VTarget = xTarget.reshaped<Eigen::RowMajor>(V.rows(), 3);
polyscope::registerSurfaceMesh("Input mesh", VTarget, F)->setEnabled(false);
polyscope::registerSurfaceMesh("Input mesh", targetV, F)->setEnabled(false);

// Run the optimization in a separate thread
thr = std::make_unique<std::thread>([&]() {
Expand Down Expand Up @@ -342,8 +338,7 @@ int main(int argc, char* argv[])
polyscope::view::flyToHomeView();

// Load input mesh to Polyscope (not displayed by default)
Eigen::MatrixXd VTarget = xTarget.reshaped<Eigen::RowMajor>(V.rows(), 3);
polyscope::registerSurfaceMesh("Input mesh", VTarget, F)->setEnabled(false);
polyscope::registerSurfaceMesh("Input mesh", targetV, F)->setEnabled(false);

// Run the optimization in a separate thread
thr = std::make_unique<std::thread>([&]() {
Expand All @@ -356,7 +351,7 @@ int main(int argc, char* argv[])
auto adjointFunc = adjointFunction(geometry, F, MrInv, theta1, E1, lambda1, lambda2, deltaLambda, thickness);

// Optimize this energy function using SGN [Zehnder et al. 2021]
sparse_gauss_newton(geometry, V, xTarget, MrInv, theta1, theta2, adjointFunc, fixedIdx, nIter, lim, wM, wL, E1,
sparse_gauss_newton(geometry, targetV, MrInv, theta1, theta2, adjointFunc, fixedIdx, nIter, lim, wM, wL, E1,
lambda1, lambda2, deltaLambda, thickness, [&](const Eigen::VectorXd& X) {
// convert X to V for visualizing individual iterations
for(int i = 0; i < V.rows(); ++i)
Expand Down Expand Up @@ -384,15 +379,14 @@ int main(int argc, char* argv[])
polyscope::getSurfaceMesh("Simulation")->addVertexScalarQuantity("theta2", theta2);

// Display distance with target mesh
Eigen::MatrixXd VTarget = xTarget.reshaped<Eigen::RowMajor>(V.rows(), 3);
Eigen::VectorXd d = (V - VTarget).cwiseProduct((V - VTarget)).rowwise().sum();
Eigen::VectorXd d = (V - targetV).cwiseProduct((V - targetV)).rowwise().sum();
d = d.array().sqrt();
std::cout << "Avg distance = "
<< 100 * d.sum() / d.size() / (VTarget.colwise().maxCoeff() - VTarget.colwise().minCoeff()).norm()
<< 100 * d.sum() / d.size() / (targetV.colwise().maxCoeff() - targetV.colwise().minCoeff()).norm()
<< "\n";
std::cout << "Max distance = "
<< 100 * d.lpNorm<Eigen::Infinity>() /
(VTarget.colwise().maxCoeff() - VTarget.colwise().minCoeff()).norm()
(targetV.colwise().maxCoeff() - targetV.colwise().minCoeff()).norm()
<< "\n";
polyscope::getSurfaceMesh("Simulation")->addVertexScalarQuantity("Distance", d)->setEnabled(true);
}
Expand Down
81 changes: 44 additions & 37 deletions src/newton.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,8 +3,8 @@
#include "functions.h"
#include "parameterization.h"
#include "simulation_utils.h"
#include "timer.h"
#include "solvers.h"
#include "timer.h"

#include <TinyAD/Utils/NewtonDecrement.hh>

Expand All @@ -17,7 +17,7 @@ void newton(Eigen::VectorXd& x,
int max_iters,
double lim,
bool verbose,
std::vector<int> fixedIdx,
const std::vector<int>& fixedIdx,
const std::function<void(const Eigen::VectorXd&)>& callback)
{
Timer timer("Newton", !verbose);
Expand Down Expand Up @@ -93,7 +93,7 @@ void newton(IntrinsicGeometryInterface& geometry,
int max_iters,
double lim,
bool verbose,
std::vector<int> fixedIdx,
const std::vector<int>& fixedIdx,
const std::function<void(const Eigen::VectorXd&)>& callback)
{
// Assemble inital x vector
Expand All @@ -108,24 +108,24 @@ void newton(IntrinsicGeometryInterface& geometry,
func.x_to_data(x, [&](Vertex v, const auto& row) { V.row(geometry.vertexIndices[v]) = row; });
}

void sparse_gauss_newton(IntrinsicGeometryInterface& geometry,
Eigen::MatrixXd& V,
const Eigen::VectorXd& xTarget,
const FaceData<Eigen::Matrix2d>& MrInv,
FaceData<double>& theta1,
VertexData<double>& theta2,
TinyAD::ScalarFunction<1, double, Eigen::Index>& adjointFunc,
std::vector<int>& fixedIdx,
int max_iters,
double lim,
double wM,
double wL,
double E1,
double lambda1,
double lambda2,
double deltaLambda,
double thickness,
const std::function<void(const Eigen::VectorXd&)>& callback)
Eigen::MatrixXd
sparse_gauss_newton(IntrinsicGeometryInterface& geometry,
const Eigen::MatrixXd& targetV,
const FaceData<Eigen::Matrix2d>& MrInv,
const FaceData<double>& theta1,
VertexData<double>& theta2,
const TinyAD::ScalarFunction<1, double, Eigen::Index>& adjointFunc,
const std::vector<int>& fixedIdx,
int max_iters,
double lim,
double wM,
double wL,
double E1,
double lambda1,
double lambda2,
double deltaLambda,
double thickness,
const std::function<void(const Eigen::VectorXd&)>& callback)
{
geometry.requireCotanLaplacian();
geometry.requireVertexLumpedMassMatrix();
Expand All @@ -135,7 +135,7 @@ void sparse_gauss_newton(IntrinsicGeometryInterface& geometry,
SurfaceMesh& mesh = geometry.mesh;

// build vector of Voronoi areas for the distance metric
Eigen::VectorXd masses = Eigen::VectorXd::Zero(V.size());
Eigen::VectorXd masses = Eigen::VectorXd::Zero(targetV.size());

// divide M by the total mesh area
double totalArea = 0;
Expand All @@ -154,12 +154,16 @@ void sparse_gauss_newton(IntrinsicGeometryInterface& geometry,
Eigen::SparseMatrix<double> L = geometry.cotanLaplacian;

// Mass matrix theta
Eigen::SparseMatrix<double> M_theta(V.rows(), V.rows());
M_theta.reserve(V.rows());
for(int i = 0; i < V.rows(); ++i)
Eigen::SparseMatrix<double> M_theta(targetV.rows(), targetV.rows());
M_theta.reserve(targetV.rows());
for(int i = 0; i < targetV.rows(); ++i)
M_theta.insert(i, i) = totalArea * masses(3 * i);

Eigen::VectorXd theta = theta2.toVector();
Eigen::VectorXd xTarget(targetV.size());
for(int i = 0; i < targetV.rows(); ++i)
for(int j = 0; j < 3; ++j)
xTarget(3 * i + j) = targetV(i, j);
Eigen::VectorXd x = xTarget;

LLTSolver adjointSolver;
Expand All @@ -176,30 +180,30 @@ void sparse_gauss_newton(IntrinsicGeometryInterface& geometry,
Eigen::SparseMatrix<double> P = projectionMatrix(fixedIdx, x.size());

// Hessian matrix H
Eigen::VectorXd X(V.size() + theta.size());
X.head(V.size()) = x;
Eigen::VectorXd X(targetV.size() + theta.size());
X.head(targetV.size()) = x;
X.tail(theta.size()) = theta;
Eigen::SparseMatrix<double> H = adjointFunc.eval_hessian(X);

// Build HGN matrix
Eigen::SparseMatrix<double> HGN = buildHGN(2 * masses, P, 2 * wM * M_theta + 2 * wL * L, H);

auto distanceGrad = [&](const Eigen::VectorXd& th) -> Eigen::VectorXd {
Eigen::VectorXd X(V.size() + th.size());
X.head(V.size()) = x;
Eigen::VectorXd X(targetV.size() + th.size());
X.head(targetV.size()) = x;
X.tail(th.size()) = th;
H = adjointFunc.eval_hessian(X);

for(int j = 0; j < V.size(); ++j)
for(int j = 0; j < targetV.size(); ++j)
H.coeffRef(j, j) += 1e-10;

Eigen::SparseMatrix<double> A = (P * H.block(0, 0, V.size(), V.size()) * P.transpose()).eval();
Eigen::SparseMatrix<double> A = (P * H.block(0, 0, targetV.size(), targetV.size()) * P.transpose()).eval();

adjointSolver.factorize(A);
if(adjointSolver.info() != Eigen::Success)
{
auto [f, g, A_proj] = adjointFunc.eval_with_hessian_proj(X);
A_proj = (P * A_proj.block(0, 0, V.size(), V.size()) * P.transpose()).eval();
A_proj = (P * A_proj.block(0, 0, targetV.size(), targetV.size()) * P.transpose()).eval();

A = 0.9 * A + 0.1 * A_proj;
adjointSolver.factorize(A);
Expand All @@ -214,7 +218,7 @@ void sparse_gauss_newton(IntrinsicGeometryInterface& geometry,

dir = P.transpose() * dir;

return -2 * H.block(V.size(), 0, th.size(), V.size()) * dir + 2 * wM * M_theta * th + 2 * wL * L * th;
return -2 * H.block(targetV.size(), 0, th.size(), targetV.size()) * dir + 2 * wM * M_theta * th + 2 * wL * L * th;
};

double energy = distance(theta);
Expand Down Expand Up @@ -242,7 +246,7 @@ void sparse_gauss_newton(IntrinsicGeometryInterface& geometry,
if(solver.info() != Eigen::Success)
{
std::cout << "Solver error\n";
return;
return targetV;
}

Eigen::VectorXd d = solver.solve(b);
Expand Down Expand Up @@ -271,10 +275,13 @@ void sparse_gauss_newton(IntrinsicGeometryInterface& geometry,

std::cout << "Final energy: " << distance(theta) << "\n";

for(int i = 0; i < V.rows(); ++i)
Eigen::MatrixXd V(targetV.rows(), 3);
for(int i = 0; i < targetV.rows(); ++i)
for(int j = 0; j < 3; ++j)
V(i, j) = x(3 * i + j);

theta2.fromVector(theta);
return V;
}

template void newton<>(IntrinsicGeometryInterface&,
Expand All @@ -283,7 +290,7 @@ template void newton<>(IntrinsicGeometryInterface&,
int,
double,
bool,
std::vector<int>,
const std::vector<int>&,
const std::function<void(const Eigen::VectorXd&)>&);

template void newton<>(IntrinsicGeometryInterface&,
Expand All @@ -292,5 +299,5 @@ template void newton<>(IntrinsicGeometryInterface&,
int,
double,
bool,
std::vector<int>,
const std::vector<int>&,
const std::function<void(const Eigen::VectorXd&)>&);
40 changes: 20 additions & 20 deletions src/newton.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ void newton(
int max_iters = 1000,
double lim = 1e-6,
bool verbose = true,
std::vector<int> fixedIdx = {},
const std::vector<int>& fixedIdx = {},
const std::function<void(const Eigen::VectorXd&)>& callBack = [](const auto&) {});

template <class Func>
Expand All @@ -26,24 +26,24 @@ void newton(
int max_iters,
double lim,
bool verbose = true,
std::vector<int> fixedIdx = {},
const std::vector<int>& fixedIdx = {},
const std::function<void(const Eigen::VectorXd&)>& callBack = [](const auto&) {});

void sparse_gauss_newton(geometrycentral::surface::IntrinsicGeometryInterface& geometry,
Eigen::MatrixXd& V,
const Eigen::VectorXd& xTarget,
const geometrycentral::surface::FaceData<Eigen::Matrix2d>& MrInv,
geometrycentral::surface::FaceData<double>& theta1,
geometrycentral::surface::VertexData<double>& theta2,
TinyAD::ScalarFunction<1, double, Eigen::Index>& adjointFunc,
std::vector<int>& fixedIdx,
int max_iters,
double lim,
double wM,
double wL,
double E1,
double lambda1,
double lambda2,
double deltaLambda,
double thickness,
const std::function<void(const Eigen::VectorXd&)>& callback = [](const auto&) {});
Eigen::MatrixXd sparse_gauss_newton(
geometrycentral::surface::IntrinsicGeometryInterface& geometry,
const Eigen::MatrixXd& targetV,
const geometrycentral::surface::FaceData<Eigen::Matrix2d>& MrInv,
const geometrycentral::surface::FaceData<double>& theta1,
geometrycentral::surface::VertexData<double>& theta2,
const TinyAD::ScalarFunction<1, double, Eigen::Index>& adjointFunc,
const std::vector<int>& fixedIdx,
int max_iters,
double lim,
double wM,
double wL,
double E1,
double lambda1,
double lambda2,
double deltaLambda,
double thickness,
const std::function<void(const Eigen::VectorXd&)>& callback = [](const auto&) {});

0 comments on commit b543f08

Please sign in to comment.