Move to github.com/motiejus/wm

main
Motiejus Jakštys 2021-05-19 23:15:44 +03:00 committed by Motiejus Jakštys
parent 72d025c1f8
commit 6d9c9d9267
27 changed files with 1 additions and 4415 deletions

View File

@ -1,13 +0,0 @@
FROM debian:bullseye
RUN apt-get update && \
DEBIAN_FRONTEND=noninteractive apt-get install -y \
texlive-full poppler-utils \
osm2pgsql postgresql-client \
python3-pygments python3-geopandas \
latexmk make \
docker.io
COPY layer2img.py /tmp/layer2img.py
RUN python3 /tmp/layer2img.py -o /tmp/foo.pdf && \
rm /tmp/layer2img.py /tmp/foo.pdf

View File

@ -1,386 +0,0 @@
OSM ?= lithuania-latest.osm.pbf
RIVERFILTER = Visinčia|Šalčia|Nemunas
SLIDES = slides-2021-03-29.pdf
GDB10LT ?= $(wildcard GDB10LT-static-*.zip)
# Max figure size (in meters) is when it's width is TEXTWIDTH_CM on scale 1:25k
SCALEDWIDTH = $(shell awk '/^TEXTWIDTH_CM/{print 25000/100*$$3}' layer2img.py)
##############################################################################
# These variables have to come before first use due to how macros are expanded
##############################################################################
NON_ARCHIVABLES = notes.txt referatui.txt slides-2021-03-29.txt
ARCHIVABLES = $(filter-out $(NON_ARCHIVABLES),$(shell git ls-files .))
LISTINGS = aggregate-rivers.sql wm.sql extract-and-generate
FIGURES = \
test-figures \
fig8-definition-of-a-bend \
fig8-elimination-gen1 \
fig8-elimination-gen2 \
fig8-elimination-gen3 \
fig5-gentle-inflection-before \
fig5-gentle-inflection-after \
inflection-1-gentle-inflection-before \
inflection-1-gentle-inflection-after \
fig6-selfcrossing \
selfcrossing-1 \
isolated-1-exaggerated
RIVERS = \
salvis-25k \
salvis-2x50k \
salvis-250k-10x \
salvis-grpk250-2x \
salvis-dp64-2x50k \
salvis-vw64-2x50k \
salvis-dpchaikin64-2x50k \
salvis-vwchaikin64-2x50k \
salvis-overlaid-dpchaikin64-2x50k \
salvis-overlaid-vwchaikin64-2x50k \
salvis-wm220-10x \
salvis-wm220-2x \
salvis-wm-overlaid-250k-zoom \
salvis-wm220
################################################################################
# FIGURES
################################################################################
test-figures_1SELECT = wm_figures
fig8-definition-of-a-bend_1SELECT = wm_debug where name='fig8' AND stage='afigures' AND gen=1
fig8-definition-of-a-bend_2SELECT = wm_debug where name='fig8' AND stage='bbends-polygon' AND gen=1
fig8-definition-of-a-bend_3SELECT = wm_debug where name='fig8' AND stage='bbends-polygon' AND gen=1
fig8-definition-of-a-bend_3LINESTYLE = dotted
fig8-elimination-gen1_1SELECT = wm_debug where name='fig8' AND stage='afigures' AND gen=1
fig8-elimination-gen1_2SELECT = wm_debug where name='fig8' AND stage='bbends-polygon' AND gen=1
fig8-elimination-gen1_3SELECT = wm_debug where name='fig8' AND stage='bbends-polygon' AND gen=1
fig8-elimination-gen1_3LINESTYLE = dotted
fig8-elimination-gen2_1SELECT = wm_debug where name='fig8' AND stage='afigures' AND gen=2
fig8-elimination-gen2_2SELECT = wm_debug where name='fig8' AND stage='bbends-polygon' AND gen=2
fig8-elimination-gen2_3SELECT = wm_debug where name='fig8' AND stage='bbends-polygon' AND gen=2
fig8-elimination-gen2_3LINESTYLE = dotted
fig8-elimination-gen3_1SELECT = wm_debug where name='fig8' AND stage='bbends' AND gen=3
fig8-elimination-gen3_2SELECT = wm_debug where name='fig8' AND stage='bbends-polygon' AND gen=3
fig8-elimination-gen3_3SELECT = wm_debug where name='fig8' AND stage='bbends-polygon' AND gen=3
fig8-elimination-gen3_3LINESTYLE = dotted
fig5-gentle-inflection-before_WITHDIV = 2
fig5-gentle-inflection-before_1SELECT = wm_debug where name='fig5' AND stage='afigures' AND gen=1
fig5-gentle-inflection-before_2SELECT = wm_debug where name='fig5' AND stage='bbends-polygon' AND gen=1
fig5-gentle-inflection-before_3SELECT = wm_debug where name='fig5' AND stage='bbends-polygon' AND gen=1
fig5-gentle-inflection-before_3LINESTYLE = dotted
fig5-gentle-inflection-after_WITHDIV = 2
fig5-gentle-inflection-after_1SELECT = wm_debug where name='fig5' AND stage='cinflections' AND gen=1
fig5-gentle-inflection-after_2SELECT = wm_debug where name='fig5' AND stage='cinflections-polygon' AND gen=1
fig5-gentle-inflection-after_3SELECT = wm_debug where name='fig5' AND stage='cinflections-polygon' AND gen=1
fig5-gentle-inflection-after_3LINESTYLE = dotted
inflection-1-gentle-inflection-before_WIDTHDIV = 2
inflection-1-gentle-inflection-before_1SELECT = wm_debug where name='inflection-1' AND stage='afigures' AND gen=1
inflection-1-gentle-inflection-before_2SELECT = wm_debug where name='inflection-1' AND stage='bbends-polygon' AND gen=1
inflection-1-gentle-inflection-before_3SELECT = wm_debug where name='inflection-1' AND stage='bbends-polygon' AND gen=1
inflection-1-gentle-inflection-before_3LINESTYLE = dotted
inflection-1-gentle-inflection-after_WIDTHDIV = 2
inflection-1-gentle-inflection-after_1SELECT = wm_debug where name='inflection-1' AND stage='cinflections' AND gen=1
inflection-1-gentle-inflection-after_2SELECT = wm_debug where name='inflection-1' AND stage='cinflections-polygon' AND gen=1
inflection-1-gentle-inflection-after_3SELECT = wm_debug where name='inflection-1' AND stage='cinflections-polygon' AND gen=1
inflection-1-gentle-inflection-after_3LINESTYLE = dotted
fig6-selfcrossing_WIDTHDIV = 2
fig6-selfcrossing_1SELECT = wm_debug where name='fig6' AND stage='afigures' AND gen=1
fig6-selfcrossing_1LINESTYLE = dotted
fig6-selfcrossing_2SELECT = wm_debug where name='fig6' AND stage='dcrossings' AND gen=1
fig6-selfcrossing_3SELECT = wm_visuals where name='fig6-baseline'
fig6-selfcrossing_3COLOR = orange
selfcrossing-1_WIDTHDIV = 2
selfcrossing-1_1SELECT = wm_debug where name='selfcrossing-1' AND stage='afigures' AND gen=1
selfcrossing-1_1LINESTYLE = dotted
selfcrossing-1_2SELECT = wm_debug where name='selfcrossing-1' AND stage='dcrossings' AND gen=1
selfcrossing-1_3SELECT = wm_visuals where name='selfcrossing-1-baseline'
selfcrossing-1_3COLOR = orange
isolated-1-exaggerated_WIDTHDIV = 2
isolated-1-exaggerated_1SELECT = wm_debug where name='isolated-1' AND stage='afigures' AND gen=2
isolated-1-exaggerated_2SELECT = wm_debug where name='isolated-1' AND stage='afigures' AND gen=1
isolated-1-exaggerated_1COLOR = orange
################################################################################
# 250K
################################################################################
salvis-wm220-250k-2x_1SELECT = wm_visuals where name='salvis-wm220'
salvis-wm220-250k-2x_WIDTHDIV = 2
salvis-wm220-250k-10x_1SELECT = wm_visuals where name='salvis-wm220'
salvis-wm220-250k-10x_WIDTHDIV = 10
salvis-250k-10x_1SELECT = wm_visuals where name='salvis-grpk10'
salvis-250k-10x_WIDTHDIV = 10
salvis-wm-overlaid-250k-zoom_1SELECT = wm_visuals where name='salvis-wm220'
salvis-wm-overlaid-250k-zoom_2SELECT = wm_visuals where name='salvis-grpk10'
salvis-wm-overlaid-250k-zoom_1COLOR = orange
salvis-grpk250-2x_1SELECT = wm_visuals where name='salvis-grpk250'
salvis-grpk250-2x_WIDTHDIV = 2
################################################################################
# 50K
################################################################################
label_wm75 = Wang--Müller 1:\numprint{50000}
label_wm220 = Wang--Müller 1:\numprint{250000}
label_vw64 = Visvalingam--Whyatt
label_dp64 = Douglas \& Peucker
label_grpk10 = GRPK 1:\numprint{10000}
label_grpk50 = GRPK 1:\numprint{50000}
label_vwchaikin64 = $(label_vw64) and Chaikin
label_dpchaikin64 = $(label_dp64) and Chaikin
legend_ = lower left
legend_tr = lower right
legend_tl = lower center
define wm_vwdp50k
RIVERS += salvis-$(1)-$(2)-$(3)-$(4)x50k$(5)
salvis-$(1)-$(2)-$(3)-$(4)x50k$(5)_1SELECT = wm_visuals where name='salvis-$(1)'
salvis-$(1)-$(2)-$(3)-$(4)x50k$(5)_1COLOR = orange
salvis-$(1)-$(2)-$(3)-$(4)x50k$(5)_1LABEL = $(label_$(1))
$(if $(2),
salvis-$(1)-$(2)-$(3)-$(4)x50k$(5)_2SELECT = wm_visuals where name='salvis-$(2)'
salvis-$(1)-$(2)-$(3)-$(4)x50k$(5)_2COLOR = green
salvis-$(1)-$(2)-$(3)-$(4)x50k$(5)_2LABEL = $(label_$(2))
,)
$(if $(3),
salvis-$(1)-$(2)-$(3)-$(4)x50k$(5)_3SELECT = wm_visuals where name='salvis-$(3)'
salvis-$(1)-$(2)-$(3)-$(4)x50k$(5)_3LINESTYLE = $(6)
salvis-$(1)-$(2)-$(3)-$(4)x50k$(5)_3LABEL = $(label_$(3))
,)
salvis-$(1)-$(2)-$(3)-$(4)x50k$(5)_WIDTHDIV = $(4)
salvis-$(1)-$(2)-$(3)-$(4)x50k$(5)_QUADRANT = $(5)
salvis-$(1)-$(2)-$(3)-$(4)x50k$(5)_LEGEND = $(legend_$(5))
endef
wm_vwdp50kblack = $(call wm_vwdp50k,$(1),$(2),$(3),$(4),$(5))
wm_vwdp50kdotted = $(call wm_vwdp50k,$(1),$(2),$(3),$(4),$(5),dotted)
$(foreach x,vw64 dp64 vwchaikin64 dpchaikin64,\
$(eval $(call wm_vwdp50kdotted,wm75,$(x),grpk10,1,)) \
)
$(eval $(call wm_vwdp50kblack,wm75,grpk50,grpk10,1))
$(eval $(call wm_vwdp50kblack,wm75,grpk50,grpk10,1,tr))
$(eval $(call wm_vwdp50kblack,wm75,grpk50,grpk10,1,tl))
$(eval $(call wm_vwdp50kblack,wm75,,grpk10,1))
$(eval $(call wm_vwdp50kblack,wm75,,grpk10,1,tr))
$(eval $(call wm_vwdp50kblack,wm75,,grpk10,1,tl))
salvis-25k_1SELECT = wm_visuals where name='salvis-grpk10'
salvis-25k_WIDTHDIV = 1
salvis-2x50k_1SELECT = wm_visuals where name='salvis-grpk10'
salvis-2x50k_WIDTHDIV = 2
salvis-dp64-2x50k_1SELECT = wm_visuals where name='salvis-dp64'
salvis-dp64-2x50k_WIDTHDIV = 2
salvis-vw64-2x50k_1SELECT = wm_visuals where name='salvis-vw64'
salvis-vw64-2x50k_WIDTHDIV = 2
salvis-dpchaikin64-2x50k_2SELECT = wm_visuals where name='salvis-dpchaikin64'
salvis-dpchaikin64-2x50k_WIDTHDIV = 2
salvis-vwchaikin64-2x50k_2SELECT = wm_visuals where name='salvis-vwchaikin64'
salvis-vwchaikin64-2x50k_WIDTHDIV = 2
salvis-overlaid-dpchaikin64-2x50k_1SELECT = wm_visuals where name='salvis-dpchaikin64'
salvis-overlaid-dpchaikin64-2x50k_2SELECT = wm_visuals where name='salvis-grpk10'
salvis-overlaid-dpchaikin64-2x50k_1COLOR = orange
salvis-overlaid-dpchaikin64-2x50k_WIDTHDIV = 2
salvis-overlaid-dpchaikin64-2x50k_QUADRANT = tl
salvis-overlaid-vwchaikin64-2x50k_1SELECT = wm_visuals where name='salvis-vwchaikin64'
salvis-overlaid-vwchaikin64-2x50k_2SELECT = wm_visuals where name='salvis-grpk10'
salvis-overlaid-vwchaikin64-2x50k_1COLOR = orange
salvis-overlaid-vwchaikin64-2x50k_WIDTHDIV = 2
salvis-overlaid-vwchaikin64-2x50k_QUADRANT = tl
salvis-wm220_1SELECT = wm_visuals where name='salvis-wm220'
salvis-wm220_WIDTHDIV = 2
define FIG_template
$(1).pdf: layer2img.py Makefile $(2)
python3 ./layer2img.py --outfile=$(1).pdf \
$$(if $$($(1)_LEGEND),--legend="$$($(1)_LEGEND)") \
$$(if $$($(1)_WIDTHDIV),--widthdiv=$$($(1)_WIDTHDIV)) \
$$(if $$($(1)_QUADRANT),--quadrant=$$($(1)_QUADRANT)) \
$$(foreach i,1 2 3, \
$$(if $$($(1)_$$(i)LABEL),--g$$(i)-label="$$($(1)_$$(i)LABEL)") \
$$(if $$($(1)_$$(i)COLOR),--g$$(i)-color="$$($(1)_$$(i)COLOR)") \
$$(if $$($(1)_$$(i)SELECT),--g$$(i)-select="$$($(1)_$$(i)SELECT)") \
$$(if $$($(1)_$$(i)LINESTYLE),--g$$(i)-linestyle="$$($(1)_$$(i)LINESTYLE)") \
)
endef
$(foreach fig,$(FIGURES),$(eval $(call FIG_template,$(fig),.faux_test)))
$(foreach fig,$(RIVERS), $(eval $(call FIG_template,$(fig),.faux_visuals)))
#################################
# The thesis, publishable version
#################################
mj-msc-full.pdf: mj-msc.pdf version.inc.tex $(ARCHIVABLES) ## Thesis for publishing
cp $< .tmp-$@
for f in $^; do \
if [ "$$f" = "$<" ]; then continue; fi; \
pdfattach .tmp-$@ $$f .tmp2-$@; \
mv .tmp2-$@ .tmp-$@; \
done
mv .tmp-$@ $@
###############################
# Auxiliary targets for humans
###############################
.PHONY: test
test: .faux_test ## Unit tests (fast)
.PHONY: visuals
visuals: .faux_visuals # Generate visuals for paper (fast)
.PHONY: test-rivers
test-rivers: .faux_test-rivers ## Rivers tests (slow)
.PHONY: slides
slides: $(SLIDES)
.PHONY: refresh-rivers
refresh-rivers: refresh-rivers-10.sql refresh-rivers-50.sql refresh-rivers-250.sql ## Refresh river data from national datasets
###########################
# The report, quick version
###########################
mj-msc.pdf: mj-msc.tex version.inc.tex vars.inc.tex bib.bib \
$(LISTINGS) $(addsuffix .pdf,$(FIGURES)) $(addsuffix .pdf,$(RIVERS))
latexmk -shell-escape -pdf $<
############################
# Report's test dependencies
############################
.PHONY: allfigs
allfigs: $(addsuffix .pdf,$(FIGURES)) $(addsuffix .pdf,$(RIVERS))
.faux_db_pre: db init.sql
bash db start
bash db -f init.sql
touch $@
.faux_db: rivers-10.sql rivers-50.sql rivers-250.sql
bash db $(addprefix -f ,$^)
touch $@
.faux_db: .EXTRA_PREREQS = .faux_db_pre
.faux_test: test.sql wm.sql .faux_db
bash db -f $<
touch $@
.faux_visuals: visuals.sql .faux_test
bash db -v scaledwidth=$(SCALEDWIDTH) -f $<
touch $@
.faux_test-rivers: test-rivers.sql wm.sql Makefile .faux_db
bash db -f $<
touch $@
################################
# Report's non-test dependencies
################################
REF = $(shell git describe --abbrev=12 --always --dirty)
version.inc.tex: Makefile $(shell git rev-parse --git-dir 2>/dev/null)
TZ=UTC date '+\gdef\VCDescribe{%F (revision $(REF))}%' > $@
vars.inc.tex: vars.awk wm.sql Makefile
awk -f $< wm.sql
###############
# Misc commands
###############
slides-2021-03-29.pdf: slides-2021-03-29.txt
pandoc -t beamer -i $< -o $@
dump-debug_wm.sql.xz:
docker exec -ti wm-mj pg_dump -Uosm osm -t wm_devug | xz -v > $@
release.zip: mj-msc.tex mj-msc.bbl version.inc.tex vars.inc.tex \
$(addsuffix .pdf,$(FIGURES)) $(addsuffix .pdf,$(RIVERS)) \
$(shell git ls-files .)
-rm $@
mkdir -p .tmp; touch .tmp/editorial-version
zip $@ $^
zip $@ -j .tmp/editorial-version
mj-msc.bbl: mj-msc.tex bib.bib
biber mj-msc
mj-msc-gray.pdf: mj-msc.pdf
gs \
-sOutputFile=$@ \
-sDEVICE=pdfwrite \
-sColorConversionStrategy=Gray \
-dProcessColorModel=/DeviceGray \
-dCompatibilityLevel=1.4 \
-dNOPAUSE \
-dBATCH \
$<
.PHONY: clean
clean: ## Clean the current working directory
-bash db stop
-rm -r .faux_test .faux_aggregate-rivers .faux_test-rivers .faux_visuals \
.faux_db .faux_db_pre version.inc.tex vars.inc.tex version.aux \
version.fdb_latexmk _minted-mj-msc .tmp \
$(shell git ls-files -o mj-msc*) \
$(addsuffix .pdf,$(FIGURES)) \
$(addsuffix .pdf,$(RIVERS)) \
$(SLIDES)
.PHONY: clean-tables
clean-tables: ## Remove tables created during unit or rivers tests
bash db -c '\dt wm_*' | awk '/_/{print "drop table "$$3";"}' | bash db -f -
-rm .faux_db
.PHONY: help
help: ## Print this help message
@awk -F':.*?## ' '/^[a-z0-9.-]*: *.*## */{printf "%-18s %s\n",$$1,$$2}' \
$(MAKEFILE_LIST)
.PHONY: wc
wc: mj-msc.pdf ## Character and page count
@pdftotext $< - | \
awk '/\yReferences\y/{exit}; {print}' | \
tr -d '[:space:]' | wc -c | \
awk '{printf("Chars: %d, pages: %.1f\n", $$1, $$1/1500)}'
define refresh_rivers_template
.PHONY: refresh-$(1)
refresh-$(1): aggregate-rivers.sql gdr2pgsql .faux_db_pre
@if [ ! -f "$$($(2))" ]; then \
echo "ERROR: $(2)-static-*.zip not found. Run env $(2)=<...>"; \
exit 1; \
fi
./gdr2pgsql "$$($(2))" "$(3)" "$(RIVERFILTER)" "$(1)"
endef
$(eval $(call rivers_template,rivers-10.sql,GDB10LT,wm_rivers))
$(eval $(call rivers_template,rivers-50.sql,grpk50LT,wm_rivers_50))
$(eval $(call rivers_template,rivers-250.sql,grpk250LT,wm_rivers_250))

