vtk3DSurfaceSTLWidget.cxx

Go to the documentation of this file.

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

  Program:   wxMaracas
  Module:    $RCSfile: vtk3DSurfaceSTLWidget.cxx,v $
  Language:  C++
  Date:      $Date: 2009/05/14 13:54:57 $
  Version:   $Revision: 1.1 $

  Copyright: (c) 2002, 2003
  License:
  
     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notice for more information.

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

#include "vtkClosePolyData.h"

#include <vtkTriangleFilter.h>
#include <vtkSTLWriter.h>
#include <vtkInteractorStyleSwitch.h>
#include <vtkPointPicker.h>
#include <vtkCamera.h>
#include <vtkPointData.h>
#include <vtkFloatArray.h>
//For STL
#include <vtkCellLocator.h>
#include <vtkInteractorStyleTrackballCamera.h>
#include <vtkCellArray.h>
#include <vtkStripper.h>
//For rubber stuff:
#include <vtkImplicitSelectionLoop.h>
#include <vtkExtractPolyDataGeometry.h>
#include <vtkPolyDataConnectivityFilter.h>
#include <vtkPolyDataNormals.h>

#include "vtk3DSurfaceSTLWidget.h"

#include "volume.hxx"
#include "marInterface.h"
#include "matrix.h"



//----------------------------------------------------------------------------
BEGIN_EVENT_TABLE( vtk3DSurfaceSTLWidget, wxVTKRenderWindowInteractor )
//EED    EVT_LEFT_DCLICK( vtk3DSurfaceSTLWidget::OnLeftDClick )
END_EVENT_TABLE( );

//----------------------------------------------------------------------------
vtk3DSurfaceSTLWidget::vtk3DSurfaceSTLWidget(
        wxWindow* parent,
        wxWindowID id,
        const wxPoint& pos,
        const wxSize& size,
        long style,
        const wxString& name)
        : wxVTKRenderWindowInteractor( parent, id, pos, size, style, name ){


        //make current interactor style be trackbal camera, asked by radiologists.
//      vtkInteractorStyleSwitch *iss = dynamic_cast<vtkInteractorStyleSwitch*>(this->GetInteractorStyle());
//      iss->SetCurrentStyleToTrackballCamera();
  

  vtkInteractorStyle3DMaracas *interactorStyle3DMaracas = vtkInteractorStyle3DMaracas::New(); 
  interactorStyle3DMaracas->SetInteractor (this);
  this->SetInteractorStyle(interactorStyle3DMaracas);


        _pRenderer  = vtkRenderer::New( );
        InitialSphere = 0;
    _centralLine       = NULL;
    _centralLineMapper = NULL;
    _centralLineActor  = NULL;
    _axesMapper = NULL;
    _axesActor = NULL;
    for(int i=0; i<4; i++)
    {
        _spheres[ i ]           = NULL;
        _spheresMapper[ i ] = NULL;
        _spheresActor[ i ]      = NULL;
    }
        _axesActor                      = NULL;
        _axesMapper                     = NULL;
        _psc                            = NULL;
        _stlExternalVessel      = _stlInternalVessel    = NULL;
        _surfMapper                     = NULL;
        _actorExternalVessel= _actorVessel                      = NULL;
        _iasc                           = NULL;
}
//----------------------------------------------------------------------------
vtk3DSurfaceSTLWidget::~vtk3DSurfaceSTLWidget()
{
  //good luck:
  if( _outActor )           _outActor           ->Delete( );
  if( _outMapper )          _outMapper          ->Delete( );
  if( _outLine )            _outLine            ->Delete( );
  if( _surfActor )          _surfActor          ->Delete( );
  if( _surfMapper )         _surfMapper         ->Delete( );
  if( _mCubes )             _mCubes                     ->Delete( );
  if( _centralLine )        _centralLine        ->Delete( );
  if( _centralLineMapper )  _centralLineMapper->Delete( );
  if( _centralLineActor )   _centralLineActor ->Delete( );
  if( _axesActor )          _axesActor          ->Delete();
  if( _axesMapper )         _axesMapper         ->Delete();
  if( _actorVessel )        _actorVessel        ->Delete();
  if( _actorExternalVessel )_actorExternalVessel->Delete();
//  if( _stlExternalVessel )  _stlExternalVessel->Delete();
//  if( _stlInternalVessel )  _stlInternalVessel->Delete();
  if( _psc )                _psc                        ->Delete();
  if(_iasc)                 _iasc                       ->Delete();
  for(int i=0; i<4; i++)
  {
    if( _spheres[ i ] )       _spheres[ i ]                     ->Delete( );
    if( _spheresMapper[ i ] ) _spheresMapper[ i ]       ->Delete( );
    if( _spheresActor[ i ] )  _spheresActor[ i ]        ->Delete( );
  }
  if( _pRenderer )          _pRenderer          ->Delete( );
}

