marAxis.cpp

Go to the documentation of this file.

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

  Program:   wxMaracas
  Module:    $RCSfile: marAxis.cpp,v $
  Language:  C++
  Date:      $Date: 2009/05/14 13:55:07 $
  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.

=========================================================================*/
#include <vtkPlaneSource.h>
#include <vtkKochanekSpline.h>
#include <vtkPolyDataConnectivityFilter.h>         
#include <vtkContourFilter.h>
#include <vtkProbeFilter.h>


#include <vtkImageChangeInformation.h>
#include <vtkImageShiftScale.h> 
#include <vtkDoubleArray.h> 


#include "marAxis.h"
#include "marVector.h"

#include "ExtractionAxe.h"
#include "FonctionsGrales.h"

#include <vtkCellArray.h>
#include <vtkCell.h>
#include <vtkCleanPolyData.h>
#include <vtkIdList.h>
#include <vtkMath.h>
#include <vtkMatrix4x4.h> 
#include <vtkPointData.h>
#include <vtkPolyData.h>
#include <vtkPlane.h>
#include <vtkSetGet.h>
#include <vtkStripper.h>
#include <vtkStructuredPoints.h>
#include <vtkTransform.h> 

// ----------------------------------------------------------------------------
marAxis::marAxis( marParameters* p ) : 
       marObject( p ), 
           kCurve( 3, INDX_count - 3 ),
      _healthySlice( -1 ), _startQuant( -1 ),
      _finishQuant( -1 ), _actualQuant( -1 )

{
        _totalAxisLenght=-1;
        _subAxisLenght=-1;
        _referenceArea=-1;
        _referenceAverDiam=-1;


    _allData = NULL;

        _healthySliceStart=-1;
        _healthySliceEnd=-1;
    reset( );
        calibration = false;
}

// ----------------------------------------------------------------------------
void marAxis::addAxisPoint( double* p )
{
    double pt[ INDX_count ];

    memcpy( pt, p, sizeof( POINTAXE ) );
    pt[ INDX_SIGNALVALUE ] = 0;
    addControlPoint( pt, pt + 3 );
}

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

double* marAxis::getNormal( unsigned int slice ) {

    double* ret = new double[ 3 ];

    if ( slice == 0 ) {
      ret[0] = _points[ slice + 1 ][0] - _points[ slice ][0];
      ret[1] = _points[ slice + 1 ][1] - _points[ slice ][1];
      ret[2] = _points[ slice + 1 ][2] - _points[ slice ][2];
    }
    else if ( slice == _points.size( ) - 1 ) {
      ret[0] = _points[ slice - 1 ][0] - _points[ slice ][0];
      ret[1] = _points[ slice - 1 ][1] - _points[ slice ][1];
      ret[2] = _points[ slice - 1 ][2] - _points[ slice ][2];
    }
    else if ( slice > 0 && slice < _points.size( ) - 1 ) {
      ret[0] = _points[ slice - 1 ][0] - _points[ slice + 1 ][0];
      ret[1] = _points[ slice - 1 ][1] - _points[ slice + 1 ][1];
      ret[2] = _points[ slice - 1 ][2] - _points[ slice + 1 ][2];
    }

    return( ret );
}

// ----------------------------------------------------------------------------
void marAxis::changeAxisResolution( )
{/*
   int i;
   int nOut = ( int )ceil(
   getParameters( )->
   getDoubleParam( marParameters::e_axis_discret_step )
   * ( double )( GetNumberOfPoints( ) )
   );

   i = GetNumberOfPoints( );
   double d = getParameters( )->
   getDoubleParam( marParameters::e_axis_discret_step );

   SetResolution( nOut );
   for( i = 0; i < _contours.size( ); i++ )
   delete _contours[ i ];
   _contours.clear( );
   for( i = 0; i < nOut; i++ )
   _contours.push_back( new marContour( getParameters( ) ) );
 */
}

// ----------------------------------------------------------------------------
void marAxis::calculateSignal( kVolume* vol )
{
    /* PSIGNAL_FLOAT PutSignalValueFloat(PSIGNAL_FLOAT signal, int pos, float value)
    {
      signal[pos] = value;
      return signal;
    }*/
  ushort maxlocal;
  double *p;

  double rayon;
  unsigned int i, mask_size = getParameters( )-> getIntParam( marParameters::e_mask_size );

  while( _signal.size( ) > 0 ) _signal.pop_back( );
  // Intensity signal

  PSIGNAL_USHORT signal = ( PSIGNAL_USHORT ) IdSigAlloc( _points.size( ), SIG_USHORT );
  for( i = 0; i < _points.size( ); i++ ) {
    p = _points[i];
    rayon = RayonLocal(
                      ( int ) ( p[ 0 ] ),
                      ( int ) ( p[ 1 ] ),
                      ( int ) ( p[ 2 ] ),
                      _points_disc, getNumberOfControlPoints( ));

    maxlocal = vol->GetMaxIntSphere2( p, rayon );
    PutSignalValue( signal, i, maxlocal );
  } // rof

  PSIGNAL_FLOAT mask_signal = ( PSIGNAL_FLOAT )IdSigAlloc( mask_size, SIG_FLOAT );

  for( i = 0; i < mask_size; i++ )
    PutSignalValueFloat( mask_signal, i, 1 );

  PSIGNAL_USHORT signalfilt = ( PSIGNAL_USHORT )IdSigAlloc( _points.size( ), SIG_USHORT );
  signalfilt = MonSigConvolve( signal, mask_signal, signalfilt, 1.0 / ( double )mask_size, 0 );
  for( i = 0; i < _points.size( ); i++ ){
    _signal.push_back( GetSignalValue( signalfilt, i ) );
  }
  IdSigFree( signal );
  IdSigFree( mask_signal );
  IdSigFree( signalfilt );
}

// ----------------------------------------------------------------------------
void marAxis::start( )
{
    /* TODO
       int swap, n = GetNumberOfPoints( );

       swap = GSL_MIN( _startQuant, _finishQuant );
       _finishQuant = GSL_MAX( _startQuant, _finishQuant );
       _startQuant = swap;

       _startQuant = ( _startQuant >= 0 )? _startQuant: 0;
       _finishQuant = ( _finishQuant >= n )? _finishQuant: n - 1;
       _actualQuant = _startQuant;
    */
}

// ----------------------------------------------------------------------------
void marAxis::cut( int slice, bool up )
{ // TODO
}


