Index: include/bn.h
===================================================================
--- include/bn.h	(Revision 59252)
+++ include/bn.h	(Arbeitskopie)
@@ -4630,11 +4630,29 @@
  * The array mustn't contain duplicated points.
  *
  * @return 0 if calculation was successful
- * @return 1 if calculation failed, e.g. because one in-parameter is a NULL-pointer
+ * @return 1 if calculation failed, e.g. because one parameter is a NULL-pointer
  */
 BN_EXPORT extern int bn_polygon_centroid(point_t *cent, size_t npts, const point_t *pts);
 
 
+/**
+ * calculate for an array of plane_eqs, which build a polyhedron, the
+ * point_t's for each face.
+ *
+ * @param[out] npts Array, which stores for every face the number of
+ * point_ts, added to pts. Needs to be allocated with npts[neqs] already.
+ * @param[out] pts 2D-array which stores the point_ts for every
+ * face. The array needs to be allocated with pts[neqs][neqs-1] already.
+ * @param[in] neqs Number of plane_ts, stored in eqs
+ * @param[in] eqs Array, that contains the plane equations, which
+ * build the polyhedron
+ *
+ * @return 0 if calculation was successful
+ * @return 1 if calculation failed, e.g. because one parameter is a NULL-Pointer
+ */
+BN_EXPORT extern int bn_polygon_mk_pts_planes(size_t *npts, point_t **pts, size_t neqs, const plane_t *eqs);
+
+
 #endif /* BN_H */
 
 /** @} */
Index: src/libbn/polygon.c
===================================================================
--- src/libbn/polygon.c	(Revision 59252)
+++ src/libbn/polygon.c	(Arbeitskopie)
@@ -131,6 +131,42 @@
 }
 
 
+int
+bn_polygon_mk_pts_planes(size_t *npts, point_t **pts, size_t neqs, const plane_t *eqs)
+{
+    size_t i, j, k, l;
+    if (!npts || !pts || neqs < 4 || !eqs)
+    	return 1;
+    /* find all vertices */
+    for (i = 0; i < neqs - 2; i++) {
+	for (j = i + 1; j < neqs - 1; j++) {
+	    for (k = j + 1; k < neqs; k++) {
+		point_t pt;
+		int keep_point = 1;
+		if (bn_mkpoint_3planes(pt, eqs[i], eqs[j], eqs[k]) < 0)
+		    continue;
+		/* discard pt if it is outside the polyhedron */
+		for (l = 0; l < neqs; l++) {
+		    if (l == i || l == j || l == k)
+			continue;
+		    if (DIST_PT_PLANE(pt, eqs[l]) > BN_TOL_DIST) {
+			keep_point = 0;
+			break;
+		    }
+		}
+		/* found a good point, add it to each of the intersecting faces */
+		if (keep_point) {
+		    VMOVE(pts[i][npts[i]], (pt)); npts[i]++;
+		    VMOVE(pts[j][npts[j]], (pt)); npts[j]++;
+		    VMOVE(pts[k][npts[k]], (pt)); npts[k]++;
+		}
+	    }
+	}
+   }
+    return 0;
+}
+
+
 /*
  * Local Variables:
  * mode: C
Index: src/libged/analyze.c
===================================================================
--- src/libged/analyze.c	(Revision 59252)
+++ src/libged/analyze.c	(Arbeitskopie)
@@ -923,12 +923,18 @@
 HIDDEN void
 analyze_arbn(struct ged *gedp, const struct rt_db_internal *ip)
 {
-    size_t i, j, k, l;
+    size_t i;
     fastf_t tot_vol = 0.0, tot_area = 0.0;
     table_t table;
     struct poly_face *faces;
     struct bu_vls tmpstr = BU_VLS_INIT_ZERO;
     struct rt_arbn_internal *aip = (struct rt_arbn_internal *)ip->idb_ptr;
+    size_t *npts;
+    point_t **tmp_pts;
+    plane_t *eqs;
+    tmp_pts = (point_t **)bu_calloc(aip->neqn, sizeof(point_t *), "analyze_arbn: tmp_pts");
+    npts = (size_t *)bu_calloc(aip->neqn, sizeof(size_t), "analyze_arbn: npts");
+    eqs = (plane_t *)bu_calloc(aip->neqn, sizeof(plane_t), "analyze_arbn: eqs");
 
     /* allocate array of face structs */
     faces = (struct poly_face *)bu_calloc(aip->neqn, sizeof(struct poly_face), "analyze_arbn: faces");
@@ -937,45 +943,22 @@
 	VUNITIZE(faces[i].plane_eqn);
 	/* allocate array of pt structs, max number of verts per faces = (# of faces) - 1 */
 	faces[i].pts = (point_t *)bu_calloc(aip->neqn - 1, sizeof(point_t), "analyze_arbn: pts");
