I added your suggestion to get this (finally):
def get_centroid(pts)
total_area = 0
total_centroids = Geom::Vector3d.new(0,0,0)
third = Geom::Transformation.scaling(1.0 / 3.0)
npts = pts.length
vec1 = Geom::Vector3d.new(pts[1].x - pts[0].x, pts[1].y - pts[0].y, pts[1].z - pts[0].z)
vec2 = Geom::Vector3d.new(pts[2].x - pts[0].x, pts[2].y - pts[0].y, pts[2].z - pts[0].z)
ref_sense = vec1.cross vec2
for i in 0...(npts-2)
vec1 = Geom::Vector3d.new(pts[i+1].x - pts[0].x, pts[i+1].y - pts[0].y, pts[i+1].z - pts[0].z)
vec2 = Geom::Vector3d.new(pts[i+2].x - pts[0].x, pts[i+2].y - pts[0].y, pts[i+2].z - pts[0].z)
vec = vec1.cross vec2
area = vec.length / 2.0
if(ref_sense.dot(vec) < 0)
area *= -1.0
end
total_area += area
centroid = (vec1 + vec2).transform(third)
t = Geom::Transformation.scaling(area)
total_centroids += centroid.transform(t)
end
c = Geom::Transformation.scaling(1.0 / total_area)
total_centroids.transform!(c) + Geom::Vector3d.new(pts[0].x,pts[0].y,pts[0].z)
end
This should work for convex and concave polygons.