Fraxinus  18.10
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"
25 
26 #include "cxTime.h"
27 #include "cxTypeConversions.h"
28 #include "cxLogger.h"
29 #include "cxDataReaderWriter.h"
31 #include "cxDoubleProperty.h"
32 #include "cxContourFilter.h"
33 #include "cxDataLocations.h"
35 #include "vtkForwardDeclarations.h"
37 #include "cxVisServices.h"
38 #include "cxUtilHelpers.h"
39 #include "FAST/Algorithms/LungSegmentation/LungSegmentation.hpp"
40 #include "FAST/Algorithms/AirwaySegmentation/AirwaySegmentation.hpp"
41 #include "FAST/Algorithms/CenterlineExtraction/CenterlineExtraction.hpp"
42 #include "FAST/Importers/ImageFileImporter.hpp"
43 #include "FAST/Exporters/VTKImageExporter.hpp"
44 #include "FAST/Exporters/VTKMeshExporter.hpp"
45 #include "FAST/Data/Segmentation.hpp"
46 #include "FAST/SceneGraph.hpp"
47 
48 namespace cx {
49 
51  FilterImpl(services),
52  mDefaultStraightCLTubesOption(false)
53 {
54  fast::Reporter::setGlobalReportMethod(fast::Reporter::COUT);
55  //Need to create OpenGL context of fast in main thread, this is done in the constructor of DeviceManger
56  fast::ImageFileImporter::pointer importer = fast::ImageFileImporter::New();
57  Q_UNUSED(importer)
58 }
59 
60 
62 }
63 
64 QString AirwaysFilter::getName() const
65 {
66  return "Airway Segmentation Filter";
67 }
68 
69 QString AirwaysFilter::getType() const
70 {
71  return "airways_filter";
72 }
73 
74 QString AirwaysFilter::getHelp() const
75 {
76  return "<html>"
77  "<h3>Airway Segmentation.</h3>"
78  "<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>"
79  "</html>";
80 }
81 
83 {
84  return "_centerline";
85 }
86 
88 {
89  return "_straight";
90 }
91 
93 {
94  return "_tubes";
95 }
96 
97 Vector3D AirwaysFilter::getSeedPointFromTool(SpaceProviderPtr spaceProvider, DataPtr data)
98 {
99  // Retrieve position of tooltip and use it as seed point
100  Vector3D point = spaceProvider->getActiveToolTipPoint(
101  spaceProvider->getD(data));
102 
103  // Have to multiply by the inverse of the spacing to get the voxel position
104  ImagePtr image = boost::dynamic_pointer_cast<Image>(data);
105  double spacingX, spacingY, spacingZ;
106  image->getBaseVtkImageData()->GetSpacing(spacingX, spacingY, spacingZ);
107  point(0) = point(0) * (1.0 / spacingX);
108  point(1) = point(1) * (1.0 / spacingY);
109  point(2) = point(2) * (1.0 / spacingZ);
110 
111  std::cout << "the selected seed point is: " << point(0) << " " << point(1)
112  << " " << point(2) << "\n";
113 
114  return point;
115 }
116 
117 int * getImageSize(DataPtr inputImage)
118 {
119  ImagePtr image = boost::dynamic_pointer_cast<Image>(inputImage);
120  return image->getBaseVtkImageData()->GetDimensions();
121 }
122 
123 bool AirwaysFilter::isSeedPointInsideImage(Vector3D seedPoint, DataPtr image)
124 {
125  int * size = getImageSize(image);
126  std::cout << "size of image is: " << size[0] << " " << size[1] << " "
127  << size[2] << "\n";
128  int x = (int) seedPoint(0);
129  int y = (int) seedPoint(1);
130  int z = (int) seedPoint(2);
131  bool result = x >= 0 && y >= 0 && z >= 0 && x < size[0] && y < size[1]
132  && z < size[2];
133  return result;
134 }
135 
137 {
138  DataPtr inputImage = mInputTypes[0].get()->getData();
139  if (!inputImage)
140  {
141  CX_LOG_ERROR() << "No input data selected";
142  return false;
143  }
144 
145  if (inputImage->getType() != "image")
146  {
147  CX_LOG_ERROR() << "Input data has to be an image";
148  return false;
149  }
150 
151  std::string filename = (patientService()->getActivePatientFolder()
152  + "/" + inputImage->getFilename()).toStdString();
153 
154  // only check seed point inside image if use seed point is checked
155  bool useManualSeedPoint = getManualSeedPointOption(mOptions)->getValue();
156  if(useManualSeedPoint)
157  {
158  seedPoint = getSeedPointFromTool(mServices->spaceProvider(), inputImage);
159  if(!isSeedPointInsideImage(seedPoint, inputImage)) {
160  CX_LOG_ERROR() << "Seed point is not inside image. Use cursor to set seed point inside trachea in the CT image.";
161  return false;
162  }
163  }
164  mInputImage = patientService()->getData<Image>(inputImage->getUid());
165 
166  return true;
167 }
168 
170 {
171  CX_LOG_INFO() << "EXECUTING AIRWAYS FILTER";
172  // Check if pre process went ok:
173  if(!mInputImage)
174  return false;
175 
176  QString q_filename = "";
177  QString activePatienFolder = patientService()->getActivePatientFolder();
178  QString inputImageFileName = mInputImage->getFilename();
179  if(!activePatienFolder.isEmpty())
180  q_filename = activePatienFolder+"/"+inputImageFileName;
181  else
182  q_filename = inputImageFileName;
183 
184  std::string filename = q_filename.toStdString();
185  try {
186  fast::Config::getTestDataPath(); // needed for initialization
187  QString cacheDir = cx::DataLocations::getCachePath();
188  fast::Config::setKernelBinaryPath(cacheDir.toStdString());
189  QString kernelDir = cx::DataLocations::findConfigFolder("/FAST", FAST_SOURCE_DIR);
190  fast::Config::setKernelSourcePath(kernelDir.toStdString());
191 
192  // Import image data from disk
193  fast::ImageFileImporter::pointer importer = fast::ImageFileImporter::New();
194  importer->setFilename(filename);
195 
196  // Need to know the data type
197  importer->update();
198  fast::Image::pointer image = importer->getOutputData<fast::Image>();
199 
200  // Do segmentation
201  fast::Segmentation::pointer segmentationData;
202  bool doLungSegmentation = getLungSegmentationOption(mOptions)->getValue();
203  bool useManualSeedPoint = getManualSeedPointOption(mOptions)->getValue();
204  try {
205  if(doLungSegmentation) {
206  fast::LungSegmentation::pointer segmentation = fast::LungSegmentation::New();
207  if(useManualSeedPoint) {
208  CX_LOG_INFO() << "Using seed point: " << seedPoint.transpose();
209  segmentation->setAirwaySeedPoint(seedPoint(0), seedPoint(1), seedPoint(2));
210  }
211  segmentation->setInputConnection(importer->getOutputPort());
212  segmentation->update();
213  segmentationData = segmentation->getOutputData<fast::Segmentation>(1);
214 
215  // Convert fast segmentation data to VTK data which CX can use (Airways)
216  vtkSmartPointer<fast::VTKImageExporter> vtkExporter = fast::VTKImageExporter::New();
217  vtkExporter->setInputConnection(segmentation->getOutputPort(1));
218  vtkExporter->Update();
219  mAirwaySegmentationOutput = vtkExporter->GetOutput();
220 
221  // Convert fast segmentation data to VTK data which CX can use (Lungs)
222  vtkSmartPointer<fast::VTKImageExporter> vtkExporter2 = fast::VTKImageExporter::New();
223  vtkExporter2->setInputConnection(segmentation->getOutputPort(0));
224  vtkExporter2->Update();
225  mLungSegmentationOutput = vtkExporter2->GetOutput();
226  } else {
227 
228  fast::AirwaySegmentation::pointer segmentation = fast::AirwaySegmentation::New();
229  if(useManualSeedPoint) {
230  CX_LOG_INFO() << "Using seed point: " << seedPoint.transpose();
231  segmentation->setSeedPoint(seedPoint(0), seedPoint(1), seedPoint(2));
232  }
233  segmentation->setInputConnection(importer->getOutputPort());
234  segmentation->update();
235  segmentationData = segmentation->getOutputData<fast::Segmentation>(0);
236 
237  // Convert fast segmentation data to VTK data which CX can use
238  vtkSmartPointer<fast::VTKImageExporter> vtkExporter = fast::VTKImageExporter::New();
239  vtkExporter->setInputConnection(segmentation->getOutputPort());
240  vtkExporter->Update();
241  mAirwaySegmentationOutput = vtkExporter->GetOutput();
242  }
243  } catch(fast::Exception & e)
244  {
245  CX_LOG_ERROR() << "The airways filter failed: \n"
246  << e.what();
247  if(!useManualSeedPoint)
248  CX_LOG_ERROR() << "Try to set the seed point manually.";
249  return false;
250  }
251 
252  CX_LOG_SUCCESS() << "FINISHED AIRWAY SEGMENTATION";
253 
254 
255  // Get the transformation of the segmentation
256  Eigen::Affine3f T = fast::SceneGraph::getEigenAffineTransformationFromData(segmentationData);
257  mTransformation.matrix() = T.matrix().cast<double>(); // cast to double
258 
259  // Extract centerline
260  fast::CenterlineExtraction::pointer centerline = fast::CenterlineExtraction::New();
261  centerline->setInputData(segmentationData);
262 
263  // Get centerline
264  vtkSmartPointer<fast::VTKMeshExporter> vtkCenterlineExporter = fast::VTKMeshExporter::New();
265  vtkCenterlineExporter->setInputConnection(centerline->getOutputPort());
266  mCenterlineOutput = vtkCenterlineExporter->GetOutput();
267  vtkCenterlineExporter->Update();
268 
269  } catch(fast::Exception& e) {
270  std::string error = e.what();
271  reportError("fast::Exception: "+qstring_cast(error));
272 
273  return false;
274  } catch(cl::Error& e) {
275  reportError("cl::Error:"+qstring_cast(e.what()));
276 
277  return false;
278  } catch (std::exception& e){
279  reportError("std::exception:"+qstring_cast(e.what()));
280 
281  return false;
282  } catch (...){
283  reportError("Airway segmentation algorithm threw a unknown exception.");
284 
285  return false;
286  }
287  return true;
288 }
289 
291 {
292  if(!mAirwaySegmentationOutput)
293  return false;
294 
295  std::cout << "POST PROCESS" << std::endl;
296 
297  // Make contour of segmented volume
298  double threshold = 1;
300  mAirwaySegmentationOutput,
301  threshold,
302  false, // reduce resolution
303  true, // smoothing
304  true, // keep topology
305  0 // target decimation
306  );
307  //outputSegmentation->get_rMd_History()->setRegistration(rMd_i);
308  //patientService()->insertData(outputSegmentation);
309 
310  // Add contour internally to cx
312  patientService(),
313  rawContour,
314  mInputImage,
315  QColor("green")
316  );
317  contour->get_rMd_History()->setRegistration(mTransformation);
318 
319  // Set output
320  mOutputTypes[1]->setValue(contour->getUid());
321 
322  if(getLungSegmentationOption(mOptions)->getValue()) {
324  mLungSegmentationOutput,
325  threshold,
326  false, // reduce resolution
327  true, // smoothing
328  true, // keep topology
329  0 // target decimation
330  );
331  //outputSegmentation->get_rMd_History()->setRegistration(rMd_i);
332  //patientService()->insertData(outputSegmentation);
333 
334  // Add contour internally to cx
335  QColor color("red");
336  color.setAlpha(100);
338  patientService(),
339  rawContour,
340  mInputImage,
341  color
342  );
343  contour->get_rMd_History()->setRegistration(mTransformation);
344 
345  // Set output
346  mOutputTypes[2]->setValue(contour->getUid());
347  }
348 
349  // Centerline
350  QString uid = mInputImage->getUid() + AirwaysFilter::getNameSuffixCenterline() + "%1";
351  QString name = mInputImage->getName() + AirwaysFilter::getNameSuffixCenterline() + "%1";
352  MeshPtr centerline = patientService()->createSpecificData<Mesh>(uid, name);
353  centerline->setVtkPolyData(mCenterlineOutput);
354  centerline->get_rMd_History()->setParentSpace(mInputImage->getUid());
355  centerline->get_rMd_History()->setRegistration(mTransformation);
356  patientService()->insertData(centerline);
357  mOutputTypes[0]->setValue(centerline->getUid());
358 
359  this->createAirwaysFromCenterline();
360 
361  return true;
362 }
363 
364 
365 void AirwaysFilter::createAirwaysFromCenterline()
366 {
368 
369  airwaysFromCLPtr->processCenterline(mCenterlineOutput);
370 
371  // Create the mesh object from the airway walls
372  QString uidMesh = mInputImage->getUid() + AirwaysFilter::getNameSuffixTubes() + "%1";
373  QString nameMesh = mInputImage->getName()+ AirwaysFilter::getNameSuffixTubes() + "%1";
374  MeshPtr airwayWalls = patientService()->createSpecificData<Mesh>(uidMesh, nameMesh);
375  airwayWalls->setVtkPolyData(airwaysFromCLPtr->generateTubes());
376  airwayWalls->get_rMd_History()->setParentSpace(mInputImage->getUid());
377  airwayWalls->get_rMd_History()->setRegistration(mTransformation);
378  airwayWalls->setColor(QColor(253, 173, 136, 255));
379  patientService()->insertData(airwayWalls);
380  mOutputTypes[4]->setValue(airwayWalls->getUid());
381 
382 
383  //insert filtered centerline from airwaysFromCenterline
384  QString uidCenterline = mInputImage->getUid() + AirwaysFilter::getNameSuffixTubes() + AirwaysFilter::getNameSuffixCenterline() + "%1";
385  QString nameCenterline = mInputImage->getName() + AirwaysFilter::getNameSuffixTubes() + AirwaysFilter::getNameSuffixCenterline() + "%1";
386  MeshPtr centerline = patientService()->createSpecificData<Mesh>(uidCenterline, nameCenterline);
387  centerline->setVtkPolyData(airwaysFromCLPtr->getVTKPoints());
388  centerline->get_rMd_History()->setParentSpace(mInputImage->getUid());
389  centerline->get_rMd_History()->setRegistration(mTransformation);
390  patientService()->insertData(centerline);
391  mOutputTypes[3]->setValue(centerline->getUid());
392 }
393 
394 void AirwaysFilter::setDefaultStraightCLTubesOption(bool defaultStraightCLTubesOption)
395 {
396  mDefaultStraightCLTubesOption = defaultStraightCLTubesOption;
397 }
398 
400 {
401  mOptionsAdapters.push_back(this->getManualSeedPointOption(mOptions));
402  mOptionsAdapters.push_back(this->getLungSegmentationOption(mOptions));
403 }
404 
406 {
408 
410  temp->setValueName("Input");
411  temp->setHelp("Select input to run airway segmentation on.");
412  mInputTypes.push_back(temp);
413 }
414 
416 {
417  StringPropertySelectMeshPtr tempMeshStringAdapter;
418  std::vector<std::pair<QString, QString>> valueHelpPairs;
419  valueHelpPairs.push_back(std::make_pair(tr("Airway Centerline"), tr("Generated centerline mesh (vtk-format).")));
420  valueHelpPairs.push_back(std::make_pair(tr("Airway Segmentation"), tr("Generated surface of the airway segmentation volume.")));
421  valueHelpPairs.push_back(std::make_pair(tr("Lung Segmentation"), tr("Generated surface of the lung segmentation volume.")));
422  valueHelpPairs.push_back(std::make_pair(tr("Straight Airway Centerline"), tr("A centerline with straight lines between the branch points.")));
423  valueHelpPairs.push_back(std::make_pair(tr("Straight Airway Tubes"), tr("Tubes based on the straight centerline")));
424 
425  foreach(auto pair, valueHelpPairs)
426  {
427  tempMeshStringAdapter = StringPropertySelectMesh::New(patientService());
428  tempMeshStringAdapter->setValueName(pair.first);
429  tempMeshStringAdapter->setHelp(pair.second);
430  mOutputTypes.push_back(tempMeshStringAdapter);
431  }
432 }
433 
434 
435 BoolPropertyPtr AirwaysFilter::getManualSeedPointOption(QDomElement root)
436 {
437  BoolPropertyPtr retval =
438  BoolProperty::initialize("Use manual seed point",
439  "",
440  "If the automatic seed point detection algorithm fails you can use cursor to set the seed point "
441  "inside trachea of the patient. "
442  "Then tick this checkbox to use the manual seed point in the airways filter.",
443  false, root);
444  return retval;
445 
446 }
447 
448 BoolPropertyPtr AirwaysFilter::getLungSegmentationOption(QDomElement root)
449 {
450  BoolPropertyPtr retval =
451  BoolProperty::initialize("Lung segmentation",
452  "",
453  "Selecting this option will also segment the two lung sacs",
454  false, root);
455  return retval;
456 
457 }
458 
459 } /* namespace cx */
460 
boost::shared_ptr< class SpaceProvider > SpaceProviderPtr
QString qstring_cast(const T &val)
boost::shared_ptr< AirwaysFromCenterline > AirwaysFromCenterlinePtr
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
#define CX_LOG_INFO
Definition: cxLogger.h:96
boost::shared_ptr< class Image > ImagePtr
Definition: cxDicomWidget.h:27
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
static QString getNameSuffixCenterline()
std::vector< SelectDataStringPropertyBasePtr > mOutputTypes
Definition: cxFilterImpl.h:74
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
static QString getNameSuffixTubes()
boost::shared_ptr< class StringPropertySelectMesh > StringPropertySelectMeshPtr
virtual QString getName() const
virtual QString getHelp() const
Namespace for all CustusX production code.