New version of the Mixed Integer Quadrangulation framework that relies only on IGL and comiso
This commit is contained in:
parent
fb2ad3f323
commit
17f1f52bae
|
|
@ -0,0 +1,327 @@
|
|||
/****************************************************************************
|
||||
* VCGLib o o *
|
||||
* Visual and Computer Graphics Library o o *
|
||||
* _ O _ *
|
||||
* Copyright(C) 2014 \/)\/ *
|
||||
* Visual Computing Lab /\/| *
|
||||
* ISTI - Italian National Research Council | *
|
||||
* \ *
|
||||
* All rights reserved. *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
|
||||
#ifndef __MIQ_PARAMETRIZATION_H
|
||||
#define __MIQ_PARAMETRIZATION_H
|
||||
|
||||
//igl stuff
|
||||
#include <igl/cross_field_missmatch.h>
|
||||
#include <igl/line_field_missmatch.h>
|
||||
#include <igl/comb_line_field.h>
|
||||
#include <igl/cut_mesh_from_singularities.h>
|
||||
#include <igl/find_cross_field_singularities.h>
|
||||
#include <igl/compute_frame_field_bisectors.h>
|
||||
#include <igl/comiso/miq.h>
|
||||
#include <vcg/complex/algorithms/parametrization/uv_utils.h>
|
||||
#include <vcg/complex/algorithms/mesh_to_matrix.h>
|
||||
|
||||
namespace vcg {
|
||||
namespace tri {
|
||||
template < class MeshType> // Classe templatata su Tipo Mesh
|
||||
class MiQParametrizer
|
||||
{
|
||||
|
||||
typedef typename MeshType::CoordType CoordType;
|
||||
typedef typename MeshType::VertexType VertexType;
|
||||
typedef typename MeshType::FaceType FaceType;
|
||||
typedef typename MeshType::ScalarType ScalarType;
|
||||
typedef typename MeshType::VertexType PolyVertexType;
|
||||
|
||||
public:
|
||||
struct MIQParameters
|
||||
{
|
||||
//the gradient of the parametrization 1 is the bb diagonal, as big is the gradient as small are the quads
|
||||
double gradient;
|
||||
//do the rounding or not across cuts... set to true to get a quadrangulation
|
||||
bool doRound;
|
||||
//do the round at once for each stiffness iteration or do it gradually.. gradually is more stable but much mor slow
|
||||
bool directRound;
|
||||
//the stiffness increment for ach iteration
|
||||
double stiffness;
|
||||
//the maximum number ofstiffness iteration to avoid folds
|
||||
int stiffness_iter;
|
||||
//local iteration to round integer variables ofr each stiffness round
|
||||
int local_iter;
|
||||
//this bool multiply the gradiens U or V separately times a constant to make hexagons
|
||||
bool hexaLine;
|
||||
//the threshold of normal dot product to consider a crease
|
||||
double crease_thr;
|
||||
//number of 90° rotation independence (4 for quad meshing, 2 to obtain a quad meshing for hexagonalization)
|
||||
int Ndir;
|
||||
//round or not the singularities
|
||||
bool round_singularities;
|
||||
//use the crease edges as feature or not
|
||||
bool crease_as_feature;
|
||||
|
||||
MIQParameters()
|
||||
{
|
||||
gradient=80;
|
||||
doRound=true;
|
||||
directRound=true;
|
||||
round_singularities=true;
|
||||
crease_as_feature=false;
|
||||
stiffness=5;
|
||||
stiffness_iter=10;
|
||||
local_iter=5;
|
||||
Ndir=4;
|
||||
crease_thr=0.2;
|
||||
hexaLine=false;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
static void GetFeatureLines(MeshType &trimesh,
|
||||
std::vector<std::vector<int> > &feature_lines)
|
||||
{
|
||||
feature_lines.clear();
|
||||
|
||||
for (size_t i=0;i<trimesh.face.size();i++)
|
||||
{
|
||||
for (int j=0;j<3;j++)
|
||||
{
|
||||
//if (!trimesh.face[i].IsB(j))continue;
|
||||
if (!trimesh.face[i].IsCrease(j))continue;
|
||||
|
||||
feature_lines.push_back(std::vector<int>());
|
||||
feature_lines.back().push_back(i);
|
||||
feature_lines.back().push_back(j);
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
private:
|
||||
|
||||
|
||||
|
||||
static void CrossFieldParam(MeshType &trimesh,
|
||||
MIQParameters &MiqP)
|
||||
{
|
||||
|
||||
Eigen::MatrixXi F;
|
||||
Eigen::MatrixXd V;
|
||||
vcg::tri::MeshToMatrix<MeshType>::GetTriMeshData(trimesh,F,V);
|
||||
|
||||
//then get the principal directions
|
||||
Eigen::MatrixXd X1,X2;
|
||||
X1=Eigen::MatrixXd(trimesh.FN(), 3);
|
||||
for (size_t i=0;i<trimesh.face.size();i++)
|
||||
{
|
||||
CoordType Dir1=trimesh.face[i].PD1();
|
||||
Dir1.Normalize();
|
||||
for (int j=0;j<3;j++)
|
||||
{
|
||||
X1(i,j)=Dir1[j];
|
||||
}
|
||||
}
|
||||
|
||||
Eigen::MatrixXd B1,B2,B3;
|
||||
igl::local_basis(V,F,B1,B2,B3);
|
||||
X2 = igl::rotate_vectors(X1, Eigen::VectorXd::Constant(1,M_PI/2), B1, B2);
|
||||
|
||||
Eigen::MatrixXd UV;
|
||||
Eigen::MatrixXi FUV;
|
||||
|
||||
//ScalarType gradsize=trimesh.bbox.Diag()*2;
|
||||
|
||||
std::vector<std::vector<int> > hard_features;
|
||||
|
||||
if (MiqP.crease_as_feature)
|
||||
GetFeatureLines(trimesh,hard_features);
|
||||
|
||||
std::vector<int> extra_round;
|
||||
|
||||
igl::miq(V,F,X1,X2,UV,FUV,MiqP.gradient,MiqP.stiffness,MiqP.directRound,
|
||||
MiqP.stiffness_iter,MiqP.local_iter,MiqP.doRound,MiqP.round_singularities,
|
||||
extra_round,hard_features);
|
||||
|
||||
// then copy UV
|
||||
for (size_t i=0;i<trimesh.face.size();i++)
|
||||
{
|
||||
for (int j=0;j<3;j++)
|
||||
{
|
||||
int index=FUV(i,j);
|
||||
trimesh.face[i].WT(j).P()[0]=UV(index,0);
|
||||
trimesh.face[i].WT(j).P()[1]=UV(index,1);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
static void LineFieldParam(MeshType &trimesh,
|
||||
MIQParameters &MiqP)
|
||||
{
|
||||
|
||||
Eigen::MatrixXi F;
|
||||
Eigen::MatrixXd V;
|
||||
vcg::tri::MeshToMatrix<MeshType>::GetTriMeshData(trimesh,F,V);
|
||||
|
||||
//then get the principal directions
|
||||
Eigen::MatrixXd X1,X2;
|
||||
X1=Eigen::MatrixXd(trimesh.FN(), 3);
|
||||
for (size_t i=0;i<trimesh.face.size();i++)
|
||||
{
|
||||
CoordType Dir1=trimesh.face[i].PD1();
|
||||
Dir1.Normalize();
|
||||
for (int j=0;j<3;j++)
|
||||
{
|
||||
X1(i,j)=Dir1[j];
|
||||
}
|
||||
}
|
||||
|
||||
Eigen::MatrixXd B1,B2,B3;
|
||||
igl::local_basis(V,F,B1,B2,B3);
|
||||
X2 = igl::rotate_vectors(X1, Eigen::VectorXd::Constant(1,M_PI/2), B1, B2);
|
||||
|
||||
// Bisector field
|
||||
Eigen::MatrixXd BIS1, BIS2;
|
||||
|
||||
// Combed bisector
|
||||
Eigen::MatrixXd BIS1_combed, BIS2_combed;
|
||||
|
||||
// Per-corner, integer mismatches
|
||||
Eigen::MatrixXi MMatch;
|
||||
|
||||
// Field singularities
|
||||
Eigen::VectorXi isSingularity, singularityIndex;
|
||||
|
||||
// Per corner seams
|
||||
Eigen::MatrixXi Seams;
|
||||
|
||||
// Combed field
|
||||
Eigen::MatrixXd X1_combed, X2_combed;
|
||||
|
||||
|
||||
// Global parametrization (with seams)
|
||||
Eigen::MatrixXd UV_seams;
|
||||
Eigen::MatrixXi FUV_seams;
|
||||
|
||||
// Global parametrization
|
||||
Eigen::MatrixXd UV;
|
||||
Eigen::MatrixXi FUV;
|
||||
|
||||
// Always work on the bisectors, it is more general
|
||||
igl::compute_frame_field_bisectors(V, F, X1, X2, BIS1, BIS2);
|
||||
|
||||
// Comb the field, implicitly defining the seams
|
||||
igl::comb_line_field(V, F, BIS1, BIS1_combed);
|
||||
igl::local_basis(V,F,B1,B2,B3);
|
||||
BIS2_combed = igl::rotate_vectors(BIS1_combed, Eigen::VectorXd::Constant(1,M_PI/2), B1, B2);
|
||||
|
||||
// Find the integer mismatches
|
||||
igl::line_field_missmatch(V, F, BIS1_combed, true, MMatch);
|
||||
|
||||
// Find the singularities
|
||||
igl::find_cross_field_singularities(V, F, MMatch, isSingularity, singularityIndex);
|
||||
|
||||
// Cut the mesh, duplicating all vertices on the seams
|
||||
igl::cut_mesh_from_singularities(V, F, MMatch, isSingularity, singularityIndex, Seams);
|
||||
|
||||
// Comb the frame-field accordingly
|
||||
igl::comb_frame_field(V, F, X1, X2, BIS1_combed, BIS2_combed, X1_combed, X2_combed);
|
||||
|
||||
std::vector<std::vector<int> > hard_features;
|
||||
|
||||
|
||||
std::vector<int> extra_round;
|
||||
|
||||
if (MiqP.crease_as_feature)
|
||||
GetFeatureLines(trimesh,hard_features);
|
||||
|
||||
//scale gradient if needed
|
||||
ScalarType sqrt3=1.732050807568877;
|
||||
|
||||
ScalarType GradX=0.5;
|
||||
ScalarType GradY=1;
|
||||
|
||||
if (MiqP.hexaLine)
|
||||
{
|
||||
for (int i=0;i<X1_combed.rows();i++)
|
||||
{
|
||||
X1_combed(i)*=GradX;//*ScaleFact;
|
||||
X2_combed(i)*=GradY;//*ScaleFact;
|
||||
}
|
||||
}
|
||||
|
||||
igl::miq(V,F,X1_combed,X2_combed,BIS1_combed,BIS2_combed,
|
||||
MMatch,isSingularity,singularityIndex,Seams,
|
||||
UV,FUV,MiqP.gradient,MiqP.stiffness,MiqP.directRound,
|
||||
MiqP.stiffness_iter,MiqP.local_iter,MiqP.doRound,MiqP.round_singularities,extra_round,hard_features);
|
||||
|
||||
// then copy UV
|
||||
for (size_t i=0;i<trimesh.face.size();i++)
|
||||
{
|
||||
for (int j=0;j<3;j++)
|
||||
{
|
||||
int index=FUV(i,j);
|
||||
trimesh.face[i].WT(j).P()[0]=UV(index,0);//*4;
|
||||
trimesh.face[i].WT(j).P()[1]=UV(index,1);//*2;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
public:
|
||||
|
||||
static void SetCreases(MeshType & mesh,
|
||||
const ScalarType &thr=0.2,
|
||||
bool setBorder=true)
|
||||
{
|
||||
for (size_t i=0;i<mesh.face.size();i++)
|
||||
for (int j=0;j<mesh.face[i].VN();j++)
|
||||
{
|
||||
FaceType *f0=&mesh.face[i];
|
||||
f0->ClearCrease(j);
|
||||
}
|
||||
|
||||
|
||||
for (size_t i=0;i<mesh.face.size();i++)
|
||||
for (int j=0;j<mesh.face[i].VN();j++)
|
||||
{
|
||||
FaceType *f0=&mesh.face[i];
|
||||
FaceType *f1=f0->FFp(j);
|
||||
|
||||
if (f0==f1){f0->SetCrease(j);continue;}
|
||||
|
||||
CoordType N0=f0->N();
|
||||
CoordType N1=f1->N();
|
||||
if ((N0*N1)>thr)continue;
|
||||
f0->SetCrease(j);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
static void MIQParametrize(MeshType &trimesh,
|
||||
MIQParameters &MiqP)
|
||||
{
|
||||
if (MiqP.crease_as_feature)
|
||||
SetCreases(trimesh,MiqP.crease_thr);
|
||||
|
||||
if (MiqP.Ndir==4)
|
||||
CrossFieldParam(trimesh,MiqP);
|
||||
else
|
||||
LineFieldParam(trimesh,MiqP);
|
||||
}
|
||||
};
|
||||
|
||||
} // end namespace tri
|
||||
} // end namespace vcg
|
||||
#endif // VORO_CLUSTER_H
|
||||
|
|
@ -0,0 +1,487 @@
|
|||
/****************************************************************************
|
||||
* VCGLib o o *
|
||||
* Visual and Computer Graphics Library o o *
|
||||
* _ O _ *
|
||||
* Copyright(C) 2014 \/)\/ *
|
||||
* Visual Computing Lab /\/| *
|
||||
* ISTI - Italian National Research Council | *
|
||||
* \ *
|
||||
* All rights reserved. *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
|
||||
#ifndef SMOOTHER_FIELD_H
|
||||
#define SMOOTHER_FIELD_H
|
||||
|
||||
//eigen stuff
|
||||
#include <eigenlib/Eigen/Sparse>
|
||||
|
||||
//vcg stuff
|
||||
#include <vcg/complex/algorithms/update/color.h>
|
||||
#include <vcg/complex/algorithms/update/quality.h>
|
||||
#include <vcg/complex/algorithms/parametrization/tangent_field_operators.h>
|
||||
#include <vcg/complex/algorithms/mesh_to_matrix.h>
|
||||
|
||||
//igl related stuff
|
||||
#include <igl/n_polyvector.h>
|
||||
#include <igl/principal_curvature.h>
|
||||
#include <igl/igl_inline.h>
|
||||
#include <igl/comiso/nrosy.h>
|
||||
|
||||
namespace vcg {
|
||||
namespace tri {
|
||||
|
||||
enum SmoothMethod{SMMiq,SMNPoly};
|
||||
|
||||
template <class MeshType>
|
||||
class FieldSmoother
|
||||
{
|
||||
typedef typename MeshType::FaceType FaceType;
|
||||
typedef typename MeshType::VertexType VertexType;
|
||||
typedef typename MeshType::ScalarType ScalarType;
|
||||
typedef typename MeshType::CoordType CoordType;
|
||||
|
||||
|
||||
static void InitQualityByAnisotropyDir(MeshType &mesh)
|
||||
{
|
||||
std::vector<ScalarType> QVal;
|
||||
for (size_t i=0;i<mesh.vert.size();i++)
|
||||
{
|
||||
ScalarType N1=fabs(mesh.vert[i].K1());
|
||||
ScalarType N2=fabs(mesh.vert[i].K2());
|
||||
QVal.push_back(N1);
|
||||
QVal.push_back(N2);
|
||||
}
|
||||
|
||||
std::sort(QVal.begin(),QVal.end());
|
||||
int percUp=int(floor(QVal.size()*0.95+0.5));
|
||||
ScalarType trimUp=QVal[percUp];
|
||||
|
||||
for (size_t i=0;i<mesh.vert.size();i++)
|
||||
{
|
||||
ScalarType N1=(mesh.vert[i].K1());
|
||||
ScalarType N2=(mesh.vert[i].K2());
|
||||
|
||||
ScalarType NMax=std::max(N1,N2)/trimUp;
|
||||
ScalarType NMin=std::min(N1,N2)/trimUp;
|
||||
|
||||
if (NMax>1)NMax=1;
|
||||
if (NMax<-1)NMax=-1;
|
||||
|
||||
if (NMin>1)NMin=1;
|
||||
if (NMin<-1)NMin=-1;
|
||||
|
||||
ScalarType CurvAni=(NMax-NMin)/2;
|
||||
mesh.vert[i].Q()=CurvAni;
|
||||
}
|
||||
vcg::tri::UpdateQuality<MeshType>::FaceFromVertex(mesh);
|
||||
}
|
||||
|
||||
static void SetEdgeDirection(FaceType *f,int edge)
|
||||
{
|
||||
CoordType dir=f->P0(edge)-f->P1(edge);
|
||||
dir.Normalize();
|
||||
ScalarType prod1=fabs(dir*f->PD1());
|
||||
ScalarType prod2=fabs(dir*f->PD2());
|
||||
if (prod1>prod2)
|
||||
{
|
||||
f->PD1()=dir;
|
||||
f->PD2()=f->N()^dir;
|
||||
}else
|
||||
{
|
||||
f->PD2()=dir;
|
||||
f->PD1()=f->N()^dir;
|
||||
}
|
||||
}
|
||||
|
||||
static void AddSharpEdgesConstraints(MeshType & mesh,
|
||||
const ScalarType &thr=0.2)
|
||||
{
|
||||
for (size_t i=0;i<mesh.face.size();i++)
|
||||
for (int j=0;j<mesh.face[i].VN();j++)
|
||||
{
|
||||
FaceType *f0=&mesh.face[i];
|
||||
FaceType *f1=f0->FFp(j);
|
||||
if (f0==f1)continue;
|
||||
CoordType N0=f0->N();
|
||||
CoordType N1=f1->N();
|
||||
if ((N0*N1)>thr)continue;
|
||||
SetEdgeDirection(f0,j);
|
||||
f0->SetS();
|
||||
}
|
||||
}
|
||||
|
||||
static void AddBorderConstraints(MeshType & mesh)
|
||||
{
|
||||
for (size_t i=0;i<mesh.face.size();i++)
|
||||
for (int j=0;j<mesh.face[i].VN();j++)
|
||||
{
|
||||
FaceType *f0=&mesh.face[i];
|
||||
FaceType *f1=f0->FFp(j);
|
||||
assert(f1!=NULL);
|
||||
if (f0!=f1)continue;
|
||||
SetEdgeDirection(f0,j);
|
||||
f0->SetS();
|
||||
}
|
||||
}
|
||||
|
||||
static void AddCurvatureConstraints(MeshType & mesh,const ScalarType &thr=0.5)
|
||||
{
|
||||
for (size_t i=0;i<mesh.face.size();i++)
|
||||
if (mesh.face[i].Q()>thr)mesh.face[i].SetS();
|
||||
}
|
||||
|
||||
//hard constraints have selected face
|
||||
static void CollectHardConstraints( MeshType & mesh,
|
||||
Eigen::VectorXi &HardI,
|
||||
Eigen::MatrixXd &HardD,
|
||||
SmoothMethod SMethod,
|
||||
int Ndir)
|
||||
{
|
||||
//count number of hard constraints
|
||||
int numS=vcg::tri::UpdateSelection<MeshType>::FaceCount(mesh);
|
||||
HardI=Eigen::MatrixXi(numS,1);
|
||||
if ((Ndir==2)||(SMethod==SMMiq))
|
||||
HardD=Eigen::MatrixXd(numS,3);
|
||||
else
|
||||
HardD=Eigen::MatrixXd(numS,6);
|
||||
//then update them
|
||||
int curr_index=0;
|
||||
for (size_t i=0;i<mesh.face.size();i++)
|
||||
{
|
||||
if (!mesh.face[i].IsS())continue;
|
||||
|
||||
CoordType dir=mesh.face[i].PD1();
|
||||
dir.Normalize();
|
||||
|
||||
HardI(curr_index,0)=i;
|
||||
|
||||
HardD(curr_index,0)=dir.X();
|
||||
HardD(curr_index,1)=dir.Y();
|
||||
HardD(curr_index,2)=dir.Z();
|
||||
if ((Ndir==4)&&(SMethod==SMNPoly))
|
||||
{
|
||||
dir=mesh.face[i].PD2();
|
||||
HardD(curr_index,3)=dir.X();
|
||||
HardD(curr_index,4)=dir.Y();
|
||||
HardD(curr_index,5)=dir.Z();
|
||||
}
|
||||
curr_index++;
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
//hard constraints have selected face
|
||||
static void CollectSoftConstraints( MeshType & mesh,
|
||||
Eigen::VectorXi &SoftI,
|
||||
Eigen::MatrixXd &SoftD,
|
||||
Eigen::VectorXd &SoftW)
|
||||
{
|
||||
//count number of soft constraints
|
||||
int numS=vcg::tri::UpdateSelection<MeshType>::FaceCount(mesh);
|
||||
numS=mesh.fn-numS;
|
||||
//allocate eigen matrix
|
||||
SoftI=Eigen::MatrixXi(numS,1);
|
||||
SoftD=Eigen::MatrixXd(numS,3);
|
||||
SoftW=Eigen::MatrixXd(numS,1);
|
||||
|
||||
//then update them
|
||||
int curr_index=0;
|
||||
for (size_t i=0;i<mesh.face.size();i++)
|
||||
{
|
||||
if (mesh.face[i].IsS())continue;
|
||||
|
||||
CoordType dir=mesh.face[i].PD1();
|
||||
dir.Normalize();
|
||||
|
||||
SoftI(curr_index,0)=i;
|
||||
|
||||
SoftD(curr_index,0)=dir.X();
|
||||
SoftD(curr_index,1)=dir.Y();
|
||||
SoftD(curr_index,2)=dir.Z();
|
||||
|
||||
SoftW(curr_index,0)=mesh.face[i].Q();
|
||||
|
||||
curr_index++;
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
static void SmoothMIQ(MeshType &mesh,
|
||||
Eigen::VectorXi &HardI, //hard constraints index
|
||||
Eigen::MatrixXd &HardD, //hard directions
|
||||
Eigen::VectorXi &SoftI, //soft constraints
|
||||
Eigen::MatrixXd &SoftD, //weight of soft constraints
|
||||
Eigen::VectorXd &SoftW, //soft directions
|
||||
ScalarType alpha_soft,
|
||||
int Ndir)
|
||||
{
|
||||
|
||||
assert((Ndir==2)||(Ndir==4));
|
||||
Eigen::MatrixXi F;
|
||||
Eigen::MatrixXd V;
|
||||
|
||||
MeshToMatrix<MeshType>::GetTriMeshData(mesh,F,V);
|
||||
|
||||
Eigen::MatrixXd output_field;
|
||||
Eigen::VectorXd output_sing;
|
||||
|
||||
igl::nrosy(V,F,HardI,HardD,SoftI,SoftW,SoftD,Ndir,alpha_soft,output_field,output_sing);
|
||||
|
||||
//finally update the principal directions
|
||||
for (size_t i=0;i<mesh.face.size();i++)
|
||||
{
|
||||
CoordType dir1;
|
||||
dir1[0]=output_field(i,0);
|
||||
dir1[1]=output_field(i,1);
|
||||
dir1[2]=output_field(i,2);
|
||||
|
||||
dir1.Normalize();
|
||||
CoordType dir2=mesh.face[i].N()^dir1;
|
||||
dir2.Normalize();
|
||||
|
||||
ScalarType Norm1=mesh.face[i].PD1().Norm();
|
||||
ScalarType Norm2=mesh.face[i].PD2().Norm();
|
||||
|
||||
mesh.face[i].PD1()=dir1*Norm1;
|
||||
mesh.face[i].PD2()=dir2*Norm2;
|
||||
}
|
||||
}
|
||||
|
||||
static void SmoothNPoly(MeshType &mesh,
|
||||
Eigen::VectorXi &HardI, //hard constraints index
|
||||
Eigen::MatrixXd &HardD, //hard directions
|
||||
int Ndir)
|
||||
{
|
||||
assert((Ndir==2)||(Ndir==4));
|
||||
|
||||
Eigen::MatrixXi F;
|
||||
Eigen::MatrixXd V;
|
||||
|
||||
MeshToMatrix<MeshType>::GetTriMeshData(mesh,F,V);
|
||||
|
||||
Eigen::MatrixXd output_field;
|
||||
Eigen::VectorXd output_sing;
|
||||
|
||||
igl::n_polyvector(V,F,HardI,HardD,output_field);
|
||||
|
||||
//finally update the principal directions
|
||||
for (size_t i=0;i<mesh.face.size();i++)
|
||||
{
|
||||
CoordType dir1;
|
||||
dir1[0]=output_field(i,0);
|
||||
dir1[1]=output_field(i,1);
|
||||
dir1[2]=output_field(i,2);
|
||||
|
||||
dir1.Normalize();
|
||||
CoordType dir2=mesh.face[i].N()^dir1;
|
||||
dir2.Normalize();
|
||||
|
||||
ScalarType Norm1=mesh.face[i].PD1().Norm();
|
||||
ScalarType Norm2=mesh.face[i].PD2().Norm();
|
||||
|
||||
mesh.face[i].PD1()=dir1*Norm1;
|
||||
mesh.face[i].PD2()=dir2*Norm2;
|
||||
}
|
||||
}
|
||||
|
||||
static void PickRandomDir(MeshType &mesh,
|
||||
int &indexF,
|
||||
CoordType &Dir)
|
||||
{
|
||||
indexF=rand()%mesh.fn;
|
||||
FaceType *currF=&mesh.face[indexF];
|
||||
CoordType dirN=currF->N();
|
||||
dirN.Normalize();
|
||||
Dir=CoordType(1,0,0);
|
||||
if (fabs(Dir*dirN)>0.9)
|
||||
Dir=CoordType(0,1,0);
|
||||
if (fabs(Dir*dirN)>0.9)
|
||||
Dir=CoordType(0,0,1);
|
||||
|
||||
Dir=dirN^Dir;
|
||||
Dir.Normalize();
|
||||
}
|
||||
|
||||
public:
|
||||
|
||||
struct SmoothParam
|
||||
{
|
||||
//the 90° rotation independence while smoothing the direction field
|
||||
int Ndir;
|
||||
//the weight of curvature if doing the smoothing keeping the field close to the original one
|
||||
ScalarType alpha_curv;
|
||||
//align the field to border or not
|
||||
bool align_borders;
|
||||
//threshold to consider some edge as sharp feature and to use as hard constraint (0, not use)
|
||||
ScalarType sharp_thr;
|
||||
//threshold to consider some edge as high curvature anisotropyand to use as hard constraint (0, not use)
|
||||
ScalarType curv_thr;
|
||||
//the method used to smooth MIQ or "Designing N-PolyVector Fields with Complex Polynomials"
|
||||
SmoothMethod SmoothM;
|
||||
//the number of faces of the ring used ot esteem the curvature
|
||||
int curvRing;
|
||||
|
||||
SmoothParam()
|
||||
{
|
||||
Ndir=4;
|
||||
curvRing=2;
|
||||
alpha_curv=0.0;
|
||||
|
||||
align_borders=false;
|
||||
|
||||
SmoothM=SMMiq;
|
||||
|
||||
sharp_thr=0.0;
|
||||
curv_thr=0.4;
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
static void SelectConstraints(MeshType &mesh,SmoothParam &SParam)
|
||||
{
|
||||
//clear all selected faces
|
||||
vcg::tri::UpdateFlags<MeshType>::FaceClearS(mesh);
|
||||
|
||||
//add curvature hard constraints
|
||||
//ScalarType Ratio=mesh.bbox.Diag()*0.01;
|
||||
|
||||
if (SParam.curv_thr>0)
|
||||
AddCurvatureConstraints(mesh,SParam.curv_thr);///Ratio);
|
||||
|
||||
//add alignment to sharp features
|
||||
if (SParam.sharp_thr>0)
|
||||
AddSharpEdgesConstraints(mesh,SParam.sharp_thr);
|
||||
|
||||
//add border constraints
|
||||
if (SParam.align_borders)
|
||||
AddBorderConstraints(mesh);
|
||||
}
|
||||
|
||||
static void GloballyOrient(MeshType &mesh)
|
||||
{
|
||||
vcg::tri::CrossField<MeshType>::MakeDirectionFaceCoherent(mesh,true);
|
||||
}
|
||||
|
||||
static void InitByCurvature(MeshType & mesh,int Nring)
|
||||
{
|
||||
|
||||
tri::RequirePerVertexCurvatureDir(mesh);
|
||||
|
||||
Eigen::MatrixXi F;
|
||||
Eigen::MatrixXd V;
|
||||
|
||||
Eigen::MatrixXd PD1,PD2,PV1,PV2;
|
||||
MeshToMatrix<MeshType>::GetTriMeshData(mesh,F,V);
|
||||
igl::principal_curvature(V,F,PD1,PD2,PV1,PV2,Nring);
|
||||
|
||||
//then copy curvature per vertex
|
||||
for (size_t i=0;i<mesh.vert.size();i++)
|
||||
{
|
||||
mesh.vert[i].PD1()=CoordType(PD1(i,0),PD1(i,1),PD1(i,2));
|
||||
mesh.vert[i].PD2()=CoordType(PD2(i,0),PD2(i,1),PD2(i,2));
|
||||
mesh.vert[i].K1()=PV1(i,0);
|
||||
mesh.vert[i].K2()=PV2(i,0);
|
||||
}
|
||||
vcg::tri::CrossField<MeshType>::SetFaceCrossVectorFromVert(mesh);
|
||||
InitQualityByAnisotropyDir(mesh);
|
||||
}
|
||||
|
||||
static void SmoothDirections(MeshType &mesh,
|
||||
int Ndir,
|
||||
SmoothMethod SMethod=SMNPoly,
|
||||
bool HardAsS=true,
|
||||
ScalarType alphaSoft=0)
|
||||
{
|
||||
|
||||
Eigen::VectorXi HardI; //hard constraints
|
||||
Eigen::MatrixXd HardD; //hard directions
|
||||
Eigen::VectorXi SoftI; //soft constraints
|
||||
Eigen::VectorXd SoftW; //weight of soft constraints
|
||||
Eigen::MatrixXd SoftD; //soft directions
|
||||
|
||||
if (HardAsS)
|
||||
CollectHardConstraints(mesh,HardI,HardD,SMethod,Ndir);
|
||||
|
||||
//collect soft constraints , miw only one that allows for soft constraints
|
||||
if ((alphaSoft>0)&&(SMethod==SMMiq))
|
||||
CollectSoftConstraints(mesh,SoftI,SoftD,SoftW);
|
||||
|
||||
//add some hard constraints if are not present
|
||||
int numC=3;
|
||||
if ((SoftI.rows()==0)&&(HardI.rows()==0))
|
||||
{
|
||||
printf("Add Forced Constraints \n");
|
||||
fflush(stdout);
|
||||
HardI=Eigen::MatrixXi(numC,1);
|
||||
|
||||
if ((Ndir==4)&&(SMethod==SMNPoly))
|
||||
HardD=Eigen::MatrixXd(numC,6);
|
||||
else
|
||||
HardD=Eigen::MatrixXd(numC,3);
|
||||
|
||||
for (int i=0;i<numC;i++)
|
||||
{
|
||||
CoordType Dir;
|
||||
int indexF;
|
||||
PickRandomDir(mesh,indexF,Dir);
|
||||
|
||||
HardI(i,0)=indexF;
|
||||
HardD(i,0)=Dir.X();
|
||||
HardD(i,1)=Dir.Y();
|
||||
HardD(i,2)=Dir.Z();
|
||||
|
||||
if ((Ndir==4)&&(SMethod==SMNPoly))
|
||||
{
|
||||
CoordType Dir1=mesh.face[indexF].N()^Dir;
|
||||
Dir1.Normalize();
|
||||
HardD(i,3)=Dir1.X();
|
||||
HardD(i,4)=Dir1.Y();
|
||||
HardD(i,5)=Dir1.Z();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//finally smooth
|
||||
if (SMethod==SMMiq)
|
||||
SmoothMIQ(mesh,HardI,HardD,SoftI,SoftD,SoftW,alphaSoft,Ndir);
|
||||
else
|
||||
{
|
||||
assert(SMethod==SMNPoly);
|
||||
SmoothNPoly(mesh,HardI,HardD,Ndir);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
static void SmoothDirections(MeshType &mesh,SmoothParam SParam)
|
||||
{
|
||||
//for the moment only cross and line field
|
||||
|
||||
//initialize direction by curvature if needed
|
||||
if ((SParam.alpha_curv>0)||
|
||||
(SParam.sharp_thr>0)||
|
||||
(SParam.curv_thr>0))
|
||||
InitByCurvature(mesh,SParam.curvRing);
|
||||
|
||||
SelectConstraints(mesh,SParam);
|
||||
//then do the actual smooth
|
||||
SmoothDirections(mesh,SParam.Ndir,SParam.SmoothM,true,SParam.alpha_curv);
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
} // end namespace tri
|
||||
} // end namespace vcg
|
||||
#endif // SMOOTHER_FIELD_H
|
||||
Loading…
Reference in New Issue