//----------------------------------------------------------------------------
void vtk3DSurfaceSTLWidget::SetSurfaceColor( float red, float green, float blue )
{
  _surfActor->GetProperty()->SetColor(red, green, blue );
  _pRenderWindow->Render();
}
//----------------------------------------------------------------------------
void vtk3DSurfaceSTLWidget::SetSurfaceVisibility( bool visible )
{
  _surfActor->SetVisibility( visible );
  _pRenderWindow->Render();
}
//----------------------------------------------------------------------------
void vtk3DSurfaceSTLWidget::SetSTLSurfaceVisibility( bool intvisible , bool extvisible)
{
  _actorVessel->SetVisibility( intvisible );
  _actorExternalVessel->SetVisibility( extvisible );
  _pRenderWindow->Render();
}
//----------------------------------------------------------------------------
void vtk3DSurfaceSTLWidget::SetSurfaceIsoValue( int isoval )
{
  _mCubes->SetValue(0, isoval);
  _pRenderWindow->Render();
}
//----------------------------------------------------------------------------
void vtk3DSurfaceSTLWidget::SetSurfaceOpacity( int opaval )
{
  //There is no way in Win32 to specify a slider different than 0->100
  _surfActor->GetProperty()->SetOpacity( (float)opaval/100 );
  _pRenderWindow->Render();
}
//----------------------------------------------------------------------------
void vtk3DSurfaceSTLWidget::GetSphereCenter( double center[3] )
{
        _spheres[ 3 ]->GetCenter( center );
}
//----------------------------------------------------------------------------
void vtk3DSurfaceSTLWidget::SetAxis( vtkPolyData *axis )
{
        _axesMapper = vtkPolyDataMapper::New( );
        _axesMapper->SetInput( axis );
        axis->Delete();

        _axesActor = vtkActor::New( );
        _axesActor->SetMapper( _axesMapper );
        _axesActor->GetProperty()->SetColor( 1, 0, 0 );
        _pRenderer->AddActor( _axesActor );

        _pRenderWindow->Render( );
}
//----------------------------------------------------------------------------
void vtk3DSurfaceSTLWidget::RemoveAxis( )
{
  if (_axesActor)
  {
      _pRenderer->RemoveActor(_axesActor);
      _axesActor->Delete();
      _axesActor = NULL;
      _pRenderWindow->Render( );
        }
}
//----------------------------------------------------------------------------
void vtk3DSurfaceSTLWidget::SetCuttingMode( bool mode )
{
  if( mode )
    {
      //start rubber/cutter mode
      if(!_iasc) 
        {
        _iasc = vtkInteractorStyleCutter::New();
        }
      this->SetInteractorStyle( _iasc );
    }
  else
    {
      if( _iasc )
      {
        _iasc->VisibilityOff(); //turn off the vtkInteractorStyleCutter
        vtkImplicitSelectionLoop *loop = vtkImplicitSelectionLoop::New();
        loop->SetLoop( _iasc->GetLoopPoints() );
        loop->SetNormal( _iasc->GetDirection() );

        vtkTriangleFilter *tf = vtkTriangleFilter::New();
        tf->SetInput( _stlInternalVessel );

        vtkExtractPolyDataGeometry *extract = vtkExtractPolyDataGeometry::New ();
        extract->SetInput ( tf->GetOutput());
        extract->SetImplicitFunction ( loop );
        extract->ExtractInsideOff ();
        extract->ExtractBoundaryCellsOff ();
        tf->Delete();
        loop->Delete();

        //#connectivity filter to keep only the largest part
        vtkPolyDataConnectivityFilter *connect = vtkPolyDataConnectivityFilter::New();
        connect->SetInput( extract->GetOutput() );
        connect->SetExtractionModeToLargestRegion ();
        extract->Delete();

        vtkClosePolyData *close = vtkClosePolyData::New();
        close->SetInput( connect->GetOutput() );
        connect->Delete();

        vtkStripper *strip = vtkStripper::New();
        strip->SetInput( close->GetOutput() );
        close->Delete();

        //If vtkClosePolyData was well written we wouldn't had to recompute the normals
        vtkPolyDataNormals *normals = vtkPolyDataNormals::New();
        normals->SetInput( strip->GetOutput() );
        strip->Delete();

        /*vtkPolyDataWriter *polywriter = vtkPolyDataWriter::New();
        polywriter->SetInput( normals->GetOutput() );
        polywriter->SetFileName( "cutter.vtk" );
        polywriter->SetFileTypeToBinary();
        polywriter->Write();
        polywriter->Delete();*/

        normals->Update(); //important before a ShallowCopy Or DeepCopy
        _stlInternalVessel->DeepCopy( normals->GetOutput() ); //Use DeepCopy when using ShallowCopy afterwards

        tf->SetInput( _stlExternalVessel );
        normals->Update();  //important before a ShallowCopy
        _stlExternalVessel->ShallowCopy( normals->GetOutput() );
        _pRenderWindow->Render( );
      }
      //make current interactor style be trackbal camera, asked by radiologists.
      vtkInteractorStyleTrackballCamera *istc = vtkInteractorStyleTrackballCamera::New();
      this->SetInteractorStyle( istc );
      istc->Delete();
    }

}
//----------------------------------------------------------------------------
void vtk3DSurfaceSTLWidget::ExportSurfaceAsSTL( const char* fileprefix )
{
   std::string prefix;

   //Now save for real in STL format
   //Don't forget to go through a triangle filter
   vtkTriangleFilter *triangle = vtkTriangleFilter::New();
   triangle->SetInput( _stlInternalVessel );
   //In difficult case (such as BESSON) mesh can be more than one piece:
   vtkPolyDataConnectivityFilter *pdcf = vtkPolyDataConnectivityFilter::New();
   pdcf->SetInput( triangle->GetOutput() );
   //In difficult case (such as BESSON) mesh can still be open...close it !
   vtkClosePolyData *close = vtkClosePolyData::New();
   close->SetInput( pdcf->GetOutput() );
   
   vtkSTLWriter *internal_mold_stl = vtkSTLWriter::New();
   internal_mold_stl->SetInput( close->GetOutput() );
   prefix = fileprefix;
   prefix += "-internal.stl";
   internal_mold_stl->SetFileName( prefix.c_str() );
   internal_mold_stl->SetFileTypeToBinary();
   internal_mold_stl->Write();
   internal_mold_stl->Delete();
   
   //convert to triangles first:
   triangle->SetInput( _stlExternalVessel );
   
   vtkSTLWriter *external_mold_stl = vtkSTLWriter::New();
   external_mold_stl->SetInput( close->GetOutput() );
   prefix = fileprefix;
   prefix += "-external.stl";
   external_mold_stl->SetFileName( prefix.c_str() );
   external_mold_stl->SetFileTypeToBinary();
   external_mold_stl->Write();
   external_mold_stl->Delete();
   
   triangle->Delete();
   close->Delete();
   pdcf->Delete();
}
//----------------------------------------------------------------------------
void vtk3DSurfaceSTLWidget::ShowMARACASData( marInterface* mar ){
    int whd[3];
    double minmax[2];

    _mar = mar;
    _marImageData = _mar->_experiment->getDynData( )->getVolume( )->castVtk();
    _marImageData->GetDimensions( whd );
    _width  = whd[0];
    _height = whd[1];
    _depth  = whd[2];

    _marImageData->GetScalarRange( minmax );
    
    //Outline v2:
    // An outline provides context around the data.
    _outLine    = vtkOutlineFilter::New( );
    _outLine->SetInput( _marImageData );
    _outMapper  = vtkPolyDataMapper::New( );
    _outMapper->SetInput( _outLine->GetOutput( ) );
    _outActor   = vtkActor::New( );
    _outActor->SetMapper( _outMapper );
    _outActor->GetProperty( )->SetColor( 0.7, 0.0, 0.9 );
    //_outActor->PickableOff( );
    _pRenderer->AddActor( _outActor );

    // Surface
    _mCubes = vtkMarchingCubes::New( );
    _surfMapper = vtkPolyDataMapper::New( );
    _surfActor = vtkActor::New( );

    _mCubes->SetInput( _marImageData );
    _mCubes->SetValue( 0, minmax[1] / 4 );

    vtkStripper *stripper = vtkStripper::New();
    stripper->SetInput( _mCubes->GetOutput( ) );

        _surfMapper->SetInput( stripper->GetOutput() );
    _surfMapper->ScalarVisibilityOff( );
    stripper->Delete();

    _surfActor->SetMapper( _surfMapper );
    _surfActor->PickableOff( );
    _surfActor->GetProperty( )->SetColor( 0.9803, 0.9215, 0.8392 );
    _surfActor->GetProperty( )->SetOpacity( 0.5 );
    _pRenderer->AddActor( _surfActor );

    // 1. ParallelProjectionOn should be set after AddActor (otherwise call vtkRenderer::ResetCameraClippingRange
    // 2. ParallelProjectionOn is *necessary* for the vtkImplicitSelectionLoop
    // otherwise this give a cone instead of a cylinder cutting.
    _pRenderer->GetActiveCamera()->ParallelProjectionOn();

    this->ConfigureVTK( );
}

