package algorithmrepository.contouring;

import oneLinersForAlgoTests.algoAsciiMatrixFile;
import oneLinersForAlgoTests.algoOneLiners;

/** Given data on a regular grid, this provides various functions
 * to follow contours of given levels.
 * 
 * Original algorithm from Scilab.
 * 
 * @author oliford <codes(at)oliford.co.uk>
 *
 */

public class ContouringTriangularScilab {
    
    int dbgI0, dbgJ0, dbgIB ,dbgJB;
    boolean dbgMBC;
    double dbgCVAL;
    
    /** Temporary for debugging */
    public void epicFail(String trap){
        
        String outPath = System.getProperty("java.io.tmpdir") + "/contourFail";
        System.err.println("**** CONTOURING FAILURE *****");
        System.err.println("Please send the directory '"+outPath+"' to oliford");
        System.err.println("If you want to continue anyway, comment out ONLY THE CALL TO epicFail() WHICH FIRED THIS TRAP.");
       
        algoAsciiMatrixFile.mustWrite(outPath+"/f.txt",new double[][]{f},false);
        algoAsciiMatrixFile.mustWrite(outPath+"/x.txt",new double[][]{x},false);
        algoAsciiMatrixFile.mustWrite(outPath+"/y.txt",new double[][]{y},false);
        algoAsciiMatrixFile.mustWrite(outPath+"/done.txt", algoOneLiners.intToDouble(new int[][]{done}), false);
        
        algoOneLiners.TextToFile(outPath+"/info.txt",
            "trap="+trap+"\n"+
            "nx="+nx+"\n"+
            "ny="+ny+"\n"+
            "i0="+dbgI0+"\n"+
            "j0="+dbgJ0+"\n"+
            "ib="+dbgIB+"\n"+
            "jb="+dbgJB+"\n"+
            "MBC="+dbgMBC+"\n"+
            "CVAL="+dbgCVAL+"\n" );
        
        //throw new RuntimeException("Contouring failure: '"+trap+"'");
    }
	
	/** The data grid - setup by constructor but f can be updated/modified later */
	protected int 		nx,ny; 	/* grid size (number of SQUARES = number of VERTICIES - 1)*/
	private int 		done[]; 	/* done flag (ny x nx)*/
	protected double 		f[]; 	/* function values at verticies (ny x nx)*/
	protected double		x[];	/* x axis value (nx) */
	protected double		y[];	/* y axis values (ny) */
	
	private ContourCallback callbacks;
	
	/* Segment IDs and return codes from findExitSegment() (internal only) */
	private static final int SEGMENT_BELOW		= 1;
	private static final int SEGMENT_LEFT		= 2;
	private static final int SEGMENT_ABOVE		= 3;
	private static final int SEGMENT_RIGHT		= 4;
	
	
	/* Return codes from contourFromSquare() (external and internal) */
	public static final int RETURN_PATHOLOGICAL_CASE = -3; //really silly cases where the contour hits a grid node and we get lost
	public static final int RETURN_HIT_GRID_BOUNDARY = -2; 
	public static final int RETURN_ABORTED		= -1;
	public static final int RETURN_CLOSED		= 0;
	
	/* done flags (bitfield)*/
	private final int HORZ_DONE	= (1 << 0);
	private final int VERT_DONE	= (1 << 1);

	
	/* Functions from macros in c */
	/* macros for 2d matrix accesses (assumes struct contData *g is present local)*/
	private final int DONE_FLAG(int x, int y){
		return done[(y)*nx + (x)];
	}
	
	private final boolean IS_VERT_DONE(int x, int y){	return ((done[(y)*nx + (x)] & VERT_DONE) != 0); }	
	private final boolean IS_HORZ_DONE(int x, int y){	return ((done[(y)*nx + (x)] & HORZ_DONE) != 0); }
	private final void SET_VERT_DONE(int x, int y){ done[(y)*nx + (x)] |= VERT_DONE; }
	private final void SET_HORZ_DONE(int x, int y){ done[(y)*nx + (x)] |= HORZ_DONE; }
		
	private final double FVAL(int x, int y){
		return f[(y)*nx + (x)];
	}
	
