wm/wm.sql

181 lines
5.6 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:45 +03:00
drop function if exists detect_bends;
2021-05-19 22:57:45 +03:00
-- detect_bends detects bends using the inflection angles. It does not do corrections.
2021-05-19 22:57:45 +03:00
create or replace function detect_bends(line geometry, OUT bends geometry[]) as $$
2021-05-19 22:57:45 +03:00
declare
pi real;
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;
2021-05-19 22:57:45 +03:00
begin
pi = radians(180);
2021-05-19 22:57:45 +03:00
-- the last vertex is iterated over twice, because the algorithm uses 3 vertices
-- to calculate the angle between them.
--
-- 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:45 +03:00
for p in (
2021-05-19 22:57:45 +03:00
(select geom from st_dumppoints(line) order by path[1] asc)
2021-05-19 22:57:45 +03:00
union all
2021-05-19 22:57:45 +03:00
(select geom from st_dumppoints(line) order by path[1] desc limit 1)
2021-05-19 22:57:45 +03:00
) loop
2021-05-19 22:57:45 +03:00
p3 = p2;
p2 = p1;
p1 = p;
2021-05-19 22:57:45 +03:00
if p3 is null then
continue;
end if;
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: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
create or replace function fix_gentle_inflections(INOUT bends geometry[]) as $$
2021-05-19 22:57:45 +03:00
declare
len int4;
2021-05-19 22:57:45 +03:00
bends1 geometry[];
bends2 geometry[];
2021-05-19 22:57:45 +03:00
begin
len = array_length(bends, 1);
2021-05-19 22:57:45 +03:00
bends1 = fix_gentle_inflections1(bends);
for i in 1..len loop
bends2[i] = st_reverse(bends1[len-i+1]);
end loop;
bends2 = fix_gentle_inflections1(bends2);
for i in 1..len loop
bends[i] = st_reverse(bends2[len-i+1]);
end loop;
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.
create or replace function fix_gentle_inflections1(INOUT bends geometry[]) as $$
2021-05-19 22:57:45 +03:00
declare
2021-05-19 22:57:45 +03:00
pi real;
2021-05-19 22:57:45 +03:00
small_angle real;
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
pi = radians(180);
2021-05-19 22:57:45 +03:00
-- the threshold when the angle is still "small", so gentle inflections can
-- be joined
small_angle := radians(30);
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:45 +03:00
-- if inflection angle between ptail[1:3] "large", stop processing this bend
if abs(st_angle(phead[1], phead[2], phead[3]) - pi) > small_angle then
--raise notice 'quitting % because angle between % % %', st_astext(ptail), st_astext(phead[1]), st_astext(phead[2]), st_astext(phead[3]);
exit;
end if;
2021-05-19 22:57:45 +03:00
2021-05-19 22:57:45 +03:00
-- distance from head's first vertex should be larger than from second vertex
if st_distance(ptail, phead[2]) < st_distance(ptail, phead[3]) then
--raise notice 'quitting % because distance', st_astext(ptail);
exit;
end if;
2021-05-19 22:57:45 +03:00
2021-05-19 22:57:45 +03:00
-- detected a gentle inflection. Move head of the tail to the tail of head
--raise notice 'fixing a gentle inflection of angle %', degrees(abs(st_angle(phead[1], phead[2], phead[3]) - pi));
2021-05-19 22:57:45 +03:00
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;
2021-05-19 22:57:45 +03:00
-- https://wiki.postgresql.org/wiki/Array_reverse
create or replace function array_reverse(anyarray) returns anyarray as $$
select array(
select $1[i]
from generate_series(
array_lower($1,1),
array_upper($1,1)
) as s(i)
order by i desc
);
$$ language 'sql' strict immutable;