ST_RedistributeVertices

ST_RedistributeVertices.sql

CREATE OR REPLACE FUNCTION ST_RedistributeVertices(
  geom geometry,
  distance double precision,
  null_for_small boolean DEFAULT false)
RETURNS geometry AS $body$

DECLARE
  perimeter numeric;
  num_vert integer;

BEGIN
  IF distance < 0 THEN
    RAISE EXCEPTION 'distance must be greater than zero, found %', distance;
  END IF;

  IF GeometryType(geom) ~ '^LINESTRING' THEN
    num_vert := round(ST_Length(geom)::numeric / distance::numeric);
    -- RAISE NOTICE 'a %, b %, a/b %, round(a/b) %', ST_Length(geom), distance, ST_Length(geom) / distance, round(ST_Length(geom)::numeric / distance::numeric);
    IF num_vert = 0 THEN
      IF null_for_small THEN
        RETURN NULL;
      END IF;
      num_vert := 1;
    END IF;
    -- RAISE NOTICE 'num_vert %', num_vert;

    RETURN ST_MakeLine(array_agg(
      ST_Line_Interpolate_Point(geom, n::numeric / num_vert::numeric)
      )) FROM generate_series(0, num_vert) AS n;

  ELSIF GeometryType(geom) ~ '^POLYGON' AND ST_NumInteriorRings(geom) = 0 THEN
    perimeter := ST_Length(ST_ExteriorRing(geom));
    num_vert := round(perimeter / distance::numeric);
    -- RAISE NOTICE 'poly num_vert %', num_vert;
    IF num_vert < 3 THEN
      IF null_for_small THEN
        RETURN NULL;
      END IF;
      -- RAISE NOTICE 'distance1 %', distance;
      distance := perimeter / 3.0;
      num_vert := round(perimeter / distance::numeric);
      -- RAISE NOTICE 'a %, b %, a/b %, round(a/b) %', perimeter, distance, perimeter / distance, round(perimeter::numeric / distance::numeric);
    END IF;

    RETURN ST_MakePolygon(ST_RedistributeVertices(ST_ExteriorRing(geom), distance, null_for_small));

  ELSIF GeometryType(geom) ~ '^POLYGON' AND ST_NumInteriorRings(geom) >= 1 THEN
    
    RETURN ST_MakePolygon(ST_ExteriorRing(ST_RedistributeVertices(ST_MakePolygon(ST_ExteriorRing(geom)), distance, null_for_small)),
                          array_agg(interior))
      FROM (
        SELECT ST_ExteriorRing(ST_RedistributeVertices(ST_MakePolygon(ST_InteriorRingN(geom, n)), distance, null_for_small)) AS interior
        FROM generate_series(1, ST_NumInteriorRings(geom)) AS n
      ) AS f
      WHERE interior NOTNULL;

  ELSIF GeometryType(geom) ~ '^(MULTI|COLLECTION)' THEN

    RETURN ST_Multi(ST_Collect(g))
      FROM (SELECT ST_RedistributeVertices((ST_Dump(geom)).geom, distance, null_for_small) AS g) AS f;

  ELSE
    RAISE EXCEPTION 'Unhandled geometry %', GeometryType(geom);
  END IF;
END;
$body$ LANGUAGE plpgsql IMMUTABLE STRICT;



-- Tests

WITH

  -- Test geometries
  line10 AS (
    SELECT 'LINESTRING(0 0, 0 10)'::geometry AS geom),

  quad AS (
    SELECT 'POLYGON((30 10, 40 40, 20 40, 10 20, 30 10))'::geometry AS geom),

  equilateral_triangle AS (
    SELECT ST_MakePolygon(ST_MakeLine(
      ARRAY[ST_MakePoint(0, 0),ST_MakePoint(0.5, sqrt(3)/2),ST_MakePoint(1, 0),ST_MakePoint(0, 0)])) AS geom),

  right_triangle AS (
    SELECT 'POLYGON((0 0, 0 1, 1 0, 0 0))'::geometry AS geom),

  quad_with_hole AS (
    SELECT 'POLYGON ((35 10, 45 45, 15 40, 10 20, 35 10),(20 30, 35 35, 30 20, 20 30))'::geometry AS geom)

  -- LineString tests
  SELECT 'same'::text AS name,
    'LINESTRING(0 0, 0 10)'::text AS expected,
    ST_AsText(ST_RedistributeVertices(geom, 10))::text AS geom
    FROM line10

  UNION ALL SELECT 'refined by two',
    'LINESTRING(0 0,0 5,0 10)',
    ST_AsText(ST_RedistributeVertices(geom, 6))
    FROM line10

  UNION ALL SELECT '*same',
    'LINESTRING(0 0, 0 10)' AS expected,
    ST_AsText(ST_RedistributeVertices(geom, 20, true))
    FROM line10

  UNION ALL SELECT 'null', NULL,
    ST_AsText(ST_RedistributeVertices(geom, 20.1, true))
    FROM line10

  -- Polygon tests
  UNION ALL SELECT 'equal distance',
    '',
    ST_AsText(ST_RedistributeVertices(geom, (2 + sqrt(2)) / 3))
    FROM right_triangle

  UNION ALL SELECT 'division by zero?',
    '',
    ST_AsText(ST_RedistributeVertices(geom, 2.0/3.0))
    FROM right_triangle

  UNION ALL SELECT 'quad to tri',
    ST_NumPoints(ST_ExteriorRing(geom)) || '->' || ST_NumPoints(ST_ExteriorRing(ST_RedistributeVertices(geom, 30))),
    ST_AsText(ST_RedistributeVertices(geom, 30))
    FROM quad

  UNION ALL SELECT 'null',
    NULL,
    ST_AsText(ST_RedistributeVertices(geom, 40, true))
    FROM quad

  UNION ALL SELECT 'modify distance',
    '',
    ST_AsText(ST_RedistributeVertices(geom, 40))
    FROM quad

  UNION ALL SELECT 'same',
    'POLYGON((0 0,0.5 0.866025403784439,1 0,0 0))',
    ST_AsText(ST_RedistributeVertices(geom, 1))
    FROM equilateral_triangle

  UNION ALL SELECT 'double resolution',
    'POLYGON((0 0,0.25 0.433012701892219,0.5 0.866025403784439,0.75 0.433012701892219,1 0,0.5 0,0 0))',
    ST_AsText(ST_RedistributeVertices(geom, 0.5))
    FROM equilateral_triangle

  UNION ALL SELECT 'same',
    'same',
    ST_AsText(ST_RedistributeVertices(geom, 2))
    FROM equilateral_triangle

  UNION ALL SELECT 'null',
    NULL,
    ST_AsText(ST_RedistributeVertices(geom, 2, true))
    FROM equilateral_triangle

  -- Polygon with holes
  UNION ALL SELECT 'null',
    NULL,
    ST_AsText(ST_RedistributeVertices(geom, 2, true))
    FROM quad_with_hole