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