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.
 * 
 * Uses the square grid algorithm and a simple Alex Meakins inspired case lookup table.
 * 
 * 
 * @author oliford <codes(at)oliford.co.uk>
 *
 */

public class Contouring {
   
	/** 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) */
	public static final int EDGE_TOP      = 0;
    public static final int EDGE_RIGHT      = 1;
    public static final int EDGE_BOTTOM		= 2;
    public static final int EDGE_LEFT		= 3;
	
	
	/* Return codes from contourFromSquare() (external and internal) */
	public static final int RETURN_FAILED = -3;
	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 TOP_DONE = (1 << 0);
    private final int LEFT_DONE = (1 << 1);

	private final boolean IS_LEFT_DONE(int x, int y){  return ((done[(y)*nx + (x)] & LEFT_DONE) != 0); }   
    private final boolean IS_TOP_DONE(int x, int y){   return ((done[(y)*nx + (x)] & TOP_DONE) != 0); }
    private final void SET_LEFT_DONE(int x, int y){ done[(y)*nx + (x)] |= LEFT_DONE; }
    private final void SET_TOP_DONE(int x, int y){ done[(y)*nx + (x)] |= TOP_DONE; }
    
    public interface GridFunction {
        double eval(int i, int j);
    };
    
    public static class FlatArrayGrid implements GridFunction {
        public double data[];
        int nx;
        public FlatArrayGrid(double data[], int nx) { this.data=data; this.nx = nx; }
        @Override
        public double eval(int i, int j) { return data[(j)*nx + (i)]; }
    }
    
    protected GridFunction f;
    
	protected final double FVAL(int x, int y){
		//
	    return f.eval(x, y);
	}

	private static final int X = -1;
	   
	/** Old leaving edge table and logic - has the very occasional failure */
	private static final int leavingEdgeTableOld[][] = new int[][]{
        // t   r   b   l   <-- entering edge
        {  X,  X,  X,  X }, /* 0:       - - / - -           */
        {  X,  2,  1,  X }, /* 1:       - - / - +           */
        {  X,  X,  3,  2 }, /* 2:       - - / + -           */
        {  X,  3,  X,  1 }, /* 3:       - - / + +           */
        {  1,  0,  X,  X }, /* 4:       - + / - -           */
        {  2,  X,  0,  X }, /* 5:       - + / - +           */
        {  1,  0,  3,  2 }, /* 6:(6a)   - + / + - center-   */
        {  3,  X,  X,  0 }, /* 7:       - + / - +           */
        {  3,  X,  X,  0 }, /* 8:       + - / - -           */
        {  3,  2,  1,  0 }, /* 9:(9a)   + - / - + center-   */
        {  2,  X,  0,  X }, /* 10:      + - / + -           */
        {  1,  0,  X,  X }, /* 11:      + - / + +           */
        {  X,  3,  X,  1 }, /* 12:      + + / - -           */
        {  X,  X,  3,  2 }, /* 13:      + + / - +           */
        {  X,  2,  1,  X }, /* 14:      + + / + -           */
        {  X,  X,  X,  X }, /* 15:      + + / + +           */
        
        {  3,  2,  1,  0 }, /* 16:(6b)  - + / + - center+   */
        {  1,  0,  3,  2 }, /* 17:(9b)  + - / - + center+   */
    };
    