// ----------------------------------------------------------------------------
void marAxis::doSpline( )
{
  unsigned int i;
  unsigned int nop, nip;
  float length=0, t;
  float p1[ 3 ], p2[ 3 ];
  double p[marAxis::INDX_count];

  vtkKochanekSpline* spX = vtkKochanekSpline::New( );
  vtkKochanekSpline* spY = vtkKochanekSpline::New( );
  vtkKochanekSpline* spZ = vtkKochanekSpline::New( );

  for( i = 0, length = 0.0; i < _controlPoints.size( ); i++ ) {
    getControlPoint(i, p, p+3);
    p2[ 0 ] = (float) p[marAxis::INDX_X];
    p2[ 1 ] = (float) p[marAxis::INDX_Y];
    p2[ 2 ] = (float) p[marAxis::INDX_Z];
    spX->AddPoint( (float) i, p2[ 0 ] );
    spY->AddPoint( (float) i, p2[ 1 ] );
    spZ->AddPoint( (float) i, p2[ 2 ] );
    if( i > 0 )
      length += (float) (sqrt( (p2[0]-p1[0])*(p2[0]-p1[0]) +
                               (p2[1]-p1[1])*(p2[1]-p1[1]) +
                               (p2[2]-p1[2])*(p2[2]-p1[2]) ) );
    p1[ 0 ] = p2[ 0 ]; p1[ 1 ] = p2[ 1 ]; p1[ 2 ] = p2[ 2 ];
  } // rof

  nop = ( unsigned int ) ( getParameters( )->getDoubleParam(
                                           marParameters::e_axis_discret_step ) * length );
  nip = _controlPoints.size( );
  _healthySlice = -1;

  while( _points.size( ) > 0 ) {
    delete _points[ _points.size( ) - 1 ];
    _points.pop_back( );
  } // fwhile

  while( _contours.size( ) > 0 ) {
    delete _contours[ _contours.size( ) - 1 ];
    _contours.pop_back( );
  } // fwhile

  for( i = 0; i < nop; i++ ) {
    t = (float) (( ( double ) nip - 1.0 ) / ( ( double ) nop - 1.0 ) * i);
    double* np = new double[ 3 ];
    np[ 0 ] =  0;
    np[ 1 ] =  0;
    np[ 2 ] =  0;
    np[ 0 ] = (double) (spX->Evaluate( t ));
    np[ 1 ] = (double) (spY->Evaluate( t ));
    np[ 2 ] = (double) (spZ->Evaluate( t ));
    _points.push_back( np );
    _contours.push_back( new marContour( i, getParameters( ) ) );

  } // rof

  spX->Delete( );
  spY->Delete( );
  spZ->Delete( );

}



// ----------------------------------------------------------------------------
void marAxis::sliceVolumeAxis( kVolume* vol, bool forceCnt )
{
/*
  int nCnts = ( int ) getNumberOfSplinePoints( );
  int sizeIma = getParameters( )->getSizeIma( );
  double dimIma = getParameters( )->getDimIma( );
  double voxSize = getParameters( )->getVoxelSize( );

  if( forceCnt )
    for( unsigned int i = 0; i < _contours.size( ); i++ )
      delete _contours[ i ];
    _contours.clear( );

    for( unsigned int i = 0; i < _slices.size( ); i++ ) {
      delete _slices[ i ];
      _3Dslices[ i ]->Delete( );
      _quantificationImages[ i ]->Delete( );
    }

  _slices.clear( );
  _3Dslices.clear( );
  _quantificationImages.clear( );


  // "Cutter" initialization
  vtkStructuredPoints* stPoints = vtkStructuredPoints::New( );
  // Perform "cut"
  double *p, *n;
  for( int k = 0; k < nCnts; k++ ) {
    //s = ( double )k / ( double )( nCnts - 1 );
    p = _points[k];
    n = getNormal( k );
    vtkPlaneSource* pSource = vtkPlaneSource::New( );
    pSource->SetOrigin( p[ 0 ], p[ 1 ], p[ 2 ] );
    pSource->SetPoint1( p[ 0 ] + dimIma - 1.0, p[ 1 ], p[ 2 ] );
    pSource->SetPoint2( p[ 0 ], p[ 1 ], p[ 2 ] + dimIma - 1.0 );
    pSource->SetResolution( sizeIma - 1, sizeIma - 1 );
    pSource->Update( );
    pSource->SetCenter( p[ 0 ], p[ 1 ], p[ 2 ] );
    pSource->SetNormal( n[ 0 ], n[ 1 ], n[ 2 ] );
    pSource->Update( );

    _3Dslices.push_back( vtkProbeFilter::New( ) );
    _3Dslices[ k ]->SetInput( ( vtkDataSet* )pSource->GetOutput( ) );
    _3Dslices[ k ]->SetSource( vol->castVtk( ) );
    _3Dslices[ k ]->Update( );
    pSource->Delete( );

    stPoints->GetPointData( )->
                SetScalars( _3Dslices[ k ]->GetOutput( )->
                GetPointData( )->GetScalars( ) );
    stPoints->SetDimensions( sizeIma, sizeIma, 1 );
    stPoints->SetScalarType( vol->castVtk( )->GetDataObjectType( ) );
    stPoints->Update( );
    if( forceCnt )
      {
                marContour* tmp_marContour = new marContour( k, getParameters( ) ) ;
                tmp_marContour->calculateVariables();
                _contours.push_back( tmp_marContour );
        }
    _slices.push_back( new kVolume( kVolume::UCHAR, 1, 1, 1 ) );

    *( _slices[ k ] ) = ( vtkImageData* )stPoints;  //copy
    _quantificationImages.push_back( (_slices[ k ])->castVtk( ) );
  } // rof

  stPoints->Delete( );
*/
}



// ----------------------------------------------------------------------------
int marAxis::getNumberOfContours( )
{
    return( _contours.size( ) );
}

// ----------------------------------------------------------------------------
int marAxis::getNumberOfSplinePoints( )
{
    return( _points.size( ) );
}

