source: ThirdParty/mpqc_open/src/lib/math/isosurf/surfst.cc@ e9f307

Candidate_v1.7.0 stable
Last change on this file since e9f307 was 860145, checked in by Frederik Heber <heber@…>, 9 years ago

Merge commit '0b990dfaa8c6007a996d030163a25f7f5fc8a7e7' as 'ThirdParty/mpqc_open'
  • Property mode set to 100644
File size: 18.2 KB
Line 
1//
2// surfst.cc
3//
4// Copyright (C) 1996 Limit Point Systems, Inc.
5//
6// Author: Curtis Janssen <cljanss@limitpt.com>
7// Maintainer: LPS
8//
9// This file is part of the SC Toolkit.
10//
11// The SC Toolkit is free software; you can redistribute it and/or modify
12// it under the terms of the GNU Library General Public License as published by
13// the Free Software Foundation; either version 2, or (at your option)
14// any later version.
15//
16// The SC Toolkit is distributed in the hope that it will be useful,
17// but WITHOUT ANY WARRANTY; without even the implied warranty of
18// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19// GNU Library General Public License for more details.
20//
21// You should have received a copy of the GNU Library General Public License
22// along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24//
25// The U.S. Government is granted a limited license as per AL 91-7.
26//
27
28#include <stdio.h>
29#include <math.h>
30
31#include <util/misc/formio.h>
32
33#include <math/scmat/matrix.h>
34#include <math/isosurf/surf.h>
35
36using namespace std;
37using namespace sc;
38
39#ifndef WRITE_OOGL // this is useful for debugging this routine
40#define WRITE_OOGL 1
41#endif
42
43#if WRITE_OOGL
44#include <util/render/oogl.h>
45#include <util/render/polygons.h>
46#endif
47
48void
49TriangulatedSurface::remove_slender_triangles(
50 int remove_slender, double height_cutoff,
51 int remove_small, double area_cutoff,
52 const Ref<Volume> &vol, double isoval)
53{
54 int i,j,k;
55 std::set<Ref<Triangle> >::iterator it,jt,kt;
56 std::set<Ref<Edge> >::iterator ie,je,ke;
57 std::set<Ref<Vertex> >::iterator iv,jv;
58
59 int surface_was_completed = _completed_surface;
60
61 if (!_completed_surface) {
62 complete_ref_arrays();
63 }
64 else {
65 clear_int_arrays();
66 }
67
68 _have_values = 0;
69 _values.clear();
70
71 if (_verbose) {
72 ExEnv::outn() << "TriangulatedSurface::remove_slender_triangles:" << endl
73 << "initial: ";
74 topology_info();
75 }
76
77#if WRITE_OOGL
78 if (_debug) {
79 render(new OOGLRender("surfstinit.oogl"));
80 }
81#endif
82
83 int deleted_edges_length;
84 do {
85 std::set<Ref<Triangle> > deleted_triangles;
86 std::set<Ref<Edge> > deleted_edges;
87 std::set<Ref<Vertex> > deleted_vertices;
88
89 std::set<Ref<Vertex> > vertices_of_deleted_triangles;
90
91 std::set<Ref<Triangle> > new_triangles;
92 std::set<Ref<Edge> > new_edges;
93 std::set<Ref<Vertex> > new_vertices;
94
95 // a vertex to set-of-connected-triangles map
96 std::map<Ref<Vertex>,std::set<Ref<Triangle> > > connected_triangle_map;
97 for (it = _triangles.begin(); it != _triangles.end(); it++) {
98 Ref<Triangle> tri = *it;
99 for (j = 0; j<3; j++) {
100 Ref<Vertex> v = tri->vertex(j);
101 connected_triangle_map[v].insert(tri);
102 }
103 }
104
105 // a vertex to set-of-connected-edges map
106 std::map<Ref<Vertex>,std::set<Ref<Edge> > > connected_edge_map;
107 for (ie = _edges.begin(); ie != _edges.end(); ie++) {
108 Ref<Edge> e = *ie;
109 for (j = 0; j<2; j++) {
110 Ref<Vertex> v = e->vertex(j);
111 connected_edge_map[v].insert(e);
112 }
113 }
114
115 for (it = _triangles.begin(); it != _triangles.end(); it++) {
116 Ref<Triangle> tri = *it;
117 // find the heights of the vertices in tri
118 double l[3], l2[3], h[3];
119 for (j=0; j<3; j++) {
120 l[j] = tri->edge(j)->straight_length();
121 if (l[j] <= 0.0) {
122 ExEnv::errn() << "TriangulatedSurface::"
123 << "remove_slender_triangles: bad edge length"
124 << endl;
125 abort();
126 }
127 l2[j] = l[j]*l[j];
128 }
129 double y = 2.0*(l2[0]*l2[1]+l2[0]*l2[2]+l2[1]*l2[2])
130 - l2[0]*l2[0] - l2[1]*l2[1] - l2[2]*l2[2];
131 if (y < 0.0) y = 0.0;
132 double x = 0.5*sqrt(y);
133 for (j=0; j<3; j++) h[j] = x/l[j];
134
135 // find the shortest height
136 int hmin;
137 if (h[0] < h[1]) hmin = 0;
138 else hmin = 1;
139 if (h[2] < h[hmin]) hmin = 2;
140
141 // see if the shortest height is below the cutoff
142 if (remove_slender && h[hmin] < height_cutoff) {
143 // find the vertex that gets eliminated
144 Ref<Vertex> vertex;
145 for (j=0; j<3; j++) {
146 if (tri->vertex(j) != tri->edge(hmin)->vertex(0)
147 &&tri->vertex(j) != tri->edge(hmin)->vertex(1)) {
148 vertex = tri->vertex(j);
149 break;
150 }
151 }
152 std::set<Ref<Triangle> > connected_triangles;
153 connected_triangles = connected_triangle_map[vertex];
154
155 // if one of the connected triangles has a vertex
156 // in a deleted triangle, save this one until the
157 // next pass
158 int skip = 0;
159 for (jt = connected_triangles.begin();
160 jt != connected_triangles.end();
161 jt++) {
162 Ref<Triangle> tri = *jt;
163 for (j=0; j<3; j++) {
164 Ref<Vertex> v = tri->vertex(j);
165 if (vertices_of_deleted_triangles.find(v)
166 != vertices_of_deleted_triangles.end()) {
167 skip = 1;
168 break;
169 }
170 }
171 if (skip) break;
172 }
173 if (skip) continue;
174
175 // find all of the edges contained in the connected triangles
176 std::set<Ref<Edge> > all_edges;
177 for (jt = connected_triangles.begin();
178 jt != connected_triangles.end();
179 jt++) {
180 Ref<Triangle> ctri = *jt;
181 Ref<Edge> e0 = ctri->edge(0);
182 Ref<Edge> e1 = ctri->edge(1);
183 Ref<Edge> e2 = ctri->edge(2);
184 all_edges.insert(e0);
185 all_edges.insert(e1);
186 all_edges.insert(e2);
187 }
188 // find all of the edges connected to the deleted vertex
189 // (including the short edge)
190 std::set<Ref<Edge> > connected_edges;
191 connected_edges = connected_edge_map[vertex];
192 // find the edges forming the perimeter of the deleted triangles
193 // (these are used to form the new triangles)
194 std::set<Ref<Edge> > perimeter_edges;
195 perimeter_edges = all_edges;
196 erase_elements_by_value(perimeter_edges,
197 connected_edges.begin(),
198 connected_edges.end());
199
200 // If deleting this point causes a flattened piece of
201 // surface, reject it. This is tested by checking for
202 // triangles that have all vertices contained in the set
203 // of vertices connected to the deleted vertex.
204 std::set<Ref<Vertex> > connected_vertices;
205 for (je = perimeter_edges.begin(); je != perimeter_edges.end();
206 je++) {
207 Ref<Edge> e = *je;
208 Ref<Vertex> v0 = e->vertex(0);
209 Ref<Vertex> v1 = e->vertex(1);
210 connected_vertices.insert(v0);
211 connected_vertices.insert(v1);
212 }
213 std::set<Ref<Triangle> > triangles_connected_to_perimeter;
214 for (jv = connected_vertices.begin();
215 jv != connected_vertices.end();
216 jv++) {
217 triangles_connected_to_perimeter.insert(
218 connected_triangle_map[*jv].begin(),
219 connected_triangle_map[*jv].end());
220 }
221 for (jt = triangles_connected_to_perimeter.begin();
222 jt != triangles_connected_to_perimeter.end();
223 jt++) {
224 Ref<Triangle> t = *jt;
225 Ref<Vertex> v0 = t->vertex(0);
226 Ref<Vertex> v1 = t->vertex(1);
227 Ref<Vertex> v2 = t->vertex(2);
228 if (connected_vertices.find(v0)!=connected_vertices.end()
229 &&connected_vertices.find(v1)!=connected_vertices.end()
230 &&connected_vertices.find(v2)!=connected_vertices.end())
231 {
232 skip = 1;
233 break;
234 }
235 }
236 if (skip) {
237 continue;
238 }
239
240 deleted_triangles.insert(connected_triangles.begin(),
241 connected_triangles.end());
242 deleted_vertices.insert(vertex);
243 deleted_edges.insert(connected_edges.begin(),
244 connected_edges.end());
245
246 for (jt = deleted_triangles.begin();
247 jt != deleted_triangles.end();
248 jt++) {
249 Ref<Triangle> t = *jt;
250 for (j=0; j<2; j++) {
251 Ref<Vertex> v = t->vertex(j);
252 vertices_of_deleted_triangles.insert(v);
253 }
254 }
255
256 // find a new point that replaces the deleted vertex
257 // (for now use one of the original, since it must lie on the
258 // analytic surface)
259 Ref<Vertex> replacement_vertex;
260 Ref<Edge> short_edge;
261 if (hmin==0) {
262 if (l[1] < l[2]) short_edge = tri->edge(1);
263 else short_edge = tri->edge(2);
264 }
265 else if (hmin==1) {
266 if (l[0] < l[2]) short_edge = tri->edge(0);
267 else short_edge = tri->edge(2);
268 }
269 else {
270 if (l[0] < l[1]) short_edge = tri->edge(0);
271 else short_edge = tri->edge(1);
272 }
273 if (short_edge->vertex(0) == vertex) {
274 replacement_vertex = short_edge->vertex(1);
275 }
276 else {
277 replacement_vertex = short_edge->vertex(0);
278 }
279 new_vertices.insert(replacement_vertex);
280 // for each vertex on the perimeter form a new edge to the
281 // replacement vertex (unless the replacement vertex
282 // is equal to the perimeter vertex)
283 std::map<Ref<Vertex>,Ref<Edge> > new_edge_map;
284 for (je = perimeter_edges.begin(); je != perimeter_edges.end();
285 je++) {
286 Ref<Edge> e = *je;
287 for (k = 0; k<2; k++) {
288 Ref<Vertex> v = e->vertex(k);
289 if (v == replacement_vertex) continue;
290 if (new_edge_map.find(v) == new_edge_map.end()) {
291 Ref<Edge> new_e;
292 // if the edge already exists then use the
293 // existing edge
294 for (ke = perimeter_edges.begin();
295 ke != perimeter_edges.end();
296 ke++) {
297 Ref<Edge> tmp = *ke;
298 if ((tmp->vertex(0) == replacement_vertex
299 &&tmp->vertex(1) == v)
300 ||(tmp->vertex(1) == replacement_vertex
301 &&tmp->vertex(0) == v)) {
302 new_e = tmp;
303 break;
304 }
305 }
306 if (ke == perimeter_edges.end()) {
307 new_e = newEdge(replacement_vertex,
308 v);
309 }
310 new_edge_map[v] = new_e;
311 new_edges.insert(new_e);
312 }
313 }
314 }
315 // for each edge on the perimeter form a new triangle with the
316 // replacement vertex (unless the edge contains the replacement
317 // vertex)
318 for (je = perimeter_edges.begin(); je != perimeter_edges.end();
319 je++) {
320 Ref<Edge> e1 = *je;
321 Ref<Vertex> v0 = e1->vertex(0);
322 Ref<Vertex> v1 = e1->vertex(1);
323 Ref<Edge> e2 = new_edge_map[v0];
324 Ref<Edge> e3 = new_edge_map[v1];
325 if (v0 == replacement_vertex
326 || v1 == replacement_vertex) continue;
327 // Compute the correct orientation of e1 within the new
328 // triangle, by finding the orientation within the old
329 // triangle.
330 int orientation = 0;
331 for (kt = connected_triangles.begin();
332 kt != connected_triangles.end();
333 kt++) {
334 if ((*kt)->contains(e1)) {
335 orientation
336 = (*kt)->orientation(e1);
337 break;
338 }
339 }
340 Ref<Triangle> newtri(newTriangle(e1,e2,e3,orientation));
341 new_triangles.insert(newtri);
342 }
343 }
344 }
345
346#if WRITE_OOGL
347 if (_debug) {
348 char filename[100];
349 static int pass = 0;
350 sprintf(filename, "surfst%04d.oogl", pass);
351 ExEnv::outn() << scprintf("PASS = %04d\n", pass);
352 Ref<Render> render = new OOGLRender(filename);
353 Ref<RenderedPolygons> poly = new RenderedPolygons;
354 poly->initialize(_vertices.size(), _triangles.size(),
355 RenderedPolygons::Vertex);
356 // the number of triangles and edges touching a vertex
357 int *n_triangle = new int[_vertices.size()];
358 int *n_edge = new int[_vertices.size()];
359 memset(n_triangle,0,sizeof(int)*_vertices.size());
360 memset(n_edge,0,sizeof(int)*_vertices.size());
361 std::set<Ref<Triangle> >::iterator it;
362 std::set<Ref<Edge> >::iterator ie;
363 std::set<Ref<Vertex> >::iterator iv;
364 std::map<Ref<Vertex>, int> vertex_to_index;
365 int i = 0;
366 for (iv = _vertices.begin(); iv != _vertices.end(); iv++, i++) {
367 Ref<Vertex> v = *iv;
368 vertex_to_index[v] = i;
369 poly->set_vertex(i,
370 v->point()[0],
371 v->point()[1],
372 v->point()[2]);
373 if (deleted_vertices.find(v) != deleted_vertices.end()) {
374 poly->set_vertex_rgb(i, 1.0, 0.0, 0.0);
375 }
376 else {
377 poly->set_vertex_rgb(i, 0.3, 0.3, 0.3);
378 }
379 }
380 i = 0;
381 for (it = _triangles.begin(); it != _triangles.end(); it++, i++) {
382 Ref<Triangle> t = *it;
383 int i0 = vertex_to_index[t->vertex(0)];
384 int i1 = vertex_to_index[t->vertex(1)];
385 int i2 = vertex_to_index[t->vertex(2)];
386 n_triangle[i0]++;
387 n_triangle[i1]++;
388 n_triangle[i2]++;
389 poly->set_face(i,i0,i1,i2);
390 }
391 for (ie = _edges.begin(); ie != _edges.end(); ie++, i++) {
392 Ref<Edge> e = *ie;
393 int i0 = vertex_to_index[e->vertex(0)];
394 int i1 = vertex_to_index[e->vertex(1)];
395 n_edge[i0]++;
396 n_edge[i1]++;
397 }
398 i = 0;
399 for (iv = _vertices.begin(); iv != _vertices.end(); iv++, i++) {
400 Ref<Vertex> v = *iv;
401 if (n_triangle[i] != n_edge[i]) {
402 ExEnv::outn() << "found bad vertex"
403 << " nedge = " << n_edge[i]
404 << " ntriangle = " << n_triangle[i]
405 << endl;
406 if (deleted_vertices.find(v) != deleted_vertices.end()) {
407 poly->set_vertex_rgb(i, 1.0, 1.0, 0.0);
408 }
409 else {
410 poly->set_vertex_rgb(i, 0.0, 1.0, 0.0);
411 }
412 }
413 }
414 render->render(poly.pointer());
415 pass++;
416 delete[] n_triangle;
417 delete[] n_edge;
418 }
419#endif
420
421 erase_elements_by_value(_triangles,
422 deleted_triangles.begin(),
423 deleted_triangles.end());
424 erase_elements_by_value(_edges,
425 deleted_edges.begin(),
426 deleted_edges.end());
427 erase_elements_by_value(_vertices,
428 deleted_vertices.begin(),
429 deleted_vertices.end());
430
431 _triangles.insert(new_triangles.begin(), new_triangles.end());
432 _edges.insert(new_edges.begin(), new_edges.end());
433 _vertices.insert(new_vertices.begin(), new_vertices.end());
434
435 if (_verbose) {
436 ExEnv::outn() << "intermediate: ";
437 topology_info();
438 }
439
440 deleted_edges_length = deleted_edges.size();
441 } while(deleted_edges_length != 0);
442
443 // fix the index maps
444 _vertex_to_index.clear();
445 _edge_to_index.clear();
446 _triangle_to_index.clear();
447
448 _index_to_vertex.clear();
449 _index_to_edge.clear();
450 _index_to_triangle.clear();
451
452 _index_to_vertex.resize(_vertices.size());
453 for (i=0, iv = _vertices.begin(); iv != _vertices.end(); i++, iv++) {
454 _vertex_to_index[*iv] = i;
455 _index_to_vertex[i] = *iv;
456 }
457
458 _index_to_edge.resize(_edges.size());
459 for (i=0, ie = _edges.begin(); ie != _edges.end(); i++, ie++) {
460 _edge_to_index[*ie] = i;
461 _index_to_edge[i] = *ie;
462 }
463
464 _index_to_triangle.resize(_triangles.size());
465 for (i=0, it = _triangles.begin(); it != _triangles.end(); i++, it++) {
466 _triangle_to_index[*it] = i;
467 _index_to_triangle[i] = *it;
468 }
469
470 // restore the int arrays if they were there to begin with
471 if (surface_was_completed) {
472 complete_int_arrays();
473 _completed_surface = 1;
474 }
475
476 if (_verbose) {
477 ExEnv::outn() << "final: ";
478 topology_info();
479 }
480}
481
482
483/////////////////////////////////////////////////////////////////////////////
484
485// Local Variables:
486// mode: c++
487// c-file-style: "CLJ"
488// End:
Note: See TracBrowser for help on using the repository browser.