/*
 * Copyright (C) 2007 askywhale
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License (version 2)
 *  as published by the Free Software Foundation.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 */

package com.askywhale.rings;

import java.util.Iterator;

import net.sf.j3d.threeD.raytracer.engine.*;
import net.sf.j3d.threeD.raytracer.shaders.*;
import net.sf.j3d.util.graphics.RGB;
import net.sf.j3d.util.Vector;

/**
 * 
 * Rings Raytracer related class
 * 
 * Isosurface class ; extend it with a getValue(Vector) that
 * return a double around 0 ; thresold is 0 and a well dimensionned
 * isosurface should be close of a sphere of a radius of 1
 * (i.e. value = sqrt(x*x+y*y+z*z)-1 )
 * positive : front side
 * 
 * 0 is positive, the limit is the point between positive and 
 * negative numbers
 * 
 * There are 3 (understandable) quality parameters, ajust them
 * at purpose (bigger is better) 
 * 
 * @author askywhale
 *
 */
public abstract class IsoSurface extends AbstractSurface {
    protected final int DEFAULT_Q_MAX_INTERSECT = 2;
    protected final double DEFAULT_Q_MAX_RELATIVE_SIZE = 2;
    protected final double DEFAULT_Q_PRECISION = 100;
    protected final boolean DEFAULT_DISCRETE_SHADE_MODE = false;
    
    protected int qMaxIntersect;
    protected double qMaxRelativeSize;
    protected double qPrecision;
    protected double qNormalPrecision;    
    protected boolean discreteShadeMode;
    
    public IsoSurface(Vector location, double size) {
        super(location, size);
        qMaxIntersect = DEFAULT_Q_MAX_INTERSECT;
        qMaxRelativeSize = DEFAULT_Q_MAX_RELATIVE_SIZE;
        qPrecision = DEFAULT_Q_PRECISION;
        discreteShadeMode = DEFAULT_DISCRETE_SHADE_MODE;
        calculateQNormalPrecision();
    }
    
    public IsoSurface(Vector location, double size, RGB color) {
        super(location, size, color);
        qMaxIntersect = DEFAULT_Q_MAX_INTERSECT;
        qMaxRelativeSize = DEFAULT_Q_MAX_RELATIVE_SIZE;
        qPrecision = DEFAULT_Q_PRECISION;
        calculateQNormalPrecision();
    }
    
    public Vector getNormalAt(Vector point) {    	
    	if(discreteShadeMode)
    		return getNormalAtDiscrete(point);
    	else
    		return getNormalAtContinue(point);
    }
    
    public Vector getNormalAtDiscrete(Vector point) {
//    	the best normal is the smaller one with all (v2-v1)<>0
//    	we start with a very small one, then grow and ajust,
//      making an average of smaller ones
    	Vector rpoint = super.getTransform(true).getInverse().transformAsLocation(point);
    	final int maxSteps = (int)(1+qNormalPrecision/30);    	
    	double mlt = Math.pow(100,1./maxSteps);
    	double quickMlt = 10*mlt;
    	double refD = (qMaxIntersect+rpoint.length())/2;
    	double baseD = refD*1e-15/qNormalPrecision; 
    	double d = baseD;
    	Vector vect = new Vector(0,0,0);
    	double baseValue = getValue(rpoint);
    	Vector rebasePoint = null;
    	double rebaseValue = 1e100;
    	boolean notAllEqualFound = false;
    	for(int i=0;i<maxSteps;i++) {
//    		double _optSqrtD = Math.sqrt(d);
    		boolean allEquals = true;
    		for(int x=-1;x<=1;x++) {
    			for(int y=-1;y<=1;y++) {
    				for(int z=-1;z<=1;z++) {
    					if(x!=0 || y!=0 || z!=0) { 
	   						Vector nVect = rpoint.add(new Vector(x*d,y*d,z*d));
	   						double value = getValue(nVect);
	   						double dValue = value - baseValue;
	   						if(dValue!=0) {
	   							allEquals = false;
		   						if(value<rebaseValue && value >= 0 && baseValue<0) {
		   							rebasePoint = nVect;
		   							rebaseValue = value;
		   						} 
		   						vect = vect.add(new Vector(
		    						x*d*dValue,y*d*dValue,z*d*dValue
//		   							x*_optSqrtD*dValue,y*_optSqrtD*dValue,z*_optSqrtD*dValue		   								
		    					));
	   						}
    					}
    				}
    			}
    		}
			
			//search only from a positive value (=>basepoints all on same side)
			if(rebasePoint!=null && baseValue<0) {
				rpoint = rebasePoint;
				d = baseD;
				vect = new Vector(0,0,0);
				//System.out.println("rebasing ("+allEqNb+","+mlt+") from "+baseValue+" to "+rebaseValue);
				baseValue = rebaseValue;
				i = 0;				
			} 
			else
				if(allEquals)
					i--;			
			if(!allEquals && !notAllEqualFound) {
				notAllEqualFound = true;
				d /= quickMlt;
			}
			if(notAllEqualFound)
				d *= mlt;
			else
				d *= quickMlt;	
    	}
    	Vector normal = vect; 
    	normal = normal.divide(normal.length());
    	normal = getTransform(true).transformAsNormal(normal).
        	multiply(getSize());
    	return normal;
    }
    
