// Aqsis
// Copyright (c) 1997 - 2001, Paul C. Gregory
//
// Contact: pgregory@aqsis.com
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public
// License as published by the Free Software Foundation; either
// version 2 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
// General Public License for more details.
//
// You should have received a copy of the GNU 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

/**
 * Copyright (C) 2006-2007  NTT DATA CORPORATION
 * 
 * Version: 1.0.0 2007/04/01
 *  
 */
package net.cellcomputing.himawari.library;

import java.io.BufferedWriter;
import java.io.FileOutputStream;
import java.io.IOException;
import java.io.OutputStreamWriter;

import net.cellcomputing.himawari.accessory.STLArray;
import net.cellcomputing.himawari.accessory.STLVector;
import net.cellcomputing.himawari.accessory.primitive.p_String;
import net.cellcomputing.himawari.accessory.primitive.p_float;
import net.cellcomputing.himawari.accessory.primitive.p_int;
import net.cellcomputing.himawari.library.types.CqColor;
import net.cellcomputing.himawari.library.types.CqMatrix;
import net.cellcomputing.himawari.library.types.CqVector2D;
import net.cellcomputing.himawari.library.types.CqVector3D;
import net.cellcomputing.himawari.library.types.CqVector4D;
import static net.cellcomputing.himawari.library.RiGlobal.QGetRenderContext;
import static net.cellcomputing.himawari.library.RiGlobal.USES;
import static  net.cellcomputing.himawari.library.RiGlobal.BilinearEvaluate;
import static net.cellcomputing.himawari.library.EqIntIndex.*;
import static net.cellcomputing.himawari.library.EqEnvVars.*;
import static net.cellcomputing.himawari.library.CqStats.*;
import static net.cellcomputing.himawari.library.EqVariableClass.*;
import static net.cellcomputing.himawari.library.EqVariableType.*;
import static net.cellcomputing.himawari.library.EqSplitDir.*;
import static net.cellcomputing.himawari.library.types.PublicFunctions.*;
import static java.lang.Math.*;
import static java.lang.Float.*;

/**
 * RenderMan NURBS surface.
 * Renderman NURBS サーフェース。
 * 
 * @author NTT DATA Corporation
 */
