package org.jzy3d.maths;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.util.Vector;

import org.apache.log4j.Logger;
import org.jzy3d.colors.Color;
import org.jzy3d.colors.IColorMappable;
import org.jzy3d.colors.colormaps.IColorMap;
import org.jzy3d.plot3d.primitives.Point;
import org.jzy3d.plot3d.primitives.Polygon;



/** The GridLoader utility offers a way to convert coordinates mode.
 * When input data is made of lists of x, y, and z points, and 
 * that these lists actually contain data that regularly lie
 * on a squared grid (orthonormal grid), the GridLoader allows 
 * a smart retrieving of data by providing:
 * 

 * 
an 1D array of sorted X coordinates (without doublon).
 * 
an 1D array of sorted Y coordinates (without doublon).
 * 
an 2D array of Z coordinates syncronized with arrays X and Y.
 * 
 * 

 * If some data are missing for representing this grid, missing data 
 * will be considered as NaN. The actual management of missing values
 * is thus left to the object that loads these data.
 * 

 * In the following example, the o represent a NaN Z value, meaning that
 * there was no (x[i],y[i],z[i]) triplet for representing it correctly:
 * 

 * 
 * 11111111 y

 * 11111111  

 * 11o11111  

 * 11111111  

 *      

 * x         

 *           

 * 
 * 

 *
 * If some points are missing in the input list to generate a complete
 * grid, the output grid will contain Z=NaN values for missing coordinates.
 *
 * @author Martin Pernollet
 */
@Deprecated
public class GridLoader implements IColorMappable{
    static Logger logger = Logger.getLogger(GridLoader.class);

	/** Initialize a GridLoader by parsing the square grid directly.
	 * When the grid is parsed, the GridLoader may return either:
	 * 

	 * 
a set of polygon describing a ring with smooth contour
	 * 

	 */
	public GridLoader(float x[], float y[], float z[]){
		bbox = new BoundingBox3d();
		setData(x, y, z);
	}
	
	/**********************************************************************/

	
	/** Return a BoundingBox indicating the bounds of the loaded data.
	 * @return a bounding box.
	 */
	public BoundingBox3d getBoundingBox(){
		return bbox;
	}
	
	/************************************************************************************************/

	/** Set the array of data. X,Y, and Z are arrays that must implicitely
	 * represent an orthonormal grid. 
	 * 
	 * If some data are missing for representing this grid, missing data 
	 * will be considered as NaN. The actual management of missing values
	 * is thus left to the object that loads these data.
	 * 
	 * In the following example, the o represent a NaN Z value, meaning that
	 * there was no (x[i],y[i],z[i]) triplet for representing it correctly:
	 * 

	 * 
	 * 11111111 y

	 * 11111111  

	 * 11o11111  

	 * 11111111  

	 *      

	 * x         

	 *           

	 * 
	 * 
 
	 * 
	 * While parsing input data, the minimal and maximal Z values are computed, 
	 * in order to have the loader ready for computing polygons color with a 
	 * colormap.
	 * 
	 * @param x list of x coordinates
	 * @param y list of y coordinates
	 * @param z list of z coordinates
	 * @throws an IllegalArgumentException if x, y , and z have not the same size
	 */
	protected void setData(float x[], float y[], float z[]){
		if(x.length!=y.length || x.length!=z.length)
			throw new IllegalArgumentException("GridLoader: x, y, and z arrays must agree in length.");
		
		// Initialize loading
		bbox.reset();
		this.x = unique(x);
		this.y = unique(y);
		this.z = new float[this.x.length][this.y.length];
		
		for(int i=0; i getSquarePolygons(IColorMap cmap, Color colorFactor){
		Vector polygons = new Vector();
		bbox.reset();
		
		boolean[] isIn; 

		for(int xi=1; xi getExtrapolatedRingPolygons(float ringMax, IColorMap cmap, Color colorFactor){
		//backup current coords and extrapolate
		float[]   xbackup = x;
		float[]   ybackup = y;
		float[][] zbackup = z;
		
		// compute required extrapolation
		float step = x[1] - x[0];
		int  nstep = x.length;
		
		int ENLARGE = 2;		
		int required = (int)( Math.ceil( (ringMax*2 - step*nstep) / step )  );
		required = required<0 ? ENLARGE:required+ENLARGE;		
		
		if(required>0)
			extrapolate(required);
		
		List polygons = getInterpolatedRingPolygons(0, ringMax, cmap, colorFactor);
		
		// get back to previous grid
		x = xbackup;
		y = ybackup;
		z = zbackup;
		
		return polygons;
	}
	
	/** Add extrapolated points on the grid. If the grid is too small for extrapolation, the arrays
	 * are maximized */
	public void extrapolate(int n){
		float[]   xnew = new float[x.length+n*2];
		float[]   ynew = new float[y.length+n*2];
		float[][] znew = new float[x.length+n*2][y.length+n*2];
		
		// assume x and y grid are allready sorted and create new grids
		float xmin = x[0];
		float xmax = x[x.length-1];
		float xgap = x[1]-x[0];		
		float ymin = y[0];
		float ymax = y[y.length-1];
		float ygap = y[1]-y[0];
		
		for(int i=0; i=n && i= x.length+n) // fill after
				xnew[i] = xmax + (i-(x.length+n)+1) * xgap;
			
			// --- y grid ---
			for(int j=0; j=n && j< (y.length+n) ){ // copy content
					ynew[j]    = y[j-n];
					
					// copy z grid
					if(i>=n && i= (y.length+n) ){ // fill after
					ynew[j]    = ymax + (j-(y.length+n)+1) * ygap;
					znew[i][j] = Float.NaN;
				}
			}
		}
		
		// extrapolation
		float olddiameter = xgap*(x.length    )/2;
		float newdiameter = xgap*(x.length-1+n*2)/2;
		olddiameter *= olddiameter;
		newdiameter *= newdiameter;
		
		int xmiddle = (xnew.length-1)/2; // assume it is an uneven grid
		int ymiddle = (ynew.length-1)/2; // assume it is an uneven grid		
		
		// start from center, and add extrapolated values iteratively on each quadrant
		for(int i=xmiddle; i= olddiameter){
					
					int xopp = i - 2*(i-xmiddle);
					int yopp = j - 2*(j-ymiddle);
					
					znew[i]   [j]    = getExtrapolatedZ(znew, i, j); // right up quadrant
					znew[xopp][j]    = getExtrapolatedZ(znew, xopp, j); // left  up
					znew[i]   [yopp] = getExtrapolatedZ(znew, i, yopp); // right down
					znew[xopp][yopp] = getExtrapolatedZ(znew, xopp, yopp); // left  down
				}
				
				// ignore values standing outside desired diameter
				else //if(sqrad > newdiameter)
					znew[i][j] = Float.NaN;	
			}
		}
		
		// store result
		x = xnew;
		y = ynew;
		z = znew;
	}
	
	private float getExtrapolatedZ(float[][] grid, int currentXi, int currentYi){
		int left   = currentXi-1>0              ? currentXi-1 : currentXi;
		int right  = currentXi+10              ? currentYi-1 : currentYi;
		int up     = currentYi+10)
			return cumval/nval;
		else
			return Float.NaN;
	}

	
	/**************************************************************************************/
	
