1 #include "modButterfly.h"
2 #include "vcg/simplex/face/component.h"
3 #include <math.h>
4
5 /****************************************************************************
6 * Rgb Triangulations Plugin *
7 * *
8 * Author: Daniele Panozzo (daniele.panozzo@gmail.com) *
9 * Copyright(C) 2007 *
10 * DISI - Department of Computer Science *
11 * University of Genova *
12 * *
13 * All rights reserved. *
14 * *
15 * This program is free software; you can redistribute it and/or modify *
16 * it under the terms of the GNU General Public License as published by *
17 * the Free Software Foundation; either version 2 of the License, or *
18 * (at your option) any later version. *
19 * *
20 * This program is distributed in the hope that it will be useful, *
21 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
23 * GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
24 * for more details. *
25 ****************************************************************************/
26
27 namespace rgbt
28 {
29
30
init(TriMeshType & m,RgbInfo & info)31 void ModButterfly::init(TriMeshType& m, RgbInfo& info)
32 {
33 for (unsigned int i = 0; i < m.vert.size(); ++i)
34 {
35 //std::cerr << "index:" << i << std::endl;
36 VertexPointer v = & m.vert[i];
37 if (!v->IsD() && v->VFp())
38 {
39
40 RgbTriangleC t = RgbTriangleC(m,info,v->VFp()->Index());
41 int ti = v->VFi();
42 RgbVertexC vr = t.V(ti);
43 assert(&vr.vert() == v);
44
45 vector<RgbVertexC> vv;
46 vv.reserve(10);
47 assert(vv.size() == 0);
48 RgbPrimitives::VV(vr,vv);
49 assert(vv.size() >= 2);
50 vr.setBaseArity(vv.size());
51
52 }
53 }
54
55 }
56
findHalfStencil(RgbVertexC & v,Pos & pos,vector<RgbVertexC> & stencil)57 void ModButterfly::findHalfStencil(RgbVertexC& v, Pos& pos, vector<RgbVertexC>& stencil)
58 {
59 assert(v.vp() == pos.v);
60
61 Pos p;
62 p = pos;
63 p.FlipV();
64 RgbVertexC v2 = RgbVertexC(*v.m,*v.rgbInfo,p.v);
65
66 int maxlevel = v.getLevel();
67 if (maxlevel < v2.getLevel())
68 maxlevel = v2.getLevel();
69
70 p = pos;
71 rotate(v,p,2);
72 RgbVertexC vtemp = move(v,p,maxlevel);
73 stencil.push_back(vtemp);
74
75 #ifndef NDEBUG
76 Pos ptemp;
77 ptemp = pos;
78 ptemp.FlipV();
79 rotate(v2,ptemp,2);
80 RgbVertexC vtemp2 = move(v2,ptemp,maxlevel);
81 assert(vtemp.index == vtemp2.index);
82 #endif
83
84 p = pos;
85 rotate(v,p,4);
86 stencil.push_back(move(v,p,maxlevel));
87
88 p = pos;
89 p.FlipV();
90 rotate(v2,p,4);
91 stencil.push_back(move(v2,p,maxlevel));
92 }
93
doSplit(RgbTriangleC & fp,int EdgeIndex,int level,TopologicalOpC & to,vector<FacePointer> * vfp)94 bool ModButterfly::doSplit(RgbTriangleC& fp, int EdgeIndex, int level, TopologicalOpC& to , vector<FacePointer> *vfp)
95 {
96 RgbVertexC v1ext;
97 RgbVertexC v2ext;
98
99 assert(EdgeIndex >= 0 && EdgeIndex <= 2);
100
101 RgbVertexC v1 = fp.V(EdgeIndex);
102 RgbVertexC v2 = fp.V((EdgeIndex+1)%3);
103 bool isBorder = fp.getEdgeIsBorder(EdgeIndex);
104
105 vector<RgbVertexC> stencil1;
106 stencil1.reserve(3);
107 vector<RgbVertexC> stencil2;
108 stencil2.reserve(3);
109
110 Point p;
111
112 if (!isBorder)
113 {
114 bool v1r = (baseArity(v1) != 6) || (v1.getIsBorder());
115 bool v2r = (baseArity(v2) != 6) || (v2.getIsBorder());
116
117 int maxlevel = v1.getLevel();
118 if (maxlevel < v2.getLevel())
119 maxlevel = v2.getLevel();
120
121 if (v1r && v2r) // both the vertices are extraordinary
122 {
123 RgbPrimitives::splitGreenEdgeIfNeeded(v1,maxlevel + 1,to);
124 RgbPrimitives::splitRedEdgeIfNeeded(v1,maxlevel + 1,to);
125 RgbPrimitives::splitGreenEdgeIfNeeded(v2,maxlevel + 1,to);
126 RgbPrimitives::splitRedEdgeIfNeeded(v2,maxlevel + 1,to);
127
128 assert(RgbPrimitives::IsValidEdge(v1,v2,&fp,&EdgeIndex));
129 RgbPrimitives::IsValidEdge(v1,v2,&fp,&EdgeIndex);
130
131 Pos pos1 = Pos(fp.face(),EdgeIndex);
132 Pos pos2 = pos1;
133 pos2.FlipF();
134 pos2.FlipV();
135
136 p = computeExtraordinary(v1,pos1);
137 p += computeExtraordinary(v2,pos2);
138 p /= 2.0;
139 }
140 else if (v1r)
141 {
142 RgbPrimitives::splitGreenEdgeIfNeeded(v1,maxlevel + 1,to);
143 RgbPrimitives::splitRedEdgeIfNeeded(v1,maxlevel + 1,to);
144
145 assert(RgbPrimitives::IsValidEdge(v1,v2,&fp,&EdgeIndex));
146 RgbPrimitives::IsValidEdge(v1,v2,&fp,&EdgeIndex);
147
148 Pos pos1 = Pos(fp.face(),EdgeIndex);
149 p = computeExtraordinary(v1,pos1);
150 }
151 else if (v2r)
152 {
153 RgbPrimitives::splitGreenEdgeIfNeeded(v2,maxlevel + 1,to);
154 RgbPrimitives::splitRedEdgeIfNeeded(v2,maxlevel + 1,to);
155
156 assert(RgbPrimitives::IsValidEdge(v1,v2,&fp,&EdgeIndex));
157 RgbPrimitives::IsValidEdge(v1,v2,&fp,&EdgeIndex);
158
159 Pos pos2 = Pos(fp.face(),EdgeIndex);
160
161 pos2.FlipV();
162 p = computeExtraordinary(v2,pos2);
163 }
164 else
165 {
166 // Additional split
167
168 RgbPrimitives::splitGreenEdgeIfNeeded(v1,maxlevel+1,to);
169 RgbPrimitives::splitGreenEdgeIfNeeded(v2,maxlevel+1,to);
170 RgbPrimitives::splitRedEdgeIfNeeded(v1,maxlevel+1,to);
171 RgbPrimitives::splitRedEdgeIfNeeded(v2,maxlevel+1,to);
172 RgbPrimitives::splitGreenEdgeIfNeeded(v1,maxlevel+1,to);
173 RgbPrimitives::splitGreenEdgeIfNeeded(v2,maxlevel+1,to);
174 RgbPrimitives::splitRedEdgeIfNeeded(v1,maxlevel+1,to);
175 RgbPrimitives::splitRedEdgeIfNeeded(v2,maxlevel+1,to);
176
177
178 if (!RgbPrimitives::IsValidEdge(v1,v2,&fp,&EdgeIndex))
179 {
180 return false; // split already done
181 }
182 // --------
183
184 Pos pos1 = Pos(fp.face(),EdgeIndex);
185 Pos pos2 = pos1;
186 pos2.FlipF();
187 pos2.FlipV();
188
189 findHalfStencil(v1,pos1,stencil1);
190 findHalfStencil(v2,pos2,stencil2);
191
192 assert(stencil1.size() == 3);
193 assert(stencil2.size() == 3);
194
195 Point acc;
196 Point pt;
197
198 pt = v1.getCoord();
199 pt *= 1.0/2.0;
200 acc = pt;
201
202 pt = v2.getCoord();
203 pt *= 1.0/2.0;
204 acc += pt;
205
206 pt = stencil1[0].getCoord();
207 pt *= 1.0/8.0;
208 acc += pt;
209
210 pt = stencil2[0].getCoord();
211 pt *= 1.0/8.0;
212 acc += pt;
213
214 pt = stencil1[1].getCoord();
215 pt *= 1.0/16.0;
216 acc -= pt;
217
218 pt = stencil1[2].getCoord();
219 pt *= 1.0/16.0;
220 acc -= pt;
221
222 pt = stencil2[1].getCoord();
223 pt *= 1.0/16.0;
224 acc -= pt;
225
226 pt = stencil2[2].getCoord();
227 pt *= 1.0/16.0;
228 acc -= pt;
229
230 p = acc;
231 }
232 }
233 else
234 {
235 // Boundary
236 // Additional split
237 int maxlevel = v1.getLevel();
238 if (maxlevel < v2.getLevel())
239 maxlevel = v2.getLevel();
240
241 RgbPrimitives::splitGreenEdgeIfNeeded(v1,maxlevel+1,to);
242 RgbPrimitives::splitGreenEdgeIfNeeded(v2,maxlevel+1,to);
243
244 if (!RgbPrimitives::IsValidEdge(v1,v2,&fp,&EdgeIndex))
245 {
246 return false; // split already done
247 }
248 // --------
249
250 Pos pos1 = Pos(fp.face(),EdgeIndex);
251 Pos pos2 = pos1;
252 pos2.FlipV();
253
254 rotateUntilBorder(v1,pos1);
255 rotateUntilBorder(v2,pos2);
256
257 v1ext = move(v1,pos1,maxlevel);
258 v2ext = move(v2,pos2,maxlevel);
259
260 assert(v1ext.getIsBorder());
261 assert(v2ext.getIsBorder());
262
263 Point acc;
264 Point pt;
265
266 pt = v1.getCoord();
267 pt *= 9.0/16.0;
268 acc = pt;
269
270 pt = v2.getCoord();
271 pt *= 9.0/16.0;
272 acc += pt;
273
274 pt = v1ext.getCoord();
275 pt *= 1.0/16.0;
276 acc -= pt;
277
278 pt = v2ext.getCoord();
279 pt *= 1.0/16.0;
280 acc -= pt;
281
282 p = acc;
283 }
284
285 VertexPointer vp;
286 vector<VertexPointer> vtemp;
287
288 if (!isBorder)
289 to.doSplit(fp.face(),EdgeIndex, p, vfp, &vtemp);
290 else
291 to.doSplitBoundary(fp.face(),EdgeIndex, p, vfp, &vtemp);
292
293 vp = vtemp[0];
294
295 RgbVertexC vNew = RgbVertexC(*(fp.m),*(fp.rgbInfo),vp);
296 vNew.resetInfo();
297 vNew.setLevel(level);
298 vNew.setCoord(p);
299 vNew.setIsNew(true);
300 vNew.setIsBorder(isBorder);
301
302 vector<RgbVertexC> vv;
303 RgbPrimitives::VV(vNew,vv,false);
304
305 RgbPrimitives::updateNormal(vNew);
306 for(vector<RgbVertexC>::iterator i = vv.begin(); i != vv.end(); i++)
307 RgbPrimitives::updateNormal(*i);
308
309 return true;
310 }
311
doCollapse(RgbTriangleC & fp,int EdgeIndex,TopologicalOpC & to,Point3<ScalarType> * p,vector<FacePointer> * vfp)312 void ModButterfly::doCollapse(RgbTriangleC& fp, int EdgeIndex, TopologicalOpC& to, Point3<ScalarType> *p, vector<FacePointer> *vfp)
313 {
314 //assert(!fp.V(EdgeIndex).getIsNew());
315 if (!fp.getEdgeIsBorder(EdgeIndex))
316 to.doCollapse(fp.face(),EdgeIndex, p, vfp);
317 else
318 to.doCollapseBoundary(fp.face(),EdgeIndex, p, vfp);
319
320 #ifndef NDEBUG
321 std::cerr << fp.m->vn << std::endl;
322 #endif
323
324 }
325
rotate(RgbVertexC & v,Pos & pos,int angle)326 void ModButterfly::rotate(RgbVertexC& v, Pos& pos,int angle)
327 {
328 //std::cerr << v.index << " requested: " << angle << " ----- ";
329 assert(v.vp() == pos.v);
330 while (angle > 0)
331 {
332 RgbTriangleC t = RgbTriangleC(v.m,v.rgbInfo,pos.F()->Index());
333 int realz = -1;
334 for (int i = 0; i<3; i++)
335 if (t.V(i).vp() == pos.v)
336 realz = i;
337 assert(realz != -1);
338
339 angle -= t.getAngle(realz);
340 //std::cerr << "| " << t.getAngle(realz);
341 pos.FlipE();
342 pos.FlipF();
343 }
344 //std::cerr << " ----- result : " << angle << std::endl;
345 //assert(angle == 0);
346 }
347
move(RgbVertexC & v,Pos & pos,int level)348 ModButterfly::RgbVertexC ModButterfly::move(RgbVertexC& v, Pos& pos, int level)
349 {
350 assert(v.vp() == pos.v);
351 assert(level >= v.getLevel());
352 int vl = level;
353 RgbVertexC lastv;
354 do
355 {
356 pos.FlipV();
357 lastv = RgbVertexC(*v.m,*v.rgbInfo,pos.v);
358 if (lastv.getLevel() <= vl)
359 break; // no need to continue, in fact a rotation of exactly 6 can be impossible in some cases
360 rotate(lastv,pos,6);
361 pos.FlipF();
362
363 assert(lastv.getLevel() > vl);
364 } while (lastv.getLevel() > vl);
365
366 return RgbVertexC(*v.m,*v.rgbInfo,pos.v);
367 }
368
rotateUntilBorder(RgbVertexC & v,Pos & pos)369 void ModButterfly::rotateUntilBorder(RgbVertexC& v, Pos& pos)
370 {
371 assert(v.vp() == pos.v);
372 do
373 {
374 pos.FlipE();
375 pos.FlipF();
376 } while (!pos.IsBorder());
377 assert(v.vp() == pos.v);
378 }
379
computeExtraordinary(RgbVertexC & v,Pos & pos)380 ModButterfly::Point ModButterfly::computeExtraordinary(RgbVertexC& v, Pos& pos)
381 {
382 int rank;
383 if (v.getLevel() == 0)
384 rank = baseArity(v);
385 else
386 rank = 4;
387
388 vector<Point> points;
389 vector<double> coeff;
390 computeExtraordinaryPattern(coeff,rank);
391 points.reserve(rank+1);
392 points.push_back(v.getCoord());
393
394 if (!v.getIsBorder())
395 {
396 int angle = rank * 2;
397 do
398 {
399 pos.FlipV();
400 points.push_back(pos.v->P());
401 assert(pos.v != v.vp());
402 pos.FlipV();
403 rotate(v,pos,2);
404 angle -= 2;
405 assert (angle >= 0);
406 }
407 while (angle != 0);
408 }
409 else
410 {
411 Pos pos2 = pos;
412 // Border
413 int angle = (rank-1) * 2;
414 do
415 {
416 pos.FlipV();
417 points.push_back(pos.v->P());
418 assert(pos.v != v.vp());
419 pos.FlipV();
420 assert(pos.v == v.vp());
421 rotate(v,pos,2);
422 angle -= 2;
423 assert (angle >= 0);
424 }
425 while ((!pos.IsBorder()) && (angle != 0));
426
427 pos.FlipV();
428 points.push_back(pos.v->P());
429 assert(pos.v != v.vp());
430 pos.FlipV();
431
432 assert(angle >= 0);
433 vector<Point> pointsreverse;
434
435 pos2.FlipF();
436 rotate(v,pos2,2);
437 angle -= 2;
438 if ((!pos2.IsBorder()) && (angle != 0))
439 {
440 do
441 {
442 pos2.FlipV();
443 pointsreverse.push_back(pos2.v->P());
444 assert(pos2.v != v.vp());
445 pos2.FlipV();
446 assert(pos2.v == v.vp());
447 rotate(v,pos2,2);
448 angle -= 2;
449 assert (angle >= 0);
450 }
451 while ((!pos2.IsBorder()) && (angle != 0));
452 }
453 assert(angle == 0);
454
455 pos.FlipV();
456 points.push_back(pos.v->P());
457 assert(pos.v != v.vp());
458 pos.FlipV();
459
460 for(int i=pointsreverse.size()-1;i>= 0; i--)
461 points.push_back(pointsreverse[i]);
462
463 }
464
465 assert(points.size() == (unsigned int)(rank + 1));
466 assert(points.size() == coeff.size());
467
468 Point acc = Point(0.0,0.0,0.0);
469
470 for (int i=0; i<rank+1; i++)
471 {
472 Point tmp = points[i] * coeff[i];
473 acc += tmp;
474 }
475
476 return acc;
477 }
478
computeExtraordinaryPattern(vector<double> & pattern,int k)479 void ModButterfly::computeExtraordinaryPattern(vector<double>& pattern, int k)
480 {
481 assert(k >= 3);
482 pattern.clear();
483 pattern.reserve(k+1);
484 pattern.push_back(3.0/4.0);
485
486 switch (k) {
487 case 3:
488 pattern.push_back(5.0/12.0);
489 pattern.push_back(-(1.0/12.0));
490 pattern.push_back(-(1.0/12.0));
491 break;
492 case 4:
493 pattern.push_back(3.0/8.0);
494 pattern.push_back(0);
495 pattern.push_back(-(1.0/8.0));
496 pattern.push_back(0);
497 break;
498 default:
499 for (int i = 0; i < k; i++)
500 {
501 double kd = k;
502 double id = i;
503 double v = (1.0/kd)*((1.0/4.0) + cos((2.0*id*M_PI)/kd) + ((1.0/2.0)*cos((4.0*id*M_PI)/kd)));
504 pattern.push_back(v);
505 }
506 break;
507 }
508 }
509
baseArity(RgbVertexC & v)510 int ModButterfly::baseArity(RgbVertexC& v)
511 {
512 int arity;
513 if (v.getLevel() > 0)
514 {
515 arity = 6;
516 }
517 else
518 {
519 arity = v.getBaseArity();
520 }
521 return arity;
522 }
523
524
525 }
526