lib/maracasVisuLib/src/kernel/marAxis.cpp

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