aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/iemmatrix_sources.c1
-rw-r--r--src/iemmatrix_sources.h1
-rw-r--r--src/mtx_qhull.c167
-rw-r--r--src/mtx_qhull/list.c294
-rw-r--r--src/mtx_qhull/list.h38
-rwxr-xr-xsrc/mtx_qhull/mtx_xprod.pd32
-rw-r--r--src/mtx_qhull/test.obj72
-rw-r--r--src/mtx_qhull/test_list.c94
-rw-r--r--src/mtx_qhull/test_write_conv_hull_obj.pd201
-rw-r--r--src/mtx_qhull/vectors.c126
-rw-r--r--src/mtx_qhull/vectors.h40
-rw-r--r--src/mtx_qhull/zhull.c482
-rw-r--r--src/mtx_qhull/zhull.h25
13 files changed, 1573 insertions, 0 deletions
diff --git a/src/iemmatrix_sources.c b/src/iemmatrix_sources.c
index 11bb4d4..0c739b6 100644
--- a/src/iemmatrix_sources.c
+++ b/src/iemmatrix_sources.c
@@ -72,6 +72,7 @@ void iemmatrix_sources_setup(void)
iemtx_powtodb_setup(); /* mtx_powtodb.c */
iemtx_print_setup(); /* mtx_print.c */
iemtx_prod_setup(); /* mtx_prod.c */
+ iemtx_qhull_setup(); /* mtx_qhull.c */
iemtx_qr_setup(); /* mtx_qr.c */
iemtx_rand_setup(); /* mtx_rand.c */
iemtx_repmat_setup(); /* mtx_repmat.c */
diff --git a/src/iemmatrix_sources.h b/src/iemmatrix_sources.h
index 0e24529..872be3a 100644
--- a/src/iemmatrix_sources.h
+++ b/src/iemmatrix_sources.h
@@ -70,6 +70,7 @@ void iemtx_pow_setup(void); /* mtx_pow.c */
void iemtx_powtodb_setup(void); /* mtx_powtodb.c */
void iemtx_print_setup(void); /* mtx_print.c */
void iemtx_prod_setup(void); /* mtx_prod.c */
+void iemtx_qhull_setup(void); /* mtx_qhull.c */
void iemtx_qr_setup(void); /* mtx_qr.c */
void iemtx_rand_setup(void); /* mtx_rand.c */
void iemtx_repmat_setup(void); /* mtx_repmat.c */
diff --git a/src/mtx_qhull.c b/src/mtx_qhull.c
new file mode 100644
index 0000000..700821f
--- /dev/null
+++ b/src/mtx_qhull.c
@@ -0,0 +1,167 @@
+/*
+ * iemmatrix
+ *
+ * objects for manipulating simple matrices
+ * own qhull algorithm implementation
+ *
+ * Copyright (c) 2012, Franz Zotter
+ * IEM, Graz, Austria
+ *
+ * For information on usage and redistribution, and for a DISCLAIMER OF ALL
+ * WARRANTIES, see the file, "LICENSE.txt," in this distribution.
+ *
+ */
+
+
+#include "iemmatrix.h"
+#include <stdlib.h>
+#include "mtx_qhull/list.h"
+#include "mtx_qhull/vectors.h"
+#include "mtx_qhull/zhull.h"
+#include "mtx_qhull/list.c"
+#include "mtx_qhull/vectors.c"
+#include "mtx_qhull/zhull.c"
+
+static t_class *mtx_qhull_class;
+
+typedef struct MTXQhull_ MTXQhull;
+struct MTXQhull_ {
+ t_object x_obj;
+ t_outlet *outl;
+ t_atom *list;
+ size_t size;
+ size_t hull_size;
+ int iter;
+ zhull_t *zh;
+};
+
+static void deleteMTXQhull(MTXQhull *xo) {
+ if (xo->zh!=0) {
+ free(xo->zh);
+ xo->zh=0;
+ }
+ if (xo->list!=0) {
+ free(xo->list);
+ xo->list=0;
+ xo->size=0;
+ }
+}
+
+void mTXQhullSetIterations(MTXQhull *xo, t_float f) {
+ xo->iter=(int)f;
+ if (xo->iter<0)
+ xo->iter=0;
+}
+
+static void *newMTXQhull(t_symbol *s, int argc, t_atom *argv) {
+ int nmax;
+ MTXQhull *xo = (MTXQhull *) pd_new (mtx_qhull_class);
+ xo->outl = outlet_new (&xo->x_obj, gensym("matrix"));
+ xo->zh=0;
+ xo->list=0;
+ xo->size=0;
+ xo->iter=0;
+ if (argc>0)
+ mTXQhullSetIterations(xo,atom_getfloat(argv));
+ return ((void *) xo);
+}
+
+
+static void mTXQhullMatrix(MTXQhull *xo, t_symbol *s,
+ int argc, t_atom *argv) {
+ int i;
+ int rows=atom_getint(argv++);
+ int columns=atom_getint(argv++);
+ int size=rows*columns;
+ int in_size=argc-2;
+ float *x;
+ float *y;
+ float *z;
+
+ /* size check */
+ if (!size)
+ post("mtx_qhull: invalid dimensions");
+ else if (in_size<size)
+ post("mtx_qhull: sparse matrix not yet supported: use \"mtx_check\"");
+ else if ((rows<4)||(columns!=3))
+ post("mtx_qhull: requires an L x 3 matrix with at least L>=4");
+ else {
+ xo->zh = (zhull_t*)malloc(sizeof(zhull_t));
+ x=(float*)malloc(rows*sizeof(float));
+ y=(float*)malloc(rows*sizeof(float));
+ z=(float*)malloc(rows*sizeof(float));
+ if ((x!=0)&&(y!=0)&&(z!=0)&&(xo->zh!=0)) {
+ for (i=0; i<rows; i++) {
+ x[i]=atom_getfloat(argv++);
+ y[i]=atom_getfloat(argv++);
+ z[i]=atom_getfloat(argv++);
+ }
+ *(xo->zh)=zhullInitPoints(x,y,z,rows);
+ calculateZHull(xo->zh,xo->iter);
+ free(x);
+ free(y);
+ free(z);
+ xo->hull_size=getLength(xo->zh->facets);
+ if (xo->list==0) {
+ xo->list = (t_atom*)malloc((xo->hull_size*3+2)*sizeof(t_atom));
+ } else {
+ xo->list = (t_atom*)realloc(xo->list,(xo->hull_size*3+2)*sizeof(t_atom));
+ }
+ if(xo->list!=0) {
+ xo->size=(xo->hull_size*3+2);
+ SETFLOAT(xo->list,(float)xo->hull_size);
+ SETFLOAT(xo->list+1,(float)3);
+ for (i=0; i<xo->hull_size; i++) {
+ SETFLOAT(xo->list+2+3*i, (float)getEntry(
+ getFacetByIndex(xo->zh->facets,i)->corners,0)+1);
+ SETFLOAT(xo->list+3+3*i, (float)getEntry(
+ getFacetByIndex(xo->zh->facets,i)->corners,1)+1);
+ SETFLOAT(xo->list+4+3*i, (float)getEntry(
+ getFacetByIndex(xo->zh->facets,i)->corners,2)+1);
+ }
+ outlet_anything(xo->outl, gensym("matrix"),
+ xo->size, xo->list);
+ freeZhull(xo->zh);
+ free(xo->zh);
+ xo->zh=0;
+ }
+ else {
+ post("mtx_qhull: memory problem, no operation!");
+ xo->size=0;
+ freeZhull(xo->zh);
+ free(xo->zh);
+ xo->zh=0;
+ }
+ }
+ else {
+ if(x!=0)
+ free(x);
+ if(y!=0)
+ free(y);
+ if(z!=0)
+ free(z);
+ if(xo->zh!=0)
+ free(xo->zh);
+ xo->zh=0;
+ post("mtx_qhull: memory error, no operation!");
+ }
+
+ }
+}
+
+void mtx_qhull_setup (void) {
+ mtx_qhull_class = class_new (
+ gensym("mtx_qhull"),
+ (t_newmethod) newMTXQhull,
+ (t_method) deleteMTXQhull,
+ sizeof(MTXQhull),
+ CLASS_DEFAULT, A_GIMME, 0);
+ class_addmethod(mtx_qhull_class, (t_method) mTXQhullMatrix, gensym("matrix"), A_GIMME, 0);
+ class_addfloat(mtx_qhull_class, (t_method) mTXQhullSetIterations);
+}
+
+void iemtx_qhull_setup(void){
+ mtx_qhull_setup();
+}
+
+
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);
generated by cgit v1.2.1 (git 2.18.0) at 2024-12-25 18:33:46 +0000