package algorithmrepository; /** * Provides fast lookup capability for which of a list of polygons is at a given point * by dividing the full space into a grid and recording polygons intersecting grid cells. * * Probably only best for: * Relatively dense polygons per unit space, otherwise * some sort of search tree will be as fast with less memory * Low number of points per polygon (n < 20), otherwise line-line * intersections probably take longer than helpful. * * @author oliford * * Copyright Oliver Ford (c) 2010 * */ public class PolygonGridLookup { public static final int POLY_OUTSIDE_CELL = 1; public static final int POLY_INTERSECT_CELL = 2; public static final int POLY_ENCLOSES_CELL = 3; public static final int POLY_ENCLOSED_BY_CELL = 4; /** Amount by which to increment polyLookup array when more space is needed */ public static final int LOOKUP_REALLOC_INCREMENT = 10; /** polygons[polyIdx][pointIdx][x/y] */ private double polygons[][][]; /** polyBoundingBoxes [polyIdx][xMin, xMax, yMin, yMax] */ private double polyBoundingBoxes[][]; private double gridX0, gridX1, gridY0, gridY1; private int gridNX, gridNY; private double gridDX, gridDY; /** polyLookup[xIdx*ny+yIdx][0=polyCount, 2*polyInGridIdx+1=POLY_xxx, 2*polyIdx+2=polyIdx] * boxes are inclusive of anything touching any edge * polyLookup[i][j*2+0] == POLY_END_OF_LIST signifies end of poly list */ private int polyLookup[][]; /** Creates grid bounded by max/min of polys and adds those polys * * @param polygons[polyIdx][pointIdx][x/y] */ public PolygonGridLookup(double polygons[][][], int nx, int ny) { gridX0 = Double.POSITIVE_INFINITY; gridX1 = Double.NEGATIVE_INFINITY; gridY0 = Double.POSITIVE_INFINITY; gridY1 = Double.NEGATIVE_INFINITY; for(int i=0; i < polygons.length; i++){ for(int j=0; j < polygons[i].length; j++){ if(polygons[i][j][0] < gridX0) gridX0 = polygons[i][j][0]; if(polygons[i][j][0] > gridX1) gridX1 = polygons[i][j][0]; if(polygons[i][j][1] < gridY0) gridY0 = polygons[i][j][1]; if(polygons[i][j][1] > gridY1) gridY1 = polygons[i][j][1]; } } this.gridDX = (gridX1 - gridX0) / nx; this.gridDY = (gridY1 - gridY0) / ny; this.gridNX = nx; this.gridNY = ny; polyLookup = new int[gridNX*gridNY][]; addPolygons(polygons); } public PolygonGridLookup(double x0, double x1, double y0, double y1, int nx, int ny) { this.gridX0 = x0; this.gridX1 = x1; this.gridY0 = y0; this.gridY1 = y1; this.gridNX = nx; this.gridNY = ny; this.gridDX = (x1 - x0) / nx; this.gridDY = (y1 - y0) / ny; polyLookup = new int[gridNX*gridNY][]; } /** Add the given polygons to the lookup object. /** @param polygons[polyIdx][pointIdx][x/y] */ public void addPolygons(double polygons[][][]){ if(this.polygons != null){ //add polys to poly tracking array first, then to grid double oldPolys[][][] = this.polygons; this.polygons = new double[oldPolys.length + polygons.length][][]; System.arraycopy(oldPolys, 0, this.polygons, 0, oldPolys.length); System.arraycopy(polygons, 0, this.polygons, oldPolys.length, polygons.length); double oldPolyBBoxes[][] = this.polyBoundingBoxes; this.polyBoundingBoxes = new double[oldPolyBBoxes.length + polygons.length][]; System.arraycopy(oldPolyBBoxes, 0, this.polyBoundingBoxes, 0, oldPolyBBoxes.length); calcPolyBoundingBoxes(oldPolyBBoxes.length); addPolysToGrid(oldPolys.length); }else{ this.polygons = polygons; this.polyBoundingBoxes = new double[polygons.length][]; calcPolyBoundingBoxes(0); addPolysToGrid(0); } } /** Calculates bounding boxes of polygons, starting from i0 */ private void calcPolyBoundingBoxes(int i0){ for(int i=i0; i < polyBoundingBoxes.length; i++){ polyBoundingBoxes[i] = new double[] { Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY, // x Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY // y }; for(int j=0; j < polygons[i].length; j++) { if( polygons[i][j][0] < polyBoundingBoxes[i][0]) polyBoundingBoxes[i][0] = polygons[i][j][0]; // minX if( polygons[i][j][0] > polyBoundingBoxes[i][1]) polyBoundingBoxes[i][1] = polygons[i][j][0]; // maxX if( polygons[i][j][1] < polyBoundingBoxes[i][2]) polyBoundingBoxes[i][2] = polygons[i][j][1]; // minY if( polygons[i][j][1] > polyBoundingBoxes[i][3]) polyBoundingBoxes[i][3] = polygons[i][j][1]; // maxY } } } /** Add the polygons in the internal polygons list to appropriate grid cells */ private void addPolysToGrid(int polyIdx0){ for(int polyIdx=polyIdx0; polyIdx < polygons.length; polyIdx++){ //calc grid cell index ranges (inclusive) that this poly might be inside int ix0 = Math.max((int)((polyBoundingBoxes[polyIdx][0] - gridX0) / gridDX), 0); int ix1 = Math.min((int)((polyBoundingBoxes[polyIdx][1] - gridX0) / gridDX), gridNX-1); int iy0 = Math.max((int)((polyBoundingBoxes[polyIdx][2] - gridY0) / gridDY), 0); int iy1 = Math.min((int)((polyBoundingBoxes[polyIdx][3] - gridY0) / gridDY), gridNY-1); //scan those cells for(int ix = ix0; ix <= ix1; ix++){ for(int iy = iy0; iy <= iy1; iy++){ int status = checkPolyCellContact(polyIdx, ix, iy); if(status != POLY_OUTSIDE_CELL) addPolyToGridCell(polyIdx, ix, iy, status); } } } } /** Adds the given polygon to the given cell for the given reason */ private void addPolyToGridCell(int polyIdx, int ix, int iy, int polyCellStatus){ int cellIdx = ix*gridNY + iy; int nPolysInCell; if(polyLookup[cellIdx] == null){ polyLookup[cellIdx] = new int[LOOKUP_REALLOC_INCREMENT * 2 + 1]; nPolysInCell = 0; polyLookup[cellIdx][0] = nPolysInCell; }else{ nPolysInCell = polyLookup[cellIdx][0]; int maxPolys = (polyLookup[cellIdx].length - 1) / 2; if(nPolysInCell >= maxPolys) { //need to reallocate int newPolys[] = new int[(maxPolys + LOOKUP_REALLOC_INCREMENT)*2 + 1]; System.arraycopy(polyLookup[cellIdx], 0, newPolys, 0, nPolysInCell*2 + 1); polyLookup[cellIdx] = newPolys; } } polyLookup[cellIdx][nPolysInCell*2 + 1] = polyCellStatus; polyLookup[cellIdx][nPolysInCell*2 + 2] = polyIdx; polyLookup[cellIdx][0]++; } /** Checks to see if, and if so how, polygon[polyIdx] is of interest to * box (x0 + xi*dx, y0 + yi*dy)-(x0 + xi*dx+dx, y0 + yi*dy+dy) * @param polyIdx * @param xi * @param yi * @return */ private int checkPolyCellContact(int polyIdx, int ix, int iy){ int pointIdx; //calc grid cell boundaries double x0 = gridX0 + ix * gridDX; double x1 = gridX0 + (ix+1) * gridDX; double y0 = gridY0 + iy * gridDY; double y1 = gridY0 + (iy+1) * gridDY; //the easy one first: is poly entirely inside cell (in which case its bbox will be) if(polyBoundingBoxes[polyIdx][0] >= x0 && polyBoundingBoxes[polyIdx][1] <= x1 && polyBoundingBoxes[polyIdx][2] >= y0 && polyBoundingBoxes[polyIdx][3] <= y1) return POLY_ENCLOSED_BY_CELL; //next we need to see if any edge of the cell hits any edge of the poly double x3 = polygons[polyIdx][polygons[polyIdx].length-1][0]; double y3 = polygons[polyIdx][polygons[polyIdx].length-1][1]; for(pointIdx=0; pointIdx < polygons[polyIdx].length; pointIdx++){ double x4 = polygons[polyIdx][pointIdx][0]; double y4 = polygons[polyIdx][pointIdx][1]; if(intersects(x0,y0, x1,y0, x3,y3, x4,y4, null)) return POLY_INTERSECT_CELL;; if(intersects(x1,y0, x1,y1, x3,y3, x4,y4, null)) return POLY_INTERSECT_CELL;; if(intersects(x1,y1, x0,y1, x3,y3, x4,y4, null)) return POLY_INTERSECT_CELL;; if(intersects(x0,y1, x0,y0, x3,y3, x4,y4, null)) return POLY_INTERSECT_CELL;; x3 = x4; y3 = y4; } //now either cell is entirely outside or entirely inside points of cell are //so we can just test a single corner of the cell being inside the poly return isPointInPoly(polyIdx, x0, y0) ? POLY_ENCLOSES_CELL : POLY_OUTSIDE_CELL; } /** Test if the line (x1,y1)-(x2,y2) intersects the line (x3,y3)-(x4,y4) */ private static final boolean intersects(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4, double hitCoord[]){ double ua,ub; ua = ((x4-x3)*(y1-y3) - (y4-y3)*(x1-x3)) / (((y4-y3)*(x2-x1)) - ((x4-x3)*(y2-y1))); ub = ((x2-x1)*(y3-y1) - (y2-y1)*(x3-x1)) / (((y2-y1)*(x4-x3)) - ((x2-x1)*(y4-y3))); if(hitCoord != null){ hitCoord[0] = x1 + ua * (x2 - x1); hitCoord[1] = y1 + ua * (y2 - y1); } return (ua >= 0 && ua <= 1 && ub >= 0 && ub <= 1); } private boolean isPointInPoly(int polyIdx, double x, double y) { double x3 = polygons[polyIdx][polygons[polyIdx].length-1][0]; double y3 = polygons[polyIdx][polygons[polyIdx].length-1][1]; int intersectCount = 0; //do a single line test to off the left of the poly's bbox // (so that we don't get any numerics at the bbox corners) double infX = polyBoundingBoxes[polyIdx][0] - (polyBoundingBoxes[polyIdx][1] - polyBoundingBoxes[polyIdx][0]); double infY = (polyBoundingBoxes[polyIdx][1] + polyBoundingBoxes[polyIdx][1]) / 2; for(int pointIdx=0; pointIdx < polygons[polyIdx].length; pointIdx++){ double x4 = polygons[polyIdx][pointIdx][0]; double y4 = polygons[polyIdx][pointIdx][1]; if(intersects(x, y, infX, infY, x3,y3, x4,y4, null)) intersectCount++; x3 = x4; y3 = y4; } return ((intersectCount % 2) != 0); } /** @return The index into the initial polygons array of the poly covering this point, or -1 if none found */ public int getPolygonAtPointBruteForce(double x, double y) { for(int i=0; i < polygons.length; i++){ //otherwise (POLY_INTERSECT_CELL or POLY_ENCLOSED_BY_CELL), we need todo full test if(isPointInPoly(i, x, y)) return i; } return -1; } /** @return The index into the initial polygons array of the poly covering this point, or -1 if none found */ public int getPolygonAtPoint(double x, double y) { if(x < gridX0 || x > gridX1 || y < gridY0 || y > gridY1) return -1; int ix = (int)((x - gridX0) / gridDX); int iy = (int)((y - gridY0) / gridDY); int cellIdx = ix * gridNY + iy; //nice easy case: no polys intersecting or in cell if(polyLookup[cellIdx] == null || polyLookup[cellIdx][0] <= 0) return -1; int nPolysInCell = polyLookup[cellIdx][0]; for(int i=0; i < nPolysInCell; i++){ int polyCellStatus = polyLookup[cellIdx][i*2 + 1]; int polyIdx = polyLookup[cellIdx][i*2 + 2]; //cell is entirely inside poly, so point is definitely in it if(polyCellStatus == POLY_ENCLOSES_CELL) return polyIdx; //otherwise (POLY_INTERSECT_CELL or POLY_ENCLOSED_BY_CELL), we need todo full test if(isPointInPoly(polyIdx, x, y)) return polyIdx; } return -1; } /** Returns the number of polygons in each grid cell */ public int[][] getPolysInGridCount() { int ret[][] = new int[gridNX][gridNY]; for(int ix=0; ix < gridNX; ix++) for(int iy=0; iy < gridNY; iy++){ int cellIdx = ix*gridNY+iy; ret[ix][iy] = polyLookup[cellIdx] == null ? 0 : polyLookup[cellIdx][0]; } return ret; } }