diff --git a/IV/wm.sql b/IV/wm.sql index 58d9a38..0c962c6 100644 --- a/IV/wm.sql +++ b/IV/wm.sql @@ -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; <> - 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;