wm/wm.sql

\set ON_ERROR_STOP on
SET plpgsql.extra_errors TO 'all';

-- detect_bends detects bends using the inflection angles. No corrections.
drop function if exists detect_bends;
create function detect_bends(line geometry, OUT bends geometry[]) as $$
declare
  pi real;
  p geometry;
  p1 geometry;
  p2 geometry;
  p3 geometry;
  bend geometry;
  prev_sign int4;
  cur_sign int4;
begin
  pi = radians(180);

  -- the last vertex is iterated over twice, because the algorithm uses 3
  -- vertices to calculate the angle between them.
  --
  -- Given 3 vertices p1, p2, p3:
  --
  --          p1___ ...
  --           /
  -- ... _____/
  --     p3   p2
  --
  -- When looping over the line, p1 will be head (lead) vertex, p2 will be the
  -- measured angle, and p3 will be trailing. The line that will be added to
  -- the bend will always be [p3,p2].
  -- So once the p1 becomes the last vertex, the loop terminates, and the
  -- [p2,p1] line will not have a chance to be added. So the loop adds the last
  -- vertex twice, so it has a chance to become p2, and be added to the bend.
  --
  for p in (
    (select geom from st_dumppoints(line) order by path[1] asc)
    union all
    (select geom from st_dumppoints(line) order by path[1] desc limit 1)
  ) loop
    p3 = p2;
    p2 = p1;
    p1 = p;
    continue when p3 is null;

    cur_sign = sign(pi - st_angle(p1, p2, p2, p3));

    if bend is null then
      bend = st_makeline(p3, p2);
    else
      bend = st_linemerge(st_union(bend, st_makeline(p3, p2)));
    end if;

    if prev_sign + cur_sign = 0 then
      if bend is not null then
        bends = bends || bend;
      end if;
      bend = st_makeline(p3, p2);
    end if;
    prev_sign = cur_sign;
  end loop;

  -- the last line may be lost if there is no "final" inflection angle. Add it.
  if (select count(1) >= 2 from st_dumppoints(bend)) then
    bends = bends || bend;
  end if;
end
$$ language plpgsql;

-- fix_gentle_inflections moves bend endpoints following "Gentle Inflection at
-- End of a Bend" section.
--
-- The text does not specify how many vertices can be "adjusted"; it can
-- equally be one or many. This function is adjusting many, as long as the
-- commulative inflection angle small (see variable below).
--
-- The implementation could be significantly optimized to avoid `st_reverse`
-- and array reversals, trading for complexity in fix_gentle_inflections1.
create or replace function fix_gentle_inflections(INOUT bends geometry[]) as $$
declare
  len int4;
  bends1 geometry[];
begin
  len = array_length(bends, 1);

  bends = fix_gentle_inflections1(bends);
  for i in 1..len loop
    bends1[i] = st_reverse(bends[len-i+1]);
  end loop;
  bends1 = fix_gentle_inflections1(bends1);

  for i in 1..len loop
    bends[i] = st_reverse(bends1[len-i+1]);
  end loop;
end
$$ language plpgsql;

-- fix_gentle_inflections1 fixes gentle inflections of an array of lines in
-- one direction. This is an implementation detail of fix_gentle_inflections.
drop function if exists fix_gentle_inflections1;
create function fix_gentle_inflections1(INOUT bends geometry[]) as $$
declare
  pi real;
  small_angle real;
  ptail geometry; -- tail point of tail bend
  phead geometry[]; -- 3 tail points of head bend
  i int4; -- bends[i] is the current head
