From 291c783983b64f4b0ae4f6f67a016233193efa8c Mon Sep 17 00:00:00 2001 From: Franz Zotter Date: Wed, 29 Aug 2012 06:58:22 +0000 Subject: cleaned up code for mtx_qhull svn path=/trunk/externals/iem/iemmatrix/; revision=16183 --- src/mtx_qhull.c | 21 ++-- src/mtx_qhull/entry.h | 14 +++ src/mtx_qhull/list.c | 54 +++++---- src/mtx_qhull/list.h | 14 +++ src/mtx_qhull/test_write_conv_hull_obj.pd | 176 +++++++++++++++++++----------- src/mtx_qhull/vectors.c | 91 ++++++++++----- src/mtx_qhull/vectors.h | 60 +++++----- src/mtx_qhull/zhull.c | 95 ++++++++++++---- src/mtx_qhull/zhull.h | 22 +++- 9 files changed, 377 insertions(+), 170 deletions(-) (limited to 'src') diff --git a/src/mtx_qhull.c b/src/mtx_qhull.c index 34a4891..9dbd05c 100644 --- a/src/mtx_qhull.c +++ b/src/mtx_qhull.c @@ -23,10 +23,11 @@ typedef struct MTXQhull_ MTXQhull; struct MTXQhull_ { t_object x_obj; t_outlet *outl; + t_outlet *outl_fl; t_atom *list; size_t size; size_t hull_size; - int iter; +// int iter; zhull_t *zh; }; @@ -42,22 +43,23 @@ static void deleteMTXQhull(MTXQhull *xo) { } } -void mTXQhullSetIterations(MTXQhull *xo, t_float f) { +/*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->outl_fl = outlet_new (&xo->x_obj, gensym("float")); xo->zh=0; xo->list=0; xo->size=0; - xo->iter=0; - if (argc>0) - mTXQhullSetIterations(xo,atom_getfloat(argv)); +// xo->iter=0; +// if (argc>0) +// mTXQhullSetIterations(xo,atom_getfloat(argv)); return ((void *) xo); } @@ -92,7 +94,8 @@ static void mTXQhullMatrix(MTXQhull *xo, t_symbol *s, z[i]=atom_getfloat(argv++); } *(xo->zh)=zhullInitPoints(x,y,z,rows); - calculateZHull(xo->zh,xo->iter); + i=calculateZHull(xo->zh); + outlet_float(xo->outl_fl, (float)i); free(x); free(y); free(z); @@ -141,7 +144,7 @@ static void mTXQhullMatrix(MTXQhull *xo, t_symbol *s, } } -void mtx_qhull_setup (void) { +static void mtx_qhull_setup (void) { mtx_qhull_class = class_new ( gensym("mtx_qhull"), (t_newmethod) newMTXQhull, @@ -149,7 +152,7 @@ void mtx_qhull_setup (void) { 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); +// class_addfloat(mtx_qhull_class, (t_method) mTXQhullSetIterations); } void iemtx_qhull_setup(void){ diff --git a/src/mtx_qhull/entry.h b/src/mtx_qhull/entry.h index bc22f35..d85da9a 100644 --- a/src/mtx_qhull/entry.h +++ b/src/mtx_qhull/entry.h @@ -4,6 +4,20 @@ #include #include + +/* + * variable entry types for + * list operations in zhull + * + * Copyright (c) 2012, IOhannes zmoelnig, + * with friendly help from + * IEM, Graz, Austria + * + * + */ + + + typedef size_t index_t; typedef union { diff --git a/src/mtx_qhull/list.c b/src/mtx_qhull/list.c index 6f3b416..cf2b797 100644 --- a/src/mtx_qhull/list.c +++ b/src/mtx_qhull/list.c @@ -2,6 +2,20 @@ #include #include +/* + * list operations for zhull + * + * Copyright (c) 2012, Franz Zotter, + * with friendly help from + * IOhannes zmoelnig + * for variable entry types + * in entry.h + * IEM, Graz, Austria + * + * + */ + + // memory things: list_t emptyList(void) { list_t generated_list; @@ -16,7 +30,6 @@ list_t allocateList(const size_t length) { 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; @@ -26,7 +39,7 @@ void reallocateList(list_t *list, const size_t length) { entry_t *old_entries = list->entries; if (length>0) { - if (list->length==0) + if (getLength(*list)==0) *list = allocateList(length); else { if (list->length != length) @@ -39,18 +52,11 @@ void reallocateList(list_t *list, } 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; @@ -72,24 +78,27 @@ 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) { +void setEntry(const list_t list, const index_t index, + const entry_t entry) { if ((index>=0)&&(index=getLength(list1)) { for (i=0; i(size_t)i) { + if (getLength(*list)>i) { setEntry(*list,i,entry); } } @@ -213,13 +222,10 @@ void removeValueFromList(list_t *list, const entry_t entry) { void appendListToList(list_t *list1, const list_t list2) { index_t i,j; - size_t siz_old = getLength(*list1); - siz_old=list1->length; + const size_t siz_old = getLength(*list1); reallocateList(list1, getLength(*list1) + getLength(list2)); - if (getLength(*list1)>siz_old) { - for (i=siz_old, j=0; ilength; i++, j++) - setEntry(*list1,i,getEntry(list2,j)); - } + for (i=siz_old, j=0; i0)&&(points!=0)) { + if (getNumPoints(*points)==0) + *points=allocatePoints(num_points); + else { + points->v = (vector_t *) realloc(points->v,sizeof(vector_t)*num_points); + if (points->v!=0) + points->num_points=num_points; + else + points->num_points=0; + } + if (points->v!=0) + points->num_points = num_points; + } + else + freePoints(points); +} + +void appendPoints(points_t *points, + const float *x, const float *y, + const float *z, const size_t num_points) { + const size_t n=getNumPoints(*points); + size_t i,j; + reallocatePoints(points,getNumPoints(*points)+num_points); + for (i=n,j=0; iv[i] = initVector(x[j],y[j],z[j]); + } +} + +vector_t crossProduct (const vector_t v1, const 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]; @@ -66,7 +107,7 @@ vector_t crossProduct (vector_t v1, vector_t v2) { return cp; } -float innerProduct (vector_t v1, vector_t v2) { +float innerProduct (const vector_t v1, const vector_t v2) { return v1.c[0]*v2.c[0] + v1.c[1]*v2.c[1] + v1.c[2]*v2.c[2]; } @@ -76,24 +117,24 @@ float distancePointPoint(const vector_t a, return lengthVector(d); } -float distancePointPlane (vector_t point, plane_t plane) { +float distancePointPlane (const vector_t point, const plane_t plane) { return innerProduct(point, plane.normal) - innerProduct(plane.point, plane.normal); } -float distancePointLine (vector_t point, line_t line) { +float distancePointLine (const vector_t point, const line_t line) { return lengthVector(crossProduct(line.direction, subtractVectors(point,line.point))); } -float distancePointLineOnPlane (vector_t point, - line_t line, plane_t plane) { +float distancePointLineOnPlane (vector_t const point, + const line_t line, const 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 addVectors(const vector_t v1, const vector_t v2) { vector_t v3; v3.c[0]=v1.c[0]+v2.c[0]; v3.c[1]=v1.c[1]+v2.c[1]; @@ -101,7 +142,7 @@ vector_t addVectors(vector_t v1, vector_t v2) { return v3; } -vector_t subtractVectors(vector_t v1, vector_t v2) { +vector_t subtractVectors(const vector_t v1, const vector_t v2) { vector_t v3; v3.c[0]=v1.c[0]-v2.c[0]; v3.c[1]=v1.c[1]-v2.c[1]; @@ -109,7 +150,7 @@ vector_t subtractVectors(vector_t v1, vector_t v2) { return v3; } -vector_t scaleVector(vector_t v1, float f) { +vector_t scaleVector(vector_t v1, const float f) { vector_t v2; v2.c[0]=f*v1.c[0]; v2.c[1]=f*v1.c[1]; @@ -227,18 +268,18 @@ line_t lineFromListedPoints(const points_t points, const list_t list) { return l; } -void printVector(vector_t v) { +void printVector(const vector_t v) { printf("[%5.2f,%5.2f,%5.2f], ", v.c[0],v.c[1],v.c[2]); } -void printPlane(plane_t p) { +void printPlane(const plane_t p) { printf("n="); printVector(p.normal); printf(", p="); printVector(p.point); } -void printLine(line_t l) { +void printLine(const line_t l) { printf("d="); printVector(l.direction); printf(", p="); diff --git a/src/mtx_qhull/vectors.h b/src/mtx_qhull/vectors.h index 06a302a..3267df9 100644 --- a/src/mtx_qhull/vectors.h +++ b/src/mtx_qhull/vectors.h @@ -4,6 +4,16 @@ #include #include #include + +/* + * vector operations for zhull + * + * Copyright (c) 2012, Franz Zotter + * IEM, Graz, Austria + * + * + */ + typedef struct vec_ { float c[3]; } vector_t; @@ -24,43 +34,43 @@ typedef struct line_ { } line_t; -vector_t initVector (float x, float y, float z); -float lengthVector(vector_t v); -vector_t normalizeVector(vector_t v); +vector_t initVector (const float x, const float y, const float z); +float lengthVector(const vector_t v); +vector_t normalizeVector(const vector_t v); -plane_t initPlane (vector_t normal, vector_t point); -line_t initLine (vector_t direction, vector_t point); +plane_t initPlane (vector_t normal, const vector_t point); +line_t initLine (vector_t direction, const vector_t point); points_t initPoints (const float *x, const float *y, const float *z, - size_t num_points); + const 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); +vector_t crossProduct (const vector_t v1, const vector_t v2); +float innerProduct (const vector_t v1, const 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); +float distancePointPlane (const vector_t point, const plane_t plane); +float distancePointLine (const vector_t point, const line_t line); +float distancePointLineOnPlane (const vector_t point, + const line_t line, const plane_t plane); +vector_t addVectors(const vector_t v1, const vector_t v2); +vector_t subtractVectors(const vector_t v1, const vector_t v2); +vector_t scaleVector(vector_t v1, const float f); /*vector_t averagePoints(points_t points); vector_t normalOfPoints(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); +plane_t planeFromThreePoints (const vector_t a, const vector_t b, const vector_t c); +line_t lineFromTwoPoints (const vector_t a, const 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); +vector_t averageListedPoints(const points_t points, const list_t list); +vector_t normalOfListedPoints(const points_t points, const list_t list); +vector_t directionOfListedPoints(const points_t points, const list_t list); +plane_t planeFromListedPoints(const points_t points, const list_t list); +line_t lineFromListedPoints(const points_t points, const list_t list); -void printPlane(plane_t p); -void printLine(line_t l); -void printVector(vector_t v); +void printPlane(const plane_t p); +void printLine(const line_t l); +void printVector(const vector_t v); #endif /* QHULL_VECTOR_H */ diff --git a/src/mtx_qhull/zhull.c b/src/mtx_qhull/zhull.c index 6454289..0c4409e 100644 --- a/src/mtx_qhull/zhull.c +++ b/src/mtx_qhull/zhull.c @@ -1,7 +1,30 @@ #include "zhull.h" +/* + * zhull + * + * own qhull algorithm implementation + * + * Copyright (c) 2012, Franz Zotter + * with friendly help from + * IOhannes zmoelnig + * IEM, Graz, Austria + * + * own Implementation after the QHULL algorithm + * that is documented in + * Barber, C.B., Dobkin, D.P., and Huhdanpaa, H.T., + * "The Quickhull algorithm for convex hulls," ACM Trans. + * on Mathematical Software, 22(4):469-483, Dec 1996, + * http://www.qhull.org + * + */ + + + + + /* facets, memory things */ -void freeFacet(facet_t *facet) { +static void freeFacet(facet_t *facet) { /* printf("deleting facet %li\n",facet); printList(facet->corners); printList(facet->outsideset); @@ -13,7 +36,7 @@ void freeFacet(facet_t *facet) { freeList(&(facet->neighbors)); } -list_t appendNewFacets(zhull_t * const zh, const size_t num_facets) { +static list_t appendNewFacets(zhull_t * const zh, const size_t num_facets) { facet_t *new_facet; index_t i; entry_t e0=entry_makeIndex(0); @@ -50,7 +73,7 @@ list_t appendNewFacets(zhull_t * const zh, const size_t num_facets) { } */ -facet_t *getFacetByIndex(const list_t facets, +static facet_t *getFacetByIndex(const list_t facets, const index_t index) { entry_t e=getEntry(facets,index); return ((facet_t*)entry_getPointer(&e)); @@ -66,11 +89,11 @@ index_t getTriangleCorner(const zhull_t * const zh, return 0; } -facet_t *getFacetByPointer(const entry_t e) { +static facet_t *getFacetByPointer(const entry_t e) { return (facet_t*) entry_getPointer(&e); } -void getHorizonEdgeByIndex(index_t *corners, +static void getHorizonEdgeByIndex(index_t *corners, const list_t horizon_fcts, const list_t horizon_fcts_edges, const list_t other_horizon_edges, const index_t index) { index_t i=(index+getLength(horizon_fcts_edges)) @@ -98,18 +121,17 @@ void getHorizonEdgeByIndex(index_t *corners, } -void removeFacetByPointer(zhull_t * const zh, +static void removeFacetByPointer(zhull_t * const zh, facet_t * const pointer) { removeValueFromList(&(zh->facets), entry_makePointer(pointer)); removeValueFromList(&(zh->facets_with_outsidepoints), entry_makePointer(pointer)); removeValueFromList(&(zh->facets_with_insidepoints), entry_makePointer(pointer)); -#warning FIXXXME - freeFacet((facet_t*)pointer); + freeFacet(pointer); } -void removeFacetByPointerList(zhull_t * const zh, +static void removeFacetByPointerList(zhull_t * const zh, const list_t pointers) { index_t i; for (i=0; ifacets, index)); } -void removeFacetByIndexList(zhull_t * const zh, +static void removeFacetByIndexList(zhull_t * const zh, const list_t indices) { facet_t *f; index_t i; @@ -135,7 +157,7 @@ void removeFacetByIndexList(zhull_t * const zh, } } -void freeFacets(zhull_t * const zh) { +static void freeFacets(zhull_t * const zh) { int i; facet_t *f; if (getLength(zh->facets)>0) { @@ -166,7 +188,7 @@ void freeZhull(zhull_t *zh) { -line_t getLine(const zhull_t * const zh, +static line_t getLine(const zhull_t * const zh, const facet_t * const f, const index_t corner) { vector_t a, b; @@ -189,7 +211,7 @@ zhull_t zhullInitPoints(const float *x, const float *y, -void dividePointsBetweenNewFacets ( +static void dividePointsBetweenNewFacets ( zhull_t * const zh, const list_t assoc, const list_t new_facets) { index_t i,j; @@ -236,7 +258,7 @@ void dividePointsBetweenNewFacets ( freeList(&point_inside_facet_list); } -void zhullInitialFacets(zhull_t *zh) { +static void zhullInitialFacets(zhull_t *zh) { list_t assoc = emptyList(); list_t new_facets = emptyList(); index_t i; @@ -269,6 +291,7 @@ void zhullInitialFacets(zhull_t *zh) { } while (idx[2]pts)); if (idx[2]pts)) { getFacetByIndex(new_facets,0)->corners = list; + appendListToList(&(zh->used_pts),list); list=initListIndex(idx,3); reverseList(&list); getFacetByIndex(new_facets,1)->corners = list; @@ -297,7 +320,7 @@ void zhullInitialFacets(zhull_t *zh) { } } -void printHorizonEdges(list_t *horizon_fcts, +static void printHorizonEdges(list_t *horizon_fcts, list_t *horizon_fcts_edges, list_t *other_horizon_edges) { index_t i; @@ -318,7 +341,7 @@ void printHorizonEdges(list_t *horizon_fcts, } -void sortHorizonEdges(list_t *horizon_fcts, +static void sortHorizonEdges(list_t *horizon_fcts, list_t *horizon_fcts_edges, list_t *other_horizon_edges) { index_t i,j; @@ -372,7 +395,7 @@ void sortHorizonEdges(list_t *horizon_fcts, } } -void removeVisibleFacetsGetHorizonAndAvailablePoints( +static void removeVisibleFacetsGetHorizonAndAvailablePoints( zhull_t * const zh, index_t point_index, facet_t *f, list_t *horizon_fcts, @@ -487,7 +510,7 @@ void removeVisibleFacetsGetHorizonAndAvailablePoints( } } -void initNewFacets(zhull_t *zh, index_t point_index, +static void initNewFacets(zhull_t *zh, index_t point_index, list_t new_facets, list_t horizon_fcts, list_t horizon_fcts_edges, list_t other_horizon_edges) { index_t i,j; @@ -572,12 +595,13 @@ void initNewFacets(zhull_t *zh, index_t point_index, } } -void makePyramidFacetsToHorizon(zhull_t *zh, index_t point_index, +static void makePyramidFacetsToHorizon(zhull_t *zh, index_t point_index, list_t horizon_fcts, list_t horizon_fcts_edges, list_t other_horizon_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,other_horizon_edges); + appendToList(&(zh->used_pts),entry_makeIndex(point_index)); /* printf("available points: "); printList(avail_points); printf("new facets : "); @@ -586,12 +610,35 @@ void makePyramidFacetsToHorizon(zhull_t *zh, index_t point_index, freeList(&new_facets); } -void calculateZHull(zhull_t *zh,int maxit) { +static appendExteriorPoints(zhull_t *zh) { + index_t i; + vector_t center = initVector(0.0f,0.0f,0.0f); + list_t facet_delete_list=emptyList(); + facet_t *f; + center=averageListedPoints(zh->pts,zh->used_pts); + printf("central point\n"); + printVector(center); + printf("\n"); + for (i=0; ifacets); i++) { + f=getFacetByIndex(zh->facets,i); + printf("distance of plane %d, d=%5.2f\n",i, + distancePointPlane(center,f->plane)); + if (distancePointPlane(center,f->plane)>-0.5f) { + appendToList(&facet_delete_list,entry_makePointer(f)); + } + } + printList(facet_delete_list); + removeFacetByPointerList(zh,facet_delete_list); + freeList(&facet_delete_list); +} + +int calculateZHull(zhull_t *zh) { index_t fli=0; index_t pi; facet_t *f; list_t outsideset; int cnt=0; + int maxit=getNumPoints(zh->pts); list_t horizon_fcts=emptyList(); list_t horizon_fcts_edges=emptyList(); list_t other_horizon_edges=emptyList(); @@ -599,8 +646,8 @@ void calculateZHull(zhull_t *zh,int maxit) { entry_t e; - if (maxit>MAXIT) - maxit=MAXIT; +// if (maxit>MAXIT) +// maxit=MAXIT; if (getNumPoints(zh->pts)!=0){ zhullInitialFacets(zh); //printZhull(zh); @@ -635,7 +682,9 @@ void calculateZHull(zhull_t *zh,int maxit) { freeList(&available_points); fli++; } +// appendExteriorPoints(zh); } + return cnt; } void printZhull(const zhull_t * const zh) { diff --git a/src/mtx_qhull/zhull.h b/src/mtx_qhull/zhull.h index fc0ff14..cc40f47 100644 --- a/src/mtx_qhull/zhull.h +++ b/src/mtx_qhull/zhull.h @@ -10,7 +10,24 @@ #define TOL_INSIDEPOINT -1e-7 #define TOL_DEGENERATE 1e-6 #define MAXIT 1000000 - +/* + * zhull + * + * own qhull algorithm implementation + * + * Copyright (c) 2012, Franz Zotter + * with friendly help from + * IOhannes zmoelnig + * IEM, Graz, Austria + * + * own Implementation after the QHULL algorithm + * that is documented in + * Barber, C.B., Dobkin, D.P., and Huhdanpaa, H.T., + * "The Quickhull algorithm for convex hulls," ACM Trans. + * on Mathematical Software, 22(4):469-483, Dec 1996, + * http://www.qhull.org + * + */ typedef struct facet_ { plane_t plane; @@ -24,12 +41,13 @@ typedef struct facet_ { typedef struct zhull_ { points_t pts; + list_t used_pts; list_t facets; list_t facets_with_outsidepoints; list_t facets_with_insidepoints; } zhull_t; -void calculateZHull(zhull_t *zh,int maxit); +int calculateZHull(zhull_t *zh); index_t getTriangleCorner(const zhull_t * const zh, const index_t triangle_idx, const index_t corner_idx); -- cgit v1.2.1