//----------------------------------------------------------------------------
void vtk3DSurfaceSTLWidget::ShowMARACASDataAndAxe( marInterface* mar ){
  marAxis* temp;
  vtkPolyData* allData;


  this->ShowMARACASData( mar );

  _pRenderer->SetBackground( 0.75, 0.75, 0.75 );

  // Axis
  _mar->_experiment->setAxis( 0 );
  temp = _mar->_experiment->getAxis( );  // ??? getActualAxis ??
  allData = temp->Draw( );
  this->SetAxis( allData );
}

//----------------------------------------------------------------------------
void vtk3DSurfaceSTLWidget::SetInitialPoint( double* pickPoint, double* cameraPos )
{
    gtm::TVector< double > pO( 3 ), pF( 3 ), pp( 3 ), cp( 3 );
    gtm::TVector< double > xc( 3 );
    gtm::TVector< double > x1( 3 ), n1( 3 );
    gtm::TVector< double > x2( 3 ), n2( 3 );
    gtm::TVector< double > tmp( 3 );
    double fac;
    bool success = true;

    if( _centralLineActor )
    {
        _pRenderer->RemoveActor( _centralLineActor );
        _centralLineActor->Delete( );
    } // fi
    if( _centralLineMapper ) _centralLineMapper->Delete( );
    if( _centralLine )       _centralLine->Delete( );

    _centralLine       = NULL;
    _centralLineMapper = NULL;
    _centralLineActor  = NULL;

    for(int i=0; i<4; i++)
    {
      if( _spheresActor[ i ] )
      {
        _pRenderer->RemoveActor( _spheresActor[ i ] );
        _spheresActor[ i ]->Delete( );
      } // fi
      if( _spheresMapper[ i ] ) _spheresMapper[ i ]->Delete( );
      if( _spheres[ i ] )       _spheres[ i ]->Delete( );
      _spheres[ i ] = NULL;
      _spheresMapper[ i ] = NULL;
      _spheresActor[ i ] = NULL;
    } //rof

    fac = GTM_MAX( _width, _height );
    fac = 2 * GTM_MAX( fac, _depth );
    pp( 0 ) = pickPoint[ 0 ]; pp( 1 ) = pickPoint[ 1 ]; pp( 2 ) = pickPoint[ 2 ];
    cp( 0 ) = cameraPos[ 0 ]; cp( 1 ) = cameraPos[ 1 ]; cp( 2 ) = cameraPos[ 2 ];

    if( this->FindCubePointsFromPoints(
      pO.GetAnsiRef( ), pF.GetAnsiRef( ),
      pp.GetAnsiRef( ), cp.GetAnsiRef( )
    ) ) {

        if( this->GetPointAndNormalIntersection(
            x1.GetAnsiRef( ), n1.GetAnsiRef( ),
            pO.GetAnsiRef( ), pF.GetAnsiRef( )
    ) ) {

            if( this->GetPointAndNormalIntersection(
          x2.GetAnsiRef( ), n2.GetAnsiRef( ),
          ( x1 - n1 ).GetAnsiRef( ), ( x1 - ( n1 * fac ) ).GetAnsiRef( )
        ) ) {

                xc = ( x2 + x1 ) * 0.5;
                _mar->_experiment->setStartPoint( (int)xc( 0 ), (int)xc( 1 ), (int)xc( 2 ) );
      
      for(int i=0; i<4; i++)
      {
        _spheres[ i ] = vtkSphereSource::New( );
        _spheresMapper[ i ] = vtkPolyDataMapper::New( );
        _spheresMapper[ i ]->SetInput( _spheres[ i ]->GetOutput( ) );
        _spheresActor[ i ] = vtkActor::New( );
        _spheresActor[ i ]->SetMapper( _spheresMapper[ i ] );
        _spheresActor[ i ]->PickableOff( );
        _pRenderer->AddActor( _spheresActor[ i ] );
      }

                _spheres[ 0 ]->SetCenter( x1( 0 ), x1( 1 ), x1( 2 ) );
                _spheres[ 1 ]->SetCenter( x2( 0 ), x2( 1 ), x2( 2 ) );
                _spheres[ 2 ]->SetCenter( xc( 0 ), xc( 1 ), xc( 2 ) );
                _spheres[ 3 ]->SetCenter( xc( 0 ), xc( 1 ), xc( 2 ) );

                _spheres[ 0 ]->SetRadius( 0.5 );
                _spheres[ 1 ]->SetRadius( 0.5 );
                _spheres[ 2 ]->SetRadius( 0.5 );
                _spheres[ 3 ]->SetRadius( ( xc - x1 ).GetNorm( ) );

                _spheresActor[ 0 ]->GetProperty( )->SetColor( 1.0, 0.0, 0.0 );
                _spheresActor[ 1 ]->GetProperty( )->SetColor( 0.0, 1.0, 0.0 );
                _spheresActor[ 2 ]->GetProperty( )->SetColor( 0.0, 0.0, 1.0 );
                _spheresActor[ 3 ]->GetProperty( )->SetColor( 1.0, 0.0, 0.0 );
                _spheresActor[ 3 ]->GetProperty( )->SetOpacity( 0.3 );

                vtkPoints* points = vtkPoints::New( );
                points->InsertNextPoint( x1.GetAnsiRef( ) );
                points->InsertNextPoint( x2.GetAnsiRef( ) );

                vtkCellArray* array = vtkCellArray::New( );
                array->InsertNextCell( 2 );
                array->InsertCellPoint( 0 );
                array->InsertCellPoint( 1 );

                _centralLine = vtkPolyData::New( );
                _centralLine->SetPoints( points );
                _centralLine->SetLines( array );
      points->Delete();
      array->Delete();

                _centralLineMapper = vtkPolyDataMapper::New( );
                _centralLineMapper->SetInput( _centralLine );

                _centralLineActor = vtkActor::New( );
                _centralLineActor->SetMapper( _centralLineMapper );
                _centralLineActor->GetProperty( )->SetColor( 1.0, 1.0, 1.0 );
                _centralLineActor->PickableOff( );
                _pRenderer->AddActor( _centralLineActor );

            } // fi

        } else success = false;

    } else success = false;

    // Show a message, if any.
    if( !success )
        wxMessageBox(
        _T("Error: Initial point can't be set.\nPlease choose another point."),
        _T("Error"), wxOK | wxCENTRE | wxICON_HAND, this
        );

    // Finish
    _pRenderWindow->Render( );
    InitialSphere = success;
}