begin
  pi = radians(180);
  -- the threshold when the angle is still "small", so gentle inflections can
  -- be joined
  small_angle := radians(30);

  for i in 2..array_length(bends, 1) loop
    -- Predicate: two bends will always share an edge. Assuming (A,B,C,D,E,F)
    -- is a bend:
    --           C________D
    --           /        \
    -- \________/          \_______/
    -- A       B           E       F
    --
    -- Then edges (A,B) and (E,F) are shared with the neighboring bends.
    --
    --
    -- Assume this curve (figure `inflection-1`), going clockwise from A:
    --
    --    \______B
    --    A      `-------. C
    --                   |
    --    G___ F         |
    --    /   `-----.____+ D
    --              E
    --
    -- After processing the curve following the definition of a bend, the bend
    -- [A-E] would be detected. Assuming inflection point E and F are "small",
    -- the bend needs to be extended by two edges to [A,G].
    select geom from st_dumppoints(bends[i-1])
      order by path[1] asc limit 1 into ptail;

    while true loop
      -- copy last 3 points of bends[i-1] (tail) to ptail
      select array(
        select geom from st_dumppoints(bends[i]) order by path[1] asc limit 3
      ) into phead;

      -- if the bend got too short, stop processing it
      exit when array_length(phead, 1) < 3;

      -- inflection angle between ptail[1:3] is "large", stop processing
      exit when abs(st_angle(phead[1], phead[2], phead[3]) - pi) > small_angle;

      -- distance from head's 1st vertex should be larger than from 2nd vertex
      exit when st_distance(ptail, phead[2]) < st_distance(ptail, phead[3]);

      -- Detected a gentle inflection.
      -- Move head of the tail to the tail of head
      bends[i] = st_removepoint(bends[i], 0);
      bends[i-1] = st_addpoint(bends[i-1], phead[3]);
    end loop;

  end loop;
end
$$ language plpgsql;

-- self_crossing eliminates self-crossing from the bends, following the
-- article's section "Self-line Crossing When Cutting a Bend".
drop function if exists self_crossing;
create function self_crossing(
  INOUT bends geometry[],
  OUT mutated boolean
) as $$
declare
  i int4;
  j int4;
  prev_length int4;
  pi real;
  angle real;
  p0 geometry;
  p1 geometry;
  p2 geometry;
  p3 geometry;
  a geometry;
  b geometry;
  bend geometry;
  multi geometry;
begin
  pi = radians(180);
  mutated = false;

  -- go through the bends and find one where sum of inflection angle is >180
  for i in 1..array_length(bends, 1) loop
    angle = 0;
    p1 = null;
    p2 = null;
    p3 = null;
    for p0 in (
      select geom from st_dumppoints(bends[i]) order by path[1] asc
    ) loop
      p3 = p2;
      p2 = p1;
      p1 = p0;
      continue when p3 is null;
      angle = angle + abs(pi - st_angle(p1, p2, p3));
    end loop;

    continue when abs(angle) <= pi;

    -- sum of inflection angles for this bend is >180, so it may be
    -- self-crossing. now try to find another bend in this line that
    -- crosses an imaginary line of end-vertices
    p0 = st_pointn(bends[i], 1);
    p1 = st_pointn(bends[i], -1);

    -- go through each bend in the given line, and see if has a potential to
    -- cross bends[i]. optimization: we care only about bends which beginning
    -- and end start at different sides of the plane, separated by endpoints
    -- p0 and p1.
    j = 0;
    while j < array_length(bends, 1) loop
      j = j + 1;
      continue when i = j;

      p2 = st_pointn(bends[j], 1);
      p3 = st_pointn(bends[j], -1);

      -- do end vertices of bend[i] cross bend[j]?
      a = st_pointn(bends[i], 1);
      b = st_pointn(bends[i], -1);
      multi = st_split(bends[j], st_makeline(a, b));
      continue when st_numgeometries(multi) = 1;
      continue when st_numgeometries(multi) = 2 and
        (st_contains(bends[j], a) or st_contains(bends[j], b));

      -- stars are aligned, we are changing the bend
      mutated = true;

      -- Sincere apologies to someone who will need to debug the block below.
      -- To understand it, I suggest you take a pencil and paper, draw a
      -- self-crossing bend (fig6 from the article works well), and figure out
      -- what happens here, by hand.
      prev_length = array_length(bends, 1);
      if j < i then
        -- remove first vertex of the following bend, because the last
        -- segment is always duplicated with the i-th bend.
        bends[i+1] = st_removepoint(bends[i+1], 0);
        bends[j] = st_geometryn(multi, 1);
        bends[j] = st_setpoint(
          bends[j], st_npoints(bends[j])-1,
          st_pointn(bends[i], st_npoints(bends[i]))
        );
        bends = bends[1:j] || bends[i+1:prev_length];
        j = i;
      else
        -- remove last vertex of the previous bend, because the last
        -- segment is duplicated with the i'th bend.
        bends[i-1] = st_removepoint(bends[i-1], st_npoints(bends[i-1])-1);
        bends[i] = st_makeline(
          st_pointn(bends[i], 1),
          st_removepoint(st_geometryn(multi, st_numgeometries(multi)), 0)
        );
        bends = bends[1:i] || bends[j+1:prev_length];
      end if;
      j = j - prev_length + array_length(bends, 1);
    end loop;
  end loop;