/*
// EED Esto toca limpiarlo o reimplementarlo en varias partes o algo
// ----------------------------------------------------------------------------
vtkPolyData* marAxis::setContour( int i, int x, int y,  std::vector< double* >* points, marContour* c )
{
        if( c != NULL ) 
        {
                c->setS( c->getS( ) );
                _contours[ i ]->copyFrom( *c );
        } 
        else 
        {
                // PS -> //    int numberOfPoints, numberOfCells; // PS
                // PS -> //    int id; // PS
                // PS -> //    gslobj_vector p( 3 ); // PS
                // PS -> //    double new_p[3]; // PS
                // PS -> //    float pf[ 3 ]; // PS
                vtkImageData* imagedata;
                
                // Contour parameters
                // PS ->     int typ = getParameters( )->getIntParam( marParameters::e_algorithm_type );
                double thr = getParameters( )->getDoubleParam( marParameters::e_threshold_isocontour );
                
                // PS ->        double thr = ( typ == marParameters::ISOCONTOURS )?
                // PS ->       getParameters( )->
                // PS ->       getDoubleParam( marParameters::e_threshold_isocontour ):
                // PS ->     ( typ == marParameters::SNAKE ) ?
                // PS ->       getParameters( )->
                // PS ->       getDoubleParam( marParameters::e_threshold_snake_isocontour ):
                // PS ->     30;
                
                
                int sizeIma = getParameters( )->getSizeIma( );
                double localint = getSignal( i );
                double sigma = getParameters( )->getDoubleParam( marParameters::e_sigma );
                double threshold = localint * thr / 100.0;
                int dimIma    = ( int ) ( ( double ) sizeIma /getParameters( )->getDoubleParam( marParameters::e_scale ) );
                int vmin = ( int )threshold;
                int vmax = ( int )threshold;
                
                x = ( x != -1 )? x: sizeIma / 2;
                y = ( y != -1 )? y: sizeIma / 2;
                
                // Isocontour
                imagedata = (vtkImageData*) ((getSlice( i , NULL ))->castVtk( ));
                float *range = imagedata->GetScalarRange();
                
                //vtkKitwareContourFilter* cntVTK = vtkKitwareContourFilter::New( );
                vtkContourFilter* cntVTK = vtkContourFilter::New( );
                cntVTK->SetInput( imagedata );
                cntVTK->SetNumberOfContours( 1 );
                //cntVTK->SetValue( 0, vmin );
                cntVTK->SetValue( 0, (range[1]*thr/100) );
                //cntVTK->SetValue( 1, vmax );
                cntVTK->Update( );
                
                vtkCleanPolyData* cpd = vtkCleanPolyData::New( );
                cpd->SetInput( cntVTK->GetOutput( ) );
                cpd->ConvertLinesToPointsOff( );
                cpd->Update( );
                
                vtkPolyDataConnectivityFilter* conn = vtkPolyDataConnectivityFilter::New( );
                //conn->SetMaxRecursionDepth( 3000 );
                
                // PS ->     //conn->SetInput( cntVTK->GetOutput( ) ); PS
                conn->SetInput( cpd->GetOutput( ) );
                
                conn->SetClosestPoint( x, y, 0 );
                conn->SetExtractionModeToClosestPointRegion( );
                conn->Update( );
                
                vtkCleanPolyData* cpd2 = vtkCleanPolyData::New( );
                cpd2->SetInput( conn->GetOutput( ) );
                cpd2->Update();
                
                vtkPolyData* polyDataResult = cpd2->GetOutput() ;
                
                double ta = _contours[ i ]->getS( );
                
                _contours[ i ]->reset( );
                
                // PS ->        gslobj_vector o( 3 ), n( 3 ), y( 3 ), vr( 3 );
                
                int id;
                int numberOfCells;
                int numberOfPoints,numberOfPoints2;
                vtkCell* cell;
                vtkIdList* pids;
                marVector p( 3 );
// PS ->                double new_p[3];
                float pf[ 3 ];

                numberOfCells = polyDataResult->GetNumberOfCells( );
        cell = polyDataResult->GetCell( 0 );
                pids = cell->GetPointIds( );
                numberOfPoints = pids->GetNumberOfIds( );

                numberOfPoints2=polyDataResult->GetNumberOfPoints();
                
                marVector o(3),n(3),y(3),vr(3);
                double cost, sint;
                y( 0 ) = y( 2 ) = 0.0;
                y( 1 ) = 1.0;
                getPoint( ( double* )o, ta );
                getTangent( ( double* )n, ta );
                vr = y.cross( n ).normalize( );
                sint = y.cross( n ).norm2( );
                cost = y.dot( n );
        
                for( int j = 0; j < numberOfPoints; j++ ) 
                {
                        id = pids->GetId( j );
                        polyDataResult->GetPoint( id, pf );
                        p( 0 ) = pf[ 0 ];
                        p( 1 ) = 0.0;
                        p( 2 ) = pf[ 1 ];
                        
                        // 3D transformation
                        p = ( p * cost ) + ( vr * ( vr.dot( p ) * ( 1.0 - cost ) ) ) -
                                ( p.cross( vr ) * sint ) + o;
                        
                // PS ->                double p_cross_vr[3];
                // PS ->                vtkMath::Cross( p, vr, p_cross_vr);
                // PS ->                for(int k=0; k<3; k++) 
                // PS ->                {
                // PS ->                        new_p[i] = (new_p[k] * cost) + ( vr(k) * ( vtkMath::Dot(vr, p) * ( 1.0 - cost ) ) ) -
                // PS ->                                ( p_cross_vr[k] * sint ) + o(i);
                // PS ->                }
                        
                        //TODO: check if p and new_p are the same...
                        _contours[ i ]->addContourPoint( p );
                } // rof
                
                // PS ->        int nbPoints=polyDataResult->GetNumberOfPoints();
                // PS ->        
                // PS ->        float* vtkContourPoint;
                // PS ->        vtkContourPoint= new float[3];
                // PS ->        double* marContourPoint;
                // PS ->        marContourPoint= new double[3];
                // PS -> 
                // PS ->        for (int j=0;j< nbPoints;j++)
                // PS ->        {
                // PS ->                vtkContourPoint=polyDataResult->GetPoint(j);
                // PS ->                marContourPoint[0]=(double)vtkContourPoint[0];
                // PS ->                marContourPoint[1]=(double)vtkContourPoint[1];
                // PS ->                marContourPoint[2]=(double)vtkContourPoint[2];
                // PS ->                _contours[ i ]->addContourPoint( marContourPoint );
                // PS ->        }
                marContour * tutu=_contours[ i ];
                _contours[ i ]->getArea();
                _contours[ i ]->getPerimeter();
                return( polyDataResult );
                //return( cpd->GetOutput( ) );
                //return( cntVTK->GetOutput( ) );
                
                 
//EED1             vtkStripper* isoStrips = vtkStripper::New( );
//EED1             isoStrips->SetInput( cpd2->GetOutput( ) );
//EED1             isoStrips->Update( );
//EED1             
//EED1                   vtkPolyData* data = isoStrips->GetOutput( );
//EED1                   data->Update( );
//EED1                   vtkCell* cell;
//EED1                   vtkIdList* pids;
//EED1                   
//EED1                     numberOfCells = data->GetNumberOfCells( );
//EED1                     cell = data->GetCell( 0 );
//EED1                     pids = cell->GetPointIds( );
//EED1                     numberOfPoints = pids->GetNumberOfIds( );

                
//EED2          _contours[ i ]->reset( );
//EED2          if( typ == marParameters::ISOCONTOURS ) {
//EED2          gslobj_vector o( 3 ), n( 3 ), y( 3 ), vr( 3 );
//EED2          double cost, sint;
//EED2          y( 0 ) = y( 2 ) = 0.0;
//EED2          y( 1 ) = 1.0;
//EED2          getPoint( ( double* )o, ta );
//EED2          getTangent( ( double* )n, ta );
//EED2          vr = y.cross( n ).normalize( );
//EED2          sint = y.cross( n ).norm2( );
//EED2          cost = y.dot( n );
//EED2          
//EED2            for( int j = 0; j < numberOfPoints; j++ ) {
//EED2            id = pids->GetId( j );
//EED2            data->GetPoint( id, pf );
//EED2            p( 0 ) = pf[ 0 ];
//EED2            p( 1 ) = 0.0;
//EED2            p( 2 ) = pf[ 1 ];
//EED2            
//EED2                  // 3D transformation
//EED2                  p = ( p * cost ) + ( vr * ( vr.dot( p ) * ( 1.0 - cost ) ) ) -
//EED2          ( p.cross( vr ) * sint ) + o;
//EED2                  
//EED2                    double p_cross_vr[3];
//EED2                    vtkMath::Cross( p, vr, p_cross_vr);
//EED2                    for(int ii=0; ii<3; ii++) {
//EED2                    new_p[i] = (new_p[ii] * cost) + ( vr(ii) * ( vtkMath::Dot(vr, p) * ( 1.0 - cost ) ) ) -
//EED2                    ( p_cross_vr[ii] * sint ) + o(i);
//EED2                    }
//EED2                    
//EED2                          //TODO: check if p and new_p are the same...
//EED2                          _contours[ i ]->addContourPoint( p );
//EED2                          } // rof
//EED2                          } else if( typ == marParameters::SNAKE || typ == marParameters::DERICHE ) {
//EED2                          // ONLY Isocontours at the moment.
//EED2                          } // fi
//EED2                          
//EED2                            cntVTK->Delete( );
//EED2                            cpd->Delete( );
//EED2                            conn->Delete( );
//EED2                            cpd2->Delete( );
//EED2          isoStrips->Delete( );

  } // fi
}

*/