    public Vector getNormalAtContinue(Vector point) {
//the best normal is the smaller one with all (v2-v1)<>0
//we start with a very small one, then grow while v2=v1, going
//everywhere around and making a sum of gradiants
    	Vector rpoint = super.getTransform(true).getInverse().transformAsLocation(point);
       	final int maxSteps = (int)qNormalPrecision;
    	final double minRel = .01/qNormalPrecision;
       	final double maxRel = 2*qMaxRelativeSize;    	
    	double d = rpoint.length()*minRel;
    	double dx = 0;
    	double dy = 0;
    	double dz = 0;
    	double mlt = Math.pow(maxRel/minRel,1/(double)maxSteps); 
    	boolean oneOrMoreIsNull = true;
		boolean twoOrMoreIsNull = true; 
		boolean allIsNull = true;
    	int nb = 0;
    	int nbNotAllIsNull = 0;
    	int nbNotTwoIsNull = 0;
		
    	while(oneOrMoreIsNull && 
    			nbNotTwoIsNull < Math.pow(maxSteps,.3) &&
    			nbNotAllIsNull < Math.sqrt(maxSteps) &&
    			nb < maxSteps) {
    		dx = getValue(rpoint.add(new Vector(d,0,0)))-
				getValue(rpoint.add(new Vector(-d,0,0)));
    		dy = getValue(rpoint.add(new Vector(0,d,0)))-
				getValue(rpoint.add(new Vector(0,-d,0)));
    		dz = getValue(rpoint.add(new Vector(0,0,d)))-
				getValue(rpoint.add(new Vector(0,0,-d)));
    		oneOrMoreIsNull = dx*dy*dz==0;
    		twoOrMoreIsNull = dx*dy==0 && dx*dz==0 && dy*dz==0; 
    		allIsNull = dx == 0 && dy == 0 && dz == 0;
    		d *= mlt;
    		if(!allIsNull)
    			nbNotAllIsNull ++;
    		else
    			nbNotAllIsNull = 0;
    		if(!twoOrMoreIsNull)
    			nbNotTwoIsNull ++;
    		else
    			nbNotTwoIsNull = 0;
    		nb ++;
    	}
    	Vector normal = new Vector(dx,dy,dz);
    	normal = new Vector(dx/normal.length(),
    			dy/normal.length(),dz/normal.length());
        normal = getTransform(true).transformAsNormal(normal).
        	multiply(getSize());
/*      Vector spNormal = point.subtract(getLocation());    	
    	Vector rpoint = point.subtract(getLocation());        
        Vector normal2 = point.subtract(super.getLocation());
        normal2 = super.getTransform(true).transformAsNormal(normal2);
System.out.println(normal2.length()/normal.length());*/        
        return normal; 
    }
    
    public boolean intersect(Ray ray) {
//        return intersectAt(ray).getIntersections().length!=0;
    	if(!sphereIntersect(ray))
    		return false;
    	Intersection intersection = intersectAt(ray); 
    	for(int i = 0;i<intersection.getIntersections().length;i++)
    		if(intersection.getIntersections()[i]!=Double.NaN)
    			return true;
    	return false;
    }
    
