00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #ifndef FITMAPS_H
00025 #define FITMAPS_H
00026
00027 #include <vcg/math/histogram.h>
00028
00029 #include <vcg/simplex/face/jumping_pos.h>
00030 #include <vcg/complex/algorithms/update/flag.h>
00031 #include <vcg/complex/algorithms/update/normal.h>
00032 #include <vcg/complex/algorithms/update/curvature.h>
00033 #include <vcg/complex/algorithms/update/topology.h>
00034 #include <vcg/complex/algorithms/update/bounding.h>
00035 #include "vcg/complex/algorithms/update/curvature_fitting.h"
00036
00037 #include <eigenlib/Eigen/Core>
00038 #include <eigenlib/Eigen/QR>
00039 #include <eigenlib/Eigen/LU>
00040 #include <eigenlib/Eigen/SVD>
00041
00042 #include <vcg/complex/algorithms/nring.h>
00043
00044 #include <vcg/complex/algorithms/smooth.h>
00045
00046
00047 using namespace Eigen;
00048
00049 namespace vcg { namespace tri {
00050
00051 template<class MeshType>
00052 class Fitmaps
00053 {
00054 public:
00055
00056 typedef typename MeshType::FaceType FaceType;
00057 typedef typename MeshType::VertexType VertexType;
00058 typedef typename MeshType::ScalarType ScalarType;
00059 typedef typename MeshType::FaceIterator FaceIterator;
00060 typedef typename MeshType::VertexIterator VertexIterator;
00061 typedef typename MeshType::CoordType CoordType;
00062
00063 typedef vcg::tri::Nring<MeshType> RingWalker;
00064
00065
00066 class Bicubic
00067 {
00068 public:
00069
00070 Bicubic() {};
00071
00072 Bicubic(vector<double>& input)
00073 {
00074 data = input;
00075 if (input.size() != 16)
00076 {
00077 assert(0);
00078 }
00079 }
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100 double& a() { return data[0];}
00101 double& b() { return data[1];}
00102 double& c() { return data[2];}
00103 double& d() { return data[3];}
00104 double& e() { return data[4];}
00105 double& f() { return data[5];}
00106 double& g() { return data[6];}
00107 double& h() { return data[7];}
00108 double& i() { return data[8];}
00109 double& l() { return data[9];}
00110 double& m() { return data[10];}
00111 double& n() { return data[11];}
00112 double& o() { return data[12];}
00113 double& p() { return data[13];}
00114 double& q() { return data[14];}
00115 double& r() { return data[15];}
00116
00117 vector<double> data;
00118
00119 double evaluate(double u, double v)
00120 {
00121
00122 return
00123 a() * u*u*u * v*v*v +
00124 b() * u*u*u * v*v +
00125 c() * u*u*u * v +
00126 d() * u*u*u +
00127 e() * u*u * v*v*v +
00128 f() * u*u * v*v +
00129 g() * u*u * v +
00130 h() * u*u +
00131 i() * u * v*v*v +
00132 l() * u * v*v +
00133 m() * u * v +
00134 n() * u +
00135 o() * v*v*v +
00136 p() * v*v +
00137 q() * v +
00138 r();
00139 }
00140
00141
00142 double distanceRMS(std::vector<CoordType>& VV)
00143 {
00144 double error = 0;
00145 for(typename std::vector<CoordType>::iterator it = VV.begin(); it != VV.end(); ++it)
00146 {
00147 double u = it->X();
00148 double v = it->Y();
00149 double n = it->Z();
00150
00151 double temp = evaluate(u,v);
00152 error += (n - temp)*(n - temp);
00153 }
00154 error /= (double) VV.size();
00155 return sqrt(error);
00156 }
00157
00158 static Bicubic fit(std::vector<CoordType> VV)
00159 {
00160 assert(VV.size() >= 16);
00161 Eigen::MatrixXf A(VV.size(),16);
00162 Eigen::MatrixXf b(VV.size(),1);
00163 Eigen::MatrixXf sol(16,1);
00164
00165 for(unsigned int c=0; c < VV.size(); ++c)
00166 {
00167 double u = VV[c].X();
00168 double v = VV[c].Y();
00169 double n = VV[c].Z();
00170
00171 A(c,0) = u*u*u * v*v*v;
00172 A(c,1) = u*u*u * v*v;
00173 A(c,2) = u*u*u * v;
00174 A(c,3) = u*u*u;
00175 A(c,4) = u*u * v*v*v;
00176 A(c,5) = u*u * v*v;
00177 A(c,6) = u*u * v;
00178 A(c,7) = u*u;
00179 A(c,8) = u * v*v*v;
00180 A(c,9) = u * v*v;
00181 A(c,10) = u * v;
00182 A(c,11) = u;
00183 A(c,12) = v*v*v;
00184 A(c,13) = v*v;
00185 A(c,14) = v;
00186 A(c,15) = 1;
00187
00188
00189 b[c] = n;
00190 }
00191
00192 A.svd().solve(b, &sol);
00193
00194 vector<double> r(16);
00195
00196 for (int i=0; i < 16; ++i)
00197 r.at(i) = sol[i];
00198
00199 return Bicubic(r);
00200 }
00201 };
00202
00203 Fitmaps()
00204 {}
00205
00206 class radSorter
00207 {
00208 public:
00209
00210 radSorter(VertexType* v)
00211 {
00212 origin = v;
00213 }
00214
00215 VertexType* origin;
00216
00217 bool operator() (VertexType* v1, VertexType* v2)
00218 {
00219 return (v1->P() - origin->P()).SquaredNorm() < (v2->P() - origin->P()).SquaredNorm();
00220 }
00221 };
00222
00223 float getMeanCurvature(VertexType* vp)
00224 {
00225 return (vp->K1() + vp->K2())/2.0;
00226 }
00227
00228 static bool fitBicubicPoints(VertexType* v, std::vector<CoordType>& ref, Bicubic& ret, std::vector<CoordType>& points, std::vector<VertexType*>& ring)
00229 {
00230 points.clear();
00231
00232 if (ring.size() < 16)
00233 {
00234 return false;
00235 }
00236
00237 typename std::vector<VertexType*>::iterator b = ring.begin();
00238 typename std::vector<VertexType*>::iterator e = ring.end();
00239
00240 while(b != e)
00241 {
00242 CoordType vT = (*b)->P() - v->P();
00243
00244 double x = vT * ref[0];
00245 double y = vT * ref[1];
00246 double z = vT * ref[2];
00247
00248 points.push_back(CoordType(x,y,z));
00249 ++b;
00250 }
00251 ret = Bicubic::fit(points);
00252 return true;
00253 }
00254
00255 static double AverageEdgeLenght(MeshType& m)
00256 {
00257 double doubleA = 0;
00258 for (FaceIterator fi = m.face.begin(); fi!=m.face.end(); fi++) if (!fi->IsD()) {
00259 doubleA+=vcg::DoubleArea(*fi);
00260 }
00261 int nquads = m.fn / 2;
00262 return sqrt( doubleA/(2*nquads) );
00263 }
00264
00265 static void computeMFitmap(MeshType& m, float perc, int ringMax = 50)
00266 {
00267 vcg::tri::UpdateCurvatureFitting<MeshType>::computeCurvature(m);
00268 vcg::tri::UpdateNormal<MeshType>::PerVertexAngleWeighted(m);
00269
00270 vcg::tri::UpdateTopology<MeshType>::FaceFace(m);
00271 vcg::tri::UpdateTopology<MeshType>::VertexFace(m);
00272
00273 vcg::tri::UpdateBounding<MeshType>::Box(m);
00274
00275 int countTemp = 0;
00276
00277 RingWalker::clearFlags(&m);
00278
00279 for(VertexIterator it=m.vert.begin(); it!=m.vert.end();++it)
00280 {
00281 if ((countTemp++ % 100) == 0)
00282 cerr << countTemp << "/" << m.vert.size() << endl;
00283
00284 RingWalker rw(&*it,&m);
00285
00286 CoordType nor = it->N();
00287
00288 float okfaces = 0;
00289 float flipfaces = 0;
00290
00291 int count = 0;
00292 do
00293 {
00294 count++;
00295 rw.expand();
00296 for(unsigned i=0; i<rw.lastF.size();++i)
00297 {
00298 CoordType vet1 = nor;
00299 CoordType vet2 = rw.lastF[i]->N();
00300
00301 vet1.Normalize();
00302 vet2.Normalize();
00303
00304
00305 double scal = vet1 * vet2;
00306 if ((scal) > 0)
00307 okfaces += (vcg::DoubleArea(*rw.lastF[i]));
00308 else
00309 flipfaces += (vcg::DoubleArea(*rw.lastF[i]));
00310 }
00311 } while ((((flipfaces)/(okfaces + flipfaces))*100.0 < perc) && (count < ringMax));
00312
00313 std::sort(rw.lastV.begin(),rw.lastV.end(),radSorter(&*it));
00314
00315 it->Q() = ((*rw.lastV.begin())->P() - it->P()).Norm();
00316 rw.clear();
00317
00318 }
00319
00320 vcg::tri::Smooth<MeshType>::VertexQualityLaplacian(m,2);
00321 }
00322
00323 static vector<VertexType*> gatherNeighsSurface(VertexType* vt, float sigma, MeshType& m)
00324 {
00325 vector<VertexType*> current;
00326
00327 RingWalker rw(vt,&m);
00328
00329 bool exit = false;
00330
00331 do
00332 {
00333 rw.expand();
00334
00335 exit = true;
00336
00337 for(typename vector<VertexType*>::iterator it = rw.lastV.begin(); it != rw.lastV.end(); ++it)
00338 {
00339 if (((*it)->P() - vt->P()).Norm() < sigma)
00340 {
00341 current.push_back(*it);
00342 exit = false;
00343 }
00344 }
00345
00346 } while (!exit);
00347
00348 rw.clear();
00349 return current;
00350 }
00351
00352 static void computeSFitmap(MeshType& m)
00353 {
00354
00355 vcg::tri::UpdateCurvatureFitting<MeshType>::computeCurvature(m);
00356 vcg::tri::UpdateNormal<MeshType>::PerVertexAngleWeighted(m);
00357
00358 vcg::tri::UpdateTopology<MeshType>::FaceFace(m);
00359 vcg::tri::UpdateTopology<MeshType>::VertexFace(m);
00360
00361
00362 vcg::tri::UpdateBounding<MeshType>::Box(m);
00363
00364 int countTemp = 0;
00365
00366 double e = AverageEdgeLenght(m);
00367
00368 int iteraz = 5;
00369
00370 for(VertexIterator it=m.vert.begin(); it!=m.vert.end();++it)
00371 {
00372 if ((countTemp++ % 100) == 0)
00373 cerr << countTemp << "/" << m.vert.size() << endl;
00374
00375 vector<float> oneX;
00376
00377
00378 for (int iteration = 0; iteration<iteraz; ++iteration)
00379 {
00380 oneX.push_back((iteration+1)*(e));
00381 }
00382
00383 std::vector<CoordType> ref(3);
00384 ref[0] = it->PD1();
00385 ref[1] = it->PD2();
00386 ref[2] = it->PD1() ^ it->PD2();
00387
00388 ref[0].Normalize();
00389 ref[1].Normalize();
00390 ref[2].Normalize();
00391
00392 Bicubic b;
00393
00394 RingWalker::clearFlags(&m);
00395
00396 std::vector<VertexType*> pointsGlobal = gatherNeighsSurface(&*it,oneX.at(iteraz-1),m);
00397
00398 vector<float> onedimensional;
00399
00400 for (int iteration = 0; iteration<iteraz; ++iteration)
00401 {
00402 std::vector<VertexType*> points;
00403
00404 std::vector<CoordType> projected;
00405 for (typename std::vector<VertexType*>::iterator it2 = pointsGlobal.begin(); it2 != pointsGlobal.end(); ++it2)
00406 {
00407 if (((*it).P() - (*it2)->P()).Norm() < oneX.at(iteration))
00408 points.push_back(*it2);
00409 }
00410
00411 std::vector<VertexType*>& pointsFitting = points;
00412
00413
00414 if (!fitBicubicPoints(&*it, ref, b, projected,pointsFitting))
00415 {
00416 onedimensional.push_back(0);
00417 }
00418 else
00419 {
00420 onedimensional.push_back(b.distanceRMS(projected));
00421 }
00422
00423 }
00424
00425
00426
00427 Eigen::MatrixXf Am(onedimensional.size(),1);
00428 Eigen::MatrixXf bm(onedimensional.size(),1);
00429 Eigen::MatrixXf sol(1,1);
00430
00431 for(unsigned int c=0; c < onedimensional.size(); ++c)
00432 {
00433 double x = oneX.at(c);
00434
00435 Am(c,0) = pow(x,4);
00436 bm[c] = onedimensional[c];
00437 }
00438
00439 Am.svd().solve(bm, &sol);
00440
00441 it->Q() = pow((double)sol[0],0.25);
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471 }
00472
00473 vcg::tri::Smooth<MeshType>::VertexQualityLaplacian(m,1);
00474 }
00475
00476
00477 ~Fitmaps(){};
00478
00479 };
00480
00481 }}
00482
00483 #endif // FITMAPS_H