mirror of
https://github.com/lucaspalomodevelop/meshlab.git
synced 2026-03-16 17:44:36 +00:00
pre-conversion to full explicit
This commit is contained in:
parent
4da7beea6a
commit
706bf71e88
@ -4,7 +4,6 @@
|
||||
#include "vcg/complex/trimesh/update/curvature.h"
|
||||
// #include "vcg/complex/trimesh/update/curvature_fitting.h"
|
||||
|
||||
|
||||
using namespace vcg;
|
||||
|
||||
//---------------------------------------------------------------------------------------//
|
||||
@ -30,8 +29,8 @@ void Balloon::init( int gridsize, int gridpad ){
|
||||
// Remember that rays will take a step JUST in their direction, so if they lie exactly
|
||||
// on the bbox, they would go outside. The code below corrects this from happening.
|
||||
Box3f enlargedbb = cloud.bbox;
|
||||
// float del = .99*vol.getDelta(); // almost +1 voxel in each direction
|
||||
float del = .50*vol.getDelta(); // ADEBUG: almost to debug correspondences
|
||||
float del = .99*vol.getDelta(); // almost +1 voxel in each direction
|
||||
// float del = .50*vol.getDelta(); // ADEBUG: almost to debug correspondences
|
||||
Point3f offset( del,del,del );
|
||||
enlargedbb.Offset( offset );
|
||||
vol.initField( enlargedbb ); // init volumetric field with the bounding box
|
||||
@ -55,9 +54,9 @@ void Balloon::init( int gridsize, int gridpad ){
|
||||
bool Balloon::initializeField(){
|
||||
//--- Setup the interpolation system
|
||||
// if we fail, signal the failure to the caller
|
||||
float OMEGA = 1e8;
|
||||
float OMEGA = 1e7; // 1e8
|
||||
LAPLACIAN type = COTANGENT; // COMBINATORIAL
|
||||
bool op_succeed = finterp.Init( &surf, type );
|
||||
bool op_succeed = finterp.Init( &surf, 1, type );
|
||||
if( !op_succeed ){
|
||||
finterp.ColorizeIllConditioned( type );
|
||||
return false;
|
||||
@ -137,7 +136,7 @@ bool Balloon::initializeField(){
|
||||
tri::UpdateQuality<CMeshO>::FaceConstant(surf, 0);
|
||||
std::vector<float> tot_w( surf.fn, 0 );
|
||||
|
||||
// We clear the ray-triangles correspondence information + distance information
|
||||
// We clear the pokingRay-triangle correspondence information + distance information
|
||||
// to get ready for the next step
|
||||
gridAccell.clearCorrespondences();
|
||||
|
||||
@ -160,30 +159,39 @@ bool Balloon::initializeField(){
|
||||
// with him if and only if this face is the closest to it. Note that we study
|
||||
// the values of t,u,v directly as a t<0 is accepted as intersecting.
|
||||
for(unsigned int i=0; i<prays.size(); i++){
|
||||
vcg::IntersectionRayTriangle<float>(prays[i]->ray, f.P(0), f.P(1), f.P(2), t, u, v);
|
||||
|
||||
Line3<float> line(prays[i]->ray.Origin(), prays[i]->ray.Direction());
|
||||
|
||||
// If the ray falls within the domain of the face
|
||||
if( u>=0 && u<=1 && v>=0 && v<=1 ){
|
||||
if( IntersectionLineTriangle(line, f.P(0), f.P(1), f.P(2), t, u, v) ){
|
||||
|
||||
// DEBUG
|
||||
// f.C() = t>0 ? Color4b(0,0,255,255) : f.C() = Color4b(255,0,0,255);
|
||||
|
||||
// DEBUG
|
||||
// if( prays[i]->f!=NULL && fabs(t)<fabs(prays[i]->t) ){
|
||||
// prays[i]->f->C() = Color4b(255,0,0,255);
|
||||
// qDebug() << "replaced: " << tri::Index(surf,prays[i]->f) << " from t: " << prays[i]->t << " to " << tri::Index(surf,f) << t;
|
||||
// }
|
||||
|
||||
// If no face was associated with this ray or this face is closer
|
||||
// than the one that I stored previously
|
||||
if( prays[i]->f==NULL || fabs(prays[i]->t)<fabs(t) ){
|
||||
if( prays[i]->f==NULL || fabs(t)<fabs(prays[i]->t) ){
|
||||
prays[i]->f=&f;
|
||||
prays[i]->t=t;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Now we scan through the rays, we visit the "best" corresponding face and we
|
||||
// set a constraint on this face. Also we modify the color of the face so that
|
||||
// an approximation can be visualized
|
||||
for(unsigned int i=0; i<gridAccell.rays.size(); i++){
|
||||
// Retrieve the corresponding face and signed distance
|
||||
Ray3f& ray = gridAccell.rays[i].ray;
|
||||
// Removing degenerate triangles could cause a ray to fail the intersection test!!
|
||||
// The use of vcg::Simplify to remove degenerate triangles could cause a ray to fail the intersection test!!
|
||||
if( gridAccell.rays[i].f == NULL){
|
||||
qDebug() << "warning: ray #" << i << "has provided NULL intersection";
|
||||
// << toString(ray.Origin()) << " " << toString(ray.Direction());
|
||||
qDebug() << "warning: ray #" << i << "has provided NULL intersection"; // << toString(ray.Origin()) << " " << toString(ray.Direction());
|
||||
continue;
|
||||
}
|
||||
CFaceO& f = *(gridAccell.rays[i].f);
|
||||
@ -194,7 +202,10 @@ bool Balloon::initializeField(){
|
||||
tot_w[ tri::Index(surf,f) ]++;
|
||||
f.Q() += t; // normalize with tot_w after
|
||||
f.SetS(); // enable the face for coloring
|
||||
|
||||
|
||||
// DEBUG
|
||||
// f.C() = Color4b(0,255,0,255);
|
||||
|
||||
// I was lazy and didn't store the u,v... we need to recompute them
|
||||
vcg::IntersectionRayTriangle<float>(ray, f.P(0), f.P(1), f.P(2), t, u, v);
|
||||
|
||||
@ -202,6 +213,7 @@ bool Balloon::initializeField(){
|
||||
finterp.AddConstraint( tri::Index(surf,f.V(0)), OMEGA*(1-u-v), t );
|
||||
finterp.AddConstraint( tri::Index(surf,f.V(1)), OMEGA*(u), t );
|
||||
finterp.AddConstraint( tri::Index(surf,f.V(2)), OMEGA*(v), t );
|
||||
assert( u>=0 && u<=1 && v>=0 && v<=1 );
|
||||
}
|
||||
|
||||
//--- Normalize in case there is more than 1 ray per-face
|
||||
@ -209,8 +221,7 @@ bool Balloon::initializeField(){
|
||||
if( tot_w[ tri::Index(surf,*fi) ] > 0 )
|
||||
fi->Q() /= tot_w[ tri::Index(surf,*fi) ];
|
||||
}
|
||||
//--- Transfer the average distance stored in face quality to a color
|
||||
// and do it only for the selection (true)
|
||||
//--- Transfer the average distance stored in face quality to a color and do it only for the selection (true)
|
||||
tri::UpdateColor<CMeshO>::FaceQualityRamp(surf, true);
|
||||
}
|
||||
|
||||
@ -221,7 +232,7 @@ void Balloon::interpolateField(){
|
||||
surf.vert.QualityEnabled = true;
|
||||
|
||||
//--- Interpolate the field
|
||||
finterp.Solve();
|
||||
finterp.SolveInQuality();
|
||||
|
||||
//--- Transfer vertex quality to surface
|
||||
rm &= ~SURF_FCOLOR; // disable face colors
|
||||
@ -231,39 +242,54 @@ void Balloon::interpolateField(){
|
||||
tri::UpdateColor<CMeshO>::VertexQualityRamp(surf,H.Percentile(0.0f),H.Percentile(1.0f));
|
||||
}
|
||||
void Balloon::computeCurvature(){
|
||||
// Enable curvature supports, How do I avoid a
|
||||
// double computation of topology here?
|
||||
surf.vert.EnableCurvature();
|
||||
surf.vert.EnableVFAdjacency();
|
||||
surf.face.EnableVFAdjacency();
|
||||
surf.face.EnableFFAdjacency();
|
||||
vcg::tri::UpdateTopology<CMeshO>::VertexFace( surf );
|
||||
vcg::tri::UpdateTopology<CMeshO>::FaceFace( surf );
|
||||
|
||||
// tri::UpdateCurvatureFitting<CMeshO>::computeCurvature( surf );
|
||||
|
||||
//--- Compute curvature and its bounds
|
||||
tri::UpdateCurvature<CMeshO>::MeanAndGaussian( surf );
|
||||
float absmax = -FLT_MAX;
|
||||
for(CMeshO::VertexIterator vi = surf.vert.begin(); vi != surf.vert.end(); ++vi){
|
||||
float cabs = fabs((*vi).Kh());
|
||||
absmax = (cabs>absmax) ? cabs : absmax;
|
||||
}
|
||||
|
||||
//--- Enable color coding
|
||||
rm &= ~SURF_FCOLOR;
|
||||
rm |= SURF_VCOLOR;
|
||||
|
||||
//--- Map curvature to two color ranges:
|
||||
// Blue => Yellow: negative values
|
||||
// Yellow => Red: positive values
|
||||
typedef unsigned char CT;
|
||||
for(CMeshO::VertexIterator vi = surf.vert.begin(); vi != surf.vert.end(); ++vi){
|
||||
if( (*vi).Kh() < 0 )
|
||||
(*vi).C().lerp(Color4<CT>::Yellow, Color4<CT>::Blue, fabs((*vi).Kh())/absmax );
|
||||
else
|
||||
(*vi).C().lerp(Color4<CT>::Yellow, Color4<CT>::Red, (*vi).Kh()/absmax);
|
||||
}
|
||||
#if FALSE
|
||||
float OMEGA = 1e2; // interpolation coefficient
|
||||
sinterp.Init( &surf, 3, COTANGENT );
|
||||
for( CMeshO::VertexIterator vi=surf.vert.begin();vi!=surf.vert.end();vi++ )
|
||||
sinterp.AddConstraint3( tri::Index(surf,*vi), OMEGA, (*vi).P()[0], (*vi).P()[1], (*vi).P()[2] );
|
||||
FieldInterpolator::XBType* coords[3];
|
||||
sinterp.Solve(coords);
|
||||
for( CMeshO::VertexIterator vi=surf.vert.begin();vi!=surf.vert.end();vi++ ){
|
||||
int vIdx = tri::Index(surf,*vi);
|
||||
(*vi).P()[0] = (*coords[0])(vIdx);
|
||||
(*vi).P()[1] = (*coords[1])(vIdx);
|
||||
(*vi).P()[2] = (*coords[2])(vIdx);
|
||||
}
|
||||
#else
|
||||
// Enable curvature supports, How do I avoid a
|
||||
// double computation of topology here?
|
||||
surf.vert.EnableCurvature();
|
||||
surf.vert.EnableVFAdjacency();
|
||||
surf.face.EnableVFAdjacency();
|
||||
surf.face.EnableFFAdjacency();
|
||||
vcg::tri::UpdateTopology<CMeshO>::VertexFace( surf );
|
||||
vcg::tri::UpdateTopology<CMeshO>::FaceFace( surf );
|
||||
|
||||
// tri::UpdateCurvatureFitting<CMeshO>::computeCurvature( surf );
|
||||
|
||||
//--- Compute curvature and its bounds
|
||||
tri::UpdateCurvature<CMeshO>::MeanAndGaussian( surf );
|
||||
float absmax = -FLT_MAX;
|
||||
for(CMeshO::VertexIterator vi = surf.vert.begin(); vi != surf.vert.end(); ++vi){
|
||||
float cabs = fabs((*vi).Kh());
|
||||
absmax = (cabs>absmax) ? cabs : absmax;
|
||||
}
|
||||
|
||||
//--- Enable color coding
|
||||
rm &= ~SURF_FCOLOR;
|
||||
rm |= SURF_VCOLOR;
|
||||
|
||||
//--- Map curvature to two color ranges:
|
||||
// Blue => Yellow: negative values
|
||||
// Yellow => Red: positive values
|
||||
typedef unsigned char CT;
|
||||
for(CMeshO::VertexIterator vi = surf.vert.begin(); vi != surf.vert.end(); ++vi){
|
||||
if( (*vi).Kh() < 0 )
|
||||
(*vi).C().lerp(Color4<CT>::Yellow, Color4<CT>::Blue, fabs((*vi).Kh())/absmax );
|
||||
else
|
||||
(*vi).C().lerp(Color4<CT>::Yellow, Color4<CT>::Red, (*vi).Kh()/absmax);
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
// HP: a correspondence has already been executed once!
|
||||
@ -363,10 +389,10 @@ void Balloon::evolveBalloon(){
|
||||
}
|
||||
// if we don't have computed the distance field, we don't really know how to
|
||||
// modulate the laplacian accordingly...
|
||||
#ifdef ENABLE_CURVATURE_IN_EVOLUTION
|
||||
if( surf.vert.CurvatureEnabled && surf.vert.QualityEnabled ){
|
||||
v.sfield += .1*k3*k2;
|
||||
}
|
||||
#ifdef TRUE // ENABLE_CURVATURE_IN_EVOLUTION
|
||||
// if( surf.vert.CurvatureEnabled && surf.vert.QualityEnabled ){
|
||||
// v.sfield += .1*k3*k2;
|
||||
// }
|
||||
else if( surf.vert.CurvatureEnabled ){
|
||||
v.sfield += .1*k3;
|
||||
}
|
||||
|
||||
@ -38,8 +38,8 @@ public:
|
||||
GridAccell gridAccell;
|
||||
/// Scalar field interpolator (one constraint per poking ray)
|
||||
FieldInterpolator finterp;
|
||||
/// Smootfield interpolator (one constraint per vertex)
|
||||
// FieldInterpolator sinterp;
|
||||
/// Smoothing interpolator (one constraint per vertex)
|
||||
FieldInterpolator sinterp;
|
||||
|
||||
/// Defines the rendering style of the ballon
|
||||
RenderModes rm;
|
||||
|
||||
@ -10,7 +10,7 @@ using namespace vcg;
|
||||
// of a SQUARE matrix with the appropriate mesh coefficients. Constraints on vertices
|
||||
// are added by the function AddConstraint
|
||||
void test_cholmod2();
|
||||
bool FieldInterpolator::Init(CMeshO* mesh, LAPLACIAN laptype){
|
||||
bool FieldInterpolator::Init(CMeshO* mesh, int numVars, LAPLACIAN laptype){
|
||||
// Testing function
|
||||
// test_cholmod2(); exit(0);
|
||||
|
||||
@ -25,12 +25,25 @@ bool FieldInterpolator::Init(CMeshO* mesh, LAPLACIAN laptype){
|
||||
// Store mesh pointer & master OMEGA
|
||||
this->mesh = mesh;
|
||||
|
||||
// Allocate memory: square matrix to feed cholesky
|
||||
A_dyn = new DynamicSparseMatrix<FieldType>(mesh->vn,mesh->vn);
|
||||
xb = new Matrix<FieldType, Dynamic, 1>(mesh->vn);
|
||||
xb->setZero(mesh->vn);
|
||||
|
||||
// Fill in the entries of A_dyn and xb
|
||||
// Allocate memory
|
||||
if( A != NULL ) delete A;
|
||||
A = new AType(mesh->vn,mesh->vn);
|
||||
assert( numVars == 1 || numVars == 3 );
|
||||
if( numVars == 1 ){
|
||||
if( xb0 != NULL ) delete xb0;
|
||||
xb0 = new XBType(mesh->vn);
|
||||
xb0->setZero(mesh->vn);
|
||||
}
|
||||
else if( numVars == 3 ){
|
||||
if( xb0 != NULL ) delete xb0;
|
||||
xb0 = new XBType(mesh->vn); xb0->setZero(mesh->vn);
|
||||
if( xb1 != NULL ) delete xb1;
|
||||
xb1 = new XBType(mesh->vn); xb1->setZero(mesh->vn);
|
||||
if( xb2 != NULL ) delete xb2;
|
||||
xb2 = new XBType(mesh->vn); xb2->setZero(mesh->vn);
|
||||
}
|
||||
|
||||
// Fill in the entries of A and xb
|
||||
if( laptype == COMBINATORIAL ){
|
||||
// Iterate throgh every vertice
|
||||
CVertexO* v;
|
||||
@ -51,14 +64,14 @@ bool FieldInterpolator::Init(CMeshO* mesh, LAPLACIAN laptype){
|
||||
// Off-diagonal entries (symmetric)
|
||||
// Only lower diagonal
|
||||
if(v_i>vNeigh_i)
|
||||
A_dyn->coeffRef(v_i,vNeigh_i) = -1;
|
||||
A->coeffRef(v_i,vNeigh_i) = -1;
|
||||
else
|
||||
A_dyn->coeffRef(vNeigh_i,v_i) = -1;
|
||||
A->coeffRef(vNeigh_i,v_i) = -1;
|
||||
numNeighs++;
|
||||
}
|
||||
while(vNeigh != firstV);
|
||||
// Diagonal entries
|
||||
A_dyn->coeffRef(v_i,v_i) = numNeighs;
|
||||
A->coeffRef(v_i,v_i) = numNeighs;
|
||||
// Was I doing it right?
|
||||
assert(pos.NumberOfIncidentVertices() == numNeighs);
|
||||
}
|
||||
@ -132,16 +145,21 @@ bool FieldInterpolator::Init(CMeshO* mesh, LAPLACIAN laptype){
|
||||
totcoeff += cotcoeff;
|
||||
// Off-diagonal entries (symmetric) Only lower diagonal
|
||||
if(vI_i > vJ_i)
|
||||
A_dyn->coeffRef(vI_i,vJ_i) = -cotcoeff;
|
||||
A->coeffRef(vI_i,vJ_i) = -cotcoeff;
|
||||
else
|
||||
A_dyn->coeffRef(vJ_i,vI_i) = -cotcoeff;
|
||||
A->coeffRef(vJ_i,vI_i) = -cotcoeff;
|
||||
}
|
||||
while(vJ != firstV);
|
||||
// Diagonal entries
|
||||
A_dyn->coeffRef(vI_i,vI_i) = totcoeff;
|
||||
A->coeffRef(vI_i,vI_i) = totcoeff;
|
||||
}
|
||||
}
|
||||
|
||||
// Bi-laplacian
|
||||
#ifdef USEBILAP
|
||||
AType Atr = A->transpose();
|
||||
(*A) = Atr * (*A);
|
||||
#endif
|
||||
return true;
|
||||
}
|
||||
void FieldInterpolator::ColorizeIllConditioned(LAPLACIAN laptype){
|
||||
@ -185,21 +203,28 @@ void FieldInterpolator::ColorizeIllConditioned(LAPLACIAN laptype){
|
||||
}
|
||||
tri::UpdateColor<CMeshO>::VertexQualityRamp(*mesh);
|
||||
}
|
||||
void FieldInterpolator::AddConstraint( unsigned int vidx, FieldType omega, FieldType val){
|
||||
void FieldInterpolator::AddConstraint( unsigned int vidx, FieldType omega, FieldType val ){
|
||||
assert( !math::IsNAN(val) );
|
||||
assert( omega != 0 );
|
||||
A_dyn->coeffRef(vidx,vidx) += omega;
|
||||
xb->coeffRef(vidx) += omega*val;
|
||||
A->coeffRef(vidx,vidx) += omega;
|
||||
xb0->coeffRef(vidx) += omega*val;
|
||||
}
|
||||
void FieldInterpolator::Solve(){
|
||||
void FieldInterpolator::AddConstraint3( unsigned int vidx, FieldType omega, FieldType val0, FieldType val1, FieldType val2 ){
|
||||
assert( xb1 != NULL && xb2 != NULL );
|
||||
A->coeffRef(vidx,vidx) += omega;
|
||||
(*xb0)(vidx) += omega*val0;
|
||||
(*xb1)(vidx) += omega*val1;
|
||||
(*xb2)(vidx) += omega*val2;
|
||||
}
|
||||
void FieldInterpolator::SolveInQuality(){
|
||||
//--- DEBUG: Damp the matrix on file, so that matlab can read it and compute solution
|
||||
// You need to remove the first few lines from the matrix.txt before being able to import
|
||||
// the matrix successfully!!
|
||||
#if 0
|
||||
ofstream myfile;
|
||||
myfile.open ("/Users/ata2/workspace/workspace_vase/cholmod_verify_posdefi/matrix.txt");
|
||||
myfile << *A_dyn;
|
||||
cout << *A_dyn;
|
||||
myfile << *A;
|
||||
cout << *A;
|
||||
myfile.close();
|
||||
myfile.open("/Users/ata2/workspace/workspace_vase/cholmod_verify_posdefi/vector.txt");
|
||||
myfile << *xb;
|
||||
@ -208,14 +233,14 @@ void FieldInterpolator::Solve(){
|
||||
#endif
|
||||
|
||||
//--- Construct the cholesky factorization
|
||||
SparseLLT<SparseMatrix<FieldType> ,Cholmod> llt(*A_dyn);
|
||||
SparseLLT<SparseMatrix<FieldType> ,Cholmod> llt(*A);
|
||||
//--- Solve the linear system
|
||||
llt.solveInPlace(*xb);
|
||||
llt.solveInPlace(*xb0);
|
||||
//-- Copy the results back in Q() field of surface
|
||||
// cout << "X: ";
|
||||
// bool hasnan = false;
|
||||
for(unsigned int vi=0; vi<mesh->vert.size(); vi++){
|
||||
(mesh->vert[vi]).Q() = (*xb)[vi];
|
||||
(mesh->vert[vi]).Q() = (*xb0)[vi];
|
||||
// hasnan = hasnan || math::IsNAN( (*xb)[vi] );
|
||||
// cout << (*xb)[vi] << " ";
|
||||
}
|
||||
@ -228,11 +253,12 @@ void FieldInterpolator::Solve(){
|
||||
infile >> (mesh->vert[vi]).Q();
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
void FieldInterpolator::Solve( Matrix<FieldType, Dynamic, 1>* sols[3] ){
|
||||
SparseLLT<SparseMatrix<FieldType> ,Cholmod> llt(*A);
|
||||
llt.solveInPlace(*xb0); sols[0] = xb0;
|
||||
llt.solveInPlace(*xb1); sols[1] = xb1;
|
||||
llt.solveInPlace(*xb2); sols[2] = xb2;
|
||||
}
|
||||
|
||||
/*
|
||||
|
||||
@ -249,36 +275,36 @@ void test_cholmod2(){
|
||||
typedef double FieldType;
|
||||
|
||||
// These are how I defined them
|
||||
Eigen::DynamicSparseMatrix<FieldType>* A_dyn;
|
||||
A_dyn = new DynamicSparseMatrix<FieldType>(6,6);
|
||||
Eigen::DynamicSparseMatrix<FieldType>* A;
|
||||
A = new DynamicSparseMatrix<FieldType>(6,6);
|
||||
// Equivalent to VectorXd or VectorXf (http://eigen.tuxfamily.org/dox/TutorialCore.html#TutorialCoreMatrixTypes)
|
||||
Eigen::Matrix<FieldType, Eigen::Dynamic, 1>* xb;
|
||||
xb = new Matrix<FieldType, Dynamic, 1>(6);
|
||||
xb->setZero(6);
|
||||
|
||||
// MOD
|
||||
A.insert(0,0) = 4; A_dyn->coeffRef(0,0) = 4;
|
||||
A.insert(1,0) = -1; A_dyn->coeffRef(1,0) = -1;
|
||||
A.insert(2,0) = -1; A_dyn->coeffRef(2,0) = -1;
|
||||
A.insert(4,0) = -1; A_dyn->coeffRef(4,0) = -1;
|
||||
A.insert(5,0) = -1; A_dyn->coeffRef(5,0) = -1;
|
||||
A.insert(1,1) = 104; A_dyn->coeffRef(1,1) = 104;
|
||||
A.insert(2,1) = -1; A_dyn->coeffRef(2,1) = -1;
|
||||
A.insert(3,1) = -1; A_dyn->coeffRef(3,1) = -1;
|
||||
A.insert(4,1) = -1; A_dyn->coeffRef(4,1) = -1;
|
||||
A.insert(2,2) = 4; A_dyn->coeffRef(2,2) = 4;
|
||||
A.insert(3,2) = -1; A_dyn->coeffRef(3,2) = -1;
|
||||
A.insert(5,2) = -1; A_dyn->coeffRef(5,2) = -1;
|
||||
A.insert(3,3) = 4; A_dyn->coeffRef(3,3) = 4;
|
||||
A.insert(4,3) = -1; A_dyn->coeffRef(4,3) = -1;
|
||||
A.insert(5,3) = -1; A_dyn->coeffRef(5,3) = -1;
|
||||
A.insert(4,4) = 4; A_dyn->coeffRef(4,4) = 4;
|
||||
A.insert(5,4) = -1; A_dyn->coeffRef(5,4) = -1;
|
||||
A.insert(5,5) = 104; A_dyn->coeffRef(5,5) = 104;
|
||||
A.insert(0,0) = 4; A->coeffRef(0,0) = 4;
|
||||
A.insert(1,0) = -1; A->coeffRef(1,0) = -1;
|
||||
A.insert(2,0) = -1; A->coeffRef(2,0) = -1;
|
||||
A.insert(4,0) = -1; A->coeffRef(4,0) = -1;
|
||||
A.insert(5,0) = -1; A->coeffRef(5,0) = -1;
|
||||
A.insert(1,1) = 104; A->coeffRef(1,1) = 104;
|
||||
A.insert(2,1) = -1; A->coeffRef(2,1) = -1;
|
||||
A.insert(3,1) = -1; A->coeffRef(3,1) = -1;
|
||||
A.insert(4,1) = -1; A->coeffRef(4,1) = -1;
|
||||
A.insert(2,2) = 4; A->coeffRef(2,2) = 4;
|
||||
A.insert(3,2) = -1; A->coeffRef(3,2) = -1;
|
||||
A.insert(5,2) = -1; A->coeffRef(5,2) = -1;
|
||||
A.insert(3,3) = 4; A->coeffRef(3,3) = 4;
|
||||
A.insert(4,3) = -1; A->coeffRef(4,3) = -1;
|
||||
A.insert(5,3) = -1; A->coeffRef(5,3) = -1;
|
||||
A.insert(4,4) = 4; A->coeffRef(4,4) = 4;
|
||||
A.insert(5,4) = -1; A->coeffRef(5,4) = -1;
|
||||
A.insert(5,5) = 104; A->coeffRef(5,5) = 104;
|
||||
A.finalize();
|
||||
|
||||
// MOD
|
||||
std::cerr << MatrixXd(*A_dyn) << "\n\n";
|
||||
std::cerr << MatrixXd(*A) << "\n\n";
|
||||
|
||||
// Create a full copy of the matrix
|
||||
MatrixXd AFULL(A);
|
||||
@ -306,11 +332,11 @@ void test_cholmod2(){
|
||||
std::cout << "SPAR: " << xs.transpose() << endl;
|
||||
|
||||
// Using my way of declaring
|
||||
SparseLLT<SparseMatrix<FieldType> ,Cholmod> llt2(*A_dyn);
|
||||
SparseLLT<SparseMatrix<FieldType> ,Cholmod> llt2(*A);
|
||||
llt2.solveInPlace(*xb);
|
||||
std::cout << "SPAR2: " << (*xb).transpose() << endl;
|
||||
|
||||
delete A_dyn;
|
||||
delete A;
|
||||
delete xb;
|
||||
}
|
||||
|
||||
|
||||
@ -22,38 +22,55 @@ namespace vcg{
|
||||
/// The type of laplacian weights the solver will use
|
||||
enum LAPLACIAN {COMBINATORIAL, COTANGENT};
|
||||
|
||||
/// CHOLMOD can only work on the DOUBLE type. NaN otherwise
|
||||
typedef double FieldType;
|
||||
|
||||
/// Interpolates a scalar field on a surface according to laplacian criterions
|
||||
class FieldInterpolator{
|
||||
public:
|
||||
/// CHOLMOD can only work on the DOUBLE type. NaN otherwise
|
||||
typedef double FieldType;
|
||||
typedef Eigen::DynamicSparseMatrix<FieldType> AType;
|
||||
typedef Eigen::Matrix<FieldType, Eigen::Dynamic, 1> XBType;
|
||||
|
||||
private:
|
||||
/// The domain over which the interpolation is defined
|
||||
CMeshO* mesh;
|
||||
/// The EIGEN matrices/vectors involved in the interpolation
|
||||
Eigen::DynamicSparseMatrix<FieldType>* A_dyn;
|
||||
Eigen::Matrix<FieldType, Eigen::Dynamic, 1>* xb;
|
||||
/// The EIGEN "A" square matrix to be involved in the interpolation
|
||||
AType* A;
|
||||
/// The EIGEN "b" coefficient vector/matrix to be involved
|
||||
XBType* xb0;
|
||||
XBType* xb1;
|
||||
XBType* xb2;
|
||||
|
||||
public:
|
||||
/// Refer to Init
|
||||
FieldInterpolator(){
|
||||
A_dyn = 0;
|
||||
xb = 0;
|
||||
A = NULL;
|
||||
xb0 = NULL;
|
||||
xb1 = NULL;
|
||||
xb2 = NULL;
|
||||
}
|
||||
/// Free the memory stored for sparse matrixes
|
||||
~FieldInterpolator(){
|
||||
if( A_dyn!=0 ) delete A_dyn;
|
||||
if( xb!=0 ) delete xb;
|
||||
if( A != NULL ) delete A;
|
||||
if( xb0 != NULL ) delete xb0;
|
||||
if( xb1 != NULL ) delete xb1;
|
||||
if( xb2 != NULL ) delete xb2;
|
||||
}
|
||||
/// Creates the sparse laplacian matrixes used by the solver using the specified on the domain of mesh
|
||||
/// @param numVars number of solves to be executed on the same domain/factorization
|
||||
/// @return false if initialization has failed for bad-conditioned triangles
|
||||
bool Init(CMeshO* mesh, LAPLACIAN laptype=COMBINATORIAL);
|
||||
bool Init(CMeshO* mesh, int numVars=1, LAPLACIAN laptype=COMBINATORIAL);
|
||||
/// Assigns color to vertexes and faces highlighting ill-conditioned elements according to the selected laplacian
|
||||
void ColorizeIllConditioned(LAPLACIAN laptype);
|
||||
/// Adds a constraint to vertex vidx
|
||||
void AddConstraint( unsigned int vidx, FieldType omega, FieldType val);
|
||||
/// Adds a constraint to vertex vidx for the pos-th set of data to be solved
|
||||
void AddConstraint3( unsigned int vidx, FieldType omega, FieldType val0, FieldType val1, FieldType val2 );
|
||||
|
||||
/// Solves the interpolation problem by Cholesky returning the solution in the Q() vertex property
|
||||
void Solve();
|
||||
void SolveInQuality();
|
||||
/// Solves the interpolation problem by Cholesky returning the solution in passed vector
|
||||
/// The returned solution vector is valid up to the next time Init is called or the object is destroyed
|
||||
void Solve( XBType* sols[3] );
|
||||
};
|
||||
|
||||
} // vcg::
|
||||
|
||||
@ -105,7 +105,7 @@ void VaseWidget::on_slice_offset_sliderMoved(int){
|
||||
void VaseWidget::on_iterationButton_released(){
|
||||
bool op_succeed = false;
|
||||
for(int i=0; i<ui.numItersSpin->value(); i++){
|
||||
gla->log->Logf(GLLogStream::FILTER, "\n----- began iteration %d -----", balloon->numiterscompleted);
|
||||
qDebug("\n----- began iteration %d -----", balloon->numiterscompleted);
|
||||
op_succeed = balloon->initializeField(); if( !op_succeed ) break;
|
||||
balloon->interpolateField(); if( !op_succeed ) break;
|
||||
balloon->computeCurvature(); if( !op_succeed ) break;
|
||||
|
||||
@ -147,7 +147,7 @@ void MyVolume::initField( const vcg::Box3f& inbbox ){
|
||||
else
|
||||
grid.Val(i,j,k) = +1;
|
||||
}
|
||||
else if( false){
|
||||
else if( true ){
|
||||
// Radius in object space
|
||||
double r = 1;
|
||||
Point3f p;
|
||||
@ -241,8 +241,7 @@ void MyVolume::isosurface( CMeshO& mesh, float offset ){
|
||||
mesh.vert.EnableVFAdjacency();
|
||||
mesh.face.EnableVFAdjacency();
|
||||
tri::UpdateTopology<CMeshO>::VertexFace( mesh );
|
||||
//tri::MCSimplify( mesh, getDelta()/4 );
|
||||
tri::MCSimplify( mesh, 0.0);
|
||||
tri::Simplify( mesh, getDelta()/10 ); // was 4 then 20
|
||||
|
||||
//--- The simplify operation removed some vertices
|
||||
tri::Allocator<CMeshO>::CompactVertexVector( mesh );
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user