	/* the coordinates between [xi,xj] along one axis for which the value of f is zCont */ 
	private final double F_INTERCEPT(double zCont, double fi, double xi, double fj, double xj){
			return ((xi) + ((zCont)-(fi)) * ((xj)-(xi)) / ((fj)-(fi)));
	}

	/* check for boundary points */
	private final boolean	IS_ON_BOUNDARY(int i, int j) {
		return ( (j)==0 || (i)==0 || (j)==(ny-1) || (i)==(nx-1) );
	}

	/* 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); }
	/* *********************************** */
	
	
	
	/**
	 *  Traces a contour starting from place it is known to be. The contour
	 *  will be traced in an abitary direction until it is aborted by the callback,
	 *  hits the grid boundary or closes on itself.
	 *  
	 *  NOTE: If it hits a boundary the rest of the contour that goes in the other 
	 *  direction from the start position will not be included.
	 *
	 *  The contour is known to be crossing the segment between 
	 *  the points (ib,jb) and (i0,j0).
	 *  
	 *  Cont: value of f along the contour
	 *   
	 *  mightBeClosed: Is it possible that this might close on itself?
	 *	(i.e did we start from inside the domain)
	 */
	int contourFromSquare(int i0, int j0, int ib, int jb, boolean mightBeClosed, double contourValue)
	{
	    dbgI0 = i0; dbgJ0 = j0; dbgIB = ib; dbgJB = jb; dbgMBC = mightBeClosed; dbgCVAL = contourValue;
	    
		int i=i0, j=j0;
		int ip, jp, im, jm, segment=0;
		boolean wflag, onBoundary[] = new boolean[1];
		int seg0 = RETURN_ABORTED;
		jp = j+1; ip = i+1;
		jm = j-1; im = i-1;

		/* we look at how the segment is starting */
		if( jb == jm)
			seg0 = SEGMENT_BELOW; /* is vertical top(jb)->bottom(j) */		    
		else  { 
			if( ib == im )
				seg0 = SEGMENT_LEFT ; /* is horizontal left(ib)->right(i) */
			else {
				if( jb == jp )
					seg0 = SEGMENT_ABOVE ; /* is vertical bottom(jb)->top(j) */
				else  if ( ib == ip )
					seg0 = SEGMENT_RIGHT; /* is horizontal right(ib)->left(i) */
			}
		}

		/* look at direction of the test segment (which the contour crosses) 
		 * and mark the initial point */
		switch (seg0) {
			case SEGMENT_BELOW: /* test line is vertical (i,j-1)-->(i,j) */
				if ( IS_VERT_DONE(i,jm))
					return RETURN_ABORTED;

				segment=SEGMENT_BELOW ; /* so the segment is below (i,j) */
				if(!callbacks.contourPoint(/*0,*/ contourValue,
						x[i], 
						F_INTERCEPT(contourValue, FVAL(i,jm), y[jm], FVAL(i,j), y[j])  ))
					return RETURN_ABORTED;
				break;
			case SEGMENT_LEFT:  /* test line is horizontal (i-1,j)-->(i,j) */
				if ( IS_HORZ_DONE(im,j) )
					return RETURN_ABORTED;

				segment=SEGMENT_LEFT ; /* so the segment is horizontal left */
				if(!callbacks.contourPoint( /*0,*/ contourValue,
						F_INTERCEPT(contourValue, FVAL(im,j), x[im], FVAL(i,j), x[i]),
						y[j]  ))
					return RETURN_ABORTED;
				break ; 
			case SEGMENT_ABOVE: /* test line is vertical (i,j+1)-->(i,j) */
				if ( IS_VERT_DONE(i,j) )
					return RETURN_ABORTED;

				segment=SEGMENT_ABOVE ; /* the segment is vertical up */
				if(!callbacks.contourPoint( /*0,*/contourValue,
						x[i],
						F_INTERCEPT(contourValue, FVAL(i,j), y[j], FVAL(i,jp), y[jp])  ))
					return RETURN_ABORTED;
				break ;
			case SEGMENT_RIGHT: /* test line is horizontal (i+1,j)-->(i,j) */
				if ( IS_HORZ_DONE(i,j))
					return RETURN_ABORTED;

				segment=SEGMENT_RIGHT ; /* segment is the horizontal right */
				if(!callbacks.contourPoint( /*0,*/ contourValue,
						F_INTERCEPT(contourValue, FVAL(i,j), x[i], FVAL(ip,j), x[ip]),
						y[j]  ))
					return RETURN_ABORTED;
				break;
			default:
				throw new RuntimeException("Unrecognised initial segment type '"+seg0+"'.");
		}

		wflag=true;
		do{
			jp=j+1; ip= i+1; jm=j-1;im=i-1;

			switch (segment) {
				case SEGMENT_BELOW: /* test line is vertical (i,j-1)-->(i,j) */					
					SET_VERT_DONE(i,jm);
					segment = findExitSegment(i, ip, ip, i, j, j, jm, jm, segment, contourValue, onBoundary);					
					/* calculate the new item, segment gives direction to explore segment */
					switch (segment) {
						case SEGMENT_BELOW: i=ip ; break;
						case SEGMENT_LEFT: i=ip; j=jm; break;						 
					}
					break;
				case SEGMENT_LEFT: /* test line is horizontal (i-1,j)-->(i,j) */
					SET_HORZ_DONE(im,j);
					segment = findExitSegment(i, i, im, im, j, jm, jm, j, segment, contourValue, onBoundary);
					switch (segment){ 
						case SEGMENT_LEFT: j = jm; break;
						case SEGMENT_ABOVE : i=im;j=jm; break;
					}
					break;
				case SEGMENT_ABOVE: /* test line is vertical (i,j+1)-->(i,j) */
					SET_VERT_DONE(i,j);
					segment = findExitSegment(i, im, im, i, j, j, jp, jp, segment, contourValue, onBoundary);
					switch (segment){ 
						case SEGMENT_ABOVE : i=im; break;
						case SEGMENT_RIGHT: i=im; j=jp; break;
					}
					break;
				case SEGMENT_RIGHT: /* test line is horizontal (i+1,j)-->(i,j) */
					SET_HORZ_DONE(i,j);
					segment = findExitSegment(i, i, ip, ip, j, jp, jp, j, segment, contourValue, onBoundary);
					switch (segment){
						case SEGMENT_RIGHT :j=jp; break;
						case SEGMENT_BELOW :i=ip; j=jp; break;
					}
					break;
				default:
					throw new RuntimeException("Unrecognised segment type '"+segment+"'.");
			}

			/** new segment is on the boundary **/
			if (onBoundary[0]){
				switch (segment) {
					case SEGMENT_BELOW: SET_VERT_DONE(i,(j-1)); break ; 
					case SEGMENT_LEFT: SET_HORZ_DONE(i-1,j);  break ; 
					case SEGMENT_ABOVE: SET_VERT_DONE(i,j); break ; 
					case SEGMENT_RIGHT: SET_HORZ_DONE(i,j); break ; 
				}
				/** we must quit the while loop **/
				wflag = false;
				
			}

			/** If we might be getting a closed surface
			 * we need to test the done flags to make sure we don't
			 * go back over ourselves
			 */
			if ( mightBeClosed && i==i0 && j==j0 && segment == seg0)
				segment = RETURN_CLOSED;
			
		}while(segment > 0 && wflag);
		
		if( (segment != RETURN_ABORTED) && (segment != RETURN_PATHOLOGICAL_CASE) && onBoundary[0])
			segment = RETURN_HIT_GRID_BOUNDARY;

		if(!callbacks.contourEnd(contourValue, segment))
			return RETURN_ABORTED;
		
		return segment;
	}

