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

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