convex hull algorithm added

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

10 files changed, 1404 insertions, 0 deletions

diff --git a/src/mtx_qhull/list.c b/src/mtx_qhull/list.c
new file mode 100644
index 0000000..07ce4f2
--- /dev/null
+++ b/src/mtx_qhull/list.c
@@ -0,0 +1,294 @@
+
+// memory things:
+list_t emptyList(void) {
+    list_t generated_list;
+    generated_list.length=0;
+    generated_list.entries=0;
+    return generated_list;
+}
+
+list_t allocateList(const size_t length) {
+    list_t generated_list = emptyList();
+    if (length>0) {
+        generated_list.entries= (entry_t*) malloc(length*sizeof(entry_t));
+        if (generated_list.entries!=0) {
+            generated_list.length=length;
+            //          printf("got list %li\n",generated_list.entries);
+        } 
+    }
+    return generated_list;
+}
+
+void reallocateList(list_t *list,
+        const size_t length) {
+    entry_t *old_entries = list->entries;
+    if (length>0) {
+        if (list->length==0) 
+            *list = allocateList(length);
+        else {
+            if (list->length != length)
+                list->entries = (entry_t*) realloc(list->entries,length*sizeof(entry_t));
+            if (list->entries!=0) 
+                list->length=length;
+            else 
+                *list=emptyList();
+        }
+    } 
+    else 
+        freeList(list);
+    /*    if ((list->entries!=old_entries)&&(old_entries!=0)) {
+          if (list->entries!=0)
+          printf("moved %li by realloc to %li\n",old_entries,list->entries);
+          else
+          printf("freed %li by realloc\n", old_entries);
+          }*/
+}
+
+
+void freeList(list_t *list) {
+    if (list->entries!=0) {
+        //        printf("deleting list %li\n",list->entries);
+        free(list->entries);
+    }
+    list->entries=0;
+    list->length=0;
+}
+
+// programming interface:
+
+size_t getLength(const list_t list) {
+    return list.length;
+}
+
+entry_t getEntry(const list_t list, const index_t index) {
+    if ((index>=0)&&(index<getLength(list)))
+        return list.entries[index];
+    else
+        return 0;
+}
+
+void setEntry(const list_t list, const index_t index, const entry_t entry) {
+    if ((index>=0)&&(index<getLength(list)))
+        list.entries[index]=entry;
+}
+
+
+list_t initList(entry_t *entries, const size_t length) {
+    index_t i;
+    list_t l = allocateList(length);
+    if (l.entries!=0) 
+        for (i=0; i<(index_t)length; i++) 
+            setEntry(l,i,entries[i]);
+    return l;
+}
+
+list_t initListFromTo(const entry_t start, const entry_t stop) {
+    index_t i;
+    size_t length;
+    entry_t c;
+    int incr;
+    if (stop>=start) {
+        length=(size_t) (stop-start+1);
+        incr=1;
+    } else {
+        length=(size_t) (start-stop+1);
+        incr=-1;
+    }
+    list_t l = allocateList(length);
+    if (l.entries!=0) 
+        for (i=0,c=start; i<length; i++, c+=incr) 
+            setEntry(l,i,c);
+    return l;
+}
+
+list_t initConstantList(const entry_t c, const size_t length){
+    entry_t i;
+    list_t l = allocateList(length);
+    if (l.entries!=0) 
+        for (i=0; i<length; i++)
+            setEntry(l,i,c);
+    return l;
+}
+
+
+list_t duplicateList(const list_t list_in) {
+    entry_t i;
+    list_t list_out=emptyList();
+    list_out = allocateList(getLength(list_in));
+    for (i=0; i<getLength(list_out); i++) 
+        setEntry(list_out,i,getEntry(list_in,i));
+    return list_out;
+}
+
+list_t mergeLists(const list_t list1, const list_t list2) {
+    list_t list_out;
+    entry_t i,j;
+    list_out = allocateList(getLength(list1)+ getLength(list2));
+    if (list_out.entries!=0) {
+        for (i=0; i<getLength(list1); i++) 
+            setEntry(list_out,i,getEntry(list1,i));
+        for (j=0; i<getLength(list_out); i++, j++) 
+            setEntry(list_out,i,getEntry(list2,j));
+    }
+    return list_out;
+}
+
+list_t getSubList(const list_t list, const list_t indices) {
+    index_t i;
+    list_t new_list = allocateList(getLength(indices));
+    for (i=0; i<getLength(new_list); i++) {
+        setEntry(new_list,i,getEntry(list,getEntry(indices,i)));
+    }
+}
+
+list_t getSubListFromTo(const list_t list, const index_t start, 
+        const index_t stop)  {
+    list_t new_list=emptyList();
+    index_t i,j;
+    int incr;
+    if ((start>0)&&(stop>0)&&(start<getLength(list))&&(stop<getLength(list))) {
+        if (start>stop) {
+            incr=-1;
+            new_list=allocateList(start-stop+1);
+        } else {
+            incr=1;
+            new_list=allocateList(start-stop+1);
+        }
+        for (j=start,i=0; i<getLength(new_list); i++, j+=incr) {
+            setEntry(new_list,i,getEntry(list,j));
+        }
+    }
+    return new_list;
+}
+
+void appendToList(list_t *list, const entry_t entry) {
+    const index_t i=(index_t)getLength(*list);
+    reallocateList(list,getLength(*list)+1);
+    if (getLength(*list)>(size_t)i) {
+        setEntry(*list,i,entry);
+    }
+}
+
+void removeEntryFromList(list_t *list, const index_t index) {
+    index_t i,j;
+    for (i=j=0; i<getLength(*list); i++) {
+        if (i!=index) 
+            setEntry(*list,j++,getEntry(*list,i));
+    }
+    reallocateList(list,j);
+}
+
+
+void removeValueFromList(list_t *list, const entry_t entry) {
+    index_t i,j;
+    for (j=i=0; i<getLength(*list); i++) {
+        if (getEntry(*list,i)!=entry) 
+            setEntry(*list, j++, getEntry(*list,i));
+    }
+    reallocateList(list,j);
+}
+
+void appendListToList(list_t *list1, const list_t list2) {
+    index_t i,j;
+    size_t siz_old = getLength(*list1);
+    siz_old=list1->length;
+    reallocateList(list1, getLength(*list1) + getLength(list2));
+    if (getLength(*list1)>siz_old) {
+        for (i=siz_old, j=0; i<list1->length; i++, j++) 
+            setEntry(*list1,i,getEntry(list2,j));
+    }
+}
+
+void removeEntryListFromList(list_t *list, const list_t indices) {
+    index_t i,j;
+    for (i=j=0; i<getLength(*list); i++)  
+        if (notInList(i,indices)) 
+            setEntry(*list, j++, getEntry(*list,i));
+    reallocateList(list,j);
+}
+
+void removeValueListFromList(list_t *list, const list_t excl_list) {
+    index_t i,j,k;
+    int keep;
+    for (j=i=0; i<getLength(*list); i++) {
+        keep=1;
+        for (k=0; k<getLength(excl_list); k++) 
+            keep=(keep)&&(getEntry(*list,i)!=getEntry(excl_list,k));
+        if (keep)
+            setEntry(*list, j++, getEntry(*list,i));
+    }
+    reallocateList(list,j);
+}
+
+void reverseList(list_t * const list) {
+    index_t i,j;
+    entry_t v;
+    const cnt = getLength(*list)/ 2;
+    if (cnt>0)
+        for (i=0, j=getLength(*list)-1; i<cnt; i++, j--) {
+            v=getEntry(*list,i);
+            setEntry(*list,i,getEntry(*list,j));
+            setEntry(*list,j,v);
+        }
+}
+
+int inList(const entry_t entry, const list_t list) {
+    index_t i;
+    for (i=0; i<getLength(list); i++) 
+        if (getEntry(list,i)==entry)
+            return 1;
+    return 0;
+}
+
+entry_t findValueInList(const entry_t entry, const list_t list) {
+    index_t i;
+    for (i=0; i<getLength(list); i++) 
+        if (entry==getEntry(list,i))
+            return i;
+    return getLength(list);  
+}
+
+void uniquefyListEntries(list_t *list)  {
+    index_t i,j,k;
+    int keep;
+    k=0;
+    for (j=0; j<getLength(*list); j++) {
+        keep=1;
+        for (i=0; (i<k)&&(keep); i++)  
+            keep = (keep)&&(list->entries[j]!=list->entries[i]);
+        if (keep) {
+            list->entries[i++]=list->entries[j];
+            k++;
+        }
+    }
+    reallocateList(list, k);
+}
+
+list_t findValueListInList(const list_t value_list,
+        const list_t list) {
+    list_t l=emptyList();
+    index_t i,j;
+    for (i=0; i<getLength(value_list); i++) 
+        appendToList(&l,findValueInList(
+                    getEntry(value_list,i),list));
+    return l;
+}
+
+int notInList(const entry_t entry, const list_t list) {
+    index_t i;
+    for (i=0; i<getLength(list); i++) 
+        if (getEntry(list,i)==entry)
+            return 0;
+    return 1;
+}
+
+void printList(list_t const list) {
+    entry_t i;
+    printf("[list]_%d=[",list.length);
+    if (getLength(list)>0)
+        printf("%d",getEntry(list,0));
+    for (i=1; i<list.length; i++) 
+        printf(", %d",getEntry(list,i));
+    printf("]\n");
+}
+
diff --git a/src/mtx_qhull/list.h b/src/mtx_qhull/list.h
new file mode 100644
index 0000000..699f9ed
--- /dev/null
+++ b/src/mtx_qhull/list.h
@@ -0,0 +1,38 @@
+
+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);
+
+// programming interface:
+size_t getLength(const list_t list);
+entry_t getEntry(const list_t list, const index_t index);
+void setEntry(const list_t list, const index_t index, const entry_t entry);
+list_t initList(entry_t *entries, const size_t length);
+list_t initListFromTo(const entry_t start, const entry_t stop);
+list_t initConstantList(const entry_t c, const size_t length);
+list_t duplicateList(const list_t list_in);
+list_t mergeLists(const list_t list1, const list_t list2);
+list_t getSubList(const list_t list, const list_t indices);
+list_t getSubListFromTo(const list_t list, const index_t start, 
+      const index_t stop);
+void appendToList(list_t *list, const entry_t entry);
+void removeValueFromList(list_t *list, const entry_t entry);
+void removeEntryFromList(list_t *list, const index_t index);
+void appendListToList(list_t *list1, const list_t list2);
+void removeValueListFromList(list_t *list, const list_t excl_list);
+void removeEntryListFromList(list_t *list, const list_t indices);
+void reverseList(list_t * const list);
+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 printList(const list_t list);
+
diff --git a/src/mtx_qhull/mtx_xprod.pd b/src/mtx_qhull/mtx_xprod.pd
new file mode 100755
index 0000000..08ac92c
--- /dev/null
+++ b/src/mtx_qhull/mtx_xprod.pd
@@ -0,0 +1,32 @@
+#N canvas 502 62 631 420 12;