//----------------------------------------------------------------------------
void vtk3DSurfaceSTLWidget::SetInitialPoint( double* point )
{
}
//----------------------------------------------------------------------------
void vtk3DSurfaceSTLWidget::ConfigureVTK( )
{
  // get render window
  _pRenderWindow =  this->GetRenderWindow(  );

  // connect renderer and render window and configure render window
  _pRenderWindow->AddRenderer( _pRenderer );

  // configure renderer
  _pRenderer->SetBackground( 0.0, 0.0, 0.0 );
  //_pRenderer->SetBackground( 1, 1, 0);  //FIXME
  _pRenderer->GetActiveCamera( )->Zoom( 1.0 );
  _pRenderer->GetActiveCamera( )->SetClippingRange( 1, 1000 );
}


//----------------------------------------------------------------------------
bool vtk3DSurfaceSTLWidget::IntersectPlaneWithLine( double* p, double* x1, double* x2, double* x3, double* x4, double* x5 )
{
    gtm::TVector< double > vx1( 3 ), vx2( 3 ), vx3( 3 );
    gtm::TVector< double > vx4( 3 ), vx5( 3 ), vx6( p, 3, false );
    gtm::TMatrix< double > mU( 4, 4 ), mD( 4, 4 );
    double t;
    bool ret;

    vx1 = x1;
    vx2 = x2;
    vx3 = x3;
    vx4 = x4;
    vx5 = x5;
    vx6 = vx5 - vx4;

    mD( 0, 0 ) = mU( 0, 0 ) = 1;
    mD( 0, 1 ) = mU( 0, 1 ) = vx1( 0 );
    mD( 0, 2 ) = mU( 0, 2 ) = vx1( 1 );
    mD( 0, 3 ) = mU( 0, 3 ) = vx1( 2 );
    mD( 1, 0 ) = mU( 1, 0 ) = 1;
    mD( 1, 1 ) = mU( 1, 1 ) = vx2( 0 );
    mD( 1, 2 ) = mU( 1, 2 ) = vx2( 1 );
    mD( 1, 3 ) = mU( 1, 3 ) = vx2( 2 );
    mD( 2, 0 ) = mU( 2, 0 ) = 1;
    mD( 2, 1 ) = mU( 2, 1 ) = vx3( 0 );
    mD( 2, 2 ) = mU( 2, 2 ) = vx3( 1 );
    mD( 2, 3 ) = mU( 2, 3 ) = vx3( 2 );
    mU( 3, 0 ) = 1;
    mU( 3, 1 ) = vx4( 0 );
    mU( 3, 2 ) = vx4( 1 );
    mU( 3, 3 ) = vx4( 2 );
    mD( 3, 0 ) = 0;
    mD( 3, 1 ) = vx6( 0 );
    mD( 3, 2 ) = vx6( 1 );
    mD( 3, 3 ) = vx6( 2 );

    ret = ( mD.Det( ) != 0 );
    if( ret ) {

        t = mU.Det( ) / mD.Det( );
        vx6 = ( ( vx4 - vx5 ) * t ) + vx4;

    } // fi
    return( ret );

}

