wm/wm.sql

489 lines
15 KiB
MySQL
Raw Normal View History

2021-05-19 22:57:45 +03:00
\set ON_ERROR_STOP on
SET plpgsql.extra_errors TO 'all';
2021-05-19 22:57:46 +03:00
-- detect_bends detects bends using the inflection angles. No corrections.
2021-05-19 22:57:45 +03:00
drop function if exists detect_bends;
2021-05-19 22:57:47 +03:00
create function detect_bends(
line geometry,
dbgname text default null,
dbgstagenum integer default null,
OUT bends geometry[]
) as $$
2021-05-19 22:57:45 +03:00
declare
2021-05-19 22:57:46 +03:00
pi constant real default radians(180);
2021-05-19 22:57:45 +03:00
p geometry;
p1 geometry;
p2 geometry;
p3 geometry;
2021-05-19 22:57:45 +03:00
bend geometry;
2021-05-19 22:57:45 +03:00
prev_sign int4;
cur_sign int4;
l_type text;
2021-05-19 22:57:45 +03:00
begin
2021-05-19 22:57:47 +03:00
l_type = st_geometrytype(line);
if l_type != 'ST_LineString' then
raise 'This function works with ST_LineString, got %', l_type;
end if;
2021-05-19 22:57:46 +03:00
-- The last vertex is iterated over twice, because the algorithm uses 3
2021-05-19 22:57:46 +03:00
-- vertices to calculate the angle between them.
2021-05-19 22:57:45 +03:00
--
-- Given 3 vertices p1, p2, p3:
--
2021-05-19 22:57:45 +03:00
-- p1___ ...
-- /
-- ... _____/
-- p3 p2
2021-05-19 22:57:45 +03:00
--
2021-05-19 22:57:45 +03:00
-- 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].
2021-05-19 22:57:45 +03:00
-- 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.
2021-05-19 22:57:47 +03:00
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
2021-05-19 22:57:45 +03:00
p3 = p2;
p2 = p1;
p1 = p;
2021-05-19 22:57:45 +03:00
continue when p3 is null;
2021-05-19 22:57:45 +03:00
cur_sign = sign(pi - st_angle(p1, p2, p2, p3));
2021-05-19 22:57:45 +03:00
2021-05-19 22:57:45 +03:00
if bend is null then
bend = st_makeline(p3, p2);
else
bend = st_linemerge(st_union(bend, st_makeline(p3, p2)));
end if;
2021-05-19 22:57:45 +03:00
if prev_sign + cur_sign = 0 then
if bend is not null then
2021-05-19 22:57:45 +03:00
bends = bends || bend;
2021-05-19 22:57:45 +03:00
end if;
bend = st_makeline(p3, p2);
end if;
prev_sign = cur_sign;
2021-05-19 22:57:45 +03:00
end loop;
2021-05-19 22:57:45 +03:00
2021-05-19 22:57:45 +03:00
-- 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
2021-05-19 22:57:45 +03:00
bends = bends || bend;
2021-05-19 22:57:45 +03:00
end if;
2021-05-19 22:57:47 +03:00
if dbgname is not null then
for i in 1..array_length(bends, 1) loop
insert into wm_debug(stage, name, gen, nbend, way) values(
'bbends',
dbgname,
dbgstagenum,
i,
bends[i]
);
insert into wm_debug(stage, name, gen, nbend, way) values(
'bbends-polygon',
dbgname,
dbgstagenum,
i,
st_makepolygon(st_addpoint(bends[i], st_startpoint(bends[i])))
);
end loop;
end if;
2021-05-19 22:57:45 +03:00
end
$$ language plpgsql;
2021-05-19 22:57:45 +03:00
2021-05-19 22:57:45 +03:00
-- fix_gentle_inflections moves bend endpoints following "Gentle Inflection at
-- End of a Bend" section.
2021-05-19 22:57:45 +03:00
--
-- 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
2021-05-19 22:57:45 +03:00
-- commulative inflection angle small (see variable below).
2021-05-19 22:57:45 +03:00
--
-- The implementation could be significantly optimized to avoid `st_reverse`
-- and array reversals, trading for complexity in fix_gentle_inflections1.
2021-05-19 22:57:47 +03:00
create or replace function fix_gentle_inflections(
INOUT bends geometry[],
dbgname text default null,
dbgstagenum integer default null
) as $$
2021-05-19 22:57:45 +03:00
declare
len int4;
2021-05-19 22:57:45 +03:00
bends1 geometry[];
begin
len = array_length(bends, 1);
2021-05-19 22:57:45 +03:00
bends = fix_gentle_inflections1(bends);
for i in 1..len loop
2021-05-19 22:57:45 +03:00
bends1[i] = st_reverse(bends[len-i+1]);
end loop;
2021-05-19 22:57:45 +03:00
bends1 = fix_gentle_inflections1(bends1);
for i in 1..len loop
2021-05-19 22:57:45 +03:00
bends[i] = st_reverse(bends1[len-i+1]);
end loop;
2021-05-19 22:57:47 +03:00
if dbgname is not null then
for i in 1..array_length(bends, 1) loop
insert into wm_debug(stage, name, gen, nbend, way) values(
'cinflections',
dbgname,
dbgstagenum,
i,
bends[i]
);
insert into wm_debug(stage, name, gen, nbend, way) values(
'cinflections-polygon',
dbgname,
dbgstagenum,
i,
st_makepolygon(st_addpoint(bends[i], st_startpoint(bends[i])))
);
end loop;
end if;
2021-05-19 22:57:45 +03:00
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.
2021-05-19 22:57:46 +03:00
drop function if exists fix_gentle_inflections1;
create function fix_gentle_inflections1(INOUT bends geometry[]) as $$
2021-05-19 22:57:45 +03:00
declare
2021-05-19 22:57:46 +03:00
pi constant real default radians(180);
2021-05-19 22:57:46 +03:00
-- the threshold when the angle is still "small", so gentle inflections can
-- be joined
2021-05-19 22:57:47 +03:00
small_angle constant real default radians(45);
2021-05-19 22:57:45 +03:00
ptail geometry; -- tail point of tail bend
phead geometry[]; -- 3 tail points of head bend
2021-05-19 22:57:45 +03:00
i int4; -- bends[i] is the current head
2021-05-19 22:57:45 +03:00
begin
2021-05-19 22:57:45 +03:00
for i in 2..array_length(bends, 1) loop
2021-05-19 22:57:45 +03:00
-- 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.
--
--
2021-05-19 22:57:45 +03:00
-- Assume this curve (figure `inflection-1`), going clockwise from A:
2021-05-19 22:57:45 +03:00
--
2021-05-19 22:57:45 +03:00
-- \______B
-- A `-------. C
-- |
-- G___ F |
-- / `-----.____+ D
-- E
2021-05-19 22:57:45 +03:00
--
-- After processing the curve following the definition of a bend, the bend
2021-05-19 22:57:45 +03:00
-- [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].
2021-05-19 22:57:45 +03:00
select geom from st_dumppoints(bends[i-1])
order by path[1] asc limit 1 into ptail;
2021-05-19 22:57:45 +03:00
while true loop
-- copy last 3 points of bends[i-1] (tail) to ptail
2021-05-19 22:57:45 +03:00
select array(
select geom from st_dumppoints(bends[i]) order by path[1] asc limit 3
) into phead;
2021-05-19 22:57:45 +03:00
2021-05-19 22:57:45 +03:00
-- if the bend got too short, stop processing it
exit when array_length(phead, 1) < 3;
2021-05-19 22:57:46 +03:00
-- inflection angle between ptail[1:3] is "large", stop processing
2021-05-19 22:57:45 +03:00
exit when abs(st_angle(phead[1], phead[2], phead[3]) - pi) > small_angle;
2021-05-19 22:57:45 +03:00
2021-05-19 22:57:46 +03:00
-- distance from head's 1st vertex should be larger than from 2nd vertex
2021-05-19 22:57:45 +03:00
exit when st_distance(ptail, phead[2]) < st_distance(ptail, phead[3]);
2021-05-19 22:57:45 +03:00
2021-05-19 22:57:46 +03:00
-- Detected a gentle inflection.
-- Move head of the tail to the tail of head
2021-05-19 22:57:45 +03:00
bends[i] = st_removepoint(bends[i], 0);
bends[i-1] = st_addpoint(bends[i-1], phead[3]);
2021-05-19 22:57:45 +03:00
end loop;
end loop;
2021-05-19 22:57:45 +03:00
end
$$ language plpgsql;
-- self_crossing eliminates self-crossing from the bends, following the
-- article's section "Self-line Crossing When Cutting a Bend".
2021-05-19 22:57:46 +03:00
drop function if exists self_crossing;
create function self_crossing(
INOUT bends geometry[],
OUT mutated boolean
) as $$
declare
2021-05-19 22:57:46 +03:00
pi constant real default radians(180);
2021-05-19 22:57:45 +03:00
i int4;
j int4;
2021-05-19 22:57:45 +03:00
prev_length int4;
a geometry;
b geometry;
2021-05-19 22:57:45 +03:00
multi geometry;
begin
mutated = false;
2021-05-19 22:57:45 +03:00
-- go through the bends and find one where sum of inflection angle is >180
2021-05-19 22:57:45 +03:00
for i in 1..array_length(bends, 1) loop
continue when abs(inflection_angle(bends[i])) <= pi;
2021-05-19 22:57:45 +03:00
-- 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
2021-05-19 22:57:45 +03:00
2021-05-19 22:57:45 +03:00
-- 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
-- of the vertex.
2021-05-19 22:57:45 +03:00
j = 0;
while j < array_length(bends, 1) loop
j = j + 1;
2021-05-19 22:57:45 +03:00
continue when i = j;
2021-05-19 22:57:45 +03:00
2021-05-19 22:57:45 +03:00
-- do end vertices of bend[i] cross bend[j]?
2021-05-19 22:57:45 +03:00
a = st_pointn(bends[i], 1);
b = st_pointn(bends[i], -1);
multi = st_split(bends[j], st_makeline(a, b));
2021-05-19 22:57:45 +03:00
continue when st_numgeometries(multi) = 1;
2021-05-19 22:57:46 +03:00
continue when st_numgeometries(multi) = 2 and
(st_contains(bends[j], a) or st_contains(bends[j], b));
2021-05-19 22:57:45 +03:00
-- stars are aligned, we are changing the bend
mutated = true;
2021-05-19 22:57:46 +03:00
-- 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.
2021-05-19 22:57:45 +03:00
prev_length = array_length(bends, 1);
if j < i then
2021-05-19 22:57:46 +03:00
-- 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);
2021-05-19 22:57:46 +03:00
bends[j] = st_setpoint(
2021-05-19 22:57:46 +03:00
bends[j],
st_npoints(bends[j])-1,
2021-05-19 22:57:46 +03:00
st_pointn(bends[i], st_npoints(bends[i]))
);
2021-05-19 22:57:45 +03:00
bends = bends[1:j] || bends[i+1:prev_length];
j = i;
else
2021-05-19 22:57:46 +03:00
-- 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)
);
2021-05-19 22:57:45 +03:00
bends = bends[1:i] || bends[j+1:prev_length];
end if;
2021-05-19 22:57:45 +03:00
j = j - prev_length + array_length(bends, 1);
2021-05-19 22:57:45 +03:00
end loop;
2021-05-19 22:57:46 +03:00
end loop;
end
$$ language plpgsql;
2021-05-19 22:57:45 +03:00
drop function if exists inflection_angle;
create function inflection_angle (IN bend geometry, OUT angle real) as $$
declare
2021-05-19 22:57:46 +03:00
pi constant real default radians(180);
p0 geometry;
p1 geometry;
p2 geometry;
p3 geometry;
begin
angle = 0;
2021-05-19 22:57:47 +03:00
for p0 in select geom from st_dumppoints(bend) 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;
end
$$ language plpgsql;
2021-05-19 22:57:46 +03:00
drop function if exists bend_attrs;
2021-05-19 22:57:47 +03:00
drop function if exists isolated_bends;
2021-05-19 22:57:46 +03:00
drop type if exists t_bend_attrs;
create type t_bend_attrs as (
bend geometry,
2021-05-19 22:57:47 +03:00
area real,
cmp real,
adjsize real,
2021-05-19 22:57:47 +03:00
baselinelength real,
2021-05-19 22:57:47 +03:00
curvature real,
isolated boolean
2021-05-19 22:57:46 +03:00
);
2021-05-19 22:57:46 +03:00
create function bend_attrs(bends geometry[], dbgname text default null) returns setof t_bend_attrs as $$
2021-05-19 22:57:46 +03:00
declare
2021-05-19 22:57:47 +03:00
fourpi constant real default 4*radians(180);
2021-05-19 22:57:46 +03:00
i int4;
2021-05-19 22:57:46 +03:00
polygon geometry;
2021-05-19 22:57:46 +03:00
bend geometry;
res t_bend_attrs;
2021-05-19 22:57:46 +03:00
begin
for i in 1..array_length(bends, 1) loop
2021-05-19 22:57:46 +03:00
bend = bends[i];
2021-05-19 22:57:47 +03:00
res = null;
res.bend = bend;
2021-05-19 22:57:47 +03:00
res.area = 0;
res.cmp = 0;
res.adjsize = 0;
2021-05-19 22:57:47 +03:00
res.baselinelength = st_distance(st_startpoint(bend), st_endpoint(bend));
2021-05-19 22:57:47 +03:00
res.curvature = inflection_angle(bend) / st_length(bend);
res.isolated = false;
2021-05-19 22:57:47 +03:00
if st_numpoints(bend) >= 3 then
polygon = st_makepolygon(st_addpoint(bend, st_startpoint(bend)));
2021-05-19 22:57:46 +03:00
-- 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 u: A = (u^2)/(4*pi)
-- 4. divide P by A: cmp = P/A = P/((u^2)*4*pi) = 4*pi*P/u^2
2021-05-19 22:57:47 +03:00
res.area = st_area(polygon);
res.cmp = fourpi*res.area/(st_perimeter(polygon)^2);
if res.cmp > 0 then
res.adjsize = (res.area*(0.75/res.cmp));
end if;
end if;
2021-05-19 22:57:47 +03:00
2021-05-19 22:57:46 +03:00
if dbgname is not null then
2021-05-19 22:57:47 +03:00
insert into wm_debug (stage, name, nbend, way, props) values(
2021-05-19 22:57:46 +03:00
'ebendattrs',
dbgname,
2021-05-19 22:57:46 +03:00
i,
2021-05-19 22:57:46 +03:00
bend,
2021-05-19 22:57:47 +03:00
jsonb_build_object(
2021-05-19 22:57:46 +03:00
'area', res.area,
'cmp', res.cmp,
'adjsize', res.adjsize,
2021-05-19 22:57:47 +03:00
'baselinelength', res.baselinelength,
2021-05-19 22:57:47 +03:00
'curvature', res.curvature
2021-05-19 22:57:46 +03:00
)
);
2021-05-19 22:57:46 +03:00
end if;
2021-05-19 22:57:46 +03:00
return next res;
2021-05-19 22:57:46 +03:00
end loop;
end;
$$ language plpgsql;
2021-05-19 22:57:47 +03:00
create function isolated_bends(INOUT bendattrs t_bend_attrs[], dbgname text default null) as $$
declare
2021-05-19 22:57:47 +03:00
isolation_threshold constant real default 0.25; -- if neighbor's curvatures are within, it's isolated
2021-05-19 22:57:47 +03:00
this real;
2021-05-19 22:57:47 +03:00
skip_next bool;
2021-05-19 22:57:47 +03:00
res t_bend_attrs;
i int4;
begin
2021-05-19 22:57:47 +03:00
for i in 2..array_length(bendattrs, 1)-1 loop
res = bendattrs[i];
if skip_next then
skip_next = false;
2021-05-19 22:57:47 +03:00
else
2021-05-19 22:57:47 +03:00
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;
end if;
2021-05-19 22:57:47 +03:00
end if;
2021-05-19 22:57:47 +03:00
if dbgname is not null then
insert into wm_debug (stage, name, nbend, way, props) values(
'fisolated_bends',
dbgname,
2021-05-19 22:57:47 +03:00
i,
2021-05-19 22:57:47 +03:00
res.bend,
jsonb_build_object(
'isolated', res.isolated
)
);
end if;
2021-05-19 22:57:47 +03:00
end loop;
end
$$ language plpgsql;
2021-05-19 22:57:46 +03:00
-- 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;
2021-05-19 22:57:46 +03:00
create function ST_SimplifyWM(geom geometry, dbgname text default null) returns geometry as $$
2021-05-19 22:57:46 +03:00
declare
stagenum integer;
2021-05-19 22:57:46 +03:00
i integer;
2021-05-19 22:57:46 +03:00
line geometry;
2021-05-19 22:57:46 +03:00
lines geometry[];
2021-05-19 22:57:46 +03:00
bends geometry[];
2021-05-19 22:57:47 +03:00
bend_attrs t_bend_attrs[];
2021-05-19 22:57:46 +03:00
mutated boolean;
2021-05-19 22:57:46 +03:00
l_type text;
begin
l_type = st_geometrytype(geom);
if l_type = 'ST_LineString' then
2021-05-19 22:57:46 +03:00
lines = array[geom];
2021-05-19 22:57:46 +03:00
elseif l_type = 'ST_MultiLineString' then
2021-05-19 22:57:46 +03:00
lines = array((select a.geom from st_dump(geom) a order by path[1] asc));
2021-05-19 22:57:46 +03:00
else
raise 'Unknown geometry type %', l_type;
end if;
2021-05-19 22:57:46 +03:00
for i in 1..array_length(lines, 1) loop
2021-05-19 22:57:46 +03:00
mutated = true;
stagenum = 1;
2021-05-19 22:57:46 +03:00
while mutated loop
2021-05-19 22:57:46 +03:00
if dbgname is not null then
2021-05-19 22:57:47 +03:00
insert into wm_debug (stage, name, gen, nbend, way) values(
2021-05-19 22:57:46 +03:00
'afigures',
dbgname,
stagenum,
2021-05-19 22:57:46 +03:00
i,
2021-05-19 22:57:46 +03:00
lines[i]
2021-05-19 22:57:46 +03:00
);
end if;
2021-05-19 22:57:47 +03:00
bends = detect_bends(lines[i], dbgname, stagenum);
2021-05-19 22:57:47 +03:00
bends = fix_gentle_inflections(bends, dbgname, stagenum);
2021-05-19 22:57:46 +03:00
2021-05-19 22:57:46 +03:00
select * from self_crossing(bends) into bends, mutated;
2021-05-19 22:57:46 +03:00
2021-05-19 22:57:46 +03:00
if dbgname is not null then
2021-05-19 22:57:47 +03:00
insert into wm_debug(stage, name, gen, nbend, way) values(
2021-05-19 22:57:46 +03:00
'dcrossings',
dbgname,
stagenum,
2021-05-19 22:57:46 +03:00
generate_subscripts(bends, 1),
2021-05-19 22:57:46 +03:00
unnest(bends)
);
end if;
2021-05-19 22:57:46 +03:00
if mutated then
lines[i] = st_linemerge(st_union(bends));
stagenum = stagenum + 1;
2021-05-19 22:57:46 +03:00
continue;
end if;
2021-05-19 22:57:46 +03:00
2021-05-19 22:57:46 +03:00
-- self-crossing mutations are done, calculate bend properties
2021-05-19 22:57:47 +03:00
bend_attrs = array((select bend_attrs(bends, dbgname)));
2021-05-19 22:57:47 +03:00
2021-05-19 22:57:47 +03:00
perform isolated_bends(bend_attrs, dbgname);
2021-05-19 22:57:46 +03:00
end loop;
2021-05-19 22:57:46 +03:00
2021-05-19 22:57:46 +03:00
end loop;
if l_type = 'ST_LineString' then
2021-05-19 22:57:46 +03:00
return st_linemerge(st_union(lines));
2021-05-19 22:57:46 +03:00
elseif l_type = 'ST_MultiLineString' then
2021-05-19 22:57:46 +03:00
return st_union(lines);
2021-05-19 22:57:46 +03:00
end if;
end
$$ language plpgsql;