+#X obj 73 53 inlet;

+#X obj 473 54 inlet;

+#X obj 73 81 t a a;

+#X obj 73 202 mtx_.*;

+#X obj 197 201 mtx_.*;

+#X obj 73 239 mtx_-;

+#X obj 73 287 outlet;

+#X text 180 34 mtx: cross product;

+#X text 291 206 matrix1: 3xN;

+#X text 291 227 matrix2: 3xN;

+#X text 291 252 ============;

+#X text 299 273 matrix: 3xN;

+#X text 407 206 [a1 a2 a3 ...];

+#X text 408 227 [b1 b2 b3 ...];

+#X text 410 272 [a1 x b1 \, a2 x b2 \, ...];

+#X obj 73 114 mtx_scroll -1;

+#X obj 197 114 mtx_scroll 1;

+#X obj 345 112 mtx_scroll 1;

+#X obj 473 113 mtx_scroll -1;

+#X connect 0 0 2 0;

+#X connect 1 0 17 0;

+#X connect 1 0 18 0;

+#X connect 2 0 15 0;

+#X connect 2 1 16 0;

+#X connect 3 0 5 0;

+#X connect 4 0 5 1;

+#X connect 5 0 6 0;

+#X connect 15 0 3 0;

+#X connect 16 0 4 0;

