1 /* quickhull, Copyright (c) 2016 Karim Naaji, karim.naaji@gmail.com
2 https://github.com/karimnaaji/3d-quickhull
7 Copyright (c) 2016 Karim Naaji, karim.naaji@gmail.com
9 Permission is hereby granted, free of charge, to any person obtaining a copy
10 of this software and associated documentation files (the "Software"), to deal
11 in the Software without restriction, including without limitation the rights
12 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
13 copies of the Software, and to permit persons to whom the Software is
14 furnished to do so, subject to the following conditions:
16 The above copyright notice and this permission notice shall be included in
17 all copies or substantial portions of the Software.
19 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
20 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
22 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
24 OUT OF OR IN CONNECTION WITH THE
27 [1] http://box2d.org/files/GDC2014/DirkGregorius_ImplementingQuickHull.pdf
28 [2] http://www.cs.smith.edu/~orourke/books/compgeom.html
29 [3] http://www.flipcode.com/archives/The_Half-Edge_Data_Structure.shtml
30 [4] http://doc.cgal.org/latest/HalfedgeDS/index.html
31 [5] http://thomasdiewald.com/blog/?p=1888
32 [6] https://fgiesen.wordpress.com/2012/02/21/half-edge-based-mesh-representations-theory/
35 #define QUICKHULL_DEBUG // Only if assertions need to be checked
36 #include "quickhull.h"
39 - 25-Feb-2018: jwz: adapted for xscreensaver
40 - 1.0.1 (2016-11-01): Various improvements over epsilon issues and
42 Debug functionalities to test final results dynamically
43 API to export hull meshes in OBJ files
44 - 1.0 (2016-09-10): Initial
47 - use float* from public interface
53 #endif /* HAVE_CONFIG_H */
55 extern const char *progname;
57 #include "quickhull.h"
59 #include "screenhackI.h" /* for jwxyz_abort */
61 #include <math.h> /* sqrt & fabs */
62 #include <stdio.h> /* FILE */
63 #include <string.h> /* memcpy */
65 /* Quickhull helpers, define your own if needed */
66 #ifndef QUICKHULL_HELPERS
67 #include <stdlib.h> /* malloc, free, realloc */
68 #define QUICKHULL_HELPERS 1
69 #define QH_MALLOC(T, N) ((T*) malloc(N * sizeof(T)))
70 #define QH_REALLOC(T, P, N) ((T*)realloc(P, sizeof(T) * N))
71 #define QH_FREE(T) free(T)
72 #define QH_SWAP(T, A, B) { T tmp = B; B = A; A = tmp; }
73 #ifdef QUICKHULL_DEBUG
74 #define QH_ASSERT(STMT) if (!(STMT)) { *(int *)0 = 0; }
75 #define QH_LOG(FMT, ...) printf(FMT, ## __VA_ARGS__)
77 #define QH_ASSERT(STMT)
78 #define QH_LOG(FMT, ...)
79 #endif /* QUICKHULL_DEBUG */
80 #endif /* QUICKHULL_HELPERS */
83 #define QH_FLT_MAX 1e+37F
87 #define QH_FLT_EPS 1E-5F
90 #ifndef QH_VERTEX_SET_SIZE
91 #define QH_VERTEX_SET_SIZE 128
94 typedef long qh_index_t;
96 typedef struct qh_half_edge {
97 qh_index_t opposite_he; /* index of the opposite half edge */
98 qh_index_t next_he; /* index of the next half edge */
99 qh_index_t previous_he; /* index of the previous half edge */
100 qh_index_t he; /* index of the current half edge */
101 qh_index_t to_vertex; /* index of the next vertex */
102 qh_index_t adjacent_face; /* index of the ajacent face */
105 typedef struct qh_index_set {
108 unsigned int capacity;
111 typedef struct qh_face {
114 qh_vertex_t centroid;
121 typedef struct qh_index_stack {
126 typedef struct qh_context {
128 qh_half_edge_t* edges;
129 qh_vertex_t* vertices;
130 qh_vertex_t centroid;
131 qh_index_stack_t facestack;
132 qh_index_stack_t scratch;
133 qh_index_stack_t horizonedges;
134 qh_index_stack_t newhorizonedges;
137 unsigned int nvertices;
140 #ifdef QUICKHULL_DEBUG
141 unsigned int maxfaces;
142 unsigned int maxedges;
147 qh__find_6eps(qh_vertex_t* vertices, unsigned int nvertices, qh_index_t* eps)
149 qh_vertex_t* ptr = vertices;
151 double minxy = +QH_FLT_MAX;
152 double minxz = +QH_FLT_MAX;
153 double minyz = +QH_FLT_MAX;
155 double maxxy = -QH_FLT_MAX;
156 double maxxz = -QH_FLT_MAX;
157 double maxyz = -QH_FLT_MAX;
160 for (i = 0; i < 6; ++i) {
164 for (i = 0; i < nvertices; ++i) {
165 if (ptr->z < minxy) {
169 if (ptr->y < minxz) {
173 if (ptr->x < minyz) {
177 if (ptr->z > maxxy) {
181 if (ptr->y > maxxz) {
185 if (ptr->x > maxyz) {
194 qh__vertex_segment_length2(qh_vertex_t* p, qh_vertex_t* a, qh_vertex_t* b)
196 double dx = b->x - a->x;
197 double dy = b->y - a->y;
198 double dz = b->z - a->z;
200 double d = dx * dx + dy * dy + dz * dz;
207 double t = ((p->x - a->x) * dx +
209 (p->z - a->z) * dz) / d;
226 return dx * dx + dy * dy + dz * dz;
230 qh__vec3_sub(qh_vec3_t* a, qh_vec3_t* b)
238 qh__vec3_add(qh_vec3_t* a, qh_vec3_t* b)
246 qh__vec3_multiply(qh_vec3_t* a, double v)
254 qh__vertex_equals_epsilon(qh_vertex_t* a, qh_vertex_t* b, double epsilon)
256 return fabs(a->x - b->x) <= epsilon &&
257 fabs(a->y - b->y) <= epsilon &&
258 fabs(a->z - b->z) <= epsilon;
262 qh__vec3_length2(qh_vec3_t* v)
264 return v->x * v->x + v->y * v->y + v->z * v->z;
268 qh__vec3_dot(qh_vec3_t* v1, qh_vec3_t* v2)
270 return v1->x * v2->x + v1->y * v2->y + v1->z * v2->z;
274 qh__vec3_normalize(qh_vec3_t* v)
276 qh__vec3_multiply(v, 1.f / sqrt(qh__vec3_length2(v)));
280 qh__find_2dps_6eps(qh_vertex_t* vertices, qh_index_t* eps,
284 double max = -QH_FLT_MAX;
286 for (i = 0; i < 6; ++i) {
287 for (j = 0; j < 6; ++j) {
295 d = vertices[eps[i]];
296 qh__vec3_sub(&d, &vertices[eps[j]]);
297 d2 = qh__vec3_length2(&d);
309 qh__vec3_cross(qh_vec3_t* v1, qh_vec3_t* v2)
313 cross.x = v1->y * v2->z - v1->z * v2->y;
314 cross.y = v1->z * v2->x - v1->x * v2->z;
315 cross.z = v1->x * v2->y - v1->y * v2->x;
321 qh__face_centroid(qh_index_t vertices[3], qh_context_t* context)
323 qh_vertex_t centroid;
326 centroid.x = centroid.y = centroid.z = 0.0;
327 for (i = 0; i < 3; ++i) {
328 qh__vec3_add(¢roid, context->vertices + vertices[i]);
331 qh__vec3_multiply(¢roid, 1.0 / 3.0);
337 qh__dist_point_plane(qh_vertex_t* v, qh_vec3_t* normal, double sdist)
339 return fabs(qh__vec3_dot(v, normal) - sdist);
343 qh__init_half_edge(qh_half_edge_t* half_edge) {
344 half_edge->adjacent_face = -1;
346 half_edge->next_he = -1;
347 half_edge->opposite_he = -1;
348 half_edge->to_vertex = -1;
349 half_edge->previous_he = -1;
352 static qh_half_edge_t*
353 qh__next_edge(qh_context_t* context)
355 qh_half_edge_t* edge = context->edges + context->nedges;
356 qh__init_half_edge(edge);
358 edge->he = context->nedges;
361 QH_ASSERT(context->nedges < context->maxedges);
367 qh__next_face(qh_context_t* context)
369 qh_face_t* face = context->faces + context->nfaces;
371 face->face = context->nfaces;
372 face->iset.indices = NULL;
373 context->valid[context->nfaces] = 1;
376 QH_ASSERT(context->nfaces < context->maxfaces);
382 qh__edge_vec3(qh_half_edge_t* edge, qh_context_t* context)
384 qh_half_edge_t prevhe = context->edges[edge->previous_he];
387 v0 = context->vertices[prevhe.to_vertex];
388 v1 = context->vertices[edge->to_vertex];
390 qh__vec3_sub(&v1, &v0);
391 qh__vec3_normalize(&v1);
397 qh__face_init(qh_face_t* face, qh_index_t vertices[3],
398 qh_context_t* context)
400 qh_half_edge_t* e0 = qh__next_edge(context);
401 qh_half_edge_t* e1 = qh__next_edge(context);
402 qh_half_edge_t* e2 = qh__next_edge(context);
404 qh_vertex_t centroid, normal;
406 e2->to_vertex = vertices[0];
407 e0->to_vertex = vertices[1];
408 e1->to_vertex = vertices[2];
410 e0->next_he = e1->he;
411 e2->previous_he = e1->he;
412 face->edges[1] = e1->he;
414 e1->next_he = e2->he;
415 e0->previous_he = e2->he;
416 face->edges[2] = e2->he;
417 v1 = qh__edge_vec3(e2, context);
419 e2->next_he = e0->he;
420 e1->previous_he = e0->he;
421 face->edges[0] = e0->he;
422 v0 = qh__edge_vec3(e0, context);
424 e2->adjacent_face = face->face;
425 e1->adjacent_face = face->face;
426 e0->adjacent_face = face->face;
428 qh__vec3_multiply(&v1, -1.f);
429 normal = qh__vec3_cross(&v0, &v1);
431 qh__vec3_normalize(&normal);
432 centroid = qh__face_centroid(vertices, context);
433 face->centroid = centroid;
434 face->sdist = qh__vec3_dot(&normal, ¢roid);
435 face->normal = normal;
436 face->iset.indices = QH_MALLOC(qh_index_t, QH_VERTEX_SET_SIZE);
437 face->iset.capacity = QH_VERTEX_SET_SIZE;
439 face->visitededges = 0;
443 qh__tetrahedron_basis(qh_context_t* context, qh_index_t vertices[3])
446 int i, j = 0, k = 0, l = 0;
447 double max = -QH_FLT_MAX;
449 qh__find_6eps(context->vertices, context->nvertices, eps);
450 qh__find_2dps_6eps(context->vertices, eps, &j, &k);
452 for (i = 0; i < 6; ++i) {
455 if (i == j || i == k) {
459 d2 = qh__vertex_segment_length2(context->vertices + eps[i],
460 context->vertices + eps[j],
461 context->vertices + eps[k]);
469 vertices[0] = eps[j];
470 vertices[1] = eps[k];
471 vertices[2] = eps[l];
475 qh__push_stack(qh_index_stack_t* stack, qh_index_t index)
477 stack->begin[stack->size] = index;
482 qh__pop_stack(qh_index_stack_t* stack)
486 if (stack->size > 0) {
487 top = stack->begin[stack->size - 1];
495 qh__furthest_point_from_plane(qh_context_t* context,
502 double max = -QH_FLT_MAX;
504 for (i = 0; i < nindices; ++i) {
505 qh_index_t index = indices ? *(indices + i) : i;
506 double dist = qh__dist_point_plane(context->vertices + index,
519 qh__face_can_see_vertex(qh_face_t* face, qh_vertex_t* v)
523 qh__vec3_sub(&tov, &face->centroid);
524 return qh__vec3_dot(&tov, &face->normal) > 0;
528 qh__face_can_see_vertex_epsilon(qh_context_t* context, qh_face_t* face,
529 qh_vertex_t* v, double epsilon)
534 qh__vec3_sub(&tov, &face->centroid);
535 dot = qh__vec3_dot(&tov, &face->normal);
542 if (dot <= epsilon && dot >= 0) {
543 qh_vec3_t n = face->normal;
545 /* allow epsilon degeneration along the face normal */
546 qh__vec3_multiply(&n, epsilon);
557 qh__assert_half_edge(qh_half_edge_t* edge, qh_context_t* context)
559 QH_ASSERT(edge->opposite_he != -1);
560 QH_ASSERT(edge->he != -1);
561 QH_ASSERT(edge->adjacent_face != -1);
562 QH_ASSERT(edge->next_he != -1);
563 QH_ASSERT(edge->previous_he != -1);
564 QH_ASSERT(edge->to_vertex != -1);
565 QH_ASSERT(context->edges[edge->opposite_he].to_vertex != edge->to_vertex);
569 qh__assert_face(qh_face_t* face, qh_context_t* context)
573 for (i = 0; i < 3; ++i) {
574 qh__assert_half_edge(context->edges + face->edges[i], context);
577 QH_ASSERT(context->valid[face->face]);
580 #ifdef QUICKHULL_DEBUG
583 qh__log_face(qh_context_t* context, qh_face_t const* face) {
584 QH_LOG("Face %ld:\n", face->face);
585 for (int i = 0; i < 3; ++i) {
586 qh_half_edge_t edge = context->edges[face->edges[i]];
587 QH_LOG("\te%d %ld\n", i, edge.he);
588 QH_LOG("\t\te%d.opposite_he %ld\n", i, edge.opposite_he);
589 QH_LOG("\t\te%d.next_he %ld\n", i, edge.next_he);
590 QH_LOG("\t\te%d.previous_he %ld\n", i, edge.previous_he);
591 QH_LOG("\t\te%d.to_vertex %ld\n", i, edge.to_vertex);
592 QH_LOG("\t\te%d.adjacent_face %ld\n", i, edge.adjacent_face);
594 QH_LOG("\tnormal %f %f %f\n", face->normal.x, face->normal.y,
596 QH_LOG("\tsdist %f\n", face->sdist);
597 QH_LOG("\tcentroid %f %f %f\n", face->centroid.x, face->centroid.y,
602 qh__test_hull(qh_context_t* context, double epsilon, int testiset)
604 unsigned int i, j, k;
606 for (i = 0; i < context->nvertices; ++i) {
607 qh_index_t vindex = i;
610 for (j = 0; j < context->nfaces; ++j) {
611 if (!context->valid[j]) {
614 qh_face_t* face = context->faces + j;
616 qh_half_edge_t* e0 = context->edges + face->edges[0];
617 qh_half_edge_t* e1 = context->edges + face->edges[1];
618 qh_half_edge_t* e2 = context->edges + face->edges[2];
620 if (e0->to_vertex == vindex ||
621 e1->to_vertex == vindex ||
622 e2->to_vertex == vindex) {
628 for (k = 0; k < face->iset.size; ++k) {
629 if (vindex == face->iset.indices[k]) {
640 for (j = 0; j < context->nfaces; ++j) {
641 if (!context->valid[j]) {
644 qh_face_t* face = context->faces + j;
646 qh_vertex_t vertex = context->vertices[vindex];
647 qh__vec3_sub(&vertex, &face->centroid);
648 if (qh__vec3_dot(&face->normal, &vertex) > epsilon) {
649 #ifdef QUICKHULL_DEBUG
650 qh__log_face(context, face);
659 #endif /* QUICKHULL_DEBUG */
663 #ifdef QUICKHULL_DEBUG
664 qh__build_hull(qh_context_t* context, double epsilon, unsigned int step,
665 unsigned int* failurestep)
667 qh__build_hull(qh_context_t* context, double epsilon)
670 qh_index_t topface = qh__pop_stack(&context->facestack);
673 #ifdef QUICKHULL_DEBUG
674 unsigned int iteration = 0;
677 while (topface != -1) {
678 qh_face_t* face = context->faces + topface;
679 qh_index_t fvi, apex;
683 #ifdef QUICKHULL_DEBUG
684 if (!context->valid[topface] || face->iset.size == 0 || iteration == step)
686 if (!context->valid[topface] || face->iset.size == 0)
689 topface = qh__pop_stack(&context->facestack);
693 #ifdef QUICKHULL_DEBUG
694 if (failurestep != NULL && !qh__test_hull(context, epsilon, 1)) {
695 if (*failurestep == 0) {
696 *failurestep = iteration;
704 fvi = qh__furthest_point_from_plane(context, face->iset.indices,
705 face->iset.size, &face->normal, face->sdist);
706 fv = context->vertices + *(face->iset.indices + fvi);
708 qh__assert_face(face, context);
710 /* Reset visited flag for faces */
712 for (i = 0; i < context->nfaces; ++i) {
713 context->faces[i].visitededges = 0;
717 /* Find horizon edge */
719 qh_index_t tovisit = topface;
720 qh_face_t* facetovisit = context->faces + tovisit;
722 /* Release scratch */
723 context->scratch.size = 0;
725 while (tovisit != -1) {
726 if (facetovisit->visitededges >= 3) {
727 context->valid[tovisit] = 0;
728 tovisit = qh__pop_stack(&context->scratch);
729 facetovisit = context->faces + tovisit;
731 qh_index_t edgeindex = facetovisit->edges[facetovisit->visitededges];
732 qh_half_edge_t* edge;
733 qh_half_edge_t* oppedge;
736 facetovisit->visitededges++;
738 edge = context->edges + edgeindex;
739 oppedge = context->edges + edge->opposite_he;
740 adjface = context->faces + oppedge->adjacent_face;
742 if (!context->valid[oppedge->adjacent_face]) { continue; }
744 qh__assert_half_edge(oppedge, context);
745 qh__assert_half_edge(edge, context);
746 qh__assert_face(adjface, context);
748 if (!qh__face_can_see_vertex(adjface, fv)) {
749 qh__push_stack(&context->horizonedges, edge->he);
751 context->valid[tovisit] = 0;
752 qh__push_stack(&context->scratch, adjface->face);
758 apex = face->iset.indices[fvi];
760 /* Sort horizon edges in CCW order */
762 qh_vertex_t triangle[3];
764 qh_vec3_t v0, v1, toapex;
767 for (i = 0; i < context->horizonedges.size; ++i) {
768 qh_index_t he0 = context->horizonedges.begin[i];
769 qh_index_t he0vert = context->edges[he0].to_vertex;
770 qh_index_t phe0 = context->edges[he0].previous_he;
771 qh_index_t phe0vert = context->edges[phe0].to_vertex;
773 for (j = i + 2; j < context->horizonedges.size; ++j) {
774 qh_index_t he1 = context->horizonedges.begin[j];
775 qh_index_t he1vert = context->edges[he1].to_vertex;
776 qh_index_t phe1 = context->edges[he1].previous_he;
777 qh_index_t phe1vert = context->edges[phe1].to_vertex;
779 if (phe1vert == he0vert || phe0vert == he1vert) {
780 QH_SWAP(qh_index_t, context->horizonedges.begin[j],
781 context->horizonedges.begin[i + 1]);
788 context->vertices[context->edges[he0].to_vertex];
793 /* Detect first triangle face ordering */
797 qh__vec3_sub(&v0, &triangle[1]);
798 qh__vec3_sub(&v1, &triangle[1]);
800 n = qh__vec3_cross(&v0, &v1);
802 /* Get the vector to the apex */
803 toapex = triangle[0];
805 qh__vec3_sub(&toapex, context->vertices + apex);
807 reversed = qh__vec3_dot(&n, &toapex) < 0.f;
811 /* Create new faces */
813 qh_index_t top = qh__pop_stack(&context->horizonedges);
814 qh_index_t last = qh__pop_stack(&context->horizonedges);
815 qh_index_t first = top;
818 QH_ASSERT(context->newhorizonedges.size == 0);
820 /* Release scratch */
821 context->scratch.size = 0;
824 qh_half_edge_t* prevhe;
825 qh_half_edge_t* nexthe;
826 qh_half_edge_t* oppedge;
827 /* qh_vec3_t normal; */
828 /* qh_vertex_t fcentroid; */
837 prevhe = context->edges + last;
838 nexthe = context->edges + top;
841 QH_SWAP(qh_half_edge_t*, prevhe, nexthe);
844 verts[0] = prevhe->to_vertex;
845 verts[1] = nexthe->to_vertex;
848 context->valid[nexthe->adjacent_face] = 0;
850 oppedge = context->edges + nexthe->opposite_he;
851 newface = qh__next_face(context);
853 qh__face_init(newface, verts, context);
855 oppedge->opposite_he = context->edges[newface->edges[0]].he;
856 context->edges[newface->edges[0]].opposite_he = oppedge->he;
858 qh__push_stack(&context->scratch, newface->face);
859 qh__push_stack(&context->newhorizonedges, newface->edges[0]);
862 last = qh__pop_stack(&context->horizonedges);
866 /* Attach point sets to newly created faces */
868 for (k = 0; k < context->nfaces; ++k) {
869 qh_face_t* f = context->faces + k;
871 if (context->valid[k] || f->iset.size == 0) {
875 if (f->visitededges == 3) {
876 context->valid[k] = 0;
879 for (i = 0; i < f->iset.size; ++i) {
880 qh_index_t vertex = f->iset.indices[i];
881 /* qh_vertex_t* v = context->vertices + vertex; */
882 qh_face_t* dface = NULL;
884 for (j = 0; j < context->scratch.size; ++j) {
885 qh_face_t* newface = context->faces + context->scratch.begin[j];
886 qh_half_edge_t* e0 = context->edges + newface->edges[0];
887 qh_half_edge_t* e1 = context->edges + newface->edges[1];
888 qh_half_edge_t* e2 = context->edges + newface->edges[2];
889 /* qh_vertex_t cv; */
891 if (e0->to_vertex == vertex ||
892 e1->to_vertex == vertex ||
893 e2->to_vertex == vertex) {
897 if (qh__face_can_see_vertex_epsilon(context, newface,
898 context->vertices + vertex,
906 if (dface->iset.size + 1 >= dface->iset.capacity) {
907 dface->iset.capacity *= 2;
908 dface->iset.indices = QH_REALLOC(qh_index_t,
909 dface->iset.indices, dface->iset.capacity);
912 dface->iset.indices[dface->iset.size++] = vertex;
920 /* Link new faces together */
922 for (i = 0; i < context->newhorizonedges.size; ++i) {
923 qh_index_t phe0, nhe1;
929 ii = (i == 0) ? context->newhorizonedges.size - 1 : (i-1);
931 ii = (i+1) % context->newhorizonedges.size;
934 phe0 = context->edges[context->newhorizonedges.begin[i]].previous_he;
935 nhe1 = context->edges[context->newhorizonedges.begin[ii]].next_he;
937 he0 = context->edges + phe0;
938 he1 = context->edges + nhe1;
940 QH_ASSERT(he1->to_vertex == apex);
941 QH_ASSERT(he0->opposite_he == -1);
942 QH_ASSERT(he1->opposite_he == -1);
944 he0->opposite_he = he1->he;
945 he1->opposite_he = he0->he;
948 context->newhorizonedges.size = 0;
951 /* Push new face to stack */
953 for (i = 0; i < context->scratch.size; ++i) {
954 qh_face_t* face = context->faces + context->scratch.begin[i];
956 if (face->iset.size > 0) {
957 qh__push_stack(&context->facestack, face->face);
961 /* Release scratch */
962 context->scratch.size = 0;
965 topface = qh__pop_stack(&context->facestack);
967 /* TODO: push all non-valid faces for reuse in face stack memory pool */
971 #ifdef QUICKHULL_FILES
973 qh_mesh_export(qh_mesh_t const* mesh, char const* filename)
975 FILE* objfile = fopen(filename, "wt");
976 fprintf(objfile, "o\n");
978 for (int i = 0; i < mesh->nvertices; ++i) {
979 qh_vertex_t v = mesh->vertices[i];
980 fprintf(objfile, "v %f %f %f\n", v.x, v.y, v.z);
983 for (int i = 0; i < mesh->nnormals; ++i) {
984 qh_vec3_t n = mesh->normals[i];
985 fprintf(objfile, "vn %f %f %f\n", n.x, n.y, n.z);
988 for (int i = 0, j = 0; i < mesh->nindices; i += 3, j++) {
989 fprintf(objfile, "f %u/%u %u/%u %u/%u\n",
990 mesh->indices[i+0] + 1, mesh->normalindices[j] + 1,
991 mesh->indices[i+1] + 1, mesh->normalindices[j] + 1,
992 mesh->indices[i+2] + 1, mesh->normalindices[j] + 1);
997 #endif /* QUICKHULL_FILES */
1000 qh__build_tetrahedron(qh_context_t* context, double epsilon)
1003 qh_index_t vertices[3];
1006 qh_vertex_t normal, centroid, vapex, tcentroid;
1008 /* Get the initial tetrahedron basis (first face) */
1009 qh__tetrahedron_basis(context, &vertices[0]);
1011 /* Find apex from the tetrahedron basis */
1016 v0 = context->vertices[vertices[1]];
1017 v1 = context->vertices[vertices[2]];
1019 qh__vec3_sub(&v0, context->vertices + vertices[0]);
1020 qh__vec3_sub(&v1, context->vertices + vertices[0]);
1022 normal = qh__vec3_cross(&v0, &v1);
1023 qh__vec3_normalize(&normal);
1025 centroid = qh__face_centroid(vertices, context);
1026 sdist = qh__vec3_dot(&normal, ¢roid);
1028 apex = qh__furthest_point_from_plane(context, NULL,
1029 context->nvertices, &normal, sdist);
1030 vapex = context->vertices[apex];
1032 qh__vec3_sub(&vapex, ¢roid);
1034 /* Whether the face is looking towards the apex */
1035 if (qh__vec3_dot(&vapex, &normal) > 0) {
1036 QH_SWAP(qh_index_t, vertices[1], vertices[2]);
1040 faces = qh__next_face(context);
1041 qh__face_init(&faces[0], vertices, context);
1043 /* Build faces from the tetrahedron basis to the apex */
1045 qh_index_t facevertices[3];
1046 for (i = 0; i < 3; ++i) {
1047 qh_half_edge_t* edge = context->edges + faces[0].edges[i];
1048 qh_half_edge_t prevedge = context->edges[edge->previous_he];
1049 qh_face_t* face = faces+i+1;
1052 facevertices[0] = edge->to_vertex;
1053 facevertices[1] = prevedge.to_vertex;
1054 facevertices[2] = apex;
1056 qh__next_face(context);
1057 qh__face_init(face, facevertices, context);
1059 e0 = context->edges + faces[i+1].edges[0];
1060 edge->opposite_he = e0->he;
1061 e0->opposite_he = edge->he;
1065 /* Attach half edges to faces tied to the apex */
1067 for (i = 0; i < 3; ++i) {
1069 qh_face_t* nextface;
1076 nextface = faces+j+1;
1078 e1 = context->edges + face->edges[1];
1079 e2 = context->edges + nextface->edges[2];
1081 QH_ASSERT(e1->opposite_he == -1);
1082 QH_ASSERT(e2->opposite_he == -1);
1084 e1->opposite_he = e2->he;
1085 e2->opposite_he = e1->he;
1087 qh__assert_half_edge(e1, context);
1088 qh__assert_half_edge(e2, context);
1092 /* Create initial point set; every point is */
1093 /* attached to the first face it can see */
1095 for (i = 0; i < context->nvertices; ++i) {
1096 /* qh_vertex_t* v; */
1097 qh_face_t* dface = NULL;
1099 if (vertices[0] == i || vertices[1] == i || vertices[2] == i) {
1103 for (j = 0; j < 4; ++j) {
1104 if (qh__face_can_see_vertex_epsilon(context, context->faces + j,
1105 context->vertices + i, epsilon)) {
1106 dface = context->faces + j;
1115 for (j = 0; j < 3; ++j) {
1116 qh_half_edge_t* e = context->edges + dface->edges[j];
1117 if (i == e->to_vertex) {
1123 if (!valid) { continue; }
1125 if (dface->iset.size + 1 >= dface->iset.capacity) {
1126 dface->iset.capacity *= 2;
1127 dface->iset.indices = QH_REALLOC(qh_index_t,
1128 dface->iset.indices, dface->iset.capacity);
1131 dface->iset.indices[dface->iset.size++] = i;
1136 /* Add initial tetrahedron faces to the face stack */
1137 tcentroid.x = tcentroid.y = tcentroid.z = 0.0;
1138 for (i = 0; i < 4; ++i) {
1139 context->valid[i] = 1;
1140 qh__assert_face(context->faces + i, context);
1141 qh__push_stack(&context->facestack, i);
1142 qh__vec3_add(&tcentroid, &context->faces[i].centroid);
1145 /* Assign the tetrahedron centroid */
1146 qh__vec3_multiply(&tcentroid, 0.25);
1147 context->centroid = tcentroid;
1149 QH_ASSERT(context->nedges == context->nfaces * 3);
1150 QH_ASSERT(context->nfaces == 4);
1151 QH_ASSERT(context->facestack.size == 4);
1157 qh__remove_vertex_duplicates(qh_context_t* context, double epsilon)
1160 for (i = 0; i < context->nvertices; ++i) {
1161 qh_vertex_t* v = context->vertices + i;
1162 if (v->x == 0) v->x = 0;
1163 if (v->y == 0) v->y = 0;
1164 if (v->z == 0) v->z = 0;
1165 for (j = i + 1; j < context->nvertices; ++j) {
1166 if (qh__vertex_equals_epsilon(context->vertices + i,
1167 context->vertices + j, epsilon))
1169 for (k = j; k < context->nvertices - 1; ++k) {
1170 context->vertices[k] = context->vertices[k+1];
1172 context->nvertices--;
1179 qh__init_context(qh_context_t* context, qh_vertex_t const* vertices,
1180 unsigned int nvertices)
1183 /* size_t nedges = 3 * nvertices - 6; */
1184 /* size_t nfaces = 2 * nvertices - 4; */
1185 unsigned int nfaces = nvertices * (nvertices - 1);
1186 unsigned int nedges = nfaces * 3;
1188 context->edges = QH_MALLOC(qh_half_edge_t, nedges);
1189 context->faces = QH_MALLOC(qh_face_t, nfaces);
1190 context->facestack.begin = QH_MALLOC(qh_index_t, nfaces);
1191 context->scratch.begin = QH_MALLOC(qh_index_t, nfaces);
1192 context->horizonedges.begin = QH_MALLOC(qh_index_t, nedges);
1193 context->newhorizonedges.begin = QH_MALLOC(qh_index_t, nedges);
1194 context->valid = QH_MALLOC(char, nfaces);
1196 if (!(context->edges &&
1198 context->facestack.begin &&
1199 context->scratch.begin &&
1200 context->horizonedges.begin &&
1201 context->newhorizonedges.begin &&
1204 jwxyz_abort ("%s: out of memory", progname);
1206 fprintf (stderr, "%s: out of memory\n", progname);
1212 context->vertices = QH_MALLOC(qh_vertex_t, nvertices);
1213 memcpy(context->vertices, vertices, sizeof(qh_vertex_t) * nvertices);
1215 context->nvertices = nvertices;
1216 context->nedges = 0;
1217 context->nfaces = 0;
1218 context->facestack.size = 0;
1219 context->scratch.size = 0;
1220 context->horizonedges.size = 0;
1221 context->newhorizonedges.size = 0;
1223 #ifdef QUICKHULL_DEBUG
1224 context->maxfaces = nfaces;
1225 context->maxedges = nedges;
1230 qh__free_context(qh_context_t* context)
1234 for (i = 0; i < context->nfaces; ++i) {
1235 QH_FREE(context->faces[i].iset.indices);
1236 context->faces[i].iset.size = 0;
1239 context->nvertices = 0;
1240 context->nfaces = 0;
1242 QH_FREE(context->edges);
1244 QH_FREE(context->faces);
1245 QH_FREE(context->facestack.begin);
1246 QH_FREE(context->scratch.begin);
1247 QH_FREE(context->horizonedges.begin);
1248 QH_FREE(context->newhorizonedges.begin);
1249 QH_FREE(context->vertices);
1250 QH_FREE(context->valid);
1254 qh_free_mesh(qh_mesh_t mesh)
1256 QH_FREE(mesh.vertices);
1257 QH_FREE(mesh.indices);
1258 QH_FREE(mesh.normalindices);
1259 QH_FREE(mesh.normals);
1263 qh__compute_epsilon(qh_vertex_t const* vertices, unsigned int nvertices)
1268 double maxxi = -QH_FLT_MAX;
1269 double maxyi = -QH_FLT_MAX;
1271 for (i = 0; i < nvertices; ++i) {
1272 double fxi = fabs(vertices[i].x);
1273 double fyi = fabs(vertices[i].y);
1283 epsilon = 2 * (maxxi + maxyi) * QH_FLT_EPS;
1289 qh_quickhull3d(qh_vertex_t const* vertices, unsigned int nvertices)
1292 qh_context_t context;
1293 /* unsigned int* indices; */
1294 unsigned int nfaces = 0, i, index, nindices;
1297 epsilon = qh__compute_epsilon(vertices, nvertices);
1299 qh__init_context(&context, vertices, nvertices);
1301 qh__remove_vertex_duplicates(&context, epsilon);
1303 /* Build the initial tetrahedron */
1304 qh__build_tetrahedron(&context, epsilon);
1306 /* Build the convex hull */
1307 #ifdef QUICKHULL_DEBUG
1308 qh__build_hull(&context, epsilon, -1, NULL);
1310 qh__build_hull(&context, epsilon);
1313 /* QH_ASSERT(qh__test_hull(&context, epsilon)); */
1315 for (i = 0; i < context.nfaces; ++i) {
1316 if (context.valid[i]) { nfaces++; }
1319 nindices = nfaces * 3;
1321 m.normals = QH_MALLOC(qh_vertex_t, nfaces);
1322 m.normalindices = QH_MALLOC(unsigned int, nfaces);
1323 m.vertices = QH_MALLOC(qh_vertex_t, nindices);
1324 m.indices = QH_MALLOC(unsigned int, nindices);
1325 m.nindices = nindices;
1326 m.nnormals = nfaces;
1331 for (i = 0; i < context.nfaces; ++i) {
1332 if (!context.valid[i]) { continue; }
1333 m.normals[index] = context.faces[i].normal;
1338 for (i = 0; i < context.nfaces; ++i) {
1339 if (!context.valid[i]) { continue; }
1340 m.normalindices[index] = index;
1345 for (i = 0; i < context.nfaces; ++i) {
1346 if (!context.valid[i]) { continue; }
1347 m.indices[index+0] = index+0;
1348 m.indices[index+1] = index+1;
1349 m.indices[index+2] = index+2;
1353 for (i = 0; i < context.nfaces; ++i) {
1354 if (!context.valid[i]) { continue; }
1355 qh_half_edge_t e0 = context.edges[context.faces[i].edges[0]];
1356 qh_half_edge_t e1 = context.edges[context.faces[i].edges[1]];
1357 qh_half_edge_t e2 = context.edges[context.faces[i].edges[2]];
1359 m.vertices[m.nvertices++] = context.vertices[e0.to_vertex];
1360 m.vertices[m.nvertices++] = context.vertices[e1.to_vertex];
1361 m.vertices[m.nvertices++] = context.vertices[e2.to_vertex];
1365 qh__free_context(&context);