diff --git a/vcg/complex/trimesh/refine_loop.h b/vcg/complex/trimesh/refine_loop.h index 81fb1ef3..e689edbd 100644 --- a/vcg/complex/trimesh/refine_loop.h +++ b/vcg/complex/trimesh/refine_loop.h @@ -20,11 +20,18 @@ * for more details. * * * ****************************************************************************/ +/*! \file refine_loop.h + * + * \brief This file contain code for Loop's subdivision scheme for triangular + * mesh and some similar method. + * + */ #ifndef __VCGLIB_REFINE_LOOP #define __VCGLIB_REFINE_LOOP -#include +#include +#include #include #include #include @@ -53,37 +60,344 @@ d5------d1------d2 -> d5--e5--d1--e2--d2 l--M--r */ -// Nuovi punti (e.g. midpoint), ossia odd vertices -// -template -struct OddPointLoop : public std::unary_function , typename MESH_TYPE::CoordType> -{ - void operator()(typename MESH_TYPE::VertexType &nv, face::Pos ep) { +/*! + * \brief Weight class for classical Loop's scheme. + * + * See Zorin, D. & Schröeder, P.: Subdivision for modeling and animation. Proc. ACM SIGGRAPH [Courses], 2000 + */ +template +struct LoopWeight { + inline SCALAR_TYPE beta(int k) { + return (k>3)?(5.0/8.0 - std::pow((3.0/8.0 + std::cos(2.0*M_PI/SCALAR_TYPE(k))/4.0),2))/SCALAR_TYPE(k):3.0/16.0; + } + + inline SCALAR_TYPE incidentRegular(int) { + return 3.0/8.0; + } + inline SCALAR_TYPE incidentIrregular(int) { + return 3.0/8.0; + } + inline SCALAR_TYPE opposite(int) { + return 1.0/8.0; + } +}; +/*! + * \brief Modified Loop's weight to optimise continuity. + * + * See Barthe, L. & Kobbelt, L.: Subdivision scheme tuning around extraordinary vertices. Computer Aided Geometric Design, 2004, 21, 561-583 + */ +template +struct RegularLoopWeight { + inline SCALAR_TYPE beta(int k) { + static SCALAR_TYPE bkPolar[] = { + .32517, + .49954, + .59549, + .625, + .63873, + .64643, + .65127, + .67358, + .68678, + .69908 + }; + + return (k<=12)?(1.0-bkPolar[k-3])/k:LoopWeight().beta(k); + } + + inline SCALAR_TYPE incidentRegular(int k) { + return 1.0 - incidentIrregular(k) - opposite(k)*2; + } + inline SCALAR_TYPE incidentIrregular(int k) { + static SCALAR_TYPE bkPolar[] = { + .15658, + .25029, + .34547, + .375, + .38877, + .39644, + .40132, + .42198, + .43423, + .44579 + }; + + return (k<=12)?bkPolar[k-3]:LoopWeight().incidentIrregular(k); + } + inline SCALAR_TYPE opposite(int k) { + static SCALAR_TYPE bkPolar[] = { + .14427, + .12524, + .11182, + .125, + .14771, + .1768, + .21092, + .20354, + .20505, + .19828 + }; + + return (k<=12)?bkPolar[k-3]:LoopWeight().opposite(k); + } +}; + +template +struct ContinuityLoopWeight { + inline SCALAR_TYPE beta(int k) { + static SCALAR_TYPE bkPolar[] = { + .32517, + .50033, + .59464, + .625, + .63903, + .67821, + .6866, + .69248, + .69678, + .70014 + }; + + return (k<=12)?(1.0-bkPolar[k-3])/k:LoopWeight().beta(k); + } + + inline SCALAR_TYPE incidentRegular(int k) { + return 1.0 - incidentIrregular(k) - opposite(k)*2; + } + inline SCALAR_TYPE incidentIrregular(int k) { + static SCALAR_TYPE bkPolar[] = { + .15658, + .26721, + .33539, + .375, + .36909, + .25579, + .2521, + .24926, + .24706, + .2452 + }; + + return (k<=12)?bkPolar[k-3]:LoopWeight().incidentIrregular(k); + } + inline SCALAR_TYPE opposite(int k) { + static SCALAR_TYPE bkPolar[] = { + .14427, + .12495, + .11252, + .125, + .14673, + .16074, + .18939, + .2222, + .25894, + .29934 + }; + + return (k<=12)?bkPolar[k-3]:LoopWeight().opposite(k); + } +}; + +// Centroid and LS3Projection classes may be pettre placed in an other file. (which one ?) + +/*! + * \brief Allow to compute classical Loop subdivision surface with the same code than LS3. + */ +template +struct Centroid { + typedef typename MESH_TYPE::ScalarType Scalar; + typedef typename MESH_TYPE::CoordType Coord; + typedef LSCALAR_TYPE LScalar; + typedef vcg::Point3 LVector; + + LVector sumP; + LScalar sumW; + + Centroid() { reset(); } + inline void reset() { + sumP.SetZero(); + sumW = 0.; + } + inline void addVertex(const typename MESH_TYPE::VertexType &v, LScalar w) { + LVector p(v.cP().X(), v.cP().Y(), v.cP().Z()); + LVector n(v.cN().X(), v.cN().Y(), v.cN().Z()); + + sumP += p * w; + sumW += w; + } + inline void project(typename MESH_TYPE::VertexType &v) const { + LVector position = sumP / sumW; + v.P() = Coord(position.X(), position.Y(), position.Z()); + } +}; + +/*! Implementation of the APSS projection for the LS3 scheme. + * + * See Gael Guennebaud and Marcel Germann and Markus Gross + * Dynamic sampling and rendering of algebraic point set surfaces. + * Computer Graphics Forum (Proceedings of Eurographics 2008), 2008, 27, 653-662 + * and Simon Boye and Gael Guennebaud and Christophe Schlick + * Least squares subdivision surfaces + * Computer Graphics Forum, 2010 + */ +template +struct LS3Projection { + typedef typename MESH_TYPE::ScalarType Scalar; + typedef typename MESH_TYPE::CoordType Coord; + typedef LSCALAR_TYPE LScalar; + typedef vcg::Point3 LVector; + + Scalar beta; + + LVector sumP; + LVector sumN; + LScalar sumDotPN; + LScalar sumDotPP; + LScalar sumW; + + inline LS3Projection(Scalar beta = 1.0) : beta(beta) { reset(); } + inline void reset() { + sumP.SetZero(); + sumN.SetZero(); + sumDotPN = 0.; + sumDotPP = 0.; + sumW = 0.; + } + inline void addVertex(const typename MESH_TYPE::VertexType &v, LScalar w) { + LVector p(v.cP().X(), v.cP().Y(), v.cP().Z()); + LVector n(v.cN().X(), v.cN().Y(), v.cN().Z()); + + sumP += p * w; + sumN += n * w; + sumDotPN += w * n.dot(p); + sumDotPP += w * vcg::SquaredNorm(p); + sumW += w; + } + void project(typename MESH_TYPE::VertexType &v) const { + LScalar invSumW = Scalar(1)/sumW; + LScalar aux4 = beta * LScalar(0.5) * + (sumDotPN - invSumW*sumP.dot(sumN)) + /(sumDotPP - invSumW*vcg::SquaredNorm(sumP)); + LVector uLinear = (sumN-sumP*(Scalar(2)*aux4))*invSumW; + LScalar uConstant = -invSumW*(uLinear.dot(sumP) + sumDotPP*aux4); + LScalar uQuad = aux4; + LVector orig = sumP*invSumW; + + // finalize + LVector position; + LVector normal; + if (fabs(uQuad)>1e-7) + { + LScalar b = 1./uQuad; + LVector center = uLinear*(-0.5*b); + LScalar radius = sqrt( vcg::SquaredNorm(center) - b*uConstant ); + + normal = orig - center; + normal.Normalize(); + position = center + normal * radius; + + normal = uLinear + position * (LScalar(2) * uQuad); + normal.Normalize(); + } + else if (uQuad==0.) + { + LScalar s = LScalar(1)/vcg::Norm(uLinear); + assert(!vcg::math::IsNAN(s) && "normal should not have zero len!"); + uLinear *= s; + uConstant *= s; + + normal = uLinear; + position = orig - uLinear * (orig.dot(uLinear) + uConstant); + } + else + { + // normalize the gradient + LScalar f = 1./sqrt(vcg::SquaredNorm(uLinear) - Scalar(4)*uConstant*uQuad); + uConstant *= f; + uLinear *= f; + uQuad *= f; + + // Newton iterations + LVector grad; + LVector dir = uLinear+orig*(2.*uQuad); + LScalar ilg = 1./vcg::Norm(dir); + dir *= ilg; + LScalar ad = uConstant + uLinear.dot(orig) + uQuad * vcg::SquaredNorm(orig); + LScalar delta = -ad*std::min(ilg,1.); + LVector p = orig + dir*delta; + for (int i=0 ; i<2 ; ++i) + { + grad = uLinear+p*(2.*uQuad); + ilg = 1./vcg::Norm(grad); + delta = -(uConstant + uLinear.dot(p) + uQuad * vcg::SquaredNorm(p))*std::min(ilg,1.); + p += dir*delta; + } + position = p; + + normal = uLinear + position * (Scalar(2) * uQuad); + normal.Normalize(); + } + + v.P() = Coord(position.X(), position.Y(), position.Z()); + v.N() = Coord(normal.X(), normal.Y(), normal.Z()); + } +}; + +template, class WEIGHT_TYPE=LoopWeight > +struct OddPointLoopGeneric : public std::unary_function , typename MESH_TYPE::VertexType> +{ + typedef METHOD_TYPE Projection; + typedef WEIGHT_TYPE Weight; + typedef typename MESH_TYPE::template PerVertexAttributeHandle ValenceAttr; + + Projection proj; + Weight weight; + ValenceAttr *valence; + + inline OddPointLoopGeneric(Projection proj = Projection(), Weight weight = Weight()) : + proj(proj), weight(weight), valence(0) {} + + void operator()(typename MESH_TYPE::VertexType &nv, face::Pos ep) { + proj.reset(); + face::Pos he(ep.f,ep.z,ep.f->V(ep.z)); - typename MESH_TYPE::CoordType *l,*r,*u,*d; - l = &he.v->P(); + typename MESH_TYPE::VertexType *l,*r,*u,*d; + l = he.v; he.FlipV(); - r = &he.v->P(); + r = he.v; if( MESH_TYPE::HasPerVertexColor()) nv.C().lerp(ep.f->V(ep.z)->C(),ep.f->V1(ep.z)->C(),.5f); if (he.IsBorder()) { - nv.P() = ((*l)*0.5 + (*r)*0.5); - + proj.addVertex(*l, 0.5); + proj.addVertex(*r, 0.5); + proj.project(nv); } else { he.FlipE(); he.FlipV(); - u = &he.v->P(); + u = he.v; he.FlipV(); he.FlipE(); - assert(&he.v->P()== r); // back to r + assert(he.v == r); // back to r he.FlipF(); he.FlipE(); he.FlipV(); - d = &he.v->P(); + d = he.v; - // abbiamo i punti l,r,u e d per ottenere M in maniera pesata - - nv.P()=((*l)*(3.0/8.0)+(*r)*(3.0/8.0)+(*d)*(1.0/8.0)+(*u)*(1.0/8.0)); + if(valence && ((*valence)[l]==6 || (*valence)[r]==6)) { + int extra = ((*valence)[l]==6)?(*valence)[r]:(*valence)[l]; + proj.addVertex(*l, ((*valence)[l]==6)?weight.incidentRegular(extra):weight.incidentIrregular(extra)); + proj.addVertex(*r, ((*valence)[r]==6)?weight.incidentRegular(extra):weight.incidentIrregular(extra)); + proj.addVertex(*u, weight.opposite(extra)); + proj.addVertex(*d, weight.opposite(extra)); + } + // unhandled case that append only at first subdivision step: use Loop's weights + else { + proj.addVertex(*l, 3.0/8.0); + proj.addVertex(*r, 3.0/8.0); + proj.addVertex(*u, 1.0/8.0); + proj.addVertex(*d, 1.0/8.0); + } + proj.project(nv); } } @@ -102,58 +416,78 @@ struct OddPointLoop : public std::unary_functionvalence = valence; + } }; -// vecchi punti, ossia even vertices -template -struct EvenPointLoop : public std::unary_function , typename MESH_TYPE::CoordType> +template, class WEIGHT_TYPE=LoopWeight > +struct EvenPointLoopGeneric : public std::unary_function , typename MESH_TYPE::VertexType> { + typedef METHOD_TYPE Projection; + typedef WEIGHT_TYPE Weight; + typedef typename MESH_TYPE::template PerVertexAttributeHandle ValenceAttr; - void operator()(typename MESH_TYPE::CoordType &nP, face::Pos ep) { - + Projection proj; + Weight weight; + ValenceAttr *valence; + + inline EvenPointLoopGeneric(Projection proj = Projection(), Weight weight = Weight()) : + proj(proj), weight(weight), valence(0) {} + + void operator()(typename MESH_TYPE::VertexType &nv, face::Pos ep) { + proj.reset(); + face::Pos he(ep.f,ep.z,ep.f->V(ep.z)); - typename MESH_TYPE::CoordType *r, *l, *curr; - curr = &he.v->P(); + typename MESH_TYPE::VertexType *r, *l, *curr; + curr = he.v; + face::Pos heStart = he; - if (he.IsBorder()) {//half edge di bordo - he.FlipV(); - r = &he.v->P(); - he.FlipV(); - assert(&he.v->P()== curr); // back to curr - he.NextB(); - if (&he.v->P() == curr) - he.FlipV(); - l = &he.v->P(); - nP = ( *(curr) * (3.0)/(4.0) + (*l)*(1.0/8.0) + (*r)*(1.0/8.0)); - } - else { - // compute valence of this vertex - int k = 0; - face::Pos heStart = he; - std::vector otherVertVec; - if(he.v->IsB())return ; - do { - he.FlipV(); - otherVertVec.push_back(he.v->P()); - he.FlipV(); - he.FlipE(); he.FlipF(); - k++; - } while(he.f!=heStart.f || he.z!=heStart.z || he.v!=heStart.v); - // while(he != heStart); - - float beta = 3.0 / 16.0; - if(k > 3 ) - beta = (1.0/(float)k) * (5.0/8.0 - pow((3.0/8.0 + 0.25 * cos(2*M_PI/k)),2)); - - *curr = *curr * (1 - k * beta) ; - typename std::vector::iterator iter; - for (iter = otherVertVec.begin(); - iter != otherVertVec.end(); - ++iter) { - *curr = *curr + (*iter) * beta; + // compute valence of this vertex or find a border + int k = 0; + do { + he.NextE(); + k++; + } while(!he.IsBorder() && he != heStart); + if (he.IsBorder()) { // Border rule + // consider valence of borders as if they are half+1 of an inner vertex. not perfect, but better than nothing. + if(valence) { + k = 0; + do { + he.NextE(); + k++; + } while(!he.IsBorder()); + (*valence)[he.V()] = std::max(2*(k-1), 3); +// (*valence)[he.V()] = 6; } - nP = *curr; + + he.FlipV(); + r = he.v; + he.FlipV(); + he.NextB(); + l = he.v; + + proj.addVertex(*curr, 3.0/4.0); + proj.addVertex(*l, 1.0/8.0); + proj.addVertex(*r, 1.0/8.0); + proj.project(nv); + } + else { // Inner rule +// assert(!he.v->IsB()); border flag no longer updated (useless) + if(valence) + (*valence)[he.V()] = k; + + typename MESH_TYPE::ScalarType beta = weight.beta(k); + + proj.addVertex(*curr, 1.0 - (typename MESH_TYPE::ScalarType)(k) * beta); + do { + proj.addVertex(*he.VFlip(), beta); + he.NextE(); + } while(he != heStart); + + proj.project(nv); } } // end of operator() @@ -179,14 +513,20 @@ struct EvenPointLoop : public std::unary_functionvalence = valence; + } }; -template struct EvenParam { - CoordType sum; - bool border; - int k; -} ; +template +struct OddPointLoop : OddPointLoopGeneric > +{ +}; +template +struct EvenPointLoop : EvenPointLoopGeneric > +{ +}; template bool RefineOddEven(MESH_TYPE &m, ODD_VERT odd, EVEN_VERT even,float length, @@ -196,41 +536,44 @@ bool RefineOddEven(MESH_TYPE &m, ODD_VERT odd, EVEN_VERT even,float length, return RefineOddEvenE(m, odd, even, ep, RefineSelected, cbOdd, cbEven); } +/*! + * \brief Perform diadic subdivision using given rules for odd and even vertices. + */ template bool RefineOddEvenE(MESH_TYPE &m, ODD_VERT odd, EVEN_VERT even, PREDICATE edgePred, bool RefineSelected=false, CallBackPos *cbOdd = 0, CallBackPos *cbEven = 0) { + typedef typename MESH_TYPE::template PerVertexAttributeHandle ValenceAttr; - // n = numero di vertici iniziali - int n = m.vn; - - // refine dei vertici odd, crea dei nuovi vertici in coda - RefineE< MESH_TYPE,OddPointLoop > (m, odd, edgePred, RefineSelected, cbOdd); // momentaneamente le callback sono identiche, almeno cbOdd deve essere passata cbEven = cbOdd; - vcg::tri::UpdateFlags::FaceBorderFromFF(m); - // aggiorno i flag perche' IsB funzioni - vcg::tri::UpdateFlags::VertexBorderFromFace (m); - //vcg::tri::UpdateColor::VertexBorderFlag(m); - - // marco i vertici even [ i primi n ] come visitati + // to mark visited vertices int evenFlag = MESH_TYPE::VertexType::NewBitFlag(); - for (int i = 0; i < n ; i++ ) { - m.vert[i].SetUserBit(evenFlag); + for (int i = 0; i < m.vn ; i++ ) { + m.vert[i].ClearUserBit(evenFlag); } - int j = 0; // di texture per wedge (uno per ogni edge) + ValenceAttr valence = vcg::tri::Allocator::AddPerVertexAttribute(m); + odd.setValenceAttr(&valence); + even.setValenceAttr(&valence); + + // store updated vertices + std::vector updatedList(m.vn, false); + std::vector newEven(m.vn); + typename MESH_TYPE::VertexIterator vi; typename MESH_TYPE::FaceIterator fi; - for (fi = m.face.begin(); fi != m.face.end(); fi++) if(!(*fi).IsD()){ //itero facce + for (fi = m.face.begin(); fi != m.face.end(); fi++) if(!(*fi).IsD() && (!RefineSelected || (*fi).IsS())){ //itero facce for (int i = 0; i < 3; i++) { //itero vert - if ( (*fi).V(i)->IsUserBit(evenFlag) && ! (*fi).V(i)->IsD() ) { - if (RefineSelected && !(*fi).V(i)->IsS() ) - break; + if ( !(*fi).V(i)->IsUserBit(evenFlag) && ! (*fi).V(i)->IsD() ) { + (*fi).V(i)->SetUserBit(evenFlag); + // use face selection, not vertex selection, to be coherent with RefineE + //if (RefineSelected && !(*fi).V(i)->IsS() ) + // break; face::Posaux (&(*fi),i); if( MESH_TYPE::HasPerVertexColor() ) { (*fi).V(i)->C().lerp((*fi).V0(i)->C() , (*fi).V1(i)->C(),0.5f); @@ -240,16 +583,34 @@ bool RefineOddEvenE(MESH_TYPE &m, ODD_VERT odd, EVEN_VERT even, PREDICATE edgePr (*cbEven)(int(100.0f * (float)j / (float)m.fn),"Refining"); j++; } - even((*fi).V(i)->P(), aux); + int index = tri::Index(m, (*fi).V(i)); + updatedList[index] = true; + even(newEven[index], aux); } } } + + MESH_TYPE::VertexType::DeleteBitFlag(evenFlag); + // refine dei vertici odd, crea dei nuovi vertici in coda + RefineE< MESH_TYPE, ODD_VERT > (m, odd, edgePred, RefineSelected, cbOdd); + + typename std::vector::iterator nei; + for (nei=newEven.begin(); nei!=newEven.end(); ++nei) { + if(updatedList[nei-newEven.begin()]) { + m.vert[nei-newEven.begin()].ImportData(*nei); + assert(m.vert[nei-newEven.begin()].N() == nei->N()); + } + } + + odd.setValenceAttr(0); + even.setValenceAttr(0); + + vcg::tri::Allocator::DeletePerVertexAttribute(m, valence); + return true; } - - } // namespace tri } // namespace vcg