Skip to content

Commit 9d88330

Browse files
mlivesumlivesu
committed
fixed ARAP (solver caching when using hard constraints missing))
1 parent acda13a commit 9d88330

File tree

2 files changed

+92
-204
lines changed

2 files changed

+92
-204
lines changed

‎include/cinolib/ARAP.cpp‎

Lines changed: 81 additions & 181 deletions
Original file line numberDiff line numberDiff line change
@@ -35,6 +35,7 @@
3535
*********************************************************************************/
3636
#include <cinolib/ARAP.h>
3737
#include <cinolib/parallel_for.h>
38+
#include <cinolib/laplacian.h>
3839

3940
namespace cinolib
4041
{
@@ -43,236 +44,137 @@ template<class M, class V, class E, class P>
4344
CINO_INLINE
4445
void ARAP(AbstractMesh<M,V,E,P> & m, ARAP_data & data)
4546
{
47+
/////////////////////////////////////////////////////////////////////////
48+
4649
auto init = [&]()
4750
{
4851
assert(m.mesh_type()==TRIMESH || m.mesh_type()==TETMESH);
4952

50-
data.init = false; // don't init next time
51-
data.xyz_loc.resize(m.num_polys()*m.verts_per_poly(0));
53+
data.init = false;
54+
55+
data.R.resize(m.num_verts());
56+
data.xyz_ref = m.vector_verts();
5257

53-
// per edge weights
5458
data.w.resize(m.num_edges());
5559
for(uint eid=0; eid<m.num_edges(); ++eid)
5660
{
5761
data.w.at(eid) = m.edge_weight(eid,data.w_type);
5862
}
5963

60-
// compute a map between matrix columns and mesh vertices
61-
// if hard constraints are used, boundary conditions will
62-
// map to -1, meaning that they do not correspond to any
63-
// column in the matrix
64-
data.col_map.resize(m.num_verts(),0);
65-
if(!data.use_soft_constraints)
66-
{
67-
for(const auto & bc : data.bcs)
68-
{
69-
data.col_map.at(bc.first) = -1;
70-
}
71-
}
72-
uint fresh_id = 0;
73-
for(uint vid=0; vid<m.num_verts(); ++vid)
74-
{
75-
if(data.col_map.at(vid)==0)
76-
{
77-
data.col_map.at(vid) = fresh_id++;
78-
}
79-
}
80-
81-
// Compute the Laplacian matrix and pre-factorize it
82-
typedef Eigen::Triplet<double> Entry;
83-
std::vector<Entry> entries;
84-
uint size = m.num_verts() - uint(data.bcs.size());
85-
if(data.use_soft_constraints)
64+
if(data.hard_constrain_handles)
8665
{
87-
size = m.num_verts() + uint(data.bcs.size());
88-
data.W.resize(size);
66+
data.A = -laplacian(m,data.w_type,1);
8967
}
90-
for(uint vid=0; vid<m.num_verts(); ++vid)
91-
{
92-
int col = data.col_map.at(vid);
93-
if(col==-1) continue; // skip, hard BC
94-
if(data.use_soft_constraints) data.W[col] = data.w_laplace;
95-
for(uint nbr : m.adj_v2v(vid))
96-
{
97-
int col_nbr = data.col_map.at(nbr);
98-
int eid = m.edge_id(vid,nbr);
99-
assert(eid>=0);
100-
entries.push_back(Entry(col, col, data.w.at(eid)));
101-
if(col_nbr==-1) continue; // skip, hard BC
102-
entries.push_back(Entry(col, col_nbr, -data.w.at(eid)));
103-
}
104-
}
105-
if(data.use_soft_constraints)
106-
{
107-
// models equation => x_bc = bc_value * w_constr
108-
uint new_row = m.num_verts();
109-
for(auto bc : data.bcs)
110-
{
111-
entries.push_back(Entry(new_row,bc.first,1));
112-
data.W[new_row] = data.w_constr;
113-
++new_row;
114-
}
115-
}
116-
data.A = Eigen::SparseMatrix<double>(size, (data.use_soft_constraints) ? m.num_verts() : size);
117-
data.A.setFromTriplets(entries.begin(), entries.end());
118-
if(data.use_soft_constraints) data.cache.derived().compute(data.A.transpose()*data.W.asDiagonal()*data.A);
119-
else data.cache.derived().compute(data.A);
120-
121-
if(data.warm_start_with_laplacian)
68+
else
12269
{
123-
Eigen::VectorXd rhs_x = Eigen::VectorXd::Zero(size);
124-
Eigen::VectorXd rhs_y = Eigen::VectorXd::Zero(size);
125-
Eigen::VectorXd rhs_z = Eigen::VectorXd::Zero(size);
70+
uint nv = m.num_verts();
71+
uint nh = data.handles.size();
72+
std::vector<Eigen::Triplet<double>> entries;
12673
for(uint vid=0; vid<m.num_verts(); ++vid)
12774
{
128-
int col = data.col_map.at(vid);
129-
if(col==-1) continue; // skip BC
130-
for(uint eid : m.adj_v2e(vid))
75+
double diag = 0;
76+
for(uint nbr : m.adj_v2v(vid))
13177
{
132-
uint nbr = m.vert_opposite_to(eid,vid);
133-
rhs_x[col] += data.w.at(eid) * (m.vert(vid).x() - m.vert(nbr).x());
134-
rhs_y[col] += data.w.at(eid) * (m.vert(vid).y() - m.vert(nbr).y());
135-
rhs_z[col] += data.w.at(eid) * (m.vert(vid).z() - m.vert(nbr).z());
136-
// move the contribution of BCs to the RHS
137-
if(data.col_map.at(nbr)==-1)
138-
{
139-
vec3d p = data.bcs.at(nbr);
140-
rhs_x[col] += data.w.at(eid) * p.x();
141-
rhs_y[col] += data.w.at(eid) * p.y();
142-
rhs_z[col] += data.w.at(eid) * p.z();
143-
}
78+
int eid = m.edge_id(vid,nbr);
79+
entries.emplace_back(vid,nbr,-data.w.at(eid));
80+
diag += data.w.at(eid);
14481
}
82+
entries.emplace_back(vid,vid,diag);
14583
}
146-
if(data.use_soft_constraints)
84+
uint off = 0;
85+
for(uint vid : data.handles)
14786
{
148-
auto x = data.cache.solve(data.A.transpose()*data.W.asDiagonal()*rhs_x).eval();
149-
auto y = data.cache.solve(data.A.transpose()*data.W.asDiagonal()*rhs_y).eval();
150-
auto z = data.cache.solve(data.A.transpose()*data.W.asDiagonal()*rhs_z).eval();
151-
data.xyz_out.resize(m.num_verts());
152-
for(uint vid=0; vid<m.num_verts(); ++vid)
153-
{
154-
data.xyz_out[vid] = vec3d(x[vid],y[vid],z[vid]);
155-
}
156-
}
157-
else
158-
{
159-
auto x = data.cache.solve(rhs_x).eval();
160-
auto y = data.cache.solve(rhs_y).eval();
161-
auto z = data.cache.solve(rhs_z).eval();
162-
data.xyz_out.resize(m.num_verts());
163-
for(uint vid=0; vid<m.num_verts(); ++vid)
164-
{
165-
int col = data.col_map[vid];
166-
if(col>=0) data.xyz_out[vid] = vec3d(x[col],y[col],z[col]);
167-
}
168-
for(const auto & bc : data.bcs)
169-
{
170-
data.xyz_out[bc.first] = bc.second;
171-
}
87+
entries.emplace_back(nv+off,vid,1.0);
88+
off++;
17289
}
173-
}
174-
else
175-
{
176-
data.xyz_out = m.vector_verts();
177-
for(const auto & bc : data.bcs) data.xyz_out.at(bc.first) = bc.second;
90+
data.A.resize(nv+nh,nv);
91+
data.A.setFromTriplets(entries.begin(), entries.end());
92+
data.cache.derived().compute(data.A.transpose()*data.A);
17893
}
17994
};
18095

96+
/////////////////////////////////////////////////////////////////////////
97+
18198
auto local_step = [&]()
18299
{
183-
PARALLEL_FOR(0, m.num_polys(), 1000, [&](uint pid)
100+
PARALLEL_FOR(0, m.num_verts(), 1000, [&](uint vid)
184101
{
185-
mat3d cov = mat3d::ZERO();
186-
for(uint eid : m.adj_p2e(pid))
102+
mat3d cov = mat3d::ZERO();
103+
for(uint nbr : m.adj_v2v(vid))
187104
{
188-
uint v0 = m.edge_vert_id(eid,0);
189-
uint v1 = m.edge_vert_id(eid,1);
190-
vec3d e_cur = data.xyz_out.at(v0) - data.xyz_out.at(v1);
191-
vec3d e_ref = m.vert(v0) - m.vert(v1);
105+
int eid = m.edge_id(vid,nbr);
106+
vec3d e_ref = (data.xyz_ref.at(vid) - data.xyz_ref.at(nbr));
107+
vec3d e_cur = (m.vert(vid) - m.vert(nbr));
108+
192109
cov += data.w.at(eid) * (e_cur * e_ref.transpose());
193110
}
194111

195-
// find closest rotation and store rotated point
196-
uint off = pid*m.verts_per_poly(pid);
197-
mat3d rot = cov.closest_orthogonal_matrix(true);
198-
for(uint i=0; i<m.verts_per_poly(pid); ++i)
199-
{
200-
data.xyz_loc.at(off+i) = rot * m.poly_vert(pid,i);
201-
}
112+
mat3d UU,VV;
113+
vec3d S;
114+
cov.SVD(UU,S,VV);
115+
mat3d I = mat3d::DIAG(1);
116+
I(2,2) = (UU*VV.transpose()).det();
117+
data.R.at(vid) = UU*I*VV.transpose();
202118
});
203119
};
204120

121+
/////////////////////////////////////////////////////////////////////////
122+
205123
auto global_step = [&]()
206124
{
207-
uint size = (data.use_soft_constraints) ? m.num_verts() + uint(data.bcs.size())
208-
: m.num_verts() - uint(data.bcs.size());
209-
Eigen::VectorXd rhs_x = Eigen::VectorXd::Zero(size);
210-
Eigen::VectorXd rhs_y = Eigen::VectorXd::Zero(size);
211-
Eigen::VectorXd rhs_z = Eigen::VectorXd::Zero(size);
125+
uint nv = m.num_verts();
126+
uint nh = (data.hard_constrain_handles) ? 0 : data.handles.size();
127+
128+
Eigen::VectorXd rhs_x(nv+nh);
129+
Eigen::VectorXd rhs_y(nv+nh);
130+
Eigen::VectorXd rhs_z(nv+nh);
131+
212132
for(uint vid=0; vid<m.num_verts(); ++vid)
213133
{
214-
int col = data.col_map.at(vid);
215-
if(col==-1) continue; // skip, vert is BC
134+
vec3d rhs(0,0,0);
216135
for(uint nbr : m.adj_v2v(vid))
217136
{
218-
int eid = m.edge_id(vid,nbr);
219-
double w = 1.0/m.adj_e2p(eid).size();
220-
for(uint pid : m.adj_e2p(eid))
221-
{
222-
uint i = m.poly_vert_offset(pid,vid);
223-
uint j = m.poly_vert_offset(pid,nbr);
224-
vec3d Re = data.xyz_loc.at(pid*m.verts_per_poly(pid)+i) - data.xyz_loc.at(pid*m.verts_per_poly(pid)+j);
225-
rhs_x[col] += w * data.w.at(eid) * Re[0];
226-
rhs_y[col] += w * data.w.at(eid) * Re[1];
227-
rhs_z[col] += w * data.w.at(eid) * Re[2];
228-
}
229-
// if nbr is a hard BC sum its contibution to the Laplacian matrix to the rhs
230-
if(data.col_map.at(nbr)==-1)
231-
{
232-
vec3d p = data.bcs.at(nbr);
233-
rhs_x[col] += data.w.at(eid) * p.x();
234-
rhs_y[col] += data.w.at(eid) * p.y();
235-
rhs_z[col] += data.w.at(eid) * p.z();
236-
}
137+
mat3d Ravg = (data.R.at(vid)+data.R.at(nbr))/2.0;
138+
vec3d e = (data.xyz_ref.at(vid) - data.xyz_ref.at(nbr));
139+
int eid = m.edge_id(vid,nbr);
140+
141+
rhs += data.w.at(eid) * Ravg * e;
237142
}
143+
rhs_x[vid] = rhs.x();
144+
rhs_y[vid] = rhs.y();
145+
rhs_z[vid] = rhs.z();
238146
}
239-
if(data.use_soft_constraints)
147+
148+
Eigen::VectorXd x,y,z;
149+
if(data.hard_constrain_handles)
240150
{
241-
// models rhs of equation => x_bc = bc_value
242-
uint new_row = m.num_verts();
243-
for(auto bc : data.bcs)
244-
{
245-
rhs_x[new_row] = bc.second.x();
246-
rhs_y[new_row] = bc.second.y();
247-
rhs_z[new_row] = bc.second.z();
248-
++new_row;
249-
}
250-
auto x = data.cache.solve(data.A.transpose()*data.W.asDiagonal()*rhs_x).eval();
251-
auto y = data.cache.solve(data.A.transpose()*data.W.asDiagonal()*rhs_y).eval();
252-
auto z = data.cache.solve(data.A.transpose()*data.W.asDiagonal()*rhs_z).eval();
253-
for(uint vid=0; vid<m.num_verts(); ++vid)
254-
{
255-
data.xyz_out[vid] = vec3d(x[vid],y[vid],z[vid]);
256-
}
151+
// TODO: use caching also with hard constraints!
152+
solve_square_system_with_bc(data.A, rhs_x, x, data.handles_x);
153+
solve_square_system_with_bc(data.A, rhs_y, y, data.handles_y);
154+
solve_square_system_with_bc(data.A, rhs_z, z, data.handles_z);
257155
}
258156
else
259157
{
260-
auto x = data.cache.solve(rhs_x).eval();
261-
auto y = data.cache.solve(rhs_y).eval();
262-
auto z = data.cache.solve(rhs_z).eval();
263-
for(uint vid=0; vid<m.num_verts(); ++vid)
158+
uint off = 0;
159+
for(uint vid : data.handles)
264160
{
265-
int col = data.col_map[vid];
266-
if(col>=0) data.xyz_out[vid] = vec3d(x[col],y[col],z[col]);
267-
}
268-
for(const auto & bc : data.bcs)
269-
{
270-
data.xyz_out[bc.first] = bc.second;
161+
rhs_x[nv+off] = data.handles_x.at(vid);
162+
rhs_y[nv+off] = data.handles_y.at(vid);
163+
rhs_z[nv+off] = data.handles_z.at(vid);
164+
++off;
271165
}
166+
x = data.cache.solve(data.A.transpose()*rhs_x).eval();
167+
y = data.cache.solve(data.A.transpose()*rhs_y).eval();
168+
z = data.cache.solve(data.A.transpose()*rhs_z).eval();
169+
}
170+
171+
for(uint vid=0; vid<m.num_verts(); ++vid)
172+
{
173+
m.vert(vid) = vec3d(x[vid],y[vid],z[vid]);
272174
}
273175
};
274176

275-
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
177+
//////////////////////////////////////////////////////////////////////////
276178

277179
if(data.init) init();
278180

@@ -281,8 +183,6 @@ void ARAP(AbstractMesh<M,V,E,P> & m, ARAP_data & data)
281183
local_step();
282184
global_step();
283185
}
284-
285-
m.vector_verts() = data.xyz_out;
286186
m.update_normals();
287187
}
288188

‎include/cinolib/ARAP.h‎

Lines changed: 11 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -36,7 +36,7 @@
3636
#ifndef CINO_ARAP_H
3737
#define CINO_ARAP_H
3838

39-
#include <cinolib/meshes/trimesh.h>
39+
#include <cinolib/meshes/abstract_mesh.h>
4040
#include <cinolib/linear_solvers.h>
4141

4242
namespace cinolib
@@ -56,33 +56,21 @@ struct ARAP_data
5656
uint n_iters = 4;
5757
bool init = true; // initialize just once (useful for multiple calls, e.g. to make more iterations)
5858

59-
std::map<uint,vec3d> bcs; // boundary conditions
60-
std::vector<int> col_map; // maps vertex ids to matrix columns
59+
std::vector<mat3d> R; // local (per vertex) rotation matrices
60+
std::vector<vec3d> xyz_ref; // reference (original) vertex positions
6161

62-
std::vector<vec3d> xyz_out; // current solution (will be the output, eventually)
63-
std::vector<vec3d> xyz_loc; // per element targets (100% rigid)
64-
std::vector<double> w; // edge weights { UNIFORM, COTANGENT }
65-
int w_type = UNIFORM; // WARNING: cot weights seem rather unstable on volume meshes in interactive deformations
62+
std::vector<double> w; // edge weights { UNIFORM, COTANGENT }
63+
int w_type = UNIFORM; // WARNING: cot weights seem rather unstable on volume meshes in interactive deformations
6664

6765
Eigen::SimplicialLLT<Eigen::SparseMatrix<double>> cache; // factorized matrix
68-
69-
// In my experience replacing hard with soft constraints works
70-
// much better for interactive shape deformation (no artifacts
71-
// due to C-0 continuity around boundary conditions).
72-
// On the other hand, things like simplicial volume maps require
73-
// the use of hard constraints to ensure boundary conformity.
74-
bool use_soft_constraints = true;
75-
double w_constr = 100.0; // weight for soft constraints
76-
double w_laplace = 1.0; // weight for the laplacian component of the matrix
77-
Eigen::VectorXd W; // diagonal matrix of weights
7866
Eigen::SparseMatrix<double> A; // a copy of the matrix (to be pre-multiplied to the rhs to form the normal equations)
7967

80-
// if true (default), the warm start will be the minimizer of
81-
// | Lx - delta |^2
82-
// where delta are the differential mesh coordinates. If false,
83-
// the warm start will simply be a copy of the input mesh, with
84-
// constrained vertices moved onto their prescribed position
85-
bool warm_start_with_laplacian = true;
68+
// deformation handles, separated for x,y,z coords to make solver call easier
69+
bool hard_constrain_handles = true;
70+
std::vector<uint> handles;
71+
std::map<uint,double> handles_x;
72+
std::map<uint,double> handles_y;
73+
std::map<uint,double> handles_z;
8674
};
8775

8876
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

0 commit comments

Comments
 (0)