From 2c7be0989a49f2a871f8b3ebe1e6c5bfa34ce7b4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Motiejus=20Jak=C5=A1tys?= <motiejus@jakstys.lt>
Date: Wed, 19 May 2021 22:57:50 +0300
Subject: [PATCH] split bendattrs to bends and attrs

this makes code easier to reason about, and some lines shorter (thus
removing splitting requirement).
---
 wm.sql | 286 ++++++++++++++++++++++-----------------------------------
 1 file changed, 108 insertions(+), 178 deletions(-)

diff --git a/wm.sql b/wm.sql
index 58d9a38..0c962c6 100644
--- a/wm.sql
+++ b/wm.sql
@@ -81,16 +81,14 @@ begin
       dbgpolygon = null;
       if st_npoints(bends[i]) >= 3 then
           dbgpolygon = st_makepolygon(
-            st_addpoint(bends[i],
-              st_startpoint(bends[i]))
+            st_addpoint(bends[i], st_startpoint(bends[i]))
           );
       end if;
       insert into wm_debug(stage, name, gen, nbend, way) values(
         'bbends-polygon', dbgname, dbggen, i, dbgpolygon);
     end loop;
   end if;
-end
-$$ language plpgsql;
+end $$ language plpgsql;
 
 -- wm_fix_gentle_inflections moves bend endpoints following "Gentle Inflection
 -- at End of a Bend" section.
@@ -141,8 +139,7 @@ begin
         'cinflections-polygon', dbgname, dbggen, i, dbgpolygon);
     end loop;
   end if;
-end
-$$ language plpgsql;
+end $$ language plpgsql;
 
 -- wm_fix_gentle_inflections1 fixes gentle inflections of an array of lines in
 -- one direction. An implementation detail of wm_fix_gentle_inflections.
@@ -204,8 +201,7 @@ begin
     end loop;
 
   end loop;
-end
-$$ language plpgsql;
+end $$ language plpgsql;
 
 -- wm_if_selfcross returns whether baseline of bendi crosses bendj.
 -- If it doesn't, returns a null geometry.
@@ -235,9 +231,7 @@ begin
   end if;
 
   return multi;
-end
-$$ language plpgsql;
-
+end $$ language plpgsql;
 
 -- wm_self_crossing eliminates self-crossing from the bends, following the
 -- article's section "Self-line Crossing When Cutting a Bend".
@@ -305,8 +299,7 @@ begin
       unnest(bends)
     );
   end if;
-end
-$$ language plpgsql;
+end $$ language plpgsql;
 
 drop function if exists wm_inflection_angle;
 create function wm_inflection_angle (IN bend geometry, OUT angle real) as $$
@@ -324,17 +317,13 @@ begin
     continue when p3 is null;
     angle = angle + abs(pi() - st_angle(p1, p2, p3));
   end loop;
-end
-$$ language plpgsql;
+end $$ language plpgsql;
 
 drop function if exists wm_bend_attrs;
-drop function if exists wm_isolated_bends;
 drop function if exists wm_elimination;
 drop function if exists wm_exaggeration;
-drop type if exists wm_t_bend_attrs;
-create type wm_t_bend_attrs as (
-  bend geometry,
-  area real,
+drop type if exists wm_t_attrs;
+create type wm_t_attrs as (
   adjsize real,
   baselinelength real,
   curvature real,
@@ -344,40 +333,66 @@ create function wm_bend_attrs(
   bends geometry[],
   dbgname text default null,
   dbggen integer default null
-) returns setof wm_t_bend_attrs as $$
+) returns wm_t_attrs[] as $$
 declare
-  cmp float;
-  i int4;
-  polygon geometry;
+  isolation_threshold constant real default 0.5;
+  attrs wm_t_attrs[];
+  attr wm_t_attrs;
   bend geometry;
-  res wm_t_bend_attrs;
+  i int4;
+  needs_curvature real;
+  skip_next boolean;
+  dbglastid integer;
 begin
   for i in 1..array_length(bends, 1) loop
     bend = bends[i];
-    res = null;
-    res.bend = bend;
-    res.adjsize = 0;
-    res.baselinelength = st_distance(st_startpoint(bend), st_endpoint(bend));
-    res.curvature = wm_inflection_angle(bend) / st_length(bend);
-    res.isolated = false;
+    attr.adjsize = 0;
+    attr.baselinelength = st_distance(st_startpoint(bend), st_endpoint(bend));
+    attr.curvature = wm_inflection_angle(bend) / st_length(bend);
+    attr.isolated = false;
     if st_numpoints(bend) >= 3 then
-      res.adjsize = wm_adjsize(bend);
+      attr.adjsize = wm_adjsize(bend);
     end if;
+    attrs[i] = attr;
+  end loop;
 
+  for i in 1..array_length(attrs, 1) loop
     if dbgname is not null then
       insert into wm_debug (stage, name, gen, nbend, way, props) values(
         'ebendattrs', dbgname, dbggen, i, bend,
         jsonb_build_object(
-          'adjsize', res.adjsize,
-          'baselinelength', res.baselinelength,
-          'curvature', res.curvature
+          'adjsize', attrs[i].adjsize,
+          'baselinelength', attrs[i].baselinelength,
+          'curvature', attrs[i].curvature,
+          'isolated', false
         )
-      );
+      ) returning id into dbglastid;
+    end if;
+
+    -- first and last bends can never be isolated by definition
+    if skip_next or i = 1 or i = array_length(attrs, 1) then
+      skip_next = false;
+      continue;
+    end if;
+
+    needs_curvature = attrs[i].curvature * isolation_threshold;
+    if attrs[i-1].curvature < needs_curvature and
+       attrs[i+1].curvature < needs_curvature then
+      attr = attrs[i];
+      attr.isolated = true;
+      attrs[i] = attr;
+      skip_next = true;
+
+      if dbgname is not null then
+        update wm_debug
+        set props=props || jsonb_build_object('isolated', true)
+        where id=dbglastid;
+      end if;
     end if;
