appli/CTArteries/CTArteries.cxx

   1 // =========================================================================
   2 // @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
   3 // =========================================================================
   4
   5 #include "CTArteries.h"
   6 #include "ui_CTArteries.h"
   7 #include "ui_Parameters.h"
   8 #include "FileDialog.h"
   9
  10 #include <set>
  11
  12 #include <itkImageFileReader.h>
  13 #include <itkImageSeriesReader.h>
  14 #include <itkQuadEdgeMesh.h>
  15
  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>
  24
  25 #include <QtGui>
  26
  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>
  35
  36 #include "algorithms/RandomWalkSegmentation.h"
  37
  38 // -------------------------------------------------------------------------
  39 #define QT_ACTION_CONN( x )                                             \
  40   this->connect(                                                        \
  41     this->m_UI->a##x, SIGNAL( triggered( ) ), this, SLOT( _s##x( ) )    \
  42     )
  43
  44 // -------------------------------------------------------------------------
  45 class ContourCallback
  46   : public vtkCommand
  47 {
  48 public:
  49   static ContourCallback* New( )
  50     {
  51       return( new ContourCallback( ) );
  52     }
  53   virtual void Execute(
  54     vtkObject* caller, unsigned long eId, void* data
  55     ) override
  56     {
  57       ivq::VTK::ImageActor* actor =
  58         ivq::VTK::ImageActor::SafeDownCast( caller );
  59       if( eId == vtkCommand::InteractionEvent && actor != NULL )
  60       {
  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( ) );
  64
  65         double abnds[ 6 ];
  66         actor->GetActorBounds( abnds );
  67         unsigned int nSlice = actor->GetSliceNumber( );
  68         for( unsigned int i = 0; i < 16; ++i )
  69         {
  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 ]
  74             );
  75
  76         } // rof
  77         this->Window->m_Spline->SetProjectionPosition(
  78           (
  79             actor->GetBounds( )[ 4 ] +
  80             actor->GetBounds( )[ 5 ]
  81             ) / double( 2 )
  82           );
  83
  84
  85         this->Window->m_ContourActor =
  86           this->Window->m_ContoursActors[ actor->GetSliceNumber( ) ];
  87         viewer->GetRenderer( )->AddViewProp( this->Window->m_ContourActor->GetActor( ) );
  88         viewer->Render( );
  89
  90       } // fi
  91     }
  92
  93 protected:
  94   ContourCallback( )
  95     : vtkCommand( )
  96     {
  97     }
  98   virtual ~ContourCallback( )
  99     {
 100     }
 101
 102 private:
 103   // Purposely not implemented
 104   ContourCallback( const ContourCallback& other );
 105   ContourCallback& operator=( const ContourCallback& other );
 106 public:
 107   CTArteries* Window;
 108 };
 109
 110 // -------------------------------------------------------------------------
 111 CTArteries::
 112 CTArteries( int argc, char* argv[], QWidget* parent )
 113   : Superclass( parent ),
 114     m_UI( new Ui::CTArteries ),
 115     m_UIParameters( new Ui::Parameters )
 116 {
 117   this->m_UI->setupUi( this );
 118   this->m_DlgParameters = new QDialog( );
 119   this->m_UIParameters->setupUi( this->m_DlgParameters );
 120
 121   // Connect signals with slots
 122   QT_ACTION_CONN( Open );
 123   QT_ACTION_CONN( Config );
 124   QT_ACTION_CONN( Process );
 125   QT_ACTION_CONN( MarkReference );
 126   QT_ACTION_CONN( MarkStenosis );
 127   QT_ACTION_CONN( SaveResults );
 128   QT_ACTION_CONN( UpdateContour );
 129
 130   // Load log
 131   if( argc == 3 )
 132   {
 133     std::string cmd( argv[ 1 ] );
 134     std::string fname( argv[ 2 ] );
 135     if( cmd == "-log" )
 136       this->_ExecuteLog( fname );
 137
 138   } // fi
 139 }
 140
 141 // -------------------------------------------------------------------------
 142 CTArteries::
 143 ~CTArteries( )
 144 {
 145   delete this->m_UI;
 146   delete this->m_UIParameters;
 147   delete this->m_DlgParameters;
 148 }
 149
 150 // -------------------------------------------------------------------------
 151 template< class _TStrings >
 152 void CTArteries::
 153 _openImage( const _TStrings& fnames )
 154 {
 155   // Save log
 156   std::ofstream str_log( "CTArteries.log" );
 157   str_log << fnames.size( ) << std::endl;
 158   for( const std::string& s: fnames )
 159     str_log << s << std::endl;
 160   str_log.close( );
 161
 162   // Create appropriate reader
 163   itk::ImageSource< TImage >::Pointer reader;
 164   if( fnames.size( ) == 1 )
 165   {
 166     itk::ImageFileReader< TImage >::Pointer r =
 167       itk::ImageFileReader< TImage >::New( );
 168     r->SetFileName( *( fnames.begin( ) ) );
 169     reader = r;
 170   }
 171   else if( fnames.size( ) > 1 )
 172   {
 173     itk::ImageSeriesReader< TImage >::Pointer r =
 174       itk::ImageSeriesReader< TImage >::New( );
 175     typename _TStrings::const_iterator i;
 176     for( i = fnames.begin( ); i != fnames.end( ); ++i )
 177       r->AddFileName( *i );
 178     reader = r;
 179
 180   } // fi
 181
 182   // Execute reader
 183   if( reader.IsNull( ) )
 184     return;
 185   try
 186   {
 187     reader->Update( );
 188     this->m_Image = reader->GetOutput( );
 189     this->m_Image->DisconnectPipeline( );
 190   }
 191   catch( std::exception& err )
 192   {
 193     this->m_Image = NULL;
 194     QMessageBox::critical( this, "Error opening image", err.what( ) );
 195
 196   } // yrt
 197   this->_showInputImage( );
 198 }
 199
 200 // -------------------------------------------------------------------------
 201 void CTArteries::
 202 _openDicom( const std::string& dirname )
 203 {
 204   ivq::Qt::DicomSeriesSelectorDialog dlg( this );
 205   dlg.setStartDir( dirname );
 206   if( dlg.exec( ) == QDialog::Accepted )
 207     this->_openImage( *( dlg.selectedFilenames( ) ) );
 208 }
 209
 210 // -------------------------------------------------------------------------
 211 void CTArteries::
 212 _showInputImage( )
 213 {
 214   if( this->m_Image.IsNotNull( ) )
 215   {
 216     this->m_VTKImage = TVTKImage::New( );
 217     this->m_VTKImage->SetInput( this->m_Image );
 218     this->m_VTKImage->Update( );
 219
 220     ivq::VTK::MPRViewers* viewers = this->m_UI->MPR->GetViewers( );
 221     viewers->SetInputData( this->m_VTKImage->GetOutput( ) );
 222     viewers->SetColorWindow( 1000 );
 223     viewers->SetColorLevel( 100 );
 224
 225     ivq::VTK::ImageViewer* view = viewers->GetView( 2 );
 226     this->m_Seeds =
 227       vtkSmartPointer< ivq::VTK::SeedWidgetOverImageActor >::New( );
 228     this->m_Seeds->SetActor( view->GetImageActor( ) );
 229     this->m_Seeds->SetInteractor( view->GetRenderWindow( )->GetInteractor( ) );
 230
 231     viewers->Render( );
 232     viewers->ResetCameras( );
 233     this->m_Seeds->EnabledOn( );
 234     viewers->Render( );
 235   }
 236   else
 237   {
 238     this->m_VTKImage = NULL;
 239
 240   } // fi
 241 }
 242
 243 // -------------------------------------------------------------------------
 244 void CTArteries::
 245 _process( const std::vector< TImage::PointType >& seeds )
 246 {
 247   // Save log
 248   std::ofstream str_log( "CTArteries.log", std::ios_base::app );
 249   str_log << this->m_UIParameters->Beta->value( ) << std::endl;
 250   str_log << this->m_UIParameters->Sigma->value( ) << std::endl;
 251   str_log << this->m_UIParameters->Radius->value( ) << std::endl;
 252   str_log << seeds.size( ) << std::endl;
 253   for( TImage::PointType seed: seeds )
 254     str_log << seed << std::endl;
 255   str_log.close( );
 256
 257   // Create algorithm
 258   typedef RandomWalkSegmentation< TImage, TScalarImage > _TSegmentation;
 259   _TSegmentation::Pointer seg = _TSegmentation::New( );
 260   seg->SetInput( this->m_Image );
 261   seg->SetBeta( this->m_UIParameters->Beta->value( ) );
 262   seg->SetSigma( this->m_UIParameters->Sigma->value( ) );
 263   seg->SetRadius( this->m_UIParameters->Radius->value( ) );
 264   seg->SetStartSeed( seeds[ 0 ] );
 265   seg->SetEndSeed( seeds[ 1 ] );
 266   for( TImage::PointType seed: seeds )
 267     seg->AddSeed( seed );
 268   try
 269   {
 270     seg->DebugOn( );
 271     seg->Update( );
 272     seg->DebugOff( );
 273   }
 274   catch( std::exception& err )
 275   {
 276     QMessageBox::critical(
 277       this,
 278       QMessageBox::tr( "Error reading image" ),
 279       QMessageBox::tr( err.what( ) )
 280       );
 281     return;
 282
 283   } // yrt
 284   this->m_Seeds->EnabledOff( );
 285
 286   this->m_Segmentation = seg->GetOutput( );
 287   this->m_Axis = seg->GetOutputAxis( );
 288   this->m_Segmentation->DisconnectPipeline( );
 289   this->m_Axis->DisconnectPipeline( );
 290 }
 291
 292 // -------------------------------------------------------------------------
 293 void CTArteries::
 294 _prepareQuantification( )
 295 {
 296   this->m_Contours.clear( );
 297   this->m_ContoursActors.clear( );
 298   this->m_Fourier.clear( );
 299
 300   TImage::RegionType ext_region = this->m_CPRImage->GetRequestedRegion( );
 301   for( long z = 0; z < ext_region.GetSize( )[ 2 ]; ++z )
 302   {
 303     TImage::IndexType slc_idx = ext_region.GetIndex( );
 304     TImage::SizeType slc_size = ext_region.GetSize( );
 305     slc_idx[ 2 ] += z;
 306     slc_size[ 2 ] = 1;
 307     TImage::RegionType slc_region;
 308     slc_region.SetIndex( slc_idx );
 309     slc_region.SetSize( slc_size );
 310
 311     typedef itk::RegionOfInterestImageFilter< TScalarImage, TScalarImage > _TROI;
 312     _TROI::Pointer roi = _TROI::New( );
 313     roi->SetInput( this->m_CPRSegmentation );
 314     roi->SetRegionOfInterest( slc_region );
 315
 316     /* TODO
 317        slc_idx[ 2 ] = 1;
 318        slc_size[ 2 ] = 0;
 319        slc_region.SetIndex( slc_idx );
 320        slc_region.SetSize( slc_size );
 321        typedef itk::Image< TBinPixel, 2 > _TBinSlice;
 322        typedef itk::ExtractImageFilter< TBinImage, _TBinSlice > _TCollapse;
 323        _TCollapse::Pointer collapse = _TCollapse::New( );
 324        collapse->SetInput( roi->GetOutput( ) );
 325        collapse->SetDirectionCollapseToIdentity( );
 326        collapse->SetExtractionRegion( slc_region );
 327        collapse->UpdateLargestPossibleRegion( );
 328     */
 329     typedef TScalarImage _TScalarSlice;
 330
 331     typedef itk::ImageToVTKImageFilter< _TScalarSlice > _TVTKScalarSlice;
 332     _TVTKScalarSlice::Pointer vtk_slice = _TVTKScalarSlice::New( );
 333     vtk_slice->SetInput( roi->GetOutput( ) );
 334     vtk_slice->Update( );
 335
 336     double r[ 2 ];
 337     vtk_slice->GetOutput( )->GetScalarRange( r );
 338
 339     vtkSmartPointer< vtkMarchingSquares > ms =
 340       vtkSmartPointer< vtkMarchingSquares >::New( );
 341     ms->SetInputData( vtk_slice->GetOutput( ) );
 342     ms->SetValue( 0, 0 );
 343
 344     vtkSmartPointer< vtkPolyDataConnectivityFilter > conn =
 345       vtkSmartPointer< vtkPolyDataConnectivityFilter >::New( );
 346     conn->SetInputConnection( ms->GetOutputPort( ) );
 347     conn->SetExtractionModeToClosestPointRegion( );
 348     conn->SetClosestPoint( 0, 0, 0 );
 349     conn->Update( );
 350
 351     vtkSmartPointer< vtkPolyData > cnt = vtkSmartPointer< vtkPolyData >::New( );
 352     cnt->DeepCopy( conn->GetOutput( ) );
 353     this->m_Contours.push_back( cnt );
 354
 355     typedef itk::QuadEdgeMesh< double, 2 > _TQEMesh;
 356     typedef _TQEMesh::PointType _T2DPoint;
 357     typedef std::vector< _T2DPoint > _T2DPoints;
 358     _T2DPoints points;
 359     _TQEMesh::Pointer mesh = _TQEMesh::New( );
 360     for( unsigned int i = 0; i < cnt->GetNumberOfPoints( ); ++i )
 361     {
 362       _T2DPoint pnt;
 363       double p[ 3 ];
 364       cnt->GetPoint( i, p );
 365       pnt[ 0 ] = p[ 0 ];
 366       pnt[ 1 ] = p[ 1 ];
 367       mesh->AddPoint( pnt );
 368
 369     } // rof
 370     vtkCellArray* lines = cnt->GetLines( );
 371     lines->InitTraversal( );
 372     vtkIdType* ids = new vtkIdType[ 2 ];
 373     vtkIdType npts;
 374     while( lines->GetNextCell( npts, ids ) != 0 )
 375       mesh->AddEdge( ids[ 0 ], ids[ 1 ] );
 376     delete ids;
 377     _TQEMesh::QEPrimal* edge = mesh->GetEdge( );
 378     if( edge != NULL )
 379     {
 380       mesh->AddFace( edge );
 381       edge = mesh->GetEdge( );
 382       for( auto eIt = edge->BeginGeomLnext( ); eIt != edge->EndGeomLnext( ); ++eIt )
 383         points.push_back( mesh->GetPoint( *eIt ) );
 384
 385       TFourier f( points.begin( ), points.end( ), 6 );
 386       f.SetOrderingToCounterClockWise( );
 387       this->m_Fourier.push_back( f );
 388     }
 389     else
 390     {
 391       TFourier f( 6 );
 392       f[ 1 ] = 0.1;
 393       this->m_Fourier.push_back( f );
 394
 395     } // fi
 396
 397     vtkSmartPointer< ivq::VTK::PolyDataActor > actor =
 398       vtkSmartPointer< ivq::VTK::PolyDataActor >::New( );
 399     actor->SetInputData( cnt );
 400     actor->GetMapper( )->ScalarVisibilityOff( );
 401     actor->Update( );
 402     actor->GetProperty( )->SetColor( 1, 0, 0 );
 403     actor->GetProperty( )->SetLineWidth( 2 );
 404     this->m_ContoursActors.push_back( actor );
 405
 406   } // rof
 407
 408   ivq::VTK::ImageViewer* slice_view = this->m_UI->Slice->GetViewer( );
 409   vtkSmartPointer< ContourCallback > cbk =
 410     vtkSmartPointer< ContourCallback >::New( );
 411   cbk->Window = this;
 412   slice_view->GetImageActor( )->AddObserver( vtkCommand::InteractionEvent, cbk );
 413   this->m_ContourActor = NULL;
 414
 415   this->m_Spline = vtkSmartPointer< vtkSplineWidget >::New( );
 416   this->m_Spline->SetCurrentRenderer( slice_view->GetRenderer( ) );
 417   this->m_Spline->SetDefaultRenderer( slice_view->GetRenderer( ) );
 418   this->m_Spline->SetInputData( slice_view->GetInput( ) );
 419   this->m_Spline->SetProp3D( slice_view->GetImageActor( ) );
 420   this->m_Spline->SetInteractor( slice_view->GetRenderWindow( )->GetInteractor( ) );
 421   double bnds[ 6 ];
 422   slice_view->GetInput( )->GetBounds( bnds );
 423   this->m_Spline->PlaceWidget(
 424     bnds[ 0 ], bnds[ 1 ],
 425     bnds[ 2 ], bnds[ 3 ],
 426     bnds[ 4 ], bnds[ 5 ]
 427     );
 428   this->m_Spline->ProjectToPlaneOn( );
 429   this->m_Spline->SetProjectionNormalToZAxes( );
 430   this->m_Spline->SetProjectionPosition(
 431     (
 432       slice_view->GetImageActor( )->GetBounds( )[ 4 ] +
 433       slice_view->GetImageActor( )->GetBounds( )[ 5 ]
 434       ) / double( 2 )
 435     );
 436   this->m_Spline->SetHandleSize( 0.005 );
 437   this->m_Spline->ClosedOn( );
 438   this->m_Spline->SetNumberOfHandles( 16 );
 439   this->m_Spline->EnabledOn( );
 440
 441   QVector< double > xdata, dmindata, dmaxdata, areadata;
 442   double data[ 3 ];
 443   for( unsigned int i = 0; i < this->m_Fourier.size( ); ++i )
 444   {
 445     double a, b, t, p;
 446     this->m_Fourier[ i ].GetEllipse( 1, a, b, t, p );
 447     double area = this->m_Fourier[ i ].GetArea( );
 448
 449     xdata.push_back( double( i ) );
 450     dmindata.push_back( b );
 451     dmaxdata.push_back( a );
 452     areadata.push_back( area );
 453
 454   } // rof
 455   this->m_UI->Plot->addGraph( this->m_UI->Plot->xAxis, this->m_UI->Plot->yAxis );
 456   this->m_UI->Plot->graph( 0 )->setPen( QPen( QColor( 255, 0, 0 ) ) );
 457   this->m_UI->Plot->graph( 0 )->setName( "Minimum diameter (mm)" );
 458   this->m_UI->Plot->graph( 0 )->setData( xdata, dmindata );
 459
 460   this->m_UI->Plot->addGraph( this->m_UI->Plot->xAxis, this->m_UI->Plot->yAxis );
 461   this->m_UI->Plot->graph( 1 )->setPen( QPen( QColor( 0, 255, 0 ) ) );
 462   this->m_UI->Plot->graph( 1 )->setName( "Maximum diameter (mm)" );
 463   this->m_UI->Plot->graph( 1 )->setData( xdata, dmaxdata );
 464
 465   this->m_UI->Plot->addGraph( this->m_UI->Plot->xAxis, this->m_UI->Plot->yAxis2 );
 466   this->m_UI->Plot->graph( 2 )->setPen( QPen( QColor( 0, 0, 255 ) ) );
 467   this->m_UI->Plot->graph( 2 )->setName( "Area (mm2)" );
 468   this->m_UI->Plot->graph( 2 )->setData( xdata, areadata );
 469
 470   this->m_UI->Plot->xAxis->setLabel( "Position" );
 471   this->m_UI->Plot->yAxis->setLabel( "Diameters (mm)" );
 472   this->m_UI->Plot->yAxis2->setLabel( "Areas (mm2)" );
 473
 474   this->m_UI->Plot->legend->setVisible( true );
 475   this->m_UI->Plot->xAxis->setVisible( true );
 476   this->m_UI->Plot->yAxis->setVisible( true );
 477   this->m_UI->Plot->yAxis2->setVisible( true );
 478
 479   this->m_UI->Plot->rescaleAxes( );
 480   this->m_UI->Plot->replot( );
 481   this->m_StenosisText = new QCPItemText( this->m_UI->Plot );
 482   this->m_StenosisText->setText( "Stenosis: NaN" );
 483   this->m_UI->Plot->replot( );
 484 }
 485
 486 // -------------------------------------------------------------------------
 487 void CTArteries::
 488 _showProcessResults( )
 489 {
 490   if( this->m_Segmentation.IsNull( ) || this->m_Axis.IsNull( ) )
 491     return;
 492
 493   this->m_VTKSegmentation = TVTKScalarImage::New( );
 494   this->m_VTKSegmentation->SetInput( this->m_Segmentation );
 495   this->m_VTKSegmentation->Update( );
 496
 497   // Show surface
 498   vtkImageData* seg = this->m_VTKSegmentation->GetOutput( );
 499   double r[ 2 ];
 500   seg->GetScalarRange( r );
 501   vtkSmartPointer< vtkMarchingCubes > mc =
 502     vtkSmartPointer< vtkMarchingCubes >::New( );
 503   mc->SetInputData( seg );
 504   mc->SetValue( 0, 0 );
 505   mc->Update( );
 506   this->m_Surface = vtkSmartPointer< ivq::VTK::PolyDataActor >::New( );
 507   this->m_Surface->SetInputConnection( mc->GetOutputPort( ) );
 508   this->m_Surface->Update( );
 509
 510   // Prepare curve
 511   this->m_Curve = TCurve::New( );
 512   for( unsigned int i = 0; i < this->m_Axis->GetSize( ); ++i )
 513     this->m_Curve->AddPoint( this->m_Axis->GetPoint( i ) );
 514   this->m_Curve->Smooth( 2 );
 515
 516   // Compute CPRs
 517   this->_CPR(
 518     this->m_CPRImage, this->m_VTKCPRImage,
 519     this->m_Image, this->m_Curve
 520     );
 521   this->_CPR(
 522     this->m_CPRSegmentation, this->m_VTKCPRSegmentation,
 523     this->m_Segmentation, this->m_Curve
 524     );
 525
 526   // Show all data
 527   vtkRenderer* renderer =
 528     this->m_UI->MPR->GetRendererWidget( )->GetRenderer( );
 529   renderer->AddViewProp( this->m_Surface->GetActor( ) );
 530   renderer->Render( );
 531   renderer->ResetCamera( );
 532   this->m_UI->MPR->GetRendererWidget( )->GetRenderWindow( )->Render( );
 533
 534   ivq::VTK::ImageViewer* cpr_view = this->m_UI->CPR->GetViewer( );
 535   ivq::VTK::ImageViewer* slice_view = this->m_UI->Slice->GetViewer( );
 536
 537   cpr_view->SetInputData( this->m_VTKCPRImage->GetOutput( ) );
 538   cpr_view->SetColorWindow( 1000 );
 539   cpr_view->SetColorLevel( 100 );
 540   cpr_view->SetSliceOrientationToXZ( );
 541   cpr_view->Render( );
 542
 543   slice_view->SetInputData( this->m_VTKCPRImage->GetOutput( ) );
 544   slice_view->SetColorWindow( 1000 );
 545   slice_view->SetColorLevel( 100 );
 546   slice_view->SetSliceOrientationToXY( );
 547   slice_view->Render( );
 548
 549   this->_prepareQuantification( );
 550   this->m_StenosisSlice = -1;
 551   this->m_ReferenceSlice = -1;
 552   this->m_StenosisLine = NULL;
 553   this->m_ReferenceLine = NULL;
 554 }
 555
 556 // -------------------------------------------------------------------------
 557 void CTArteries::
 558 _showStenosis( )
 559 {
 560   if( this->m_StenosisSlice >= 0 )
 561   {
 562     if( this->m_StenosisLine == NULL )
 563       this->m_StenosisLine = new QCPItemLine( this->m_UI->Plot );
 564     this->m_StenosisLine->start->setCoords( this->m_StenosisSlice, 0 );
 565     this->m_StenosisLine->end->setCoords(
 566       this->m_StenosisSlice,
 567       this->m_UI->Plot->yAxis->range( ).size( )
 568       );
 569     this->m_UI->Plot->replot( );
 570
 571   } // fi
 572   if( this->m_ReferenceSlice >= 0 )
 573   {
 574     if( this->m_ReferenceLine == NULL )
 575       this->m_ReferenceLine = new QCPItemLine( this->m_UI->Plot );
 576     this->m_ReferenceLine->start->setCoords( this->m_ReferenceSlice, 0 );
 577     this->m_ReferenceLine->end->setCoords(
 578       this->m_ReferenceSlice,
 579       this->m_UI->Plot->yAxis->range( ).size( )
 580       );
 581     this->m_UI->Plot->replot( );
 582
 583   } // fi
 584   if( this->m_StenosisSlice < 0 || this->m_ReferenceSlice < 0 )
 585     return;
 586
 587   double sa, sb, st, sp;
 588   this->m_Fourier[ this->m_StenosisSlice ].GetEllipse( 1, sa, sb, st, sp );
 589   double sarea = this->m_Fourier[ this->m_StenosisSlice ].GetArea( );
 590   double ra, rb, rt, rp;
 591   this->m_Fourier[ this->m_ReferenceSlice ].GetEllipse( 1, ra, rb, rt, rp );
 592   double rarea = this->m_Fourier[ this->m_ReferenceSlice ].GetArea( );
 593   double qmin = 100.0 * ( 1.0 - ( sb / rb ) );
 594   double qmax = 100.0 * ( 1.0 - ( sa / ra ) );
 595   double qarea = 100.0 * ( 1.0 - ( sarea / rarea ) );
 596
 597   std::stringstream s;
 598   s << "Stenosis: " << std::endl
 599     << qmin << "% (min diameter)" << std::endl
 600     << qmax << "% (max diameter)" << std::endl
 601     << qarea << "% (area)";
 602   this->m_StenosisText->setText( s.str( ).c_str( ) );
 603   this->m_StenosisText->position->setCoords(
 604     ( this->m_ReferenceSlice + this->m_StenosisSlice ) * 0.5,
 605     this->m_UI->Plot->yAxis->range( ).center( )
 606     );
 607   this->m_UI->Plot->replot( );
 608 }
 609
 610 // -------------------------------------------------------------------------
 611 void CTArteries::
 612 _ExecuteLog( const std::string& fname )
 613 {
 614   std::ifstream str_log( fname.c_str( ) );
 615   if( str_log )
 616   {
 617     str_log.close( );
 618
 619   } // fi
 620 }
 621
 622 // -------------------------------------------------------------------------
 623 void CTArteries::
 624 _sOpen( )
 625 {
 626   FileDialog dlg( this );
 627   dlg.setNameFilter(
 628     "Medical images (*.mhd *.dcm);;Other images (*.png *.jpg *.bmp);;\
 629 All files (*)"
 630     );
 631   if( dlg.exec( ) == QDialog::Accepted )
 632   {
 633     std::set< std::string > filenames, dirnames;
 634     QStringList names = dlg.selectedFiles( );
 635     for( QString name: names )
 636     {
 637       QFileInfo info( name );
 638       if( info.isDir( ) )
 639         dirnames.insert( name.toStdString( ) );
 640       else
 641         filenames.insert( name.toStdString( ) );
 642
 643     } // rof
 644
 645     if( dirnames.size( ) == 1 && filenames.size( ) == 0 )
 646       this->_openDicom( *( dirnames.begin( ) ) );
 647     else if( dirnames.size( ) == 0 && filenames.size( ) > 0 )
 648       this->_openImage( filenames );
 649     else
 650       QMessageBox::critical(
 651         this,
 652         "Error opening image",
 653         "Directory - files mixed up, don't know what to do."
 654         );
 655
 656   } // fi
 657 }
 658
 659 // -------------------------------------------------------------------------
 660 void CTArteries::
 661 _sConfig( )
 662 {
 663   this->m_DlgParameters->exec( );
 664 }
 665
 666 // -------------------------------------------------------------------------
 667 void CTArteries::
 668 _sProcess( )
 669 {
 670   // Get seeds
 671   typedef ivq::VTK::SeedWidgetOverImageActor::TSeeds _TSeeds;
 672   if( this->m_Image.IsNull( ) || this->m_Seeds.GetPointer( ) == NULL )
 673   {
 674     QMessageBox::critical( this, "Error processing", "No valid input image." );
 675     return;
 676
 677   } // fi
 678   std::vector< TImage::PointType > seeds;
 679   for( _TSeeds::value_type sValue: this->m_Seeds->GetSeeds( ) )
 680   {
 681     for( unsigned int i = 0; i < sValue.second.size( ); i += 3 )
 682     {
 683       TImage::PointType pnt;
 684       pnt[ 0 ] = sValue.second[ i + 0 ];
 685       pnt[ 1 ] = sValue.second[ i + 1 ];
 686       pnt[ 2 ] = sValue.second[ i + 2 ];
 687       seeds.push_back( pnt );
 688
 689     } // rof
 690
 691   } // rof
 692   if( seeds.size( ) < 2 )
 693   {
 694     QMessageBox::critical( this, "Error processing", "Not enough seeds." );
 695     return;
 696
 697   } // fi
 698   this->_process( seeds );
 699   this->_showProcessResults( );
 700 }
 701
 702 // -------------------------------------------------------------------------
 703 void CTArteries::
 704 _sMarkReference( )
 705 {
 706   if( this->m_Fourier.size( ) == 0 )
 707     return;
 708   ivq::VTK::ImageViewer* slice_view = this->m_UI->Slice->GetViewer( );
 709   this->m_ReferenceSlice = slice_view->GetImageActor( )->GetSliceNumber( );
 710   this->_showStenosis( );
 711 }
 712
 713 // -------------------------------------------------------------------------
 714 void CTArteries::
 715 _sMarkStenosis( )
 716 {
 717   if( this->m_Fourier.size( ) == 0 )
 718     return;
 719   ivq::VTK::ImageViewer* slice_view = this->m_UI->Slice->GetViewer( );
 720   this->m_StenosisSlice = slice_view->GetImageActor( )->GetSliceNumber( );
 721   this->_showStenosis( );
 722 }
 723
 724 // -------------------------------------------------------------------------
 725 void CTArteries::
 726 _sSaveResults( )
 727 {
 728   std::string fname =
 729     QFileDialog::getSaveFileName(
 730       this,
 731       "Select one or more files to open",
 732       ".",
 733       "CSV files (*.csv);;All files (*)"
 734       ).toStdString( );
 735   if( fname == "" )
 736     return;
 737
 738   std::ofstream out( fname.c_str( ) );
 739   out << "slice;min_diameter;max_diameter;area" << std::endl;
 740   for( unsigned int i = 0; i < this->m_Fourier.size( ); ++i )
 741   {
 742     double a, b, t, p;
 743     this->m_Fourier[ i ].GetEllipse( 1, a, b, t, p );
 744     double area = this->m_Fourier[ i ].GetArea( );
 745     out << i << ";" << b << ";" << a << ";" << area << std::endl;
 746
 747   } // rof
 748   out << "############################" << std::endl;
 749   typedef ivq::VTK::SeedWidgetOverImageActor::TSeeds _TSeeds;
 750   const _TSeeds& seeds = this->m_Seeds->GetSeeds( );
 751   for( _TSeeds::value_type seed: seeds )
 752   {
 753     for( unsigned int i = 0; i < seed.second.size( ); i += 3 )
 754     {
 755       out
 756         << seed.second[ i + 0 ] << " "
 757         << seed.second[ i + 1 ] << " "
 758         << seed.second[ i + 2 ] << std::endl;
 759
 760     } // rof
 761
 762   } // rof
 763   out << "############################" << std::endl;
 764   out.close( );
 765 }
 766
 767 // -------------------------------------------------------------------------
 768 void CTArteries::
 769 _sUpdateContour( )
 770 {
 771   if( this->m_ContourActor == NULL )
 772     return;
 773   vtkSmartPointer< vtkPolyData > new_pd =
 774     vtkSmartPointer< vtkPolyData >::New( );
 775   this->m_Spline->GetPolyData( new_pd );
 776
 777   ivq::VTK::ImageViewer* slice_view = this->m_UI->Slice->GetViewer( );
 778   ivq::VTK::ImageActor* slice_actor = slice_view->GetImageActor( );
 779   this->m_ContoursActors[ slice_actor->GetSliceNumber( ) ]->
 780     SetInputData( new_pd );
 781   slice_view->Render( );
 782
 783
 784   typedef itk::QuadEdgeMesh< double, 2 > _TQEMesh;
 785   typedef _TQEMesh::PointType _T2DPoint;
 786   typedef std::vector< _T2DPoint > _T2DPoints;
 787   _T2DPoints points;
 788   _TQEMesh::Pointer mesh = _TQEMesh::New( );
 789   for( unsigned int i = 0; i < new_pd->GetNumberOfPoints( ); ++i )
 790   {
 791     _T2DPoint pnt;
 792     double p[ 3 ];
 793     new_pd->GetPoint( i, p );
 794     pnt[ 0 ] = p[ 0 ];
 795     pnt[ 1 ] = p[ 1 ];
 796     mesh->AddPoint( pnt );
 797
 798   } // rof
 799   vtkCellArray* lines = new_pd->GetLines( );
 800   lines->InitTraversal( );
 801   vtkIdType* ids = new vtkIdType[ 2 ];
 802   vtkIdType npts;
 803   while( lines->GetNextCell( npts, ids ) != 0 )
 804     mesh->AddEdge( ids[ 0 ], ids[ 1 ] );
 805   delete ids;
 806   _TQEMesh::QEPrimal* edge = mesh->GetEdge( );
 807   if( edge != NULL )
 808   {
 809     mesh->AddFace( edge );
 810     edge = mesh->GetEdge( );
 811     for( auto eIt = edge->BeginGeomLnext( ); eIt != edge->EndGeomLnext( ); ++eIt )
 812       points.push_back( mesh->GetPoint( *eIt ) );
 813
 814     TFourier f( points.begin( ), points.end( ), 6 );
 815     f.SetOrderingToCounterClockWise( );
 816     this->m_Fourier[ slice_actor->GetSliceNumber( ) ] = f;
 817
 818   } // fi
 819 }
 820
 821 // -------------------------------------------------------------------------
 822 template< class _TImagePtr, class _TVTKImagePtr, class _TCurvePtr >
 823 void CTArteries::
 824 _CPR(
 825   _TImagePtr& output, _TVTKImagePtr& vtk_output,
 826   const _TImagePtr& input, const _TCurvePtr& curve
 827   )
 828 {
 829   typedef typename _TImagePtr::ObjectType _TImage;
 830   typedef typename _TCurvePtr::ObjectType _TCurve;
 831   typedef ivq::ITK::CPRImageFilter< _TImage, _TCurve > _TFilter;
 832
 833   typename _TFilter::Pointer cpr = _TFilter::New( );
 834   cpr->SetInput( input );
 835   cpr->SetInputCurve( curve );
 836   cpr->SetSliceRadius( 40 );
 837   cpr->Update( );
 838   output = cpr->GetOutput( );
 839   output->DisconnectPipeline( );
 840
 841   vtk_output = _TVTKImagePtr::ObjectType::New( );
 842   vtk_output->SetInput( output );
 843   vtk_output->Update( );
 844 }
 845
 846 // -------------------------------------------------------------------------
 847 #include <QApplication>
 848
 849 int main( int argc, char* argv[] )
 850 {
 851   QApplication a( argc, argv );
 852   CTArteries w( argc, argv );
 853   w.show( );
 854   return( a.exec( ) );
 855 }
 856
 857 // eof - $RCSfile$