Instrument Neutral Distributed Interface INDI  1.9.5
ConvexHull.cpp
Go to the documentation of this file.
1 
9 // This C++ code is based on the c code described below
10 // it was ported to c++ and modded to work on on coordinate types
11 // other than integers by Roger James in December 2013
12 /*
13 This code is described in "Computational Geometry in C" (Second Edition),
14 Chapter 4. It is not written to be comprehensible without the
15 explanation in that book.
16 
17 Input: 3n integer coordinates for the points.
18 Output: the 3D convex hull, in postscript with embedded comments
19  showing the vertices and faces.
20 
21 Compile: gcc -o chull chull.c (or simply: make)
22 
23 Written by Joseph O'Rourke, with contributions by
24  Kristy Anderson, John Kutcher, Catherine Schevon, Susan Weller.
25 Last modified: May 2000
26 Questions to orourke@cs.smith.edu.
27 
28 --------------------------------------------------------------------
29 This code is Copyright 2000 by Joseph O'Rourke. It may be freely
30 redistributed in its entirety provided that this copyright notice is
31 not removed.
32 --------------------------------------------------------------------
33 */
34 
35 #include "ConvexHull.h"
36 
37 #include <iostream>
38 #include <map>
39 
40 namespace INDI
41 {
42 namespace AlignmentSubsystem
43 {
44 ConvexHull::ConvexHull() : vertices(nullptr), edges(nullptr), faces(nullptr), debug(false), check(false),
45  ScaleFactor(SAFE-1)
46 {
47 }
48 
50 {
51  tFace f;
52  tEdge e, temp;
53  int vol;
54  bool vis = false;
55 
56  if (debug)
57  {
58  std::cerr << "AddOne: starting to add v" << p->vnum << ".\n";
59  //PrintOut( vertices );
60  }
61 
62  /* Mark faces visible from p. */
63  f = faces;
64  do
65  {
66  vol = VolumeSign(f, p);
67  if (debug)
68  std::cerr << "faddr: " << std::hex << f << " paddr: " << p << " Vol = " << std::dec << vol << '\n';
69  if (vol < 0)
70  {
71  f->visible = VISIBLE;
72  vis = true;
73  }
74  f = f->next;
75  } while (f != faces);
76 
77  /* If no faces are visible from p, then p is inside the hull. */
78  if (!vis)
79  {
80  p->onhull = !ONHULL;
81  return false;
82  }
83 
84  /* Mark edges in interior of visible region for deletion.
85  Erect a newface based on each border edge. */
86  e = edges;
87  do
88  {
89  temp = e->next;
90  if (e->adjface[0]->visible && e->adjface[1]->visible)
91  /* e interior: mark for deletion. */
92  e->delete_it = REMOVED;
93  else if (e->adjface[0]->visible || e->adjface[1]->visible)
94  /* e border: make a new face. */
95  e->newface = MakeConeFace(e, p);
96  e = temp;
97  } while (e != edges);
98  return true;
99 }
100 
102 {
103  int i;
104  tFace fstart;
105  tEdge e;
106  tVertex v;
107  bool error = false;
108 
109  fstart = faces;
110  if (faces)
111  do
112  {
113  for (i = 0; i < 3; ++i)
114  {
115  v = faces->vertex[i];
116  e = faces->edge[i];
117  if (v != e->endpts[0] && v != e->endpts[1])
118  {
119  error = true;
120  std::cerr << "CheckEndpts: Error!\n";
121  std::cerr << " addr: " << std::hex << faces << ':';
122  std::cerr << " edges:";
123  std::cerr << "(" << e->endpts[0]->vnum << "," << e->endpts[1]->vnum << ")";
124  std::cerr << "\n";
125  }
126  }
127  faces = faces->next;
128  } while (faces != fstart);
129 
130  if (error)
131  std::cerr << "Checks: ERROR found and reported above.\n";
132  else
133  std::cerr << "Checks: All endpts of all edges of all faces check.\n";
134 }
135 
136 void ConvexHull::CheckEuler(int V, int E, int F)
137 {
138  if (check)
139  std::cerr << "Checks: V, E, F = " << V << ' ' << E << ' ' << F << ":\t";
140 
141  if ((V - E + F) != 2)
142  std::cerr << "Checks: V-E+F != 2\n";
143  else if (check)
144  std::cerr << "V-E+F = 2\t";
145 
146  if (F != (2 * V - 4))
147  std::cerr << "Checks: F=" << F << " != 2V-4=" << 2 * V - 4 << "; V=" << V << '\n';
148  else if (check)
149  std::cerr << "F = 2V-4\t";
150 
151  if ((2 * E) != (3 * F))
152  std::cerr << "Checks: 2E=" << 2 * E << " != 3F=" << 3 * F << "; E=" << E << ", F=" << F << '\n';
153  else if (check)
154  std::cerr << "2E = 3F\n";
155 }
156 
158 {
159  tVertex v;
160  tEdge e;
161  tFace f;
162  int V = 0, E = 0, F = 0;
163 
164  Consistency();
165  Convexity();
166  if ((v = vertices) != nullptr)
167  do
168  {
169  if (v->mark)
170  V++;
171  v = v->next;
172  } while (v != vertices);
173 
174  if ((e = edges) != nullptr)
175  do
176  {
177  E++;
178  e = e->next;
179  } while (e != edges);
180  if ((f = faces) != nullptr)
181  do
182  {
183  F++;
184  f = f->next;
185  } while (f != faces);
186  CheckEuler(V, E, F);
187  CheckEndpts();
188 }
189 
191 {
192  tEdge e; /* Primary index into edge list. */
193  tEdge t; /* Temporary edge pointer. */
194 
195  /* Integrate the newface's into the data structure. */
196  /* Check every edge. */
197  e = edges;
198  do
199  {
200  if (e->newface)
201  {
202  if (e->adjface[0]->visible)
203  e->adjface[0] = e->newface;
204  else
205  e->adjface[1] = e->newface;
206  e->newface = nullptr;
207  }
208  e = e->next;
209  } while (e != edges);
210 
211  /* Delete any edges marked for deletion. */
212  while (edges && edges->delete_it)
213  {
214  e = edges;
215  remove<tEdge>(edges, e);
216  }
217  e = edges->next;
218  do
219  {
220  if (e->delete_it)
221  {
222  t = e;
223  e = e->next;
224  remove<tEdge>(edges, t);
225  }
226  else
227  e = e->next;
228  } while (e != edges);
229 }
230 
232 {
233  tFace f; /* Primary pointer into face list. */
234  tFace t; /* Temporary pointer, for deleting. */
235 
236  while (faces && faces->visible)
237  {
238  f = faces;
239  remove<tFace>(faces, f);
240  }
241  f = faces->next;
242  do
243  {
244  if (f->visible)
245  {
246  t = f;
247  f = f->next;
248  remove<tFace>(faces, t);
249  }
250  else
251  f = f->next;
252  } while (f != faces);
253 }
254 
256 {
257  CleanEdges();
258  CleanFaces();
259  CleanVertices(pvnext);
260 }
261 
263 {
264  tEdge e;
265  tVertex v, t;
266 
267  /* Mark all vertices incident to some undeleted edge as on the hull. */
268  e = edges;
269  do
270  {
271  e->endpts[0]->onhull = e->endpts[1]->onhull = ONHULL;
272  e = e->next;
273  } while (e != edges);
274 
275  /* Delete all vertices that have been processed but
276  are not on the hull. */
277  while (vertices && vertices->mark && !vertices->onhull)
278  {
279  /* If about to delete vnext, advance it first. */
280  v = vertices;
281  if (v == *pvnext)
282  *pvnext = v->next;
283  remove<tVertex>(vertices, v);
284  }
285  v = vertices->next;
286  do
287  {
288  if (v->mark && !v->onhull)
289  {
290  t = v;
291  v = v->next;
292  if (t == *pvnext)
293  *pvnext = t->next;
294  remove<tVertex>(vertices, t);
295  }
296  else
297  v = v->next;
298  } while (v != vertices);
299 
300  /* Reset flags. */
301  v = vertices;
302  do
303  {
304  v->duplicate = nullptr;
305  v->onhull = !ONHULL;
306  v = v->next;
307  } while (v != vertices);
308 }
309 
311 {
312  return (c->v[Z] - a->v[Z]) * (b->v[Y] - a->v[Y]) - (b->v[Z] - a->v[Z]) * (c->v[Y] - a->v[Y]) == 0 &&
313  (b->v[Z] - a->v[Z]) * (c->v[X] - a->v[X]) - (b->v[X] - a->v[X]) * (c->v[Z] - a->v[Z]) == 0 &&
314  (b->v[X] - a->v[X]) * (c->v[Y] - a->v[Y]) - (b->v[Y] - a->v[Y]) * (c->v[X] - a->v[X]) == 0;
315 }
316 
318 {
319  tEdge e;
320  int i, j;
321 
322  e = edges;
323 
324  do
325  {
326  /* find index of endpoint[0] in adjacent face[0] */
327  for (i = 0; e->adjface[0]->vertex[i] != e->endpts[0]; ++i)
328  ;
329 
330  /* find index of endpoint[0] in adjacent face[1] */
331  for (j = 0; e->adjface[1]->vertex[j] != e->endpts[0]; ++j)
332  ;
333 
334  /* check if the endpoints occur in opposite order */
335  if (!(e->adjface[0]->vertex[(i + 1) % 3] == e->adjface[1]->vertex[(j + 2) % 3] ||
336  e->adjface[0]->vertex[(i + 2) % 3] == e->adjface[1]->vertex[(j + 1) % 3]))
337  break;
338  e = e->next;
339  } while (e != edges);
340 
341  if (e != edges)
342  std::cerr << "Checks: edges are NOT consistent.\n";
343  else
344  std::cerr << "Checks: edges consistent.\n";
345 }
346 
348 {
349  tVertex v, vnext;
350 
351  v = vertices;
352  do
353  {
354  vnext = v->next;
355  if (!v->mark)
356  {
357  v->mark = PROCESSED;
358  AddOne(v);
359  CleanUp(&vnext); /* Pass down vnext in case it gets deleted. */
360 
361  if (check)
362  {
363  std::cerr << "ConstructHull: After Add of " << v->vnum << " & Cleanup:\n";
364  Checks();
365  }
366  // if ( debug )
367  // PrintOut( v );
368  }
369  v = vnext;
370  } while (v != vertices);
371 }
372 
374 {
375  tFace f;
376  tVertex v;
377  int vol;
378 
379  f = faces;
380 
381  do
382  {
383  v = vertices;
384  do
385  {
386  if (v->mark)
387  {
388  vol = VolumeSign(f, v);
389  if (vol < 0)
390  break;
391  }
392  v = v->next;
393  } while (v != vertices);
394 
395  f = f->next;
396 
397  } while (f != faces);
398 
399  if (f != faces)
400  std::cerr << "Checks: NOT convex.\n";
401  else if (check)
402  std::cerr << "Checks: convex.\n";
403 }
404 
406 {
407  tVertex v0, v1, v2, v3;
408  tFace f0, f1 = nullptr;
409  int vol;
410 
411  /* Find 3 noncollinear points. */
412  v0 = vertices;
413  while (Collinear(v0, v0->next, v0->next->next))
414  if ((v0 = v0->next) == vertices)
415  {
416  std::cout << "DoubleTriangle: All points are Collinear!\n";
417  exit(0);
418  }
419  v1 = v0->next;
420  v2 = v1->next;
421 
422  /* Mark the vertices as processed. */
423  v0->mark = PROCESSED;
424  v1->mark = PROCESSED;
425  v2->mark = PROCESSED;
426 
427  /* Create the two "twin" faces. */
428  f0 = MakeFace(v0, v1, v2, f1);
429  f1 = MakeFace(v2, v1, v0, f0);
430 
431  /* Link adjacent face fields. */
432  f0->edge[0]->adjface[1] = f1;
433  f0->edge[1]->adjface[1] = f1;
434  f0->edge[2]->adjface[1] = f1;
435  f1->edge[0]->adjface[1] = f0;
436  f1->edge[1]->adjface[1] = f0;
437  f1->edge[2]->adjface[1] = f0;
438 
439  /* Find a fourth, noncoplanar point to form tetrahedron. */
440  v3 = v2->next;
441  vol = VolumeSign(f0, v3);
442  while (!vol)
443  {
444  if ((v3 = v3->next) == v0)
445  {
446  std::cout << "DoubleTriangle: All points are coplanar!\n";
447  exit(0);
448  }
449  vol = VolumeSign(f0, v3);
450  }
451 
452  /* Insure that v3 will be the first added. */
453  vertices = v3;
454  if (debug)
455  {
456  std::cerr << "DoubleTriangle: finished. Head repositioned at v3.\n";
457  //PrintOut( vertices );
458  }
459 }
460 
462 {
463  tFace f = faces;
464  tEdge newEdge;
465  int i, j;
466 
467  do
468  {
469  for (i = 0; i < 3; i++)
470  {
471  if (!(((f->edge[i]->endpts[0] == f->vertex[i]) && (f->edge[i]->endpts[1] == f->vertex[(i + 1) % 3])) ||
472  ((f->edge[i]->endpts[1] == f->vertex[i]) && (f->edge[i]->endpts[0] == f->vertex[(i + 1) % 3]))))
473  {
474  /* Change the order of the edges on the face: */
475  for (j = 0; j < 3; j++)
476  {
477  /* find the edge that should be there */
478  if (((f->edge[j]->endpts[0] == f->vertex[i]) &&
479  (f->edge[j]->endpts[1] == f->vertex[(i + 1) % 3])) ||
480  ((f->edge[j]->endpts[1] == f->vertex[i]) && (f->edge[j]->endpts[0] == f->vertex[(i + 1) % 3])))
481  {
482  /* Swap it with the one erroneously put into its place: */
483  if (debug)
484  std::cerr << "Making a swap in EdgeOrderOnFaces: F(" << f->vertex[0]->vnum << ','
485  << f->vertex[1]->vnum << ',' << f->vertex[2]->vnum << "): e[" << i << "] and e[" << j
486  << "]\n";
487  newEdge = f->edge[i];
488  f->edge[i] = f->edge[j];
489  f->edge[j] = newEdge;
490  }
491  }
492  }
493  }
494  f = f->next;
495  } while (f != faces);
496 }
497 
499 {
500  tFace fv; /* The visible face adjacent to e */
501  int i; /* Index of e->endpoint[0] in fv. */
502  tEdge s = nullptr; /* Temporary, for swapping */
503 
504  if (e->adjface[0]->visible)
505  fv = e->adjface[0];
506  else
507  fv = e->adjface[1];
508 
509  /* Set vertex[0] & [1] of f to have the same orientation
510  as do the corresponding vertices of fv. */
511  for (i = 0; fv->vertex[i] != e->endpts[0]; ++i)
512  ;
513  /* Orient f the same as fv. */
514  if (fv->vertex[(i + 1) % 3] != e->endpts[1])
515  {
516  f->vertex[0] = e->endpts[1];
517  f->vertex[1] = e->endpts[0];
518  }
519  else
520  {
521  f->vertex[0] = e->endpts[0];
522  f->vertex[1] = e->endpts[1];
523  swap<tEdge>(s, f->edge[1], f->edge[2]);
524  }
525  /* This swap is tricky. e is edge[0]. edge[1] is based on endpt[0],
526  edge[2] on endpt[1]. So if e is oriented "forwards," we
527  need to move edge[1] to follow [0], because it precedes. */
528 
529  f->vertex[2] = p;
530 }
531 
533 {
534  tEdge new_edge[2];
535  tFace new_face;
536  int i, j;
537 
538  /* Make two new edges (if don't already exist). */
539  for (i = 0; i < 2; ++i)
540  /* If the edge exists, copy it into new_edge. */
541  if (!(new_edge[i] = e->endpts[i]->duplicate))
542  {
543  /* Otherwise (duplicate is nullptr), MakeNullEdge. */
544  new_edge[i] = MakeNullEdge();
545  new_edge[i]->endpts[0] = e->endpts[i];
546  new_edge[i]->endpts[1] = p;
547  e->endpts[i]->duplicate = new_edge[i];
548  }
549 
550  /* Make the new face. */
551  new_face = MakeNullFace();
552  new_face->edge[0] = e;
553  new_face->edge[1] = new_edge[0];
554  new_face->edge[2] = new_edge[1];
555  MakeCcw(new_face, e, p);
556 
557  /* Set the adjacent face pointers. */
558  for (i = 0; i < 2; ++i)
559  for (j = 0; j < 2; ++j)
560  /* Only one nullptr link should be set to new_face. */
561  if (!new_edge[i]->adjface[j])
562  {
563  new_edge[i]->adjface[j] = new_face;
564  break;
565  }
566 
567  return new_face;
568 }
569 
571 {
572  tFace f;
573  tEdge e0, e1, e2;
574 
575  /* Create edges of the initial triangle. */
576  if (!fold)
577  {
578  e0 = MakeNullEdge();
579  e1 = MakeNullEdge();
580  e2 = MakeNullEdge();
581  }
582  else /* Copy from fold, in reverse order. */
583  {
584  e0 = fold->edge[2];
585  e1 = fold->edge[1];
586  e2 = fold->edge[0];
587  }
588  e0->endpts[0] = v0;
589  e0->endpts[1] = v1;
590  e1->endpts[0] = v1;
591  e1->endpts[1] = v2;
592  e2->endpts[0] = v2;
593  e2->endpts[1] = v0;
594 
595  /* Create face for triangle. */
596  f = MakeNullFace();
597  f->edge[0] = e0;
598  f->edge[1] = e1;
599  f->edge[2] = e2;
600  f->vertex[0] = v0;
601  f->vertex[1] = v1;
602  f->vertex[2] = v2;
603 
604  /* Link edges to face. */
605  e0->adjface[0] = e1->adjface[0] = e2->adjface[0] = f;
606 
607  return f;
608 }
609 
610 void ConvexHull::MakeNewVertex(double x, double y, double z, int VertexId)
611 {
612  tVertex v;
613 
614  v = MakeNullVertex();
615  v->v[X] = x * ScaleFactor;
616  v->v[Y] = y * ScaleFactor;
617  v->v[Z] = z * ScaleFactor;
618  v->vnum = VertexId;
619  if ((std::abs(x) > SAFE) || (std::abs(y) > SAFE) || (std::abs(z) > SAFE))
620  {
621  std::cout << "Coordinate of vertex below might be too large: run with -d flag\n";
622  PrintPoint(v);
623  }
624 }
625 
627 {
628  tEdge e;
629 
630  e = new tsEdge;
631  e->adjface[0] = e->adjface[1] = e->newface = nullptr;
632  e->endpts[0] = e->endpts[1] = nullptr;
633  e->delete_it = !REMOVED;
634  add<tEdge>(edges, e);
635  return e;
636 }
637 
639 {
640  tFace f;
641 
642  f = new tsFace;
643  for (int i = 0; i < 3; ++i)
644  {
645  f->edge[i] = nullptr;
646  f->vertex[i] = nullptr;
647  }
648  f->visible = !VISIBLE;
649  add<tFace>(faces, f);
650  return f;
651 }
652 
654 {
655  tVertex v;
656 
657  v = new tsVertex;
658  v->duplicate = nullptr;
659  v->onhull = !ONHULL;
660  v->mark = !PROCESSED;
661  add<tVertex>(vertices, v);
662 
663  return v;
664 }
665 
667 {
668  /* Pointers to vertices, edges, faces. */
669  tVertex v;
670  tEdge e;
671  tFace f;
672  int xmin, ymin, xmax, ymax;
673  int a[3], b[3]; /* used to compute normal vector */
674  /* Counters for Euler's formula. */
675  int V = 0, E = 0, F = 0;
676  /* Note: lowercase==pointer, uppercase==counter. */
677 
678  /*-- find X min & max --*/
679  v = vertices;
680  xmin = xmax = v->v[X];
681  do
682  {
683  if (v->v[X] > xmax)
684  xmax = v->v[X];
685  else if (v->v[X] < xmin)
686  xmin = v->v[X];
687  v = v->next;
688  } while (v != vertices);
689 
690  /*-- find Y min & max --*/
691  v = vertices;
692  ymin = ymax = v->v[Y];
693  do
694  {
695  if (v->v[Y] > ymax)
696  ymax = v->v[Y];
697  else if (v->v[Y] < ymin)
698  ymin = v->v[Y];
699  v = v->next;
700  } while (v != vertices);
701 
702  /* PostScript header */
703  std::cout << "%!PS\n";
704  std::cout << "%%BoundingBox: " << xmin << ' ' << ymin << ' ' << xmax << ' ' << ymax << '\n';
705  std::cout << ".00 .00 setlinewidth\n";
706  std::cout << -xmin + 72 << ' ' << -ymin + 72 << " translate\n";
707  /* The +72 shifts the figure one inch from the lower left corner */
708 
709  /* Vertices. */
710  v = vertices;
711  do
712  {
713  if (v->mark)
714  V++;
715  v = v->next;
716  } while (v != vertices);
717  std::cout << "\n%% Vertices:\tV = " << V << '\n';
718  std::cout << "%% index:\t\tx\ty\tz\n";
719  do
720  {
721  std::cout << "%% " << v->vnum << ":\t" << v->v[X] << '\t' << v->v[Y] << '\t' << v->v[Z] << '\n';
722  v = v->next;
723  } while (v != vertices);
724 
725  /* Faces. */
726  /* visible faces are printed as PS output */
727  f = faces;
728  do
729  {
730  ++F;
731  f = f->next;
732  } while (f != faces);
733  std::cout << "\n%% Faces:\tF = " << F << '\n';
734  std::cout << "%% Visible faces only: \n";
735  do
736  {
737  /* Print face only if it is visible: if normal vector >= 0 */
738  // RFJ This code is 2d so what is calculated below
739  // is actually the perp product or signed area.
740  SubVec(f->vertex[1]->v, f->vertex[0]->v, a);
741  SubVec(f->vertex[2]->v, f->vertex[1]->v, b);
742  if ((a[0] * b[1] - a[1] * b[0]) >= 0)
743  {
744  std::cout << "%% vnums: " << f->vertex[0]->vnum << " " << f->vertex[1]->vnum << " " << f->vertex[2]->vnum
745  << '\n';
746  std::cout << "newpath\n";
747  std::cout << f->vertex[0]->v[X] << '\t' << f->vertex[0]->v[Y] << "\tmoveto\n";
748  std::cout << f->vertex[1]->v[X] << '\t' << f->vertex[1]->v[Y] << "\tlineto\n";
749  std::cout << f->vertex[2]->v[X] << '\t' << f->vertex[2]->v[Y] << "\tlineto\n";
750  std::cout << "closepath stroke\n\n";
751  }
752  f = f->next;
753  } while (f != faces);
754 
755  /* prints a list of all faces */
756  std::cout << "%% List of all faces: \n";
757  std::cout << "%%\tv0\tv1\tv2\t(vertex indices)\n";
758  do
759  {
760  std::cout << "%%\t" << f->vertex[0]->vnum << '\t' << f->vertex[1]->vnum << '\t' << f->vertex[2]->vnum << '\n';
761  f = f->next;
762  } while (f != faces);
763 
764  /* Edges. */
765  e = edges;
766  do
767  {
768  E++;
769  e = e->next;
770  } while (e != edges);
771  std::cout << "\n%% Edges:\tE = " << E << '\n';
772  /* Edges not printed out (but easily added). */
773 
774  std::cout << "\nshowpage\n\n";
775 
776  check = true;
777  CheckEuler(V, E, F);
778 }
779 
780 void ConvexHull::PrintEdges(std::ofstream &Ofile)
781 {
782  tEdge temp;
783  int i;
784 
785  temp = edges;
786  Ofile << "Edge List\n";
787  if (edges)
788  do
789  {
790  Ofile << " addr: " << std::hex << edges << '\t';
791  Ofile << "adj: ";
792  for (i = 0; i < 2; ++i)
793  Ofile << edges->adjface[i] << ' ';
794  Ofile << " endpts:" << std::dec;
795  for (i = 0; i < 2; ++i)
796  Ofile << edges->endpts[i]->vnum << ' ';
797  Ofile << " del:" << edges->delete_it << '\n';
798  edges = edges->next;
799  } while (edges != temp);
800 }
801 
802 void ConvexHull::PrintFaces(std::ofstream &Ofile)
803 {
804  int i;
805  tFace temp;
806 
807  temp = faces;
808  Ofile << "Face List\n";
809  if (faces)
810  do
811  {
812  Ofile << " addr: " << std::hex << faces << " ";
813  Ofile << " edges:" << std::hex;
814  for (i = 0; i < 3; ++i)
815  Ofile << faces->edge[i] << ' ';
816  Ofile << " vert:" << std::dec;
817  for (i = 0; i < 3; ++i)
818  Ofile << ' ' << faces->vertex[i]->vnum;
819  Ofile << " vis: " << faces->visible << '\n';
820  faces = faces->next;
821  } while (faces != temp);
822 }
823 
824 void ConvexHull::PrintObj(const char *FileName)
825 {
826  tVertex v;
827  tFace f;
828  std::map<int, int> vnumToOffsetMap;
829  int a[3], b[3]; /* used to compute normal vector */
830  double c[3], length;
831  std::ofstream Ofile;
832 
833  Ofile.open(FileName, std::ios_base::out | std::ios_base::trunc);
834 
835  Ofile << "# obj file written by chull\n";
836  Ofile << "mtllib chull.mtl\n";
837  Ofile << "g Object001\n";
838  Ofile << "s 1\n";
839  Ofile << "usemtl default\n";
840 
841  // The convex hull code removes vertices from the list that at are on the hull
842  // so the vertices list might have missing vnums. So I need to construct a map of
843  // vnums to new vertex array indices.
844  int Offset = 1;
845  v = vertices;
846  do
847  {
848  vnumToOffsetMap[v->vnum] = Offset;
849  Ofile << "v " << v->v[X] << ' ' << v->v[Y] << ' ' << v->v[Z] << '\n';
850  Offset++;
851  v = v->next;
852  } while (v != vertices);
853 
854  // normals
855  f = faces;
856  do
857  {
858  // get two tangent vectors
859  SubVec(f->vertex[1]->v, f->vertex[0]->v, a);
860  // SubVec( f->vertex[2]->v, f->vertex[1]->v, b );
861  SubVec(f->vertex[2]->v, f->vertex[0]->v, b);
862  // cross product for the normal
863  c[0] = a[1] * b[2] - a[2] * b[1];
864  c[1] = a[2] * b[0] - a[0] * b[2];
865  c[2] = a[0] * b[1] - a[1] * b[0];
866  // normalise
867  length = sqrt((c[0] * c[0]) + (c[1] * c[1]) + (c[2] * c[2]));
868  c[0] = c[0] / length;
869  c[1] = c[1] / length;
870  c[2] = c[2] / length;
871  Ofile << "vn " << c[0] << ' ' << c[1] << ' ' << c[2] << '\n';
872  f = f->next;
873  } while (f != faces);
874 
875  // Faces
876  int i = 1;
877  f = faces;
878  do
879  {
880  Ofile << "f " << vnumToOffsetMap[f->vertex[0]->vnum] << "//" << i << ' ' << vnumToOffsetMap[f->vertex[1]->vnum]
881  << "//" << i << ' ' << vnumToOffsetMap[f->vertex[2]->vnum] << "//" << i << '\n';
882  i++;
883  f = f->next;
884  } while (f != faces);
885 
886  Ofile.close();
887 
888  Ofile.open("chull.mtl", std::ios_base::out | std::ios_base::trunc);
889 
890  Ofile << "newmtl default\n";
891  Ofile << "Ka 0.2 0 0\n";
892  Ofile << "Kd 0.8 0 0\n";
893  Ofile << "illum 1\n";
894 
895  Ofile.close();
896 }
897 
898 void ConvexHull::PrintOut(const char *FileName, tVertex v)
899 {
900  std::ofstream Ofile;
901  Ofile.open(FileName, std::ios_base::out | std::ios_base::trunc);
902 
903  Ofile << "\nHead vertex " << v->vnum << " = " << std::hex << v << " :\n";
904 
905  PrintVertices(Ofile);
906  PrintEdges(Ofile);
907  PrintFaces(Ofile);
908 
909  Ofile.close();
910 }
911 
913 {
914  for (int i = 0; i < 3; i++)
915  std::cout << '\t' << p->v[i];
916 
917  std::cout << '\n';
918 }
919 
920 void ConvexHull::PrintVertices(std::ofstream &Ofile)
921 {
922  tVertex temp;
923 
924  temp = vertices;
925  Ofile << "Vertex List\n";
926  if (vertices)
927  do
928  {
929  Ofile << " addr " << std::hex << vertices << "\t";
930  Ofile << " vnum " << std::dec << vertices->vnum;
931  Ofile << '(' << vertices->v[X] << ',' << vertices->v[Y] << ',' << vertices->v[Z] << ')';
932  Ofile << " active:" << vertices->onhull;
933  Ofile << " dup:" << std::hex << vertices->duplicate;
934  Ofile << " mark:" << std::dec << vertices->mark << '\n';
936  } while (vertices != temp);
937 }
938 
940 {
941  tVertex v;
942  int x, y, z;
943  int vnum = 0;
944 
945  while (!(std::cin.eof() || std::cin.fail()))
946  {
947  std::cin >> x >> y >> z;
948  v = MakeNullVertex();
949  v->v[X] = x;
950  v->v[Y] = y;
951  v->v[Z] = z;
952  v->vnum = vnum++;
953  if ((abs(x) > SAFE) || (abs(y) > SAFE) || (abs(z) > SAFE))
954  {
955  std::cout << "Coordinate of vertex below might be too large: run with -d flag\n";
956  PrintPoint(v);
957  }
958  }
959 }
960 
962 {
963  tVertex CurrentVertex = vertices;
964  tEdge CurrentEdge = edges;
965  tFace CurrentFace = faces;
966 
967  if (nullptr != CurrentVertex)
968  {
969  do
970  {
971  tVertex TempVertex = CurrentVertex;
972  CurrentVertex = CurrentVertex->next;
973  delete TempVertex;
974  } while (CurrentVertex != vertices);
975  vertices = nullptr;
976  }
977 
978  if (nullptr != CurrentEdge)
979  {
980  do
981  {
982  tEdge TempEdge = CurrentEdge;
983  CurrentEdge = CurrentEdge->next;
984  delete TempEdge;
985  } while (CurrentEdge != edges);
986  edges = nullptr;
987  }
988 
989  if (nullptr != CurrentFace)
990  {
991  do
992  {
993  tFace TempFace = CurrentFace;
994  CurrentFace = CurrentFace->next;
995  delete TempFace;
996  } while (CurrentFace != faces);
997  faces = nullptr;
998  }
999 
1000  debug = false;
1001  check = false;
1002 }
1003 
1004 void ConvexHull::SubVec(int a[3], int b[3], int c[3])
1005 {
1006  int i;
1007 
1008  for (i = 0; i < 3; i++) // RFJ
1009  //for( i=0; i < 2; i++ )
1010  c[i] = a[i] - b[i];
1011 }
1012 
1014 {
1015  int vol;
1016  int ax, ay, az, bx, by, bz, cx, cy, cz;
1017 
1018  ax = f->vertex[0]->v[X] - p->v[X];
1019  ay = f->vertex[0]->v[Y] - p->v[Y];
1020  az = f->vertex[0]->v[Z] - p->v[Z];
1021  bx = f->vertex[1]->v[X] - p->v[X];
1022  by = f->vertex[1]->v[Y] - p->v[Y];
1023  bz = f->vertex[1]->v[Z] - p->v[Z];
1024  cx = f->vertex[2]->v[X] - p->v[X];
1025  cy = f->vertex[2]->v[Y] - p->v[Y];
1026  cz = f->vertex[2]->v[Z] - p->v[Z];
1027 
1028  vol = (ax * (by * cz - bz * cy) + ay * (bz * cx - bx * cz) + az * (bx * cy - by * cx));
1029 
1030  return vol;
1031 }
1032 
1034 {
1035  double vol;
1036  int voli;
1037  double ax, ay, az, bx, by, bz, cx, cy, cz;
1038 
1039  ax = f->vertex[0]->v[X] - p->v[X];
1040  ay = f->vertex[0]->v[Y] - p->v[Y];
1041  az = f->vertex[0]->v[Z] - p->v[Z];
1042  bx = f->vertex[1]->v[X] - p->v[X];
1043  by = f->vertex[1]->v[Y] - p->v[Y];
1044  bz = f->vertex[1]->v[Z] - p->v[Z];
1045  cx = f->vertex[2]->v[X] - p->v[X];
1046  cy = f->vertex[2]->v[Y] - p->v[Y];
1047  cz = f->vertex[2]->v[Z] - p->v[Z];
1048 
1049  vol = ax * (by * cz - bz * cy) + ay * (bz * cx - bx * cz) + az * (bx * cy - by * cx);
1050 
1051  if (debug)
1052  {
1053  /* Compute the volume using integers for comparison. */
1054  voli = Volumei(f, p);
1055  std::cerr << "Face=" << std::hex << f << "; Vertex=" << std::dec << p->vnum << ": vol(int) = " << voli
1056  << ", vol(double) = " << vol << "\n";
1057  }
1058 
1059  /* The volume should be an integer. */
1060  if (vol > 0.5)
1061  return 1;
1062  else if (vol < -0.5)
1063  return -1;
1064  else
1065  return 0;
1066 }
1067 
1068 } // namespace AlignmentSubsystem
1069 } // namespace INDI
INDI::AlignmentSubsystem::ConvexHull::Y
@ Y
Definition: ConvexHull.h:68
INDI::AlignmentSubsystem::ConvexHull::SAFE
static const int SAFE
Definition: ConvexHull.h:162
INDI::AlignmentSubsystem::ConvexHull::tFaceStructure
Definition: ConvexHull.h:146
INDI::AlignmentSubsystem::ConvexHull::Collinear
bool Collinear(tVertex a, tVertex b, tVertex c)
Collinear checks to see if the three points given are collinear, by checking to see if each element o...
Definition: ConvexHull.cpp:310
INDI::AlignmentSubsystem::ConvexHull::ReadVertices
void ReadVertices()
ReadVertices: Reads in the vertices, and links them into a circular list with MakeNullVertex....
Definition: ConvexHull.cpp:939
INDI::AlignmentSubsystem::ConvexHull::ConstructHull
void ConstructHull()
ConstructHull adds the vertices to the hull one at a time. The hull vertices are those in the list ma...
Definition: ConvexHull.cpp:347
INDI::AlignmentSubsystem::ConvexHull::PrintObj
void PrintObj(const char *FileName="chull.obj")
Outputs the faces in Lightwave obj format for 3d viewing. The files chull.obj and chull....
Definition: ConvexHull.cpp:824
INDI::AlignmentSubsystem::ConvexHull::Print
void Print()
Print: Prints out the vertices and the faces. Uses the vnum indices corresponding to the order in whi...
Definition: ConvexHull.cpp:666
INDI::AlignmentSubsystem::ConvexHull::SubVec
void SubVec(int a[3], int b[3], int c[3])
SubVec: Computes a - b and puts it into c.
Definition: ConvexHull.cpp:1004
INDI::AlignmentSubsystem::ConvexHull::Z
@ Z
Definition: ConvexHull.h:69
INDI::AlignmentSubsystem::ConvexHull::MakeNewVertex
void MakeNewVertex(double x, double y, double z, int VertexId)
Makes a vertex from the supplied information and adds it to the vertices list.
Definition: ConvexHull.cpp:610
INDI::AlignmentSubsystem::ConvexHull::MakeFace
tFace MakeFace(tVertex v0, tVertex v1, tVertex v2, tFace f)
MakeFace creates a new face structure from three vertices (in ccw order). It returns a pointer to the...
Definition: ConvexHull.cpp:570
INDI::AlignmentSubsystem::ConvexHull::MakeCcw
void MakeCcw(tFace f, tEdge e, tVertex p)
MakeCcw puts the vertices in the face structure in counterclock wise order. We want to store the vert...
Definition: ConvexHull.cpp:498
INDI::AlignmentSubsystem::ConvexHull::CheckEndpts
void CheckEndpts()
Checks that, for each face, for each i={0,1,2}, the [i]th vertex of that face is either the [0]th or ...
Definition: ConvexHull.cpp:101
INDI::AlignmentSubsystem::ConvexHull::MakeConeFace
tFace MakeConeFace(tEdge e, tVertex p)
MakeConeFace makes a new face and two new edges between the edge and the point that are passed to it....
Definition: ConvexHull.cpp:532
INDI::AlignmentSubsystem::ConvexHull::tVertexStructure::v
int v[3]
Definition: ConvexHull.h:129
INDI::AlignmentSubsystem::ConvexHull::tEdgeStructure
Definition: ConvexHull.h:137
INDI::AlignmentSubsystem::ConvexHull::Volumei
int Volumei(tFace f, tVertex p)
Volumei returns the volume of the tetrahedron determined by f and p.
Definition: ConvexHull.cpp:1013
INDI::AlignmentSubsystem::ConvexHull::tVertexStructure::duplicate
tEdge duplicate
Definition: ConvexHull.h:131
INDI::AlignmentSubsystem::ConvexHull::PROCESSED
static const bool PROCESSED
Definition: ConvexHull.h:161
INDI::AlignmentSubsystem::ConvexHull::MakeNullEdge
tEdge MakeNullEdge()
MakeNullEdge creates a new cell and initializes all pointers to NULL and sets all flags to off....
Definition: ConvexHull.cpp:626
INDI::AlignmentSubsystem::ConvexHull::edges
tEdge edges
Definition: ConvexHull.h:165
INDI::AlignmentSubsystem::ConvexHull::tVertexStructure
Definition: ConvexHull.h:127
INDI::AlignmentSubsystem::ConvexHull::ONHULL
static const bool ONHULL
Definition: ConvexHull.h:158
INDI::AlignmentSubsystem::ConvexHull::vertices
tVertex vertices
Definition: ConvexHull.h:164
INDI::AlignmentSubsystem::ConvexHull::tEdgeStructure::next
tEdge next
Definition: ConvexHull.h:143
INDI::AlignmentSubsystem::ConvexHull::ConvexHull
ConvexHull()
Definition: ConvexHull.cpp:44
INDI::AlignmentSubsystem::ConvexHull::Reset
void Reset()
Frees the vertices edges and faces lists and resets the debug and check flags.
Definition: ConvexHull.cpp:961
INDI::AlignmentSubsystem::ConvexHull::AddOne
bool AddOne(tVertex p)
AddOne is passed a vertex. It first determines all faces visible from that point. If none are visible...
Definition: ConvexHull.cpp:49
INDI::AlignmentSubsystem::ConvexHull::VolumeSign
int VolumeSign(tFace f, tVertex p)
VolumeSign returns the sign of the volume of the tetrahedron determined by f and p....
Definition: ConvexHull.cpp:1033
INDI::AlignmentSubsystem::ConvexHull::tVertexStructure::vnum
int vnum
Definition: ConvexHull.h:130
INDI::AlignmentSubsystem::ConvexHull::CleanUp
void CleanUp(tVertex *pvnext)
CleanUp goes through each data structure list and clears all flags and NULLs out some pointers....
Definition: ConvexHull.cpp:255
INDI::AlignmentSubsystem::ConvexHull::MakeNullFace
tFace MakeNullFace()
MakeNullFace creates a new face structure and initializes all of its flags to NULL and sets all the f...
Definition: ConvexHull.cpp:638
dec
double dec
Definition: ieqprolegacydriver.cpp:44
INDI::AlignmentSubsystem::ConvexHull::EdgeOrderOnFaces
void EdgeOrderOnFaces()
EdgeOrderOnFaces: puts e0 between v0 and v1, e1 between v1 and v2, e2 between v2 and v0 on each face....
Definition: ConvexHull.cpp:461
INDI::AlignmentSubsystem::ConvexHull::CleanVertices
void CleanVertices(tVertex *pvnext)
CleanVertices runs through the vertex list and deletes the vertices that are marked as processed but ...
Definition: ConvexHull.cpp:262
INDI::AlignmentSubsystem::ConvexHull::tVertexStructure::next
tVertex next
Definition: ConvexHull.h:134
INDI::AlignmentSubsystem::ConvexHull::X
@ X
Definition: ConvexHull.h:67
INDI::AlignmentSubsystem::ConvexHull::CleanFaces
void CleanFaces()
CleanFaces runs through the face list and deletes any face marked visible.
Definition: ConvexHull.cpp:231
INDI::AlignmentSubsystem::ConvexHull::Checks
void Checks()
Checks the consistency of the hull and prints the results to the standard error output.
Definition: ConvexHull.cpp:157
INDI::AlignmentSubsystem::ConvexHull::tsEdge
struct tEdgeStructure tsEdge
Definition: ConvexHull.h:121
INDI::AlignmentSubsystem::ConvexHull::PrintPoint
void PrintPoint(tVertex p)
Prints a single vertex to the standard output.
Definition: ConvexHull.cpp:912
INDI::AlignmentSubsystem::ConvexHull::VISIBLE
static const bool VISIBLE
Definition: ConvexHull.h:160
INDI::AlignmentSubsystem::ConvexHull::tEdgeStructure::endpts
tVertex endpts[2]
Definition: ConvexHull.h:140
INDI::AlignmentSubsystem::ConvexHull::tEdgeStructure::adjface
tFace adjface[2]
Definition: ConvexHull.h:139
INDI::AlignmentSubsystem::ConvexHull::DoubleTriangle
void DoubleTriangle()
DoubleTriangle builds the initial double triangle. It first finds 3 noncollinear points and makes two...
Definition: ConvexHull.cpp:405
INDI::AlignmentSubsystem::ConvexHull::CheckEuler
void CheckEuler(int V, int E, int F)
CheckEuler checks Euler's relation, as well as its implications when all faces are known to be triang...
Definition: ConvexHull.cpp:136
INDI::AlignmentSubsystem::ConvexHull::PrintEdges
void PrintEdges(std::ofstream &Ofile)
Prints the edges Ofile.
Definition: ConvexHull.cpp:780
INDI::AlignmentSubsystem::ConvexHull::PrintVertices
void PrintVertices(std::ofstream &Ofile)
Prints vertices to Ofile.
Definition: ConvexHull.cpp:920
INDI::AlignmentSubsystem::ConvexHull::faces
tFace faces
Definition: ConvexHull.h:166
INDI::AlignmentSubsystem::ConvexHull::tVertexStructure::mark
bool mark
Definition: ConvexHull.h:133
INDI::AlignmentSubsystem::ConvexHull::tFaceStructure::edge
tEdge edge[3]
Definition: ConvexHull.h:150
INDI::AlignmentSubsystem::ConvexHull::Consistency
void Consistency()
Consistency runs through the edge list and checks that all adjacent faces have their endpoints in opp...
Definition: ConvexHull.cpp:317
INDI::AlignmentSubsystem::ConvexHull::PrintFaces
void PrintFaces(std::ofstream &Ofile)
Prints the faces to Ofile.
Definition: ConvexHull.cpp:802
INDI::AlignmentSubsystem::ConvexHull::tFaceStructure::next
tFace next
Definition: ConvexHull.h:153
INDI::AlignmentSubsystem::ConvexHull::tVertexStructure::onhull
bool onhull
Definition: ConvexHull.h:132
INDI
Namespace to encapsulate INDI client, drivers, and mediator classes.
Definition: AlignmentSubsystemForClients.cpp:11
INDI::AlignmentSubsystem::ConvexHull::CleanEdges
void CleanEdges()
CleanEdges runs through the edge list and cleans up the structure. If there is a newface then it will...
Definition: ConvexHull.cpp:190
INDI::AlignmentSubsystem::ConvexHull::tEdgeStructure::delete_it
bool delete_it
Definition: ConvexHull.h:142
INDI::AlignmentSubsystem::ConvexHull::Convexity
void Convexity()
Convexity checks that the volume between every face and every point is negative. This shows that each...
Definition: ConvexHull.cpp:373
INDI::AlignmentSubsystem::ConvexHull::PrintOut
void PrintOut(const char *FileName, tVertex v)
Prints vertices, edges and faces to the standard error output.
Definition: ConvexHull.cpp:898
INDI::AlignmentSubsystem::ConvexHull::tEdgeStructure::newface
tFace newface
Definition: ConvexHull.h:141
INDI::AlignmentSubsystem::ConvexHull::tsFace
struct tFaceStructure tsFace
Definition: ConvexHull.h:124
INDI::AlignmentSubsystem::ConvexHull::REMOVED
static const bool REMOVED
Definition: ConvexHull.h:159
INDI::AlignmentSubsystem::ConvexHull::MakeNullVertex
tVertex MakeNullVertex()
MakeNullVertex: Makes a vertex, nulls out fields.
Definition: ConvexHull.cpp:653
ConvexHull.h
INDI::AlignmentSubsystem::ConvexHull::tFaceStructure::vertex
tVertex vertex[3]
Definition: ConvexHull.h:151
INDI::AlignmentSubsystem::ConvexHull::tFaceStructure::visible
bool visible
Definition: ConvexHull.h:152
INDI::AlignmentSubsystem::ConvexHull::tsVertex
struct tVertexStructure tsVertex
Definition: ConvexHull.h:118