From 04be907cbba161af7eda231f8a6a1856fc3a6042 Mon Sep 17 00:00:00 2001 From: Franz Zotter Date: Sat, 25 Aug 2012 22:51:30 +0000 Subject: improved deterministic behavior of [mtx_qhull], however planar polygons do not fully work yet. svn path=/trunk/externals/iem/iemmatrix/; revision=16172 --- src/mtx_qhull.c | 13 +- src/mtx_qhull/list.c | 1 - src/mtx_qhull/list.h | 5 +- src/mtx_qhull/test.obj | 608 ++++++++++++++++++++++++++---- src/mtx_qhull/test_list.c | 2 +- src/mtx_qhull/test_write_conv_hull_obj.pd | 310 ++++++++------- src/mtx_qhull/vectors.c | 139 ++++++- src/mtx_qhull/vectors.h | 36 +- src/mtx_qhull/zhull.c | 323 +++++++++------- src/mtx_qhull/zhull.h | 17 +- 10 files changed, 1069 insertions(+), 385 deletions(-) (limited to 'src') diff --git a/src/mtx_qhull.c b/src/mtx_qhull.c index 7223b4b..34a4891 100644 --- a/src/mtx_qhull.c +++ b/src/mtx_qhull.c @@ -15,8 +15,6 @@ #include "iemmatrix.h" #include -#include "mtx_qhull/list.h" -#include "mtx_qhull/vectors.h" #include "mtx_qhull/zhull.h" static t_class *mtx_qhull_class; @@ -109,14 +107,9 @@ static void mTXQhullMatrix(MTXQhull *xo, t_symbol *s, SETFLOAT(xo->list,(float)xo->hull_size); SETFLOAT(xo->list+1,(float)3); for (i=0; ihull_size; i++) { - t_atom*ap=xo->list+2+3*i; - float f=(float)getEntry(getFacetByIndex(xo->zh->facets,i)->corners,0)+1; - SETFLOAT(xo->list+2+3*i, (float)getEntry( - getFacetByIndex(xo->zh->facets,i)->corners,0)+1); - SETFLOAT(xo->list+3+3*i, (float)getEntry( - getFacetByIndex(xo->zh->facets,i)->corners,1)+1); - SETFLOAT(xo->list+4+3*i, (float)getEntry( - getFacetByIndex(xo->zh->facets,i)->corners,2)+1); + SETFLOAT(xo->list+2+3*i, (float)getTriangleCorner(xo->zh,i,0)+1); + SETFLOAT(xo->list+3+3*i, (float)getTriangleCorner(xo->zh,i,1)+1); + SETFLOAT(xo->list+4+3*i, (float)getTriangleCorner(xo->zh,i,2)+1); } outlet_anything(xo->outl, gensym("matrix"), xo->size, xo->list); diff --git a/src/mtx_qhull/list.c b/src/mtx_qhull/list.c index a464ae9..1f8f58b 100644 --- a/src/mtx_qhull/list.c +++ b/src/mtx_qhull/list.c @@ -1,4 +1,3 @@ - #include "list.h" #include #include diff --git a/src/mtx_qhull/list.h b/src/mtx_qhull/list.h index dfccede..1a98256 100644 --- a/src/mtx_qhull/list.h +++ b/src/mtx_qhull/list.h @@ -1,16 +1,18 @@ #ifndef QHULL_LIST_H #define QHULL_LIST_H +#include +#include #include typedef long int entry_t; typedef long int index_t; - typedef struct list_ { entry_t *entries; size_t length; } list_t; + // memory things: list_t emptyList(void); void freeList(list_t *list); @@ -38,6 +40,7 @@ int inList(const entry_t entry, const list_t list); int notInList(const entry_t entry, const list_t list); list_t findValueListInList(const list_t value_list, const list_t list); entry_t findValueInList(const entry_t entry, const list_t list); +void uniquefyListEntries(list_t *list); void printList(const list_t list); #endif /* QHULL_LIST_H */ diff --git a/src/mtx_qhull/test.obj b/src/mtx_qhull/test.obj index 4b4e602..27d5b09 100644 --- a/src/mtx_qhull/test.obj +++ b/src/mtx_qhull/test.obj @@ -1,72 +1,536 @@ -v 1.5 0 0 -v 1.37304 0.603839 0.0119395 -v 1.00025 1.11771 0.0145916 -v 0.448097 1.4314 0.0170812 -v -0.339153 1.46107 0.0157868 -v -1.12324 0.994132 0.00513365 -v -1.49996 0.00517674 0.00970638 -v -1.1198 -0.99802 0.00258554 -v -0.284272 -1.47273 0.0159703 -v 0.518517 -1.40748 0.0113897 -v 1.06359 -1.05766 0.0118678 -v 1.39616 -0.548269 0.0118545 -v 1.21634 0.507813 0.715999 -v 0.496716 1.22144 0.715096 -v -0.543141 1.2091 0.702188 -v -1.21884 0.496413 0.719732 -v -1.20562 -0.527837 0.719625 -v -0.520718 -1.21172 0.714562 -v 0.548338 -1.19922 0.714976 -v 1.22223 -0.511408 0.703294 -v 0.558944 0.5949 1.25844 -v -0.560745 0.592992 1.25854 -v -0.566714 -0.600261 1.25241 -v 0.581676 -0.550909 1.26813 -f 24 23 22 -f 24 22 21 -f 19 23 24 -f 24 20 19 -f 23 19 18 -f 23 18 17 -f 16 22 23 -f 23 17 16 -f 22 16 15 -f 14 21 22 -f 22 15 14 -f 24 21 13 -f 13 20 24 -f 21 14 13 -f 11 19 20 -f 20 12 11 -f 19 11 10 -f 9 18 19 -f 19 10 9 -f 8 17 18 -f 18 9 8 -f 8 9 10 -f 7 16 17 -f 17 8 7 -f 6 15 16 -f 16 7 6 -f 6 7 8 -f 5 14 15 -f 15 6 5 -f 14 5 4 -f 4 5 6 -f 3 13 14 -f 3 5 6 -f 14 4 3 -f 3 4 6 -f 13 3 2 -f 20 13 1 -f 1 12 20 -f 1 11 12 -f 1 10 11 -f 1 8 10 -f 1 6 8 -f 1 5 6 -f 1 4 6 -f 1 3 6 -f 1 3 4 -f 13 2 1 -f 1 2 3 +v -1.29956 0.335546 0.669732 +v 0.335546 0.669732 -1.29956 +v 1.29956 -0.335546 0.669732 +v 0.669732 -1.29956 0.335546 +v -0.335546 0.669732 1.29956 +v 0.669732 1.29956 -0.335546 +v -0.335546 -0.669732 -1.29956 +v -0.669732 1.29956 0.335546 +v 0.335546 -0.669732 1.29956 +v -0.669732 -1.29956 -0.335546 +v -1.29956 -0.335546 -0.669732 +v 1.29956 0.335546 -0.669732 +v -1.21027 -0.692637 -0.552728 +v -0.692637 -0.552728 -1.21027 +v 1.21027 0.692637 -0.552728 +v -0.552728 -1.21027 -0.692637 +v 0.692637 -0.552728 1.21027 +v -0.552728 1.21027 0.692637 +v 0.692637 0.552728 -1.21027 +v 0.552728 1.21027 -0.692637 +v -0.692637 0.552728 1.21027 +v 0.552728 -1.21027 0.692637 +v -1.21027 0.692637 0.552728 +v 1.21027 -0.692637 0.552728 +v -0.202263 -0.060033 1.48509 +v -0.060033 1.48509 -0.202263 +v 0.202263 0.060033 1.48509 +v 1.48509 -0.202263 -0.060033 +v 0.060033 1.48509 0.202263 +v 1.48509 0.202263 0.060033 +v 0.060033 -1.48509 -0.202263 +v -1.48509 0.202263 -0.060033 +v -0.060033 -1.48509 0.202263 +v -1.48509 -0.202263 0.060033 +v -0.202263 0.060033 -1.48509 +v 0.202263 -0.060033 -1.48509 +v 0.074691 -0.419607 -1.43818 +v -0.419607 -1.43818 0.074691 +v -0.074691 0.419607 -1.43818 +v -1.43818 0.074691 -0.419607 +v 0.419607 -1.43818 -0.074691 +v -1.43818 -0.074691 0.419607 +v 0.419607 1.43818 0.074691 +v 1.43818 -0.074691 -0.419607 +v -0.419607 1.43818 -0.074691 +v 1.43818 0.074691 0.419607 +v 0.074691 0.419607 1.43818 +v -0.074691 -0.419607 1.43818 +v 0.308206 -0.289353 1.4392 +v -0.289353 1.4392 0.308206 +v -0.308206 0.289353 1.4392 +v 1.4392 0.308206 -0.289353 +v 0.289353 1.4392 -0.308206 +v 1.4392 -0.308206 0.289353 +v 0.289353 -1.4392 0.308206 +v -1.4392 -0.308206 -0.289353 +v -0.289353 -1.4392 -0.308206 +v -1.4392 0.308206 0.289353 +v 0.308206 0.289353 -1.4392 +v -0.308206 -0.289353 -1.4392 +v -0.418357 1.15801 -0.856735 +v 1.15801 -0.856735 -0.418357 +v 0.418357 -1.15801 -0.856735 +v -0.856735 -0.418357 1.15801 +v -1.15801 -0.856735 0.418357 +v -0.856735 0.418357 -1.15801 +v -1.15801 0.856735 -0.418357 +v 0.856735 0.418357 1.15801 +v 1.15801 0.856735 0.418357 +v 0.856735 -0.418357 -1.15801 +v -0.418357 -1.15801 0.856735 +v 0.418357 1.15801 0.856735 +v 1.36855 -0.589797 -0.170943 +v -0.589797 -0.170943 1.36855 +v -1.36855 0.589797 -0.170943 +v -0.170943 1.36855 -0.589797 +v 0.589797 -0.170943 -1.36855 +v -0.170943 -1.36855 0.589797 +v 0.589797 0.170943 1.36855 +v 0.170943 -1.36855 -0.589797 +v -0.589797 0.170943 -1.36855 +v 0.170943 1.36855 0.589797 +v 1.36855 0.589797 0.170943 +v -1.36855 -0.589797 0.170943 +v 1.27299 -0.019365 -0.793166 +v -0.019365 -0.793166 1.27299 +v -1.27299 0.019365 -0.793166 +v -0.793166 1.27299 -0.019365 +v 0.019365 -0.793166 -1.27299 +v -0.793166 -1.27299 0.019365 +v 0.019365 0.793166 1.27299 +v 0.793166 -1.27299 -0.019365 +v -0.019365 0.793166 -1.27299 +v 0.793166 1.27299 0.019365 +v 1.27299 0.019365 0.793166 +v -1.27299 -0.019365 0.793166 +v -1.04378 0.316747 1.02966 +v 0.316747 1.02966 -1.04378 +v 1.04378 -0.316747 1.02966 +v 1.02966 -1.04378 0.316747 +v -0.316747 1.02966 1.04378 +v 1.02966 1.04378 -0.316747 +v -0.316747 -1.02966 -1.04378 +v -1.02966 1.04378 0.316747 +v 0.316747 -1.02966 1.04378 +v -1.02966 -1.04378 -0.316747 +v -1.04378 -0.316747 -1.02966 +v 1.04378 0.316747 -1.02966 +v -0.392577 -0.872445 1.1553 +v -0.872445 1.1553 -0.392577 +v 0.392577 0.872445 1.1553 +v 1.1553 -0.392577 -0.872445 +v 0.872445 1.1553 0.392577 +v 1.1553 0.392577 0.872445 +v 0.872445 -1.1553 -0.392577 +v -1.1553 0.392577 -0.872445 +v -0.872445 -1.1553 0.392577 +v -1.1553 -0.392577 0.872445 +v -0.392577 0.872445 -1.1553 +v 0.392577 -0.872445 -1.1553 +v 0.759204 1.05149 0.753644 +v 1.05149 0.753644 0.759204 +v -0.759204 -1.05149 0.753644 +v 0.753644 0.759204 1.05149 +v -1.05149 0.753644 -0.759204 +v 0.753644 -0.759204 -1.05149 +v -1.05149 -0.753644 0.759204 +v -0.753644 -0.759204 1.05149 +v 1.05149 -0.753644 -0.759204 +v -0.753644 0.759204 -1.05149 +v 0.759204 -1.05149 -0.753644 +v -0.759204 1.05149 -0.753644 +v -0.661122 0.903363 0.998425 +v 0.903363 0.998425 -0.661122 +v 0.661122 -0.903363 0.998425 +v 0.998425 -0.661122 0.903363 +v -0.903363 0.998425 0.661122 +v 0.998425 0.661122 -0.903363 +v -0.903363 -0.998425 -0.661122 +v -0.998425 0.661122 0.903363 +v 0.903363 -0.998425 0.661122 +v -0.998425 -0.661122 -0.903363 +v -0.661122 -0.903363 -0.998425 +v 0.661122 0.903363 -0.998425 +v -1.33539 -0.440278 0.522396 +v -0.440278 0.522396 -1.33539 +v 1.33539 0.440278 0.522396 +v 0.522396 -1.33539 -0.440278 +v 0.440278 0.522396 1.33539 +v 0.522396 1.33539 0.440278 +v 0.440278 -0.522396 -1.33539 +v -0.522396 1.33539 -0.440278 +v -0.440278 -0.522396 1.33539 +v -0.522396 -1.33539 0.440278 +v -1.33539 0.440278 -0.522396 +v 1.33539 -0.440278 -0.522396 +v 0.992958 0.047085 1.12331 +v 0.047085 1.12331 0.992958 +v -0.992958 -0.047085 1.12331 +v 1.12331 0.992958 0.047085 +v -0.047085 1.12331 -0.992958 +v 1.12331 -0.992958 -0.047085 +v -0.047085 -1.12331 0.992958 +v -1.12331 -0.992958 0.047085 +v 0.047085 -1.12331 -0.992958 +v -1.12331 0.992958 -0.047085 +v 0.992958 -0.047085 -1.12331 +v -0.992958 0.047085 -1.12331 +v -0.188599 -1.31655 -0.69364 +v -1.31655 -0.69364 -0.188599 +v 0.188599 1.31655 -0.69364 +v -0.69364 -0.188599 -1.31655 +v 1.31655 -0.69364 0.188599 +v -0.69364 0.188599 1.31655 +v 1.31655 0.69364 -0.188599 +v 0.69364 0.188599 -1.31655 +v -1.31655 0.69364 0.188599 +v 0.69364 -0.188599 1.31655 +v -0.188599 1.31655 0.69364 +v 0.188599 -1.31655 0.69364 +f 71 78 163 +f 78 71 154 +f 57 38 10 +f 86 153 109 +f 163 86 109 +f 71 163 109 +f 57 10 16 +f 55 4 22 +f 4 141 22 +f 107 172 14 +f 55 78 33 +f 154 38 33 +f 78 154 33 +f 154 117 90 +f 38 154 90 +f 10 38 90 +f 154 71 123 +f 117 154 123 +f 78 55 180 +f 163 78 180 +f 55 22 180 +f 129 62 131 +f 126 129 131 +f 63 126 131 +f 143 14 7 +f 74 159 64 +f 153 74 64 +f 118 127 64 +f 159 118 64 +f 146 35 81 +f 53 20 171 +f 53 171 76 +f 53 76 26 +f 26 45 29 +f 5 47 91 +f 101 5 91 +f 158 101 91 +f 29 45 50 +f 57 16 169 +f 126 63 120 +f 151 126 120 +f 89 151 120 +f 89 7 37 +f 151 89 37 +f 162 62 73 +f 109 153 128 +f 71 109 128 +f 123 71 128 +f 127 123 128 +f 64 127 128 +f 153 64 128 +f 122 147 69 +f 150 94 43 +f 29 150 43 +f 26 29 43 +f 53 26 43 +f 86 9 48 +f 153 86 48 +f 93 146 119 +f 159 74 174 +f 101 18 133 +f 5 101 133 +f 21 5 133 +f 174 74 51 +f 21 174 51 +f 5 21 51 +f 47 5 51 +f 117 123 65 +f 123 127 65 +f 16 10 139 +f 143 16 139 +f 35 37 60 +f 81 35 60 +f 172 81 60 +f 14 172 60 +f 7 14 60 +f 37 7 60 +f 94 102 6 +f 43 94 6 +f 53 43 6 +f 20 53 6 +f 94 150 113 +f 27 79 149 +f 47 27 149 +f 91 47 149 +f 127 118 145 +f 65 127 145 +f 10 90 106 +f 139 10 106 +f 143 139 142 +f 14 143 142 +f 107 14 142 +f 122 69 121 +f 124 122 121 +f 72 124 121 +f 150 72 121 +f 113 150 121 +f 69 113 121 +f 172 107 168 +f 81 172 168 +f 119 146 130 +f 61 119 130 +f 61 130 132 +f 104 88 166 +f 177 104 166 +f 88 110 166 +f 110 132 67 +f 166 110 67 +f 86 163 105 +f 9 86 105 +f 180 22 105 +f 163 180 105 +f 118 159 96 +f 145 118 96 +f 62 129 156 +f 73 62 156 +f 17 136 99 +f 136 3 99 +f 9 17 49 +f 48 9 49 +f 27 48 49 +f 79 27 49 +f 76 171 161 +f 61 76 161 +f 119 61 161 +f 93 119 161 +f 72 150 82 +f 158 72 82 +f 29 50 82 +f 150 29 82 +f 57 169 80 +f 104 177 23 +f 90 117 164 +f 106 90 164 +f 170 106 164 +f 117 65 164 +f 147 30 83 +f 69 147 83 +f 65 145 84 +f 164 65 84 +f 170 164 84 +f 145 34 84 +f 146 81 66 +f 130 146 66 +f 168 116 66 +f 81 168 66 +f 107 142 11 +f 116 155 125 +f 66 116 125 +f 130 66 125 +f 132 130 125 +f 67 132 125 +f 155 67 125 +f 73 156 28 +f 79 49 178 +f 157 79 178 +f 99 157 178 +f 17 99 178 +f 49 17 178 +f 158 82 179 +f 101 158 179 +f 18 101 179 +f 50 18 179 +f 82 50 179 +f 63 131 148 +f 80 63 148 +f 83 30 52 +f 175 83 52 +f 104 23 137 +f 133 18 137 +f 84 34 56 +f 170 84 56 +f 11 142 13 +f 56 11 13 +f 170 56 13 +f 106 170 13 +f 139 106 13 +f 142 139 13 +f 156 85 44 +f 28 156 44 +f 30 28 44 +f 52 30 44 +f 62 162 115 +f 131 62 115 +f 148 131 115 +f 92 148 115 +f 162 92 115 +f 2 59 39 +f 93 2 39 +f 146 93 39 +f 35 146 39 +f 1 96 97 +f 174 21 97 +f 159 174 97 +f 96 159 97 +f 89 120 165 +f 80 169 165 +f 63 80 165 +f 120 63 165 +f 162 73 173 +f 122 124 68 +f 79 157 68 +f 149 79 68 +f 124 149 68 +f 17 9 135 +f 136 17 135 +f 141 136 135 +f 22 141 135 +f 105 22 135 +f 9 105 135 +f 57 80 31 +f 38 57 31 +f 33 38 31 +f 80 148 31 +f 168 107 87 +f 116 168 87 +f 155 116 87 +f 107 11 87 +f 23 1 140 +f 137 23 140 +f 133 137 140 +f 21 133 140 +f 97 21 140 +f 1 97 140 +f 73 28 54 +f 173 73 54 +f 155 32 75 +f 67 155 75 +f 166 67 75 +f 177 166 75 +f 96 1 42 +f 145 96 42 +f 34 145 42 +f 32 34 42 +f 2 144 19 +f 59 2 19 +f 144 138 19 +f 155 87 40 +f 32 155 40 +f 34 32 40 +f 56 34 40 +f 11 56 40 +f 87 11 40 +f 85 167 108 +f 19 138 108 +f 88 104 8 +f 45 88 8 +f 50 45 8 +f 18 50 8 +f 137 18 8 +f 104 137 8 +f 102 94 160 +f 175 102 160 +f 83 175 160 +f 69 83 160 +f 113 69 160 +f 94 113 160 +f 91 149 111 +f 158 91 111 +f 72 158 111 +f 124 72 111 +f 149 124 111 +f 2 93 98 +f 144 2 98 +f 20 144 98 +f 171 20 98 +f 161 171 98 +f 93 161 98 +f 175 52 15 +f 102 175 15 +f 147 122 114 +f 68 157 114 +f 122 68 114 +f 55 33 41 +f 4 55 41 +f 92 4 41 +f 148 92 41 +f 31 148 41 +f 33 31 41 +f 77 167 70 +f 151 77 70 +f 126 151 70 +f 129 126 70 +f 102 15 134 +f 6 102 134 +f 20 6 134 +f 144 20 134 +f 138 144 134 +f 15 138 134 +f 147 114 95 +f 99 3 95 +f 157 99 95 +f 114 157 95 +f 77 151 36 +f 59 77 36 +f 39 59 36 +f 35 39 36 +f 37 35 36 +f 151 37 36 +f 76 61 152 +f 26 76 152 +f 45 26 152 +f 88 45 152 +f 110 88 152 +f 132 110 152 +f 61 132 152 +f 162 173 100 +f 92 162 100 +f 4 92 100 +f 141 4 100 +f 77 59 176 +f 167 77 176 +f 108 167 176 +f 19 108 176 +f 59 19 176 +f 129 70 112 +f 156 129 112 +f 85 156 112 +f 167 85 112 +f 70 167 112 +f 141 100 24 +f 136 141 24 +f 3 136 24 +f 54 3 24 +f 173 54 24 +f 100 173 24 +f 7 89 103 +f 143 7 103 +f 16 143 103 +f 169 16 103 +f 165 169 103 +f 89 165 103 +f 108 138 12 +f 85 108 12 +f 44 85 12 +f 52 44 12 +f 15 52 12 +f 138 15 12 +f 74 153 25 +f 51 74 25 +f 47 51 25 +f 27 47 25 +f 48 27 25 +f 153 48 25 +f 147 95 46 +f 30 147 46 +f 28 30 46 +f 54 28 46 +f 3 54 46 +f 95 3 46 +f 1 23 58 +f 42 1 58 +f 32 42 58 +f 75 32 58 +f 177 75 58 +f 23 177 58 diff --git a/src/mtx_qhull/test_list.c b/src/mtx_qhull/test_list.c index 4ddddde..9ac5924 100644 --- a/src/mtx_qhull/test_list.c +++ b/src/mtx_qhull/test_list.c @@ -6,7 +6,7 @@ int main(char **argv, int argc) { list_t l1=emptyList(); list_t l2=emptyList(); list_t l3=emptyList(); - int x[]={0, 2, 4, 6}; + entry_t x[]={0, 2, 4, 6}; printf("empty list:\n"); printList(l1); diff --git a/src/mtx_qhull/test_write_conv_hull_obj.pd b/src/mtx_qhull/test_write_conv_hull_obj.pd index d27370d..f73dff1 100644 --- a/src/mtx_qhull/test_write_conv_hull_obj.pd +++ b/src/mtx_qhull/test_write_conv_hull_obj.pd @@ -1,6 +1,4 @@ -#N canvas 534 0 850 668 10; -#X obj 595 205 convhull 0.001; -#X obj 595 176 t a a; +#N canvas 368 486 824 447 10; #N canvas 69 11 539 331 writeobject 0; #X obj 163 270 textfile; #X obj 163 161 mtx; @@ -48,7 +46,7 @@ #X connect 17 0 0 0; #X connect 18 0 11 0; #X connect 19 0 12 0; -#X restore 116 490 pd writeobject; +#X restore 217 444 pd writeobject; #N canvas 0 22 450 300 gemwin 0; #X obj 132 136 gemwin; #X obj 67 89 outlet; @@ -71,131 +69,187 @@ #X connect 7 0 0 0; #X connect 8 0 0 0; #X connect 9 0 8 0; -#X restore 556 423 pd gemwin; -#X msg 556 404 destroy; -#X obj 668 404 gemhead; -#X obj 668 423 world_light; -#X obj 362 -1 gemhead; -#X obj 393 151 model; -#X floatatom 377 79 5 0 0 0 - - -; -#X floatatom 420 80 5 0 0 0 - - -; -#X floatatom 460 80 5 0 0 0 - - -; -#X msg 17 196 open test.obj; -#X msg 500 97 rescale 0; -#X obj 364 38 color 0.9 0.9 0.9; -#X obj 595 227 extra_point; -#X obj 209 467 mtx_* 1.5; -#X obj 363 18 ortho; -#X obj 488 -3 tgl 15 0 empty empty empty 17 7 0 10 -262144 -1 -1 0 +#X restore 376 128 pd gemwin; +#X msg 376 109 destroy; +#X obj 488 109 gemhead; +#X obj 494 129 world_light; +#X obj -11 319 gemhead; +#X obj 79 452 model; +#X floatatom 6 382 5 0 0 0 - - -; +#X floatatom 49 383 5 0 0 0 - - -; +#X floatatom 89 383 5 0 0 0 - - -; +#X msg 33 333 open test.obj; +#X msg 129 400 rescale 0; +#X obj -9 358 color 0.9 0.9 0.9; +#X obj 335 426 mtx_* 1.5; +#X obj -10 338 ortho; +#X obj 137 309 tgl 15 0 empty empty empty 17 7 0 10 -262144 -1 -1 0 1; -#X obj 485 38 f; -#X obj 540 38 % 360; -#X obj 480 17 metro 50; -#X msg 360 176 smooth 0; -#X obj 595 266 convhull 0.001; -#X obj 595 248 t a a; -#X obj 85 411 mtx; -#X obj 84 446 mtx_slice; -#X msg 123 425 1 1 \$1 end; -#X floatatom 18 250 5 0 0 0 - - -; -#X obj 221 440 mtx; -#X obj -7 284 t b b f b; -#X obj 125 394 mtx_size; -#X obj 101 375 t a a; -#X msg 184 396 set \$1; -#X obj 355 101 rotateXYZ -74 0 0; -#X obj 513 38 + 3; -#X obj 360 128 t b a; -#X obj 113 332 convhull 0.001; -#X obj 159 229 t a a; -#X msg 105 72 4 40; -#X msg 109 92 5 67; -#X msg 103 50 3 26; -#X msg 101 137 read designsN/N\$1_\$2.mtx \, bang; -#X obj 60 135 mtx; -#X obj 183 280 demux; -#X obj 210 213 tgl 15 0 empty empty empty 17 7 0 10 -262144 -1 -1 1 +#X obj 112 358 f; +#X obj 167 358 % 360; +#X obj 129 329 metro 50; +#X msg 9 452 smooth 0; +#X obj 213 386 mtx; +#X obj 212 421 mtx_slice; +#X msg 251 400 1 1 \$1 end; +#X floatatom 73 276 5 0 0 0 - - -; +#X obj 334 402 mtx; +#X obj 113 277 t b b f b; +#X obj 253 369 mtx_size; +#X obj 229 350 t a a; +#X msg 312 371 set \$1; +#X obj -16 404 rotateXYZ -74 0 0; +#X obj 140 358 + 3; +#X obj 46 429 t b a; +#X obj 335 317 convhull 0.001; +#X obj 231 271 t a a; +#X msg 151 83 4 40; +#X msg 188 66 5 67; +#X msg 149 67 3 26; +#X msg 43 124 read designsN/N\$1_\$2.mtx \, bang; +#X obj 41 205 mtx; +#X obj 236 294 demux; +#X obj 123 182 tgl 15 0 empty empty empty 17 7 0 10 -262144 -1 -1 0 1; -#X obj 249 237 t b f; -#X obj 60 165 t b a b; -#X msg 111 111 9 180; -#X msg 97 1 1.5 6; -#X msg 104 27 2 14; -#X msg 181 43 2.5 12; -#X floatatom 323 280 5 0 0 0 - - -; -#X obj 327 301 t b f; -#X obj 416 257 mtx_rand; -#X msg 427 235 200 3; -#X obj 245 336 mtx_qhull 500; -#X msg 175 77 read designsN/dode.mtx \, bang; -#X msg 221 167 read Rls.mtx \, bang; -#X connect 0 0 15 0; -#X connect 1 0 0 0; -#X connect 1 1 15 1; +#X obj 71 180 t b f; +#X msg 190 85 9 180; +#X msg 48 70 1.5 6; +#X msg 97 66 2 14; +#X msg 46 88 2.5 12; +#X floatatom 396 171 5 0 0 0 - - -; +#X obj 378 199 t b f; +#X obj -18 163 mtx_rand; +#X msg -7 141 200 3; +#X obj 237 319 mtx_qhull 500; +#X msg -123 202 read Rls.mtx \, bang; +#X msg 68 150 matrix 4 3 0 0 0 1 0 0 0 1 0 1 1 0; +#X msg 100 85 2.5 20; +#X obj 107 239 mtx_size; +#X obj 42 235 t b a a; +#X obj 175 238 * 2; +#X obj 207 238 - 4; +#X floatatom 236 237 5 0 0 0 - - -; +#X obj 396 64 mtx_print; +#X msg 278 37 matrix 8 3 0 0 0 0 0 1 0 1 0 0 1 1 1 0 0 1 0 1 1 1 0 +1 1 1; +#X obj 293 74 mtx_roll -1; +#X obj 388 256 mtx_print; +#N canvas 0 0 624 397 makepolygon 0; +#X obj 75 114 mtx_:; +#X obj 80 44 t f f; +#X obj 83 70 - 1; +#X msg 83 90 0 \$1; +#X obj 89 139 mtx_./ 1; +#X obj 87 187 mtx_* 3.14159; +#X obj 92 163 mtx_* 2; +#X obj 94 208 t a a a; +#X obj 150 210 mtx_size; +#X msg 218 208 1 \$1; +#X obj 258 208 mtx_zeros; +#X obj 103 238 mtx_cos; +#X obj 136 265 mtx_sin; +#X obj 143 341 mtx_transpose; +#X obj 147 363 outlet; +#X obj 114 17 inlet; +#X obj 279 340 mtx_print; +#X obj 135 317 mtx_concat col; +#X obj 173 291 mtx_concat col; +#X msg 68 19 5; +#X connect 0 0 4 0; +#X connect 1 0 2 0; +#X connect 1 1 4 1; +#X connect 2 0 3 0; +#X connect 3 0 0 0; +#X connect 4 0 6 0; +#X connect 5 0 7 0; +#X connect 6 0 5 0; +#X connect 7 0 11 0; +#X connect 7 1 12 0; +#X connect 7 2 8 0; +#X connect 8 1 9 0; +#X connect 9 0 10 0; +#X connect 10 0 18 1; +#X connect 11 0 17 0; +#X connect 12 0 18 0; +#X connect 13 0 14 0; +#X connect 13 0 16 0; +#X connect 15 0 1 0; +#X connect 17 0 13 0; +#X connect 18 0 17 1; +#X connect 19 0 1 0; +#X restore 251 176 pd makepolygon; +#X msg 323 150 5; +#X connect 1 0 2 0; +#X connect 2 0 1 0; #X connect 3 0 4 0; -#X connect 4 0 3 0; -#X connect 5 0 6 0; -#X connect 7 0 17 0; -#X connect 9 0 34 1; -#X connect 10 0 34 2; -#X connect 11 0 34 3; -#X connect 12 0 36 0; -#X connect 13 0 8 0; -#X connect 14 0 34 0; -#X connect 15 0 24 0; -#X connect 16 0 2 1; -#X connect 17 0 14 0; -#X connect 18 0 21 0; -#X connect 19 0 11 0; -#X connect 19 0 35 0; -#X connect 20 0 19 1; -#X connect 21 0 19 0; -#X connect 22 0 8 0; -#X connect 23 0 32 0; -#X connect 24 0 23 0; -#X connect 24 1 29 0; -#X connect 25 0 26 0; -#X connect 26 0 2 0; -#X connect 27 0 26 1; -#X connect 28 0 30 0; -#X connect 29 0 16 0; -#X connect 30 0 12 0; -#X connect 30 1 25 0; -#X connect 30 2 27 0; -#X connect 30 3 29 0; -#X connect 31 0 27 0; -#X connect 31 0 33 0; -#X connect 32 0 25 0; -#X connect 32 1 31 0; -#X connect 33 0 28 0; -#X connect 34 0 8 0; -#X connect 35 0 20 0; -#X connect 36 0 22 0; -#X connect 36 1 8 0; -#X connect 37 0 32 0; -#X connect 38 0 44 0; -#X connect 38 1 29 0; -#X connect 39 0 42 0; -#X connect 40 0 42 0; -#X connect 41 0 42 0; -#X connect 42 0 43 0; -#X connect 43 0 47 0; +#X connect 5 0 14 0; +#X connect 7 0 29 1; +#X connect 8 0 29 2; +#X connect 9 0 29 3; +#X connect 10 0 31 0; +#X connect 11 0 6 0; +#X connect 12 0 29 0; +#X connect 13 0 0 1; +#X connect 14 0 12 0; +#X connect 15 0 18 0; +#X connect 16 0 9 0; +#X connect 16 0 30 0; +#X connect 17 0 16 1; +#X connect 18 0 16 0; +#X connect 19 0 6 0; +#X connect 20 0 21 0; +#X connect 21 0 0 0; +#X connect 22 0 21 1; +#X connect 23 0 25 0; +#X connect 24 0 13 0; +#X connect 25 0 10 0; +#X connect 25 1 20 0; +#X connect 25 2 22 0; +#X connect 25 3 24 0; +#X connect 26 0 22 0; +#X connect 26 0 28 0; +#X connect 27 0 20 0; +#X connect 27 1 26 0; +#X connect 28 0 23 0; +#X connect 29 0 6 0; +#X connect 30 0 17 0; +#X connect 31 0 19 0; +#X connect 31 1 6 0; +#X connect 32 0 27 0; +#X connect 33 0 39 0; +#X connect 33 1 24 0; +#X connect 34 0 37 0; +#X connect 35 0 37 0; +#X connect 36 0 37 0; +#X connect 37 0 38 0; +#X connect 38 0 55 0; +#X connect 39 0 50 0; +#X connect 39 1 32 0; +#X connect 40 0 41 0; +#X connect 41 0 38 0; +#X connect 41 1 39 1; +#X connect 42 0 37 0; +#X connect 43 0 37 0; #X connect 44 0 37 0; -#X connect 44 1 56 0; -#X connect 45 0 46 0; -#X connect 46 0 43 0; -#X connect 46 1 44 1; -#X connect 47 0 12 0; -#X connect 47 1 38 0; -#X connect 48 0 42 0; -#X connect 49 0 42 0; -#X connect 50 0 42 0; -#X connect 51 0 42 0; -#X connect 52 0 53 0; -#X connect 53 0 43 0; -#X connect 53 1 56 0; -#X connect 54 0 43 0; -#X connect 55 0 54 0; -#X connect 56 0 32 0; -#X connect 57 0 43 0; -#X connect 58 0 43 0; +#X connect 45 0 37 0; +#X connect 46 0 47 0; +#X connect 47 0 38 0; +#X connect 47 1 50 0; +#X connect 48 0 38 0; +#X connect 49 0 48 0; +#X connect 50 0 27 0; +#X connect 50 0 62 0; +#X connect 51 0 38 0; +#X connect 52 0 38 0; +#X connect 53 0 37 0; +#X connect 54 0 56 0; +#X connect 55 0 10 0; +#X connect 55 1 33 0; +#X connect 55 2 54 0; +#X connect 56 0 57 0; +#X connect 57 0 58 0; +#X connect 60 0 61 0; +#X connect 61 0 38 0; +#X connect 61 0 59 0; +#X connect 63 0 38 0; +#X connect 64 0 63 0; diff --git a/src/mtx_qhull/vectors.c b/src/mtx_qhull/vectors.c index 6cc3dd8..d9d3bd0 100644 --- a/src/mtx_qhull/vectors.c +++ b/src/mtx_qhull/vectors.c @@ -1,7 +1,4 @@ #include "vectors.h" -#include -#include -#include vector_t initVector (float x, float y, float z) { @@ -9,13 +6,16 @@ vector_t initVector (float x, float y, float z) { return vec; } +float lengthVector(vector_t v) { + return sqrtf(v.c[0]*v.c[0]+v.c[1]*v.c[1]+v.c[2]*v.c[2]); +} vector_t normalizeVector(vector_t v) { - float r=sqrtf(v.c[0]*v.c[0]+v.c[1]*v.c[1]+v.c[2]*v.c[2]); + float r=lengthVector(v); v.c[0]/=r; v.c[1]/=r; v.c[2]/=r; return v; -}; +} plane_t initPlane (vector_t normal, vector_t point) { plane_t plane; @@ -24,6 +24,13 @@ plane_t initPlane (vector_t normal, vector_t point) { return plane; } +line_t initLine (vector_t direction, vector_t point) { + line_t line; + line.point = point; + line.direction = normalizeVector(direction); + return line; +} + points_t initPoints (const float *x, const float *y, const float *z, size_t num_points) { @@ -63,11 +70,29 @@ float innerProduct (vector_t v1, vector_t v2) { return v1.c[0]*v2.c[0] + v1.c[1]*v2.c[1] + v1.c[2]*v2.c[2]; } +float distancePointPoint(const vector_t a, + const vector_t b) { + vector_t d=subtractVectors(b,a); + return lengthVector(d); +} + float distancePointPlane (vector_t point, plane_t plane) { return innerProduct(point, plane.normal) - innerProduct(plane.point, plane.normal); } +float distancePointLine (vector_t point, line_t line) { + return lengthVector(crossProduct(line.direction, + subtractVectors(point,line.point))); +} + +float distancePointLineOnPlane (vector_t point, + line_t line, plane_t plane) { + vector_t normal_in_plane = normalizeVector(crossProduct( + line.direction, plane.normal)); + return innerProduct(subtractVectors(point,line.point),normal_in_plane); +} + vector_t addVectors(vector_t v1, vector_t v2) { vector_t v3; v3.c[0]=v1.c[0]+v2.c[0]; @@ -92,31 +117,113 @@ vector_t scaleVector(vector_t v1, float f) { return v2; } +/* vector_t averagePoints(points_t points) { vector_t m = initVector(0.0f, 0.0f, 0.0f); size_t i; - for (i=0; i=0)&&(index +#include +#include typedef struct vec_ { float c[3]; } vector_t; @@ -18,31 +18,49 @@ typedef struct plane_ { vector_t point; } plane_t; +typedef struct line_ { + vector_t direction; + vector_t point; +} line_t; + + vector_t initVector (float x, float y, float z); +float lengthVector(vector_t v); vector_t normalizeVector(vector_t v); plane_t initPlane (vector_t normal, vector_t point); +line_t initLine (vector_t direction, vector_t point); points_t initPoints (const float *x, const float *y, const float *z, size_t num_points); void freePoints (points_t *points); +size_t getNumPoints(const points_t points); +vector_t getPoint(const points_t points, const index_t index); vector_t crossProduct (vector_t v1, vector_t v2); float innerProduct (vector_t v1, vector_t v2); +float distancePointPoint(const vector_t a, const vector_t b); float distancePointPlane (vector_t point, plane_t plane); +float distancePointLine (vector_t point, line_t line); +float distancePointLineOnPlane (vector_t point, + line_t line, plane_t plane); vector_t addVectors(vector_t v1, vector_t v2); vector_t subtractVectors(vector_t v1, vector_t v2); vector_t scaleVector(vector_t v1, float f); -vector_t averagePoints(points_t points); +/*vector_t averagePoints(points_t points); vector_t normalOfPoints(points_t points); -plane_t planeFromPoints(points_t points); +plane_t planeFromPoints(points_t points);*/ +plane_t planeFromThreePoints (vector_t a, vector_t b, vector_t c); +line_t lineFromTwoPoints (vector_t a, vector_t b); vector_t averageListedPoints(points_t points, list_t list); vector_t normalOfListedPoints(points_t points, list_t list); +vector_t directionOfListedPoints(points_t points, list_t list); plane_t planeFromListedPoints(points_t points, list_t list); +line_t lineFromListedPoints(points_t points, list_t list); void printPlane(plane_t p); +void printLine(line_t l); void printVector(vector_t v); - -#endif /* QHULL_VECTORS_H */ +#endif /* QHULL_VECTOR_H */ diff --git a/src/mtx_qhull/zhull.c b/src/mtx_qhull/zhull.c index 58f636e..5730160 100644 --- a/src/mtx_qhull/zhull.c +++ b/src/mtx_qhull/zhull.c @@ -1,52 +1,12 @@ #include "zhull.h" - -#include -#include - -#define TOL 1e-7 -#define MAXIT 1000000 - - - -vector_t averageListedPoints(const points_t points, const list_t list) { - index_t i; - vector_t m = initVector(0.0f, 0.0f, 0.0f); - for (i=0; icorners); - printList(facet->outsideset); - printList(facet->insideset); - printList(facet->neighbors);*/ + /* printf("deleting facet %li\n",facet); + printList(facet->corners); + printList(facet->outsideset); + printList(facet->insideset); + printList(facet->neighbors);*/ freeList(&(facet->corners)); freeList(&(facet->outsideset)); freeList(&(facet->insideset)); @@ -67,9 +27,9 @@ list_t appendNewFacets(zhull_t * const zh, const size_t num_facets) { new_facet->insideset=emptyList(); new_facet->maxdistance=0; new_facet->farthest_outside_point=0; -// printf("%li, %li\n", new_facet, (entry_t)new_facet); + // printf("%li, %li\n", new_facet, (entry_t)new_facet); setEntry(new_facets,i,(entry_t)new_facet); -// printf("created facet %li\n",new_facet); + // printf("created facet %li\n",new_facet); } appendListToList(&(zh->facets),new_facets); return new_facets; @@ -91,6 +51,15 @@ facet_t *getFacetByIndex(const list_t facets, return (facet_t*) getEntry(facets,index); } +index_t getTriangleCorner(const zhull_t * const zh, + const index_t triangle_idx, + const index_t corner_idx) { + if (triangle_idxfacets)) + return getEntry(getFacetByIndex(zh->facets,triangle_idx)->corners,corner_idx); + else + return 0; +} + facet_t *getFacetByPointer(const entry_t ptr) { return (facet_t*) ptr; } @@ -153,7 +122,7 @@ void freeFacets(zhull_t * const zh) { if (getLength(zh->facets)>0) { for (i=0; ifacets); i++) { f=getFacetByIndex(zh->facets,i); -// printf("deleting facet %li\n",i); + // printf("deleting facet %li\n",i); freeFacet(f); } freeList(&(zh->facets)); @@ -162,7 +131,7 @@ void freeFacets(zhull_t * const zh) { void freeZhull(zhull_t *zh) { if (zh!=0) { -// printf("free zhull\n"); + // printf("free zhull\n"); freeFacets(zh); freeList(&(zh->facets_with_insidepoints)); freeList(&(zh->facets_with_outsidepoints)); @@ -176,9 +145,15 @@ void freeZhull(zhull_t *zh) { // interface -vector_t getPoint(const zhull_t * const zh, - const index_t index) { - return (vector_t) zh->pts.v[index]; + + +line_t getLine(const zhull_t * const zh, + const facet_t * const f, + const index_t corner) { + vector_t a, b; + a=getPoint(zh->pts,getEntry(f->corners, (corner)%getLength(f->corners))); + b=getPoint(zh->pts,getEntry(f->corners,(corner+1)%getLength(f->corners))); + return lineFromTwoPoints(a,b); } zhull_t zhullInitPoints(const float *x, const float *y, @@ -195,68 +170,91 @@ zhull_t zhullInitPoints(const float *x, const float *y, void dividePointsBetweenNewFacets ( zhull_t * const zh, const list_t assoc, - const list_t new_facets, const list_t horizon_fcts) { + const list_t new_facets) { index_t i,j; - list_t new_and_horizon_fcts=mergeLists(new_facets,horizon_fcts); - printList(new_and_horizon_fcts); - uniquefyListEntries(&new_and_horizon_fcts); - printList(new_and_horizon_fcts); + facet_t *f; + list_t point_inside_facet_list=emptyList(); float d; for (i=0; iplane); - // printf("distance %5.2f to %d\n",d,findValueInList((entry_t)getFacetByIndex(new_and_horizon_fcts,j),zh->facets)); - /* - if ((d>=-TOL) && (d<=TOL)) { - appendToList(&(getFacetByIndex(new_facets,j)->insideset), - getEntry(assoc,i)); - if (notInList(getEntry(new_facets,j),zh->facets_with_insidepoints)) - appendToList(&(zh->facets_with_insidepoints), - getEntry(new_facets,j)); - break; - } - else */ - if (d>=-TOL) { - appendToList(&(getFacetByIndex(new_and_horizon_fcts,j)->outsideset), - getEntry(assoc,i)); - // printf("appended "); - if (notInList(getEntry(new_and_horizon_fcts,j),zh->facets_with_outsidepoints)) - appendToList(&(zh->facets_with_outsidepoints),getEntry(new_and_horizon_fcts,j)); - if (getFacetByIndex(new_and_horizon_fcts,j)->maxdistancemaxdistance=d; - getFacetByIndex(new_and_horizon_fcts,j)->farthest_outside_point=getEntry(assoc,i); - } + for (j=0; jpts,getEntry(assoc,i)), + f->plane); + if (d>=TOL_OUTSIDEPOINT) break; + else if (d>=TOL_INSIDEPOINT) { + appendToList(&point_inside_facet_list,(entry_t)f); + } + } + if (d>=TOL_OUTSIDEPOINT) { + appendToList(&(f->outsideset), getEntry(assoc,i)); + if (notInList((entry_t)f,zh->facets_with_outsidepoints)) + appendToList(&(zh->facets_with_outsidepoints),(entry_t)f); + if (f->maxdistancemaxdistance=d; + f->farthest_outside_point=getEntry(assoc,i); + } + } + else { + if (getLength(point_inside_facet_list)>0) { + for (j=0; jinsideset),getEntry(assoc,i)); + } + appendListToList(&(zh->facets_with_insidepoints),point_inside_facet_list); + uniquefyListEntries(&(zh->facets_with_insidepoints)); } } } + freeList(&point_inside_facet_list); } void zhullInitialFacets(zhull_t *zh) { list_t assoc = emptyList(); list_t new_facets = emptyList(); index_t i; - if (zh->pts.num_points >= 3) { - assoc = initListFromTo(0,zh->pts.num_points-1); + entry_t idx[3]={0,1,2}; + list_t list=initList(idx,3); + if (getNumPoints(zh->pts)>= 3) { + assoc = initListFromTo(0,getNumPoints(zh->pts)-1); new_facets = appendNewFacets(zh,2); if (getLength(new_facets)==2) { - getFacetByIndex(new_facets,0)->corners = initListFromTo(0,2); - getFacetByIndex(new_facets,1)->corners = initListFromTo(2,0); - getFacetByIndex(new_facets,0)->neighbors = - initConstantList((entry_t)getFacetByIndex(new_facets,1),3); - getFacetByIndex(new_facets,1)->neighbors = - initConstantList((entry_t)getFacetByIndex(new_facets,0),3); - for (i=0; i<2; i++) { - getFacetByIndex(new_facets,i)->plane = - planeFromListedPoints(zh->pts, getFacetByIndex(new_facets,i)->corners); - getFacetByIndex(new_facets,i)->outsideset = emptyList(); - getFacetByIndex(new_facets,i)->insideset = emptyList(); - getFacetByIndex(new_facets,i)->maxdistance = 0.0f; - removeValueListFromList(&assoc,getFacetByIndex(new_facets,i)->corners); + do { // circumvent coincident points + if (distancePointPoint(getPoint(zh->pts,idx[0]),getPoint(zh->pts,idx[1]))>TOL_DEGENERATE) + break; + else + idx[1]++; + } while (idx[1]pts)); + if (idx[1]pts)) { + do { // circumvent degenerate triangles + list=initList(idx,3); + if (lengthVector(normalOfListedPoints(zh->pts,list))>TOL_DEGENERATE) + break; + else { + idx[2]++; + freeList(&list); + } + } while (idx[2]pts)); + if (idx[2]pts)) { + getFacetByIndex(new_facets,0)->corners = list; + list=initList(idx,3); + reverseList(&list); + getFacetByIndex(new_facets,1)->corners = list; + getFacetByIndex(new_facets,0)->neighbors = + initConstantList((entry_t)getFacetByIndex(new_facets,1),3); + getFacetByIndex(new_facets,1)->neighbors = + initConstantList((entry_t)getFacetByIndex(new_facets,0),3); + for (i=0; i<2; i++) { + getFacetByIndex(new_facets,i)->plane = + planeFromListedPoints(zh->pts, getFacetByIndex(new_facets,i)->corners); + getFacetByIndex(new_facets,i)->outsideset = emptyList(); + getFacetByIndex(new_facets,i)->insideset = emptyList(); + getFacetByIndex(new_facets,i)->maxdistance = 0.0f; + removeValueListFromList(&assoc,getFacetByIndex(new_facets,i)->corners); + } + dividePointsBetweenNewFacets(zh, assoc, new_facets); + } } - dividePointsBetweenNewFacets(zh, assoc, new_facets, emptyList()); } freeList(&new_facets); freeList(&assoc); @@ -308,59 +306,78 @@ void sortHorizonEdges(list_t *horizon_fcts, void removeVisibleFacetsGetHorizonAndAvailablePoints( zhull_t * const zh, index_t point_index, - facet_t *facet, + facet_t *f, list_t *horizon_fcts, list_t *horizon_fcts_edges, list_t *avail_points) { index_t i,j,k; - facet_t *f, *n; + facet_t *n; float d; list_t visible_fcts = emptyList(); - list_t indices_for_printing = emptyList(); + list_t fcts_to_visit = emptyList(); + list_t fcts_visited = emptyList(); *avail_points = emptyList(); *horizon_fcts = emptyList(); *horizon_fcts_edges = emptyList(); - appendToList(&visible_fcts, (entry_t)facet); - for(i=0;ioutsideset); - appendListToList(avail_points, f->insideset); - for (j=0; jneighbors); j++) { - n=getFacetByIndex(f->neighbors,j); - if (notInList((entry_t) n, visible_fcts)) { - d=distancePointPlane(getPoint(zh,point_index), - n->plane); - if (d>-TOL) { - if (notInList((entry_t)n,visible_fcts)) - appendToList(&visible_fcts,(entry_t)n); - } + i=0; + appendToList(&fcts_to_visit,(entry_t)f); + if (getLength(fcts_to_visit)==0) + return; + do { + f=getFacetByIndex(fcts_to_visit,0); + d=distancePointPlane(getPoint(zh->pts,point_index),f->plane); + if (d>=TOL_OUTSIDEPOINT) { // visible facet + appendToList(&visible_fcts,(entry_t)f); + appendListToList(avail_points, f->outsideset); + for (j=0; jneighbors); j++) { + n=getFacetByIndex(f->neighbors,j); + d=distancePointPlane(getPoint(zh->pts,point_index),n->plane); + if (d>=TOL_INSIDEPOINT) + appendToList(&fcts_to_visit,(entry_t)n); else { - appendToList(horizon_fcts,(entry_t)n); + // invisible neighbor: horizon edge found k=getEntry(f->corners,(j+1)%getLength(f->corners)); - // printf("searching for corner with %d in\n",k); - // printList(n->corners); k=findValueInList(k,n->corners); - // printf("found %d\n",k); + appendToList(horizon_fcts,(entry_t)n); appendToList(horizon_fcts_edges,k); - } + } + } + } + else if (d>=TOL_INSIDEPOINT) { // coincident facet + for (j=0; jneighbors); j++) { + // is point outside polygon? searching edges. + d=distancePointLineOnPlane(getPoint(zh->pts,point_index), + getLine(zh,f,j),f->plane); + if (d>=TOL_OUTSIDEPOINT) { // point outside edge: register horizon edge + appendToList(horizon_fcts,(entry_t)f); + appendToList(horizon_fcts_edges,j); + // exceptional: initial facet with opposing normal + if (innerProduct(f->plane.normal, getFacetByIndex(f->neighbors,j)->plane.normal)< + -1.0f+TOL_OUTSIDEPOINT) + appendToList(&fcts_to_visit,(entry_t) + getFacetByIndex(f->neighbors,j)); + } // otherwise point is discared } } - } - printf("removing facets"); - indices_for_printing=findValueListInList(visible_fcts, zh->facets); - printList(indices_for_printing); + appendToList(&fcts_visited,(entry_t)f); + removeValueListFromList(&fcts_to_visit,fcts_visited); + } while (getLength(fcts_to_visit)>0); + +// printf("removing facets"); +// indices_for_printing=findValueListInList(visible_fcts, zh->facets); +// printList(indices_for_printing); +// freeList(&indices_for_printing); removeFacetByPointerList(zh,visible_fcts); freeList(&visible_fcts); - // printHorizonEdges(horizon_fcts, horizon_fcts_edges); + freeList(&fcts_to_visit); sortHorizonEdges(horizon_fcts, horizon_fcts_edges); - printHorizonEdges(horizon_fcts, horizon_fcts_edges); +// printHorizonEdges(horizon_fcts, horizon_fcts_edges); } - -void initNewFacets(zhull_t *zh, int point_index, +void initNewFacets(zhull_t *zh, index_t point_index, list_t new_facets, list_t horizon_fcts, list_t horizon_fcts_edges) { index_t i,j; @@ -368,8 +385,6 @@ void initNewFacets(zhull_t *zh, int point_index, array[0]=point_index; for (i=0; icorners=initList(array,3); - // printList(getFacetByIndex(new_facets,i)->corners); // neighbors // previous new neighbor @@ -391,6 +405,11 @@ void initNewFacets(zhull_t *zh, int point_index, getFacetByIndex(new_facets,i)->neighbors= initList(array,3); + removeValueFromList(&(getFacetByIndex(horizon_fcts,i)->insideset),point_index); + if (getLength(getFacetByIndex(horizon_fcts,i)->insideset)==0) { + removeValueFromList(&(zh->facets_with_insidepoints),(entry_t)getFacetByIndex(horizon_fcts,i)); + } + // registering at old neighbor: setEntry(getFacetByIndex(horizon_fcts,i)->neighbors, getEntry(horizon_fcts_edges,i), getEntry(new_facets,i)); @@ -399,6 +418,9 @@ void initNewFacets(zhull_t *zh, int point_index, getFacetByIndex(new_facets,i)->plane = planeFromListedPoints(zh->pts, getFacetByIndex(new_facets,i)->corners); + +// printf("new facet %d ",findValueInList(getEntry(new_facets,i),zh->facets)); +// printFacet(zh, getFacetByIndex(new_facets,i)); } } @@ -406,12 +428,13 @@ void makePyramidFacetsToHorizon(zhull_t *zh, index_t point_index, list_t horizon_fcts, list_t horizon_fcts_edges, list_t avail_points) { list_t new_facets = appendNewFacets(zh, getLength(horizon_fcts_edges)); +// printf("making new pyramid of %d facets\n",getLength(horizon_fcts_edges)); initNewFacets(zh,point_index,new_facets,horizon_fcts,horizon_fcts_edges); // printf("available points: "); // printList(avail_points); // printf("new facets : "); // printList(new_facets); - dividePointsBetweenNewFacets(zh, avail_points, new_facets, horizon_fcts); + dividePointsBetweenNewFacets(zh, avail_points, new_facets); freeList(&new_facets); } @@ -426,22 +449,32 @@ void calculateZHull(zhull_t *zh,int maxit) { list_t available_points=emptyList(); if (maxit>MAXIT) maxit=MAXIT; - if (zh->pts.num_points!=0){ + if (getNumPoints(zh->pts)!=0){ zhullInitialFacets(zh); - printZhull(zh); - while((getLength(zh->facets_with_outsidepoints)>0)&&(cnt++facets_with_outsidepoints); - f=getFacetByIndex(zh->facets_with_outsidepoints,fli); - pi=f->farthest_outside_point; +// printZhull(zh); + while(((getLength(zh->facets_with_insidepoints)>0) + ||(getLength(zh->facets_with_outsidepoints)>0)) + &&(cnt++facets_with_insidepoints)>0) { + fli%=getLength(zh->facets_with_insidepoints); + f=getFacetByIndex(zh->facets_with_insidepoints,fli); + pi=getEntry(f->insideset,0); +// printf("insidepoint\n"); +// printList(zh->facets_with_insidepoints); + } + else { + fli%=getLength(zh->facets_with_outsidepoints); + f=getFacetByIndex(zh->facets_with_outsidepoints,fli); + pi=f->farthest_outside_point; + } removeVisibleFacetsGetHorizonAndAvailablePoints(zh,pi,f, &horizon_fcts, &horizon_fcts_edges, &available_points); removeValueFromList(&available_points, pi); makePyramidFacetsToHorizon(zh,pi,horizon_fcts,horizon_fcts_edges, available_points); - printZhull(zh); - +// printZhull(zh); freeList(&horizon_fcts); freeList(&horizon_fcts_edges); freeList(&available_points); @@ -453,13 +486,15 @@ void calculateZHull(zhull_t *zh,int maxit) { void printZhull(const zhull_t * const zh) { index_t fi; list_t indices = emptyList(); - printf("zhull from %d points\n", zh->pts.num_points); + printf("zhull from %d points\n", getNumPoints(zh->pts)); printf("facets with outsidepoints: "); indices=findValueListInList(zh->facets_with_outsidepoints,zh->facets); printList(indices); freeList(&indices); -// printf("facets with insidepoints: "); -// printList(zh->facets_with_insidepoints); + printf("facets with insidepoints: "); + indices=findValueListInList(zh->facets_with_insidepoints,zh->facets); + printList(indices); + freeList(&indices); printf("zhull has %d facets\n", getLength(zh->facets)); for (fi=0; fifacets); fi++) { printf("facet %d<%d>: ",fi,getFacetByIndex(zh->facets,fi)); diff --git a/src/mtx_qhull/zhull.h b/src/mtx_qhull/zhull.h index 2a20768..fc0ff14 100644 --- a/src/mtx_qhull/zhull.h +++ b/src/mtx_qhull/zhull.h @@ -1,9 +1,16 @@ #ifndef QHULL_ZHULL_H #define QHULL_ZHULL_H -#include -#include "list.h" #include "vectors.h" +#include "list.h" +#include +#include +#include +#define TOL_OUTSIDEPOINT 1e-7 +#define TOL_INSIDEPOINT -1e-7 +#define TOL_DEGENERATE 1e-6 +#define MAXIT 1000000 + typedef struct facet_ { plane_t plane; @@ -23,13 +30,13 @@ typedef struct zhull_ { } zhull_t; void calculateZHull(zhull_t *zh,int maxit); +index_t getTriangleCorner(const zhull_t * const zh, + const index_t triangle_idx, + const index_t corner_idx); void printZhull(const zhull_t * const zh); void freeZhull(zhull_t *zh); zhull_t zhullInitPoints(const float *x, const float *y, const float *z, const size_t num_points); void printFacet(const zhull_t * const zh, const facet_t * const f); - -facet_t *getFacetByIndex(const list_t facets, const index_t index); - #endif /* QHULL_ZHULL_H */ -- cgit v1.2.1