/*
* Copyright (C) 2010.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License version 3 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.
*/
package test.display;
import java.util.ArrayList;
import java.util.List;
import java.util.Locale;
import java.util.TreeMap;
import java.util.Map.Entry;
/**
* Stand alone program to display the DEM file. This stores a
* digital elevation model used to display shaded surface and profiles of tracks.
*
* @author Ronny Klier
* @author Gerd Petermann
*
*/
public class DemDisplay
extends CommonDisplay
{
int countErr =
0;
int countGood =
0;
int countSorry =
0;
private enum EncType
{
HYBRID, LEN
}
private enum WrapType
{
WRAP_0, WRAP_1, WRAP_2
}
private enum CalcType
{
CALC_STD, CALC_PLATEAU_ZERO, CALC_PLATEAU_NON_ZERO
}
private static final int[] plateauUnit =
{ 1,
1,
1,
1,
2,
2,
2,
2,
4,
4,
4,
4,
8,
8,
8,
8,
16,
16,
32,
32,
64,
64,
128 };
private static final int[] plateauBinBits =
{ 0,
0,
0,
1,
1,
1,
1,
2,
2,
2,
2,
3,
3,
3,
3,
4,
4,
5,
5,
6,
6,
7,
8 };
private static int evaluateData
(int oldsum,
int elemcount,
int newdata,
int region
) {
switch (region
) {
case 0:
return -
1 - oldsum - elemcount
;
case 1:
return 2 * (newdata + elemcount
) +
3;
case 2:
return 2 * newdata -
1;
case 3:
return 2 * (newdata - elemcount
) -
5;
default:
return 1 - oldsum + elemcount
;
}
}
private static int getEvaluateDataRegion
(int oldsum,
int elemcount,
int newdata
) {
if (elemcount
< 63) {
if (newdata
< -
2 -
((oldsum +
3 * elemcount
) >> 1)) {
return 0;
} else if (newdata
< -
((oldsum + elemcount
) >> 1)) {
return 1;
} else if (newdata
< 2 -
((oldsum - elemcount
) >> 1)) {
return 2;
} else if (newdata
< 4 -
((oldsum -
3 * elemcount
) >> 1)) {
return 3;
} else {
return 4;
}
} else {
if (newdata
< -
2 -
((oldsum +
3 * elemcount
) >> 1)) {
return 0;
} else if (newdata
< -
((oldsum + elemcount
) >> 1) -
1) {
// special case in if !
return 1;
} else if (newdata
< 2 -
((oldsum - elemcount
) >> 1)) {
return 2;
} else if (newdata
< 4 -
((oldsum -
3 * elemcount
) >> 1)) {
return 3;
} else {
return 4;
}
}
}
private static int getStartHUnit
(int maxHeight
) {
if (maxHeight
< 0x9f
)
return 1;
else if (maxHeight
< 0x11f
)
return 2;
else if (maxHeight
< 0x21f
)
return 4;
else if (maxHeight
< 0x41f
)
return 8;
else if (maxHeight
< 0x81f
)
return 16;
else if (maxHeight
< 0x101f
)
return 32;
else if (maxHeight
< 0x201f
)
return 64;
else if (maxHeight
< 0x401f
)
return 128;
return 256;
}
/**
*
* @param hu
* @return
*/
private static int normalizeHUnit
(int hu
) {
if (hu
> 0) {
return Integer.
highestOneBit(hu
);
}
return 0;
}
private class ValDecoder
{
private EncType encType
;
private WrapType wrapType
;
private int sumH
;
private int sumL
;
private int elemCount
;
private int hunit
;
private final CalcType type
;
private final int unitDelta
;
// used for debugging
private final int maxHeight
;
private int followerDDiff
;
private int shrinkFactor
;
public ValDecoder
(CalcType type,
int maxHeight,
int shrinkFactor
) {
this.
type = type
;
this.
maxHeight = maxHeight
;
this.
shrinkFactor = shrinkFactor
;
unitDelta =
Math.
max(0, maxHeight - 0x5f
) / 0x40
; // TODO what does it mean?
encType = EncType.
HYBRID;
wrapType = WrapType.
WRAP_0;
hunit = getStartHUnit
(maxHeight
);
}
public EncType getEncType
() {
return encType
;
}
private int decode1
(int rawVal
) {
int v
;
if (wrapType == WrapType.
WRAP_0)
v = rawVal
;
else if (wrapType == WrapType.
WRAP_1)
v =
1 - rawVal
;
else {// if (wrapType == WrapType.WRAP02)
v = -rawVal
;
}
return v
;
}
private int processVal
(int delta1
) {
if (type == CalcType.
CALC_STD) {
// calculate threshold sum hybrid
sumH += delta1
> 0 ? delta1 : -delta1
;
if (sumH + unitDelta +
1 >= 0xffff
)
sumH -= 0x10000
;
// calculate threshold sum for length encoding
int evalRegion = -
1;
int workData = delta1
;
if (elemCount ==
63) {
evalRegion = getEvaluateDataRegion
(sumL, elemCount, delta1
);
boolean datagerade = delta1
% 2 ==
0;
boolean sumL1 =
(sumL -
1) % 4 ==
0;
switch (evalRegion
) {
case 0:
case 2:
case 4:
if ((sumL1
&& !datagerade
) ||
(!sumL1
&& datagerade
)) {
workData++
;
}
break;
case 1:
workData++
;
if ((sumL1
&& !datagerade
) ||
(!sumL1
&& datagerade
)) {
workData++
;
}
break;
case 3:
if ((sumL1
&& datagerade
) ||
(!sumL1
&& !datagerade
)) {
workData--
;
}
break;
}
}
if (evalRegion
< 0)
evalRegion = getEvaluateDataRegion
(sumL, elemCount, workData
);
int eval = evaluateData
(sumL, elemCount, workData, evalRegion
);
sumL += eval
;
// now update elem counter
elemCount++
;
if (elemCount ==
64) {
elemCount =
32;
sumH =
((sumH - unitDelta
) >> 1) -
1;
sumL /=
2;
if (sumL
% 2 !=
0) {
sumL++
;
}
}
// calculate new hunit
hunit = normalizeHUnit
((unitDelta + sumH +
1) /
(elemCount +
1));
// finally determine encode type for next value
wrapType = WrapType.
WRAP_0;
if (hunit
> 0) {
encType = EncType.
HYBRID;
} else {
encType = EncType.
LEN;
if (sumL
> 0 && shrinkFactor
<=
1) {
wrapType = WrapType.
WRAP_1;
}
}
} else if (type == CalcType.
CALC_PLATEAU_ZERO) {
// calculate threshold sum hybrid
sumH += delta1
> 0 ? delta1 :
1 - delta1
; // different to standard
if (sumH + unitDelta +
1 >= 0xffff
)
sumH -= 0x10000
;
// calculate threshold sum for length encoding
sumL += delta1
<=
0 ? -
1 :
1;
elemCount++
;
if (elemCount ==
64) {
elemCount =
32;
sumH =
((sumH - unitDelta
) >> 1) -
1;
sumL /=
2;
if (sumL
% 2 !=
0) {
sumL++
;
}
}
// calculate new hunit
hunit = normalizeHUnit
((unitDelta + sumH +
1 - elemCount /
2) /
(elemCount +
1));
// finally determine encode type for next value
wrapType = WrapType.
WRAP_0;
if (hunit
> 0) {
encType = EncType.
HYBRID;
} else {
encType = EncType.
LEN;
if (sumL
>=
0)
wrapType = WrapType.
WRAP_1;
}
if (delta1
<=
0)
--delta1
;
} else {
assert type == CalcType.
CALC_PLATEAU_NON_ZERO;
// calculate threshold sum hybrid
sumH += delta1
< 0 ? -delta1 : delta1
; // simple absolute sum
if (sumH + unitDelta +
1 >= 0xffff
)
sumH -= 0x10000
;
// calculate threshold sum for length encoding
sumL += delta1
<=
0 ? -
1 :
1;
elemCount++
;
if (elemCount ==
64) {
elemCount =
32;
sumH =
((sumH - unitDelta
) >> 1) -
1;
sumL /=
2;
if (sumL
% 2 !=
0) {
sumL--
; // different to CALC_PLATEAU_ZERO !
}
}
// calculate new hunit
hunit = normalizeHUnit
((unitDelta + sumH +
1) /
(elemCount +
1));
// finally determine encode type for next value
wrapType = WrapType.
WRAP_0;
if (hunit
> 0) {
encType = EncType.
HYBRID;
} else {
encType = EncType.
LEN;
if (sumL
<=
0)
wrapType = WrapType.
WRAP_2;
}
if (followerDDiff
> 0) {
delta1 = -delta1
;
}
}
return delta1
;
}
public int getHUnit
() {
return hunit
;
}
public void setDDiff
(int followerDDiff
) {
this.
followerDDiff = followerDDiff
;
}
}
private class DemTile
{
private DemSection section
;
private int tileNumberLat
;
private int tileNumberLon
;
private int height
;
private int width
;
private int offset
; // offset from section.dataOffset2
private int dataLength
; // number of bytes for compressed elevation data
private int baseHeight
; // base or minimum height in this tile
private int maxDeltaHeight
; // maximum height above baseHeight (elevation?)
private int encDeltaHeight
;
private byte encodingType
; // seems to determine how the highest values are displayed
private int bitPos
;
private int bytePos
;
private int lastTablePos =
0; // current position in plateau tables
private StringBuilder bitString =
new StringBuilder();
DemTile
(DemSection parent,
int numLat,
int numLon
) {
section = parent
;
tileNumberLat = numLat
;
tileNumberLon = numLon
;
width =
(numLon +
1 == parent.
tilesLon) ? parent.
nonStdWidth : parent.
pointsPerLon;
height =
(numLat +
1 == parent.
tilesLat) ? parent.
nonStdHeight : parent.
pointsPerLat;
}
void readHeader
(int baseOffset
) {
Displayer d =
new Displayer
(reader
);
d.
setTitle("DEM Tile Header");
DisplayItem item = d.
item();
item.
addText("Zoom level %1$s Tile row %2$s column %3$s", section.
zoomLevel, tileNumberLat, tileNumberLon
);
reader.
position(section.
dataOffset +
(section.
tilesLon * tileNumberLat + tileNumberLon
) * section.
tileDescSize);
switch (section.
offsetSize) {
case 1:
offset = d.
byteValue("data section: at offset %#08x") & 0xff
;
break;
case 2:
offset = d.
charValue("data section: at offset %#08x") & 0xffff
;
break;
case 3:
offset = d.
int3Value("data section: at offset %#08x") & 0xffffff
;
break;
case 4:
offset = d.
intValue("data section: at offset %#08x");
break;
}
assert offset
>=
0;
offset += baseOffset
;
if (1 == section.
baseSize) {
baseHeight = d.
sByteValue("minimum (base) height: %d");
} else {
baseHeight = d.
shortValue("minimum (base) height: %d");
}
if (1 == section.
differenceSize) {
maxDeltaHeight = d.
byteValue("max height delta: %d") & 0xff
;
} else {
maxDeltaHeight = d.
shortValue("max height delta: %d");
}
if (section.
hasExtra) {
encodingType = d.
byteValue("encoding type: %d");
}
d.
print(outStream
);
assert maxDeltaHeight
>=
0;
return;
}
void setDataLength
(int nextOffset
) {
if (maxDeltaHeight ==
0)
return;
dataLength = nextOffset - offset
;
if (dataLength
< 0)
dataLength =
0;
}
public void decodeBitStreamData
(int lonOff,
int latOff
) {
// if (section.shrinkFactor > 1)
for (int loop =
0; loop
< 2; loop++
) {
try {
decodeBitStreamData0
(lonOff, latOff
);
countGood++
;
break;
} catch (Exception e
) {
System.
out.
println(e
);
System.
out.
println("Error: shrink=" + section.
shrinkFactor +
", maxDeltaHeight=" + maxDeltaHeight
+
", width=" + width +
", height=" + height
);
if (loop ==
0)
countErr++
;
}
}
}
public void decodeBitStreamData0
(int lonOff,
int latOff
) {
reader.
position(offset
);
Displayer d =
new Displayer
(reader
);
d.
setTitle(String.
format("Bitstream data for tile row %1$s column %2$s at zoom level %3$s", tileNumberLat, tileNumberLon, section.
zoomLevel));
d.
rawValue(dataLength,
"DEM bitstream data");
boolean tryRead =
true;
if (dataLength ==
0) {
DisplayItem item = d.
item();
item.
addText("no data");
tryRead =
false;
}
encDeltaHeight = maxDeltaHeight / section.
shrinkFactor;
int rest = maxDeltaHeight
% section.
shrinkFactor;
d.
print(outStream
);
if (!tryRead
)
return;
bytePos =
0;
bitPos =
0;
lastTablePos =
0;
bitString.
setLength(0);
reader.
position(offset
);
byte[] bytes =
new byte[dataLength
];
for (int i =
0; i
< dataLength
; i ++
) {
bytes
[i
] = reader.
get();
}
int maxZeroBits= getMaxLengthZeroBits
(section.
shrinkFactor, maxDeltaHeight
);
// assert 22 == getMaxLengthZeroBits(3, 117);
if (rest
> 0)
encDeltaHeight++
;
System.
out.
println(String.
format("mdh=%d mh=%d sf=%d rest=%d bh=%d maxZero=%d bsLen=%d",maxDeltaHeight, encDeltaHeight, section.
shrinkFactor,rest, baseHeight, maxZeroBits, dataLength
));
if (section.
shrinkFactor > 1 && rest
!=
0) {
System.
out.
println("sorry, don't know how to interpret data with shrink factor > 1 and rest > 0");
++countSorry
;
return;
}
// silent = (maxDeltaHeight == 211 && section.shrinkFactor == 3) == false;
// silent = (section.shrinkFactor == 3) == false;
boolean readStd =
false;
ValDecoder decoderStd =
new ValDecoder
(CalcType.
CALC_STD, encDeltaHeight, section.
shrinkFactor);
ValDecoder decoderPF0 =
new ValDecoder
(CalcType.
CALC_PLATEAU_ZERO, encDeltaHeight, section.
shrinkFactor);
ValDecoder decoderPFN0 =
new ValDecoder
(CalcType.
CALC_PLATEAU_NON_ZERO, encDeltaHeight, section.
shrinkFactor);
int n =
0;
int m =
0;
int numVals = height
* width
;
List<Integer> heights =
new ArrayList<>(numVals
);
// we always start with a plateau
int pLen =
0;
CalcType ct = CalcType.
CALC_PLATEAU_ZERO;
ValDecoder currDecoder = decoderPF0
;
while (true) {
if (readStd
)
currDecoder = decoderStd
;
else {
pLen = decodePlateauLen
(n, bytes
);
if (!silent
)
reportBits
("plateau len " + pLen
);
if (pLen
> 0) {
int v = getHeight
(n-
1, m, heights
);
for (int i =
0; i
< pLen
; i++
) {
addHeight
(v, heights
);
check
(heights, section.
zoomLevel, latOff, lonOff
);
}
n = heights.
size();
m = n / width
;
n = n - width
* m
;
if (m
>= height
)
break;
if (n ==
0)
continue;
}
int followerDDiff =
0;
if (n
> 0) {
followerDDiff = getHeight
(n, m-
1, heights
) - getHeight
(n-
1, m, heights
);
}
currDecoder =
(followerDDiff ==
0) ? decoderPF0: decoderPFN0
;
currDecoder.
setDDiff(followerDDiff
);
}
ct = currDecoder.
type;
int numZeroBits = countZeroBits
(bytes
);
int delta
;
int critNumZero = maxZeroBits -
(currDecoder == decoderStd
? 0 :
1 + plateauBinBits
[lastTablePos
]);
if (numZeroBits
<= critNumZero
&& section.
shrinkFactor > 1 && numZeroBits
> maxZeroBits -
8)
System.
out.
println("len??? " + currDecoder.
type +
" " + numZeroBits +
" " + maxZeroBits +
" " + critNumZero
);
if (numZeroBits
> critNumZero
) {
if (section.
shrinkFactor > 1)
System.
out.
println("big bin " + currDecoder.
type +
" " + numZeroBits +
" " + maxZeroBits +
" " + critNumZero
);
delta = decodeBigBin
(bytes, getBigBinBits
(encDeltaHeight
));
} else if (currDecoder.
getEncType() == EncType.
HYBRID) {
delta = decodeHybrid
(bytes, numZeroBits, currDecoder.
getHUnit());
} else {
delta = decodeLen
(numZeroBits
);
}
int delta1 = currDecoder.
decode1(delta
);
int delta2 = currDecoder.
processVal(delta1
);
int v
;
int hUp = getHeight
(n, m-
1, heights
);
// calculate predicted value and add delta from file
if (ct == CalcType.
CALC_STD) {
int predict
;
int hLeft = getHeight
(n-
1, m, heights
);
int hUpLeft = getHeight
(n-
1, m-
1, heights
);
// calculate incline of upper values
int hdiffUp = hUp- hUpLeft
;
// we expect same incline as in the upper row
predict = hLeft + hdiffUp
;
if (predict
> encDeltaHeight
) {
predict = encDeltaHeight
;
} else if (predict
< 0)
predict =
0;
if (hUp- hLeft
> 0)
v = predict - delta2
;
else
v = predict + delta2
;
} else {
// plateau follower: predicted value is simply upper value
v = hUp + delta2
;
}
// unwrap
if (v
< 0) {
v += encDeltaHeight +
1;
} else if (v
> encDeltaHeight
) {
v -= encDeltaHeight +
1;
}
if (v
< 0 || v
> encDeltaHeight
)
throw new RuntimeException("invalid v:" + v
);
if (!silent
)
reportBits
((n+m
*width
) +
"(" + n +
"," + m +
") raw= " + delta1 +
" -> d= " + v
);
addHeight
(v, heights
);
check
(heights, section.
zoomLevel, latOff, lonOff
);
if (bytePos
> this.
dataLength)
break;
// calc pos in matrix
n = heights.
size();
m = n / width
;
n = n - width
* m
;
if (m == height
)
break;
// determine how to interpret next bits
int hl = getHeight
(n-
1,m,heights
);
int dDiff =
0;
if (m
> 0) {
if (n
> 0) {
int hu = getHeight
(n,m-
1,heights
);
dDiff = hu- hl
;
}
} else {
dDiff = -hl
;
}
readStd =
(dDiff
!=
0);
}
if (bytePos +
1 < dataLength
) {
System.
out.
println(("Decoded sequence: " + heights
));
throw new RuntimeException("did not read all avalaible bytes, remaining: " +
(dataLength - bytePos+
1) );
} else if (heights.
size() != numVals
) {
throw new RuntimeException("did not find all expected values");
} else
System.
out.
println("ok, decoded expected " + numVals +
" values for shrink="+section.
shrinkFactor +
", maxDeltaHeight=" + maxDeltaHeight +
", width=" + width +
", height=" + height
);
return;
}
private void reportBits
(String string
) {
bitString.
append(' ').
append(string
);
System.
out.
println(bitString
);
bitString.
setLength(0);
}
/**
* Some plausibility checks
* @param heights
* @param zoomLevel
* @param lonOff
* @param latOff
*/
private void check
(List<Integer> heights,
int zoomLevel,
int latOff,
int lonOff
) {
if (true && heights.
size() > 0 ) {
int pos = heights.
size() -
1;
int x = pos
% width
;
int y = pos / width
;
int v = getHeight
(x,y,heights
);
boolean undef = isUndef
(v
);
// plausibility checks
int lat = section.
top - latOff - y
* section.
pointsDistanceLat;
int lon = section.
left + lonOff + x
* section.
pointsDistanceLon;
if (!silent
) {
int realHeight = v
* section.
shrinkFactor + baseHeight
;
boolean limited =
false;
if (v == encDeltaHeight
) {
if (realHeight
> baseHeight + maxDeltaHeight
) {
realHeight = baseHeight + maxDeltaHeight
;
limited =
true;
}
}
StringBuilder bs =
new StringBuilder(64);
bs.
append("h(").
append(x
).
append(',').
append(y
).
append("):");
if (undef
) {
bs.
append("undef ");
} else {
realHeight =
Math.
round(useMeters
? realHeight : realHeight
* 0.3048f
);
bs.
append(realHeight
).
append(" m ");
if (limited
) {
bs.
append("(limitied) ");
}
}
bs.
append(pos2String
(lat, lon
));
System.
out.
println(bs
);
}
// if (heights.size() == 1)
// return;
// if (!undef & !silent && zoomLevel == 0 && section.pointsDistanceLat < 10000) {
// int maxDiff = 3;
// int left = getHeight(x-1,y,heights);
// if (!isUndef(left)) {
// int hDiff = (v - left) ;
// if (Math.abs(hDiff) > maxDiff && x > 0) {
// System.out.println("large hDiff : (" + x + "," + y+"):" + hDiff + " " + pos2String(lat, lon) + " " + this);
// return;
// }
// }
// int up = getHeight(x,y-1,heights);
// if (!isUndef(up)) {
// int vDiff = (v - up);
// if (Math.abs(vDiff) > maxDiff && y > 0) {
// System.out.println("large vDiff: (" + x + "," + y+"):" + vDiff + " " + pos2String(lat, lon) + " " + this);
// return;
// }
// }
// }
}
}
/**
* Evaluate if the height means "undefined"
* @param h the height
* @return true if height means undefined, else false
*/
private boolean isUndef
(int h
) {
// the number of 1 bits in the field encodingType seems to give the number of heights that should be considered invalid?
if (encodingType ==
0)
return false;
if (h == encDeltaHeight
&& (encodingType ==
1 || encodingType ==
2 || encodingType ==
4))
return true;
if (h
> encDeltaHeight -
1 && (encodingType ==
3 || encodingType ==
5 || encodingType ==
6)) {
// TODO: not sure about this case
return true;
}
return false;
}
@
Override
public String toString
() {
return "DemTile [tileNumberLat=" + tileNumberLat +
", tileNumberLon=" + tileNumberLon +
", height=" + height
+
", width=" + width +
", offset=" + offset +
", dataLength=" + dataLength +
", baseHeight="
+ baseHeight +
", maxDeltaHeight=" + maxDeltaHeight +
"]";
}
private void addHeight
(int v,
List<Integer> heights
) {
heights.
add(v
);
}
private int getHeight
(int x,
int y,
List<Integer> heights
) {
if (y
< 0)
return 0; // virtual 1st line
if (x
< 0) {
// virtual column
return y ==
0 ? 0 : heights.
get((y -
1) * width
);
}
int pos = x + y
* width
;
if (heights.
size() > pos
)
return heights.
get(pos
);
return 0; // hmmm ?
}
private int decodePlateauLen
(int col,
byte[] bytes
) {
int currCol = col
;
int plateauLen =
0;
int unit =
1;
int tablePos =
0;
if (lastTablePos
> 0)
tablePos = lastTablePos
;
int rowLen = currCol
;
while (true) {
boolean b = consumeBit
(bytes
);
if (!b
)
break;
unit = plateauUnit
[tablePos
];
rowLen += unit
;
tablePos++
;
if (rowLen
>= width
)
break;
}
if (tablePos
> 0 && rowLen
!= width
)
tablePos--
;
if (rowLen
< width
) {
int binBitsToRead = plateauBinBits
[tablePos
];
if (binBitsToRead
> 0) {
int binVal = decodeBin
(bytes, binBitsToRead
);
rowLen += binVal
;
}
}
lastTablePos = tablePos
;
if (rowLen
< width
) {
plateauLen += rowLen - currCol
;
if (plateauLen
> 0 && section.
shrinkFactor > 1) {
long dd =
4;
}
return plateauLen
;
}
if (rowLen
> width
) {
long dd =
4;
}
plateauLen += width - currCol
;
return plateauLen
;
}
/**
* Decode a value given as number of zero bits.
* @param numZeroBits
* @return the decoded value
*/
private int decodeLen
(int numZeroBits
) {
// different meaning of odd / even number of bits
int v0 =
((numZeroBits
& 1) !=
0) ? (numZeroBits +
1) >> 1 : -
(numZeroBits
>> 1);
return v0
;
}
private int countZeroBits
(byte[] bytes
) {
int s =
0;
// read 0-bits
while (!consumeBit
(bytes
)) {
s++
;
}
return s
;
}
private int decodeHybrid
(byte[] bytes,
int numZeroBits,
int hunit
) {
assert hunit
> 0;
int numBits =
Integer.
numberOfTrailingZeros(hunit
);
int binVal = decodeBin
(bytes, numBits
);
// read sign bit
boolean isPositive = consumeBit
(bytes
);
int num = numZeroBits
* hunit + binVal
;
num = isPositive
? num +
1 : -num
;
return num
;
}
private int decodeBin
(byte[] bytes,
int numBits
) {
if (numBits ==
0)
return 0;
int binVal =
0;
int meaning =
1 << numBits
;
for (int i =
0; i
< numBits
; i++
) {
meaning
>>=
1;
if (consumeBit
(bytes
))
binVal |= meaning
;
}
return binVal
;
}
private int decodeBigBin
(byte[] bytes,
int numBits
) {
int v = decodeBin
(bytes, numBits-
1);
v++
;
boolean sign = consumeBit
(bytes
);
return sign
? -v : v
;
}
private boolean consumeBit
(byte[] bytes
) {
if (bitPos
> 7) {
bitPos =
0;
bytePos++
;
}
if (bytePos
>= bytes.
length)
throw new RuntimeException("tried to read nonexisting byte");
int b = bytes
[bytePos
] & 0xff
;
b
>>=
7 - bitPos++
;
boolean rc =
(b
& 1) ==
1;
if (!silent
) {
bitString.
append(rc
? '1' :
'.');
}
return rc
;
}
}
/**
* Calculate the longest sequence of zero bits that is considered as a valid number.
* The number depends on the field maxDeltaHeight in the header structure and
* the shrink factor.
* @param shrink
* @param maxRealDelta
* @return the number of bits
*/
private static int getMaxLengthZeroBits
(int shrink,
int maxRealDelta
) {
if (shrink ==
3 && maxRealDelta
<=
127) {
if (maxRealDelta
< 69)
return 21;
if (maxRealDelta
< 96)
return 23;
return 22;
}
int m = getMaxLengthZeroBitsUnchecked
(shrink, maxRealDelta
);
if (m
< 9)
return 9;
return m
;
}
/**
* Calculate the longest sequence of zero bits that is considered as a valid number.
* The number depends on the field maxDeltaHeight in the header structure and
* the shrink factor.
* @param shrink
* @param maxRealDelta
* @return the number of bits, sometimes wrong, see getMaxLengthZeroBits
*/
private static int getMaxLengthZeroBitsUnchecked
(int shrink,
int maxRealDelta
) {
if (shrink ==
1) {
int ld = int_ld
(maxRealDelta
);
return (maxRealDelta
< 256) ? ld +
15 : ld
* 3 +
1;
}
if (maxRealDelta
> 0) {
TreeMap<Integer,
Integer> tmp =
new TreeMap<>();
int s = int_ld
((shrink -
1) /
2);
for (int i =
0; i
< 16; i++
) {
int v =
3 * (i +
1) + s
;
int k =
(int) Math.
pow(2, i
);
tmp.
put(k, v
);
}
for (int i =
1; i
< 16; i++
) {
int k = shrink
* ((int) Math.
pow(2, i
) -
1) +
1;
if (k
>=
65536)
break;
if (tmp.
containsKey(k
)) {
tmp.
put(k, tmp.
get(k
) -
1);
} else
tmp.
put(k, -
1);
}
Entry
<Integer,
Integer> e
;
Integer hit = tmp.
get(maxRealDelta
);
if (hit
!=
null) {
if (hit
> 0)
return hit
;
e = tmp.
lowerEntry(maxRealDelta
);
return e.
getValue() -
1;
}
e = tmp.
lowerEntry(maxRealDelta
);
if (e.
getValue() > 0)
return e.
getValue();
e = tmp.
lowerEntry(e.
getKey());
return e.
getValue() -
1;
}
throw new RuntimeException("did not find max MaxLengthZeroBits for " + shrink +
" " + maxRealDelta
);
}
private static int int_ld
(int v
) {
assert v
> 0;
// return Integer.numberOfTrailingZeros(Integer.highestOneBit(v));
return (int) Math.
floor(Math.
log(v
) /
Math.
log(2));
}
static int getBigBinBits
(int maxHeight
) {
int bits
;
if (maxHeight
<16384) {
int n =
Integer.
highestOneBit(maxHeight
);
bits =
Integer.
numberOfTrailingZeros(n
) +
1;
} else
bits =
15;
return bits
;
// if (maxHeight < 2)
// return 1;
// else if (maxHeight < 4)
// return 2;
// else if (maxHeight < 8)
// return 3;
// else if (maxHeight < 16)
// return 4;
// else if (maxHeight < 32)
// return 5;
// else if (maxHeight < 64)
// return 6;
// else if (maxHeight < 128)
// return 7;
// else if (maxHeight < 256)
// return 8;
// else if (maxHeight < 512)
// return 9;
// else if (maxHeight < 1024)
// return 10;
// else if (maxHeight < 2048)
// return 11;
// else if (maxHeight < 4096)
// return 12;
// else if (maxHeight < 8192)
// return 13;
// else if (maxHeight < 16384)
// return 14;
// else
// return 15;
}
private class DemSection
{
private int zoomLevel
;
private int pointsPerLat
;
private int pointsPerLon
;
private int nonStdHeight
;
private int nonStdWidth
;
private int shrinkFactor
;
private int tilesLat
;
private int tilesLon
;
private byte offsetSize
;
private byte baseSize
;
private byte differenceSize
;
private boolean hasExtra
;
private short tileDescSize
;
private int dataOffset
;
private int dataOffset2
;
private int pointsDistanceLat
;
private int pointsDistanceLon
;
private int top
;
private int left
;
private int maxRealHeight
;
private DemTile
[] tiles
;
protected void print
() {
readHeader
();
}
private void readHeader
() {
Displayer d =
new Displayer
(reader
);
d.
setTitle("DEM Sector Header");
d.
byteValue("unknown byte %d"); //0x00
zoomLevel = d.
byteValue("zoom level %d"); //0x01
pointsPerLat = d.
intValue("Points per tile (latitude) %d"); //0x02 - 0x05
pointsPerLon = d.
intValue("Points per tile (longitude) %d");//0x06 - 0x09
nonStdHeight =
1 + d.
intValue("height -1 of non-std tile at bottom: %d "); //0x0a - 0x0d
nonStdWidth =
1 + d.
intValue("width - 1 of non-std tile at right: %d"); //0x0e - 0x11
int shrinkVal = d.
shortValue("shrink value"); //0x12 - 0x13
shrinkFactor =
2 * shrinkVal +
1;
d.
item().
addText("2*value + 1 = shrink factor => " + shrinkFactor
);
tilesLon =
1 + d.
intValue("Tiles -1 per longitude: %d"); //0x14 - 0x17
tilesLat =
1 + d.
intValue("Tiles -1 per latitude: %d"); //0x18 - 0x1B
short recordDesc = d.
shortValue("Flags for tile description"); //0x1C - 0x1D
offsetSize =
(byte)((recordDesc
& 3) +
1);
baseSize =
(byte)(((recordDesc
& 4) >> 2) +
1);
differenceSize =
(byte) (((recordDesc
& 8) >> 3) +
1);
hasExtra =
(recordDesc
& 0x10
) !=
0;
DisplayItem item = d.
item();
item.
addText("Data offset has %s bytes", offsetSize
);
item.
addText("Base height has %s bytes", baseSize
);
item.
addText("Height difference has %s bytes", differenceSize
);
item.
addText("Extra flag: %d",
(recordDesc
& 0x10
) >> 4);
tileDescSize = d.
shortValue("Size of tile description: %d");//0x1E - 0x1F
dataOffset = d.
intValue("offset of first tile description: %#08x");// 0x20 - 0x23
dataOffset2 = d.
intValue("offset of first tiles data: %#08x"); //0x24 - 0x27
left = d.
intValue("western bound DEM units: %d "); //0x28-0x2B
top = d.
intValue("northern bound DEM units: %d"); //0x2C - 0x2F
d.
item().
addText("Upper left in degrees lat,lon: %s,%s",
DemDisplay.
toDegrees(top
),
DemDisplay.
toDegrees(left
));
pointsDistanceLat = d.
intValue("DEM units between points (latitude): %d"); //0x30 - 0x33
pointsDistanceLon = d.
intValue("DEM units between points (longitude): %d"); //0x34 - 0x37
int right = left +
((tilesLon-
1) * pointsPerLon + nonStdWidth
) * pointsDistanceLon
;
int bottom = top -
((tilesLat-
1) * pointsPerLat + nonStdHeight
) * pointsDistanceLat
;
d.
item().
addText("Calc. eastern bound DEM units: " + right
);
d.
item().
addText("Calc. southern bound DEM units: " + bottom
);
d.
item().
addText("Calc. lower right corner in degrees lat,lon: %s,%s",
DemDisplay.
toDegrees(bottom
),
DemDisplay.
toDegrees(right
));
d.
shortValue("minimum (base) height: %d"); //0x38 - 0x39
maxRealHeight = d.
shortValue("maximum height: %d"); //0x3A - 0x3B
d.
print(outStream
);
tiles =
new DemTile
[tilesLat
* tilesLon
];
DemTile lastDataTile =
null;
// now read the tile description
for (int i =
0; i
< tiles.
length; ++i
) {
tiles
[i
] =
new DemTile
(this, i / tilesLon, i
% tilesLon
);
tiles
[i
].
readHeader(dataOffset2
);
if (i
!=
0 && lastDataTile
!=
null) {
lastDataTile.
setDataLength(tiles
[i
].
offset);
}
if (tiles
[i
].
maxDeltaHeight !=
0)
lastDataTile = tiles
[i
];
// data length of last tile has to be set after reading next section
}
}
private void setLastDataLength
(int nextSectionOffset
) {
for (int i = tiles.
length -
1; i
>=
0; i--
) {
if (tiles
[i
].
maxDeltaHeight !=
0) {
tiles
[i
].
setDataLength(nextSectionOffset
);
break;
}
}
}
public void decodeBitstreams
() {
Displayer d =
new Displayer
(reader
);
d.
setTitle(String.
format("Bitstream data of zoom level %s", zoomLevel
));
d.
print(outStream
);
int lonOff =
0;
int latOff =
0;
for (int i =
0; i
< tiles.
length; ++i
) {
// tiles[i].printRawData(reader);
lonOff = tiles
[i
].
tileNumberLon * pointsPerLon
* pointsDistanceLon
;
latOff = tiles
[i
].
tileNumberLat * pointsPerLat
* pointsDistanceLat
;
tiles
[i
].
decodeBitStreamData(lonOff, latOff
);
}
}
}
private DemSection
[] demSections
;
private int sectionOffset
;
private short sectionSize
;
private boolean useMeters
;
private boolean silent =
true;
@
Override
protected void print
() {
reader.
position(reader.
getGMPOffset());
readCommonHeader
();
readFileHeader
();
}
private void readFileHeader
() {
Displayer d =
new Displayer
(reader
);
d.
setTitle("DEM Header");
useMeters =
1 !=
(d.
intValue("Flags ?") & 1);
DisplayItem item = d.
item();
if (useMeters
) {
item.
addText("elevation is given in meters");
} else {
item.
addText("elevation is given in feet");
}
demSections =
new DemSection
[d.
shortValue("Number of zoom levels/sectors")];
d.
intValue("unknown integer");
sectionSize = d.
shortValue("size of sector headers");
sectionOffset = d.
intValue("offset of first sector header");
if (getHeaderLen
() > 0x25
) {
d.
intValue("unknown integer at 0x25");
}
d.
print(outStream
);
reader.
position(sectionOffset
);
for (int i =
0; i
< demSections.
length; i++
) {
reader.
position(sectionOffset + sectionSize
* i
);
demSections
[i
] =
new DemSection
();
demSections
[i
].
print();
if (i
!=
0) {
// data of last tile of previous section ends before first tile of this section
demSections
[i-
1].
setLastDataLength(demSections
[i
].
dataOffset);
}
if (i +
1 == demSections.
length) {
// data of last tile of last section ends before begin of section description
demSections
[i
].
setLastDataLength(sectionOffset
);
}
}
// now iterate a second time over section and tiles and print out raw data
for (int i =
0; i
< demSections.
length; ++i
) {
demSections
[i
].
decodeBitstreams();
}
}
private static String toDegrees
(int val
) {
return String.
format(Locale.
ENGLISH,
"%.7f", val
* (45.0 /
(1<<29)));
}
private static String pos2String
(int lat,
int lon
) {
return "at " + toDegrees
(lat
) +
" " + toDegrees
(lon
);
}
/**
* Helper code to print table of max zero bit lengths.
*/
private static void printMaxZeroTable
() {
for (int s =
1; s
<=
15; s+=
2) {
int lastM = -
1;
for (int i = s
; i
<=
Short.
MAX_VALUE; i += s
) {
int m = getMaxLengthZeroBits
(s, i
);
if (m
!= lastM
)
System.
out.
println(String.
format("s=%d h=%d h2=%d m=%d", s ,i,i/s,m
));
lastM = m
;
}
}
}
public static void main
(String[] args
) {
if (args.
length < 1) {
System.
err.
println("Usage: DemdDsplay [--verbose] <filename> ");
System.
exit(1);
}
DemDisplay td =
new DemDisplay
();
if ("--verbose".
equals(args
[0]))
td.
silent =
false;
td.
display(args
[args.
length-
1],
"DEM");
System.
out.
println("number of errors " + td.
countErr +
" and " + td.
countGood +
" were OK, " + td.
countSorry +
" are not yet understood");
}
}