improved deterministic behavior of [mtx_qhull], however planar polygons

do not fully work yet.

svn path=/trunk/externals/iem/iemmatrix/; revision=16172

9 files changed, 1066 insertions, 375 deletions

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 <stdlib.h>
 #include <stdio.h>
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 <stdio.h>
+#include <stdlib.h>
 
 #include <sys/types.h>
 
 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 <math.h>
-#include <stdlib.h>
-#include <stdio.h>
 
 
 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<points.num_points; i++)
-        m=addVectors(points.v[i],m);
-    m=scaleVector(m,1.0f/((float)points.num_points));
+    for (i=0; i<getNumPoints(points); i++)
+        m=addVectors(getPoint(points,m));
+    m=scaleVector(m,1.0f/((float)getNumPoints(points)));
     return m;
 }
-
 vector_t normalOfPoints(points_t points) {
     vector_t m = averagePoints(points);
     vector_t n = initVector(0.0f, 0.0f, 0.0f);
     size_t i;
-    for (i=0; i<points.num_points; i++) {
-        n=addVectors(crossProduct(points.v[i],m),n);
+    for (i=0; i<getNumPoints(points); i++) {
+        n=addVectors(crossProduct(getPoint(points,i),m),n);
     }
     return n;
 }
-
 plane_t planeFromPoints(points_t points) {
     vector_t p=averagePoints(points);
     vector_t n=normalOfPoints(points);
     plane_t pl=initPlane(n,p);
     return pl;
 }
+*/
+
+plane_t planeFromThreePoints(const vector_t a, 
+        const vector_t b, const vector_t c) {
+    vector_t ab = subtractVectors(b,a);
+    vector_t ac = subtractVectors(c,a);
+    vector_t n = normalizeVector(crossProduct(ab,ac));
+    return initPlane(n,a);
+}
+
+line_t lineFromTwoPoints(const vector_t a, 
+        const vector_t b) {
+    vector_t direction = subtractVectors(b,a);
+    return initLine(direction,a);
+}
+
+size_t getNumPoints(const points_t points) {
+    return points.num_points;
+}
+
+vector_t getPoint(const points_t points, 
+        const index_t index) {
+    if ((index>=0)&&(index<getNumPoints(points)))
+        return points.v[index];
+    return initVector(0.0f,0.0f,0.0f);
+}
+
+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; i<getLength(list); i++) {
+        m=addVectors(getPoint(points,getEntry(list,i)),m);
+    }
+    m=scaleVector(m,1.0f/((float)getLength(list)));
+    return m;
+}
+
+vector_t normalOfListedPoints(const points_t points, const list_t list) {
+    vector_t m = averageListedPoints(points,list);
+    vector_t n = initVector(0.0f, 0.0f, 0.0f);
+    vector_t d1,d2,c;
+    index_t i;
+    for (i=1; i<=getLength(list); i++) {
+        d1=subtractVectors(getPoint(points,getEntry(list,i-1)),m);
+        d2=subtractVectors(getPoint(points,getEntry(list,i%getLength(list))),m);
+        c=crossProduct(d1,d2);
+        n=addVectors(c,n);
+    }
+    return n;
+}
+
+vector_t directionOfListedPoints(const points_t points, const list_t list) {
+    vector_t m = averageListedPoints(points,list);
+    vector_t dir = initVector(0.0f, 0.0f, 0.0f);
+    vector_t d,c;
+    index_t i;
+    for (i=1; i<getLength(list); i++) {
+        d=subtractVectors(getPoint(points,getEntry(list,i-1)),
+                getPoint(points,getEntry(list,i)));
+        dir=addVectors(d,dir);
+    }
+    return dir;
+}
+
+plane_t planeFromListedPoints(const points_t points, const list_t list) {
+    vector_t p=averageListedPoints(points,list);
+    vector_t n=normalOfListedPoints(points,list);
+    plane_t pl=initPlane(n,p);
+    return pl;
+}
+
+
+
+line_t lineFromListedPoints(const points_t points, const list_t list) {
+    vector_t p=averageListedPoints(points,list);
+    vector_t n=directionOfListedPoints(points,list);
+    line_t l=initLine(n,p);
+    return l;
+}
+
+void printVector(vector_t v) {
+   printf("[%5.2f,%5.2f,%5.2f], ", v.c[0],v.c[1],v.c[2]);
+}
 
 void printPlane(plane_t p) {
    printf("n=");
@@ -125,7 +232,11 @@ void printPlane(plane_t p) {
    printVector(p.point);
 }
 
-void printVector(vector_t v) {
-   printf("[%5.2f,%5.2f,%5.2f], ", v.c[0],v.c[1],v.c[2]);
+void printLine(line_t l) {
+   printf("d=");
+   printVector(l.direction);
+   printf(", p=");
+   printVector(l.point);
 }
 
+
diff --git a/src/mtx_qhull/vectors.h b/src/mtx_qhull/vectors.h
index fadb370..06a302a 100644
--- a/src/mtx_qhull/vectors.h
+++ b/src/mtx_qhull/vectors.h
@@ -1,9 +1,9 @@
-#ifndef QHULL_VECTORS_H
-#define QHULL_VECTORS_H
-
-#include <sys/types.h>
+#ifndef QHULL_VECTOR_H
+#define QHULL_VECTOR_H 
 #include "list.h"
-
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
 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 <stdlib.h>
-#include <stdio.h>
-
-#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; i<getLength(list); i++) {
-        m=addVectors(points.v[getEntry(list,i)],m);
-    }
-    m=scaleVector(m,1.0f/((float)getLength(list)));
-    return m;
-}
-
-vector_t normalOfListedPoints(const points_t points, const list_t list) {
-    vector_t m = averageListedPoints(points,list);
-    vector_t n = initVector(0.0f, 0.0f, 0.0f);
-    vector_t d1,d2,c;
-    index_t i;
-    for (i=1; i<getLength(list); i++) {
-        d1=subtractVectors(points.v[getEntry(list,i-1)],m);
-        d2=subtractVectors(points.v[getEntry(list,i)],m);
-        c=crossProduct(d1,d2);
-        n=addVectors(c,n);
-    }
-    return n;
-}
-
-plane_t planeFromListedPoints(const points_t points, const list_t list) {
-    vector_t p=averageListedPoints(points,list);
-    vector_t n=normalOfListedPoints(points,list);
-    plane_t pl=initPlane(n,p);
-    return pl;
-}
-
 /* facets, memory things */
 
 void freeFacet(facet_t *facet) {
-/*    printf("deleting facet %li\n",facet);
-    printList(facet->corners);
-    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_idx<getLength(zh->facets)) 
+        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; i<getLength(zh->facets); 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; i<getLength(assoc); i++) {
-        printf("p %d\n", getEntry(assoc,i));
-        for (j=0; j<getLength(new_and_horizon_fcts); j++) {
-            d=distancePointPlane(getPoint(zh,getEntry(assoc,i)),
-                    getFacetByIndex(new_and_horizon_fcts,j)->plane); 
- //           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)->maxdistance<d) {
-                    getFacetByIndex(new_and_horizon_fcts,j)->maxdistance=d;
-                    getFacetByIndex(new_and_horizon_fcts,j)->farthest_outside_point=getEntry(assoc,i);
-                }
+        for (j=0; j<getLength(new_facets); j++) {
+            f=getFacetByIndex(new_facets,j);
+            d=distancePointPlane(getPoint(zh->pts,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->maxdistance<d) {
+                f->maxdistance=d;
+                f->farthest_outside_point=getEntry(assoc,i);
+            }
+        }
+        else {
+            if (getLength(point_inside_facet_list)>0) {
+                for (j=0; j<getLength(point_inside_facet_list); j++) {
+                    f=getFacetByIndex(point_inside_facet_list,j);
+                    appendToList(&(f->insideset),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]<getNumPoints(zh->pts));
+            if (idx[1]<getNumPoints(zh->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]<getNumPoints(zh->pts));
+                if (idx[2]<getNumPoints(zh->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;i<getLength(visible_fcts);i++) {
-        f=getFacetByIndex(visible_fcts,i);
-        appendListToList(avail_points, f->outsideset);
-        appendListToList(avail_points, f->insideset);
-        for (j=0; j<getLength(f->neighbors); 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; j<getLength(f->neighbors); 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; j<getLength(f->neighbors); 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; i<getLength(new_facets); i++) {
-        //        printf("adding facet %d\n",getEntry(new_facets,i));
-
         // corners
         getHorizonEdgeByIndex(array,horizon_fcts, horizon_fcts_edges,i);
         array[2]=array[1];
@@ -377,7 +392,6 @@ void initNewFacets(zhull_t *zh, int point_index,
         array[0]=array[2];
         array[2]=point_index;
         getFacetByIndex(new_facets,i)->corners=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++<maxit))  {
-            printf("//////////////// ITERATION %d ///////\n",cnt);
-            fli%=getLength(zh->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++<maxit))  {
+//            printf("//////////////// ITERATION %d ///////\n",cnt);
+            if (getLength(zh->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; fi<getLength(zh->facets); 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 <sys/types.h>
-#include "list.h"
 #include "vectors.h"
+#include "list.h"
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#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 */