lib/maracasVisuLib/src/kernel/include/marAxis.cpp

   1 /*# ---------------------------------------------------------------------
   2 #
   3 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
   4 #                        pour la Sant�)
   5 # Authors : Eduardo Davila, Frederic Cervenansky, Claire Mouton
   6 # Previous Authors : Laurent Guigues, Jean-Pierre Roux
   7 # CreaTools website : www.creatis.insa-lyon.fr/site/fr/creatools_accueil
   8 #
   9 #  This software is governed by the CeCILL-B license under French law and
  10 #  abiding by the rules of distribution of free software. You can  use,
  11 #  modify and/ or redistribute the software under the terms of the CeCILL-B
  12 #  license as circulated by CEA, CNRS and INRIA at the following URL
  13 #  http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  14 #  or in the file LICENSE.txt.
  15 #
  16 #  As a counterpart to the access to the source code and  rights to copy,
  17 #  modify and redistribute granted by the license, users are provided only
  18 #  with a limited warranty  and the software's author,  the holder of the
  19 #  economic rights,  and the successive licensors  have only  limited
  20 #  liability.
  21 #
  22 #  The fact that you are presently reading this means that you have had
  23 #  knowledge of the CeCILL-B license and that you accept its terms.
  24 # ------------------------------------------------------------------------ */
  25
  26 /*=========================================================================
  27
  28   Program:   wxMaracas
  29   Module:    $RCSfile: marAxis.cpp,v $
  30   Language:  C++
  31   Date:      $Date: 2012/11/15 14:15:31 $
  32   Version:   $Revision: 1.2 $
  33
  34   Copyright: (c) 2002, 2003
  35   License:
  36
  37      This software is distributed WITHOUT ANY WARRANTY; without even
  38      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  39      PURPOSE.  See the above copyright notice for more information.
  40
  41 =========================================================================*/
  42 #include <vtkPlaneSource.h>
  43 #include <vtkKochanekSpline.h>
  44 #include <vtkPolyDataConnectivityFilter.h>
  45 #include <vtkContourFilter.h>
  46 #include <vtkProbeFilter.h>
  47
  48
  49 #include <vtkImageChangeInformation.h>
  50 #include <vtkImageShiftScale.h>
  51 #include <vtkDoubleArray.h>
  52
  53
  54 #include "marAxis.h"
  55 #include "marVector.h"
  56
  57 #include "ExtractionAxe.h"
  58 #include "FonctionsGrales.h"
  59
  60 #include <vtkCellArray.h>
  61 #include <vtkCell.h>
  62 #include <vtkCleanPolyData.h>
  63 #include <vtkIdList.h>
  64 #include <vtkMath.h>
  65 #include <vtkMatrix4x4.h>
  66 #include <vtkPointData.h>
  67 #include <vtkPolyData.h>
  68 #include <vtkPlane.h>
  69 #include <vtkSetGet.h>
  70 #include <vtkStripper.h>
  71 #include <vtkStructuredPoints.h>
  72 #include <vtkTransform.h>
  73
  74 // ----------------------------------------------------------------------------
  75 marAxis::marAxis( marParameters* p ) :
  76        marObject( p ),
  77            kCurve( 3, INDX_count - 3 ),
  78       _healthySlice( -1 ), _startQuant( -1 ),
  79       _finishQuant( -1 ), _actualQuant( -1 )
  80
  81 {
  82         _totalAxisLenght=-1;
  83         _subAxisLenght=-1;
  84         _referenceArea=-1;
  85         _referenceAverDiam=-1;
  86
  87
  88     _allData = NULL;
  89
  90         _healthySliceStart=-1;
  91         _healthySliceEnd=-1;
  92     reset( );
  93         calibration = false;
  94 }
  95
  96 // ----------------------------------------------------------------------------
  97 void marAxis::addAxisPoint( double* p )
  98 {
  99     double pt[ INDX_count ];
 100
 101     memcpy( pt, p, sizeof( POINTAXE ) );
 102     pt[ INDX_SIGNALVALUE ] = 0;
 103     addControlPoint( pt, pt + 3 );
 104 }
 105
 106 // ----------------------------------------------------------------------------
 107
 108 double* marAxis::getNormal( unsigned int slice ) {
 109
 110     double* ret = new double[ 3 ];
 111
 112     if ( slice == 0 ) {
 113       ret[0] = _points[ slice + 1 ][0] - _points[ slice ][0];
 114       ret[1] = _points[ slice + 1 ][1] - _points[ slice ][1];
 115       ret[2] = _points[ slice + 1 ][2] - _points[ slice ][2];
 116     }
 117     else if ( slice == _points.size( ) - 1 ) {
 118       ret[0] = _points[ slice - 1 ][0] - _points[ slice ][0];
 119       ret[1] = _points[ slice - 1 ][1] - _points[ slice ][1];
 120       ret[2] = _points[ slice - 1 ][2] - _points[ slice ][2];
 121     }
 122     else if ( slice > 0 && slice < _points.size( ) - 1 ) {
 123       ret[0] = _points[ slice - 1 ][0] - _points[ slice + 1 ][0];
 124       ret[1] = _points[ slice - 1 ][1] - _points[ slice + 1 ][1];
 125       ret[2] = _points[ slice - 1 ][2] - _points[ slice + 1 ][2];
 126     }
 127
 128     return( ret );
 129 }
 130
 131 // ----------------------------------------------------------------------------
 132 void marAxis::changeAxisResolution( )
 133 {/*
 134    int i;
 135    int nOut = ( int )ceil(
 136    getParameters( )->
 137    getDoubleParam( marParameters::e_axis_discret_step )
 138    * ( double )( GetNumberOfPoints( ) )
 139    );
 140
 141    i = GetNumberOfPoints( );
 142    double d = getParameters( )->
 143    getDoubleParam( marParameters::e_axis_discret_step );
 144
 145    SetResolution( nOut );
 146    for( i = 0; i < _contours.size( ); i++ )
 147    delete _contours[ i ];
 148    _contours.clear( );
 149    for( i = 0; i < nOut; i++ )
 150    _contours.push_back( new marContour( getParameters( ) ) );
 151  */
 152 }
 153
 154 // ----------------------------------------------------------------------------
 155 void marAxis::calculateSignal( kVolume* vol )
 156 {
 157    /**
 158     * PSIGNAL_USHORT <-> unsigned short *
 159     */
 160     /* PSIGNAL_FLOAT PutSignalValueFloat(PSIGNAL_FLOAT signal, int pos, float value)
 161     {
 162       signal[pos] = value;
 163       return signal;
 164     }*/
 165   ushort maxlocal;
 166   double *p;
 167
 168   double rayon;
 169   unsigned int i, mask_size = getParameters( )-> getIntParam( marParameters::e_mask_size );
 170
 171   while( _signal.size( ) > 0 ) _signal.pop_back( );
 172   // Intensity signal
 173
 174   PSIGNAL_USHORT signal = ( PSIGNAL_USHORT ) IdSigAlloc( _points.size( ), SIG_USHORT );
 175   for( i = 0; i < _points.size( ); i++ ) {
 176     p = _points[i];
 177     rayon = RayonLocal(
 178                       ( int ) ( p[ 0 ] ),
 179                       ( int ) ( p[ 1 ] ),
 180                       ( int ) ( p[ 2 ] ),
 181                       _points_disc, getNumberOfControlPoints( ));
 182
 183     maxlocal = vol->GetMaxIntSphere2( p, rayon );
 184     PutSignalValue( signal, i, maxlocal );
 185   } // rof
 186
 187   PSIGNAL_FLOAT mask_signal = ( PSIGNAL_FLOAT )IdSigAlloc( mask_size, SIG_FLOAT );
 188
 189   for( i = 0; i < mask_size; i++ )
 190     PutSignalValueFloat( mask_signal, i, 1 );
 191
 192   PSIGNAL_USHORT signalfilt = ( PSIGNAL_USHORT )IdSigAlloc( _points.size( ), SIG_USHORT );
 193   signalfilt = MonSigConvolve( signal, mask_signal, signalfilt, 1.0 / ( double )mask_size, 0 );
 194   for( i = 0; i < _points.size( ); i++ ){
 195     _signal.push_back( GetSignalValue( signalfilt, i ) );
 196   }
 197   IdSigFree( signal );
 198   IdSigFree( mask_signal );
 199   IdSigFree( signalfilt );
 200 }
 201
 202 // ----------------------------------------------------------------------------
 203 void marAxis::start( )
 204 {
 205     /* TODO
 206        int swap, n = GetNumberOfPoints( );
 207
 208        swap = GSL_MIN( _startQuant, _finishQuant );
 209        _finishQuant = GSL_MAX( _startQuant, _finishQuant );
 210        _startQuant = swap;
 211
 212        _startQuant = ( _startQuant >= 0 )? _startQuant: 0;
 213        _finishQuant = ( _finishQuant >= n )? _finishQuant: n - 1;
 214        _actualQuant = _startQuant;
 215     */
 216 }
 217
 218 // ----------------------------------------------------------------------------
 219 void marAxis::cut( int slice, bool up )
 220 { // TODO
 221 }
 222
 223
 224 // ----------------------------------------------------------------------------
 225 void marAxis::doSpline( )
 226 {
 227   unsigned int i;
 228   unsigned int nop, nip;
 229   float length=0, t;
 230   float p1[ 3 ], p2[ 3 ];
 231   double p[marAxis::INDX_count];
 232
 233   vtkKochanekSpline* spX = vtkKochanekSpline::New( );
 234   vtkKochanekSpline* spY = vtkKochanekSpline::New( );
 235   vtkKochanekSpline* spZ = vtkKochanekSpline::New( );
 236
 237   for( i = 0, length = 0.0; i < _controlPoints.size( ); i++ ) {
 238     getControlPoint(i, p, p+3);
 239     p2[ 0 ] = (float) p[marAxis::INDX_X];
 240     p2[ 1 ] = (float) p[marAxis::INDX_Y];
 241     p2[ 2 ] = (float) p[marAxis::INDX_Z];
 242     spX->AddPoint( (float) i, p2[ 0 ] );
 243     spY->AddPoint( (float) i, p2[ 1 ] );
 244     spZ->AddPoint( (float) i, p2[ 2 ] );
 245     if( i > 0 )
 246       length += (float) (sqrt( (p2[0]-p1[0])*(p2[0]-p1[0]) +
 247                                (p2[1]-p1[1])*(p2[1]-p1[1]) +
 248                                (p2[2]-p1[2])*(p2[2]-p1[2]) ) );
 249     p1[ 0 ] = p2[ 0 ]; p1[ 1 ] = p2[ 1 ]; p1[ 2 ] = p2[ 2 ];
 250   } // rof
 251
 252   nop = ( unsigned int ) ( getParameters( )->getDoubleParam(
 253                                            marParameters::e_axis_discret_step ) * length );
 254   nip = _controlPoints.size( );
 255   _healthySlice = -1;
 256
 257   while( _points.size( ) > 0 ) {
 258     delete _points[ _points.size( ) - 1 ];
 259     _points.pop_back( );
 260   } // fwhile
 261
 262   while( _contours.size( ) > 0 ) {
 263     delete _contours[ _contours.size( ) - 1 ];
 264     _contours.pop_back( );
 265   } // fwhile
 266
 267   for( i = 0; i < nop; i++ ) {
 268     t = (float) (( ( double ) nip - 1.0 ) / ( ( double ) nop - 1.0 ) * i);
 269     double* np = new double[ 3 ];
 270     np[ 0 ] =  0;
 271     np[ 1 ] =  0;
 272     np[ 2 ] =  0;
 273     np[ 0 ] = (double) (spX->Evaluate( t ));
 274     np[ 1 ] = (double) (spY->Evaluate( t ));
 275     np[ 2 ] = (double) (spZ->Evaluate( t ));
 276     _points.push_back( np );
 277     _contours.push_back( new marContour( i, getParameters( ) ) );
 278
 279   } // rof
 280
 281   spX->Delete( );
 282   spY->Delete( );
 283   spZ->Delete( );
 284
 285 }
 286
 287
 288
 289 // ----------------------------------------------------------------------------
 290 void marAxis::sliceVolumeAxis( kVolume* vol, bool forceCnt )
 291 {
 292 /*
 293   int nCnts = ( int ) getNumberOfSplinePoints( );
 294   int sizeIma = getParameters( )->getSizeIma( );
 295   double dimIma = getParameters( )->getDimIma( );
 296   double voxSize = getParameters( )->getVoxelSize( );
 297
 298   if( forceCnt )
 299     for( unsigned int i = 0; i < _contours.size( ); i++ )
 300       delete _contours[ i ];
 301     _contours.clear( );
 302
 303     for( unsigned int i = 0; i < _slices.size( ); i++ ) {
 304       delete _slices[ i ];
 305       _3Dslices[ i ]->Delete( );
 306       _quantificationImages[ i ]->Delete( );
 307     }
 308
 309   _slices.clear( );
 310   _3Dslices.clear( );
 311   _quantificationImages.clear( );
 312
 313
 314   // "Cutter" initialization
 315   vtkStructuredPoints* stPoints = vtkStructuredPoints::New( );
 316   // Perform "cut"
 317   double *p, *n;
 318   for( int k = 0; k < nCnts; k++ ) {
 319     //s = ( double )k / ( double )( nCnts - 1 );
 320     p = _points[k];
 321     n = getNormal( k );
 322     vtkPlaneSource* pSource = vtkPlaneSource::New( );
 323     pSource->SetOrigin( p[ 0 ], p[ 1 ], p[ 2 ] );
 324     pSource->SetPoint1( p[ 0 ] + dimIma - 1.0, p[ 1 ], p[ 2 ] );
 325     pSource->SetPoint2( p[ 0 ], p[ 1 ], p[ 2 ] + dimIma - 1.0 );
 326     pSource->SetResolution( sizeIma - 1, sizeIma - 1 );
 327     pSource->Update( );
 328     pSource->SetCenter( p[ 0 ], p[ 1 ], p[ 2 ] );
 329     pSource->SetNormal( n[ 0 ], n[ 1 ], n[ 2 ] );
 330     pSource->Update( );
 331
 332     _3Dslices.push_back( vtkProbeFilter::New( ) );
 333     _3Dslices[ k ]->SetInput( ( vtkDataSet* )pSource->GetOutput( ) );
 334     _3Dslices[ k ]->SetSource( vol->castVtk( ) );
 335     _3Dslices[ k ]->Update( );
 336     pSource->Delete( );
 337
 338     stPoints->GetPointData( )->
 339                 SetScalars( _3Dslices[ k ]->GetOutput( )->
 340                 GetPointData( )->GetScalars( ) );
 341     stPoints->SetDimensions( sizeIma, sizeIma, 1 );
 342     stPoints->SetScalarType( vol->castVtk( )->GetDataObjectType( ) );
 343     stPoints->Update( );
 344     if( forceCnt )
 345       {
 346                 marContour* tmp_marContour = new marContour( k, getParameters( ) ) ;
 347                 tmp_marContour->calculateVariables();
 348                 _contours.push_back( tmp_marContour );
 349         }
 350     _slices.push_back( new kVolume( kVolume::UCHAR, 1, 1, 1 ) );
 351
 352     *( _slices[ k ] ) = ( vtkImageData* )stPoints;  //copy
 353     _quantificationImages.push_back( (_slices[ k ])->castVtk( ) );
 354   } // rof
 355
 356   stPoints->Delete( );
 357 */
 358 }
 359
 360
 361
 362 // ----------------------------------------------------------------------------
 363 int marAxis::getNumberOfContours( )
 364 {
 365     return( _contours.size( ) );
 366 }
 367
 368 // ----------------------------------------------------------------------------
 369 int marAxis::getNumberOfSplinePoints( )
 370 {
 371     return( _points.size( ) );
 372 }
 373
 374 /*
 375 // EED Esto toca limpiarlo o reimplementarlo en varias partes o algo
 376 // ----------------------------------------------------------------------------
 377 vtkPolyData* marAxis::setContour( int i, int x, int y,  std::vector< double* >* points, marContour* c )
 378 {
 379         if( c != NULL )
 380         {
 381                 c->setS( c->getS( ) );
 382                 _contours[ i ]->copyFrom( *c );
 383         }
 384         else
 385         {
 386                 // PS -> //    int numberOfPoints, numberOfCells; // PS
 387                 // PS -> //    int id; // PS
 388                 // PS -> //    gslobj_vector p( 3 ); // PS
 389                 // PS -> //    double new_p[3]; // PS
 390                 // PS -> //    float pf[ 3 ]; // PS
 391                 vtkImageData* imagedata;
 392
 393                 // Contour parameters
 394                 // PS ->     int typ = getParameters( )->getIntParam( marParameters::e_algorithm_type );
 395                 double thr = getParameters( )->getDoubleParam( marParameters::e_threshold_isocontour );
 396
 397                 // PS ->        double thr = ( typ == marParameters::ISOCONTOURS )?
 398                 // PS ->       getParameters( )->
 399                 // PS ->       getDoubleParam( marParameters::e_threshold_isocontour ):
 400                 // PS ->     ( typ == marParameters::SNAKE ) ?
 401                 // PS ->       getParameters( )->
 402                 // PS ->       getDoubleParam( marParameters::e_threshold_snake_isocontour ):
 403                 // PS ->     30;
 404
 405
 406                 int sizeIma = getParameters( )->getSizeIma( );
 407                 double localint = getSignal( i );
 408                 double sigma = getParameters( )->getDoubleParam( marParameters::e_sigma );
 409                 double threshold = localint * thr / 100.0;
 410                 int dimIma    = ( int ) ( ( double ) sizeIma /getParameters( )->getDoubleParam( marParameters::e_scale ) );
 411                 int vmin = ( int )threshold;
 412                 int vmax = ( int )threshold;
 413
 414                 x = ( x != -1 )? x: sizeIma / 2;
 415                 y = ( y != -1 )? y: sizeIma / 2;
 416
 417                 // Isocontour
 418                 imagedata = (vtkImageData*) ((getSlice( i , NULL ))->castVtk( ));
 419                 float *range = imagedata->GetScalarRange();
 420
 421                 //vtkKitwareContourFilter* cntVTK = vtkKitwareContourFilter::New( );
 422                 vtkContourFilter* cntVTK = vtkContourFilter::New( );
 423                 cntVTK->SetInput( imagedata );
 424                 cntVTK->SetNumberOfContours( 1 );
 425                 //cntVTK->SetValue( 0, vmin );
 426                 cntVTK->SetValue( 0, (range[1]*thr/100) );
 427                 //cntVTK->SetValue( 1, vmax );
 428                 cntVTK->Update( );
 429
 430                 vtkCleanPolyData* cpd = vtkCleanPolyData::New( );
 431                 cpd->SetInput( cntVTK->GetOutput( ) );
 432                 cpd->ConvertLinesToPointsOff( );
 433                 cpd->Update( );
 434
 435                 vtkPolyDataConnectivityFilter* conn = vtkPolyDataConnectivityFilter::New( );
 436                 //conn->SetMaxRecursionDepth( 3000 );
 437
 438                 // PS ->     //conn->SetInput( cntVTK->GetOutput( ) ); PS
 439                 conn->SetInput( cpd->GetOutput( ) );
 440
 441                 conn->SetClosestPoint( x, y, 0 );
 442                 conn->SetExtractionModeToClosestPointRegion( );
 443                 conn->Update( );
 444
 445                 vtkCleanPolyData* cpd2 = vtkCleanPolyData::New( );
 446                 cpd2->SetInput( conn->GetOutput( ) );
 447                 cpd2->Update();
 448
 449                 vtkPolyData* polyDataResult = cpd2->GetOutput() ;
 450
 451                 double ta = _contours[ i ]->getS( );
 452
 453                 _contours[ i ]->reset( );
 454
 455                 // PS ->        gslobj_vector o( 3 ), n( 3 ), y( 3 ), vr( 3 );
 456
 457                 int id;
 458                 int numberOfCells;
 459                 int numberOfPoints,numberOfPoints2;
 460                 vtkCell* cell;
 461                 vtkIdList* pids;
 462                 marVector p( 3 );
 463 // PS ->                double new_p[3];
 464                 float pf[ 3 ];
 465
 466                 numberOfCells = polyDataResult->GetNumberOfCells( );
 467         cell = polyDataResult->GetCell( 0 );
 468                 pids = cell->GetPointIds( );
 469                 numberOfPoints = pids->GetNumberOfIds( );
 470
 471                 numberOfPoints2=polyDataResult->GetNumberOfPoints();
 472
 473                 marVector o(3),n(3),y(3),vr(3);
 474                 double cost, sint;
 475                 y( 0 ) = y( 2 ) = 0.0;
 476                 y( 1 ) = 1.0;
 477                 getPoint( ( double* )o, ta );
 478                 getTangent( ( double* )n, ta );
 479                 vr = y.cross( n ).normalize( );
 480                 sint = y.cross( n ).norm2( );
 481                 cost = y.dot( n );
 482
 483                 for( int j = 0; j < numberOfPoints; j++ )
 484                 {
 485                         id = pids->GetId( j );
 486                         polyDataResult->GetPoint( id, pf );
 487                         p( 0 ) = pf[ 0 ];
 488                         p( 1 ) = 0.0;
 489                         p( 2 ) = pf[ 1 ];
 490
 491                         // 3D transformation
 492                         p = ( p * cost ) + ( vr * ( vr.dot( p ) * ( 1.0 - cost ) ) ) -
 493                                 ( p.cross( vr ) * sint ) + o;
 494
 495                 // PS ->                double p_cross_vr[3];
 496                 // PS ->                vtkMath::Cross( p, vr, p_cross_vr);
 497                 // PS ->                for(int k=0; k<3; k++)
 498                 // PS ->                {
 499                 // PS ->                        new_p[i] = (new_p[k] * cost) + ( vr(k) * ( vtkMath::Dot(vr, p) * ( 1.0 - cost ) ) ) -
 500                 // PS ->                                ( p_cross_vr[k] * sint ) + o(i);
 501                 // PS ->                }
 502
 503                         //TODO: check if p and new_p are the same...
 504                         _contours[ i ]->addContourPoint( p );
 505                 } // rof
 506
 507                 // PS ->        int nbPoints=polyDataResult->GetNumberOfPoints();
 508                 // PS ->
 509                 // PS ->        float* vtkContourPoint;
 510                 // PS ->        vtkContourPoint= new float[3];
 511                 // PS ->        double* marContourPoint;
 512                 // PS ->        marContourPoint= new double[3];
 513                 // PS ->
 514                 // PS ->        for (int j=0;j< nbPoints;j++)
 515                 // PS ->        {
 516                 // PS ->                vtkContourPoint=polyDataResult->GetPoint(j);
 517                 // PS ->                marContourPoint[0]=(double)vtkContourPoint[0];
 518                 // PS ->                marContourPoint[1]=(double)vtkContourPoint[1];
 519                 // PS ->                marContourPoint[2]=(double)vtkContourPoint[2];
 520                 // PS ->                _contours[ i ]->addContourPoint( marContourPoint );
 521                 // PS ->        }
 522                 marContour * tutu=_contours[ i ];
 523                 _contours[ i ]->getArea();
 524                 _contours[ i ]->getPerimeter();
 525                 return( polyDataResult );
 526                 //return( cpd->GetOutput( ) );
 527                 //return( cntVTK->GetOutput( ) );
 528
 529
 530 //EED1             vtkStripper* isoStrips = vtkStripper::New( );
 531 //EED1             isoStrips->SetInput( cpd2->GetOutput( ) );
 532 //EED1             isoStrips->Update( );
 533 //EED1
 534 //EED1                   vtkPolyData* data = isoStrips->GetOutput( );
 535 //EED1                   data->Update( );
 536 //EED1                   vtkCell* cell;
 537 //EED1                   vtkIdList* pids;
 538 //EED1
 539 //EED1                     numberOfCells = data->GetNumberOfCells( );
 540 //EED1                     cell = data->GetCell( 0 );
 541 //EED1                     pids = cell->GetPointIds( );
 542 //EED1                     numberOfPoints = pids->GetNumberOfIds( );
 543
 544
 545 //EED2          _contours[ i ]->reset( );
 546 //EED2          if( typ == marParameters::ISOCONTOURS ) {
 547 //EED2          gslobj_vector o( 3 ), n( 3 ), y( 3 ), vr( 3 );
 548 //EED2          double cost, sint;
 549 //EED2          y( 0 ) = y( 2 ) = 0.0;
 550 //EED2          y( 1 ) = 1.0;
 551 //EED2          getPoint( ( double* )o, ta );
 552 //EED2          getTangent( ( double* )n, ta );
 553 //EED2          vr = y.cross( n ).normalize( );
 554 //EED2          sint = y.cross( n ).norm2( );
 555 //EED2          cost = y.dot( n );
 556 //EED2
 557 //EED2            for( int j = 0; j < numberOfPoints; j++ ) {
 558 //EED2            id = pids->GetId( j );
 559 //EED2            data->GetPoint( id, pf );
 560 //EED2            p( 0 ) = pf[ 0 ];
 561 //EED2            p( 1 ) = 0.0;
 562 //EED2            p( 2 ) = pf[ 1 ];
 563 //EED2
 564 //EED2                  // 3D transformation
 565 //EED2                  p = ( p * cost ) + ( vr * ( vr.dot( p ) * ( 1.0 - cost ) ) ) -
 566 //EED2          ( p.cross( vr ) * sint ) + o;
 567 //EED2
 568 //EED2                    double p_cross_vr[3];
 569 //EED2                    vtkMath::Cross( p, vr, p_cross_vr);
 570 //EED2                    for(int ii=0; ii<3; ii++) {
 571 //EED2                    new_p[i] = (new_p[ii] * cost) + ( vr(ii) * ( vtkMath::Dot(vr, p) * ( 1.0 - cost ) ) ) -
 572 //EED2                    ( p_cross_vr[ii] * sint ) + o(i);
 573 //EED2                    }
 574 //EED2
 575 //EED2                          //TODO: check if p and new_p are the same...
 576 //EED2                          _contours[ i ]->addContourPoint( p );
 577 //EED2                          } // rof
 578 //EED2                          } else if( typ == marParameters::SNAKE || typ == marParameters::DERICHE ) {
 579 //EED2                          // ONLY Isocontours at the moment.
 580 //EED2                          } // fi
 581 //EED2
 582 //EED2                            cntVTK->Delete( );
 583 //EED2                            cpd->Delete( );
 584 //EED2                            conn->Delete( );
 585 //EED2                            cpd2->Delete( );
 586 //EED2          isoStrips->Delete( );
 587
 588   } // fi
 589 }
 590
 591 */
 592
 593
 594
 595
 596 // ----------------------------------------------------------------------------
 597 void marAxis::createContour(int i, kVolume* vol)
 598 {
 599
 600         int j,id;
 601         int numberOfPoints,numberOfCells;
 602         vtkCell* cell;
 603         vtkIdList* pids;
 604         double p[ 3 ];
 605         vtkPolyData* vtkPolyData_2DContour;
 606         double x,y;
 607
 608         int sizeIma                             =       getParameters( )->getSizeIma( );
 609         double dimIma                   =       getParameters( )->getDimIma( );
 610
 611 // EED 24 Dec 2007
 612   dimIma=sizeIma/3;
 613
 614         _contours[i]                    =       new marContour( i, getParameters( ) ) ;
 615         vtkPolyData_2DContour   =       get2Dcontour(i,vol);
 616
 617         numberOfCells                   =       vtkPolyData_2DContour->GetNumberOfCells( );
 618     cell                                        =       vtkPolyData_2DContour->GetCell( 0 );
 619         pids                                    =       cell->GetPointIds( );
 620         numberOfPoints                  =       pids->GetNumberOfIds( );
 621
 622         for( j = 0; j < numberOfPoints; j++ )
 623         {
 624                 id = pids->GetId( j );
 625                 vtkPolyData_2DContour->GetPoint( id, p);
 626                 x=p[0]-64.0;
 627                 y=p[1]-64.0;
 628                 x=x * dimIma / ( float ) sizeIma;
 629                 y=y * dimIma / ( float ) sizeIma;
 630                 _contours[ i ]->addContourPoint( x , y );
 631         }
 632
 633         _contours[i]->do_spline();
 634         _contours[i]->calculateVariables();
 635
 636 }
 637
 638 // ----------------------------------------------------------------------------
 639 void marAxis::EraseContour(int i){
 640         if (_3Dcontour[i]!=NULL){
 641                 _3Dcontour[i]->Delete();
 642                 _3Dcontour[i]=NULL;
 643         }
 644
 645         if (_2DDiameterMax[i]!=NULL){
 646                 _2DDiameterMax[i]->Delete();
 647                 _2DDiameterMax[i]=NULL;
 648         }
 649
 650         if (_2DDiameterMin[i]!=NULL){
 651                 _2DDiameterMin[i]->Delete();
 652                 _2DDiameterMin[i]=NULL;
 653         }
 654
 655         if (_2Dcontours[i]!=NULL){
 656                 _2Dcontours[i]->Delete();
 657                 _2Dcontours[i]=NULL;
 658         }
 659
 660         if (_contours[i]!=NULL){
 661                 delete _contours[i];
 662                 _contours[i]=NULL;
 663         }
 664 }
 665
 666 // ----------------------------------------------------------------------------
 667 void marAxis::replaceContour2D(int i,int size,double *vx,double *vy){
 668
 669         EraseContour(i);
 670
 671
 672         vtkPoints *_pts = vtkPoints::New();
 673         _pts->SetNumberOfPoints(size);
 674         int j;
 675
 676         for (j=0 ; j<size ; j++){
 677                 _pts->SetPoint(j,       vx[j]   , vy[j] , 0 );
 678         }
 679 //      _pts->SetPoint(0,       vx[0]   , vy[0] , 0 );
 680
 681         vtkCellArray *lines = vtkCellArray::New();
 682         lines->InsertNextCell( size );
 683         for ( j=0 ; j<size+1 ; j++ ){
 684                 lines->InsertCellPoint(j % size );
 685         }
 686
 687         vtkPolyData *_pd = vtkPolyData::New();
 688         _pd->SetPoints( _pts );
 689         _pd->SetLines( lines );
 690         lines->Delete();  //do not delete lines ??
 691         _pts->Delete();
 692
 693         _2Dcontours[i]=_pd;
 694         createContour(i,NULL);
 695 }
 696 // ----------------------------------------------------------------------------
 697 void marAxis::createSlice(int i , kVolume* vol){
 698         int k=i;
 699         int sizeIma = getParameters( )->getSizeIma( );
 700
 701         vtkStructuredPoints* stPoints = vtkStructuredPoints::New( );
 702         double *p, *n;
 703     p = _points[k];
 704     n = getNormal( k );
 705
 706     stPoints->GetPointData( )->SetScalars(  get3DSlice(k,vol)->GetOutput()->GetPointData()->GetScalars()  );
 707     stPoints->SetDimensions( sizeIma, sizeIma, 1 );
 708
 709
 710     stPoints->SetScalarType( vol->castVtk( )->GetDataObjectType( ) );
 711         stPoints->SetScalarTypeToShort ();
 712     stPoints->Update();
 713
 714     vtkImageChangeInformation* change = vtkImageChangeInformation::New();
 715     change->SetInput( stPoints );
 716     change->Update();    //important
 717
 718
 719     if (_slices[k]!=NULL) { delete _slices[k]; }
 720     _slices[k] = new kVolume( change->GetOutput() ) ;
 721
 722         stPoints->Delete( );
 723         change->Delete();
 724 }
 725 // ----------------------------------------------------------------------------
 726 void marAxis::create3DSlice(int i , kVolume* vol )
 727 {
 728
 729   int k=i;
 730   int sizeIma   = getParameters( )->getSizeIma( );
 731   double dimIma = getParameters( )->getDimIma( );
 732
 733   dimIma=sizeIma/3;
 734
 735
 736 //  double voxSize = getParameters( )->getVoxelSize( );
 737
 738 //  if( forceCnt )
 739 //    for( unsigned int i = 0; i < _contours.size( ); i++ )
 740 //      delete _contours[ i ];
 741 //    _contours.clear( );
 742
 743 //    for( unsigned int i = 0; i < _slices.size( ); i++ ) {
 744 //      delete _slices[ i ];
 745 //      _3Dslices[ i ]->Delete( );
 746 //      _quantificationImages[ i ]->Delete( );
 747 //    }
 748
 749 //  _slices.clear( );
 750 // _3Dslices.clear( );
 751 //  _quantificationImages.clear( );
 752
 753   // "Cutter" initialization
 754 //  vtkStructuredPoints* stPoints = vtkStructuredPoints::New( );
 755
 756   // Perform "cut"
 757   double *p, *n;
 758 //  for( int k = 0; k < nCnts; k++ ) {
 759     //s = ( double )k / ( double )( nCnts - 1 );
 760     p = _points[k];
 761     n = getNormal( k );
 762
 763     vtkPlaneSource* pSource = vtkPlaneSource::New( );
 764     pSource->SetOrigin( p[ 0 ], p[ 1 ], p[ 2 ] );
 765     pSource->SetPoint1( p[ 0 ] + dimIma - 1.0, p[ 1 ], p[ 2 ] );
 766     pSource->SetPoint2( p[ 0 ], p[ 1 ], p[ 2 ] + dimIma - 1.0 );
 767     pSource->SetResolution( sizeIma - 1, sizeIma - 1 );
 768     pSource->Update( );
 769     pSource->SetCenter( p[ 0 ], p[ 1 ], p[ 2 ] );
 770     pSource->SetNormal( n[ 0 ], n[ 1 ], n[ 2 ] );
 771     pSource->Update( );
 772
 773         if (_3Dslices[ k ]!=NULL){ _3Dslices[ k ]->Delete(); }
 774     _3Dslices[ k ] = vtkProbeFilter::New( ) ;
 775     _3Dslices[ k ]->SetInput( ( vtkDataSet* )pSource->GetOutput( ) );
 776     _3Dslices[ k ]->SetSource( vol->castVtk( ) );
 777     _3Dslices[ k ]->Update( );
 778     pSource->Delete( );
 779
 780 //    stPoints->GetPointData( )->
 781 //              SetScalars( _3Dslices[ k ]->GetOutput( )->
 782 //              GetPointData( )->GetScalars( ) );
 783 //    stPoints->SetDimensions( sizeIma, sizeIma, 1 );
 784 //    stPoints->SetScalarType( vol->castVtk( )->GetDataObjectType( ) );
 785 //    stPoints->Update( );
 786 //
 787 //    if( forceCnt )
 788 //      {
 789 //              marContour* tmp_marContour = new marContour( k, getParameters( ) ) ;
 790 //              tmp_marContour->calculateVariables();
 791 //              _contours.push_back( tmp_marContour );
 792 //      }
 793 //    _slices.push_back( new kVolume( kVolume::UCHAR, 1, 1, 1 ) );
 794 //    *( _slices[ k ] ) = ( vtkImageData* )stPoints;  //copy
 795 //    _quantificationImages.push_back( (_slices[ k ])->castVtk( ) );
 796 //  } // rof
 797
 798 //  stPoints->Delete( );
 799
 800 //      _3Dslices[i]=NULL;
 801 }
 802
 803
 804
 805
 806 // ----------------------------------------------------------------------------
 807 void marAxis::create3Dcontour(int i, kVolume* vol)
 808 {
 809         vtkPoints *_vtkPoints;
 810         double *o;
 811         double *n;
 812         vtkMatrix4x4* mat;
 813         vtkTransform* trans;
 814         double nn[3];
 815         double c[3],p1[3],p2[3];
 816         double nA,nB,nC;
 817
 818         _vtkPoints                                              = vtkPoints::New();
 819         o                                                               = _points[i];
 820         n                                                               = getNormal( i );
 821         mat                                                             = vtkMatrix4x4::New();
 822         trans                                                   = vtkTransform::New();
 823
 824
 825         nn[0]=n[0];  nn[1]=n[1];  nn[2]=n[2];
 826         nC=sqrt( nn[0]*nn[0] + nn[1]*nn[1] + nn[2]*nn[2] );
 827         nn[0]=nn[0]/nC;  nn[1]=nn[1]/nC;  nn[2]=nn[2]/nC;
 828
 829     vtkPlaneSource* pSource = vtkPlaneSource::New( );
 830     pSource->SetOrigin( 0, 0    , 0                             );
 831     pSource->SetPoint1( 1, 0    , 0                             );
 832     pSource->SetPoint2( 0, 0    , 1.0   );
 833     pSource->SetCenter( 0,0,0 );
 834     pSource->SetNormal( nn[ 0 ], nn[ 1 ], nn[ 2 ] );
 835     pSource->Update( );
 836 //    pSource->Update( );
 837         pSource->GetOrigin(c);
 838         pSource->GetPoint1(p1);
 839         pSource->GetPoint2(p2);
 840         pSource->Delete();
 841         p1[0]=p1[0]-c[0];  p1[1]=p1[1]-c[1];  p1[2]=p1[2]-c[2];
 842         p2[0]=p2[0]-c[0];  p2[1]=p2[1]-c[1];  p2[2]=p2[2]-c[2];
 843         nA=sqrt( p1[0]*p1[0] + p1[1]*p1[1] + p1[2]*p1[2] );
 844         nB=sqrt( p2[0]*p2[0] + p2[1]*p2[1] + p2[2]*p2[2] );
 845         p1[0]=p1[0]/nA;  p1[1]=p1[1]/nA;  p1[2]=p1[2]/nA;
 846         p2[0]=p2[0]/nB;  p2[1]=p2[1]/nB;  p2[2]=p2[2]/nB;
 847
 848     mat->SetElement (0,0, nn[0]);
 849     mat->SetElement (1,0, nn[1]);
 850     mat->SetElement (2,0, nn[2]);
 851     mat->SetElement (3,0, 0);
 852     mat->SetElement (0,1, p2[0]);
 853     mat->SetElement (1,1, p2[1]);
 854     mat->SetElement (2,1, p2[2]);
 855     mat->SetElement (3,1, 0);
 856         mat->SetElement (0,2, p1[0]);
 857     mat->SetElement (1,2, p1[1]);
 858     mat->SetElement (2,2, p1[2]);
 859     mat->SetElement (3,2, 0);
 860         mat->SetElement (0,3, 0);
 861     mat->SetElement (1,3, 0);
 862     mat->SetElement (2,3, 0);
 863     mat->SetElement (3,3, 1);
 864
 865 //      double deter=mat->Determinant(mat);
 866         trans->Identity();
 867         trans->SetMatrix(mat);
 868         float ang;
 869         ang=-90;
 870         trans->RotateWXYZ  ( ang,0,1,0);
 871
 872 // EED Borrame
 873 //      double scale = 1;
 874 //      trans->Scale  ( scale , scale, scale);
 875
 876         trans->Update();
 877
 878         vtkMatrix4x4 *m=vtkMatrix4x4::New();
 879         trans->GetMatrix  (  m  );
 880
 881
 882         int j,numberOfPoints;
 883         marContour* contour2D   = getContour(i,vol);
 884         numberOfPoints                  = contour2D->getNumberOfPoints( );
 885         double fpA[3];
 886         double pA[3],ppp[3];
 887
 888
 889         for( j = 0; j <= numberOfPoints; j++ ) {
 890                 contour2D->getPoint( fpA,j % numberOfPoints);
 891
 892                 pA[0] = fpA[0] + 0.2;  // correction EED
 893                 pA[1] = fpA[1] + 0.2;  // correction EED
 894
 895                 pA[2] = 0;
 896
 897 // Why this does not works...(release version)????
 898 //              trans->TransformPoint(pA,ppp);
 899 // or
 900 //              trans->TransformVector(pA,ppp);
 901 // so..
 902                 ppp[0] = m->GetElement(0,0)*pA[0] + m->GetElement(0,1)*pA[1] + m->GetElement(0,2)*pA[2] + m->GetElement(0,3);
 903                 ppp[1] = m->GetElement(1,0)*pA[0] + m->GetElement(1,1)*pA[1] + m->GetElement(1,2)*pA[2] + m->GetElement(1,3);
 904                 ppp[2] = m->GetElement(2,0)*pA[0] + m->GetElement(2,1)*pA[1] + m->GetElement(2,2)*pA[2] + m->GetElement(2,3);
 905
 906 //              ppp[0]=ppp[0]+o[0]-c[0]; ppp[1]=ppp[1]+o[1]-c[1]; ppp[2]=ppp[2]+o[2]-c[2];
 907                 ppp[0]=ppp[0]+o[0]; ppp[1]=ppp[1]+o[1]; ppp[2]=ppp[2]+o[2];
 908 //              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;;
 909                 _vtkPoints->InsertNextPoint( (float)ppp[0],(float)ppp[1],(float)ppp[2] );
 910         }
 911
 912         m->Delete();
 913         mat->Delete();
 914         trans->Delete();
 915         if (_3Dcontour[i]!=NULL){ _3Dcontour[i]->Delete(); }
 916         _3Dcontour[i]=_vtkPoints;
 917
 918 }
 919 // ----------------------------------------------------------------------------
 920 void marAxis::createSliceImage(int i  , kVolume* vol )
 921 {
 922   _quantificationImages[i] = getSlice(i,vol)->castVtk( ) ;
 923
 924 }
 925 // ----------------------------------------------------------------------------
 926 void marAxis::create2Dcontour(int i , kVolume* vol )
 927 {
 928         vtkImageData* imagedata;
 929
 930         double  thr                     = getParameters( )->getDoubleParam( marParameters::e_threshold_isocontour );
 931         int             sizeIma         = getParameters( )->getSizeIma( );
 932         double  localint        = getSignal( i );
 933         double  sigma           = getParameters( )->getDoubleParam( marParameters::e_sigma );
 934         double  threshold       = localint * thr / 100.0;
 935
 936         int             dimIma          = ( int ) ( ( double ) sizeIma /getParameters( )->getDoubleParam( marParameters::e_scale ) );
 937         int             vmin            = ( int )threshold;
 938         int             vmax            = ( int )threshold;
 939
 940         int             x                       = -1;
 941         int             y                       = -1;
 942
 943         x = ( x != -1 )? x: sizeIma / 2;
 944         y = ( y != -1 )? y: sizeIma / 2;
 945
 946         // Isocontour
 947         imagedata = (vtkImageData*) ((getSlice( i , vol ))->castVtk( ));
 948         double *range = imagedata->GetScalarRange();
 949
 950         //vtkKitwareContourFilter* cntVTK = vtkKitwareContourFilter::New( );
 951         vtkContourFilter* cntVTK = vtkContourFilter::New( );
 952         cntVTK->SetInput( imagedata );
 953         cntVTK->SetNumberOfContours( 1 );
 954         //cntVTK->SetValue( 0, vmin );
 955         cntVTK->SetValue( 0, (range[1]*thr/100) );
 956         //cntVTK->SetValue( 1, vmax );
 957         cntVTK->Update( );
 958
 959         vtkCleanPolyData* cpd = vtkCleanPolyData::New( );
 960         cpd->SetInput( cntVTK->GetOutput( ) );
 961         cpd->ConvertLinesToPointsOff( );
 962         cpd->Update( );
 963
 964         vtkPolyDataConnectivityFilter* conn = vtkPolyDataConnectivityFilter::New( );
 965         //conn->SetMaxRecursionDepth( 3000 );
 966
 967         // PS ->     //conn->SetInput( cntVTK->GetOutput( ) ); PS
 968         conn->SetInput( cpd->GetOutput( ) );
 969
 970         conn->SetClosestPoint( x, y, 0 );
 971         conn->SetExtractionModeToClosestPointRegion( );
 972         conn->Update( );
 973
 974         vtkCleanPolyData* cpd2 = vtkCleanPolyData::New( );
 975         cpd2->SetInput( conn->GetOutput( ) );
 976         cpd2->Update();
 977
 978         vtkStripper* vtkstripper = vtkStripper::New( );
 979         vtkstripper->SetInput( cpd2->GetOutput() );
 980         vtkstripper->Update();
 981
 982         vtkPolyData* polyDataResult =  vtkstripper->GetOutput();
 983     polyDataResult->Update( );
 984
 985         cntVTK          -> Delete();
 986         cpd2            -> Delete();
 987         cpd                     -> Delete();
 988         conn            -> Delete();
 989
 990 /*
 991 // EED 27 Oct 2007
 992         double *p;
 993         int ii,size=polyDataResult->GetNumberOfPoints();
 994 printf("EED marAxis::create2Dcontour size %d\n",size);
 995         for (ii=0;ii<size;ii++)
 996         {
 997                 p=polyDataResult->GetPoint(ii);
 998 printf("EED marAxis::create2Dcontour %d >> %f %f %f\n",ii,p[0], p[1], p[2]);
 999                 p[2]=-100;
1000         }
1001
1002         for (ii=0;ii<size;ii++)
1003         {
1004                 p=polyDataResult->GetPoint(ii);
1005 printf("EED marAxis::create2Dcontour xxxxx %d >> %f %f %f\n",ii,p[0], p[1], p[2]);
1006         }
1007 */
1008
1009         if ( _2Dcontours[i]!=NULL ) { _2Dcontours[i]->Delete(); }
1010         _2Dcontours[i]=polyDataResult;
1011 }
1012 // ----------------------------------------------------------------------------
1013 void marAxis::create2DDiameterMin(      int i , kVolume* vol ){
1014         double p1[2],p2[2];
1015         marContour *marcontour=getContour(i,vol);
1016         marcontour->getMinimumLine(p1,p2);
1017         vtkPoints *_vtkPoints = vtkPoints::New();
1018         int sizeIma = getParameters( )->getSizeIma( );
1019
1020         double dimIma = getParameters( )->getDimIma( );
1021
1022 // EED 24 Dec 2007
1023   dimIma=sizeIma/3;
1024
1025         double factor=( float ) sizeIma / dimIma;
1026         p1[0]=p1[0]*factor+sizeIma/2;
1027         p1[1]=p1[1]*factor+sizeIma/2;
1028         p2[0]=p2[0]*factor+sizeIma/2;
1029         p2[1]=p2[1]*factor+sizeIma/2;
1030         _vtkPoints->InsertNextPoint( p1[0] , p1[1] , 1 );
1031         _vtkPoints->InsertNextPoint( p2[0] , p2[1] , 1 );
1032         if ( _2DDiameterMin[i]!=NULL ) { _2DDiameterMin[i]->Delete(); }
1033         _2DDiameterMin[i] = _vtkPoints;
1034 }
1035 // ----------------------------------------------------------------------------
1036 void marAxis::create2DDiameterMax(      int i , kVolume* vol ){
1037         double p1[2],p2[2];
1038         marContour *marcontour=getContour(i,vol);
1039         marcontour->getMaximumLine(p1,p2);
1040         vtkPoints *_vtkPoints = vtkPoints::New();
1041         int sizeIma = getParameters( )->getSizeIma( );
1042
1043         double dimIma = getParameters( )->getDimIma( );
1044 // EED 24 Dec 2007
1045 dimIma=sizeIma/3;
1046
1047         double factor=( float ) sizeIma / dimIma;
1048         p1[0]=p1[0]*factor+sizeIma/2;
1049         p1[1]=p1[1]*factor+sizeIma/2;
1050         p2[0]=p2[0]*factor+sizeIma/2;
1051         p2[1]=p2[1]*factor+sizeIma/2;
1052         _vtkPoints->InsertNextPoint( p1[0] , p1[1] , 1 );
1053         _vtkPoints->InsertNextPoint( p2[0] , p2[1] , 1 );
1054         if ( _2DDiameterMax[i]!=NULL ) { _2DDiameterMax[i]->Delete(); }
1055         _2DDiameterMax[i] = _vtkPoints;
1056 }
1057 // ----------------------------------------------------------------------------
1058 void marAxis::createEmptyVectors()
1059 {
1060         int nCnts = ( int ) getNumberOfSplinePoints( );
1061         clearAllVectors();
1062         int i;
1063 //      int xxxx=getNumberOfContours();
1064         for (i=0;i<nCnts;i++){
1065                 _contours.push_back( NULL );
1066                 _slices.push_back( NULL );
1067                 _3Dslices.push_back( NULL );
1068                 _3Dcontour.push_back( NULL );
1069                 _quantificationImages.push_back( NULL );
1070                 _2Dcontours.push_back( NULL );
1071                 _2DDiameterMin.push_back( NULL );
1072                 _2DDiameterMax.push_back( NULL );
1073         }
1074 }
1075 // ----------------------------------------------------------------------------
1076 void marAxis::clearAllVectors()
1077 {
1078         unsigned int i;
1079
1080         for (i=0;i<_contours.size();i++){
1081                 if (_contours[i]!=NULL) {                                       //  DATA-MODEL-2D
1082                         delete _contours[i];
1083                         _contours[i]=NULL;
1084                 }
1085         }
1086         _contours.clear();
1087
1088         for (i=0;i<_slices.size();i++){                                 //  DATA-MODEL-Voxel XxYx1
1089                 if (_slices[i]!=NULL) {
1090                         delete _slices[i];
1091                         _slices[i]=NULL;
1092                 }
1093         }
1094         _slices.clear();
1095
1096
1097         for (i=0;i<_3Dslices.size();i++){                                       //  VISUALISATION_VTK 3D
1098                 if (_3Dslices[i]!=NULL) {
1099                         _3Dslices[i]->Delete();
1100                         _3Dslices[i]=NULL;
1101                 }
1102         }
1103         _3Dslices.clear();
1104
1105         for (i=0;i<_3Dcontour.size();i++){                              //  VISUALISATION_VTK 3D
1106                 if (_3Dcontour[i]!=NULL) {
1107                         _3Dcontour[i]->Delete();
1108                         _3Dcontour[i]=NULL;
1109                 }
1110         }
1111         _3Dcontour.clear();
1112
1113         for (i=0;i<_quantificationImages.size();i++){   //  VISUALISATION_VTK 2D
1114                 if (_quantificationImages[i]!=NULL) {
1115 //                      _quantificationImages[i]->Delete();
1116 //                      _quantificationImages[i]=NULL;
1117                 }
1118         }
1119         _quantificationImages.clear();
1120
1121         for (i=0;i<_2Dcontours.size();i++){                             //  VISUALISATION_VTK 2D
1122                 if (_2Dcontours[i]!=NULL) {
1123 // EED ???  This object was allready erased bye the VTK pipeline
1124 //                      _2Dcontours[i]->Delete();
1125                         _2Dcontours[i]=NULL;
1126                 }
1127         }
1128         _2Dcontours.clear();
1129
1130         for (i=0;i<_2DDiameterMin.size();i++){                  //  VISUALISATION_VTK 2D
1131                 if (_2DDiameterMin[i]!=NULL) {
1132                         _2DDiameterMin[i]->Delete() ;
1133                         _2DDiameterMin[i]=NULL;
1134                 }
1135         }
1136         _2DDiameterMin.clear();
1137
1138         for (i=0;i<_2DDiameterMax.size();i++){                  //  VISUALISATION_VTK 2D
1139                 if (_2DDiameterMax[i]!=NULL) {
1140                         _2DDiameterMax[i]->Delete();
1141                         _2DDiameterMax[i]=NULL;
1142                 }
1143         }
1144         _2DDiameterMax.clear();
1145
1146 }
1147 // ----------------------------------------------------------------------------
1148 void marAxis::eraseContourVectorsContent()
1149 {
1150         unsigned int i;
1151         for (i=0;i<_contours.size();i++){
1152                 if (_contours[i]!=NULL) {                                       //  DATA-MODEL-2D
1153                         delete _contours[i];
1154                         _contours[i]=NULL;
1155                 }
1156         }
1157         for (i=0;i<_2Dcontours.size();i++){                             //  VISUALISATION_VTK 2D
1158                 if (_2Dcontours[i]!=NULL) {
1159                         _2Dcontours[i]->Delete();
1160                         _2Dcontours[i]=NULL;
1161                 }
1162         }
1163         for (i=0;i<_2DDiameterMin.size();i++){                          //  VISUALISATION_VTK 2D
1164                 if (_2DDiameterMin[i]!=NULL) {
1165                         _2DDiameterMin[i]->Delete();
1166                         _2DDiameterMin[i]=NULL;
1167                 }
1168         }
1169         for (i=0;i<_2DDiameterMax.size();i++){                          //  VISUALISATION_VTK 2D
1170                 if (_2DDiameterMax[i]!=NULL) {
1171                         _2DDiameterMax[i]->Delete();
1172                         _2DDiameterMax[i]=NULL;
1173                 }
1174         }
1175         for (i=0;i<_3Dcontour.size();i++){                              //  VISUALISATION_VTK 3D
1176                 if (_3Dcontour[i]!=NULL) {
1177                         _3Dcontour[i]->Delete();
1178                         _3Dcontour[i]=NULL;
1179                 }
1180         }
1181 }
1182
1183 // ----------------------------------------------------------------------------
1184 marContour* marAxis::getContour(int i, kVolume* vol)                            // DATA-MODEL 2D
1185 {
1186         if (_contours[i]==NULL){ createContour(i, vol); }
1187         return _contours[i];
1188 }
1189 // ----------------------------------------------------------------------------
1190 kVolume* marAxis::getSlice(int i, kVolume* vol)                                         // DATA-MODEL-Voxel XxYx1
1191 {
1192         if (_slices[i]==NULL){ createSlice( i , vol ); }
1193         return _slices[i];
1194 }
1195 // ----------------------------------------------------------------------------
1196
1197 bool marAxis::if3DcontourExist(int i)
1198 {
1199         bool result=true;
1200         if (_3Dcontour[i]==NULL)
1201         {
1202                 result=false;
1203         }
1204         return result;
1205 }
1206
1207 // ----------------------------------------------------------------------------
1208 void marAxis::Save3Dcontour(FILE *ff,int id)
1209 {
1210         int i,size;
1211         double point[3];
1212         if (_3Dcontour[id]!=NULL)
1213         {
1214 //Ramiro                fprintf(ff,"contour_id: %d \n", id);
1215                 size = _3Dcontour[id]->GetNumberOfPoints();
1216 //Ramiro                fprintf(ff,"numberOfPoints: %d \n", size);
1217             for ( i=0 ; i<size  ;i++ ){
1218                 _3Dcontour[id]->GetPoint( i , point );
1219                         fprintf(ff,"%f %f %f %d\n", point[0], point[1],point[2],id);
1220             }
1221         }
1222 }
1223 // ----------------------------------------------------------------------------
1224 void marAxis::SaveExisting3DContours(FILE *ff)
1225 {
1226         if (ff!=NULL)
1227         {
1228                 // counting existing 3Dcontours
1229                 int acum=0;
1230                 int i,size=_3Dcontour.size();
1231                 for (i=0; i<size; i++)
1232                 {
1233                         if (_3Dcontour[i]!=NULL)
1234                         {
1235                                 acum++;
1236                         }
1237                 }
1238 //Ramiro                fprintf(ff, "NumberOf3DContours: %d \n", acum);
1239                 // saving existing 3Dcontours
1240                 for (i=0; i<size; i++)
1241                 {
1242                         Save3Dcontour(ff,i);
1243                 }
1244         }
1245 }
1246
1247 // ----------------------------------------------------------------------------
1248 vtkProbeFilter* marAxis::get3DSlice(int i, kVolume* vol){                       // VISUALISATION-VTK 3D
1249         if (_3Dslices[i]==NULL){ create3DSlice(i,vol); }
1250         return _3Dslices[i];
1251 }
1252 // ----------------------------------------------------------------------------
1253 vtkPoints* marAxis::get3Dcontour(int i, kVolume* vol){                          // VISUALISATION-VTK 3D
1254         if (_3Dcontour[i]==NULL){ create3Dcontour(i, vol); }
1255         return _3Dcontour[i];
1256 }
1257 // ----------------------------------------------------------------------------
1258 vtkImageData* marAxis::getSliceImage(int i, kVolume* vol){                      // VISUALISATION-VTK 2D
1259         if (_quantificationImages[i]==NULL){ createSliceImage(i,vol); }
1260         return _quantificationImages[i];
1261 }
1262 // ----------------------------------------------------------------------------
1263 vtkPolyData* marAxis::get2Dcontour(int i , kVolume* vol){                       // VISUALISATION-VTK 2D
1264         if (_2Dcontours[i]==NULL){ create2Dcontour(i,vol); }
1265         return _2Dcontours[i];
1266 }
1267 // ----------------------------------------------------------------------------
1268 vtkPoints* marAxis::get2DDiameterMin(int i , kVolume* vol){             // VISUALISATION-VTK 2D
1269         if (_2DDiameterMin[i]==NULL){ create2DDiameterMin(i,vol); }
1270         return _2DDiameterMin[i];
1271 }
1272 // ----------------------------------------------------------------------------
1273 vtkPoints* marAxis::get2DDiameterMax(int i , kVolume* vol){             // VISUALISATION-VTK 2D
1274         if (_2DDiameterMax[i]==NULL){ create2DDiameterMax(i,vol); }
1275         return _2DDiameterMax[i];
1276 }
1277 // ----------------------------------------------------------------------------
1278 double marAxis::getAxisLenght(int pIni,int pEnd){
1279
1280     if (pIni>pEnd){
1281                 int tmp=pIni;
1282                 pIni=pEnd;
1283                 pEnd=tmp;
1284         }
1285
1286         marVector* pO = new marVector(2);
1287         marVector* pF = new marVector(2);
1288         double L;
1289     unsigned int j;
1290     for( j = pIni, L = 0.0; j < pEnd-1; j++ ) {
1291                 (*pO)=_points[j];
1292                 (*pF)=_points[j+1];
1293                 (*pF)=(*pF)-(*pO);
1294                 L+=pF->norm2();
1295     } // rof
1296         delete pO;
1297         delete pF;
1298         return L;
1299 }
1300 // ----------------------------------------------------------------------------
1301 void marAxis::calculateTotalAxisLenght(){
1302         _totalAxisLenght = getAxisLenght( 0 , _points.size( ) );
1303 }
1304 // ----------------------------------------------------------------------------
1305 double marAxis::getTotalLength(){
1306         if (_totalAxisLenght==-1){
1307                 calculateTotalAxisLenght();
1308         }
1309         return _totalAxisLenght;
1310 }
1311
1312 // ----------------------------------------------------------------------------
1313 void marAxis::calculateSubAxisLength(){
1314         _subAxisLenght = getAxisLenght( _startQuant , _finishQuant );
1315 }
1316 // ----------------------------------------------------------------------------
1317 double marAxis::getSubAxisLength(){
1318         if (_subAxisLenght==-1){
1319                 calculateSubAxisLength();
1320         }
1321         return _subAxisLenght;
1322 }
1323 // ----------------------------------------------------------------------------
1324 double marAxis::getAverageArea(int pIni, int pEnd, kVolume* vol){
1325         marContour      *marcontourHealthy;
1326         int                     ihealthySlice,itemp;
1327         double          acumArea        =       0;
1328         for (ihealthySlice = pIni ; ihealthySlice <= pEnd ; ihealthySlice++ ){
1329                 itemp=ihealthySlice;
1330                 if (itemp<0)                            { itemp=0;                                      }
1331                 if (itemp>=_points.size( ))     { itemp=_points.size( )-1;      }
1332                 marcontourHealthy = getContour( itemp , vol );
1333                 acumArea = acumArea + marcontourHealthy->getArea();
1334         }
1335         return acumArea / (pEnd - pIni + 1);
1336 }
1337 // ----------------------------------------------------------------------------
1338 void marAxis::calculateReferenceArea(kVolume* vol){
1339         _referenceArea = getAverageArea(_healthySliceStart,_healthySliceEnd,  vol);
1340 }
1341 // ----------------------------------------------------------------------------
1342 double marAxis::getReferenceArea(kVolume* vol){
1343         if (_referenceArea==-1){
1344                 calculateReferenceArea(vol);
1345         }
1346         return _referenceArea;
1347 }
1348 // ----------------------------------------------------------------------------
1349 void marAxis::calculateReferenceAverDiam(kVolume* vol){
1350         marContour *marcontourHealthy;
1351         int ihealthySlice;
1352         double acumDiam=0;
1353         for (ihealthySlice=_healthySliceStart; ihealthySlice<=_healthySliceEnd; ihealthySlice++){
1354                 marcontourHealthy = getContour( ihealthySlice , vol );
1355                 acumDiam = acumDiam + marcontourHealthy->getAverageDiameter();
1356         }
1357         _referenceAverDiam = acumDiam / (_healthySliceEnd-_healthySliceStart+1);
1358 }
1359 // ----------------------------------------------------------------------------
1360 double marAxis::getReferenceAverDiam(kVolume* vol){
1361         if (_referenceAverDiam==-1){
1362                 calculateReferenceAverDiam(vol);
1363         }
1364         return _referenceAverDiam;
1365 }
1366
1367 // ----------------------------------------------------------------------------
1368 void marAxis::reset( )
1369 {
1370 //EED Borrame
1371 //    unsigned int i;
1372     kCurve::reset( );
1373
1374     _description                = "";
1375     _healthySlice               = -1;
1376     _startQuant                 = -1;
1377     _finishQuant                = -1;
1378     _actualQuant                = -1;
1379
1380         _subAxisLenght          = -1;
1381         _referenceArea          = -1;
1382         _referenceAverDiam      = -1;
1383
1384         Delete( );
1385
1386 // EED Borrame
1387 /*
1388     for( i = 0; i < _contours.size( ); i++ )
1389       delete _contours[ i ];
1390     _contours.clear( );
1391     for( i = 0; i < _slices.size( ); i++ )
1392         delete _slices[ i ];
1393     _slices.clear( );
1394 */
1395 }
1396
1397 // ----------------------------------------------------------------------------
1398 void marAxis::copyFrom( const marObject& from )
1399 { // TODO
1400 }
1401
1402 // ----------------------------------------------------------------------------
1403 bool marAxis::save( std::ofstream& os )
1404 {
1405     double p[ INDX_count ];
1406     unsigned int i;
1407
1408     i = _description.length( );
1409     os.write( ( const char* )&i, sizeof( int ) );
1410     os.write( ( char* )_description.c_str( ), i * sizeof( char ) );
1411
1412     i = getNumberOfControlPoints( );
1413     os.write( ( const char* )&i, sizeof( int ) );
1414     for( i = 0; i < getNumberOfControlPoints( ); i++ ) {
1415
1416         memcpy( p, _controlPoints[ i ], 3 * sizeof( double ) );
1417         memcpy( p + 3, _controlState[ i ],
1418                 ( INDX_count - 3 ) * sizeof( double ) );
1419         os.write( ( const char* )p, INDX_count * sizeof( double ) );
1420
1421     } // rof
1422     i = getNumberOfContours( );
1423     os.write( ( const char* )&i, sizeof( int ) );
1424     for( i = 0; i < getNumberOfContours( ); i++ )
1425       _contours[ i ]->save( os );
1426
1427     return( true );
1428 }
1429
1430 // ----------------------------------------------------------------------------
1431 bool marAxis::load( std::ifstream& is )
1432 {
1433     double p[ INDX_count ];
1434     int i, n;
1435
1436     reset( );
1437
1438     is.read( ( char* )&n, sizeof( int ) );
1439     _description.resize( n );
1440     is.read( ( char* )_description.c_str( ), n * sizeof( char ) );
1441
1442     is.read( ( char* )&n, sizeof( int ) );
1443     for( i = 0; i < n; i++ ) {
1444
1445         is.read( ( char* )p, INDX_count * sizeof( double ) );
1446         addControlPoint( p, p + 3 );
1447
1448     } // rof
1449     is.read( ( char* )&n, sizeof( int ) );
1450     for( i = 0; i < n; i++ ) {
1451
1452       _contours.push_back( new marContour( 0, getParameters( ) ) );
1453       _contours[ i ]->load( is );
1454
1455     } // rof
1456     return( true );
1457 }
1458 // ----------------------------------------------------------------------------
1459 vtkPolyData* marAxis::Draw( )
1460 {
1461   unsigned int i, j;
1462   double *p;
1463
1464   vtkPoints* allPoints = vtkPoints::New( );
1465   vtkCellArray* allTopology = vtkCellArray::New( );
1466
1467   allTopology->InsertNextCell( _points.size( ) );
1468   for( i = 0, j=0; i < _points.size( ); i++, j++ ) {
1469     p = _points[i];
1470     allPoints->InsertNextPoint( p[ 0 ],  p[ 1 ], p[ 2 ] );
1471     allTopology->InsertCellPoint( j );
1472   } // rof
1473
1474   _allData = vtkPolyData::New( );
1475   _allData->SetPoints( allPoints );
1476   _allData->SetLines( allTopology );
1477   allPoints->Delete();
1478   allTopology->Delete();
1479
1480   return ( _allData );
1481 }
1482 // ----------------------------------------------------------------------------
1483 vtkPolyData *marAxis::GetAxisData( )
1484 {
1485         unsigned int i;
1486         double point[ 3 ];
1487         double radio;
1488         double p[marAxis::INDX_count];
1489
1490
1491         vtkPoints                       *allPoints              = vtkPoints::New( );
1492         vtkCellArray            *allTopology    = vtkCellArray::New( );
1493         vtkDoubleArray          *allRadios              = vtkDoubleArray::New( );
1494         allRadios->SetName("radio");
1495
1496         for( i = 0; i < _controlPoints.size( ); i++ ) {
1497                 getControlPoint(i, p, p+3);
1498                 point[ 0 ]      =  p[marAxis::INDX_X];
1499                 point[ 1 ]      =  p[marAxis::INDX_Y];
1500                 point[ 2 ]      =  p[marAxis::INDX_Z];
1501                 radio           =  p[marAxis::INDX_RAYON]*2.0;
1502                 allPoints    -> InsertPoint( i, point[ 0 ] ,  point[ 1 ] , point[ 2 ]  );  //para saber exactamante el indice que se le asigno
1503
1504                 allRadios       -> InsertValue(i,radio);
1505
1506     if(i>0){
1507        allTopology->InsertNextCell(2);
1508        allTopology->InsertCellPoint(i);
1509        allTopology->InsertCellPoint(i-1);
1510       }
1511
1512         } // rof
1513
1514         _allData = vtkPolyData::New( );
1515         _allData        -> SetPoints( allPoints );
1516         _allData        -> SetLines( allTopology );
1517         _allData        -> GetPointData()->SetScalars(allRadios);
1518
1519
1520         return ( _allData );
1521 }
1522 // ----------------------------------------------------------------------------
1523 void marAxis::Delete( ){
1524         clearAllVectors();
1525
1526         int i;
1527         for (i=0;i<_points.size();i++){
1528                 if (_points[i]!=NULL) {
1529                         delete _points[i];
1530                         _points[i]=NULL;
1531                 }
1532         }
1533         _points.clear();
1534         _signal.clear();
1535
1536
1537         if(_allData) {
1538                 _allData->Delete();
1539                 _allData = NULL;
1540         }
1541
1542 //  reset();
1543 }
1544 // ----------------------------------------------------------------------------
1545 void marAxis::setHealthySlice( int hsS, int hs, int hsE ){
1546           _healthySliceStart    = hsS;
1547           _healthySlice                 = hs;
1548           _healthySliceEnd              = hsE;
1549           _referenceArea=-1;
1550           _referenceAverDiam=-1;
1551   };
1552 // ----------------------------------------------------------------------------
1553 void marAxis::setStartQuant( int sq ){
1554         _startQuant             =       sq;
1555         _subAxisLenght  =       -1;
1556 };
1557 // ----------------------------------------------------------------------------
1558 void marAxis::setFinishQuant( int fq ){
1559         _finishQuant    =       fq;
1560         _subAxisLenght  =       -1;
1561 };
1562
1563 // ----------------------------------------------------------------------------
1564 void marAxis::AddPointToList(double x, double y, double z, int signal)
1565 {
1566         double *p;
1567         p=(double *)malloc(sizeof(double)*3);
1568         p[0]=x;
1569         p[1]=y;
1570         p[2]=z;
1571         _points.push_back( p );
1572         _signal.push_back( signal );
1573 }
1574
1575
1576
1577 // eof - axis.cxx