diff --git a/pom.xml b/pom.xml index 519a379..a93ef5e 100644 --- a/pom.xml +++ b/pom.xml @@ -25,6 +25,11 @@ junit 4.1 + + com.vividsolutions + jts + 1.13 + diff --git a/src/main/java/de/rondiplomatico/nds/NDSBBox.java b/src/main/java/de/rondiplomatico/nds/NDSBBox.java index 044a83a..82c8b71 100644 --- a/src/main/java/de/rondiplomatico/nds/NDSBBox.java +++ b/src/main/java/de/rondiplomatico/nds/NDSBBox.java @@ -30,6 +30,14 @@ public class NDSBBox { private final int south; private final int west; + // constructor to get global bounding box in NDS coordinates + public NDSBBox() { + north = 90; + east = 180; + south = -90; + west = -180; + } + /** * * Gets the south west corner of the bounding box diff --git a/src/main/java/de/rondiplomatico/nds/NDSCoordinate.java b/src/main/java/de/rondiplomatico/nds/NDSCoordinate.java index 343468f..eae0140 100644 --- a/src/main/java/de/rondiplomatico/nds/NDSCoordinate.java +++ b/src/main/java/de/rondiplomatico/nds/NDSCoordinate.java @@ -96,6 +96,21 @@ public NDSCoordinate(long ndsMortonCoordinates) { latitude = lat; longitude = lon; } + + /* + * Compute midpoint of this and given NDSCoordinate + * + * @param p + * the given coordinate + */ + public NDSCoordinate getMidpoint(NDSCoordinate p) { + double lat = p.toWGS84().getLatitude(); + double lon = p.toWGS84().getLongitude(); + WGS84Coordinate lonlat = toWGS84(); + double longitude = lonlat.getLongitude(); + double latitude = lonlat.getLatitude(); + return new NDSCoordinate(longitude + 0.5 * (lon - longitude), latitude + 0.5 * (lat - latitude)); + } private void verify(int lon, int lat) { if (lat < MIN_LATITUDE || MAX_LATITUDE < lat) { diff --git a/src/main/java/de/rondiplomatico/nds/NDSEnvelope.java b/src/main/java/de/rondiplomatico/nds/NDSEnvelope.java new file mode 100644 index 0000000..eb1a834 --- /dev/null +++ b/src/main/java/de/rondiplomatico/nds/NDSEnvelope.java @@ -0,0 +1,86 @@ +package de.rondiplomatico.nds; + +import com.vividsolutions.jts.geom.Coordinate; +import com.vividsolutions.jts.geom.OctagonalEnvelope; + +/** + * Class to get a vividsolutions envelope of coordinates. + * + * No warranties for correctness, use at own risk. + * + * @author Andreas Hessenthaler + * @since 13.02.2020 + */ +public class NDSEnvelope { + + private static OctagonalEnvelope vividEnvelope = new OctagonalEnvelope(); + + public NDSEnvelope(double[][] polygonCoordinates) { + for(int i = 0; i < polygonCoordinates.length; i++) { + Coordinate vividCoord = new Coordinate(polygonCoordinates[i][0], polygonCoordinates[i][1]); + vividEnvelope.expandToInclude(vividCoord); + } + } + + public NDSEnvelope(OctagonalEnvelope e) { + vividEnvelope = e; + } + + public OctagonalEnvelope getEnvelope() { + return vividEnvelope; + } + + public NDSCoordinate getSouthWest() { + return new NDSCoordinate(vividEnvelope.getMinX(), vividEnvelope.getMinY()); + } + + public NDSCoordinate getSouthEast() { + return new NDSCoordinate(vividEnvelope.getMaxX(), vividEnvelope.getMinY()); + } + + public NDSCoordinate getNorthEast() { + return new NDSCoordinate(vividEnvelope.getMinX(), vividEnvelope.getMaxY()); + } + + public NDSCoordinate getNorthWest() { + return new NDSCoordinate(vividEnvelope.getMaxX(), vividEnvelope.getMaxY()); + } + + public int[] getMasterTileInfo(int maxLevels) { + NDSCoordinate point0 = getSouthWest(); + NDSCoordinate point1 = getSouthEast(); + NDSCoordinate point2 = getNorthEast(); + NDSCoordinate point3 = getNorthWest(); + return getMasterTileInfo(point0, point1, point2, point3, maxLevels); + } + + public int[] getMasterTileInfo(NDSCoordinate point0, NDSCoordinate point1, NDSCoordinate point2, NDSCoordinate point3, int maxLevels) { + int masterTileLevel = -1; + int masterTileID = -1; + for (int li = 0; li < maxLevels; li++) { + NDSTile currTile0 = new NDSTile(li, point0); + NDSTile currTile1 = new NDSTile(li, point1); + NDSTile currTile2 = new NDSTile(li, point2); + NDSTile currTile3 = new NDSTile(li, point3); + int currTileID0 = currTile0.getTileNumber(); + int currTileID1 = currTile1.getTileNumber(); + int currTileID2 = currTile2.getTileNumber(); + int currTileID3 = currTile3.getTileNumber(); + boolean singleTileID = (currTileID0 == currTileID1) && (currTileID0 == currTileID2) && (currTileID0 == currTileID3); + // if at least one tile ID is different, we discard the tile IDs and keep our previously detected tile ID (i.e. on the previous level) + if(!singleTileID) { + break; + } + // store tile info + masterTileLevel = li; + masterTileID = currTileID0; + } + // check if valid result + if (masterTileID == -1) { + System.out.println(">>>ERROR: Invalid master tile ID."); + System.exit(1); + } + int[] masterTileInfo = {masterTileLevel, masterTileID}; + return masterTileInfo; + } +} \ No newline at end of file diff --git a/src/main/java/de/rondiplomatico/nds/NDSTile.java b/src/main/java/de/rondiplomatico/nds/NDSTile.java index b3e1c2a..56a444b 100644 --- a/src/main/java/de/rondiplomatico/nds/NDSTile.java +++ b/src/main/java/de/rondiplomatico/nds/NDSTile.java @@ -1,5 +1,17 @@ package de.rondiplomatico.nds; +import java.awt.image.BufferedImage; +import java.io.File; +import java.io.IOException; +import java.util.ArrayList; +import java.util.Arrays; +import java.util.Collections; +import java.util.HashMap; +import java.util.List; +import java.util.Map; + +import javax.imageio.ImageIO; + import lombok.EqualsAndHashCode; import lombok.Getter; import lombok.ToString; @@ -11,6 +23,7 @@ * No warranties for correctness, use at own risk. * * @author Daniel Wirtz + * @author Andreas Hessenthaler * @since 20.02.2020 */ @Getter @@ -108,6 +121,256 @@ public NDSTile(int level, WGS84Coordinate coord) { this(level, new NDSCoordinate(coord.getLongitude(), coord.getLatitude())); } + /** + * Get tile indices in x / y from id and level. + */ + public int[] getTileXY() { + // variables for storing x/y-tile indices + int tileX, tileY; + tileX = tileY = -1; + // get tile size, assuming that there are 2x1 tiles on level 0 + double tileSizeX, tileSizeY; + tileSizeX = 360.0 / (double) (Math.pow(2, level+1)); + tileSizeY = 180.0 / (double) (Math.pow(2, level )); + // get tile center and southwest corner + NDSCoordinate c = getCenter(); + double lon = c.toWGS84().getLongitude() - 0.5 * tileSizeX; + double lat = c.toWGS84().getLatitude() - 0.5 * tileSizeY; + // compute tile indices + tileX = (int) Math.round((lon + 180.0) / tileSizeX); + tileY = (int) Math.round((lat + 90.0) / tileSizeY); + return new int[] {tileX, tileY}; + } + + /** + * Get tile indices in x / y from given id and level. + * + * @param level + * the level + * @param nr + * the tile number + */ + public int[] getTileXYfromTileNumber(int level, int nr) { + NDSTile tile = new NDSTile(level, nr); + return tile.getTileXY(); + } + + /** + * Get tile number from indices in x / y and level. + * + * @param level + * the level + * @param tileX + * the tile x-coordinate + * @param tileY + * the tile y-coordinate + */ + public int getTileNumberfromTileXY(int level, int tileX, int tileY) { + // get tile size, assuming that there are 2x1 tiles on level 0 + double tileSizeX = 360.0 / (double) (Math.pow(2, level+1)); + double tileSizeY = 180.0 / (double) (Math.pow(2, level )); + // get center + double clon = tileX * tileSizeX + 0.5 * tileSizeX - 180.0; + double clat = tileY * tileSizeY + 0.5 * tileSizeY - 90.0; + NDSCoordinate c = new NDSCoordinate(clon, clat); + // get tile + NDSTile t = new NDSTile(level, c); + // return tile number + return t.getTileNumber(); + } + + /** + * Get tile quadkey from tile indices in x / y and level. + * + * @param level + * the level + * @param tileX + * the tile x-coordinate + * @param tileY + * the tile y-coordinate + */ + public String getTileQuadkeyFromTileXY(int level, int tileX, int tileY) { + // setup object for building quadkey + StringBuilder quadkey = new StringBuilder(); + // loop over levels + for (int i = level+1; i > 0; i--) { + char digit = '0'; + int mask = 1 << (i - 1); + if ((tileX & mask) != 0) { + digit++; + } + if ((tileY & mask) != 0) { + digit++; + digit++; + } + // add to quadkey + quadkey.append(digit); + } + // return quadkey + return quadkey.toString(); + } + + /** + * Get tile indices in x / y from tile quadkey. + * + * @param quadkey + * the quadkey + */ + public int[] getTileXYfromTileQuadkey(String quadkey) { + int[] tileXY = new int[2]; + // 0-based level index + int level = quadkey.length(); + for (int i = level; i > 0; i--) { + int mask = 1 << (i - 1); + switch (quadkey.charAt(level - i)) { + case '0': + break; + case '1': + tileXY[0] |= mask; + break; + case '2': + tileXY[1] |= mask; + break; + case '3': + tileXY[0] |= mask; + tileXY[1] |= mask; + break; + default: + System.out.println(">>>ERROR:"); + throw new IllegalArgumentException("Invalid quadkey digit sequence."); + } + } + return tileXY; + } + + /** + * Print map with labeled tiles. + * + * @param level + * the level + * @param tiles + * the relevant tiles + * @param printType + * the print type (e.g., png or text) + * @param fileName + * the output filename + */ + public void printMap(int level, int[] tiles, String printType, String fileName) { + // get number of tiles in x/y + int numTilesX = (int) Math.pow(2, level+1); + int numTilesY = (int) Math.pow(2, level ); + // print labeled cells + if (printType.equalsIgnoreCase("text")) { + if (level > 4) { + System.out.println(">>>WARNING: Writing to console won't give readable output. Skipping.."); + return; + } + String header = new String(new char[(numTilesX+1)*2-1]).replace('\0', '-'); + for (int idxY = numTilesY-1; idxY >= 0; idxY--) { + System.out.println(header); + for (int idxX = 0; idxX < numTilesX; idxX++) { + int tileID = getTileNumberfromTileXY(level, idxX, idxY); + // print label + // note: if check is quite inefficient.. + if (Arrays.stream(tiles).anyMatch(i -> i == tileID)) { + System.out.print("|X"); + } else { + System.out.print("| "); + } + } + System.out.println("|"); + } + System.out.println(header); + } else if (printType.equalsIgnoreCase("png")) { + if (level > 13) { + System.out.println(">>>WARNING: Required heap space is "+((long)numTilesX*(long)numTilesY*(long)4/(long)1024/(long)1024)+" MB."); + System.out.println(">>>WARNING: Writing to buffered image may fail due to limited heap space. Skipping.."); + return; + } + // create image using default value for type int: 0 + BufferedImage image = new BufferedImage(numTilesX, numTilesY, BufferedImage.TYPE_INT_RGB); + for (int idx = 0; idx < tiles.length; idx++) { + int[] tileXY = getTileXYfromTileNumber(level, tiles[idx]); + image.setRGB(tileXY[0], numTilesY-1-tileXY[1], -16000); + } + File ImageFile = new File(System.getProperty("user.home"), fileName+".png"); + try { + ImageIO.write(image, "png", ImageFile); + } catch (IOException e) { + e.printStackTrace(); + } + } + return; + } + + /** + * Compute a hash map from tile numbers. + * + * @param level + * the level + * @param tileNumbers + * the tileNumbers + * @return + */ + public Map> tileNumbersToHM(int level, int[] tileNumbers){ + // get the master tile number 0 + NDSTile masterTile = new NDSTile(0, 0); + // create hashmap of (key, value) pairs, where key is tileY and val is tileX + // note: val may contain multiple values + Map> tileHM = new HashMap>(); + for (int ti = 0; ti < tileNumbers.length; ti++) { + int[] tileXY = getTileXYfromTileNumber(level, tileNumbers[ti]); + int key = tileXY[1]; + int newVal = tileXY[0]; + // if key already exists, add value to sorted list; else create new list + if (tileHM.containsKey(key)) { + List prevList = tileHM.get(key); + prevList.add(newVal); + Collections.sort(prevList); + tileHM.put(key, prevList); + } else { + List newList = new ArrayList(); + newList.add(newVal); + tileHM.put(key, newList); + } + } + return tileHM; + } + + /** + * Compute tile numbers from hash map. + * + * @param level + * the level + * @param hm + * the hash map + * @return + */ + public int[] hmToTileNumbers(int level, Map> hm) { + NDSTile masterTile = new NDSTile(0, 0); + int numVals = getNumberOfValuesHM(hm); + int[] filledTileIDs = new int[numVals]; + int idx = 0; + for (Map.Entry> entry : hm.entrySet()) { + int key = entry.getKey(); + List currList = entry.getValue(); + for (int idx0 = 0; idx0 < currList.size(); idx0++) { + filledTileIDs[idx] = getTileNumberfromTileXY(level, currList.get(idx0), key); + idx++; + } + } + return filledTileIDs; + } + + + private int getNumberOfValuesHM(Map> hm) { + int numVals = 0; + for (Map.Entry> entry : hm.entrySet()) { + numVals = numVals + entry.getValue().size(); + } + return numVals; + } + /** * Checks if the current Tile contains a certain coordinate. * @@ -155,6 +418,41 @@ public NDSCoordinate getCenter() { } return center; } + + /* + * Returns the corners of the tile + * + * @param masterTile + * the masterTile + */ + public NDSCoordinate[] getCorners() { + NDSBBox bb = getBBox(); + return new NDSCoordinate[] {bb.southWest(), bb.southEast(), bb.northEast(), bb.northWest()}; + } + + /* + * Return the child tile numbers + */ + public int[] getChildTileNumbers() { + int id0 = (tileNumber << 2); + return new int[] {id0, id0+1, id0+3, id0+2}; + } + + public int getChildTileNumberSouthWest() { + return (tileNumber << 2); + } + + public int getChildTileNumberSouthEast() { + return (tileNumber << 2) + 1; + } + + public int getChildTileNumberNorthEast() { + return (tileNumber << 2) + 3; + } + + public int getChildTileNumberNorthWest() { + return (tileNumber << 2) + 2; + } /** * Creates a bounding box for the current tile. diff --git a/src/main/java/de/rondiplomatico/nds/floodFill.java b/src/main/java/de/rondiplomatico/nds/floodFill.java new file mode 100644 index 0000000..bf89250 --- /dev/null +++ b/src/main/java/de/rondiplomatico/nds/floodFill.java @@ -0,0 +1,271 @@ +package de.rondiplomatico.nds; + +import java.awt.Point; +import java.util.LinkedList; +import java.util.Queue; + +/** + * Utility class for flood fill algorithm applied to marked tiles deduced from polygon points. + * + * @author Andreas Hessenthaler + * @since 13.02.2020 + */ +public class floodFill { + int[][] ff; + int minX = Integer.MAX_VALUE; + int maxX = Integer.MIN_VALUE; + int minY = Integer.MAX_VALUE; + int maxY = Integer.MIN_VALUE; + int dimX = -1; + int dimY = -1; + int bg = 0; + int bound = 1; + int fill = bound; + int startAtX = 0; + int startAtY = 0; + + public floodFill(int[][] f, int a, int b, int c, int d) { + ff = f; + minX = a; + maxX = b; + minY = c; + maxY = d; + dimX = maxX - minX + 1; + dimY = maxY - minY + 1; + setStartingPoint(); + } + + /* + * Four-way recursive flood fill algorithm + * + * @param x + * the tile x coordinate + * @param y + * the tile y coordinate + */ + public void recursiveFloodFill(int x, int y) { + if ((x < 0) || (y < 0) || (x > dimX-1 || (y > dimY-1))) { + return; + } + if (ff[x][y] == bg) { + ff[x][y] = fill; + recursiveFloodFill(x+1, y ); // right + recursiveFloodFill(x, y+1); // top + recursiveFloodFill(x-1, y ); // left + recursiveFloodFill(x, y-1); // bottom + } + return; + } + + /* + * Eight-way recursive flood fill algorithm + * + * @param x + * the tile x coordinate + * @param y + * the tile y coordinate + */ + public void recursiveFloodFill8(int x, int y) { + if ((x < 0) || (y < 0) || (x > dimX-1 || (y > dimY-1))) { + return; + } + if (ff[x][y] == bg) { + ff[x][y] = fill; + recursiveFloodFill(x+1, y ); // right + recursiveFloodFill(x, y+1); // top + recursiveFloodFill(x-1, y ); // left + recursiveFloodFill(x, y-1); // bottom + recursiveFloodFill(x+1, y+1); // top right + recursiveFloodFill(x-1, y+1); // top left + recursiveFloodFill(x+1, y-1); // bottom right + recursiveFloodFill(x-1, y-1); // bottom left + } + return; + } + + /* + * Iterative four-way flood fill algorithm - avoids stack overflow issues + * + */ + public void iterativeFloodFill() { + Queue q = new LinkedList(); + q.add(new Point(startAtX, startAtY)); + while (!q.isEmpty()) { + Point p = q.remove(); + if ((p.x < 0) || (p.y < 0) || (p.x > dimX-1 || (p.y > dimY-1))) { + continue; + } + if (ff[p.x][p.y] == bg) { + ff[p.x][p.y] = fill; + q.add(new Point(p.x+1, p.y )); // right + q.add(new Point(p.x, p.y+1)); // top + q.add(new Point(p.x-1, p.y )); // left + q.add(new Point(p.x, p.y-1)); // bottom + } + } + } + + /* + * Iterative eight-way flood fill algorithm - avoids stack overflow issues + * + */ + public void iterativeFloodFill8() { + Queue q = new LinkedList(); + q.add(new Point(startAtX, startAtY)); + while (!q.isEmpty()) { + Point p = q.remove(); + if ((p.x < 0) || (p.y < 0) || (p.x > dimX-1 || (p.y > dimY-1))) { + continue; + } + if (ff[p.x][p.y] == bg) { + ff[p.x][p.y] = fill; + q.add(new Point(p.x+1, p.y )); // right + q.add(new Point(p.x, p.y+1)); // top + q.add(new Point(p.x-1, p.y )); // left + q.add(new Point(p.x, p.y-1)); // bottom + q.add(new Point(p.x+1, p.y+1)); // top right + q.add(new Point(p.x-1, p.y+1)); // top left + q.add(new Point(p.x+1, p.y-1)); // bottom right + q.add(new Point(p.x-1, p.y-1)); // bottom left + } + } + } + + /* + * Safety net for filling left over lines after applying standard flood fill algorithm + * + */ + public void fillHoles() { + // first, find vertical lines + boolean vline = false; + boolean vlineEnded = false; + int currStartY = -1; + int currStopY = -1; + for (int col = 1; col < dimY-1; col++) { + vline = false; + vlineEnded = false; + for (int row = 1; row < dimX-1; row++) { + // minimum requirement: left and right neighbors are not background + if ((ff[row][col-1] != bg) && (ff[row][col+1] != bg)) { + if ((ff[row][col] == bg) + && (ff[row-1][col] != bg)) { + // beginning of vertical line + currStartY = row; + vline = true; + } else if (vline && (ff[row][col] == bg)) { + // vertical line continues + } else if (vline && (ff[row][col] != bg)) { + // end of vertical line + currStopY = row; + vlineEnded = true; + } + } else if (vline && ((ff[row][col+1] != bg) || (ff[row][col-1] != bg)) && (ff[row][col] != bg)) { + // end of horizontal line + currStopY = row; + vlineEnded = true; + } else { + vline = false; + vlineEnded = false; + } + if (vline && vlineEnded) { + // fill line and reset booleans + System.out.println(">>>>>>INFO: Filling vertical line."); + for (int idx = currStartY; idx < currStopY; idx++) { + ff[idx][col] = fill; + } + vline = false; + vlineEnded = false; + } + } + } + // second, find horizontal lines + boolean hline = false; + boolean hlineEnded = false; + for (int row = 1; row < dimX-1; row++) { + hline = false; + hlineEnded = false; + for (int col = 1; col < dimY-1; col++) { + // minimum requirement: top and bottom neighbors are not background + if ((ff[row+1][col] != bg) && (ff[row-1][col] != bg)) { + if ((ff[row][col] == bg) + && (ff[row][col-1] != bg)) { + // beginning of horizontal line + currStartY = col; + hline = true; + } else if (hline && (ff[row][col] == bg)) { + // vertical line continues + } else if (hline && (ff[row][col] != bg)) { + // end of vertical line + currStopY = col; + hlineEnded = true; + } + } else if (hline && ((ff[row+1][col] != bg) || (ff[row-1][col] != bg)) && (ff[row][col] != bg)) { + // end of vertical line + currStopY = col; + hlineEnded = true; + } else { + hline = false; + hlineEnded = false; + } + if (hline && hlineEnded) { + // fill line and reset booleans + System.out.println(">>>>>>INFO: Filling horizontal line."); + for (int idx = currStartY; idx < currStopY; idx++) { + ff[row][idx] = fill; + } + hline = false; + hlineEnded = false; + } + } + } + + return; + } + + /* + * Set initial point for flood fill algorithm + * + * @todo need to do more extensive checks if we can always find an initial point if we only check the right half of the box + */ + public void setStartingPoint() { + startAtX = (int) (0.5 * dimX); + startAtY = 0; + // we're starting on background, skip that.. + while (ff[startAtX][startAtY] == bg) { + startAtY++; + } + // we're continuing on an edge, skip that.. + while (ff[startAtX][startAtY] == bound) { + startAtY++; + } + for (int row = startAtY; row < dimY; row++, startAtY++) { + int rowSum = 0; + for (int col = startAtX; col < dimX; col++) { + if (ff[col][row-1] != ff[col][row]) { + rowSum = rowSum + ff[col][row]; + } + } + if ((rowSum % 2) == 1) { + // we've found a point that's inside + break; + } + } + return; + } + + /* + * Get initial point for flood fill algorithm + * + */ + public int[] getStartingPoint() { + return new int[] {startAtX, startAtY}; + } + + /* + * Get flood filled image + * + */ + public int[][] getFloodFill(){ + return ff; + } +} diff --git a/src/main/java/de/rondiplomatico/nds/polygon.java b/src/main/java/de/rondiplomatico/nds/polygon.java new file mode 100644 index 0000000..8865b73 --- /dev/null +++ b/src/main/java/de/rondiplomatico/nds/polygon.java @@ -0,0 +1,264 @@ +package de.rondiplomatico.nds; + +import java.util.ArrayList; +import java.util.HashMap; +import java.util.List; +import java.util.Map; +import java.util.stream.IntStream; + +public class polygon { + + /** + * Main method to test selecting all tiles covered by a sample polygon. + * + * @param args + */ + public static void main(String[] args) { + + long t0 = System.currentTimeMillis(); + + // get some random polygon data for testing purposes +// double[][] polygonCoordinates = new double[][] { +// {10.5, 45.9}, +// {13.0, 50.3}, +// {15.0, 47.0}, +// {13.4, 70.0}, +// {10.5, 45.9} +// }; + + // get some random polygon data for testing purposes +// double[][] polygonCoordinates = new double[][] { +// {10.5, 45.9}, +// {13.0, 63.3}, +// {15.0, 47.0}, +// {13.4, 70.0}, +// {10.5, 45.9} +// }; + + // get some random polygon data for testing purposes +// double[][] polygonCoordinates = new double[][] { +// {10.5, 45.9}, +// {13.0, 50.3}, +// {15.0, 47.0}, +// {17.0, 50.3}, +// {20.0, 47.0}, +// {13.4, 60.0}, +// {13.4, 70.0}, +// {10.5, 45.9} +// }; + + // get some random polygon data approximating Germany + // https://www.mapsofworld.com/lat_long/germany-lat-long.html + double[][] polygonCoordinates = new double[][] { + {10.5, 45.9}, + {13.0, 45.9}, + {14.0, 49.0}, + {12.0, 50.0}, + {15.0, 51.0}, + {15.0, 54.0}, + {13.5, 54.5}, + {11.0, 54.0}, + {10.0, 55.0}, + { 8.5, 55.0}, + { 9.0, 54.0}, + { 7.0, 53.5}, + { 6.0, 52.0}, + { 6.1, 50.0}, + { 8.0, 49.0}, + { 7.5, 47.5}, + {10.5, 45.9} + }; + + // number of levels in map hierarchy + int maxLevels = 15; + // get a bounding octagonal envelope (defaults to quadrilateral in 2D case) + // for sample data corresponding to a polygon (e.g. borders of a country) + NDSEnvelope envelope = new NDSEnvelope(polygonCoordinates); + // get corresponding bounding tile ID, i.e. find level where all polygon points are on the same tile + int[] masterTileInfo = envelope.getMasterTileInfo(maxLevels); + int masterTileLevel = masterTileInfo[0]; + int masterTileID = masterTileInfo[1]; + // store the master tile + NDSTile masterTile = new NDSTile(masterTileLevel, masterTileID); + // get all tiles covered by the polygon on the tstLevel + int tstLevel = 8; // 13 + // refine polygon coordinates to avoid edges that are crossing a tile + // set refinement factor to -1 to adaptively refine + int numSamples = -1; + double[][] polygonCoordinatesRef = polygon.refinePolygon(tstLevel, polygonCoordinates, numSamples); + // let's grab all tiles with polygon points + int[] uniqueTileIDs = polygon.getUniqueTileNumbersOnLevel(tstLevel, polygonCoordinatesRef); + // dump to image file for debugging + masterTile.printMap(tstLevel, uniqueTileIDs, "png", "map"); + // get x/y tile indices for tiles with polygon points + // key is tileY, val is tileX (note: val may contain multiple tileX indices) + Map> tileHM = masterTile.tileNumbersToHM(tstLevel, uniqueTileIDs); + // dump to image file for debugging + Map> filledTileHM = polygon.mapFillPolygon(tileHM, uniqueTileIDs, "flood-fill-safe"); + int[] filledTileIDs = masterTile.hmToTileNumbers(tstLevel, filledTileHM); + masterTile.printMap(tstLevel, filledTileIDs, "png", "mapFilled"); + + long t1 = System.currentTimeMillis(); + System.out.println("\n>>>INFO: Program finished in "+(t1-t0)+" ms."); + } + + /** + * Applies a flood-fill algorithm to a polygon + * + * @param hmi + * the hash map input hmi + * @param tileNumbers + * the tileNumbers + * @param method + * the flood-fill algorithm (flood-fill, row-fill) + * @return + */ + public static Map> mapFillPolygon(Map> hmi, int[] tileNumbers, String method){ + Map> hmo = new HashMap>(); + if (method.contains("flood-fill")) { + // first get a minimum 'bounding tile box' + int minX = Integer.MAX_VALUE; + int maxX = Integer.MIN_VALUE; + int minY = Integer.MAX_VALUE; + int maxY = Integer.MIN_VALUE; + // loop over (key, values) pairs to get boundary + for (Map.Entry> entry : hmi.entrySet()) { + int key = entry.getKey(); + List currList = entry.getValue(); + // we have a sorted list, so we can just get the first and last element for the comparison + minX = Math.min(currList.get(0), minX); + maxX = Math.max(currList.get(currList.size()-1), maxX); + // key is the y-index + minY = Math.min(key, minY); + maxY = Math.max(key, maxY); + } + // get dimensions + int dimX = maxX - minX + 1; + int dimY = maxY - minY + 1; + // allocate int[][] for flood fill algorithm; default value 0 + int[][] tmp = new int[dimX][dimY]; + // loop over (key, values) pairs to label boundary + for (Map.Entry> entry : hmi.entrySet()) { + int key = entry.getKey() - minY; + List currList = entry.getValue(); + // loop over values to label boundary + for (int idx0 = 0; idx0 < currList.size(); idx0++) { + int val = currList.get(idx0) - minX; + tmp[val][key] = 1; + } + } + // perform flood fill from midpoint + + floodFill ff = new floodFill(tmp, minX, maxX, minY, maxY); + ff.iterativeFloodFill(); + tmp = ff.getFloodFill(); + if (method.equalsIgnoreCase("flood-fill-safe")) { + ff.fillHoles(); + } + // loop over y + for (int key = 0; key < dimY; key++) { + List currList = new ArrayList(); + // loop over x + for (int x = 0; x < dimX; x++) { + // get all elements that are on boundary or inside + if (tmp[x][key] > 0) { + currList.add(x+minX); + } + } + hmo.put(key+minY, currList); + } + } else if (method.contains("row-fill")) { + System.out.println(">>>WARNING: Line-filling may not respect boundaries of concave shape. Use flood fill algorithm instead."); + for (Map.Entry> entry : hmi.entrySet()) { + int key = entry.getKey(); + List currList = entry.getValue(); + for (int idx0 = 0; idx0 < currList.size()-1; idx0++) { + int currVal = currList.get(idx0); + int nextVal = currList.get(idx0+1); + if (currVal+1 < nextVal) { + currList.add(idx0+1, currVal+1); + } + } + hmo.put(key, currList); + } + } else { + System.out.println(">>>ERROR: Unknown method "+method); + System.exit(1); + } + return hmo; + } + + /** + * Refine polygon data by subsampling. + * + * @param level + * the level + * @param polygonCoordinates + * the polygon coordinates + * @param numSamples + * number of samples between each polygon edge (-1 requests adaptive sampling) + * @return polygonCoordinatesRef + * the refined polygon coordinates + */ + public static double[][] refinePolygon(int level, double[][] polygonCoordinates, int numSamples) { + if (numSamples == 0) { + System.out.println(">>>INFO: Polygon refinement off. Consider switching to adaptive refinement."); + return polygonCoordinates; + } + if (numSamples == -1) { + double maxDist = 0.0; + double tileSizeX = 360.0 / (double) (Math.pow(2, level+1)); + for (int idx = 0; idx < polygonCoordinates.length-1; idx++) { + double p0x = polygonCoordinates[idx ][0]; + double p0y = polygonCoordinates[idx ][1]; + double p1x = polygonCoordinates[idx+1][0]; + double p1y = polygonCoordinates[idx+1][1]; + double dist = Math.sqrt(Math.pow(p1x-p0x, 2.0) + Math.pow(p1y-p0y, 2.0)); + maxDist = Math.max(dist, maxDist); + } + double maxTargetDist = tileSizeX * 0.4; + if (maxDist <= maxTargetDist) { + numSamples = 1; + } else { + numSamples = (int) Math.ceil(maxDist / maxTargetDist); + } + System.out.println(">>>INFO: Setting adaptive refinement to "+numSamples); + } + int numCoord = polygonCoordinates.length; + int numCoordRef = (numCoord - 1) * (numSamples + 1) + 1; + double[][] polygonCoordinatesRef = new double[numCoordRef][2]; + // sample between all points except for the last and first because those coincide + for (int idx = 0; idx < numCoord-1; idx++) { + double x0 = polygonCoordinates[idx ][0]; + double x1 = polygonCoordinates[idx+1][0]; + double y0 = polygonCoordinates[idx ][1]; + double y1 = polygonCoordinates[idx+1][1]; + for (int jdx = 0; jdx < numSamples+1; jdx++) { + polygonCoordinatesRef[idx*(numSamples+1)+jdx][0] = x0 + (x1 - x0) * (double)jdx / (double)(numSamples + 1); + polygonCoordinatesRef[idx*(numSamples+1)+jdx][1] = y0 + (y1 - y0) * (double)jdx / (double)(numSamples + 1); + } + } + // final point to close polygon + polygonCoordinatesRef[numCoordRef-1][0] = polygonCoordinates[numCoord-1][0]; + polygonCoordinatesRef[numCoordRef-1][1] = polygonCoordinates[numCoord-1][1]; + // return refined polygon + return polygonCoordinatesRef; + } + + public static int[] getTileNumbersOnLevel(int level, double[][] polygonCoordinates){ + int numCoord = polygonCoordinates.length; + int[] tileNumbers = new int[numCoord]; + // loop over all coordinates to get map tile IDs + for (int idx = 0; idx < numCoord; idx++) { + NDSCoordinate currCoord = new NDSCoordinate(polygonCoordinates[idx][0], polygonCoordinates[idx][1]); + NDSTile currTile = new NDSTile(level, currCoord); + tileNumbers[idx] = currTile.getTileNumber(); + } + return tileNumbers; + } + + public static int[] getUniqueTileNumbersOnLevel(int level, double[][] polygonCoordinates){ + int[] tileNumbers = getTileNumbersOnLevel(level, polygonCoordinates); + return IntStream.of(tileNumbers).distinct().toArray(); + } +} diff --git a/src/test/java/de/rondiplomatico/nds/NDSTileTest.java b/src/test/java/de/rondiplomatico/nds/NDSTileTest.java index 9d58ad6..0f78f21 100644 --- a/src/test/java/de/rondiplomatico/nds/NDSTileTest.java +++ b/src/test/java/de/rondiplomatico/nds/NDSTileTest.java @@ -6,6 +6,8 @@ import org.junit.Test; +import com.vividsolutions.jts.geom.OctagonalEnvelope; + import de.rondiplomatico.nds.NDSBBox; import de.rondiplomatico.nds.NDSCoordinate; import de.rondiplomatico.nds.NDSTile; @@ -298,6 +300,222 @@ public void testCenterWorks() { } } + @Test + public void testCornersWorks() { + // Tile 0 + NDSTile t = new NDSTile(0, 0); + NDSCoordinate c = t.getCenter(); + NDSCoordinate[] corners = t.getCorners(); + NDSCoordinate sw = corners[0]; + NDSCoordinate se = corners[1]; + NDSCoordinate ne = corners[2]; + NDSCoordinate nw = corners[3]; + System.out.println(sw.toWGS84()); + System.out.println(se.toWGS84()); + System.out.println(nw.toWGS84()); + System.out.println(ne.toWGS84()); + assertEquals( 0.0, sw.toWGS84().getLongitude(), eps); + assertEquals(180.0, se.toWGS84().getLongitude(), eps); + assertEquals( 0.0, nw.toWGS84().getLongitude(), eps); + assertEquals(180.0, ne.toWGS84().getLongitude(), eps); + assertEquals(-90.0, sw.toWGS84().getLatitude(), eps); + assertEquals(-90.0, se.toWGS84().getLatitude(), eps); + assertEquals( 90.0, nw.toWGS84().getLatitude(), eps); + assertEquals( 90.0, ne.toWGS84().getLatitude(), eps); + // Tile 1 + t = new NDSTile(0, 1); + c = t.getCenter(); + corners = t.getCorners(); + sw = corners[0]; + se = corners[1]; + ne = corners[2]; + nw = corners[3]; + System.out.println(sw.toWGS84()); + System.out.println(se.toWGS84()); + System.out.println(nw.toWGS84()); + System.out.println(ne.toWGS84()); + assertEquals(-180.0, sw.toWGS84().getLongitude(), eps); + assertEquals( 0.0, se.toWGS84().getLongitude(), eps); + assertEquals(-180.0, nw.toWGS84().getLongitude(), eps); + assertEquals( 0.0, ne.toWGS84().getLongitude(), eps); + assertEquals(-90.0, sw.toWGS84().getLatitude(), eps); + assertEquals(-90.0, se.toWGS84().getLatitude(), eps); + assertEquals( 90.0, nw.toWGS84().getLatitude(), eps); + assertEquals( 90.0, ne.toWGS84().getLatitude(), eps); + } + + @Test + public void testMidpointsWorks() { + // Tile 0 + NDSTile t = new NDSTile(0, 0); + NDSCoordinate c = t.getCenter(); + NDSCoordinate[] corners = t.getCorners(); + NDSCoordinate sw = corners[0]; + NDSCoordinate se = corners[1]; + NDSCoordinate ne = corners[2]; + NDSCoordinate nw = corners[3]; + NDSCoordinate csw = sw.getMidpoint(c); + NDSCoordinate cse = se.getMidpoint(c); + NDSCoordinate cne = ne.getMidpoint(c); + NDSCoordinate cnw = nw.getMidpoint(c); + assertEquals( 45.0, csw.toWGS84().getLongitude(), eps); + assertEquals(135.0, cse.toWGS84().getLongitude(), eps); + assertEquals( 45.0, cnw.toWGS84().getLongitude(), eps); + assertEquals(135.0, cne.toWGS84().getLongitude(), eps); + assertEquals(-45.0, csw.toWGS84().getLatitude(), eps); + assertEquals(-45.0, cse.toWGS84().getLatitude(), eps); + assertEquals( 45.0, cnw.toWGS84().getLatitude(), eps); + assertEquals( 45.0, cne.toWGS84().getLatitude(), eps); + // Tile 1 + t = new NDSTile(0, 1); + c = t.getCenter(); + corners = t.getCorners(); + sw = corners[0]; + se = corners[1]; + ne = corners[2]; + nw = corners[3]; + csw = sw.getMidpoint(c); + cse = se.getMidpoint(c); + cne = ne.getMidpoint(c); + cnw = nw.getMidpoint(c); + assertEquals(-135.0, csw.toWGS84().getLongitude(), eps); + assertEquals( -45.0, cse.toWGS84().getLongitude(), eps); + assertEquals(-135.0, cnw.toWGS84().getLongitude(), eps); + assertEquals( -45.0, cne.toWGS84().getLongitude(), eps); + assertEquals( -45.0, csw.toWGS84().getLatitude(), eps); + assertEquals( -45.0, cse.toWGS84().getLatitude(), eps); + assertEquals( 45.0, cnw.toWGS84().getLatitude(), eps); + assertEquals( 45.0, cne.toWGS84().getLatitude(), eps); + } + + @Test + public void testChildTileNumbersWorks() { + // get tile and its center + int masterTileLevel = 3; + NDSTile t = new NDSTile(masterTileLevel, 8); + NDSCoordinate c = t.getCenter(); + // store the master tile corners; ordering is: sw se ne nw + NDSCoordinate[] corners = t.getCorners(); + NDSCoordinate sw = corners[0]; + NDSCoordinate se = corners[1]; + NDSCoordinate ne = corners[2]; + NDSCoordinate nw = corners[3]; + // first, let's create four points that we know are in the child tiles + int childTileLevel = masterTileLevel + 1; + NDSCoordinate csw = sw.getMidpoint(c); + NDSCoordinate cse = se.getMidpoint(c); + NDSCoordinate cne = ne.getMidpoint(c); + NDSCoordinate cnw = nw.getMidpoint(c); + // get child tiles based on master level for double-checking + NDSTile tsw = new NDSTile(childTileLevel, csw); + NDSTile tse = new NDSTile(childTileLevel, cse); + NDSTile tne = new NDSTile(childTileLevel, cne); + NDSTile tnw = new NDSTile(childTileLevel, cnw); + int swNumber = t.getChildTileNumberSouthWest(); + int seNumber = t.getChildTileNumberSouthEast(); + int neNumber = t.getChildTileNumberNorthEast(); + int nwNumber = t.getChildTileNumberNorthWest(); + // now perform check + assertEquals(tsw.getTileNumber(), swNumber); + assertEquals(tse.getTileNumber(), seNumber); + assertEquals(tne.getTileNumber(), neNumber); + assertEquals(tnw.getTileNumber(), nwNumber); + } + + @Test + public void testQuadkeyWorks() { + // rough polygon data approximating Germany + double[][] polygonCoordinates = new double[][] { + {10.5, 45.9}, + {13.0, 45.9}, + {14.0, 49.0}, + {12.0, 50.0}, + {15.0, 51.0}, + {15.0, 54.0}, + {13.5, 54.5}, + {11.0, 54.0}, + {10.0, 55.0}, + { 8.5, 55.0}, + { 9.0, 54.0}, + { 7.0, 53.5}, + { 6.0, 52.0}, + { 6.1, 50.0}, + { 8.0, 49.0}, + { 7.5, 47.5}, + {10.5, 45.9} + }; + // number of levels in map hierarchy + int maxLevels = 15; + // get a bounding octagonal envelope (defaults to quadrilateral in 2D case) + // for sample data corresponding to a polygon (e.g. borders of a country) + NDSEnvelope envelope = new NDSEnvelope(polygonCoordinates); + // get corresponding bounding tile ID, i.e. find level where all polygon points are on the same tile + int[] masterTileInfo = envelope.getMasterTileInfo(maxLevels); + int masterTileLevel = masterTileInfo[0]; + int masterTileID = masterTileInfo[1]; + // store the master tile + NDSTile masterTile = new NDSTile(masterTileLevel, masterTileID); + // set some random tile number and level for testing + int tstTileNumber = 7; + int tstTileLevel = 3; + int[] tstTileXY = masterTile.getTileXYfromTileID(tstTileLevel, tstTileNumber); + int tstTileX = tstTileXY[0]; + int tstTileY = tstTileXY[1]; + assertEquals(11, tstTileX); + assertEquals( 5, tstTileY); + String tstTileQuadkey = masterTile.getTileQuadkeyFromTileXY(tstTileLevel, tstTileX, tstTileY); + assertEquals("1213", tstTileQuadkey); + assertEquals("00313102310", masterTile.getTileQuadkeyFromTileXY(10, 486, 332)); + assertEquals("01202102332221212", masterTile.getTileQuadkeyFromTileXY(16, 35210, 21493)); + // inverse direction + tstTileXY = masterTile.getTileXYfromTileQuadkey(tstTileQuadkey); + assertEquals(tstTileX, tstTileXY[0]); + assertEquals(tstTileY, tstTileXY[1]); + assertEquals(tstTileNumber, masterTile.getTileIDfromTileXY(tstTileLevel, tstTileX, tstTileY)); + } + + @Test + public void testEnvelopeWorks() { + // rough polygon data approximating Germany + double[][] polygonCoordinates = new double[][] { + {10.5, 45.9}, + {13.0, 45.9}, + {14.0, 49.0}, + {12.0, 50.0}, + {15.0, 51.0}, + {15.0, 54.0}, + {13.5, 54.5}, + {11.0, 54.0}, + {10.0, 55.0}, + { 8.5, 55.0}, + { 9.0, 54.0}, + { 7.0, 53.5}, + { 6.0, 52.0}, + { 6.1, 50.0}, + { 8.0, 49.0}, + { 7.5, 47.5}, + {10.5, 45.9} + }; + NDSEnvelope envelope = new NDSEnvelope(polygonCoordinates); + OctagonalEnvelope vividEnvelope = envelope.getEnvelope(); + // check envelope boundaries + assertEquals( 6.0, vividEnvelope.getMinX(), eps); + assertEquals(15.0, vividEnvelope.getMaxX(), eps); + assertEquals(45.9, vividEnvelope.getMinY(), eps); + assertEquals(55.0, vividEnvelope.getMaxY(), eps); + // check bounding box corners + assertEquals(581131592357515410L, envelope.getSouthWest().getMortonCode()); + assertEquals(595825689965249734L, envelope.getSouthEast().getMortonCode()); + assertEquals(592806849050488888L, envelope.getNorthEast().getMortonCode()); + assertEquals(607500946658223212L, envelope.getNorthWest().getMortonCode()); + // check master tile details + int[] info = envelope.getMasterTileInfo(15); + int level = info[0]; + int number = info[1]; + assertEquals(3, level); + assertEquals(8, number); + } + @Test public void testConstructor() { try {