// 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 net.cellcomputing.himawari.accessory.STLVector;
import net.cellcomputing.himawari.accessory.primitive.p_float;
import net.cellcomputing.himawari.library.types.CqVector2D;
import net.cellcomputing.himawari.library.types.CqVector3D;

/**
 * Trim curve, based on NURBS surface, with control to restrict to 2d.
 * 	2dに制限するコントロールがあるNURBSの表面に基づいたトリム曲線。
 * 
 * @author NTT DATA Corporation
 */
public strictfp class CqTrimCurve {

	/** default constructer */
	public CqTrimCurve(){}
	
	/**　copy constructer */
	public CqTrimCurve( CqTrimCurve From )
	{
		this.assignment( From );
	}
	
	/** assignment */
	public void assignment( CqTrimCurve From )
	{
		int i;
		this.m_aKnots.resize( From.cKnots() );
		this.m_aVerts.resize( From.cVerts() );
		for( i = 0; i < From.cKnots(); i++ ) this.m_aKnots.get( i ).value = From.m_aKnots.get( i ).value;
		for( i = 0; i < From.cVerts(); i++ ) this.m_aVerts.get( i ).assignment( From.m_aVerts.get( i )) ;
		this.m_cVerts = From.m_cVerts;
		this.m_Order = From.m_Order;
	}
	
	/** destructer */
	public void destract(){}

	
    /**  NURBS表面のu方向の向きを取得する。 */
	public int	Order() 
    {
        return ( m_Order );
    }
	
    /**  NURBS表面のv方向の向きを取得する。 */
	public int	Degree() 
    {
        return ( m_Order -1 );
    }
	
    /**  u方向のコントロールポイントの数を取得する。*/
	public int	cVerts()
    {
        return ( m_cVerts );
    }
	
    /**  Knotのベクタのu方向の長さを取得する。 */
	public int	cKnots()
    {
        return ( m_cVerts + m_Order );
    }
	
    /**  Knotのベクタのu方向のへの参照を返す。*/
	public STLVector<p_float> aKnots()
    {
        return ( m_aKnots );
    }

    /** 
     * uとvの指定されたindexのコントロールポイントを取得する。
     * 
     * @param u　u方向のuのIndex 
     * @return 4Dの均質の制御点の参照。
     */
    public CqVector3D CP( final int u )
    {
        return ( m_aVerts.get( u ) );
    }
  
    /** 
     * NURBS構造を初期化して、規定寸法のNURBSカーブを取ってください。
     * 
     * @param Order 要求されたOrder
     * @param cVerts 必要な制御点カウント。
     */
    public void	Init( int Order, int cVerts )
    {
        int cKnots = cVerts + Order;
        m_aKnots.resize( cKnots );
        m_aVerts.resize( cVerts );
        m_Order = Order;
        m_cVerts = cVerts;
    }

    /**
     * NURBSのスパンを見つける。
     * 
     * @param u
     * @return　スパン
     */
    public int	FindSpan( float u )
    {
        if ( u >= m_aKnots.get( m_cVerts ).value )
            return ( m_cVerts -1 );
        if ( u <= m_aKnots.get( Degree() ).value )
            return ( Degree() );

        int low = 0;
        int high = m_cVerts + 1;
        int mid = ( low + high ) / 2;

        while ( u < m_aKnots.get( mid ).value || u >= m_aKnots.get( mid + 1 ).value )
        {
            if ( u < m_aKnots.get( mid ).value )
                high = mid;
            else
                low = mid;
            mid = ( low + high ) / 2;
        }
        return ( mid ); 
    }
    
    /**
     * 基本機能
     * 与えられた点データから、NURBS曲線を生成する。
     * 
     * @param u
     * @param span
     * @param BasisVals
     */
    public void	BasisFunctions( float u, int span, float[] BasisVals )
    {
        int r, s, i;
        double omega;

        BasisVals[ 0 ] = 1.0f;
        for ( r = 2; r <= (int)( m_Order ); r++ )
        {
            i = span - r + 1;
            BasisVals[ r - 1 ] = 0.0f;
            for ( s = r - 2; s >= 0; s-- )
            {
                i++;
                if ( i < 0 )
                    omega = 0;
                else
                    omega = ( u - m_aKnots.get( i ).value ) / ( m_aKnots.get( i + r - 1 ).value - m_aKnots.get( i ).value );

                BasisVals[ s + 1 ] = (float)(BasisVals[ s + 1 ] + ( 1 - omega ) * BasisVals[ s ] );
                BasisVals[ s ] = (float)(omega * BasisVals[ s ] ) ;
            }
        }
    }
    
    /**
     * 評価機能と思われます。
     * @param u　方向性を決める。
     */
    public CqVector2D	Evaluate( float u )
    {
    	float basis[] = new float[m_Order];
        CqVector3D r = new CqVector3D( 0, 0, 0 );

        // Evaluate non-uniform basis functions (and derivatives)

        int span = FindSpan( u );
        int first = span - m_Order + 1;
        BasisFunctions( u, span, basis );

        // Weight control points against the basis functions
        int j;
        for ( j = 0; j < m_Order; j++ )
        {
            int rj = m_Order - 1 - j;
            
            float tmp = basis[ rj ];
            CqVector3D cp = CP( j + first );

            r.x =  r.x + cp.x * tmp ;
            r.y =  r.y + cp.y * tmp ;
            r.z =  r.z + cp.z * tmp ;
        }

        return ( new CqVector2D( r.x / r.z, r.y / r.z ) );
    }
    
    /**
     * 
     * 
     * 
     * @param u
     * @param r
     * @return
     */
    public int	InsertKnot( float u, int r )
    {
        // Work on a copy.
        CqTrimCurve nS = new CqTrimCurve( this );

        // Compute k and s      u = [ u_k , u_k+1)  with u_k having multiplicity s
        int k = m_aKnots.size() - 1, s = 0;
        int i, j;
        int p = Degree();

        if ( u < m_aKnots.get( Degree() ).value || u > m_aKnots.get( m_cVerts ).value )
            return ( 0 );

        int size = (int)( m_aKnots.size() );
        for ( i = 0; i < size; i++ )
        {
            if ( m_aKnots.get( i ).value > u )
            {
                k = i - 1;
                break;
            }
        }

        if ( u <= m_aKnots.get( k ).value )
        {
            s = 1;
            for ( i = k; i > 0; i-- )
            {
                if ( m_aKnots.get( i ).value <= m_aKnots.get( i - 1 ).value )
                    s++;
                else
                    break;
            }
        }
        else
            s = 0;

        if ( ( r + s ) > p + 1 )
            r = p + 1 - s;

        if ( r <= 0 )
            return ( 0 );

        nS.Init( m_Order, m_cVerts + r );

        // Load new knot vector
        for ( i = 0;i <= k;i++ ) 
        	nS.m_aKnots.get( i ).value = m_aKnots.get( i ).value;
        for ( i = 1;i <= r;i++ ) 
        	nS.m_aKnots.get( k + i ).value = u ;
        
        size = (int)( m_aKnots.size() );
        for ( i = k + 1;i < size; i++ )
            nS.m_aKnots.get( i + r ).value = m_aKnots.get( i ).value ;

		// Save unaltered control points 
        CqVector3D R[] = new CqVector3D[ p + 1 ];
        // Insert control points as required on each row.
        for ( i = 0; i <= k - p; i++ ) 
        	nS.CP( i ).assignment( CP( i ) );
        size = (int)( m_cVerts );
        for ( i = k - s; i < size; i++ )
            nS.CP( i + r ).assignment( CP( i ) );
        for ( i = 0; i <= p - s; i++ ) 
        	R [ i ] = CP( k - p + i ) ;

        // Insert the knot r times
        int L = 0 ;
        float alpha;
        for ( j = 1; j <= r; j++ )
        {
            L = k - p + j;
            int i1;
            int limit = p - j - s;
            for ( i = 0;i <= limit;i++ )
            {
                i1 = i + 1;
                alpha = ( u - m_aKnots.get( L + i ).value ) / ( m_aKnots.get( k + i1 ).value - m_aKnots.get( L + i ).value );
                R[ i ] = ( new CqVector3D( (float)(alpha * R[ i1 ].x + ( 1.0 - alpha ) * R[ i ].x),
                                     (float)(alpha * R[ i1 ].y + ( 1.0 - alpha ) * R[ i ].y),
                                     (float)(alpha * R[ i1 ].z + ( 1.0 - alpha ) * R[ i ].z ) ));
            }
            nS.CP( L ).assignment( R[ 0 ] );
            if ( p - j - s > 0 )
                nS.CP( k + r - j - s ).assignment( R[ p - j - s ] );
        }

        // Load remaining control points
        for ( i = L + 1; i < k - s; i++ )
            nS.CP( i ).assignment( R[ i - L ]);

        this.assignment( nS );

        return ( r );
    }
    
    /**
     * Ensure a nonperiodic (clamped) knot vector 
     * by inserting U[p] and U[m-p] multiple times.
     * 
     */
    public void	Clamp()
    {
        int n1 = InsertKnot( m_aKnots.get( Degree() ).value, Degree() );
        int n2 = InsertKnot( m_aKnots.get( m_cVerts ).value, Degree() );

        // Now trim unnecessary knots and control points
        if ( n1 != 0 || n2 != 0 )
        {
            CqTrimCurve nS = new CqTrimCurve( this );
            m_aKnots.resize( m_aKnots.size() - n1 - n2 );
            m_aVerts.resize( m_cVerts - n1 - n2 );
            m_cVerts -= n1 + n2;
            
            int i;
            int size = nS.m_aKnots.size() - n2;
            
            for ( i = n1; i < (int)( size ); i++ )
                m_aKnots.get( i - n1 ).value = nS.m_aKnots.get( i ).value ;
            size = nS.m_cVerts - n2;
            
            for ( i = n1; i < (int)( size ); i++ )
                CP( i - n1 ).assignment( nS.CP( i ) );
       
        }
    }

    protected STLVector<p_float>	m_aKnots = new STLVector<p_float>( p_float.class );	///< Knot vector.
    protected int	m_Order;	///< Surface order.
    protected int	m_cVerts;	///< Control point.
    protected STLVector<CqVector3D>	m_aVerts = new STLVector<CqVector3D>( CqVector3D.class );	///< Nurbs control points.

}
