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