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

  Program:   wxMaracas
  Module:    $RCSfile: vtk3DSurfaceSTLWidget.cxx,v $
  Language:  C++
  Date:      $Date: 2008/10/31 16:32:41 $
  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"



/**
 * Again wxVTK is an hybrid class, double click was done using wxWindows, simply
 * because VTK doesn't provide such facility (as opposed to wheel support that is
 * supposed to be merged in VTK trunk sooner or later:
 *
 * http://public.kitware.com/pipermail/vtkusers/2003-August/019548.html
 * http://www.creatis.insa-lyon.fr/~malaterre/vtk/wheel.patch
 */

//----------------------------------------------------------------------------
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)
{
  ///\todo : STL visibility: both internal and external mold:
  _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
    {
   /**
    * Copy / paste from vtkPolyDataCutter
    */
      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 );

}
//----------------------------------------------------------------------------
 /**
  * What we are doing here is from a LFV marAxis to turn it into a vtkPolyData
  * Therefore it should be a lot more easier to integrate in a real vtk class
  * for filtering : vtkMaracasImageToSTL
  */
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( );
}
//----------------------------------------------------------------------------