1 #ifndef _OPT_PATCHES
2 #define _OPT_PATCHES
3
4 #include <algorithm>
5 #include <vcg/complex/complex.h>
6
7
8 template <class MeshType>
9 class PatchesOptimizer
10 {
11 typedef typename MeshType::VertexType VertexType;
12 typedef typename MeshType::FaceType FaceType;
13 typedef typename MeshType::CoordType CoordType;
14 typedef typename MeshType::ScalarType ScalarType;
15
16 public:
17 struct minInfoUV
18 {
19 public:
20
21 VertexType* to_optimize;
22 std::vector<VertexType*> Hres_vert;
23 MeshType *parametrized_domain;
24 MeshType *base_domain;
25 MeshType hres_mesh;
26 };
27
28 ScalarType averageArea;
29 ScalarType averageLength;
30 MeshType &base_mesh;
31 MeshType &final_mesh;
32 int global_mark;
33 vcg::SimpleTempData<typename MeshType::VertContainer,int> markers;
34
35 ///energy for equilararity and equiareal minimization
Equi_energy(float * p,float * x,int,int,void * data)36 static void Equi_energy(float *p, float *x, int /*m*/, int/* n*/, void *data)
37 {
38
39 /*const float MaxVal=10000.f;*/
40 minInfoUV &inf = *(minInfoUV *)data;
41
42 ///assign coordinate to central vertex
43 inf.to_optimize->T().U()=p[0];
44 inf.to_optimize->T().V()=p[1];
45
46 ///control that the parametrization is non folded
47 std::vector<FaceType*> folded;
48 bool b=NonFolded<MeshType>(*inf.parametrized_domain,folded);
49 if (!b)
50 {
51 x[0]=std::numeric_limits<float>::max();
52 x[1]=std::numeric_limits<float>::max();
53 return;
54 }
55
56 ////set rest positions for survived vertex
57 ///get the non border one that is the one survived
58 CoordType val;
59 bool found0,found1;
60
61 ///update 3d position of central vertex
62 found0=GetCoordFromUV<MeshType>(inf.hres_mesh,inf.to_optimize->T().U(),inf.to_optimize->T().V(),val,true);
63 if (!found0)
64 found1=GetCoordFromUV<MeshType>(*inf.parametrized_domain,inf.to_optimize->T().U(),inf.to_optimize->T().V(),val,true);
65
66 //assert ((found0)||(found1));
67 if ((found0)||(found1))
68 inf.to_optimize->RPos=val;
69
70 ///clear assigned vertices
71 for (unsigned int i=0;i<inf.parametrized_domain->face.size();i++)
72 inf.parametrized_domain->face[i].vertices_bary.resize(0);
73
74 ///update alphabeta from UV to calculate edge_length and area
75 bool inside=true;
76 for (unsigned int i=0;i<inf.Hres_vert.size();i++)
77 {
78 VertexType *test=inf.Hres_vert[i];
79 ScalarType u=test->T().U();
80 ScalarType v=test->T().V();
81 CoordType bary;
82 int index;
83 inside &=GetBaryFaceFromUV(*inf.parametrized_domain,u,v,bary,index);
84 FaceType* chosen;
85 if (!inside)///ack
86 {
87 chosen=test->father;
88 bary=test->Bary;
89 }
90 else
91 {
92 chosen=&inf.parametrized_domain->face[index];
93 }
94 chosen->vertices_bary.push_back(std::pair<VertexType*,vcg::Point3f>(test,bary));
95 test->father=chosen;
96 assert(!chosen->IsD());
97 test->Bary=bary;
98 }
99
100 if (!inside)///ack
101 {
102 x[0]=std::numeric_limits<float>::max();
103 x[1]=std::numeric_limits<float>::max();
104 return;
105 }
106
107 ScalarType maxEdge=0;
108 ScalarType minEdge=std::numeric_limits<float>::max();
109 ScalarType maxArea=0;
110 ScalarType minArea=std::numeric_limits<float>::max();
111
112 ///find minimum and maximum of estimated area
113 for (unsigned int i=0;i<inf.parametrized_domain->face.size();i++)
114 {
115 ScalarType area=EstimateAreaByParam<FaceType>(&inf.parametrized_domain->face[i]);
116 if (area<minArea)
117 minArea=area;
118 if (area>maxArea)
119 maxArea=area;
120 }
121
122 ///find minimum and maximum of edges
123 for (unsigned int i=0;i<inf.parametrized_domain->vert.size();i++)
124 {
125 VertexType *v0=&inf.parametrized_domain->vert[i];
126 VertexType *v1=inf.to_optimize;
127 if (v0!=v1)
128 {
129 std::vector<typename MeshType::FaceType*> on_edge,faces1,faces2;
130 getSharedFace<MeshType>(v0,v1,on_edge,faces1,faces2);
131 FaceType* edgeF[2];
132 edgeF[0]=on_edge[0];
133 edgeF[1]=on_edge[1];
134 ScalarType length=EstimateLengthByParam<MeshType>(v0,v1,edgeF);
135 if (length<minEdge)
136 minEdge=length;
137 if (length>maxEdge)
138 maxEdge=length;
139 }
140 }
141 //if ((minArea<=0)||(minEdge<=0))
142
143 /*assert(minArea>0);
144 assert(minEdge>0);*/
145 if (minArea==0)
146 minArea=(ScalarType)0.00001;
147 if (minEdge==0)
148 minEdge=(ScalarType)0.00001;
149
150 x[0]=((maxArea/minArea)*(ScalarType)2.0);
151 x[1]=pow(maxEdge/minEdge,(ScalarType)2.0);
152
153 }
154
LengthPath(VertexType * v0,VertexType * v1)155 static ScalarType LengthPath(VertexType *v0,VertexType *v1)
156 {
157 std::vector<FaceType*> on_edge,faces1,faces2;
158 getSharedFace<MeshType>(v0,v1,on_edge,faces1,faces2);
159 FaceType* edgeF[2];
160 edgeF[0]=on_edge[0];
161 edgeF[1]=on_edge[1];
162 ScalarType length=EstimateLengthByParam<FaceType>(v0,v1,edgeF);
163 return length;
164 }
165
166 /////return the priority of vertex processing
167 //ScalarType Priority(VertexType *v)
168 //{
169 // std::vector<typename MeshType::VertexType*> star;
170 // getVertexStar<MeshType>(v,star);
171 // ScalarType prior=0;
172 // for (int i=0;i<star.size();i++)
173 // {
174 // VertexType *v1=star[i];
175 // ScalarType length=LengthPath(v,v1);//EstimateLengthByParam<FaceType>(v0,v1,edgeF);
176 // prior+=pow((length-averageLength),(ScalarType)2);
177 // }
178 // std::vector<VertexType*> vertices;
179 // std::vector<FaceType*> faces;
180 // vertices.push_back(v);
181 // getSharedFace<MeshType>(vertices,faces);
182 // for (int i=0;i<faces.size();i++)
183 // prior+=pow((EstimateAreaByParam<FaceType>(faces[i])-averageArea),2);
184
185 // return prior;
186 //}
187
188 ///return the priority of vertex processing
Priority(VertexType * v)189 ScalarType Priority(VertexType *v)
190 {
191 std::vector<typename MeshType::VertexType*> star;
192 getVertexStar<MeshType>(v,star);
193 ScalarType priorL=0,priorA=0;
194 std::vector<ScalarType> Lengths,Areas;
195 Lengths.resize(star.size());
196
197 std::vector<VertexType*> vertices;
198 std::vector<FaceType*> faces;
199 vertices.push_back(v);
200 getSharedFace<MeshType>(vertices,faces);
201 Areas.resize(faces.size());
202
203 ScalarType aveL=0;
204 ScalarType aveA=0;
205 for (unsigned int i=0;i<star.size();i++)
206 {
207 VertexType *v1=star[i];
208 ScalarType length=LengthPath(v,v1);//EstimateLengthByParam<FaceType>(v0,v1,edgeF);
209 Lengths[i]=length;
210 aveL+=length;
211 }
212 aveL/=(ScalarType)star.size();
213
214 for (unsigned int i=0;i<faces.size();i++)
215 {
216 Areas[i]=EstimateAreaByParam<FaceType>(faces[i]);
217 aveA+=Areas[i];
218 }
219 aveA/=(ScalarType)faces.size();
220
221 for (unsigned int i=0;i<Lengths.size();i++)
222 priorL+=pow((Lengths[i]-aveL),(ScalarType)2);
223
224 for (unsigned int i=0;i<Areas.size();i++)
225 priorA+=pow((Areas[i]-aveA),(ScalarType)2);
226
227 return (pow(priorL,(ScalarType)2)/2.0+priorA);
228 }
229
FindVarianceLenghtArea(MeshType & base_mesh,const ScalarType & averageLength,const ScalarType & averageArea,ScalarType & varianceL,ScalarType & varianceA)230 static void FindVarianceLenghtArea(MeshType &base_mesh,
231 const ScalarType &averageLength,
232 const ScalarType &averageArea,
233 ScalarType &varianceL,
234 ScalarType &varianceA)
235 {
236 typename MeshType::FaceIterator Fi;
237 varianceA=0;
238 varianceL=0;
239 int num_edge=0;
240 for (Fi=base_mesh.face.begin();Fi!=base_mesh.face.end();Fi++)
241 {
242 ScalarType area=EstimateAreaByParam<FaceType>(&(*Fi));
243 varianceA+=pow((area-averageArea),(ScalarType)2.0);
244 for (int i=0;i<3;i++)
245 {
246 VertexType *v0=(*Fi).V(i);
247 VertexType *v1=(*Fi).V((i+1)%3);
248 /* std::vector<FaceType*> on_edge,faces1,faces2;
249 getSharedFace<MeshType>(v0,v1,on_edge,faces1,faces2);
250 FaceType* edgeF[2];
251 edgeF[0]=on_edge[0];
252 edgeF[1]=on_edge[1];*/
253 if (v0>v1)
254 {
255 ScalarType length=LengthPath(v0,v1);//EstimateLengthByParam<FaceType>(v0,v1,edgeF);
256 varianceL+=pow((length-averageLength),(ScalarType)2);
257 num_edge++;
258 }
259 }
260 }
261 varianceL=sqrt(varianceL/(ScalarType)num_edge);
262 varianceA=sqrt(varianceA/(ScalarType)base_mesh.fn);
263 }
264
265 ///optimize UV of central vertex
OptimizeUV(VertexType * center,MeshType & base_domain)266 static void OptimizeUV(VertexType *center,MeshType &base_domain)
267 {
268 ///parametrize base domain star and subvertices
269 ParametrizeStarEquilateral<MeshType>(center,true);
270
271 ///get incident faces
272 std::vector<FaceType*> faces;
273 std::vector<VertexType*> vertices;
274 vertices.push_back(center);
275 getSharedFace<MeshType>(vertices,faces);
276 MeshType domain;
277
278 ///get Hres Vertices
279 std::vector<VertexType*> Hres_vert;
280 getHresVertex<typename MeshType::FaceType>(faces,Hres_vert);
281
282 ///make a copy of base mesh
283 std::vector<FaceType*> ordFaces;
284 CreateMeshVertexStar<MeshType>(vertices,ordFaces,domain);
285 assert(ordFaces.size()==domain.face.size());
286 assert(ordFaces.size()==faces.size());
287 /* assert(Test(ordFaces,faces));
288 assert(Test1(ordFaces,domain.face));
289 assert(Test1(faces,domain.face));*/
290 /*assert(Test2(domain,Hres_vert.size()));*/
291
292 UpdateTopologies<MeshType>(&domain);
293
294 ///minimization
295 minInfoUV Minf;
296 ///setting parameters for minimization
297 //Minf.base_domain=base_domain;
298 Minf.parametrized_domain=&domain;
299 //Minf.base_domain=&base_mesh;
300 Minf.Hres_vert=std::vector<VertexType*>(Hres_vert.begin(),Hres_vert.end());
301
302 ///create a copy of hres mesh
303 std::vector<VertexType*> OrderedVertices;
304 std::vector<FaceType*> OrderedFaces;
305 CopyMeshFromVertices<MeshType>(Hres_vert,OrderedVertices,OrderedFaces,Minf.hres_mesh);
306
307 ///get the vertex to optimize position
308 int i=0;
309 while (domain.vert[i].IsB()) i++;
310 Minf.to_optimize=&domain.vert[i];
311
312 ///texture value of central vertex
313 ///that should be optimized
314 float *p=new float [2];
315 p[0]=0;
316 p[1]=0;
317
318 ///allocate vector of output
319 float *x=new float [2];
320 x[0]=0;
321 x[1]=0;
322
323 float opts[LM_OPTS_SZ], info[LM_INFO_SZ];
324 opts[0]=(float)LM_INIT_MU;
325 opts[1]=(float)1E-15;
326 opts[2]=(float)1E-15;
327 opts[3]=(float)1E-20;
328 opts[4]=(float)LM_DIFF_DELTA;
329
330 /*int num=*/slevmar_dif(Equi_energy,p,x,2,2,1000,opts,info,NULL,NULL,&Minf);
331
332
333 ///copy back values
334
335 //clear old values
336 for (unsigned int i=0;i<ordFaces.size();i++)
337 ordFaces[i]->vertices_bary.resize(0);
338
339
340 //reassing
341 int num=0;
342 for (unsigned int i=0;i<domain.face.size();i++)
343 {
344 for (unsigned int j=0;j<domain.face[i].vertices_bary.size();j++)
345 {
346 VertexType *vert=domain.face[i].vertices_bary[j].first;
347 CoordType bary=domain.face[i].vertices_bary[j].second;
348 ordFaces[i]->vertices_bary.push_back(std::pair<VertexType*,vcg::Point3f>(vert,bary));
349 /*vert->father=ordFaces[i];
350 assert(!ordFaces[i]->IsD());
351 vert->Bary=bary;*/
352 AssingFather(*vert,ordFaces[i],bary,base_domain);
353 num++;
354 }
355 }
356 //assert(num==Minf.Hres_vert.size());
357 if (size_t(num)!=Minf.Hres_vert.size())
358 {
359 printf("num0 %d \n",num);
360 printf("num1 %d \n",int(Minf.Hres_vert.size()));
361 }
362
363 center->RPos=Minf.to_optimize->RPos;
364 delete[] x;
365 delete[] p;
366 }
367
368 struct Elem
369 {
370 public:
371 VertexType *center;
372 ScalarType priority;
373 int t_mark;
374
375 //bool operator <(Elem* e){return (priority<e.priority);}
376 inline bool operator <(const Elem& e) const
377 {
378 return (priority<e.priority);
379 //return (locModPtr->Priority() < h.locModPtr->Priority());
380 }
381
ElemElem382 Elem(VertexType *_center,ScalarType _priority,int _t_mark)
383 {
384 center=_center;
385 priority=_priority;
386 t_mark=_t_mark;
387 }
388 };
389
390 std::vector<Elem> Operations;
391
392
Execute(VertexType * center)393 void Execute(VertexType *center)
394 {
395 OptimizeUV(center,base_mesh);
396 std::vector<typename MeshType::VertexType*> neigh;
397 getVertexStar<MeshType>(center,neigh);
398
399 //push back neighbors and update mark
400 global_mark++;
401 ///update mark neighbors
402 for(unsigned int i=0;i<neigh.size();i++)
403 markers[neigh[i]]=global_mark;
404
405 ///push_back neighbors to the heap
406 for(unsigned int i=0;i<neigh.size();i++)
407 {
408 Operations.push_back(Elem(neigh[i],Priority(neigh[i]),global_mark));
409 std::push_heap( Operations.begin(), Operations.end());
410 }
411 }
412
413
414
OptimizePatches()415 void OptimizePatches()
416 {
417 global_mark=0;
418 markers.Init(global_mark);
419
420 Operations.clear();
421
422 const ScalarType sqrtsrt3=(ScalarType)1.31607401;
423
424 averageArea=Area(final_mesh)/((ScalarType)base_mesh.fn*(ScalarType)2.0);
425 averageLength=(ScalarType)2.0*sqrt(averageArea)/sqrtsrt3;
426
427 ScalarType varianceL,varianceA;
428 FindVarianceLenghtArea(base_mesh,averageLength,averageArea,varianceL,varianceA);
429
430 #ifndef _MESHLAB
431 printf("Variance length:%f\n",varianceL*100.f/averageLength);
432 printf("Variance area:%f\n",varianceA*100.f/averageArea);
433 #endif
434 //Initialize heap
435 for (unsigned int i=0;i<base_mesh.vert.size();i++)
436 {
437 VertexType *v=&base_mesh.vert[i];
438 Operations.push_back(Elem(v,Priority(v),global_mark));
439 }
440 std::make_heap(Operations.begin(),Operations.end());
441 float gap=(ScalarType)0.05;
442 int n_oper=0;
443 ///start Optimization
444 ScalarType varianceL0=varianceL;
445 ScalarType varianceA0=varianceA;
446 ScalarType varianceL1;
447 ScalarType varianceA1;
448 bool continue_opt=true;
449 while (continue_opt)
450 {
451 int temp_oper=0;
452 while (temp_oper<20)
453 {
454 std::pop_heap(Operations.begin(),Operations.end());
455 VertexType* oper=Operations.back().center;
456 int t_mark=Operations.back().t_mark;
457 Operations.pop_back();
458 if (markers[oper]<=t_mark)
459 {
460 Execute(oper);
461 n_oper++;
462 temp_oper++;
463 }
464 }
465 FindVarianceLenghtArea(base_mesh,averageLength,averageArea,varianceL1,varianceA1);
466 ScalarType percL=(varianceL0-varianceL1)*100/averageLength;
467 ScalarType percA=(varianceA0-varianceA1)*100/averageArea;
468 ScalarType curr_gap=percL+percA;
469 #ifndef _MESHLAB
470 printf("gap:%f\n",curr_gap);
471 #endif
472 varianceL0=varianceL1;
473 varianceA0=varianceA1;
474 continue_opt=curr_gap>gap;
475 }
476 FindVarianceLenghtArea(base_mesh,averageLength,averageArea,varianceL,varianceA);
477
478 #ifndef _MESHLAB
479 printf("Num Oper:%i\n",n_oper);
480 printf("Variance length:%f\n",varianceL*100.f/averageLength);
481 printf("Variance area:%f\n",varianceA*100.f/averageArea);
482 #endif
483
484 }
485
PatchesOptimizer(MeshType & _base_mesh,MeshType & _final_mesh)486 PatchesOptimizer(MeshType &_base_mesh,MeshType &_final_mesh):base_mesh(_base_mesh),final_mesh(_final_mesh),markers(_base_mesh.vert){}
487
HresMesh()488 static MeshType* &HresMesh()
489 {
490 static MeshType* mesh;
491 return mesh;
492 }
493
BaseMesh()494 static MeshType* &BaseMesh()
495 {
496 static MeshType* mesh;
497 return mesh;
498 }
499
500 };
501 #endif
502