    public Intersection intersectAt(Ray ray) {
    	if(!sphereIntersect(ray)) {
    		return new Intersection(ray, this, new double[0]);
    	}
   
    	Intersection si = sphereIntersectAt(ray);
        Ray newRay = new Ray(
                super.getTransform(true).getInverse().transformAsLocation(ray.getOrigin()),
                super.getTransform(true).getInverse().transformAsOffset(ray.getDirection()));

        int nb = 0;
        double tsearch[] = new double[qMaxIntersect];        
        for(int i=0;i<qMaxIntersect;i++)
        	tsearch[i] = Double.NaN;
        
        double addN = 1/qPrecision/newRay.getDirection().length();
        //double addN = .3/qPrecision/newRay.getDirection().length();
        
        Vector addVect = newRay.getDirection().multiply(addN);
     
        double start = Math.max(0, si.getIntersections()[1])-addN;
        double end = (Math.max(0, si.getIntersections()[0])+addN);

        double n = start;
        Vector pos = newRay.getOrigin().add(newRay.getDirection().multiply(start));
        double v = getValue(pos);

        while(n<=end && nb<qMaxIntersect) {
        	double oldV = v, oldN = n;
        	Vector oldPos = pos;
        	pos = pos.add(addVect);
        	n += addN;
        	v = getValue(pos);        	
        	if((oldV<0) ^ (v<0)) { //yay! a xor!
        		//should be monotoneous        	
        		double nPrc = .5;
        		double prc = .5;
            	double minPositPrc = 0;
            	double minPositV = oldV;
            	if(oldV<0) {
            		minPositPrc = 1;
                	minPositV = v;
            	}
               	for(int i=0;i<1+qPrecision*.2;i++) {            	
         			nPrc *= .5;
        			double vi = getValue(getIntermediateVect(oldPos, pos, prc));
        			if(vi>=0 && vi<=minPositV) {
        				minPositPrc = prc;
        				minPositV = vi;
        			}
        			if((oldV<0) ^ (vi<0))
        				prc -= nPrc;
        			else
        				prc += nPrc;
        		}
        		//tsearch[nb] = getIntermediateZero(
        		//	new PosValue(n1,v),new PosValue(n2,nv));
        		tsearch[nb] = oldN+minPositPrc*addN;
        		nb++;
        	}
        }
        double t[] = new double[nb];
        for(int i=0;i<nb;i++)
        	t[nb-i-1] = tsearch[i];
        return new Intersection(ray, this, t);
        //return si;
    }
    protected Vector getIntermediateVect(Vector p1, Vector p2, double prc) {
    	return p2.subtract(p1).multiply(prc).add(p1);
    }
    
    protected abstract double getValue(Vector point);


	public int getQMaxIntersect() {
		return qMaxIntersect;
	}	
    /**
     * @param maxIntersect max intersections foundable for a ray ; 
     * default to 256 (expand freely)
     */
	public void setQMaxIntersect(int maxIntersect) {
		qMaxIntersect = maxIntersect;
	}

	public double getQMaxRelativeSize() {
		return qMaxRelativeSize;
	}
	/**
	 * @param maxRelativeSize limit relative to size ; default 2, 1 
	 * for a well prepared surface, can be bigger (if you have spikes
	 * on your surface)  
	 */
	public void setQMaxRelativeSize(double maxRelativeSize) {
		qMaxRelativeSize = maxRelativeSize;
	}

	public double getQPrecision() {
		return qPrecision;
	}
	/**
	 * @param precision step or intersection search (rel. to size) ; 
	 * increase for better limits ; default to 100 - 1..1000 
	 */
	public void setQPrecision(double precision) {		
		qPrecision = precision;
		calculateQNormalPrecision();
	}


	public boolean isDiscreteShadeMode() {
		return discreteShadeMode;
	}

	/**
	 * 
	 * @param discreteShadeMode true for discrete shading mode
	 * set to true if your the gradients of your surface follow
	 * a discrete function (e.g. : the space inside your surface
	 * has the value -1 and outside 1, without gradient) 
	 * Default : false
	 */
	public void setDiscreteShadeMode(boolean discreteShadeMode) {
		this.discreteShadeMode = discreteShadeMode;
	}

