package algorithmrepository.contouring; import java.util.Collections; import java.util.Comparator; import java.util.Vector; import algorithmrepository.DynamicDoubleArray; /** Provides contouring of data on a regular grid within an arbitrary shaped closed boundary * * Currently this is mostly useful starting inside the boundary * and immeadiately aborting when you hit it. * * The general 'return all contours' isn't that great the moment as * contours which are 'closed' inside the grid boundary but not closed inside * the inner boundary will come back in bits. */ public class BoundedContouring extends Contouring { private double x0, x1, y0, y1; private double dx, dy; private double fFOC = Double.NaN; private int FOCStartGridPoint[] = null; /** The boundary split into elements within the contouring triangular grid */ private double boundaryCache[][]; public OptimisedBoundaryHandler boundary; private BoundaryCheckCallback boundCheckCallback; /* Same sign or 0 if either is 1 - ??? */ protected final boolean NOT_SAME_SIGN(double v1, double v2){ return ( (Double.isNaN(v1) || Double.isNaN(v2)) ? false : ((v1>=0) ? ((v2<0)?true:false) : ((v2>=0)?true:false)) ); } /* Is i odd? */ protected final boolean ISODD(int i){ return ( (i % 2) != 0); } /* *********************************** */ /** Creates object for contouring within a boundary * @param callbacks ContourCollector or other contouting callback * @param x0, x1, nx Grid definition (regular) * @param y0, y1, ny Grid definition (regular) * @param f a Contouring.GridFunction * @param boundX List of X coords of boundary * @param boundY List of Y coords of boundary */ public BoundedContouring(ContourCallback callbacks, double x0, double x1, int nx, double y0, double y1, int ny, GridFunction f, double boundX[], double boundY[]) { super(null, x0, x1, nx, y0, y1, ny, f); initBounded(callbacks, x0, x1, nx, y0, y1, ny, boundX, boundY); } /** Creates object for contouring within a boundary * @param callbacks ContourCollector or other contouting callback * @param x0, x1, nx Grid definition (regular) * @param y0, y1, ny Grid definition (regular) * @param f f[yIdx*nx + xIdx] - Function to contour (flat 2D array) * @param boundX List of X coords of boundary * @param boundY List of Y coords of boundary */ public BoundedContouring(ContourCallback callbacks, double x0, double x1, int nx, double y0, double y1, int ny, double f[], double boundX[], double boundY[]) { super(null, x0, x1, nx, y0, y1, ny, f); initBounded(callbacks, x0, x1, nx, y0, y1, ny, boundX, boundY); } /** Creates object for contouring within a boundary * @param callbacks ContourCollector or other contouting callback * @param x0, x1, nx Grid definition (regular) * @param y0, y1, ny Grid definition (regular) * @param f f[yIdx][xIdx] - Function to contour * @param boundX List of X coords of boundary * @param boundY List of Y coords of boundary */ public BoundedContouring(ContourCallback callbacks, double x0, double x1, int nx, double y0, double y1, int ny, double f[][], double boundX[], double boundY[]) { super(null, x0, x1, nx, y0, y1, ny, f); initBounded(callbacks, x0, x1, nx, y0, y1, ny, boundX, boundY); } public void initBounded(ContourCallback callbacks, double x0, double x1, int nx, double y0, double y1, int ny, double boundX[], double boundY[]){ this.x0 = x0; this.x1 = x1; this.y0 = y0; this.y1 = y1; this.dx = (x1 - x0) / ((double)nx - 1.0); this.dy = (y1 - y0) / ((double)ny - 1.0); boundary = new OptimisedBoundaryHandler(boundX, boundY, x0, x1, nx, y0, y1, ny); boundCheckCallback = new BoundaryCheckCallback(callbacks, boundary); super.setCallbacks(boundCheckCallback); //splitBoundary(); } /** Traces all contours on the whole grid but splits them at the boundary */ public final void getAllGridContours(double f[]) { //flag that we want it to split the contour and continue at the boundary boundCheckCallback.abortOnBoundary = false; for (int i=0 ; i < f.length; i++) { /* for each contour */ clearFlags(); getGridBoundaryContours(f[i]); getGridRemainingContours(f[i]); } } public final void setCallbacks(ContourCallback callbacks){ /* Don't set ours, set the inner ones */ boundCheckCallback.innerCallback = callbacks; } /** Gets contours starting from the grid boundary */ public void getGridBoundaryContours(double Cont) { //this is just the normal one super.getBoundaryContours(Cont); } /** Gets the rest of all the contours in the grid after getGridBoundaryContours() */ public void getGridRemainingContours(double Cont) { //this is just the normal one super.getRemainingContours(Cont); } @Override /** Gets contours starting from the boundary */ public void getBoundaryContours(double Cont) { if(true) throw new RuntimeException("Not yet implemented."); //walk the split boundary (last point == first) for(int k=0; k < (boundaryCache.length-1); k++){ //element ends double xA = boundaryCache[k][0], xB = boundaryCache[k+1][0]; double yA = boundaryCache[k][1], yB = boundaryCache[k+1][1]; int iA = (int)boundaryCache[k][3], jA = (int)boundaryCache[k][4]; int iB = (int)boundaryCache[k+1][3], jB = (int)boundaryCache[k+1][4]; //get function val at each end of the element using the cached interpolation coefficients double fA = boundaryCache[k][5] * FVAL(jA, iA) + boundaryCache[k][6] * FVAL((jA+1), iA) + boundaryCache[k][7] * FVAL(jA, (iA+1)) + boundaryCache[k][8] * FVAL((jA+1), (iA+1)) ; double fB = boundaryCache[k+1][5] * FVAL(jB, iB) + boundaryCache[k+1][6] * FVAL((jB+1), iB) + boundaryCache[k+1][7] * FVAL(jB, (iB+1)) + boundaryCache[k+1][8] * FVAL((jB+1), (iB+1)) ; if (NOT_SAME_SIGN ( fA - Cont, fB - Cont) ){ //so we know the contour hits this segment. double l = (Cont - fA) / (fB - fA); //normalised dist along element double x = xA + l * (xB - xA); //position of contact with contour double y = yA + l * (yB - yA); int i = (int)((x - x0)/dx); int j = (int)((y - y0)/dy); double fx = x - (x0 + j * dx); //fractional position of contact with contour double fy = y - (y0 + i * dy); //We need to fire the callback, (ours as well) with that point so it starts ON the boundary boundCheckCallback.contourPoint(Cont, x, y); //value at centre of square double fc = ( FVAL(j, i) + FVAL(j+1, i) + FVAL(j, i+1) + FVAL(j+1, i+1)) / 4.0; /* int triang = whichTriangle(fx, fy); switch(triang){ case 0: //Top if (NOT_SAME_SIGN ( f[i*nx+j] - Cont, f[i*nx+(j+1)] - Cont) ){ //Hits top //We can just fire this } break; case 1: //Left break; case 2: //Right break; case 3: //Bottom break; } */ } } } public double splitBX[]; public double splitBY[]; /** Take the boundary and sllit each segment so no element covers more than 1 triangle * of the contouring grid. */ public void splitBoundary(){ throw new RuntimeException("Needs replacing with square grid version"); /*Vector allPoints = new Vector(); / Trace boundary (open) contours / double bx[] = boundary.getBoundX(); double by[] = boundary.getBoundY(); for(int k=0; k < bx.length-1; k++){ //for each boundary element //get the grid rectangle enclosing this boundary segment int box[] = boundary.getBoxEnclosing(bx[k],by[k],bx[k+1],by[k+1]); int i0 = box[0], i1 = box[1]; int j0 = box[2], j1 = box[3]; Vector points = new Vector(); /* for each grid cell in the box, test for intersections with these edges: * *---0---* * |\ / * | 4 1 * | \ / * 3 X * | / \ * | 2 5 * |/ \ * * * / double hitCoord[] = new double[9]; //x, y and l = dist along boundary element, i, j, a, b, c, d //a,b,c,d are the coefficient of the function value four corners to get the //function value at this coordinate, given the triangular grid //add the start point (not the end one, it'll be the start of the next boundary element //work out which triangle the point is in int i = (int)( (by[k] - y0)/dy ); //grid square that point is in int j = (int)( (bx[k] - x0)/dx ); double fdy = (by[k] - (y0 + i * dy)) / dy; //fractional position in square double fdx = (bx[k] - (x0 + j * dx)) / dx; //fractional position in square double cf[] = getInterpolationCoefficients(fdx, fdy); points.add(new double[]{ bx[k], by[k], 0, i, j, cf[0], cf[1], cf[2], cf[3] }); for(i=0; i <= ny; i++){ double gridY0 = y0 + i * dy; double gridY1 = gridY0 + dy; double gridYC = gridY0 + dy/2; for( j=0; j <= nx; j++){ double gridX0 = x0 + j * dx; double gridX1 = gridX0 + dx; double gridXC = gridX0 + dx/2; //coordinates in order as shown on diagram double sx[] = new double[]{ gridX0, gridX1, gridXC, gridX0, gridX0, gridXC, gridX1 }; double sy[] = new double[]{ gridY0, gridY0, gridYC, gridY1, gridY0, gridYC, gridY1 }; //test intersections with each segmenet 0-5 for(int c=0; c<6; c++){ if( OptimisedBoundaryHandler.intersects(sx[c], sy[c], sx[c+1], sy[c+1], bx[k], by[k], bx[k+1], by[k+1], hitCoord) ){ double lx = (hitCoord[0] - bx[k]); double ly = (hitCoord[1] - by[k]); hitCoord[2] = lx*lx + ly*ly; hitCoord[3] = i; hitCoord[4] = j; //also cache the coefficients of the grid square corners to //later work out the value of f at this point cf = getInterpolationCoefficients( (hitCoord[0] - gridX0)/dx, (hitCoord[1] - gridY0)/dy ); System.arraycopy(cf, 0, hitCoord, 5, 4); points.add(hitCoord.clone()); //add to points array } } } } //now sort the points class PtCompare implements Comparator{ public int compare(double[] o1, double[] o2) { return (int)Math.signum(o1[2] - o2[2]); } }; Collections.sort(points, new PtCompare() ); //and add them in that order to the main points array allPoints.addAll(points); } //add the very last point to complete the loop (it will be the same as the first, so we can use that info) allPoints.add(allPoints.get(0).clone()); boundaryCache = new double[allPoints.size()][]; for(int i=0;i gridMax){ gridMax = FVAL(j, i); maxI = i; maxJ = j; } if(FVAL(j, i) < gridMin){ gridMin = FVAL(j, i); minI = i; minJ = j; } } double slX0, slY0; if(fromGridMax){ slX0 = x[maxJ]; slY0 = y[maxI]; }else{ slX0 = x[minJ]; slY0 = y[minI]; } return getLastClosedContour(slX0, slY0, nBisects, targetDf, passDownAttempts); } /** Starting from the specified point (slX0, slY0), * searches the line from there in the +x direction to the boundary * for the highest f value (assuming it monotonically increases, * that gives a closed contour within the boundary. Once a bracket * LCC < f < FOC is found, the contour LCC is re-run with the * exterior callbacks enabled. * * @param targetDf Required seperation (f_FOC - f_LCC) * @param passDownAttempts If true all attempts at finding the LCC are done with the * exterior callbacks enabled. * * @return bracket double[]{ f of last known closed contour, f of first known open contour } */ public double[] getLastClosedContour(double slX0, double slY0, int nBisects, double targetDf, boolean passDownAttempts){ double fLastClosedContour = Double.NaN; double fFirstOpenContour = Double.NaN; // Last closed and first open contours /* define the line from min/max point to boundary */ double slX1, slY1; double hitCoord[] = new double[2]; if( !boundary.isCrossing(slX0, slY0, x[nx-1], slY0, hitCoord) ) throw new RuntimeException("No boundary found between start point and right grid edge!!"); slX1 = hitCoord[0]; slY1 = hitCoord[1]; // System.out.println("Search line: ("+slX0+","+slY0+")-("+slX1+","+slY1+")"); if( ((slY1 - slY0) * (slY1 - slY0)) > 1e-24){ throw new RuntimeException("sanity check failed: Search line should be at const Y. Y0 = " + slY0 + ", Y1 = " + slY1); } /* Setup and enable search line intersection counting in countouring callbacks */ boundCheckCallback.searchLine = new double[][]{ { slX0, slX1 }, { slY0, slY1 } }; /* Disable anything further down the line (and store it)*/ ContourCallback storedCallback = boundCheckCallback.innerCallback; if(!passDownAttempts) boundCheckCallback.innerCallback = null; boolean storedAbort = boundCheckCallback.abortOnBoundary; boundCheckCallback.abortOnBoundary = true; double curY = slY0; int i = (int)((curY - y[0]) / dy); double l0 = 0, l1 = 1; //normalised positions along search line of current search region int jLCC = 0; for(int k=0; k < nBisects; k++){ /* clc new position to search from */ double l = l0 + (l1 - l0) / 2.0; double curX = slX0 + l * (slX1 - slX0); int j = (int)((curX - x[0]) / dx); /* calc grid cell this is in */ double jp = (curX - x[0] - (j * dx)) / dx; /* calc normalised position within the cell */ double curF = (1.0-jp)*FVAL(j, i) + jp*FVAL(j+1, i); /* and linear interpolate to get psi of the point */ boundCheckCallback.sliCount = 0; /* reset the line intersrection count */ clearFlags(); //int ret = contourFromSquare(j+1, i, j, i, true, curF); int ret = contourFromEdge(j, i-1, EDGE_BOTTOM, curF); /* contourPoint() will count the SLI point which we start from twice (at beginning and at end), so drop one * sometimes this doesn't work and it only hits the starting SLI point once * This should be ok in the enclosing case but will break the even/odd test if the surface * does intersect us multiple times. * FIXME: Think about this a bit * */ if(boundCheckCallback.sliCount >= 2) boundCheckCallback.sliCount--; //else // System.err.println("WARNING: Contouring find LCC only hit the search line "+boundCheckCallback.sliCount+" times. This is definately odd."); //System.out.println("k="+k+", j="+j+", x="+curX+", curF="+curF+", sliCount="+boundCheckCallback.sliCount+", ret="+ret); /* if the contour is closed AND passes the search line an odd number of times it must * enclose the magnetic axis, so we want to look further out */ if( ret == RETURN_FAILED){ //mumble mumble System.err.println("Contouring failed in LCFS Finder. Nugdging"); l1 = l1 + 0.01 * (l1 - l0); continue; } if( ret == RETURN_CLOSED ){ if( ISODD(boundCheckCallback.sliCount) ){ l0 = l; fLastClosedContour = curF; jLCC = j; } else { System.err.println("Search line: ("+slX0+","+slY0+")-("+slX1+","+slY1+")"); throw new ArithmeticException("Fatal: Found a closed surface (f="+curF+") that doesn't enclose the "+ "search line ends - it hit the line " + boundCheckCallback.sliCount + " times. Don't know where to go next!!"); } } else { /* otherwise we hit the boundary (or the grid boundary)- so it's open and we need to look further in */ l1 = l; fFirstOpenContour = curF; this.FOCStartGridPoint = new int[]{ i, j }; } /* Check if we're reached target f_LCC - f_FCC */ if(!Double.isNaN(fLastClosedContour) && !Double.isNaN(fFirstOpenContour)){ double df = (fLastClosedContour - fFirstOpenContour); if( (df*df) < (targetDf*targetDf)) break; } } /* Restore the inner callback (if any) */ boundCheckCallback.innerCallback = storedCallback; /* If there is one, run the LCC through it */ if(!Double.isNaN(fLastClosedContour) && storedCallback != null){ boundCheckCallback.sliCount = 0; /* reset the line intersrection count */ clearFlags(); int ret = contourFromEdge(jLCC, i, EDGE_TOP, fLastClosedContour); } this.fFOC = fFirstOpenContour; // Save for calls to getFirstOpenContour /* disable search line intersection counting and restore the abort status*/ boundCheckCallback.abortOnBoundary = storedAbort; boundCheckCallback.searchLine = null; //System.out.println("LCC Finder terminated: LCC = " + fLastClosedContour + ", FOC = "+fFirstOpenContour); return new double[]{ fLastClosedContour, fFirstOpenContour }; } /** Traces the first open contour from the point it was found on the * FCC search line in getLastClosedContour. * * This will actually (probably) return 2 contours, each starting * from the point on the search line and ending on the boundary. */ public void traceFirstOpenContourComplete(){ if(Double.isNaN(fFOC)) throw new RuntimeException("No first open contour found, have you run getLastClosedContour()??"); boundCheckCallback.abortOnBoundary = true; int i = FOCStartGridPoint[0]; int j = FOCStartGridPoint[1]; //Start from the known start point in the original direction clearFlags(); int ret = contourFromEdge(j, i, EDGE_TOP, fFOC); //and do it again in the other direction clearFlags(); ret = contourFromEdge(j, i-1, EDGE_BOTTOM, fFOC); } }