View File

@ -1,80 +1,4 @@
WangMüller algorithm in PostGIS
--------------------------------
This is a work-in-progress implementation following "Line generalization based
on analysis of shape characteristics" by Wang and Müller, 1998.
Structure
---------
There will be 2 deliverables from this folder:
- `wm.sql`, the implementation.
- paper `mj-msc-full.pdf`, a MSc thesis, explaining it.
It contains a few supporting files, notably:
- `tests.sql` synthetic unit tests.
- `test-rivers.sql` tests with real rivers.
- `Makefile` glues everything together.
- `layer2img.py` converts a PostGIS layer to an embeddable image.
- `aggregate-rivers.sql` combines multiple river objects (linestrings or
multilinestrings) to a single one.
- `init.sql` initializes PostGIS database for running the tests.
- `rivers-*.sql` are national dataset snapshots of rivers (`Makefile`
contains code to update them).
- ... and a few more files necessary to build the paper.
Running
-------
`make help` lists the select commands for humans. As of writing:
```
# make help
mj-msc-full.pdf Thesis for publishing
test Unit tests (fast)
test-rivers Rivers tests (slow)
clean Clean the current working directory
clean-tables Remove tables created during unit or rivers tests
help Print this help message
wc Character and page count
refresh-rivers Refresh rivers-*.sql from Open Street Maps
```
To execute the algorithm, run:
- `make test` for tests with synthetic data.
- `make test-rivers` for tests with real rivers. You may adjust the rivers and
data source (e.g. use a different country instead of Lithuania) by changing
the `Makefile` and the test files. Left as an exercise for the reader.
Building the paper (pdf)
------------------------
```
# make -j$(nproc) mj-msc-full.pdf
```
`mj-msc.tex` results in `mj-msc-full.pdf`, which will be at some point
published. It needs quite a few dependencies, including a functioning Docker
environment, postgresql client, geopandas, pygments, osm2pgsql, poppler, and a
"quite extensive" LaTeX installation. This was tested on Debian 11.
Contributing
------------
Please reach out to me before contributing. As of writing, this is not ready to
accept broader contributions, TODO:
- [x] Elimination operator.
- [x] Exaggeration operator.
- [ ] Combination operator.
- [ ] CI (unlikely to happen).
- [x] Known bug in `wm_self_crossing`: the program crashes with a river in
Lithuania.
License
-------
GPL 2.0 or later, same as PostGIS.
Moved to [github.com/motiejus/wm](https://github.com/motiejus/wm).

View File

@ -1,49 +0,0 @@
/* Aggregates rivers by name and proximity. */
drop function if exists aggregate_rivers;
create function aggregate_rivers() returns table(
id integer,
name text,
way geometry
) as $$
declare
c record;
cc record;
changed boolean;
begin
while (select count(1) from wm_rivers_tmp) > 0 loop
select * from wm_rivers_tmp limit 1 into c;
delete from wm_rivers_tmp a where a.id = c.id;
changed = true;
while changed loop
changed = false;
for cc in (
select * from wm_rivers_tmp a where
a.name = c.name and
st_dwithin(a.way, c.way, 500)
) loop
c.way = st_linemerge(st_union(c.way, cc.way));
delete from wm_rivers_tmp a where a.id = cc.id;
changed = true;
end loop;
end loop; -- while changed
return query select c.id, c.name, c.way;
end loop; -- count(1) from wm_rivers_tmp > 0
return;
end
$$ language plpgsql;
drop index if exists wm_rivers_tmp_id;
drop index if exists wm_rivers_tmp_gix;
drop table if exists wm_rivers_tmp;
create temporary table wm_rivers_tmp (id serial, name text, way geometry);
create index wm_rivers_tmp_id on wm_rivers_tmp(id);
create index wm_rivers_tmp_gix on wm_rivers_tmp using gist(way) include(name);
insert into wm_rivers_tmp (name, way)
select p.vardas as name, p.shape as way from :srctable p;
drop table if exists :dsttable;
create table :dsttable as (
select * from aggregate_rivers() where st_length(way) >= 50000
);
drop table wm_rivers_tmp;

Binary file not shown.

Before

Width:  |  Height:  |  Size: 13 KiB

View File

@ -1,262 +0,0 @@
@article{wang1998line,
title={Line generalization based on analysis of shape characteristics},
author={Wang, Zeshen and M{\"u}ller, Jean-Claude},
journal={Cartography and Geographic Information Systems},
volume={25},
number={1},
pages={3--15},
year={1998},
publisher={Taylor \& Francis}
}
@article{Kolanowski_2018,
title={Cartographic Line Generalization Based on Radius of Curvature Analysis},
volume={7},
ISSN={2220-9964},
url={http://dx.doi.org/10.3390/ijgi7120477},
DOI={10.3390/ijgi7120477},
number={12},
journal={ISPRS International Journal of Geo-Information},
publisher={MDPI AG},
author={Kolanowski, Bogdan and Augustyniak, Jacek and Latos, Dorota},
year={2018},
month={12},
pages={477}
}
@article{visvalingam1993line,
title={Line generalisation by repeated elimination of points},
author={Visvalingam, Maheswari and Whyatt, James D},
journal={The cartographic journal},
volume={30},
number={1},
pages={46--51},
year={1993},
publisher={Taylor \& Francis}
}
@article{muller1991generalization,
title={Generalization of spatial databases},
author={Muller, Jean-Claude},
journal={Geographical information systems},
volume={1},
pages={457--475},
year={1991},
publisher={John Wiley and Sons}
}
@article{miuller1995generalization,
title={Generalization-state of the art and issues},
author={Miuller, JC and Weibel, R and Lagrange, J and {\"E}alge, F},
journal={GIS and Generalisation: Methodology and Practice},
pages={3--17},
year={1995}
}
@inproceedings{mcmaster1992generalization,
title={Generalization in digital cartography},
author={McMaster, Robert Brainerd and Shea, K Stuart},
year={1992},
organization={Association of American Geographers Washington, DC}
}
@article{douglas1973algorithms,
title={Algorithms for the reduction of the number of points required to represent a digitized line or its caricature},
author={Douglas, David H and Peucker, Thomas K},
journal={Cartographica: the international journal for geographic information and geovisualization},
volume={10},
number={2},
pages={112--122},
year={1973},
publisher={University of Toronto Press}
}
% Algorithms for generalization, not reaching satisfactory results
@inproceedings{monmonier1986toward,
title={Toward a practicable model of cartographic generalisation.},
author={Monmonier, Mark},
booktitle={Auto Carto London. Proc. conference, 1986. Vol. 2},
pages={257--266},
year={1986},
organization={distributed Royal Institution of Chartered Surveyors}
}
@inproceedings{mcmaster1993spatial,
title={A spatial-object level organization of transformations for cartographic generalization},
author={McMaster, RB and Barnett, Leone},
booktitle={AUTOCARTO-CONFERENCE-},
pages={386--386},
year={1993},
organization={Citeseer}
}
@inproceedings{jiang2003line,
title={Line simplification using self-organizing maps},
author={Jiang, Bin and Nakos, Byron},
booktitle={Proceedings of the ISPRS Workshop on Spatial Analysis and Decision Making, Hong Kong, China},
pages={3--5},
year={2003}
}
@article{dyken2009simultaneous,
title={Simultaneous curve simplification},
author={Dyken, Christopher and D{\ae}hlen, Morten and Sevaldrud, Thomas},
journal={Journal of geographical systems},
volume={11},
number={3},
pages={273--289},
year={2009},
publisher={Springer}
}
@article{mustafa2006dynamic,
title={Dynamic simplification and visualization of large maps},
author={Mustafa, Nabil and Krishnan, Shankar and Varadhan, Gokul and Venkatasubramanian, Suresh},
journal={International Journal of Geographical Information Science},
volume={20},
number={3},
pages={273--302},
year={2006},
publisher={Taylor \& Francis}
}
@article{nollenburg2008morphing,
title={Morphing polylines: A step towards continuous generalization},
author={N{\"o}llenburg, Martin and Merrick, Damian and Wolff, Alexander and Benkert, Marc},
journal={Computers, Environment and Urban Systems},
volume={32},
number={4},
pages={248--260},
year={2008},
publisher={Elsevier}
}
@inproceedings{stanislawski2012automated,
title={Automated metric assessment of line simplification in humid landscapes},
author={Stanislawski, Lawrence V and Raposo, Paulo and Howard, Michael and Buttenfield, Barbara P},
booktitle={Proceedings of the AutoCarto},
year={2012}
}
% LIKELY UNNEEDED
@book{buttenfield1991map,
title={Map Generalization: Making rules for knowledge representation},
author={Buttenfield, Barbara Pfeil and McMaster, Robert Brainerd},
year={1991},
publisher={Longman Scientific \& Technical London}
}
@article{chaikin1974algorithm,
title={An algorithm for high-speed curve generation},
author={Chaikin, George Merrill},
journal={Computer graphics and image processing},
volume={3},
number={4},
pages={346--349},
year={1974},
publisher={Elsevier}
}
@article{knuth1976big,
title={Big omicron and big omega and big theta},
author={Knuth, Donald E},
journal={ACM Sigact News},
volume={8},
number={2},
pages={18--24},
year={1976},
publisher={ACM New York, NY, USA}
}
@book{bachmann1894analytische,
title={Die analytische zahlentheorie},
author={Bachmann, Paul},
volume={2},
year={1894},
publisher={Teubner}
}
@article{landau1911,
title={Handbuch der Lehre von der Verteilung der Primzahlen},
journal={Monatshefte f{\"u}r Mathematik und Physik},
year={1911},
month={12},
day={01},
volume={22},
number={1},
pages={A26-A26},
issn={1436-5081},
doi={10.1007/BF01742852},
}
@online{mappingunits,
author={Aileen Buckley},
title={Guidelines for minimum size for text and symbols on maps},
date={2008-01-16},
url={https://www.esri.com/arcgis-blog/products/product/mapping/guidelines-for-minimum-size-for-text-and-symbols-on-maps/},
organization={Esri},
urldate={2021-05-03},
}
@online{cartoucheMinimalDimensions,
author={CartouCHe},
title={Cartographic Design for Screen Maps},
subtitle={Minimum Dimensions},
date={2012-01-26},
url={http://www.e-cartouche.ch/content_reg/cartouche/cartdesign/en/html/GenRules_learningObject3.html},
urldate={2021-05-03},
}
@online{epsg3857,
author={MapTiler Team},
title={WGS 84/Pseudo-Mercator},
url={https://epsg.io/3857},
urldate={2021-05-03},
}
@online{postgis311,
author={PostGIS Team},
title={PostGIS 3.1.1},
url={https://postgis.net/2021/01/28/postgis-3.1.1/},
urldate={2021-05-12},
}
@online{postgisref,
author={PostGIS Team},
title={PostGIS Reference},
url={https://postgis.net/docs/reference.html},
urldate={2021-05-12},
}
@online{wmsql,
author={Motiejus Jakštys},
title={Wang--M{\"u}ller implementation in PostGIS},
url={https://github.com/motiejus/wm},
urldate={2021-05-19},
}
@online{openstreetmap,
author={OpenStreetMap contributors},
title={Project that creates and distributes free world's geographic data},
url={https://www.openstreetmap.org},
urldate={2021-05-15},
}
@online{nzt,
author={Nacionalinė Žemės Tarnyba Prie Žemės Ūkio Ministerijos},
title={Atviri Duomenys},
url={http://nzt.lt/go.php/lit/Atviri-duomenys},
urldate={2021-05-15},
}
@online{openmapwm,
author={Tomas Straupis},
title={Test harness for Wang--M{\"u}ller algorithm},
url={https://dev.openmap.lt/webgl/wm.html},
urldate={2021-05-15},
}
@article{devangleserrorbends,
author={Gökgöz, Türkay and Sen, Alper and Memduhoğlu, Abdulkadir and Hacar, Müslüm},
year={2015},
month={10},
pages={2185-2204},
title={A New Algorithm for Cartographic Simplification of Streams and Lakes Using Deviation Angles and Error Bands},
volume={4},
journal={ISPRS International Journal of Geo-Information},
doi={10.3390/ijgi4042185}
}

49
IV/db
View File

@ -1,49 +0,0 @@
#!/bin/bash
set -euo pipefail
name=wm-mj
_psql() {
env \
PGPASSWORD=osm \
PGHOST=127.0.0.1 \
PGUSER=osm \
PGDATABASE=osm \
psql "$@"
}
_wait_for_postgres() {
>&2 echo -n "Waiting for postgres"
for _ in $(seq 240); do
if _psql -qc '\q' 2>/dev/null; then
>&2 echo " up"
exit 0
fi
>&2 echo -n .
sleep 1
done
>&2 echo " down"
exit 1
}
case ${1:-} in
start)
_psql -qc '\q' 2>/dev/null && exit 0
docker run -d --rm \
--net=host \
-e POSTGRES_DBNAME=osm \
-e POSTGRES_USER=osm \
-e POSTGRES_PASSWORD=osm \
--name "$name" \
postgis/postgis:13-3.1-alpine \
-c log_statement=all \
-c listen_addresses=127.0.0.1
_wait_for_postgres
;;
stop)
docker stop "$name"
;;
*)
_psql "$@"
;;
esac

View File

@ -1,11 +0,0 @@
#!/bin/bash -eu
s=${1:-mj-msc-full.pdf}
d=$(mktemp -d)
f=mj-msc.pdf
l="$d/make.log"
echo "Extracting $s to workdir $d/"; pdfdetach -saveall -o "$d" "$s"
echo "Logs in $l ..."; make -j "$(nproc)" -C "$d" "$f" &> "$l" || {
echo "Failed to generate. $l extract:"; tail -20 "$l"; exit 1
}
echo "Opening $d/$f ..."; xdg-open "$d/$f"
echo "$d/$f was closed. Removing $d"; rm -r "$d"

View File

@ -1,34 +0,0 @@
#!/bin/bash
set -euo pipefail
dbzip=$(realpath -s "$1")
tbl=$2
filter=$3
outfile=$4
if [[ "$dbzip" =~ " " ]]; then
echo "ERROR: $dbzip contains spaces"
exit 1
fi
set -x
gdbname=$(unzip -Z1 "$dbzip" | awk -Fgdb '/.gdb/{print $1"gdb";exit}')
now=$(TZ=UTC date +"%FT%TZ")
here=$(whoami)@$(hostname -f)
mkdir -p ".tmp"
ogr2ogr -skipfailures -f PGDump /dev/stdout -t_srs epsg:3857 \
"/vsizip/$dbzip/$gdbname" -nln "src_$tbl" hidro_l | \
awk "!/^INSERT/{print}; /^INSERT/&&/${filter}/{print;next}" | \
bash ./db | \
grep -v '^INSERT 0 1'
bash db -f aggregate-rivers.sql -v "srctable=src_$tbl" -v "dsttable=$tbl"
(
echo "-- Generated at $now on $here";
echo "-- Rivers: $filter";
docker exec wm-mj pg_dump --clean -Uosm osm -t "$tbl" | tr -d '\r'
) > ".tmp/$outfile"
mv ".tmp/$outfile" "$outfile"

View File

@ -1,23 +0,0 @@
#!/bin/bash
set -euo pipefail
# Prefix the 'make <...>' with this script to build the artifact in an isolated
# container. This means host dependencies can only be Docker and a shell (to
# run this script).
#
# Usage:
# ./in-container make help
# ./in-container make -j mj-msc-full.pdf
# ...
NAME=wm-mj-build
if [[ -z "$(docker images -q --filter "reference=$NAME")" ]]; then
docker build -t "$NAME" .
fi
exec docker run -ti --rm \
--net=host \
-v /var/run/docker.sock:/var/run/docker.sock \
-v $(git rev-parse --show-toplevel):/x \
-w /x/$(basename ${PWD}) \
"$NAME" "$@"

View File

@ -1,33 +0,0 @@
-- This file initializes tables for unit and river tests.
-- ST_SimplifyWM, when dbgname is non-empty, expects `wm_debug` table to be
-- created.
-- to preview this somewhat conveniently in QGIS:
-- stage || '_' || name || ' gen:' || coalesce(gen, 'Ø') || ' nbend:' || lpad(nbend, 4, '0')
drop table if exists wm_debug;
create table wm_debug(
id serial,
stage text not null,
name text not null,
gen bigint not null,
nbend bigint,
way geometry,
props jsonb
);
drop table if exists wm_manual;
create table wm_manual (
id serial,
name text,
way geometry,
props jsonb
);
-- Run ST_SimplifyWM in debug mode, so `wm_debug` is populated. That table
-- is used for geometric assertions later in the file.
drop table if exists wm_demo;
create table wm_demo (name text, i bigint, way geometry);
-- wm_visuals holds visual aids for the paper.
drop table if exists wm_visuals;
create table wm_visuals (name text, way geometry);

View File

@ -1,133 +0,0 @@
#!/usr/bin/python3
"""
Convert PostGIS geometries to an image. To scale.
Accepts a few geometry fine-tuning parameters.
"""
import argparse
import geopandas
import psycopg2
from matplotlib import rc
import matplotlib.pyplot as plt
CMAP = 'tab20c' # 'Set3' # is nice too
PSQL_CREDS = "host=127.0.0.1 dbname=osm user=osm password=osm"
COLORS = {
'black': '#000000',
'green': '#1b9e77',
'orange': '#d95f02',
'purple': '#7570b3',
}
# see `NOTICE` in the LaTeX document; this is the width of the main text block.
TEXTWIDTH_CM = 12.12364
QUADRANTS = {'tr':1, 'tl':2, 'bl':3, 'br':4}
def color(string):
return COLORS[string if string else 'black']
def inch(cm):
return cm / 2.54
def parse_args():
kwcolor = {'type': color, 'default': 'black'}
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('--g1-select')
parser.add_argument('--g1-linestyle')
parser.add_argument('--g1-label')
parser.add_argument('--g1-color', **kwcolor)
parser.add_argument('--g2-select')
parser.add_argument('--g2-linestyle')
parser.add_argument('--g2-label')
parser.add_argument('--g2-color', **kwcolor)
parser.add_argument('--g3-select')
parser.add_argument('--g3-linestyle')
parser.add_argument('--g3-label')
parser.add_argument('--g3-color', **kwcolor)
parser.add_argument('--legend',
help="Legend location, following matplotlib rules", default='best')
parser.add_argument('--widthdiv', default=1, type=float,
help="Divide the width by this number "
"(useful when two images are laid horizontally "
"in the resulting file")
parser.add_argument('--quadrant', choices=QUADRANTS.keys(),
help="Image is comprised of 4 quadrants. This variable, "
"when non-empty, will clip and return the requested quadrant")
parser.add_argument('--outfile', metavar='<file>',
help="If unset, displayed on the screen")
return parser.parse_args()
def read_layer(select, width, maybe_quadrant):
if not select:
return
way = "way"
if maybe_quadrant:
way = "wm_quadrant(way, {})".format(QUADRANTS[maybe_quadrant])
conn = psycopg2.connect(PSQL_CREDS)
sql = "SELECT {way} as way1 FROM {select}".format(way=way, select=select)
return geopandas.read_postgis(sql, con=conn, geom_col='way1')
def plot_args(geom, color, maybe_linestyle, maybe_label):
if geom is None:
return
# polygons either have fillings or lines
if geom.geom_type[0] == 'Polygon':
if maybe_linestyle:
return {
'edgecolor': 'black',
'linestyle': maybe_linestyle,
'color': (0, 0, 0, 0),
}
else:
return {'cmap': CMAP, 'alpha': .25}
r = {'color': color}
if maybe_linestyle == 'invisible':
r['color'] = (0, 0, 0, 0)
elif maybe_linestyle:
r['linestyle'] = maybe_linestyle
if maybe_label:
r['label'] = '\\normalfont %s' % maybe_label
return r
def main():
args = parse_args()
width = TEXTWIDTH_CM / args.widthdiv
g1 = read_layer(args.g1_select, width, args.quadrant)
g2 = read_layer(args.g2_select, width, args.quadrant)
g3 = read_layer(args.g3_select, width, args.quadrant)
c1 = plot_args(g1, args.g1_color, args.g1_linestyle, args.g1_label)
c2 = plot_args(g2, args.g2_color, args.g2_linestyle, args.g2_label)
c3 = plot_args(g3, args.g3_color, args.g3_linestyle, args.g3_label)
rc('text', usetex=True)
rc('text.latex', preamble='\\usepackage{numprint}\n')
fig, ax = plt.subplots(constrained_layout=True)
fig.set_figwidth(inch(width))
g1 is not None and g1.plot(ax=ax, linewidth=.75, **c1)
g2 is not None and g2.plot(ax=ax, linewidth=.75, **c2)
g3 is not None and g3.plot(ax=ax, linewidth=.75, **c3)
ax.legend(loc=args.legend, frameon=False)
ax.axis('off')
ax.margins(0, 0)
if args.outfile:
fig.savefig(args.outfile, bbox_inches='tight', dpi=600)
else:
plt.show()
if __name__ == '__main__':
main()

File diff suppressed because it is too large Load Diff

View File

@ -1,117 +0,0 @@
Self-line crossing when cutting a bend
--------------------------------------
The self-line-crossing may happen after a few bends have been skipped. E.g.
ends of A<->B cross the line, but "swallow" a few more in between:
,______
/ \
|___A | \ \
\ | B\ | __
\ | | | / \
/ | | |___,---,___/A |
/ | \_________________|
\ |
\ | \ \
/ / B\ | _ __
----/ / | | / \ / \
/ ,____/ | |___/ \___/A |
/ B| \_________________|
|
If a bend with 180+ deg sum of inflection angles is found, its line between
inflection angles (AB in our examples) must be crossed with all the other bends
to detect a possible line-crossing. This is O(N*M), where N is the total number
of line segments, and M is the number of qualifying bends. It is expensive.
Also, there is another way to remove self-crossing, without removing most of
the bend. E.g. from:
\ /
B\ | __
| | / \
| |____/A |
\__________|
Instead of:
\ /
\/ A'
B
To:
\ \_
B\ `-,_.__
| A' \
| |
\__________|
But perhaps it doesn't look quite as natural. I will trust the original
article to do the right thing here and remove the bend altogether.
ALSO: the bends should be iterated from different directions:
for i := 0; i < len(bends); i++ {
for j := 0; j < i; j++ {
...
}
for j := len(bends); j > i; j-- {
...
}
}
So if there are multiple bends between the baseline, they will be cut
correctly.
The Context of a Bend
---------------------
Similar bends:
> For example, if bend 1 has four unit areas and bend 2 has six unit areas, the
> average size is five units, and the normalized areas of bends 1 and 2 are
> 4/5=0.8 and 6/5=1.2, respectively.
My comment: everything until this sentence is clear. However, "unit areas" is
misleading: there is little reason to normalize areas, but leave the distances
intact (if we'd like to normalize areas, it would make sense to square-root
them).
Removing that removes changes the meaning of the sentence that **euclidean
distance** is normalized (the composite of the bend properties), rather than
a single component.
Offered structure
-----------------
- Introduction
- Previous research overview
- Methodology and methodics
- Results
- Conclusions
- Literature review
- Appendix
for 2021-04-19
--------------
- literatūros šaltinių analizė
- literatūros šaltinių priskyrimas atskiroms magistro darbo struktūrinėms dalims
analizės uždaviniai:
- galutinis problemos formulavimas
- darbo tikslo formulavimas
- uždavinių formulavimas
- aktualumo
- naujumo
- pritaikomumo formulavimas
Angl.:
- trūksta literatūros apžvalgos: pervadinti šiuolaikinius sprendimus į tai
- mažiausiai 2 poskyriai skyriuje.
- techninė dalis -- į rezultatus.
- "eksperimento rezultatai" eina į "darbo rezultatus".

Binary file not shown.

Before

Width:  |  Height:  |  Size: 65 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 198 KiB

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

View File

@ -1,90 +0,0 @@
---
title: Wang & Müller linijų generalizacijos algoritmo įgyvendinimas upių pavyzdžiu
author: Motiejus Jakštys
description: |
Prezentacija apie magistrinio darbo progresą vertinimo grupei
date: |
2021-03-29 \
![](../misc/Logo_vilniaus_universitetas.png){width=2in}
lang: lt-LT
header-includes: |
\definecolor{vulightgrey}{RGB}{220,220,220}
\definecolor{vudarkgrey}{RGB}{65,65,65}
\definecolor{vupurple}{RGB}{123,0,63}
\definecolor{darkgreen}{RGB}{32,96,32}
\setbeamercolor{title}{fg=vupurple}
\setbeamercolor{frametitle}{fg=vupurple}
\setbeamercolor{item}{fg=vupurple}
\setbeamercolor{normal text}{fg=vudarkgrey}
...
# Problema
- Dabartiniai atvirai prieinami algoritmai
- Visvalingam-Whyatt
- Douglas-Peucker
- Problematika
- nepritaikyti gamtiniams objektams: upėms, krantų linijoms
- prarandami raiškūs gamtinių objektų elementai: vingiai, kilpos
- Alternatyvos
- egzistuoja tik teorinės
- nėra implementacijos
- net teorinėms -- neapibrėžti parametrai
# Aktualumas
- Smulkesnių mastelių žemėlapių georeferenciniai duomenys kuriami taikant
automatizuoto generalizavimo metodus (pvz., GDR10LT generaluzuojant į
GDR50LT, GDR250LT).
- Nėra atvirai prieinamo sprendimo gamtinių linijinių objektų generalizacijai
# Etapai: Wang & Müller techninio įgyvendinimo
::: nonincremental
- Įlinkių atpažinimas
- Silpnas pasislinkimas įlinkio gale
- Įlinkio kirtimas kitu įlinkiu
- Patikrinimas su Nemunu, Šalčia ir Visinčia
- Įlinkių matavimai:
- dydis ir forma
- izoliuoti, panašūs
- Operatoriai:
- Specifinio įlinkio pašalinimo
- Dviejų įlinkių kombinavimo į vieną
- Izoliuotų įlinkių padidinimo
- Pritaikymas visoms Lietuvos upėms
:::
# Etapai: analizės ir darbo rašymo
::: nonincremental
- Esamų algoritmų apžvalga
- Problemos aprašymas
- Pasirinkto generalizacijos sprendimo apibūdinimas
- Sukurtos implementacijos aptarimas
- privalumai
- trūkumai
- Nauda praktikoje
- Tolimesnių tyrimų kryptys
:::
# Naujumas
Šis darbas pasiūlys iki šiol neegzistuojantį atvirą sprendimą gamtinių
linijinių objektų generalizacijai (Lietuvoje -- GDR50 ir GDR250 analogus)
# Tolimesnis pritaikymas
- Sukurta algoritmo techninė realizacija atviro kodo GIS bus pradžia tolesniam jo
pritaikymui kitiems gamtiniams objektams: ežerų kranto linija, jūrų,
kontinentų, miškų ribos.
- Kadangi rezultatas bus techninė realizacija, o ne rašinėlis, tai suteikia
toliau galimybę vystyti sprendimą.
# Klausimai
Ačiū už dėmesį.

View File

@ -1,20 +0,0 @@
\set ON_ERROR_STOP on
SET plpgsql.extra_errors TO 'all';
-- This fails with real rivers since dcf4c02307baeece51470a961a113a8fad68fad5
-- (adding GDB10LT data). The same rivers from OpenStreetMaps work.
-- There seems to be a bug in wm_exaggeration.
do $$
declare
npoints bigint;
secs bigint;
begin
select * from ST_SimplifyWM_Estimate((select st_union(way) from wm_rivers)) into npoints, secs;
raise notice 'Total points: %', npoints;
raise notice 'Expected duration: %s (+-%s)', ceil(secs), floor(secs*.5);
end $$ language plpgsql;
delete from wm_debug where name in (select distinct name from wm_rivers);
delete from wm_demo where name in (select distinct name from wm_rivers);
insert into wm_demo (name, way) select name, ST_SimplifyWM(way, 75, null, name) from wm_rivers;

View File

@ -1,198 +0,0 @@
\i wm.sql
-- https://stackoverflow.com/questions/19982373/which-tools-libraries-do-you-use-to-unit-test-your-pl-pgsql
create or replace function assert_equals(expected anyelement, actual anyelement) returns void as $$
begin
if expected = actual or (expected is null and actual is null) then
--do nothing
else
raise exception 'Assertion Error. Expected <%> but was <%>', expected, actual;
end if;
end $$ language plpgsql;
drop function if exists dbg_geomsummary;
create function dbg_geomsummary(geoms geometry[]) returns void as $$
declare i int4;
begin
raise notice 'len: %', array_length(geoms, 1);
for i in 1..array_length(geoms, 1) loop
raise notice '% %: %', st_geometrytype(geoms[i]), lpad(i::text, 2, '0'), st_astext(geoms[i]);
end loop;
end
$$ language plpgsql;
drop table if exists wm_figures;
create table wm_figures (name text, way geometry);
-- add fig8.gpkg to postgis:
-- ogr2ogr -update -f PostgreSQL PG:"host=127.0.0.1 user=osm password=osm dbname=osm" fig8.gpkg
-- to "normalize" a new line when it's in `f`:
-- select st_astext(st_snaptogrid(st_transscale(geometry, -19.5, .016, 4000, 4000), 1)) from f;
insert into wm_figures (name, way) values
('fig3', 'LINESTRING(0 0,12 0,13 4,20 2,20 0,32 0,33 10,38 16,43 15,44 10,44 0,60 0)'::geometry),
('fig3-1', 'LINESTRING(0 0,12 0,13 4,20 2,20 0,32 0,33 10,38 16,43 15,44 10,44 0)'::geometry),
('fig5', 'LINESTRING(0 39,19 52,27 77,26 104,41 115,49 115,65 103,65 75,53 45,63 15,91 0)'::geometry),
('fig6', 'LINESTRING(84 47,91 59,114 64,122 80,116 92,110 93,106 106,117 118,136 107,135 76,120 45,125 39,141 39,147 32)'::geometry),
('fig8', 'LINESTRING(173 12,174 10,180 8,186 8,186 13,191 11,189 6,201 5,203 11,216 16,216 6,222 6,229 3,236 2,239 6,243 8,248 6)'::geometry),
('inflection-1', 'LINESTRING(110 24,114 20,133 20,145 15,145 0,136 8,123 10,114 10,111 2)'::geometry),
('multi-island', 'MULTILINESTRING((-15 10,-10 10,-5 11,0 11,5 11,10 10,11 9,13 10,15 9),(-5 11,-2 15,0 16,2 15,5 11))'::geometry),
('selfcrossing-1','LINESTRING(-27 180,-20 166,-21 142,-18 136,55 136,55 136,71 145,44 165,37 146,22 145,14 164,11 164,3 146,-12 146,-13 176,-18 184)'::geometry),
('isolated-1', 'LINESTRING(-50 103,-48 102,-30 103,-31 105,-31 107,-27 107,-26 103,-6 103,-4 104)'::geometry),
('isolated-2', 'LINESTRING(250 100,246 104,234 105,230 106,225 101,224 93,217 78,206 69)'::geometry);
insert into wm_figures (name, way) values ('fig6-rev',ST_Reverse(
ST_Translate((select way from wm_figures where name='fig6'), 60, 0)));
insert into wm_figures (name, way) values ('fig6-combi', ST_Union(
ST_Translate((select way from wm_figures where name='fig6'), 0, 90),
ST_Translate((select way from wm_figures where name='fig6'), 80, 90)
));
insert into wm_figures (name, way) values ('selfcrossing-1-rev',ST_Reverse(
ST_Translate((select way from wm_figures where name='selfcrossing-1'), 0, 60)));
-- 3395 is now "reserved" for figures.
update wm_figures set way=st_setsrid(way, 3395);
delete from wm_debug where name in (select distinct name from wm_figures);
delete from wm_demo where name in (select distinct name from wm_figures);
insert into wm_demo (name, way) select name, ST_SimplifyWM(way, .1, null, name) from wm_figures where name not in ('fig8', 'isolated-1');
insert into wm_demo (name, way) select name, ST_SimplifyWM(way, 14, null, name) from wm_figures where name in ('fig8', 'isolated-1', 'isolated-2');
drop function if exists wm_debug_get;
create function wm_debug_get( _stage text, _name text, OUT ways geometry[]) as $$
declare
begin
ways = array((select way from wm_debug where stage=_stage and name=_name order by id));
end $$ language plpgsql;
do $$
declare
vbends geometry[];
begin
vbends = wm_debug_get('bbends', 'fig3');
perform assert_equals(5, array_length(vbends, 1));
perform assert_equals('LINESTRING(0 0,12 0,13 4)', st_astext(vbends[1]));
perform assert_equals('LINESTRING(12 0,13 4,20 2,20 0)', st_astext(vbends[2]));
perform assert_equals('LINESTRING(20 2,20 0,32 0,33 10)', st_astext(vbends[3]));
perform assert_equals('LINESTRING(32 0,33 10,38 16,43 15,44 10,44 0)', st_astext(vbends[4]));
perform assert_equals(4, array_length(wm_detect_bends((select way from wm_figures where name='fig3-1')), 1));
select wm_detect_bends((select way from wm_figures where name='fig5')) into vbends;
perform assert_equals(3, array_length(vbends, 1));
end $$ language plpgsql;
do $$
declare
vbends geometry[];
vinflections geometry[];
begin
vinflections = wm_debug_get('cinflections', 'fig5');
perform assert_equals('LINESTRING(0 39,19 52,27 77)', st_astext(vinflections[1]));
perform assert_equals('LINESTRING(19 52,27 77,26 104,41 115,49 115,65 103,65 75,53 45)', st_astext(vinflections[2]));
perform assert_equals('LINESTRING(65 75,53 45,63 15,91 0)', st_astext(vinflections[3]));
-- inflections-1, the example in fix_gentle_inflections docstring
select array((select way from wm_debug where name='inflection-1' and stage='bbends')) into vbends;
select array((select way from wm_debug where name='inflection-1' and stage='cinflections')) into vinflections;
perform assert_equals(vbends[1], vinflections[1]); -- unchanged
perform assert_equals('LINESTRING(114 20,133 20,145 15,145 0,136 8,123 10,114 10)', st_astext(vinflections[2]));
perform assert_equals('LINESTRING(123 10,114 10,111 2)', st_astext(vinflections[3]));
end $$ language plpgsql;
do $$
declare
fig6 constant text default 'LINESTRING(84 47,91 59,114 64,120 45,125 39,141 39,147 32)';
selfcrossing1 constant text default 'LINESTRING(-27 180,-20 166,-13 176,-18 184)';
vcrossings geometry[];
mutated boolean;
begin
select * from wm_self_crossing(wm_debug_get('cinflections', 'fig6')) into vcrossings, mutated;
perform assert_equals(true, mutated);
perform assert_equals(
fig6,
(select st_astext(
st_linemerge(st_union(way))
) from (select unnest(vcrossings) way) a)
);
select * from wm_self_crossing(wm_debug_get('cinflections', 'fig6-rev')) into vcrossings, mutated;
perform assert_equals(true, mutated);
perform assert_equals(
fig6,
(select st_astext(
st_translate(st_reverse(st_linemerge(st_union(way))), -60, 0)
) from (select unnest(vcrossings) way) a)
);
select * from wm_self_crossing(wm_debug_get('cinflections', 'fig6-combi')) into vcrossings, mutated;
perform assert_equals(true, mutated);
perform assert_equals(
'MULTILINESTRING((84 137,91 149,114 154,120 135,125 129,141 129,147 122),(164 137,171 149,194 154,200 135,205 129,221 129,227 122))',
(select st_astext(
st_linemerge(st_union(way))
) from (select unnest(vcrossings) way) a)
);
select * from wm_self_crossing(wm_debug_get('cinflections', 'selfcrossing-1')) into vcrossings, mutated;
perform assert_equals(true, mutated);
perform assert_equals(
selfcrossing1,
(select st_astext(
st_linemerge(st_union(way))
) from (select unnest(vcrossings) way) a)
);
select * from wm_self_crossing(wm_debug_get('cinflections', 'selfcrossing-1-rev')) into vcrossings, mutated;
perform assert_equals(true, mutated);
perform assert_equals(
selfcrossing1,
(select st_astext(
st_translate(st_reverse(st_linemerge(st_union(way))), 0, -60)
) from (select unnest(vcrossings) way) a)
);
end $$ language plpgsql;
-- verifying bends in fig8 are eliminated like explained in the WM paper
do $$
declare
fig8gen2 constant text default 'LINESTRING(173 12,174 10,180 8,186 8,189 6,201 5,203 11,216 16,216 6,229 3,236 2,239 6,243 8,248 6)';
fig8gen3 constant text default 'LINESTRING(173 12,174 10,180 8,189 6,201 5,203 11,216 16,216 6,229 3,236 2,239 6,243 8,248 6)';
eliminations geometry[];
begin
eliminations = wm_debug_get('afigures', 'fig8');
perform assert_equals(fig8gen2, st_astext(eliminations[2]));
perform assert_equals(fig8gen3, st_astext(eliminations[3]));
end $$ language plpgsql;
-- testing wm_exaggerate_bend2 in isolation
do $$
declare
fig3b2 geometry;
bend geometry;
size float;
begin
select way from wm_debug where name='fig3' and stage='bbends' and gen=1 and nbend=2 into fig3b2;
size = wm_adjsize(fig3b2);
bend = wm_exaggerate_bend2(fig3b2, size, 50.);
perform assert_equals('ST_LineString', st_geometrytype(bend));
insert into wm_debug(stage, name, gen, nbend, way) values('manual', 'fig3', 1, 1, bend);
end $$ language plpgsql;
-- misc visuals
do $$
declare fig6b1 geometry;
declare fig6b2 geometry;
declare sclong geometry;
declare scshort geometry;
begin
select way from wm_debug where name='fig6' and stage='bbends' and gen=1 into fig6b1 limit 1 offset 0;
select way from wm_debug where name='fig6' and stage='bbends' and gen=1 into fig6b2 limit 1 offset 2;
insert into wm_visuals (name, way) values('fig6-baseline', st_makeline(st_startpoint(fig6b2), st_endpoint(fig6b2)));
insert into wm_visuals (name, way) values('fig6-newline', st_makeline(st_endpoint(fig6b1), st_endpoint(fig6b2)));
select way from wm_debug where name='selfcrossing-1' and stage='bbends' and gen=1 into sclong limit 1 offset 1;
select way from wm_debug where name='selfcrossing-1' and stage='bbends' and gen=1 into scshort limit 1 offset 4;
insert into wm_visuals (name, way) values('selfcrossing-1-baseline', st_makeline(st_startpoint(sclong), st_endpoint(sclong)));
insert into wm_visuals (name, way) values('selfcrossing-1-newline', st_makeline(st_startpoint(sclong), st_endpoint(scshort)));
end $$ language plpgsql;

View File

@ -1,26 +0,0 @@
#!/usr/bin/awk -f
BEGIN { FS="[(); ]" }
/small_angle constant real default radians/ {
x1 += 1;
d1 = sprintf("\\newcommand{\\smallAngle}{$%d^\\circ$}",$8);
}
/isolation_threshold constant real default / {
x2 += 1;
d2 = sprintf("\\newcommand{\\isolationThreshold}{%.1f}",$7);
}
/scale2 constant float default / {
x3 += 1;
d3 = sprintf("\\newcommand{\\exaggerationEnthusiasm}{%.1f}",$7);
}
END{
if(x1 == 1 && x2 == 1 && x3 == 1) {
print d1 > "vars.inc.tex"
print d2 >> "vars.inc.tex"
print d3 >> "vars.inc.tex"
} else {
exit 1
}
}

View File

@ -1,146 +0,0 @@
\set ON_ERROR_STOP on
SET plpgsql.extra_errors TO 'all';
-- wm_bbox clips a geometry by a bounding box around a given object,
-- matching dimensions of A-class paper (1 by sqrt(2)).
drop function if exists wm_bbox;
create function wm_bbox(
center geometry,
scaledwidth float
) returns geometry as $$
declare
halfX float;
halfY float;
begin
halfX = scaledwidth / 2;
halfY = halfX * sqrt(2);
return st_envelope(
st_union(
st_translate(center, halfX, halfY),
st_translate(center, -halfX, -halfY)
)
);
end $$ language plpgsql;
-- wm_quadrant divides the given geometry to 4 rectangles
-- and returns the requested quadrant following cartesian
-- convention:
-- +----------+
-- | II | I |
--- +----------+
-- | III | IV |
-- +-----+----+
-- matching dimensions of A-class paper (1 by sqrt(2).
drop function if exists wm_quadrant;
create function wm_quadrant(
geom geometry,
quadrant integer
) returns geometry as $$
declare
xmin float;
xmax float;
ymin float;
ymax float;
begin
xmin = st_xmin(geom);
xmax = st_xmax(geom);
ymin = st_ymin(geom);
ymax = st_ymax(geom);
if quadrant = 1 or quadrant = 2 then
ymin = (ymin + ymax)/2;
else
ymax = (ymin + ymax)/2;
end if;
if quadrant = 2 or quadrant = 3 then
xmax = (xmin + xmax)/2;
else
xmin = (xmin + xmax)/2;
end if;
return st_intersection(
geom,
st_makeenvelope(xmin, ymin, xmax, ymax, st_srid(geom))
);
end $$ language plpgsql;
drop function if exists wm_salvisbbox;
create function wm_salvisbbox(
geom geometry,
scaledwidth float
) returns geometry as $$
declare
ret geometry;
begin
with multismall as (
select st_intersection(
st_union(geom),
wm_bbox(
st_closestpoint(
(select way from wm_rivers where name='Šalčia'),
(select way from wm_rivers where name='Visinčia')
),
scaledwidth
)
) ways
)
-- protecting against very small bends that were cut
-- in the corner of the picture
select st_union(a.geom)
from st_dump((select ways from multismall)) a
where st_length(a.geom) >= 100
into ret;
return ret;
end $$ language plpgsql;
delete from wm_debug where name like 'salvis%';
delete from wm_visuals where name like 'salvis%';
insert into wm_visuals(name, way) values
('salvis-grpk10', (
wm_salvisbbox(
(select st_union(way) from wm_rivers where name in ('Šalčia', 'Visinčia')),
:scaledwidth
)
)),
('salvis-grpk50', (
wm_salvisbbox(
(select st_union(way) from wm_rivers_50 where name in ('Šalčia', 'Visinčia')),
:scaledwidth
)
)),
('salvis-grpk250', (
wm_salvisbbox(
(select st_union(way) from wm_rivers_250 where name in ('Šalčia', 'Visinčia')),
:scaledwidth
)
));
do $$
declare
i integer;
geom1 geometry;
geom2 geometry;
geom3 geometry;
begin
foreach i in array array[16, 32, 64, 256] loop
geom1 = st_simplify((select way from wm_visuals where name='salvis-grpk10'), i);
geom2 = st_simplifyvw((select way from wm_visuals where name='salvis-grpk10'), i*i);
insert into wm_visuals(name, way) values
('salvis-dp' || i, geom1),
('salvis-dpchaikin' || i, st_chaikinsmoothing(geom1, 5)),
('salvis-vw' || i, geom2),
('salvis-vwchaikin' || i, st_chaikinsmoothing(geom2, 5));
end loop;
-- more than 220 doesn't work, because there is an exaggerated bend near
-- Šalčia-Visinčia crossing, and it "exaggerates" to the
-- other river.
foreach i in array array[75, 220] loop
geom3 = st_simplifywm((select way from wm_visuals where name='salvis-grpk10'), i, 50, 'salvis-wm' || i);
insert into wm_visuals(name, way) values
('salvis-wm' || i, geom3);
end loop;
end $$ language plpgsql;

BIN
IV/vu.pdf

Binary file not shown.

Binary file not shown.

Before

Width:  |  Height:  |  Size: 85 KiB

850
IV/wm.sql
View File

@ -1,850 +0,0 @@
\set ON_ERROR_STOP on
SET plpgsql.extra_errors TO 'all';
-- wm_detect_bends detects bends using the inflection angles. No corrections.
drop function if exists wm_detect_bends;
create function wm_detect_bends(
line geometry,
dbgname text default null,
dbggen integer default null,
OUT bends geometry[]
) as $$
declare
p geometry;
p1 geometry;
p2 geometry;
p3 geometry;
bend geometry;
prev_sign int4;
cur_sign int4;
l_type text;
dbgpolygon geometry;
begin
l_type = st_geometrytype(line);
if l_type != 'ST_LineString' then
raise 'This function works with ST_LineString, got %', l_type;
end if;
-- 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:
--
-- p1___ ...
-- /
-- ... _____/
-- p3 p2
--
-- 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].
-- 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.
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
p3 = p2;
p2 = p1;
p1 = p;
continue when p3 is null;
cur_sign = sign(pi() - st_angle(p1, p2, p2, p3));
if bend is null then
bend = st_makeline(p3, p2);
else
bend = st_linemerge(st_union(bend, st_makeline(p3, p2)));
end if;
if prev_sign + cur_sign = 0 then
if bend is not null then
bends = bends || bend;
end if;
bend = st_makeline(p3, p2);
end if;
prev_sign = cur_sign;
end loop;
-- 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
bends = bends || bend;
end if;
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, dbggen, i, bends[i]);
dbgpolygon = null;
if st_npoints(bends[i]) >= 3 then
dbgpolygon = st_makepolygon(
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;
-- wm_fix_gentle_inflections moves bend endpoints following "Gentle Inflection
-- at End of a Bend" section.
--
-- 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
-- cumulative inflection angle is small (see variable below).
--
-- The implementation could be significantly optimized to avoid `st_reverse`
-- and array reversals, trading for complexity in wm_fix_gentle_inflections1.
drop function if exists wm_fix_gentle_inflections;
create function wm_fix_gentle_inflections(
INOUT bends geometry[],
dbgname text default null,
dbggen integer default null
) as $$
declare
len int4;
bends1 geometry[];
dbgpolygon geometry;
begin
len = array_length(bends, 1);
bends = wm_fix_gentle_inflections1(bends);
for i in 1..len loop
bends1[i] = st_reverse(bends[len-i+1]);
end loop;
bends1 = wm_fix_gentle_inflections1(bends1);
for i in 1..len loop
bends[i] = st_reverse(bends1[len-i+1]);
end loop;
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, dbggen, i, bends[i]);
dbgpolygon = null;
if st_npoints(bends[i]) >= 3 then
dbgpolygon = st_makepolygon(
st_addpoint(bends[i],
st_startpoint(bends[i]))
);
end if;
insert into wm_debug(stage, name, gen, nbend, way) values(
'cinflections-polygon', dbgname, dbggen, i, dbgpolygon);
end loop;
end if;
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.
drop function if exists wm_fix_gentle_inflections1;
create function wm_fix_gentle_inflections1(INOUT bends geometry[]) as $$
declare
-- the threshold when the angle is still "small", so gentle inflections can
-- be joined
small_angle constant real default radians(45);
ptail geometry; -- tail point of tail bend
phead geometry[]; -- 3 tail points of head bend
i int4; -- bends[i] is the current head
begin
for i in 2..array_length(bends, 1) loop
-- 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.
--
--
-- Assume this curve (figure `inflection-1`), going clockwise from A:
--
-- \______B
-- A `-------. C
-- |
-- G___ F |
-- / `-----.____+ D
-- E
--
-- After processing the curve following the definition of a bend, the bend
-- [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].
select geom from st_dumppoints(bends[i-1])
order by path[1] asc limit 1 into ptail;
while true loop
-- copy last 3 points of bends[i-1] (tail) to ptail
select array(
select geom from st_dumppoints(bends[i]) order by path[1] asc limit 3
) into phead;
-- if the bend got too short, stop processing it
exit when array_length(phead, 1) < 3;
-- inflection angle between ptail[1:3] is "large", stop processing
exit when abs(st_angle(phead[1], phead[2], phead[3]) - pi()) > small_angle;
-- distance from head's 1st vertex should be larger than from 2nd vertex
exit when st_distance(ptail, phead[2]) < st_distance(ptail, phead[3]);
-- Between two bends, bend with smaller baseline wins when two
-- neighboring bends can have gentle inflections. This is a heuristic
-- that can be safely removed, but in practice has shown to avoid
-- creating some very bendy lines.
exit when st_distance(st_pointn(bends[i], 1), st_pointn(bends[i], -1)) <
st_distance(st_pointn(bends[i-1], 1), st_pointn(bends[i-1], -1));
-- Detected a gentle inflection.
-- Move head of the tail to the tail of head
bends[i] = st_removepoint(bends[i], 0);
bends[i-1] = st_addpoint(bends[i-1], phead[3]);
end loop;
end loop;
end $$ language plpgsql;
-- wm_if_selfcross returns whether baseline of bendi crosses bendj.
-- If it doesn't, returns a null geometry.
-- Otherwise, it will return the baseline split into a few parts where it
-- crosses bendj.
drop function if exists wm_if_selfcross;
create function wm_if_selfcross(
bendi geometry,
bendj geometry
) returns geometry as $$
declare
a geometry;
b geometry;
multi geometry;
begin
a = st_pointn(bendi, 1);
b = st_pointn(bendi, -1);
multi = st_split(bendj, st_makeline(a, b));
if st_numgeometries(multi) = 1 then
return null;
end if;
if st_numgeometries(multi) = 2 and
(st_contains(bendj, a) or st_contains(bendj, b)) then
return null;
end if;
return multi;
end $$ language plpgsql;
-- wm_self_crossing eliminates self-crossing from the bends, following the
-- article's section "Self-line Crossing When Cutting a Bend".
drop function if exists wm_self_crossing;
create function wm_self_crossing(
INOUT bends geometry[],
dbgname text default null,
dbggen integer default null,
OUT mutated boolean
) as $$
declare
i int4;
j int4;
multi geometry;
begin
mutated = false;
<<bendloop>>
for i in 1..array_length(bends, 1) loop
continue when abs(wm_inflection_angle(bends[i])) <= pi();
-- 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
-- Go through each bend in the given line, and see if has a potential to
-- cross bends[i]. The line-cut process is different when i<j and i>j;
-- therefore there are two loops, one for each case.
for j in 1..i-1 loop
multi = wm_if_selfcross(bends[i], bends[j]);
continue when multi is null;
mutated = true;
-- 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);
bends[j] = st_setpoint(
bends[j],
st_npoints(bends[j])-1,
st_pointn(bends[i], st_npoints(bends[i]))
);
bends = bends[1:j] || bends[i+1:];
continue bendloop;
end loop;
for j in reverse array_length(bends, 1)..i+1 loop
multi = wm_if_selfcross(bends[i], bends[j]);
continue when multi is null;
mutated = true;
-- 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)
);
bends = bends[1:i] || bends[j+1:];
continue bendloop;
end loop;
end loop;
if dbgname is not null then
insert into wm_debug(stage, name, gen, nbend, way) values(
'dcrossings', dbgname, dbggen, generate_subscripts(bends, 1),
unnest(bends)
);
end if;
end $$ language plpgsql;
drop function if exists wm_inflection_angle;
create function wm_inflection_angle (IN bend geometry, OUT angle real) as $$
declare
p0 geometry;
p1 geometry;
p2 geometry;
p3 geometry;
begin
angle = 0;
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;
drop function if exists wm_bend_attrs;
drop function if exists wm_elimination;
drop function if exists wm_exaggeration;
drop type if exists wm_t_attrs;
create type wm_t_attrs as (
adjsize real,
baselinelength real,
curvature real,
isolated boolean
);
create function wm_bend_attrs(
bends geometry[],
dbgname text default null,
dbggen integer default null
) returns wm_t_attrs[] as $$
declare
isolation_threshold constant real default 0.5;
attrs wm_t_attrs[];
attr wm_t_attrs;
bend geometry;
i int4;
needs_curvature real;
skip_next boolean;
dbglastid integer;
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;
if st_numpoints(bend) >= 3 then
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, bends[i],
jsonb_build_object(
'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
-- invariant: two bends that touch cannot be isolated.
if st_npoints(bends[i]) > 3 then
skip_next = false;
end if;
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;
end loop;
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
drop function if exists wm_st_split;
create function wm_st_split(
input geometry,
blade geometry
) returns geometry as $$
declare
type1 text;
type2 text;
begin
type1 = st_geometrytype(input);
type2 = st_geometrytype(blade);
if not (type1 = 'ST_LineString' and
type2 = 'ST_Point') then
raise 'Arguments must be LineString and Point, got: % and %', type1, type2;
end if;
return st_split(st_snap(input, blade, 0.00000001), blade);
end $$ language plpgsql;
-- wm_exaggerate_bend2 is the second version of bend exaggeration. Uses
-- non-linear interpolation by point azimuth. Slower, but produces nicer
-- exaggerated geometries.
drop function if exists wm_exaggerate_bend2;
create function wm_exaggerate_bend2(
INOUT bend geometry,
size float,
desired_size float
) as $$
declare
scale2 constant float default 1.2; -- exaggeration enthusiasm
midpoint geometry; -- midpoint of the baseline
points geometry[];
startazimuth float;
azimuth float;
diffazimuth float;
point geometry;
sss float;
protect int = 10;
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);
startazimuth = st_azimuth(midpoint, st_pointn(bend, 1));
while (size < desired_size) and (protect > 0) loop
protect = protect - 1;
for i in 2..st_npoints(bend)-1 loop
point = st_pointn(bend, i);
azimuth = st_azimuth(midpoint, point);
diffazimuth = degrees(azimuth - startazimuth);
if diffazimuth > 180 then
diffazimuth = diffazimuth - 360;
elseif diffazimuth < -180 then
diffazimuth = diffazimuth + 360;
end if;
diffazimuth = abs(diffazimuth);
if diffazimuth > 90 then
diffazimuth = 180 - diffazimuth;
end if;
sss = ((scale2-1) * (diffazimuth / 90)^0.5);
point = st_transform(
st_project(
st_transform(point, 4326)::geography,
st_distance(midpoint, point) * sss, azimuth)::geometry,
st_srid(midpoint)
);
bend = st_setpoint(bend, i-1, point);
end loop;
size = wm_adjsize(bend);
end loop;
end $$ language plpgsql;
-- wm_exaggerate_bend exaggerates a given bend. Uses naive linear
-- interpolation. Faster than wm_exaggerate_bend2, but result visually looks
-- worse.
drop function if exists wm_exaggerate_bend;
create function wm_exaggerate_bend(
INOUT bend geometry,
size float,
desired_size float
) as $$
declare
scale constant float default 1.2; -- exaggeration enthusiasm
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);
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;
create function wm_adjsize(bend geometry, OUT adjsize float) as $$
declare
polygon geometry;
area float;
cmp float;
begin
adjsize = 0;
polygon = st_makepolygon(st_addpoint(bend, st_startpoint(bend)));
-- 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/(4pi)
-- 4. divide P by A: cmp = P/A = P/(u^2/(4pi)) = 4pi*P/u^2
area = st_area(polygon);
cmp = 4*pi()*area/(st_perimeter(polygon)^2);
if cmp > 0 then
adjsize = (area*(0.75/cmp));
end if;
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[],
dhalfcircle float,
intersect_patience integer,
dbgname text default null,
dbggen integer default null,
OUT mutated boolean
) as $$
declare
desired_size constant float default pi()*(dhalfcircle^2)/8;
bend geometry;
tmpint geometry;
i integer;
n integer;
last_id integer;
begin
mutated = false;
<<bendloop>>
for i in 1..array_length(attrs, 1) loop
if attrs[i].isolated and attrs[i].adjsize < desired_size then
bend = wm_exaggerate_bend2(bends[i], attrs[i].adjsize, desired_size);
-- 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(bends[i-1]);
if n > 3 then
continue when st_intersects(bend,
st_removepoint(st_removepoint(bends[i-1], n-1), n-2));
end if;
if n > 2 then
tmpint = st_intersection(bend, st_removepoint(bends[i-1], n-1));
continue when st_npoints(tmpint) > 1;
end if;
n = st_npoints(bends[i+1]);
if n > 3 then
continue when st_intersects(bend,
st_removepoint(st_removepoint(bends[i+1], 0), 0));
end if;
if n > 2 then
tmpint = st_intersection(bend, st_removepoint(bends[i+1], 0));
continue when st_npoints(tmpint) > 1;
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(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 = bends[i+n];
if st_npoints(tmpint) > 2 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(bends[i+n-1]) = 3 then
tmpint = st_removepoint(tmpint, 0);
end if;
end if;
continue bendloop when st_intersects(bend, tmpint);
end loop;
-- No intersections within intersect_patience, mutate bend!
mutated = true;
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.
bends[i-1] = st_addpoint(
st_removepoint(bends[i-1], st_npoints(bends[i-1])-1),
st_pointn(bends[i], 1),
-1
);
bends[i+1] = st_addpoint(
st_removepoint(bends[i+1], 0),
st_pointn(bends[i], st_npoints(bends[i])-1),
0
);
if dbgname is not null then
insert into wm_debug (stage, name, gen, nbend, way) values(
'gexaggeration', dbgname, dbggen, i, bends[i]);
end if;
end if;
end loop;
end $$ language plpgsql;
create function wm_elimination(
INOUT bends geometry[],
attrs wm_t_attrs[],
dhalfcircle float,
dbgname text default null,
dbggen integer default null,
OUT mutated boolean
) as $$
declare
desired_size constant float default pi()*(dhalfcircle^2)/8;
leftsize float;
rightsize float;
i int4;
begin
mutated = false;
i = 1;
while i < array_length(attrs, 1)-1 loop
i = i + 1;
continue when attrs[i].adjsize = 0;
continue when attrs[i].adjsize > desired_size;
if i = 2 then
leftsize = attrs[i].adjsize + 1;
else
leftsize = attrs[i-1].adjsize;
end if;
if i = array_length(attrs, 1)-1 then
rightsize = attrs[i].adjsize + 1;
else
rightsize = attrs[i+1].adjsize;
end if;
continue when attrs[i].adjsize >= leftsize;
continue when attrs[i].adjsize >= rightsize;
-- Local minimum. Elminate bend!
mutated = true;
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
bends[i-1] = st_addpoint(
st_removepoint(bends[i-1], st_npoints(bends[i-1])-1),
st_pointn(bends[i], 1),
-1
);
bends[i+1] = st_addpoint(
st_removepoint(bends[i+1], 0),
st_pointn(bends[i], st_npoints(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
insert into wm_debug(stage, name, gen, nbend, way) values(
'helimination',
dbgname,
dbggen,
generate_subscripts(bends, 1),
unnest(bends)
);
end if;
end $$ language plpgsql;
drop function if exists ST_SimplifyWM_Estimate;
create function ST_SimplifyWM_Estimate(
geom geometry,
OUT npoints bigint,
OUT secs bigint
) as $$
declare
lines geometry[];
l_type text;
begin
l_type = st_geometrytype(geom);
if l_type = 'ST_LineString' then
lines = array[geom];
elseif l_type = 'ST_MultiLineString' then
lines = array((select a.geom from st_dump(geom) a order by path[1] asc));
else
raise 'Unknown geometry type %', l_type;
end if;
npoints = 0;
for i in 1..array_length(lines, 1) loop
npoints = npoints + st_numpoints(lines[i]);
end loop;
secs = npoints / 33;
end $$ language plpgsql;
-- ST_SimplifyWM simplifies a given geometry using Wang & Müller's
-- "Line Generalization Based on Analysis of Shape Characteristics" algorithm,
-- 1998.
-- Input parameters:
-- - geom: ST_LineString or ST_MultiLineString: the geometry to be simplified
-- - dhalfcircle: the diameter of a half-circle, whose area is an approximate
-- threshold for small bend elimination. If bend's area is larger than that,
-- the bend will be left alone.
drop function if exists ST_SimplifyWM;
create function ST_SimplifyWM(
geom geometry,
dhalfcircle float,
intersect_patience integer default 10,
dbgname text default null
) returns geometry as $$
declare
gen integer;
i integer;
j integer;
line geometry;
lines geometry[];
bends geometry[];
attrs wm_t_attrs[];
mutated boolean;
l_type text;
begin
if intersect_patience is null then
intersect_patience = 10;
end if;
l_type = st_geometrytype(geom);
if l_type = 'ST_LineString' then
lines = array[geom];
elseif l_type = 'ST_MultiLineString' then
lines = array((select a.geom from st_dump(geom) a order by path[1] asc));
else
raise 'Unknown geometry type %', l_type;
end if;
<<lineloop>>
for i in 1..array_length(lines, 1) loop
mutated = true;
gen = 1;
while mutated loop
if dbgname is not null then
insert into wm_debug (stage, name, gen, nbend, way) values(
'afigures', dbgname, gen, i, lines[i]);
end if;
bends = wm_detect_bends(lines[i], dbgname, gen);
bends = wm_fix_gentle_inflections(bends, dbgname, gen);
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;
end if;
-- TODO: wm_combination
if not mutated then
select * from wm_elimination(bends, attrs,
dhalfcircle, dbgname, gen) into bends, mutated;
end if;
if mutated then
lines[i] = st_linemerge(st_union(bends));
if st_geometrytype(lines[i]) != 'ST_LineString' then
-- 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.',
gen, st_geometrytype(lines[i]), dbgname;
--exit lineloop;
end if;
gen = gen + 1;
continue;
end if;
end loop;
end loop;
if l_type = 'ST_LineString' then
return st_linemerge(st_union(lines));
elseif l_type = 'ST_MultiLineString' then
return st_union(lines);
end if;
end $$ language plpgsql;