Fix thread-safety bug in geom.c

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,27 +2223,40 @@ 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);
     if (b < slist[ns]) break;
     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);
 }

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))))
 
 ; ****************************************************************
-