/*=========================================================================

 Program:   wxMaracas
 Module:    $RCSfile: marContour.cpp,v $
 Language:  C++
 Date:      $Date: 2008/10/31 16:32:54 $
 Version:   $Revision: 1.1 $
 
  Copyright: (c) 2002, 2003
  License:
  
   This software is distributed WITHOUT ANY WARRANTY; without even 
   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
   PURPOSE.  See the above copyright notice for more information.
   
=========================================================================*/



extern "C" {
#include"idcnt.h"
};

//#include <gsl/gsl_blas.h> //PS
#include "marVector.h"
#include <math.h> //PS
#include "marMathConst.h" //PS

#include "marContour.h"
#include <vtkPolyLine.h>
#include <vtkKochanekSpline.h> 

#include "load_snake.h"

// -------------------------------------------------------------------------
marContour::marContour( double s, marParameters* p )
: marObject( p ),
_s( s )
{
}

// -------------------------------------------------------------------------

void marContour::addContourPoint(  double x, double y  )
{
    double* np = new double[ 2 ];
	
//    memcpy( np, p, 2 * sizeof( double ) );
	np[0]=x;
	np[1]=y;
    _points.push_back( np );
}

// -------------------------------------------------------------------------
void marContour::resample( )
{ // TODO
}

// -------------------------------------------------------------------------
bool marContour::isCalculated( )
{
    return( _points.size( ) > 0 );
}

// -------------------------------------------------------------------------
int marContour::getNumberOfPoints( )
{
    return( _points.size( ) );
}

// -------------------------------------------------------------------------
void marContour::getPoint( double* p, int i )
{
    memcpy( p, _points[ i ], 2 * sizeof( double ) );
}
// -------------------------------------------------------------------------
void marContour::calculateVariables()
{

	 _area					=	-1;
	 _perimeter				=	-1;
	 _diameterFromArea		=	-1;
	 _diameterFromPerimeter	=	-1;
	 _minimumDiameter		=	-1;
	 _maximumDiameter		=	-1;
	 _averageDiameter		=	-1;


	calculateArea();
	calculatePerimeter();
	calculateDiameterFromArea();	
	calculateDiameterFromPerimeter();
	calculateMinimumMaximumDiameter();
	calculateAverageDiameter();	

}
/* METHOD DESCRIPTION ****************************************************
 *                                                                       *
 * marContour::do_spline (method)                                        *
 *                                                                       *
 * DESCRIPTION : Calculates a spline function in order to add points to  *
 *               actual points.                                          *
 *                                                                       *
 * SYNTAX : ctn.do_spline( )                                             *
 *                                                                       *
 ******************************************************* END DESCRIPTION */
void marContour::do_spline( ) {
    unsigned int i, to;
    vtkKochanekSpline* CntSplineX = vtkKochanekSpline::New( );
    vtkKochanekSpline* CntSplineY = vtkKochanekSpline::New( );

//EED	temPoint< float, 2 >* tmp;
	double *tmp;

//EED    float t;
    unsigned int nps = 1000; // TODO: 1000?????
    CntSplineX->SetDefaultTension( 0 ); CntSplineX->SetDefaultBias( 0 ); CntSplineX->SetDefaultContinuity( 0 );
    CntSplineY->SetDefaultTension( 0 ); CntSplineY->SetDefaultBias( 0 ); CntSplineY->SetDefaultContinuity( 0 );
    for( i = 0; i < _points.size( ); i++ ) {
        CntSplineX->AddPoint( i, ( _points[ i ] )[ 0 ] );
        CntSplineY->AddPoint( i, ( _points[ i ] )[ 1 ] );
    } //  rof

    if( i > 0 ) {
        CntSplineX->AddPoint( i, ( _points[ 0 ] )[ 0 ] );
        CntSplineY->AddPoint( i, ( _points[ 0 ] )[ 1 ] );
    } // fi
    to = _points.size( );
    while( _no_spline_points.size( ) > 0 ) {
        delete _no_spline_points[ _no_spline_points.size( ) - 1 ];
        _no_spline_points.pop_back( );
    } // fwhile
    for( i = 0; i < _points.size( ); i++ ) {
//EED         tmp = new temPoint< float, 2 >;
		tmp = (double*)malloc ( 2*sizeof(double) );
         ( tmp )[ 0 ] = ( _points[ i ] )[ 0 ];
         ( tmp )[ 1 ] = ( _points[ i ] )[ 1 ];
        _no_spline_points.push_back( tmp );
    } // rof

//EED    this->clean( );
//EED    for( i = 0; i < nps; i++ ) {
//EED        tmp = new temPoint< float, 2 >;
//EED        t = ( float ) ( to - 1 ) / ( float ) ( nps - 1 ) * ( float ) i;
//EED        ( *tmp )[ 0 ] = CntSplineX->Evaluate( t );
//EED        ( *tmp )[ 1 ] = CntSplineY->Evaluate( t );
//EED        _points.push_back( tmp );
//EED    } // rof
//EED    _updated_measures = false;

	CntSplineX->Delete();
	CntSplineY->Delete();


}