	/**
	 * 
	 * @param ray entering ray
	 * @return do intersect
	 */
	private boolean sphereIntersect(Ray ray) {
        Ray newRay = new Ray(
                super.getTransform(true).getInverse().
                	transformAsLocation(ray.getOrigin())/*.
                	divide(qMaxRelativeSize)*/,
                super.getTransform(true).getInverse().
                	transformAsOffset(ray.getDirection()));
       
        Vector a = newRay.getOrigin();
        Vector d = newRay.getDirection();
        double b = d.dotProduct(a);
        double c = a.dotProduct(a);
        
        double discriminant = (b * b) - (d.dotProduct(d)) * (c - qMaxRelativeSize);
        
        return discriminant >= 0;
    }

    public Intersection sphereIntersectAt(Ray ray) {
        Ray newRay = new Ray(
                super.getTransform(true).getInverse().
                	transformAsLocation(ray.getOrigin())/*.
                	divide(qMaxRelativeSize)*/,
                super.getTransform(true).getInverse().
                	transformAsOffset(ray.getDirection()));
        
        Vector a = newRay.getOrigin();
        Vector d = newRay.getDirection();
        double b = d.dotProduct(a);
        double c = a.dotProduct(a);
        
        double discriminant = (b * b) - (d.dotProduct(d)) * (c - qMaxRelativeSize);
        double discriminantSqrt = Math.sqrt(discriminant);
        
        double t[] = new double[2];
        
        t[0] = (-b + discriminantSqrt) / (d.dotProduct(d));
        t[1] = (-b - discriminantSqrt) / (d.dotProduct(d));
        return new Intersection(ray, this, t);
    }	

    //todo: verify this implementation
    public double getIndexOfRefraction(Vector p) {
//		if (getValue(p.subtract(getLocation()))<0) {
        if(getValue(super.getTransform(true).getInverse().
        		transformAsLocation(p))<0)
			return super.getIndexOfRefraction();
		else
			return 1.0;
	}    
    
    /**
     * calculate a special precision quality index used
     * by normal : higher if some reflection/refraction are
     * needed
     */
    private void calculateQNormalPrecision() {    	
    	double maxReflectivity = 0;
    	Iterator it = getShaderSet().iterator();
    	while(it.hasNext()) {
    		Shader s = (Shader)it.next();
    		if(s instanceof ReflectionShader)
    			maxReflectivity = Math.max(maxReflectivity, 
    					((ReflectionShader)s).getReflectivity());
    		if(s instanceof RefractionShader)
    			maxReflectivity = Math.max(maxReflectivity, 
    					((ReflectionShader)s).getReflectivity());    		
    	}   		
    	qNormalPrecision = qPrecision * (1+maxReflectivity*2); 
    }
/*    
    private class PosValue {
    	public double pos, value;
    	public PosValue(double p,double v){ pos = p; value = v;}
    }
*/
    /**
     * 
     * @param ray : NORMALIZED ray
     * @return 
     */
    /*private double getNormalV(Ray ray) {
    	return ray.getOrigin().
    		dotProduct(ray.getDirection())/
    		ray.getDirection().length();
    }
    */
    /**
     * 
     * @param ray
     * @param pv : 4 pos+values
     * @return
     */
    /*private PosValue[] predicMinAndMax(Ray ray, PosValue[] pv) {
    	//v=f(p) 
    	//f(x) = a*x3+b*x2+cx+d cubic interp. from 4PVs
    	
    	//f'(x) = (a/3)*x2+(b/2)*x+c
    	
    	//3 cases : 2, 1 or 0 solutions for f'(x)=0, x between points 0 & 1
    	
    	//if 2 : return min then max
    	
    	//else : return 0 or 1 PV
    }*/
/*    
    private double getIntermediateZero(PosValue v1, PosValue v2) {
//    	y = ax+b
    	double a = (v2.value-v1.value)/(v2.pos-v1.pos);
		double b = v2.value - v2.pos*a;
		return -b/a;
    }
*/
    
}
