diff --git a/utils/geom.c b/utils/geom.c index dddf0c3..4b847d2 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -1007,8 +1007,6 @@ static number compute_dot_cross(vector3 a, vector3 b, vector3 c) { Requires that geometry_lattice global has been initialized, etcetera. */ void geom_get_bounding_box(geometric_object o, geom_box *box) { - geom_fix_object_ptr(&o); - /* initialize to empty box at the center of the object: */ box->low = box->high = o.center; @@ -2148,7 +2146,7 @@ static int dcmp(const void *pd1, const void *pd2) { /* the intersection s-values sorted in ascending order. */ /* the return value is the number of intersections. */ /******************************************************************/ -int intersect_line_with_prism(prism *prsm, vector3 pc, vector3 dc, double *slist) { +int intersect_line_with_prism(prism *prsm, vector3 pc, vector3 dc, double *slist, int slist_len) { vector3 pp = prism_coordinate_c2p(prsm, pc); vector3 dp = prism_vector_c2p(prsm, dc); vector3 *vps_bottom = prsm->vertices_p.items; @@ -2183,7 +2181,7 @@ int intersect_line_with_prism(prism *prsm, vector3 pc, vector3 dc, double *slist vector3 tus = matrix3x3_vector3_mult(matrix3x3_inverse(M), RHS); if (tus.x < -tus_tolerance || tus.x > 1+tus_tolerance || tus.y < -tus_tolerance || tus.y > 1+tus_tolerance) continue; double s = tus.z; - slist[num_intersections++] = s; + if (num_intersections++ < slist_len) slist[num_intersections-1] = s; } // identify intersections with prism ceiling and floor faces @@ -2194,9 +2192,12 @@ int intersect_line_with_prism(prism *prsm, vector3 pc, vector3 dc, double *slist double s = (z0p - pp.z) / dp.z; vector3 *vps = LowerUpper ? vps_top : vps_bottom; if (!node_in_polygon(pp.x + s * dp.x, pp.y + s * dp.y, vps, num_vertices)) continue; - slist[num_intersections++] = s; + if (num_intersections++ < slist_len) slist[num_intersections-1] = s; } + // caller needs a bigger slist array + if (num_intersections > slist_len) return num_intersections; + qsort((void *)slist, num_intersections, sizeof(double), dcmp); // if num_intersections is zero then just return that if (num_intersections == 0) return num_intersections; @@ -2222,20 +2223,30 @@ int intersect_line_with_prism(prism *prsm, vector3 pc, vector3 dc, double *slist /***************************************************************/ /***************************************************************/ double intersect_line_segment_with_prism(prism *prsm, vector3 pc, vector3 dc, double a, double b) { - double *slist = prsm->workspace.items; - int num_intersections = intersect_line_with_prism(prsm, pc, dc, slist); + int num_vertices = prsm->vertices_p.num_items; + double *slist; + int num_intersections, num_intersections_alloc; + + double slist_stack[1024]; // usually, a small stack-allocated array will be enough + slist = slist_stack; + num_intersections = num_intersections_alloc = intersect_line_with_prism(prsm, pc, dc, slist, 1024); + if (num_intersections_alloc > 1024) { // slow fallback + slist = (double *) malloc(num_intersections_alloc * sizeof(double)); + num_intersections = intersect_line_with_prism(prsm, pc, dc, slist, num_intersections_alloc); + CHECK(num_intersections <= num_intersections_alloc, "bug in intersect_line_with_prism"); + } // na=smallest index such that slist[na] > a int na = -1; int ns; + double ds = 0.0; for (ns = 0; na == -1 && ns < num_intersections; ns++) if (slist[ns] > a) na = ns; - if (na == -1) return 0.0; + if (na == -1) goto done; int inside = ((na % 2) == 0 ? 0 : 1); double last_s = a; - double ds = 0.0; for (ns = na; ns < num_intersections; ns++) { double this_s = fmin(b, slist[ns]); if (inside) ds += (this_s - last_s); @@ -2243,6 +2254,9 @@ double intersect_line_segment_with_prism(prism *prsm, vector3 pc, vector3 dc, do inside = (1 - inside); last_s = this_s; } + +done: + if (num_intersections_alloc > 1024) free(slist); return ds > 0.0 ? ds : 0.0; } @@ -2822,11 +2836,6 @@ void init_prism(geometric_object *o) { for (nv = 0; nv < num_vertices; nv++) { prsm->vertices_top.items[nv] = prism_coordinate_p2c(prsm, prsm->vertices_top_p.items[nv]); } - - // workspace is an internally-stored double-valued array of length num_vertices+2 - // that is used by some geometry routines - prsm->workspace.num_items = num_vertices + 2; - prsm->workspace.items = (double *)malloc((num_vertices + 2) * sizeof(double)); } /* like init_prism, but works with an already-initialied prism */ @@ -2838,7 +2847,6 @@ void reinit_prism(geometric_object *o) { free(prsm->top_polygon_diff_vectors_p.items); free(prsm->top_polygon_diff_vectors_scaled_p.items); free(prsm->vertices_top.items); - free(prsm->workspace.items); init_prism(o); } diff --git a/utils/geom.scm b/utils/geom.scm index 6b03a94..e2a9284 100644 --- a/utils/geom.scm +++ b/utils/geom.scm @@ -135,7 +135,6 @@ (define-property vertices_top_p '() (make-list-type 'vector3)) (define-property vertices_top '() (make-list-type 'vector3)) (define-property centroid (vector3 0 0 0) 'vector3) - (define-property workspace '() (make-list-type 'number)) (define-property m_c2p identity_matrix 'matrix3x3) (define-property m_p2c identity_matrix 'matrix3x3)) @@ -355,4 +354,3 @@ (reciprocal->cartesian v)))) ; **************************************************************** -