PostGIS shape descriptor, "as rectangle" metrics vector

pgis_shape_descriptor.sql

CREATE FUNCTION shapedescr_sizes(
    -- Shape-descriptor "as rectangle" for geometry description by sizes. 
    gbase geometry,              -- input
    -- p_seqs integer DEFAULT 8,    -- deprecated? for st_buffer(g,w,p_seqs) or point-buffer inference 
    -- p_shape varchar DEFAULT '',  -- will be endcap indicator
    p_decplacesof_zero integer DEFAULT 6, -- precision of zero when rounding delta
    p_dwmin float DEFAULT 99999999.0,     -- change to ex. 0.0001, if to use.
    p_deltaimpact float DEFAULT 9999.0    -- internal (maximized by probability of negative delta)
)  RETURNS float[] AS $f$
  DECLARE
    ret float[];
    dw float;
    b float;
    L_estim float;
    H_estim float;
    aorig float;
    gaux geometry;
    g1 geometry;
    A0 float;
    A1 float;
    c float;
    delta float;
    per float;
    errcod float;
  BEGIN
    errcod=0.0;
    IF gbase IS NULL OR NOT(ST_IsClosed(gbase)) THEN
        errcod=1;                  -- ERROR1 (die)
        RAISE EXCEPTION 'error %: invalid input geometry',errcod;
    END IF;
    A0 := ST_Area(gbase);
    per := st_perimeter(gbase);
    dw := sqrt(A0)/p_deltaimpact;
    IF dw>p_dwmin THEN dw:=p_dwmin; END IF;
    g1 = ST_Buffer(gbase,dw);
    A1 = ST_area(g1);
    IF A0>A1 THEN 
        errcod=10;                 -- ERROR2 (die)
        RAISE EXCEPTION 'error %: invalid buffer/geometry with A0=% g.t. A1=%',errcod,A0,A1;
    END IF;
    IF (A1-A0)>1.001*dw*per THEN 
        gaux := ST_Buffer(g1,-dw);  -- closing operation.
        A0 = ST_Area(gaux);         -- changed area
        per := ST_Perimeter(gaux);  -- changed
        errcod:=errcod + 0.1;       -- Warning3
    END IF;
    C := 2.0*dw;
    b := -(A1-A0)/C+C;
    delta := b^2-4.0*A0;
    IF delta<0.0 AND round(delta,p_decplacesof_zero)<=0.0 THEN
           delta=0.0; -- for regular shapes like the square
           errcod:=errcod + 0.01;  -- Warning2
    END IF;
    IF delta<0.0 THEN
        L_estim := NULL; 
        H_estim := NULL;
        errcod:=errcod+100;        -- ERROR3
    ELSE
        L_estim := (-b + sqrt(delta))/2.0; 
        H_estim := (-b - sqrt(delta))/2.0;
    END IF;
    IF abs(A0-L_estim*H_estim)>0.001 THEN 
        errcod:=errcod + 0.001;    -- Warning1
    END IF;
    ret := array[L_estim,H_estim,a0,per,dw,errcod];
    return ret;
  END
$f$ LANGUAGE plpgsql IMMUTABLE;


CREATE or replace FUNCTION shapedescr_sizes_tr(
    -- Translator for shapedescr_sizes(). Uses ROUND(float).
    float[],               -- shapedescr_sizes() returned vector
    integer DEFAULT 0,     -- general round parameter
    integer DEFAULT 3      -- number of "decimal warnings" 
) RETURNS varchar[] -- length, height, area, perimeter, dw, radius, err_message
AS $f$
   SELECT array[
      round(L,$2+1)::varchar, round(H,$2+1)::varchar, round(area,$2)::varchar,
      round(perim,$2+1)::varchar, round(dw,$2+3)::varchar,
      round(sqrt(L*H/pi()),$2+1)::varchar,  -- radius for "shape as circle" 
      CASE WHEN errcod>3.0 THEN 'ERROR '  ||round(errcod-$3) 
           WHEN errcod>0.0  THEN 'WARNING '||round(errcod)||CASE
                   WHEN round(10^(-errcod),$3-1)!=($1)[6] THEN '.'|| ($1)[6]*10^$3
                   ELSE ''
		END
           ELSE '' 
      END::varchar
   ]
   FROM (
        SELECT ($1)[1] as L, ($1)[2] as H, ($1)[3] as area, ($1)[4] as perim, 
               ($1)[5] as dw, log(($1)[6]*10.0^($3+1)+1.0) as errcod
   ) as t; 
$f$ LANGUAGE SQL IMMUTABLE;


CREATE FUNCTION shapedescr_sizes_tr(geometry, integer DEFAULT 0, integer DEFAULT 3) 
RETURNS varchar[] AS $f$
   SELECT shapedescr_sizes_tr(shapedescr_sizes($1),$2,$3);
$f$ LANGUAGE SQL IMMUTABLE;


-- visual check:
SELECT shapedescr_sizes_tr(ST_MakePolygon(ST_GeomFromText('LINESTRING(75 30,77 29,78 29.5, 75 30)'))); -- 0.6

-- simulate areas and its central lines, and calculates estimatives.
SELECT tipo, round("L",1) as "L", 2*w as "H", 
       round(descr[1],1) as "L_estim", round(descr[2],1) as "H_estim",
       round(ST_Area(g0)/"L",1) as "HbyL_Area", round(0.5*ST_Perimeter(g0)-"L",1) as "HbyL_perim" 