	/** Load data standing on an orthonormal grid.
	 * 

	 * Each input point (i.e. the association of x[i], y[j], z[i][j]) will be
	 * represented by a polygon centered on this point. The default coordinates
	 * of this polygon will be:
	 * 

	 * 
x[i-1], y[j+1], z[i-1][j+1]
	 * 
x[i-1], y[j-1], z[i-1][j-1]
	 * 
x[i+1], y[j-1], z[i+1][j-1]
	 * 
x[i+1], y[j+1], z[i+1][j+1]
	 * 
	 * There are thus three types of polygons:
	 * 

	 * 
those that stand completely inside the ringMin and ringMax radius and
	 * that have the previous coordinates.
	 * 
those that stand completely outside the ringMin and ringMax radius and
	 * that won't be added to the list of polygons.
	 * 
those that have some points in and some points out of the ringMin and 
	 * ringMax radius. These polygons are recomputed so that "out" points are replaced
	 * by two points that make the smooth contour. According to the number of "out"
	 * points, the modified polygon will gather 3, 4, or 5 points.
	 * 

	 * As a consequence, it is suggested to provide data ranging outside of ringMin
	 * and ringMax, in order to be sure to have a perfect round surface.
	 * 
	 * @param ringMin the minimum radius of this ring.
	 * @param ringMax the maximum radius of this ring.
	 * @param colorFactor a weighting factor for the color returned by the Colormap.
	 */
	public List getInterpolatedRingPolygons(float ringMin, float ringMax, IColorMap cmap, Color colorFactor){
		List polygons = new ArrayList();
		bbox.reset();
		
		boolean[] isIn; 

		for(int xi=1; xi= minRadius;
		isIn[1] = !Float.isNaN(p[1].xyz.z) && radius[1] < maxRadius && radius[1] >= minRadius;
		isIn[2] = !Float.isNaN(p[2].xyz.z) && radius[2] < maxRadius && radius[2] >= minRadius;
		isIn[3] = !Float.isNaN(p[3].xyz.z) && radius[3] < maxRadius && radius[3] >= minRadius;
		
		return isIn;
	}
	
	private float radius2d(Point p){
		return (float)Math.sqrt(p.xyz.x*p.xyz.x + p.xyz.y*p.xyz.y);
	}
	
	/**
	 * Return a point that is the intersection between a segment and a circle
	 * @throws ArithmeticException if points do not stand on an squared (orthonormal) grid.
	 */
	private Point findPoint(Point p1, Point p2, float ringRadius){		
		// We know that the seeked point is on a horizontal or vertial line
		float x3, y3, z3;
		float w1, w2;
		double alpha;
		
		// We know x3 and radius and seek y3, using intermediate alpha
		if(p1.xyz.x==p2.xyz.x){
			x3 = p1.xyz.x;
			alpha = Math.acos(x3/ringRadius);
			
			if(p1.xyz.y<0 && p2.xyz.y<0)
				y3 = -(float)Math.sin(alpha)*ringRadius;
			else if(p1.xyz.y>0 && p2.xyz.y>0)
				y3 =  (float)Math.sin(alpha)*ringRadius;
			else if(p1.xyz.y==-p2.xyz.y)
				y3 = 0; // ne peut pas arriver
			else
				throw new ArithmeticException(("no alignement between p1(" +p1.xyz.x + "," + p1.xyz.y + "," + p1.xyz.z + ") and p2(" +p2.xyz.x + "," + p2.xyz.y + "," + p2.xyz.z + ")"));
			
			// and now get z3
			if(!Float.isNaN(p1.xyz.z) && Float.isNaN(p2.xyz.z))
				z3 = p1.xyz.z;
			else if(Float.isNaN(p1.xyz.z) && !Float.isNaN(p2.xyz.z))
				z3 = p2.xyz.z;
			else if(!Float.isNaN(p1.xyz.z) && !Float.isNaN(p2.xyz.z)){					
				w2 = (float)( Math.sqrt( (x3      -p1.xyz.x)*(x3      -p1.xyz.x) + (y3      -p1.xyz.y)*(y3      -p1.xyz.y) )
				            / Math.sqrt( (p2.xyz.x-p1.xyz.x)*(p2.xyz.x-p1.xyz.x) + (p2.xyz.y-p1.xyz.y)*(p2.xyz.y-p1.xyz.y) ));
				w1 = 1-w2;
				z3 = w1*p1.xyz.z+w2*p2.xyz.z;
			}
			else
				throw new ArithmeticException(("can't compute z3 with p1(" +p1.xyz.x + "," + p1.xyz.y + ") and p2(" +p2.xyz.x + "," + p2.xyz.y + ")"));
		}
		// We know y3 and radius and seek x3, using intermediate alpha
		else if(p1.xyz.y==p2.xyz.y){
			y3 = p1.xyz.y;
			alpha = Math.asin(y3/ringRadius);
			
			if(p1.xyz.x<0 && p2.xyz.x<0)
				x3 = -(float)Math.cos(alpha)*ringRadius;
			else if(p1.xyz.x>0 && p2.xyz.x>0)
				x3 =  (float)Math.cos(alpha)*ringRadius;
			else if(p1.xyz.x==-p2.xyz.x)
				x3 = 0; // ne peut pas arriver
			else
				throw new ArithmeticException(("no alignement between p1(" +p1.xyz.x + "," + p1.xyz.y + "," + p1.xyz.z + ") and p2(" +p2.xyz.x + "," + p2.xyz.y + "," + p2.xyz.z + ")"));
			
			// and now get z3
			if(!Float.isNaN(p1.xyz.z) && Float.isNaN(p2.xyz.z))
				z3 = p1.xyz.z;
			else if(Float.isNaN(p1.xyz.z) && !Float.isNaN(p2.xyz.z))
				z3 = p2.xyz.z;
			else if(!Float.isNaN(p1.xyz.z) && !Float.isNaN(p2.xyz.z)){		
				w2 = (float)( Math.sqrt( (x3      -p1.xyz.x)*(x3      -p1.xyz.x) + (y3      -p1.xyz.y)*(y3      -p1.xyz.y) )
			                / Math.sqrt( (p2.xyz.x-p1.xyz.x)*(p2.xyz.x-p1.xyz.x) + (p2.xyz.y-p1.xyz.y)*(p2.xyz.y-p1.xyz.y) ));
				w1 = 1-w2;
				z3 = w1*p1.xyz.z+w2*p2.xyz.z;
			}
			else
				throw new ArithmeticException(("can't compute z3 with p1(" +p1.xyz.x + "," + p1.xyz.y + ") and p2(" +p2.xyz.x + "," + p2.xyz.y + ")"));
		}
		else
			throw new ArithmeticException(("no alignement between p1(" +p1.xyz.x + "," + p1.xyz.y + ") and p2(" +p2.xyz.x + "," + p2.xyz.y + ")"));
		
		return new Point(new Coord3d(x3, y3, z3));
	}
	
	/**********************************************************************/
	
	@Override
    public double getMin(){
		return colorZmin;
	}
	
	@Override
    public double getMax(){
		return colorZmax;
	}
	
	@Override
    public void setMin(double zmin){
		colorZmin = (float)zmin;
	}
	
	@Override
    public void setMax(double zmax){
		colorZmax = (float)zmax;
	}
		
	/**********************************************************************/
	
	private float colorZmin = Float.NaN;
	private float colorZmax = Float.NaN;
	
	public float x[];
	public float y[];
	public float z[][];
	
	private BoundingBox3d bbox;
	
	private int   findxi;
	private int   findyj;
}