manolas
/
vcglib
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

smoothing v0.01 SMOOTH TETRA VERT IN VERTEXCOORDLAPLACIAN

This commit is contained in:
T.Alderighi 2018-05-08 13:16:45 +02:00
parent 25339718d9
commit 85501c3992
1 changed files with 1259 additions and 1182 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

205
vcg/complex/algorithms/smooth.h
View File

@ -21,7 +21,6 @@
* *
****************************************************************************/
#ifndef __VCGLIB__SMOOTH
#define __VCGLIB__SMOOTH
@ -31,7 +30,6 @@
#include <vcg/complex/algorithms/closest.h>
#include <vcg/space/index/kdtree/kdtree.h>
namespace vcg
{
namespace tri
@ -71,7 +69,6 @@ public:
// normal n with a speed equal to the mean curvature
void VertexCoordLaplacianCurvatureFlow(MeshType & /*m*/, int /*step*/, ScalarType /*delta*/)
{
}
// Another Laplacian smoothing variant,
@ -90,7 +87,8 @@ static void VertexCoordLaplacianAngleWeighted(MeshType &m, int step, ScalarType
for (vi = m.vert.begin(); vi != m.vert.end(); ++vi)
TD[*vi] = lpz;
ScalarType a[3];
for(fi=m.face.begin();fi!=m.face.end();++fi)if(!(*fi).IsD())
for (fi = m.face.begin(); fi != m.face.end(); ++fi)
if (!(*fi).IsD())
{
CoordType mp = ((*fi).V(0)->P() + (*fi).V(1)->P() + (*fi).V(2)->P()) / 3.0;
CoordType e0 = ((*fi).V(0)->P() - (*fi).V(1)->P()).Normalize();
@ -102,7 +100,8 @@ static void VertexCoordLaplacianAngleWeighted(MeshType &m, int step, ScalarType
a[2] = AngleN(-e2, e1);
//assert(fabs(M_PI -a[0] -a[1] -a[2])<0.0000001);
for(int j=0;j<3;++j){
for (int j = 0; j < 3; ++j)
{
CoordType dir = (mp - (*fi).V(j)->P()).Normalize();
TD[(*fi).V(j)].PntSum += dir * a[j];
TD[(*fi).V(j)].LenSum += a[j]; // well, it should be named angleSum
@ -111,7 +110,6 @@ static void VertexCoordLaplacianAngleWeighted(MeshType &m, int step, ScalarType
for (vi = m.vert.begin(); vi != m.vert.end(); ++vi)
if (!(*vi).IsD() && TD[*vi].LenSum > 0)
(*vi).P() = (*vi).P() + (TD[*vi].PntSum / TD[*vi].LenSum) * delta;
}
};
@ -138,9 +136,11 @@ static void VertexCoordScaleDependentLaplacian_Fujiwara(MeshType &m, int step, S
for (vi = m.vert.begin(); vi != m.vert.end(); ++vi)
TD[*vi] = lpz;
for(fi=m.face.begin();fi!=m.face.end();++fi)if(!(*fi).IsD())
for (fi = m.face.begin(); fi != m.face.end(); ++fi)
if (!(*fi).IsD())
for (int j = 0; j < 3; ++j)
if(!(*fi).IsB(j)) {
if (!(*fi).IsB(j))
{
CoordType edge = (*fi).V1(j)->P() - (*fi).V(j)->P();
ScalarType len = Norm(edge);
edge /= len;
@ -150,18 +150,20 @@ static void VertexCoordScaleDependentLaplacian_Fujiwara(MeshType &m, int step, S
TD[(*fi).V1(j)].LenSum += len;
}
for(fi=m.face.begin();fi!=m.face.end();++fi)if(!(*fi).IsD())
for (fi = m.face.begin(); fi != m.face.end(); ++fi)
if (!(*fi).IsD())
for (int j = 0; j < 3; ++j)
// se l'edge j e' di bordo si riazzera tutto e si riparte
if((*fi).IsB(j)) {
if ((*fi).IsB(j))
{
TD[(*fi).V(j)].PntSum = CoordType(0, 0, 0);
TD[(*fi).V1(j)].PntSum = CoordType(0, 0, 0);
TD[(*fi).V(j)].LenSum = 0;
TD[(*fi).V1(j)].LenSum = 0;
}
for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD())
for (fi = m.face.begin(); fi != m.face.end(); ++fi)
if (!(*fi).IsD())
for (int j = 0; j < 3; ++j)
if ((*fi).IsB(j))
{
@ -188,7 +190,6 @@ static void VertexCoordScaleDependentLaplacian_Fujiwara(MeshType &m, int step, S
}
};
class LaplacianInfo
{
public:
@ -212,6 +213,73 @@ public:
static void AccumulateLaplacianInfo(MeshType &m, SimpleTempData<typename MeshType::VertContainer, LaplacianInfo> &TD, bool cotangentFlag = false)
{
float weight = 1.0f;
//if we are applying to a tetrahedral mesh:
ForEachTetra(m, [&](MeshType::TetraType &t) {
for (int i = 0; i < 4; ++i)
if (!t.IsB(i))
{
VertexPointer v0, v1, v2;
v0 = t.V(Tetra::VofF(i, 0));
v1 = t.V(Tetra::VofF(i, 1));
v2 = t.V(Tetra::VofF(i, 2));
TD[v0].sum += v1->P() * weight;
TD[v0].sum += v2->P() * weight;
TD[v0].cnt += 2 * weight;
TD[v1].sum += v0->P() * weight;
TD[v1].sum += v2->P() * weight;
TD[v1].cnt += 2 * weight;
TD[v2].sum += v0->P() * weight;
TD[v2].sum += v1->P() * weight;
TD[v2].cnt += 2 * weight;
}
});
ForEachTetra(m, [&](MeshType::TetraType &t) {
for (int i = 0; i < 4; ++i)
if (t.IsB(i))
{
VertexPointer v0, v1, v2;
v0 = t.V(Tetra::VofF(i, 0));
v1 = t.V(Tetra::VofF(i, 1));
v2 = t.V(Tetra::VofF(i, 2));
TD[v0].sum = v0->P();
TD[v1].sum = v1->P();
TD[v2].sum = v2->P();
TD[v0].cnt = 1;
TD[v1].cnt = 1;
TD[v2].cnt = 1;
}
});
ForEachTetra(m, [&](MeshType::TetraType &t) {
for (int i = 0; i < 4; ++i)
if (t.IsB(i))
{
VertexPointer v0, v1, v2;
v0 = t.V(Tetra::VofF(i, 0));
v1 = t.V(Tetra::VofF(i, 1));
v2 = t.V(Tetra::VofF(i, 2));
TD[v0].sum += v1->P() * weight;
TD[v0].sum += v2->P() * weight;
TD[v0].cnt += 2 * weight;
TD[v1].sum += v0->P() * weight;
TD[v1].sum += v2->P() * weight;
TD[v1].cnt += 2 * weight;
TD[v2].sum += v0->P() * weight;
TD[v2].sum += v1->P() * weight;
TD[v2].cnt += 2 * weight;
}
});
FaceIterator fi;
for (fi = m.face.begin(); fi != m.face.end(); ++fi)
{
@ -219,7 +287,8 @@ static void AccumulateLaplacianInfo(MeshType &m, SimpleTempData<typename MeshTyp
for (int j = 0; j < 3; ++j)
if (!(*fi).IsB(j))
{
if(cotangentFlag) {
if (cotangentFlag)
{
float angle = Angle(fi->P1(j) - fi->P2(j), fi->P0(j) - fi->P2(j));
weight = tan((M_PI * 0.5) - angle);
}
@ -265,7 +334,8 @@ static void VertexCoordLaplacian(MeshType &m, int step, bool SmoothSelected=fals
SimpleTempData<typename MeshType::VertContainer, LaplacianInfo> TD(m.vert, lpz);
for (int i = 0; i < step; ++i)
{
if(cb)cb(100*i/step, "Classic Laplacian Smoothing");
if (cb)
cb(100 * i / step, "Classic Laplacian Smoothing");
TD.Init(lpz);
AccumulateLaplacianInfo(m, TD, cotangentWeight);
for (auto vi = m.vert.begin(); vi != m.vert.end(); ++vi)
@ -284,7 +354,8 @@ static void VertexCoordPlanarLaplacian(MeshType &m, int step, float AngleThrRad
SimpleTempData<typename MeshType::VertContainer, LaplacianInfo> TD(m.vert, lpz);
for (int i = 0; i < step; ++i)
{
if(cb)cb(100*i/step, "Planar Laplacian Smoothing");
if (cb)
cb(100 * i / step, "Planar Laplacian Smoothing");
TD.Init(lpz);
AccumulateLaplacianInfo(m, TD);
// First normalize the AccumulateLaplacianInfo
@ -295,18 +366,24 @@ static void VertexCoordPlanarLaplacian(MeshType &m, int step, float AngleThrRad
TD[*vi].sum = ((*vi).P() + TD[*vi].sum) / (TD[*vi].cnt + 1);
}
for(auto fi=m.face.begin();fi!=m.face.end();++fi){
if(!(*fi).IsD()){
for (int j = 0; j < 3; ++j) {
for (auto fi = m.face.begin(); fi != m.face.end(); ++fi)
{
if (!(*fi).IsD())
{
for (int j = 0; j < 3; ++j)
{
if (Angle(Normal(TD[(*fi).V0(j)].sum, (*fi).P1(j), (*fi).P2(j)),
Normal((*fi).P0(j), (*fi).P1(j), (*fi).P2(j))) > AngleThrRad)
TD[(*fi).V0(j)].sum = (*fi).P0(j);
}
}
}
for(auto fi=m.face.begin();fi!=m.face.end();++fi){
if(!(*fi).IsD()){
for (int j = 0; j < 3; ++j) {
for (auto fi = m.face.begin(); fi != m.face.end(); ++fi)
{
if (!(*fi).IsD())
{
for (int j = 0; j < 3; ++j)
{
if (Angle(Normal(TD[(*fi).V0(j)].sum, TD[(*fi).V1(j)].sum, (*fi).P2(j)),
Normal((*fi).P0(j), (*fi).P1(j), (*fi).P2(j))) > AngleThrRad)
{
@ -364,7 +441,6 @@ So if
* lambda == 0.5, kpb==0.01 -> mu = 1/(0.01 - 2) = -0.502
*/
static void VertexCoordTaubin(MeshType &m, int step, float lambda, float mu, bool SmoothSelected = false, vcg::CallBackPos *cb = 0)
{
LaplacianInfo lpz(CoordType(0, 0, 0), 0);
@ -372,7 +448,8 @@ static void VertexCoordTaubin(MeshType &m, int step, float lambda, float mu, boo
VertexIterator vi;
for (int i = 0; i < step; ++i)
{
if(cb) cb(100*i/step, "Taubin Smoothing");
if (cb)
cb(100 * i / step, "Taubin Smoothing");
TD.Init(lpz);
AccumulateLaplacianInfo(m, TD);
for (vi = m.vert.begin(); vi != m.vert.end(); ++vi)
@ -408,7 +485,8 @@ static void VertexQualityTaubin(MeshType &m, int step, float lambda, float mu, b
for (size_t i = 0; i < m.vert.size(); i++)
OldQ[i] = m.vert[i].Q();
if(cb) cb(100*i/step, "Taubin Smoothing");
if (cb)
cb(100 * i / step, "Taubin Smoothing");
VertexQualityLaplacian(m, 1, SmoothSelected);
@ -436,7 +514,6 @@ static void VertexQualityTaubin(MeshType &m, int step, float lambda, float mu, b
} // end for step
}
static void VertexCoordLaplacianQuality(MeshType &m, int step, bool SmoothSelected = false)
{
LaplacianInfo lpz;
@ -481,7 +558,8 @@ static void VertexCoordLaplacianHC(MeshType &m, int step, bool SmoothSelected=fa
SimpleTempData<typename MeshType::VertContainer, HCSmoothInfo> TD(m.vert, lpz);
// First Loop compute the laplacian
FaceIterator fi;
for(fi=m.face.begin();fi!=m.face.end();++fi)if(!(*fi).IsD())
for (fi = m.face.begin(); fi != m.face.end(); ++fi)
if (!(*fi).IsD())
{
for (int j = 0; j < 3; ++j)
{
@ -500,11 +578,13 @@ static void VertexCoordLaplacianHC(MeshType &m, int step, bool SmoothSelected=fa
}
}
VertexIterator vi;
for(vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
for (vi = m.vert.begin(); vi != m.vert.end(); ++vi)
if (!(*vi).IsD())
TD[*vi].sum /= (float)TD[*vi].cnt;
// Second Loop compute average difference
for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD())
for (fi = m.face.begin(); fi != m.face.end(); ++fi)
if (!(*fi).IsD())
{
for (int j = 0; j < 3; ++j)
{
@ -530,7 +610,6 @@ static void VertexCoordLaplacianHC(MeshType &m, int step, bool SmoothSelected=fa
// Laplacian smooth of the quality.
class ColorSmoothInfo
{
public:
@ -544,12 +623,16 @@ public:
static void VertexColorLaplacian(MeshType &m, int step, bool SmoothSelected = false, vcg::CallBackPos *cb = 0)
{
ColorSmoothInfo csi;
csi.r=0; csi.g=0; csi.b=0; csi.cnt=0;
csi.r = 0;
csi.g = 0;
csi.b = 0;
csi.cnt = 0;
SimpleTempData<typename MeshType::VertContainer, ColorSmoothInfo> TD(m.vert, csi);
for (int i = 0; i < step; ++i)
{
if(cb) cb(100*i/step, "Vertex Color Laplacian Smoothing");
if (cb)
cb(100 * i / step, "Vertex Color Laplacian Smoothing");
VertexIterator vi;
for (vi = m.vert.begin(); vi != m.vert.end(); ++vi)
TD[*vi] = csi;
@ -619,12 +702,16 @@ static void VertexColorLaplacian(MeshType &m, int step, bool SmoothSelected=fals
static void FaceColorLaplacian(MeshType &m, int step, bool SmoothSelected = false, vcg::CallBackPos *cb = 0)
{
ColorSmoothInfo csi;
csi.r=0; csi.g=0; csi.b=0; csi.cnt=0;
csi.r = 0;
csi.g = 0;
csi.b = 0;
csi.cnt = 0;
SimpleTempData<typename MeshType::FaceContainer, ColorSmoothInfo> TD(m.face, csi);
for (int i = 0; i < step; ++i)
{
if(cb) cb(100*i/step, "Face Color Laplacian Smoothing");
if (cb)
cb(100 * i / step, "Face Color Laplacian Smoothing");
FaceIterator fi;
for (fi = m.face.begin(); fi != m.face.end(); ++fi)
TD[*fi] = csi;
@ -705,7 +792,6 @@ static void PointCloudQualityMedian(MeshType &m, int medianSize=8)
for (int i = 0; i < m.vn; ++i)
m.vert[i].Q() = newQVec[i];
}
static void VertexQualityLaplacian(MeshType &m, int step = 1, bool SmoothSelected = false)
@ -817,7 +903,6 @@ static void VertexNormalLaplacian(MeshType &m, int step,bool SmoothSelected=fals
if (!SmoothSelected || (*vi).IsS())
(*vi).N() = TD[*vi].sum / TD[*vi].cnt;
}
};
// Smooth solo lungo la direzione di vista
@ -877,15 +962,14 @@ static void VertexCoordViewDepth(MeshType &m,
if (!(*vi).IsD() && TD[*vi].cnt > 0)
{
CoordType np = TD[*vi].sum / TD[*vi].cnt;
CoordType d = (*vi).cP() - viewpoint; d.Normalize();
CoordType d = (*vi).cP() - viewpoint;
d.Normalize();
ScalarType s = d.dot(np - (*vi).cP());
(*vi).P() += d * (s * alpha);
}
}
}
/****************************************************************************************************************/
/****************************************************************************************************************/
// Paso Double Smoothing
@ -905,7 +989,6 @@ public:
CoordType np;
};
class PDFaceInfo
{
public:
@ -923,7 +1006,6 @@ public:
// This is the Normal Smoothing approach of Shen and Berner
// Fuzzy Vector Median-Based Surface Smoothing TVCG 2004
void FaceNormalFuzzyVectorSB(MeshType &m,
SimpleTempData<typename MeshType::FaceContainer, PDFaceInfo> &TD,
ScalarType sigma)
@ -962,7 +1044,8 @@ void FaceNormalFuzzyVectorSB(MeshType &m,
ScalarType ang = AngleN(ep.f->N(), (*fi).N());
mm += ep.f->N() * exp((-sigma) * ang * ang / dd);
}
else mm+=ep.f->N();
else
mm += ep.f->N();
(*ep.f).SetV();
}
++ep;
@ -986,7 +1069,8 @@ static void FaceNormalLaplacianVF(MeshType &m)
tri::UpdateNormal<MeshType>::NormalizePerFaceByArea(m);
for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD())
for (FaceIterator fi = m.face.begin(); fi != m.face.end(); ++fi)
if (!(*fi).IsD())
{
// 1) Clear all the visited flag of faces that are vertex-adjacent to fi
for (int i = 0; i < 3; ++i)
@ -1026,7 +1110,6 @@ static void FaceNormalLaplacianVF(MeshType &m)
// Normals are normalized:
// FF adjacency is present.
static void FaceNormalLaplacianFF(MeshType &m, int step = 1, bool SmoothSelected = false)
{
PDFaceInfo lpzf(CoordType(0, 0, 0));
@ -1037,7 +1120,8 @@ static void FaceNormalLaplacianFF(MeshType &m, int step=1, bool SmoothSelected=f
tri::UpdateNormal<MeshType>::NormalizePerFaceByArea(m);
for (int iStep = 0; iStep < step; ++iStep)
{
for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD())
for (fi = m.face.begin(); fi != m.face.end(); ++fi)
if (!(*fi).IsD())
{
CoordType normalSum = (*fi).N();
@ -1054,7 +1138,6 @@ static void FaceNormalLaplacianFF(MeshType &m, int step=1, bool SmoothSelected=f
}
}
/***************************************************************************/
// Paso Doble Step 1 compute the smoothed normals
/***************************************************************************/
@ -1068,7 +1151,6 @@ static void FaceNormalLaplacianFF(MeshType &m, int step=1, bool SmoothSelected=f
// sigma == 1 Nothing is averaged.
// Only within the specified range are averaged toghether. The averagin is weighted with the
static void FaceNormalAngleThreshold(MeshType &m,
SimpleTempData<typename MeshType::FaceContainer, PDFaceInfo> &TD,
ScalarType sigma)
@ -1169,14 +1251,12 @@ A(...) (2-2nm) =
static CoordType FaceErrorGrad(CoordType &p, CoordType &p0, CoordType &p1, CoordType &m)
{
return TriAreaGradient(p,p0,p1) *2.0f
- CrossProdGradient(p,p0,p1,m) *2.0f ;
return TriAreaGradient(p, p0, p1) * 2.0f - CrossProdGradient(p, p0, p1, m) * 2.0f;
}
/***************************************************************************/
// Paso Doble Step 2 Fitta la mesh a un dato insieme di normali
/***************************************************************************/
static void FitMesh(MeshType &m,
SimpleTempData<typename MeshType::VertContainer, PDVertInfo> &TDV,
SimpleTempData<typename MeshType::FaceContainer, PDFaceInfo> &TDF,
@ -1200,12 +1280,9 @@ static void FitMesh(MeshType &m,
for (vi = m.vert.begin(); vi != m.vert.end(); ++vi)
(*vi).P() = TDV[*vi].np;
}
/****************************************************************************************************************/
static void FastFitMesh(MeshType &m,
SimpleTempData<typename MeshType::VertContainer, PDVertInfo> &TDV,
//SimpleTempData<typename MeshType::FaceContainer, PDFaceInfo> &TDF,
@ -1232,7 +1309,8 @@ static void FastFitMesh(MeshType &m,
if (OnlySelected)
{
for (vi = m.vert.begin(); vi != m.vert.end(); ++vi)
if((*vi).IsS()) (*vi).P()=TDV[*vi].np;
if ((*vi).IsS())
(*vi).P() = TDV[*vi].np;
}
else
{
@ -1241,8 +1319,6 @@ static void FastFitMesh(MeshType &m,
}
}
// The sigma parameter affect the normal smoothing step
static void VertexCoordPasoDoble(MeshType &m, int NormalSmoothStep, typename MeshType::ScalarType Sigma = 0, int FitStep = 50, bool SmoothSelected = false)
@ -1264,14 +1340,15 @@ static void VertexCoordPasoDoble(MeshType &m, int NormalSmoothStep, typename Mes
FastFitMesh(m, TDV, SmoothSelected);
}
static void VertexNormalPointCloud(MeshType &m, int neighborNum, int iterNum, KdTree<ScalarType> *tp = 0)
{
SimpleTempData<typename MeshType::VertContainer, CoordType> TD(m.vert, CoordType(0, 0, 0));
VertexConstDataWrapper<MeshType> ww(m);
KdTree<ScalarType> *tree = 0;
if(tp==0) tree = new KdTree<ScalarType>(ww);
else tree=tp;
if (tp == 0)
tree = new KdTree<ScalarType>(ww);
else
tree = tp;
typename KdTree<ScalarType>::PriorityQueue nq;
// tree->setMaxNofNeighbors(neighborNum);
@ -1298,7 +1375,8 @@ static void VertexNormalPointCloud(MeshType &m, int neighborNum, int iterNum, Kd
tri::UpdateNormal<MeshType>::NormalizePerVertex(m);
}
if(tp==0) delete tree;
if (tp == 0)
delete tree;
}
//! Laplacian smoothing with a reprojection on a target surface.
@ -1313,7 +1391,6 @@ static void VertexCoordLaplacianReproject(MeshType& m, GRID& grid, MeshTypeTri&
}
}
template <class GRID, class MeshTypeTri>
static void VertexCoordLaplacianReproject(MeshType &m, GRID &grid, MeshTypeTri &gridmesh, typename MeshType::VertexType *vp)
{
@ -1334,7 +1411,6 @@ static void VertexCoordLaplacianReproject(MeshType& m, GRID& grid, MeshTypeTri&
// move vertex
vp->P() = ct;
vector<FacePointer> faces2 = HalfEdgeTopology<MeshType>::get_incident_faces(vp);
// estimate normal
@ -1349,7 +1425,6 @@ static void VertexCoordLaplacianReproject(MeshType& m, GRID& grid, MeshTypeTri&
avgn += vcg::Normal<typename MeshType::CoordType>(vertices[0]->cP(), vertices[1]->cP(), vertices[2]->cP());
avgn += vcg::Normal<typename MeshType::CoordType>(vertices[2]->cP(), vertices[3]->cP(), vertices[0]->cP());
}
avgn.Normalize();
@ -1384,7 +1459,9 @@ static void VertexCoordLaplacianReproject(MeshType& m, GRID& grid, MeshTypeTri&
faces.push_back(fr);
}
if (fr) if (fr->N()*raydir<0) fr=0; // discard: inverse normal;
if (fr)
if (fr->N() * raydir < 0)
fr = 0; // discard: inverse normal;
typename MeshType::CoordType newPos;
if (minDists.size() == 0)
{