]> Creatis software - creaMaracasVisu.git/blob - lib/maracasVisuLib/src/interface/wxWindows/widgets/include/vtk3DSurfaceSTLWidget.cxx
*** empty log message ***
[creaMaracasVisu.git] / 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 //----------------------------------------------------------------------------