/*******************************************************************************
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * 
 * This library 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 Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 * USA
 * 
 * Contact Info:
 * 	Bruce Donald
 * 	Duke University
 * 	Department of Computer Science
 * 	Levine Science Research Center (LSRC)
 * 	Durham
 * 	NC 27708-0129 
 * 	USA
 * 	brd@cs.duke.edu
 * 
 * Copyright (C) 2011 Jeffrey W. Martin and Bruce R. Donald
 * 
 * <signature of Bruce Donald>, April 2011
 * Bruce Donald, Professor of Computer Science
 ******************************************************************************/


package edu.duke.donaldLab.jdshot.grid.base;

import edu.duke.donaldLab.jdshot.search.BoundContext;
import edu.duke.donaldLab.jdshot.search.SubunitType;
import edu.duke.donaldLab.share.geom.AxisAlignedBox;
import edu.duke.donaldLab.share.geom.OrientedBox;
import edu.duke.donaldLab.share.geom.Vector3;
import edu.duke.donaldLab.share.math.CompareReal;
import edu.duke.donaldLab.share.math.Distribution;
import edu.duke.donaldLab.share.math.Matrix3;
import edu.duke.donaldLab.share.protein.SubunitOrder;

public class CellCheckerBase
{
	/**************************
	 *   Functions
	 **************************/
	
	protected void getOrientedBox( OrientedBox orientedBox, BoundContext boundContext, Vector3 k )
	{
		// allocate memory
		Vector3 u = new Vector3();
		Matrix3 rot = new Matrix3();
		Vector3 d = new Vector3();
		Vector3 axis = new Vector3();
		
		// k is the point we're thinking about rotating
		double kmagSq = k.getSquaredLength();
		
		// if k is at the origin, just return a degenerate box
		if( CompareReal.eq( 0.0, kmagSq ) )
		{
			orientedBox.setDegenerate();
			return;
		}
		
		double kmag = Math.sqrt( kmagSq );
		
		// DEBUG
		//System.out.println( "k = " + k );
		//System.out.println( "kmag = " + k.getLength() );
		//System.out.println( "e = " + boundContext.epsilon );
		
		// calculate u - a rotated version of k and the center of the ball
		u.set( k );
		Matrix3.getRotation( rot, boundContext.a, boundContext.alpha );
		rot.multiply( u );
		
		// Jeff: 06/18/2009 - this is faster, but won't work on the Cn side
		//assert( boundContext.alpha == Math.PI ) : "This code assumes alpha = pi";
		//double rb = 2.0 * boundContext.epsilon * ( Math.abs( k.getDot( boundContext.a ) ) + kmag );
		
		// calculate rb - the radius of the ball at u
		double mcos = 1.0 - Math.cos( boundContext.alpha );
		double sin = Math.sin( boundContext.alpha );
		double rb =
			boundContext.epsilon * (
				Math.sqrt(
					sin * sin * kmagSq
					+ mcos * mcos * k.getDot( boundContext.a ) * k.getDot( boundContext.a )
				)
				+ Math.abs( 1.0 - Math.cos( boundContext.alpha ) ) * kmag
			);
		double rbSq = rb * rb;
		
		// calculate rs (hehe, yeah) - the radius of the sphere at the origin
		double rs = kmag;
		double rsSq = kmagSq;
		
		// DEBUG
		//System.out.println( "rb = " + rb );
		//System.out.println( "rs = " + rs );
		
		// calculate the dimensions of the oriented box
		double xhalfwidth = 0.0;
		if( rb >= 2.0 * rs )
		{
			xhalfwidth = rs;
		}
		else
		{
			xhalfwidth = rbSq / 4.0 / rs;
		}
		
		double yzhalfwidth = 0.0;
		if( rb >= rs )
		{
			yzhalfwidth = rs;
		}
		else
		{
			yzhalfwidth = rb / rs * Math.sqrt( 4.0 * rsSq - rbSq ) / 2.0;
		}
		
		// calculate d - the unit vector pointing from u to the origin (and the center of the sphere)
		d.set( u );
		d.negate();
		d.normalize();
		
		// let the x axis of the box be d
		orientedBox.x.set( d );
		
		// for the box's y axis, find a unit axis whose dot product with d is closest to 0
		// then the box's y axis is just the cross product of that unit axis and the box's x axis
		Vector3.getUnitX( axis );
		double xdot = Math.abs( d.getDot( axis ) );
		Vector3.getUnitY( axis );
		double ydot = Math.abs( d.getDot( axis ) );
		Vector3.getUnitZ( axis );
		double zdot = Math.abs( d.getDot( axis ) );
		if( xdot <= ydot && xdot <= zdot )
		{
			Vector3.getUnitX( axis );
			orientedBox.x.getCross( orientedBox.y, axis );
		}
		else if( ydot <= xdot && ydot <= zdot )
		{
			Vector3.getUnitY( axis );
			orientedBox.x.getCross( orientedBox.y, axis );
		}
		else
		{
			Vector3.getUnitZ( axis );
			orientedBox.x.getCross( orientedBox.y, axis );
		}
		
		// the box's z axis is just x cross y
		orientedBox.x.getCross( orientedBox.z, orientedBox.y );
		
		// these vectors should already be normalized, but precision screws us up
		orientedBox.x.normalize();
		orientedBox.y.normalize();
		orientedBox.z.normalize();
		
		// now, get the min and max points for the box
		d.scale( xhalfwidth );
		d.add( u );
		orientedBox.getPointInBoxSpace( d );
		orientedBox.min.set(
			d.x - xhalfwidth,
			d.y - yzhalfwidth,
			d.z - yzhalfwidth
		);
		orientedBox.max.set(
			d.x + xhalfwidth,
			d.y + yzhalfwidth,
			d.z + yzhalfwidth
		);
	}
	
