1 // =========================================================================
2 // @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
3 // =========================================================================
5 #include "CTArteries.h"
6 #include "ui_CTArteries.h"
7 #include "ui_Parameters.h"
8 #include "FileDialog.h"
12 #include <itkImageFileReader.h>
13 #include <itkImageSeriesReader.h>
14 #include <itkQuadEdgeMesh.h>
16 #include <vtkMarchingCubes.h>
17 #include <vtkMarchingSquares.h>
18 #include <vtkPolyDataConnectivityFilter.h>
19 #include <vtkPolyDataMapper.h>
20 #include <vtkProperty.h>
21 #include <vtkRenderer.h>
22 #include <vtkRenderWindow.h>
23 #include <vtkSplineWidget.h>
27 #include <ivq/ITK/CPRImageFilter.h>
28 #include <ivq/VTK/ImageActor.h>
29 #include <ivq/VTK/ImageViewer.h>
30 #include <ivq/VTK/MPRViewers.h>
31 #include <ivq/VTK/PolyDataActor.h>
32 #include <ivq/VTK/SeedWidgetOverImageActor.h>
33 #include <ivq/Qt/DicomSeriesSelectorDialog.h>
34 #include <ivq/Qt/RendererWidget.h>
36 #include "algorithms/RandomWalkSegmentation.h"
38 // -------------------------------------------------------------------------
39 #define QT_ACTION_CONN( x ) \
41 this->m_UI->a##x, SIGNAL( triggered( ) ), this, SLOT( _s##x( ) ) \
44 // -------------------------------------------------------------------------
49 static ContourCallback* New( )
51 return( new ContourCallback( ) );
54 vtkObject* caller, unsigned long eId, void* data
57 ivq::VTK::ImageActor* actor =
58 ivq::VTK::ImageActor::SafeDownCast( caller );
59 if( eId == vtkCommand::InteractionEvent && actor != NULL )
61 ivq::VTK::ImageViewer* viewer = this->Window->m_UI->Slice->GetViewer( );
62 if( this->Window->m_ContourActor != NULL )
63 viewer->GetRenderer( )->RemoveViewProp( this->Window->m_ContourActor->GetActor( ) );
66 actor->GetActorBounds( abnds );
67 unsigned int nSlice = actor->GetSliceNumber( );
68 for( unsigned int i = 0; i < 16; ++i )
70 double w = double( 6.28 ) * double( i ) / double( 16 );
71 CTArteries::TFourier::TPoint v = this->Window->m_Fourier[ nSlice ]( w );
72 this->Window->m_Spline->SetHandlePosition(
73 i, v[ 0 ], v[ 1 ], abnds[ 4 ]
77 this->Window->m_Spline->SetProjectionPosition(
79 actor->GetBounds( )[ 4 ] +
80 actor->GetBounds( )[ 5 ]
85 this->Window->m_ContourActor =
86 this->Window->m_ContoursActors[ actor->GetSliceNumber( ) ];
87 viewer->GetRenderer( )->AddViewProp( this->Window->m_ContourActor->GetActor( ) );
98 virtual ~ContourCallback( )
103 // Purposely not implemented
104 ContourCallback( const ContourCallback& other );
105 ContourCallback& operator=( const ContourCallback& other );
110 // -------------------------------------------------------------------------
112 CTArteries( int argc, char* argv[], QWidget* parent )
113 : Superclass( parent ),
114 m_UI( new Ui::CTArteries ),
115 m_UIParameters( new Ui::Parameters )
117 this->m_UI->setupUi( this );
118 this->m_DlgParameters = new QDialog( );
119 this->m_UIParameters->setupUi( this->m_DlgParameters );
121 // Connect signals with slots
122 QT_ACTION_CONN( Open );
123 QT_ACTION_CONN( Config );
124 QT_ACTION_CONN( Process );
129 std::string cmd( argv[ 1 ] );
130 std::string fname( argv[ 2 ] );
132 this->_ExecuteLog( fname );
137 // -------------------------------------------------------------------------
142 delete this->m_UIParameters;
143 delete this->m_DlgParameters;
146 // -------------------------------------------------------------------------
147 template< class _TStrings >
149 _openImage( const _TStrings& fnames )
152 std::ofstream str_log( "CTArteries.log" );
153 str_log << fnames.size( ) << std::endl;
154 for( const std::string& s: fnames )
155 str_log << s << std::endl;
158 // Create appropriate reader
159 itk::ImageSource< TImage >::Pointer reader;
160 if( fnames.size( ) == 1 )
162 itk::ImageFileReader< TImage >::Pointer r =
163 itk::ImageFileReader< TImage >::New( );
164 r->SetFileName( *( fnames.begin( ) ) );
167 else if( fnames.size( ) > 1 )
169 itk::ImageSeriesReader< TImage >::Pointer r =
170 itk::ImageSeriesReader< TImage >::New( );
171 typename _TStrings::const_iterator i;
172 for( i = fnames.begin( ); i != fnames.end( ); ++i )
173 r->AddFileName( *i );
179 if( reader.IsNull( ) )
184 this->m_Image = reader->GetOutput( );
185 this->m_Image->DisconnectPipeline( );
187 catch( std::exception& err )
189 this->m_Image = NULL;
190 QMessageBox::critical( this, "Error opening image", err.what( ) );
193 this->_showInputImage( );
196 // -------------------------------------------------------------------------
198 _openDicom( const std::string& dirname )
200 ivq::Qt::DicomSeriesSelectorDialog dlg( this );
201 dlg.setStartDir( dirname );
202 if( dlg.exec( ) == QDialog::Accepted )
203 this->_openImage( *( dlg.selectedFilenames( ) ) );
206 // -------------------------------------------------------------------------
210 if( this->m_Image.IsNotNull( ) )
212 this->m_VTKImage = TVTKImage::New( );
213 this->m_VTKImage->SetInput( this->m_Image );
214 this->m_VTKImage->Update( );
216 ivq::VTK::MPRViewers* viewers = this->m_UI->MPR->GetViewers( );
217 viewers->SetInputData( this->m_VTKImage->GetOutput( ) );
218 viewers->SetColorWindow( 1000 );
219 viewers->SetColorLevel( 100 );
221 ivq::VTK::ImageViewer* view = viewers->GetView( 2 );
223 vtkSmartPointer< ivq::VTK::SeedWidgetOverImageActor >::New( );
224 this->m_Seeds->SetActor( view->GetImageActor( ) );
225 this->m_Seeds->SetInteractor( view->GetRenderWindow( )->GetInteractor( ) );
228 viewers->ResetCameras( );
229 this->m_Seeds->EnabledOn( );
234 this->m_VTKImage = NULL;
239 // -------------------------------------------------------------------------
241 _process( const std::vector< TImage::PointType >& seeds )
244 std::ofstream str_log( "CTArteries.log", std::ios_base::app );
245 str_log << this->m_UIParameters->Beta->value( ) << std::endl;
246 str_log << this->m_UIParameters->Sigma->value( ) << std::endl;
247 str_log << this->m_UIParameters->Radius->value( ) << std::endl;
248 str_log << seeds.size( ) << std::endl;
249 for( TImage::PointType seed: seeds )
250 str_log << seed << std::endl;
254 typedef RandomWalkSegmentation< TImage, TScalarImage > _TSegmentation;
255 _TSegmentation::Pointer seg = _TSegmentation::New( );
256 seg->SetInput( this->m_Image );
257 seg->SetBeta( this->m_UIParameters->Beta->value( ) );
258 seg->SetSigma( this->m_UIParameters->Sigma->value( ) );
259 seg->SetRadius( this->m_UIParameters->Radius->value( ) );
260 seg->SetStartSeed( seeds[ 0 ] );
261 seg->SetEndSeed( seeds[ 1 ] );
262 for( TImage::PointType seed: seeds )
263 seg->AddSeed( seed );
270 catch( std::exception& err )
272 QMessageBox::critical(
274 QMessageBox::tr( "Error reading image" ),
275 QMessageBox::tr( err.what( ) )
280 this->m_Seeds->EnabledOff( );
282 this->m_Segmentation = seg->GetOutput( );
283 this->m_Axis = seg->GetOutputAxis( );
284 this->m_Segmentation->DisconnectPipeline( );
285 this->m_Axis->DisconnectPipeline( );
288 // -------------------------------------------------------------------------
290 _prepareQuantification( )
292 this->m_Contours.clear( );
293 this->m_ContoursActors.clear( );
294 this->m_Fourier.clear( );
296 TImage::RegionType ext_region = this->m_CPRImage->GetRequestedRegion( );
297 for( long z = 0; z < ext_region.GetSize( )[ 2 ]; ++z )
299 TImage::IndexType slc_idx = ext_region.GetIndex( );
300 TImage::SizeType slc_size = ext_region.GetSize( );
303 TImage::RegionType slc_region;
304 slc_region.SetIndex( slc_idx );
305 slc_region.SetSize( slc_size );
307 typedef itk::RegionOfInterestImageFilter< TScalarImage, TScalarImage > _TROI;
308 _TROI::Pointer roi = _TROI::New( );
309 roi->SetInput( this->m_CPRSegmentation );
310 roi->SetRegionOfInterest( slc_region );
315 slc_region.SetIndex( slc_idx );
316 slc_region.SetSize( slc_size );
317 typedef itk::Image< TBinPixel, 2 > _TBinSlice;
318 typedef itk::ExtractImageFilter< TBinImage, _TBinSlice > _TCollapse;
319 _TCollapse::Pointer collapse = _TCollapse::New( );
320 collapse->SetInput( roi->GetOutput( ) );
321 collapse->SetDirectionCollapseToIdentity( );
322 collapse->SetExtractionRegion( slc_region );
323 collapse->UpdateLargestPossibleRegion( );
325 typedef TScalarImage _TScalarSlice;
327 typedef itk::ImageToVTKImageFilter< _TScalarSlice > _TVTKScalarSlice;
328 _TVTKScalarSlice::Pointer vtk_slice = _TVTKScalarSlice::New( );
329 vtk_slice->SetInput( roi->GetOutput( ) );
330 vtk_slice->Update( );
333 vtk_slice->GetOutput( )->GetScalarRange( r );
335 vtkSmartPointer< vtkMarchingSquares > ms =
336 vtkSmartPointer< vtkMarchingSquares >::New( );
337 ms->SetInputData( vtk_slice->GetOutput( ) );
338 ms->SetValue( 0, 0 );
340 vtkSmartPointer< vtkPolyDataConnectivityFilter > conn =
341 vtkSmartPointer< vtkPolyDataConnectivityFilter >::New( );
342 conn->SetInputConnection( ms->GetOutputPort( ) );
343 conn->SetExtractionModeToClosestPointRegion( );
344 conn->SetClosestPoint( 0, 0, 0 );
347 vtkSmartPointer< vtkPolyData > cnt = vtkSmartPointer< vtkPolyData >::New( );
348 cnt->DeepCopy( conn->GetOutput( ) );
349 this->m_Contours.push_back( cnt );
351 typedef itk::QuadEdgeMesh< double, 2 > _TQEMesh;
352 typedef _TQEMesh::PointType _T2DPoint;
353 typedef std::vector< _T2DPoint > _T2DPoints;
355 _TQEMesh::Pointer mesh = _TQEMesh::New( );
356 for( unsigned int i = 0; i < cnt->GetNumberOfPoints( ); ++i )
360 cnt->GetPoint( i, p );
363 mesh->AddPoint( pnt );
366 vtkCellArray* lines = cnt->GetLines( );
367 lines->InitTraversal( );
368 vtkIdType* ids = new vtkIdType[ 2 ];
370 while( lines->GetNextCell( npts, ids ) != 0 )
371 mesh->AddEdge( ids[ 0 ], ids[ 1 ] );
373 _TQEMesh::QEPrimal* edge = mesh->GetEdge( );
376 mesh->AddFace( edge );
377 edge = mesh->GetEdge( );
378 for( auto eIt = edge->BeginGeomLnext( ); eIt != edge->EndGeomLnext( ); ++eIt )
379 points.push_back( mesh->GetPoint( *eIt ) );
383 TFourier f( points.begin( ), points.end( ), 6 );
384 f.SetOrderingToCounterClockWise( );
385 this->m_Fourier.push_back( f );
387 vtkSmartPointer< ivq::VTK::PolyDataActor > actor =
388 vtkSmartPointer< ivq::VTK::PolyDataActor >::New( );
389 actor->SetInputData( cnt );
390 actor->GetMapper( )->ScalarVisibilityOff( );
392 actor->GetProperty( )->SetColor( 1, 0, 0 );
393 actor->GetProperty( )->SetLineWidth( 2 );
394 this->m_ContoursActors.push_back( actor );
398 ivq::VTK::ImageViewer* slice_view = this->m_UI->Slice->GetViewer( );
399 vtkSmartPointer< ContourCallback > cbk =
400 vtkSmartPointer< ContourCallback >::New( );
402 slice_view->GetImageActor( )->AddObserver( vtkCommand::InteractionEvent, cbk );
403 this->m_ContourActor = NULL;
405 this->m_Spline = vtkSmartPointer< vtkSplineWidget >::New( );
406 this->m_Spline->SetCurrentRenderer( slice_view->GetRenderer( ) );
407 this->m_Spline->SetDefaultRenderer( slice_view->GetRenderer( ) );
408 this->m_Spline->SetInputData( slice_view->GetInput( ) );
409 this->m_Spline->SetProp3D( slice_view->GetImageActor( ) );
410 this->m_Spline->SetInteractor( slice_view->GetRenderWindow( )->GetInteractor( ) );
412 slice_view->GetInput( )->GetBounds( bnds );
413 this->m_Spline->PlaceWidget(
414 bnds[ 0 ], bnds[ 1 ],
415 bnds[ 2 ], bnds[ 3 ],
418 this->m_Spline->ProjectToPlaneOn( );
419 this->m_Spline->SetProjectionNormalToZAxes( );
420 this->m_Spline->SetProjectionPosition(
422 slice_view->GetImageActor( )->GetBounds( )[ 4 ] +
423 slice_view->GetImageActor( )->GetBounds( )[ 5 ]
426 this->m_Spline->SetHandleSize( 0.005 );
427 this->m_Spline->ClosedOn( );
428 this->m_Spline->SetNumberOfHandles( 16 );
429 this->m_Spline->EnabledOn( );
431 QVector< double > xdata, dmindata, dmaxdata, areadata;
433 for( unsigned int i = 0; i < this->m_Fourier.size( ); ++i )
436 this->m_Fourier[ i ].GetEllipse( 1, a, b, t, p );
437 double area = this->m_Fourier[ i ].GetArea( );
439 xdata.push_back( double( i ) );
440 dmindata.push_back( b );
441 dmaxdata.push_back( a );
442 areadata.push_back( area );
445 this->m_UI->Plot->addGraph( this->m_UI->Plot->xAxis, this->m_UI->Plot->yAxis );
446 this->m_UI->Plot->graph( 0 )->setPen( QPen( QColor( 255, 0, 0 ) ) );
447 this->m_UI->Plot->graph( 0 )->setName( "Minimum diameter (mm)" );
448 this->m_UI->Plot->graph( 0 )->setData( xdata, dmindata );
450 this->m_UI->Plot->addGraph( this->m_UI->Plot->xAxis, this->m_UI->Plot->yAxis );
451 this->m_UI->Plot->graph( 1 )->setPen( QPen( QColor( 0, 255, 0 ) ) );
452 this->m_UI->Plot->graph( 1 )->setName( "Maximum diameter (mm)" );
453 this->m_UI->Plot->graph( 1 )->setData( xdata, dmaxdata );
455 this->m_UI->Plot->addGraph( this->m_UI->Plot->xAxis, this->m_UI->Plot->yAxis2 );
456 this->m_UI->Plot->graph( 2 )->setPen( QPen( QColor( 0, 0, 255 ) ) );
457 this->m_UI->Plot->graph( 2 )->setName( "Area (mm2)" );
458 this->m_UI->Plot->graph( 2 )->setData( xdata, areadata );
460 this->m_UI->Plot->xAxis->setLabel( "Position" );
461 this->m_UI->Plot->yAxis->setLabel( "Diameters (mm)" );
462 this->m_UI->Plot->yAxis2->setLabel( "Areas (mm2)" );
464 this->m_UI->Plot->legend->setVisible( true );
465 this->m_UI->Plot->xAxis->setVisible( true );
466 this->m_UI->Plot->yAxis->setVisible( true );
467 this->m_UI->Plot->yAxis2->setVisible( true );
469 this->m_UI->Plot->rescaleAxes( );
470 this->m_UI->Plot->replot( );
473 // -------------------------------------------------------------------------
475 _showProcessResults( )
477 if( this->m_Segmentation.IsNull( ) || this->m_Axis.IsNull( ) )
480 this->m_VTKSegmentation = TVTKScalarImage::New( );
481 this->m_VTKSegmentation->SetInput( this->m_Segmentation );
482 this->m_VTKSegmentation->Update( );
485 vtkImageData* seg = this->m_VTKSegmentation->GetOutput( );
487 seg->GetScalarRange( r );
488 vtkSmartPointer< vtkMarchingCubes > mc =
489 vtkSmartPointer< vtkMarchingCubes >::New( );
490 mc->SetInputData( seg );
491 mc->SetValue( 0, 0 );
493 this->m_Surface = vtkSmartPointer< ivq::VTK::PolyDataActor >::New( );
494 this->m_Surface->SetInputConnection( mc->GetOutputPort( ) );
495 this->m_Surface->Update( );
498 this->m_Curve = TCurve::New( );
499 for( unsigned int i = 0; i < this->m_Axis->GetSize( ); ++i )
500 this->m_Curve->AddPoint( this->m_Axis->GetPoint( i ) );
501 this->m_Curve->Smooth( 2 );
505 this->m_CPRImage, this->m_VTKCPRImage,
506 this->m_Image, this->m_Curve
509 this->m_CPRSegmentation, this->m_VTKCPRSegmentation,
510 this->m_Segmentation, this->m_Curve
514 vtkRenderer* renderer =
515 this->m_UI->MPR->GetRendererWidget( )->GetRenderer( );
516 renderer->AddViewProp( this->m_Surface->GetActor( ) );
518 renderer->ResetCamera( );
519 this->m_UI->MPR->GetRendererWidget( )->GetRenderWindow( )->Render( );
521 ivq::VTK::ImageViewer* cpr_view = this->m_UI->CPR->GetViewer( );
522 ivq::VTK::ImageViewer* slice_view = this->m_UI->Slice->GetViewer( );
524 cpr_view->SetInputData( this->m_VTKCPRImage->GetOutput( ) );
525 cpr_view->SetColorWindow( 1000 );
526 cpr_view->SetColorLevel( 100 );
527 cpr_view->SetSliceOrientationToXZ( );
530 slice_view->SetInputData( this->m_VTKCPRImage->GetOutput( ) );
531 slice_view->SetColorWindow( 1000 );
532 slice_view->SetColorLevel( 100 );
533 slice_view->SetSliceOrientationToXY( );
534 slice_view->Render( );
536 this->_prepareQuantification( );
539 // -------------------------------------------------------------------------
541 _ExecuteLog( const std::string& fname )
543 std::ifstream str_log( fname.c_str( ) );
551 // -------------------------------------------------------------------------
555 FileDialog dlg( this );
557 "Medical images (*.mhd *.dcm);;Other images (*.png *.jpg *.bmp);;\
560 if( dlg.exec( ) == QDialog::Accepted )
562 std::set< std::string > filenames, dirnames;
563 QStringList names = dlg.selectedFiles( );
564 for( QString name: names )
566 QFileInfo info( name );
568 dirnames.insert( name.toStdString( ) );
570 filenames.insert( name.toStdString( ) );
574 if( dirnames.size( ) == 1 && filenames.size( ) == 0 )
575 this->_openDicom( *( dirnames.begin( ) ) );
576 else if( dirnames.size( ) == 0 && filenames.size( ) > 0 )
577 this->_openImage( filenames );
579 QMessageBox::critical(
581 "Error opening image",
582 "Directory - files mixed up, don't know what to do."
588 // -------------------------------------------------------------------------
592 this->m_DlgParameters->exec( );
595 // -------------------------------------------------------------------------
600 typedef ivq::VTK::SeedWidgetOverImageActor::TSeeds _TSeeds;
601 if( this->m_Image.IsNull( ) || this->m_Seeds.GetPointer( ) == NULL )
603 QMessageBox::critical( this, "Error processing", "No valid input image." );
607 std::vector< TImage::PointType > seeds;
608 for( _TSeeds::value_type sValue: this->m_Seeds->GetSeeds( ) )
610 for( unsigned int i = 0; i < sValue.second.size( ); i += 3 )
612 TImage::PointType pnt;
613 pnt[ 0 ] = sValue.second[ i + 0 ];
614 pnt[ 1 ] = sValue.second[ i + 1 ];
615 pnt[ 2 ] = sValue.second[ i + 2 ];
616 seeds.push_back( pnt );
621 if( seeds.size( ) < 2 )
623 QMessageBox::critical( this, "Error processing", "Not enough seeds." );
627 this->_process( seeds );
628 this->_showProcessResults( );
631 // -------------------------------------------------------------------------
632 template< class _TImagePtr, class _TVTKImagePtr, class _TCurvePtr >
635 _TImagePtr& output, _TVTKImagePtr& vtk_output,
636 const _TImagePtr& input, const _TCurvePtr& curve
639 typedef typename _TImagePtr::ObjectType _TImage;
640 typedef typename _TCurvePtr::ObjectType _TCurve;
641 typedef ivq::ITK::CPRImageFilter< _TImage, _TCurve > _TFilter;
643 typename _TFilter::Pointer cpr = _TFilter::New( );
644 cpr->SetInput( input );
645 cpr->SetInputCurve( curve );
646 cpr->SetSliceRadius( 40 );
648 output = cpr->GetOutput( );
649 output->DisconnectPipeline( );
651 vtk_output = _TVTKImagePtr::ObjectType::New( );
652 vtk_output->SetInput( output );
653 vtk_output->Update( );
656 // -------------------------------------------------------------------------
657 #include <QApplication>
659 int main( int argc, char* argv[] )
661 QApplication a( argc, argv );
662 CTArteries w( argc, argv );