commit 48a5dbdc0134ea6a208ca22e1a46b3e21fde9ee5 (tree)
parent a6ec2169acb8353e70d7f42dfbbe23fb50aabbf3
Author: Motiejus Jakštys <motiejus@uber.com>
Date: Fri, 7 May 2021 13:21:06 +0300
bring back working wm.sql
Diffstat:
| M | IV/wm.sql | | | 468 | +++++++++++++++++++++++++++++++++++++++++++------------------------------------ |
1 file changed, 257 insertions(+), 211 deletions(-)
diff --git a/IV/wm.sql b/IV/wm.sql
@@ -81,14 +81,16 @@ 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.
@@ -139,7 +141,8 @@ 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.
@@ -201,7 +204,8 @@ 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.
@@ -231,7 +235,9 @@ 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".
@@ -299,7 +305,8 @@ 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 $$
@@ -317,14 +324,17 @@ 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 function if exists wm_st_intersects_neighbors;
-drop type if exists wm_t_attrs;
-create type wm_t_attrs as (
+drop type if exists wm_t_bend_attrs;
+create type wm_t_bend_attrs as (
+ bend geometry,
+ area real,
adjsize real,
baselinelength real,
curvature real,
@@ -334,66 +344,40 @@ create function wm_bend_attrs(
bends geometry[],
dbgname text default null,
dbggen integer default null
-) returns wm_t_attrs[] as $$
+) returns setof wm_t_bend_attrs as $$
declare
- isolation_threshold constant real default 0.5;
- attrs wm_t_attrs[];
- attr wm_t_attrs;
- bend geometry;
+ cmp float;
i int4;
- needs_curvature real;
- skip_next boolean;
- dbglastid integer;
+ polygon geometry;
+ bend geometry;
+ res wm_t_bend_attrs;
begin
for i in 1..array_length(bends, 1) loop
bend = bends[i];
- 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;
+ 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;
if st_numpoints(bend) >= 3 then
- attr.adjsize = wm_adjsize(bend);
+ res.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', attrs[i].adjsize,
- 'baselinelength', attrs[i].baselinelength,
- 'curvature', attrs[i].curvature,
- 'isolated', false
+ 'adjsize', res.adjsize,
+ 'baselinelength', res.baselinelength,
+ 'curvature', res.curvature
)
- ) 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;
-
- return attrs;
-end $$ language plpgsql;
+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
@@ -417,44 +401,55 @@ 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(INOUT bend geometry) as $$
+create function wm_exaggerate_bend(
+ INOUT bend geometry,
+ size float,
+ desired_size float
+) as $$
declare
- scale constant float default 1.2; -- exaggeration enthusiasm
+ scale constant float default 2; -- per-step scaling factor
midpoint geometry; -- midpoint of the baseline
splitbend geometry; -- bend split across its half
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);
- 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:
- -- 1. draw a line between midpoint and the point on the bend.
- -- 2. multiply the line length by M. Midpoint stays intact.
- -- 3. the new set of lines form a new bend.
- -- Uses linear interpolation; can be updated to gaussian or similar;
- -- then interpolate manually instead of relying on st_addmeasure.
- bendm = st_collect(
- st_addmeasure(st_geometryn(splitbend, 1), 1, scale),
- st_addmeasure(st_geometryn(splitbend, 2), scale, 1)
- );
-
- points = array((
- select st_scale(
- st_makepoint(st_x(geom), st_y(geom)),
- st_makepoint(st_m(geom), st_m(geom)),
- midpoint
- )
- from st_dumppoints(bendm)
- order by path[1], path[2]
- ));
-
- bend = st_setsrid(st_makeline(points), st_srid(bend));
-end $$ language plpgsql;
+ 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:
+ -- 1. draw a line between midpoint and the point on the bend.
+ -- 2. multiply the line length by M. Midpoint stays intact.
+ -- 3. the new set of lines form a new bend.
+ -- Uses linear interpolation; can be updated to gaussian or similar;
+ -- then interpolate manually instead of relying on st_addmeasure.
+ bendm = st_collect(
+ st_addmeasure(st_geometryn(splitbend, 1), 1, scale),
+ st_addmeasure(st_geometryn(splitbend, 2), scale, 1)
+ );
+
+ points = array((
+ select st_scale(
+ st_makepoint(st_x(geom), st_y(geom)),
+ st_makepoint(st_m(geom), st_m(geom)),
+ midpoint
+ )
+ from st_dumppoints(bendm)
+ order by path[1], path[2]
+ ));
+
+ bend = st_setsrid(st_makeline(points), st_srid(bend));
+ size = wm_adjsize(bend);
+ end loop;
+end
+$$ language plpgsql;
-- wm_adjsize calculates adjusted size for a polygon. Can return 0.
drop function if exists wm_adjsize;
@@ -479,62 +474,12 @@ begin
if cmp > 0 then
adjsize = (area*(0.75/cmp));
end if;
-end $$ language plpgsql;
-
-drop function if exists wm_st_intersects_neighbors;
-create function wm_st_intersects_neighbors(
- bends geometry[],
- 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(bends[i-1]);
- if n > 3 and st_intersects(bend,
- st_removepoint(st_removepoint(bends[i-1], n-1), n-2)) then
- return true;
- end if;
-
- n = st_npoints(bends[i+1]);
- if n > 3 and st_intersects(bend,
- st_removepoint(st_removepoint(bends[i+1], 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(bends, 1);
-
- -- More special handling: if the neigbhoring bend has 3 vertices, the
- -- neighbor's neighbor may just touch the bend; in this
- -- case, the nearest vertex should be removed before comparing.
- neighbor = bends[i+n];
- if st_npoints(neighbor) > 2 then
- if n = -2 and st_npoints(bends[i+n+1]) = 3 then
- neighbor = st_removepoint(neighbor, st_npoints(neighbor)-1);
- elsif n = 2 and st_npoints(bends[i+n-1]) = 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;
+end
+$$ language plpgsql;
-- wm_exaggeration is the Exaggeration Operator described in the WM paper.
create function wm_exaggeration(
- INOUT bends geometry[],
- attrs wm_t_attrs[],
+ INOUT bendattrs wm_t_bend_attrs[],
dhalfcircle float,
intersect_patience integer,
dbgname text default null,
@@ -543,58 +488,87 @@ create function wm_exaggeration(
) as $$
declare
desired_size constant float default pi()*(dhalfcircle^2)/8;
- bend geometry;
- neighbor geometry;
- size float;
+ tmpbendattr wm_t_bend_attrs;
+ tmpint geometry;
i integer;
- this_mutated boolean;
+ n integer;
+ last_id integer;
begin
mutated = false;
- i = 0;
- while i < array_length(attrs, 1) loop
- i = i + 1;
+ <<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
+ );
- continue when not attrs[i].isolated;
- this_mutated = false;
+ -- does tmpbendattrs.bend intersect with the previous or next
+ -- intersect_patience bends? If they do, abort exaggeration for this one.
- -- 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.
- size = attrs[i].adjsize;
- while size < desired_size loop
- bend = wm_exaggerate_bend(bends[i]);
- exit when wm_st_intersects_neighbors(bends, bend, i, intersect_patience);
- this_mutated = true;
- bends[i] = bend;
- size = wm_adjsize(bend);
- end loop;
- continue when not this_mutated;
+ -- 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;
- 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.
- bend = st_removepoint(bends[i-1], st_npoints(bends[i-1])-1);
- bends[i-1] = bend;
- bend = st_removepoint(bends[i+1], 0);
- bends[i+1] = bend;
+ 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;
- if dbgname is not null then
- insert into wm_debug (stage, name, gen, nbend, way) values(
- 'gexaggeration', dbgname, dbggen, i, bends[i]);
- 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;
- -- Next bend was modified, so `isolated` is not valid for neighbor's
- -- neighbor. Skip over.
- i = i + 2;
+ continue bendloop when st_intersects(tmpbendattr.bend, tmpint);
+ end loop;
+
+ -- No intersections within intersect_patience, mutate bend!
+ mutated = true;
+ bendattrs[i] = tmpbendattr;
+
+ -- 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;
+
+ if dbgname is not null then
+ insert into wm_debug (stage, name, gen, nbend, way) values(
+ 'gexaggeration', dbgname, dbggen, i, bendattrs[i].bend
+ );
+ end if;
+ end if;
end loop;
end $$ language plpgsql;
create function wm_elimination(
- INOUT bends geometry[],
- attrs wm_t_attrs[],
+ INOUT bendattrs wm_t_bend_attrs[],
dhalfcircle float,
dbgname text default null,
dbggen integer default null,
@@ -605,54 +579,119 @@ 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(attrs, 1)-1 loop
+ while i < array_length(bendattrs, 1)-1 loop
i = i + 1;
- continue when attrs[i].adjsize = 0;
- continue when attrs[i].adjsize > desired_size;
+ continue when bendattrs[i].adjsize = 0;
+ continue when bendattrs[i].adjsize > desired_size;
if i = 2 then
- leftsize = attrs[i].adjsize + 1;
+ leftsize = bendattrs[i].adjsize + 1;
else
- leftsize = attrs[i-1].adjsize;
+ leftsize = bendattrs[i-1].adjsize;
end if;
- if i = array_length(attrs, 1)-1 then
- rightsize = attrs[i].adjsize + 1;
+ if i = array_length(bendattrs, 1)-1 then
+ rightsize = bendattrs[i].adjsize + 1;
else
- rightsize = attrs[i+1].adjsize;
+ rightsize = bendattrs[i+1].adjsize;
end if;
- continue when attrs[i].adjsize >= leftsize;
- continue when attrs[i].adjsize >= rightsize;
+ continue when bendattrs[i].adjsize >= leftsize;
+ continue when bendattrs[i].adjsize >= rightsize;
-- Local minimum. Elminate bend!
mutated = true;
- bends[i] = st_makeline(st_pointn(bends[i], 1), st_pointn(bends[i], -1));
-
+ tmpbendattr.bend = st_makeline(
+ st_pointn(bendattrs[i].bend, 1),
+ st_pointn(bendattrs[i].bend, -1)
+ );
+ bendattrs[i] = tmpbendattr;
-- remove last vertex of the previous bend and
-- first vertex of the next bend, because bends always
-- share a line segment together
- bends[i-1] = st_removepoint(bends[i-1], st_npoints(bends[i-1])-1);
- bends[i+1] = st_removepoint(bends[i+1], 0);
+ 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;
-- 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(bends, 1),
- unnest(bends)
+ generate_subscripts(dbgbends, 1),
+ unnest(dbgbends)
);
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(
@@ -678,7 +717,8 @@ 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,
@@ -702,7 +742,7 @@ declare
line geometry;
lines geometry[];
bends geometry[];
- attrs wm_t_attrs[];
+ bendattrs wm_t_bend_attrs[];
mutated boolean;
l_type text;
begin
@@ -733,37 +773,42 @@ begin
select * from wm_self_crossing(bends, dbgname, gen) into bends, mutated;
- 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;
+ if mutated then
+ lines[i] = st_linemerge(st_union(bends));
+ gen = gen + 1;
+ continue;
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(
- bends, attrs, dhalfcircle, dbgname, gen) into bends, mutated;
+ bendattrs, dhalfcircle, dbgname, gen) into bendattrs, 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
- raise notice 'Got % (in %) instead of ST_LineString. '
+ -- For manual debugging:
+ --insert into wm_manual(name, way)
+ --select 'non-linestring-' || a.path[1], a.geom
+ --from st_dump(lines[i]) a
+ --order by a.path[1];
+ raise 'Got % (in %) instead of ST_LineString. '
'Does the exaggerated bend intersect with the line? '
'If so, try increasing intersect_patience.',
st_geometrytype(lines[i]), dbgname;
- -- For manual debugging, usually when wm_exaggeration returns
- -- ST_MultiLineString. Uncomment the code below and change `raise`
- -- to `raise notice` above.
- insert into wm_manual(name, way)
- select 'non-linestring-' || a.path[1], a.geom
- from st_dump(lines[i]) a
- order by a.path[1];
- exit lineloop;
+ --exit lineloop;
end if;
gen = gen + 1;
continue;
@@ -776,4 +821,5 @@ begin
elseif l_type = 'ST_MultiLineString' then
return st_union(lines);
end if;
-end $$ language plpgsql;
+end
+$$ language plpgsql;
