Line data Source code
1 : /*
2 : * Copyright (c) 2014-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Sebastian Reiter
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #include <fstream>
34 : #include <queue>
35 : #include "lib_grid/lg_base.h"
36 : #include "common/profiler/profiler.h"
37 : #include "simple_grid.h"
38 : #include "edge_length_adjustment_extended.h"
39 : #include "lib_grid/refinement/regular_refinement.h"
40 : #include "common/node_tree/node_tree.h"
41 : #include "lib_grid/algorithms/trees/octree.h"
42 : #include "lib_grid/callbacks/callbacks.h"
43 : #include "../trees/octree.h"
44 :
45 : using namespace std;
46 :
47 : namespace ug
48 : {
49 :
50 :
51 0 : AdjustEdgeLengthDesc::
52 0 : AdjustEdgeLengthDesc() :
53 0 : minEdgeLen(1),
54 0 : maxEdgeLen(2),
55 0 : approximation(0.9),
56 0 : triQuality(0.9),
57 : // adaptivity(0.5),
58 0 : projectPoints(true)
59 0 : {}
60 :
61 :
62 : /// only for debugging purposes!!!
63 : /** Output value pairs to gnuplot...
64 : * \{ */
65 : //#define EDGE_LENGTH_ADJUSTMENT__GPLOT_ENABLED
66 : #ifdef EDGE_LENGTH_ADJUSTMENT__GPLOT_ENABLED
67 : typedef vector<pair<number, number> > GnuplotData;
68 : static GnuplotData gplotLengthFac;
69 : static GnuplotData gplotMinCurvature;
70 : static GnuplotData gplotAverageCurvature;
71 :
72 : void WriteGnuplotData(const char* filename, const GnuplotData& data)
73 : {
74 : ofstream out(filename);
75 : if(!out)
76 : return;
77 :
78 : for(size_t i = 0; i < data.size(); ++i)
79 : out << data[i].first << " " << data[i].second << endl;
80 :
81 : out.close();
82 : }
83 :
84 : #define GPLOTPOINT(dataName, x, y) dataName.push_back(make_pair<number, number>((x), (y)));
85 : #define GPLOTSAVE() {WriteGnuplotData("length_fac.gplot", gplotLengthFac);\
86 : WriteGnuplotData("min_curvature.gplot", gplotMinCurvature);\
87 : WriteGnuplotData("average_curvature.gplot", gplotAverageCurvature);}
88 : #else
89 : // do nothing if EDGE_LENGTH_ADJUSTMENT__GPLOT_ENABLED is false
90 : #define GPLOTPOINT(dataName, x, y)
91 : #define GPLOTSAVE()
92 : #endif
93 : /** \} */
94 :
95 :
96 : /*
97 : vector3 PNTrianglePos(const vector3& p0, const vector3& p1, const vector3& p2,
98 : const vector3& n0, const vector3& n1, const vector3& n2);
99 :
100 : vector3 PNTriangleNorm(const vector3& p0, const vector3& p1, const vector3& p2,
101 : const vector3& n0, const vector3& n1, const vector3& n2);
102 :
103 : vector3 PNCTrianglePos(const vector3& p0, const vector3& p1, const vector3& p2,
104 : const vector3& n0, const vector3& n1, const vector3& n2,
105 : const vector3& cn0, const vector3& cn1, const vector3& cn2);
106 :
107 : vector3 PNCTriangleNorm(const vector3& p0, const vector3& p1, const vector3& p2,
108 : const vector3& n0, const vector3& n1, const vector3& n2,
109 : const vector3& cn0, const vector3& cn1, const vector3& cn2);
110 :
111 :
112 : vector3 pos(number bc0, number bc1);
113 : vector3 norm(number bc0, number bc1);
114 : };
115 :
116 : class ProjectedVertex{
117 : Vertex* vertex();
118 : vector3 vertex_position();
119 : vector3 vertex_normal();
120 : vector3 surface_normal();
121 : vector3 surface_position();
122 : size_t num_barycentric_coords();
123 : number barycentric_coord(size_t index);
124 : };
125 :
126 : class SurfaceRepresentation{
127 : public:
128 :
129 : };
130 : */
131 :
132 : ////////////////////////////////////////////////////////////////////////
133 0 : static void AssignFixedVertices(Grid& grid, SubsetHandler& shMarks)
134 : {
135 0 : grid.begin_marking();
136 :
137 : // mark all vertices that are not regular crease-vertices as fixed
138 : for(EdgeIterator iter = shMarks.begin<Edge>(REM_CREASE);
139 0 : iter != shMarks.end<Edge>(REM_CREASE); ++iter)
140 : {
141 : Edge* e = *iter;
142 0 : for(size_t i = 0; i < 2; ++i){
143 0 : Vertex* vrt = e->vertex(i);
144 0 : if(!grid.is_marked(vrt)){
145 : grid.mark(vrt);
146 : int counter = 0;
147 0 : for(Grid::AssociatedEdgeIterator nbIter = grid.associated_edges_begin(vrt);
148 0 : nbIter != grid.associated_edges_end(vrt); ++nbIter)
149 : {
150 0 : if(shMarks.get_subset_index(*nbIter) != REM_NONE)
151 0 : ++counter;
152 : }
153 :
154 0 : if(counter != 2)
155 0 : shMarks.assign_subset(vrt, REM_FIXED);
156 : }
157 : }
158 : }
159 :
160 0 : grid.end_marking();
161 :
162 : // mark all vertices that lie on a fixed edge as fixed vertex
163 : for(EdgeIterator iter = shMarks.begin<Edge>(REM_FIXED);
164 0 : iter != shMarks.end<Edge>(REM_FIXED); ++iter)
165 : {
166 : Edge* e = *iter;
167 0 : shMarks.assign_subset(e->vertex(0), REM_FIXED);
168 0 : shMarks.assign_subset(e->vertex(1), REM_FIXED);
169 : }
170 0 : }
171 :
172 0 : static void AssignCreaseVertices(Grid& grid, SubsetHandler& shMarks)
173 : {
174 : // mark all vertices that lie on a crease and which are not fixed
175 : // as crease vertices.
176 0 : if(shMarks.num_subsets() <= REM_CREASE)
177 : return;
178 :
179 : for(EdgeIterator iter = shMarks.begin<Edge>(REM_CREASE);
180 0 : iter != shMarks.end<Edge>(REM_CREASE); ++iter)
181 : {
182 : Edge* e = *iter;
183 0 : for(uint i = 0; i < 2; ++i)
184 0 : if(shMarks.get_subset_index(e->vertex(i)) != REM_FIXED)
185 0 : shMarks.assign_subset(e->vertex(i), REM_CREASE);
186 : }
187 : }
188 :
189 :
190 : ////////////////////////////////////////////////////////////////////////
191 : template <class TVertexPositionAccessor>
192 : static number CalculateNormalDot(TriangleDescriptor& td1, TriangleDescriptor& td2,
193 : TVertexPositionAccessor& aaPos)
194 : {
195 : vector3 n1;
196 : CalculateTriangleNormal(n1, aaPos[td1.vertex(0)],
197 : aaPos[td1.vertex(1)], aaPos[td1.vertex(2)]);
198 : vector3 n2;
199 : CalculateTriangleNormal(n2, aaPos[td2.vertex(0)],
200 : aaPos[td2.vertex(1)], aaPos[td2.vertex(2)]);
201 : return VecDot(n1, n2);
202 : }
203 :
204 :
205 : // ////////////////////////////////////////////////////////////////////////
206 : // // CalculateCurvature
207 : // template <class TAAPosVRT>
208 : // number CalculateMinCurvature(Grid& grid, SubsetHandler& shMarks,
209 : // Vertex* vrt, TAAPosVRT& aaPos)
210 : // {
211 : // //TODO: check whether static vNormals brings any benefits.
212 : // //TODO: special cases for crease vertices
213 :
214 : // // face normals
215 : // static vector<vector3> vNormals;
216 : // // vertex normal (mean face normal)
217 : // vector3 n(0, 0, 0);
218 :
219 : // // calculate the normals of associated faces
220 : // vNormals.clear();
221 : // Grid::AssociatedFaceIterator iterEnd = grid.associated_faces_end(vrt);
222 : // for(Grid::AssociatedFaceIterator iter = grid.associated_faces_begin(vrt);
223 : // iter != iterEnd; ++iter)
224 : // {
225 : // vector3 nTmp;
226 : // CalculateNormal(nTmp, *iter, aaPos);
227 : // vNormals.push_back(nTmp);
228 : // VecAdd(n, n, nTmp);
229 : // }
230 :
231 : // // the vertex normal
232 : // VecNormalize(n, n);
233 :
234 : // // get the min dot-product of the vertex normal with associated faces-normals
235 : // number minDot = 1;
236 : // for(size_t i = 0; i < vNormals.size(); ++i)
237 : // minDot = std::min(minDot, VecDot(n, vNormals[i]));
238 :
239 : // //todo: Think about converting to radiants.
240 : // //minDot = acos(minDot) / (0.5 * PI);
241 : // //...
242 : // // done
243 : // GPLOTPOINT(gplotMinCurvature, 0, minDot);
244 : // return minDot;
245 : // }
246 :
247 : // ////////////////////////////////////////////////////////////////////////
248 : // template <class TAAPosVRT>
249 : // number CalculateAverageCurvature(Grid& grid, SubsetHandler& shMarks,
250 : // Edge* e, TAAPosVRT& aaPos)
251 : // {
252 : // number avCurv = 0.5 * (CalculateMinCurvature(grid, shMarks, e->vertex(0), aaPos)
253 : // + CalculateMinCurvature(grid, shMarks, e->vertex(1), aaPos));
254 : // GPLOTPOINT(gplotAverageCurvature, 0.5, avCurv);
255 : // return avCurv;
256 : // }
257 :
258 : // ////////////////////////////////////////////////////////////////////////
259 : // template <class TAAPosVRT>
260 : // number CalculateLengthFac(Grid& grid, SubsetHandler& shMarks,
261 : // Edge* e, TAAPosVRT& aaPos)
262 : // {
263 : // number lenFac = CalculateAverageCurvature(grid, shMarks, e, aaPos);
264 : // lenFac = (lenFac - 0.95) / 0.05;
265 : // lenFac = max(number(0.25), lenFac);
266 : // GPLOTPOINT(gplotLengthFac, 0.5, lenFac);
267 : // return lenFac;
268 : // }
269 :
270 : // ////////////////////////////////////////////////////////////////////////
271 : // ////////////////////////////////////////////////////////////////////////
272 : // template <class TAAPosVRT, class TAANormVRT, class TAAIntVRT>
273 : // bool TrySwap(const AdjustEdgeLengthDesc& desc, Grid& grid, Edge* e,
274 : // TAAPosVRT& aaPos, TAANormVRT& aaNorm,
275 : // TAAIntVRT& aaInt, SubsetHandler* pshMarks = NULL,
276 : // EdgeSelector* pCandidates = NULL)
277 : // {
278 : // // swaps are neither allowed for crease edges nor for fixed edges
279 : // if(pshMarks){
280 : // if(pshMarks->get_subset_index(e) == REM_FIXED ||
281 : // pshMarks->get_subset_index(e) == REM_CREASE)
282 : // return false;
283 : // }
284 :
285 : // // get the associated faces. we need two of them
286 : // Face* f[2];
287 : // if(GetAssociatedFaces(f, grid, e, 2) != 2)
288 : // return false;
289 :
290 : // // make sure both faces are triangles
291 : // if((f[0]->num_vertices() != 3) || (f[1]->num_vertices() != 3))
292 : // return false;
293 :
294 : // // create a simple grid
295 : // SimpleGrid sg;
296 : // if(!ObtainSimpleGrid(sg, grid, e->vertex(0), e->vertex(1), 0,
297 : // aaPos, aaNorm, aaInt))
298 : // {
299 : // LOG("ObtainSimpleGrid failed. ignoring edge...\n");
300 : // return false;
301 : // }
302 :
303 : // // calculate geometric-approximation-degree and triangle quality
304 : // number approxDeg = GeometricApproximationDegree(sg);
305 : // number shapeDeg = ShapeQualityDegree(sg);
306 : // // number smoothDeg = VecDot(sg.triangleNormals[0], sg.triangleNormals[1]);
307 :
308 : // // this is a new test. the idea is that each edge should be orthogonal to
309 : // // the normals of its endpoints - at least in a perfectly smooth surface.
310 :
311 : // // number approxDeg;
312 : // // number newApproxDeg;
313 : // // {
314 : // // vector3 v;
315 : // // VecSubtract(v, sg.vertices[1], sg.vertices[0]);
316 : // // VecNormalize(v, v);
317 : // // approxDeg = 1. - 0.5*(fabs(VecDot(v, sg.vertexNormals[0])) +
318 : // // fabs(VecDot(v, sg.vertexNormals[1])));
319 :
320 : // // VecSubtract(v, sg.vertices[3], sg.vertices[2]);
321 : // // VecNormalize(v, v);
322 : // // newApproxDeg = 1. - 0.5*(fabs(VecDot(v, sg.vertexNormals[2])) +
323 : // // fabs(VecDot(v, sg.vertexNormals[3])));
324 : // // }
325 :
326 : // // perform a swap on the simple grid
327 : // if(!SwapEdge(sg))
328 : // {
329 : // LOG("swap edge failed...\n");
330 : // return false;
331 : // }
332 :
333 : // // calculate new geometric-approximation-degree and triangle quality
334 : // number newApproxDeg = GeometricApproximationDegree(sg);
335 : // number newShapeDeg = ShapeQualityDegree(sg);
336 : // // number newSmoothDeg = VecDot(sg.triangleNormals[0], sg.triangleNormals[1]);
337 :
338 : // // // neither the shapeDeg nor the approxDeg may get too bad.
339 : // // if((newApproxDeg < 0.5 * approxDeg) || (newShapeDeg < 0.5 * shapeDeg))
340 : // // return false;
341 :
342 : // // // make sure that the swap does not destroy the smoothness
343 : // // if(newSmoothDeg < 0.1 * smoothDeg)
344 : // // return false;
345 :
346 : // // if(0.2 * (newApproxDeg - approxDeg) + 0.8 * (newShapeDeg - shapeDeg) > 0)
347 : // //if(newShapeDeg > shapeDeg)
348 : // //if(newApproxDeg > approxDeg)
349 :
350 :
351 : // if((newApproxDeg > desc.approximation * approxDeg) ||
352 : // (newApproxDeg >= desc.approximation && newShapeDeg > shapeDeg))
353 : // {
354 : // // swap the edge
355 : // e = SwapEdge(grid, e);
356 :
357 : // if(e){
358 : // // swap was successful
359 : // // if pCandidates was specified then add new candidates
360 : // if(pCandidates){
361 : // for(int i = 0; i < 2; ++i)
362 : // pCandidates->select(grid.associated_edges_begin(e->vertex(i)),
363 : // grid.associated_edges_end(e->vertex(i)));
364 : // // e was selected but is not really a candidate
365 : // pCandidates->deselect(e);
366 : // }
367 :
368 : // return true;
369 : // }
370 : // }
371 :
372 : // return false;
373 : // }
374 :
375 : // ////////////////////////////////////////////////////////////////////////
376 : // template <class TAAPosVRT, class TAANormVRT, class TAAIntVRT>
377 : // bool PerformSwaps(const AdjustEdgeLengthDesc& desc, Grid& grid,
378 : // SubsetHandler& shMarks, EdgeSelector& esel,
379 : // TAAPosVRT& aaPos, TAANormVRT& aaNorm, TAAIntVRT& aaInt)
380 : // {
381 : // PROFILE_FUNC();
382 : // LOG(" performing swaps\n");
383 : // int numSwaps = 0;
384 : // int maxNumSwaps = esel.num<Edge>() * 2;
385 :
386 : // while(!esel.empty())
387 : // {
388 : // Edge* e = *esel.begin<Edge>();
389 : // esel.deselect(e);
390 :
391 : // if(TrySwap(desc, grid, e, aaPos, aaNorm, aaInt, &shMarks, &esel)){
392 : // ++numSwaps;
393 : // if(numSwaps > maxNumSwaps){
394 : // UG_LOG(" aborting since maxNumSwaps was reached...");
395 : // esel.clear();
396 : // }
397 : // }
398 : // }
399 : // LOG(" swaps performed: " << numSwaps << endl);
400 :
401 : // return true;
402 : // }
403 :
404 : /** returns the resulting vertex or NULL, if no collapse was performed.*/
405 : template <class TAAPosVRT, class TAANormVRT, class TAAIntVRT>
406 0 : static Vertex* TryCollapse(const AdjustEdgeLengthDesc& desc, Grid& grid, Edge* e,
407 : TAAPosVRT& aaPos, TAANormVRT& aaNorm,
408 : TAAIntVRT& aaInt, SubsetHandler* pshMarks = NULL,
409 : EdgeSelector* pCandidates = NULL)
410 : {
411 0 : if(pshMarks)
412 : {
413 : SubsetHandler& shMarks = *pshMarks;
414 : // collapses are not allowed for fixed edges
415 0 : if(shMarks.get_subset_index(e) == REM_FIXED)
416 : return NULL;
417 :
418 : // if both endpoints of are fixed vertices then
419 : // we may not collapse
420 : int vrtSI[2];
421 0 : vrtSI[0] = shMarks.get_subset_index(e->vertex(0));
422 0 : vrtSI[1] = shMarks.get_subset_index(e->vertex(1));
423 :
424 0 : if((vrtSI[0] == REM_FIXED) && (vrtSI[1] == REM_FIXED))
425 : return NULL;
426 :
427 : // if both endpoints are somehow marked, e has to be a
428 : // crease edge
429 0 : if((vrtSI[0] != REM_NONE) && (vrtSI[1] != REM_NONE)
430 0 : && (shMarks.get_subset_index(e) != REM_CREASE))
431 : return NULL;
432 : }
433 :
434 : // check whether the edge can be collapsed
435 0 : if(EdgeCollapseIsValid(grid, e))
436 : {
437 : // test the collapse using a simple-grid
438 : SimpleGrid sgSrc;
439 0 : if(!ObtainSimpleGrid(sgSrc, grid, e->vertex(0), e->vertex(1), 1,
440 : aaPos, aaNorm, aaInt))
441 : {
442 : LOG("ObtainSimpleGrid failed. ignoring edge...\n");
443 : return NULL;
444 : }
445 :
446 : // calculate geometric-approximation-degree and triangle quality
447 0 : number approxDeg = GeometricApproximationDegree(sgSrc);
448 0 : number minApproxDeg = desc.approximation * approxDeg;
449 0 : number shapeDeg = ShapeQualityDegree(sgSrc);
450 :
451 : // perform a collapse on the simple grid
452 : SimpleGrid sgDest;
453 0 : if(!ObtainSimpleGrid_CollapseEdge(sgDest, grid, e,
454 : 1, aaPos, aaNorm, aaInt))
455 : {
456 : LOG("collapse edge failed...\n");
457 : return NULL;
458 : }
459 :
460 : // get the positions of the old endpoints
461 : static const int numTestPositions = 3;
462 : int newInd = 0;
463 0 : vector3 v[numTestPositions];
464 0 : v[0] = aaPos[e->vertex(0)];
465 0 : v[1] = aaPos[e->vertex(1)];
466 : v[2] = sgDest.vertices[newInd];
467 :
468 0 : vector3 n[numTestPositions];
469 0 : n[0] = aaNorm[e->vertex(0)];
470 0 : n[1] = aaNorm[e->vertex(1)];
471 : n[2] = sgDest.vertexNormals[newInd];
472 :
473 : // we'll compare 3 approximation degrees and three shape degrees
474 : number newApproxDeg[numTestPositions];
475 : number newShapeDeg[numTestPositions];
476 :
477 : // check which position is the best
478 : int bestIndex = -1;
479 : // the vertex subset index is used to support marks (crease and fixed vertices)
480 : int vrtSI[2];
481 : vrtSI[0] = vrtSI[1] = REM_NONE;
482 :
483 0 : if(pshMarks)
484 : {
485 0 : vrtSI[0] = pshMarks->get_subset_index(e->vertex(0));
486 0 : vrtSI[1] = pshMarks->get_subset_index(e->vertex(1));
487 :
488 0 : if((vrtSI[0] == REM_FIXED) || ((vrtSI[0] != REM_NONE) && (vrtSI[1] == REM_NONE))){
489 : bestIndex = 0;
490 : sgDest.vertices[newInd] = v[0];
491 : sgDest.vertexNormals[newInd] = n[0];
492 0 : newApproxDeg[0] = GeometricApproximationDegree(sgDest);
493 0 : newShapeDeg[0] = ShapeQualityDegree(sgDest);
494 : }
495 0 : else if((vrtSI[1] == REM_FIXED) || ((vrtSI[1] != REM_NONE) && (vrtSI[0] == REM_NONE))){
496 : bestIndex = 1;
497 : sgDest.vertices[newInd] = v[1];
498 : sgDest.vertexNormals[newInd] = n[1];
499 0 : newApproxDeg[1] = GeometricApproximationDegree(sgDest);
500 0 : newShapeDeg[1] = ShapeQualityDegree(sgDest);
501 : }
502 : }
503 :
504 : if(bestIndex == -1){
505 : // check all three approximation degrees
506 0 : for(int i = 0; i < numTestPositions; ++i){
507 : // we'll compute all qualities with the averaged normal
508 : sgDest.vertices[newInd] = v[i];
509 : sgDest.vertexNormals[newInd] = n[i];
510 0 : CalculateTriangleNormals(sgDest);
511 0 : newApproxDeg[i] = GeometricApproximationDegree(sgDest);
512 0 : newShapeDeg[i] = ShapeQualityDegree(sgDest);
513 : }
514 : // get the best one
515 : bestIndex = -1;
516 :
517 0 : for(int i = 0; i < numTestPositions; ++i){
518 0 : if(newApproxDeg[i] >= minApproxDeg){
519 0 : if(bestIndex == -1)
520 : bestIndex = i;
521 0 : else if(newShapeDeg[i] > newShapeDeg[bestIndex])
522 : bestIndex = i;
523 : }
524 : }
525 :
526 : // for(int i = 1; i < numTestPositions; ++i){
527 : // if(newShapeDeg[i] > newShapeDeg[bestIndex])
528 : // bestIndex = i;
529 : // }
530 : }
531 :
532 0 : if(bestIndex == -1)
533 : return NULL;
534 :
535 : // if the shape-degree of the collapsed region is too bad, we'll skip the collapse
536 0 : if(newShapeDeg[bestIndex] < desc.triQuality * shapeDeg)
537 : return NULL;
538 : /*
539 : // the approximation degree is only interesting if both endpoints of the
540 : // edge are regular surface vertices
541 : bool regularNeighbourhood = IsRegularSurfaceVertex(grid, e->vertex(0)) &&
542 : IsRegularSurfaceVertex(grid, e->vertex(1));
543 : */
544 :
545 : // if the best approximation degree is not too bad, we'll perform the collapse
546 0 : if(/*!regularNeighbourhood || */(newApproxDeg[bestIndex] >= minApproxDeg))
547 : {
548 : // pick one of the endpoints to be the one that resides
549 : // This has to be done with care, since the residing vertex
550 : // determines which edges will be removed and which will remain
551 : // (this is important, since we want crease edges to remain!).
552 : // If only one of the endpoints lies on a crease or is fixed, the
553 : // decision is easy (the marked one has to remain).
554 : // If both are marked, we have to investigate the triangles which
555 : // are adjacent to the collapsed edge (quadrilaterals don't make
556 : // problems here).
557 : // - If all three corner vertices are marked we have to distinguish:
558 : // (a) all three edges are marked -> abort
559 : // (b) two edges are marked -> the vertex connecting both remains.
560 : // (c) one edge is marked -> the chosen vertex doesn't affect edge-marks.
561 : // - If only the two are marked they don't affect the edge-marks.
562 : //
563 : // Even more care has to be taken if a fixed vertex is involved.
564 : // If case (b) applies and the creas-vertex is has to remain, it
565 : // has to be moved to the position of the fixed-vertex and has to
566 : // be marked fixed itself.
567 :
568 : // choose the vertex that shall remain.
569 0 : Vertex* vrt = e->vertex(0);
570 :
571 0 : if(vrtSI[0] != REM_FIXED && vrtSI[1] != REM_NONE)
572 : {
573 0 : vrt = e->vertex(1);
574 : }
575 :
576 0 : if(vrtSI[0] != REM_NONE && vrtSI[1] != REM_NONE){
577 : // both are marked. Things are getting complicated now!
578 : // get adjacent faces of e
579 : vector<Face*> vFaces;
580 0 : vFaces.reserve(2);
581 0 : CollectFaces(vFaces, grid, e);
582 :
583 : vector<Edge*> vEdges;
584 0 : vEdges.reserve(4);
585 0 : for(size_t i = 0; i < vFaces.size(); ++i){
586 0 : Face* f = vFaces[i];
587 : // only triangles are interesting
588 0 : if(f->num_edges() == 3){
589 : // get the number of marked edges
590 0 : CollectEdges(vEdges, grid, f);
591 : // count the marked edges
592 : int numMarked = 0;
593 0 : for(size_t j = 0; j < vEdges.size(); ++j){
594 : // note that pshMarks exists since vrtSI != REM_NONE
595 0 : if(pshMarks->get_subset_index(vEdges[j]) != REM_NONE)
596 0 : ++numMarked;
597 : }
598 :
599 0 : if(numMarked == 3){
600 : return NULL; // case (a) applies
601 : }
602 0 : else if(numMarked == 2){
603 : // check which of the vrts is connected to two
604 : // marked edges of the triangle
605 0 : for(size_t j = 0; j < 2; ++j){
606 : int numMarked = 0;
607 0 : for(size_t k = 0; k < vEdges.size(); ++k){
608 0 : if(pshMarks->get_subset_index(vEdges[k]) != REM_NONE){
609 0 : if(EdgeContains(vEdges[k], e->vertex(j)))
610 0 : ++numMarked;
611 : }
612 : }
613 : // if numMarked == 2 we found it
614 0 : if(numMarked == 2){
615 : // the connected edge has to be marked as a crease
616 0 : Edge* ce = GetConnectedEdge(grid, e->vertex(j), f);
617 0 : if(ce)
618 0 : pshMarks->assign_subset(ce, REM_CREASE);
619 : // we're done. break
620 : break;
621 : }
622 : }
623 : }
624 : else{
625 : }
626 : }
627 : }
628 0 : }
629 :
630 : // collapse the edge
631 0 : CollapseEdge(grid, e, vrt);
632 :
633 : // assign best position
634 : aaPos[vrt] = v[bestIndex];
635 : // assign the normal
636 : aaNorm[vrt] = n[bestIndex];
637 : /*
638 : if(pCandidates){
639 : //TODO: all edges that belong to associated faces are new candidates.
640 : // associated edges of vrt are possible new candidates
641 : pCandidates->select(grid.associated_edges_begin(vrt),
642 : grid.associated_edges_end(vrt));
643 : }*/
644 :
645 0 : return vrt;
646 : }
647 0 : }
648 :
649 : return NULL;
650 : }
651 :
652 :
653 : ////////////////////////////////////////////////////////////////////////
654 : template <class TAAPosVRT, class TAANormVRT, class TAAIntVRT>
655 0 : static bool PerformCollapses(const AdjustEdgeLengthDesc& desc, Grid& grid,
656 : SubsetHandler& shMarks, EdgeSelector& esel,
657 : TAAPosVRT& aaPos, TAANormVRT& aaNorm, TAAIntVRT& aaInt)
658 : {
659 : PROFILE_FUNC();
660 : LOG(" performing collapses\n");
661 : vector<Edge*> edges;
662 : int numCollapses = 0;
663 : // compare squares
664 0 : number minEdgeLenSq = sq(desc.minEdgeLen);
665 :
666 0 : while(!esel.empty())
667 : {
668 : Edge* e = *esel.begin<Edge>();
669 0 : esel.deselect(e);
670 :
671 : // the higher the curvature the larger minEdgeLen (lenFac is the inverse).
672 : // number dot = max(SMALL, VecDot(aaNorm[e->vertex(0)], aaNorm[e->vertex(1)]));
673 : // number lenFacSq = sq(desc.adaptivity * dot + (1. - desc.adaptivity));
674 : // adaptivity doesn't work at the moment...
675 : number lenFacSq = 1;
676 :
677 : // check whether the edge is short enough
678 0 : if(lenFacSq * VecDistanceSq(aaPos[e->vertex(0)], aaPos[e->vertex(1)]) < minEdgeLenSq)
679 : {
680 0 : Vertex* vrt = TryCollapse(desc, grid, e, aaPos, aaNorm, aaInt, &shMarks, &esel);
681 0 : if(vrt){
682 0 : ++numCollapses;
683 :
684 : // we'll deselect associated edges of vrt, to avoid cascade-collapses
685 : CollectAssociated(edges, grid, vrt);
686 : esel.deselect(edges.begin(), edges.end());
687 : }
688 : }
689 : }
690 0 : LOG(" collapses performed: " << numCollapses << endl);
691 :
692 0 : return true;
693 0 : }
694 :
695 : ////////////////////////////////////////////////////////////////////////
696 : template <class TAAPosVRT, class TAANormVRT>
697 0 : static bool TrySplit(Grid& grid, Edge* e, TAAPosVRT& aaPos, TAANormVRT& aaNorm,
698 : EdgeSelector* pCandidates = NULL, SubsetHandler* pshMarks = NULL)
699 : {
700 : // bool bCreaseEdge = false;
701 : // splits are not allowed for fixed edges
702 0 : if(pshMarks){
703 0 : if(pshMarks->get_subset_index(e) == REM_FIXED)
704 : return false;
705 : // else if(pshMarks->get_subset_index(e) == REM_CREASE)
706 : // bCreaseEdge = true;
707 : }
708 :
709 : // get the center of the edges
710 0 : vector3 vCenter = CalculateCenter(e, aaPos);
711 :
712 : // the new normal
713 : vector3 n;
714 0 : VecAdd(n, aaNorm[e->vertex(0)], aaNorm[e->vertex(1)]);
715 0 : VecNormalize(n, n);
716 :
717 : // if(bCreaseEdge){
718 : // // interpolating the normal can cause severe problems at highly
719 : // // irregular vertices or if one vertecs lies on a very
720 : // // sharp edge (the normals of the endpoints thus point
721 : // // in different directions.)
722 : // // VecAdd(n, aaNorm[e->vertex(0)], aaNorm[e->vertex(1)]);
723 : // // VecNormalize(n, n);
724 : // CalculateNormal(n, grid, e, aaPos);
725 : // }
726 :
727 : // split the edge
728 : RegularVertex* vrt = SplitEdge<RegularVertex>(grid, e, false);
729 :
730 : // assign the new position
731 : aaPos[vrt] = vCenter;
732 :
733 : // assign the new normal. calculate it if required
734 : // if(!bCreaseEdge)
735 : // CalculateVertexNormal(n, grid, vrt, aaPos);
736 :
737 : aaNorm[vrt] = n;
738 :
739 : // associated edges of vrt are candidates again
740 : // if(pCandidates)
741 : // pCandidates->select(grid.associated_edges_begin(vrt),
742 : // grid.associated_edges_end(vrt));
743 :
744 0 : return true;
745 : }
746 :
747 : ////////////////////////////////////////////////////////////////////////
748 : template <class TAAPosVRT, class TAANormVRT>
749 0 : static bool PerformSplits(const AdjustEdgeLengthDesc& desc, Grid& grid,
750 : SubsetHandler& shMarks, EdgeSelector& esel,
751 : TAAPosVRT& aaPos, TAANormVRT& aaNorm)
752 : {
753 : PROFILE_FUNC();
754 : // compare squares
755 0 : number maxEdgeLenSq = sq(desc.maxEdgeLen);
756 :
757 : LOG(" performing splits\n");
758 : int numSplits = 0;
759 :
760 0 : while(!esel.empty())
761 : {
762 : // get an edge
763 : Edge* e = *esel.begin<Edge>();
764 0 : esel.deselect(e);
765 :
766 : // the higher the curvature the smaller maxEdgeLen (lenFac is the inverse).
767 : // adaptivity doesn't work currently...
768 : // number dot = max(SMALL, VecDot(aaNorm[e->vertex(0)], aaNorm[e->vertex(1)]));
769 : // number lenFacSq = sq(desc.adaptivity * dot + (1. - desc.adaptivity));
770 : number lenFacSq = 1;
771 :
772 : // check whether the edge should be splitted
773 0 : if(VecDistanceSq(aaPos[e->vertex(0)], aaPos[e->vertex(1)]) > lenFacSq * maxEdgeLenSq)
774 : {
775 0 : if(TrySplit(grid, e, aaPos, aaNorm, &esel, &shMarks))
776 0 : ++numSplits;
777 : }
778 : }
779 :
780 0 : LOG(" splits performed: " << numSplits << endl);
781 0 : return true;
782 : }
783 :
784 : ////////////////////////////////////////////////////////////////////////
785 : // relocate point by smoothing
786 : // static void RelocatePointBySmoothing(vector3& vOut, const vector3&v,
787 : // std::vector<vector3>& vNodes,
788 : // size_t numIterations, number stepSize,
789 : // std::vector<number>* pvWeights = NULL)
790 : // {
791 : // if(vNodes.empty())
792 : // return;
793 :
794 : // //TODO: incorporate weights
795 : // // iterate through all nodes, calculate the direction
796 : // // from v to each node, and add the scaled sum to v.
797 : // vector3 t, vOld;
798 : // vOut = v;
799 : // stepSize /= (number)vNodes.size();
800 :
801 : // for(size_t j = 0; j < numIterations; ++j){
802 : // vOld = vOut;
803 : // for(size_t i = 0; i < vNodes.size(); ++i){
804 : // VecSubtract(t, vNodes[i], vOld);
805 : // VecScale(t, t, stepSize);
806 : // VecAdd(vOut, vOut, t);
807 : // }
808 : // }
809 : // }
810 :
811 : ////////////////////////////////////////////////////////////////////////
812 : // FixBadTriangles
813 : // template <class TAAPosVRT, class TAANormVRT>
814 : // static bool FixBadTriangles(Grid& grid, SubsetHandler& shMarks, EdgeSelector& esel,
815 : // TAAPosVRT& aaPos, TAANormVRT& aaNorm,
816 : // number qualityThreshold)
817 : // {
818 : // LOG(" fixing bad triangles... ");
819 : // // iterate through marked edges and check whether adjacent triangles are
820 : // // badly shaped. If this is the case, then split the non-marked
821 : // // edges of the triangle.
822 : // // store newly inserted vertices and smooth them at the end of the algo
823 : // vector<Vertex*> vNewVrts;
824 : // vector<Face*> vFaces;
825 : // vector<Edge*> vEdges;
826 :
827 : // // we wont assign this selector to the grid until it is clear that
828 : // // we'll need it.
829 : // Selector sel;
830 :
831 :
832 : // EdgeIterator iter = esel.begin<Edge>();
833 : // while(iter != esel.end<Edge>())
834 : // {
835 : // // store edge and increase iterator immediatly
836 : // Edge* e = *iter;
837 : // ++iter;
838 :
839 : // // get the adjacent faces
840 : // CollectFaces(vFaces, grid, e);
841 :
842 : // // check whether one of them is degenerated.
843 : // for(size_t i = 0; i < vFaces.size(); ++i){
844 : // Face* f = vFaces[i];
845 : // number q = FaceQuality(f, aaPos);
846 : // // if the quality is smaller than threshold, we have to
847 : // // do something
848 : // if(q < qualityThreshold){
849 : // // make sure that the selector is connected to the grid
850 : // if(sel.grid() == NULL)
851 : // sel.assign_grid(grid);
852 :
853 : // // get associated edges and mark them for refinement
854 : // CollectEdges(vEdges, grid, f);
855 :
856 : // for(size_t j = 0; j < vEdges.size(); ++j){
857 : // if(shMarks.get_subset_index(vEdges[j]) == REM_NONE)
858 : // sel.select(vEdges[j]);
859 : // }
860 : // }
861 : // }
862 : // }
863 :
864 : // // if edges have been selected, we'll now call refine.
865 : // if(sel.grid() != NULL){
866 : // if(sel.num<Edge>() > 0){
867 : // if(Refine(grid, sel)){
868 : // LOG(sel.num<Vertex>() << " new vertices... ");
869 : // // retriangulate surface
870 : // if(grid.num<Quadrilateral>() > 0)
871 : // Triangulate(grid, grid.begin<Quadrilateral>(), grid.end<Quadrilateral>());
872 :
873 : // // calculate normals, then
874 : // // smooth new vertices (all vertices selected in sel).
875 : // vector<Vertex*> vNeighbours;
876 : // vector<vector3> vNodes;
877 :
878 : // // calculate normals
879 : // for(VertexIterator iter = sel.begin<Vertex>();
880 : // iter != sel.end<Vertex>(); ++iter)
881 : // {
882 : // // collect neighbour nodes
883 : // Vertex* vrt = *iter;
884 : // CollectNeighbors(vNeighbours, grid, vrt);
885 :
886 : // // sum their normals and interpolate it
887 : // // make sure to only add valid normals
888 : // vector3 n(0, 0, 0);
889 : // for(size_t i = 0; i < vNeighbours.size(); ++i){
890 : // if(!sel.is_selected(vNeighbours[i]))
891 : // VecAdd(n, n, aaNorm[vNeighbours[i]]);
892 : // }
893 : // VecNormalize(aaNorm[vrt], n);
894 : // }
895 :
896 : // // repeat smoothing.
897 : // for(size_t i = 0; i < 5; ++i){
898 : // for(VertexIterator iter = sel.begin<Vertex>();
899 : // iter != sel.end<Vertex>(); ++iter)
900 : // {
901 : // Vertex* vrt = *iter;
902 :
903 : // // collect the neighbours and project them to the plane
904 : // // that is defined by vrt and its normal
905 : // vector3 v = aaPos[vrt];
906 : // vector3 n = aaNorm[vrt];
907 :
908 : // CollectNeighbors(vNeighbours, grid, vrt);
909 : // vNodes.resize(vNeighbours.size());
910 :
911 : // for(size_t j = 0; j < vNodes.size(); ++j)
912 : // ProjectPointToPlane(vNodes[j], aaPos[vNeighbours[j]], v, n);
913 :
914 : // // perform point relocation
915 : // RelocatePointBySmoothing(aaPos[vrt], v, vNodes, 5, 0.05);
916 : // }
917 : // }
918 : // }
919 : // else{
920 : // LOG("refine failed!\n");
921 : // return false;
922 : // }
923 : // }
924 : // }
925 :
926 : // LOG("done\n");
927 : // return true;
928 : // }
929 :
930 : ////////////////////////////////////////////////////////////////////////
931 : // PerformSmoothing
932 : // template <class TAAPosVRT, class TAANormVRT>
933 : // static void PerformSmoothing(Grid& grid, SubsetHandler& shMarks,
934 : // TAAPosVRT& aaPos, TAANormVRT& aaNorm,
935 : // size_t numIterations, number stepSize)
936 : // {
937 : // vector<vector3> vNodes;
938 : // vector<Vertex*> vNeighbours;
939 : // for(size_t i = 0; i < numIterations; ++i){
940 : // CalculateVertexNormals(grid, aaPos, aaNorm);
941 : // for(VertexIterator iter = grid.begin<Vertex>();
942 : // iter != grid.end<Vertex>(); ++iter)
943 : // {
944 : // Vertex* vrt = *iter;
945 : // // if the vertex is fixed then leave it where it is.
946 : // if(shMarks.get_subset_index(vrt) == REM_FIXED)
947 : // continue;
948 : // /*
949 : // if(shMarks.get_subset_index(vrt) == REM_CREASE)
950 : // continue;
951 : // */
952 : // // collect the neighbours and project them to the plane
953 : // // that is defined by vrt and its normal
954 : // vector3 v = aaPos[vrt];
955 : // vector3 n = aaNorm[vrt];
956 :
957 : // bool bProjectPointsToPlane = true;
958 :
959 : // if(shMarks.get_subset_index(vrt) == REM_CREASE){
960 : // CollectNeighbors(vNeighbours, grid, vrt, NHT_EDGE_NEIGHBORS,
961 : // IsInSubset(shMarks, REM_CREASE));
962 :
963 : // // we have to choose a special normal
964 : // if(vNeighbours.size() != 2){
965 : // UG_LOG("n"<<vNeighbours.size());
966 : // continue;
967 : // }
968 : // else{
969 : // vector3 v1, v2;
970 : // VecSubtract(v1, v, aaPos[vNeighbours[0]]);
971 : // VecSubtract(v2, v, aaPos[vNeighbours[1]]);
972 : // VecNormalize(v1, v1);
973 : // VecNormalize(v2, v2);
974 : // VecAdd(n, v1, v2);
975 : // if(VecLengthSq(n) > SMALL)
976 : // VecNormalize(n, n);
977 : // else {
978 : // // both edges have the same direction.
979 : // // don't project normals
980 : // bProjectPointsToPlane = false;
981 : // }
982 : // }
983 : // }
984 : // else{
985 : // CollectNeighbors(vNeighbours, grid, vrt);
986 : // }
987 :
988 : // vNodes.resize(vNeighbours.size());
989 :
990 : // if(bProjectPointsToPlane){
991 : // for(size_t j = 0; j < vNodes.size(); ++j)
992 : // ProjectPointToPlane(vNodes[j], aaPos[vNeighbours[j]], v, n);
993 : // }
994 : // else{
995 : // for(size_t j = 0; j < vNodes.size(); ++j)
996 : // vNodes[j] = aaPos[vNeighbours[j]];
997 : // }
998 : // // perform point relocation
999 : // RelocatePointBySmoothing(aaPos[vrt], v, vNodes, 5, stepSize);
1000 : // }
1001 : // }
1002 : // }
1003 :
1004 : /** Make sure that elements in gridOut directly correspond to
1005 : * elements in gridIn*/
1006 : template <class TGeomObj>
1007 : static void CopySelectionStatus(Selector& selOut, Grid& gridOut,
1008 : Selector& selIn, Grid& gridIn)
1009 : {
1010 : typedef typename geometry_traits<TGeomObj>::iterator GeomObjIter;
1011 : GeomObjIter endOut = gridOut.end<TGeomObj>();
1012 : GeomObjIter endIn = gridIn.end<TGeomObj>();
1013 : GeomObjIter iterOut = gridOut.begin<TGeomObj>();
1014 : GeomObjIter iterIn = gridIn.begin<TGeomObj>();
1015 :
1016 : for(; (iterOut != endOut) && (iterIn != endIn); ++iterOut, ++iterIn)
1017 : {
1018 : if(selIn.is_selected(*iterIn))
1019 : selOut.select(*iterOut);
1020 : }
1021 : }
1022 :
1023 : ////////////////////////////////////////////////////////////////////////
1024 0 : bool AdjustEdgeLength(Grid& grid, SubsetHandler& shMarks,
1025 : const AdjustEdgeLengthDesc& desc, int numIterations)
1026 : {
1027 : PROFILE_FUNC();
1028 :
1029 : //TODO: replace this by a template parameter
1030 : APosition aPos = aPosition;
1031 :
1032 : // we have to make sure that the mesh consist of triangles only,
1033 : // since the algorithm would produce bad results if not.
1034 0 : if(grid.num<Quadrilateral>() > 0)
1035 : {
1036 : UG_LOG(" INFO: grid contains quadrilaterals. Converting to triangles...\n");
1037 :
1038 : //TODO: not gridIn but a copy of gridIn (pRefGrid) should be triangulated.
1039 0 : Triangulate(grid, grid.begin<Quadrilateral>(),
1040 0 : grid.end<Quadrilateral>());
1041 : }
1042 :
1043 : // make sure that grid and pRefGrid have position-attachments
1044 0 : if(!grid.has_vertex_attachment(aPos)){
1045 : UG_LOG(" vertex-position-attachment missing in AdjustEdgeLength. Aborting.\n");
1046 : return false;
1047 : }
1048 :
1049 : // make sure that faces create associated edges
1050 0 : if(!grid.option_is_enabled(FACEOPT_AUTOGENERATE_EDGES))
1051 : {
1052 : LOG(" INFO: auto-enabling FACEOPT_AUTOGENERATE_EDGES in AdjustEdgeLength.\n");
1053 0 : grid.enable_options(FACEOPT_AUTOGENERATE_EDGES);
1054 : }
1055 :
1056 : // position attachment
1057 : Grid::VertexAttachmentAccessor<APosition> aaPos(grid, aPos);
1058 :
1059 : // we need an integer attachment (a helper for ObtainSimpleGrid)
1060 : AInt aInt;
1061 0 : grid.attach_to<Vertex>(aInt);
1062 : Grid::AttachmentAccessor<Vertex, AInt> aaInt(grid, aInt);
1063 :
1064 : // TODO: normals shouldn't be created here but instead be passed to the method.
1065 : // attach the vertex normals.
1066 : ANormal aNorm;
1067 0 : grid.attach_to<Vertex>(aNorm);
1068 : Grid::AttachmentAccessor<Vertex, ANormal> aaNorm(grid, aNorm);
1069 0 : CalculateVertexNormals(grid, aPos, aNorm);
1070 :
1071 : // assign vertex marks
1072 0 : AssignFixedVertices(grid, shMarks);
1073 0 : AssignCreaseVertices(grid, shMarks);
1074 :
1075 : // we need an selector that holds all edges that are candidates for a operation
1076 0 : EdgeSelector esel(grid);
1077 0 : esel.enable_selection_inheritance(false);
1078 :
1079 : // sort the triangles of pRefGrid into a octree to speed-up projection performance
1080 : SPOctree octree;
1081 0 : node_tree::Traverser_ProjectPoint pojectionTraverser;
1082 0 : if(desc.projectPoints){
1083 : //PROFILE_BEGIN(octree_construction);
1084 0 : octree = CreateOctree(grid, grid.begin<Triangle>(),
1085 0 : grid.end<Triangle>(),
1086 0 : 10, 30, false, aPos);
1087 :
1088 : //PROFILE_END();
1089 0 : if(!octree.valid()){
1090 : UG_LOG(" Octree creation failed in AdjustEdgeLength. Aborting.\n");
1091 : return false;
1092 : }
1093 : }
1094 :
1095 : // start the main iteration
1096 0 : for(int iteration = 0; iteration < numIterations; ++iteration)
1097 : {
1098 : // perform splits
1099 0 : esel.select(grid.begin<Edge>(), grid.end<Edge>());
1100 0 : if(!PerformSplits(desc, grid, shMarks, esel, aaPos, aaNorm))
1101 : return false;
1102 :
1103 : // // perform collapses
1104 0 : esel.select(grid.begin<Edge>(), grid.end<Edge>());
1105 0 : if(!PerformCollapses(desc, grid, shMarks, esel, aaPos, aaNorm, aaInt))
1106 : return false;
1107 :
1108 : // perform swaps
1109 : // esel.select(grid.begin<Edge>(), grid.end<Edge>());
1110 : // if(!PerformSwaps(desc, grid, shMarks, esel, aaPos, aaNorm, aaInt))
1111 : // return false;
1112 :
1113 : /*
1114 : // This is commented out, since it didn't help with the problems encountered
1115 : // in the geometries at that time.
1116 : // The algorithm however works and may prove useful in the future.
1117 : // fix bad triangles
1118 : // adjacent to crease-edges badly shaped triangles may occur,
1119 : // which have to be treated somehow.
1120 : esel.clear();
1121 : esel.select(shMarks.begin<Edge>(REM_CREASE), shMarks.end<Edge>(REM_CREASE));
1122 : esel.select(shMarks.begin<Edge>(REM_FIXED), shMarks.end<Edge>(REM_FIXED));
1123 : FixBadTriangles(grid, shMarks, esel, aaPos, aaNorm, 0.1);
1124 : */
1125 : // relocate points
1126 : /*LOG(" smoothing points...");
1127 : PerformSmoothing(grid, shMarks, aaPos, aaNorm, 10, 0.1);
1128 : LOG(" done\n");*/
1129 :
1130 : // LOG(" updating normals...\n");
1131 : // CalculateVertexNormals(grid, aPos, aNorm);
1132 :
1133 : // project points back on the surface
1134 0 : if(desc.projectPoints)
1135 : {
1136 : LOG(" projecting points...");
1137 : //PROFILE_BEGIN(projecting_points);
1138 : for(VertexIterator iter = grid.vertices_begin();
1139 0 : iter != grid.vertices_end(); ++iter)
1140 : {
1141 : //TODO: project crease vertices onto creases only! Don't project fixed vertices
1142 0 : if(shMarks.get_subset_index(*iter) != REM_FIXED){
1143 : vector3 vNew;
1144 0 : if(pojectionTraverser.project(aaPos[*iter], octree/*, &aaNorm[*iter]*/)){
1145 : aaPos[*iter] = pojectionTraverser.get_closest_point();
1146 : }
1147 : else{
1148 : LOG("f");
1149 : }
1150 : }
1151 : }
1152 : //PROFILE_END();
1153 : LOG(" done\n");
1154 : }
1155 :
1156 : // if(iteration < numIterations - 1){
1157 : // LOG(" updating normals...");
1158 : // CalculateVertexNormals(grid, aPos, aNorm);
1159 : // }
1160 : }
1161 :
1162 :
1163 : // detach
1164 : grid.detach_from<Vertex>(aInt);
1165 : grid.detach_from<Vertex>(aNorm);
1166 :
1167 : GPLOTSAVE();
1168 : return true;
1169 0 : }
1170 :
1171 :
1172 :
1173 :
1174 : /*
1175 : ////////////////////////////////////////////////////////////////////////
1176 : // This is an alternative version for PerformSplits.
1177 : // It uses the Refine method.
1178 : // It is not yet optimized for maximum speed.
1179 : // While it performs a little less splits, overall runtime of
1180 : // AdjustEdgeLength is not better than with the original
1181 : // PerformSplits method.
1182 : template <class TAAPosVRT, class TAANormVRT>
1183 : bool PerformSplits(Grid& grid, SubsetHandler& shMarks, EdgeSelector& esel,
1184 : number maxEdgeLen, TAAPosVRT& aaPos, TAANormVRT& aaNorm)
1185 : {
1186 : AInt aInt;
1187 : grid.attach_to_edges(aInt);
1188 :
1189 : // compare squares
1190 : maxEdgeLen *= maxEdgeLen;
1191 :
1192 : Selector sel(grid);
1193 : sel.enable_autoselection(true);
1194 : sel.select(esel.begin<Edge>(), esel.end<Edge>());
1195 :
1196 : LOG(" performing splits\n");
1197 : int numSplits = 0;
1198 :
1199 : while(!sel.empty()){
1200 : // deselect all vertices and faces
1201 : sel.clear_selection<Vertex>();
1202 : sel.clear_selection<Face>();
1203 : sel.clear_selection<Volume>();
1204 :
1205 : // deselect all edges that shall not be splitted
1206 : EdgeIterator iter = sel.begin<Edge>();
1207 : while(iter != sel.end<Edge>()){
1208 : Edge* e = *iter;
1209 : ++iter;
1210 :
1211 : // the higher the curvature the smaller the maxEdgeLen.
1212 : // minimal lenFac is 0.1
1213 : number lenFac = CalculateLengthFac(grid, shMarks, e, aaPos);
1214 :
1215 : // fixed edges will not be refined
1216 : if(shMarks.get_subset_index(e) == REM_FIXED)
1217 : sel.deselect(e);
1218 : else if(VecDistanceSq(aaPos[e->vertex(0)], aaPos[e->vertex(1)]) < lenFac * maxEdgeLen)
1219 : sel.deselect(e);
1220 : }
1221 :
1222 : // refine the grid
1223 : Refine(grid, sel, aInt);
1224 :
1225 : // new vertices are selected
1226 : numSplits += sel.num<Vertex>();
1227 :
1228 : // re-triangulate
1229 : Triangulate(grid, &aaPos);
1230 :
1231 : // calculate normal for new vertices
1232 : //TODO: be careful with crease edges
1233 : for(VertexIterator iter = sel.begin<Vertex>();
1234 : iter != sel.end<Vertex>(); ++iter)
1235 : {
1236 : CalculateVertexNormal(aaNorm[*iter], grid, *iter, aaPos);
1237 : }
1238 :
1239 : sel.clear_selection<Vertex>();
1240 : sel.clear_selection<Face>();
1241 : sel.clear_selection<Volume>();
1242 : }
1243 :
1244 : grid.detach_from_edges(aInt);
1245 : LOG(" splits performed: " << numSplits << endl);
1246 : return true;
1247 : }
1248 : */
1249 :
1250 : }// end of namespace
|
