CustusX  18.04
An IGT application
cxAirwaysFilterService.cpp
Go to the documentation of this file.
1 /*=========================================================================
2 This file is part of CustusX, an Image Guided Therapy Application.
3 
4 Copyright (c) SINTEF Department of Medical Technology.
5 All rights reserved.
6 
7 CustusX is released under a BSD 3-Clause license.
8 
9 See Lisence.txt (https://github.com/SINTEFMedtek/CustusX/blob/master/License.txt) for details.
10 =========================================================================*/
11 
12 #include "cxAirwaysFilterService.h"
13 
14 #include <QTimer>
15 
16 #include <vtkImageImport.h>
17 #include <vtkImageData.h>
18 #include <vtkImageShiftScale.h>
19 #include <ctkPluginContext.h>
20 #include <vtkImplicitModeller.h>
21 #include <vtkContourFilter.h>
22 #include "cxBranchList.h"
24 
25 #include "cxTime.h"
26 #include "cxTypeConversions.h"
27 #include "cxLogger.h"
28 #include "cxDataReaderWriter.h"
30 #include "cxDoubleProperty.h"
31 #include "cxContourFilter.h"
32 #include "cxDataLocations.h"
34 #include "vtkForwardDeclarations.h"
36 #include "cxVisServices.h"
37 #include "cxUtilHelpers.h"
38 #include "FAST/Algorithms/LungSegmentation/LungSegmentation.hpp"
39 #include "FAST/Algorithms/AirwaySegmentation/AirwaySegmentation.hpp"
40 #include "FAST/Algorithms/CenterlineExtraction/CenterlineExtraction.hpp"
41 #include "FAST/Importers/ImageFileImporter.hpp"
42 #include "FAST/Exporters/VTKImageExporter.hpp"
43 #include "FAST/Exporters/VTKMeshExporter.hpp"
44 #include "FAST/Data/Segmentation.hpp"
45 #include "FAST/SceneGraph.hpp"
46 
47 namespace cx {
48 
50  FilterImpl(services),
51  mDefaultStraightCLTubesOption(false)
52 {
53  fast::Reporter::setGlobalReportMethod(fast::Reporter::COUT);
54  //Need to create OpenGL context of fast in main thread, this is done in the constructor of DeviceManger
55  fast::ImageFileImporter::pointer importer = fast::ImageFileImporter::New();
56  Q_UNUSED(importer)
57 }
58 
59 
61 }
62 
63 QString AirwaysFilter::getName() const
64 {
65  return "Airway Segmentation Filter";
66 }
67 
68 QString AirwaysFilter::getType() const
69 {
70  return "airways_filter";
71 }
72 
73 QString AirwaysFilter::getHelp() const
74 {
75  return "<html>"
76  "<h3>Airway Segmentation.</h3>"
77  "<p><i>Extracts segmentation and centerline from a CT volume. If method fails, try to crop volume. </br>Algorithm written by Erik Smistad.</i></p>"
78  "</html>";
79 }
80 
82 {
83  return "_centerline";
84 }
85 
87 {
88  return "_straight";
89 }
90 
92 {
93  return "_tubes";
94 }
95 
96 Vector3D AirwaysFilter::getSeedPointFromTool(SpaceProviderPtr spaceProvider, DataPtr data)
97 {
98  // Retrieve position of tooltip and use it as seed point
99  Vector3D point = spaceProvider->getActiveToolTipPoint(
100  spaceProvider->getD(data));
101 
102  // Have to multiply by the inverse of the spacing to get the voxel position
103  ImagePtr image = boost::dynamic_pointer_cast<Image>(data);
104  double spacingX, spacingY, spacingZ;
105  image->getBaseVtkImageData()->GetSpacing(spacingX, spacingY, spacingZ);
106  point(0) = point(0) * (1.0 / spacingX);
107  point(1) = point(1) * (1.0 / spacingY);
108  point(2) = point(2) * (1.0 / spacingZ);
109 
110  std::cout << "the selected seed point is: " << point(0) << " " << point(1)
111  << " " << point(2) << "\n";
112 
113  return point;
114 }
115 
116 int * getImageSize(DataPtr inputImage)
117 {
118  ImagePtr image = boost::dynamic_pointer_cast<Image>(inputImage);
119  return image->getBaseVtkImageData()->GetDimensions();
120 }
121 
122 bool AirwaysFilter::isSeedPointInsideImage(Vector3D seedPoint, DataPtr image)
123 {
124  int * size = getImageSize(image);
125  std::cout << "size of image is: " << size[0] << " " << size[1] << " "
126  << size[2] << "\n";
127  int x = (int) seedPoint(0);
128  int y = (int) seedPoint(1);
129  int z = (int) seedPoint(2);
130  bool result = x >= 0 && y >= 0 && z >= 0 && x < size[0] && y < size[1]
131  && z < size[2];
132  return result;
133 }
134 
136 {
137  DataPtr inputImage = mInputTypes[0].get()->getData();
138  if (!inputImage)
139  {
140  CX_LOG_ERROR() << "No input data selected";
141  return false;
142  }
143 
144  if (inputImage->getType() != "image")
145  {
146  CX_LOG_ERROR() << "Input data has to be an image";
147  return false;
148  }
149 
150  std::string filename = (patientService()->getActivePatientFolder()
151  + "/" + inputImage->getFilename()).toStdString();
152 
153  // only check seed point inside image if use seed point is checked
154  bool useManualSeedPoint = getManualSeedPointOption(mOptions)->getValue();
155  if(useManualSeedPoint)
156  {
157  seedPoint = getSeedPointFromTool(mServices->spaceProvider(), inputImage);
158  if(!isSeedPointInsideImage(seedPoint, inputImage)) {
159  CX_LOG_ERROR() << "Seed point is not inside image. Use cursor to set seed point inside trachea in the CT image.";
160  return false;
161  }
162  }
163  mInputImage = patientService()->getData<Image>(inputImage->getUid());
164 
165  return true;
166 }
167 
169 {
170  CX_LOG_INFO() << "EXECUTING AIRWAYS FILTER";
171  // Check if pre process went ok:
172  if(!mInputImage)
173  return false;
174 
175  QString q_filename = "";
176  QString activePatienFolder = patientService()->getActivePatientFolder();
177  QString inputImageFileName = mInputImage->getFilename();
178  if(!activePatienFolder.isEmpty())
179  q_filename = activePatienFolder+"/"+inputImageFileName;
180  else
181  q_filename = inputImageFileName;
182 
183  std::string filename = q_filename.toStdString();
184  try {
185  fast::Config::getTestDataPath(); // needed for initialization
186  QString cacheDir = cx::DataLocations::getCachePath();
187  fast::Config::setKernelBinaryPath(cacheDir.toStdString());
188  QString kernelDir = cx::DataLocations::findConfigFolder("/FAST", FAST_SOURCE_DIR);
189  fast::Config::setKernelSourcePath(kernelDir.toStdString());
190 
191  // Import image data from disk
192  fast::ImageFileImporter::pointer importer = fast::ImageFileImporter::New();
193  importer->setFilename(filename);
194 
195  // Need to know the data type
196  importer->update();
197  fast::Image::pointer image = importer->getOutputData<fast::Image>();
198 
199  // Do segmentation
200  fast::Segmentation::pointer segmentationData;
201  bool doLungSegmentation = getLungSegmentationOption(mOptions)->getValue();
202  bool useManualSeedPoint = getManualSeedPointOption(mOptions)->getValue();
203  try {
204  if(doLungSegmentation) {
205  fast::LungSegmentation::pointer segmentation = fast::LungSegmentation::New();
206  if(useManualSeedPoint) {
207  CX_LOG_INFO() << "Using seed point: " << seedPoint.transpose();
208  segmentation->setAirwaySeedPoint(seedPoint(0), seedPoint(1), seedPoint(2));
209  }
210  segmentation->setInputConnection(importer->getOutputPort());
211  segmentation->update();
212  segmentationData = segmentation->getOutputData<fast::Segmentation>(1);
213 
214  // Convert fast segmentation data to VTK data which CX can use (Airways)
215  vtkSmartPointer<fast::VTKImageExporter> vtkExporter = fast::VTKImageExporter::New();
216  vtkExporter->setInputConnection(segmentation->getOutputPort(1));
217  vtkExporter->Update();
218  mAirwaySegmentationOutput = vtkExporter->GetOutput();
219 
220  // Convert fast segmentation data to VTK data which CX can use (Lungs)
221  vtkSmartPointer<fast::VTKImageExporter> vtkExporter2 = fast::VTKImageExporter::New();
222  vtkExporter2->setInputConnection(segmentation->getOutputPort(0));
223  vtkExporter2->Update();
224  mLungSegmentationOutput = vtkExporter2->GetOutput();
225  } else {
226 
227  fast::AirwaySegmentation::pointer segmentation = fast::AirwaySegmentation::New();
228  if(useManualSeedPoint) {
229  CX_LOG_INFO() << "Using seed point: " << seedPoint.transpose();
230  segmentation->setSeedPoint(seedPoint(0), seedPoint(1), seedPoint(2));
231  }
232  segmentation->setInputConnection(importer->getOutputPort());
233  segmentation->update();
234  segmentationData = segmentation->getOutputData<fast::Segmentation>(0);
235 
236  // Convert fast segmentation data to VTK data which CX can use
237  vtkSmartPointer<fast::VTKImageExporter> vtkExporter = fast::VTKImageExporter::New();
238  vtkExporter->setInputConnection(segmentation->getOutputPort());
239  vtkExporter->Update();
240  mAirwaySegmentationOutput = vtkExporter->GetOutput();
241  }
242  } catch(fast::Exception & e)
243  {
244  CX_LOG_ERROR() << "The airways filter failed: \n"
245  << e.what();
246  if(!useManualSeedPoint)
247  CX_LOG_ERROR() << "Try to set the seed point manually.";
248  return false;
249  }
250 
251  CX_LOG_SUCCESS() << "FINISHED AIRWAY SEGMENTATION";
252 
253 
254  // Get the transformation of the segmentation
255  Eigen::Affine3f T = fast::SceneGraph::getEigenAffineTransformationFromData(segmentationData);
256  mTransformation.matrix() = T.matrix().cast<double>(); // cast to double
257 
258  // Extract centerline
259  fast::CenterlineExtraction::pointer centerline = fast::CenterlineExtraction::New();
260  centerline->setInputData(segmentationData);
261 
262  // Get centerline
263  vtkSmartPointer<fast::VTKMeshExporter> vtkCenterlineExporter = fast::VTKMeshExporter::New();
264  vtkCenterlineExporter->setInputConnection(centerline->getOutputPort());
265  mCenterlineOutput = vtkCenterlineExporter->GetOutput();
266  vtkCenterlineExporter->Update();
267 
268  } catch(fast::Exception& e) {
269  std::string error = e.what();
270  reportError("fast::Exception: "+qstring_cast(error));
271 
272  return false;
273  } catch(cl::Error& e) {
274  reportError("cl::Error:"+qstring_cast(e.what()));
275 
276  return false;
277  } catch (std::exception& e){
278  reportError("std::exception:"+qstring_cast(e.what()));
279 
280  return false;
281  } catch (...){
282  reportError("Airway segmentation algorithm threw a unknown exception.");
283 
284  return false;
285  }
286  return true;
287 }
288 
290 {
291  if(!mAirwaySegmentationOutput)
292  return false;
293 
294  std::cout << "POST PROCESS" << std::endl;
295 
296  // Make contour of segmented volume
297  double threshold = 1;
299  mAirwaySegmentationOutput,
300  threshold,
301  false, // reduce resolution
302  true, // smoothing
303  true, // keep topology
304  0 // target decimation
305  );
306  //outputSegmentation->get_rMd_History()->setRegistration(rMd_i);
307  //patientService()->insertData(outputSegmentation);
308 
309  // Add contour internally to cx
311  patientService(),
312  rawContour,
313  mInputImage,
314  QColor("green")
315  );
316  contour->get_rMd_History()->setRegistration(mTransformation);
317 
318  // Set output
319  mOutputTypes[1]->setValue(contour->getUid());
320 
321  if(getLungSegmentationOption(mOptions)->getValue()) {
323  mLungSegmentationOutput,
324  threshold,
325  false, // reduce resolution
326  true, // smoothing
327  true, // keep topology
328  0 // target decimation
329  );
330  //outputSegmentation->get_rMd_History()->setRegistration(rMd_i);
331  //patientService()->insertData(outputSegmentation);
332 
333  // Add contour internally to cx
334  QColor color("red");
335  color.setAlpha(100);
337  patientService(),
338  rawContour,
339  mInputImage,
340  color
341  );
342  contour->get_rMd_History()->setRegistration(mTransformation);
343 
344  // Set output
345  mOutputTypes[2]->setValue(contour->getUid());
346  }
347 
348  // Centerline
349  QString uid = mInputImage->getUid() + AirwaysFilter::getNameSuffix() + "%1";
350  QString name = mInputImage->getName() + AirwaysFilter::getNameSuffix() + "%1";
351  MeshPtr centerline = patientService()->createSpecificData<Mesh>(uid, name);
352  centerline->setVtkPolyData(mCenterlineOutput);
353  centerline->get_rMd_History()->setParentSpace(mInputImage->getUid());
354  centerline->get_rMd_History()->setRegistration(mTransformation);
355  patientService()->insertData(centerline);
356  mOutputTypes[0]->setValue(centerline->getUid());
357 
358  // Straight centerline and tubes
359  if(getStraightCLTubesOption(mOptions)->getValue())
360  {
361  this->createStraightCL();
362  this->createTubes();
363  }
364 
365  return true;
366 }
367 
384 void AirwaysFilter::createTubes()
385 {
386  // Get the straight centerline to model the tubes around.
387  QString straightCLUid = mOutputTypes[3]->getValue();
388  MeshPtr straightCL = boost::dynamic_pointer_cast<Mesh>(patientService()->getData(straightCLUid));
389  if(!straightCL)
390  return;
391  vtkPolyDataPtr clPolyData = straightCL->getVtkPolyData();
392 
393  // Create the implicit modeller
394  vtkSmartPointer<vtkImplicitModeller> blobbyLogoImp =
395  vtkSmartPointer<vtkImplicitModeller>::New();
396  blobbyLogoImp->SetInputData(clPolyData);
397  blobbyLogoImp->SetMaximumDistance(0.1);
398  blobbyLogoImp->SetSampleDimensions(256, 256, 256);
399  blobbyLogoImp->SetAdjustDistance(0.1);
400 
401  // Extract an iso surface, i.e. the tube shell
402  vtkSmartPointer<vtkContourFilter> blobbyLogoIso =
403  vtkSmartPointer<vtkContourFilter>::New();
404  blobbyLogoIso->SetInputConnection(blobbyLogoImp->GetOutputPort());
405  blobbyLogoIso->SetValue(1, 1.5); //orig
406  blobbyLogoIso->Update();
407 
408  // Create the mesh object from the tube shell
409  QString uid = mInputImage->getUid() + AirwaysFilter::getNameSuffix() + AirwaysFilter::getNameSuffixStraight() + AirwaysFilter::getNameSuffixTubes() + "%1";
410  QString name = mInputImage->getName() + AirwaysFilter::getNameSuffix() + AirwaysFilter::getNameSuffixStraight() + AirwaysFilter::getNameSuffixTubes() + "%1";
411  MeshPtr centerline = patientService()->createSpecificData<Mesh>(uid, name);
412  centerline->setVtkPolyData(blobbyLogoIso->GetOutput());
413  centerline->get_rMd_History()->setParentSpace(mInputImage->getUid());
414  centerline->get_rMd_History()->setRegistration(mTransformation);
415  // The color is taken from the new Fraxinus logo. Blue is the common color for lungs/airways. Partly transparent for a nice effect in Fraxinus.
416  centerline->setColor(QColor(118, 178, 226, 200));
417  patientService()->insertData(centerline);
418  mOutputTypes[4]->setValue(centerline->getUid());
419 
420 }
421 
422 void AirwaysFilter::setDefaultStraightCLTubesOption(bool defaultStraightCLTubesOption)
423 {
424  mDefaultStraightCLTubesOption = defaultStraightCLTubesOption;
425 }
426 
427 void AirwaysFilter::createStraightCL()
428 {
429  QString uid = mInputImage->getUid() + AirwaysFilter::getNameSuffix() + AirwaysFilter::getNameSuffixStraight() + "%1";
430  QString name = mInputImage->getName() + AirwaysFilter::getNameSuffix() + AirwaysFilter::getNameSuffixStraight() + "%1";
431  MeshPtr centerline = patientService()->createSpecificData<Mesh>(uid, name);
432 
434 
435  Eigen::MatrixXd CLpoints = makeTransformedMatrix(mCenterlineOutput);
436  bl->findBranchesInCenterline(CLpoints);
437  vtkPolyDataPtr retval = bl->createVtkPolyDataFromBranches(false, true);
438 
439  centerline->setVtkPolyData(retval);
440  centerline->get_rMd_History()->setParentSpace(mInputImage->getUid());
441  centerline->get_rMd_History()->setRegistration(mTransformation);
442  patientService()->insertData(centerline);
443  mOutputTypes[3]->setValue(centerline->getUid());
444 }
445 
447 {
448  mOptionsAdapters.push_back(this->getManualSeedPointOption(mOptions));
449  mOptionsAdapters.push_back(this->getLungSegmentationOption(mOptions));
450  mOptionsAdapters.push_back(this->getStraightCLTubesOption(mOptions));
451 }
452 
454 {
456 
458  temp->setValueName("Input");
459  temp->setHelp("Select input to run airway segmentation on.");
460  mInputTypes.push_back(temp);
461 }
462 
464 {
465  StringPropertySelectMeshPtr tempMeshStringAdapter;
466  std::vector<std::pair<QString, QString>> valueHelpPairs;
467  valueHelpPairs.push_back(std::make_pair(tr("Airway Centerline"), tr("Generated centerline mesh (vtk-format).")));
468  valueHelpPairs.push_back(std::make_pair(tr("Airway Segmentation"), tr("Generated surface of the airway segmentation volume.")));
469  valueHelpPairs.push_back(std::make_pair(tr("Lung Segmentation"), tr("Generated surface of the lung segmentation volume.")));
470  valueHelpPairs.push_back(std::make_pair(tr("Straight Airway Centerline"), tr("A centerline with straight lines between the branch points.")));
471  valueHelpPairs.push_back(std::make_pair(tr("Straight Airway Tubes"), tr("Tubes based on the straight centerline")));
472 
473  foreach(auto pair, valueHelpPairs)
474  {
475  tempMeshStringAdapter = StringPropertySelectMesh::New(patientService());
476  tempMeshStringAdapter->setValueName(pair.first);
477  tempMeshStringAdapter->setHelp(pair.second);
478  mOutputTypes.push_back(tempMeshStringAdapter);
479  }
480 }
481 
482 
483 BoolPropertyPtr AirwaysFilter::getManualSeedPointOption(QDomElement root)
484 {
485  BoolPropertyPtr retval =
486  BoolProperty::initialize("Use manual seed point",
487  "",
488  "If the automatic seed point detection algorithm fails you can use cursor to set the seed point "
489  "inside trachea of the patient. "
490  "Then tick this checkbox to use the manual seed point in the airways filter.",
491  false, root);
492  return retval;
493 
494 }
495 
496 BoolPropertyPtr AirwaysFilter::getLungSegmentationOption(QDomElement root)
497 {
498  BoolPropertyPtr retval =
499  BoolProperty::initialize("Lung segmentation",
500  "",
501  "Selecting this option will also segment the two lung sacs",
502  false, root);
503  return retval;
504 
505 }
506 
507 BoolPropertyPtr AirwaysFilter::getStraightCLTubesOption(QDomElement root)
508 {
509  BoolPropertyPtr retval =
510  BoolProperty::initialize("Straight centerline and tubes",
511  "",
512  "Use this option to generate a centerline with straight branches between "
513  "the branch points. "
514  "You also get tubes based on this straight line.",
515  mDefaultStraightCLTubesOption, root);
516  return retval;
517 }
518 
519 
520 } /* namespace cx */
521 
boost::shared_ptr< class SpaceProvider > SpaceProviderPtr
QString qstring_cast(const T &val)
static BoolPropertyPtr initialize(const QString &uid, QString name, QString help, bool value, QDomNode root=QDomNode())
std::vector< SelectDataStringPropertyBasePtr > mInputTypes
Definition: cxFilterImpl.h:73
void reportError(QString msg)
Definition: cxLogger.cpp:71
A mesh data set.
Definition: cxMesh.h:45
boost::shared_ptr< class VisServices > VisServicesPtr
Definition: cxMainWindow.h:40
boost::shared_ptr< class BranchList > BranchListPtr
#define CX_LOG_INFO
Definition: cxLogger.h:96
boost::shared_ptr< class Image > ImagePtr
Definition: cxDicomWidget.h:27
virtual vtkPolyDataPtr getVtkPolyData() const
Definition: cxMesh.cpp:129
virtual vtkImageDataPtr getBaseVtkImageData()
Definition: cxImage.cpp:335
virtual void createOutputTypes()
AirwaysFilter(VisServicesPtr services)
std::vector< PropertyPtr > mOptionsAdapters
Definition: cxFilterImpl.h:75
boost::shared_ptr< class Data > DataPtr
static QString findConfigFolder(QString pathRelativeToConfigRoot, QString alternativeAbsolutePath="")
VisServicesPtr mServices
Definition: cxFilterImpl.h:82
boost::shared_ptr< class SelectDataStringPropertyBase > SelectDataStringPropertyBasePtr
PatientModelServicePtr patientService()
virtual bool execute()
static QString getCachePath()
return path to a folder that is used during execution, will be cleared at start and stop...
void setVtkPolyData(const vtkPolyDataPtr &polyData)
Definition: cxMesh.cpp:91
#define CX_LOG_ERROR
Definition: cxLogger.h:99
A volumetric data set.
Definition: cxImage.h:45
int * getImageSize(DataPtr inputImage)
#define CX_LOG_SUCCESS
Definition: cxLogger.h:97
QDomElement mOptions
Definition: cxFilterImpl.h:76
virtual QString getType() const
vtkSmartPointer< vtkPolyData > vtkPolyDataPtr
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:42
std::vector< SelectDataStringPropertyBasePtr > mOutputTypes
Definition: cxFilterImpl.h:74
static QString getNameSuffix()
static QString getNameSuffixStraight()
void setDefaultStraightCLTubesOption(bool defaultStraightCLTubesOption)
static StringPropertySelectMeshPtr New(PatientModelServicePtr patientModelService)
virtual bool postProcess()
static StringPropertySelectImagePtr New(PatientModelServicePtr patientModelService)
boost::shared_ptr< class BoolProperty > BoolPropertyPtr
boost::shared_ptr< class Mesh > MeshPtr
Eigen::MatrixXd makeTransformedMatrix(vtkPolyDataPtr linesPolyData, Transform3D rMd)
makeTransformedMatrix This method takes an vtkpolydata as input, runs it through a transform and retu...
static QString getNameSuffixTubes()
boost::shared_ptr< class StringPropertySelectMesh > StringPropertySelectMeshPtr
virtual QString getName() const
virtual QString getHelp() const
Namespace for all CustusX production code.