	/** 
	 *   findExitSegment : This function is called with four points defined clockwise
	 * The contour is known to pass between (i1,j1) and (i4,j4)
	 * Mark the points it passes through the triangles of this cell
	 * and return which segment is leaves through
	 *  
	 *//* -end javadoc-
	 *
	 * (i1,j1) *------------* (i2,j2)
	 *         |\          /|
	 *         | \   __   / |
	 *         |  \ /  \ /  |
	 *         |   /    \   |
	 *         |  / \  / \  |
	 *         | /   \/   \_______-->
	 * _________/    /\     |
	 *         |    /  \    |
	 *         |   /    \   |
	 *         |  /      \  |
	 *         | /        \ |
	 *         |/          \|
	 * (i4,j4) *------------* (i3,j3)
	 *
	 *-----------------------------------------------------------------------*/
	private final int  findExitSegment(int i1, int i2, int i3, int i4,
						int j1, int j2, int j3, int j4,
						int segment, double Cont, boolean onBoundary[]){
		
		double phi1, phi2, phi3, phi4, xav, yav, phiav;
		int squareFlipped,i;
		
		/* difference between f at each vertex, and target f */
		phi1 = FVAL(i1, j1) - Cont;
		phi2 = FVAL(i2, j2) - Cont;
		phi3 = FVAL(i3, j3) - Cont;
		phi4 = FVAL(i4, j4) - Cont;
		squareFlipped = 0;
		onBoundary[0] = false;
		
		//cases where the corners are 0 may often brake the algorithm
		//Simplest way to deal with these is to just give up and go
		//and do a very slighltly higher or lower contour
		if(phi1 == 0 || phi2 == 0 || phi3 == 0 || phi4 == 0)
		    return RETURN_PATHOLOGICAL_CASE;

		/* the point at the center of the rectangle (x/y average)*/
		xav = (x[i1] + x[i3]) / 2.0 ; 
		yav = (y[j1] + y[j3]) / 2.0 ; 
		phiav = (phi1 + phi2 + phi3 + phi4) / 4.0;	/* average diff at all 4 verticies */
		
		//check it actually passes through the entry edge
		if( !NOT_SAME_SIGN(phi1, phi4) ){
		    System.err.println("Not entering left edge trapped on square:\n"+
		        "f("+i1+","+j1+") = "+phi1+"\n" +
		        "f("+i2+","+j2+") = "+phi2+"\n" +
		        "f("+i3+","+j3+") = "+phi3+"\n" +
		        "f("+i4+","+j4+") = "+phi4+"\n" +
                "phiav = " + phiav);
	            
		    if(phi1 == 0 || phi4 == 0)
		        System.err.println("Contouring detected pathological case phi1 = "+phi1+", phi4 = "+phi4+". Continuing with hope...");
		    else
		        epicFail("Not entering first edge");
		    return RETURN_PATHOLOGICAL_CASE;
		}

		if ( Double.isNaN(phiav) )
			return RETURN_ABORTED;

		/* does it pass between centre and 4th? (i.e anticlockwise) */
		if ( NOT_SAME_SIGN(phiav, phi4) ) 
		{
			int iTmp, jTmp; 
			double phiTmp;

			squareFlipped = 1 ; 
			/* flip points vertically to make it clockwise */
			iTmp = i4; i4 = i1; i1 = iTmp; /* swap points 1 & 4 */
			jTmp = j4; j4 = j1; j1 = jTmp;

			iTmp = i3; i3 = i2; i2 = iTmp; /* swap points 2 & 3*/
			jTmp = j3; j3 = j2; j2 = jTmp;

			phiTmp = phi1; phi1 = phi4; phi4= phiTmp;
			phiTmp = phi3; phi3 = phi2; phi2= phiTmp;
		}

		/* we store a new item */
		if(!callbacks.contourPoint(/*1,*/ Cont,
				F_INTERCEPT(0.0, phi1, x[i1], phiav, xav),
				F_INTERCEPT(0.0, phi1, y[j1], phiav, yav) ))
			return RETURN_ABORTED;

		/*
		* Run through each triangle in clockwise dir until we find which edge the contour leaves through
		*
		*/
		for (i=0;  ; i++){ /* i= 0,1,2 */
			int iTmp, jTmp;
			double phiTmp;

			/** if it leaves through the square through this edge, exit the loop **/
			if ( NOT_SAME_SIGN(phi1, phi2) )
				break;

			if  ( phiav != 0.0 ) { /*ignoring the case where it goes exactly through the sqaure's centre */
				/* report the point where is crosses into the next triangle */
				if(!callbacks.contourPoint(/*1,*/ Cont,
					F_INTERCEPT(0.0, phi2, x[i2], phiav, xav),
					F_INTERCEPT(0.0, phi2, y[j2], phiav, yav)))
				return RETURN_ABORTED;
			} 
			
			if (i >= 3) { //we ought to find an exit route, otherwise we have lost the plot
			    System.err.println("I>=3 trapped on square:\n"+
	                "f("+i1+","+j1+") = "+phi1+"\n" +
	                "f("+i2+","+j2+") = "+phi2+"\n" +
	                "f("+i3+","+j3+") = "+phi3+"\n" +
	                "f("+i4+","+j4+") = "+phi4+"\n" +
	                "phiav = " + phiav);
	            
			    epicFail("I>=3 trap");
			    return RETURN_PATHOLOGICAL_CASE;
			}

			/** rotate points anticlockwise so we look next at the next section in the clockwise dir **/
			iTmp = i1; i1 = i2; i2 = i3; i3 = i4; i4 = iTmp;
			jTmp = j1; j1 = j2; j2 = j3; j3 = j4; j4 = jTmp;
			phiTmp = phi1; phi1 = phi2; phi2 = phi3; phi3 = phi4; phi4 = phiTmp;
		}

		/* report the point where it leaves the square */
		if(!callbacks.contourPoint(/*1,*/ Cont,
				F_INTERCEPT(0.0, phi1, x[i1], phi2, x[i2]),
				F_INTERCEPT(0.0, phi1, y[j1], phi2, y[j2]) ))
			return RETURN_ABORTED;

		/* if the new segment is on an edge, report the end */
		if ( IS_ON_BOUNDARY(i1, j1) && IS_ON_BOUNDARY(i2, j2))
			onBoundary[0] = true;

		/* i now tells us the edge number it left through relative to the entry edge.
		 * segment tells us which segment it entered through
		 * so we return the segment number its leaving through
		 */

		/* but deal with the case where we flipped the grid */
		if ( squareFlipped == 1 && !ISODD(i) )
			i = i+2;

		return ( 1 + ( ( i + segment + 2) % 4));
	}
	