public strictfp class CqSurfaceNURBS extends CqSurface {
	
	public static final int  NURBS_H_INCLUDED = 1;
	
//	nttdata 高速化
	//protected  STLVector<p_float> m_auKnots = new STLVector<p_float>( p_float.class );	///< Knot vector for the u direction.
	//protected  STLVector<p_float> m_avKnots = new STLVector<p_float>( p_float.class );	///< Knot vector for the v direction.
	protected  STLArray<p_float> m_auKnots = new STLArray<p_float>( p_float.class );	///< Knot vector for the u direction.
	protected  STLArray<p_float> m_avKnots = new STLArray<p_float>( p_float.class );	///< Knot vector for the v direction.
	
	protected  int	 m_uOrder;	///< Surface order in the u direction. ( Unsigned int => int )
	protected  int	 m_vOrder;	///< Surface order in the v direction. ( Unsigned int => int )
	protected  int   m_cuVerts;	///< Control point count in the u direction. ( Unsigned int => int )
	protected  int   m_cvVerts;	///< Control point count in the v direction. ( Unsigned int => int )
	protected  float m_umin;	///< Minimum value of u over surface.
	protected  float m_umax;	///< Maximum value of u over surface.
	protected  float m_vmin;	///< Minimum value of v over surface.
	protected  float m_vmax;	///< Maximum value of v over surface.
	protected  boolean	m_fPatchMesh;	///< Flag indicating this is an unsubdivided mesh.
	protected  CqTrimLoopArray	m_TrimLoops = new CqTrimLoopArray();	///< Local trim curves, prepared for this surface.

	
	
	/**
	 * デフォルトコンストラクタ。<BR>初期化を行う。
	 */
	@SuppressWarnings("unchecked")
	public CqSurfaceNURBS() {
		
		super() ;
		m_uOrder = 0 ; 
		m_vOrder = 0 ;
		m_cuVerts = 0 ;
		m_cvVerts = 0 ;
		m_umin = 0.0f ;
		m_umax = 1.0f ;
		m_vmin = 0.0f ;
		m_vmax = 1.0f ;
		m_fPatchMesh = false ;
		
		m_TrimLoops.aLoops().clear();
		//nttdata change
		//m_TrimLoops.aLoops().addAll( (STLVector<CqTrimLoop>)((CqAttributes) pAttributes()).TrimLoops().aLoops());
		for(CqTrimLoop tmp :(STLVector<CqTrimLoop>)((CqAttributes) pAttributes()).TrimLoops().aLoops())
		{
			m_TrimLoops.aLoops().add( new CqTrimLoop(tmp) );
		}
		STATS_INC( GPR_nurbs );
	}
	
	
	/**
	 * コピーコンストラクタ
	 * 
	 * @param From	コピー元
	 */
	public CqSurfaceNURBS( final CqSurfaceNURBS From )
	{
		this.assignment( From );
		STATS_INC( GPR_nurbs );
	}
	
	public 	void destruct()
	{
		//実装なし　
	}
	
	public  void Torus()
	{
		//	実装なし　
	}
	
	/**
	 * Get the order of the NURBS surface in the u direction.
	 * 
	 * @return	m_uOrder
	 */
	public int	uOrder()
	{
		return ( m_uOrder );
	}
	
	/**
	 * Get the order of the NURBS surface in the v direction.
	 * 
	 * @return	m_vOrder
	 */
	
	public int	vOrder() 
	{
		return ( m_vOrder );
	}
	
	/**
	 * Get the degree of the NURBS surface in the u direction.
	 * 
	 * @return
	 */
	public int	uDegree() 
	{
		return ( m_uOrder -1 );
	}
	
	/**
	 * Get the degree of the NURBS surface in the v direction.
	 * 
	 * @return
	 */
	public int	vDegree() 
	{
		return ( m_vOrder -1 );
	}
	
	/**
	 * Get the number of control points in the u direction.
	 * 
	 * @return	m_cuVerts
	 */
	public int	cuVerts() 
	{
		return ( m_cuVerts );
	}
	
	/**
	 * 
	 * Get the number of control points in the v direction.
	 * @return	m_cvVerts
	 */
	public int	cvVerts() 
	{
		return ( m_cvVerts );
	}
	
	/**
	 * Get the length of the knot vector for the u direction.
	 * 
	 * @return
	 */
	public int	cuKnots() 
	{
		return ( m_cuVerts + m_uOrder );
	}
	
	/**
	 * Get the length of the knot vector for the v direction.
	 * 
	 * @return
	 */
	public int	cvKnots() 
	{
		return ( m_cvVerts + m_vOrder );
	}
	
	/**
	 * Get the minimum u value of the surface.
	 * 
	 * @return	m_umin
	 */
	public float umin() 
	{
		return ( m_umin );
	}
	
	/**
	 * Set the minimum u value of the surface.
	 * 
	 * @param umin
	 */
	public void	Setumin( float umin )
	{
		m_umin = umin;
	}
	
	/**
	 * Get the minimum v value of the surface.
	 * 
	 * @return	m_vmin
	 */
	public float vmin() 
	{
		return ( m_vmin );
	}
	
	/**
	 * Set the minimum v value of the surface.
	 * 
	 * @param vmin
	 */
	public void	Setvmin( float vmin )
	{
		m_vmin = vmin;
	}
	
	/**
	 * Get the maximum u value of the surface.
	 * 
	 * @return
	 */
	public float umax() 
	{
		return ( m_umax );
	}
	
	/** 
	 * Set the maximum u value of the surface.
	 *
	 * @param	umax
	 */
	public void	Setumax( float umax )
	{
		m_umax = umax;
	}
	
	/**
	 * Get the maximum v value of the surface.
	 * 
	 * @return	m_vmax
	 */
	public float vmax() 
	{
		return ( m_vmax );
	}
	
	/**
	 * Set the maximum v value of the surface.
	 * 
	 * @param vmax
	 */
	public void	Setvmax( float vmax )
	{
		m_vmax = vmax;
	}
	
	/**
	 * Get a reference to the knot vector for the u direction.
	 * 
	 * @return	m_auKnots
	 */
	//nttdata 高速化
//	public STLVector<p_float> auKnots()
//	{
//		return ( m_auKnots );
//	}
//	
//	/**
//	 * Get a reference to the knot vector for the v direction.
//	 * 
//	 * @return	m_avKnots
//	 */
//	public STLVector<p_float> avKnots()
//	{
//		return ( m_avKnots );
//	}
	public STLArray<p_float> auKnots()
	{
		return ( m_auKnots );
	}
	
	public STLArray<p_float> avKnots()
	{
		return ( m_avKnots );
	}
	
	/**
	 * Determine how many segments in u for this surface patch.
	 * 
	 * @return
	 */
	public int cuSegments() 
	{
		return ( 1 + m_cuVerts - m_uOrder );
	}
	
	/**
	 * Determine how many segments in v for this surface patch.
	 * 
	 * @return
	 */
	public int cvSegments() 
	{
		return ( 1 + m_cvVerts - m_vOrder );
	}
	
	/**
	 * Determine whether this surface is a mesh, <BR>
	 * and needs to be split into segments before continuing.
	 * 
	 * @return	m_fPatchMesh
	 */
	public boolean fPatchMesh() 
	{
		return ( m_fPatchMesh );
	}
	
	/**
	 * Mark this mesh as being part of a mesh or not.
	 * 同メソッドをオーバーロード<BR>
	 * 初期値　:　true　を与える。
	 */
	public void SetfPatchMesh()
	{
		SetfPatchMesh( true );
	}
	
	/**
	 * Mark this mesh as being part of a mesh or not.
	 * 
	 * @param fPatchMesh
	 */
	public void SetfPatchMesh( boolean fPatchMesh )
	{
		m_fPatchMesh = fPatchMesh;
	}
	
	/**
	 * 代入演算子のオーバーロード
	 * 
	 * 
	 * @param From	代入元
	 */
	@SuppressWarnings("unchecked")
	public void assignment( final CqSurfaceNURBS From )
	{
		// Use the CqSurface assignment operator.
		super.assignment( From );
		
		// Initialise the NURBS surface.
		Init( From.m_uOrder, From.m_vOrder, From.m_cuVerts, From.m_cvVerts );
		
		m_umin = From.m_umin;
		m_umax = From.m_umax;
		m_vmin = From.m_vmin;
		m_vmax = From.m_vmax;
		
		m_fPatchMesh = From.m_fPatchMesh;
		
		int i;
		for ( i = From.m_auKnots.size() - 1; i >= 0; i-- )
			m_auKnots.get( i ).value = From.m_auKnots.get( i ).value;
		
		for ( i = From.m_avKnots.size() - 1; i >= 0; i-- )
			m_avKnots.get( i ).value = From.m_avKnots.get( i ).value;
		
		m_TrimLoops.aLoops().clear();
		m_TrimLoops.aLoops().addAll( (STLVector<CqTrimLoop>)((CqAttributes) pAttributes()).TrimLoops().aLoops() );
	}
	
	/**
	 * 比較演算子のオーバーロード
	 * 
	 * 
	 * @return	比較結果
	 */
	public int compareTo( final CqSurfaceNURBS from )
	{
		if ( from.m_cuVerts != m_cuVerts || from.m_cvVerts != m_cvVerts )
			return ( 0 );
		
		if ( from.m_uOrder != m_uOrder || from.m_vOrder != m_vOrder )
			return ( 0 );
		
		int i;
		for ( i = P().Size() - 1; i >= 0; i-- ) {
			if ( !((CqVector4D)P().pValue_get( i , 0 )).equals((CqVector4D)from.P().pValue_get( i , 0 ) ))
				return ( 0 );
		}
		
		for ( i = m_auKnots.size() - 1; i >= 0; i-- ){
//			if ( m_auKnots[ i ] != from.m_auKnots[ i ] )
			if ( m_auKnots.get( i ).value != from.m_auKnots.get( i ).value )
				return ( 0 );
		}
		
		for ( i = m_avKnots.size() - 1; i >= 0; i-- ){
//			if ( m_avKnots[ i ] != from.m_avKnots[ i ] )
			if ( m_avKnots.get( i ).value != from.m_avKnots.get( i ).value )
				return ( 0 );
		}
		return ( 1 );
	}
	
	/** 
	 * Get the control point at the specified u,v index.
	 * @param u Index in the u direction.
	 * @param v Index in the v direction.
	 * @return Reference to the 4D homogenous control point.
	 */
	public final CqVector4D CP( final int u, final int v )
	{
		return ( CqVector4D )( P().pValue_get(( v * m_cuVerts )+ u , 0 ) );
	}
	
//	
//	/**
//	*  Get the control point at the specified u,v index.
//	* @param u Index in the u direction.
//	* @param v Index in the v direction.
//	* @return Reference to the 4D homogenous control point.
//	*/
//	public final CqVector4D CP( final int u, int v )
//	{
//	return ( CqVector4D )( P().pValue_get( ( v * m_cuVerts ) + u )[0] );
//	}
	
	/** 
	 * Initialise the NURBS structures to take a NURBS surfafe of the specified dimensions.
	 * @param uOrder	 The required order in the u direction.
	 * @param vOrder	 The required order in the v direction.
	 * @param cuVerts	 The required control point count in the u direction.
	 * @param cvVerts	 The required control point count in the v direction.
	 */
	public void	Init( int uOrder, int vOrder, int cuVerts, int cvVerts )
	{
		int uKnots = cuVerts + uOrder;
		int vKnots = cvVerts + vOrder;
		m_auKnots.resize( uKnots );
		m_avKnots.resize( vKnots );
		m_uOrder = uOrder;
		m_vOrder = vOrder;
		m_cuVerts = cuVerts;
		m_cvVerts = cvVerts;
	}
	
	/**
	 * Find the span in the U knot vector containing 
	 * the specified parameter value.
	 * 
	 * @param u
	 * @return
	 */
	public int	FindSpanU( float u )
	{
		if ( u >= m_auKnots.get( m_cuVerts ).value )
			return ( m_cuVerts - 1 );
		if ( u <= m_auKnots.get( uDegree() ).value )
			return ( uDegree() );
		
		int low = 0;
		int high = m_cuVerts + 1;
		int mid = ( low + high ) / 2;
		
		while ( u < m_auKnots.get( mid ).value || u >= m_auKnots.get( mid + 1 ).value )
		{
			if ( u < m_auKnots.get( mid ).value )
				high = mid;
			else
				low = mid;
			mid = ( low + high ) / 2;
		}
		return ( mid );
	}
	
	/** 
	 * Find the span in the V knot vector containing 
	 * the specified parameter value.
	 * 
	 * @param v
	 * @return
	 */
	public int	FindSpanV( float v )
	{
		if ( v >= m_avKnots.get( m_cvVerts ).value )
			return ( m_cvVerts -1 );
		if ( v <= m_avKnots.get( vDegree()).value )
			return ( vDegree() );
		
		int low = 0;
		int high = m_cvVerts + 1;
		int mid = ( low + high ) / 2;
		
		while ( v < m_avKnots.get( mid ).value || v >= m_avKnots.get( mid + 1 ).value )
		{
			if ( v < m_avKnots.get( mid ).value )
				high = mid;
			else
				low = mid;
			mid = ( low + high ) / 2;
		}
		return ( mid );
	}
	
	/**
	 * Return the basis functions for the specified parameter value.
	 * 
	 * @param u
	 * @param span
	 * @param aKnots
	 * @param k
	 * @param BasisVals
	 */
	//nttdata　高速化
//	public void BasisFunctions( float u, int i, STLVector<p_float> U, int k, STLVector<p_float> N )
//	{
//		int j, r;
//		float saved, temp;
//		
//		float[] left = new float[ k ];
//		float[] right = new float[ k ];
//		
//		N.get( 0 ).value = 1.0f ;
//		for ( j = 1; j <= k - 1; j++ )
//		{
//			left[ j ] = u - U.get( i + 1 - j ).value ;
//			right[ j ] = U.get( i + j ).value - u ;
//			saved = 0.0f;
//			for ( r = 0; r < j; r++ )
//			{
//				temp = N.get( r ).value / ( right[ r + 1 ] + left[ j - r ] );
//				N.get( r ).value = saved + right[ r + 1 ] * temp ;
//				saved = left[ j - r ] * temp;
//			}
//			N.get( j ).value = saved;
//		}
//	}
	public void BasisFunctions( float u, int i, STLArray<p_float> U, int k, STLArray<p_float> N )
	{
		int j, r;
		float saved, temp;
		
		float[] left = new float[ k ];
		float[] right = new float[ k ];
		
		N.get( 0 ).value = 1.0f ;
		for ( j = 1; j <= k - 1; j++ )
		{
			left[ j ] = u - U.get( i + 1 - j ).value ;
			right[ j ] = U.get( i + j ).value - u ;
			saved = 0.0f;
			for ( r = 0; r < j; r++ )
			{
				temp = N.get( r ).value / ( right[ r + 1 ] + left[ j - r ] );
				N.get( r ).value = saved + right[ r + 1 ] * temp ;
				saved = left[ j - r ] * temp;
			}
			N.get( j ).value = saved;
		}
	}
	
	/**
	 * Return the basis functions for the specified parameter value.
	 * 
	 * 
	 * @param u
	 * @param i
	 */
	public void DersBasisFunctions( float u, int i, STLArray<p_float> U, int k, int n, STLVector<STLVector<p_float>> ders )
	{
		int j, r;
		float saved, temp;
		
		float left[] = new float[ k ];
		float right[] = new float[ k ];
		STLVector<STLVector<p_float>> ndu = new STLVector<STLVector<p_float>>( 2, p_float.class, k );
		STLVector<STLVector<p_float>> a = new STLVector<STLVector<p_float>>( 2, p_float.class, 2 );
		
		for ( j = 0; j < k; j++ ) 
			ndu.get( j ).resize( k );
		ders.resize( n + 1 );
		for ( j = 0; j < n + 1; j++ ) 
			ders.get( j ).resize( k );
		a.get( 0 ).resize( k );
		a.get( 1 ).resize( k );
		
		int p = k - 1;
		
		ndu.get( 0 ).get( 0 ).value = 1.0f ;
		for ( j = 1; j <= p; j++ )
		{
			left[ j ] = u - U.get( i + 1 - j ).value ;
			right[ j ] = U.get( i + j ).value - u ;
			saved = 0.0f;
			for ( r = 0; r < j; r++ )
			{
				ndu.get( j ).get( r ).value = right[ r + 1 ] + left[ j - r ] ;
				temp = ndu.get( r ).get( j - 1 ).value / ndu.get( j ).get( r ).value;
				
				ndu.get( r ).get( j ).value = saved + right[ r + 1 ] * temp ;
				saved = left[ j - r ] * temp;
			}
			ndu.get( j ).get( j ).value = saved ;
		}
		
		// Load the basis functions
		for ( j = 0; j <= p; j++ )
			ders.get( 0 ).get( j ).value =  ndu.get( j ).get( p ).value ;
		
		// Compute the derivatives.
		for ( r = 0; r <= p; r ++ )
		{
			// Alternate rows in array a.
			int s1 = 0;
			int s2 = 1;
			a.get( 0 ).get( 0 ).value = 1.0f ;
			int j1, j2;
			
			// Loop to compute the kth derivative
			for ( k = 1; k <= n; k++ )
			{
				float d = 0.0f;
				int rk = r - k;
				int pk = p - k;
				if ( r >= k )
				{
					a.get( s2 ).get( 0 ).value = a.get( s1 ).get( 0 ).value / ndu.get( pk + 1 ).get( rk ).value ;
					d = a.get( s2 ).get( 0 ).value * ndu.get( rk ).get( pk ).value ;
				}
				if ( rk >= -1 )
					j1 = 1;
				else
					j1 = -rk;
				
				if ( r - 1 <= pk )
					j2 = k - 1;
				else
					j2 = p - r;
				
				for ( j = j1; j <= j2; j++ )
				{
					a.get( s2 ).get( j ).value = ( a.get( s1 ).get( j ).value - a.get( s1 ).get( j - 1 ).value ) / ndu.get( pk + 1 ).get( rk + j ).value ;
					d += a.get( s2 ).get( j ).value * ndu.get( rk + j ).get( pk ).value;
				}
				
				if ( r <= pk )
				{
					a.get( s2 ).get( k ).value = -a.get( s1 ).get( k - 1 ).value / ndu.get( pk + 1 ).get( r ).value ;
					d += a.get( s2 ).get( k ).value * ndu.get( r ).get( pk ).value;
				}
				ders.get( k ).get( r ).value = d ;
				
				// Switch rows.
				j = s1;
				s1 = s2;
				s2 = j;
			}
		}
		// Multiply through the correct factors.
		r = p; 
		for ( k = 1; k <= n; k++ )
		{
			for ( j = 0; j <= p; j++ )
				ders.get( k ).get( j ).value *= r;
			r *= ( p - k );
		}
	}
	//nttdata 高速化
//	public void DersBasisFunctions( float u, int i, STLVector<p_float> U, int k, int n, STLVector<STLVector<p_float>> ders )
//	{
//		int j, r;
//		float saved, temp;
//		
//		float left[] = new float[ k ];
//		float right[] = new float[ k ];
//		STLVector<STLVector<p_float>> ndu = new STLVector<STLVector<p_float>>( 2, p_float.class, k );
//		STLVector<STLVector<p_float>> a = new STLVector<STLVector<p_float>>( 2, p_float.class, 2 );
//		
//		for ( j = 0; j < k; j++ ) 
//			ndu.get( j ).resize( k );
//		ders.resize( n + 1 );
//		for ( j = 0; j < n + 1; j++ ) 
//			ders.get( j ).resize( k );
//		a.get( 0 ).resize( k );
//		a.get( 1 ).resize( k );
//		
//		int p = k - 1;
//		
//		ndu.get( 0 ).get( 0 ).value = 1.0f ;
//		for ( j = 1; j <= p; j++ )
//		{
//			left[ j ] = u - U.get( i + 1 - j ).value ;
//			right[ j ] = U.get( i + j ).value - u ;
//			saved = 0.0f;
//			for ( r = 0; r < j; r++ )
//			{
//				ndu.get( j ).get( r ).value = right[ r + 1 ] + left[ j - r ] ;
//				temp = ndu.get( r ).get( j - 1 ).value / ndu.get( j ).get( r ).value;
//				
//				ndu.get( r ).get( j ).value = saved + right[ r + 1 ] * temp ;
//				saved = left[ j - r ] * temp;
//			}
//			ndu.get( j ).get( j ).value = saved ;
//		}
//		
//		// Load the basis functions
//		for ( j = 0; j <= p; j++ )
//			ders.get( 0 ).get( j ).value =  ndu.get( j ).get( p ).value ;
//		
//		// Compute the derivatives.
//		for ( r = 0; r <= p; r ++ )
//		{
//			// Alternate rows in array a.
//			int s1 = 0;
//			int s2 = 1;
//			a.get( 0 ).get( 0 ).value = 1.0f ;
//			int j1, j2;
//			
//			// Loop to compute the kth derivative
//			for ( k = 1; k <= n; k++ )
//			{
//				float d = 0.0f;
//				int rk = r - k;
//				int pk = p - k;
//				if ( r >= k )
//				{
//					a.get( s2 ).get( 0 ).value = a.get( s1 ).get( 0 ).value / ndu.get( pk + 1 ).get( rk ).value ;
//					d = a.get( s2 ).get( 0 ).value * ndu.get( rk ).get( pk ).value ;
//				}
//				if ( rk >= -1 )
//					j1 = 1;
//				else
//					j1 = -rk;
//				
//				if ( r - 1 <= pk )
//					j2 = k - 1;
//				else
//					j2 = p - r;
//				
//				for ( j = j1; j <= j2; j++ )
//				{
//					a.get( s2 ).get( j ).value = ( a.get( s1 ).get( j ).value - a.get( s1 ).get( j - 1 ).value ) / ndu.get( pk + 1 ).get( rk + j ).value ;
//					d += a.get( s2 ).get( j ).value * ndu.get( rk + j ).get( pk ).value;
//				}
//				
//				if ( r <= pk )
//				{
//					a.get( s2 ).get( k ).value = -a.get( s1 ).get( k - 1 ).value / ndu.get( pk + 1 ).get( r ).value ;
//					d += a.get( s2 ).get( k ).value * ndu.get( r ).get( pk ).value;
//				}
//				ders.get( k ).get( r ).value = d ;
//				
//				// Switch rows.
//				j = s1;
//				s1 = s2;
//				s2 = j;
//			}
//		}
//		// Multiply through the correct factors.
//		r = p; 
//		for ( k = 1; k <= n; k++ )
//		{
//			for ( j = 0; j <= p; j++ )
//				ders.get( k ).get( j ).value *= r;
//			r *= ( p - k );
//		}
//	}
//	
	/**
	 * Evaluate the nurbs surface at parameter values u,v.
	 * 
	 * 
	 * @param <T>
	 * @param <SLT>
	 * @param u
	 * @param v
	 * @param pParam
	 * @return
	 */ 
	@SuppressWarnings("unchecked")
	public <T, SLT> T Evaluate( float u, float v, CqParameterTyped<T, SLT> pParam )
	{
		//nttdata 高速化
//		STLVector<p_float> Nu = new STLVector<p_float>( p_float.class, m_uOrder );
//		STLVector<p_float> Nv = new STLVector<p_float>( p_float.class, m_vOrder );
		STLArray<p_float> Nu = new STLArray<p_float>( p_float.class, m_uOrder );
		STLArray<p_float> Nv = new STLArray<p_float>( p_float.class, m_vOrder );
		
		
		/* Evaluate non-uniform basis functions (and derivatives) */
		
		int uspan = FindSpanU( u );
		BasisFunctions( u, uspan, m_auKnots, m_uOrder, Nu );
		int vspan = FindSpanV( v );
		BasisFunctions( v, vspan, m_avKnots, m_vOrder, Nv );
		int uind = uspan - uDegree();
		
		T S = (T)pParam.pValue_get( 0 , 0 );
		int l, k;
		if( S instanceof p_float )
		{ 
			S = (T) new p_float(); 
			float temp = 0;
			for ( l = 0; l <= vDegree(); l++ ){
				int vind = vspan - vDegree() + l;
				for ( k = 0; k <= uDegree(); k++ )
					temp = temp + Nu.get( k ).value * ((p_float)pParam.pValue_get( vind * m_cuVerts + uind + k ,0)).value;
				((p_float)S).value =(((p_float)S).value + Nv.get( 0 ).value * temp );
			}	
		}else if( S instanceof p_int )
		{
			S = (T) new p_int();
			for ( l = 0; l <= vDegree(); l++ ){
				int temp = 0;
				int vind = vspan - vDegree() + l;
				for ( k = 0; k <= uDegree(); k++ )
					temp = temp + (int)Nu.get( k ).value * ((p_int)pParam.pValue_get(( vind * m_cuVerts ) + uind + k ,0)).value ;
				((p_int)S).value = (((p_int)S).value + (int)Nv.get( l ).value * temp );
			} 
		}else if( S instanceof CqVector3D )
		{
			S = (T) new CqVector3D();
			for ( l = 0; l <= vDegree(); l++ ){
				CqVector3D temp = new CqVector3D();
				int vind = vspan - vDegree() + l;
				for ( k = 0; k <= uDegree(); k++ )
					temp.assignment( temp.add( ((CqVector3D)pParam.pValue_get(( vind * m_cuVerts ) + uind + k ,0)).mul( Nu.get( k ).value )));
				((CqVector3D)S).assignment( ((CqVector3D)S).add( temp.mul( Nv.get( l ).value )) );
			}
		}else if( S instanceof CqVector4D )
		{
			S = (T) new CqVector4D();
			for ( l = 0; l <= vDegree(); l++ ){
				CqVector4D temp = new CqVector4D();
				int vind = vspan - vDegree() + l;
				//nttdata 高速化
//				for ( k = 0; k <= uDegree(); k++ ) {
//					temp.assignment( temp.add( ((CqVector4D)pParam.pValue_get(( vind * m_cuVerts ) + uind + k ,0)).mul( Nu.get( k ).value )));
//				}
//				    ((CqVector4D)S).assignment( ((CqVector4D)S).add( temp.mul( Nv.get( l ).value )));
				CqVector4D temp2 = new CqVector4D();
				for ( k = 0; k <= uDegree(); k++ ) {
					temp2.assignment((CqVector4D)pParam.pValue_get(( vind * m_cuVerts ) + uind + k ,0));
					temp.assignAdd( temp2.assignMul( Nu.get( k ).value ));
				}
					((CqVector4D)S).assignAdd( temp.assignMul( Nv.get( l ).value ));
						
			}
		}else if( S instanceof CqColor )
		{
			S = (T) new CqColor();
			for ( l = 0; l <= vDegree(); l++ ){
				CqColor temp = new CqColor();
				int vind = vspan - vDegree() + l;
				for ( k = 0; k <= uDegree(); k++ )
					temp.assignment( temp.add( ((CqColor)pParam.pValue_get(( vind * m_cuVerts ) + uind + k ,0)).mul( Nu.get( k ).value )));
				((CqColor)S).assignment( ((CqColor)S).add( temp.mul( Nv.get( l ).value ) ));
			}
		}else if( S instanceof CqMatrix )
		{
			S = (T) new CqMatrix();
			for ( l = 0; l <= vDegree(); l++ ){
				CqMatrix temp = new CqMatrix();
				int vind = vspan - vDegree() + l;
				for ( k = 0; k <= uDegree(); k++ )
					temp.assignment( temp.add( ((CqMatrix)pParam.pValue_get(( vind * m_cuVerts ) + uind + k ,0)).mulInv( Nu.get( k ).value )));
				((CqMatrix)S).assignment( ((CqMatrix)S).add( temp.mulInv( Nv.get( l ).value )) );
			}
		}else if( S instanceof p_String )
		{
			S = (T) new p_String();
			for ( l = 0; l <= vDegree(); l++ ){
				p_String temp = new p_String();
				int vind = vspan - vDegree() + l;
				for ( k = 0; k <= uDegree(); k++ )
					temp.value = temp.value + ((p_String)pParam.pValue_get(( vind * m_cuVerts ) + uind + k ,0)).value;
//				temp.value = temp.value + Nu.get( k ).value * ((p_String)pParam.pValue_get(( vind * m_cuVerts ) + uind + k ,0)).value )));
				((p_String)S).value = ((p_String)S).value + temp.value ;
//				((p_String)S).value = ((p_String)S).value + Nv.get( l ).value * temp.value ;
			}
		}
		return ( S );
	}
	
	/**
	 *  Evaluate the nurbs surface at parameter values u,v.
	 * 
	 * 
	 * @param u
	 * @param v
	 * @return
	 */
	public CqVector4D	EvaluateWithNormal( float u, float v, CqVector4D P )
	{
		CqVector4D N;
		int d = 1;	// Default to 1st order derivatives.
		int k, l, s, r;
		int p = uDegree();
		int q = vDegree();
		
		STLVector<STLVector<CqVector4D>> SKL = new STLVector<STLVector<CqVector4D>>( 2, CqVector4D.class, d + 1 );
		for ( k = 0; k <= d; k++ ) 
			SKL.get( k ).resize( d + 1 );
		STLVector<STLVector<p_float>> Nu = new STLVector<STLVector<p_float>>( 2, p_float.class );
		STLVector<STLVector<p_float>> Nv = new STLVector<STLVector<p_float>>( 2, p_float.class );
		STLVector<CqVector4D> temp = new STLVector<CqVector4D>( CqVector4D.class, q + 1 );
		
		int du = min( d, p );
		for ( k = p + 1; k <= d; k++ ) // 初期化はいらない気がする。nttdata
			for ( l = 0; l <= d - k; l++ )
				SKL.get( k ).set( l, new CqVector4D( 0.0f, 0.0f, 0.0f, 1.0f ));
		
		int dv = min( d, q );
		for ( l = q + 1; l <= d; l++ )
			for ( k = 0; k <= d - l; k++ )
				SKL.get( k ).set( l, new CqVector4D( 0.0f, 0.0f, 0.0f, 1.0f ));
		
		int uspan = FindSpanU( u );
		DersBasisFunctions( u, uspan, m_auKnots, m_uOrder, du, Nu );
		int vspan = FindSpanV( v );
		DersBasisFunctions( v, vspan, m_avKnots, m_vOrder, dv, Nv );
		
		for ( k = 0; k <= du; k++ )
		{
			for ( s = 0; s <= q; s++ )
			{
//				temp.set( s, new CqVector4D( 0.0f, 0.0f, 0.0f, 1.0f ));
				for ( r = 0; r <= p; r++ )
					temp.get( s ).assignment( temp.get( s ).add( CP( uspan - p + r, vspan - q + s ).mul( Nu.get( k ).get( r ).value )));
			}
			int dd = min( d - k, dv );
			for ( l = 0; l <= dd; l++ )
			{
//				SKL.get( k ).set( l, new CqVector4D( 0.0f, 0.0f, 0.0f, 1.0f ));
				for ( s = 0; s <= q; s++ )
					SKL.get( k ).get( l ).assignment( SKL.get( k ).get( l ).add( temp.get( s ).mul( Nv.get( l ).get( s ).value )));
			}
		}
		N = SKL.get( 1 ).get( 0 ).mod( SKL.get( 0 ).get( 1 ));
		N.Unit();
		
		P.assignment( SKL.get( 0 ).get( 0 ));
		
		return ( N );
	}
	
	/**
	 * Split this NURBS surface into two subsurfaces along u or v depending on dirflag (TRUE=u)
	 * 
	 * @param	nrbA
	 * @param	nrbB
	 * @param	dirflag
	 */
	public  void	SplitNURBS( CqSurfaceNURBS nrbA, CqSurfaceNURBS nrbB, boolean dirflag )
	{
		//nttdata 高速化
		//STLVector<p_float> aKnots = ( dirflag ) ? m_auKnots : m_avKnots;
		STLArray<p_float> aKnots = ( dirflag ) ? m_auKnots : m_avKnots;
		int Order = ( dirflag ) ? m_uOrder : m_vOrder;
		
		// Add a multiplicty k knot to the knot vector in the direction
		// specified by dirflag, and refine the surface.  This creates two
		// adjacent surfaces with c0 discontinuity at the seam.
		int extra = 0;
		int last = ( dirflag ) ? ( m_cuVerts + m_uOrder - 1 ) : ( m_cvVerts + m_vOrder - 1 );
		//    int middex=last/2;
		//    float midVal=aKnots[middex];
		float midVal = ( aKnots.get( 0 ).value + aKnots.get( last ).value ) / 2;
		int middex = ( dirflag ) ? FindSpanU( midVal ) : FindSpanV( midVal );
		
		// Search forward and backward to see if multiple knot is already there
		int i = 0;
		int same = 0;
		if ( aKnots.get( middex ).value == midVal )
		{
			i = middex + 1 ;
			same = 1;
			while ( ( i < last ) && ( aKnots.get( i ).value == midVal ) )
			{
				i++;
				same++;
			}
			
			i = middex - 1;
			while ( ( i > 0L ) && ( aKnots.get( i ).value == midVal ) )
			{
				i--;
				middex--;	// middex is start of multiple knot
				same++;
			}
		}
		
		if ( i <= 0 )           	    // No knot in middle, must create it
		{
			midVal = ( aKnots.get( 0 ).value + aKnots.get( last ).value ) / 2.0f ;
			middex = 0;
			while ( aKnots.get( middex + 1 ).value < midVal )
				middex++;
			same = 0;
		}
		
		extra = Order - same;
		//STLVector<p_float> anewKnots = new STLVector<p_float>( p_float.class, extra );
		STLArray<p_float> anewKnots = new STLArray<p_float>( p_float.class, extra );
		
		if ( same < Order )           	    // Must add knots
		{
			for ( i = 0; i < extra; i++ ) 
				anewKnots.get( i ).value = midVal ;
		}
		
		int SplitPoint = ( extra < Order ) ? middex - 1 : middex;
		if ( dirflag ) 
			RefineKnotU( anewKnots );
		else	
			RefineKnotV( anewKnots );
		
		// Build the two child surfaces, and copy the data from the refined
		// version of the parent (tmp) into the two children
		
		// First half
		nrbA.Init( m_uOrder, m_vOrder, ( dirflag ) ? SplitPoint + 1 : m_cuVerts, ( dirflag ) ? m_cvVerts : SplitPoint + 1 );
		int j;
		
		for ( CqParameter iUP : m_aUserParams )
		{
			if ( iUP.Class().getValue() == class_vertex )
			{
				CqParameter pNewA = iUP.CloneType( iUP.strName(), iUP.Count() );
				pNewA.SetSize( nrbA.cuVerts() * nrbA.cvVerts() );
				for ( i = 0; i < nrbA.m_cvVerts; i++ )
					for ( j = 0; j < nrbA.m_cuVerts; j++ )
						pNewA.SetValue( iUP, ( i * nrbA.cuVerts() ) + j, ( i * m_cuVerts ) + j );
				nrbA.AddPrimitiveVariable( pNewA );
			}
		}
		
		for ( i = 0; i < nrbA.m_uOrder + nrbA.m_cuVerts; i++ )
			nrbA.m_auKnots.get( i ).value = m_auKnots.get( i ).value;
		for ( i = 0; i < nrbA.m_vOrder + nrbA.m_cvVerts; i++ )
			nrbA.m_avKnots.get( i ).value = m_avKnots.get( i ).value;
		
		// Second half
		SplitPoint++;
		nrbB.Init( m_uOrder, m_vOrder, ( dirflag ) ? m_cuVerts - SplitPoint : m_cuVerts, ( dirflag ) ? m_cvVerts : m_cvVerts - SplitPoint );
		for ( CqParameter iUP : m_aUserParams )
		{
			if ( iUP.Class().getValue() == class_vertex )
			{
				CqParameter  pNewB = iUP.CloneType( iUP.strName(), iUP.Count() );
				pNewB.SetSize( nrbB.cuVerts() * nrbB.cvVerts() );
				for ( i = 0; i < nrbB.m_cvVerts; i++ )
				{
					for ( j = 0; j < nrbB.m_cuVerts; j++ )
					{
						int iSrc = ( dirflag ) ? i : i + SplitPoint;
						iSrc *= m_cuVerts;
						iSrc += ( dirflag ) ? j + SplitPoint : j;
						pNewB.SetValue( iUP, ( i * nrbB.cuVerts() + j ), iSrc );
					}
				}
				nrbB.AddPrimitiveVariable( pNewB );
			}
		}

		for ( i = 0; i < nrbB.m_uOrder + nrbB.m_cuVerts; i++ ) 
			nrbB.m_auKnots.get( i ).value = m_auKnots.get(( dirflag ) ? ( i + SplitPoint ) : i ).value ;
		for ( i = 0; i < nrbB.m_vOrder + nrbB.m_cvVerts; i++ )
			nrbB.m_avKnots.get( i ).value = m_avKnots.get(( dirflag ) ? i : ( i + SplitPoint )).value ;
		
	}
	
	/**
	 * Split the NURBS surface into B-Spline (sub) surfaces
	 * 
	 * @param	Array
	 */
	//nttdata 高速化
	//public  void SubdivideSegments( STLVector<CqSurfaceNURBS> S )
	public  void SubdivideSegments( STLArray<CqSurfaceNURBS> S )
	{
		int uSplits = cuSegments();
		int vSplits = cvSegments();
		
		// Resize the array to hold the aplit surfaces.
		S.resize( uSplits * vSplits );
		
		int iu, iv;
		
		// An array to hold the split points in u and v, fill in the first one for us.
//		STLVector<p_int> uSplitPoint = new STLVector<p_int>( p_int.class ,uSplits + 1 );
//		STLVector<p_int> vSplitPoint = new STLVector<p_int>( p_int.class ,vSplits + 1 );
		int uSplitPoint[] = new int[ uSplits + 1 ];
		int vSplitPoint[] = new int[ vSplits + 1 ];
		uSplitPoint[ 0 ] = 0 ;
		vSplitPoint[ 0 ] = 0 ;
		
		// Refine the knot vectors as appropriate to generate the required split points in u
		for ( iu = 1; iu < uSplits; iu++ )
		{
			float su = ( (float) iu / (float) uSplits )
			* ( m_auKnots.get( m_cuVerts ).value - m_auKnots.get( m_uOrder - 1 ).value )
			+ m_auKnots.get( m_uOrder - 1 ).value;
			
			int extra = 0;
			int last = m_cuVerts + m_uOrder - 1;
			float midVal = su;
			int middex = FindSpanU( midVal );
			
			// Search forward and backward to see if multiple knot is already there
			int i = 0;
			int same = 0;
			if ( auKnots().get( middex ).value == midVal )
			{
				i = middex + 1;
				same = 1;
				while ( ( i < last ) && ( auKnots().get( i ).value == midVal ) ){
					i++;
					same++;
				}
				
				i = middex - 1;
				while ( ( i > 0L ) && ( auKnots().get( i ).value == midVal ) ){
					i--;
					middex--;	// middex is start of multiple knot
					same++;
				}
			}
			
			if ( i <= 0L )           	    // No knot in middle, must create it
			{
				middex = 0;
				while ( auKnots().get( middex + 1 ).value < midVal )
					middex++;
				same = 0;
			}
			
			extra = m_uOrder - same;
			
			//nttdata 高速化
			//STLVector<p_float> anewKnots = new STLVector<p_float>( p_float.class, extra );
			STLArray<p_float> anewKnots = new STLArray<p_float>( p_float.class, extra );
			
			if ( same < m_uOrder )           	    // Must add knots
			{
				for ( i = 0; i < extra; i++ )
					anewKnots.get( i ).value =  midVal ;
			}
			
			uSplitPoint[ iu ] = ( extra < m_uOrder ) ? middex - 1 : middex;
			RefineKnotU( anewKnots );
		}
		
		// Refine the knot vectors as appropriate to generate the required split points in v
		for ( iv = 1; iv < vSplits; iv++ )
		{
			float sv = ( (float) iv / (float) vSplits )
			* ( m_avKnots.get( m_cvVerts ).value - m_avKnots.get( m_vOrder - 1 ).value )
			+ m_avKnots.get( m_vOrder - 1 ).value;
			
			int extra = 0;
			int last = m_cvVerts + m_vOrder - 1;
			float midVal = sv;
			int middex = FindSpanV( midVal );
			// Search forward and backward to see if multiple knot is already there
			int i = 0;
			int same = 0;
			if ( avKnots().get( middex ).value == midVal )
			{
				i = middex + 1;
				same = 1;
				while ( ( i < last ) && ( avKnots().get( i ).value == midVal ) ){
					i++;
					same++;
				}
				
				i = middex - 1;
				while ( ( i > 0 ) && ( avKnots().get( i ).value == midVal ) ){
					i--;
					middex--;	// middex is start of multiple knot
					same++;
				}
			}
			
			if ( i <= 0L )           	    // No knot in middle, must create it
			{
				middex = 0;
				while ( avKnots().get( middex + 1 ).value < midVal )
					middex++;
				same = 0;
			}
			
			extra = m_vOrder - same;
			//nttdata 高速化
			//STLVector<p_float> anewKnots = new STLVector<p_float>( p_float.class, extra );
			STLArray<p_float> anewKnots = new STLArray<p_float>( p_float.class, extra );
			
			if ( same < m_vOrder )           	    // Must add knots
			{
				for ( i = 0; i < extra; i++ )
					anewKnots.get( i ).value = midVal ;
			}
			
			vSplitPoint[ iv ] = ( extra < m_vOrder ) ? middex - 1 : middex;
			RefineKnotV( anewKnots );
		}
		
		// Fill in the end points for the last split.
		uSplitPoint[ uSplits ] = m_cuVerts - 1;
		vSplitPoint[ vSplits ] = m_cvVerts - 1;
		
		// Now go over the surface, generating the new patches at the split points in the arrays.
		int uPatch, vPatch;
		// Initialise the offset for the first segment.
		int vOffset = 0;
		for ( vPatch = 0; vPatch < vSplits; vPatch++ )
		{
			// Initialise the offset for the first segment.
			int uOffset = 0;
			// Get the end of the next segment in v.
			int vEnd = vSplitPoint[ vPatch + 1 ];
			
			// Loop across u rows, filling points and knot vectors.
			for ( uPatch = 0; uPatch < uSplits; uPatch++ )
			{
				int uEnd = uSplitPoint[ uPatch + 1 ];
				
				// The index of the patch we are working on.
				int iS = ( vPatch * uSplits ) + uPatch;
				S.set( iS, new CqSurfaceNURBS() );
				S.get( iS ).SetfPatchMesh( false );
				// Initialise it to the same orders as us, with the calculated control point densities.
				S.get( iS ).Init( m_uOrder, m_vOrder, ( uEnd + 1 ) - uOffset, ( vEnd + 1 ) - vOffset );
				
				// Copy any 'vertex' class user primitive variables.
				int iPu, iPv;
				for ( CqParameter iUP : aUserParams() )
				{
					if ( iUP.Class().getValue() == class_vertex )
					{
						CqParameter pNewUP = iUP.CloneType( iUP.strName(), iUP.Count());
						pNewUP.SetSize( S.get( iS ).cVertex() );
						
						for ( iPv = 0; iPv <= vEnd - vOffset; iPv++ )
						{
							int iPIndex = ( ( vOffset + iPv ) * m_cuVerts ) + uOffset;
							for ( iPu = 0; iPu <= uEnd - uOffset; iPu++ )
							{
								int iSP = ( iPv * S.get( iS ).cuVerts() ) + iPu;
								pNewUP.SetValue( iUP, iSP, iPIndex++ );
							}
						}
						S.get( iS ).AddPrimitiveVariable( pNewUP );
					}
				}
				
				// Copy the knot vectors
				int iuK, ivK;
				for ( iuK = 0; iuK < S.get( iS ).uOrder() + S.get( iS ).cuVerts(); iuK++ )
					S.get( iS ).auKnots().get( iuK ).value = auKnots().get( uOffset + iuK ).value;
				for ( ivK = 0; ivK < S.get( iS ).vOrder() + S.get( iS ).cvVerts(); ivK++ )
					S.get( iS ).avKnots().get( ivK ).value = avKnots().get( vOffset + ivK ).value;
				
				// Set the offset to just after the end of this segment.
				uOffset = uEnd + 1;
			}
			// Set the offset to just after the end of this segment.
			vOffset = vEnd + 1;
		}
		
		// Now setup any 'varying', 'uniform' or 'constant' class variables on the segments.
		int irow, icol;
		int nuSegs = uSplits;
		int nvSegs = vSplits;
		for ( icol = 0; icol < nvSegs; icol++ )
		{
			for ( irow = 0; irow < nuSegs; irow++ )
			{
				int iPatch = ( icol * nuSegs ) + irow;
				int iA = ( icol * ( nuSegs + 1 ) ) + irow;
				int iB = ( icol * ( nuSegs + 1 ) ) + irow + 1;
				int iC = ( ( icol + 1 ) * ( nuSegs + 1 ) ) + irow;
				int iD = ( ( icol + 1 ) * ( nuSegs + 1 ) ) + irow + 1;
				
				for ( CqParameter iUP : aUserParams() )
				{
					if ( iUP.Class().getValue() == class_varying )
					{
						CqParameter  pNewUP = iUP.CloneType( iUP.strName(), iUP.Count() );
						pNewUP.SetSize( 4 );
						pNewUP.SetValue( iUP, 0, iA );
						pNewUP.SetValue( iUP, 1, iB );
						pNewUP.SetValue( iUP, 2, iC );
						pNewUP.SetValue( iUP, 3, iD );
						
						S.get( iPatch ).AddPrimitiveVariable( pNewUP );
					}
					else if ( iUP.Class().getValue() == class_uniform )
					{
						CqParameter  pNewUP = iUP.CloneType( iUP.strName(), iUP.Count() );
						pNewUP.SetSize( 1 );
						pNewUP.SetValue( iUP, 0, iPatch );
						S.get( iPatch ).AddPrimitiveVariable( pNewUP );
					}
					else if ( iUP.Class().getValue() == class_constant )
					{
						CqParameter pNewUP = iUP.Clone( );
						S.get( iPatch ).AddPrimitiveVariable( pNewUP );
					}
				}
			}
		}
	}
	
	
	/**
	 * Insert the specified knots into the U knot vector, and refine the control points accordingly.
	 * 
	 * 
	 * @param X
	 */
	@SuppressWarnings("unchecked")
	//nttdata 高速化
	//public  void RefineKnotU( final STLVector<p_float> X )
	public  void RefineKnotU( final STLArray<p_float> X )
	{
		if ( X.size() <= 0 )
			return ;
		
		int n = m_cuVerts - 1;
		int p = uDegree();
		int m = n + p + 1;
		int a, b;
		int r = X.size() - 1;
		int j, row;
		
		// Find the insertion points for the start and end of the knot vector to be inserted.
		a = FindSpanU( X.get( 0 ).value );
		b = FindSpanU( X.get( r ).value );
		++b;
		
		int i = b + p - 1;
		int k = b + p + r;
		m_cuVerts = r + 1 + n + 1;
		
		//nttdata 高速化
		//STLVector<p_float>	auHold = new STLVector<p_float>( p_float.class, m_auKnots.size() );
		STLArray<p_float>	auHold = new STLArray<p_float>( p_float.class, m_auKnots.size() );
		for( j = 0; j < m_auKnots.size(); j++ )	
			auHold.get( j ).value = m_auKnots.get( j ).value ;
		
		m_auKnots.resize( m_cuVerts + m_uOrder );
		
		// Copy the knot values up to the first insertion point.
		for ( j = 0; j <= a; j++ )
			m_auKnots.get( j ).value = auHold.get( j ).value;
		
		// Copy the knot values from the second insertion point to the end.
		for ( j = b + p; j <= m; j++ )
			m_auKnots.get( j + r + 1 ).value = auHold.get( j ).value;
		
		for ( j = r; j >= 0; j-- ){
			while ( X.get( j ).value <= m_auKnots.get( i ).value && i > a )
				m_auKnots.get( k-- ).value = auHold.get( i-- ).value;
			m_auKnots.get( k-- ).value = X.get( j ).value;
		}
		
		i = b + p - 1;
		k = b + p + r;
		
		// Now process all the 'vertex' class variables.
		for ( CqParameter iUP : m_aUserParams )
		{
			if ( iUP.Class().getValue() == class_vertex )
			{
				i = b + p - 1;
				k = b + p + r;
				
				CqParameter pHold = iUP.Clone();
				iUP.SetSize( m_cuVerts * m_cvVerts );
				
				// Copy the control points from the original
				for ( row = 0; row < m_cvVerts; row++ )
				{
					// Copy the CPs up to the first insertion point minus the degree (this is the number of control points in that section).
					int rowoff = ( row * m_cuVerts );
					for ( j = 0; j <= a - p ; j++ ) 
						iUP.SetValue( pHold, rowoff + j, ( row * ( n + 1 ) ) + j );
					// Copy the CPs beyond the second insertion point to the end.
					for ( j = b - 1; j <= n; j++ )
						iUP.SetValue( pHold, rowoff + j + r + 1, ( row * ( n + 1 ) ) + j );
				}
				
				for ( j = r; j >= 0; j-- )
				{
					while ( X.get( j ).value <= m_auKnots.get( i ).value && i > a )
					{
						for ( row = 0; row < m_cvVerts; row++ )
							iUP.SetValue( pHold, ( row * m_cuVerts ) + k - p - 1, ( row * ( n + 1 ) ) + i - p - 1 );
						--k;
						--i;
					}
					
					for ( row = 0; row < m_cvVerts; row++ )
						iUP.SetValue( iUP, ( row * m_cuVerts ) + k - p - 1, ( row * m_cuVerts ) + k - p );
					
					int l;
					for ( l = 1; l <= p ; l++ )
					{
						int ind = k - p + l;
						float alpha = m_auKnots.get( k + l ).value - X.get( j ).value;
						if ( alpha == 0.0f )
						{
							for ( row = 0; row < m_cvVerts; row++ )
								iUP.SetValue( iUP, ( row * m_cuVerts ) + ind - 1 , ( row * m_cuVerts ) + ind );
						}
						else
						{
							alpha /= m_auKnots.get( k + l ).value - auHold.get( i - p + l ).value;
							// Make sure index is OK.
							switch ( iUP.Type().getValue() )
							{
							case type_float:
							{
								CqParameterTyped<p_float, p_float> pTParam = ( CqParameterTyped<p_float, p_float> ) iUP ;
								for ( row = 0; row < m_cvVerts; row++ )
									((p_float)pTParam.pValue_get( ( row * m_cuVerts ) + ind - 1 , 0 )).value = alpha * ((p_float)pTParam.pValue_get(0, ( row * m_cuVerts ) + ind - 1 )).value + ( 1.0f - alpha ) * ((p_float) pTParam.pValue_get(0,(row * m_cuVerts) + ind)).value;
								break;
							}
							case type_integer:
							{
								CqParameterTyped<p_int, p_float> pTParam = ( CqParameterTyped<p_int, p_float> ) iUP ;
								for ( row = 0; row < m_cvVerts; row++ )
									((p_int)pTParam.pValue_get( ( row * m_cuVerts ) + ind - 1 , 0 )).value =
										(int)( alpha * ((p_int)pTParam.pValue_get(0,(row * m_cuVerts ) + ind - 1)).value + ( 1.0f - alpha ) * ((p_int)pTParam.pValue_get(0,(row * m_cuVerts ) + ind )).value );
								break;
							}
							case type_point:
							case type_normal:
							case type_vector:
							{
								CqParameterTyped<CqVector3D, CqVector3D> pTParam = ( CqParameterTyped<CqVector3D, CqVector3D> ) iUP;
								for ( row = 0; row < m_cvVerts; row++ )
									((CqVector3D)pTParam.pValue_get(( row * m_cuVerts ) + ind - 1, 0 )).assignment( ( ((CqVector3D)pTParam.pValue_get(0,(row * m_cuVerts) + ind - 1)).mul( alpha ) ).add( ((CqVector3D)pTParam.pValue_get(0,(row * m_cuVerts) + ind)).mul( 1.0f - alpha ) ));
								break;
							}
							case type_hpoint:
							{
								CqParameterTyped<CqVector4D, CqVector3D> pTParam = ( CqParameterTyped<CqVector4D, CqVector3D> ) iUP;
								for ( row = 0; row < m_cvVerts; row++ )
								{
									CqVector4D cp = new CqVector4D( 
											alpha * ((CqVector4D)pTParam.pValue_get(0,( row * m_cuVerts )+ ind - 1)).x + (1.0f - alpha) * ((CqVector4D)pTParam.pValue_get(0,(row * m_cuVerts)+ ind )).x,
											alpha * ((CqVector4D)pTParam.pValue_get(0,( row * m_cuVerts )+ ind - 1)).y + (1.0f - alpha) * ((CqVector4D)pTParam.pValue_get(0,(row * m_cuVerts)+ ind )).y,
											alpha * ((CqVector4D)pTParam.pValue_get(0,( row * m_cuVerts )+ ind - 1)).z + (1.0f - alpha) * ((CqVector4D)pTParam.pValue_get(0,(row * m_cuVerts)+ ind )).z,
											alpha * ((CqVector4D)pTParam.pValue_get(0,( row * m_cuVerts )+ ind - 1)).w + (1.0f - alpha) * ((CqVector4D)pTParam.pValue_get(0,(row * m_cuVerts)+ ind )).w );
									
									((CqVector4D)pTParam.pValue_get(( row * m_cuVerts ) + ind - 1 , 0 )).assignment( cp );
								}
								break;
							}
							case type_color:
							{
								CqParameterTyped<CqColor, CqColor> pTParam = ( CqParameterTyped<CqColor, CqColor> ) iUP;
								for ( row = 0; row < m_cvVerts; row++ )
									((CqColor)pTParam.pValue_get(( row * m_cuVerts ) + ind - 1 , 0 )).assignment(( ((CqColor)pTParam.pValue_get(0,( row * m_cuVerts ) + ind - 1)).mul( alpha )).add( ((CqColor)pTParam.pValue_get(0,( row * m_cuVerts ) + ind)).mul( 1.0f - alpha )));
								break;
							}
							case type_string: 
							{
								CqParameterTyped<p_String, p_String> pTParam = ( CqParameterTyped<p_String, p_String> ) iUP;
								for ( row = 0; row < m_cvVerts; row++ )
									((p_String)pTParam.pValue_get(( row * m_cuVerts ) + ind - 1, 0 )).value = ((p_String)pTParam.pValue_get(0,( row * m_cuVerts ) + ind - 1)).value + ((p_String)pTParam.pValue_get(0,( row * m_cuVerts ) + ind)).value;
//								((p_String)pTParam.pValue_get(( row * m_cuVerts ) + ind - 1, 0 )).value = alpha * ((p_String)pTParam.pValue_get(,( row * m_cuVerts ) + ind - 1)).value + ( 1.0f - alpha ) * ((p_String)pTParam.pValue_get(,( row * m_cuVerts ) + ind)).value;
								break;
							}
							case type_matrix:
							{
								CqParameterTyped<CqMatrix, CqMatrix> pTParam = ( CqParameterTyped<CqMatrix, CqMatrix> ) iUP;
								for ( row = 0; row < m_cvVerts; row++ )
									((CqMatrix)pTParam.pValue_get(( row * m_cuVerts ) + ind - 1, 0 )).assignment( (((CqMatrix)pTParam.pValue_get(0,( row * m_cuVerts ) + ind - 1)).mulInv( alpha )).add( ((CqMatrix)pTParam.pValue_get(0,( row * m_cuVerts ) + ind)).mulInv( 1.0f - alpha ) ));
								break;
							}
							default:
							{
								// left blank to avoid compiler warnings about unhandled types
								break;
							}
							}
						}
					}
					--k;
				}
				pHold = null;
			}
		}
	}
	
	/**
	 * Insert the specified knots into the V knot vector,
	 * and refine the control points accordingly.
	 * 
	 * @param X
	 */
	@SuppressWarnings("unchecked")
	//nttdata 高速化
	//public  void RefineKnotV( final STLVector<p_float> X )
	public  void RefineKnotV( final STLArray<p_float> X )
	{
		if ( X.size() <= 0 )
			return ;
		
		int n = m_cvVerts - 1;
		int p = vDegree();
		int m = n + p + 1;
		int a, b;
		int r = X.size() - 1;
		int j, col;
		
		a = FindSpanV( X.get( 0 ).value ) ;
		b = FindSpanV( X.get( r ).value ) ;
		++b;
		
		int i = b + p - 1;
		int k = b + p + r;
		m_cvVerts = r + 1 + n + 1;
		
		//nttdata 高速化
		//STLVector<p_float>	avHold = new STLVector<p_float>( p_float.class, m_avKnots.size() );
		STLArray<p_float>	avHold = new STLArray<p_float>( p_float.class, m_avKnots.size() );
		for( j = 0; j < m_avKnots.size(); j++ )	avHold.get( j ).value = m_avKnots.get( j ).value;
		
		m_avKnots.resize( m_cvVerts + m_vOrder );
		
		for ( j = 0; j <= a; j++ )
			m_avKnots.get( j ).value = avHold.get( j ).value;
		
		for ( j = b + p; j <= m; j++ )
			m_avKnots.get( j + r + 1 ).value = avHold.get( j ).value;
		
		for ( j = r; j >= 0 ; j-- ){
			while ( X.get( j ).value <= m_avKnots.get( i ).value && i > a )
				m_avKnots.get( k-- ).value = avHold.get( i-- ).value;
			m_avKnots.get( k-- ).value = X.get( j ).value;
		}
		
		i = b + p - 1;
		k = b + p + r;
		
		// Now process all the 'vertex' class variables.
		for ( CqParameter iUP : m_aUserParams ){
			if ( iUP.Class().getValue() == class_vertex )
			{
				i = b + p - 1;
				k = b + p + r;
				
				CqParameter pHold = iUP.Clone();
				iUP.SetSize( m_cuVerts * m_cvVerts );
				
				for ( col = 0; col <m_cuVerts; col++ ){
					for ( j = 0; j <= a - p; j++ ) iUP.SetValue( pHold, ( j * m_cuVerts ) + col, ( j * m_cuVerts ) + col );
					for ( j = b - 1; j <= n; j++ ) iUP.SetValue( pHold, ( ( j + r + 1 ) * m_cuVerts ) + col, ( j * m_cuVerts ) + col );
				}
				
				
				for ( j = r; j >= 0 ; j-- )
				{
					while ( X.get( j ).value <= m_avKnots.get( i ).value && i > a )
					{
						for ( col = 0; col <m_cuVerts; col++ )
							iUP.SetValue( pHold, ( ( k - p - 1 ) * m_cuVerts ) + col, ( ( i - p - 1 ) * m_cuVerts ) + col );
						--k;
						--i;
					}
					for ( col = 0; col < m_cuVerts; col++ )
						iUP.SetValue( iUP, ( ( k - p - 1 ) * m_cuVerts ) + col, ( ( k - p ) * m_cuVerts ) + col );
					
					int l;
					for ( l = 1; l <= p; l++ )
					{
						int ind = k - p + l;
						float alpha = m_avKnots.get( k + l ).value - X.get( j ).value;
						if ( alpha == 0.0f )
						{
							for ( col = 0; col <m_cuVerts; col++ )
								iUP.SetValue( iUP, ( ( ind - 1 ) * m_cuVerts ) + col, ( ind * m_cuVerts ) + col );
						}
						else
						{
							alpha /= m_avKnots.get( k + l ).value - avHold.get( i - p + l ).value;
							switch ( iUP.Type().getValue() )
							{
							case type_float:
							{
								CqParameterTyped<p_float, p_float> pTParam = ( CqParameterTyped<p_float, p_float> ) iUP;
								for ( col = 0; col <m_cuVerts ; col++ )
									((p_float)pTParam.pValue_get( 0 ,(( ind - 1 ) * m_cuVerts) + col)).value = alpha * ((p_float)pTParam.pValue_get( 0 ,(( ind - 1 ) * m_cuVerts ) + col)).value + ( 1.0f - alpha ) * ((p_float)pTParam.pValue_get( 0 ,( ind * m_cuVerts ) + col)).value;
								break;
							}
							case type_integer:
							{
								CqParameterTyped<p_int, p_float> pTParam = ( CqParameterTyped<p_int, p_float> ) iUP;
								for ( col = 0; col <m_cuVerts; col++ )
									((p_int)pTParam.pValue_get( 0 ,(( ind - 1 ) * m_cuVerts ) + col)).value = 
										(int)( alpha * ((p_int)pTParam.pValue_get( 0 ,(( ind - 1 ) * m_cuVerts ) + col)).value + ( 1.0f - alpha ) * ((p_int)pTParam.pValue_get( 0 ,(ind * m_cuVerts) + col)).value );
								break;
							}
							case type_point:
							case type_normal:
							case type_vector:
							{
								CqParameterTyped<CqVector3D, CqVector3D> pTParam = ( CqParameterTyped<CqVector3D, CqVector3D> ) iUP;
								for ( col = 0; col <m_cuVerts; col++ )
									((CqVector3D)pTParam.pValue_get( 0 ,(( ind - 1 ) * m_cuVerts ) + col)).assignment(( ((CqVector3D)pTParam.pValue_get( 0 ,(( ind - 1 )* m_cuVerts) + col)).mul( alpha ) ).add( ((CqVector3D)pTParam.pValue_get( 0 ,(ind * m_cuVerts) + col)).mul( 1.0f - alpha ) ));
								break;
							}
							case type_hpoint:
							{
								CqParameterTyped<CqVector4D, CqVector3D> pTParam = ( CqParameterTyped<CqVector4D, CqVector3D> ) iUP;
								for ( col = 0; col <m_cuVerts; col++ )
								{
									CqVector4D cp = new CqVector4D (
											alpha * ((CqVector4D)pTParam.pValue_get( 0 , ( ( ind - 1 ) * m_cuVerts ) + col )).x + ( 1.0f - alpha ) * ((CqVector4D)pTParam.pValue_get( 0 , ( ind * m_cuVerts ) + col )).x,
											alpha * ((CqVector4D)pTParam.pValue_get( 0 , ( ( ind - 1 ) * m_cuVerts ) + col )).y + ( 1.0f - alpha ) * ((CqVector4D)pTParam.pValue_get( 0 , ( ind * m_cuVerts ) + col )).y,
											alpha * ((CqVector4D)pTParam.pValue_get( 0 , ( ( ind - 1 ) * m_cuVerts ) + col )).z + ( 1.0f - alpha ) * ((CqVector4D)pTParam.pValue_get( 0 , ( ind * m_cuVerts ) + col )).z,
											alpha * ((CqVector4D)pTParam.pValue_get( 0 , ( ( ind - 1 ) * m_cuVerts ) + col )).w + ( 1.0f - alpha ) * ((CqVector4D)pTParam.pValue_get( 0 , ( ind * m_cuVerts ) + col )).w );

									((CqVector4D)pTParam.pValue_get( 0 , (( ind - 1 ) * m_cuVerts ) + col )).assignment( cp );
								}
								break;
							}
							
							case type_color:
							{
								CqParameterTyped<CqColor, CqColor> pTParam = ( CqParameterTyped<CqColor, CqColor> ) iUP;
								for ( col = 0; col < m_cuVerts; col++ )
									((CqColor)pTParam.pValue_get( 0 ,(( ind - 1 ) * m_cuVerts) + col)).assignment(( ((CqColor)pTParam.pValue_get( 0 ,(( ind - 1 ) * m_cuVerts) + col)).mul( alpha )).add( ((CqColor)pTParam.pValue_get( 0 ,(ind * m_cuVerts) + col)).mul(  1.0f - alpha )));
								break;
							}
							case type_string:
							{	
								CqParameterTyped<p_String, p_String> pTParam = ( CqParameterTyped<p_String, p_String> ) iUP;
								for ( col = 0; col < m_cuVerts; col++ )
									((p_String)pTParam.pValue_get( 0 ,(( ind - 1 ) * m_cuVerts) + col)).value = ((p_String)pTParam.pValue_get( 0 ,(( ind - 1 ) * m_cuVerts) + col)).value + ((p_String)pTParam.pValue_get( 0 ,(ind * m_cuVerts) + col)).value ;
//								((p_String)pTParam.pValue_get( 0 ,(( ind - 1 ) * m_cuVerts) + col)).value = alpha * ((p_String)pTParam.pValue_get( 0 ,(( ind - 1 ) * m_cuVerts) + col)).value + ( 1.0f - alpha ) * ((p_String)pTParam.pValue_get( 0 ,(ind * m_cuVerts) + col)).value ;
								break;
							}
							case type_matrix:
							{
								CqParameterTyped<CqMatrix, CqMatrix> pTParam = ( CqParameterTyped<CqMatrix, CqMatrix> ) iUP;
								for ( col = 0; col < m_cuVerts; col++ )
									((CqMatrix)pTParam.pValue_get( 0 ,(( ind - 1 ) * m_cuVerts ) + col)).assignment(( ((CqMatrix)pTParam.pValue_get( 0 ,(( ind - 1 )* m_cuVerts ) + col)).mulInv( alpha )).add( ((CqMatrix)pTParam.pValue_get( 0 ,(ind* m_cuVerts) + col)).mulInv( 1.0f - alpha ) ));
								break;
							}
							default:
							{
								// left blank to avoid compiler warnings about unhandled types
								break;
							}
							}
						}
					}
					--k;
				}
				pHold = null;
			}
		}
	}
	
	/**
	 * Insert the specified knot into the U knot vector, 
	 * and refine the control points accordingly.
	 * 
	 * @param u
	 * @param r
	 * @return	The number of new knots created.
	 */
	@SuppressWarnings("unchecked")
	public  int	InsertKnotU( float u, int r )
	{
		int n = m_cuVerts;
		int k = m_auKnots.size() - 1;
		int s = 0;
		int i, j;
		int p = uDegree();
		
		// If the specified u value falls outside the current range, then fail.
		if ( u < m_auKnots.get( uDegree()).value || u > m_auKnots.get( m_cuVerts ).value )
			return ( 0 );
		
		// Calculate k as the index of the last knot value <= the specified value.
		for ( i = 0; i < m_auKnots.size(); i++ ){
			if ( m_auKnots.get( i ).value > u ){
				k = i - 1;
				break;
			}
		}
		
		// Calculate the number of knots at the insertion point with the same value as the specified knot.
		if ( u <= m_auKnots.get( k ).value ){
			s = 1;
			for ( i = k; i > uDegree(); i-- ){
				if ( m_auKnots.get( i ).value <= m_auKnots.get( i - 1 ).value )
					s++;
				else
					break;
			}
		}
		else
			s = 0;
		
		// Adjust the number of insertions to take into account the number of knots with that value already in the vector.
		if ( ( r + s ) > p + 1 )
			r = p + 1 - s;
		
		// If this means we don't have to do anything then exit.
		if ( r <= 0 )
			return ( 0 );
		
		m_cuVerts = m_cuVerts + r;
		m_auKnots.resize( m_cuVerts + m_uOrder );
		
		//nttdata 高速化
		//STLVector<p_float> auHold = new STLVector<p_float>( p_float.class, m_auKnots.size() ); 
		STLArray<p_float> auHold = new STLArray<p_float>( p_float.class, m_auKnots.size() );
		for( i = 0; i < m_auKnots.size(); i++ )	
			auHold.get( i ).value = m_auKnots.get( i ).value;
		
		// Load new knot vector
		// Copy up to the insertion point.
		for ( i = 0; i <= k;i++ )
			m_auKnots.get( i ).value = auHold.get( i ).value; 
		
		// Add the specified value 'r' times at the insertion point.
		for ( i = 1; i <= r;i++ )
			m_auKnots.get( k + i ).value = u;
		
		// Copy after the insertion point up to the end.
		// 2005.02.14.Start.change
		//  for ( i = k + 1; i < static_cast<int>( m_auKnots.size() ); i++ )
		for ( i = k + 1; i + r < m_auKnots.size(); i++ )
			// 2005.02.14.End  .change
			m_auKnots.get( i + r ).value = auHold.get( i ).value;
		
		// Now process all the 'vertex' class variables.
		for ( CqParameter iUP : m_aUserParams )
		{
			if ( iUP.Class().getValue() == class_vertex )
			{
				CqParameter pHold = iUP.Clone();
				iUP.SetSize( m_cuVerts * m_cvVerts );
				
				// Save unaltered control points
				CqParameter R = iUP.CloneType( "R" );
				R.SetSize( p + 1 );
				
				// Insert control points as required on each row.
				int row;
				for ( row = 0; row < m_cvVerts; row++ )
				{
					// First copy the first set of control points up to the insertion point minus the degree
					for ( i = 0; i <= k - p; i++ )
						// Qw[i][row] = Pw[i][row]
						iUP.SetValue( pHold, ( row * m_cuVerts ) + i, ( row * n ) + i );
					for ( i = k - s; i < m_cuVerts - r; i++ )
						// Qw[i+r][row] = Pw[i][row]
						iUP.SetValue( pHold, ( row * m_cuVerts ) + i + r, ( row * n ) + i );
					for ( i = 0; i <= p - s; i++ )
						// Rw[i] = Pw[k-p+i][row]
						R.SetValue( pHold, i, ( row * n ) + k - p + i );
					
					// Insert the knot r times
					int L = 0 ;
					float alpha;
					for ( j = 1; j <= r; j++ )
					{
						L = k - p + j;
						for ( i = 0;i <= p - j - s;i++ )
						{
							alpha = ( u - auHold.get( L + i ).value ) / ( auHold.get( i + k + 1 ).value - auHold.get( L + i ).value );
							
							switch ( iUP.Type().getValue() )
							{
							case type_float:
							{
								CqParameterTyped<p_float, p_float> pTR = ( CqParameterTyped<p_float, p_float> ) R ;
								((p_float)pTR.pValue_get( i ,0)).value = alpha * ((p_float)pTR.pValue_get( i + 1 ,0)).value + ( 1.0f - alpha ) * ((p_float)pTR.pValue_get( i ,0)).value ;
								break;
							}
							case type_integer:
							{
								CqParameterTyped<p_int, p_float> pTR = ( CqParameterTyped<p_int, p_float> ) R ;
								((p_int)pTR.pValue_get( i ,0)).value = 
									(int)( alpha * ((p_int)pTR.pValue_get( i + 1 ,0)).value + ( 1.0 - alpha ) * ((p_int)pTR.pValue_get( i ,0)).value );
								break;
							}
							case type_point:
							case type_normal:
							case type_vector:
							{
								CqParameterTyped<CqVector3D, CqVector3D> pTR = ( CqParameterTyped<CqVector3D, CqVector3D> ) R ;
								((CqVector3D)pTR.pValue_get( i ,0)).assignment( (((CqVector3D)pTR.pValue_get( i + 1 ,0)).mul( alpha )).add( ((CqVector3D)pTR.pValue_get( i ,0)).mul( 1.0f - alpha )));
								break;
							}
							case type_hpoint:
							{
								CqParameterTyped<CqVector4D, CqVector3D> pTR = ( CqParameterTyped<CqVector4D, CqVector3D> ) R ;
								CqVector4D cp = new CqVector4D ( 
										alpha * ((CqVector4D)pTR.pValue_get( i + 1 ,0)).x + ( 1.0f - alpha ) * ((CqVector4D)pTR.pValue_get( i ,0)).x,
										alpha * ((CqVector4D)pTR.pValue_get( i + 1 ,0)).y + ( 1.0f - alpha ) * ((CqVector4D)pTR.pValue_get( i ,0)).y,
										alpha * ((CqVector4D)pTR.pValue_get( i + 1 ,0)).z + ( 1.0f - alpha ) * ((CqVector4D)pTR.pValue_get( i ,0)).z,
										alpha * ((CqVector4D)pTR.pValue_get( i + 1 ,0)).w + ( 1.0f - alpha ) * ((CqVector4D)pTR.pValue_get( i ,0)).w 
										);
								((CqVector4D)pTR.pValue_get( i ,0)).assignment( cp );
								break;
							}
							case type_color:
							{
								CqParameterTyped<CqColor, CqColor> pTR = ( CqParameterTyped<CqColor, CqColor> ) R ;
								((CqColor)pTR.pValue_get( i ,0)).assignment( (((CqColor)pTR.pValue_get( i + 1 ,0)).mul( alpha )).add( ((CqColor)pTR.pValue_get( i ,0)).mul( 1.0f - alpha )) );
								break;
							}
							//							case type_string:
							//							{
							//								CqParameterTyped<CqString, CqString>* pTR = static_cast<CqParameterTyped<CqString, CqString>*>(R);
							//								(*pTR->pValue_get( i )) = alpha * (*pTR->pValue_get( i + 1 )) + ( 1.0 - alpha ) * (*pTR->pValue_get( i ));
							//								break;
							//							}
							case type_matrix:
							{
								CqParameterTyped<CqMatrix, CqMatrix> pTR = ( CqParameterTyped<CqMatrix, CqMatrix> ) R ;
								((CqMatrix)pTR.pValue_get( i ,0)).assignment( (((CqMatrix)pTR.pValue_get( i + 1 ,0)).mulInv( alpha )).add( ((CqMatrix)pTR.pValue_get( i ,0)).mulInv( 1.0f - alpha )) );
								break;
							}
							default:
							{
								// left blank to avoid compiler warnings about unhandled types in the switch
								break;
							}
							}
						}
						// Qw[L][row] = Rw[0]
						iUP.SetValue( R, (row * m_cuVerts)+ L, 0 );
						if ( p - j - s > 0 )
							// Qw[k+r-j-s][row] = Rw[p-j-s]
							iUP.SetValue( R, (row * m_cuVerts) + k + r - j - s, p - j - s );
					}
					// Load remaining control points
					for ( i = L + 1; i < k - s; i++ )
						// Qw[i][row] = Rw[i-L]
						iUP.SetValue( R, ( row * m_cuVerts ) + i, i - L );
				}
				R = null;
				pHold = null;
			}
		}
		return ( r );
	}
	
	
	/**
	 * Insert the specified knot into the V knot vector, 
	 * and refine the control points accordingly.
	 * 
	 * @param v
	 * @param r
	 * @return	The number of new knots created.
	 */
	@SuppressWarnings("unchecked")
	public  int	InsertKnotV( float v, int r )
	{
		// Compute k and s      v = [ v_k , v_k+1)  with v_k having multiplicity s
		int m = m_cvVerts;
		int k = m_avKnots.size() - 1, s = 0;
		int i, j;
		int p = vDegree();
		
		if ( v < m_avKnots.get( vDegree()).value || v > m_avKnots.get( m_cvVerts ).value )
			return ( 0 );
		
		int size = m_avKnots.size();
		for ( i = 0; i < size; i++ ){
			if ( m_avKnots.get( i ).value > v ) {
				k = i - 1;
				break;
			}
		}
		
		if ( v <= m_avKnots.get( k ).value ) {
			s = 1;
			for ( i = k; i > vDegree(); i-- ){
				if ( m_avKnots.get( i ).value <= m_avKnots.get( i - 1 ).value )
					s++;
				else
					break;
			}
		}
		else
			s = 0;
		
		if ( ( r + s ) > p + 1 )
			r = p + 1 - s;
		
		if ( r <= 0 )
			return ( 0 );
		
		// Work on a copy.
		m_cvVerts = m_cvVerts + r;
		m_avKnots.resize( m_cvVerts + m_vOrder );
		
		//nttdata 高速化
		//STLVector<p_float>	avHold = new STLVector<p_float>( p_float.class, m_avKnots.size() );
		STLArray<p_float>	avHold = new STLArray<p_float>( p_float.class, m_avKnots.size() );
		for( i = 0; i < m_avKnots.size(); i++ ) 
			avHold.get( i ).value = m_avKnots.get( i ).value;
		
		// Load new knot vector
		for ( i = 0; i <= k;i++ )
			m_avKnots.get( i ).value = avHold.get( i ).value;
		for ( i = 1; i <= r;i++ )
			m_avKnots.get( k + i ).value = v;
		size = m_avKnots.size();
		// 2005.02.14.Start.change
		// for ( i = k + 1;i < size; i++ )
		for ( i = k + 1; i + r < size; i++ )
			// 2005.02.14.End  .change
			m_avKnots.get( i + r ).value = avHold.get( i ).value;
		
		// Now process all the 'vertex' class variables.
		for ( CqParameter iUP : m_aUserParams )
		{
			if ( iUP.Class().getValue() == class_vertex )
			{
				CqParameter pHold = iUP.Clone();
				iUP.SetSize( m_cuVerts * m_cvVerts );
				
				// Save unaltered control points
				CqParameter R = iUP.CloneType( "R" );
				R.SetSize( p + 1 );
				
				// Insert control points as required on each row.
				for ( int col = 0; col < m_cuVerts; col++ )
				{
					for ( i = 0; i <= k - p; i++ )
						// Qw[col][i] = Pw[col][i]
						iUP.SetValue( pHold, (i * m_cuVerts)+ col, (i * m_cuVerts)+ col );
					for ( i = k - s; i < m ; i++ )
						// Qw[col][i+r] = Pw[col][i]
						iUP.SetValue( pHold, (( i + r )* m_cuVerts ) + col, (i * m_cuVerts)+ col );
					for ( i = 0; i <= p - s; i++ )
						// Rw[i] = Pw[col][k-p+i]
						R.SetValue( pHold, i, (( k - p + i )* m_cuVerts ) + col );
					
					// Insert the knot r times
					int L = 0 ;
					float alpha;
					for ( j = 1; j <= r; j++ )
					{
						L = k - p + j;
						for ( i = 0;i <= p - j - s;i++ )
						{
							alpha = ( v - avHold.get( L + i ).value ) / ( avHold.get( i + k + 1 ).value - avHold.get( L + i ).value );
							
							switch ( iUP.Type().getValue() )
							{
							case type_float:
							{
								CqParameterTyped<p_float, p_float> pTR = ( CqParameterTyped<p_float, p_float> ) R ;
								((p_float) pTR.pValue_get( i ,0)).value = alpha * ((p_float) pTR.pValue_get( i + 1 ,0)).value + ( 1.0f - alpha ) * ((p_float) pTR.pValue_get( i ,0)).value;
								break;
							}
							case type_integer:
							{
								CqParameterTyped<p_int, p_float> pTR = ( CqParameterTyped<p_int, p_float> ) R ;
								((p_int) pTR.pValue_get( i ,0)).value = (int)(
										alpha * ((p_int) pTR.pValue_get( i + 1 ,0)).value + ( 1.0f - alpha ) * ((p_int) pTR.pValue_get( i ,0)).value );
								break;
							}
							case type_point:
							case type_normal:
							case type_vector:
							{
								CqParameterTyped<CqVector3D, CqVector3D> pTR = ( CqParameterTyped<CqVector3D, CqVector3D> ) R ;
								((CqVector3D) pTR.pValue_get( i ,0)).assignment( (((CqVector3D)pTR.pValue_get( i + 1 ,0)).mul( alpha ) ).add( ((CqVector3D) pTR.pValue_get( i ,0)).mul( 1.0f - alpha ) ));
								break;
							}
							
							case type_hpoint:
							{
								CqParameterTyped<CqVector4D, CqVector3D> pTR = ( CqParameterTyped<CqVector4D, CqVector3D> ) R ;
								CqVector4D cp = new CqVector4D (
										alpha * ((CqVector4D) pTR.pValue_get( i + 1 ,0)).x + ( 1.0f - alpha ) * ((CqVector4D) pTR.pValue_get( i ,0)).x,
										alpha * ((CqVector4D) pTR.pValue_get( i + 1 ,0)).y + ( 1.0f - alpha ) * ((CqVector4D) pTR.pValue_get( i ,0)).y,
										alpha * ((CqVector4D) pTR.pValue_get( i + 1 ,0)).z + ( 1.0f - alpha ) * ((CqVector4D) pTR.pValue_get( i ,0)).z,
										alpha * ((CqVector4D) pTR.pValue_get( i + 1 ,0)).w + ( 1.0f - alpha ) * ((CqVector4D) pTR.pValue_get( i ,0)).w 
								);
								((CqVector4D) pTR.pValue_get( i ,0)).assignment( cp );
								break;
							}
							
							case type_color:
							{
								CqParameterTyped<CqColor, CqColor> pTR = ( CqParameterTyped<CqColor, CqColor> ) R ;
								((CqColor) pTR.pValue_get( i ,0)).assignment( (((CqColor) pTR.pValue_get( i + 1 ,0)).mul(alpha)).add( ((CqColor) pTR.pValue_get( i ,0)).mul( 1.0f - alpha )) );
								break;
							}
							
							//							case type_string:
							//							{
							//								CqParameterTyped<CqString, CqString>* pTR = static_cast<CqParameterTyped<CqString, CqString>*>(R);
							//								(*pTR.pValue_get( i )) = alpha * (*pTR.pValue_get( i + 1 )) + ( 1.0 - alpha ) * (*pTR.pValue_get( i ));
							//								break;
							//							}
							
							case type_matrix:
							{
								CqParameterTyped<CqMatrix, CqMatrix> pTR = ( CqParameterTyped<CqMatrix, CqMatrix> ) R ;
								((CqMatrix) pTR.pValue_get( i ,0)).assignment( (((CqMatrix) pTR.pValue_get( i + 1 ,0)).mulInv( alpha )).add( ((CqMatrix) pTR.pValue_get( i ,0)).mulInv( 1.0f - alpha )) );
								break;
							}
							
							default:
							{
								// left blank to avoid compiler warnings about unhandled types
								break;
							}
							}
						}
						// Qw[col][L] = Rw[0]
						iUP.SetValue( R, (L * m_cuVerts) + col, 0 );
						if ( p - j - s > 0 )
							// Qw[col][k+r-j-s] = Rw[p-j-s]
							iUP.SetValue( R, (( k + r - j - s ) * m_cuVerts ) + col, p - j - s );
					}
					// Load remaining control points
					for ( i = L + 1; i < k - s; i++ )
						// Qw[col][i] = Rw[i-L]
						iUP.SetValue( R, ( i * m_cuVerts ) + col, i - L );
				}
				R = null;
				pHold = null;
			}
		}
		return ( r );
	}
	
	/**
	 * Ensure a nonperiodic (clamped) knot vector by inserting U[p] and U[m-p] multiple times.
	 * 
	 * 
	 */
	public  void	ClampU()
	{
		float u1 = m_auKnots.get( uDegree() ).value; 
		float u2 = m_auKnots.get( m_cuVerts ).value; 
		
		int n1 = InsertKnotU( u1, uDegree() );
		int n2 = InsertKnotU( u2, uDegree() );
		int i = 0;
		
		// Now trim unnecessary knots and control points
		if ( n1 != 0 || n2 != 0 )
		{
			//nttdata 高速化
			//STLVector<p_float> auHold = new STLVector<p_float>( p_float.class, m_auKnots.size() );
			STLArray<p_float> auHold = new STLArray<p_float>( p_float.class, m_auKnots.size() );
			for( i = 0; i < m_auKnots.size(); i++ ) 
				auHold.get( i ).value = m_auKnots.get( i ).value;
			
			m_auKnots.resize( m_auKnots.size() - n1 - n2 );
			
			for ( i = n1; i < auHold.size() - n2; i++ )
				m_auKnots.get( i - n1 ).value = auHold.get( i ).value;
			
			int n = m_cuVerts;
			m_cuVerts -= n1 + n2;
			
			// Now process all the 'vertex' class variables.
			for ( CqParameter iUP : m_aUserParams )
			{
				if ( iUP.Class().getValue() == class_vertex )
				{
					CqParameter pHold = iUP.Clone();
					iUP.SetSize( m_cuVerts * m_cvVerts );
					
					int row;
					for ( row = 0; row < m_cvVerts; row++ )
					{
						for ( i = n1; i < n - n2; i++ )
							iUP.SetValue( pHold, row * m_cuVerts + i - n1, row * n + i );
					}
					pHold = null;
				}
			}
		}
	}
	
	/**
	 * Ensure a nonperiodic (clamped) knot vector by inserting V[p] and V[m-p] multiple times.
	 * 
	 * 
	 */
	public  void	ClampV()
	{
		float v1 = m_avKnots.get( vDegree() ).value;
		float v2 = m_avKnots.get( m_cvVerts ).value;
		
		int n1 = InsertKnotV( v1, vDegree() );
		int n2 = InsertKnotV( v2, vDegree() );
		
		int i = 0;  
		// Now trim unnecessary knots and control points
		if ( n1 != 0 || n2 != 0 )
		{
			//nttdata 高速化
			//STLVector<p_float> avHold = new STLVector<p_float>( p_float.class, m_avKnots.size() );
			STLArray<p_float> avHold = new STLArray<p_float>( p_float.class, m_avKnots.size() );
			for( i = 0; i < m_avKnots.size(); i++ )	
				avHold.get( i ).value = m_avKnots.get( i ).value;
			
			m_avKnots.resize( m_avKnots.size() - n1 - n2 );
			
			for ( i = n1; i < avHold.size() - n2; i++ )
				m_avKnots.get( i - n1 ).value = avHold.get( i ).value;
			
			int n = m_cvVerts;
			m_cvVerts -= n1 + n2;
			
			// Now process all the 'vertex' class variables.
			for ( CqParameter iUP : m_aUserParams )
			{
				if ( iUP.Class().getValue() == class_vertex )
				{
					CqParameter pHold = iUP.Clone();
					iUP.SetSize(  m_cvVerts * m_cuVerts );
					
					for ( int col = 0; col < m_cuVerts; col++ )
					{
						for ( i = n1; i < n - n2; i++ )
							iUP.SetValue( pHold, ((i - n1) * m_cuVerts)+ col, (i * m_cuVerts) + col );
					}
					pHold = null;
				}
			}
		}
	}
	
	/**
	 * Clamp the surface to ensure the knot vectors are 0-1 in each direction.
	 * 
	 * 
	 */
	public  void	Clamp()
	{
		ClampU(); ClampV();
	}
	
	/**
	 * 
	 * 
	 * 
	 */
	public  void	OutputMesh()
	{
		int Granularity = 30;  // Controls the number of steps in u and v
		
		//nttdata 高速化
		//STLVector<CqSurfaceNURBS> S = new STLVector<CqSurfaceNURBS>( CqSurfaceNURBS.class );
		STLArray<CqSurfaceNURBS> S = new STLArray<CqSurfaceNURBS>( CqSurfaceNURBS.class );
		S.add( this );
		
		try {
			// Save the grid as a .raw file.
			BufferedWriter bw = new BufferedWriter( new OutputStreamWriter( new FileOutputStream( "NURBS.RAW" )));
			
			int s;
			for ( s = 0; s < S.size(); s++ )
			{
				bw.write( "Surface_" + s + "\n" );
				STLVector<STLVector<CqVector3D>> aaPoints = new STLVector<STLVector<CqVector3D>>( 2, CqVector3D.class, Granularity + 1 );
				
				for ( int p = 0; p <= Granularity; p++ ) //二次元配列
					aaPoints.get( p ).resize( Granularity + 1 );
				
				
				// Compute points on curve
				int i;
				for ( i = 0; i <= Granularity; i++ )
				{
					float v = ((float) i / (float)Granularity )
					* ( S.get( s ).m_avKnots.get( S.get( s ).m_cvVerts ).value - S.get( s ).m_avKnots.get( S.get( s ).m_vOrder - 1 ).value )
					+ S.get( s ).m_avKnots.get( S.get( s ).m_vOrder - 1 ).value;
					
					for ( int j = 0; j <= Granularity; j++ ) {
						float u = ((float) j / (float)Granularity )
						* ( S.get( s ).m_auKnots.get( S.get( s ).m_cuVerts ).value - S.get( s ).m_auKnots.get( S.get( s ).m_uOrder - 1 ).value )
						+ S.get( s ).m_auKnots.get( S.get( s ).m_uOrder - 1 ).value;
						
						aaPoints.get( i ).get( j ).assignment( S.get( s ).Evaluate( u, v, P()) );
					}
				}
				
				for ( i = 0; i < Granularity; i++ ){
					for ( int j = 0; j < Granularity; j++ )
					{
						bw.write( aaPoints.get( i ).get( j ).x +" "+ aaPoints.get( i ).get( j ).y +" "+ aaPoints.get( i ).get( j ).z +" "+
								  aaPoints.get( i + 1 ).get( j + 1 ).x +" "+ aaPoints.get( i + 1 ).get( j + 1 ).y +" "+ aaPoints.get( i + 1 ).get( j + 1 ).z +" "+
								  aaPoints.get( i + 1 ).get( j ).x +" "+ aaPoints.get( i + 1 ).get( j ).y +" "+ aaPoints.get( i + 1 ).get( j ).z + "\n"
						);
						bw.write( aaPoints.get( i ).get( j ).x +" "+ aaPoints.get( i ).get( j ).y +" "+ aaPoints.get( i ).get( j ).z +" "+ 
								  aaPoints.get( i ).get( j + 1 ).x +" "+ aaPoints.get( i ).get( j + 1 ).y +" "+ aaPoints.get( i ).get( j + 1 ).z +" "+ 
								  aaPoints.get( i + 1 ).get( j + 1 ).x +" "+ aaPoints.get( i + 1 ).get( j + 1 ).y +" "+ aaPoints.get( i + 1 ).get( j + 1 ).z + "\n"
						);
					}
				}
			}
			bw.close();
		}catch ( IOException ioe ) { ioe.printStackTrace(); }
	}
	
	/**
	 * 
	 * 
	 * 
	 * @param name
	 * @param index
	 */
	public  void	AppendMesh( final String name, int index )
	{
		int Granularity = 10;  // Controls the number of steps in u and v
		
		try {
			// Save the grid as a .raw file.
			BufferedWriter bw = new BufferedWriter( new OutputStreamWriter( new FileOutputStream( name, true )));
			bw.write( "Surface_" + index + "\n" );
			
			STLVector<STLVector<CqVector3D>> aaPoints = new STLVector<STLVector<CqVector3D>>( 2, CqVector3D.class, Granularity + 1 );
			
			for ( int p = 0; p <= Granularity; p++ ) 
				aaPoints.get( p ).resize( Granularity + 1 );
			
			// Compute points on curve
			int i;
			for ( i = 0; i <= Granularity; i++ )
			{
				float v = ((float) i / (float) Granularity )
				* ( m_avKnots.get( m_cvVerts ).value - m_avKnots.get( m_vOrder - 1 ).value )
				+ m_avKnots.get( m_vOrder - 1 ).value;
				
				for ( int j = 0; j <= Granularity; j++ ){
					float u = ((float) j / (float) Granularity )
					* ( m_auKnots.get( m_cuVerts ).value - m_auKnots.get( m_uOrder - 1 ).value )
					+ m_auKnots.get( m_uOrder - 1 ).value;
					
					aaPoints.get( i ).get( j ).assignment( Evaluate( u, v, P()) );
				}
			}
			
			for ( i = 0; i < Granularity; i++ ){
				for ( int j = 0; j < Granularity; j++ )
				{
					bw.write( aaPoints.get( i ).get( j ).x +" "+ aaPoints.get( i ).get( j ).y +" "+ aaPoints.get( i ).get( j ).z +" "+
							aaPoints.get( i + 1 ).get( j + 1 ).x +" "+ aaPoints.get( i + 1 ).get( j + 1 ).y +" "+ aaPoints.get( i + 1 ).get( j + 1 ).z +" "+
							aaPoints.get( i + 1 ).get( j ).x +" "+ aaPoints.get( i + 1 ).get( j ).y +" "+ aaPoints.get( i + 1 ).get( j ).z + "\n"
					);
					bw.write( aaPoints.get( i ).get( j ).x +" "+ aaPoints.get( i ).get( j ).y +" "+ aaPoints.get( i ).get( j ).z +
							aaPoints.get( i ).get( j + 1 ).x +" "+ aaPoints.get( i ).get( j + 1 ).y +" "+ aaPoints.get( i ).get( j + 1 ).z + 
							aaPoints.get( i + 1 ).get( j + 1 ).x +" "+ aaPoints.get( i + 1 ).get( j + 1 ).y +" "+ aaPoints.get( i + 1 ).get( j + 1 ).z + "\n"
					);
				}
			}
			bw.close();
		}catch( IOException ioe ){ ioe.printStackTrace(); }
	}
	
	/**
	 * 
	 * 
	 * 
	 * @param name
	 */
	public  void Output( final String name )
	{
		int i;
		try {
			BufferedWriter bos = new BufferedWriter( new OutputStreamWriter( new FileOutputStream( name ))); 
			bos.write( "NuPatch " + m_cuVerts + " " + m_uOrder + "[\n" );
			
			for ( i = 0; i < m_auKnots.size(); i++ )
				bos.write( m_auKnots.get( i ).value + " \n" );
			bos.write( 0.0f + " " + 1.0f + " " + m_cvVerts + " " + m_vOrder + "[\n" );
			
			for ( i = 0; i < m_avKnots.size(); i++ )
				bos.write( m_avKnots.get( i ).value + " \n" );
			bos.write( "]\n" + 0.0f + " " + 1.0f + " " + "\"Pw\" [\n" );
			
			for( i = 0; i < P().Size(); i++ )
				bos.write( ((CqVector4D)P().pValue_get( i , 0 )).toString() );
			
			bos.close();
		}catch( IOException ioe ){ ioe.printStackTrace(); }
	}
	
	public CqTrimLoopArray TrimLoops()
	{
		return ( m_TrimLoops );
	}
	
	//****************************************************
	//***  Function from CqSurface
	//****************************************************
	
	/**
	 * Subdivide a bicubic patch in the u direction, return the left side.
	 * 
	 * 
	 * @param pnrbA
	 * @param pnrbB
	 */
	public  void uSubdivide( CqSurfaceNURBS pnrbA, CqSurfaceNURBS pnrbB )
	{
		pnrbA = new CqSurfaceNURBS();
		pnrbB = new CqSurfaceNURBS();
		
		SplitNURBS( pnrbA, pnrbB, true );
		
		uSubdivideUserParameters( pnrbA, pnrbB );
	}
	
	/**
	 * Subdivide a bicubic patch in the v direction, return the top side.
	 * 
	 * 
	 * @param pnrbA
	 * @param pnrbB
	 */
	public  void vSubdivide( CqSurfaceNURBS pnrbA, CqSurfaceNURBS pnrbB )
	{
		pnrbA = new CqSurfaceNURBS();
		pnrbB = new CqSurfaceNURBS();
		
		SplitNURBS( pnrbA, pnrbB, false );
		
		vSubdivideUserParameters( pnrbA, pnrbB );
	}
	
	/**
	 *  Dice the patch into a mesh of micropolygons.
	 * 
	 * 
	 * @param pParameter
	 * @param uDiceSize
	 * @param vDiceSize
	 */
	@SuppressWarnings("unchecked")
	public  void NaturalDice( CqParameter pParameter, int uDiceSize, int vDiceSize, IqShaderData pData )
	{
		//CqVector4D vec1 = new CqVector4D();
		int iv;
		for ( iv = 0; iv <= vDiceSize; iv++ )
		{
			float sv = ((float) iv / (float) vDiceSize )
			* ( m_avKnots.get( m_cvVerts ).value - m_avKnots.get( m_vOrder - 1 ).value )
			+ m_avKnots.get( m_vOrder - 1 ).value;
			int iu;
			for ( iu = 0; iu <= uDiceSize; iu++ )
			{
				int igrid = ( iv * ( uDiceSize + 1 ) ) + iu;
				float su = ((float)iu / (float) uDiceSize )
				* ( m_auKnots.get( m_cuVerts ).value - m_auKnots.get( m_uOrder - 1 ).value )
				+ m_auKnots.get( m_uOrder - 1 ).value;
				
				switch ( pParameter.Type().getValue() )
				{
				case type_float:
				{
					CqParameterTyped<p_float, p_float> pTParam = ( CqParameterTyped<p_float, p_float> ) pParameter ;
					pData.SetValue( ((p_float)Evaluate( su, sv, pTParam )).value, igrid );
					break;
				}
				case type_integer:
				{
					CqParameterTyped<p_int, p_float> pTParam = ( CqParameterTyped<p_int, p_float> ) pParameter ;
					pData.SetValue( ((p_int)Evaluate( su, sv, pTParam )).value, igrid );
					break;
				}
				case type_point:
				case type_normal:
				case type_vector:
				{
					CqParameterTyped<CqVector3D, CqVector3D> pTParam = ( CqParameterTyped<CqVector3D, CqVector3D> ) pParameter ;
					pData.SetValue( (CqVector3D)Evaluate( su, sv, pTParam ), igrid );
					break;
				}
				case type_hpoint: 
				{
					CqParameterTyped<CqVector4D, CqVector3D> pTParam = ( CqParameterTyped<CqVector4D, CqVector3D> ) pParameter ;
					pData.SetValue( new CqVector3D( Evaluate( su, sv, pTParam )), igrid );
					break;
				}
				case type_color:
				{
					CqParameterTyped<CqColor, CqColor> pTParam = ( CqParameterTyped<CqColor, CqColor> ) pParameter ;
					pData.SetValue( (CqColor)Evaluate( su, sv, pTParam ), igrid );
					break;
				}
				case type_string:
				{
					CqParameterTyped<p_String, p_String> pTParam = ( CqParameterTyped<p_String, p_String> ) pParameter ;
					pData.SetValue( ((p_String)Evaluate( su, sv, pTParam )).value, igrid );
					break;
				}
				case type_matrix:
				{
					CqParameterTyped<CqMatrix, CqMatrix> pTParam = ( CqParameterTyped<CqMatrix, CqMatrix> ) pParameter ;
					pData.SetValue( (CqMatrix)Evaluate( su, sv, pTParam ), igrid );
					break;
				}
				default:
				{
					// left blank to avoid compiler warnings about unhandled types
					break;
				}
				}
			}
		}
	}
	
	
	//virtual TqBool		CanGenerateNormals() const	{ return( TqTrue ); }
	
	
	/**
	 * Generate the vertex normals if not specified.
	 * 
	 * 
	 * @param uDiceSize
	 * @param vDiceSize
	 */
	public 	void	GenerateGeometricNormals( int uDiceSize, int vDiceSize, IqShaderData pNormals )
	{
		// Get the handedness of the coordinate system (at the time of creation) and
		// the coordinate system specified, to check for normal flipping.
		assert( P() != null );
		
		boolean CSO = pTransform().GetHandedness( pTransform().Time( 0 ));
		boolean O = pAttributes().GetIntegerAttribute( "System", "Orientation" ) [ 0 ] != 0;
		
		CqVector3D N = new CqVector3D();
		CqVector4D P = new CqVector4D();
		
		for ( int iv = 0; iv <= vDiceSize; iv++ )
		{
			float sv = ( (float)iv / (float) vDiceSize )
			* ( m_avKnots.get( m_cvVerts ).value - m_avKnots.get( m_vOrder - 1 ).value )
			+ m_avKnots.get( m_vOrder - 1 ).value;
			for ( int iu = 0; iu <= uDiceSize; iu++ )
			{
				float su = ( (float) iu / (float) uDiceSize )
				* ( m_auKnots.get( m_cuVerts ).value - m_auKnots.get( m_uOrder - 1 ).value )
				+ m_auKnots.get( m_uOrder - 1 ).value;
				int igrid = ( iv * ( uDiceSize + 1 )) + iu;
				N.assignment( EvaluateWithNormal( su, sv, P )); 
				N = ( O == CSO ) ? N : N.mul( -1f );
				pNormals.SetNormal( N, igrid );
				
				///TODO This would be more efficient if we can store the P here as well, instead of calculating it twice.
//				pP->SetPoint( P, igrid );
			}
		}
	}
	
	/**
	 * Return the boundary extents in object space of the surface patch
	 * @see net.cellcomputing.himawari.library.IqSurface#Bound()
	 */
	public 	CqBound	Bound() 
	{
		// Get the boundary in camera space.
		CqVector3D	vecA = new CqVector3D( MAX_VALUE, MAX_VALUE, MAX_VALUE );
		CqVector3D	vecB = new CqVector3D( -MAX_VALUE, -MAX_VALUE, -MAX_VALUE );
		
		int cVerts = m_cuVerts * m_cvVerts;
		for ( int i = 0; i < cVerts; i++ )
		{
			CqVector3D	vecV = new CqVector3D( (CqVector4D)P().pValue_get( i , 0 ) );
			if ( vecV.x < vecA.x ) vecA.x( vecV.x );
			if ( vecV.y < vecA.y ) vecA.y( vecV.y );
			if ( vecV.x > vecB.x ) vecB.x( vecV.x );
			if ( vecV.y > vecB.y ) vecB.y( vecV.y );
			if ( vecV.z < vecA.z ) vecA.z( vecV.z );
			if ( vecV.z > vecB.z ) vecB.z( vecV.z );
		}
		CqBound	B = new CqBound();
		B.vecMin().assignment( vecA );
		B.vecMax().assignment( vecB );
		
		return ( AdjustBoundForTransformationMotion( B ) );
	}
	
	/**
	 * Split the patch into smaller patches.
	 * 
	 * 
	 * @return
	 */
	public int Split( STLVector<CqBasicSurface> aSplits )
	{
		int cSplits = 0;
		
		if ( fPatchMesh() && ( cuSegments() > 1 || cvSegments() > 1 ) )
		{
			//nttdata 高速化
			//STLVector<CqSurfaceNURBS> S = new STLVector<CqSurfaceNURBS>( CqSurfaceNURBS.class );
			STLArray<CqSurfaceNURBS> S = new STLArray<CqSurfaceNURBS>( CqSurfaceNURBS.class );
			SubdivideSegments( S );
			int i;
			for ( i = 0; i < S.size(); i++ )
			{
				S.get( i ).SetSurfaceParameters( this );
				S.get( i ).TrimLoops().aLoops().addAll( TrimLoops().aLoops() );
				S.get( i ).m_fDiceable = true;
				S.get( i ).m_SplitDir.setValue( m_SplitDir.getValue());
				S.get( i ).m_EyeSplitCount = m_EyeSplitCount;
				//ADDREF( S[ i ] );
				aSplits.add( S.get( i ) );
			}
			return ( i );
		}
		
		// Split the surface in u or v
		CqSurfaceNURBS pNew1 = new CqSurfaceNURBS() ;
		CqSurfaceNURBS pNew2 = new CqSurfaceNURBS() ;
		
		// If this primitive is being split because it spans the e and hither planes, then
		// we should split in both directions to ensure we overcome the crossing.
		SplitNURBS( pNew1, pNew2, ( m_SplitDir.getValue() == SplitDir_U || !m_fDiceable ));
		
		for ( CqParameter iUP : m_aUserParams )
		{
			if ( iUP.Class().getValue() != class_vertex )
			{
				CqParameter pNewA = iUP.Clone();
				CqParameter pNewB = iUP.Clone();
				iUP.Subdivide( pNewA, pNewB, SplitDir() == SplitDir_U, this );
				pNew1.AddPrimitiveVariable( pNewA );
				pNew2.AddPrimitiveVariable( pNewB );
			}
		}
		
		pNew1.SetSurfaceParameters( this );
		pNew2.SetSurfaceParameters( this );
		pNew1.TrimLoops().aLoops().clear();
		pNew2.TrimLoops().aLoops().clear();
		pNew1.TrimLoops().aLoops().addAll( TrimLoops().aLoops() );
		pNew2.TrimLoops().aLoops().addAll( TrimLoops().aLoops() );
		pNew1.m_fDiceable = true;
		pNew2.m_fDiceable = true;
		pNew1.m_SplitDir.setValue( ( m_SplitDir.getValue() == SplitDir_U )? SplitDir_V : SplitDir_U );
		pNew2.m_SplitDir.setValue( ( m_SplitDir.getValue() == SplitDir_U )? SplitDir_V : SplitDir_U );
		pNew1.m_EyeSplitCount = m_EyeSplitCount;
		pNew2.m_EyeSplitCount = m_EyeSplitCount;
		pNew1.SetfPatchMesh( false );
		pNew2.SetfPatchMesh( false );
		
		aSplits.add( pNew1 );
		aSplits.add( pNew2 );
		
		cSplits = 2;
		
		if ( !m_fDiceable )
		{
			STLVector<CqBasicSurface> aSplits0 = new STLVector<CqBasicSurface>( CqBasicSurface.class );
			STLVector<CqBasicSurface> aSplits1 = new STLVector<CqBasicSurface>( CqBasicSurface.class );
			
			cSplits = aSplits.get( 0 ).Split( aSplits0 );
			cSplits += aSplits.get( 1 ).Split( aSplits1 );
			
			aSplits.clear();
			aSplits.addAll( aSplits0 );
			aSplits0.clear();
			aSplits.addAll( aSplits1 );
		}
		
		return ( cSplits );
	}
	
	/**
	 * Return whether or not the patch is diceable
	 * @see net.cellcomputing.himawari.library.CqBasicSurface#Diceable()
	 */
	public  boolean	Diceable()
	{
		assert( null != P() );
		
		// If the cull check showed that the primitive cannot be diced due to crossing the e and hither planes,
		// then we can return immediately.
		if ( !m_fDiceable )
			return ( false );
		
		// Otherwise we should continue to try to find the most advantageous split direction, OR the dice size.
		// Convert the control hull to raster space.
		CqVector2D[] avecHull = new CqVector2D[ m_cuVerts * m_cvVerts ];
		
		float gs = 16.0f;
		final float[] poptGridSize = QGetRenderContext().optCurrent().GetFloatOption( "System", "SqrtGridSize" );
		if( poptGridSize != null )
			gs = poptGridSize[0];
		
		float gridsize = 1.0f;
		if ( gs >= 1.0 ) 
			gridsize = gs * gs;
		
		float ShadingRateSqrt = pAttributes().GetFloatAttribute( "System", "ShadingRateSqrt" ) [ 0 ];
		final CqMatrix matCtoR = QGetRenderContext().matSpaceToSpace( "camera", "raster", new CqMatrix(), pTransform().matObjectToWorld(pTransform().Time(0)), QGetRenderContext().Time() );
		
		for ( int i = 0; i < m_cuVerts*m_cvVerts; i++ )
		{
			//nttdata 高速化
			//CqVector3D vT = new CqVector3D( (CqVector4D)P().pValue_get( i )[0] );
			CqVector3D vT = new CqVector3D( (CqVector4D)P().pValue_get( i , 0) );
			//vT.h(1.0f);
			vT = matCtoR.multiply( vT );
			avecHull[ i ] = new CqVector2D( vT );
		}
		
		// Now work out the longest continuous line in raster space for u and v.
		float uLen = 0f;
		float vLen = 0f;
		float MaxuLen = 0f;
		float MaxvLen = 0f;
		int v, u;
		
		for ( v = 0; v < m_cvVerts; v++ ){
			for ( u = 0; u < m_cuVerts - 1; u++ )
				uLen += ( avecHull[( v * m_cuVerts ) + u + 1].sub( avecHull[( v * m_cuVerts ) + u ]) ).Magnitude();
			
			if ( uLen > MaxuLen ) MaxuLen = uLen;
			uLen = 0;
		}
		for ( u = 0; u < m_cuVerts; u++ ) {
			for ( v = 0; v < m_cvVerts - 1; v++ )
				vLen += ( avecHull[(( v + 1 ) * m_cuVerts ) + u].sub( avecHull[( v * m_cuVerts ) + u ]) ).Magnitude();
					
			if ( vLen > MaxvLen ) MaxvLen = vLen;
			vLen = 0;
		}
		
		if ( MaxvLen > gridsize || MaxuLen > gridsize )
		{
			m_SplitDir.setValue(( MaxuLen > MaxvLen ) ? SplitDir_U : SplitDir_V );
			avecHull = null;
			return ( false );
		}
		
		if( ShadingRateSqrt>0 ){
			MaxuLen /= ShadingRateSqrt;
			MaxvLen /= ShadingRateSqrt;
		}
		else{
			MaxuLen = 0;
			MaxvLen = 0;
		}
		
		m_uDiceSize = max( round( MaxuLen ), 1 );
		m_vDiceSize = max( round( MaxvLen ), 1 );
		
		// Ensure power of 2 to avoid cracking
		final int[] binary = pAttributes().GetIntegerAttribute( "dice", "binary" );
		if ( binary != null && binary[0] != 0 )
		{
			m_uDiceSize = ( int )CEIL_POW2( m_uDiceSize );
			m_vDiceSize = ( int )CEIL_POW2( m_vDiceSize );
		}
		
		if ( MaxuLen < Float_h.FLT_EPSILON || MaxvLen < Float_h.FLT_EPSILON  ) 
		{
			m_fDiscard = true;
			avecHull = null;
			return ( false );
		}
		
		avecHull = null;
		m_SplitDir.setValue(( MaxuLen > MaxvLen ) ? SplitDir_U : SplitDir_V );
		
		if ( m_uDiceSize > gs) return false;
		if ( m_vDiceSize > gs) return false;
		
		return ( true );
	}
	
	/**
	 * Determine whether the passed surface is valid to be used as a
	 * frame in motion blur for this surface.
	 * 
	 * @return
	 */
	public boolean IsMotionBlurMatch( CqBasicSurface pSurf )
	{
		return( false );
	}
	
	/**
	 * 
	 * 
	 * 
	 * @param bUseDef_st
	 */
	public void SetDefaultPrimitiveVariables( boolean bUseDef_st )
	{
		int bUses = Uses();
		
		if ( USES( bUses, EnvVars_u ) )
		{
			AddPrimitiveVariable( new CqParameterTypedVarying<p_float, p_float>( "u", new EqVariableType( type_float ), p_float.class, p_float.class ));
			u().SetSize( cVarying() );
			
			float uinc = ( m_umax - m_umin ) / ( cuSegments() ); 
			 
			int c, r;
			int i = 0;
			for ( c = 0; c < cvSegments() + 1; c++ )
			{
				float uval = m_umin;
				for ( r = 0; r < cuSegments() + 1; r++ )
				{
					//nttdata 高速化
					//((p_float)u().pValue_get() [ i++ ]).value = uval;
					((p_float)u().pValue_get(0, i++ )).value = uval;
					uval += uinc;
				}
			}
		}
		
		if ( USES( bUses, EnvVars_v ) )
		{
			AddPrimitiveVariable( new CqParameterTypedVarying<p_float, p_float>( "v", new EqVariableType( type_float ), p_float.class, p_float.class ));
			v().SetSize( cVarying() );
			
			float vinc = ( m_vmax - m_vmin ) / ( cvSegments() );
			float vval = m_vmin;
			
			int c, r;
			int i = 0;
			for ( c = 0; c < cvSegments() + 1; c++ )
			{
				for ( r = 0; r < cuSegments() + 1; r++ )
					//nttdata 高速化
					//((p_float)v().pValue_get() [ i++ ] ).value = vval;
					((p_float)v().pValue_get(0, i++ ) ).value = vval;
				vval += vinc;
			}
		}
		
		final float[] pTC = pAttributes().GetFloatAttribute( "System", "TextureCoordinates" );
		CqVector2D st1 = new CqVector2D( pTC[ 0 ], pTC[ 1 ] );
		CqVector2D st2 = new CqVector2D( pTC[ 2 ], pTC[ 3 ] );
		CqVector2D st3 = new CqVector2D( pTC[ 4 ], pTC[ 5 ] );
		CqVector2D st4 = new CqVector2D( pTC[ 6 ], pTC[ 7 ] );
		
		if ( USES( bUses, EnvVars_s ) && !bHasVar(EnvVars_s) && bUseDef_st )
		{
			AddPrimitiveVariable( new CqParameterTypedVarying<p_float, p_float>( "s", new EqVariableType( type_float ), p_float.class, p_float.class ));
			s().SetSize( cVarying() );
			
			int c,r;
			int i = 0;
			for ( c = 0; c <= cvSegments(); c++ )
			{
				float v = ( 1.0f / ( cvSegments() ) ) * c;
				for ( r = 0; r <= cuSegments(); r++ )
				{
					float u = ( 1.0f / ( cuSegments() ) ) * r;
					((p_float) s().pValue_get( 0 , i++ )).value = ((p_float)BilinearEvaluate( new p_float(st1.x), new p_float(st2.x), new p_float(st3.x), new p_float(st4.x), u, v , p_float.class )).value;
				}
			}
		}
		
		if ( USES( bUses, EnvVars_t ) && !bHasVar(EnvVars_t) && bUseDef_st )
		{
			AddPrimitiveVariable( new CqParameterTypedVarying<p_float, p_float>( "t", new EqVariableType( type_float ), p_float.class, p_float.class ));
			t().SetSize( cVarying() );
			
			int c, r;
			int i = 0;
			for ( c = 0; c <= cvSegments(); c++ )
			{
				float v = ( 1.0f / ( cvSegments() ) ) * c;
				for ( r = 0; r <= cuSegments(); r++ )
				{
					float u = ( 1.0f / ( cuSegments() ) ) * r;
					((p_float) t().pValue_get( 0 , i++ )).value = ((p_float)BilinearEvaluate( new p_float(st1.y), new p_float(st2.y), new p_float(st3.y), new p_float(st4.y), u, v, p_float.class )).value;
				}
			}
		}
	}
	
	/* (non-Javadoc)
	 * @see net.cellcomputing.himawari.library.IqSurface#cUniform()
	 */
	public int cUniform() 
	{
		int nuSegments = ( 1 + m_cuVerts - m_uOrder );
		int nvSegments = ( 1 + m_cvVerts - m_vOrder );
		return ( nuSegments * nvSegments );
	}
	/* (non-Javadoc)
	 * @see net.cellcomputing.himawari.library.IqSurface#cVarying()
	 */
	public 	int cVarying()
	{
		int nuSegments = ( 1 + m_cuVerts - m_uOrder );
		int nvSegments = ( 1 + m_cvVerts - m_vOrder );
		return ( ( nuSegments + 1 ) * ( nvSegments + 1 ) );
	}
	/* (non-Javadoc)
	 * @see net.cellcomputing.himawari.library.IqSurface#cVertex()
	 */
	public 	int cVertex() 
	{
		return ( m_cuVerts * m_cvVerts );
	}
	/* (non-Javadoc)
	 * @see net.cellcomputing.himawari.library.IqSurface#cFaceVarying()
	 */
	public int cFaceVarying()
	{
		//TODO Must work out what this value should be.
		return ( 1 );
	}
	
	/* (non-Javadoc)
	 * @see net.cellcomputing.himawari.library.CqSurface#bCanBeTrimmed()
	 */
	public final boolean bCanBeTrimmed() 
	{
		return ( true );
	}
	
	/* (non-Javadoc)
	 * @see net.cellcomputing.himawari.library.CqSurface#bIsPointTrimmed(net.cellcomputing.himawari.library.types.CqVector2D)
	 */
	public final boolean bIsPointTrimmed( final CqVector2D p ) 
	{
		return  m_TrimLoops.TrimPoint( p );
	}
	
	/* (non-Javadoc)
	 * @see net.cellcomputing.himawari.library.CqSurface#bIsLineIntersecting(net.cellcomputing.himawari.library.types.CqVector2D, net.cellcomputing.himawari.library.types.CqVector2D)
	 */
	public final boolean bIsLineIntersecting( final CqVector2D v1, final CqVector2D v2 ) 
	{
		return ( m_TrimLoops.LineIntersects( v1, v2 ) );
	}
	
	/**
	 * Determine the segment count for the specified trim curve to make each segment the appropriate size
	 *  for the current shading rate.
	 * @see net.cellcomputing.himawari.library.CqSurface#TrimDecimation(net.cellcomputing.himawari.library.CqTrimCurve)
	 */
	public 	int	TrimDecimation( final CqTrimCurve Curve )
	{
		float Len = 0;
		float MaxLen = 0;
		int cSegments = 0;
		CqMatrix matCtoR = QGetRenderContext().matSpaceToSpace( "camera", "raster", new CqMatrix(), pTransform().matObjectToWorld( pTransform().Time( 0 )), QGetRenderContext().Time() );
		
		int iTrimCurvePoint;
		for ( iTrimCurvePoint = 0; iTrimCurvePoint < Curve.cVerts() - 1; iTrimCurvePoint++ )
		{
			// Get the u,v of the current point.
			CqVector3D vecCP = new CqVector3D();
			vecCP = Curve.CP( iTrimCurvePoint );
			float u = vecCP.x;
			float v = vecCP.y;
			
			// Get the u,v of the next point.
			vecCP = Curve.CP( iTrimCurvePoint + 1 );
			float u2 = vecCP.x;
			float v2 = vecCP.y;
			
			CqVector3D vecP = new CqVector3D( Evaluate( u, v, P() ));
			vecP = matCtoR.multiply( vecP ); 
			CqVector3D vecP2 = new CqVector3D( Evaluate( u2, v2, P() ));
			vecP2 = matCtoR.multiply( vecP2 ); 
			
			Len = ( vecP2.sub( vecP )).Magnitude();
			if ( Len > MaxLen ) MaxLen = Len;
			cSegments++;
		}
		float ShadingRateSqrt = pAttributes().GetFloatAttribute( "System", "ShadingRateSqrt" ) [ 0 ];
		MaxLen /= ShadingRateSqrt;
		
		int SplitCount = (int)max( MaxLen, 1 );
		
		return ( SplitCount * cSegments );
	}
	
	/* (non-Javadoc)
	 * @see net.cellcomputing.himawari.library.CqBasicSurface#PrepareTrimCurve()
	 */
	public 	void	PrepareTrimCurve()
	{
		m_TrimLoops.Prepare( this );
	}
	
}

