/*
	ConRec.as
	
	Created on 05 march 2011
	
	Copyright ( c ) 2011 by Jochem van der Spek
 	
	This software is copyrighted by Jochem van der Spek. This code is base on the work of
	Paul D. Bourke CONREC.F routine and Nicholas Yue's C++ implementation
	
	The authors hereby grant permission to use, copy, and distribute this
	software and its documentation for any purpose, provided that existing
	copyright notices are retained in all copies and that this notice is included
	verbatim in any distributions. Additionally, the authors grant permission to
	modify this software and its documentation for any purpose, provided that
	such modifications are not distributed without the explicit consent of the
	authors and that existing copyright notices are retained in all copies. Some
	of the algorithms implemented by this software are patented, observe all
	applicable patent law.
	
	IN NO EVENT SHALL THE AUTHORS OR DISTRIBUTORS BE LIABLE TO ANY PARTY FOR
	DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT
	OF THE USE OF THIS SOFTWARE, ITS DOCUMENTATION, OR ANY DERIVATIVES THEREOF,
	EVEN IF THE AUTHORS HAVE BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
	
	THE AUTHORS AND DISTRIBUTORS SPECIFICALLY DISCLAIM ANY WARRANTIES, INCLUDING,
	BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
	PARTICULAR PURPOSE, AND NON - INFRINGEMENT.  THIS SOFTWARE IS PROVIDED ON AN
	"AS IS" BASIS, AND THE AUTHORS AND DISTRIBUTORS HAVE NO OBLIGATION TO PROVIDE
	MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.

	Ported to Actionscript3 from the C++  code by Nicholas Yue ( see above copyright notice ).
	@see http://paulbourke.net/papers/conrec for full description
	of code and original C++  source.
	
	@author  Jochem van der Spek http://www.jvanderspek.com
	@version 1.0 

	ConRec is a contouring subroutine for rectangularily spaced data 
		
	It emits calls to a line drawing subroutine supplied by the user
	which draws a contour map corresponding to real*4data on a randomly
	spaced rectangular grid. The coordinates emitted are in the same
	units given in the x(  ) and y(  ) arrays.

	Any number of contour levels may be specified but they must be in order of increasing value.
	@param d   -  matrix of data to contour
	@param ilb,iub,jlb,jub  -  index bounds of data matrix
				The following two, one dimensional arrays ( x and y ) contain the horizontal and
	            vertical coordinates of each sample points.
	@param x   -  data matrix column coordinates
	@param y   -  data matrix row coordinates
	@param nc  -  number of contour levels
	@param z   -  contour levels in increasing order.

	here's an example usage routine, the parameters for the constructor are kept identical
	to the original C++ from which I derived my class, which is a bit awkward, but makes 
	for clear code to compare to the original code.
	
	function createHeightMap( ):Sprite {
		var heightMap:Vector.< Vector.< Number > > = new Vector.< Vector.< Number > >();
		var xLower:int = 0;
		var xUpper:int = 10;
		var yLower:int = 0;
		var yUpper:int = 10;
		var xCoords:Vector.< Number > = new Vector.< Number >();
		var yCoords:Vector.< Number > = new Vector.< Number >();
		var numZCoords:int = 5;
		var zCoords:Vector.< Number > = new Vector.< Number >();
		for( var i:int = 0; i <= xUpper; i++ ){
			xCoords[ i ] = i * 10.0;
			yCoords[ i ] = i * 10.0;
			heightMap[ i ] = new Vector.<Number>();
			for( var j:int = 0; j <= yUpper; j++ ){
				heightMap[ i ][ j ] = Math.random();
			}
		}
		for( i = 0; i < numZCoords; i++ ){
			zCoords[ i ] = 0.1 * i;
		}
	
		return( new ConRec( heightMap, xLower, xUpper, yLower, yUpper, xCoords, yCoords, numZCoords, zCoords ) );
	} 

*/