	/** Traces a contour starting from the position (x0,y0) */
	public int contourFromPoint(double x0, double y0, double Cont){
	    int i0=0, j0=0;
	    for(int i=0; i<(ny-1); i++)
	        if(y[i] < y0 && y[i+1] > y0){
	            i0 = i;
	            break;
	        }
	    
	    for(int j=0; j<(nx-1); j++)
            if(y[j] < x0 && y[j+1] > x0){
                j0 = j;
                break;
            }
	    
	    return contourFromPoint(i0, j0, x0, y0, Cont);
	}
	
	/** Traces a contour from the point (x0,y0) which must be in square (j0,i0)-(j0+1,i0+1) */
	public int contourFromPoint(int i0, int j0, double x0, double y0, double Cont){
	    
	    callbacks.contourPoint(/*0,*/Cont, x0, y0);
	    
	    double fx = (x0 - x[j0]) / (x[j0+1] - x[j0]);
        double fy = (y0 - y[i0]) / (y[i0+1] - y[i0]);
        
        if(true)throw new RuntimeException("Not yet implemented");
        	    
	    switch(whichTriangle(fx, fy)){
        case 0: //Top
            break;
        case 1: //Left
            break;
        case 2: //Right
            break;
        case 3: //Bottom
            break;
	    }
	     
	    
	    return 0;
	}
	    
	
	/** Clears cell done flags ready for next contour */
	protected final void clearFlags() {
		done = new int[ny*nx];
	}
	