// ----------------------------------------------------------------------------
void marAxis::createContour(int i, kVolume* vol)
{

        int j,id;
        int numberOfPoints,numberOfCells;
        vtkCell* cell;
        vtkIdList* pids;
        double p[ 3 ];
        vtkPolyData* vtkPolyData_2DContour;
        double x,y;

        int sizeIma                             =       getParameters( )->getSizeIma( );
        double dimIma                   =       getParameters( )->getDimIma( );

// EED 24 Dec 2007
  dimIma=sizeIma/3;

        _contours[i]                    =       new marContour( i, getParameters( ) ) ;
        vtkPolyData_2DContour   =       get2Dcontour(i,vol);

        numberOfCells                   =       vtkPolyData_2DContour->GetNumberOfCells( );
    cell                                        =       vtkPolyData_2DContour->GetCell( 0 );
        pids                                    =       cell->GetPointIds( );
        numberOfPoints                  =       pids->GetNumberOfIds( );

        for( j = 0; j < numberOfPoints; j++ ) 
        {
                id = pids->GetId( j );
                vtkPolyData_2DContour->GetPoint( id, p);
                x=p[0]-64.0;
                y=p[1]-64.0;
                x=x * dimIma / ( float ) sizeIma;
                y=y * dimIma / ( float ) sizeIma;
                _contours[ i ]->addContourPoint( x , y );
        }

        _contours[i]->do_spline();
        _contours[i]->calculateVariables();

}

// ----------------------------------------------------------------------------
void marAxis::EraseContour(int i){
        if (_3Dcontour[i]!=NULL){
                _3Dcontour[i]->Delete();
                _3Dcontour[i]=NULL;
        }

        if (_2DDiameterMax[i]!=NULL){
                _2DDiameterMax[i]->Delete();
                _2DDiameterMax[i]=NULL;
        }

        if (_2DDiameterMin[i]!=NULL){
                _2DDiameterMin[i]->Delete();
                _2DDiameterMin[i]=NULL;
        }

        if (_2Dcontours[i]!=NULL){
                _2Dcontours[i]->Delete();
                _2Dcontours[i]=NULL;
        }

        if (_contours[i]!=NULL){
                delete _contours[i];
                _contours[i]=NULL;
        }
}

// ----------------------------------------------------------------------------
void marAxis::replaceContour2D(int i,int size,double *vx,double *vy){

        EraseContour(i);


        vtkPoints *_pts = vtkPoints::New();
        _pts->SetNumberOfPoints(size);
        int j;

        for (j=0 ; j<size ; j++){
                _pts->SetPoint(j,       vx[j]   , vy[j] , 0 );
        }
//      _pts->SetPoint(0,       vx[0]   , vy[0] , 0 );

        vtkCellArray *lines = vtkCellArray::New();
        lines->InsertNextCell( size );
        for ( j=0 ; j<size+1 ; j++ ){
                lines->InsertCellPoint(j % size );
        }

        vtkPolyData *_pd = vtkPolyData::New();
        _pd->SetPoints( _pts );
        _pd->SetLines( lines );
        lines->Delete();  //do not delete lines ??
        _pts->Delete();  

        _2Dcontours[i]=_pd;
        createContour(i,NULL);
}
// ----------------------------------------------------------------------------
void marAxis::createSlice(int i , kVolume* vol){
        int k=i;
        int sizeIma = getParameters( )->getSizeIma( );

        vtkStructuredPoints* stPoints = vtkStructuredPoints::New( );
        double *p, *n;
    p = _points[k];
    n = getNormal( k );

    stPoints->GetPointData( )->SetScalars(  get3DSlice(k,vol)->GetOutput()->GetPointData()->GetScalars()  );
    stPoints->SetDimensions( sizeIma, sizeIma, 1 );


    stPoints->SetScalarType( vol->castVtk( )->GetDataObjectType( ) );
        stPoints->SetScalarTypeToShort ();
    stPoints->Update();
        
    vtkImageChangeInformation* change = vtkImageChangeInformation::New();
    change->SetInput( stPoints );  
    change->Update();    //important


    if (_slices[k]!=NULL) { delete _slices[k]; }
    _slices[k] = new kVolume( change->GetOutput() ) ;

        stPoints->Delete( );
        change->Delete();
}
// ----------------------------------------------------------------------------
void marAxis::create3DSlice(int i , kVolume* vol )
{

  int k=i;
  int sizeIma   = getParameters( )->getSizeIma( );
  double dimIma = getParameters( )->getDimIma( );

  dimIma=sizeIma/3;


//  double voxSize = getParameters( )->getVoxelSize( );

//  if( forceCnt )
//    for( unsigned int i = 0; i < _contours.size( ); i++ )
//      delete _contours[ i ];
//    _contours.clear( );

//    for( unsigned int i = 0; i < _slices.size( ); i++ ) {
//      delete _slices[ i ];
//      _3Dslices[ i ]->Delete( );
//      _quantificationImages[ i ]->Delete( );
//    }

//  _slices.clear( );
// _3Dslices.clear( );
//  _quantificationImages.clear( );

  // "Cutter" initialization
//  vtkStructuredPoints* stPoints = vtkStructuredPoints::New( );

  // Perform "cut"
  double *p, *n;
//  for( int k = 0; k < nCnts; k++ ) {
    //s = ( double )k / ( double )( nCnts - 1 );
    p = _points[k];
    n = getNormal( k );

    vtkPlaneSource* pSource = vtkPlaneSource::New( );
    pSource->SetOrigin( p[ 0 ], p[ 1 ], p[ 2 ] );
    pSource->SetPoint1( p[ 0 ] + dimIma - 1.0, p[ 1 ], p[ 2 ] );
    pSource->SetPoint2( p[ 0 ], p[ 1 ], p[ 2 ] + dimIma - 1.0 );
    pSource->SetResolution( sizeIma - 1, sizeIma - 1 );
    pSource->Update( );
    pSource->SetCenter( p[ 0 ], p[ 1 ], p[ 2 ] );
    pSource->SetNormal( n[ 0 ], n[ 1 ], n[ 2 ] );
    pSource->Update( );

        if (_3Dslices[ k ]!=NULL){ _3Dslices[ k ]->Delete(); }
    _3Dslices[ k ] = vtkProbeFilter::New( ) ;
    _3Dslices[ k ]->SetInput( ( vtkDataSet* )pSource->GetOutput( ) );
    _3Dslices[ k ]->SetSource( vol->castVtk( ) );
    _3Dslices[ k ]->Update( );
    pSource->Delete( );

//    stPoints->GetPointData( )->
//              SetScalars( _3Dslices[ k ]->GetOutput( )->
//              GetPointData( )->GetScalars( ) );
//    stPoints->SetDimensions( sizeIma, sizeIma, 1 );
//    stPoints->SetScalarType( vol->castVtk( )->GetDataObjectType( ) );
//    stPoints->Update( );
//
//    if( forceCnt )
//      {
//              marContour* tmp_marContour = new marContour( k, getParameters( ) ) ;
//              tmp_marContour->calculateVariables();
//              _contours.push_back( tmp_marContour );
//      }
//    _slices.push_back( new kVolume( kVolume::UCHAR, 1, 1, 1 ) );
//    *( _slices[ k ] ) = ( vtkImageData* )stPoints;  //copy
//    _quantificationImages.push_back( (_slices[ k ])->castVtk( ) );
//  } // rof

//  stPoints->Delete( );

//      _3Dslices[i]=NULL;
}




