motiejus/wm
1
Fork 0
Code Releases Activity

split bendattrs to bends and attrs

Browse Source 
this makes code easier to reason about, and some lines shorter (thus
removing splitting requirement).
This commit is contained in:
Motiejus Jakštys 2021-05-19 22:57:50 +03:00 committed by Motiejus Jakštys
parent 5783dd3f3e
commit 2c7be0989a
1 changed files with 108 additions and 178 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

286
wm.sql
View File

@ -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;