Branch data Line data Source code
1 : : /*
2 : : * QslimDecimation.cpp
3 : : *
4 : : * Created on: Mar 10, 2010
5 : : * Author: iulian
6 : : */
7 : :
8 : : #include <assert.h>
9 : : #include "QslimDecimation.hpp"
10 : :
11 : : // for proximity searches
12 : : #include "moab/AdaptiveKDTree.hpp"
13 : : #include "moab/ReadUtilIface.hpp"
14 : :
15 : : #include "Mat4.h"
16 : : #include "defs.h"
17 : : #include "quadrics.h"
18 : : #include <time.h>
19 : : #include <map>
20 : : #include <limits>
21 : :
22 : : // those are used in model navigation/simplification
23 : : #include "primitives.h"
24 : :
25 : : // this is the global thing, that everybody will use
26 : : moab::Interface * mb;
27 : : moab::Tag uniqIDtag; // this will be used to mark vertices moab::EntityHandles
28 : : moab::Tag validTag;
29 : :
30 : : moab::Tag costTag; // simplification induces an error cost at each vertex
31 : : // try to keep adding the cost, to see if it is spreading nicely
32 : :
33 : : // this will be used to store plane data for each triangle, as 4 doubles
34 : : // will be updated only if needed ( -m option == opts.will_preserve_mesh_quality)
35 : : moab::Tag planeDataTag;
36 : :
37 : 78 : moab::Range verts; // original list of vertices, that are part of the original triangles
38 : 78 : moab::Range triangles;
39 : 78 : moab::Range edgs;
40 : 78 : QslimOptions opts; // external
41 : : moab::EntityHandle iniSet;
42 : :
43 : 227992 : int uniqID(moab::EntityHandle v) {
44 : : int val;
45 [ + - ]: 227992 : moab::ErrorCode rval = mb->tag_get_data(uniqIDtag, &v, 1, &val);
46 : : (void) rval;
47 [ - + ]: 227992 : assert(rval==moab::MB_SUCCESS);
48 : 227992 : return val;
49 : : }
50 : : // the vertices are not deleted anymore, just invalidated
51 : : // the edges are deleted, though, and triangles
52 : 44424 : int ehIsValid(moab::EntityHandle v) {
53 : : unsigned char val;
54 [ + - ]: 44424 : moab::ErrorCode rval = mb->tag_get_data(validTag, &v, 1, &val);
55 : : (void) rval;
56 [ - + ]: 44424 : assert(rval==moab::MB_SUCCESS);
57 : 44424 : return (int) val;
58 : : }
59 : :
60 : : // include here the main classes used for decimation
61 : :
62 : : #include "Heap.hpp"
63 : : // prox grid is used for proximity grid only
64 : : //#include "ProxGrid.h"
65 : :
66 : : class pair_info: public Heapable {
67 : : public:
68 : : moab::EntityHandle v0, v1; // Vertex *v0, *v1;
69 : :
70 : : Vec3 candidate;
71 : : double cost;
72 : :
73 : : //pair_info ( Vertex *a,Vertex *b ) { v0=a; v1=b; cost=HUGE; }
74 : 30400 : pair_info(moab::EntityHandle a, moab::EntityHandle b) {
75 : 15200 : v0 = a;
76 : 15200 : v1 = b;
77 : 15200 : cost = std::numeric_limits<double>::max();
78 : 15200 : }
79 : :
80 : 1088 : bool isValid() {
81 [ + - ][ + - ]: 1088 : return ehIsValid(v0) && ehIsValid(v1);
82 : : }// :v0->isValid() && v1->isValid(); }
83 : : };
84 : :
85 : : typedef buffer<pair_info *> pair_buffer;
86 : :
87 : 10404 : class vert_info {
88 : : public:
89 : :
90 : : pair_buffer pairs;
91 : :
92 : : Mat4 Q;
93 : : double norm;
94 : :
95 : 5202 : vert_info() :
96 [ + - ]: 5202 : Q(Mat4::zero) {
97 [ + - ]: 5202 : pairs.init(2);
98 : 5202 : norm = 0.0;
99 : 5202 : }
100 : : };
101 : :
102 : : // these are
103 : : static Heap *heap;
104 : 78 : static array<vert_info> vinfo;
105 : : static double proximity_limit; // distance threshold squared
106 : :
107 : : int validFaceCount;
108 : : int validVertCount;
109 : :
110 : : ////////////////////////////////////////////////////////////////////////
111 : : //
112 : : // Low-level routines for manipulating pairs
113 : : //
114 : :
115 : 227992 : static inline vert_info& vertex_info(moab::EntityHandle v)//Vertex *v )
116 : : {
117 : : // moab::EntityHandle should return an integer tag with an
118 : : // index in the big array of vert_info
119 : : // something like: return tag
120 : : // for the time being, we can return the simple id...
121 : 227992 : return vinfo(uniqID(v));
122 : : }
123 : :
124 : : static
125 : 2404 : bool check_for_pair(moab::EntityHandle v0, moab::EntityHandle v1)//Vertex *v0, Vertex *v1 )
126 : : {
127 : 2404 : const pair_buffer& pairs = vertex_info(v0).pairs;
128 : :
129 [ + + ]: 12926 : for (int i = 0; i < pairs.length(); i++) {
130 [ + + ][ + + ]: 11524 : if (pairs(i)->v0 == v1 || pairs(i)->v1 == v1)
[ + + ]
131 : 1002 : return true;
132 : : }
133 : :
134 : 1402 : return false;
135 : : }
136 : :
137 : 15200 : static pair_info *new_pair(moab::EntityHandle v0, moab::EntityHandle v1)// Vertex *v0, Vertex *v1 )
138 : : {
139 [ + - ]: 15200 : vert_info& v0_info = vertex_info(v0);
140 [ + - ]: 15200 : vert_info& v1_info = vertex_info(v1);
141 : :
142 [ + - ][ + - ]: 15200 : pair_info *pair = new pair_info(v0, v1);
143 [ + - ]: 15200 : v0_info.pairs.add(pair);
144 [ + - ]: 15200 : v1_info.pairs.add(pair);
145 : :
146 : 15200 : return pair;
147 : : }
148 : :
149 : : static
150 : 1546 : void delete_pair(pair_info *pair) {
151 : 1546 : vert_info& v0_info = vertex_info(pair->v0);
152 : 1546 : vert_info& v1_info = vertex_info(pair->v1);
153 : :
154 : 1546 : v0_info.pairs.remove(v0_info.pairs.find(pair));
155 : 1546 : v1_info.pairs.remove(v1_info.pairs.find(pair));
156 : :
157 [ + + ]: 1546 : if (pair->isInHeap())
158 : 1002 : heap->kill(pair->getHeapPos());
159 : :
160 : 1546 : delete pair;
161 : 1546 : }
162 : :
163 : : ////////////////////////////////////////////////////////////////////////
164 : : //
165 : : // The actual guts of the algorithm:
166 : : //
167 : : // - pair_is_valid
168 : : // - compute_pair_info
169 : : // - do_contract
170 : : //
171 : : /*
172 : : static
173 : : bool pair_is_valid(moab::EntityHandle u, moab::EntityHandle v)// Vertex *u, Vertex *v )
174 : : {
175 : : //
176 : : Vec3 vu = getVec3FromMBVertex(mb, u);
177 : : Vec3 vv = getVec3FromMBVertex(mb, v);
178 : : return norm2(vu - vv) < proximity_limit;
179 : : //return norm2 ( *u - *v ) < proximity_limit;
180 : : }*/
181 : :
182 : : static
183 : 556758 : int predict_face(moab::EntityHandle tria, moab::EntityHandle v1,
184 : : moab::EntityHandle v2, /*Face& F, Vertex *v1, Vertex *v2,*/
185 : : Vec3& vnew, Vec3& f1, Vec3& f2, Vec3& f3) {
186 : 556758 : int nmapped = 0;
187 : : const moab::EntityHandle * conn;
188 : : int num_nodes;
189 [ + - ]: 556758 : moab::ErrorCode rval = mb->get_connectivity(tria, conn, num_nodes);
190 : : (void) rval;
191 [ + - ][ - + ]: 556758 : assert(3==num_nodes && rval == moab::MB_SUCCESS);
192 [ + + ][ + + ]: 556758 : if (conn[0] == v1 || conn[0] == v2) {
193 [ + - ]: 220630 : f1 = vnew;
194 : 220630 : nmapped++;
195 : : } else
196 [ + - ][ + - ]: 336128 : f1 = getVec3FromMBVertex(mb, conn[0]);
197 : :
198 [ + + ][ + + ]: 556758 : if (conn[1] == v1 || conn[1] == v2) {
199 [ + - ]: 219620 : f2 = vnew;
200 : 219620 : nmapped++;
201 : : } else
202 [ + - ][ + - ]: 337138 : f2 = getVec3FromMBVertex(mb, conn[1]);
203 : :
204 [ + + ][ + + ]: 556758 : if (conn[2] == v1 || conn[2] == v2) {
205 [ + - ]: 227508 : f3 = vnew;
206 : 227508 : nmapped++;
207 : : } else
208 [ + - ][ + - ]: 329250 : f3 = getVec3FromMBVertex(mb, conn[2]);
209 : :
210 : : // find vertices in face tria
211 : : /*
212 : : if ( F.vertex ( 0 ) == v1 || F.vertex ( 0 ) == v2 )
213 : : { f1 = vnew; nmapped++; }
214 : : else f1 = *F.vertex ( 0 );
215 : :
216 : : if ( F.vertex ( 1 ) == v1 || F.vertex ( 1 ) == v2 )
217 : : { f2 = vnew; nmapped++; }
218 : : else f2 = *F.vertex ( 1 );
219 : :
220 : : if ( F.vertex ( 2 ) == v1 || F.vertex ( 2 ) == v2 )
221 : : { f3 = vnew; nmapped++; }
222 : : else f3 = *F.vertex ( 2 );
223 : : */
224 : 556758 : return nmapped;
225 : : }
226 : :
227 : : #define MESH_INVERSION_PENALTY 1e9
228 : :
229 : : static
230 : 60616 : double pair_mesh_positivity(/* Model& M,*/moab::EntityHandle v1,
231 : : moab::EntityHandle v2, /*Vertex *v1, Vertex *v2,*/
232 : : Vec3& vnew) {
233 [ + - ]: 60616 : std::vector<moab::EntityHandle> changed;
234 : :
235 : : // : construct the list of faces influenced by the
236 : : // moving of vertices v1 and v2 into vnew
237 : : //M.contractionRegion ( v1, v2, changed );
238 [ + - ]: 60616 : moab::ErrorCode rval = contractionRegion(mb, v1, v2, changed);
239 [ - + ]: 60616 : if (rval != moab::MB_SUCCESS) {
240 [ # # ]: 0 : std::cout << "error in getting adjacency information vs: "
241 [ # # ][ # # ]: 0 : << mb->id_from_handle(v1) << " " << mb->id_from_handle(v2) << "\n";
[ # # ][ # # ]
[ # # ][ # # ]
242 : : }
243 : :
244 : : // double Nsum = 0;
245 [ - + ][ # # ]: 60616 : if (opts.logfile && opts.selected_output & OUTPUT_CONTRACTIONS)
246 [ # # ][ # # ]: 0 : *opts.logfile << " positivity for v1, v2: " << mb->id_from_handle(v1)
[ # # ]
247 [ # # ][ # # ]: 0 : << " " << mb->id_from_handle(v2) << std::endl;
[ # # ][ # # ]
248 : :
249 [ + + ]: 609338 : for (unsigned int i = 0; i < changed.size(); i++) {
250 : : //Face& F = *changed ( i );
251 [ + - ]: 556758 : moab::EntityHandle F = changed[i];
252 [ + - ][ + - ]: 556758 : Vec3 f1, f2, f3;
[ + - ]
253 : :
254 [ + - ]: 556758 : int nmapped = predict_face(F, v1, v2, vnew, f1, f2, f3);
255 : :
256 : : //
257 : : // Only consider non-degenerate faces
258 [ + + ]: 556758 : if (nmapped < 2) {
259 : : //Plane Pnew ( f1, f2, f3 );
260 : : #if 0
261 : : Vec3 normalNew = Pnew.normal();
262 : : if ( normalNew[Z] < positivityMin )
263 : : positivityMin=normalNew[Z]; // Z direction!!!
264 : : if (opts.logfile && opts.selected_output&OUTPUT_CONTRACTIONS )
265 : : *opts.logfile << "Triangle " << mb->id_from_handle(F)
266 : : << " nmapped " << nmapped << std::endl;
267 : : if (opts.logfile && positivityMin<=0 && opts.selected_output&OUTPUT_CONTRACTIONS )
268 : : *opts.logfile << "Triangle " << mb->id_from_handle(F)
269 : : << " normal Z:" << normalNew[Z] << std::endl;
270 : : #endif
271 [ + - ][ + - ]: 445758 : double positiv = (f2[0] - f1[0]) * (f3[1] - f1[1]) - (f2[1] - f1[1])
[ + - ][ + - ]
[ + - ][ + - ]
272 [ + - ][ + - ]: 445758 : * (f3[0] - f1[0]);
273 [ + + ]: 445758 : if (positiv <= 0) {
274 [ - + ][ # # ]: 8036 : if (opts.logfile && opts.selected_output & OUTPUT_CONTRACTIONS)
275 [ # # ][ # # ]: 0 : *opts.logfile << "Triangle " << mb->id_from_handle(F) << " nmapped "
[ # # ][ # # ]
276 [ # # ][ # # ]: 0 : << nmapped << " orient: " << positiv << std::endl;
[ # # ][ # # ]
277 : 445758 : return MESH_INVERSION_PENALTY * 10;
278 : : }
279 : : }
280 : : }
281 : :
282 : : //return (-Nmin) * MESH_INVERSION_PENALTY;
283 : :
284 : 60616 : return 0.0;
285 : : }
286 : :
287 : : // will make sure that no hole would be closed
288 : : /*
289 : : * will accomplish that by looking at the edges connected to both vertices
290 : : * if there would be edges that would merge without a connecting triangle,
291 : : * it would increase the cost (penalize it harshly)
292 : : */
293 : 0 : static double pair_mesh_topology(moab::EntityHandle v1, moab::EntityHandle v2)
294 : : {
295 : : // first, find nodes v3 that are connected to both v1 and v2 by edges
296 : : // if there is no triangle that is formed with v1, v2, v3, it means it would
297 : : // collapse a hole
298 [ # # ][ # # ]: 0 : std::vector<moab::EntityHandle> edges1, edges2;
299 [ # # ][ # # ]: 0 : moab::Range nodes1, nodes2;
300 [ # # ]: 0 : moab::ErrorCode rval = mb->get_adjacencies(&v1, 1, 1, false, edges1);
301 [ # # ]: 0 : assert (moab::MB_SUCCESS == rval);
302 [ # # ]: 0 : filterValid(mb, edges1);
303 [ # # ]: 0 : rval = mb->get_adjacencies(&v2, 1, 1, false, edges2);
304 [ # # ]: 0 : assert (moab::MB_SUCCESS == rval);
305 [ # # ]: 0 : filterValid(mb, edges2);
306 [ # # ][ # # ]: 0 : rval = mb->get_connectivity(&edges1[0], edges1.size(), nodes1);
307 [ # # ]: 0 : assert (moab::MB_SUCCESS == rval);
308 [ # # ][ # # ]: 0 : rval = mb->get_connectivity(&edges2[0], edges2.size(), nodes2);
309 [ # # ]: 0 : assert (moab::MB_SUCCESS == rval);
310 [ # # ]: 0 : moab::Range commonV=intersect (nodes1, nodes2);
311 [ # # ]: 0 : moab::Range v12;
312 [ # # ][ # # ]: 0 : v12.insert(v1); v12.insert(v2);
313 [ # # ]: 0 : moab::Range leftOver = subtract(commonV, v12);
314 [ # # ][ # # ]: 0 : for (moab::Range::iterator it = leftOver.begin(); it!=leftOver.end(); it++)
[ # # ][ # # ]
[ # # ]
315 : : {
316 [ # # ]: 0 : moab::EntityHandle thirdNode = *it;
317 [ # # ][ # # ]: 0 : if (!ehIsValid(thirdNode))
318 : 0 : continue;
319 : : // the moment we find a triple v1, v2, thirdNode that is not a triangle, return penalty
320 [ # # ]: 0 : moab::Range triple=v12;
321 [ # # ]: 0 : triple.insert(thirdNode);
322 [ # # ][ # # ]: 0 : moab::Range connTris;
323 [ # # ]: 0 : rval = mb->get_adjacencies(triple, 2, false, connTris );
324 [ # # ]: 0 : if (moab::MB_SUCCESS == rval)
325 : : {
326 [ # # ][ # # ]: 0 : if (connTris.empty())
327 : 0 : return MESH_INVERSION_PENALTY; // this means that there are no
328 : : /// triangles connected to the 3 nodes
329 : : // so it would close a hole/tunnel
330 : :
331 : 0 : int noValidTris = 0;
332 [ # # ][ # # ]: 0 : for (moab::Range::iterator it2 = connTris.begin(); it2!=connTris.end(); it2++)
[ # # ][ # # ]
[ # # ]
333 : : {
334 [ # # ][ # # ]: 0 : if (ehIsValid(*it2))
[ # # ]
335 : : {
336 : 0 : noValidTris++;
337 : 0 : break;
338 : : }
339 : : }
340 [ # # ]: 0 : if (0==noValidTris)
341 [ # # ]: 0 : return MESH_INVERSION_PENALTY; // this means that there are no
342 : : // triangles connected to the 3 nodes
343 : : // so it would close a hole/tunnel
344 : :
345 : : }
346 : :
347 : 0 : }
348 : 0 : return 0.; // no penalty
349 : : }
350 : : static
351 : 0 : double pair_mesh_penalty( /*Model& M, Vertex *v1, Vertex *v2,*/
352 : : moab::EntityHandle v1, moab::EntityHandle v2, Vec3& vnew) {
353 [ # # ]: 0 : std::vector<moab::EntityHandle> changed;
354 : :
355 : : // construct the list of faces influenced by the
356 : : // moving of vertices v1 and v2 into vnew
357 : : //M.contractionRegion ( v1, v2, changed );
358 [ # # ]: 0 : moab::ErrorCode rval = contractionRegion(mb, v1, v2, changed);
359 [ # # ]: 0 : if (rval != moab::MB_SUCCESS) {
360 [ # # ]: 0 : std::cout << "error in getting adjacency information vs: "
361 [ # # ][ # # ]: 0 : << mb->id_from_handle(v1) << " " << mb->id_from_handle(v2) << "\n";
[ # # ][ # # ]
[ # # ][ # # ]
362 : : }
363 : :
364 : : // double Nsum = 0;
365 : 0 : double Nmin = 0;
366 : :
367 [ # # ]: 0 : for (unsigned int i = 0; i < changed.size(); i++) {
368 : : //Face& F = *changed ( i );
369 [ # # ]: 0 : moab::EntityHandle F = changed[i];
370 [ # # ][ # # ]: 0 : Vec3 f1, f2, f3;
[ # # ]
371 : :
372 [ # # ]: 0 : int nmapped = predict_face(F, v1, v2, vnew, f1, f2, f3);
373 : : //
374 : : // Only consider non-degenerate faces
375 [ # # ]: 0 : if (nmapped < 2) {
376 [ # # ]: 0 : Plane Pnew(f1, f2, f3);
377 [ # # ]: 0 : Plane p = trianglePlane(mb, F);
378 [ # # ][ # # ]: 0 : double delta = Pnew.normal() * p.normal(); // Pnew.normal() * F.plane().normal();
[ # # ]
379 : :
380 [ # # ]: 0 : if (Nmin > delta)
381 : 0 : Nmin = delta;
382 : : }
383 : : }
384 : :
385 : : //return (-Nmin) * MESH_INVERSION_PENALTY;
386 [ # # ]: 0 : if (Nmin < 0.0)
387 : 0 : return MESH_INVERSION_PENALTY;
388 : : else
389 : 0 : return 0.0;
390 : : }
391 : :
392 : : static
393 : 60616 : void compute_pair_info(pair_info *pair /* Model * pM0,*/) {
394 : 60616 : moab::EntityHandle v0 = pair->v0;
395 : 60616 : moab::EntityHandle v1 = pair->v1;
396 : :
397 : : // Vertex *v0 = pair->v0;
398 : : // Vertex *v1 = pair->v1;
399 : :
400 [ + - ]: 60616 : vert_info& v0_info = vertex_info(v0);
401 [ + - ]: 60616 : vert_info& v1_info = vertex_info(v1);
402 : :
403 [ + - ]: 60616 : Mat4 Q = v0_info.Q + v1_info.Q;
404 : 60616 : double norm = v0_info.norm + v1_info.norm;
405 : :
406 [ + - ]: 60616 : pair->cost = quadrix_pair_target(Q, v0, v1, pair->candidate);
407 : :
408 [ - + ]: 60616 : if (opts.will_weight_by_area)
409 : 0 : pair->cost /= norm;
410 : :
411 [ - + ]: 60616 : if (opts.will_preserve_mesh_quality)
412 [ # # ]: 0 : pair->cost += pair_mesh_penalty(/* *pM0,*/v0, v1, pair->candidate);
413 : :
414 [ + - ]: 60616 : if (opts.height_fields)
415 [ + - ]: 60616 : pair->cost += pair_mesh_positivity(/* *pM0, */v0, v1, pair->candidate);
416 : :
417 [ - + ]: 60616 : if (opts.topology)
418 [ # # ]: 0 : pair->cost += pair_mesh_topology(v0, v1);
419 : :
420 [ - + ][ # # ]: 60616 : if (opts.logfile && opts.selected_output & OUTPUT_CONTRACTIONS)
421 : :
422 [ # # ][ # # ]: 0 : *opts.logfile << "pair ( v" << mb->id_from_handle(v0) << " v"
[ # # ][ # # ]
423 [ # # ][ # # ]: 0 : << mb->id_from_handle(v1) << " ) cost: " << -pair->cost << std::endl;
[ # # ][ # # ]
[ # # ]
424 : : //
425 : : // NOTICE: In the heap we use the negative cost. That's because
426 : : // the heap is implemented as a MAX heap.
427 : : //
428 [ + - ][ + + ]: 60616 : if (pair->isInHeap()) {
429 [ + - ]: 45416 : heap->update(pair, -pair->cost);
430 : : } else {
431 [ + - ]: 15200 : heap->insert(pair, -pair->cost);
432 : : }
433 [ - + ][ # # ]: 60616 : if (opts.logfile && opts.selected_output & OUTPUT_COST) {
434 [ # # ]: 0 : heap_node *top = heap->top();
435 : 0 : pair_info *pairTop = (pair_info *) top->obj;
436 [ # # ][ # # ]: 0 : *opts.logfile << " i/u pair (" << uniqID(pair->v0) << "," << uniqID(
[ # # ][ # # ]
437 [ # # ][ # # ]: 0 : pair->v1) << ")=" << pair->cost << " min : (" << uniqID(pairTop->v0)
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
438 [ # # ][ # # ]: 0 : << "," << uniqID(pairTop->v1) << ") " << pairTop->cost << std::endl;
[ # # ][ # # ]
[ # # ][ # # ]
439 : : }
440 : 60616 : }
441 : 544 : void recomputeChangedPairsCost(std::vector<moab::EntityHandle> & changed, moab::EntityHandle v0) {
442 : : //
443 [ + + ]: 5032 : for (unsigned int i = 0; i < changed.size(); i++) {
444 : :
445 [ + - ]: 4488 : moab::EntityHandle F = changed[i];
446 : : const moab::EntityHandle * conn;
447 : : int num_nodes;
448 [ + - ]: 4488 : mb->get_connectivity(F, conn, num_nodes);
449 : : //if (!F->isValid())
450 : : // continue;
451 : : // recompute the pair that is not connected to vertex v0
452 : : // loop over all the vertices of F that are not v0, and recompute the costs
453 : : // of all the pairs associated that do not contain v0
454 : : // we do not have to recreate or delete any pair, we just recompute what we have
455 : : // some will be recomputed 2 times, but it is OK
456 [ + + ]: 17952 : for (int k = 0; k < 3; k++) {
457 : 13464 : moab::EntityHandle v = conn[k];
458 [ + + ]: 13464 : if (v == v0)
459 : 4488 : continue;
460 [ + - ]: 8976 : vert_info & v_info = vertex_info(v);
461 [ + - ][ + + ]: 58558 : for (int j = 0; j < v_info.pairs.length(); j++) {
462 [ + - ]: 49582 : pair_info *p = v_info.pairs(j);
463 [ + + ][ + + ]: 49582 : if (p->v0 == v0 || p->v1 == v0)
464 : 8976 : continue; // do not recompute cost of pairs already computed
465 [ - + ][ # # ]: 40606 : if (opts.logfile && (opts.selected_output & OUTPUT_COST))
466 [ # # ][ # # ]: 0 : *opts.logfile << "recompute cost of pair (v" << uniqID(p->v0) + 1
[ # # ]
467 [ # # ][ # # ]: 0 : << " v" << uniqID(p->v1) + 1 << ")" << std::endl;
[ # # ][ # # ]
[ # # ]
468 [ + - ]: 40606 : compute_pair_info(p);
469 : : }
470 : :
471 : : }
472 : :
473 : : }
474 : 544 : }
475 : :
476 : : static
477 : 544 : void do_contract(pair_info *pair) {
478 : :
479 : 544 : moab::EntityHandle v0 = pair->v0;
480 : 544 : moab::EntityHandle v1 = pair->v1;
481 : : // cost of contraction is accumulated at v0
482 : 544 : double costToContract = pair->cost;
483 [ + - ]: 544 : vert_info& v0_info = vertex_info(v0);
484 [ + - ]: 544 : vert_info& v1_info = vertex_info(v1);
485 [ + - ]: 544 : Vec3 vnew = pair->candidate;
486 [ - + ][ # # ]: 544 : if (opts.logfile && (opts.selected_output & OUTPUT_COST)) {
487 [ # # ][ # # ]: 0 : *opts.logfile << "---- do contract v0:" << uniqID(v0) + 1 << " v1:"
[ # # ][ # # ]
488 [ # # ][ # # ]: 0 : << uniqID(v1) + 1 << std::endl;
[ # # ]
489 : : }
490 : : int i;
491 : :
492 : : //
493 : : // Make v0 be the new vertex
494 [ + - ]: 544 : v0_info.Q += v1_info.Q;
495 : 544 : v0_info.norm += v1_info.norm;
496 : :
497 : : //
498 : : // Perform the actual contraction
499 [ + - ]: 544 : std::vector<moab::EntityHandle> changed;
500 [ + - ]: 544 : moab::ErrorCode rval1 = contract(mb, v0, v1, vnew, changed);
501 : : (void) rval1;
502 : : // this is a list of triangles still connected to v0 (are they valid? probably)
503 [ - + ]: 544 : if (opts.will_preserve_mesh_quality) {
504 : : // recompute normals only in this case, because they are not needed otherwise
505 : 0 : int size = changed.size();
506 [ # # ]: 0 : for (int i = 0; i < size; i++) {
507 [ # # ][ # # ]: 0 : computeTrianglePlane(mb, changed[i]);
508 : : }
509 : : }
510 [ - + ][ # # ]: 544 : assert (moab::MB_SUCCESS == rval1 || moab::MB_MULTIPLE_ENTITIES_FOUND == rval1);
511 : :
512 : : #ifdef SUPPORT_VCOLOR
513 : : //
514 : : // If the vertices are colored, color the new vertex
515 : : // using the average of the old colors.
516 : : v0->props->color += v1->props->color;
517 : : v0->props->color /= 2;
518 : : #endif
519 : :
520 : : //
521 : : // Remove the pair that we just contracted
522 [ + - ]: 544 : delete_pair(pair);
523 : : //
524 : : // Recalculate pairs associated with v0
525 [ + - ][ + + ]: 3952 : for (i = 0; i < v0_info.pairs.length(); i++) {
526 [ + - ]: 3408 : pair_info *p = v0_info.pairs(i);
527 [ + - ]: 3408 : compute_pair_info(p);
528 : : }
529 : :
530 : : //
531 : : // Process pairs associated with now dead vertex
532 : :
533 [ + + ][ + - ]: 544 : static pair_buffer condemned(6); // collect condemned pairs for execution
[ + - ][ # # ]
534 [ + - ]: 544 : condemned.reset();
535 : :
536 [ + - ][ + + ]: 2948 : for (i = 0; i < v1_info.pairs.length(); i++) {
537 [ + - ]: 2404 : pair_info *p = v1_info.pairs(i);
538 : :
539 : 2404 : moab::EntityHandle u = 0L;
540 [ + + ]: 2404 : if (p->v0 == v1)
541 : 1390 : u = p->v1;
542 [ + - ]: 1014 : else if (p->v1 == v1)
543 : 1014 : u = p->v0;
544 : : else
545 [ # # ][ # # ]: 0 : std::cerr << "YOW! This is a bogus pair." << std::endl;
546 : :
547 [ + - ][ + + ]: 2404 : if (!check_for_pair(u, v0)) {
548 : 1402 : p->v0 = v0;
549 : 1402 : p->v1 = u;
550 [ + - ]: 1402 : v0_info.pairs.add(p);
551 [ + - ]: 1402 : compute_pair_info(p);
552 : : } else
553 [ + - ]: 1002 : condemned.add(p);
554 : : }
555 : :
556 [ + - ][ + + ]: 1546 : for (i = 0; i < condemned.length(); i++)
557 : : // Do you have any last requests?
558 [ + - ][ + - ]: 1002 : delete_pair(condemned(i));
559 : : // only now you can delete the vertex v1 from database
560 : : // moab::ErrorCode rval = mb->delete_entities(&v1, 1);
561 : : // no, it is better to invalidate the vertex, do not delete it
562 : : // maybe we will delete at the end all that are invalid ??
563 : 544 : int invalid = 0;
564 [ + - ]: 544 : moab::ErrorCode rval = mb->tag_set_data(validTag, &v1, 1, &invalid);
565 : : (void) rval;
566 [ - + ]: 544 : assert (moab::MB_SUCCESS == rval);
567 : :
568 [ - + ]: 544 : if (opts.plotCost) {
569 : 0 : double cost_at_v0 = 0; // maybe it is already set before
570 [ # # ]: 0 : rval = mb->tag_get_data(costTag, &v0, 1, &cost_at_v0);
571 : 0 : cost_at_v0 += costToContract;
572 [ # # ]: 0 : rval = mb->tag_set_data(costTag, &v0, 1, &cost_at_v0);
573 : : }
574 : :
575 [ + - ]: 544 : v1_info.pairs.reset(); // safety precaution
576 [ + - ]: 544 : recomputeChangedPairsCost(changed, v0);
577 : :
578 : 544 : }
579 : :
580 : : ////////////////////////////////////////////////////////////////////////
581 : : //
582 : : // External interface: setup and single step iteration
583 : : //
584 : :
585 : 0 : bool decimate_quadric(moab::EntityHandle v, Mat4& Q) {
586 [ # # ]: 0 : if (vinfo.length() > 0) {
587 : 0 : vert_info & vinf = vertex_info(v);
588 : 0 : Q = vinf.Q;
589 : 0 : return true;
590 : : } else
591 : 0 : return false;
592 : : }
593 : :
594 : : // it is assumed it is mb, moab interface
595 : 544 : void decimate_contract() {
596 : : heap_node *top;
597 : : pair_info *pair;
598 : :
599 : : for (;;) {
600 : 544 : top = heap->extract();
601 [ - + ]: 544 : if (!top)
602 : 0 : return;
603 : 544 : pair = (pair_info *) top->obj;
604 : :
605 : : //
606 : : // This may or may not be necessary. I'm just not quite
607 : : // willing to assume that all the junk has been removed from the
608 : : // heap.
609 [ + - ]: 544 : if (pair->isValid())
610 : 544 : break;
611 : :
612 : 0 : delete_pair(pair);
613 : : }
614 : :
615 [ - + ][ # # ]: 544 : if (opts.logfile && (opts.selected_output & OUTPUT_COST))
616 : 0 : *opts.logfile << "#$cost " << validFaceCount << " before contract: "
617 : 0 : << pair->cost << std::endl;
618 : :
619 : 544 : do_contract(pair);
620 : :
621 [ - + ][ # # ]: 544 : if (opts.logfile && (opts.selected_output & OUTPUT_COST))
622 : 0 : *opts.logfile << "#$cost " << validFaceCount << std::endl;
623 : :
624 : 544 : validVertCount--; // Attempt to maintain valid vertex information
625 : : }
626 : :
627 : 0 : double decimate_error(moab::EntityHandle v) {
628 : 0 : vert_info& info = vertex_info(v);
629 : :
630 [ # # ]: 0 : double err = quadrix_evaluate_vertex(v, info.Q);
631 : :
632 [ # # ]: 0 : if (opts.will_weight_by_area)
633 : 0 : err /= info.norm;
634 : :
635 : 0 : return err;
636 : : }
637 : :
638 : 544 : double decimate_min_error() {
639 : : heap_node *top;
640 : : pair_info *pair;
641 : :
642 : : for (;;) {
643 : 544 : top = heap->top();
644 [ - + ]: 544 : if (!top)
645 : 0 : return -1.0;
646 : 544 : pair = (pair_info *) top->obj;
647 : :
648 [ + - ]: 544 : if (pair->isValid())
649 : 544 : break;
650 : :
651 : 0 : top = heap->extract();
652 : 0 : delete_pair(pair);
653 : : }
654 : :
655 : 544 : return pair->cost;
656 : : }
657 : : #if 0
658 : : double decimate_max_error ( )
659 : : {
660 : : real max_err = 0;
661 : :
662 : : for ( int i=0; i<m.vertCount(); i++ )
663 : : if ( m.vertex ( i )->isValid() )
664 : : {
665 : : max_err = MAX ( max_err, decimate_error ( m.vertex ( i ) ) );
666 : : }
667 : :
668 : : return max_err;
669 : : }
670 : : #endif
671 : : namespace MeshKit {
672 : :
673 : 2 : QslimDecimation::QslimDecimation(moab::Interface * mbi,
674 : 2 : moab::EntityHandle root_set) {
675 : 2 : _mb = mbi;
676 : 2 : iniSet = root_set;// it is not necessarily the root set; this is external; bad design
677 : 2 : }
678 : :
679 : 0 : QslimDecimation::~QslimDecimation() {
680 [ # # ]: 0 : }
681 : 2 : int QslimDecimation::decimate(QslimOptions & iOpts, moab::Range & oRange) {
682 : : // opts is external
683 : 2 : opts = iOpts;
684 : :
685 : 2 : mb = _mb; // (reinterpret_cast<MBiMesh *> (m_mesh))->mbImpl;
686 : : // need to get all the triangles from the set
687 : : // also all the edges, and all vertices
688 : : // not a good design here; mb is extern !
689 : : //
690 [ - + ]: 2 : if (NULL == mb)
691 : 0 : return 1;// error
692 : : //moab::EntityHandle mbSet = reinterpret_cast<moab::EntityHandle>(m_InitialSet);
693 : :
694 : 2 : clock_t start_time = clock();
695 : 2 : int err = Init();
696 [ - + ]: 2 : if (err) {
697 : 0 : std::cerr << "Error in initialization of decimation. Abort\n";
698 : 0 : return 1;
699 : : }
700 : 2 : clock_t init_time = clock();
701 : 2 : std::cout << " Initialization " << (double) (init_time - start_time)
702 : 2 : / CLOCKS_PER_SEC << " s.\n";
703 : :
704 : 2 : int faces_reduction = validFaceCount - opts.face_target;
705 : 2 : int counter = 0, interval = 0;
706 : 2 : clock_t currentTime = init_time;
707 [ + + ][ + - ]: 546 : while (validFaceCount > opts.face_target && decimate_min_error()
[ + + ]
708 : 544 : < opts.error_tolerance) {
709 : 544 : int initf = validFaceCount;
710 : : // big routine
711 : 544 : decimate_contract();
712 : 544 : counter += (initf - validFaceCount);
713 [ + + ]: 544 : if (counter > faces_reduction / opts.timingIntervals) {
714 : : // print some time stats
715 : 18 : clock_t p10_time = clock();
716 : 18 : std::cout << " " << ++interval << "/" << opts.timingIntervals
717 : 18 : << " reduce to " << validFaceCount << " faces in "
718 : 18 : << (double) (p10_time - currentTime) / CLOCKS_PER_SEC << " s, total:"
719 : 18 : << (double) (p10_time - init_time) / CLOCKS_PER_SEC << " s.\n";
720 : 18 : counter = 0;
721 : 18 : currentTime = p10_time;
722 : : }
723 : : }
724 : :
725 : 2 : clock_t finish_time = clock();
726 : 2 : std::cout << " Decimation: " << (double) (finish_time - init_time)
727 : 2 : / CLOCKS_PER_SEC << " s.\n";
728 : :
729 [ - + ]: 2 : if (opts.create_range) {
730 : : // the remaining nodes and triangles are copied in the range
731 : : // they are put in the old set, too
732 : : // maybe this has to change
733 : : // count first valid vertices and triangles
734 [ # # ]: 0 : moab::Range::iterator it;
735 [ # # ]: 0 : std::vector<moab::EntityHandle> validVerts;
736 [ # # ]: 0 : std::vector<moab::EntityHandle> validTris;
737 [ # # ][ # # ]: 0 : for (it = triangles.begin(); it != triangles.end(); it++) {
[ # # ][ # # ]
[ # # ]
738 [ # # ][ # # ]: 0 : if (ehIsValid(*it))
[ # # ]
739 [ # # ][ # # ]: 0 : validTris.push_back(*it);
740 : : }
741 [ # # ][ # # ]: 0 : for (it = verts.begin(); it != verts.end(); it++) {
[ # # ][ # # ]
[ # # ]
742 [ # # ][ # # ]: 0 : if (ehIsValid(*it))
[ # # ]
743 [ # # ][ # # ]: 0 : validVerts.push_back(*it);
744 : : }
745 : : //
746 [ # # ]: 0 : std::vector<double> coords;
747 : 0 : int numNodes = (int) validVerts.size();
748 : 0 : int numTriangles = (int) validTris.size();
749 [ # # ]: 0 : coords.resize(3 * numNodes);
750 : :
751 [ # # ]: 0 : moab::ErrorCode rval = mb->get_coords(&(validVerts[0]), numNodes,
752 [ # # ][ # # ]: 0 : &(coords[0]));
753 : : (void) rval;
754 [ # # ]: 0 : assert(moab::MB_SUCCESS==rval);
755 : :
756 : : // create the new vertices, at the same coordinates
757 : : // put those verts in the range that is output
758 : :
759 : : // to make sure, the range is cleared
760 [ # # ]: 0 : oRange.clear();
761 [ # # ][ # # ]: 0 : rval = mb->create_vertices(&coords[0], numNodes, oRange);
762 [ # # ]: 0 : assert(moab::MB_SUCCESS==rval);
763 : : // the first element in the range must be the start of the new vertex sequence
764 [ # # ]: 0 : std::map<moab::EntityHandle, moab::EntityHandle> mapping; // this will be from old
765 : : // to new vertices, for connectivity
766 [ # # ]: 0 : for (int i = 0; i < numNodes; i++) {
767 [ # # ][ # # ]: 0 : mapping[validVerts[i]] = oRange[i];
[ # # ]
768 : : }
769 : :
770 : : //get the query interface, which we will use to create the triangles directly
771 : : moab::ReadUtilIface *iface;
772 [ # # ]: 0 : rval = mb -> query_interface(iface);// use the new query interface
773 [ # # ]: 0 : assert(moab::MB_SUCCESS==rval);
774 : :
775 : : //create the triangles, get a direct ptr to connectivity
776 : : moab::EntityHandle starth, *connect;
777 : : rval = iface -> get_element_connect(numTriangles, 3, moab::MBTRI, 1,
778 [ # # ]: 0 : starth, connect);
779 [ # # ]: 0 : assert(moab::MB_SUCCESS==rval);
780 : : // get first the connectivity of the old triangles
781 : : const moab::EntityHandle * conn3;
782 [ # # ]: 0 : for (int i = 0; i < numTriangles; i++) {
783 : : int num_nodes;
784 : : // get each valid triangle one by one, and set the new connectivity directly
785 [ # # ][ # # ]: 0 : rval = mb->get_connectivity(validTris[i], conn3, num_nodes);
786 [ # # ][ # # ]: 0 : assert( (moab::MB_SUCCESS==rval) && (num_nodes==3));
787 : : // directly modify the connect array in database
788 [ # # ]: 0 : for (int j = 0; j < 3; j++)
789 [ # # ]: 0 : connect[j] = mapping[conn3[j]];
790 : :
791 : : // update adjacencies
792 : : // not very smart...; we would like to update once and for all triangles
793 : : // not in a loop
794 [ # # ]: 0 : rval = iface -> update_adjacencies(starth+i, 1, 3, connect);
795 [ # # ]: 0 : assert(moab::MB_SUCCESS==rval);
796 : :
797 : 0 : connect += 3; // advance
798 : : }
799 : :
800 : :
801 : :
802 : : // clear completely the initial set, after deleting all elements from it...
803 : : // ok, we are done, commit to ME ?
804 [ # # ]: 0 : rval = mb->delete_entities(triangles);
805 [ # # ]: 0 : assert(moab::MB_SUCCESS==rval);
806 [ # # ]: 0 : rval = mb->delete_entities(edgs);
807 [ # # ]: 0 : assert(moab::MB_SUCCESS==rval);
808 [ # # ]: 0 : rval = mb->delete_entities(verts);
809 [ # # ]: 0 : assert(moab::MB_SUCCESS==rval);
810 : : // remove everything from the initial set, because we will add the new triangles
811 [ # # ]: 0 : mb->remove_entities(iniSet, triangles);
812 [ # # ]: 0 : mb->remove_entities(iniSet, verts);
813 [ # # ]: 0 : mb->remove_entities(iniSet, edgs);
814 : :
815 : : //add triangles to output range (for the MESelection)
816 [ # # ]: 0 : oRange.insert(starth, starth + numTriangles - 1);
817 : : // add all entities from the range to the initial set, now
818 [ # # ]: 0 : rval = mb->add_entities(iniSet, oRange);
819 [ # # ]: 0 : assert(moab::MB_SUCCESS==rval);
820 : : //me->commit_mesh(mit->second, COMPLETE_MESH);
821 : : // end copy
822 : : } else {
823 [ + - ]: 2 : moab::Range::const_reverse_iterator rit;
824 [ + - ]: 2 : if (opts.useDelayedDeletion) {
825 : :
826 : : // put in a range triangles to delete
827 [ + - ]: 2 : moab::Range delRange;
828 : : // delete triangles and edges that are invalid
829 [ + - ][ + - ]: 10002 : for (rit = triangles.rbegin(); rit != triangles.rend(); rit++) {
[ + - ][ + - ]
[ + + ]
830 [ + - ]: 10000 : moab::EntityHandle tr = *rit;
831 : : // check the validity
832 [ + - ][ + + ]: 10000 : if (ehIsValid(tr))
833 : 8998 : continue;
834 [ + - ]: 1002 : mb->delete_entities(&tr, 1);
835 [ + - ]: 1002 : delRange.insert(tr);
836 : : }
837 [ + - ]: 2 : mb->remove_entities(iniSet, delRange);
838 : : // maybe we should delete all edges, but for now, we will keep them
839 [ + - ][ + - ]: 15202 : for (rit = edgs.rbegin(); rit != edgs.rend(); rit++) {
[ + - ][ + - ]
[ + + ]
840 [ + - ]: 15200 : moab::EntityHandle v = *rit;
841 : : // check the validity
842 [ + - ][ + + ]: 15200 : if (ehIsValid(v))
843 : 13654 : continue;
844 [ + - ]: 1546 : mb->delete_entities(&v, 1);
845 : 2 : }
846 : :
847 : : }
848 : : // delete them one by one
849 [ + - ][ + - ]: 5204 : for (rit = verts.rbegin(); rit != verts.rend(); rit++) {
[ + - ][ + - ]
[ + + ]
850 [ + - ]: 5202 : moab::EntityHandle v = *rit;
851 : : // check the validity
852 [ + - ][ + + ]: 5202 : if (ehIsValid(v))
853 : 4658 : continue;
854 [ + - ]: 544 : mb->delete_entities(&v, 1);
855 : : }
856 : : }
857 : 2 : clock_t delete_vTime = clock();
858 : 2 : std::cout << " Delete Vertices: " << (double) (delete_vTime - finish_time)
859 : 2 : / CLOCKS_PER_SEC << " s.\n";
860 : : // we need to delete the tags we created; they are artificial
861 : : // list of tags to delete:
862 : : // moab::Tag uniqIDtag; // this will be used to mark vertices moab::EntityHandles
863 : : // moab::Tag validTag;
864 : : // moab::Tag planeDataTag;
865 : :
866 : : // moab::Tag costTag; // simplification induces an error cost at each vertex
867 : : // try to keep adding the cost, to see if it is spreading nicely
868 : :
869 : : // keep only the cost, the rest are artificial
870 : 2 : mb->tag_delete(uniqIDtag);
871 : 2 : mb->tag_delete(validTag);
872 : 2 : mb->tag_delete(planeDataTag);
873 : : //
874 : :
875 : 2 : return 0;
876 : : }
877 : :
878 : 2 : int QslimDecimation::Init() {
879 : : int i;
880 : : unsigned int j;
881 : :
882 : : //moab::EntityHandle * set = reinterpret_cast<moab::EntityHandle *> (&_InitialSet);
883 : : moab::ErrorCode rval = mb->get_entities_by_type(iniSet, moab::MBTRI,
884 [ + - ]: 2 : triangles);
885 [ + - ]: 2 : validFaceCount = triangles.size();// this gets reduced every time we simplify the model
886 : :
887 : : // store the normals/planes computed at each triangle
888 : : // we may need just the normals, but compute planes, it is about the same job
889 : 2 : double defPlane[] = { 0., 0., 1., 0. };
890 : : rval = mb->tag_get_handle("PlaneTriangleData", 4, moab::MB_TYPE_DOUBLE,
891 [ + - ]: 2 : planeDataTag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT, &defPlane);
892 : :
893 : : // compute the triangle plane and store it, for each triangle
894 [ + - ][ + - ]: 10002 : for (moab::Range::iterator itr = triangles.begin(); itr != triangles.end(); itr++) {
[ + - ][ + - ]
[ + + ]
895 : : // this index i will be the index in the vinfo array
896 [ + - ]: 10000 : moab::EntityHandle tri = *itr;
897 [ + - ]: 10000 : computeTrianglePlane(mb, tri);
898 : : // setting the data for the tag/triangle is done in the compute
899 : : //rval = mb->tag_set_data(planeDataTag, &tri, 1, &plane);
900 : : }
901 : :
902 : : // create all the edges if not existing
903 [ + - ]: 2 : mb->get_adjacencies(triangles, 1, true, edgs, moab::Interface::UNION);
904 : :
905 : : // moab::Range verts;// the vertices are always there, they do not need to be created
906 [ + - ]: 2 : mb->get_adjacencies(triangles, 0, true, verts, moab::Interface::UNION);
907 [ + - ]: 2 : int numNodes = verts.size();
908 : 2 : validVertCount = numNodes; // this will be kept
909 [ + - ]: 2 : vinfo.init(numNodes);
910 : : // set a unique integer tag with the position in vinfo array
911 : : // this will be used instead of v->uniqID in the vinfo array
912 : 2 : int def_data = -1;
913 : :
914 : : rval = mb->tag_get_handle("uniqID", 1, moab::MB_TYPE_INTEGER, uniqIDtag,
915 [ + - ]: 2 : moab::MB_TAG_DENSE | moab::MB_TAG_EXCL, &def_data);
916 [ - + ]: 2 : if (moab::MB_SUCCESS != rval)
917 : 0 : return 1;
918 : :
919 : 2 : unsigned char def_data_bit = 1;// valid by default
920 : :
921 : : rval = mb->tag_get_handle("valid", 1, moab::MB_TYPE_BIT, validTag,
922 [ + - ]: 2 : moab::MB_TAG_EXCL | moab::MB_TAG_BIT, &def_data_bit);
923 [ - + ]: 2 : if (moab::MB_SUCCESS != rval)
924 : 0 : return 1;
925 : :
926 [ - + ]: 2 : if (opts.plotCost) {
927 : 0 : double cost_default = 0.;
928 : :
929 : : rval = mb->tag_get_handle("costTAG", 1, moab::MB_TYPE_DOUBLE, costTag,
930 [ # # ]: 0 : moab::MB_TAG_DENSE | moab::MB_TAG_CREAT, &cost_default);
931 [ # # ]: 0 : if (moab::MB_SUCCESS != rval)
932 : 0 : return 1;
933 : : }
934 : :
935 : : // set tag for each vertex; this will not be changed during simplification
936 : 2 : i = 0; // for index
937 [ + - ][ + - ]: 5204 : for (moab::Range::iterator it = verts.begin(); it != verts.end(); it++, i++) {
[ + - ][ + - ]
[ + + ]
938 : : // this index i will be the index in the vinfo array
939 [ + - ]: 5202 : moab::EntityHandle v = *it;
940 [ + - ]: 5202 : rval = mb->tag_set_data(uniqIDtag, &v, 1, &i);
941 : : // the default valid data is 1; set it to 0 only to "mark" the vertex invalid for future deletion
942 : : // is it really necessary if we put the default value as 1 anyway
943 : :
944 : : //rval = mb->tag_set_data(validTag, &v, 1, &def_data_bit);
945 : : }
946 [ + - ]: 2 : moab::Range::iterator it;
947 [ - + ][ # # ]: 2 : if (opts.logfile && opts.selected_output & OUTPUT_MODEL_DEFN) {
948 [ # # ][ # # ]: 0 : for (it = verts.begin(); it != verts.end(); it++) {
[ # # ][ # # ]
[ # # ]
949 [ # # ]: 0 : moab::EntityHandle v = *it;
950 : : double coords[3];
951 [ # # ]: 0 : rval = mb->get_coords(&v, 1, coords);
952 [ # # ][ # # ]: 0 : *opts.logfile << "v: " << uniqID(v) << " " << mb->id_from_handle(v)
[ # # ][ # # ]
[ # # ][ # # ]
953 [ # # ][ # # ]: 0 : << " " << coords[0] << " " << coords[1] << " " << coords[2]
[ # # ][ # # ]
[ # # ][ # # ]
954 [ # # ]: 0 : << std::endl;
955 : : }
956 : : }
957 [ + - ][ + - ]: 2 : std::cout << " Decimate: Distributing shape constraints." << std::endl;
958 : :
959 [ - + ]: 2 : if (opts.will_use_vertex_constraint)
960 [ # # ][ # # ]: 0 : for (it = verts.begin(); it != verts.end(); it++) {
[ # # ][ # # ]
[ # # ]
961 [ # # ]: 0 : moab::EntityHandle v = *it;
962 [ # # ][ # # ]: 0 : vertex_info(v).Q = quadrix_vertex_constraint(v);
[ # # ]
963 : : }
964 : :
965 [ + - ]: 2 : if (opts.will_use_plane_constraint) {
966 [ + - ][ + - ]: 10002 : for (it = triangles.begin(); it != triangles.end(); it++) {
[ + - ][ + - ]
[ + + ]
967 : :
968 [ + - ]: 10000 : moab::EntityHandle tr = *it;
969 [ + - ]: 10000 : Mat4 Q = quadrix_plane_constraint(tr);
970 : 10000 : double norm = 0.0;
971 : :
972 [ - + ]: 10000 : if (opts.will_weight_by_area) {
973 : 0 : norm = 1; // triangle area : m_model->face ( i )->area();
974 [ # # ]: 0 : Q *= norm;
975 : : }
976 : : const moab::EntityHandle * conn;
977 : : int num_nodes;
978 [ + - ]: 10000 : rval = mb->get_connectivity(tr, conn, num_nodes);
979 [ + + ]: 40000 : for (j = 0; j < 3; j++) {
980 [ + - ]: 30000 : vert_info& vj_info = vertex_info(conn[j]);
981 [ + - ]: 30000 : vj_info.Q += Q;
982 [ + - ]: 30000 : vertex_info(conn[j]).norm += norm;
983 : :
984 : : }
985 : : }
986 : : }
987 : : // just define (one uniqTag for a triangle, see what is happening)
988 [ + - ]: 2 : moab::EntityHandle tr1 = triangles[0];
989 [ + - ]: 2 : rval = mb->tag_set_data(uniqIDtag, &tr1, 1, &i);// just some value
990 : :
991 [ + - ]: 2 : if (opts.will_constrain_boundaries) {
992 [ + - ]: 2 : std::cout << " Decimate: Accumulating discontinuity constraints."
993 [ + - ]: 2 : << std::endl;
994 [ + - ][ + - ]: 15202 : for (it = edgs.begin(); it != edgs.end(); it++) {
[ + - ][ + - ]
[ + + ]
995 [ + - ]: 15200 : moab::EntityHandle edg = *it;
996 [ + - ][ + + ]: 15200 : if (is_border(edg)) {
997 : : const moab::EntityHandle * conn;
998 : : int num_nodes;
999 [ + - ]: 400 : rval = mb->get_connectivity(edg, conn, num_nodes);
1000 [ - + ]: 400 : if (moab::MB_SUCCESS != rval)
1001 : 0 : return 1;// fail
1002 [ + - ]: 400 : Mat4 B = quadrix_discontinuity_constraint(edg);
1003 : 400 : double norm = 0.0;
1004 : :
1005 [ - + ]: 400 : if (opts.will_weight_by_area) {
1006 [ # # ]: 0 : Vec3 ve1 = getVec3FromMBVertex(mb, conn[0]);
1007 [ # # ]: 0 : Vec3 ve2 = getVec3FromMBVertex(mb, conn[1]);
1008 [ # # ][ # # ]: 0 : norm = norm2(ve1 - ve2);
1009 [ # # ]: 0 : B *= norm;
1010 : : }
1011 : :
1012 [ + - ]: 400 : B *= opts.boundary_constraint_weight;
1013 [ + - ]: 400 : vert_info& v0_info = vertex_info(conn[0]);
1014 [ + - ]: 400 : vert_info& v1_info = vertex_info(conn[1]);
1015 [ + - ]: 400 : v0_info.Q += B;
1016 : 400 : v0_info.norm += norm;
1017 [ + - ]: 400 : v1_info.Q += B;
1018 : 400 : v1_info.norm += norm;
1019 : : }
1020 : : }
1021 : : }
1022 : :
1023 [ + - ][ + - ]: 2 : std::cout << " Decimate: Allocating heap." << std::endl;
1024 [ + - ][ + - ]: 2 : heap = new Heap(edgs.size());
[ + - ]
1025 : :
1026 : 2 : int pair_count = 0;
1027 : :
1028 [ + - ][ + - ]: 2 : std::cout << " Decimate: Collecting pairs [edges]." << std::endl;
1029 [ + - ][ + - ]: 15202 : for (it = edgs.begin(); it != edgs.end(); it++) {
[ + - ][ + - ]
[ + + ]
1030 [ + - ]: 15200 : moab::EntityHandle edg = *it;
1031 : : const moab::EntityHandle * conn;
1032 : : int num_nodes;
1033 [ + - ]: 15200 : rval = mb->get_connectivity(edg, conn, num_nodes);
1034 [ - + ]: 15200 : if (moab::MB_SUCCESS != rval)
1035 : 0 : return 1;// fail
1036 [ + - ]: 15200 : pair_info *pair = new_pair(conn[0], conn[1]);
1037 [ + - ]: 15200 : compute_pair_info(pair);
1038 : 15200 : pair_count++;
1039 : : }
1040 : :
1041 [ - + ]: 2 : if (opts.pair_selection_tolerance < 0) {
1042 : 0 : opts.pair_selection_tolerance = 1;//m_model->bounds.radius * 0.05;
1043 [ # # ]: 0 : std::cout << " Decimate: Auto-limiting at 5% of model radius."
1044 [ # # ]: 0 : << std::endl;
1045 : : }
1046 : : proximity_limit = opts.pair_selection_tolerance
1047 : 2 : * opts.pair_selection_tolerance;
1048 [ - + ]: 2 : if (proximity_limit > 0) {
1049 [ # # ]: 0 : std::cout << " Decimate: Collecting pairs [limit="
1050 [ # # ][ # # ]: 0 : << opts.pair_selection_tolerance << "]." << std::endl;
[ # # ]
1051 : : // use adaptive kd tree to find proximity vertices
1052 : 0 : moab::EntityHandle tree_root = 0;
1053 [ # # ]: 0 : moab::AdaptiveKDTree kd(mb);
1054 [ # # ]: 0 : rval = kd.build_tree(verts, &tree_root);
1055 [ # # ]: 0 : if (rval != moab::MB_SUCCESS) {
1056 [ # # ][ # # ]: 0 : std::cout << "Can't build tree for vertices" << std::endl;
1057 : 0 : return 1;
1058 : : }
1059 : :
1060 [ # # ][ # # ]: 0 : for (it = verts.begin(); it != verts.end(); it++) {
[ # # ][ # # ]
[ # # ][ # # ]
1061 [ # # ]: 0 : moab::Range closeVertices;
1062 [ # # ]: 0 : closeVertices.clear();
1063 [ # # ]: 0 : moab::EntityHandle v = *it;
1064 : : double coords[3];
1065 [ # # ]: 0 : mb->get_coords(&v, 1, coords);
1066 : : //moab::CartVect v1(coords);
1067 [ # # ]: 0 : std::vector<moab::EntityHandle> leaves; // sets of vertices close by
1068 : : kd.distance_search( coords,
1069 [ # # ]: 0 : opts.pair_selection_tolerance, leaves);
1070 : : // add to the list of close vertices
1071 [ # # ]: 0 : for (j = 0; j < leaves.size(); j++) {
1072 [ # # ]: 0 : rval = mb->get_entities_by_type(leaves[j], moab::MBVERTEX,
1073 [ # # ]: 0 : closeVertices);// add to the list
1074 : : }
1075 : :
1076 [ # # ][ # # ]: 0 : for (moab::Range::iterator it2 = closeVertices.begin(); it2
[ # # ][ # # ]
1077 [ # # ]: 0 : != closeVertices.end(); it2++) {
1078 [ # # ]: 0 : moab::EntityHandle vclose = *it2;
1079 [ # # ]: 0 : if (v == vclose)
1080 : 0 : continue;
1081 : : double coords2[3];
1082 [ # # ]: 0 : mb->get_coords(&vclose, 1, coords2);
1083 : :
1084 : : //moab::CartVect v2(coords2);
1085 : 0 : double dd = (coords[0] - coords2[0]) * (coords[0] - coords2[0])
1086 : 0 : + (coords[1] - coords2[1]) * (coords[1] - coords2[1]) + (coords[2]
1087 : 0 : - coords2[2]) * (coords[2] - coords2[2]);
1088 : :
1089 [ # # ]: 0 : if (dd > proximity_limit)
1090 : 0 : continue;
1091 : :
1092 : : /*
1093 : : #ifdef SAFETY
1094 : : assert ( pair_is_valid ( v1,v2 ) );
1095 : : #endif
1096 : : */
1097 [ # # ][ # # ]: 0 : if (!check_for_pair(v, vclose)) {
1098 [ # # ]: 0 : pair_info *pair = new_pair(v, vclose);
1099 [ # # ]: 0 : compute_pair_info(pair);
1100 : 0 : pair_count++;
1101 : : }
1102 : : }
1103 : :
1104 : 0 : }
1105 : : } else
1106 [ + - ][ + - ]: 2 : std::cout << " Decimate: Ignoring non-edge pairs [limit=0]." << std::endl;
1107 : :
1108 [ + - ][ + - ]: 2 : std::cout << " Decimate: Designated " << pair_count << " pairs."
[ + - ]
1109 [ + - ]: 2 : << std::endl;
1110 : :
1111 : 2 : return 0;// no error
1112 : : }
1113 : :
1114 [ + - ][ + - ]: 312 : } // namespace MeshKit
|