// ----------------------------------------------------------------------------
void marAxis::create3Dcontour(int i, kVolume* vol)
{       
        vtkPoints *_vtkPoints;
        double *o;
        double *n;
        vtkMatrix4x4* mat;
        vtkTransform* trans;
        double nn[3];
        double c[3],p1[3],p2[3];
        double nA,nB,nC;

        _vtkPoints                                              = vtkPoints::New();
        o                                                               = _points[i];
        n                                                               = getNormal( i );
        mat                                                             = vtkMatrix4x4::New();
        trans                                                   = vtkTransform::New();


        nn[0]=n[0];  nn[1]=n[1];  nn[2]=n[2];
        nC=sqrt( nn[0]*nn[0] + nn[1]*nn[1] + nn[2]*nn[2] );
        nn[0]=nn[0]/nC;  nn[1]=nn[1]/nC;  nn[2]=nn[2]/nC;

    vtkPlaneSource* pSource = vtkPlaneSource::New( );
    pSource->SetOrigin( 0, 0    , 0                             );
    pSource->SetPoint1( 1, 0    , 0                             );
    pSource->SetPoint2( 0, 0    , 1.0   );
    pSource->SetCenter( 0,0,0 );
    pSource->SetNormal( nn[ 0 ], nn[ 1 ], nn[ 2 ] );
    pSource->Update( );
//    pSource->Update( );
        pSource->GetOrigin(c);
        pSource->GetPoint1(p1);
        pSource->GetPoint2(p2);
        pSource->Delete();
        p1[0]=p1[0]-c[0];  p1[1]=p1[1]-c[1];  p1[2]=p1[2]-c[2];
        p2[0]=p2[0]-c[0];  p2[1]=p2[1]-c[1];  p2[2]=p2[2]-c[2];
        nA=sqrt( p1[0]*p1[0] + p1[1]*p1[1] + p1[2]*p1[2] );
        nB=sqrt( p2[0]*p2[0] + p2[1]*p2[1] + p2[2]*p2[2] );
        p1[0]=p1[0]/nA;  p1[1]=p1[1]/nA;  p1[2]=p1[2]/nA;
        p2[0]=p2[0]/nB;  p2[1]=p2[1]/nB;  p2[2]=p2[2]/nB;

    mat->SetElement (0,0, nn[0]);
    mat->SetElement (1,0, nn[1]);
    mat->SetElement (2,0, nn[2]);
    mat->SetElement (3,0, 0);
    mat->SetElement (0,1, p2[0]);
    mat->SetElement (1,1, p2[1]);
    mat->SetElement (2,1, p2[2]);
    mat->SetElement (3,1, 0);
        mat->SetElement (0,2, p1[0]);
    mat->SetElement (1,2, p1[1]);
    mat->SetElement (2,2, p1[2]);
    mat->SetElement (3,2, 0);
        mat->SetElement (0,3, 0);
    mat->SetElement (1,3, 0);
    mat->SetElement (2,3, 0);
    mat->SetElement (3,3, 1);

//      double deter=mat->Determinant(mat);
        trans->Identity();
        trans->SetMatrix(mat);
        float ang;
        ang=-90;
        trans->RotateWXYZ  ( ang,0,1,0); 

// EED Borrame
//      double scale = 1; 
//      trans->Scale  ( scale , scale, scale); 

        trans->Update();

        vtkMatrix4x4 *m=vtkMatrix4x4::New();
        trans->GetMatrix  (  m  );
    

        int j,numberOfPoints;
        marContour* contour2D   = getContour(i,vol);
        numberOfPoints                  = contour2D->getNumberOfPoints( );
        double fpA[3];
        double pA[3],ppp[3];


        for( j = 0; j <= numberOfPoints; j++ ) {
                contour2D->getPoint( fpA,j % numberOfPoints);

                pA[0] = fpA[0] + 0.2;  // correction EED 
                pA[1] = fpA[1] + 0.2;  // correction EED

                pA[2] = 0;

// Why this does not works...(release version)????
//              trans->TransformPoint(pA,ppp); 
// or
//              trans->TransformVector(pA,ppp);
// so..
                ppp[0] = m->GetElement(0,0)*pA[0] + m->GetElement(0,1)*pA[1] + m->GetElement(0,2)*pA[2] + m->GetElement(0,3);
                ppp[1] = m->GetElement(1,0)*pA[0] + m->GetElement(1,1)*pA[1] + m->GetElement(1,2)*pA[2] + m->GetElement(1,3);
                ppp[2] = m->GetElement(2,0)*pA[0] + m->GetElement(2,1)*pA[1] + m->GetElement(2,2)*pA[2] + m->GetElement(2,3);

//              ppp[0]=ppp[0]+o[0]-c[0]; ppp[1]=ppp[1]+o[1]-c[1]; ppp[2]=ppp[2]+o[2]-c[2];  
                ppp[0]=ppp[0]+o[0]; ppp[1]=ppp[1]+o[1]; ppp[2]=ppp[2]+o[2];  
//              ppp[0]=ppp[0]+o[0]-c[0]*0.781; ppp[1]=ppp[1]+o[1]-c[1]*0.781;; ppp[2]=ppp[2]+o[2]-c[2]*0.781;;  
                _vtkPoints->InsertNextPoint( (float)ppp[0],(float)ppp[1],(float)ppp[2] );
        }

        m->Delete();
        mat->Delete();
        trans->Delete();
        if (_3Dcontour[i]!=NULL){ _3Dcontour[i]->Delete(); }
        _3Dcontour[i]=_vtkPoints;       

}
// ----------------------------------------------------------------------------
void marAxis::createSliceImage(int i  , kVolume* vol )
{
  _quantificationImages[i] = getSlice(i,vol)->castVtk( ) ;

}
// ----------------------------------------------------------------------------
void marAxis::create2Dcontour(int i , kVolume* vol )
{
        vtkImageData* imagedata;
                
        double  thr                     = getParameters( )->getDoubleParam( marParameters::e_threshold_isocontour );
        int             sizeIma         = getParameters( )->getSizeIma( );
        double  localint        = getSignal( i );
        double  sigma           = getParameters( )->getDoubleParam( marParameters::e_sigma );
        double  threshold       = localint * thr / 100.0;

        int             dimIma          = ( int ) ( ( double ) sizeIma /getParameters( )->getDoubleParam( marParameters::e_scale ) );
        int             vmin            = ( int )threshold;
        int             vmax            = ( int )threshold;

        int             x                       = -1;
        int             y                       = -1;

        x = ( x != -1 )? x: sizeIma / 2;
        y = ( y != -1 )? y: sizeIma / 2;
        
        // Isocontour
        imagedata = (vtkImageData*) ((getSlice( i , vol ))->castVtk( ));
        double *range = imagedata->GetScalarRange();    

        //vtkKitwareContourFilter* cntVTK = vtkKitwareContourFilter::New( );
        vtkContourFilter* cntVTK = vtkContourFilter::New( );
        cntVTK->SetInput( imagedata );
        cntVTK->SetNumberOfContours( 1 );
        //cntVTK->SetValue( 0, vmin );
        cntVTK->SetValue( 0, (range[1]*thr/100) );
        //cntVTK->SetValue( 1, vmax );
        cntVTK->Update( );
                
        vtkCleanPolyData* cpd = vtkCleanPolyData::New( );
        cpd->SetInput( cntVTK->GetOutput( ) );
        cpd->ConvertLinesToPointsOff( );
        cpd->Update( );
                
        vtkPolyDataConnectivityFilter* conn = vtkPolyDataConnectivityFilter::New( );
        //conn->SetMaxRecursionDepth( 3000 );
                
        // PS ->     //conn->SetInput( cntVTK->GetOutput( ) ); PS
        conn->SetInput( cpd->GetOutput( ) );
                
        conn->SetClosestPoint( x, y, 0 );
        conn->SetExtractionModeToClosestPointRegion( );
        conn->Update( );
                
        vtkCleanPolyData* cpd2 = vtkCleanPolyData::New( );
        cpd2->SetInput( conn->GetOutput( ) );
        cpd2->Update();

        vtkStripper* vtkstripper = vtkStripper::New( );
        vtkstripper->SetInput( cpd2->GetOutput() );
        vtkstripper->Update();

        vtkPolyData* polyDataResult =  vtkstripper->GetOutput();
    polyDataResult->Update( );
        
        cntVTK          -> Delete();
        cpd2            -> Delete();
        cpd                     -> Delete();
        conn            -> Delete();

/*
// EED 27 Oct 2007
        double *p;
        int ii,size=polyDataResult->GetNumberOfPoints();
printf("EED marAxis::create2Dcontour size %d\n",size);
        for (ii=0;ii<size;ii++)
        {
                p=polyDataResult->GetPoint(ii);
printf("EED marAxis::create2Dcontour %d >> %f %f %f\n",ii,p[0], p[1], p[2]);
                p[2]=-100;      
        }

        for (ii=0;ii<size;ii++)
        {
                p=polyDataResult->GetPoint(ii);
printf("EED marAxis::create2Dcontour xxxxx %d >> %f %f %f\n",ii,p[0], p[1], p[2]);
        }
*/

        if ( _2Dcontours[i]!=NULL ) { _2Dcontours[i]->Delete(); }
        _2Dcontours[i]=polyDataResult;
}
// ----------------------------------------------------------------------------
void marAxis::create2DDiameterMin(      int i , kVolume* vol ){
        double p1[2],p2[2];
        marContour *marcontour=getContour(i,vol);
        marcontour->getMinimumLine(p1,p2);
        vtkPoints *_vtkPoints = vtkPoints::New();
        int sizeIma = getParameters( )->getSizeIma( );

        double dimIma = getParameters( )->getDimIma( );

// EED 24 Dec 2007
  dimIma=sizeIma/3;

        double factor=( float ) sizeIma / dimIma;
        p1[0]=p1[0]*factor+sizeIma/2;
        p1[1]=p1[1]*factor+sizeIma/2;
        p2[0]=p2[0]*factor+sizeIma/2;
        p2[1]=p2[1]*factor+sizeIma/2;
        _vtkPoints->InsertNextPoint( p1[0] , p1[1] , 1 );
        _vtkPoints->InsertNextPoint( p2[0] , p2[1] , 1 );
        if ( _2DDiameterMin[i]!=NULL ) { _2DDiameterMin[i]->Delete(); }
        _2DDiameterMin[i] = _vtkPoints;
}
// ----------------------------------------------------------------------------
void marAxis::create2DDiameterMax(      int i , kVolume* vol ){
        double p1[2],p2[2];
        marContour *marcontour=getContour(i,vol);
        marcontour->getMaximumLine(p1,p2);
        vtkPoints *_vtkPoints = vtkPoints::New();
        int sizeIma = getParameters( )->getSizeIma( );

        double dimIma = getParameters( )->getDimIma( );
// EED 24 Dec 2007
dimIma=sizeIma/3;

        double factor=( float ) sizeIma / dimIma;
        p1[0]=p1[0]*factor+sizeIma/2;
        p1[1]=p1[1]*factor+sizeIma/2;
        p2[0]=p2[0]*factor+sizeIma/2;
        p2[1]=p2[1]*factor+sizeIma/2;
        _vtkPoints->InsertNextPoint( p1[0] , p1[1] , 1 );
        _vtkPoints->InsertNextPoint( p2[0] , p2[1] , 1 );
        if ( _2DDiameterMax[i]!=NULL ) { _2DDiameterMax[i]->Delete(); }
        _2DDiameterMax[i] = _vtkPoints;
}
// ----------------------------------------------------------------------------
void marAxis::createEmptyVectors()
{
        int nCnts = ( int ) getNumberOfSplinePoints( );
        clearAllVectors();
        int i;
//      int xxxx=getNumberOfContours();
        for (i=0;i<nCnts;i++){
                _contours.push_back( NULL );
                _slices.push_back( NULL );
                _3Dslices.push_back( NULL );
                _3Dcontour.push_back( NULL );
                _quantificationImages.push_back( NULL );
                _2Dcontours.push_back( NULL );
                _2DDiameterMin.push_back( NULL );
                _2DDiameterMax.push_back( NULL );
        }
}
// ----------------------------------------------------------------------------
void marAxis::clearAllVectors()
{
        unsigned int i;
        
        for (i=0;i<_contours.size();i++){
                if (_contours[i]!=NULL) {                                       //  DATA-MODEL-2D
                        delete _contours[i];
                        _contours[i]=NULL;
                }
        }
        _contours.clear();

        for (i=0;i<_slices.size();i++){                                 //  DATA-MODEL-Voxel XxYx1
                if (_slices[i]!=NULL) {                         
                        delete _slices[i];
                        _slices[i]=NULL;
                }
        }
        _slices.clear();


        for (i=0;i<_3Dslices.size();i++){                                       //  VISUALISATION_VTK 3D
                if (_3Dslices[i]!=NULL) {
                        _3Dslices[i]->Delete();
                        _3Dslices[i]=NULL;
                }
        }
        _3Dslices.clear();

        for (i=0;i<_3Dcontour.size();i++){                              //  VISUALISATION_VTK 3D
                if (_3Dcontour[i]!=NULL) {                              
                        _3Dcontour[i]->Delete();
                        _3Dcontour[i]=NULL;
                }
        }
        _3Dcontour.clear();

        for (i=0;i<_quantificationImages.size();i++){   //  VISUALISATION_VTK 2D
                if (_quantificationImages[i]!=NULL) {                           
//                      _quantificationImages[i]->Delete();
//                      _quantificationImages[i]=NULL;
                }
        }
        _quantificationImages.clear();

        for (i=0;i<_2Dcontours.size();i++){                             //  VISUALISATION_VTK 2D
                if (_2Dcontours[i]!=NULL) {                             
// EED ???  This object was allready erased bye the VTK pipeline
//                      _2Dcontours[i]->Delete();
                        _2Dcontours[i]=NULL;
                }
        }
        _2Dcontours.clear();

        for (i=0;i<_2DDiameterMin.size();i++){                  //  VISUALISATION_VTK 2D
                if (_2DDiameterMin[i]!=NULL) {                          
                        _2DDiameterMin[i]->Delete() ;
                        _2DDiameterMin[i]=NULL;
                }
        }
        _2DDiameterMin.clear();

        for (i=0;i<_2DDiameterMax.size();i++){                  //  VISUALISATION_VTK 2D
                if (_2DDiameterMax[i]!=NULL) {                          
                        _2DDiameterMax[i]->Delete();
                        _2DDiameterMax[i]=NULL;
                }
        }
        _2DDiameterMax.clear();

}
// ----------------------------------------------------------------------------
void marAxis::eraseContourVectorsContent()
{
        unsigned int i;
        for (i=0;i<_contours.size();i++){
                if (_contours[i]!=NULL) {                                       //  DATA-MODEL-2D
                        delete _contours[i];
                        _contours[i]=NULL;
                }
        }
        for (i=0;i<_2Dcontours.size();i++){                             //  VISUALISATION_VTK 2D
                if (_2Dcontours[i]!=NULL) {                             
                        _2Dcontours[i]->Delete();
                        _2Dcontours[i]=NULL;
                }
        }
        for (i=0;i<_2DDiameterMin.size();i++){                          //  VISUALISATION_VTK 2D
                if (_2DDiameterMin[i]!=NULL) {                          
                        _2DDiameterMin[i]->Delete();
                        _2DDiameterMin[i]=NULL;
                }
        }
        for (i=0;i<_2DDiameterMax.size();i++){                          //  VISUALISATION_VTK 2D
                if (_2DDiameterMax[i]!=NULL) {                          
                        _2DDiameterMax[i]->Delete();
                        _2DDiameterMax[i]=NULL;
                }
        }
        for (i=0;i<_3Dcontour.size();i++){                              //  VISUALISATION_VTK 3D
                if (_3Dcontour[i]!=NULL) {                              
                        _3Dcontour[i]->Delete();
                        _3Dcontour[i]=NULL;
                }
        }
}