	private static final double[] linearSpace(double x0, double x1, int nx){
	    double x[] = new double[nx];
	    double dx = (x1 - x0) / ((double)nx - 1.0);
	    for(int i=0; i<nx; i++)
	        x[i] = x0 + i * dx;
	    return x;
	}
	
	/** Creates contouring object based on regular grid of data.     
     * 
     * @param f values at grid points (flat - x is deepest dimension)
     */
    public ContouringTriangularScilab(ContourCallback callbacks,
                        double x0, double x1, int nx,
                        double y0, double y1, int ny,
                        double f[]) {
        this.callbacks = callbacks;
        x = linearSpace(x0, x1, nx);
        y = linearSpace(y0, y1, ny);        
        this.nx = nx; this.ny = ny;
        this.f = f;
        if(f != null && f.length != (nx*ny))
            throw new IllegalArgumentException("f[] should to be (x.length * y.length) long.");
        
        clearFlags();
    }
    
    /** Creates contouring object based on regular grid of data.     
     * 
     * @param f values at grid points ( is deepest dimension)
     */
    public ContouringTriangularScilab(ContourCallback callbacks,
                        double x0, double x1, int nx,
                        double y0, double y1, int ny,
                        double f[][]) {
        this.callbacks = callbacks;
        x = linearSpace(x0, x1, nx);
        y = linearSpace(y0, y1, ny);        
        this.nx = nx; this.ny = ny;
        
        setF(f);
        
        clearFlags();
    }