    /** @returns the edge through which a contour, entering the given edge of the grid cell will leave.
     * Using oliford's logic - (which has very occasional failures) */
    public final int getLeavingEdgeOld(int i, int j, int enteringEdge, double contourValue){
        boolean tlAbove = FVAL(i, j) >= contourValue;
        boolean trAbove = FVAL(i+1, j) >= contourValue;
        boolean brAbove = FVAL(i+1, j+1) >= contourValue;
        boolean blAbove = FVAL(i, j+1) >= contourValue;
        
        int tabIdx = (tlAbove?8:0) + (trAbove?4:0) + (blAbove?2:0) + (brAbove?1:0);
        
       // System.out.print("Cell ("+i+", "+j+") is case "+tabIdx+". f="+contourValue+". Entering="+enteringEdge+"... ");
        
        //cases 6 and 9 are ambiguous but decided by the average, and may become 6b(16) or 9b(17)
        if(tabIdx == 6 && ((FVAL(i, j) + FVAL(i, j+1) + FVAL(i+1, j+1) + FVAL(i+1, j)) / 4) > contourValue)
            tabIdx = 16;
        else if(tabIdx == 9 && ((FVAL(i, j) + FVAL(i, j+1) + FVAL(i+1, j+1) + FVAL(i+1, j)) / 4) > contourValue)
               tabIdx = 17;
        
        // 0 and 15 have all nodes higher or all nodes lower, which probably means with have a node identically equal to contourValue
        // we have to deal with this specially
        if(tabIdx == 0 || tabIdx == 15){
            //quick hack, just nudge the contour value toward the mean
            double fC = (FVAL(i, j) + FVAL(i, j+1) + FVAL(i+1, j+1) + FVAL(i+1, j)) / 4;
            return getLeavingEdgeOld(i, j, enteringEdge, contourValue + 1e-6 * (fC - contourValue));
        }
            
                
        int ret = leavingEdgeTableOld[tabIdx][enteringEdge];
        
        if(ret == X){
            debugDumpState(i, j, enteringEdge, contourValue);
            throw new RuntimeException("Contouring got lost. Shouldn't be entering edge "+enteringEdge+" of a square type " + tabIdx + ".\n");
        }
        
        //System.out.println("OK. Leaving="+ret);
        
        return ret;
    }
    
    /** Leaving edge table replacement, curtesy of Alex Meakins. */
    private static final int leavingEdgeMainTable[][] = new int[][]{
        // t   r   b   l   <-- entering edge
        {  X,  X,  X,  X   },  // 0
        {  3,  X,  X,  0   },  // 1
        {  1,  0,  X,  X   },  // 2
        {  X,  3,  X,  1   },  // 3
        {  X,  2,  1,  X   },  // 4
        {  4,  4,  4,  4   },  // 5 (ambiguous)
        {  2,  X,  0,  X   },  // 6
        {  X,  X,  3,  2   },  // 7
        {  X,  X,  3,  2   },  // 8
        {  2,  X,  0,  X   },  // 9
        {  5,  5,  5,  5   },  // 10 (ambiguous)
        {  X,  2,  1,  X   },  // 11
        {  X,  3,  X,  1   },  // 12
        {  1,  0,  X,  X   },  // 13
        {  3,  X,  X,  0   },  // 14
        {  X,  X,  X,  X   },  // 15
        
    };

   private static final int leavingEdgeAmbiguousTable[][] = new int[][]{
        // t   r   b   l   
        {  3,  2,  1,  0   },  // 5a, 10b
        {  1,  0,  3,  2   },  // 5b, 10a
    };

    
   /** @returns the edge through which a contour, entering the given edge of the grid cell will leave.
    * Using Alex Meakin's logic - have not yet seen any failures */
    public final int getLeavingEdge(int i, int j, int enteringEdge, double contourValue){
        boolean tlLower = FVAL(i, j) < contourValue;
        boolean trLower = FVAL(i+1, j) < contourValue;
        boolean brLower = FVAL(i+1, j+1) < contourValue;
        boolean blLower = FVAL(i, j+1) < contourValue;
        
        int tabIdx = (blLower?8:0) + (brLower?4:0) + (trLower?2:0) + (tlLower?1:0);

        int ret = leavingEdgeMainTable[tabIdx][enteringEdge];
        
        // handle ambiguous cases
        if ((ret & 4) != 0) {
            if (ret == 4) {
                // case 5
                if (((FVAL(i, j) + FVAL(i, j+1) + FVAL(i+1, j+1) + FVAL(i+1, j)) / 4) < contourValue)
                    ret = leavingEdgeAmbiguousTable[1][enteringEdge];
                else
                    ret = leavingEdgeAmbiguousTable[0][enteringEdge];
                
            }
            else {
                // case 10
                if (((FVAL(i, j) + FVAL(i, j+1) + FVAL(i+1, j+1) + FVAL(i+1, j)) / 4) < contourValue)
                    ret = leavingEdgeAmbiguousTable[0][enteringEdge];
                else
                    ret = leavingEdgeAmbiguousTable[1][enteringEdge];               
            }
        }
        
        //debugDumpState(i, j, enteringEdge, contourValue);
        if((ret <= X) || (ret > 3)){
            debugDumpState(i, j, enteringEdge, contourValue);
            throw new RuntimeException("Contouring got lost. Shouldn't be entering edge "+enteringEdge+" of a square type " + tabIdx + ".\n");
        }
        
        //System.out.println("OK. Leaving="+ret);
        
        return ret;
    }
	