//! use default package
package {
	import flash.display.Sprite;

	//! extend sprite, so we can draw the contours in the graphics layer
	public class ConRec extends Sprite {

		//! create contained vectors for the xsect and ysect functions
		private var h:Vector.< Number > = new Vector.< Number >( );
		private var sh:Vector.< int > = new Vector.< int >( );
		private var xh:Vector.< Number > = new Vector.< Number >( );
		private var yh:Vector.< Number > = new Vector.< Number >( );

		public function ConRec( d:Vector.< Vector.< Number > >, ilb:int, iub:int, jlb:int, jub:int, xv:Vector.< Number >, yv:Vector.< Number >, nc:int, zv:Vector.< Number >  ):void {
			// The indexing of im and jm should be noted as it has to start from zero
			// unlike the fortran counter part
			var 		im:Vector.< int > = new Vector.< int >();
			var 		jm:Vector.< int > = new Vector.< int >();
			// Note that castab is arranged differently from the FORTRAN code because
			// Fortran and C/C +  +  arrays are transposed of each other, in this case
			// it is more tricky as castab is in 3 dimension
			var 		castab:Vector.< Vector.< Vector.< int > > > = new Vector.< Vector.< Vector.< int > > >();
			var         m1:int;
			var         m2:int;
			var         m3:int;
			var         case_value:int;
			var			dmin:Number;
			var			dmax:Number;
			var      	x1:Number = 0.0;
			var      	x2:Number = 0.0;
			var      	y1:Number = 0.0;
			var      	y2:Number = 0.0;
			var 		i:int;
			var 		j:int;
			var 		k:int;
			var 		m:int;

			//! don't know how to initialize a vector 
			//! ( new Vector.<int>( [1,2,3,4] ) doesn't seem to work )
			im[ 0 ] = 0;
			im[ 1 ] = 1;
			im[ 2 ] = 1;
			im[ 3 ] = 0;

			jm[ 0 ] = 0;
			jm[ 1 ] = 0;
			jm[ 2 ] = 1;
			jm[ 3 ] = 1;

			h[ 0 ] = 0.0;
			h[ 1 ] = 0.0;
			h[ 2 ] = 0.0;
			h[ 3 ] = 0.0;
			h[ 4 ] = 0.0;

			sh[ 0 ] = 0;
			sh[ 1 ] = 0;
			sh[ 2 ] = 0;
			sh[ 3 ] = 0;
			sh[ 4 ] = 0;

			xh[ 0 ] = 0.0;
			xh[ 1 ] = 0.0;
			xh[ 2 ] = 0.0;
			xh[ 3 ] = 0.0;
			xh[ 4 ] = 0.0;

			yh[ 0 ] = 0.0;
			yh[ 1 ] = 0.0;
			yh[ 2 ] = 0.0;
			yh[ 3 ] = 0.0;
			yh[ 4 ] = 0.0;

			for( i = 0; i < 3; i++ ){
				castab[ i ] = new Vector.< Vector.< int > >();
				for( j = 0; j < 3; j++ ){
					castab[ i ][ j ] = new Vector.< int >( );
				}
			}
			castab[ 0 ][ 0 ][ 0 ] = 0;
			castab[ 0 ][ 0 ][ 1 ] = 0;
			castab[ 0 ][ 0 ][ 2 ] = 8;
			castab[ 0 ][ 1 ][ 0 ] = 0;
			castab[ 0 ][ 1 ][ 1 ] = 2;
			castab[ 0 ][ 1 ][ 2 ] = 5;
			castab[ 0 ][ 2 ][ 0 ] = 7;
			castab[ 0 ][ 2 ][ 1 ] = 6;
			castab[ 0 ][ 2 ][ 2 ] = 9;

			castab[ 1 ][ 0 ][ 0 ] = 0;
			castab[ 1 ][ 0 ][ 1 ] = 3;
			castab[ 1 ][ 0 ][ 2 ] = 4;
			castab[ 1 ][ 1 ][ 0 ] = 1;
			castab[ 1 ][ 1 ][ 1 ] = 3;
			castab[ 1 ][ 1 ][ 2 ] = 1;
			castab[ 1 ][ 2 ][ 0 ] = 4;
			castab[ 1 ][ 2 ][ 1 ] = 3;
			castab[ 1 ][ 2 ][ 2 ] = 0;

			castab[ 2 ][ 0 ][ 0 ] = 9;
			castab[ 2 ][ 0 ][ 1 ] = 6;
			castab[ 2 ][ 0 ][ 2 ] = 7;
			castab[ 2 ][ 1 ][ 0 ] = 5;
			castab[ 2 ][ 1 ][ 1 ] = 2;
			castab[ 2 ][ 1 ][ 2 ] = 0;
			castab[ 2 ][ 2 ][ 0 ] = 8;
			castab[ 2 ][ 2 ][ 1 ] = 0;
			castab[ 2 ][ 2 ][ 2 ] = 0;

			for( j = ( jub - 1 ); j >= jlb; j-- ) {
				for( i = ilb; i <= iub - 1; i++ ) {
					var temp1:Number;
					var temp2:Number;
					temp1 = Math.min( d[ i ][ j ], d[ i ][ j + 1 ] );
					temp2 = Math.min( d[ i + 1 ][ j ], d[ i + 1 ][ j + 1 ] );
					dmin  = Math.min( temp1, temp2 );
					temp1 = Math.max( d[ i ][ j ], d[ i ][ j + 1 ] );
					temp2 = Math.max( d[ i + 1 ][ j ], d[ i + 1 ][ j + 1 ] );
					dmax  = Math.max( temp1, temp2 );
					
					if ( dmax >= zv[ 0 ] && dmin <= zv[ nc - 1 ] ) {
						for ( k = 0; k < nc; k++ ){
							if ( zv[ k ] >= dmin && zv[ k ] <= dmax ) {
								for ( m = 4; m >= 0; m--  ) {
									if( m > 0 ) {
										// The indexing of im and jm should be noted as it has to
										// start from zero
										h[ m ] = d[ i + im[ m - 1 ] ][ j + jm[ m - 1 ] ] - zv[ k ];
										xh[ m ] = xv[ i + im[ m - 1 ] ];
										yh[ m ] = yv[ j + jm[ m - 1 ] ];
									} 
									else {
										h[ 0 ] = 0.25 * ( h[ 1 ] + h[ 2 ] + h[ 3 ] + h[ 4 ] );
										xh[ 0 ]= 0.5 * ( xv[ i ] + xv[ i + 1 ] );
										yh[ 0 ]= 0.5 * ( yv[ j ] + yv[ j + 1 ] );
									}
									if ( h[ m ] > 0.0 ) {
										sh[ m ] = 1;
									} 
									else if ( h[ m ] < 0.0 ) {
										sh[ m ] = - 1;
									} 
									else{
										sh[ m ] = 0;
									}
								}
								//
								// Note: at this stage the relative heights of the corners and the
								// centre are in the h array, and the corresponding coordinates are
								// in the xh and yh arrays. The centre of the box is indexed by 0
								// and the 4 corners by 1 to 4 as shown below.
								// Each triangle is then indexed by the parameter m, and the 3
								// vertices of each triangle are indexed by parameters m1,m2,and
								// m3.
								// It is assumed that the centre of the box is always vertex 2
								// though this isimportant only when all 3 vertices lie exactly on
								// the same contour level, in which case only the side of the box
								// is drawn.
								//
								//
								//      vertex 4  +  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  +  vertex 3
								//               | \               / |
								//               |   \    m - 3    /   |
								//               |     \       /     |
								//               |       \   /       |
								//               |  m=2    X   m=2   |       the centre is vertex 0
								//               |       /   \       |
								//               |     /       \     |
								//               |   /    m=1    \   |
								//               | /               \ |
								//      vertex 1  +  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  +  vertex 2
								//
								//
								//
								//               Scan each triangle in the box
								//
								for( m = 1; m <= 4; m++ ) {
									m1 = m;
									m2 = 0;
									if( m != 4 ){
										m3 = m + 1;
									} 
									else {
										m3 = 1;
									}
									case_value = castab[ sh[ m1 ] + 1 ][ sh[ m2 ] + 1 ][ sh[ m3 ] + 1 ];
									if ( case_value != 0 ) {
										switch ( case_value ) {
											case 1: // Line between vertices 1 and 2
												x1 = xh[ m1 ];
												y1 = yh[ m1 ];
												x2 = xh[ m2 ];
												y2 = yh[ m2 ];
												break;
											case 2: // Line between vertices 2 and 3
												x1 = xh[ m2 ];
												y1 = yh[ m2 ];
												x2 = xh[ m3 ];
												y2 = yh[ m3 ];
												break;
											case 3: // Line between vertices 3 and 1
												x1 = xh[ m3 ];
												y1 = yh[ m3 ];
												x2 = xh[ m1 ];
												y2 = yh[ m1 ];
												break;
											case 4: // Line between vertex 1 and side 2 - 3
												x1 = xh[ m1 ];
												y1 = yh[ m1 ];
												x2 = xsect( m2, m3 );
												y2 = ysect( m2, m3 );
												break;
											case 5: // Line between vertex 2 and side 3 - 1
												x1 = xh[ m2 ];
												y1 = yh[ m2 ];
												x2 = xsect( m3, m1 );
												y2 = ysect( m3, m1 );
												break;
											case 6: //  Line between vertex 3 and side 1 - 2
												x1 = xh[ m3 ];
												y1 = yh[ m3 ];
												x2 = xsect( m1, m2 );
												y2 = ysect( m1, m2 );
												break;
											case 7: // Line between sides 1 - 2 and 2 - 3
												x1 = xsect( m1, m2 );
												y1 = ysect( m1, m2 );
												x2 = xsect( m2, m3 );
												y2 = ysect( m2, m3 );
												break;
											case 8: // Line between sides 2 - 3 and 3 - 1
												x1 = xsect( m2, m3 );
												y1 = ysect( m2, m3 );
												x2 = xsect( m3, m1 );
												y2 = ysect( m3, m1 );
												break;
											case 9: // Line between sides 3 - 1 and 1 - 2
												x1 = xsect( m3, m1 );
												y1 = ysect( m3, m1 );
												x2 = xsect( m1, m2 );
												y2 = ysect( m1, m2 );
												break;
											default:
												break;
										}
										// Put your processing code here and comment out the printf
										trace(  "segment:"  +  x1  +  ","  +  y1  +  ","  +  x2  +  ","  +  y2  +  ","  +  zv[ k ]  );
										this.graphics.lineStyle( 1.0, 0xFF0000 );
										this.graphics.moveTo( x1, y1 );
										this.graphics.lineTo( x2, y2 );
									}
								}
							}
						}
					}
				}
			}
		}
		
		private function xsect( p1:int, p2:int  ):Number{
			return(  ( h[ p2 ] * xh[ p1 ] - h[ p1 ] * xh[ p2 ] ) / ( h[ p2 ] - h[ p1 ] )  );
		}
		
		private function ysect( p1:int, p2:int  ):Number{
			return(  ( h[ p2 ] * yh[ p1 ] - h[ p1 ] * yh[ p2 ] ) / ( h[ p2 ] - h[ p1 ] )  );
		}
	}
}