Adaptive bend exaggeration

IV/test.sql

@@ -190,6 +190,7 @@ declare
 begin
   select way from wm_debug where name='fig3' and stage='bbends' and gen=1 and nbend=2 into fig3b2;
   size = wm_adjsize(fig3b2);
-  bend = wm_exaggerate_bend(fig3b2, size, 50.);
+  bend = wm_exaggerate_bend(fig3b2);
+  perform assert_equals('ST_LineString', st_geometrytype(bend));
   insert into wm_debug(stage, name, gen, nbend, way) values('manual', 'fig3', 1, 1, bend);
 end $$ language plpgsql;

IV/wm.sql

@@ -331,6 +331,7 @@ 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 function if exists wm_st_intersects_neighbors;
 drop type if exists wm_t_bend_attrs;
 create type wm_t_bend_attrs as (
   bend geometry,
@@ -401,11 +402,7 @@ end $$ language plpgsql;
 
 -- wm_exaggerate_bend exaggerates a given bend. Must be a simple linestring.
 drop function if exists wm_exaggerate_bend;
-create function wm_exaggerate_bend(
+create function wm_exaggerate_bend(INOUT bend geometry) as $$
-  INOUT bend geometry,
-  size float,
-  desired_size float
-) as $$
 declare
   scale constant float default 2; -- per-step scaling factor
   midpoint geometry; -- midpoint of the baseline
@@ -413,15 +410,11 @@ declare
   bendm geometry; -- bend with coefficients to prolong the lines
   points geometry[];
 begin
-  if size = 0 then
-    raise 'invalid input: zero-area bend';
-  end if;
   midpoint = st_lineinterpolatepoint(st_makeline(
       st_pointn(bend, 1),
       st_pointn(bend, -1)
     ), .5);
 
-  while size < desired_size loop
   splitbend = wm_st_split(bend, st_lineinterpolatepoint(bend, .5));
   -- Convert bend to LINESTRINGM, where M is the fraction by how
   -- much the point will be prolonged:
@@ -446,8 +439,6 @@ begin
     ));
 
   bend = st_setsrid(st_makeline(points), st_srid(bend));
-    size = wm_adjsize(bend);
-  end loop;
 end
 $$ language plpgsql;
 
@@ -474,8 +465,57 @@ begin
   if cmp > 0 then
     adjsize = (area*(0.75/cmp));
   end if;