end
$$ language plpgsql;

drop function if exists bend_attrs;
create function bend_attrs(bends geometry[], dbg boolean default false) returns table(bend geometry, area real, cmp real, adjsize real) as $$
declare
  i int4;
  fourpi real;
  polygon geometry;
begin
  fourpi = 4*radians(180);
  for i in 1..array_length(bends, 1) loop
    bend = bends[i];
    select st_makepolygon(st_addpoint(bend, st_startpoint(bend))) into polygon;
    -- Compactness Index (cmp) is defined as "the ratio of the area of the
    -- polygon over the circle whose circumference length is the same as the
    -- length of the circumference of the polygon". I assume they meant the
    -- area of the circle. So here goes:
    -- 1. get polygon area P.
    -- 2. get polygon perimeter = u. Pretend it's our circle's circumference.
    -- 3. get A (area) of the circle from circumference u := (u^2)/(4*pi)
    -- 4. divide area by A: cmp = P/((u^2)*4*pi) = 4*pi*P/u^2
    select st_area(polygon) into area;
    select fourpi*area/(st_perimeter(polygon)^2) into cmp;
    if cmp > 0 then
      select (area*(0.75/cmp)) into adjsize;
    end if;
    if dbg then
      insert into debug_wm (name, way, props) values(
        'bend_attrs_' || i,
        polygon,
        json_build_object('area', area, 'cmp', cmp, 'adjsize', adjsize)
      );
    end if;
    return next;
  end loop;
end;
$$ language plpgsql;

-- ST_SimplifyWM simplifies a given geometry using Wang & Müller's
-- "Line Generalization Based on Analysis of Shape Characteristics" algorithm,
-- 1998.
drop function if exists ST_SimplifyWM;
create function ST_SimplifyWM(geom geometry, dbg boolean default false) returns geometry as $$
declare
  dbg_stage integer;
  i integer;
  line geometry;
  lines geometry[];
  bends geometry[];
  mutated boolean;
  l_type text;
begin
  l_type = st_geometrytype(geom);
  if l_type = 'ST_LineString' then
    lines = array[geom];
  elseif l_type = 'ST_MultiLineString' then
    lines = array((select a.geom from st_dump(geom) a order by path[1] asc));
  else
    raise 'Unknown geometry type %', l_type;
  end if;

  for i in 1..array_length(lines, 1) loop
    mutated = true;
    dbg_stage = 1;
    while mutated loop
      if dbg then
        insert into debug_wm (name, way) values(
          dbg_stage || 'afigures_' || i,
          lines[i]
        );
      end if;

      bends = detect_bends(lines[i]);

      if dbg then
        insert into debug_wm(name, way) values(
          dbg_stage || 'bbends_' || i || '_' || generate_subscripts(bends, 1),
          unnest(bends)
        );
      end if;

      bends = fix_gentle_inflections(bends);

      if dbg then
        insert into debug_wm(name, way) values(
          dbg_stage || 'cinflections' || i || '_' || generate_subscripts(bends, 1),
          unnest(bends)
        );
      end if;

      select * from self_crossing(bends) into bends, mutated;

      if dbg then
        insert into debug_wm(name, way) values(
          dbg_stage || 'dcrossings' || i || '_' || generate_subscripts(bends, 1),
          unnest(bends)
        );
      end if;

      if mutated then
        lines[i] = st_linemerge(st_union(bends));
        dbg_stage = dbg_stage + 1;
        continue;
      end if;

      -- self-crossing mutations are done, calculate bend properties
      perform bend_attrs(bends, true);
    end loop;

  end loop;

  if l_type = 'ST_LineString' then
    return st_linemerge(st_union(lines));
  elseif l_type = 'ST_MultiLineString' then
    return st_union(lines);
  end if;
end
$$ language plpgsql;