// ----------------------------------------------------------------------------
marContour* marAxis::getContour(int i, kVolume* vol)                            // DATA-MODEL 2D
{                                                               
        if (_contours[i]==NULL){ createContour(i, vol); }
        return _contours[i];
} 
// ----------------------------------------------------------------------------
kVolume* marAxis::getSlice(int i, kVolume* vol)                                         // DATA-MODEL-Voxel XxYx1
{
        if (_slices[i]==NULL){ createSlice( i , vol ); }
        return _slices[i];
}
// ----------------------------------------------------------------------------

bool marAxis::if3DcontourExist(int i)
{
        bool result=true;
        if (_3Dcontour[i]==NULL)
        {
                result=false;
        } 
        return result;
}

// ----------------------------------------------------------------------------
void marAxis::Save3Dcontour(FILE *ff,int id)
{
        int i,size;
        double point[3];
        if (_3Dcontour[id]!=NULL)
        {
//Ramiro                fprintf(ff,"contour_id: %d \n", id);
                size = _3Dcontour[id]->GetNumberOfPoints();
//Ramiro                fprintf(ff,"numberOfPoints: %d \n", size);
            for ( i=0 ; i<size  ;i++ ){
                _3Dcontour[id]->GetPoint( i , point );
                        fprintf(ff,"%f %f %f %d\n", point[0], point[1],point[2],id);
            }
        }
}
// ----------------------------------------------------------------------------
void marAxis::SaveExisting3DContours(FILE *ff)
{
        if (ff!=NULL)
        {
                // counting existing 3Dcontours
                int acum=0;
                int i,size=_3Dcontour.size();
                for (i=0; i<size; i++)
                {
                        if (_3Dcontour[i]!=NULL)
                        {
                                acum++;
                        }
                }
//Ramiro                fprintf(ff, "NumberOf3DContours: %d \n", acum);
                // saving existing 3Dcontours
                for (i=0; i<size; i++)
                {
                        Save3Dcontour(ff,i);
                }               
        }
}