	/** Creates contouring object based on grid of data.	 
	 * 
	 * @param x x coordinates of grid points.
	 * @param y y coordinates of grid points.
	 * @param f values at grid points (flat - x is deepest dimension)
	 */
	public ContouringTriangularScilab(ContourCallback callbacks, double x[], double y[], double f[]) {
	    this.callbacks = callbacks;
		this.x = x;
		this.y = y;
		this.nx = x.length;
		this.ny = y.length;
		this.f = f;
		if(f != null && f.length != (nx*ny))
			throw new IllegalArgumentException("f[] should to be (x.length * y.length) long.");
		
		clearFlags();
	}
	
	/** Creates contouring object based on grid of data.	 
	 * 
	 * @param x x coordinates of grid points.
	 * @param y y coordinates of grid points.
	 * @param f values at grid points (x is deepest dimension)
	 */
	public ContouringTriangularScilab(ContourCallback callbacks, double x[], double y[], double f[][]) {
	    this.callbacks = callbacks;
		this.x = x;
		this.y = y;
		this.nx = x.length;
		this.ny = y.length;
		
		setF(f);
		
		clearFlags();
	}
	
	public final void setF(double f[]){
	    this.f = f;
	}
	
	public final void setF(double f[][]){
	    if(f.length != ny || f[0].length != nx)
            throw new IllegalArgumentException("f[][] should to be (y.length) x (x.length)");
        
        this.f = new double[ny*nx];
        for(int i=0;i<ny;i++)
            System.arraycopy(f[i], 0, this.f, i*nx, nx);
	}
	
	public void setCallbacks(ContourCallback callbacks){ this.callbacks = callbacks; }
	
	
	/** Traces all contours */
	public final void getAllContours(double f[]) {
		

		for (int i=0 ; i < f.length; i++) { /* for each contour */
			/** done is a flag array to memorise checked parts of the grid **/
			clearFlags();
			getBoundaryContours(f[i]);
			getRemainingContours(f[i]);
		}
	}
	
	/** Trace boundary (open) contours */
	public void getBoundaryContours(double Cont) {
			int i,j,k,ib,jb;
			
			int boundary[][] = calcBoundary();

			/** all the boundary segments **/
			for (k = 1; k < boundary[0].length; k++) {
				i = boundary[0][k];
				j = boundary[1][k];
				ib = boundary[0][k-1];
				jb = boundary[1][k-1];
				if  (NOT_SAME_SIGN (	FVAL(i, j) - Cont,
							FVAL(ib, jb) - Cont) )
					contourFromSquare(i, j, ib, jb, false, Cont);
			}
	}
		
	/** Gets the rest of the contours (closed contours) after getBoundaryContours() */
	public void getRemainingContours(double Cont) {
			int i,j;

			/** inside segments **/
			for ( i = 1 ; i < nx-1; i++)
				for ( j = 1 ; j < ny-1 ; j++)
					if ( NOT_SAME_SIGN(	FVAL(i, j) - Cont,
								FVAL(i, j-1) - Cont) )
						contourFromSquare( i, j, i, j-1, true , Cont);
	}


