Index: include/bn.h
===================================================================
--- include/bn.h	(Revision 59246)
+++ include/bn.h	(Arbeitskopie)
@@ -4621,6 +4621,20 @@
 BN_EXPORT extern int bn_polygon_area(fastf_t *area, size_t npts, const point_t *pts);
 
 
+/**
+ * calculate the centroid of a non self-intersecting polygon
+ *
+ * @param[out] cent The centroid of the polygon
+ * @param[in] npts Number of point_ts, stored in pts
+ * @param[in] pts all points of the polygon, sorted counter-clockwise.
+ * 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
+ */
+BN_EXPORT extern int bn_polygon_centroid(point_t *cent, size_t npts, const point_t *pts);
+
+
 #endif /* BN_H */
 
 /** @} */
Index: src/libbn/polygon.c
===================================================================
--- src/libbn/polygon.c	(Revision 59246)
+++ src/libbn/polygon.c	(Arbeitskopie)
@@ -64,6 +64,73 @@
     *area = fabs(VDOT(plane_eqn, tot)) * 0.5;
     return 0;
 }
+
+
+int
+bn_polygon_centroid(point_t *cent, size_t npts, const point_t *pts)
+{
+    size_t i;
+    fastf_t x_0 = 0.0;
+    fastf_t x_1 = 0.0;
+    fastf_t y_0 = 0.0;
+    fastf_t y_1 = 0.0;
+    fastf_t z_0 = 0.0;
+    fastf_t z_1 = 0.0;
+    fastf_t a = 0.0;
+    fastf_t signedArea = 0.0;
+
+    if (!pts || !cent || npts < 3)
+	return 1;
+    /* Calculate Centroid projection for face for x-y-plane */
+    for (i = 0; i < npts-1; i++) {
+	x_0 = pts[i][0];
+	y_0 = pts[i][1];
+	x_1 = pts[i+1][0];
+	y_1 = pts[i+1][1];
+	a = x_0 *y_1 - x_1*y_0;
+	signedArea += a;
+	*cent[0] += (x_0 + x_1)*a;
+	*cent[1] += (y_0 + y_1)*a;
+    }
+    x_0 = pts[i][0];
+    y_0 = pts[i][1];
+    x_1 = pts[0][0];
+    y_1 = pts[0][1];
+    a = x_0 *y_1 - x_1*y_0;
+    signedArea += a;
+    *cent[0] += (x_0 + x_1)*a;
+    *cent[1] += (y_0 + y_1)*a;
+
+    signedArea *= 0.5;
+    *cent[0] /= (6.0*signedArea);
+    *cent[1] /= (6.0*signedArea);
+
+    /* calculate Centroid projection for face for x-z-plane */
+
+    signedArea = 0.0;
+    for (i = 0; i < npts-1; i++) {
+	x_0 = pts[i][0];
+	z_0 = pts[i][2];
+	x_1 = pts[i+1][0];
+	z_1 = pts[i+1][2];
+	a = x_0 *z_1 - x_1*z_0;
+	signedArea += a;
+	*cent[2] += (z_0 + z_1)*a;
+    }
+    x_0 = pts[i][0];
+    z_0 = pts[i][2];
+    x_1 = pts[0][0];
+    z_0 = pts[0][2];
+    a = x_0 *z_1 - x_1*z_0;
+    signedArea += a;
+    *cent[2] += (z_0 + z_1)*a;
+
+    signedArea *= 0.5;
+    *cent[2] /= (6.0*signedArea);
+    return 0;
+}
+
+
 /*
  * Local Variables:
  * mode: C
Index: src/librt/primitives/arbn/arbn.c
===================================================================
--- src/librt/primitives/arbn/arbn.c	(Revision 59246)
+++ src/librt/primitives/arbn/arbn.c	(Arbeitskopie)
@@ -1413,7 +1413,7 @@
 {
     struct poly_face *faces;
     struct rt_arbn_internal *aip = (struct rt_arbn_internal *)ip->idb_ptr;
-    size_t i, j;
+    size_t i;
     point_t arbit_point = VINIT_ZERO;
     fastf_t volume = 0.0;
 
@@ -1432,60 +1432,7 @@
     }
     rt_arbn_faces_area(faces, aip);
     for (i = 0; i < aip->neqn; i++) {
-	fastf_t x_0 = 0.0;
-	fastf_t x_1 = 0.0;
-        fastf_t y_0 = 0.0;
-	fastf_t y_1 = 0.0;
-	fastf_t z_0 = 0.0;
-	fastf_t z_1 = 0.0;
-	fastf_t a = 0.0;
-	fastf_t signedArea = 0.0;
-	/* Calculate Centroid projection for face for x-y-plane */
-	for (j = 0; j < faces[i].npts-1; j++) {
-	    x_0 = faces[i].pts[j][0];
-	    y_0 = faces[i].pts[j][1];
-	    x_1 = faces[i].pts[j+1][0];
-	    y_1 = faces[i].pts[j+1][1];
-	    a = x_0 *y_1 - x_1*y_0;
-	    signedArea += a;
-	    faces[i].cent[0] += (x_0 + x_1)*a;
-	    faces[i].cent[1] += (y_0 + y_1)*a;
-	}
-	x_0 = faces[i].pts[j][0];
-	y_0 = faces[i].pts[j][1];
-	x_1 = faces[i].pts[0][0];
-	y_1 = faces[i].pts[0][1];
-	a = x_0 *y_1 - x_1*y_0;
-	signedArea += a;
-	faces[i].cent[0] += (x_0 + x_1)*a;
-	faces[i].cent[1] += (y_0 + y_1)*a;
-
-	signedArea *= 0.5;
-	faces[i].cent[0] /= (6.0*signedArea);
-	faces[i].cent[1] /= (6.0*signedArea);
-
-	/* calculate Centroid projection for face for x-z-plane */
-
-	signedArea = 0.0;
-	for (j = 0; j < faces[i].npts-1; j++) {
-	    x_0 = faces[i].pts[j][0];
-	    z_0 = faces[i].pts[j][2];
-	    x_1 = faces[i].pts[j+1][0];
-	    z_1 = faces[i].pts[j+1][2];
-	    a = x_0 *z_1 - x_1*z_0;
-	    signedArea += a;
-	    faces[i].cent[2] += (z_0 + z_1)*a;
-	}
-	x_0 = faces[i].pts[j][0];
-	z_0 = faces[i].pts[j][2];
-	x_1 = faces[i].pts[0][0];
-	z_0 = faces[i].pts[0][2];
-	a = x_0 *z_1 - x_1*z_0;
-	signedArea += a;
-	faces[i].cent[2] += (z_0 + z_1)*a;
-
-	signedArea *= 0.5;
-	faces[i].cent[2] /= (6.0*signedArea);
+	bn_polygon_centroid(&faces[i].cent, faces[i].npts, (const point_t *) faces[i].pts);
 	VADD2(arbit_point, arbit_point, faces[i].cent);
 
     }