-end
+end $$ language plpgsql;
-$$ language plpgsql;
+
+drop function if exists wm_st_intersects_neighbors;
+create function wm_st_intersects_neighbors(
+  bendattrs wm_t_bend_attrs[],
+  bend geometry,
+  i integer,
+  intersect_patience integer
+) returns boolean as $$
+declare
+  n integer;
+  neighbor geometry;
+begin
+  -- 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);
+  if n > 3 and st_intersects(bend,
+    st_removepoint(st_removepoint(bendattrs[i-1].bend, n-1), n-2)) then
+    return true;
+  end if;
+
+  n = st_npoints(bendattrs[i+1].bend);
+  if n > 3 and st_intersects(bend,
+    st_removepoint(st_removepoint(bendattrs[i+1].bend, 0), 0)) then
+    return true;
+  end if;
+
+  -- Go through all the neighbors and see if the enlarged bend intersects
+  -- with any of them. If yes -- abort enlargement.
+  for n in -intersect_patience+1..intersect_patience-1 loop
+    continue when i+n < 1 or n in (-1, 0, 1);
+    continue when i+n > array_length(bendattrs, 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.
+    neighbor = bendattrs[i+n].bend;
+    if st_npoints(neighbor) > 2 then
+      if n = -2 and st_npoints(bendattrs[i+n+1].bend) = 3 then
+        neighbor = st_removepoint(neighbor, st_npoints(neighbor)-1);
+      elsif n = 2 and st_npoints(bendattrs[i+n-1].bend) = 3 then
+        neighbor = st_removepoint(neighbor, 0);
+      end if;
+    end if;
+    if st_intersects(bend, neighbor) then
+      return true;
+    end if;
+  end loop;
+
+  return false;
+end $$ language plpgsql;
 
 -- wm_exaggeration is the Exaggeration Operator described in the WM paper.
 create function wm_exaggeration(
@@ -489,64 +529,42 @@ create function wm_exaggeration(
 declare
   desired_size constant float default pi()*(dhalfcircle^2)/8;
   tmpbendattr wm_t_bend_attrs;
-  tmpint geometry;
+  neighbor geometry;
+  size float;
   i integer;
   n integer;
-  last_id integer;
+  this_mutated boolean;
 begin
   mutated = false;
   <<bendloop>>
   for i in 1..array_length(bendattrs, 1) loop
-    if bendattrs[i].isolated and bendattrs[i].adjsize < desired_size then
+    continue when not bendattrs[i].isolated;
-      tmpbendattr.bend = wm_exaggerate_bend(
+    size = bendattrs[i].adjsize;
-        bendattrs[i].bend,
+    this_mutated = false;
-        bendattrs[i].adjsize,
+
-        desired_size
+    -- keep increasing this bend until either size permits, or it hits a
+    -- neighboring bend (intersect_patience number of neighbors will be
+    -- checked). When the size is right, stick it to the line.
+    while size < desired_size loop
+      tmpbendattr.bend = wm_exaggerate_bend(bendattrs[i].bend);
+      exit when wm_st_intersects_neighbors(
+        bendattrs,
+        tmpbendattr.bend,
+        i,
+        intersect_patience
       );
 
-      -- does tmpbendattrs.bend intersect with the previous or next
+      this_mutated = true;
-      -- 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);
-      if n > 3 then
-        continue when st_intersects(tmpbendattr.bend,
-          st_removepoint(st_removepoint(bendattrs[i-1].bend, n-1), n-2));
-      end if;
-
-      n = st_npoints(bendattrs[i+1].bend);
-      if n > 3 then
-        continue when st_intersects(tmpbendattr.bend,
-          st_removepoint(st_removepoint(bendattrs[i+1].bend, 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);
-
-        -- 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;
-        if st_npoints(tmpint) > 2 then
-          if n = -2 and st_npoints(bendattrs[i+n+1].bend) = 3 then
-            tmpint = st_removepoint(tmpint, st_npoints(tmpint)-1);
-          elsif n = 2 and st_npoints(bendattrs[i+n-1].bend) = 3 then
-            tmpint = st_removepoint(tmpint, 0);
-          end if;
-        end if;
-
-        continue bendloop when st_intersects(tmpbendattr.bend, tmpint);
-      end loop;
-
-      -- No intersections within intersect_patience, mutate bend!
-      mutated = true;
       bendattrs[i] = tmpbendattr;
+      size = wm_adjsize(bendattrs[i].bend);
+    end loop;
+    continue when not this_mutated;
 
-      -- remove last vertex of the previous bend and first vertex of the next
+    mutated = true;
+
+    -- Stick the right-sized bend to the line
+    -----------------------------------------
+    -- 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.
@@ -563,7 +581,6 @@ begin
         'gexaggeration', dbgname, dbggen, i, bendattrs[i].bend
       );
     end if;
-    end if;
   end loop;
 end $$ language plpgsql;
 
@@ -799,7 +816,8 @@ begin
         end loop;
         lines[i] = st_linemerge(st_union(bends));
         if st_geometrytype(lines[i]) != 'ST_LineString' then
-          -- For manual debugging:
+          -- For manual debugging, usually when wm_exaggeration returns
+          -- ST_MultiLineString. Uncomment the code below and `raise notice`.
           --insert into wm_manual(name, way)
           --select 'non-linestring-' || a.path[1], a.geom
           --from st_dump(lines[i]) a