// -------------------------------------------------------------------------
void marContour::calculateArea( )
{
    double area;
    int i, j;
	
    // This uses Green's theorem:
    // A = 1/2 * sum( xiyi+1 - xi+1yi); pO == pN
    // A < 0 -> A = |A| (a negative value could raise because points are
    // given in clockwise order).
    for( i = 0, area = 0.0; i < _points.size( ); i++ ) {
		
        j = ( i + 1 ) % _points.size( );
		
        //  Area
        area +=
            ( _points[ i ][ 0 ] * _points[ j ][ 1 ] ) -
            ( _points[ j ][ 0 ] * _points[ i ][ 1 ] );
		

    } // rof
    area /= 2.0;
    area = fabs( area );

    this->_area = area;
}








// -------------------------------------------------------------------------
void marContour::calculatePerimeter( )
{
	// PS ->     //gsl_vector* pO = gsl_vector_alloc( 2 ); //PS
	marVector* pO = new marVector(2);
	// PS ->     //gsl_vector* pF = gsl_vector_alloc( 2 ); //PS
	marVector* pF = new marVector(2);
    double L;
    unsigned int j;
	
	int nbPoints=_points.size( );

    for( j = 0, L = 0.0; j < nbPoints; j++ ) {
		
		// PS -> 	//memcpy( pO->data, _points[ j ], 2 * sizeof( double ) ); //PS
		(*pO)=_points[j];
		// PS -> 	//memcpy( pF->data, _points[ j + 1 ], 2 * sizeof( double ) ); //PS
		(*pF)=_points[(j+1)%nbPoints];
		
		// PS -> 	// gsl_vector_sub( pF, pO );//PS
		(*pF)=(*pF)-(*pO);
		// PS -> 	//L += gsl_blas_dnrm2( pF );//PS
		L+=pF->norm2();
		
    } // rof
	/*
	// PS ->     //memcpy( pO->data, _points[ j - 1 ], 2 * sizeof( double ) ); //PS
     (*pO)=_points[j-1];
	// PS -> 	//memcpy( pF->data, _points[ 0 ], 2 * sizeof( double ) ); //PS
 	(*pF)=_points[0];
	
	
	// PS -> 	// gsl_vector_sub( pF, pO );//PS
 	(*pF)=(*pF)-(*pO);
	// PS -> 	//L += gsl_blas_dnrm2( pF );//PS
 	L+=pF->norm2();
	
	// PS ->     //gsl_vector_free( pO );//PS
	// PS ->     //gsl_vector_free( pF );//PS
	*/
	delete pO;
	delete pF;

    _perimeter=L ;
}

// -------------------------------------------------------------------------
/* METHOD DESCRIPTION ****************************************************
 *                                                                       *
 * marContour::recalcule_measures (method)                               *
 *                                                                       *
 * DESCRIPTION : Does one step to calcule measures. Don't call this      *
 *               outside this file.                                      *
 *                                                                       *
 * SYNTAX : this->recalcule_measures( )                                  *
 *                                                                       *
 ******************************************************* END DESCRIPTION */