	/** Fills in the boundary definition part of contData */
	private final int[][] calcBoundary() {
		int j,k,ci;

		int nBound = 2*(nx)+2*(ny)-3; /* Boundary size?? */

		int xBound[] = new int[nBound];
		int yBound[] = new int[nBound];
		
		/* construct the parametrisation of the boundary points */
		/*
		 *	*-->*-->*-->*
		 *	^           |
		 *	|           v
		 *	*   *   *   *
		 *	^          / 
		 *	|        /  
		 *	*<--*<--*   *
		 */
		for (int i=0; i < ny; i++){ /* left */
			xBound[i] = 0;
			yBound[i] = i;
		}
		
		for (int i=1; i <  nx; i++){ /* bottom */
			xBound[ny-1 + i] = i;
			yBound[ny-1 + i] = ny-1;
		}
		
		for (int i=ny-2; i >= 0 ; i--){ /* right */
			xBound[2*ny + nx-3 - i] = nx-1;
			yBound[2*ny + nx-3 - i] = i;
		}
		
		for (int i=nx-2; i >= 0; i--){ /* top */
			xBound[2*ny + 2*nx-4 - i] = i;
			yBound[2*ny + 2*nx-4 - i] = 0;
		}
		return new int[][]{ xBound, yBound };
	}
	
	/** Returns the coefficients of the four corners of the grid cell
	 * to give the value of f at the specified fractional position (fx,fy) 
	 * consistent with the triangular grid linear approximation that the  
	 * contouring assumes.
	 * 
	 * @return coeffs double[]{ c_TL, c_TR, c_BL, c_BR }
	 */ 
	public static double[] getInterpolationCoefficients(double fx, double fy){
	    if(fx < 0 || fx > 1 || fy < 0 || fy > 1) //FIXME: Remove this for efficiency once checked!	        
	        throw new RuntimeException("Sanity check: getInterpolationCoefficients("+fx+","+fy+") called outside (0,0)-(1,1)");
	    
	    if(fx > fy){ //Top or Right
            if(fx > (1-fy)) //Right
                return new double[]{    (1 - fx)/2,
                                        (0.5 + fx/2 - fy),
                                        (1 - fx)/2,
                                        (fy - fx/2 - 0.5) };
            else //Top
                return new double[]{    (1 - fx - fy/2),
                                        (fx - fy/2),
                                        fy/2,
                                        fy/2 }; 
        }else{ //Bottom or Left
            if(fx > (1-fy)) //Bottom
                return new double[]{    (1 - fy)/2,
                                        (1 - fy)/2,
                                        (0.5 + fy/2 - fx),
                                        (fx + fy/2 - 0.5) };     
            else //Left
                return new double[]{    (1 - fy - fx/2),
                                        fx/2,
                                        (fy - fx/2),
                                        fx/2 };
        }
	    
	    
	}
	

	
	/** Returns which triangle (0=Top, 1=Left, 2=Right, 3=Bottom)
	 * the fractional position (fx, fy) is in.
	 * 
	 */
	public static final int whichTriangle(double fx, double fy){
	    if(fx > fy){ //Top or Right
            if(fx > (1-fy))
                return 2;
            else //Top
                return 0; 
            
        }else{ //Bottom or Left
            if(fx > (1-fy)) //Bottom
                return 3;     
            else //Left
                return 1;
            
        }
           
	}
	
	/** Returns the fractional coordinates of the corners of the triangle 
	 * in which the fractional position (fx,fy) lies.
	  * @return fracCoords double[][]{ { x0, x1, x2, x0}, { y0, y1, y2, y0 } } - clockwise
	 */
	public static final double[][] getTriangleVerticies(double fx, double fy){
	    if(fx > fy){ //Top or Right
            if(fx > (1-fy)) //Right
                return new double[][]{ { 1, 1, 0.5, 1 }, { 0, 1, 0.5, 0 } };
            else //Top
                return new double[][]{ { 0, 1, 0.5, 1 }, { 0, 0, 0.5, 0 } };
            
        }else{ //Bottom or Left
            if(fx > (1-fy)) //Bottom
                return new double[][]{ { 0.5, 1, 0, 0.5 }, { 0.5, 1, 1, 0.5 } };
            else //Left
                return new double[][]{ { 0, 0.5, 0, 0 }, { 0, 0.5, 1, 0 } };
            
        }
	}
	
	
	
}