	/** Begins contouring from the given square by looking at each edge in turn.
	 * @param i
	 * @param j
	 * @param contourValue
	 * @return
	 */
	public int contourFromSquare(int i, int j, double contourValue){
	    int enteringEdge;
	    
	    for(int edge=0; edge < 3; edge++){
	        if(!isEdgeDone(i, j, edge) && isOnEdge(i, j, edge, contourValue) )
	            return contourFromEdge(i, j, edge, contourValue);
	        
	    }
	    
	    return RETURN_ABORTED;
	}
	
	public boolean isEdgeDone(int i, int j, int edge){
	    switch(edge){
            case EDGE_TOP: return IS_TOP_DONE(i, j);
            case EDGE_RIGHT: return IS_LEFT_DONE(i+1, j);
            case EDGE_BOTTOM: return IS_TOP_DONE(i, j+1);
            case EDGE_LEFT: return IS_LEFT_DONE(i, j);
	    }
	    return false;
    
	}

	/** Returns if the given function value exists on the given edge of the given cell */
	public boolean isOnEdge(int i, int j, int enteringEdge, double contourValue){
        double a,b;
        switch(enteringEdge){
            case EDGE_TOP: a = FVAL(i+1, j); b = FVAL(i, j); break;
            case EDGE_RIGHT: a = FVAL(i+1, j+1); b = FVAL(i+1, j); break;
            case EDGE_BOTTOM: a = FVAL(i+1, j+1); b =  FVAL(i, j+1); break;
            case EDGE_LEFT: a = FVAL(i, j+1); b = FVAL(i, j); break;
            default:
                throw new RuntimeException("Unrecognised edge type '"+enteringEdge+"'.");
        }  
        
        return contourValue > Math.min(a, b) && contourValue <= Math.max(a, b);
        
        /*switch(enteringEdge){
            case EDGE_TOP: f = (contourValue - FVAL(i, j)) / (FVAL(i+1, j) - FVAL(i, j)); break;
            case EDGE_RIGHT: f = (contourValue - FVAL(i+1, j)) / (FVAL(i+1, j+1) - FVAL(i+1, j)); break;
            case EDGE_BOTTOM: f = (contourValue - FVAL(i, j+1)) / (FVAL(i+1, j+1) - FVAL(i, j+1)); break;
            case EDGE_LEFT: f = (contourValue - FVAL(i, j)) / (FVAL(i, j+1) - FVAL(i, j)); break;
            default:
                throw new RuntimeException("Unrecognised edge type '"+enteringEdge+"'.");
        }
        
        return (f >= 0 && f < 1);*/
    
	}
	
	public double getCrossingPoint(int i, int j, int edge, double contourValue, double hitPos[], boolean checkOnEdge){
	    double f;
        //mark at point here
        switch(edge){
            case EDGE_TOP:
                f = (contourValue - FVAL(i, j)) / (FVAL(i+1, j) - FVAL(i, j));
                hitPos[0] = x[i] + f * (x[i+1]-x[i]);
                hitPos[1] = y[j];
                break;
            case EDGE_RIGHT:
                f = (contourValue - FVAL(i+1, j)) / (FVAL(i+1, j+1) - FVAL(i+1, j));
                hitPos[0] = x[i+1];
                hitPos[1] = y[j] + f * (y[j+1]-y[j]);
                break;
            case EDGE_BOTTOM:
                f = (contourValue - FVAL(i, j+1)) / (FVAL(i+1, j+1) - FVAL(i, j+1));
                hitPos[0] = x[i] + f * (x[i+1]-x[i]);
                hitPos[1] = y[j+1];
                break;
            case EDGE_LEFT:
                f = (contourValue - FVAL(i, j)) / (FVAL(i, j+1) - FVAL(i, j));;
                hitPos[0] = x[i];
                hitPos[1] = y[j] + f * (y[j+1]-y[j]);
                break;
            default:
                    throw new RuntimeException("Unrecognised edge type '"+edge+"'.");
        }
        
        //input check
        if(checkOnEdge && (f < 0 || f > 1)){
            debugDumpState(i, j, edge, contourValue);
            throw new IllegalArgumentException("Edge "+edge+" of cell ("+i+","+j+") doesn't contain the level f="+contourValue);
        }
        
        return f;
	}
	