void marContour::calculateMinimumMaximumDiameter( )
{

	double _minimum_diameter;
	double _maximum_diameter;
	double _p1_min[2];
	double _p2_min[2];
	double _p1_max[2];
	double _p2_max[2];
//-----

    unsigned int i, j;
    float  tmp1[2], tmp2[2];
    vtkPoints* points;
    vtkPolyLine* poly;
    vtkUnstructuredGrid* grid;
    PCONTOUR_FLOAT ctr, ctr2;

    double mom_min;
    double mom_max;

//EED    if( !_updated_measures ) {
        //  Preparation
//EED        _area = 0.0;
//EED        _perimeter = 0.0;
        _minimum_diameter = 0.0;
        _maximum_diameter = 0.0;
        _p1_min[ 0 ] = _p1_min[ 1 ] = 0.0;
        _p2_min[ 0 ] = _p2_min[ 1 ] = 0.0;
        _p1_max[ 0 ] = _p1_max[ 1 ] = 0.0;
        _p2_max[ 0 ] = _p2_max[ 1 ] = 0.0;
//EED        //  This uses Green's theorem:
//EED        //  A = 1/2 * sum( xiyi+1 - xi+1yi); pO == pN
//EED        //  A < 0 -> A = |A| (a negative value is because points are given in
//EED        //  clockwise order).
//EED        //
//EED        //  P = sum( |Pj_Pj+1| ); pO == pN
//EED        //
//EED        //  Prepares data for maximum & minimum process.
//EED        //  The scale factor must be applied.
        points = vtkPoints::New( );
        poly   = vtkPolyLine::New( );
        grid   = vtkUnstructuredGrid::New( );

        poly->GetPointIds( )->SetNumberOfIds( _points.size( ) );
        ctr2 = ( PCONTOUR_FLOAT ) IdCntAlloc( _points.size( ), CNT_FLOAT );

        for( i = 0, j = 1; i < _points.size( ); i++, j = ( j + 1 ) % _points.size( ) ) {
            tmp1[0] = _points[ i ][0];
            tmp1[1] = _points[ i ][1];
//EED            tmp1 = *_points[ i ];
            tmp2[0] = _points[ j ][0];
            tmp2[1] = _points[ j ][1];
//EED            tmp2 = *_points[ j ];
            //  VTK stuff for maximum & minimum
            poly->GetPointIds( )->SetId( i, i );
            points->InsertNextPoint( tmp1[ 0 ], tmp1[ 1 ], 1 );
            //  Point scale factor
//EED           tmp1 *= ( _parameters->get_dim_ima( ) / ( float ) _parameters->get_size_ima( ) );
//EED           tmp2 *= ( _parameters->get_dim_ima( ) / ( float ) _parameters->get_size_ima( ) );
            //  Area & perimeter
//EED            IdCntAddPointG( ctr2, tmp1[ 0 ], tmp1[ 1 ] );
//EED            _area += ( ( tmp1[ 0 ] * tmp2[ 1 ] ) - ( tmp2[ 0 ] * tmp1[ 1 ] ) );
//EED            _perimeter += tmp1 * tmp2;
        } // rof
//EED        _area *= 0.5;
//EED        _area = ( float ) fabs( ( double ) _area );
        //  TODO : Just testing: LibIDO is better than green theorem's?
        //  which is the difference?
//EED        _area = ( float ) IdCntCalculAireG( ( PCONTOUR ) ctr2 );
        //  Libido contour
        ctr = ( PCONTOUR_FLOAT ) IdCntAlloc( _no_spline_points.size( ), CNT_FLOAT );
        for( i = 0; i < _no_spline_points.size( ); i++ ) {
            IdCntAddPointG(
                ctr,
                ( _no_spline_points[ i ] )[ 0 ],
                ( _no_spline_points[ i ] )[ 1 ]
            );
        } // rof
        //  Diameters...

        if( _points.size( ) > 0 ) {
            grid->Allocate( 1, 1 );
            grid->InsertNextCell( poly->GetCellType( ), poly->GetPointIds( ) );
            grid->SetPoints( points );
            axe_max_min(
                ctr,
                _no_spline_points.size( ),
                ( int ) grid->GetCell( 0 ),
                &_minimum_diameter, &_maximum_diameter,
                &( _p1_min[ 0 ] ), &( _p1_min[ 1 ] ),
                &( _p2_min[ 0 ] ), &( _p2_min[ 1 ] ),
                &( _p1_max[ 0 ] ), &( _p1_max[ 1 ] ),
                &( _p2_max[ 0 ] ), &( _p2_max[ 1 ] ),
                &mom_min, &mom_max
            );
			_min[0]=_p1_min[ 0 ]; 
			_min[1]=_p1_min[ 1 ]; 
			_min[2]=_p2_min[ 0 ]; 
			_min[3]=_p2_min[ 1 ]; 
			_max[0]=_p1_max[ 0 ]; 
			_max[1]=_p1_max[ 1 ]; 
			_max[2]=_p2_max[ 0 ]; 
			_max[3]=_p2_max[ 1 ]; 

//EED
//			float param_dim  =(float) _parameters->get_dim_ima( );
//			float param_size =(float) _parameters->get_size_ima( );
//			float relation=param_dim / param_size
//            _minimum_diameter *= relation;
//            _maximum_diameter *= relation;

        } // fi
//EED     _updated_measures = true;
//EED   } // fi
//}
//-----
	_minimumDiameter = _minimum_diameter;
	_maximumDiameter = _maximum_diameter;

	points->Delete();
	poly->Delete();
	grid->Delete();
	
}


