/*
* Copyright (C) 2009 Christian Gawron
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License version 2 as
* published by the Free Software Foundation.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
*
* Author: Christian Gawron
* Create date: 03-Jul-2009
*/
package uk.me.parabola.mkgmap.reader.dem;
import java.io.File;
import java.io.FileNotFoundException;
import java.lang.reflect.Constructor;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import uk.me.parabola.imgfmt.ExitException;
import uk.me.parabola.imgfmt.FileSystemParam;
import uk.me.parabola.imgfmt.FormatException;
import uk.me.parabola.imgfmt.Utils;
import uk.me.parabola.imgfmt.app.Area;
import uk.me.parabola.imgfmt.app.Coord;
import uk.me.parabola.imgfmt.app.map.Map;
import uk.me.parabola.log.Logger;
import uk.me.parabola.mkgmap.build.MapBuilder;
import uk.me.parabola.mkgmap.general.LevelInfo;
import uk.me.parabola.mkgmap.general.LoadableMapDataSource;
import uk.me.parabola.mkgmap.osmstyle.StyleImpl;
import uk.me.parabola.mkgmap.osmstyle.StyledConverter;
import uk.me.parabola.mkgmap.reader.MapperBasedMapDataSource;
import uk.me.parabola.mkgmap.reader.osm.OsmConverter;
import uk.me.parabola.mkgmap.reader.osm.Style;
import uk.me.parabola.mkgmap.reader.osm.Way;
import uk.me.parabola.mkgmap.srt.SrtTextReader;
import uk.me.parabola.util.EnhancedProperties;
/**
* Create contour lines using an algorithm similar to that described in <a
* href="http://mapcontext.com/autocarto/proceedings/auto-carto-5/pdf/an-adaptive-grid-contouring-algorithm.pdf">An
* Adaptive Grid Contouring Algorithm</a> by Downing and Zoraster.
*/
public abstract class DEM
{
private static final Logger log =
Logger.
getLogger(DEM.
class);
private final static double epsilon = 1e-9
;
final protected static double delta =
1.5;
private final static int maxPoints =
200000;
private final static double minDist =
15;
private final static double maxDist =
21;
protected static int M =
1200;
protected static int N = M
;
protected static double res =
1.0 / N
;
private static int id = -
1;
protected int lat
;
protected int lon
;
protected abstract double ele
(int x,
int y
);
protected abstract void read
(int minLon,
int minLat,
int maxLon,
int maxLat
);
public static void createContours
(LoadableMapDataSource mapData, EnhancedProperties config
) {
Area bounds = mapData.
getBounds();
double minLat = Utils.
toDegrees(bounds.
getMinLat());
double minLon = Utils.
toDegrees(bounds.
getMinLong());
double maxLat = Utils.
toDegrees(bounds.
getMaxLat());
double maxLon = Utils.
toDegrees(bounds.
getMaxLong());
System.
out.
printf("bounds: %f %f %f %f\n", minLat, minLon, maxLat, maxLon
);
DEM data
;
String demType = config.
getProperty("dem-type",
"SRTM");
try {
String dataPath
;
Class<?> demClass
;
switch (demType
) {
case "ASTER":
dataPath = config.
getProperty("dem-path",
"ASTER");
demClass =
Class.
forName("uk.me.parabola.mkgmap.reader.dem.optional.GeoTiffDEM$ASTER");
break;
case "CGIAR":
dataPath = config.
getProperty("dem-path",
"CGIAR");
demClass =
Class.
forName("uk.me.parabola.mkgmap.reader.dem.optional.GeoTiffDEM$CGIAR");
break;
default:
dataPath = config.
getProperty("dem-path",
"SRTM");
demClass =
Class.
forName("uk.me.parabola.mkgmap.reader.dem.HGTDEM");
break;
}
Constructor<DEM
> constructor =
(Constructor<DEM
>) demClass.
getConstructor(String.
class,
Double.
TYPE,
Double.
TYPE,
Double.
TYPE,
Double.
TYPE);
data = constructor.
newInstance(dataPath, minLat, minLon, maxLat, maxLon
);
}
catch (Exception ex
) {
throw new ExitException
("failed to create DEM", ex
);
}
Isolines lines = data.
new Isolines
(data, minLat, minLon, maxLat, maxLon
);
int increment = config.
getProperty("dem-increment",
10);
double minHeight = lines.
getMinHeight();
double maxHeight = lines.
getMaxHeight();
int maxLevels = config.
getProperty("dem-maxlevels",
100);
while ((maxHeight - minHeight
) / increment
> maxLevels
)
increment
*=
2;
Style style = StyleImpl.
readStyle(config
);
LoadableMapDataSource dest = mapData
;
if (config.
getProperty("dem-separate-img",
false)) {
dest =
new DEMMapDataSource
(mapData, config
);
}
OsmConverter converter =
new StyledConverter
(style,
((MapperBasedMapDataSource
) dest
).
getMapper(), config
);
for (int level =
0; level
< maxHeight
; level += increment
) {
if (level
< minHeight
) continue;
// create isolines
lines.
addLevel(level
);
for (Isolines.
Isoline line : lines.
isolines) {
Way way =
new Way
(id--, line.
points);
way.
addTag("contour",
"elevation");
way.
addTag("ele",
String.
format("%d",
(int) line.
level));
converter.
convertWay(way
);
}
lines.
isolines.
clear();
}
if (config.
getProperty("dem-separate-img",
false)) {
MapBuilder builder =
new MapBuilder
();
builder.
config(config
);
// Get output directory
String DEFAULT_DIR =
".";
String fileOutputDir = config.
getProperty("output-dir", DEFAULT_DIR
);
// Test if directory exists
File outputDir =
new File(fileOutputDir
);
if (!outputDir.
exists()) {
System.
out.
println("Output directory not found. Creating directory '" + fileOutputDir +
"'");
if (!outputDir.
mkdirs()) {
System.
err.
println("Unable to create output directory! Using default directory instead");
fileOutputDir = DEFAULT_DIR
;
}
} else if (!outputDir.
isDirectory()) {
System.
err.
println("The --output-dir parameter must specify a directory. The parameter is being ignored, writing to default directory instead.");
fileOutputDir = DEFAULT_DIR
;
}
FileSystemParam params =
new FileSystemParam
();
params.
setMapDescription("contour lines");
long mapName =
Integer.
valueOf(config.
getProperty("mapname",
"63240000"));
try {
String mapname =
String.
format("%08d", mapName +
10000000);
Map map =
Map.
createMap(mapname, fileOutputDir, params, mapname, SrtTextReader.
sortForCodepage(1252));
builder.
makeMap(map, dest
);
map.
close();
}
catch (Exception ex
) {
throw new ExitException
("could not open " + mapName, ex
);
}
}
}
private int lastXi = -
1;
private int lastYi = -
1;
private final static int[][] bcInv =
{
{1,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0},
{0,
0,
0,
0,
0,
0,
0,
0,
1,
0,
0,
0,
0,
0,
0,
0},
{-
3,
0,
0,
3,
0,
0,
0,
0, -
2,
0,
0, -
1,
0,
0,
0,
0},
{2,
0,
0, -
2,
0,
0,
0,
0,
1,
0,
0,
1,
0,
0,
0,
0},
{0,
0,
0,
0,
1,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0},
{0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
1,
0,
0,
0},
{0,
0,
0,
0, -
3,
0,
0,
3,
0,
0,
0,
0, -
2,
0,
0, -
1},
{0,
0,
0,
0,
2,
0,
0, -
2,
0,
0,
0,
0,
1,
0,
0,
1},
{-
3,
3,
0,
0, -
2, -
1,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0},
{0,
0,
0,
0,
0,
0,
0,
0, -
3,
3,
0,
0, -
2, -
1,
0,
0},
{9, -
9,
9, -
9,
6,
3, -
3, -
6,
6, -
6, -
3,
3,
4,
2,
1,
2},
{-
6,
6, -
6,
6, -
4, -
2,
2,
4, -
3,
3,
3, -
3, -
2, -
1, -
1, -
2},
{2, -
2,
0,
0,
1,
1,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0},
{0,
0,
0,
0,
0,
0,
0,
0,
2, -
2,
0,
0,
1,
1,
0,
0},
{-
6,
6, -
6,
6, -
3, -
3,
3,
3, -
4,
4,
2, -
2, -
2, -
2, -
1, -
1},
{4, -
4,
4, -
4,
2,
2, -
2, -
2,
2, -
2, -
2,
2,
1,
1,
1,
1}
};
private static int lastId =
1000000000;
private static double lastX
;
private static double lastY
;
private static final int[][] off0 =
{{0,
0},
{0,
0},
{0,
1},
{1,
1}};
private static final int[][] off1 =
{{1,
0},
{0,
1},
{1,
1},
{1,
0}};
private static final int[] brd =
{1,
2,
4,
8};
private static final int[] rev =
{2,
3,
0,
1};
private static final int[][] mov =
{{0, -
1},
{-
1,
0},
{0,
1},
{1,
0}};
private final double[][] bc =
new double[4][4];
private final double[] bc_y =
new double[4];
private final double[] bc_y1 =
new double[4];
private final double[] bc_y2 =
new double[4];
private final double[] bc_y12 =
new double[4];
private final double[] bc_Coeff =
new double[16];
private final double[] bc_x =
new double[16];
private void recalculateCoefficients
(int xi,
int yi
) {
double v00 = ele
(xi, yi
);
double v0p = ele
(xi, yi +
1);
double vpp = ele
(xi +
1, yi +
1);
double vp0 = ele
(xi +
1, yi
);
double vm0 = ele
(xi -
1, yi
);
double v0m = ele
(xi, yi -
1);
double vmp = ele
(xi -
1, yi +
1);
double vpm = ele
(xi +
1, yi -
1);
double vmm = ele
(xi -
1, yi -
1);
double vmP = ele
(xi +
2, yi -
1);
double vPm = ele
(xi -
1, yi +
2);
double vP0 = ele
(xi +
2, yi
);
double v0P = ele
(xi, yi +
2);
double vPp = ele
(xi +
2, yi +
1);
double vpP = ele
(xi +
1, yi +
2);
double vPP = ele
(xi +
2, yi +
2);
bc_y
[0] = v00
;
bc_y
[1] = vp0
;
bc_y
[2] = vpp
;
bc_y
[3] = v0p
;
bc_y1
[0] =
(vp0 - vm0
) /
2;
bc_y1
[1] =
(vP0 - v00
) /
2;
bc_y1
[2] =
(vPp - v0p
) /
2;
bc_y1
[3] =
(vpp - vmp
) /
2;
bc_y2
[0] =
(v0p - v0m
) /
2;
bc_y2
[1] =
(vpp - vpm
) /
2;
bc_y2
[2] =
(vpP - vp0
) /
2;
bc_y2
[3] =
(v0P - v00
) /
2;
bc_y12
[0] =
(vpp - vpm - vmp + vmm
) /
4;
bc_y12
[0] =
(vPp - vPm - v0p + v0m
) /
4;
bc_y12
[2] =
(vPP - vP0 - v0P + v00
) /
4;
bc_y12
[0] =
(vpP - vp0 - vmP + vm0
) /
4;
int i
;
for (i =
0; i
< 4; i++
) {
bc_x
[i
] = bc_y
[i
];
bc_x
[i +
4] = bc_y1
[i
];
bc_x
[i +
8] = bc_y2
[i
];
bc_x
[i +
12] = bc_y12
[i
];
}
for (i =
0; i
< 16; i++
) {
double s =
0;
for (int k =
0; k
< 16; k++
) s += bcInv
[i
][k
] * bc_x
[k
];
bc_Coeff
[i
] = s
;
}
int l =
0;
for (i =
0; i
< 4; i++
)
for (int j =
0; j
< 4; j++
)
bc
[i
][j
] = bc_Coeff
[l++
];
}
protected double gradient
(double lat,
double lon,
double[] grad
) {
grad
[0] =
0;
grad
[1] =
0;
double x =
(lon -
this.
lon) / res
;
double y =
(lat -
this.
lat) / res
;
int xi =
(int) x
;
int yi =
(int) y
;
if (lastXi
!= xi || lastYi
!= yi
) {
log.
debug("new Cell for interpolation: %d %d", xi, yi
);
recalculateCoefficients
(xi, yi
);
lastXi = xi
;
lastYi = yi
;
}
double t = x - xi
;
double u = y - yi
;
if (xi
< 0 || xi
> N +
1 || yi
< 0 || yi
> N +
1)
throw new IndexOutOfBoundsException(String.
format("(%f, %f)->(%d, %d)", lat, lon, xi, yi
));
double val =
0;
for (int i =
3; i
>=
0; i--
) {
val = t
* val +
((bc
[i
][3] * u + bc
[i
][2]) * u + bc
[i
][1]) * u + bc
[i
][0];
grad
[0] = u
* grad
[0] +
(3 * bc
[3][i
] * t +
2 * bc
[2][i
]) * t + bc
[1][i
];
grad
[1] = t
* grad
[1] +
(3 * bc
[i
][3] * t +
2 * bc
[i
][2]) * t + bc
[i
][1];
}
return val
;
}
protected double elevation
(double lat,
double lon
) {
double x =
(lon -
this.
lon) / res
;
double y =
(lat -
this.
lat) / res
;
int xi =
(int) x
;
int yi =
(int) y
;
if (lastXi
!= xi || lastYi
!= yi
) {
log.
debug("new Cell for interpolation: %d %d", xi, yi
);
recalculateCoefficients
(xi, yi
);
lastXi = xi
;
lastYi = yi
;
}
double t = x - xi
;
double u = y - yi
;
if (xi
< 0 || xi
> N +
1 || yi
< 0 || yi
> N +
1)
throw new IndexOutOfBoundsException(String.
format("(%f, %f)->(%d, %d)", lat, lon, xi, yi
));
double val =
0;
for (int i =
3; i
>=
0; i--
) {
val = t
* val +
((bc
[i
][3] * u + bc
[i
][2]) * u + bc
[i
][1]) * u + bc
[i
][0];
}
return val
;
}
protected double elevation
(int x,
int y
) {
if (x
< 0 || x
> N || y
< 0 || y
> N
)
throw new IndexOutOfBoundsException(String.
format("elevation: %d %d", x, y
));
return ele
(x, y
);
}
class Isolines
{
final DEM data
;
final int minX
;
final int maxX
;
final int minY
;
final int maxY
;
double min
;
double max
;
final ArrayList<Isoline
> isolines =
new ArrayList<>();
class Isoline
{
final int id
;
final ArrayList<Coord
> points
;
final double level
;
private Isoline
(double level
) {
this.
level = level
;
id = lastId++
;
points =
new ArrayList<>();
}
private class Edge
implements Brent.
Function {
final double x0
;
final double y0
;
final double x1
;
final double y1
;
Edge
(double x0,
double y0,
double x1,
double y1
) {
this.
x0 = x0
;
this.
y0 = y0
;
this.
x1 = x1
;
this.
y1 = y1
;
}
public double eval
(double d
) {
return data.
elevation(x0 + d
* (x1 - x0
), y0 + d
* (y1 - y0
)) - level
;
}
}
private class FN
implements Brent.
Function {
double x0, y0
;
double dx, dy
;
public void setParameter
(double x0,
double y0,
double dx,
double dy
) {
this.
x0 = x0
;
this.
y0 = y0
;
this.
dx = dx
;
this.
dy = dy
;
}
public double eval
(double t
) {
return data.
elevation(y0 + t
* dy, x0 + t
* dx
) - level
;
}
}
private final FN fn =
new FN
();
final double[] grad =
new double[2];
final double[] px =
new double[4];
final double[] py =
new double[4];
final int[] edges =
new int[4];
boolean addCell
(Position p,
int direction
) {
log.
debug("addCell: %f %d %d %d %d", level, p.
ix, p.
iy, p.
edge, direction
);
int c =
0;
for (int k =
0; k
< 4; k++
) {
if (k == p.
edge)
continue;
int x0 = p.
ix + off0
[k
][0];
int y0 = p.
iy + off0
[k
][1];
int x1 = p.
ix + off1
[k
][0];
int y1 = p.
iy + off1
[k
][1];
double l0 = elevation
(x0, y0
) - level
;
double l1 = elevation
(x1, y1
) - level
;
if (Math.
abs(l1
) < epsilon || l0
* l1
< 0) {
edges
[c
] = k
;
Brent.
Function f =
new Edge
(data.
lat + y0
* DEM.
res, data.
lon + x0
* DEM.
res, data.
lat + y1
* DEM.
res, data.
lon + x1
* DEM.
res);
double f0 = elevation
(x0, y0
) - level
;
double delta
;
if (Math.
abs(1) < epsilon
) {
delta =
1;
} else if (Math.
abs(f0
) < epsilon
)
throw new ExitException
("implementation error!");
else
delta = Brent.
zero(f, epsilon,
1 - epsilon
);
px
[c
] = data.
lon +
(x0 + delta
* (x1 - x0
)) * DEM.
res;
py
[c
] = data.
lat +
(y0 + delta
* (y1 - y0
)) * DEM.
res;
c++
;
}
}
if (c ==
1) {
p.
edge = edges
[0];
double px0 = p.
x;
double py0 = p.
y;
p.
x = px
[0];
p.
y = py
[0];
double px1 = p.
x;
double py1 = p.
y;
double xMin = data.
lon + p.
ix * DEM.
res;
double xMax = xMin + DEM.
res;
double yMin = data.
lat + p.
iy * DEM.
res;
double yMax = yMin + DEM.
res;
refineAdaptively
(xMin, yMin, xMax, yMax, px0, py0, px1, py1, direction, maxDist
);
addPoint
(p.
x, p.
y, direction
);
p.
moveCell();
return true;
} else {
log.
debug("addCellByStepping: %d", c
);
return addCellByStepping
(p, direction, c, edges, px, py
);
}
}
private void refineAdaptively
(double xMin,
double yMin,
double xMax,
double yMax,
double x0,
double y0,
double x1,
double y1,
int direction,
double maxDist
)
{
double dist = quickDistance
(x0, y0, x1, y1
);
if (dist
> maxDist
) {
double dx = x1 - x0
;
double dy = y1 - y0
;
double xm = x0 +
0.5 * dx
;
double ym = y0 +
0.5 * dy
;
double n =
Math.
sqrt(dx
* dx + dy
* dy
);
fn.
setParameter(xm, ym, -dy / n, dx / n
);
Brent.
Function f = fn
;
double t0 = -
0.05 * res
;
double t1 =
0.05 * res
;
double f0 = f.
eval(t0
);
double f1 = f.
eval(t1
);
int count =
0;
while (f0
* f1
> 0 && count++
< 20) {
if ((count
& 1) > 0)
t0 -=
0.05 * res
;
else
t1 +=
0.05 * res
;
f0 = f.
eval(t0
);
f1 = f.
eval(t1
);
log.
debug("refine: %f %f %f %f", t0, t1, f0, f1
);
}
if (f0
* f1
< 0) {
double t = Brent.
zero(f, t0, t1
);
xm -= t
* dy
;
ym += t
* dx
;
} else {
log.
debug("refine failed: %f %f %f %f", t0, t1, f0, f1
);
return;
}
if (xm
> xMin
&& xm
< xMax
&& ym
> yMin
&& ym
< yMax
)
refineAdaptively
(xMin, yMin, xMax, yMax, x0, y0, xm, ym, direction, maxDist
* 1.1);
addPoint
(xm, ym, direction
);
if (xm
> xMin
&& xm
< xMax
&& ym
> yMin
&& ym
< yMax
)
refineAdaptively
(xMin, yMin, xMax, yMax, xm, ym, x1, y1, direction, maxDist
* 1.1);
}
}
boolean addCellByStepping
(Position p,
int direction,
int numEdges,
int[] edges,
double[] px,
double[] py
) {
log.
debug("addCellByStepping: %f %d %d %d %d", level, p.
ix, p.
iy, p.
edge, direction
);
int iMin = -
1;
double md =
5000;
for (int i =
0; i
< numEdges
; i++
) {
gradient
(p.
y, p.
x, grad
);
double dist = quickDistance
(p.
x, p.
y, px
[i
], py
[i
]);
log.
debug("distance %d: %f", i, dist
);
if (dist
< md
&& (visited
[p.
iy * (N +
1) + p.
ix] & brd
[edges
[i
]]) ==
0) {
md = dist
;
iMin = i
;
}
}
p.
edge = edges
[iMin
];
double px0 = p.
x;
double py0 = p.
y;
p.
x = px
[iMin
];
p.
y = py
[iMin
];
double px1 = p.
x;
double py1 = p.
y;
double xMin = data.
lon + p.
ix * DEM.
res;
double xMax = xMin + DEM.
res;
double yMin = data.
lat + p.
iy * DEM.
res;
double yMax = yMin + DEM.
res;
refineAdaptively
(xMin, yMin, xMax, yMax, px0, py0, px1, py1, direction, maxDist
);
addPoint
(p.
x, p.
y, direction
);
p.
moveCell();
return true;
}
private void addPoint
(double x,
double y,
int direction
) {
double dist = quickDistance
(x, y, lastX, lastY
);
log.
debug("addPoint: %f %f %f", x, y, dist
);
if (dist
> minDist
) {
if (direction
> 0)
points.
add(0,
new Coord
(y, x
));
else
points.
add(points.
size(),
new Coord
(y, x
));
lastX = x
;
lastY = y
;
}
}
private void close
() {
points.
add(points.
size(), points.
get(0));
}
}
public Isolines
(DEM data,
double minLat,
double minLon,
double maxLat,
double maxLon
) {
System.
out.
printf("init: %f %f %f %f\n", minLat, minLon, maxLat, maxLon
);
this.
data = data
;
this.
minX =
(int) ((minLon - data.
lon) / res
);
this.
minY =
(int) ((minLat - data.
lat) / res
);
this.
maxX =
(int) ((maxLon - data.
lon) / res
);
this.
maxY =
(int) ((maxLat - data.
lat) / res
);
init
();
}
private void init
() {
System.
out.
printf("init: %d %d %d %d\n", minX, minY, maxX, maxY
);
data.
read(minX -
2, minY -
2, maxX +
2, maxY +
2);
// we need some overlap for bicubic interpolation
max = -
1000;
min =
10000;
for (int i = minX
; i
< maxX
; i++
)
for (int j = minY
; j
< maxY
; j++
) {
if (data.
elevation(i, j
) < min
) min = data.
elevation(i, j
);
if (data.
elevation(i, j
) > max
) max = data.
elevation(i, j
);
}
log.
debug("min: %f, max: %f\n", min, max
);
}
double getMinHeight
() {
return min
;
}
double getMaxHeight
() {
return max
;
}
private class Edge
implements Brent.
Function {
final double x0, y0, x1, y1, level
;
Edge
(double x0,
double y0,
double x1,
double y1,
double level
) {
this.
x0 = x0
;
this.
y0 = y0
;
this.
x1 = x1
;
this.
y1 = y1
;
this.
level = level
;
}
public double eval
(double d
) {
return data.
elevation(x0 + d
* (x1 - x0
), y0 + d
* (y1 - y0
)) - level
;
}
}
private class Position {
int ix, iy
;
double x, y
;
int edge
;
Position(int ix,
int iy,
double x,
double y,
int edge
) {
this.
ix = ix
;
this.
iy = iy
;
this.
x = x
;
this.
y = y
;
this.
edge = edge
;
}
Position(Position p
) {
this.
ix = p.
ix;
this.
iy = p.
iy;
this.
x = p.
x;
this.
y = p.
y;
this.
edge = p.
edge;
}
void markEdge
() {
log.
debug("marking edge: %d %d %d %d", ix, iy, edge, brd
[edge
]);
visited
[iy
* (N +
1) + ix
] |= brd
[edge
];
}
void moveCell
() {
markEdge
();
ix += mov
[edge
][0];
iy += mov
[edge
][1];
edge = rev
[edge
];
markEdge
();
}
}
final byte[] visited =
new byte[(N +
1) * (N +
1)];
public void addLevel
(double level
) {
if (level
< min || level
> max
)
return;
System.
out.
printf("addLevel: %f\n", level
);
Arrays.
fill(visited,
(byte) 0);
for (int y = minY
; y
< maxY
; y++
) {
for (int x = minX
; x
< maxX
; x++
) {
byte v =
0;
int direction
;
// Mark the borders of the cell, represented by the four points (i, j), (i+1, j), (i, j+1), (i+1, j+1),
// which are intersected by the contour. The values are:
// 1: top
// 2: left
// 4: bottom
// 8: right
if (data.
elevation(x, y
) >= level
) {
if (data.
elevation(x +
1, y
) < level
) {
v |=
1;
}
if (data.
elevation(x, y +
1) < level
) {
v |=
2;
}
direction =
1;
} else {
if (data.
elevation(x +
1, y
) > level
) {
v |=
1;
}
if (data.
elevation(x, y +
1) > level
) {
v |=
2;
}
direction = -
1;
}
int k = -
1;
if ((v
& 1) > 0 && (visited
[y
* (N +
1) + x
] & 1) ==
0) {
k =
0;
} else if ((v
& 2) > 0 && (visited
[y
* (N +
1) + x
] & 2) ==
0) {
k =
1;
}
if (k
>=
0) {
int x0 = x + off0
[k
][0];
int y0 = y + off0
[k
][1];
int x1 = x + off1
[k
][0];
int y1 = y + off1
[k
][1];
try {
Brent.
Function f =
new Edge
(data.
lat + y0
* DEM.
res, data.
lon + x0
* DEM.
res,
data.
lat + y1
* DEM.
res, data.
lon + x1
* DEM.
res,
level
);
double f0 = elevation
(x0, y0
) - level
;
double f1 = elevation
(x1, y1
) - level
;
double delta
;
if (Math.
abs(f0
) < epsilon
) {
delta =
0;
} else if (Math.
abs(f1
) < epsilon
)
continue;
else
delta = Brent.
zero(f,
0,
1 - epsilon
);
Position p =
new Position(x, y, data.
lon +
(x0 + delta
* (x1 - x0
)) * DEM.
res, data.
lat +
(y0 + delta
* (y1 - y0
)) * DEM.
res, k
);
p.
markEdge();
isolines.
add(traceByStepping
(level, p, direction
));
}
catch (RuntimeException ex
) {
log.
debug("error: %s", ex.
toString());
ex.
printStackTrace();
return;
}
}
}
}
}
private Isoline traceByStepping
(double level,
Position p,
int direction
) {
log.
debug("traceByStepping: starting contour %f %d %d %f %f %d", level, p.
ix, p.
iy, p.
x, p.
y, p.
edge);
int n =
0;
Position startP =
new Position(p
);
boolean foundEnd =
false;
Isoline line =
new Isoline
(level
);
while (true) {
log.
debug("traceByStepping: %f %d %d %f %f %d", level, p.
ix, p.
iy, p.
x, p.
y, p.
edge);
visited
[p.
iy * (N +
1) + p.
ix] |= brd
[p.
edge];
if (n
> 0 && p.
ix == startP.
ix && p.
iy == startP.
iy && quickDistance
(p.
x, p.
y, startP.
x, startP.
y) < 5) {
log.
debug("closed curve!");
line.
close();
break;
} else if (p.
ix < minX || p.
iy < minY || p.
ix >= maxX || p.
iy >= maxY
) {
if (foundEnd
) // did we already reach one end?
{
log.
debug("second border reached!");
break;
} else {
log.
debug("border reached!");
foundEnd =
true;
n =
0;
direction
*= -
1;
p =
new Position(startP
);
p.
moveCell();
continue;
}
}
n++
;
if (!line.
addCell(p, direction
) || line.
points.
size() > maxPoints
) {
log.
debug("ending contour");
isolines.
add(line
);
return line
;
}
}
return line
;
}
}
private static double quickDistance
(double long1,
double lat1,
double long2,
double lat2
) {
double latDiff
;
if (lat1
< lat2
)
latDiff = lat2 - lat1
;
else
latDiff = lat1 - lat2
;
if (latDiff
> 90)
latDiff -=
180;
double longDiff
;
if (long1
< long2
)
longDiff = long2 - long1
;
else
longDiff = long1 - long2
;
if (longDiff
> 180)
longDiff -=
360;
// scale longDiff by cosine of average latitude
longDiff
*=
Math.
cos(Math.
PI /
180 * Math.
abs((lat1 + lat2
) /
2));
double distDegSq =
(latDiff
* latDiff
) +
(longDiff
* longDiff
);
return 40075000 * Math.
sqrt(distDegSq
) /
360;
}
private static class DEMMapDataSource
extends MapperBasedMapDataSource
implements LoadableMapDataSource
{
final LoadableMapDataSource parent
;
final List<String> copyright =
new ArrayList<>();
DEMMapDataSource
(LoadableMapDataSource parent, EnhancedProperties props
) {
this.
parent = parent
;
config
(props
);
}
public boolean isFileSupported
(String name
) {
return false;
}
public void load
(String name
)
throws FileNotFoundException, FormatException
{
throw new FormatException
("load not supported");
}
public LevelInfo
[] mapLevels
() {
return parent.
mapLevels();
}
public LevelInfo
[] overviewMapLevels
() {
return parent.
overviewMapLevels();
}
public String[] copyrightMessages
() {
return copyright.
toArray(new String[1]);
}
}
}