//----------------------------------------------------------------------------
bool vtk3DSurfaceSTLWidget::FindCubePointsFromPoints( double* pO, double* pF, double* pickPoint, double* cameraPos )
{
    gtm::TVector< double > p1( 3 ), p2( 3 ), p3( 3 ), p4( 3 ), p5( 3 ), p6( 3 );
    gtm::TVector< double > c( 3 );
    gtm::TVector< double >* swp;
    gtm::TVector< double > tpO( pO, 3, false ), tpF( pF, 3, false );
    std::vector< gtm::TVector< double >* > points;
    double x1[ 3 ], x2[ 3 ], x3[ 3 ], d1, d2;
    int i, j;

    // 1st plane intersection
    x1[ 0 ] = 0; x1[ 1 ] = _height; x1[ 2 ] =      0;
    x2[ 0 ] = 0; x2[ 1 ] =       0; x2[ 2 ] =      0;
    x3[ 0 ] = 0; x3[ 1 ] =       0; x3[ 2 ] = _depth;
    if( this->IntersectPlaneWithLine( p1.GetAnsiRef( ), x1, x2, x3, cameraPos, pickPoint ) )
        points.push_back( &p1 );

    // 2nd plane intersection
    x1[ 0 ] =      0; x1[ 1 ] = _height; x1[ 2 ] = _depth;
    x2[ 0 ] =      0; x2[ 1 ] =       0; x2[ 2 ] = _depth;
    x3[ 0 ] = _width; x3[ 1 ] =       0; x3[ 2 ] = _depth;
    if( this->IntersectPlaneWithLine( p2.GetAnsiRef( ), x1, x2, x3, cameraPos, pickPoint ) )
        points.push_back( &p2 );

    // 3rd plane intersection
    x1[ 0 ] =      0; x1[ 1 ] = 0; x1[ 2 ] = _depth;
    x2[ 0 ] =      0; x2[ 1 ] = 0; x2[ 2 ] =      0;
    x3[ 0 ] = _width; x3[ 1 ] = 0; x3[ 2 ] =      0;
    if( this->IntersectPlaneWithLine( p3.GetAnsiRef( ), x1, x2, x3, cameraPos, pickPoint ) )
        points.push_back( &p3 );

    // 4th plane intersection
    x1[ 0 ] = _width; x1[ 1 ] = _height; x1[ 2 ] = _depth;
    x2[ 0 ] = _width; x2[ 1 ] =       0; x2[ 2 ] = _depth;
    x3[ 0 ] = _width; x3[ 1 ] =       0; x3[ 2 ] =      0;
    if( this->IntersectPlaneWithLine( p4.GetAnsiRef( ), x1, x2, x3, cameraPos, pickPoint ) )
        points.push_back( &p4 );

    // 5th plane intersection
    x1[ 0 ] = _width; x1[ 1 ] =       0; x1[ 2 ] = 0;
    x2[ 0 ] =      0; x2[ 1 ] =       0; x2[ 2 ] = 0;
    x3[ 0 ] =      0; x3[ 1 ] = _height; x3[ 2 ] = 0;
    if( this->IntersectPlaneWithLine( p5.GetAnsiRef( ), x1, x2, x3, cameraPos, pickPoint ) )
        points.push_back( &p5 );

    // 6th plane intersection
    x1[ 0 ] =      0; x1[ 1 ] = _height; x1[ 2 ] =      0;
    x2[ 0 ] =      0; x2[ 1 ] = _height; x2[ 2 ] = _depth;
    x3[ 0 ] = _width; x3[ 1 ] = _height; x3[ 2 ] = _depth;
    if( this->IntersectPlaneWithLine( p6.GetAnsiRef( ), x1, x2, x3, cameraPos, pickPoint ) )
        points.push_back( &p6 );

    if( points.size( ) >= 2 ) { // Did I find at least 2 points?

        c( 0 ) = ( double )_width  / 2.0;
        c( 1 ) = ( double )_height / 2.0;
        c( 2 ) = ( double )_depth  / 2.0;

        // Sort with bubble sort. Only 30 iterations!
        for( i = 0; i < points.size( ); i++ ) {
            for( j = 0; j < points.size( ) - 1; j++ ) {

                d1 = ( c - *points[ j ] ).GetNorm( );
                d2 = ( c - *points[ j + 1 ] ).GetNorm( );
                if( d2 < d1 ) {

                    swp = points[ j ];
                    points[ j ] = points[ j + 1 ];
                    points[ j + 1 ] = swp;

                } // fi

    } // rof

        } // rof

        // Order the two points according to distance to camera.
        c = cameraPos;
        d1 = ( c - *points[ 0 ] ).GetNorm( );
        d2 = ( c - *points[ 1 ] ).GetNorm( );
        tpO = ( d1 < d2 )? *points[ 0 ]: *points[ 1 ];
        tpF = ( d1 > d2 )? *points[ 0 ]: *points[ 1 ];
        return( true );

    } else return( false );

}