	/** Begins contouring from the given edge of the given cell. */
	public int contourFromEdge(int i0, int j0, int enteringEdge, double contourValue){
	    int i = i0;
	    int j = j0;
	    double hitPos[] = new double[2];
	    
        //mark the starting edge visited
        switch (enteringEdge) {
            case EDGE_TOP: SET_TOP_DONE(i, j); break;
            case EDGE_RIGHT: SET_LEFT_DONE(i+1, j); break;
            case EDGE_BOTTOM: SET_TOP_DONE(i, j+1); break;
            case EDGE_LEFT: SET_LEFT_DONE(i, j); break;            
        }
	    
	    do{
	        if( Double.isNaN(FVAL(i, j)) || Double.isNaN(FVAL(i+1, j)) 
	                || Double.isNaN(FVAL(i+1, j+1)) || Double.isNaN(FVAL(i, j+1)) ){
	            System.err.println("WARNING: Contouring entered grid square with NaN value.");
	            callbacks.contourEnd(contourValue, RETURN_FAILED);
                return RETURN_FAILED;
	        }
	        	            
	        getCrossingPoint(i, j, enteringEdge, contourValue, hitPos, true);
	        //register the point
	        if(!callbacks.contourPoint(contourValue, hitPos[0], hitPos[1])){
	            //if they request, abort the contour and return aborted
	            callbacks.contourEnd(contourValue, RETURN_ABORTED);
                return RETURN_ABORTED;
	        }
	        
	        int leavingEdge = getLeavingEdge(i, j, enteringEdge, contourValue);
/*	        int leavingEdge2 = getLeavingEdgeOld(i, j, enteringEdge, contourValue);
	        if(leavingEdge != leavingEdge2){
	            debugDumpState(i, j, enteringEdge, contourValue);
	            System.err.println("Alex's got "+leavingEdge+", oli's got "+leavingEdge2);
	            leavingEdge = getLeavingEdge(i, j, enteringEdge, contourValue);
	        }
*/
	        //move to that cell and set the new entering edge
	        boolean closed = false;
	        switch(leavingEdge){
    	        case EDGE_TOP: closed = IS_TOP_DONE(i, j); j--; enteringEdge = EDGE_BOTTOM; break;
    	        case EDGE_RIGHT: closed = IS_LEFT_DONE(i+1, j); i++; enteringEdge = EDGE_LEFT; break;
    	        case EDGE_BOTTOM: closed = IS_TOP_DONE(i, j+1); j++; enteringEdge = EDGE_TOP; break;
    	        case EDGE_LEFT: closed = IS_LEFT_DONE(i, j); i--; enteringEdge = EDGE_RIGHT; break;            
	        }
	        
	        //mark the new edge visited
            switch (enteringEdge) {
                case EDGE_TOP: SET_TOP_DONE(i, j); break;
                case EDGE_RIGHT: SET_LEFT_DONE(i+1, j); break;
                case EDGE_BOTTOM: SET_TOP_DONE(i, j+1); break;
                case EDGE_LEFT: SET_LEFT_DONE(i, j); break;            
            }
	        
	        //check we've not left the grid
	        if(i < 0 || i >= (nx-1) || j < 0 || j >= (ny-1) ){
	            getCrossingPoint(i, j, enteringEdge, contourValue, hitPos, true);
                if(!callbacks.contourPoint(contourValue, hitPos[0], hitPos[1])){
                    callbacks.contourEnd(contourValue, RETURN_ABORTED);
                    return RETURN_ABORTED;
                }
	            if(!callbacks.contourEnd(contourValue, RETURN_HIT_GRID_BOUNDARY))
	                return RETURN_ABORTED;
	            return RETURN_HIT_GRID_BOUNDARY;
	        }
	        
	        //if this square has been visited, we've closed it or we're dealing with a split contour
	        if( closed ){
	            int retCode = (i==i0 && j==j0) ? RETURN_CLOSED : RETURN_ABORTED; 
	            getCrossingPoint(i, j, enteringEdge, contourValue, hitPos, true);
	            
	            if(!callbacks.contourPoint(contourValue, hitPos[0], hitPos[1])){
	                retCode = RETURN_ABORTED;
	            }
	            if(!callbacks.contourEnd(contourValue, retCode))
                    return RETURN_ABORTED;
	            
                return retCode;
	        }
	        
	    }while(true);
	    
	}
	