+	tmp_pts[i] = faces[i].pts;
+    	HMOVE(eqs[i], faces[i].plane_eqn);
     }
     /* allocate table rows, 1 row per plane eqn */
     table.rows = (row_t *)bu_calloc(aip->neqn, sizeof(row_t), "analyze_arbn: rows");
     table.nrows = aip->neqn;
 
-    /* find all vertices */
-    for (i = 0; i < aip->neqn - 2; i++) {
-	for (j = i + 1; j < aip->neqn - 1; j++) {
-	    for (k = j + 1; k < aip->neqn; k++) {
-		point_t pt;
-		int keep_point = 1;
-		if (bn_mkpoint_3planes(pt, aip->eqn[i], aip->eqn[j], aip->eqn[k]) < 0) {
-		    continue;
-		}
-		/* discard pt if it is outside the arbn */
-		for (l = 0; l < aip->neqn; l++) {
-		    if (l == i || l == j || l == k) {
-			continue;
-		    }
-		    if (DIST_PT_PLANE(pt, aip->eqn[l]) > gedp->ged_wdbp->wdb_tol.dist) {
-			keep_point = 0;
-			break;
-		    }
-		}
-		/* found a good point, add it to each of the intersecting faces */
-		if (keep_point) {
-		    ADD_PT(faces[i], pt);
-		    ADD_PT(faces[j], pt);
-		    ADD_PT(faces[k], pt);
-		}
-	    }
-	}
-    }
+    bn_polygon_mk_pts_planes(npts, tmp_pts, aip->neqn, (const plane_t *)eqs);
 
     for (i = 0; i < aip->neqn; i++) {
 	vect_t tmp;
 	bu_vls_sprintf(&tmpstr, "%4zu", i);
 	snprintf(faces[i].label, sizeof(faces[i].label), "%s", bu_vls_addr(&tmpstr));
 
+	faces[i].npts = npts[i];
+
 	/* calculate surface area */
 	analyze_poly_face(gedp, &faces[i], &table.rows[i]);
 	tot_area += faces[i].area;
@@ -1003,6 +986,9 @@
     }
     bu_free((char *)faces, "analyze_arbn: faces");
     bu_free((char *)table.rows, "analyze_arbn: rows");
+    bu_free((char *)tmp_pts, "analyze_arbn: tmp_pts");
+    bu_free((char *)npts, "analyze_arbn: npts");
+    bu_free((char *)eqs, "analyze_arbn: eqs");
     bu_vls_free(&tmpstr);
 }
 
Index: src/librt/primitives/arbn/arbn.c
===================================================================
--- src/librt/primitives/arbn/arbn.c	(Revision 59252)
+++ src/librt/primitives/arbn/arbn.c	(Arbeitskopie)
@@ -1348,41 +1348,29 @@
 static void
 rt_arbn_faces_area(struct poly_face* faces, struct rt_arbn_internal* aip)
 {
-    size_t i, j, l ,k;
+    size_t i;
+    size_t *npts;
+    point_t **tmp_pts;
+    plane_t *eqs;
+    tmp_pts = (point_t **)bu_calloc(aip->neqn, sizeof(point_t *), "rt_arbn_faces_area: tmp_pts");
+    npts = (size_t *)bu_calloc(aip->neqn, sizeof(size_t), "rt_arbn_faces_area: npts");
+    eqs = (plane_t *)bu_calloc(aip->neqn, sizeof(plane_t), "rt_arbn_faces_area: eqs");
+
     for (i = 0; i < aip->neqn; i++) {
     	HMOVE(faces[i].plane_eqn, aip->eqn[i]);
     	VUNITIZE(faces[i].plane_eqn);
+	tmp_pts[i] = faces[i].pts;
+    	HMOVE(eqs[i], faces[i].plane_eqn);
     }
-    /* find all vertices */
-    for (i = 0; i < aip->neqn - 2; i++) {
-	for (j = i + 1; j < aip->neqn - 1; j++) {
-	    for (k = j + 1; k < aip->neqn; k++) {
-		point_t pt;
-		int keep_point = 1;
-		if (bn_mkpoint_3planes(pt, aip->eqn[i], aip->eqn[j], aip->eqn[k]) < 0)
-		    continue;
-		/* discard pt if it is outside the arbn */
-		for (l = 0; l < aip->neqn; l++) {
-		    if (l == i || l == j || l == k)
-			continue;
-		    if (DIST_PT_PLANE(pt, aip->eqn[l]) > BN_TOL_DIST) {
-			keep_point = 0;
-			break;
-		    }
-		}
-		/* found a good point, add it to each of the intersecting faces */
-		if (keep_point) {
-		    VMOVE((faces[i]).pts[(faces[i]).npts], (pt)); (faces[i]).npts++;
-		    VMOVE((faces[j]).pts[(faces[j]).npts], (pt)); (faces[j]).npts++;
-		    VMOVE((faces[k]).pts[(faces[k]).npts], (pt)); (faces[k]).npts++;
-		}
-	    }
-	}
-    }
+    bn_polygon_mk_pts_planes(npts, tmp_pts, aip->neqn, (const plane_t *)eqs);
     for (i = 0; i < aip->neqn; i++) {
+	faces[i].npts = npts[i];
     	bu_sort(faces[i].pts, faces[i].npts, sizeof(point_t), ccw, &faces[i].plane_eqn);
     	bn_polygon_area(&faces[i].area, faces[i].npts, (const point_t *)faces[i].pts);
     }
+    bu_free((char *)tmp_pts, "rt_arbn_faces_area: tmp_pts");
+    bu_free((char *)npts, "rt_arbn_faces_area: npts");
+    bu_free((char *)eqs, "rt_arbn_faces_area: eqs");
 }
 
 