//----------------------------------------------------------------------------
bool vtk3DSurfaceSTLWidget::GetPointAndNormalIntersection( double* p, double* n, double* pO, double* pF )
{
    gtm::TVector< double > n1( 3 ), n2( 3 ), n3( 3 );
    int subId, cellId, returnVal;
    double t, pcoords[ 3 ], x[ 3 ];
    double fpO[ 3 ], fpF[ 3 ];
    double p1[ 3 ], p2[ 3 ], p3[ 3 ];
    vtkPolyData* data = _mCubes->GetOutput( );
    vtkCellLocator* locator = vtkCellLocator::New( );

    locator->SetDataSet( data );
    locator->Initialize( );
    locator->Update( );

    fpO[ 0 ] = pO[ 0 ]; fpO[ 1 ] = pO[ 1 ]; fpO[ 2 ] = pO[ 2 ];
    fpF[ 0 ] = pF[ 0 ]; fpF[ 1 ] = pF[ 1 ]; fpF[ 2 ] = pF[ 2 ];
    returnVal = locator->IntersectWithLine( fpO, fpF, 0.1, t, x, pcoords, subId, cellId );
    locator->Delete( );

    if( returnVal )
    {
        vtkCell* cell = data->GetCell( cellId );
        vtkPoints* points = cell->GetPoints( );

        data->GetPointData( )->GetNormals( )->GetTuple( cell->GetPointIds( )->GetId( 0 ), n1.GetAnsiRef( ) );
        data->GetPointData( )->GetNormals( )->GetTuple( cell->GetPointIds( )->GetId( 1 ), n2.GetAnsiRef( ) );
        data->GetPointData( )->GetNormals( )->GetTuple( cell->GetPointIds( )->GetId( 2 ), n3.GetAnsiRef( ) );
   
        n1 += n2 + n3;
        n1 *= ( 1.0 / 3.0 );
        n1.Normalize( );
        n[ 0 ] = n1( 0 ); n[ 1 ] = n1( 1 ); n[ 2 ] = n1( 2 );

        points->GetPoint( 0, p1 );
        points->GetPoint( 1, p2 );
        points->GetPoint( 2, p3 );
        this->IntersectPlaneWithLine( p, p1, p2, p3, pO, pF );
        return( true );
    } else return( false );

}
//----------------------------------------------------------------------------
vtkPolyData* vtk3DSurfaceSTLWidget::ConvertMarAxisToPolyData()
{
  double p1[marAxis::INDX_count];
  //retrieve the actual axis
  marAxis* ax = _mar->_experiment->getAxis();
  //Number of axis points
  int numPoints = ax->getNumberOfControlPoints();

  vtkPolyData *axis = vtkPolyData::New();
  vtkPoints *points = vtkPoints::New();
  vtkFloatArray *scalars = vtkFloatArray::New();
  scalars->SetNumberOfComponents (2);
  //Field are : 1. radius (=INDX_RAYON) , 2. average (=INDX_SIGNALVALUE)

  double tmp1,tmp2;
  for(int i=0; i<numPoints; i++){
    ax->getControlPoint(i, p1, p1+3);
    points->InsertPoint(i, p1);
        tmp1 = p1[marAxis::INDX_RAYON];
//      tmp2 = p1[marAxis::INDX_SIGNALVALUE];
        tmp2 = (double)ax->getSignal(i);
    scalars->InsertTuple2(i,tmp1, tmp2);

  }

  // We now assign the pieces to the vtkPolyData.
  axis->SetPoints(points);
  points->Delete();
  axis->GetPointData()->SetScalars(scalars);
  scalars->Delete();
 
  return axis;
}
//----------------------------------------------------------------------------
void vtk3DSurfaceSTLWidget::ConstructVessel( )
{
  vtkPolyData *axis = this->ConvertMarAxisToPolyData( );

  //retrieve the vtk image data

  if( !_psc) this->_psc = vtkImagePolyDataSeedConnectivity::New();
  this->_psc->SetInput( _marImageData );
  this->_psc->SetAxis( axis );
  this->_psc->Update();
  axis->Delete();

/*  vtkStripper *strip = vtkStripper::New();
  strip->SetInput( this->_psc->GetOutput() );  //GetOutput() -> inner mold
  strip->Update(); //Important before a ShallowCopy !!!
  
  if( !_stlInternalVessel) _stlInternalVessel = vtkPolyData::New();
  _stlInternalVessel->ShallowCopy( strip->GetOutput() );
  strip->Delete();*/
  _stlInternalVessel = this->_psc->GetOutput();

/*  if( !_stlExternalVessel) _stlExternalVessel = vtkPolyData::New();
  _stlExternalVessel->ShallowCopy( this->_psc->GetOuterMold() );
  //this->_psc->GetOuterMold()->Delete();*/
  _stlExternalVessel = this->_psc->GetOuterMold();

  vtkPolyDataMapper *dsm1 = vtkPolyDataMapper ::New();
  dsm1->SetInput ( _stlInternalVessel ); //fasten stuff
  dsm1->ScalarVisibilityOff();

  _actorVessel = vtkActor::New();
  _actorVessel->SetMapper (dsm1);

  //improve visibility:
  _actorVessel->GetProperty()->SetColor (1,0,0);

  vtkPolyDataMapper *dsm2 = vtkPolyDataMapper ::New();
  dsm2->SetInput ( _stlExternalVessel );
  dsm2->ScalarVisibilityOff();

  _actorExternalVessel = vtkActor::New();
  _actorExternalVessel->SetMapper (dsm2);

  //improve visibility:
  _actorExternalVessel->GetProperty()->SetRepresentationToWireframe();

  _pRenderer->AddActor( _actorVessel );
  _pRenderer->AddActor( _actorExternalVessel );
  _pRenderWindow->Render( );

  dsm1->Delete();
  dsm2->Delete();
}
//----------------------------------------------------------------------------
void vtk3DSurfaceSTLWidget::SetSTLThresholdRatio(double ratio)
{
  this->_psc->SetThresholdRatio( (double)(ratio/100) );
  _pRenderWindow->Render( );
}
//----------------------------------------------------------------------------
float vtk3DSurfaceSTLWidget::GetSTLThreshold()
{
  return  100*this->_psc->GetThresholdRatio( );
}
//----------------------------------------------------------------------------