	private void debugDumpState(int i, int j, int edge, double contourValue){
	    double fC = (FVAL(i, j) + FVAL(i+1, j) + FVAL(i+1, j+1) + FVAL(i, j+1)) / 4;
	    System.err.println("edge = " + edge + "\ncontourValue = " + contourValue + "\n" +
            "tl: f= "+ FVAL(i, j) + "\t\tdiff=" + (FVAL(i, j) - contourValue) + "\n" +
            "tr: f= "+ FVAL(i+1, j) + "\t\tdiff=" + (FVAL(i+1, j) - contourValue) + "\n" +
            "br: f= "+ FVAL(i+1, j+1) + "\t\tdiff=" + (FVAL(i+1, j+1) - contourValue) + "\n" +
            "bl: f= "+ FVAL(i, j+1) + "\t\tdiff=" + (FVAL(i, j+1) - contourValue) + "\n" + 
            "center: " + fC + "\t\tdiff=" + (fC - contourValue) + "\n");
	    
	}
		
	/** 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 Contouring(ContourCallback callbacks,double x0, double x1, int nx,double y0, double y1, int ny,GridFunction f) {
        this(callbacks, linearSpace(x0, x1, nx), linearSpace(y0, y1, ny), f);
    }
	
	/** Creates contouring object based on regular grid of data.     
     * 
     * @param f values at grid points (flat - x is deepest dimension)
     */
    public Contouring(ContourCallback callbacks,double x0, double x1, int nx,double y0, double y1, int ny,double f[]) {
        this(callbacks, linearSpace(x0, x1, nx), linearSpace(y0, y1, ny), f);
    }
    
    /** Creates contouring object based on regular grid of data.
     * @param f values at grid points ( is deepest dimension)
     */
    public Contouring(ContourCallback callbacks, double x0, double x1, int nx, double y0, double y1, int ny, double f[][]) {
        this(callbacks, linearSpace(x0, x1, nx), linearSpace(y0, y1, ny), f);        
    }

	/** 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 Contouring(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;
		if(f != null && f.length != (nx*ny))
            throw new IllegalArgumentException("f[] should to be (x.length * y.length) long.");
		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 (x is deepest dimension)
	 */
	public Contouring(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 Contouring(ContourCallback callbacks, double x[], double y[], GridFunction 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 = new FlatArrayGrid(f, nx);
	}
	
	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)");
        
        double f1d[] = new double[ny*nx];
        for(int i=0;i<ny;i++)
            System.arraycopy(f[i], 0, f1d, i*nx, nx);
        
        setF(f1d);
	}
	
	public final void setF(GridFunction f){
	    this.f = f;
	}
	
	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 contourValue) {
			
	    // top row
	    for(int i=0; i < (nx-1); i++)
	        if(!IS_TOP_DONE(i, 0) && isOnEdge(i, 0, EDGE_TOP, contourValue))
	            contourFromEdge(i, 0, EDGE_TOP, contourValue);

	    // right column
	    for(int j=0; j < (ny-1); j++)
	        if(!IS_LEFT_DONE((nx-1), j) && isOnEdge((nx-2), j, EDGE_RIGHT, contourValue))
	            contourFromEdge((nx-2), j, EDGE_RIGHT, contourValue);

	    // bottom row
	    for(int i=0; i < (nx-1); i++)
	        if(!IS_TOP_DONE(i, (ny-1)) && isOnEdge(i, (ny-2), EDGE_BOTTOM, contourValue))
	            contourFromEdge(i, (ny-2), EDGE_BOTTOM, contourValue);

	    // left column
	    for(int j=0; j < (ny-1); j++)
	        if(!IS_LEFT_DONE(0, j) && isOnEdge(0, j, EDGE_LEFT, contourValue))
	            contourFromEdge(0, j, EDGE_LEFT, contourValue);
                         
	}
		
	/** Gets the rest of the contours (closed contours) after getBoundaryContours() */
	public void getRemainingContours(double contourValue) {
			int i,j;

			/** inside segments **/
			for(i=0; i < (nx-1); i++)
				for(j=0; j < (ny-1); j++)
				    contourFromSquare(i, j, contourValue);
					
	}
	
	
}