1 /*# ---------------------------------------------------------------------
3 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
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
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.
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
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 # ------------------------------------------------------------------------ */
26 /*=========================================================================
29 Module: $RCSfile: vtk3DSurfaceSTLWidget.cxx,v $
31 Date: $Date: 2012/11/15 14:15:18 $
32 Version: $Revision: 1.2 $
34 Copyright: (c) 2002, 2003
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.
41 =========================================================================*/
43 #include "vtkClosePolyData.h"
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>
53 #include <vtkCellLocator.h>
54 #include <vtkInteractorStyleTrackballCamera.h>
55 #include <vtkCellArray.h>
56 #include <vtkStripper.h>
58 #include <vtkImplicitSelectionLoop.h>
59 #include <vtkExtractPolyDataGeometry.h>
60 #include <vtkPolyDataConnectivityFilter.h>
61 #include <vtkPolyDataNormals.h>
63 #include "vtk3DSurfaceSTLWidget.h"
66 #include "marInterface.h"
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:
76 * http://public.kitware.com/pipermail/vtkusers/2003-August/019548.html
77 * http://www.creatis.insa-lyon.fr/~malaterre/vtk/wheel.patch
80 //----------------------------------------------------------------------------
81 BEGIN_EVENT_TABLE( vtk3DSurfaceSTLWidget, wxVTKRenderWindowInteractor )
82 //EED EVT_LEFT_DCLICK( vtk3DSurfaceSTLWidget::OnLeftDClick )
85 //----------------------------------------------------------------------------
86 vtk3DSurfaceSTLWidget::vtk3DSurfaceSTLWidget(
93 : wxVTKRenderWindowInteractor( parent, id, pos, size, style, name ){
96 //make current interactor style be trackbal camera, asked by radiologists.
97 // vtkInteractorStyleSwitch *iss = dynamic_cast<vtkInteractorStyleSwitch*>(this->GetInteractorStyle());
98 // iss->SetCurrentStyleToTrackballCamera();
101 vtkInteractorStyle3DMaracas *interactorStyle3DMaracas = vtkInteractorStyle3DMaracas::New();
102 interactorStyle3DMaracas->SetInteractor (this);
103 this->SetInteractorStyle(interactorStyle3DMaracas);
106 _pRenderer = vtkRenderer::New( );
109 _centralLineMapper = NULL;
110 _centralLineActor = NULL;
113 for(int i=0; i<4; i++)
115 _spheres[ i ] = NULL;
116 _spheresMapper[ i ] = NULL;
117 _spheresActor[ i ] = NULL;
122 _stlExternalVessel = _stlInternalVessel = NULL;
124 _actorExternalVessel= _actorVessel = NULL;
127 //----------------------------------------------------------------------------
128 vtk3DSurfaceSTLWidget::~vtk3DSurfaceSTLWidget()
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++)
150 if( _spheres[ i ] ) _spheres[ i ] ->Delete( );
151 if( _spheresMapper[ i ] ) _spheresMapper[ i ] ->Delete( );
152 if( _spheresActor[ i ] ) _spheresActor[ i ] ->Delete( );
154 if( _pRenderer ) _pRenderer ->Delete( );
157 //----------------------------------------------------------------------------
158 void vtk3DSurfaceSTLWidget::SetSurfaceColor( float red, float green, float blue )
160 _surfActor->GetProperty()->SetColor(red, green, blue );
161 _pRenderWindow->Render();
163 //----------------------------------------------------------------------------
164 void vtk3DSurfaceSTLWidget::SetSurfaceVisibility( bool visible )
166 _surfActor->SetVisibility( visible );
167 _pRenderWindow->Render();
169 //----------------------------------------------------------------------------
170 void vtk3DSurfaceSTLWidget::SetSTLSurfaceVisibility( bool intvisible , bool extvisible)
172 ///\todo : STL visibility: both internal and external mold:
173 _actorVessel->SetVisibility( intvisible );
174 _actorExternalVessel->SetVisibility( extvisible );
175 _pRenderWindow->Render();
177 //----------------------------------------------------------------------------
178 void vtk3DSurfaceSTLWidget::SetSurfaceIsoValue( int isoval )
180 _mCubes->SetValue(0, isoval);
181 _pRenderWindow->Render();
183 //----------------------------------------------------------------------------
184 void vtk3DSurfaceSTLWidget::SetSurfaceOpacity( int opaval )
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();
190 //----------------------------------------------------------------------------
191 void vtk3DSurfaceSTLWidget::GetSphereCenter( double center[3] )
193 _spheres[ 3 ]->GetCenter( center );
195 //----------------------------------------------------------------------------
196 void vtk3DSurfaceSTLWidget::SetAxis( vtkPolyData *axis )
198 _axesMapper = vtkPolyDataMapper::New( );
199 _axesMapper->SetInput( axis );
202 _axesActor = vtkActor::New( );
203 _axesActor->SetMapper( _axesMapper );
204 _axesActor->GetProperty()->SetColor( 1, 0, 0 );
205 _pRenderer->AddActor( _axesActor );
207 _pRenderWindow->Render( );
209 //----------------------------------------------------------------------------
210 void vtk3DSurfaceSTLWidget::RemoveAxis( )
214 _pRenderer->RemoveActor(_axesActor);
215 _axesActor->Delete();
217 _pRenderWindow->Render( );
220 //----------------------------------------------------------------------------
221 void vtk3DSurfaceSTLWidget::SetCuttingMode( bool mode )
225 //start rubber/cutter mode
228 _iasc = vtkInteractorStyleCutter::New();
230 this->SetInteractorStyle( _iasc );
235 * Copy / paste from vtkPolyDataCutter
239 _iasc->VisibilityOff(); //turn off the vtkInteractorStyleCutter
240 vtkImplicitSelectionLoop *loop = vtkImplicitSelectionLoop::New();
241 loop->SetLoop( _iasc->GetLoopPoints() );
242 loop->SetNormal( _iasc->GetDirection() );
244 vtkTriangleFilter *tf = vtkTriangleFilter::New();
245 tf->SetInput( _stlInternalVessel );
247 vtkExtractPolyDataGeometry *extract = vtkExtractPolyDataGeometry::New ();
248 extract->SetInput ( tf->GetOutput());
249 extract->SetImplicitFunction ( loop );
250 extract->ExtractInsideOff ();
251 extract->ExtractBoundaryCellsOff ();
255 //#connectivity filter to keep only the largest part
256 vtkPolyDataConnectivityFilter *connect = vtkPolyDataConnectivityFilter::New();
257 connect->SetInput( extract->GetOutput() );
258 connect->SetExtractionModeToLargestRegion ();
261 vtkClosePolyData *close = vtkClosePolyData::New();
262 close->SetInput( connect->GetOutput() );
265 vtkStripper *strip = vtkStripper::New();
266 strip->SetInput( close->GetOutput() );
269 //If vtkClosePolyData was well written we wouldn't had to recompute the normals
270 vtkPolyDataNormals *normals = vtkPolyDataNormals::New();
271 normals->SetInput( strip->GetOutput() );
274 /*vtkPolyDataWriter *polywriter = vtkPolyDataWriter::New();
275 polywriter->SetInput( normals->GetOutput() );
276 polywriter->SetFileName( "cutter.vtk" );
277 polywriter->SetFileTypeToBinary();
279 polywriter->Delete();*/
281 normals->Update(); //important before a ShallowCopy Or DeepCopy
282 _stlInternalVessel->DeepCopy( normals->GetOutput() ); //Use DeepCopy when using ShallowCopy afterwards
284 tf->SetInput( _stlExternalVessel );
285 normals->Update(); //important before a ShallowCopy
286 _stlExternalVessel->ShallowCopy( normals->GetOutput() );
287 _pRenderWindow->Render( );
289 //make current interactor style be trackbal camera, asked by radiologists.
290 vtkInteractorStyleTrackballCamera *istc = vtkInteractorStyleTrackballCamera::New();
291 this->SetInteractorStyle( istc );
296 //----------------------------------------------------------------------------
297 void vtk3DSurfaceSTLWidget::ExportSurfaceAsSTL( const char* fileprefix )
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() );
312 vtkSTLWriter *internal_mold_stl = vtkSTLWriter::New();
313 internal_mold_stl->SetInput( close->GetOutput() );
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();
321 //convert to triangles first:
322 triangle->SetInput( _stlExternalVessel );
324 vtkSTLWriter *external_mold_stl = vtkSTLWriter::New();
325 external_mold_stl->SetInput( close->GetOutput() );
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();
337 //----------------------------------------------------------------------------
338 void vtk3DSurfaceSTLWidget::ShowMARACASData( marInterface* mar ){
343 _marImageData = _mar->_experiment->getDynData( )->getVolume( )->castVtk();
344 _marImageData->GetDimensions( whd );
349 _marImageData->GetScalarRange( minmax );
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 );
364 _mCubes = vtkMarchingCubes::New( );
365 _surfMapper = vtkPolyDataMapper::New( );
366 _surfActor = vtkActor::New( );
368 _mCubes->SetInput( _marImageData );
369 _mCubes->SetValue( 0, minmax[1] / 4 );
371 vtkStripper *stripper = vtkStripper::New();
372 stripper->SetInput( _mCubes->GetOutput( ) );
374 _surfMapper->SetInput( stripper->GetOutput() );
375 _surfMapper->ScalarVisibilityOff( );
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 );
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();
389 this->ConfigureVTK( );
392 //----------------------------------------------------------------------------
393 void vtk3DSurfaceSTLWidget::ShowMARACASDataAndAxe( marInterface* mar ){
395 vtkPolyData* allData;
398 this->ShowMARACASData( mar );
400 _pRenderer->SetBackground( 0.75, 0.75, 0.75 );
403 _mar->_experiment->setAxis( 0 );
404 temp = _mar->_experiment->getAxis( ); // ??? getActualAxis ??
405 allData = temp->Draw( );
406 this->SetAxis( allData );
409 //----------------------------------------------------------------------------
410 void vtk3DSurfaceSTLWidget::SetInitialPoint( double* pickPoint, double* cameraPos )
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 );
420 if( _centralLineActor )
422 _pRenderer->RemoveActor( _centralLineActor );
423 _centralLineActor->Delete( );
425 if( _centralLineMapper ) _centralLineMapper->Delete( );
426 if( _centralLine ) _centralLine->Delete( );
429 _centralLineMapper = NULL;
430 _centralLineActor = NULL;
432 for(int i=0; i<4; i++)
434 if( _spheresActor[ i ] )
436 _pRenderer->RemoveActor( _spheresActor[ i ] );
437 _spheresActor[ i ]->Delete( );
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;
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 ];
451 if( this->FindCubePointsFromPoints(
452 pO.GetAnsiRef( ), pF.GetAnsiRef( ),
453 pp.GetAnsiRef( ), cp.GetAnsiRef( )
456 if( this->GetPointAndNormalIntersection(
457 x1.GetAnsiRef( ), n1.GetAnsiRef( ),
458 pO.GetAnsiRef( ), pF.GetAnsiRef( )
461 if( this->GetPointAndNormalIntersection(
462 x2.GetAnsiRef( ), n2.GetAnsiRef( ),
463 ( x1 - n1 ).GetAnsiRef( ), ( x1 - ( n1 * fac ) ).GetAnsiRef( )
466 xc = ( x2 + x1 ) * 0.5;
467 _mar->_experiment->setStartPoint( (int)xc( 0 ), (int)xc( 1 ), (int)xc( 2 ) );
469 for(int i=0; i<4; i++)
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 ] );
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 ) );
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( ) );
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 );
496 vtkPoints* points = vtkPoints::New( );
497 points->InsertNextPoint( x1.GetAnsiRef( ) );
498 points->InsertNextPoint( x2.GetAnsiRef( ) );
500 vtkCellArray* array = vtkCellArray::New( );
501 array->InsertNextCell( 2 );
502 array->InsertCellPoint( 0 );
503 array->InsertCellPoint( 1 );
505 _centralLine = vtkPolyData::New( );
506 _centralLine->SetPoints( points );
507 _centralLine->SetLines( array );
511 _centralLineMapper = vtkPolyDataMapper::New( );
512 _centralLineMapper->SetInput( _centralLine );
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 );
522 } else success = false;
524 } else success = false;
526 // Show a message, if any.
529 _T("Error: Initial point can't be set.\nPlease choose another point."),
530 _T("Error"), wxOK | wxCENTRE | wxICON_HAND, this
534 _pRenderWindow->Render( );
535 InitialSphere = success;
538 //----------------------------------------------------------------------------
539 void vtk3DSurfaceSTLWidget::SetInitialPoint( double* point )
542 //----------------------------------------------------------------------------
543 void vtk3DSurfaceSTLWidget::ConfigureVTK( )
546 _pRenderWindow = this->GetRenderWindow( );
548 // connect renderer and render window and configure render window
549 _pRenderWindow->AddRenderer( _pRenderer );
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 );
559 //----------------------------------------------------------------------------
560 bool vtk3DSurfaceSTLWidget::IntersectPlaneWithLine( double* p, double* x1, double* x2, double* x3, double* x4, double* x5 )
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 );
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 );
588 mU( 3, 1 ) = vx4( 0 );
589 mU( 3, 2 ) = vx4( 1 );
590 mU( 3, 3 ) = vx4( 2 );
592 mD( 3, 1 ) = vx6( 0 );
593 mD( 3, 2 ) = vx6( 1 );
594 mD( 3, 3 ) = vx6( 2 );
596 ret = ( mD.Det( ) != 0 );
599 t = mU.Det( ) / mD.Det( );
600 vx6 = ( ( vx4 - vx5 ) * t ) + vx4;
607 //----------------------------------------------------------------------------
608 bool vtk3DSurfaceSTLWidget::FindCubePointsFromPoints( double* pO, double* pF, double* pickPoint, double* cameraPos )
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;
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 );
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 );
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 );
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 );
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 );
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 );
660 if( points.size( ) >= 2 ) { // Did I find at least 2 points?
662 c( 0 ) = ( double )_width / 2.0;
663 c( 1 ) = ( double )_height / 2.0;
664 c( 2 ) = ( double )_depth / 2.0;
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++ ) {
670 d1 = ( c - *points[ j ] ).GetNorm( );
671 d2 = ( c - *points[ j + 1 ] ).GetNorm( );
675 points[ j ] = points[ j + 1 ];
676 points[ j + 1 ] = swp;
684 // Order the two points according to distance to camera.
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 ];
692 } else return( false );
696 //----------------------------------------------------------------------------
697 bool vtk3DSurfaceSTLWidget::GetPointAndNormalIntersection( double* p, double* n, double* pO, double* pF )
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( );
707 locator->SetDataSet( data );
708 locator->Initialize( );
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 );
718 vtkCell* cell = data->GetCell( cellId );
719 vtkPoints* points = cell->GetPoints( );
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( ) );
728 n[ 0 ] = n1( 0 ); n[ 1 ] = n1( 1 ); n[ 2 ] = n1( 2 );
730 points->GetPoint( 0, p1 );
731 points->GetPoint( 1, p2 );
732 points->GetPoint( 2, p3 );
733 this->IntersectPlaneWithLine( p, p1, p2, p3, pO, pF );
735 } else return( false );
738 //----------------------------------------------------------------------------
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
744 vtkPolyData* vtk3DSurfaceSTLWidget::ConvertMarAxisToPolyData()
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();
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)
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);
769 // We now assign the pieces to the vtkPolyData.
770 axis->SetPoints(points);
772 axis->GetPointData()->SetScalars(scalars);
777 //----------------------------------------------------------------------------
778 void vtk3DSurfaceSTLWidget::ConstructVessel( )
780 vtkPolyData *axis = this->ConvertMarAxisToPolyData( );
782 //retrieve the vtk image data
784 if( !_psc) this->_psc = vtkImagePolyDataSeedConnectivity::New();
785 this->_psc->SetInput( _marImageData );
786 this->_psc->SetAxis( axis );
787 this->_psc->Update();
790 /* vtkStripper *strip = vtkStripper::New();
791 strip->SetInput( this->_psc->GetOutput() ); //GetOutput() -> inner mold
792 strip->Update(); //Important before a ShallowCopy !!!
794 if( !_stlInternalVessel) _stlInternalVessel = vtkPolyData::New();
795 _stlInternalVessel->ShallowCopy( strip->GetOutput() );
797 _stlInternalVessel = this->_psc->GetOutput();
799 /* if( !_stlExternalVessel) _stlExternalVessel = vtkPolyData::New();
800 _stlExternalVessel->ShallowCopy( this->_psc->GetOuterMold() );
801 //this->_psc->GetOuterMold()->Delete();*/
802 _stlExternalVessel = this->_psc->GetOuterMold();
804 vtkPolyDataMapper *dsm1 = vtkPolyDataMapper ::New();
805 dsm1->SetInput ( _stlInternalVessel ); //fasten stuff
806 dsm1->ScalarVisibilityOff();
808 _actorVessel = vtkActor::New();
809 _actorVessel->SetMapper (dsm1);
811 //improve visibility:
812 _actorVessel->GetProperty()->SetColor (1,0,0);
814 vtkPolyDataMapper *dsm2 = vtkPolyDataMapper ::New();
815 dsm2->SetInput ( _stlExternalVessel );
816 dsm2->ScalarVisibilityOff();
818 _actorExternalVessel = vtkActor::New();
819 _actorExternalVessel->SetMapper (dsm2);
821 //improve visibility:
822 _actorExternalVessel->GetProperty()->SetRepresentationToWireframe();
824 _pRenderer->AddActor( _actorVessel );
825 _pRenderer->AddActor( _actorExternalVessel );
826 _pRenderWindow->Render( );
831 //----------------------------------------------------------------------------
832 void vtk3DSurfaceSTLWidget::SetSTLThresholdRatio(double ratio)
834 this->_psc->SetThresholdRatio( (double)(ratio/100) );
835 _pRenderWindow->Render( );
837 //----------------------------------------------------------------------------
838 float vtk3DSurfaceSTLWidget::GetSTLThreshold()
840 return 100*this->_psc->GetThresholdRatio( );
842 //----------------------------------------------------------------------------