	protected void getSmallAxisAlignedBox( AxisAlignedBox axisAlignedBox, BoundContext boundContext, Vector3 k )
	{
		// allocate memory
		Vector3 corner = new Vector3();
		Matrix3 boxToWorld = new Matrix3();
		OrientedBox orientedBox = new OrientedBox();
		
		getOrientedBox( orientedBox, boundContext, k );
		
		// start with an aab that has corners at infinity
		axisAlignedBox.min.set( Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY );
		axisAlignedBox.max.set( Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY );
		
		// transform oriented box corners into world space
		orientedBox.getBoxToWorldTransform( boxToWorld );
		for( int i=0; i<OrientedBox.NumCorners; i++ )
		{
			orientedBox.getCornerInBoxSpace( corner, i );
			boxToWorld.multiply( corner );
			
			if( corner.x < axisAlignedBox.min.x )
			{
				axisAlignedBox.min.x = corner.x;
			}
			if( corner.x > axisAlignedBox.max.x )
			{
				axisAlignedBox.max.x = corner.x;
			}
			if( corner.y < axisAlignedBox.min.y )
			{
				axisAlignedBox.min.y = corner.y;
			}
			if( corner.y > axisAlignedBox.max.y )
			{
				axisAlignedBox.max.y = corner.y;
			}
			if( corner.z < axisAlignedBox.min.z )
			{
				axisAlignedBox.min.z = corner.z;
			}
			if( corner.z > axisAlignedBox.max.z )
			{
				axisAlignedBox.max.z = corner.z;
			}
		}
	}
	
	protected void getBoundingBox( AxisAlignedBox bound, BoundContext boundContext )
	{
		// allocate memory
		Distribution diagonalSq = new Distribution();
		Vector3 k = new Vector3();
		AxisAlignedBox axisAlignedBox = new AxisAlignedBox();
		
		// set the corners to infinity
		bound.min.set( Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY );
		bound.max.set( Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY );
		
		// for each corner in our bound context...
		for( Vector3 t : boundContext.corners )
		{
			// calculate k by subtracting t from q
			k.set( boundContext.q );
			k.subtract( t );
			
			// get the aabb
			getSmallAxisAlignedBox( axisAlignedBox, boundContext, k );
			
			// run statistics on the squared diagonal length
			diagonalSq.add( axisAlignedBox.getDiagonalSquared() );
			
			// translate the aabb by the corner
			axisAlignedBox.min.add( t );
			axisAlignedBox.max.add( t );
			
			// DEBUG
			//System.out.println( "\t" + axisAlignedBox );
			
			// update our bound
			if( axisAlignedBox.min.x < bound.min.x )
			{
				bound.min.x = axisAlignedBox.min.x;
			}
			if( axisAlignedBox.max.x > bound.max.x )
			{
				bound.max.x = axisAlignedBox.max.x;
			}
			if( axisAlignedBox.min.y < bound.min.y )
			{
				bound.min.y = axisAlignedBox.min.y;
			}
			if( axisAlignedBox.max.y > bound.max.y )
			{
				bound.max.y = axisAlignedBox.max.y;
			}
			if( axisAlignedBox.min.z < bound.min.z )
			{
				bound.min.z = axisAlignedBox.min.z;
			}
			if( axisAlignedBox.max.z > bound.max.z )
			{
				bound.max.z = axisAlignedBox.max.z;
			}
		}
		
		boundContext.avgDiagonalSq = diagonalSq.getMean();
	}
	
	// UNDONE: these functions only work for D_2 proteins!
	// UNDONE: these functions no longer work at all! We need to rethink how subunit mapping works for Dn proteins
	protected int getSubunitId( SubunitType type, SubunitOrder subunitOrder )
	{
		// TEMP
		throw new UnsupportedOperationException();
		//return subunitOrder.mapComputedToReference( type.ordinal() );
	}
	
	protected SubunitType getSubunitType( int subunitId, SubunitOrder subunitOrder )
	{
		// TEMP
		throw new UnsupportedOperationException();
		//return SubunitType.values()[subunitOrder.mapComputedToReference( subunitId )];
	}
}