Index: src/librt/primitives/nmg/nmg.c
===================================================================
--- src/librt/primitives/nmg/nmg.c	(Revision 59246)
+++ src/librt/primitives/nmg/nmg.c	(Arbeitskopie)
@@ -3226,7 +3226,7 @@
     struct bu_ptbl nmg_faces;
     fastf_t volume = 0.0;
     point_t arbit_point = VINIT_ZERO;
-    unsigned int num_faces, i, j;
+    size_t num_faces, i;
 
     *cent[0] = 0.0;
     *cent[1] = 0.0;
@@ -3246,61 +3246,7 @@
     }
     rt_nmg_faces_area(faces, s);
     for (i = 0; i < num_faces; i++) {
-	fastf_t x_0 = 0.0;
-	fastf_t x_1 = 0.0;
-	fastf_t y_0 = 0.0;
-	fastf_t y_1 = 0.0;
-	fastf_t z_0 = 0.0;
-	fastf_t z_1 = 0.0;
-	fastf_t a = 0.0;
-	fastf_t signedArea = 0.0;
-
-	/* Calculate Centroid projection for face for x-y-plane */
-	for (j = 0; j < faces[i].npts-1; j++) {
-	    x_0 = faces[i].pts[j][0];
-	    y_0 = faces[i].pts[j][1];
-	    x_1 = faces[i].pts[j+1][0];
-	    y_1 = faces[i].pts[j+1][1];
-	    a = x_0 *y_1 - x_1*y_0;
-	    signedArea += a;
-	    faces[i].cent[0] += (x_0 + x_1)*a;
-	    faces[i].cent[1] += (y_0 + y_1)*a;
-	}
-	x_0 = faces[i].pts[j][0];
-	y_0 = faces[i].pts[j][1];
-	x_1 = faces[i].pts[0][0];
-	y_1 = faces[i].pts[0][1];
-	a = x_0 *y_1 - x_1*y_0;
-	signedArea += a;
-	faces[i].cent[0] += (x_0 + x_1)*a;
-	faces[i].cent[1] += (y_0 + y_1)*a;
-
-	signedArea *= 0.5;
-	faces[i].cent[0] /= (6.0*signedArea);
-	faces[i].cent[1] /= (6.0*signedArea);
-
-	/* calculate Centroid projection for face for x-z-plane */
-
-	signedArea = 0.0;
-	for (j = 0; j < faces[i].npts-1; j++) {
-	    x_0 = faces[i].pts[j][0];
-	    z_0 = faces[i].pts[j][2];
-	    x_1 = faces[i].pts[j+1][0];
-	    z_1 = faces[i].pts[j+1][2];
-	    a = x_0 *z_1 - x_1*z_0;
-	    signedArea += a;
-	    faces[i].cent[2] += (z_0 + z_1)*a;
-	}
-	x_0 = faces[i].pts[j][0];
-	z_0 = faces[i].pts[j][2];
-	x_1 = faces[i].pts[0][0];
-	z_0 = faces[i].pts[0][2];
-	a = x_0 *z_1 - x_1*z_0;
-	signedArea += a;
-	faces[i].cent[2] += (z_0 + z_1)*a;
-
-	signedArea *= 0.5;
-	faces[i].cent[2] /= (6.0*signedArea);
+	bn_polygon_centroid(&faces[i].cent, faces[i].npts, (const point_t *) faces[i].pts);
 	VADD2(arbit_point, arbit_point, faces[i].cent);
 
     }