// -------------------------------------------------------------------------
void marContour::getMaximumLine( double* pO, double* pF )
{
    memcpy( pO, _max, 2 * sizeof( double ) );
    memcpy( pF, _max + 2, 2 * sizeof( double ) );
}

// -------------------------------------------------------------------------
void marContour::getMinimumLine( double* pO, double* pF )
{
    memcpy( pO, _min, 2 * sizeof( double ) );
    memcpy( pF, _min + 2, 2 * sizeof( double ) );
}

// -------------------------------------------------------------------------
void marContour::reset( )
{
    int i;
	
    for( i = 0; i < _points.size( ); i++ )
		delete _points[ i ];
    _points.clear( );

    for( i = 0; i < _no_spline_points.size( ); i++ )
		delete _no_spline_points[ i ];
    _no_spline_points.clear( );

    _min[ 0 ] = _min[ 1 ] = _min[ 2 ] = _min[ 3 ] = 0.0;
    _max[ 0 ] = _max[ 1 ] = _max[ 2 ] = _max[ 3 ] = 0.0;
}

// -------------------------------------------------------------------------
void marContour::copyFrom( const marObject& from )
{
    double* tmp;
    int i;
    marContour* f = &( ( marContour& )from );
	
    reset( );
    memcpy( _min, f->_min, 4 * sizeof( double ) );
    memcpy( _max, f->_max, 4 * sizeof( double ) );
    for( i = 0; i < f->_points.size( ); i++ ) {
		
		tmp = new double[ 3 ];
		memcpy( tmp, f->_points[ i ], 3 * sizeof( double ) );
		_points.push_back( tmp );
		
    } // rof
}

// -------------------------------------------------------------------------
bool marContour::save( std::ofstream& os )
{
    int i;
	
    i = _points.size( );
    os.write( ( const char* )&i, sizeof( int ) );
    for( i = 0; i < _points.size( ); i++ )
		os.write( ( const char* )_points[ i ], 3 * sizeof( double ) );
	
    return( true );
}

// -------------------------------------------------------------------------
bool marContour::load( std::ifstream& is )
{
    int i, size;
    double* tmp;
	
    reset( );
    is.read( ( char* )&size, sizeof( int ) );
    for( i = 0; i < size; i++ ) {
		
		tmp = new double[ 3 ];
		is.read( ( char* )tmp, 3 * sizeof( double ) );
		_points.push_back( tmp );
		
    } // rof
	
    return( true );
}


// ----------------------------------------------------------------------------
vtkUnstructuredGrid* marContour::Draw( )
{
	int taille;
	double p[2];
	
	vtkPoints* PointsContour;
	vtkPolyLine* aPolyLine;
	int i;
	
	taille = (int) (this->getNumberOfPoints( ) / 2);
	
	PointsContour = vtkPoints::New( );
	
	aPolyLine = vtkPolyLine::New( );
	( aPolyLine->GetPointIds( ) )->SetNumberOfIds( taille );
	for ( i = 0;  i < taille; i++) {
		( aPolyLine->GetPointIds( ) )->SetId( i , i);
		getPoint((double *) p, i*2);
		PointsContour->InsertNextPoint( p[ 0 ], p[ 1 ], 0 );
	}
	
	_allData = vtkUnstructuredGrid::New( );
	_allData->Allocate( 1 , 1 );
	_allData->InsertNextCell( aPolyLine->GetCellType( ) , 
		aPolyLine->GetPointIds( ) );
	_allData->SetPoints( PointsContour );
	
	PointsContour->Delete( );
	aPolyLine->Delete( );
	return ( _allData );
}

// ----------------------------------------------------------------------------
void marContour::Delete( )
{
	if(_allData) {
		_allData->Delete();
		_allData = NULL;
	}
	
	reset();
}


// eof - contour.cxx