// ----------------------------------------------------------------------------
vtkProbeFilter* marAxis::get3DSlice(int i, kVolume* vol){                       // VISUALISATION-VTK 3D
        if (_3Dslices[i]==NULL){ create3DSlice(i,vol); }
        return _3Dslices[i];
}
// ----------------------------------------------------------------------------
vtkPoints* marAxis::get3Dcontour(int i, kVolume* vol){                          // VISUALISATION-VTK 3D
        if (_3Dcontour[i]==NULL){ create3Dcontour(i, vol); }
        return _3Dcontour[i];
}
// ----------------------------------------------------------------------------
vtkImageData* marAxis::getSliceImage(int i, kVolume* vol){                      // VISUALISATION-VTK 2D
        if (_quantificationImages[i]==NULL){ createSliceImage(i,vol); }
        return _quantificationImages[i];
}
// ----------------------------------------------------------------------------
vtkPolyData* marAxis::get2Dcontour(int i , kVolume* vol){                       // VISUALISATION-VTK 2D
        if (_2Dcontours[i]==NULL){ create2Dcontour(i,vol); }
        return _2Dcontours[i];
}
// ----------------------------------------------------------------------------
vtkPoints* marAxis::get2DDiameterMin(int i , kVolume* vol){             // VISUALISATION-VTK 2D
        if (_2DDiameterMin[i]==NULL){ create2DDiameterMin(i,vol); }
        return _2DDiameterMin[i];
}
// ----------------------------------------------------------------------------
vtkPoints* marAxis::get2DDiameterMax(int i , kVolume* vol){             // VISUALISATION-VTK 2D
        if (_2DDiameterMax[i]==NULL){ create2DDiameterMax(i,vol); }
        return _2DDiameterMax[i];
}
// ----------------------------------------------------------------------------
double marAxis::getAxisLenght(int pIni,int pEnd){

    if (pIni>pEnd){
                int tmp=pIni;
                pIni=pEnd;
                pEnd=tmp;
        }

        marVector* pO = new marVector(2);
        marVector* pF = new marVector(2);
        double L;
    unsigned int j;
    for( j = pIni, L = 0.0; j < pEnd-1; j++ ) {
                (*pO)=_points[j];
                (*pF)=_points[j+1];
                (*pF)=(*pF)-(*pO);
                L+=pF->norm2();
    } // rof
        delete pO;
        delete pF;
        return L;
}
// ----------------------------------------------------------------------------
void marAxis::calculateTotalAxisLenght(){
        _totalAxisLenght = getAxisLenght( 0 , _points.size( ) );
}
// ----------------------------------------------------------------------------
double marAxis::getTotalLength(){
        if (_totalAxisLenght==-1){
                calculateTotalAxisLenght();
        }
        return _totalAxisLenght;
}

// ----------------------------------------------------------------------------
void marAxis::calculateSubAxisLength(){
        _subAxisLenght = getAxisLenght( _startQuant , _finishQuant );
}
// ----------------------------------------------------------------------------
double marAxis::getSubAxisLength(){
        if (_subAxisLenght==-1){
                calculateSubAxisLength();
        }
        return _subAxisLenght;
}
// ----------------------------------------------------------------------------
double marAxis::getAverageArea(int pIni, int pEnd, kVolume* vol){
        marContour      *marcontourHealthy;
        int                     ihealthySlice,itemp;
        double          acumArea        =       0;
        for (ihealthySlice = pIni ; ihealthySlice <= pEnd ; ihealthySlice++ ){
                itemp=ihealthySlice;
                if (itemp<0)                            { itemp=0;                                      }
                if (itemp>=_points.size( ))     { itemp=_points.size( )-1;      }
                marcontourHealthy = getContour( itemp , vol );
                acumArea = acumArea + marcontourHealthy->getArea();
        }
        return acumArea / (pEnd - pIni + 1);
}
// ----------------------------------------------------------------------------
void marAxis::calculateReferenceArea(kVolume* vol){
        _referenceArea = getAverageArea(_healthySliceStart,_healthySliceEnd,  vol);
}
// ----------------------------------------------------------------------------
double marAxis::getReferenceArea(kVolume* vol){
        if (_referenceArea==-1){
                calculateReferenceArea(vol);
        }
        return _referenceArea;
}
// ----------------------------------------------------------------------------
void marAxis::calculateReferenceAverDiam(kVolume* vol){
        marContour *marcontourHealthy;
        int ihealthySlice;
        double acumDiam=0;
        for (ihealthySlice=_healthySliceStart; ihealthySlice<=_healthySliceEnd; ihealthySlice++){
                marcontourHealthy = getContour( ihealthySlice , vol );
                acumDiam = acumDiam + marcontourHealthy->getAverageDiameter();
        }
        _referenceAverDiam = acumDiam / (_healthySliceEnd-_healthySliceStart+1);
}
// ----------------------------------------------------------------------------
double marAxis::getReferenceAverDiam(kVolume* vol){
        if (_referenceAverDiam==-1){
                calculateReferenceAverDiam(vol);
        }
        return _referenceAverDiam;
}

