1 /* Copyright (c) 2011 Khaled Mamou (kmamou at gmail dot com)
2 All rights reserved.
3
4
5 Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
6
7 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
8
9 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
10
11 3. The names of the contributors may not be used to endorse or promote products derived from this software without specific prior written permission.
12
13 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
14 */
15 #ifndef _CRT_SECURE_NO_WARNINGS
16 #define _CRT_SECURE_NO_WARNINGS
17 #endif //_CRT_SECURE_NO_WARNINGS
18
19 #include <sstream>
20 #include "hacdGraph.h"
21 #include "hacdHACD.h"
22 #include "hacdICHull.h"
23 #include <string.h>
24 #include <algorithm>
25 #include <iterator>
26 #include <limits>
27 #include "assert.h"
28
29 bool gCancelRequest = false;
30 namespace HACD
31 {
Concavity(ICHull & ch,std::map<long,DPoint> & distPoints)32 double HACD::Concavity(ICHull& ch, std::map<long, DPoint>& distPoints)
33 {
34 double concavity = 0.0;
35 double distance = 0.0;
36 std::map<long, DPoint>::iterator itDP(distPoints.begin());
37 std::map<long, DPoint>::iterator itDPEnd(distPoints.end());
38 for (; itDP != itDPEnd; ++itDP)
39 {
40 if (!(itDP->second).m_computed)
41 {
42 if (itDP->first >= 0)
43 {
44 distance = ch.ComputeDistance(itDP->first, m_points[itDP->first], m_normals[itDP->first], (itDP->second).m_computed, true);
45 }
46 else
47 {
48 distance = ch.ComputeDistance(itDP->first, m_facePoints[-itDP->first - 1], m_faceNormals[-itDP->first - 1], (itDP->second).m_computed, true);
49 }
50 }
51 else
52 {
53 distance = (itDP->second).m_dist;
54 }
55 if (concavity < distance)
56 {
57 concavity = distance;
58 }
59 }
60 return concavity;
61 }
62
CreateGraph()63 void HACD::CreateGraph()
64 {
65 // vertex to triangle adjacency information
66 std::vector<std::set<long> > vertexToTriangles;
67 vertexToTriangles.resize(m_nPoints);
68 for (size_t t = 0; t < m_nTriangles; ++t)
69 {
70 vertexToTriangles[m_triangles[t].X()].insert(static_cast<long>(t));
71 vertexToTriangles[m_triangles[t].Y()].insert(static_cast<long>(t));
72 vertexToTriangles[m_triangles[t].Z()].insert(static_cast<long>(t));
73 }
74
75 m_graph.Clear();
76 m_graph.Allocate(m_nTriangles, 5 * m_nTriangles);
77 unsigned long long tr1[3];
78 unsigned long long tr2[3];
79 long i1, j1, k1, i2, j2, k2;
80 long t1, t2;
81 for (size_t v = 0; v < m_nPoints; v++)
82 {
83 std::set<long>::const_iterator it1(vertexToTriangles[v].begin()), itEnd(vertexToTriangles[v].end());
84 for (; it1 != itEnd; ++it1)
85 {
86 t1 = *it1;
87 i1 = m_triangles[t1].X();
88 j1 = m_triangles[t1].Y();
89 k1 = m_triangles[t1].Z();
90 tr1[0] = GetEdgeIndex(i1, j1);
91 tr1[1] = GetEdgeIndex(j1, k1);
92 tr1[2] = GetEdgeIndex(k1, i1);
93 std::set<long>::const_iterator it2(it1);
94 for (++it2; it2 != itEnd; ++it2)
95 {
96 t2 = *it2;
97 i2 = m_triangles[t2].X();
98 j2 = m_triangles[t2].Y();
99 k2 = m_triangles[t2].Z();
100 tr2[0] = GetEdgeIndex(i2, j2);
101 tr2[1] = GetEdgeIndex(j2, k2);
102 tr2[2] = GetEdgeIndex(k2, i2);
103 int shared = 0;
104 for (int i = 0; i < 3; ++i)
105 {
106 for (int j = 0; j < 3; ++j)
107 {
108 if (tr1[i] == tr2[j])
109 {
110 shared++;
111 }
112 }
113 }
114 if (shared == 1) // two triangles are connected if they share exactly one edge
115 {
116 m_graph.AddEdge(t1, t2);
117 }
118 }
119 }
120 }
121 if (m_ccConnectDist >= 0.0)
122 {
123 m_graph.ExtractCCs();
124 if (m_graph.m_nCCs > 1)
125 {
126 std::vector<std::set<long> > cc2V;
127 cc2V.resize(m_graph.m_nCCs);
128 long cc;
129 for (size_t t = 0; t < m_nTriangles; ++t)
130 {
131 cc = m_graph.m_vertices[t].m_cc;
132 cc2V[cc].insert(m_triangles[t].X());
133 cc2V[cc].insert(m_triangles[t].Y());
134 cc2V[cc].insert(m_triangles[t].Z());
135 }
136
137 for (size_t cc1 = 0; cc1 < m_graph.m_nCCs; ++cc1)
138 {
139 for (size_t cc2 = cc1 + 1; cc2 < m_graph.m_nCCs; ++cc2)
140 {
141 std::set<long>::const_iterator itV1(cc2V[cc1].begin()), itVEnd1(cc2V[cc1].end());
142 for (; itV1 != itVEnd1; ++itV1)
143 {
144 double distC1C2 = std::numeric_limits<double>::max();
145 double dist;
146 t1 = -1;
147 t2 = -1;
148 std::set<long>::const_iterator itV2(cc2V[cc2].begin()), itVEnd2(cc2V[cc2].end());
149 for (; itV2 != itVEnd2; ++itV2)
150 {
151 dist = (m_points[*itV1] - m_points[*itV2]).GetNorm();
152 if (dist < distC1C2)
153 {
154 distC1C2 = dist;
155 t1 = *vertexToTriangles[*itV1].begin();
156
157 std::set<long>::const_iterator it2(vertexToTriangles[*itV2].begin()),
158 it2End(vertexToTriangles[*itV2].end());
159 t2 = -1;
160 for (; it2 != it2End; ++it2)
161 {
162 if (*it2 != t1)
163 {
164 t2 = *it2;
165 break;
166 }
167 }
168 }
169 }
170 if (distC1C2 <= m_ccConnectDist && t1 > 0 && t2 > 0)
171 {
172 m_graph.AddEdge(t1, t2);
173 }
174 }
175 }
176 }
177 }
178 }
179 }
InitializeDualGraph()180 void HACD::InitializeDualGraph()
181 {
182 long i, j, k;
183 Vec3<Real> u, v, w, normal;
184 delete[] m_normals;
185 m_normals = new Vec3<Real>[m_nPoints];
186 if (m_addFacesPoints)
187 {
188 delete[] m_facePoints;
189 delete[] m_faceNormals;
190 m_facePoints = new Vec3<Real>[m_nTriangles];
191 m_faceNormals = new Vec3<Real>[m_nTriangles];
192 }
193 memset(m_normals, 0, sizeof(Vec3<Real>) * m_nPoints);
194 for (unsigned long f = 0; f < m_nTriangles; f++)
195 {
196 if (m_callBack) (*m_callBack)("+ InitializeDualGraph\n", f, m_nTriangles, 0);
197
198 if (gCancelRequest)
199 return;
200
201 i = m_triangles[f].X();
202 j = m_triangles[f].Y();
203 k = m_triangles[f].Z();
204
205 m_graph.m_vertices[f].m_distPoints[i].m_distOnly = false;
206 m_graph.m_vertices[f].m_distPoints[j].m_distOnly = false;
207 m_graph.m_vertices[f].m_distPoints[k].m_distOnly = false;
208
209 ICHull* ch = new ICHull;
210 m_graph.m_vertices[f].m_convexHull = ch;
211 ch->AddPoint(m_points[i], i);
212 ch->AddPoint(m_points[j], j);
213 ch->AddPoint(m_points[k], k);
214 ch->SetDistPoints(&m_graph.m_vertices[f].m_distPoints);
215
216 u = m_points[j] - m_points[i];
217 v = m_points[k] - m_points[i];
218 w = m_points[k] - m_points[j];
219 normal = u ^ v;
220
221 m_normals[i] += normal;
222 m_normals[j] += normal;
223 m_normals[k] += normal;
224
225 m_graph.m_vertices[f].m_surf = normal.GetNorm();
226 m_graph.m_vertices[f].m_perimeter = u.GetNorm() + v.GetNorm() + w.GetNorm();
227
228 normal.Normalize();
229
230 m_graph.m_vertices[f].m_boudaryEdges.insert(GetEdgeIndex(i, j));
231 m_graph.m_vertices[f].m_boudaryEdges.insert(GetEdgeIndex(j, k));
232 m_graph.m_vertices[f].m_boudaryEdges.insert(GetEdgeIndex(k, i));
233 if (m_addFacesPoints)
234 {
235 m_faceNormals[f] = normal;
236 m_facePoints[f] = (m_points[i] + m_points[j] + m_points[k]) / 3.0;
237 m_graph.m_vertices[f].m_distPoints[-static_cast<long>(f) - 1].m_distOnly = true;
238 }
239 if (m_addExtraDistPoints)
240 { // we need a kd-tree structure to accelerate this part!
241 long i1, j1, k1;
242 Vec3<Real> u1, v1, normal1;
243 normal = -normal;
244 double distance = 0.0;
245 double distMin = 0.0;
246 size_t faceIndex = m_nTriangles;
247 Vec3<Real> seedPoint((m_points[i] + m_points[j] + m_points[k]) / 3.0);
248 long nhit = 0;
249 for (size_t f1 = 0; f1 < m_nTriangles; f1++)
250 {
251 i1 = m_triangles[f1].X();
252 j1 = m_triangles[f1].Y();
253 k1 = m_triangles[f1].Z();
254 u1 = m_points[j1] - m_points[i1];
255 v1 = m_points[k1] - m_points[i1];
256 normal1 = (u1 ^ v1);
257 if (normal * normal1 > 0.0)
258 {
259 nhit = IntersectRayTriangle(Vec3<double>(seedPoint.X(), seedPoint.Y(), seedPoint.Z()),
260 Vec3<double>(normal.X(), normal.Y(), normal.Z()),
261 Vec3<double>(m_points[i1].X(), m_points[i1].Y(), m_points[i1].Z()),
262 Vec3<double>(m_points[j1].X(), m_points[j1].Y(), m_points[j1].Z()),
263 Vec3<double>(m_points[k1].X(), m_points[k1].Y(), m_points[k1].Z()),
264 distance);
265 if ((nhit == 1) && ((distMin > distance) || (faceIndex == m_nTriangles)))
266 {
267 distMin = distance;
268 faceIndex = f1;
269 }
270 }
271 }
272 if (faceIndex < m_nTriangles)
273 {
274 i1 = m_triangles[faceIndex].X();
275 j1 = m_triangles[faceIndex].Y();
276 k1 = m_triangles[faceIndex].Z();
277 m_graph.m_vertices[f].m_distPoints[i1].m_distOnly = true;
278 m_graph.m_vertices[f].m_distPoints[j1].m_distOnly = true;
279 m_graph.m_vertices[f].m_distPoints[k1].m_distOnly = true;
280 if (m_addFacesPoints)
281 {
282 m_graph.m_vertices[f].m_distPoints[-static_cast<long>(faceIndex) - 1].m_distOnly = true;
283 }
284 }
285 }
286 }
287 for (size_t v = 0; v < m_nPoints; v++)
288 {
289 m_normals[v].Normalize();
290 }
291 }
292
NormalizeData()293 void HACD::NormalizeData()
294 {
295 if (m_nPoints == 0)
296 {
297 return;
298 }
299 m_barycenter = m_points[0];
300 Vec3<Real> min = m_points[0];
301 Vec3<Real> max = m_points[0];
302 Real x, y, z;
303 for (size_t v = 1; v < m_nPoints; v++)
304 {
305 m_barycenter += m_points[v];
306 x = m_points[v].X();
307 y = m_points[v].Y();
308 z = m_points[v].Z();
309 if (x < min.X())
310 min.X() = x;
311 else if (x > max.X())
312 max.X() = x;
313 if (y < min.Y())
314 min.Y() = y;
315 else if (y > max.Y())
316 max.Y() = y;
317 if (z < min.Z())
318 min.Z() = z;
319 else if (z > max.Z())
320 max.Z() = z;
321 }
322 m_barycenter /= static_cast<Real>(m_nPoints);
323 m_diag = (max - min).GetNorm();
324 const Real invDiag = static_cast<Real>(2.0 * m_scale / m_diag);
325 if (m_diag != 0.0)
326 {
327 for (size_t v = 0; v < m_nPoints; v++)
328 {
329 m_points[v] = (m_points[v] - m_barycenter) * invDiag;
330 }
331 }
332 }
DenormalizeData()333 void HACD::DenormalizeData()
334 {
335 if (m_nPoints == 0)
336 {
337 return;
338 }
339 if (m_diag != 0.0)
340 {
341 const Real diag = static_cast<Real>(m_diag / (2.0 * m_scale));
342 for (size_t v = 0; v < m_nPoints; v++)
343 {
344 m_points[v] = m_points[v] * diag + m_barycenter;
345 }
346 }
347 }
HACD(void)348 HACD::HACD(void)
349 {
350 m_convexHulls = 0;
351 m_triangles = 0;
352 m_points = 0;
353 m_normals = 0;
354 m_nTriangles = 0;
355 m_nPoints = 0;
356 m_nClusters = 0;
357 m_concavity = 0.0;
358 m_diag = 1.0;
359 m_barycenter = Vec3<Real>(0.0, 0.0, 0.0);
360 m_alpha = 0.1;
361 m_beta = 0.1;
362 m_nVerticesPerCH = 30;
363 m_callBack = 0;
364 m_addExtraDistPoints = false;
365 m_addNeighboursDistPoints = false;
366 m_scale = 1000.0;
367 m_partition = 0;
368 m_nMinClusters = 3;
369 m_facePoints = 0;
370 m_faceNormals = 0;
371 m_ccConnectDist = 30;
372 }
~HACD(void)373 HACD::~HACD(void)
374 {
375 delete[] m_normals;
376 delete[] m_convexHulls;
377 delete[] m_partition;
378 delete[] m_facePoints;
379 delete[] m_faceNormals;
380 }
381 int iteration = 0;
ComputeEdgeCost(size_t e)382 void HACD::ComputeEdgeCost(size_t e)
383 {
384 GraphEdge& gE = m_graph.m_edges[e];
385 long v1 = gE.m_v1;
386 long v2 = gE.m_v2;
387
388 if (m_graph.m_vertices[v2].m_distPoints.size() > m_graph.m_vertices[v1].m_distPoints.size())
389 {
390 gE.m_v1 = v2;
391 gE.m_v2 = v1;
392 //std::swap<long>(v1, v2);
393 std::swap(v1, v2);
394 }
395 GraphVertex& gV1 = m_graph.m_vertices[v1];
396 GraphVertex& gV2 = m_graph.m_vertices[v2];
397
398 // delete old convex-hull
399 delete gE.m_convexHull;
400 // create the edge's convex-hull
401 ICHull* ch = new ICHull;
402 gE.m_convexHull = ch;
403 (*ch) = (*gV1.m_convexHull);
404
405 // update distPoints
406 gE.m_distPoints = gV1.m_distPoints;
407 std::map<long, DPoint>::iterator itDP(gV2.m_distPoints.begin());
408 std::map<long, DPoint>::iterator itDPEnd(gV2.m_distPoints.end());
409 std::map<long, DPoint>::iterator itDP1;
410
411 for (; itDP != itDPEnd; ++itDP)
412 {
413 itDP1 = gE.m_distPoints.find(itDP->first);
414 if (itDP1 == gE.m_distPoints.end())
415 {
416 gE.m_distPoints[itDP->first].m_distOnly = (itDP->second).m_distOnly;
417 if (!(itDP->second).m_distOnly)
418 {
419 ch->AddPoint(m_points[itDP->first], itDP->first);
420 }
421 }
422 else
423 {
424 if ((itDP1->second).m_distOnly && !(itDP->second).m_distOnly)
425 {
426 gE.m_distPoints[itDP->first].m_distOnly = false;
427 ch->AddPoint(m_points[itDP->first], itDP->first);
428 }
429 }
430 }
431
432 ch->SetDistPoints(&gE.m_distPoints);
433 // create the convex-hull
434 while (ch->Process() == ICHullErrorInconsistent) // if we face problems when constructing the visual-hull. really ugly!!!!
435 {
436 // if (m_callBack) (*m_callBack)("\t Problem with convex-hull construction [HACD::ComputeEdgeCost]\n", 0.0, 0.0, 0);
437 ch = new ICHull;
438 CircularList<TMMVertex>& verticesCH = (gE.m_convexHull)->GetMesh().m_vertices;
439 size_t nV = verticesCH.GetSize();
440 long ptIndex = 0;
441 verticesCH.Next();
442 for (size_t v = 1; v < nV; ++v)
443 {
444 ptIndex = verticesCH.GetHead()->GetData().m_name;
445 if (ptIndex != ICHull::sc_dummyIndex /* && ptIndex < m_nPoints*/)
446 ch->AddPoint(m_points[ptIndex], ptIndex);
447 verticesCH.Next();
448 }
449 delete gE.m_convexHull;
450 gE.m_convexHull = ch;
451 }
452 double volume = 0.0;
453 double concavity = 0.0;
454 if (ch->IsFlat())
455 {
456 bool insideHull;
457 std::map<long, DPoint>::iterator itDP(gE.m_distPoints.begin());
458 std::map<long, DPoint>::iterator itDPEnd(gE.m_distPoints.end());
459 for (; itDP != itDPEnd; ++itDP)
460 {
461 if (itDP->first >= 0)
462 {
463 concavity = std::max<double>(concavity, ch->ComputeDistance(itDP->first, m_points[itDP->first], m_normals[itDP->first], insideHull, false));
464 }
465 }
466 }
467 else
468 {
469 if (m_addNeighboursDistPoints)
470 { // add distance points from adjacent clusters
471 std::set<long> eEdges;
472 std::set_union(gV1.m_edges.begin(),
473 gV1.m_edges.end(),
474 gV2.m_edges.begin(),
475 gV2.m_edges.end(),
476 std::inserter(eEdges, eEdges.begin()));
477
478 std::set<long>::const_iterator ed(eEdges.begin());
479 std::set<long>::const_iterator itEnd(eEdges.end());
480 long a, b, c;
481 for (; ed != itEnd; ++ed)
482 {
483 a = m_graph.m_edges[*ed].m_v1;
484 b = m_graph.m_edges[*ed].m_v2;
485 if (a != v2 && a != v1)
486 {
487 c = a;
488 }
489 else if (b != v2 && b != v1)
490 {
491 c = b;
492 }
493 else
494 {
495 c = -1;
496 }
497 if (c > 0)
498 {
499 GraphVertex& gVC = m_graph.m_vertices[c];
500 std::map<long, DPoint>::iterator itDP(gVC.m_distPoints.begin());
501 std::map<long, DPoint>::iterator itDPEnd(gVC.m_distPoints.end());
502 std::map<long, DPoint>::iterator itDP1;
503 for (; itDP != itDPEnd; ++itDP)
504 {
505 itDP1 = gE.m_distPoints.find(itDP->first);
506 if (itDP1 == gE.m_distPoints.end())
507 {
508 if (itDP->first >= 0 && itDP1 == gE.m_distPoints.end() && ch->IsInside(m_points[itDP->first]))
509 {
510 gE.m_distPoints[itDP->first].m_distOnly = true;
511 }
512 else if (itDP->first < 0 && ch->IsInside(m_facePoints[-itDP->first - 1]))
513 {
514 gE.m_distPoints[itDP->first].m_distOnly = true;
515 }
516 }
517 }
518 }
519 }
520 }
521 concavity = Concavity(*ch, gE.m_distPoints);
522 }
523
524 // compute boudary edges
525 double perimeter = 0.0;
526 double surf = 1.0;
527 if (m_alpha > 0.0)
528 {
529 gE.m_boudaryEdges.clear();
530 std::set_symmetric_difference(gV1.m_boudaryEdges.begin(),
531 gV1.m_boudaryEdges.end(),
532 gV2.m_boudaryEdges.begin(),
533 gV2.m_boudaryEdges.end(),
534 std::inserter(gE.m_boudaryEdges,
535 gE.m_boudaryEdges.begin()));
536
537 std::set<unsigned long long>::const_iterator itBE(gE.m_boudaryEdges.begin());
538 std::set<unsigned long long>::const_iterator itBEEnd(gE.m_boudaryEdges.end());
539 for (; itBE != itBEEnd; ++itBE)
540 {
541 perimeter += (m_points[static_cast<long>((*itBE) >> 32)] -
542 m_points[static_cast<long>((*itBE) & 0xFFFFFFFFULL)])
543 .GetNorm();
544 }
545 surf = gV1.m_surf + gV2.m_surf;
546 }
547 double ratio = perimeter * perimeter / (4.0 * sc_pi * surf);
548 gE.m_volume = (m_beta == 0.0) ? 0.0 : ch->ComputeVolume() / pow(m_scale, 3.0); // cluster's volume
549 gE.m_surf = surf; // cluster's area
550 gE.m_perimeter = perimeter; // cluster's perimeter
551 gE.m_concavity = concavity; // cluster's concavity
552 gE.m_error = static_cast<Real>(concavity + m_alpha * ratio + m_beta * volume); // cluster's priority
553 }
InitializePriorityQueue()554 bool HACD::InitializePriorityQueue()
555 {
556 m_pqueue.reserve(m_graph.m_nE + 100);
557 for (size_t e = 0; e < m_graph.m_nE; ++e)
558 {
559 ComputeEdgeCost(static_cast<long>(e));
560 m_pqueue.push(GraphEdgePriorityQueue(static_cast<long>(e), m_graph.m_edges[e].m_error));
561 }
562 return true;
563 }
Simplify()564 void HACD::Simplify()
565 {
566 long v1 = -1;
567 long v2 = -1;
568 double progressOld = -1.0;
569 double progress = 0.0;
570 double globalConcavity = 0.0;
571 char msg[1024];
572 double ptgStep = 1.0;
573 while ((globalConcavity < m_concavity) &&
574 (m_graph.GetNVertices() > m_nMinClusters) &&
575 (m_graph.GetNEdges() > 0))
576 {
577 progress = 100.0 - m_graph.GetNVertices() * 100.0 / m_nTriangles;
578 if (fabs(progress - progressOld) > ptgStep && m_callBack)
579 {
580 sprintf(msg, "%3.2f %% V = %lu \t C = %f \t \t \r", progress, static_cast<unsigned long>(m_graph.GetNVertices()), globalConcavity);
581 (*m_callBack)(msg, progress, globalConcavity, m_graph.GetNVertices());
582 progressOld = progress;
583 if (progress > 99.0)
584 {
585 ptgStep = 0.01;
586 }
587 else if (progress > 90.0)
588 {
589 ptgStep = 0.1;
590 }
591 }
592
593 GraphEdgePriorityQueue currentEdge(0, 0.0);
594 bool done = false;
595 do
596 {
597 done = false;
598 if (m_pqueue.size() == 0)
599 {
600 done = true;
601 break;
602 }
603 currentEdge = m_pqueue.top();
604 m_pqueue.pop();
605 } while (m_graph.m_edges[currentEdge.m_name].m_deleted ||
606 m_graph.m_edges[currentEdge.m_name].m_error != currentEdge.m_priority);
607
608 if (m_graph.m_edges[currentEdge.m_name].m_concavity < m_concavity && !done)
609 {
610 globalConcavity = std::max<double>(globalConcavity, m_graph.m_edges[currentEdge.m_name].m_concavity);
611 v1 = m_graph.m_edges[currentEdge.m_name].m_v1;
612 v2 = m_graph.m_edges[currentEdge.m_name].m_v2;
613 // update vertex info
614 m_graph.m_vertices[v1].m_error = m_graph.m_edges[currentEdge.m_name].m_error;
615 m_graph.m_vertices[v1].m_surf = m_graph.m_edges[currentEdge.m_name].m_surf;
616 m_graph.m_vertices[v1].m_volume = m_graph.m_edges[currentEdge.m_name].m_volume;
617 m_graph.m_vertices[v1].m_concavity = m_graph.m_edges[currentEdge.m_name].m_concavity;
618 m_graph.m_vertices[v1].m_perimeter = m_graph.m_edges[currentEdge.m_name].m_perimeter;
619 m_graph.m_vertices[v1].m_distPoints = m_graph.m_edges[currentEdge.m_name].m_distPoints;
620 (*m_graph.m_vertices[v1].m_convexHull) = (*m_graph.m_edges[currentEdge.m_name].m_convexHull);
621 (m_graph.m_vertices[v1].m_convexHull)->SetDistPoints(&(m_graph.m_vertices[v1].m_distPoints));
622 m_graph.m_vertices[v1].m_boudaryEdges = m_graph.m_edges[currentEdge.m_name].m_boudaryEdges;
623
624 // We apply the optimal ecol
625 // std::cout << "v1 " << v1 << " v2 " << v2 << std::endl;
626 m_graph.EdgeCollapse(v1, v2);
627 // recompute the adjacent edges costs
628 std::set<long>::const_iterator itE(m_graph.m_vertices[v1].m_edges.begin()),
629 itEEnd(m_graph.m_vertices[v1].m_edges.end());
630 for (; itE != itEEnd; ++itE)
631 {
632 size_t e = *itE;
633 ComputeEdgeCost(static_cast<long>(e));
634 m_pqueue.push(GraphEdgePriorityQueue(static_cast<long>(e), m_graph.m_edges[e].m_error));
635 }
636 }
637 else
638 {
639 break;
640 }
641 }
642 while (!m_pqueue.empty())
643 {
644 m_pqueue.pop();
645 }
646
647 m_cVertices.clear();
648 m_nClusters = m_graph.GetNVertices();
649 m_cVertices.reserve(m_nClusters);
650 for (size_t p = 0, v = 0; v != m_graph.m_vertices.size(); ++v)
651 {
652 if (!m_graph.m_vertices[v].m_deleted)
653 {
654 if (m_callBack)
655 {
656 char msg[1024];
657 sprintf(msg, "\t CH \t %lu \t %lf \t %lf\n", static_cast<unsigned long>(p), m_graph.m_vertices[v].m_concavity, m_graph.m_vertices[v].m_error);
658 (*m_callBack)(msg, 0.0, 0.0, m_nClusters);
659 p++;
660 }
661 m_cVertices.push_back(static_cast<long>(v));
662 }
663 }
664 if (m_callBack)
665 {
666 sprintf(msg, "# clusters = %lu \t C = %f\n", static_cast<unsigned long>(m_nClusters), globalConcavity);
667 (*m_callBack)(msg, progress, globalConcavity, m_graph.GetNVertices());
668 }
669 }
670
Compute(bool fullCH,bool exportDistPoints)671 bool HACD::Compute(bool fullCH, bool exportDistPoints)
672 {
673 gCancelRequest = false;
674
675 if (!m_points || !m_triangles || !m_nPoints || !m_nTriangles)
676 {
677 return false;
678 }
679 size_t nV = m_nTriangles;
680 if (m_callBack)
681 {
682 std::ostringstream msg;
683 msg << "+ Mesh" << std::endl;
684 msg << "\t # vertices \t" << m_nPoints << std::endl;
685 msg << "\t # triangles \t" << m_nTriangles << std::endl;
686 msg << "+ Parameters" << std::endl;
687 msg << "\t min # of clusters \t" << m_nMinClusters << std::endl;
688 msg << "\t max concavity \t" << m_concavity << std::endl;
689 msg << "\t compacity weigth \t" << m_alpha << std::endl;
690 msg << "\t volume weigth \t" << m_beta << std::endl;
691 msg << "\t # vertices per convex-hull \t" << m_nVerticesPerCH << std::endl;
692 msg << "\t scale \t" << m_scale << std::endl;
693 msg << "\t add extra distance points \t" << m_addExtraDistPoints << std::endl;
694 msg << "\t add neighbours distance points \t" << m_addNeighboursDistPoints << std::endl;
695 msg << "\t add face distance points \t" << m_addFacesPoints << std::endl;
696 msg << "\t produce full convex-hulls \t" << fullCH << std::endl;
697 msg << "\t max. distance to connect CCs \t" << m_ccConnectDist << std::endl;
698 (*m_callBack)(msg.str().c_str(), 0.0, 0.0, nV);
699 }
700 if (m_callBack) (*m_callBack)("+ Normalizing Data\n", 0.0, 0.0, nV);
701 NormalizeData();
702 if (m_callBack) (*m_callBack)("+ Creating Graph\n", 0.0, 0.0, nV);
703 CreateGraph();
704 // Compute the surfaces and perimeters of all the faces
705 if (m_callBack) (*m_callBack)("+ Initializing Dual Graph\n", 0.0, 0.0, nV);
706 if (gCancelRequest)
707 return false;
708
709 InitializeDualGraph();
710 if (m_callBack) (*m_callBack)("+ Initializing Priority Queue\n", 0.0, 0.0, nV);
711 if (gCancelRequest)
712 return false;
713
714 InitializePriorityQueue();
715 // we simplify the graph
716 if (m_callBack) (*m_callBack)("+ Simplification ...\n", 0.0, 0.0, m_nTriangles);
717 Simplify();
718 if (m_callBack) (*m_callBack)("+ Denormalizing Data\n", 0.0, 0.0, m_nClusters);
719 DenormalizeData();
720 if (m_callBack) (*m_callBack)("+ Computing final convex-hulls\n", 0.0, 0.0, m_nClusters);
721 delete[] m_convexHulls;
722 m_convexHulls = new ICHull[m_nClusters];
723 delete[] m_partition;
724 m_partition = new long[m_nTriangles];
725 for (size_t p = 0; p != m_cVertices.size(); ++p)
726 {
727 size_t v = m_cVertices[p];
728 m_partition[v] = static_cast<long>(p);
729 for (size_t a = 0; a < m_graph.m_vertices[v].m_ancestors.size(); a++)
730 {
731 m_partition[m_graph.m_vertices[v].m_ancestors[a]] = static_cast<long>(p);
732 }
733 // compute the convex-hull
734 const std::map<long, DPoint>& pointsCH = m_graph.m_vertices[v].m_distPoints;
735 std::map<long, DPoint>::const_iterator itCH(pointsCH.begin());
736 std::map<long, DPoint>::const_iterator itCHEnd(pointsCH.end());
737 for (; itCH != itCHEnd; ++itCH)
738 {
739 if (!(itCH->second).m_distOnly)
740 {
741 m_convexHulls[p].AddPoint(m_points[itCH->first], itCH->first);
742 }
743 }
744 m_convexHulls[p].SetDistPoints(&m_graph.m_vertices[v].m_distPoints);
745 if (fullCH)
746 {
747 m_convexHulls[p].Process();
748 }
749 else
750 {
751 m_convexHulls[p].Process(static_cast<unsigned long>(m_nVerticesPerCH));
752 }
753 if (exportDistPoints)
754 {
755 itCH = pointsCH.begin();
756 for (; itCH != itCHEnd; ++itCH)
757 {
758 if ((itCH->second).m_distOnly)
759 {
760 if (itCH->first >= 0)
761 {
762 m_convexHulls[p].AddPoint(m_points[itCH->first], itCH->first);
763 }
764 else
765 {
766 m_convexHulls[p].AddPoint(m_facePoints[-itCH->first - 1], itCH->first);
767 }
768 }
769 }
770 }
771 }
772 return true;
773 }
774
GetNTrianglesCH(size_t numCH) const775 size_t HACD::GetNTrianglesCH(size_t numCH) const
776 {
777 if (numCH >= m_nClusters)
778 {
779 return 0;
780 }
781 return m_convexHulls[numCH].GetMesh().GetNTriangles();
782 }
GetNPointsCH(size_t numCH) const783 size_t HACD::GetNPointsCH(size_t numCH) const
784 {
785 if (numCH >= m_nClusters)
786 {
787 return 0;
788 }
789 return m_convexHulls[numCH].GetMesh().GetNVertices();
790 }
791
GetCH(size_t numCH,Vec3<Real> * const points,Vec3<long> * const triangles)792 bool HACD::GetCH(size_t numCH, Vec3<Real>* const points, Vec3<long>* const triangles)
793 {
794 if (numCH >= m_nClusters)
795 {
796 return false;
797 }
798 m_convexHulls[numCH].GetMesh().GetIFS(points, triangles);
799 return true;
800 }
801
Save(const char * fileName,bool uniColor,long numCluster) const802 bool HACD::Save(const char* fileName, bool uniColor, long numCluster) const
803 {
804 std::ofstream fout(fileName);
805 if (fout.is_open())
806 {
807 if (m_callBack)
808 {
809 char msg[1024];
810 sprintf(msg, "Saving %s\n", fileName);
811 (*m_callBack)(msg, 0.0, 0.0, m_graph.GetNVertices());
812 }
813 Material mat;
814 if (numCluster < 0)
815 {
816 for (size_t p = 0; p != m_nClusters; ++p)
817 {
818 if (!uniColor)
819 {
820 mat.m_diffuseColor.X() = mat.m_diffuseColor.Y() = mat.m_diffuseColor.Z() = 0.0;
821 while (mat.m_diffuseColor.X() == mat.m_diffuseColor.Y() ||
822 mat.m_diffuseColor.Z() == mat.m_diffuseColor.Y() ||
823 mat.m_diffuseColor.Z() == mat.m_diffuseColor.X())
824 {
825 mat.m_diffuseColor.X() = (rand() % 100) / 100.0;
826 mat.m_diffuseColor.Y() = (rand() % 100) / 100.0;
827 mat.m_diffuseColor.Z() = (rand() % 100) / 100.0;
828 }
829 }
830 m_convexHulls[p].GetMesh().SaveVRML2(fout, mat);
831 }
832 }
833 else if (numCluster < static_cast<long>(m_cVertices.size()))
834 {
835 m_convexHulls[numCluster].GetMesh().SaveVRML2(fout, mat);
836 }
837 fout.close();
838 return true;
839 }
840 else
841 {
842 if (m_callBack)
843 {
844 char msg[1024];
845 sprintf(msg, "Error saving %s\n", fileName);
846 (*m_callBack)(msg, 0.0, 0.0, m_graph.GetNVertices());
847 }
848 return false;
849 }
850 }
851 } // namespace HACD
852