FROM (
 SELECT *, st_length(g) as "L", shapedescr_sizes(g0) AS descr
 FROM (
      SELECT *, st_buffer(g,w,tipo) as g0
      FROM (
        SELECT ST_GeomFromText('LINESTRING(50 50,150 150,160 180,200 200)') AS g,
	 unnest(array[2.0,10.0,20.0,50.0]) AS w
      ) AS t, (
          SELECT unnest(array['','endcap=square','endcap=flat join=bevel']) AS tipo
      ) AS t2
 ) as t3
) as t4;

-- tests with database:
--SELECT shapedescr_sizes_tr(the_geom) as descr FROM quadras LIMIT 100; 

test_second_method.sql

-- -- -- -- -- -- -- -- -- -- -- 
-- -- -- -- -- -- -- -- -- -- -- 
-- ADDING TESTS OF SECOND METHOD

CREATE or replace FUNCTION shapedescr_sizes_test(
    -- TESTING grow of precision: result NO.
    -- Shape-descriptor "as rectangle" for geometry description by sizes. 
    gbase geometry,              -- input
    -- p_seqs integer DEFAULT 8,    -- deprecated? for st_buffer(g,w,p_seqs) or point-buffer inference 
    -- p_shape varchar DEFAULT '',  -- will be endcap indicator
    p_decplacesof_zero integer DEFAULT 4, -- precision of zero when rounding delta
    p_dwmin float DEFAULT 99999999.0,     -- change to ex. 0.0001, if to use.
    p_deltaimpact float DEFAULT 9999.0    -- internal (maximized by probability of negative delta)
)  RETURNS float[] AS $f$
  DECLARE
    ret float[];
    dw float;
    b float;
    L_estim float;
    H_estim float;
    aorig float;
    gaux geometry;
    g1 geometry;
    A0 float;
    A1 float;
    A2 float;
    delta float;
    per float;
    errcod float;
  BEGIN
    errcod=0.0;
    IF gbase IS NULL OR NOT(ST_IsClosed(gbase)) THEN
        errcod=1;                  -- ERROR1 (die)
        RAISE EXCEPTION 'error %: invalid input geometry',errcod;
    END IF;
    A0 := ST_Area(gbase);
    per := st_perimeter(gbase); -- only for dw*per cheching, and for output 
    dw := sqrt(A0)/p_deltaimpact;
    IF dw>p_dwmin THEN dw:=p_dwmin; END IF;
    g1 = ST_Buffer(gbase,dw);
    A1 = ST_area(g1);
    IF A0>A1 THEN 
        errcod=10;                 -- ERROR2 (die)
        RAISE EXCEPTION 'error %: invalid buffer/geometry with A0=% g.t. A1=%',errcod,A0,A1;
    END IF;
    IF (A1-A0)>1.001*dw*per THEN 
        gaux := ST_Buffer(g1,-dw);  -- closing operation.
        A0 = ST_Area(gaux);         -- changed area
        per := ST_Perimeter(gaux);  -- changed
        errcod:=errcod + 0.1;       -- Warning3
        A2 = ST_Area(ST_Buffer(gaux,-dw));         -- second changed area
    ELSE
        A2 = ST_Area(ST_Buffer(gbase,-dw));
    END IF;
    -- C := 2.0*dw;
    b := -(A1-A2)/(4.0*dw); -- -(A1-A0)/C+C;
    delta := b^2-4.0*A0;
    IF delta<0.0 AND round(delta,p_decplacesof_zero)<=0.0 THEN
           delta=0.0; -- for regular shapes like the square
           errcod:=errcod + 0.01;  -- Warning2
    END IF;
    IF delta<0.0 THEN
        L_estim := NULL; 
        H_estim := NULL;
        errcod:=errcod+100;        -- ERROR3
    ELSE
        L_estim := (-b + sqrt(delta))/2.0; 
        H_estim := (-b - sqrt(delta))/2.0;
    END IF;
    IF abs(A0-L_estim*H_estim)>0.001 THEN 
        errcod:=errcod + 0.001;    -- Warning1
    END IF;
    ret := array[L_estim,H_estim,a0,per,dw,errcod];
    return ret;
  END
$f$ LANGUAGE plpgsql IMMUTABLE;


-- COMPARING METHODS 
SELECT tipo, round("L",1) as "L", 2*w as "H", 
       round(descr[1],3) as "L_estim", round(descr[2],3) as "H_estim",
       round(descr2[1],3) as "L_estim_v2", round(descr2[2],3) as "H_estim_v2"
FROM (
 SELECT *, st_length(g) as "L", shapedescr_sizes(g0) AS descr,
        shapedescr_sizes_test(g0) AS descr2
 FROM (
      SELECT *, st_buffer(g,w,tipo) as g0
      FROM (
        SELECT ST_GeomFromText('LINESTRING(50 50,150 150,160 180,200 200)') AS g,
	 unnest(array[2.0,10.0,20.0,50.0]) AS w, 0.005 as dw
      ) AS t, (
          SELECT unnest(array['','endcap=square','endcap=flat join=bevel']) AS tipo
      ) AS t2
 ) as t3
) as t4;  -- result = very little gain, NOT justify heavy shapedescr_sizes_test().