// ----------------------------------------------------------------------------
void marAxis::reset( )
{
//EED Borrame
//    unsigned int i;
    kCurve::reset( );

    _description                = "";
    _healthySlice               = -1;
    _startQuant                 = -1;
    _finishQuant                = -1;
    _actualQuant                = -1;

        _subAxisLenght          = -1;
        _referenceArea          = -1;
        _referenceAverDiam      = -1;

        Delete( );

// EED Borrame  
/*
    for( i = 0; i < _contours.size( ); i++ )
      delete _contours[ i ];
    _contours.clear( );
    for( i = 0; i < _slices.size( ); i++ )
        delete _slices[ i ];
    _slices.clear( );
*/
}

// ----------------------------------------------------------------------------
void marAxis::copyFrom( const marObject& from )
{ // TODO
}

// ----------------------------------------------------------------------------
bool marAxis::save( std::ofstream& os )
{
    double p[ INDX_count ];
    unsigned int i;

    i = _description.length( );
    os.write( ( const char* )&i, sizeof( int ) );
    os.write( ( char* )_description.c_str( ), i * sizeof( char ) );

    i = getNumberOfControlPoints( );
    os.write( ( const char* )&i, sizeof( int ) );
    for( i = 0; i < getNumberOfControlPoints( ); i++ ) {

        memcpy( p, _controlPoints[ i ], 3 * sizeof( double ) );
        memcpy( p + 3, _controlState[ i ],
                ( INDX_count - 3 ) * sizeof( double ) );
        os.write( ( const char* )p, INDX_count * sizeof( double ) );

    } // rof
    i = getNumberOfContours( );
    os.write( ( const char* )&i, sizeof( int ) );
    for( i = 0; i < getNumberOfContours( ); i++ )
      _contours[ i ]->save( os );

    return( true );
}

// ----------------------------------------------------------------------------
bool marAxis::load( std::ifstream& is )
{
    double p[ INDX_count ];
    int i, n;

    reset( );

    is.read( ( char* )&n, sizeof( int ) );
    _description.resize( n );
    is.read( ( char* )_description.c_str( ), n * sizeof( char ) );

    is.read( ( char* )&n, sizeof( int ) );
    for( i = 0; i < n; i++ ) {

        is.read( ( char* )p, INDX_count * sizeof( double ) );
        addControlPoint( p, p + 3 );

    } // rof
    is.read( ( char* )&n, sizeof( int ) );
    for( i = 0; i < n; i++ ) {

      _contours.push_back( new marContour( 0, getParameters( ) ) );
      _contours[ i ]->load( is );

    } // rof
    return( true );
}
// ----------------------------------------------------------------------------
vtkPolyData* marAxis::Draw( )
{
  unsigned int i, j;
  double *p;

  vtkPoints* allPoints = vtkPoints::New( );
  vtkCellArray* allTopology = vtkCellArray::New( );

  allTopology->InsertNextCell( _points.size( ) );
  for( i = 0, j=0; i < _points.size( ); i++, j++ ) {
    p = _points[i];
    allPoints->InsertNextPoint( p[ 0 ],  p[ 1 ], p[ 2 ] );
    allTopology->InsertCellPoint( j );
  } // rof

  _allData = vtkPolyData::New( );
  _allData->SetPoints( allPoints );
  _allData->SetLines( allTopology );
  allPoints->Delete();
  allTopology->Delete();

  return ( _allData );
}
// ----------------------------------------------------------------------------
vtkPolyData *marAxis::GetAxisData( )
{
        unsigned int i;
        double point[ 3 ];
        double radio;
        double p[marAxis::INDX_count];


        vtkPoints                       *allPoints              = vtkPoints::New( );
        vtkCellArray            *allTopology    = vtkCellArray::New( );
        vtkDoubleArray          *allRadios              = vtkDoubleArray::New( );
        allRadios->SetName("radio");

        for( i = 0; i < _controlPoints.size( ); i++ ) {
                getControlPoint(i, p, p+3);
                point[ 0 ]      =  p[marAxis::INDX_X];
                point[ 1 ]      =  p[marAxis::INDX_Y];
                point[ 2 ]      =  p[marAxis::INDX_Z];
                radio           =  p[marAxis::INDX_RAYON]*2.0;
                allPoints    -> InsertPoint( i, point[ 0 ] ,  point[ 1 ] , point[ 2 ]  );  //para saber exactamante el indice que se le asigno

                allRadios       -> InsertValue(i,radio);
 
    if(i>0){
       allTopology->InsertNextCell(2);
       allTopology->InsertCellPoint(i);
       allTopology->InsertCellPoint(i-1);
      }

        } // rof

        _allData = vtkPolyData::New( );
        _allData        -> SetPoints( allPoints );
        _allData        -> SetLines( allTopology );
        _allData        -> GetPointData()->SetScalars(allRadios);

        
        return ( _allData );
}
// ----------------------------------------------------------------------------
void marAxis::Delete( ){
        clearAllVectors();

        int i;
        for (i=0;i<_points.size();i++){         
                if (_points[i]!=NULL) {                         
                        delete _points[i];
                        _points[i]=NULL;
                }
        }
        _points.clear();
        _signal.clear();


        if(_allData) {
                _allData->Delete();
                _allData = NULL;
        }

//  reset();
}
// ----------------------------------------------------------------------------
void marAxis::setHealthySlice( int hsS, int hs, int hsE ){ 
          _healthySliceStart    = hsS;
          _healthySlice                 = hs;   
          _healthySliceEnd              = hsE;
          _referenceArea=-1;
          _referenceAverDiam=-1;
  };
// ----------------------------------------------------------------------------
void marAxis::setStartQuant( int sq ){
        _startQuant             =       sq;             
        _subAxisLenght  =       -1;
};
// ----------------------------------------------------------------------------
void marAxis::setFinishQuant( int fq ){ 
        _finishQuant    =       fq;             
        _subAxisLenght  =       -1;
};

// ----------------------------------------------------------------------------
void marAxis::AddPointToList(double x, double y, double z, int signal)
{
        double *p;
        p=(double *)malloc(sizeof(double)*3);
        p[0]=x;
        p[1]=y;
        p[2]=z;
        _points.push_back( p );
        _signal.push_back( signal );
}



// eof - axis.cxx