-    return next res;
   end loop;
-end;
-$$ language plpgsql;
+
+  return attrs;
+end $$ language plpgsql;
 
 -- sm_st_split a line by a point in a more robust way than st_split.
 -- See https://trac.osgeo.org/postgis/ticket/2192
@@ -448,8 +463,7 @@ begin
     bend = st_setsrid(st_makeline(points), st_srid(bend));
     size = wm_adjsize(bend);
   end loop;
-end
-$$ language plpgsql;
+end $$ language plpgsql;
 
 -- wm_adjsize calculates adjusted size for a polygon. Can return 0.
 drop function if exists wm_adjsize;
@@ -474,12 +488,12 @@ begin
   if cmp > 0 then
     adjsize = (area*(0.75/cmp));
   end if;
-end
-$$ language plpgsql;
+end $$ language plpgsql;
 
 -- wm_exaggeration is the Exaggeration Operator described in the WM paper.
 create function wm_exaggeration(
-  INOUT bendattrs wm_t_bend_attrs[],
+  INOUT bends geometry[],
+  attrs wm_t_attrs[],
   dhalfcircle float,
   intersect_patience integer,
   dbgname text default null,
@@ -488,7 +502,7 @@ create function wm_exaggeration(
 ) as $$
 declare
   desired_size constant float default pi()*(dhalfcircle^2)/8;
-  tmpbendattr wm_t_bend_attrs;
+  bend geometry;
   tmpint geometry;
   i integer;
   n integer;
@@ -496,79 +510,69 @@ declare
 begin
   mutated = false;
   <<bendloop>>
-  for i in 1..array_length(bendattrs, 1) loop
-    if bendattrs[i].isolated and bendattrs[i].adjsize < desired_size then
-      tmpbendattr.bend = wm_exaggerate_bend(
-        bendattrs[i].bend,
-        bendattrs[i].adjsize,
-        desired_size
-      );
+  for i in 1..array_length(attrs, 1) loop
+    if attrs[i].isolated and attrs[i].adjsize < desired_size then
+      bend = wm_exaggerate_bend(bends[i], attrs[i].adjsize, desired_size);
 
-      -- does tmpbendattrs.bend intersect with the previous or next
+      -- does bend intersect with the previous or next
       -- intersect_patience bends? If they do, abort exaggeration for this one.
 
-
       -- Do close-by bends intersect with this one? Special
       -- handling first, because 2 vertices need to be removed before checking.
-      n = st_npoints(bendattrs[i-1].bend);
+      n = st_npoints(bends[i-1]);
       if n > 3 then
-        continue when st_intersects(tmpbendattr.bend,
-          st_removepoint(st_removepoint(bendattrs[i-1].bend, n-1), n-2));
+        continue when st_intersects(bend,
+          st_removepoint(st_removepoint(bends[i-1], n-1), n-2));
       end if;
 
-      n = st_npoints(bendattrs[i+1].bend);
+      n = st_npoints(bends[i+1]);
       if n > 3 then
-        continue when st_intersects(tmpbendattr.bend,
-          st_removepoint(st_removepoint(bendattrs[i+1].bend, 0), 0));
+        continue when st_intersects(bend,
+          st_removepoint(st_removepoint(bends[i+1], 0), 0));
       end if;
 
       for n in -intersect_patience+1..intersect_patience-1 loop
         continue when n in (-1, 0, 1);
         continue when i+n < 1;
-        continue when i+n > array_length(bendattrs, 1);
+        continue when i+n > array_length(attrs, 1);
 
         -- More special handling: if the neigbhoring bend has 3 vertices, the
         -- neighbor's neighbor may just touch the tmpbendattr.bend; in this
         -- case, the nearest vertex should be removed before comparing.
-        tmpint = bendattrs[i+n].bend;
+        tmpint = bends[i+n];
         if st_npoints(tmpint) > 2 then
-          if n = -2 and st_npoints(bendattrs[i+n+1].bend) = 3 then
+          if n = -2 and st_npoints(bends[i+n+1]) = 3 then
             tmpint = st_removepoint(tmpint, st_npoints(tmpint)-1);
-          elsif n = 2 and st_npoints(bendattrs[i+n-1].bend) = 3 then
+          elsif n = 2 and st_npoints(bends[i+n-1]) = 3 then
             tmpint = st_removepoint(tmpint, 0);
           end if;
         end if;
 
-        continue bendloop when st_intersects(tmpbendattr.bend, tmpint);
+        continue bendloop when st_intersects(bend, tmpint);
       end loop;
 
       -- No intersections within intersect_patience, mutate bend!
       mutated = true;
-      bendattrs[i] = tmpbendattr;
+      bends[i] = bend;
 
       -- remove last vertex of the previous bend and first vertex of the next
       -- bend, because bends always share a line segment together this is
       -- duplicated in a few places, because PostGIS does not allow (?)
       -- mutating an array when passed to a function.
-      tmpbendattr.bend = st_removepoint(
-        bendattrs[i-1].bend,
-        st_npoints(bendattrs[i-1].bend)-1
-      );
-      bendattrs[i-1] = tmpbendattr;
-      tmpbendattr.bend = st_removepoint(bendattrs[i+1].bend, 0);
-      bendattrs[i+1] = tmpbendattr;
+      bends[i-1] = st_removepoint(bends[i-1], st_npoints(bends[i-1])-1);
+      bends[i+1] = st_removepoint(bends[i+1], 0);
 
       if dbgname is not null then
         insert into wm_debug (stage, name, gen, nbend, way) values(
-          'gexaggeration', dbgname, dbggen, i, bendattrs[i].bend
-        );
+          'gexaggeration', dbgname, dbggen, i, bends[i]);
       end if;
     end if;
   end loop;
 end $$ language plpgsql;
 
 create function wm_elimination(
-  INOUT bendattrs wm_t_bend_attrs[],
+  INOUT bends geometry[],
+  attrs wm_t_attrs[],
   dhalfcircle float,
   dbgname text default null,
   dbggen integer default null,
@@ -579,119 +583,55 @@ declare
   leftsize float;
   rightsize float;
   i int4;
-  j int4;
-  tmpbendattr wm_t_bend_attrs;
-  dbgbends geometry[];
 begin
   mutated = false;
 
   i = 1;
-  while i < array_length(bendattrs, 1)-1 loop
+  while i < array_length(attrs, 1)-1 loop
     i = i + 1;
-    continue when bendattrs[i].adjsize = 0;
-    continue when bendattrs[i].adjsize > desired_size;
+    continue when attrs[i].adjsize = 0;
+    continue when attrs[i].adjsize > desired_size;
 
     if i = 2 then
-      leftsize = bendattrs[i].adjsize + 1;
+      leftsize = attrs[i].adjsize + 1;
     else
-      leftsize = bendattrs[i-1].adjsize;
+      leftsize = attrs[i-1].adjsize;
     end if;
 
-    if i = array_length(bendattrs, 1)-1 then
-      rightsize = bendattrs[i].adjsize + 1;
+    if i = array_length(attrs, 1)-1 then
+      rightsize = attrs[i].adjsize + 1;
     else
-      rightsize = bendattrs[i+1].adjsize;
+      rightsize = attrs[i+1].adjsize;
     end if;
 
-    continue when bendattrs[i].adjsize >= leftsize;
-    continue when bendattrs[i].adjsize >= rightsize;
+    continue when attrs[i].adjsize >= leftsize;
+    continue when attrs[i].adjsize >= rightsize;
 
     -- Local minimum. Elminate bend!
     mutated = true;
-    tmpbendattr.bend = st_makeline(
-      st_pointn(bendattrs[i].bend, 1),
-      st_pointn(bendattrs[i].bend, -1)
-    );
-    bendattrs[i] = tmpbendattr;
+    bends[i] = st_makeline(st_pointn(bends[i], 1), st_pointn(bends[i], -1));
+
     -- remove last vertex of the previous bend and
     -- first vertex of the next bend, because bends always
     -- share a line segment together
-    tmpbendattr.bend = st_removepoint(
-      bendattrs[i-1].bend,
-      st_npoints(bendattrs[i-1].bend)-1
-    );
-    bendattrs[i-1] = tmpbendattr;
-    tmpbendattr.bend = st_removepoint(bendattrs[i+1].bend, 0);
-    bendattrs[i+1] = tmpbendattr;
+    bends[i-1] = st_removepoint(bends[i-1], st_npoints(bends[i-1])-1);
+    bends[i+1] = st_removepoint(bends[i+1], 0);
     -- the next bend's adjsize is now messed up; it should not be taken
     -- into consideration for other local minimas. Skip over 2.
     i = i + 2;
   end loop;
 
   if dbgname is not null then
-    for j in 1..array_length(bendattrs, 1) loop
-      dbgbends[j] = bendattrs[j].bend;
-    end loop;
-
     insert into wm_debug(stage, name, gen, nbend, way) values(
       'helimination',
       dbgname,
       dbggen,
-      generate_subscripts(dbgbends, 1),
-      unnest(dbgbends)
+      generate_subscripts(bends, 1),
+      unnest(bends)
     );
   end if;
-end
-$$ language plpgsql;
+end $$ language plpgsql;
 
-create function wm_isolated_bends(
-  INOUT bendattrs wm_t_bend_attrs[],
-  dbgname text default null,
-  dbggen integer default null
-) as $$
-declare
-  -- if neighbor's curvatures are within this fraction of the current bend
-  isolation_threshold constant real default 0.5;
-  this real;
-  skip_next bool;
-  res wm_t_bend_attrs;
-  i int4;
-  last_id integer;
-begin
-  for i in 1..array_length(bendattrs, 1) loop
-    if dbgname is not null then
-      insert into wm_debug (stage, name, gen, nbend, way, props) values(
-        'fisolated_bends', dbgname, dbggen, i, bendattrs[i].bend,
-        jsonb_build_object('isolated', false)
-      ) returning id into last_id;
-    end if;
-    -- first and last bends cannot be isolated
-    if i = 1 or i = array_length(bendattrs, 1) then
-      continue;
-    end if;
-
-    res = bendattrs[i];
-    if skip_next then
-      skip_next = false;
-    else
-      this = bendattrs[i].curvature * isolation_threshold;
-      if bendattrs[i-1].curvature < this and
-         bendattrs[i+1].curvature < this then
-        res.isolated = true;
-        bendattrs[i] = res;
-        skip_next = true;
-
-        if dbgname is not null then
-          update wm_debug
-          set props=props || jsonb_build_object('isolated', true)
-          where id=last_id;
-        end if;
-      end if;
-    end if;
-
-  end loop;
-end
-$$ language plpgsql;
 
 drop function if exists ST_SimplifyWM_Estimate;
 create function ST_SimplifyWM_Estimate(
@@ -717,8 +657,7 @@ begin
     npoints = npoints + st_numpoints(lines[i]);
   end loop;
   secs = npoints / 33;
-end
-$$ language plpgsql;
+end $$ language plpgsql;
 
 -- ST_SimplifyWM simplifies a given geometry using Wang & Müller's
 -- "Line Generalization Based on Analysis of Shape Characteristics" algorithm,
@@ -742,7 +681,7 @@ declare
   line geometry;
   lines geometry[];
   bends geometry[];
-  bendattrs wm_t_bend_attrs[];
+  attrs wm_t_attrs[];
   mutated boolean;
   l_type text;
 begin
@@ -773,31 +712,23 @@ begin
 
       select * from wm_self_crossing(bends, dbgname, gen) into bends, mutated;
 
-      if mutated then
-        lines[i] = st_linemerge(st_union(bends));
-        gen = gen + 1;
-        continue;
+      if not mutated then
+        attrs = wm_bend_attrs(bends, dbgname, gen);
+
+        select * from wm_exaggeration(bends, attrs,
+          dhalfcircle, intersect_patience, dbgname, gen) into bends, mutated;
       end if;
 
-      bendattrs = array((select wm_bend_attrs(bends, dbgname, gen)));
-      bendattrs = wm_isolated_bends(bendattrs, dbgname, gen);
-
-      select * from wm_exaggeration(
-        bendattrs, dhalfcircle, intersect_patience, dbgname, gen
-      ) into bendattrs, mutated;
-
       -- TODO: wm_combination
 
       if not mutated then
-        select * from wm_elimination(
-          bendattrs, dhalfcircle, dbgname, gen) into bendattrs, mutated;
+        select * from wm_elimination(bends, attrs,
+          dhalfcircle, dbgname, gen) into bends, mutated;
       end if;
 
       if mutated then
-        for j in 1..array_length(bendattrs, 1) loop
-          bends[j] = bendattrs[j].bend;
-        end loop;
         lines[i] = st_linemerge(st_union(bends));
+
         if st_geometrytype(lines[i]) != 'ST_LineString' then
           -- For manual debugging:
           --insert into wm_manual(name, way)
@@ -821,5 +752,4 @@ begin
   elseif l_type = 'ST_MultiLineString' then
     return st_union(lines);
   end if;
-end
-$$ language plpgsql;
+end $$ language plpgsql;