lib/maracasVisuLib/src/interface/wxWindows/widgets/include/vtk3DSurfaceSTLWidget.cxx

   1 /*=========================================================================
   2
   3   Program:   wxMaracas
   4   Module:    $RCSfile: vtk3DSurfaceSTLWidget.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2009/05/14 13:54:57 $
   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
  18 #include "vtkClosePolyData.h"
  19
  20 #include <vtkTriangleFilter.h>
  21 #include <vtkSTLWriter.h>
  22 #include <vtkInteractorStyleSwitch.h>
  23 #include <vtkPointPicker.h>
  24 #include <vtkCamera.h>
  25 #include <vtkPointData.h>
  26 #include <vtkFloatArray.h>
  27 //For STL
  28 #include <vtkCellLocator.h>
  29 #include <vtkInteractorStyleTrackballCamera.h>
  30 #include <vtkCellArray.h>
  31 #include <vtkStripper.h>
  32 //For rubber stuff:
  33 #include <vtkImplicitSelectionLoop.h>
  34 #include <vtkExtractPolyDataGeometry.h>
  35 #include <vtkPolyDataConnectivityFilter.h>
  36 #include <vtkPolyDataNormals.h>
  37
  38 #include "vtk3DSurfaceSTLWidget.h"
  39
  40 #include "volume.hxx"
  41 #include "marInterface.h"
  42 #include "matrix.h"
  43
  44
  45
  46 /**
  47  * Again wxVTK is an hybrid class, double click was done using wxWindows, simply
  48  * because VTK doesn't provide such facility (as opposed to wheel support that is
  49  * supposed to be merged in VTK trunk sooner or later:
  50  *
  51  * http://public.kitware.com/pipermail/vtkusers/2003-August/019548.html
  52  * http://www.creatis.insa-lyon.fr/~malaterre/vtk/wheel.patch
  53  */
  54
  55 //----------------------------------------------------------------------------
  56 BEGIN_EVENT_TABLE( vtk3DSurfaceSTLWidget, wxVTKRenderWindowInteractor )
  57 //EED    EVT_LEFT_DCLICK( vtk3DSurfaceSTLWidget::OnLeftDClick )
  58 END_EVENT_TABLE( );
  59
  60 //----------------------------------------------------------------------------
  61 vtk3DSurfaceSTLWidget::vtk3DSurfaceSTLWidget(
  62         wxWindow* parent,
  63         wxWindowID id,
  64         const wxPoint& pos,
  65         const wxSize& size,
  66         long style,
  67         const wxString& name)
  68         : wxVTKRenderWindowInteractor( parent, id, pos, size, style, name ){
  69
  70
  71         //make current interactor style be trackbal camera, asked by radiologists.
  72 //      vtkInteractorStyleSwitch *iss = dynamic_cast<vtkInteractorStyleSwitch*>(this->GetInteractorStyle());
  73 //      iss->SetCurrentStyleToTrackballCamera();
  74
  75
  76   vtkInteractorStyle3DMaracas *interactorStyle3DMaracas = vtkInteractorStyle3DMaracas::New();
  77   interactorStyle3DMaracas->SetInteractor (this);
  78   this->SetInteractorStyle(interactorStyle3DMaracas);
  79
  80
  81         _pRenderer  = vtkRenderer::New( );
  82         InitialSphere = 0;
  83     _centralLine       = NULL;
  84     _centralLineMapper = NULL;
  85     _centralLineActor  = NULL;
  86     _axesMapper = NULL;
  87     _axesActor = NULL;
  88     for(int i=0; i<4; i++)
  89     {
  90         _spheres[ i ]           = NULL;
  91         _spheresMapper[ i ] = NULL;
  92         _spheresActor[ i ]      = NULL;
  93     }
  94         _axesActor                      = NULL;
  95         _axesMapper                     = NULL;
  96         _psc                            = NULL;
  97         _stlExternalVessel      = _stlInternalVessel    = NULL;
  98         _surfMapper                     = NULL;
  99         _actorExternalVessel= _actorVessel                      = NULL;
 100         _iasc                           = NULL;
 101 }
 102 //----------------------------------------------------------------------------
 103 vtk3DSurfaceSTLWidget::~vtk3DSurfaceSTLWidget()
 104 {
 105   //good luck:
 106   if( _outActor )           _outActor           ->Delete( );
 107   if( _outMapper )          _outMapper          ->Delete( );
 108   if( _outLine )            _outLine            ->Delete( );
 109   if( _surfActor )          _surfActor          ->Delete( );
 110   if( _surfMapper )         _surfMapper         ->Delete( );
 111   if( _mCubes )             _mCubes                     ->Delete( );
 112   if( _centralLine )        _centralLine        ->Delete( );
 113   if( _centralLineMapper )  _centralLineMapper->Delete( );
 114   if( _centralLineActor )   _centralLineActor ->Delete( );
 115   if( _axesActor )          _axesActor          ->Delete();
 116   if( _axesMapper )         _axesMapper         ->Delete();
 117   if( _actorVessel )        _actorVessel        ->Delete();
 118   if( _actorExternalVessel )_actorExternalVessel->Delete();
 119 //  if( _stlExternalVessel )  _stlExternalVessel->Delete();
 120 //  if( _stlInternalVessel )  _stlInternalVessel->Delete();
 121   if( _psc )                _psc                        ->Delete();
 122   if(_iasc)                 _iasc                       ->Delete();
 123   for(int i=0; i<4; i++)
 124   {
 125     if( _spheres[ i ] )       _spheres[ i ]                     ->Delete( );
 126     if( _spheresMapper[ i ] ) _spheresMapper[ i ]       ->Delete( );
 127     if( _spheresActor[ i ] )  _spheresActor[ i ]        ->Delete( );
 128   }
 129   if( _pRenderer )          _pRenderer          ->Delete( );
 130 }
 131
 132 //----------------------------------------------------------------------------
 133 void vtk3DSurfaceSTLWidget::SetSurfaceColor( float red, float green, float blue )
 134 {
 135   _surfActor->GetProperty()->SetColor(red, green, blue );
 136   _pRenderWindow->Render();
 137 }
 138 //----------------------------------------------------------------------------
 139 void vtk3DSurfaceSTLWidget::SetSurfaceVisibility( bool visible )
 140 {
 141   _surfActor->SetVisibility( visible );
 142   _pRenderWindow->Render();
 143 }
 144 //----------------------------------------------------------------------------
 145 void vtk3DSurfaceSTLWidget::SetSTLSurfaceVisibility( bool intvisible , bool extvisible)
 146 {
 147   ///\todo : STL visibility: both internal and external mold:
 148   _actorVessel->SetVisibility( intvisible );
 149   _actorExternalVessel->SetVisibility( extvisible );
 150   _pRenderWindow->Render();
 151 }
 152 //----------------------------------------------------------------------------
 153 void vtk3DSurfaceSTLWidget::SetSurfaceIsoValue( int isoval )
 154 {
 155   _mCubes->SetValue(0, isoval);
 156   _pRenderWindow->Render();
 157 }
 158 //----------------------------------------------------------------------------
 159 void vtk3DSurfaceSTLWidget::SetSurfaceOpacity( int opaval )
 160 {
 161   //There is no way in Win32 to specify a slider different than 0->100
 162   _surfActor->GetProperty()->SetOpacity( (float)opaval/100 );
 163   _pRenderWindow->Render();
 164 }
 165 //----------------------------------------------------------------------------
 166 void vtk3DSurfaceSTLWidget::GetSphereCenter( double center[3] )
 167 {
 168         _spheres[ 3 ]->GetCenter( center );
 169 }
 170 //----------------------------------------------------------------------------
 171 void vtk3DSurfaceSTLWidget::SetAxis( vtkPolyData *axis )
 172 {
 173         _axesMapper = vtkPolyDataMapper::New( );
 174         _axesMapper->SetInput( axis );
 175         axis->Delete();
 176
 177         _axesActor = vtkActor::New( );
 178         _axesActor->SetMapper( _axesMapper );
 179         _axesActor->GetProperty()->SetColor( 1, 0, 0 );
 180         _pRenderer->AddActor( _axesActor );
 181
 182         _pRenderWindow->Render( );
 183 }
 184 //----------------------------------------------------------------------------
 185 void vtk3DSurfaceSTLWidget::RemoveAxis( )
 186 {
 187   if (_axesActor)
 188   {
 189       _pRenderer->RemoveActor(_axesActor);
 190       _axesActor->Delete();
 191       _axesActor = NULL;
 192       _pRenderWindow->Render( );
 193         }
 194 }
 195 //----------------------------------------------------------------------------
 196 void vtk3DSurfaceSTLWidget::SetCuttingMode( bool mode )
 197 {
 198   if( mode )
 199     {
 200       //start rubber/cutter mode
 201       if(!_iasc)
 202         {
 203         _iasc = vtkInteractorStyleCutter::New();
 204         }
 205       this->SetInteractorStyle( _iasc );
 206     }
 207   else
 208     {
 209    /**
 210     * Copy / paste from vtkPolyDataCutter
 211     */
 212       if( _iasc )
 213       {
 214         _iasc->VisibilityOff(); //turn off the vtkInteractorStyleCutter
 215         vtkImplicitSelectionLoop *loop = vtkImplicitSelectionLoop::New();
 216         loop->SetLoop( _iasc->GetLoopPoints() );
 217         loop->SetNormal( _iasc->GetDirection() );
 218
 219         vtkTriangleFilter *tf = vtkTriangleFilter::New();
 220         tf->SetInput( _stlInternalVessel );
 221
 222         vtkExtractPolyDataGeometry *extract = vtkExtractPolyDataGeometry::New ();
 223         extract->SetInput ( tf->GetOutput());
 224         extract->SetImplicitFunction ( loop );
 225         extract->ExtractInsideOff ();
 226         extract->ExtractBoundaryCellsOff ();
 227         tf->Delete();
 228         loop->Delete();
 229
 230         //#connectivity filter to keep only the largest part
 231         vtkPolyDataConnectivityFilter *connect = vtkPolyDataConnectivityFilter::New();
 232         connect->SetInput( extract->GetOutput() );
 233         connect->SetExtractionModeToLargestRegion ();
 234         extract->Delete();
 235
 236         vtkClosePolyData *close = vtkClosePolyData::New();
 237         close->SetInput( connect->GetOutput() );
 238         connect->Delete();
 239
 240         vtkStripper *strip = vtkStripper::New();
 241         strip->SetInput( close->GetOutput() );
 242         close->Delete();
 243
 244         //If vtkClosePolyData was well written we wouldn't had to recompute the normals
 245         vtkPolyDataNormals *normals = vtkPolyDataNormals::New();
 246         normals->SetInput( strip->GetOutput() );
 247         strip->Delete();
 248
 249         /*vtkPolyDataWriter *polywriter = vtkPolyDataWriter::New();
 250         polywriter->SetInput( normals->GetOutput() );
 251         polywriter->SetFileName( "cutter.vtk" );
 252         polywriter->SetFileTypeToBinary();
 253         polywriter->Write();
 254         polywriter->Delete();*/
 255
 256         normals->Update(); //important before a ShallowCopy Or DeepCopy
 257         _stlInternalVessel->DeepCopy( normals->GetOutput() ); //Use DeepCopy when using ShallowCopy afterwards
 258
 259         tf->SetInput( _stlExternalVessel );
 260         normals->Update();  //important before a ShallowCopy
 261         _stlExternalVessel->ShallowCopy( normals->GetOutput() );
 262         _pRenderWindow->Render( );
 263       }
 264       //make current interactor style be trackbal camera, asked by radiologists.
 265       vtkInteractorStyleTrackballCamera *istc = vtkInteractorStyleTrackballCamera::New();
 266       this->SetInteractorStyle( istc );
 267       istc->Delete();
 268     }
 269
 270 }
 271 //----------------------------------------------------------------------------
 272 void vtk3DSurfaceSTLWidget::ExportSurfaceAsSTL( const char* fileprefix )
 273 {
 274    std::string prefix;
 275
 276    //Now save for real in STL format
 277    //Don't forget to go through a triangle filter
 278    vtkTriangleFilter *triangle = vtkTriangleFilter::New();
 279    triangle->SetInput( _stlInternalVessel );
 280    //In difficult case (such as BESSON) mesh can be more than one piece:
 281    vtkPolyDataConnectivityFilter *pdcf = vtkPolyDataConnectivityFilter::New();
 282    pdcf->SetInput( triangle->GetOutput() );
 283    //In difficult case (such as BESSON) mesh can still be open...close it !
 284    vtkClosePolyData *close = vtkClosePolyData::New();
 285    close->SetInput( pdcf->GetOutput() );
 286
 287    vtkSTLWriter *internal_mold_stl = vtkSTLWriter::New();
 288    internal_mold_stl->SetInput( close->GetOutput() );
 289    prefix = fileprefix;
 290    prefix += "-internal.stl";
 291    internal_mold_stl->SetFileName( prefix.c_str() );
 292    internal_mold_stl->SetFileTypeToBinary();
 293    internal_mold_stl->Write();
 294    internal_mold_stl->Delete();
 295
 296    //convert to triangles first:
 297    triangle->SetInput( _stlExternalVessel );
 298
 299    vtkSTLWriter *external_mold_stl = vtkSTLWriter::New();
 300    external_mold_stl->SetInput( close->GetOutput() );
 301    prefix = fileprefix;
 302    prefix += "-external.stl";
 303    external_mold_stl->SetFileName( prefix.c_str() );
 304    external_mold_stl->SetFileTypeToBinary();
 305    external_mold_stl->Write();
 306    external_mold_stl->Delete();
 307
 308    triangle->Delete();
 309    close->Delete();
 310    pdcf->Delete();
 311 }
 312 //----------------------------------------------------------------------------
 313 void vtk3DSurfaceSTLWidget::ShowMARACASData( marInterface* mar ){
 314     int whd[3];
 315     double minmax[2];
 316
 317     _mar = mar;
 318     _marImageData = _mar->_experiment->getDynData( )->getVolume( )->castVtk();
 319     _marImageData->GetDimensions( whd );
 320     _width  = whd[0];
 321     _height = whd[1];
 322     _depth  = whd[2];
 323
 324     _marImageData->GetScalarRange( minmax );
 325
 326     //Outline v2:
 327     // An outline provides context around the data.
 328     _outLine    = vtkOutlineFilter::New( );
 329     _outLine->SetInput( _marImageData );
 330     _outMapper  = vtkPolyDataMapper::New( );
 331     _outMapper->SetInput( _outLine->GetOutput( ) );
 332     _outActor   = vtkActor::New( );
 333     _outActor->SetMapper( _outMapper );
 334     _outActor->GetProperty( )->SetColor( 0.7, 0.0, 0.9 );
 335     //_outActor->PickableOff( );
 336     _pRenderer->AddActor( _outActor );
 337
 338     // Surface
 339     _mCubes = vtkMarchingCubes::New( );
 340     _surfMapper = vtkPolyDataMapper::New( );
 341     _surfActor = vtkActor::New( );
 342
 343     _mCubes->SetInput( _marImageData );
 344     _mCubes->SetValue( 0, minmax[1] / 4 );
 345
 346     vtkStripper *stripper = vtkStripper::New();
 347     stripper->SetInput( _mCubes->GetOutput( ) );
 348
 349         _surfMapper->SetInput( stripper->GetOutput() );
 350     _surfMapper->ScalarVisibilityOff( );
 351     stripper->Delete();
 352
 353     _surfActor->SetMapper( _surfMapper );
 354     _surfActor->PickableOff( );
 355     _surfActor->GetProperty( )->SetColor( 0.9803, 0.9215, 0.8392 );
 356     _surfActor->GetProperty( )->SetOpacity( 0.5 );
 357     _pRenderer->AddActor( _surfActor );
 358
 359     // 1. ParallelProjectionOn should be set after AddActor (otherwise call vtkRenderer::ResetCameraClippingRange
 360     // 2. ParallelProjectionOn is *necessary* for the vtkImplicitSelectionLoop
 361     // otherwise this give a cone instead of a cylinder cutting.
 362     _pRenderer->GetActiveCamera()->ParallelProjectionOn();
 363
 364     this->ConfigureVTK( );
 365 }
 366
 367 //----------------------------------------------------------------------------
 368 void vtk3DSurfaceSTLWidget::ShowMARACASDataAndAxe( marInterface* mar ){
 369   marAxis* temp;
 370   vtkPolyData* allData;
 371
 372
 373   this->ShowMARACASData( mar );
 374
 375   _pRenderer->SetBackground( 0.75, 0.75, 0.75 );
 376
 377   // Axis
 378   _mar->_experiment->setAxis( 0 );
 379   temp = _mar->_experiment->getAxis( );  // ??? getActualAxis ??
 380   allData = temp->Draw( );
 381   this->SetAxis( allData );
 382 }
 383
 384 //----------------------------------------------------------------------------
 385 void vtk3DSurfaceSTLWidget::SetInitialPoint( double* pickPoint, double* cameraPos )
 386 {
 387     gtm::TVector< double > pO( 3 ), pF( 3 ), pp( 3 ), cp( 3 );
 388     gtm::TVector< double > xc( 3 );
 389     gtm::TVector< double > x1( 3 ), n1( 3 );
 390     gtm::TVector< double > x2( 3 ), n2( 3 );
 391     gtm::TVector< double > tmp( 3 );
 392     double fac;
 393     bool success = true;
 394
 395     if( _centralLineActor )
 396     {
 397         _pRenderer->RemoveActor( _centralLineActor );
 398         _centralLineActor->Delete( );
 399     } // fi
 400     if( _centralLineMapper ) _centralLineMapper->Delete( );
 401     if( _centralLine )       _centralLine->Delete( );
 402
 403     _centralLine       = NULL;
 404     _centralLineMapper = NULL;
 405     _centralLineActor  = NULL;
 406
 407     for(int i=0; i<4; i++)
 408     {
 409       if( _spheresActor[ i ] )
 410       {
 411         _pRenderer->RemoveActor( _spheresActor[ i ] );
 412         _spheresActor[ i ]->Delete( );
 413       } // fi
 414       if( _spheresMapper[ i ] ) _spheresMapper[ i ]->Delete( );
 415       if( _spheres[ i ] )       _spheres[ i ]->Delete( );
 416       _spheres[ i ] = NULL;
 417       _spheresMapper[ i ] = NULL;
 418       _spheresActor[ i ] = NULL;
 419     } //rof
 420
 421     fac = GTM_MAX( _width, _height );
 422     fac = 2 * GTM_MAX( fac, _depth );
 423     pp( 0 ) = pickPoint[ 0 ]; pp( 1 ) = pickPoint[ 1 ]; pp( 2 ) = pickPoint[ 2 ];
 424     cp( 0 ) = cameraPos[ 0 ]; cp( 1 ) = cameraPos[ 1 ]; cp( 2 ) = cameraPos[ 2 ];
 425
 426     if( this->FindCubePointsFromPoints(
 427       pO.GetAnsiRef( ), pF.GetAnsiRef( ),
 428       pp.GetAnsiRef( ), cp.GetAnsiRef( )
 429     ) ) {
 430
 431         if( this->GetPointAndNormalIntersection(
 432             x1.GetAnsiRef( ), n1.GetAnsiRef( ),
 433             pO.GetAnsiRef( ), pF.GetAnsiRef( )
 434     ) ) {
 435
 436             if( this->GetPointAndNormalIntersection(
 437           x2.GetAnsiRef( ), n2.GetAnsiRef( ),
 438           ( x1 - n1 ).GetAnsiRef( ), ( x1 - ( n1 * fac ) ).GetAnsiRef( )
 439         ) ) {
 440
 441                 xc = ( x2 + x1 ) * 0.5;
 442                 _mar->_experiment->setStartPoint( (int)xc( 0 ), (int)xc( 1 ), (int)xc( 2 ) );
 443
 444       for(int i=0; i<4; i++)
 445       {
 446         _spheres[ i ] = vtkSphereSource::New( );
 447         _spheresMapper[ i ] = vtkPolyDataMapper::New( );
 448         _spheresMapper[ i ]->SetInput( _spheres[ i ]->GetOutput( ) );
 449         _spheresActor[ i ] = vtkActor::New( );
 450         _spheresActor[ i ]->SetMapper( _spheresMapper[ i ] );
 451         _spheresActor[ i ]->PickableOff( );
 452         _pRenderer->AddActor( _spheresActor[ i ] );
 453       }
 454
 455                 _spheres[ 0 ]->SetCenter( x1( 0 ), x1( 1 ), x1( 2 ) );
 456                 _spheres[ 1 ]->SetCenter( x2( 0 ), x2( 1 ), x2( 2 ) );
 457                 _spheres[ 2 ]->SetCenter( xc( 0 ), xc( 1 ), xc( 2 ) );
 458                 _spheres[ 3 ]->SetCenter( xc( 0 ), xc( 1 ), xc( 2 ) );
 459
 460                 _spheres[ 0 ]->SetRadius( 0.5 );
 461                 _spheres[ 1 ]->SetRadius( 0.5 );
 462                 _spheres[ 2 ]->SetRadius( 0.5 );
 463                 _spheres[ 3 ]->SetRadius( ( xc - x1 ).GetNorm( ) );
 464
 465                 _spheresActor[ 0 ]->GetProperty( )->SetColor( 1.0, 0.0, 0.0 );
 466                 _spheresActor[ 1 ]->GetProperty( )->SetColor( 0.0, 1.0, 0.0 );
 467                 _spheresActor[ 2 ]->GetProperty( )->SetColor( 0.0, 0.0, 1.0 );
 468                 _spheresActor[ 3 ]->GetProperty( )->SetColor( 1.0, 0.0, 0.0 );
 469                 _spheresActor[ 3 ]->GetProperty( )->SetOpacity( 0.3 );
 470
 471                 vtkPoints* points = vtkPoints::New( );
 472                 points->InsertNextPoint( x1.GetAnsiRef( ) );
 473                 points->InsertNextPoint( x2.GetAnsiRef( ) );
 474
 475                 vtkCellArray* array = vtkCellArray::New( );
 476                 array->InsertNextCell( 2 );
 477                 array->InsertCellPoint( 0 );
 478                 array->InsertCellPoint( 1 );
 479
 480                 _centralLine = vtkPolyData::New( );
 481                 _centralLine->SetPoints( points );
 482                 _centralLine->SetLines( array );
 483       points->Delete();
 484       array->Delete();
 485
 486                 _centralLineMapper = vtkPolyDataMapper::New( );
 487                 _centralLineMapper->SetInput( _centralLine );
 488
 489                 _centralLineActor = vtkActor::New( );
 490                 _centralLineActor->SetMapper( _centralLineMapper );
 491                 _centralLineActor->GetProperty( )->SetColor( 1.0, 1.0, 1.0 );
 492                 _centralLineActor->PickableOff( );
 493                 _pRenderer->AddActor( _centralLineActor );
 494
 495             } // fi
 496
 497         } else success = false;
 498
 499     } else success = false;
 500
 501     // Show a message, if any.
 502     if( !success )
 503         wxMessageBox(
 504         _T("Error: Initial point can't be set.\nPlease choose another point."),
 505         _T("Error"), wxOK | wxCENTRE | wxICON_HAND, this
 506         );
 507
 508     // Finish
 509     _pRenderWindow->Render( );
 510     InitialSphere = success;
 511 }
 512
 513 //----------------------------------------------------------------------------
 514 void vtk3DSurfaceSTLWidget::SetInitialPoint( double* point )
 515 {
 516 }
 517 //----------------------------------------------------------------------------
 518 void vtk3DSurfaceSTLWidget::ConfigureVTK( )
 519 {
 520   // get render window
 521   _pRenderWindow =  this->GetRenderWindow(  );
 522
 523   // connect renderer and render window and configure render window
 524   _pRenderWindow->AddRenderer( _pRenderer );
 525
 526   // configure renderer
 527   _pRenderer->SetBackground( 0.0, 0.0, 0.0 );
 528   //_pRenderer->SetBackground( 1, 1, 0);  //FIXME
 529   _pRenderer->GetActiveCamera( )->Zoom( 1.0 );
 530   _pRenderer->GetActiveCamera( )->SetClippingRange( 1, 1000 );
 531 }
 532
 533
 534 //----------------------------------------------------------------------------
 535 bool vtk3DSurfaceSTLWidget::IntersectPlaneWithLine( double* p, double* x1, double* x2, double* x3, double* x4, double* x5 )
 536 {
 537     gtm::TVector< double > vx1( 3 ), vx2( 3 ), vx3( 3 );
 538     gtm::TVector< double > vx4( 3 ), vx5( 3 ), vx6( p, 3, false );
 539     gtm::TMatrix< double > mU( 4, 4 ), mD( 4, 4 );
 540     double t;
 541     bool ret;
 542
 543     vx1 = x1;
 544     vx2 = x2;
 545     vx3 = x3;
 546     vx4 = x4;
 547     vx5 = x5;
 548     vx6 = vx5 - vx4;
 549
 550     mD( 0, 0 ) = mU( 0, 0 ) = 1;
 551     mD( 0, 1 ) = mU( 0, 1 ) = vx1( 0 );
 552     mD( 0, 2 ) = mU( 0, 2 ) = vx1( 1 );
 553     mD( 0, 3 ) = mU( 0, 3 ) = vx1( 2 );
 554     mD( 1, 0 ) = mU( 1, 0 ) = 1;
 555     mD( 1, 1 ) = mU( 1, 1 ) = vx2( 0 );
 556     mD( 1, 2 ) = mU( 1, 2 ) = vx2( 1 );
 557     mD( 1, 3 ) = mU( 1, 3 ) = vx2( 2 );
 558     mD( 2, 0 ) = mU( 2, 0 ) = 1;
 559     mD( 2, 1 ) = mU( 2, 1 ) = vx3( 0 );
 560     mD( 2, 2 ) = mU( 2, 2 ) = vx3( 1 );
 561     mD( 2, 3 ) = mU( 2, 3 ) = vx3( 2 );
 562     mU( 3, 0 ) = 1;
 563     mU( 3, 1 ) = vx4( 0 );
 564     mU( 3, 2 ) = vx4( 1 );
 565     mU( 3, 3 ) = vx4( 2 );
 566     mD( 3, 0 ) = 0;
 567     mD( 3, 1 ) = vx6( 0 );
 568     mD( 3, 2 ) = vx6( 1 );
 569     mD( 3, 3 ) = vx6( 2 );
 570
 571     ret = ( mD.Det( ) != 0 );
 572     if( ret ) {
 573
 574         t = mU.Det( ) / mD.Det( );
 575         vx6 = ( ( vx4 - vx5 ) * t ) + vx4;
 576
 577     } // fi
 578     return( ret );
 579
 580 }
 581
 582 //----------------------------------------------------------------------------
 583 bool vtk3DSurfaceSTLWidget::FindCubePointsFromPoints( double* pO, double* pF, double* pickPoint, double* cameraPos )
 584 {
 585     gtm::TVector< double > p1( 3 ), p2( 3 ), p3( 3 ), p4( 3 ), p5( 3 ), p6( 3 );
 586     gtm::TVector< double > c( 3 );
 587     gtm::TVector< double >* swp;
 588     gtm::TVector< double > tpO( pO, 3, false ), tpF( pF, 3, false );
 589     std::vector< gtm::TVector< double >* > points;
 590     double x1[ 3 ], x2[ 3 ], x3[ 3 ], d1, d2;
 591     int i, j;
 592
 593     // 1st plane intersection
 594     x1[ 0 ] = 0; x1[ 1 ] = _height; x1[ 2 ] =      0;
 595     x2[ 0 ] = 0; x2[ 1 ] =       0; x2[ 2 ] =      0;
 596     x3[ 0 ] = 0; x3[ 1 ] =       0; x3[ 2 ] = _depth;
 597     if( this->IntersectPlaneWithLine( p1.GetAnsiRef( ), x1, x2, x3, cameraPos, pickPoint ) )
 598         points.push_back( &p1 );
 599
 600     // 2nd plane intersection
 601     x1[ 0 ] =      0; x1[ 1 ] = _height; x1[ 2 ] = _depth;
 602     x2[ 0 ] =      0; x2[ 1 ] =       0; x2[ 2 ] = _depth;
 603     x3[ 0 ] = _width; x3[ 1 ] =       0; x3[ 2 ] = _depth;
 604     if( this->IntersectPlaneWithLine( p2.GetAnsiRef( ), x1, x2, x3, cameraPos, pickPoint ) )
 605         points.push_back( &p2 );
 606
 607     // 3rd plane intersection
 608     x1[ 0 ] =      0; x1[ 1 ] = 0; x1[ 2 ] = _depth;
 609     x2[ 0 ] =      0; x2[ 1 ] = 0; x2[ 2 ] =      0;
 610     x3[ 0 ] = _width; x3[ 1 ] = 0; x3[ 2 ] =      0;
 611     if( this->IntersectPlaneWithLine( p3.GetAnsiRef( ), x1, x2, x3, cameraPos, pickPoint ) )
 612         points.push_back( &p3 );
 613
 614     // 4th plane intersection
 615     x1[ 0 ] = _width; x1[ 1 ] = _height; x1[ 2 ] = _depth;
 616     x2[ 0 ] = _width; x2[ 1 ] =       0; x2[ 2 ] = _depth;
 617     x3[ 0 ] = _width; x3[ 1 ] =       0; x3[ 2 ] =      0;
 618     if( this->IntersectPlaneWithLine( p4.GetAnsiRef( ), x1, x2, x3, cameraPos, pickPoint ) )
 619         points.push_back( &p4 );
 620
 621     // 5th plane intersection
 622     x1[ 0 ] = _width; x1[ 1 ] =       0; x1[ 2 ] = 0;
 623     x2[ 0 ] =      0; x2[ 1 ] =       0; x2[ 2 ] = 0;
 624     x3[ 0 ] =      0; x3[ 1 ] = _height; x3[ 2 ] = 0;
 625     if( this->IntersectPlaneWithLine( p5.GetAnsiRef( ), x1, x2, x3, cameraPos, pickPoint ) )
 626         points.push_back( &p5 );
 627
 628     // 6th plane intersection
 629     x1[ 0 ] =      0; x1[ 1 ] = _height; x1[ 2 ] =      0;
 630     x2[ 0 ] =      0; x2[ 1 ] = _height; x2[ 2 ] = _depth;
 631     x3[ 0 ] = _width; x3[ 1 ] = _height; x3[ 2 ] = _depth;
 632     if( this->IntersectPlaneWithLine( p6.GetAnsiRef( ), x1, x2, x3, cameraPos, pickPoint ) )
 633         points.push_back( &p6 );
 634
 635     if( points.size( ) >= 2 ) { // Did I find at least 2 points?
 636
 637         c( 0 ) = ( double )_width  / 2.0;
 638         c( 1 ) = ( double )_height / 2.0;
 639         c( 2 ) = ( double )_depth  / 2.0;
 640
 641         // Sort with bubble sort. Only 30 iterations!
 642         for( i = 0; i < points.size( ); i++ ) {
 643             for( j = 0; j < points.size( ) - 1; j++ ) {
 644
 645                 d1 = ( c - *points[ j ] ).GetNorm( );
 646                 d2 = ( c - *points[ j + 1 ] ).GetNorm( );
 647                 if( d2 < d1 ) {
 648
 649                     swp = points[ j ];
 650                     points[ j ] = points[ j + 1 ];
 651                     points[ j + 1 ] = swp;
 652
 653                 } // fi
 654
 655     } // rof
 656
 657         } // rof
 658
 659         // Order the two points according to distance to camera.
 660         c = cameraPos;
 661         d1 = ( c - *points[ 0 ] ).GetNorm( );
 662         d2 = ( c - *points[ 1 ] ).GetNorm( );
 663         tpO = ( d1 < d2 )? *points[ 0 ]: *points[ 1 ];
 664         tpF = ( d1 > d2 )? *points[ 0 ]: *points[ 1 ];
 665         return( true );
 666
 667     } else return( false );
 668
 669 }
 670
 671 //----------------------------------------------------------------------------
 672 bool vtk3DSurfaceSTLWidget::GetPointAndNormalIntersection( double* p, double* n, double* pO, double* pF )
 673 {
 674     gtm::TVector< double > n1( 3 ), n2( 3 ), n3( 3 );
 675     int subId, cellId, returnVal;
 676     double t, pcoords[ 3 ], x[ 3 ];
 677     double fpO[ 3 ], fpF[ 3 ];
 678     double p1[ 3 ], p2[ 3 ], p3[ 3 ];
 679     vtkPolyData* data = _mCubes->GetOutput( );
 680     vtkCellLocator* locator = vtkCellLocator::New( );
 681
 682     locator->SetDataSet( data );
 683     locator->Initialize( );
 684     locator->Update( );
 685
 686     fpO[ 0 ] = pO[ 0 ]; fpO[ 1 ] = pO[ 1 ]; fpO[ 2 ] = pO[ 2 ];
 687     fpF[ 0 ] = pF[ 0 ]; fpF[ 1 ] = pF[ 1 ]; fpF[ 2 ] = pF[ 2 ];
 688     returnVal = locator->IntersectWithLine( fpO, fpF, 0.1, t, x, pcoords, subId, cellId );
 689     locator->Delete( );
 690
 691     if( returnVal )
 692     {
 693         vtkCell* cell = data->GetCell( cellId );
 694         vtkPoints* points = cell->GetPoints( );
 695
 696         data->GetPointData( )->GetNormals( )->GetTuple( cell->GetPointIds( )->GetId( 0 ), n1.GetAnsiRef( ) );
 697         data->GetPointData( )->GetNormals( )->GetTuple( cell->GetPointIds( )->GetId( 1 ), n2.GetAnsiRef( ) );
 698         data->GetPointData( )->GetNormals( )->GetTuple( cell->GetPointIds( )->GetId( 2 ), n3.GetAnsiRef( ) );
 699
 700         n1 += n2 + n3;
 701         n1 *= ( 1.0 / 3.0 );
 702         n1.Normalize( );
 703         n[ 0 ] = n1( 0 ); n[ 1 ] = n1( 1 ); n[ 2 ] = n1( 2 );
 704
 705         points->GetPoint( 0, p1 );
 706         points->GetPoint( 1, p2 );
 707         points->GetPoint( 2, p3 );
 708         this->IntersectPlaneWithLine( p, p1, p2, p3, pO, pF );
 709         return( true );
 710     } else return( false );
 711
 712 }
 713 //----------------------------------------------------------------------------
 714  /**
 715   * What we are doing here is from a LFV marAxis to turn it into a vtkPolyData
 716   * Therefore it should be a lot more easier to integrate in a real vtk class
 717   * for filtering : vtkMaracasImageToSTL
 718   */
 719 vtkPolyData* vtk3DSurfaceSTLWidget::ConvertMarAxisToPolyData()
 720 {
 721   double p1[marAxis::INDX_count];
 722   //retrieve the actual axis
 723   marAxis* ax = _mar->_experiment->getAxis();
 724   //Number of axis points
 725   int numPoints = ax->getNumberOfControlPoints();
 726
 727   vtkPolyData *axis = vtkPolyData::New();
 728   vtkPoints *points = vtkPoints::New();
 729   vtkFloatArray *scalars = vtkFloatArray::New();
 730   scalars->SetNumberOfComponents (2);
 731   //Field are : 1. radius (=INDX_RAYON) , 2. average (=INDX_SIGNALVALUE)
 732
 733   double tmp1,tmp2;
 734   for(int i=0; i<numPoints; i++){
 735     ax->getControlPoint(i, p1, p1+3);
 736     points->InsertPoint(i, p1);
 737         tmp1 = p1[marAxis::INDX_RAYON];
 738 //      tmp2 = p1[marAxis::INDX_SIGNALVALUE];
 739         tmp2 = (double)ax->getSignal(i);
 740     scalars->InsertTuple2(i,tmp1, tmp2);
 741
 742   }
 743
 744   // We now assign the pieces to the vtkPolyData.
 745   axis->SetPoints(points);
 746   points->Delete();
 747   axis->GetPointData()->SetScalars(scalars);
 748   scalars->Delete();
 749
 750   return axis;
 751 }
 752 //----------------------------------------------------------------------------
 753 void vtk3DSurfaceSTLWidget::ConstructVessel( )
 754 {
 755   vtkPolyData *axis = this->ConvertMarAxisToPolyData( );
 756
 757   //retrieve the vtk image data
 758
 759   if( !_psc) this->_psc = vtkImagePolyDataSeedConnectivity::New();
 760   this->_psc->SetInput( _marImageData );
 761   this->_psc->SetAxis( axis );
 762   this->_psc->Update();
 763   axis->Delete();
 764
 765 /*  vtkStripper *strip = vtkStripper::New();
 766   strip->SetInput( this->_psc->GetOutput() );  //GetOutput() -> inner mold
 767   strip->Update(); //Important before a ShallowCopy !!!
 768
 769   if( !_stlInternalVessel) _stlInternalVessel = vtkPolyData::New();
 770   _stlInternalVessel->ShallowCopy( strip->GetOutput() );
 771   strip->Delete();*/
 772   _stlInternalVessel = this->_psc->GetOutput();
 773
 774 /*  if( !_stlExternalVessel) _stlExternalVessel = vtkPolyData::New();
 775   _stlExternalVessel->ShallowCopy( this->_psc->GetOuterMold() );
 776   //this->_psc->GetOuterMold()->Delete();*/
 777   _stlExternalVessel = this->_psc->GetOuterMold();
 778
 779   vtkPolyDataMapper *dsm1 = vtkPolyDataMapper ::New();
 780   dsm1->SetInput ( _stlInternalVessel ); //fasten stuff
 781   dsm1->ScalarVisibilityOff();
 782
 783   _actorVessel = vtkActor::New();
 784   _actorVessel->SetMapper (dsm1);
 785
 786   //improve visibility:
 787   _actorVessel->GetProperty()->SetColor (1,0,0);
 788
 789   vtkPolyDataMapper *dsm2 = vtkPolyDataMapper ::New();
 790   dsm2->SetInput ( _stlExternalVessel );
 791   dsm2->ScalarVisibilityOff();
 792
 793   _actorExternalVessel = vtkActor::New();
 794   _actorExternalVessel->SetMapper (dsm2);
 795
 796   //improve visibility:
 797   _actorExternalVessel->GetProperty()->SetRepresentationToWireframe();
 798
 799   _pRenderer->AddActor( _actorVessel );
 800   _pRenderer->AddActor( _actorExternalVessel );
 801   _pRenderWindow->Render( );
 802
 803   dsm1->Delete();
 804   dsm2->Delete();
 805 }
 806 //----------------------------------------------------------------------------
 807 void vtk3DSurfaceSTLWidget::SetSTLThresholdRatio(double ratio)
 808 {
 809   this->_psc->SetThresholdRatio( (double)(ratio/100) );
 810   _pRenderWindow->Render( );
 811 }
 812 //----------------------------------------------------------------------------
 813 float vtk3DSurfaceSTLWidget::GetSTLThreshold()
 814 {
 815   return  100*this->_psc->GetThresholdRatio( );
 816 }
 817 //----------------------------------------------------------------------------