1 #ifndef _ISO_TANGENTSPACE 2 #define _ISO_TANGENTSPACE 3 #include "iso_parametrization.h" 4 #include <vcg/complex/algorithms/update/curvature.h> 5 #include <vcg/complex/algorithms/update/normal.h> 6 7 class TangentSpace{ 8 typedef IsoParametrization::CoordType CoordType; 9 typedef IsoParametrization::ScalarType ScalarType; 10 11 12 vcg::SimpleTempData<typename ParamMesh::VertContainer,vcg::Matrix33<ScalarType> > *ProjMatrix; 13 14 public: 15 IsoParametrization *isoParam; 16 17 isoParamTheta(int i,vcg::Point2<ScalarType> p,vcg::Point3<ScalarType> & res)18 bool isoParamTheta(int i, vcg::Point2<ScalarType> p, vcg::Point3<ScalarType> &res) const{ 19 return isoParam->Theta(i,p,res); 20 } 21 22 // Theta(int i,const vcg::Point2<ScalarType> & UV,CoordType & pos3D)23 void Theta(int i, 24 const vcg::Point2<ScalarType> &UV, 25 CoordType &pos3D){ 26 isoParamTheta(i,UV,pos3D); 27 } 28 Theta(const int & I,const vcg::Point2<ScalarType> & alpha_beta,std::vector<ParamFace * > & face,std::vector<CoordType> & baryVal)29 bool Theta(const int &I, 30 const vcg::Point2<ScalarType> &alpha_beta, // alphaBeta 31 std::vector<ParamFace*> &face, 32 std::vector<CoordType> &baryVal) 33 { 34 int ret=isoParam->Theta(I,alpha_beta,face,baryVal); 35 return (ret!=-1); 36 } 37 38 ///initialize the sampler 39 void Init(IsoParametrization *_isoParam,ScalarType radius=(ScalarType)0.1) 40 { 41 isoParam=_isoParam; 42 ProjMatrix = new vcg::SimpleTempData<typename ParamMesh::VertContainer,vcg::Matrix33<ScalarType> > (isoParam->ParaMesh()->vert); 43 44 InitProjectionMatrix(radius); 45 vcg::tri::UpdateNormals<ParamMesh>::PerFaceNormalized(*isoParam->ParaMesh()); 46 vcg::tri::UpdateCurvature<ParamMesh>::PrincipalDirectionsNormalCycles(*isoParam->ParaMesh()); 47 } 48 49 ///given an initial position in parametric space (I0,bary0) 50 ///and a 2D vector (vect) expressed in parametric space modify the final 51 ///position (I1,bary1) and return true if everything was ok, false otherwise Sum(const int & I0,const vcg::Point2<ScalarType> & bary0,const vcg::Point2<ScalarType> & vect,int & I1,vcg::Point2<ScalarType> & bary1,int & domain)52 bool Sum(const int &I0,const vcg::Point2<ScalarType> &bary0, 53 const vcg::Point2<ScalarType> &vect, 54 int &I1,vcg::Point2<ScalarType> &bary1,int &domain) const 55 { 56 57 58 vcg::Point2<ScalarType> dest=bary0+vect; 59 //vect[0]*Xaxis + vect[1]*Yaxis; 60 ScalarType alpha=dest.X(); 61 ScalarType beta=dest.Y(); 62 ///point inside the face 63 if ((alpha>=0)&&(alpha<=1)&&(beta>=0)&&(beta<=1)&&((alpha+beta)<=1)) 64 { 65 bary1=dest; 66 I1=I0; 67 domain=0; 68 return true; 69 } 70 71 ///control edges 72 int edge=-1; 73 if ((alpha<=1)&&(beta<=1)&&((alpha+beta)>=1)) 74 edge=0; 75 else 76 if ((alpha<=0)&&(beta<=1)&&((alpha+beta)>=0)) 77 edge=1; 78 else 79 if ((alpha<=1)&&(beta<=0)&&((alpha+beta)>=0)) 80 edge=2; 81 if (edge!=-1) 82 { 83 int DiamIndex=isoParam->GetDiamond(I0,edge); 84 vcg::Point2<ScalarType> UVDiam; 85 ///transform to diamond coordinates 86 isoParam->GE1(I0,dest,DiamIndex,UVDiam); 87 ///trasform back to I,alpha,beta 88 isoParam->inv_GE1(DiamIndex,UVDiam,I1,bary1); 89 domain=1; 90 return true; 91 } 92 int star=-1; 93 ScalarType gamma=(1-alpha-beta); 94 if ((alpha>beta)&&(alpha>gamma)) 95 star=0; 96 else 97 if ((beta>alpha)&&(beta>gamma)) 98 star=1; 99 else 100 star=2; 101 102 ///get the index of star 103 int StarIndex=isoParam->GetStarIndex(I0,star); 104 vcg::Point2<ScalarType> UVHstar; 105 ///transform to UV 106 bool found=isoParam->GE0(I0,dest,StarIndex,UVHstar); 107 ///trasform back to I,alpha,beta 108 if (!found) 109 return false; 110 found=isoParam->inv_GE0(StarIndex,UVHstar,I1,bary1); 111 /*AbstractFace* f0=&isoParam->AbsMesh()->face[I0]; 112 AbstractFace* f1=&isoParam->AbsMesh()->face[I1]; 113 AbstractVertex *v0=f0->V(0); 114 AbstractVertex *v1=f0->V(1); 115 AbstractVertex *v2=f0->V(2); 116 AbstractVertex *v3=f1->V(0); 117 AbstractVertex *v4=f1->V(1); 118 AbstractVertex *v5=f1->V(2);*/ 119 domain=2; 120 return found; 121 } 122 123 ///given two positions in parametric space (I0,bary0) and (I1,bary1) 124 ///modify the 2D vector (vect) and return true if everything was ok, false otherwise Sub(const int & I0,const vcg::Point2<ScalarType> & bary0,const int & I1,const vcg::Point2<ScalarType> & bary1,vcg::Point2<ScalarType> & vect,int & num)125 bool Sub(const int &I0,const vcg::Point2<ScalarType> &bary0, 126 const int &I1,const vcg::Point2<ScalarType> &bary1, 127 vcg::Point2<ScalarType> &vect,int &num) const 128 { 129 int IndexDomain; 130 num=isoParam->InterpolationSpace(I0,I1,IndexDomain); 131 ///is a face 132 if (num==0) 133 { 134 //printf("F"); 135 assert(I0==I1); 136 vect=bary0-bary1; 137 return true; 138 } 139 else 140 ///a diamond 141 if (num==1) 142 { 143 //printf("D"); 144 ///transform in UV space 145 vcg::Point2<ScalarType> UVDiam; 146 isoParam->GE1(I1,bary1,IndexDomain,UVDiam); 147 ///then find bary coords which respect to the first face 148 vcg::Point2<ScalarType> bary2; 149 isoParam->inv_GE1_fixedI(IndexDomain,UVDiam,I0,bary2); 150 vect=bary0-bary2; 151 return true; 152 } 153 else 154 ///a star 155 if (num==2) 156 { 157 //printf("S"); 158 ///transform in UV space 159 vcg::Point2<ScalarType> UVStar; 160 isoParam->GE0(I1,bary1,IndexDomain,UVStar); 161 ///then find bary coords which respect to the first face 162 vcg::Point2<ScalarType> bary2; 163 isoParam->inv_GE0_fixedI(IndexDomain,UVStar,I0,bary2); 164 vect=bary0-bary2; 165 166 return true; 167 } 168 else 169 return false; 170 } 171 172 bool TangentDir(const int &I,const vcg::Point2<ScalarType> &bary, 173 CoordType &XAxis,CoordType &YAxis,ScalarType radius=(ScalarType)0.5) const 174 { 175 ///3d position of origin 176 //CoordType origin; 177 //isoParam->Theta(I,bary,origin); 178 179 // two axis in alpha-beta space that will be orthogonal in UV-Space 180 const ScalarType h=(ScalarType)2.0/sqrt((ScalarType)3.0); 181 vcg::Point2<ScalarType> XP=vcg::Point2<ScalarType>(1,0); 182 vcg::Point2<ScalarType> YP=vcg::Point2<ScalarType>(-0.5,h); 183 184 XP*=radius; 185 YP*=radius; 186 vcg::Point2<ScalarType> XM=-XP; 187 vcg::Point2<ScalarType> YM=-YP; 188 //Yaxis.Normalize(); 189 190 vcg::Point2<ScalarType> bary0,bary1,bary2,bary3; 191 int I0,I1,I2,I3; 192 CoordType X0,X1,Y0,Y1; 193 ///find paraCoords of four neighbors 194 int domain; 195 bool done=true; 196 done&=Sum(I,bary,XP,I0,bary0,domain); 197 done&=Sum(I,bary,XM,I1,bary1,domain); 198 done&=Sum(I,bary,YP,I2,bary2,domain); 199 done&=Sum(I,bary,YM,I3,bary3,domain); 200 if (!done) 201 return false; 202 203 //if (I0!=I1 || I1!=I2 || I2!=I3) return false; 204 205 ///get 3d position 206 isoParamTheta(I0,bary0,X0); 207 isoParamTheta(I1,bary1,X1); 208 isoParamTheta(I2,bary2,Y0); 209 isoParamTheta(I3,bary3,Y1); 210 211 212 ///average .. considering one is opposite respect to the other 213 XAxis=(X0-X1)/(ScalarType)2.0; 214 YAxis=(Y0-Y1)/(ScalarType)2.0; 215 ///final scaling 216 XAxis/=radius; 217 YAxis/=radius; 218 return true; 219 } 220 221 void Test(int Sample=100,int Ite=100) 222 { 223 int max=isoParam->AbsMesh()->face.size(); 224 for (int I=0;I<max;I++) 225 { 226 printf("\n TESTING %d \n",I); 227 for (int i=0;i<Sample;i++) 228 { 229 for (int j=0;j<Sample;j++) 230 { 231 if ((i+j)<Sample) 232 { 233 ScalarType alpha=(ScalarType)i/(ScalarType)Sample; 234 ScalarType beta=(ScalarType)j/(ScalarType)Sample; 235 assert((alpha+beta)<=1.0); 236 vcg::Point2<ScalarType> bary=vcg::Point2<ScalarType>(alpha,beta); 237 assert((bary.X()>=0)&&(bary.X()<=1)); 238 assert((bary.Y()>=0)&&(bary.Y()<=1)); 239 assert((bary.X()+bary.Y()<=1)); 240 241 for (int k=0;k<Ite;k++) 242 { 243 ScalarType d0=((ScalarType)(rand()%1000)); 244 ScalarType d1=((ScalarType)(rand()%1000)); 245 vcg::Point2<ScalarType> vect=vcg::Point2<ScalarType>(d0,d1); 246 vect.Normalize(); 247 ScalarType norm=(ScalarType)0.05;//((ScalarType)(rand()%1000))/(ScalarType)2000; 248 assert(norm<1); 249 vect*=norm; 250 int I1; 251 vcg::Point2<ScalarType> bary1; 252 vcg::Point2<ScalarType> vect1; 253 int domain0; 254 bool b1=Sum(I,bary,vect,I1,bary1,domain0); 255 assert(b1); 256 257 assert((bary1.X()>=0)&&(bary1.X()<=1)); 258 assert((bary1.Y()>=0)&&(bary1.Y()<=1)); 259 if(!((bary1.X()+bary1.Y())<=1.00001)) 260 { 261 printf("\n SUM %.4f \n",bary1.X()+bary1.Y()); 262 assert(0); 263 } 264 assert(I1<max); 265 int domain; 266 bool b2=Sub(I,bary,I1,bary1,vect1,domain); 267 assert(b2); 268 ScalarType diff=(vect1+vect).Norm(); 269 if (domain0==0) 270 assert(domain==0); 271 if ((domain0==1)&&(domain==2)) 272 assert(0); 273 /*if (domain0!=domain) 274 assert(0);*/ 275 if (diff>0.001) 276 { 277 printf("\n DIFF %.4f domain SUM %d domain SUB %d \n",diff,domain0,domain); 278 //assert(0); 279 } 280 //assert(fabs(vect1.X()-vect.X())<0.0001); 281 //assert(fabs(vect1.Y()-vect.Y())<0.0001); 282 283 } 284 } 285 } 286 } 287 } 288 } 289 290 void InitProjectionMatrix(ScalarType radius=(ScalarType)0.1) 291 { 292 for (int i=0;i<isoParam->ParaMesh()->vert.size();i++) 293 { 294 int I=isoParam->ParaMesh()->vert[i].T().N(); 295 vcg::Point2<ScalarType> bary=isoParam->ParaMesh()->vert[i].T().P(); 296 297 CoordType origin; 298 isoParamTheta(I,bary,origin); 299 CoordType XAxis,YAxis; // tangent axis in 3D Object space 300 ///get tangent directions 301 302 bool done=TangentDir(I,bary,XAxis,YAxis,(ScalarType)0.1); 303 if (!done) 304 { 305 (*ProjMatrix)[i].SetIdentity(); 306 continue; 307 } 308 309 // must compute res2d such that: res2d.X() * XAxis + res2d.Y() * YAxis + dontCare * ZAxis = vect3d 310 CoordType ZAxis = -(XAxis^YAxis).Normalize()*XAxis.Norm(); 311 312 (*ProjMatrix)[i].SetColumn(0,XAxis); 313 (*ProjMatrix)[i].SetColumn(1,YAxis); 314 (*ProjMatrix)[i].SetColumn(2,ZAxis); 315 vcg::Invert((*ProjMatrix)[i]); 316 } 317 } 318 GetCurvature(const int & I,const vcg::Point2<ScalarType> & alpha_beta,CoordType & d1,CoordType & d2,ScalarType & k1,ScalarType & k2)319 void GetCurvature(const int &I,const vcg::Point2<ScalarType> &alpha_beta, 320 CoordType &d1,CoordType &d2,ScalarType &k1,ScalarType &k2) 321 { 322 ParamFace* face=NULL; 323 CoordType baryVal; 324 isoParam->Theta(I,alpha_beta,face,baryVal); 325 ParamVertex *v0=face->V(0); 326 ParamVertex *v1=face->V(1); 327 ParamVertex *v2=face->V(2); 328 d1=v0->PD1()*baryVal.X()+v1->PD1()*baryVal.Y()+v2->PD1()*baryVal.Z(); 329 d2=v0->PD2()*baryVal.X()+v1->PD2()*baryVal.Y()+v2->PD2()*baryVal.Z(); 330 k1=v0->K1()*baryVal.X()+v1->K1()*baryVal.Y()+v2->K1()*baryVal.Z(); 331 k2=v0->K2()*baryVal.X()+v1->K2()*baryVal.Y()+v2->K2()*baryVal.Z(); 332 } 333 GetProjectionMatrix(const int & I,const vcg::Point2<ScalarType> & alpha_beta,vcg::Matrix33<ScalarType> & projMatr)334 void GetProjectionMatrix(const int &I,const vcg::Point2<ScalarType> &alpha_beta, 335 vcg::Matrix33<ScalarType> &projMatr) const 336 { 337 ParamFace* face=NULL; 338 CoordType baryVal; 339 int dom=isoParam->Theta(I,alpha_beta,face,baryVal); 340 if (dom==-1) 341 { 342 projMatr.SetIdentity(); 343 return; 344 } 345 ParamVertex *v0=face->V(0); 346 ParamVertex *v1=face->V(1); 347 ParamVertex *v2=face->V(2); 348 349 projMatr=(*ProjMatrix)[v0]*baryVal.X(); 350 projMatr+=(*ProjMatrix)[v1]*baryVal.Y(); 351 projMatr+=(*ProjMatrix)[v2]*baryVal.Z(); 352 } 353 354 Project(const int & I,const vcg::Point2<ScalarType> & bary,const CoordType & vect3d,vcg::Point2<ScalarType> & res2d)355 bool Project(const int &I,const vcg::Point2<ScalarType> &bary, 356 const CoordType &vect3d, // in object space 357 vcg::Point2<ScalarType> &res2d) const 358 { 359 vcg::Matrix33f m; 360 GetProjectionMatrix(I,bary,m); 361 362 ScalarType deltaX=vect3d*m.GetRow(0); 363 ScalarType deltaY=vect3d*m.GetRow(1); 364 365 // two axis in alpha-beta space that will be orthogonal in UV-Space 366 const ScalarType h=(ScalarType)2.0/sqrt((ScalarType)3.0); 367 vcg::Point2<ScalarType> XP=vcg::Point2<ScalarType>(1,0); 368 vcg::Point2<ScalarType> YP=vcg::Point2<ScalarType>(-0.5,h); 369 res2d = XP*deltaX + YP*deltaY; 370 //res2d=vcg::Point2<ScalarType>(deltaX,deltaY); 371 372 return true; 373 } 374 375 // WEIGHTED INTERPOLATION OF POINTS IN TANGENT SPACE 376 /// weights MUST be normalized Interpolate(const int & I0,const vcg::Point2<ScalarType> & alpha_beta0,const int & I1,const vcg::Point2<ScalarType> & alpha_beta1,ScalarType weight,int & I_res,vcg::Point2<ScalarType> & alpha_beta_res)377 bool Interpolate(const int &I0,const vcg::Point2<ScalarType> &alpha_beta0, 378 const int &I1,const vcg::Point2<ScalarType> &alpha_beta1, 379 ScalarType weight, 380 int &I_res, 381 vcg::Point2<ScalarType> &alpha_beta_res) 382 { 383 int IndexDomain; 384 int kind=isoParam->InterpolationSpace(I0,I1,IndexDomain); 385 if (kind==-1) 386 return false; 387 388 vcg::Point2<ScalarType> transformed0,transformed1; 389 390 ///interpolate in a face 391 if (kind==0) 392 { 393 isoParam->GE2(IndexDomain,alpha_beta0,transformed0); 394 isoParam->GE2(IndexDomain,alpha_beta1,transformed1); 395 } 396 else 397 ///interpolate in a diamond 398 if (kind==1) 399 { 400 isoParam->GE1(I0,alpha_beta0,IndexDomain,transformed0); 401 isoParam->GE1(I1,alpha_beta1,IndexDomain,transformed1); 402 } 403 else 404 { 405 isoParam->GE0(I0,alpha_beta0,IndexDomain,transformed0); 406 isoParam->GE0(I1,alpha_beta1,IndexDomain,transformed1); 407 } 408 409 vcg::Point2<ScalarType> UV_interp=transformed0*weight+transformed1*(1.0-weight); 410 ///FINALLY...... 411 ///transform back to alpha,beta,I 412 if (kind==0) 413 { 414 isoParam->inv_GE2(IndexDomain,UV_interp,alpha_beta_res); 415 I_res=IndexDomain; 416 } 417 else 418 if (kind==1) 419 { 420 isoParam->inv_GE1(IndexDomain,UV_interp,I_res,alpha_beta_res); 421 } 422 else 423 isoParam->inv_GE0(IndexDomain,UV_interp,I_res,alpha_beta_res); 424 return true; 425 } 426 AbstractArea()427 ScalarType AbstractArea() 428 { 429 ScalarType Cnum=sqrt(3.0)/4.0; 430 return(isoParam->AbsMesh()->fn*Cnum); 431 } 432 433 // WEIGHTED INTERPOLATION OF POINTS IN TANGENT SPACE 434 /// weights MUST be normalized 435 bool Interpolate(const std::vector<int> &I, 436 const std::vector<vcg::Point2<ScalarType> > &alpha_beta, 437 const std::vector<ScalarType> &weights, 438 int &I_res, 439 vcg::Point2<ScalarType> &alpha_beta_res, 440 int *num=NULL) 441 { 442 int size; 443 if (num==NULL) 444 size=alpha_beta.size(); 445 else 446 size=*num; 447 448 int IndexDomain; 449 int kind=isoParam->InterpolationSpace(I,IndexDomain,num); 450 if (kind==-1) 451 return false; 452 453 std::vector<vcg::Point2<ScalarType> > transformed; 454 transformed.resize(size); 455 456 ///interpolate in a face 457 if (kind==0) 458 { 459 for (int i=0;i<size;i++) 460 isoParam->GE2(IndexDomain,alpha_beta[i],transformed[i]); 461 } 462 else 463 ///interpolate in a diamond 464 if (kind==1) 465 { 466 for (int i=0;i<size;i++) 467 isoParam->GE1(I[i],alpha_beta[i],IndexDomain,transformed[i]); 468 } 469 else 470 { 471 for (int i=0;i<size;i++) 472 bool b0=isoParam->GE0(I[i],alpha_beta[i],IndexDomain,transformed[i]); 473 } 474 475 /// do the interpolation 476 vcg::Point2<ScalarType> UV_interp=vcg::Point2<ScalarType>(0,0); 477 for (int i=0;i<size;i++) 478 UV_interp+=(transformed[i]*weights[i]); 479 480 ///FINALLY...... 481 ///transform back to alpha,beta,I 482 if (kind==0) 483 { 484 isoParam->inv_GE2(IndexDomain,UV_interp,alpha_beta_res); 485 I_res=IndexDomain; 486 } 487 else 488 if (kind==1) 489 { 490 isoParam->inv_GE1(IndexDomain,UV_interp,I_res,alpha_beta_res); 491 } 492 else 493 isoParam->inv_GE0(IndexDomain,UV_interp,I_res,alpha_beta_res); 494 return true; 495 } 496 497 }; 498 #endif 499