+#X connect 17 0 3 1;

+#X connect 18 0 4 1;

diff --git a/src/mtx_qhull/test.obj b/src/mtx_qhull/test.obj
new file mode 100644
index 0000000..4b4e602
--- /dev/null
+++ b/src/mtx_qhull/test.obj
@@ -0,0 +1,72 @@
+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
diff --git a/src/mtx_qhull/test_list.c b/src/mtx_qhull/test_list.c
new file mode 100644
index 0000000..4ddddde
--- /dev/null
+++ b/src/mtx_qhull/test_list.c
@@ -0,0 +1,94 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include "list.h"
+
+int main(char **argv, int argc) {
+   list_t l1=emptyList();
+   list_t l2=emptyList();
+   list_t l3=emptyList();
+   int x[]={0, 2, 4, 6};
+
+   printf("empty list:\n");
+   printList(l1);
+   freeList(&l1);
+
+   printf("constant list with 10 1 entries:\n");
+   l1=initConstantList(1,10);
+   printList(l1);
+   freeList(&l1);
+
+   printf("list from 2 to 3:\n");
+   l1=initListFromTo(2,3);
+   printList(l1);
+   printf("list from 7 to 1:\n");
+   l2=initListFromTo(7,1);
+   printList(l2);
+   printf("duplicate list from 2 to 3:\n");
+   l3=duplicateList(l1);
+   printList(l3);
+   freeList(&l3);
+
+   printf("merge list from 1 to 7 with list from 2 to 3:\n");
+   l3=mergeLists(l1,l2);
+   printList(l3);
+
+   printf("remove entry 6 list from list:\n");
+   removeValueFromList(&l3,6);
+   printList(l3);
+
+   printf("remove entries [2, 3] from list:\n");
+   removeValueListFromList(&l3,l1);
+   printList(l3);
+
+   printf("reverse list:\n");
+   reverseList (&l3);
+   printList(l3);
+
+   printf("append entry 8 to list:\n");
+   appendToList(&l3,8);
+   printList(l3);
+
+   printf("is 8 not in list?: %d\n",notInList(8,l3));
+   printf("is 3 not in list?: %d\n",notInList(3,l3));
+
+   printf("remove item 4 from list\n");
+   removeEntryFromList(&l3,4);
+   printList(l3);
+
+   printf("remove item 1 from list\n");
+   removeEntryFromList(&l3,1);
+   printList(l3);
+
+   printf("remove item 0 from list\n");
+   removeEntryFromList(&l3,0);
+   printList(l3);
+
+   printf("...taking a longer list\n");
+   printList(l2);
+   freeList(&l3);
+   l3=initList(x,4);
+   printf("removing index list ");
+   printList(l3);
+   removeEntryListFromList(&l2,l3);
+   printList(l2);
+   freeList(&l1);
+   freeList(&l2);
+   freeList(&l3);
+
+   l1=initListFromTo(1,3);
+   l2=initListFromTo(0,5);
+   printf("...taking a longer list\n");
+   printList(l1);
+   printf("finding indices of values in list");
+   printList(l2);
+   l3=findValueListInList(l2,l1);
+   printList(l3);
+   
+
+   freeList(&l1);
+   freeList(&l2);
+   freeList(&l3);
+   
+}
+
+// gcc list.c test_list.c && ./a.out
diff --git a/src/mtx_qhull/test_write_conv_hull_obj.pd b/src/mtx_qhull/test_write_conv_hull_obj.pd
new file mode 100644
index 0000000..d27370d
--- /dev/null
+++ b/src/mtx_qhull/test_write_conv_hull_obj.pd
@@ -0,0 +1,201 @@
+#N canvas 534 0 850 668 10;
+#X obj 595 205 convhull 0.001;
+#X obj 595 176 t a a;
+#N canvas 69 11 539 331 writeobject 0;
+#X obj 163 270 textfile;
+#X obj 163 161 mtx;
+#X obj 163 123 until;
+#X obj 200 123 f;
+#X obj 226 123 + 1;
+#X msg 225 102 1;
+#X obj 163 81 t a a b;
+#X obj 163 102 mtx_size;
+#X obj 143 32 inlet;
+#X obj 373 44 inlet;
+#X msg 199 142 row \$1;
+#X obj 163 238 list trim;
+#X msg 269 196 add \$1;
+#X msg 329 174 symbol v;
+#X obj 373 65 t a b;
+#X obj 143 54 t b a b;
+#X msg 45 230 write test.obj cr;
+#X msg 399 108 clear;
+#X obj 164 218 list prepend;
+#X msg 269 174 symbol f;
+#X connect 1 0 18 0;
+#X connect 2 0 3 0;
+#X connect 3 0 4 0;
+#X connect 3 0 10 0;
+#X connect 4 0 3 1;
+#X connect 5 0 3 1;
+#X connect 6 0 7 0;
+#X connect 6 1 1 1;
+#X connect 6 2 5 0;
+#X connect 7 0 2 0;
+#X connect 8 0 15 0;
+#X connect 9 0 14 0;
+#X connect 10 0 1 0;
+#X connect 11 0 0 0;
+#X connect 12 0 18 1;
+#X connect 13 0 12 0;
+#X connect 14 0 6 0;
+#X connect 14 1 13 0;
+#X connect 14 1 17 0;
+#X connect 15 0 16 0;
+#X connect 15 1 6 0;
+#X connect 15 2 19 0;
+#X connect 16 0 0 0;
+#X connect 17 0 0 0;
+#X connect 18 0 11 0;
+#X connect 19 0 12 0;
+#X restore 116 490 pd writeobject;
+#N canvas 0 22 450 300 gemwin 0;
+#X obj 132 136 gemwin;
+#X obj 67 89 outlet;
+#X obj 67 10 inlet;
+#X obj 67 41 route create;
+#X msg 67 70 set destroy;
+#X msg 142 68 set create;
+#X msg 132 112 create \, 1 \, lighting 1;
+#X msg 298 112 destroy \, reset;
+#X msg 256 79 color 0.78 0.8 1;
+#X obj 256 57 loadbang;
+#X connect 2 0 3 0;
+#X connect 3 0 4 0;
+#X connect 3 0 6 0;
+#X connect 3 1 5 0;
+#X connect 3 1 7 0;
+#X connect 4 0 1 0;
+#X connect 5 0 1 0;
+#X connect 6 0 0 0;
+#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
+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
+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 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 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;
diff --git a/src/mtx_qhull/vectors.c b/src/mtx_qhull/vectors.c
new file mode 100644
index 0000000..909543e
--- /dev/null
+++ b/src/mtx_qhull/vectors.c
@@ -0,0 +1,126 @@
+
+vector_t initVector (float x, float y, float z) {
+    vector_t vec={x, y, z};
+    return vec;
+}
+
+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]);
+    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;
+    plane.point = point;
+    plane.normal = normalizeVector(normal);
+    return plane;
+}
+
+points_t initPoints (const float *x,
+        const float *y, const float *z, 
+        size_t num_points) {
+    points_t points;
+    size_t i;
+
+    points.v = (vector_t*) malloc (sizeof(vector_t)*num_points);
+    if (points.v!=0) {
+//        printf("created %li\n",points.v);
+        points.num_points = num_points;
+        for (i=0; i<num_points; i++) {
+            points.v[i] = initVector(x[i],y[i],z[i]);
+        }
+    }
+    return points; 
+}
+
+void freePoints (points_t *points) {
+//    printf("deleting %li\n",points->v);
+    if (points!=0) {
+        if (points->v!=0)
+            free(points->v);
+        points->v = 0;
+        points->num_points = 0;
+    }
+}
+
+vector_t crossProduct (vector_t v1, vector_t v2) {
+    vector_t cp;
+    cp.c[0]= v1.c[1]*v2.c[2]-v1.c[2]*v2.c[1];
+    cp.c[1]=-v1.c[0]*v2.c[2]+v1.c[2]*v2.c[0];
+    cp.c[2]= v1.c[0]*v2.c[1]-v1.c[1]*v2.c[0];
+    return cp;
+}
+
+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 distancePointPlane (vector_t point, plane_t plane) {
+    return innerProduct(point, plane.normal) - 
+        innerProduct(plane.point, plane.normal);
+}
+
+vector_t addVectors(vector_t v1, vector_t v2) {
+    vector_t v3;
+    v3.c[0]=v1.c[0]+v2.c[0];
+    v3.c[1]=v1.c[1]+v2.c[1];
+    v3.c[2]=v1.c[2]+v2.c[2];
+    return v3;
+}
+
+vector_t subtractVectors(vector_t v1, vector_t v2) {
+    vector_t v3;
+    v3.c[0]=v1.c[0]-v2.c[0];
+    v3.c[1]=v1.c[1]-v2.c[1];
+    v3.c[2]=v1.c[2]-v2.c[2];
+    return v3;
+}
+
+vector_t scaleVector(vector_t v1, float f) {
+    vector_t v2;
+    v2.c[0]=f*v1.c[0];
+    v2.c[1]=f*v1.c[1];
+    v2.c[2]=f*v1.c[2];
+    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));
+    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);
+    }
+    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;
+}
+
+void printPlane(plane_t p) {
+   printf("n=");
+   printVector(p.normal);
+   printf(", 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]);
+}
+
diff --git a/src/mtx_qhull/vectors.h b/src/mtx_qhull/vectors.h
new file mode 100644
index 0000000..201205a
--- /dev/null
+++ b/src/mtx_qhull/vectors.h
@@ -0,0 +1,40 @@
+typedef struct vec_ {
+    float c[3];
+} vector_t;
+
+typedef struct points_ {
+    vector_t *v;
+    size_t num_points;
+} points_t;
+
+typedef struct plane_ {
+    vector_t normal;
+    vector_t point;
+} plane_t;
+
+vector_t initVector (float x, float y, float z);
+vector_t normalizeVector(vector_t v);
+
+plane_t initPlane (vector_t normal, vector_t point);
+
+points_t initPoints (const float *x,
+        const float *y, const float *z, 
+        size_t num_points);
+void freePoints (points_t *points);
+
+vector_t crossProduct (vector_t v1, vector_t v2);
+float innerProduct (vector_t v1, vector_t v2);
+float distancePointPlane (vector_t point, 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 normalOfPoints(points_t points);
+plane_t planeFromPoints(points_t points);
+
+vector_t averageListedPoints(points_t points, list_t list);
+vector_t normalOfListedPoints(points_t points, list_t list);
+plane_t planeFromListedPoints(points_t points, list_t list);
+
+void printPlane(plane_t p);
+void printVector(vector_t v);
diff --git a/src/mtx_qhull/zhull.c b/src/mtx_qhull/zhull.c
new file mode 100644
index 0000000..d13260e
--- /dev/null
+++ b/src/mtx_qhull/zhull.c
@@ -0,0 +1,482 @@
+#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);*/
+    freeList(&(facet->corners));
+    freeList(&(facet->outsideset));
+    freeList(&(facet->insideset));
+    freeList(&(facet->neighbors));
+}
+
+list_t appendNewFacets(zhull_t * const zh, const size_t num_facets) {
+    facet_t *new_facet;
+    list_t new_facets=initConstantList(0,num_facets);
+    index_t i;
+    for (i=0; i<getLength(new_facets); i++) {
+        new_facet = (facet_t*) malloc(sizeof(facet_t));
+        if (new_facet==0)
+            break;
+        new_facet->neighbors=emptyList();
+        new_facet->corners=emptyList();
+        new_facet->outsideset=emptyList();
+        new_facet->insideset=emptyList();
+        new_facet->maxdistance=0;
+        new_facet->farthest_outside_point=0;
+//        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);
+    }
+    appendListToList(&(zh->facets),new_facets);
+    return new_facets;
+}
+
+/* facets, interface */
+/*entry_t getFacetCornerByIndex(const facet_t *facet,
+  index_t index_corner) {
+  if (facet!=0)
+  return getEntry(facet->corners, 
+  index_corner%getLength(facet->corners));
+  else
+  return 0;
+  }
+  */
+
+facet_t *getFacetByIndex(const list_t facets, 
+        const index_t index) {
+    return (facet_t*) getEntry(facets,index);
+}
+
+facet_t *getFacetByPointer(const entry_t ptr) {
+    return (facet_t*) ptr;
+}
+
+void getHorizonEdgeByIndex(index_t *corners,
+        const list_t horizon_fcts, const list_t horizon_fcts_edges, 
+        const index_t index) {
+    index_t i=(index+getLength(horizon_fcts_edges))
+        %getLength(horizon_fcts_edges);
+    facet_t *f = getFacetByIndex(horizon_fcts,i);
+    index_t j=getEntry(horizon_fcts_edges,i);
+
+    corners[0]=getEntry(f->corners, j);
+    if (getLength(f->corners)>0) 
+        corners[1]=getEntry(f->corners, (j+1) % getLength(f->corners));
+    else 
+        corners[1]=corners[0];
+}
+
+
+void removeFacetByPointer(zhull_t * const zh, 
+        facet_t * const pointer) {
+    removeValueFromList(&(zh->facets), (entry_t)pointer);
+    removeValueFromList(&(zh->facets_with_outsidepoints), 
+            (entry_t)pointer);
+    removeValueFromList(&(zh->facets_with_insidepoints), 
+            (entry_t)pointer);
+    freeFacet((facet_t*)pointer);
+} 
+
+void removeFacetByPointerList(zhull_t * const zh, 
+        const list_t pointers) {
+    index_t i;
+    for (i=0; i<getLength(pointers); i++) {
+        removeFacetByPointer(zh,
+                getFacetByIndex(pointers,i));
+    }
+} 
+
+void removeFacetByIndex(zhull_t * const zh, 
+        const index_t index) {
+    removeFacetByPointer(zh,
+            getFacetByIndex(zh->facets, index));
+} 
+
+void removeFacetByIndexList(zhull_t * const zh, 
+        const list_t indices) {
+    facet_t *f; 
+    index_t i;
+    for (i=0; i<getLength(indices); i++) {
+        f = getFacetByIndex(zh->facets, getEntry(indices,i));
+        removeFacetByPointer(zh,f);
+    }
+} 
+
+
+void freeFacets(zhull_t * const zh) {
+    int i;
+    facet_t *f;
+    if (getLength(zh->facets)>0) {
+        for (i=0; i<getLength(zh->facets); i++) {
+            f=getFacetByIndex(zh->facets,i);
+//            printf("deleting facet %li\n",i);
+            freeFacet(f);
+        }
+        freeList(&(zh->facets));
+    }
+}
+
+void freeZhull(zhull_t *zh) {
+    if (zh!=0) {
+//        printf("free zhull\n");
+        freeFacets(zh);
+        freeList(&(zh->facets_with_insidepoints));
+        freeList(&(zh->facets_with_outsidepoints));
+        freePoints(&(zh->pts));
+    }
+}
+
+
+// ***********************************
+//
+// interface
+
+
+vector_t getPoint(const zhull_t * const zh, 
+        const index_t index) {
+    return (vector_t) zh->pts.v[index];
+}
+
+zhull_t zhullInitPoints(const float *x, const float *y, 
+        const float *z, const size_t num_points) {
+    zhull_t zh;
+    zh.pts=initPoints(x,y,z,num_points);
+    zh.facets=emptyList();
+    zh.facets_with_outsidepoints=emptyList();
+    zh.facets_with_insidepoints=emptyList();
+    return zh;
+}
+
+
+
+void dividePointsBetweenNewFacets (
+        zhull_t * const zh, const list_t assoc, 
+        const list_t new_facets, const list_t horizon_fcts) {
+    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);
+    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);
+                }
+                break;
+            }
+        }
+    }
+}
+
+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);
+        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);
+            }
+            dividePointsBetweenNewFacets(zh, assoc, new_facets, emptyList());
+        }
+        freeList(&new_facets);
+        freeList(&assoc); 
+    }
+}
+
+void printHorizonEdges(list_t *horizon_fcts, 
+        list_t *horizon_fcts_edges) {
+    index_t i;
+    index_t c1[2];
+
+    printf("horizon edges: ");
+    for (i=0; i<getLength(*horizon_fcts_edges); i++) {
+        getHorizonEdgeByIndex(c1, *horizon_fcts, *horizon_fcts_edges, i);
+        printf(", %d->%d",c1[0],c1[1]);
+    }
+    printf("\n");
+}
+
+
+void sortHorizonEdges(list_t *horizon_fcts, 
+        list_t *horizon_fcts_edges) {
+    index_t i,j;
+    facet_t *fi;
+    index_t ei;
+    index_t c1[2];
+    index_t c2[2];
+
+    if (getLength(*horizon_fcts_edges)==0)
+        return;
+    for (i=0; i<getLength(*horizon_fcts_edges)-1; i++) {
+        getHorizonEdgeByIndex(c1, *horizon_fcts, *horizon_fcts_edges, i);
+        for (j=i+1; j<getLength(*horizon_fcts_edges); j++) {
+            getHorizonEdgeByIndex(c2, *horizon_fcts, *horizon_fcts_edges, j);
+            if (c1[1]==c2[0]) { //found edge continuation: swap positions
+                ei=getEntry(*horizon_fcts,j);
+                setEntry(*horizon_fcts,j,getEntry(*horizon_fcts,i+1));
+                setEntry(*horizon_fcts,i+1,ei);
+                ei=getEntry(*horizon_fcts_edges,j);
+                setEntry(*horizon_fcts_edges,j,getEntry(*horizon_fcts_edges,i+1));
+                setEntry(*horizon_fcts_edges,i+1,ei);
+                break;
+            }
+            // HIER MUSS MEHR LOGIK HINEIN: HORIZONT KANN
+            // FEHLER BEINHALTEN
+        }
+    }
+}
+
+void removeVisibleFacetsGetHorizonAndAvailablePoints(
+        zhull_t * const zh, index_t point_index, 
+        facet_t *facet,
+        list_t *horizon_fcts, 
+        list_t *horizon_fcts_edges,
+        list_t *avail_points) {
+
+    index_t i,j,k;
+    facet_t *f, *n;
+    float d;
+    list_t visible_fcts = emptyList();
+    list_t indices_for_printing = 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);
+                }
+                else {
+                    appendToList(horizon_fcts,(entry_t)n);
+                    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_edges,k);
+                }
+            }
+        }
+    }
+    printf("removing facets");
+    indices_for_printing=findValueListInList(visible_fcts, zh->facets);
+    printList(indices_for_printing);
+    removeFacetByPointerList(zh,visible_fcts);
+    freeList(&visible_fcts);
+    //    printHorizonEdges(horizon_fcts, horizon_fcts_edges);
+    sortHorizonEdges(horizon_fcts, horizon_fcts_edges);
+    printHorizonEdges(horizon_fcts, horizon_fcts_edges);
+}
+
+
+void initNewFacets(zhull_t *zh, int point_index,
+        list_t new_facets, list_t horizon_fcts,
+        list_t horizon_fcts_edges) {
+    index_t i,j;
+    entry_t array[3];
+
+    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];
+        array[1]=array[0];
+        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
+        j=(getLength(horizon_fcts)+i-1) % getLength(horizon_fcts);
+        array[1]=getEntry(new_facets,j);
+        // next new neighbor
+        j=(i+1) % getLength(horizon_fcts);
+        array[2]=getEntry(new_facets,j);
+        // old neighbor
+        array[0]=getEntry(horizon_fcts,i);
+        getFacetByIndex(new_facets,i)->neighbors=
+            initList(array,3);
+
+        // registering at old neighbor:
+        setEntry(getFacetByIndex(horizon_fcts,i)->neighbors,
+                getEntry(horizon_fcts_edges,i), getEntry(new_facets,i));
+
+        // initializing normal vectors and lists
+        getFacetByIndex(new_facets,i)->plane = 
+            planeFromListedPoints(zh->pts, 
+                    getFacetByIndex(new_facets,i)->corners);
+    }
+}
+
+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));
+    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);
+    freeList(&new_facets);
+}
+
+void calculateZHull(zhull_t *zh,int maxit) {
+    index_t fli=0;
+    index_t pi;
+    facet_t *f;
+    list_t outsideset;
+    int cnt=0;
+    list_t horizon_fcts=emptyList();
+    list_t horizon_fcts_edges=emptyList();
+    list_t available_points=emptyList();
+    if (maxit>MAXIT)
+        maxit=MAXIT;
+    if (zh->pts.num_points!=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;
+            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);
+
+            freeList(&horizon_fcts);
+            freeList(&horizon_fcts_edges);
+            freeList(&available_points);
+            fli++;
+        }
+    }
+}
+
+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("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("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));
+        printFacet(zh,getFacetByIndex(zh->facets,fi));
+    }
+}
+
+void printFacet(const zhull_t * const zh, 
+        const facet_t * const f) {
+    list_t indices=emptyList();
+    indices=findValueListInList(f->neighbors,zh->facets);
+    printf("plane: ");
+    printPlane(f->plane);
+    printf("\n");
+    printf("corners: ");
+    printList(f->corners);
+    printf("outsideset: ");
+    printList(f->outsideset);
+    printf("insideset: ");
+    printList(f->insideset);
+    printf("neighbors: ");
+    printList(indices);
+    freeList(&indices);
+    printf("pt %d with maxdist %5.2f\n",f->farthest_outside_point, f->maxdistance);
+}
+
diff --git a/src/mtx_qhull/zhull.h b/src/mtx_qhull/zhull.h
new file mode 100644
index 0000000..dd520e1
--- /dev/null
+++ b/src/mtx_qhull/zhull.h
@@ -0,0 +1,25 @@
+
+typedef struct facet_ {
+    plane_t plane;
+    list_t corners;
+    list_t outsideset;
+    list_t insideset;
+    size_t farthest_outside_point;
+    list_t neighbors;
+    float maxdistance;
+} facet_t;
+
+typedef struct zhull_ {
+    points_t pts;
+    list_t facets;
+    list_t facets_with_outsidepoints;
+    list_t facets_with_insidepoints;
+} zhull_t;
+
+void calculateZHull(zhull_t *zh,int maxit);
+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);