3D centroid in PostGIS

Question

How could I find the 3D centroid of this feature in PostGIS:

SELECT
  ST_AsText(
    ST_GeomFromText('POLYGON Z ((0 0 0, 1 0 0, 1 0 1, 1 1 1, 0 1 1, 0 0 1, 0 0 0 ))')
  ) geom;
                        geom                        
---------------------------------------------------------
 POLYGON Z ((0 0 0,1 0 0,1 0 1,1 1 1,0 1 1,0 0 1,0 0 0))
(1 row)

Unfortunately, the following query (because it's based on GEOS) returns the centroid of the 2D projection of the feature in the XY plane:

SELECT
  ST_AsText(
    ST_Centroid(
      ST_GeomFromText('POLYGON Z ((0 0 0, 1 0 0, 1 0 1, 1 1 1, 0 1 1, 0 0 1, 0 0 0 ))')
    )
  ) geom;
   geom    
----------------
 POINT(0.5 0.5)
(1 row)

What I'm searching for is actually the barycenter, aka the center of mass of all vertices of the input shape.
So the expected result should be: POINT Z (0.5 0.333333 0.666667).

A numpy equivalent in Python would simply be:

import numpy as np

A = np.array([[0, 0, 0], [1, 0, 0], [1, 0, 1], [1, 1, 1], [0, 1, 1], [0, 0, 1]])

print(A.mean(axis=0))

>: [0.5        0.33333333 0.66666667]

(There is the ST_GeometricMedian() which can handle 3D correctly, but it's not giving the true centroid according to its documentation.)

Apparently and surprisingly, there is currently no such ST_3DCentroid() like function in SFCGAL:
https://postgis.net/docs/reference.html#reference_sfcgal

So, how could I compute the centroid, aka the center of mass of a 3D feature in PostGIS?

Version info:

PostgreSQL 15.1 (Debian 15.1-1.pgdg110+1) on x86_64-pc-linux-gnu,
  compiled by gcc (Debian 10.2.1-6) 10.2.1 20210110, 64-bit

POSTGIS="3.4.0dev 3.3.0rc2-390-gc2a0b2024"
  [EXTENSION]
    PGSQL="150"
    GEOS="3.12.0dev-CAPI-1.18.0"
    SFCGAL="SFCGAL 1.4.1, CGAL 5.5.1, BOOST 1.74.0"
    PROJ="9.2.0"
    LIBXML="2.9.10"
    LIBJSON="0.15"
    LIBPROTOBUF="1.3.3"
    WAGYU="0.5.0 (Internal)"
    TOPOLOGY

Answer 1

I'm also posting here what I have found so far as an answer to keep things clearly separated.

So this is my current workaround, maybe other (future) people could have better/faster/more robust solutions:

WITH cte AS
(SELECT 
  geom(
    ST_DumpPoints(
      ST_RemovePoint(
        ST_ExteriorRing(
          ST_GeomFromText('POLYGON Z ((0 0 0, 1 0 0, 1 0 1, 1 1 1, 0 1 1, 0 0 1, 0 0 0 ))')
        ), 0 -- remove first point because, for a polygon, first point = last point by definition
      )
    )
  )
)
SELECT 
  ST_AsText(
    ST_MakePoint(
      AVG(ST_X(geom)),
      AVG(ST_Y(geom)),
      AVG(ST_Z(geom))
    )
  ) geom
FROM cte;
                       geom                        
---------------------------------------------------
 POINT Z (0.5 0.333333333333333 0.666666666666667)