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 #include "quickhull.h"
57 #include <math.h> /* sqrt & fabs */
58 #include <stdio.h> /* FILE */
59 #include <string.h> /* memcpy */
61 /* Quickhull helpers, define your own if needed */
62 #ifndef QUICKHULL_HELPERS
63 #include <stdlib.h> /* malloc, free, realloc */
64 #define QUICKHULL_HELPERS 1
65 #define QH_MALLOC(T, N) ((T*) malloc(N * sizeof(T)))
66 #define QH_REALLOC(T, P, N) ((T*)realloc(P, sizeof(T) * N))
67 #define QH_FREE(T) free(T)
68 #define QH_SWAP(T, A, B) { T tmp = B; B = A; A = tmp; }
69 #ifdef QUICKHULL_DEBUG
70 #define QH_ASSERT(STMT) if (!(STMT)) { *(int *)0 = 0; }
71 #define QH_LOG(FMT, ...) printf(FMT, ## __VA_ARGS__)
73 #define QH_ASSERT(STMT)
74 #define QH_LOG(FMT, ...)
75 #endif /* QUICKHULL_DEBUG */
76 #endif /* QUICKHULL_HELPERS */
79 #define QH_FLT_MAX 1e+37F
83 #define QH_FLT_EPS 1E-5F
86 #ifndef QH_VERTEX_SET_SIZE
87 #define QH_VERTEX_SET_SIZE 128
90 typedef long qh_index_t;
92 typedef struct qh_half_edge {
93 qh_index_t opposite_he; /* index of the opposite half edge */
94 qh_index_t next_he; /* index of the next half edge */
95 qh_index_t previous_he; /* index of the previous half edge */
96 qh_index_t he; /* index of the current half edge */
97 qh_index_t to_vertex; /* index of the next vertex */
98 qh_index_t adjacent_face; /* index of the ajacent face */
101 typedef struct qh_index_set {
104 unsigned int capacity;
107 typedef struct qh_face {
110 qh_vertex_t centroid;
117 typedef struct qh_index_stack {
122 typedef struct qh_context {
124 qh_half_edge_t* edges;
125 qh_vertex_t* vertices;
126 qh_vertex_t centroid;
127 qh_index_stack_t facestack;
128 qh_index_stack_t scratch;
129 qh_index_stack_t horizonedges;
130 qh_index_stack_t newhorizonedges;
133 unsigned int nvertices;
136 #ifdef QUICKHULL_DEBUG
137 unsigned int maxfaces;
138 unsigned int maxedges;
143 qh__find_6eps(qh_vertex_t* vertices, unsigned int nvertices, qh_index_t* eps)
145 qh_vertex_t* ptr = vertices;
147 double minxy = +QH_FLT_MAX;
148 double minxz = +QH_FLT_MAX;
149 double minyz = +QH_FLT_MAX;
151 double maxxy = -QH_FLT_MAX;
152 double maxxz = -QH_FLT_MAX;
153 double maxyz = -QH_FLT_MAX;
156 for (i = 0; i < 6; ++i) {
160 for (i = 0; i < nvertices; ++i) {
161 if (ptr->z < minxy) {
165 if (ptr->y < minxz) {
169 if (ptr->x < minyz) {
173 if (ptr->z > maxxy) {
177 if (ptr->y > maxxz) {
181 if (ptr->x > maxyz) {
190 qh__vertex_segment_length2(qh_vertex_t* p, qh_vertex_t* a, qh_vertex_t* b)
192 double dx = b->x - a->x;
193 double dy = b->y - a->y;
194 double dz = b->z - a->z;
196 double d = dx * dx + dy * dy + dz * dz;
203 double t = ((p->x - a->x) * dx +
205 (p->z - a->z) * dz) / d;
222 return dx * dx + dy * dy + dz * dz;
226 qh__vec3_sub(qh_vec3_t* a, qh_vec3_t* b)
234 qh__vec3_add(qh_vec3_t* a, qh_vec3_t* b)
242 qh__vec3_multiply(qh_vec3_t* a, double v)
250 qh__vertex_equals_epsilon(qh_vertex_t* a, qh_vertex_t* b, double epsilon)
252 return fabs(a->x - b->x) <= epsilon &&
253 fabs(a->y - b->y) <= epsilon &&
254 fabs(a->z - b->z) <= epsilon;
258 qh__vec3_length2(qh_vec3_t* v)
260 return v->x * v->x + v->y * v->y + v->z * v->z;
264 qh__vec3_dot(qh_vec3_t* v1, qh_vec3_t* v2)
266 return v1->x * v2->x + v1->y * v2->y + v1->z * v2->z;
270 qh__vec3_normalize(qh_vec3_t* v)
272 qh__vec3_multiply(v, 1.f / sqrt(qh__vec3_length2(v)));
276 qh__find_2dps_6eps(qh_vertex_t* vertices, qh_index_t* eps,
280 double max = -QH_FLT_MAX;
282 for (i = 0; i < 6; ++i) {
283 for (j = 0; j < 6; ++j) {
291 d = vertices[eps[i]];
292 qh__vec3_sub(&d, &vertices[eps[j]]);
293 d2 = qh__vec3_length2(&d);
305 qh__vec3_cross(qh_vec3_t* v1, qh_vec3_t* v2)
309 cross.x = v1->y * v2->z - v1->z * v2->y;
310 cross.y = v1->z * v2->x - v1->x * v2->z;
311 cross.z = v1->x * v2->y - v1->y * v2->x;
317 qh__face_centroid(qh_index_t vertices[3], qh_context_t* context)
319 qh_vertex_t centroid;
322 centroid.x = centroid.y = centroid.z = 0.0;
323 for (i = 0; i < 3; ++i) {
324 qh__vec3_add(¢roid, context->vertices + vertices[i]);
327 qh__vec3_multiply(¢roid, 1.0 / 3.0);
333 qh__dist_point_plane(qh_vertex_t* v, qh_vec3_t* normal, double sdist)
335 return fabs(qh__vec3_dot(v, normal) - sdist);
339 qh__init_half_edge(qh_half_edge_t* half_edge) {
340 half_edge->adjacent_face = -1;
342 half_edge->next_he = -1;
343 half_edge->opposite_he = -1;
344 half_edge->to_vertex = -1;
345 half_edge->previous_he = -1;
348 static qh_half_edge_t*
349 qh__next_edge(qh_context_t* context)
351 qh_half_edge_t* edge = context->edges + context->nedges;
352 qh__init_half_edge(edge);
354 edge->he = context->nedges;
357 QH_ASSERT(context->nedges < context->maxedges);
363 qh__next_face(qh_context_t* context)
365 qh_face_t* face = context->faces + context->nfaces;
367 face->face = context->nfaces;
368 face->iset.indices = NULL;
369 context->valid[context->nfaces] = 1;
372 QH_ASSERT(context->nfaces < context->maxfaces);
378 qh__edge_vec3(qh_half_edge_t* edge, qh_context_t* context)
380 qh_half_edge_t prevhe = context->edges[edge->previous_he];
383 v0 = context->vertices[prevhe.to_vertex];
384 v1 = context->vertices[edge->to_vertex];
386 qh__vec3_sub(&v1, &v0);
387 qh__vec3_normalize(&v1);
393 qh__face_init(qh_face_t* face, qh_index_t vertices[3],
394 qh_context_t* context)
396 qh_half_edge_t* e0 = qh__next_edge(context);
397 qh_half_edge_t* e1 = qh__next_edge(context);
398 qh_half_edge_t* e2 = qh__next_edge(context);
400 qh_vertex_t centroid, normal;
402 e2->to_vertex = vertices[0];
403 e0->to_vertex = vertices[1];
404 e1->to_vertex = vertices[2];
406 e0->next_he = e1->he;
407 e2->previous_he = e1->he;
408 face->edges[1] = e1->he;
410 e1->next_he = e2->he;
411 e0->previous_he = e2->he;
412 face->edges[2] = e2->he;
413 v1 = qh__edge_vec3(e2, context);
415 e2->next_he = e0->he;
416 e1->previous_he = e0->he;
417 face->edges[0] = e0->he;
418 v0 = qh__edge_vec3(e0, context);
420 e2->adjacent_face = face->face;
421 e1->adjacent_face = face->face;
422 e0->adjacent_face = face->face;
424 qh__vec3_multiply(&v1, -1.f);
425 normal = qh__vec3_cross(&v0, &v1);
427 qh__vec3_normalize(&normal);
428 centroid = qh__face_centroid(vertices, context);
429 face->centroid = centroid;
430 face->sdist = qh__vec3_dot(&normal, ¢roid);
431 face->normal = normal;
432 face->iset.indices = QH_MALLOC(qh_index_t, QH_VERTEX_SET_SIZE);
433 face->iset.capacity = QH_VERTEX_SET_SIZE;
435 face->visitededges = 0;
439 qh__tetrahedron_basis(qh_context_t* context, qh_index_t vertices[3])
442 int i, j = 0, k = 0, l = 0;
443 double max = -QH_FLT_MAX;
445 qh__find_6eps(context->vertices, context->nvertices, eps);
446 qh__find_2dps_6eps(context->vertices, eps, &j, &k);
448 for (i = 0; i < 6; ++i) {
451 if (i == j || i == k) {
455 d2 = qh__vertex_segment_length2(context->vertices + eps[i],
456 context->vertices + eps[j],
457 context->vertices + eps[k]);
465 vertices[0] = eps[j];
466 vertices[1] = eps[k];
467 vertices[2] = eps[l];
471 qh__push_stack(qh_index_stack_t* stack, qh_index_t index)
473 stack->begin[stack->size] = index;
478 qh__pop_stack(qh_index_stack_t* stack)
482 if (stack->size > 0) {
483 top = stack->begin[stack->size - 1];
491 qh__furthest_point_from_plane(qh_context_t* context,
498 double max = -QH_FLT_MAX;
500 for (i = 0; i < nindices; ++i) {
501 qh_index_t index = indices ? *(indices + i) : i;
502 double dist = qh__dist_point_plane(context->vertices + index,
515 qh__face_can_see_vertex(qh_face_t* face, qh_vertex_t* v)
519 qh__vec3_sub(&tov, &face->centroid);
520 return qh__vec3_dot(&tov, &face->normal) > 0;
524 qh__face_can_see_vertex_epsilon(qh_context_t* context, qh_face_t* face,
525 qh_vertex_t* v, double epsilon)
530 qh__vec3_sub(&tov, &face->centroid);
531 dot = qh__vec3_dot(&tov, &face->normal);
538 if (dot <= epsilon && dot >= 0) {
539 qh_vec3_t n = face->normal;
541 /* allow epsilon degeneration along the face normal */
542 qh__vec3_multiply(&n, epsilon);
553 qh__assert_half_edge(qh_half_edge_t* edge, qh_context_t* context)
555 QH_ASSERT(edge->opposite_he != -1);
556 QH_ASSERT(edge->he != -1);
557 QH_ASSERT(edge->adjacent_face != -1);
558 QH_ASSERT(edge->next_he != -1);
559 QH_ASSERT(edge->previous_he != -1);
560 QH_ASSERT(edge->to_vertex != -1);
561 QH_ASSERT(context->edges[edge->opposite_he].to_vertex != edge->to_vertex);
565 qh__assert_face(qh_face_t* face, qh_context_t* context)
569 for (i = 0; i < 3; ++i) {
570 qh__assert_half_edge(context->edges + face->edges[i], context);
573 QH_ASSERT(context->valid[face->face]);
576 #ifdef QUICKHULL_DEBUG
579 qh__log_face(qh_context_t* context, qh_face_t const* face) {
580 QH_LOG("Face %ld:\n", face->face);
581 for (int i = 0; i < 3; ++i) {
582 qh_half_edge_t edge = context->edges[face->edges[i]];
583 QH_LOG("\te%d %ld\n", i, edge.he);
584 QH_LOG("\t\te%d.opposite_he %ld\n", i, edge.opposite_he);
585 QH_LOG("\t\te%d.next_he %ld\n", i, edge.next_he);
586 QH_LOG("\t\te%d.previous_he %ld\n", i, edge.previous_he);
587 QH_LOG("\t\te%d.to_vertex %ld\n", i, edge.to_vertex);
588 QH_LOG("\t\te%d.adjacent_face %ld\n", i, edge.adjacent_face);
590 QH_LOG("\tnormal %f %f %f\n", face->normal.x, face->normal.y,
592 QH_LOG("\tsdist %f\n", face->sdist);
593 QH_LOG("\tcentroid %f %f %f\n", face->centroid.x, face->centroid.y,
598 qh__test_hull(qh_context_t* context, double epsilon, int testiset)
600 unsigned int i, j, k;
602 for (i = 0; i < context->nvertices; ++i) {
603 qh_index_t vindex = i;
606 for (j = 0; j < context->nfaces; ++j) {
607 if (!context->valid[j]) {
610 qh_face_t* face = context->faces + j;
612 qh_half_edge_t* e0 = context->edges + face->edges[0];
613 qh_half_edge_t* e1 = context->edges + face->edges[1];
614 qh_half_edge_t* e2 = context->edges + face->edges[2];
616 if (e0->to_vertex == vindex ||
617 e1->to_vertex == vindex ||
618 e2->to_vertex == vindex) {
624 for (k = 0; k < face->iset.size; ++k) {
625 if (vindex == face->iset.indices[k]) {
636 for (j = 0; j < context->nfaces; ++j) {
637 if (!context->valid[j]) {
640 qh_face_t* face = context->faces + j;
642 qh_vertex_t vertex = context->vertices[vindex];
643 qh__vec3_sub(&vertex, &face->centroid);
644 if (qh__vec3_dot(&face->normal, &vertex) > epsilon) {
645 #ifdef QUICKHULL_DEBUG
646 qh__log_face(context, face);
655 #endif /* QUICKHULL_DEBUG */
659 #ifdef QUICKHULL_DEBUG
660 qh__build_hull(qh_context_t* context, double epsilon, unsigned int step,
661 unsigned int* failurestep)
663 qh__build_hull(qh_context_t* context, double epsilon)
666 qh_index_t topface = qh__pop_stack(&context->facestack);
669 #ifdef QUICKHULL_DEBUG
670 unsigned int iteration = 0;
673 while (topface != -1) {
674 qh_face_t* face = context->faces + topface;
675 qh_index_t fvi, apex;
679 #ifdef QUICKHULL_DEBUG
680 if (!context->valid[topface] || face->iset.size == 0 || iteration == step)
682 if (!context->valid[topface] || face->iset.size == 0)
685 topface = qh__pop_stack(&context->facestack);
689 #ifdef QUICKHULL_DEBUG
690 if (failurestep != NULL && !qh__test_hull(context, epsilon, 1)) {
691 if (*failurestep == 0) {
692 *failurestep = iteration;
700 fvi = qh__furthest_point_from_plane(context, face->iset.indices,
701 face->iset.size, &face->normal, face->sdist);
702 fv = context->vertices + *(face->iset.indices + fvi);
704 qh__assert_face(face, context);
706 /* Reset visited flag for faces */
708 for (i = 0; i < context->nfaces; ++i) {
709 context->faces[i].visitededges = 0;
713 /* Find horizon edge */
715 qh_index_t tovisit = topface;
716 qh_face_t* facetovisit = context->faces + tovisit;
718 /* Release scratch */
719 context->scratch.size = 0;
721 while (tovisit != -1) {
722 if (facetovisit->visitededges >= 3) {
723 context->valid[tovisit] = 0;
724 tovisit = qh__pop_stack(&context->scratch);
725 facetovisit = context->faces + tovisit;
727 qh_index_t edgeindex = facetovisit->edges[facetovisit->visitededges];
728 qh_half_edge_t* edge;
729 qh_half_edge_t* oppedge;
732 facetovisit->visitededges++;
734 edge = context->edges + edgeindex;
735 oppedge = context->edges + edge->opposite_he;
736 adjface = context->faces + oppedge->adjacent_face;
738 if (!context->valid[oppedge->adjacent_face]) { continue; }
740 qh__assert_half_edge(oppedge, context);
741 qh__assert_half_edge(edge, context);
742 qh__assert_face(adjface, context);
744 if (!qh__face_can_see_vertex(adjface, fv)) {
745 qh__push_stack(&context->horizonedges, edge->he);
747 context->valid[tovisit] = 0;
748 qh__push_stack(&context->scratch, adjface->face);
754 apex = face->iset.indices[fvi];
756 /* Sort horizon edges in CCW order */
758 qh_vertex_t triangle[3];
760 qh_vec3_t v0, v1, toapex;
763 for (i = 0; i < context->horizonedges.size; ++i) {
764 qh_index_t he0 = context->horizonedges.begin[i];
765 qh_index_t he0vert = context->edges[he0].to_vertex;
766 qh_index_t phe0 = context->edges[he0].previous_he;
767 qh_index_t phe0vert = context->edges[phe0].to_vertex;
769 for (j = i + 2; j < context->horizonedges.size; ++j) {
770 qh_index_t he1 = context->horizonedges.begin[j];
771 qh_index_t he1vert = context->edges[he1].to_vertex;
772 qh_index_t phe1 = context->edges[he1].previous_he;
773 qh_index_t phe1vert = context->edges[phe1].to_vertex;
775 if (phe1vert == he0vert || phe0vert == he1vert) {
776 QH_SWAP(qh_index_t, context->horizonedges.begin[j],
777 context->horizonedges.begin[i + 1]);
784 context->vertices[context->edges[he0].to_vertex];
789 /* Detect first triangle face ordering */
793 qh__vec3_sub(&v0, &triangle[1]);
794 qh__vec3_sub(&v1, &triangle[1]);
796 n = qh__vec3_cross(&v0, &v1);
798 /* Get the vector to the apex */
799 toapex = triangle[0];
801 qh__vec3_sub(&toapex, context->vertices + apex);
803 reversed = qh__vec3_dot(&n, &toapex) < 0.f;
807 /* Create new faces */
809 qh_index_t top = qh__pop_stack(&context->horizonedges);
810 qh_index_t last = qh__pop_stack(&context->horizonedges);
811 qh_index_t first = top;
814 QH_ASSERT(context->newhorizonedges.size == 0);
816 /* Release scratch */
817 context->scratch.size = 0;
820 qh_half_edge_t* prevhe;
821 qh_half_edge_t* nexthe;
822 qh_half_edge_t* oppedge;
823 /* qh_vec3_t normal; */
824 /* qh_vertex_t fcentroid; */
833 prevhe = context->edges + last;
834 nexthe = context->edges + top;
837 QH_SWAP(qh_half_edge_t*, prevhe, nexthe);
840 verts[0] = prevhe->to_vertex;
841 verts[1] = nexthe->to_vertex;
844 context->valid[nexthe->adjacent_face] = 0;
846 oppedge = context->edges + nexthe->opposite_he;
847 newface = qh__next_face(context);
849 qh__face_init(newface, verts, context);
851 oppedge->opposite_he = context->edges[newface->edges[0]].he;
852 context->edges[newface->edges[0]].opposite_he = oppedge->he;
854 qh__push_stack(&context->scratch, newface->face);
855 qh__push_stack(&context->newhorizonedges, newface->edges[0]);
858 last = qh__pop_stack(&context->horizonedges);
862 /* Attach point sets to newly created faces */
864 for (k = 0; k < context->nfaces; ++k) {
865 qh_face_t* f = context->faces + k;
867 if (context->valid[k] || f->iset.size == 0) {
871 if (f->visitededges == 3) {
872 context->valid[k] = 0;
875 for (i = 0; i < f->iset.size; ++i) {
876 qh_index_t vertex = f->iset.indices[i];
877 /* qh_vertex_t* v = context->vertices + vertex; */
878 qh_face_t* dface = NULL;
880 for (j = 0; j < context->scratch.size; ++j) {
881 qh_face_t* newface = context->faces + context->scratch.begin[j];
882 qh_half_edge_t* e0 = context->edges + newface->edges[0];
883 qh_half_edge_t* e1 = context->edges + newface->edges[1];
884 qh_half_edge_t* e2 = context->edges + newface->edges[2];
885 /* qh_vertex_t cv; */
887 if (e0->to_vertex == vertex ||
888 e1->to_vertex == vertex ||
889 e2->to_vertex == vertex) {
893 if (qh__face_can_see_vertex_epsilon(context, newface,
894 context->vertices + vertex,
902 if (dface->iset.size + 1 >= dface->iset.capacity) {
903 dface->iset.capacity *= 2;
904 dface->iset.indices = QH_REALLOC(qh_index_t,
905 dface->iset.indices, dface->iset.capacity);
908 dface->iset.indices[dface->iset.size++] = vertex;
916 /* Link new faces together */
918 for (i = 0; i < context->newhorizonedges.size; ++i) {
919 qh_index_t phe0, nhe1;
925 ii = (i == 0) ? context->newhorizonedges.size - 1 : (i-1);
927 ii = (i+1) % context->newhorizonedges.size;
930 phe0 = context->edges[context->newhorizonedges.begin[i]].previous_he;
931 nhe1 = context->edges[context->newhorizonedges.begin[ii]].next_he;
933 he0 = context->edges + phe0;
934 he1 = context->edges + nhe1;
936 QH_ASSERT(he1->to_vertex == apex);
937 QH_ASSERT(he0->opposite_he == -1);
938 QH_ASSERT(he1->opposite_he == -1);
940 he0->opposite_he = he1->he;
941 he1->opposite_he = he0->he;
944 context->newhorizonedges.size = 0;
947 /* Push new face to stack */
949 for (i = 0; i < context->scratch.size; ++i) {
950 qh_face_t* face = context->faces + context->scratch.begin[i];
952 if (face->iset.size > 0) {
953 qh__push_stack(&context->facestack, face->face);
957 /* Release scratch */
958 context->scratch.size = 0;
961 topface = qh__pop_stack(&context->facestack);
963 /* TODO: push all non-valid faces for reuse in face stack memory pool */
967 #ifdef QUICKHULL_FILES
969 qh_mesh_export(qh_mesh_t const* mesh, char const* filename)
971 FILE* objfile = fopen(filename, "wt");
972 fprintf(objfile, "o\n");
974 for (int i = 0; i < mesh->nvertices; ++i) {
975 qh_vertex_t v = mesh->vertices[i];
976 fprintf(objfile, "v %f %f %f\n", v.x, v.y, v.z);
979 for (int i = 0; i < mesh->nnormals; ++i) {
980 qh_vec3_t n = mesh->normals[i];
981 fprintf(objfile, "vn %f %f %f\n", n.x, n.y, n.z);
984 for (int i = 0, j = 0; i < mesh->nindices; i += 3, j++) {
985 fprintf(objfile, "f %u/%u %u/%u %u/%u\n",
986 mesh->indices[i+0] + 1, mesh->normalindices[j] + 1,
987 mesh->indices[i+1] + 1, mesh->normalindices[j] + 1,
988 mesh->indices[i+2] + 1, mesh->normalindices[j] + 1);
993 #endif /* QUICKHULL_FILES */
996 qh__build_tetrahedron(qh_context_t* context, double epsilon)
999 qh_index_t vertices[3];
1002 qh_vertex_t normal, centroid, vapex, tcentroid;
1004 /* Get the initial tetrahedron basis (first face) */
1005 qh__tetrahedron_basis(context, &vertices[0]);
1007 /* Find apex from the tetrahedron basis */
1012 v0 = context->vertices[vertices[1]];
1013 v1 = context->vertices[vertices[2]];
1015 qh__vec3_sub(&v0, context->vertices + vertices[0]);
1016 qh__vec3_sub(&v1, context->vertices + vertices[0]);
1018 normal = qh__vec3_cross(&v0, &v1);
1019 qh__vec3_normalize(&normal);
1021 centroid = qh__face_centroid(vertices, context);
1022 sdist = qh__vec3_dot(&normal, ¢roid);
1024 apex = qh__furthest_point_from_plane(context, NULL,
1025 context->nvertices, &normal, sdist);
1026 vapex = context->vertices[apex];
1028 qh__vec3_sub(&vapex, ¢roid);
1030 /* Whether the face is looking towards the apex */
1031 if (qh__vec3_dot(&vapex, &normal) > 0) {
1032 QH_SWAP(qh_index_t, vertices[1], vertices[2]);
1036 faces = qh__next_face(context);
1037 qh__face_init(&faces[0], vertices, context);
1039 /* Build faces from the tetrahedron basis to the apex */
1041 qh_index_t facevertices[3];
1042 for (i = 0; i < 3; ++i) {
1043 qh_half_edge_t* edge = context->edges + faces[0].edges[i];
1044 qh_half_edge_t prevedge = context->edges[edge->previous_he];
1045 qh_face_t* face = faces+i+1;
1048 facevertices[0] = edge->to_vertex;
1049 facevertices[1] = prevedge.to_vertex;
1050 facevertices[2] = apex;
1052 qh__next_face(context);
1053 qh__face_init(face, facevertices, context);
1055 e0 = context->edges + faces[i+1].edges[0];
1056 edge->opposite_he = e0->he;
1057 e0->opposite_he = edge->he;
1061 /* Attach half edges to faces tied to the apex */
1063 for (i = 0; i < 3; ++i) {
1065 qh_face_t* nextface;
1072 nextface = faces+j+1;
1074 e1 = context->edges + face->edges[1];
1075 e2 = context->edges + nextface->edges[2];
1077 QH_ASSERT(e1->opposite_he == -1);
1078 QH_ASSERT(e2->opposite_he == -1);
1080 e1->opposite_he = e2->he;
1081 e2->opposite_he = e1->he;
1083 qh__assert_half_edge(e1, context);
1084 qh__assert_half_edge(e2, context);
1088 /* Create initial point set; every point is */
1089 /* attached to the first face it can see */
1091 for (i = 0; i < context->nvertices; ++i) {
1092 /* qh_vertex_t* v; */
1093 qh_face_t* dface = NULL;
1095 if (vertices[0] == i || vertices[1] == i || vertices[2] == i) {
1099 for (j = 0; j < 4; ++j) {
1100 if (qh__face_can_see_vertex_epsilon(context, context->faces + j,
1101 context->vertices + i, epsilon)) {
1102 dface = context->faces + j;
1111 for (j = 0; j < 3; ++j) {
1112 qh_half_edge_t* e = context->edges + dface->edges[j];
1113 if (i == e->to_vertex) {
1119 if (!valid) { continue; }
1121 if (dface->iset.size + 1 >= dface->iset.capacity) {
1122 dface->iset.capacity *= 2;
1123 dface->iset.indices = QH_REALLOC(qh_index_t,
1124 dface->iset.indices, dface->iset.capacity);
1127 dface->iset.indices[dface->iset.size++] = i;
1132 /* Add initial tetrahedron faces to the face stack */
1133 tcentroid.x = tcentroid.y = tcentroid.z = 0.0;
1134 for (i = 0; i < 4; ++i) {
1135 context->valid[i] = 1;
1136 qh__assert_face(context->faces + i, context);
1137 qh__push_stack(&context->facestack, i);
1138 qh__vec3_add(&tcentroid, &context->faces[i].centroid);
1141 /* Assign the tetrahedron centroid */
1142 qh__vec3_multiply(&tcentroid, 0.25);
1143 context->centroid = tcentroid;
1145 QH_ASSERT(context->nedges == context->nfaces * 3);
1146 QH_ASSERT(context->nfaces == 4);
1147 QH_ASSERT(context->facestack.size == 4);
1153 qh__remove_vertex_duplicates(qh_context_t* context, double epsilon)
1156 for (i = 0; i < context->nvertices; ++i) {
1157 qh_vertex_t* v = context->vertices + i;
1158 if (v->x == 0) v->x = 0;
1159 if (v->y == 0) v->y = 0;
1160 if (v->z == 0) v->z = 0;
1161 for (j = i + 1; j < context->nvertices; ++j) {
1162 if (qh__vertex_equals_epsilon(context->vertices + i,
1163 context->vertices + j, epsilon))
1165 for (k = j; k < context->nvertices - 1; ++k) {
1166 context->vertices[k] = context->vertices[k+1];
1168 context->nvertices--;
1175 qh__init_context(qh_context_t* context, qh_vertex_t const* vertices,
1176 unsigned int nvertices)
1179 /* size_t nedges = 3 * nvertices - 6; */
1180 /* size_t nfaces = 2 * nvertices - 4; */
1181 unsigned int nfaces = nvertices * (nvertices - 1);
1182 unsigned int nedges = nfaces * 3;
1184 context->edges = QH_MALLOC(qh_half_edge_t, nedges);
1185 context->faces = QH_MALLOC(qh_face_t, nfaces);
1186 context->facestack.begin = QH_MALLOC(qh_index_t, nfaces);
1187 context->scratch.begin = QH_MALLOC(qh_index_t, nfaces);
1188 context->horizonedges.begin = QH_MALLOC(qh_index_t, nedges);
1189 context->newhorizonedges.begin = QH_MALLOC(qh_index_t, nedges);
1190 context->valid = QH_MALLOC(char, nfaces);
1192 context->vertices = QH_MALLOC(qh_vertex_t, nvertices);
1193 memcpy(context->vertices, vertices, sizeof(qh_vertex_t) * nvertices);
1195 context->nvertices = nvertices;
1196 context->nedges = 0;
1197 context->nfaces = 0;
1198 context->facestack.size = 0;
1199 context->scratch.size = 0;
1200 context->horizonedges.size = 0;
1201 context->newhorizonedges.size = 0;
1203 #ifdef QUICKHULL_DEBUG
1204 context->maxfaces = nfaces;
1205 context->maxedges = nedges;
1210 qh__free_context(qh_context_t* context)
1214 for (i = 0; i < context->nfaces; ++i) {
1215 QH_FREE(context->faces[i].iset.indices);
1216 context->faces[i].iset.size = 0;
1219 context->nvertices = 0;
1220 context->nfaces = 0;
1222 QH_FREE(context->edges);
1224 QH_FREE(context->faces);
1225 QH_FREE(context->facestack.begin);
1226 QH_FREE(context->scratch.begin);
1227 QH_FREE(context->horizonedges.begin);
1228 QH_FREE(context->newhorizonedges.begin);
1229 QH_FREE(context->vertices);
1230 QH_FREE(context->valid);
1234 qh_free_mesh(qh_mesh_t mesh)
1236 QH_FREE(mesh.vertices);
1237 QH_FREE(mesh.indices);
1238 QH_FREE(mesh.normalindices);
1239 QH_FREE(mesh.normals);
1243 qh__compute_epsilon(qh_vertex_t const* vertices, unsigned int nvertices)
1248 double maxxi = -QH_FLT_MAX;
1249 double maxyi = -QH_FLT_MAX;
1251 for (i = 0; i < nvertices; ++i) {
1252 double fxi = fabs(vertices[i].x);
1253 double fyi = fabs(vertices[i].y);
1263 epsilon = 2 * (maxxi + maxyi) * QH_FLT_EPS;
1269 qh_quickhull3d(qh_vertex_t const* vertices, unsigned int nvertices)
1272 qh_context_t context;
1273 /* unsigned int* indices; */
1274 unsigned int nfaces = 0, i, index, nindices;
1277 epsilon = qh__compute_epsilon(vertices, nvertices);
1279 qh__init_context(&context, vertices, nvertices);
1281 qh__remove_vertex_duplicates(&context, epsilon);
1283 /* Build the initial tetrahedron */
1284 qh__build_tetrahedron(&context, epsilon);
1286 /* Build the convex hull */
1287 #ifdef QUICKHULL_DEBUG
1288 qh__build_hull(&context, epsilon, -1, NULL);
1290 qh__build_hull(&context, epsilon);
1293 /* QH_ASSERT(qh__test_hull(&context, epsilon)); */
1295 for (i = 0; i < context.nfaces; ++i) {
1296 if (context.valid[i]) { nfaces++; }
1299 nindices = nfaces * 3;
1301 m.normals = QH_MALLOC(qh_vertex_t, nfaces);
1302 m.normalindices = QH_MALLOC(unsigned int, nfaces);
1303 m.vertices = QH_MALLOC(qh_vertex_t, nindices);
1304 m.indices = QH_MALLOC(unsigned int, nindices);
1305 m.nindices = nindices;
1306 m.nnormals = nfaces;
1311 for (i = 0; i < context.nfaces; ++i) {
1312 if (!context.valid[i]) { continue; }
1313 m.normals[index] = context.faces[i].normal;
1318 for (i = 0; i < context.nfaces; ++i) {
1319 if (!context.valid[i]) { continue; }
1320 m.normalindices[index] = index;
1325 for (i = 0; i < context.nfaces; ++i) {
1326 if (!context.valid[i]) { continue; }
1327 m.indices[index+0] = index+0;
1328 m.indices[index+1] = index+1;
1329 m.indices[index+2] = index+2;
1333 for (i = 0; i < context.nfaces; ++i) {
1334 if (!context.valid[i]) { continue; }
1335 qh_half_edge_t e0 = context.edges[context.faces[i].edges[0]];
1336 qh_half_edge_t e1 = context.edges[context.faces[i].edges[1]];
1337 qh_half_edge_t e2 = context.edges[context.faces[i].edges[2]];
1339 m.vertices[m.nvertices++] = context.vertices[e0.to_vertex];
1340 m.vertices[m.nvertices++] = context.vertices[e1.to_vertex];
1341 m.vertices[m.nvertices++] = context.vertices[e2.to_vertex];
1345 qh__free_context(&context);