From b952d96711bc45170eae72ed3510c740b185655c Mon Sep 17 00:00:00 2001 From: nicopietroni Date: Sun, 24 May 2015 14:52:50 +0000 Subject: [PATCH] Added support for tracing separatrix from singularities --- .../parametrization/tangent_field_operators.h | 516 +++++++++++++++++- 1 file changed, 513 insertions(+), 3 deletions(-) diff --git a/vcg/complex/algorithms/parametrization/tangent_field_operators.h b/vcg/complex/algorithms/parametrization/tangent_field_operators.h index 87c36356..2365da8d 100644 --- a/vcg/complex/algorithms/parametrization/tangent_field_operators.h +++ b/vcg/complex/algorithms/parametrization/tangent_field_operators.h @@ -32,6 +32,8 @@ namespace vcg { namespace tri{ + + template < typename ScalarType > vcg::Point2 InterpolateNRosy2D(const std::vector > &V, const std::vector &W, @@ -108,17 +110,25 @@ vcg::Point3 InterpolateNRosy3D(const std::vector RotNorm=vcg::RotationMatrix(Norm[i],TargetN); + //std::cout << "Norm[i] " << Norm[i].X() << " " << Norm[i].Y() << " " << Norm[i].Z()<< std::endl; + //std::cout << "TargetN " << TargetN.X() << " " << TargetN.Y() << " " << TargetN.Z()<< std::endl<< std::flush; + CoordType rotV=RotNorm*V[i]; //assert(fabs(rotV*TargetN)<0.000001); rotV.Normalize(); + //std::cout << "rotV " << rotV.X() << " " << rotV.Y() << " " << rotV.Z()<< std::endl<< std::flush; + ///trassform to the reference frame rotV=RotFrame*rotV; //it's 2D from now on Cross2D.push_back(vcg::Point2(rotV.X(),rotV.Y())); + } + vcg::Point2 AvDir2D=InterpolateNRosy2D(Cross2D,W,N); CoordType AvDir3D=CoordType(AvDir2D.X(),AvDir2D.Y(),0); //transform back to 3D @@ -133,12 +143,512 @@ class CrossField typedef typename MeshType::VertexType VertexType; typedef typename MeshType::CoordType CoordType; typedef typename MeshType::ScalarType ScalarType; + typedef typename vcg::face::Pos PosType; + typedef typename vcg::Triangle3 TriangleType; private: + + static void SubDivideDir(const CoordType &Edge0, + const CoordType &Edge1, + std::vector &SubDEdges, + int Nsub) + { + CoordType Dir0=Edge0; + CoordType Dir1=Edge1; + ScalarType a=Edge0.Norm(); + Dir0.Normalize(); + Dir1.Normalize(); + /* + * + * + * A + * / | + * Dir1 / | + * ^ / | + * | c / | + * / / |b + * | / | + * / / | + * / | + * / | + * B____________________C ->Dir0 + * a + * + * */ + ScalarType BTotal=vcg::Angle(Dir0,Dir1); + //get angle step + ScalarType StepAngle=BTotal/(ScalarType)Nsub; + //get the other edge + CoordType Edge2=Edge1-Edge0; + + //and its direction + CoordType Dir2=Edge2; + Dir2.Normalize(); + //find angle C + ScalarType C=vcg::Angle(Dir2,-Dir0); + //safety checks + assert(BTotal<=(M_PI)); + assert(BTotal>=0); + assert(C<=(M_PI)); + assert(C>=0); + + //push the first one + SubDEdges.push_back(Edge0); + for (size_t i=1;i=0); + //find current lenght + ScalarType b=a*sin(B)/sin(A); + //then move along the direction of edge b + CoordType intervB=Dir2*b; + //finally find absolute position summing edge 0 + intervB+=Edge0; + SubDEdges.push_back(intervB); + } + //push the last one + SubDEdges.push_back(Edge1); + } + + static void FindSubDir(vcg::Triangle3 T3, + size_t Nvert, + std::vector &SubDEdges, + int Nsub) + { + CoordType P0=T3.P0(Nvert); + CoordType P1=T3.P1(Nvert); + CoordType P2=T3.P2(Nvert); + P1-=P0; + P2-=P0; + SubDivideDir(P1,P2,SubDEdges,Nsub); + for (size_t j=0;j T3, + size_t Nvert, + std::vector > &SubTris, + int Nsub) + { + std::vector SplittedPos; + FindSubDir(T3,Nvert,SplittedPos,Nsub); + CoordType P0=T3.P(Nvert); + //then create the triangles + for (size_t j=0;j0)?+1:-1);} + static void FindSubTriangles(const typename vcg::face::Pos &vPos, + std::vector &SubFaces, + std::vector &OriginalFace, + int numSub=3) + { + SubFaces.clear(); + OriginalFace.clear(); + + //not ct on border + assert(!vPos.IsBorder()); + //the vertex should be the one on the pos + assert(vPos.F()->V(vPos.E())==vPos.V()); + + // collect all faces around the star of the vertex + std::vector StarPos; + vcg::face::VFOrderedStarFF(vPos, StarPos); + + //all the infos for the strip of triangles + + VertexType* v0=vPos.V(); + CoordType P0=v0->P(); + //create the strip of triangles + for (size_t i = 0; i < StarPos.size(); ++i) + { + PosType currPos=StarPos[i]; + FaceType *CurrF=currPos.F(); + OriginalFace.push_back(CurrF); + + VertexType *v1=CurrF->V2(currPos.E()); + VertexType *v2=CurrF->V1(currPos.E()); + CoordType P1=v1->P(); + CoordType P2=v2->P(); + assert(v1!=v2); + assert(v0!=v1); + assert(v0!=v2); + + SubdivideTris(vcg::Triangle3(P0,P1,P2),0,SubFaces,numSub); + } + } + + static void InterpolateCrossSubTriangles(const std::vector &SubFaces, + const std::vector &OriginalFace, + std::vector &Dir, + std::vector &NormSubF) + { + Dir.clear(); + //check + assert(SubFaces.size()>OriginalFace.size()); + + int SubDivisionSize=SubFaces.size()/OriginalFace.size(); + assert(SubDivisionSize>=3); + assert(SubDivisionSize%2==1);//should be odd + int MiddlePos=SubDivisionSize/2+1;//the one in the middle that should preserve face's cross field + + //calculate the interpolation weights and intervals + std::vector > InterpFaces; + std::vector InterpWeights; + + //the one in the middle should spond to one of the + //two faces the rest is interpolated + for (size_t i=0;iMiddlePos) + { + IndexF0=interval; + IndexF1=(interval+1) % OriginalFace.size(); + alpha=1-(sub_int-MiddlePos)/(ScalarType)SubDivisionSize; + } + else //sub_int==MiddlePos + { + IndexF0=interval; + IndexF1=(interval+1) % OriginalFace.size(); + alpha=1; + } + assert((IndexF0>=0)&&(IndexF0=0)&&(IndexF1(F0,F1)); + InterpWeights.push_back(alpha); + } + assert(InterpFaces.size()==InterpWeights.size()); + //then calculate the interpolated cross field + for (size_t i=0;i > TangVect; + std::vector > Norms; + std::vector W; + Norms.push_back(InterpFaces[i].first->N()); + Norms.push_back(InterpFaces[i].second->N()); + CoordType Dir0=InterpFaces[i].first->PD1(); + CoordType Dir1=InterpFaces[i].second->PD1(); + TangVect.push_back(Dir0); + TangVect.push_back(Dir1); + ScalarType CurrW=InterpWeights[i]; + W.push_back(CurrW); + W.push_back(1-CurrW); + + CoordType TargetN=InterpFaces[i].first->N(); + if (CurrW<0.5) + TargetN=InterpFaces[i].second->N(); + + NormSubF.push_back(TargetN); + //CoordType TargetN=vcg::Normal(SubFaces[i].P(0),SubFaces[i].P(1),SubFaces[i].P(2)); + TargetN.Normalize(); + CoordType InterpD0=InterpolateNRosy3D(TangVect,Norms,W,4,TargetN); + //CoordType InterpD0=Dir1; + //if (CurrW>0.5)InterpD0=Dir0; + + Dir.push_back(InterpD0); + } + } + + static ScalarType turn (const CoordType &norm, const CoordType &d0, const CoordType &d1) + { + //first check if they are coplanar up to a certain delta + return ((d0 ^ d1).normalized() * norm); + } + + static void InterpolateDir(const CoordType &Dir0, + const CoordType &Dir1, + const CoordType &Sep0, + const CoordType &Sep1, + const TriangleType &t0, + const TriangleType &t1, + CoordType &Interpolated, + size_t &Face) + { + //find smallest edge + ScalarType smallestE=std::numeric_limits::max(); + for (int j=0;j<3;j++) + { + ScalarType L0=(t0.P0(j)-t0.P1(j)).Norm(); + ScalarType L1=(t1.P0(j)-t1.P1(j)).Norm(); + if (L0 rotation=vcg::RotationMatrix(N0,N1); + + //transform the first triangle + T0Rot.P(1)=rotation*T0Rot.P(1); + T0Rot.P(2)=rotation*T0Rot.P(2); + + //also rotate the directions + CoordType Dir0Rot=rotation*Dir0; + CoordType Dir1Rot=Dir1; + CoordType Sep0Rot=rotation*Sep0; + CoordType Sep1Rot=Sep1; + + //find the interpolation angles + ScalarType Angle0=vcg::Angle(Dir0Rot,Sep0Rot); + ScalarType Angle1=vcg::Angle(Dir1Rot,Sep1Rot); + assert(Angle0>=0); + assert(Angle1>=0); + assert(Angle0<=(M_PI/2)); + assert(Angle1<=(M_PI/2)); + ScalarType alpha=0.5;//Angle0/(Angle0+Angle1); + + //then interpolate the direction + //CoordType Interp=Dir0Rot*(1-alpha)+Dir1Rot*alpha; + Interpolated=Sep0Rot*(1-alpha)+Sep1Rot*alpha; + Interpolated.Normalize(); + Interpolated*=smallestE; + + //then find the triangle which falls into + CoordType bary0,bary1; + bool Inside0=vcg::InterpolationParameters(T0Rot,Interpolated,bary0); + bool Inside1=vcg::InterpolationParameters(T1Rot,Interpolated,bary1); + //assert(Inside0 || Inside1); + if (!(Inside0 || Inside1)) + { + std::cout << "Not Inside " << Interpolated.X() << "," << Interpolated.Y() << "," << Interpolated.Z() << std::endl; + std::cout << "bary0 " << bary0.X() << "," << bary0.Y() << "," << bary0.Z() << std::endl; + std::cout << "bary1 " << bary1.X() << "," << bary1.Y() << "," << bary1.Z() << std::endl; + std::cout << "Diff0 " << fabs(bary0.Norm() - 1) << std::endl; + std::cout << "Diff1 " << fabs(bary1.Norm() - 1) << std::endl; + } + + if (Inside0) + { + Interpolated=t0.P(0)*bary0.X()+t0.P(1)*bary0.Y()+t0.P(2)*bary0.Z(); + Interpolated-=Origin; + Face=0; + } + else + Face=1; + + //otherwise is already in the tangent space of t0 + Interpolated.Normalize(); + } + + static void ReduceOneDirectionField(std::vector &directions, + std::vector &faces) + { + //compute distribution and find closest pait + std::vector DirProd; + ScalarType MaxProd=-1; + int MaxInd0=-1,MaxInd1=-1; + for (size_t i=0;iMaxProd) + { + MaxProd=currP; + MaxInd0=Index0; + MaxInd1=Index1; + } + prod+=currP; + } + DirProd.push_back(prod); + } + assert(MaxInd0!=MaxInd1); + + //then find the one to be deleted + int IndexDel=MaxInd0; + if (DirProd[MaxInd1]>DirProd[MaxInd0])IndexDel=MaxInd1; + + std::vector SwapV(directions.begin(),directions.end()); + std::vector SwapF(faces.begin(),faces.end()); + + directions.clear(); + for (size_t i=0;i &vPos, + std::vector &directions, + std::vector &faces, + std::vector &WrongTris, + int expVal=-1, + int numSub=3) + { + + directions.clear(); + + //not ct on border + assert(!vPos.IsBorder()); + //the vertex should be the one on the pos + assert(vPos.F()->V(vPos.E())==vPos.V()); + + std::vector SubFaces; + std::vector Dir1,Dir2; + std::vector Norms; + std::vector OriginalFaces; + + //find subfaces + FindSubTriangles(vPos,SubFaces,OriginalFaces,numSub); + + //interpolate the cross field + InterpolateCrossSubTriangles(SubFaces,OriginalFaces,Dir1,Norms); + + //then calculate the orthogonal + for (size_t i=0;iP(); + for (size_t i = 0; i < SubFaces.size(); ++i) + { + TriangleType t0=SubFaces[i]; + TriangleType t1=SubFaces[(i+1)%SubFaces.size()]; + CoordType N0=Norms[i];//vcg::Normal(t0.P(0),t0.P(1),t0.P(2)); + CoordType N1=Norms[(i+1)%Norms.size()];//vcg::Normal(t1.P(0),t1.P(1),t1.P(2)); + N0.Normalize(); + N1.Normalize(); + + std::vector SubDEdges0; + std::vector SubDEdges1; + FindSubDir(t0,0,SubDEdges0,2); + FindSubDir(t1,0,SubDEdges1,2); + CoordType Bary0=SubDEdges0[1]; + CoordType Bary1=SubDEdges1[1]; + //then get the directions to the barycenter + Bary0-=CenterStar; + Bary1-=CenterStar; + Bary0.Normalize(); + Bary1.Normalize(); + + //then get the cross field of the 2 faces + CoordType Dir1F0=Dir1[i]; + CoordType Dir2F0=Dir2[i]; + + assert(fabs(Dir1F0*N0)<0.001); + assert(fabs(Dir1F0*Dir2F0)<0.001); + + if ((Dir1F0*Bary0)<0)Dir1F0=-Dir1F0; + if ((Dir2F0*Bary0)<0)Dir2F0=-Dir2F0; + + CoordType Dir1F1=Dir1[(i+1)%Dir1.size()]; + CoordType Dir2F1=Dir2[(i+1)%Dir2.size()]; + + assert(fabs(Dir1F1*N1)<0.001); + assert(fabs(Dir1F1*Dir2F1)<0.01); + + //find the most similar rotation of the cross field + vcg::Matrix33 rotation=vcg::RotationMatrix(N0,N1); + CoordType Dir1F0R=rotation*Dir1F0; + CoordType Dir2F0R=rotation*Dir2F0; + + //then get the closest upf to K*PI/2 rotations + Dir1F1=vcg::tri::CrossField::K_PI(Dir1F1,Dir1F0R,N1); + Dir2F1=vcg::tri::CrossField::K_PI(Dir2F1,Dir2F0R,N1); + + //then check if cross the direction of the barycenter + ScalarType versef0D1 = turn(N0, Bary0, Dir1F0); + ScalarType versef1D1 = turn(N1, Bary1, Dir1F1); + ScalarType versef0D2 = turn(N0, Bary0, Dir2F0); + ScalarType versef1D2 = turn(N1, Bary1, Dir2F1); + + if ((Bary0*Bary1)<0) + { + + WrongTris.push_back(t0); + WrongTris.push_back(t1); + } + + + CoordType InterpDir; + size_t tri_Index; + if ((versef0D1 * versef1D1) < ScalarType(0)) + { + + InterpolateDir(Dir1F0,Dir1F1,Bary0,Bary1,t0,t1,InterpDir,tri_Index); + + } + else + if ((versef0D2 * versef1D2 )< ScalarType(0)) + { + directions.push_back(InterpolateDir(Dir2F0,Dir2F1,Bary0,Bary1,t0,t1)); + } + + //retrieve original face + assert((tri_Index==0)||(tri_Index==1)); + int OrigFIndex=((i+tri_Index)%SubFaces.size())/numSub; + assert(OrigFIndex>=0); + assert(OrigFIndex CrossToUV(FaceType &f) + static vcg::Point2 CrossToUV(FaceType &f,int numD=0) { typedef typename FaceType::ScalarType ScalarType; typedef typename FaceType::CoordType CoordType; - CoordType Curv=CrossVector(f,0); + CoordType Curv=CrossVector(f,numD); Curv.Normalize(); CoordType bary3d=(f.P(0)+f.P(1)+f.P(2))/3.0;