Skip to content

Instantly share code, notes, and snippets.

@geozelot

geozelot/STCG_ApproximateMedialSpine.sql

Last active November 6, 2025 15:25
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save geozelot/096627a91cf36004f0c42510eee5d108 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save geozelot/096627a91cf36004f0c42510eee5d108 to your computer and use it in GitHub Desktop.
Download ZIP
PostgreSQL/PostGIS - Approximate the spine of the medial axis of an aeral geometry assuming maximum span across the skeletal graph
Raw
STCG_ApproximateMedialSpine.sql
/*
* Wraps ST_ApproximateMedialAxis (https://postgis.net/docs/ST_ApproximateMedialAxis.html) and attempts
* to extract the skeletal spine by assuming maximum span across the axis tree.
*
* Note:
* The current version of this function extracts the shell of each (component of each MULTI-)POLYGON
* to avoid handling cyclic graphs around inner rings, i.e. the resulting spines may cross any polygon hole.
*
*
* Two versions, targeting either the ST_ prefix or CG_ prefix functions.
*/
CREATE OR REPLACE FUNCTION public.ST_ApproximateMedialSpine(
geom GEOMETRY
)
RETURNS GEOMETRY
LANGUAGE SQL
IMMUTABLE STRICT
PARALLEL SAFE
COST 100
AS
$$
WITH RECURSIVE
reduce AS (
SELECT
component.path[1] AS seq,
ST_NumGeometries(noded_axis) AS n,
noded_axis AS geom
FROM
ST_Dump($1) AS component,
LATERAL ST_MakePolygon(ST_ExteriorRing(component.geom)) AS shell,
LATERAL ST_ApproximateMedialAxis(shell) AS medial_axis,
LATERAL ST_LineMerge(medial_axis) AS noded_axis
UNION ALL
SELECT
merge.seq,
ST_NumGeometries(merge.geom) AS n,
merge.geom
FROM
(
SELECT
eval.seq,
ST_LineMerge(ST_Collect(eval.edge_geom)) AS geom
FROM
(
SELECT
card.seq,
card.edge_id, card.edge_geom,
CASE WHEN card.is_dangle
THEN ROW_NUMBER() OVER(PARTITION BY card.seq, card.node_id, card.is_dangle ORDER BY card.edge_length DESC) - card.degree = 0
ELSE ROW_NUMBER() OVER(PARTITION BY card.seq, card.edge_id ORDER BY card.node_id) - card.degree = 0
END AS is_major
FROM
(
SELECT
conn.seq,
conn.edge_id,
conn.node_id,
conn.degree, MIN(conn.degree) OVER(PARTITION BY conn.seq, conn.edge_id) = 0 AS is_dangle,
conn.edge_length, conn.edge_geom
FROM
(
SELECT
reduce.seq,
edge.path[1] AS edge_id,
DENSE_RANK() OVER(PARTITION BY reduce.seq ORDER BY node.geom) AS node_id,
(COUNT(edge.path[1]) OVER(PARTITION BY reduce.seq, node.geom) > 1)::INT AS degree,
ST_Length(edge.geom::GEOGRAPHY) AS edge_length,
edge.geom AS edge_geom
FROM
reduce,
ST_Dump(reduce.geom) AS edge,
LATERAL ST_Boundary(edge.geom) AS nodes,
LATERAL ST_Dump(nodes) AS node
) AS conn
) AS card
) AS eval
WHERE
eval.is_major
GROUP BY
eval.seq
) AS merge
WHERE
ST_NumGeometries(merge.geom) > 0
)
SELECT
ST_LineMerge(ST_Collect(final.geom)) AS geom
FROM
(
SELECT
DENSE_RANK() OVER(PARTITION BY reduce.seq ORDER BY reduce.n DESC) = 1 AS keep_n,
ROW_NUMBER() OVER(PARTITION BY reduce.seq, reduce.n ORDER BY ST_Length(dmp.geom::GEOGRAPHY) DESC) <= 2 AS keep_len,
reduce.seq,
dmp.geom
FROM
reduce,
LATERAL ST_Dump(reduce.geom) AS dmp
WHERE
reduce.n <= 3
) AS final
WHERE
final.keep_n AND final.keep_len
;
$$
;
CREATE OR REPLACE FUNCTION public.CG_ApproximateMedialSpine(
geom GEOMETRY
)
RETURNS GEOMETRY
LANGUAGE SQL
IMMUTABLE STRICT
PARALLEL SAFE
COST 100
AS
$$
WITH RECURSIVE
reduce AS (
SELECT
component.path[1] AS seq,
ST_NumGeometries(noded_axis) AS n,
noded_axis AS geom
FROM
ST_Dump($1) AS component,
LATERAL ST_MakePolygon(ST_ExteriorRing(component.geom)) AS shell,
LATERAL CG_ApproximateMedialAxis(shell) AS medial_axis,
LATERAL ST_LineMerge(medial_axis) AS noded_axis
UNION ALL
SELECT
merge.seq,
ST_NumGeometries(merge.geom) AS n,
merge.geom
FROM
(
SELECT
eval.seq,
ST_LineMerge(ST_Collect(eval.edge_geom)) AS geom
FROM
(
SELECT
card.seq,
card.edge_id, card.edge_geom,
CASE WHEN card.is_dangle
THEN ROW_NUMBER() OVER(PARTITION BY card.seq, card.node_id, card.is_dangle ORDER BY card.edge_length DESC) - card.degree = 0
ELSE ROW_NUMBER() OVER(PARTITION BY card.seq, card.edge_id ORDER BY card.node_id) - card.degree = 0
END AS is_major
FROM
(
SELECT
conn.seq,
conn.edge_id,
conn.node_id,
conn.degree, MIN(conn.degree) OVER(PARTITION BY conn.seq, conn.edge_id) = 0 AS is_dangle,
conn.edge_length, conn.edge_geom
FROM
(
SELECT
reduce.seq,
edge.path[1] AS edge_id,
DENSE_RANK() OVER(PARTITION BY reduce.seq ORDER BY node.geom) AS node_id,
(COUNT(edge.path[1]) OVER(PARTITION BY reduce.seq, node.geom) > 1)::INT AS degree,
ST_Length(edge.geom::GEOGRAPHY) AS edge_length,
edge.geom AS edge_geom
FROM
reduce,
ST_Dump(reduce.geom) AS edge,
LATERAL ST_Boundary(edge.geom) AS nodes,
LATERAL ST_Dump(nodes) AS node
) AS conn
) AS card
) AS eval
WHERE
eval.is_major
GROUP BY
eval.seq
) AS merge
WHERE
ST_NumGeometries(merge.geom) > 0
)
SELECT
ST_LineMerge(ST_Collect(final.geom)) AS geom
FROM
(
SELECT
DENSE_RANK() OVER(PARTITION BY reduce.seq ORDER BY reduce.n DESC) = 1 AS keep_n,
ROW_NUMBER() OVER(PARTITION BY reduce.seq, reduce.n ORDER BY ST_Length(dmp.geom::GEOGRAPHY) DESC) <= 2 AS keep_len,
reduce.seq,
dmp.geom
FROM
reduce,
LATERAL ST_Dump(reduce.geom) AS dmp
WHERE
reduce.n <= 3
) AS final
WHERE
final.keep_n AND final.keep_len
;
$$
;
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment