diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/ConvertSrtmTile.java b/brouter-map-creator/src/main/java/btools/mapcreator/ConvertSrtmTile.java index 67346d2..f628152 100644 --- a/brouter-map-creator/src/main/java/btools/mapcreator/ConvertSrtmTile.java +++ b/brouter-map-creator/src/main/java/btools/mapcreator/ConvertSrtmTile.java @@ -8,6 +8,7 @@ public class ConvertSrtmTile public static int NROWS; public static int NCOLS; + public static final short SKIPDATA = -32766; // >50 degree skipped pixel public static final short NODATA2 = -32767; // bil-formats nodata public static final short NODATA = Short.MIN_VALUE; @@ -17,8 +18,6 @@ public class ConvertSrtmTile private static void readBilZip( String filename, int rowOffset, int colOffset, boolean halfCols ) throws Exception { - int fileRows = 3601; - int fileCols = 3601; ZipInputStream zis = new ZipInputStream( new BufferedInputStream( new FileInputStream( filename ) ) ); try { @@ -27,7 +26,7 @@ public class ConvertSrtmTile ZipEntry ze = zis.getNextEntry(); if ( ze.getName().endsWith( ".bil" ) ) { - readBilFromStream( zis, rowOffset, colOffset, fileRows, fileCols, halfCols ); + readBilFromStream( zis, rowOffset, colOffset, halfCols ); return; } } @@ -38,80 +37,69 @@ public class ConvertSrtmTile } } - private static void readBilFromStream( InputStream is, int rowOffset, int colOffset, int fileRows, int fileCols, boolean halfCols ) + private static void readBilFromStream( InputStream is, int rowOffset, int colOffset, boolean halfCols ) throws Exception { DataInputStream dis = new DataInputStream( new BufferedInputStream( is ) ); - for ( int ir = 0; ir < fileRows; ir++ ) + for ( int ir = 0; ir < 3601; ir++ ) { int row = rowOffset + ir; - short lastVal = 0; - boolean fillGap = false; - - for ( int ic = 0; ic < fileCols; ic++ ) + for ( int ic = 0; ic < 3601; ic++ ) { int col = colOffset + ic; - short val; + if ( ( ic % 2 ) == 1 && halfCols ) { - fillGap = true; + if ( getPixel( row, col ) == NODATA ) + { + setPixel( row, col, SKIPDATA ); + } + continue; } - else + + int i0 = dis.read(); + int i1 = dis.read(); + + if ( i0 == -1 || i1 == -1 ) + throw new RuntimeException( "unexcepted end of file reading bil entry!" ); + + short val = (short) ( ( i1 << 8 ) | i0 ); + + if ( val == NODATA2 ) { - int i0 = dis.read(); - int i1 = dis.read(); - - if ( i0 == -1 || i1 == -1 ) - throw new RuntimeException( "unexcepted end of file reading bil entry!" ); - - val = (short) ( ( i1 << 8 ) | i0 ); - - if ( val == NODATA2 ) - { - val = NODATA; - } - - if ( fillGap ) - { - setPixel( row, col - 1, val, lastVal ); - fillGap = false; - } - -if ( row == 18010 ) System.out.print( val + " " ); - - setPixel( row, col, val, val ); - lastVal = val; + val = NODATA; } + + setPixel( row, col, val ); } } } - private static void setPixel( int row, int col, short val1, short val2 ) + + private static void setPixel( int row, int col, short val ) { if ( row >= 0 && row < NROWS && col >= 0 && col < NCOLS ) { - if ( val1 != NODATA && val2 != NODATA ) - { - int val = val1 + val2; - if ( val < -32768 || val > 32767 ) - throw new IllegalArgumentException( "val1=" + val1 + " val2=" + val2 ); - - imagePixels[row * NCOLS + col] = (short) ( val ); - } + imagePixels[row * NCOLS + col] = val; } } - public static void main( String[] args ) throws Exception + private static short getPixel( int row, int col ) { - doConvert( args[0], Integer.parseInt( args[1] ), Integer.parseInt( args[2] ), args[3], null ); + if ( row >= 0 && row < NROWS && col >= 0 && col < NCOLS ) + { + return imagePixels[row * NCOLS + col]; + } + return NODATA; } - public static void doConvert( String inputDir, int lonDegreeStart, int latDegreeStart, String outputFile, SrtmRaster raster90 ) throws Exception + + public static void doConvert( String inputDir, String v1Dir, int lonDegreeStart, int latDegreeStart, String outputFile, SrtmRaster raster90 ) throws Exception { int extraBorder = 10; int datacells = 0; - int matchingdatacells = 0; + int mismatches = 0; NROWS = 5 * 3600 + 1 + 2 * extraBorder; NCOLS = 5 * 3600 + 1 + 2 * extraBorder; @@ -152,59 +140,63 @@ if ( row == 18010 ) System.out.print( val + " " ); } } - if ( raster90 != null ) - { - for ( int row90 = 0; row90 < 6001; row90++ ) - { - int crow = 3 * row90 + extraBorder; // center row of 3x3 - for ( int col90 = 0; col90 < 6001; col90++ ) - { - int ccol = 3 * col90 + extraBorder; // center col of 3x3 + boolean halfCol5 = latDegreeStart >= 50 || latDegreeStart < -50; + for ( int row90 = 0; row90 < 6001; row90++ ) + { + int crow = 3 * row90 + extraBorder; // center row of 3x3 + for ( int col90 = 0; col90 < 6001; col90++ ) + { + int ccol = 3 * col90 + extraBorder; // center col of 3x3 + + // evaluate 3x3 area + if ( raster90 != null && (!halfCol5 || (col90 % 2) == 0 ) ) + { short v90 = raster90.eval_array[row90 * 6001 + col90]; - // evaluate 3x3 area int sum = 0; int nodatas = 0; + int datas = 0; + int colstep = halfCol5 ? 2 : 1; for ( int row = crow - 1; row <= crow + 1; row++ ) { - for ( int col = ccol - 1; col <= ccol + 1; col++ ) + for ( int col = ccol - colstep; col <= ccol + colstep; col += colstep ) { short v30 = imagePixels[row * NCOLS + col]; - sum += v30; if ( v30 == NODATA ) { nodatas++; } + else if ( v30 != SKIPDATA ) + { + sum += v30; + datas++; + } } } - boolean doReplace = nodatas > 0 || v90 == NODATA; - short replaceValue = NODATA; + boolean doReplace = nodatas > 0 || v90 == NODATA || datas < 7; if ( !doReplace ) { datacells++; - int diff = sum - 9 * v90; - replaceValue= (short)(diff / 9); - doReplace = true; - if ( diff < -9 || diff > 9 ) + int diff = sum - datas * v90; + if ( diff < -4 || diff > 4 ) { -// System.out.println( "replacing due to median missmatch: sum=" + sum + " v90=" + v90 ); doReplace = true; + mismatches++; } - if ( diff > -50 && diff < 50 ) + + if ( diff > -50 && diff < 50 && ( row90 % 1200 ) != 0 && ( col90 % 1200 ) != 0 ) { - matchingdatacells++; - diffs[diff+50]++; + diffs[diff + 50]++; } } if ( doReplace ) { for ( int row = crow - 1; row <= crow + 1; row++ ) { - for ( int col = ccol - 1; col <= ccol + 1; col++ ) + for ( int col = ccol - colstep; col <= ccol + colstep; col += colstep ) { -// imagePixels[row * NCOLS + col] = v90; - imagePixels[row * NCOLS + col] = replaceValue; + imagePixels[row * NCOLS + col] = v90; } } } @@ -215,6 +207,7 @@ if ( row == 18010 ) System.out.print( val + " " ); SrtmRaster raster = new SrtmRaster(); raster.nrows = NROWS; raster.ncols = NCOLS; + raster.halfcol = halfCol5; raster.noDataValue = NODATA; raster.cellsize = 1 / 3600.; raster.xllcorner = lonDegreeStart - ( 0.5 + extraBorder ) * raster.cellsize; @@ -235,16 +228,32 @@ if ( row == 18010 ) System.out.print( val + " " ); if ( pix2.length != imagePixels.length ) throw new RuntimeException( "length mismatch!" ); - for ( int i = 0; i < pix2.length; i++ ) + // compare decoding result + for ( int row = 0; row < NROWS; row++ ) { - if ( pix2[i] != imagePixels[i] ) + int colstep = halfCol5 ? 2 : 1; + for ( int col = 0; col < NCOLS; col += colstep ) { - throw new RuntimeException( "content mismatch!" ); + int idx = row * NCOLS + col; + if ( imagePixels[idx] == SKIPDATA ) + { + continue; + } + short p2 = pix2[idx]; + if ( p2 > SKIPDATA ) + { + p2 /= 2; + } + if ( p2 != imagePixels[idx] ) + { + throw new RuntimeException( "content mismatch!" ); + } } } - for(int i=0; i<100;i++) System.out.println( "diff[" + (i-50) + "] = " + diffs[i] ); - System.out.println( "datacells=" + datacells + " mismatch%=" + 100.*(datacells-matchingdatacells)/datacells ); + for(int i=1; i<100;i++) System.out.println( "diff[" + (i-50) + "] = " + diffs[i] ); + System.out.println( "datacells=" + datacells + " mismatch%=" + (100.*mismatches)/datacells ); +btools.util.MixCoderDataOutputStream.stats(); // test( raster ); // raster.calcWeights( 50. ); // test( raster ); diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/ConvertUrlList.java b/brouter-map-creator/src/main/java/btools/mapcreator/ConvertUrlList.java index c2ba42b..b41e31a 100644 --- a/brouter-map-creator/src/main/java/btools/mapcreator/ConvertUrlList.java +++ b/brouter-map-creator/src/main/java/btools/mapcreator/ConvertUrlList.java @@ -42,11 +42,6 @@ public class ConvertUrlList int ilon_base = ( srtmLonIdx - 1 ) * 5 - 180; int ilat_base = 150 - srtmLatIdx * 5 - 90; - if ( ilon_base < 5 || ilon_base > 10 || ilat_base < 45 || ilat_base > 50 ) - { - continue; - } - SrtmRaster raster90 = null; File file90 = new File( new File( args[1] ), filename90 ); @@ -56,7 +51,7 @@ public class ConvertUrlList raster90 = new SrtmData( file90 ).getRaster(); } - ConvertSrtmTile.doConvert( args[2], ilon_base, ilat_base, filename30, raster90 ); + ConvertSrtmTile.doConvert( args[2], args[3], ilon_base, ilat_base, filename30, raster90 ); } br.close(); } diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/RasterCoder.java b/brouter-map-creator/src/main/java/btools/mapcreator/RasterCoder.java index 4471612..b7d3495 100644 --- a/brouter-map-creator/src/main/java/btools/mapcreator/RasterCoder.java +++ b/brouter-map-creator/src/main/java/btools/mapcreator/RasterCoder.java @@ -17,6 +17,7 @@ public class RasterCoder dos.writeInt(raster.ncols); dos.writeInt(raster.nrows); + dos.writeBoolean(raster.halfcol); dos.writeDouble(raster.xllcorner); dos.writeDouble(raster.yllcorner); dos.writeDouble(raster.cellsize); @@ -37,6 +38,7 @@ public class RasterCoder SrtmRaster raster = new SrtmRaster(); raster.ncols = dis.readInt(); raster.nrows = dis.readInt(); + raster.halfcol = dis.readBoolean(); raster.xllcorner = dis.readDouble(); raster.yllcorner = dis.readDouble(); raster.cellsize = dis.readDouble(); @@ -59,12 +61,26 @@ public class RasterCoder int nrows = raster.nrows; int ncols = raster.ncols; short[] pixels = raster.eval_array; + int colstep = raster.halfcol ? 2 : 1; for (int row = 0; row < nrows; row++) { - for (int col = 0; col < ncols; col++) + short lastval = Short.MIN_VALUE; // nodata + for (int col = 0; col < ncols; col += colstep ) { - mco.writeMixed(pixels[row * ncols + col]); + short val = pixels[row * ncols + col]; + if ( val == -32766 ) + { + val = lastval; // replace remaining (border) skips + } + else + { + lastval = val; + } + + // remap nodata + int code = val == Short.MIN_VALUE ? -1 : ( val < 0 ? val-1 : val ); + mco.writeMixed( code ); } } mco.flush(); @@ -76,12 +92,35 @@ public class RasterCoder int nrows = raster.nrows; int ncols = raster.ncols; short[] pixels = raster.eval_array; + int colstep = raster.halfcol ? 2 : 1; for (int row = 0; row < nrows; row++) { - for (int col = 0; col < ncols; col++) + for (int col = 0; col < ncols; col += colstep ) { - pixels[row * ncols + col] = (short) mci.readMixed(); + int code = mci.readMixed(); + + // remap nodata + int v30 = code == -1 ? Short.MIN_VALUE : ( code < 0 ? code + 1 : code ); + if ( v30 > -32766 ) + { + v30 *= 2; + } + pixels[row * ncols + col] = (short) ( v30 ); + } + if ( raster.halfcol ) + { + for (int col = 1; col < ncols-1; col += colstep ) + { + int l = (int)pixels[row * ncols + col - 1]; + int r = (int)pixels[row * ncols + col + 1]; + short v30 = Short.MIN_VALUE; // nodata + if ( l > -32766 && r > -32766 ) + { + v30 = (short)((l+r)/2); + } + pixels[row * ncols + col] = v30; + } } } } diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/SrtmData.java b/brouter-map-creator/src/main/java/btools/mapcreator/SrtmData.java index ea42f96..026db04 100644 --- a/brouter-map-creator/src/main/java/btools/mapcreator/SrtmData.java +++ b/brouter-map-creator/src/main/java/btools/mapcreator/SrtmData.java @@ -101,7 +101,7 @@ public class SrtmData { if ( negative ) n = -n; - short val = n < -250 ? Short.MIN_VALUE : (short) (n*2); + short val = n < -250 ? Short.MIN_VALUE : (short) (n); raster.eval_array[row * raster.ncols + col] = val; if ( ++col == raster.ncols ) diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/SrtmRaster.java b/brouter-map-creator/src/main/java/btools/mapcreator/SrtmRaster.java index 67f95e4..2e10e64 100644 --- a/brouter-map-creator/src/main/java/btools/mapcreator/SrtmRaster.java +++ b/brouter-map-creator/src/main/java/btools/mapcreator/SrtmRaster.java @@ -11,6 +11,7 @@ public class SrtmRaster { public int ncols; public int nrows; + public boolean halfcol; public double xllcorner; public double yllcorner; public double cellsize; @@ -50,7 +51,7 @@ public class SrtmRaster + (1.-wrow)*( wcol)*get(row ,col+1) + ( wrow)*( wcol)*get(row+1,col+1); // System.out.println( "eval=" + eval ); - return missingData ? Short.MIN_VALUE : (short)(eval*2); + return missingData ? Short.MIN_VALUE : (short)(eval*4); } private short get( int r, int c ) @@ -72,10 +73,7 @@ public class SrtmRaster double drow = (lat - yllcorner)/cellsize; int row = (int)drow; int col = (int)dcol; - if ( col < 9 || col >= ncols-9 || row < 9 || row >= nrows-9 ) - { - throw new IllegalArgumentException( "out of boundary range: col=" + col + " row=" + row ); - } + double dgx = (dcol-col)*gridSteps; double dgy = (drow-row)*gridSteps; @@ -110,6 +108,8 @@ public class SrtmRaster double m = (1.-wlat) * m0 + wlat * m1; return (short)(m * 2); } + + private ReducedMedianFilter rmf = new ReducedMedianFilter( 256 ); private double getElevation( Weights w, int row, int col ) { @@ -122,9 +122,9 @@ public class SrtmRaster int mx = nx / 2; // mean pixels int my = ny / 2; - System.out.println( "nx="+ nx + " ny=" + ny ); - - ReducedMedianFilter rmf = new ReducedMedianFilter( nx * ny ); + // System.out.println( "nx="+ nx + " ny=" + ny ); + + rmf.reset(); for( int ix = 0; ix < nx; ix ++ ) { @@ -134,7 +134,7 @@ public class SrtmRaster rmf.addSample( w.getWeight( ix, iy ), val ); } } - return missingData ? rmf.calcEdgeReducedMedian( filterCenterFraction ) : 0.; + return missingData ? 0. : rmf.calcEdgeReducedMedian( filterCenterFraction ); } @@ -190,8 +190,8 @@ public class SrtmRaster private static int gridSteps = 10; private static Weights[][][] allShiftWeights = new Weights[17][][]; - private static double filterCenterFraction = 0.5; - private static double filterDiscRadius = 2.9; // in pixels + private static double filterCenterFraction = 0.2; + private static double filterDiscRadius = 4.999; // in pixels static { diff --git a/brouter-util/src/main/java/btools/util/MixCoderDataInputStream.java b/brouter-util/src/main/java/btools/util/MixCoderDataInputStream.java index a7d911a..8e86904 100644 --- a/brouter-util/src/main/java/btools/util/MixCoderDataInputStream.java +++ b/brouter-util/src/main/java/btools/util/MixCoderDataInputStream.java @@ -12,46 +12,74 @@ import java.io.InputStream; public final class MixCoderDataInputStream extends DataInputStream { - private long lastValue; - private long repCount; + private int lastValue; + private int repCount; + private int diffshift; + + private int bm = 0x100; + private int b; public MixCoderDataInputStream( InputStream is ) { super( is ); } - public long readSigned() throws IOException - { - long v = readUnsigned(); - return ( v & 1 ) == 0 ? v >> 1 : -(v >> 1 ); - } - - public long readUnsigned() throws IOException - { - long v = 0; - int shift = 0; - for(;;) - { - long i7 = readByte() & 0xff; - v |= (( i7 & 0x7f ) << shift); - if ( ( i7 & 0x80 ) == 0 ) break; - shift += 7; - } - return v; - } - - public long readMixed() throws IOException + public int readMixed() throws IOException { if ( repCount == 0 ) { - long b = readByte() & 0xff; - long repCode = b >> 6; - long diffcode = b & 0x3f; - repCount = repCode == 0 ? readUnsigned() : repCode; - lastValue += diffcode == 0 ? readSigned() : diffcode - 32; + boolean negative = decodeBit(); + int d = decodeVarBits() + diffshift; + repCount = decodeVarBits() + 1; + lastValue += negative ? -d : d; + diffshift = 1; } repCount--; return lastValue; } + public final boolean decodeBit() throws IOException + { + if ( bm == 0x100 ) + { + bm = 1; + b = readByte(); + } + boolean value = ( ( b & bm ) != 0 ); + bm <<= 1; + return value; + } + + public final int decodeVarBits() throws IOException + { + int range = 0; + int value = 0; + while (!decodeBit()) + { + value += range + 1; + range = 2 * range + 1; + } + return value + decodeBounded( range ); + } + + + public final int decodeBounded( int max ) throws IOException + { + int value = 0; + int im = 1; // integer mask + while (( value | im ) <= max) + { + if ( bm == 0x100 ) + { + bm = 1; + b = readByte(); + } + if ( ( b & bm ) != 0 ) + value |= im; + bm <<= 1; + im <<= 1; + } + return value; + } + } diff --git a/brouter-util/src/main/java/btools/util/MixCoderDataOutputStream.java b/brouter-util/src/main/java/btools/util/MixCoderDataOutputStream.java index b4a5469..ad6d81a 100644 --- a/brouter-util/src/main/java/btools/util/MixCoderDataOutputStream.java +++ b/brouter-util/src/main/java/btools/util/MixCoderDataOutputStream.java @@ -12,54 +12,41 @@ import java.io.OutputStream; public final class MixCoderDataOutputStream extends DataOutputStream { - private long lastValue; - private long lastLastValue; - private long repCount; - private boolean doFlush; + private int lastValue; + private int lastLastValue; + private int repCount; + private int diffshift; + + private int bm = 1; // byte mask (write mode) + private int b = 0; + + public static int[] diffs = new int[100]; + public static int[] counts = new int[100]; public MixCoderDataOutputStream( OutputStream os ) { super( os ); } - public void writeSigned( long v ) throws IOException - { - writeUnsigned( v < 0 ? ( (-v) << 1 ) | 1 : v << 1 ); - } - - public void writeUnsigned( long v ) throws IOException - { - if ( v < 0 ) throw new IllegalArgumentException( "writeUnsigned: " + v ); - do - { - long i7 = v & 0x7f; - v >>= 7; - if ( v != 0 ) i7 |= 0x80; - writeByte( (byte)( i7 & 0xff ) ); - } - while( v != 0 ); - } - - public void writeMixed( long v ) throws IOException + public void writeMixed( int v ) throws IOException { if ( v != lastValue && repCount > 0 ) { - long d = lastValue - lastLastValue; + int d = lastValue - lastLastValue; lastLastValue = lastValue; - - // if diff fits within 6 bits and rep-count < 4, write a single byte - int repCode = repCount < 4 ? (int)repCount : 0; - int diffcode = (int)(d > -32 && d < 32 ? d+32 : 0); - - writeByte( (byte)( diffcode | repCode << 6 ) ); - if ( repCode == 0) + + encodeBit( d < 0 ); + if ( d < 0 ) { - writeUnsigned( repCount ); - } - if ( diffcode == 0) - { - writeSigned( d ); + d = -d; } + encodeVarBits( d-diffshift ); + encodeVarBits( repCount-1 ); + + if ( d < 100 ) diffs[d]++; + if ( repCount < 100 ) counts[repCount]++; + + diffshift = 1; repCount = 0; } lastValue = v; @@ -69,11 +56,68 @@ public final class MixCoderDataOutputStream extends DataOutputStream @Override public void flush() throws IOException { - // todo: does this keep stream consistency after flush ? - long v = lastValue; + int v = lastValue; writeMixed( v+1 ); lastValue = v; repCount = 0; + if ( bm > 1 ) + { + writeByte( (byte)b ); // flush bit-coding + } } + public final void encodeBit( boolean value ) throws IOException + { + if ( bm == 0x100 ) + { + writeByte( (byte)b ); + bm = 1; + b = 0; + } + if ( value ) + { + b |= bm; + } + bm <<= 1; + } + + public final void encodeVarBits( int value ) throws IOException + { + int range = 0; + while (value > range) + { + encodeBit( false ); + value -= range + 1; + range = 2 * range + 1; + } + encodeBit( true ); + encodeBounded( range, value ); + } + + public final void encodeBounded( int max, int value ) throws IOException + { + int im = 1; // integer mask + while (im <= max) + { + if ( bm == 0x100 ) + { + writeByte( (byte)b ); + bm = 1; + b = 0; + } + if ( ( value & im ) != 0 ) + { + b |= bm; + max -= im; + } + bm <<= 1; + im <<= 1; + } + } + + public static void stats() + { + for(int i=1; i<100;i++) System.out.println( "diff[" + i + "] = " + diffs[i] ); + for(int i=1; i<100;i++) System.out.println( "counts[" + i + "] = " + counts[i] ); + } } diff --git a/brouter-util/src/main/java/btools/util/ReducedMedianFilter.java b/brouter-util/src/main/java/btools/util/ReducedMedianFilter.java index 1b7868a..7de4f08 100644 --- a/brouter-util/src/main/java/btools/util/ReducedMedianFilter.java +++ b/brouter-util/src/main/java/btools/util/ReducedMedianFilter.java @@ -24,6 +24,14 @@ public final class ReducedMedianFilter { if ( weight > 0. ) { + for( int i=0; i