Index: src/librt/primitives/nmg/nmg.c
===================================================================
--- src/librt/primitives/nmg/nmg.c	(Revision 59252)
+++ src/librt/primitives/nmg/nmg.c	(Arbeitskopie)
@@ -3136,45 +3136,33 @@
 rt_nmg_faces_area(struct poly_face* faces, struct shell* s)
 {
     struct bu_ptbl nmg_faces;
-    unsigned int num_faces, i, j, k, l;
+    unsigned int num_faces, i;
+    size_t *npts;
+    point_t **tmp_pts;
+    plane_t *eqs;
     nmg_face_tabulate(&nmg_faces, &s->l.magic);
     num_faces = BU_PTBL_LEN(&nmg_faces);
+    tmp_pts = (point_t **)bu_calloc(num_faces, sizeof(point_t *), "rt_nmg_faces_area: tmp_pts");
+    npts = (size_t *)bu_calloc(num_faces, sizeof(size_t), "rt_nmg_faces_area: npts");
+    eqs = (plane_t *)bu_calloc(num_faces, sizeof(plane_t), "rt_nmg_faces_area: eqs");
+
     for(i = 0; i < num_faces; i++) {
 	struct face *f;
 	f = (struct face *)BU_PTBL_GET(&nmg_faces, i);
 	HMOVE(faces[i].plane_eqn, f->g.plane_p->N);
 	VUNITIZE(faces[i].plane_eqn);
+	tmp_pts[i] = faces[i].pts;
+    	HMOVE(eqs[i], faces[i].plane_eqn);
     }
-    /* find all vertices */
-    for (i = 0; i < num_faces - 2; i++) {
-	for (j = i + 1; j < num_faces - 1; j++) {
-	    for (k = j + 1; k < num_faces; k++) {
-		point_t pt;
-		int keep_point = 1;
-		if (bn_mkpoint_3planes(pt, faces[i].plane_eqn, faces[j].plane_eqn, faces[k].plane_eqn) < 0)
-		    continue;
-		/* discard pt if it is outside the arbn */
-		for (l = 0; l < num_faces; l++) {
-		    if (l == i || l == j || l == k)
-			continue;
-		    if (DIST_PT_PLANE(pt, faces[l].plane_eqn) > BN_TOL_DIST) {
-			keep_point = 0;
-			break;
-		    }
-		}
-		/* found a good point, add it to each of the intersecting faces */
-		if (keep_point) {
-		    VMOVE((faces[i]).pts[(faces[i]).npts], (pt)); (faces[i]).npts++;
-		    VMOVE((faces[j]).pts[(faces[j]).npts], (pt)); (faces[j]).npts++;
-		    VMOVE((faces[k]).pts[(faces[k]).npts], (pt)); (faces[k]).npts++;
-		}
-	    }
-	}
-    }
+    bn_polygon_mk_pts_planes(npts, tmp_pts, num_faces, (const plane_t *)eqs);
     for (i = 0; i < num_faces; i++) {
+	faces[i].npts = npts[i];
 	bu_sort(faces[i].pts, faces[i].npts, sizeof(point_t), nmg_ccw, &faces[i].plane_eqn);
 	bn_polygon_area(&faces[i].area, faces[i].npts, (const point_t *)faces[i].pts);
     }
+    bu_free((char *)tmp_pts, "rt_nmg_faces_area: tmp_pts");
+    bu_free((char *)npts, "rt_nmg_faces_area: npts");
+    bu_free((char *)eqs, "rt_nmg_faces_area: eqs");
 }
 
 
