CustusX  16.5
An IGT application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
cxLevelSetFilterService.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) 2008-2014, SINTEF Department of Medical Technology
5 All rights reserved.
6 
7 Redistribution and use in source and binary forms, with or without
8 modification, are permitted provided that the following conditions are met:
9 
10 1. Redistributions of source code must retain the above copyright notice,
11  this list of conditions and the following disclaimer.
12 
13 2. Redistributions in binary form must reproduce the above copyright notice,
14  this list of conditions and the following disclaimer in the documentation
15  and/or other materials provided with the distribution.
16 
17 3. Neither the name of the copyright holder nor the names of its contributors
18  may be used to endorse or promote products derived from this software
19  without specific prior written permission.
20 
21 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
24 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
25 FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
26 DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
27 SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
29 OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 =========================================================================*/
32 
33 #include "cxLevelSetFilterService.h"
34 
35 
36 #include <ctkPluginContext.h>
37 #include "cxTime.h"
38 #include "cxMesh.h"
39 #include "cxTypeConversions.h"
40 #include "cxReporter.h"
41 #include "cxRegistrationTransform.h"
42 #include "cxStringProperty.h"
43 #include "cxDoubleProperty.h"
44 #include "cxSelectDataStringProperty.h"
45 #include "cxCoordinateSystemHelpers.h"
46 #include "cxData.h"
47 #include "cxImage.h"
48 #include <vtkImageImport.h>
49 #include <vtkImageData.h>
50 #include <vtkImageShiftScale.h>
51 #include <vtkImageData.h>
52 #include "cxContourFilter.h"
53 #include "cxDataLocations.h"
54 
55 #include <itkBinaryMorphologicalClosingImageFilter.h>
56 #include <itkBinaryBallStructuringElement.h>
57 #include "cxAlgorithmHelpers.h"
58 #include <vtkImageCast.h>
59 
60 #include "levelSet.hpp"
61 #include "OpenCLManager.hpp"
62 #include "HelperFunctions.hpp"
63 #include "cxSpaceProvider.h"
64 #include "cxPatientModelServiceProxy.h"
65 
66 #include "level-set-segmentation-config.h"
67 #include "OulConfig.hpp"
68 #include "cxVisServices.h"
69 
70 namespace cx
71 {
72 
73 LevelSetFilter::LevelSetFilter(ctkPluginContext *pluginContext) :
74  FilterImpl(VisServices::create(pluginContext))
75 {
76 }
77 
78 QString LevelSetFilter::getName() const
79 {
80  return "Level Set Segmentation";
81 }
82 
83 QString LevelSetFilter::getType() const
84 {
85  return "LevelSetFilter";
86 }
87 
88 QString LevelSetFilter::getHelp() const
89 {
90  return "<html>"
91  "<h3>Level Set Segmentation Filter.</h3>"
92  "</html>";
93 }
94 
95 Vector3D LevelSetFilter::getSeedPointFromTool(SpaceProviderPtr spaceProvider, DataPtr data)
96 {
97  // Retrieve position of tooltip and use it as seed point
98  Vector3D point = spaceProvider->getActiveToolTipPoint(
99  spaceProvider->getD(data));
100 
101  // Have to multiply by the inverse of the spacing to get the voxel position
102  ImagePtr image = boost::dynamic_pointer_cast<Image>(data);
103  double spacingX, spacingY, spacingZ;
104  image->getBaseVtkImageData()->GetSpacing(spacingX, spacingY, spacingZ);
105  point(0) = point(0) * (1.0 / spacingX);
106  point(1) = point(1) * (1.0 / spacingY);
107  point(2) = point(2) * (1.0 / spacingZ);
108 
109  std::cout << "the selected seed point is: " << point(0) << " " << point(1)
110  << " " << point(2) << "\n";
111 
112  return point;
113 }
114 
115 int * getImageSize(DataPtr inputImage)
116 {
117  ImagePtr image = boost::dynamic_pointer_cast<Image>(inputImage);
118  return image->getBaseVtkImageData()->GetDimensions();
119 }
120 
121 bool LevelSetFilter::isSeedPointInsideImage(Vector3D seedPoint, DataPtr image)
122 {
123  int * size = getImageSize(image);
124  std::cout << "size of image is: " << size[0] << " " << size[1] << " "
125  << size[2] << "\n";
126  int x = (int) seedPoint(0);
127  int y = (int) seedPoint(1);
128  int z = (int) seedPoint(2);
129  bool result = x >= 0 && y >= 0 && z >= 0 && x < size[0] && y < size[1]
130  && z < size[2];
131  return result;
132 }
133 
134 bool LevelSetFilter::preProcess()
135 {
136  DataPtr inputImage = mInputTypes[0].get()->getData();
137  if (!inputImage)
138  {
139  std::cout << "No input data selected" << std::endl;
140  return false;
141  }
142 
143  if (inputImage->getType() != "image")
144  {
145  std::cout << "Input data has to be an image" << std::endl;
146  return false;
147  }
148 
149  filename = (patientService()->getActivePatientFolder()
150  + "/" + inputImage->getFilename()).toStdString();
151 
152  seedPoint = getSeedPointFromTool(mServices->spaceProvider(), inputImage);
153  if (!isSeedPointInsideImage(seedPoint, inputImage))
154  {
155  std::cout << "Seed point is not inside image!" << std::endl;
156  return false;
157  }
158  image = patientService()->getData<Image>(inputImage->getUid());
159 
160  return true;
161 }
162 
163 bool LevelSetFilter::execute()
164 {
165  DataPtr inputImage = mInputTypes[0].get()->getData();
166 
167  float threshold = getThresholdOption(mOptions)->getValue();
168  float epsilon = getEpsilonOption(mOptions)->getValue();
169  float alpha = getAlphaOption(mOptions)->getValue();
170 
171  std::cout << "Parameters are set to: " << threshold << " " << epsilon << " "
172  << alpha << std::endl;
173 
174  // Run level set segmentation
175  SIPL::int3 seed(seedPoint(0), seedPoint(1), seedPoint(2));
176  try
177  {
178  QString kernelDir = cx::DataLocations::findConfigFolder("/lss", KERNELS_DIR);
179  QString oulDir = cx::DataLocations::findConfigFolder("/tsf", OUL_DIR);
180 
181  SIPL::Volume<char> * result = runLevelSet(filename.c_str(), seed, 10, // seed radius
182  1000, // iterations per narrow band
183  threshold, epsilon, alpha, kernelDir.toStdString(), oulDir.toStdString());
184  SIPL::int3 size = result->getSize();
185  rawSegmentation = this->convertToVtkImageData(
186  (char *) result->getData(), size.x, size.y, size.z, image);
187  delete result;
188 
189  return true;
190  } catch (SIPL::SIPLException &e)
191  {
192  std::string error = e.what();
193  reportError(
194  "SIPL::SIPLException: " + qstring_cast(error));
195 
196  return false;
197  } catch (cl::Error &e)
198  {
199  if (e.err() == CL_MEM_OBJECT_ALLOCATION_FAILURE
200  || e.err() == CL_OUT_OF_RESOURCES)
201  {
202  reportError(
203  "cl::Error: Not enough memory on the device to process this image. ("
204  + qstring_cast(oul::getCLErrorString(e.err()))
205  + ")");
206  }
207  else
208  {
209  reportError(
210  "cl::Error: "
211  + qstring_cast(oul::getCLErrorString(e.err())));
212  }
213  return false;
214  } catch (...)
215  {
216  reportError("Unknown exception occured.");
217  return false;
218  }
219 }
220 
221 bool LevelSetFilter::postProcess()
222 {
223  if(!rawSegmentation)
224  return false;
225 
226  //add segmentation internally to cx
227  QString uidSegmentation = image->getUid() + "_seg%1";
228  QString nameSegmentation = image->getName() + "_seg%1";
229 
230  ImagePtr outputSegmentation2 = patientService()->createSpecificData<Image>(uidSegmentation, nameSegmentation);
231  outputSegmentation2->intitializeFromParentImage(image);
232  outputSegmentation2->setVtkImageData(rawSegmentation);
233 // ImagePtr outputSegmentation2 = dataManager()->createDerivedImage(
234 // rawSegmentation, uidSegmentation, nameSegmentation, image);
235  ImagePtr outputSegmentation;
236  if (!outputSegmentation2)
237  return false;
238 
239  float radius = getRadiusOption(mOptions)->getValue();
240  if (radius > 0)
241  {
242  std::cout << "Performing morphological closing on segmentation result"
243  << std::endl;
244 
245  // Convert radius in mm to radius in voxels for the structuring element
246  Eigen::Array3d spacing = image->getSpacing();
247  itk::Size<3> radiusInVoxels;
248  radiusInVoxels[0] = radius / spacing(0);
249  radiusInVoxels[1] = radius / spacing(1);
250  radiusInVoxels[2] = radius / spacing(2);
251 
252  itkImageType::ConstPointer itkImage =
253  AlgorithmHelper::getITKfromSSCImage(outputSegmentation2);
254 
255  // Create structuring element
256  typedef itk::BinaryBallStructuringElement<unsigned char, 3> StructuringElementType;
257  StructuringElementType structuringElement;
258  structuringElement.SetRadius(radiusInVoxels);
259  structuringElement.CreateStructuringElement();
260 
261  // Morphological closing
262  typedef itk::BinaryMorphologicalClosingImageFilter<itkImageType,
263  itkImageType, StructuringElementType> closingFilterType;
264  closingFilterType::Pointer closingFilter = closingFilterType::New();
265  closingFilter->SetInput(itkImage);
266  closingFilter->SetKernel(structuringElement);
267  closingFilter->SetForegroundValue(1);
268  closingFilter->Update();
269  itkImage = closingFilter->GetOutput();
270 
271  //Convert ITK to VTK
272  itkToVtkFilterType::Pointer itkToVtkFilter = itkToVtkFilterType::New();
273  itkToVtkFilter->SetInput(itkImage);
274  itkToVtkFilter->Update();
275 
276  vtkImageDataPtr rawResult = vtkImageDataPtr::New();
277  rawResult->DeepCopy(itkToVtkFilter->GetOutput());
278 
279  vtkImageCastPtr imageCast = vtkImageCastPtr::New();
280  imageCast->SetInputData(rawResult);
281  imageCast->SetOutputScalarTypeToUnsignedChar();
282  rawResult = imageCast->GetOutput();
283 
284  outputSegmentation = patientService()->createSpecificData<Image>(uidSegmentation, nameSegmentation);
285  outputSegmentation->intitializeFromParentImage(image);
286  outputSegmentation->setVtkImageData(rawResult);
287  rawSegmentation = rawResult;
288  }
289  else
290  {
291  outputSegmentation = outputSegmentation2;
292  }
293 
294  //make contour of segmented volume
295  double threshold = 1;
296  vtkPolyDataPtr rawContour = ContourFilter::execute(rawSegmentation,
297  threshold);
298  Transform3D rMd_i = image->get_rMd(); //transform from the volumes coordinate system to our reference coordinate system
299  outputSegmentation->get_rMd_History()->setRegistration(rMd_i);
300  patientService()->insertData(outputSegmentation);
301 
302  //add contour internally to cx
303  MeshPtr contour = ContourFilter::postProcess(patientService(), rawContour, image,
304  QColor("blue"));
305  contour->get_rMd_History()->setRegistration(rMd_i);
306 
307  //set output
308  mOutputTypes[0]->setValue(outputSegmentation->getUid());
309  mOutputTypes[1]->setValue(contour->getUid());
310 
311  return true;
312 }
313 
314 void LevelSetFilter::createOptions()
315 {
316  mOptionsAdapters.push_back(this->getThresholdOption(mOptions));
317  mOptionsAdapters.push_back(this->getEpsilonOption(mOptions));
318  mOptionsAdapters.push_back(this->getAlphaOption(mOptions));
319  mOptionsAdapters.push_back(this->getRadiusOption(mOptions));
320 }
321 
322 QDomElement LevelSetFilter::getmOptions()
323 {
324  return this->mOptions;
325 }
326 
327 void LevelSetFilter::createInputTypes()
328 {
329  SelectDataStringPropertyBasePtr temp;
330 
331  temp = StringPropertySelectImage::New(patientService());
332  temp->setValueName("Input");
333  temp->setHelp("Select image input for thresholding");
334  mInputTypes.push_back(temp);
335 }
336 
337 void LevelSetFilter::createOutputTypes()
338 {
339  SelectDataStringPropertyBasePtr temp;
340 
341  temp = StringPropertySelectData::New(patientService());
342  temp->setValueName("Output");
343  temp->setHelp("Output segmented binary image");
344  mOutputTypes.push_back(temp);
345 
346  temp = StringPropertySelectData::New(patientService());
347  temp->setValueName("Contour");
348  temp->setHelp("Output contour generated from thresholded binary image.");
349  mOutputTypes.push_back(temp);
350 }
351 
352 void LevelSetFilter::setActive(bool on)
353 {
354  FilterImpl::setActive(on);
355 
356  /*
357  if (!mActive)
358  RepManager::getInstance()->getThresholdPreview()->removePreview();
359  */
360 }
361 
362 vtkImageDataPtr LevelSetFilter::convertToVtkImageData(char * data, int size_x,
363  int size_y, int size_z, ImagePtr input)
364 {
365  vtkImageDataPtr retval = this->importRawImageData((void*) data, size_x,
366  size_y, size_z, input, VTK_UNSIGNED_CHAR);
367  return retval;
368 }
369 
370 //From vtkType.h (on Ubuntu 12.04)
371 //#define VTK_VOID 0
372 //#define VTK_BIT 1
373 //#define VTK_CHAR 2
374 //#define VTK_SIGNED_CHAR 15
375 //#define VTK_UNSIGNED_CHAR 3
376 //#define VTK_SHORT 4
377 //#define VTK_UNSIGNED_SHORT 5
378 //#define VTK_INT 6
379 //#define VTK_UNSIGNED_INT 7
380 //#define VTK_LONG 8
381 //#define VTK_UNSIGNED_LONG 9
382 //#define VTK_FLOAT 10
383 //#define VTK_DOUBLE 11
384 //#define VTK_ID_TYPE 12
385 vtkImageDataPtr LevelSetFilter::importRawImageData(void * data, int size_x,
386  int size_y, int size_z, ImagePtr input, int type)
387 {
388  vtkImageImportPtr imageImport = vtkImageImportPtr::New();
389 
390  imageImport->SetWholeExtent(0, size_x - 1, 0, size_y - 1, 0, size_z - 1);
391  imageImport->SetDataExtentToWholeExtent();
392  imageImport->SetDataScalarType(type);
393  imageImport->SetNumberOfScalarComponents(1);
394  imageImport->SetDataSpacing(input->getBaseVtkImageData()->GetSpacing());
395  imageImport->SetImportVoidPointer(data);
396  imageImport->Update();
397  imageImport->Modified();
398 
399  vtkImageDataPtr retval = vtkImageDataPtr::New();
400  retval->DeepCopy(imageImport->GetOutput());
401 
402  return retval;
403 }
404 
405 DoublePropertyPtr LevelSetFilter::getThresholdOption(QDomElement root)
406 {
407  DoublePropertyPtr retval = DoubleProperty::initialize(
408  "Threshold", "", "Select threshold for the segmentation", 1,
409  DoubleRange(-5000, 5000, 0.0000001), 0, root);
410  return retval;
411 
412 }
413 
414 DoublePropertyPtr LevelSetFilter::getEpsilonOption(QDomElement root)
415 {
416  DoublePropertyPtr retval = DoubleProperty::initialize("Epsilon",
417  "", "Select epsilon for the segmentation", 1,
418  DoubleRange(-5000, 5000, 0.0000001), 0, root);
419  return retval;
420 
421 }
422 
423 DoublePropertyPtr LevelSetFilter::getAlphaOption(QDomElement root)
424 {
425  DoublePropertyPtr retval = DoubleProperty::initialize("Alpha",
426  "", "Select alpha for the segmentation", 0.1,
427  DoubleRange(0, 1, 0.01), 2, root);
428  retval->setGuiRepresentation(DoubleProperty::grSLIDER);
429  return retval;
430 
431 }
432 
433 DoublePropertyPtr LevelSetFilter::getRadiusOption(QDomElement root)
434 {
435  DoublePropertyPtr retval =
436  DoubleProperty::initialize("Radius for morphological closing",
437  "",
438  "Select radius (in mm) for the morphological closing of the final result. Radius at 0 will skip this process.",
439  0.0, DoubleRange(0, 20, 0.5), 2, root);
440  retval->setGuiRepresentation(DoubleProperty::grSLIDER);
441  return retval;
442 
443 }
444 
445 } // end namespace cx
446 
cx::LevelSetFilter::getEpsilonOption
DoublePropertyPtr getEpsilonOption(QDomElement root)
Definition: cxLevelSetFilterService.cpp:414
cx::SpaceProviderPtr
boost::shared_ptr< class SpaceProvider > SpaceProviderPtr
Definition: cxLegacySingletons.h:41
qstring_cast
QString qstring_cast(const T &val)
Definition: cxTypeConversions.h:67
cx::FilterImpl::mInputTypes
std::vector< SelectDataStringPropertyBasePtr > mInputTypes
Definition: cxFilterImpl.h:94
cx::reportError
void reportError(QString msg)
Definition: cxLogger.cpp:92
cx::LevelSetFilter::isSeedPointInsideImage
static bool isSeedPointInsideImage(Vector3D, DataPtr)
Definition: cxLevelSetFilterService.cpp:121
cx::Transform3D
Transform3D Transform3D
Transform3D is a representation of an affine 3D transform.
Definition: cxLandmarkPatientRegistrationWidget.h:54
cx::StringPropertySelectData::New
static StringPropertySelectDataPtr New(PatientModelServicePtr patientModelService, QString typeRegexp=".*")
Definition: cxSelectDataStringProperty.h:122
cx::DoubleRange
Utility class for describing a bounded numeric range.
Definition: cxDoubleRange.h:53
cx::ImagePtr
boost::shared_ptr< class Image > ImagePtr
Definition: cxDicomWidget.h:48
cx::LevelSetFilter::getSeedPointFromTool
static Vector3D getSeedPointFromTool(SpaceProviderPtr spaceProvider, DataPtr image)
Definition: cxLevelSetFilterService.cpp:95
cx::Image::getBaseVtkImageData
virtual vtkImageDataPtr getBaseVtkImageData()
Definition: cxImage.cpp:356
cx::LevelSetFilter::LevelSetFilter
LevelSetFilter(ctkPluginContext *pluginContext)
Definition: cxLevelSetFilterService.cpp:73
cx::LevelSetFilter::setActive
virtual void setActive(bool on)
Definition: cxLevelSetFilterService.cpp:352
cx::LevelSetFilter::getHelp
virtual QString getHelp() const
Definition: cxLevelSetFilterService.cpp:88
cx::Image::intitializeFromParentImage
virtual void intitializeFromParentImage(ImagePtr parentImage)
Definition: cxImage.cpp:172
cx::FilterImpl::mOptionsAdapters
std::vector< PropertyPtr > mOptionsAdapters
Definition: cxFilterImpl.h:96
cx::DataPtr
boost::shared_ptr< class Data > DataPtr
Definition: cxRegistrationApplicator.h:43
cx::DataLocations::findConfigFolder
static QString findConfigFolder(QString pathRelativeToConfigRoot, QString alternativeAbsolutePath="")
Definition: cxDataLocations.cpp:260
vtkPolyDataPtr
vtkSmartPointer< class vtkPolyData > vtkPolyDataPtr
Definition: cxAirwaysFilterService.h:52
cx::FilterImpl::mServices
VisServicesPtr mServices
Definition: cxFilterImpl.h:103
cx::SelectDataStringPropertyBasePtr
boost::shared_ptr< class SelectDataStringPropertyBase > SelectDataStringPropertyBasePtr
Definition: cxMeshInfoWidget.h:48
cx::FilterImpl::patientService
PatientModelServicePtr patientService()
Definition: cxFilterImpl.cpp:51
cx::ContourFilter::execute
virtual bool execute()
Definition: cxContourFilter.cpp:215
cx::Image
A volumetric data set.
Definition: cxImage.h:66
cx::getImageSize
int * getImageSize(DataPtr inputImage)
Definition: cxLevelSetFilterService.cpp:115
cx::FilterImpl::setActive
virtual void setActive(bool on)
Definition: cxFilterImpl.cpp:100
cx::FilterImpl::mOptions
QDomElement mOptions
Definition: cxFilterImpl.h:97
cx::LevelSetFilter::getThresholdOption
DoublePropertyPtr getThresholdOption(QDomElement root)
Definition: cxLevelSetFilterService.cpp:405
vtkImageCastPtr
vtkSmartPointer< class vtkImageCast > vtkImageCastPtr
Definition: vtkForwardDeclarations.h:88
cx::spaceProvider
cxLogicManager_EXPORT SpaceProviderPtr spaceProvider()
Definition: cxLogicManager.cpp:64
cx::DoublePropertyPtr
boost::shared_ptr< class DoubleProperty > DoublePropertyPtr
Definition: cxReconstructionMethodService.h:54
cx::Vector3D
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:63
vtkImageImportPtr
vtkSmartPointer< class vtkImageImport > vtkImageImportPtr
Definition: cxVideoConnection.h:42
cx::itkImageType
itk::Image< PixelType, Dimension > itkImageType
Definition: cxAlgorithmHelpers.h:51
cx::FilterImpl::mOutputTypes
std::vector< SelectDataStringPropertyBasePtr > mOutputTypes
Definition: cxFilterImpl.h:95
cx::DoubleProperty::initialize
static DoublePropertyPtr initialize(const QString &uid, QString name, QString help, double value, DoubleRange range, int decimals, QDomNode root=QDomNode())
Definition: cxDoubleProperty.cpp:51
cx::ContourFilter::postProcess
virtual bool postProcess()
Definition: cxContourFilter.cpp:316
cx::AlgorithmHelper::getITKfromSSCImage
static itkImageType::ConstPointer getITKfromSSCImage(ImagePtr image)
Definition: cxAlgorithmHelpers.cpp:56
cx::StringPropertySelectImage::New
static StringPropertySelectImagePtr New(PatientModelServicePtr patientModelService)
Definition: cxSelectDataStringProperty.h:97
cx::LevelSetFilter::getName
virtual QString getName() const
Definition: cxLevelSetFilterService.cpp:78
cx::LevelSetFilter::getType
virtual QString getType() const
Definition: cxLevelSetFilterService.cpp:83
cx::MeshPtr
boost::shared_ptr< class Mesh > MeshPtr
Definition: cxForwardDeclarations.h:68
vtkImageDataPtr
vtkSmartPointer< class vtkImageData > vtkImageDataPtr
Definition: cxVideoConnectionWidget.h:51
cx::LevelSetFilter::getAlphaOption
DoublePropertyPtr getAlphaOption(QDomElement root)
Definition: cxLevelSetFilterService.cpp:423
cx::LevelSetFilter::getRadiusOption
DoublePropertyPtr getRadiusOption(QDomElement root)
Definition: cxLevelSetFilterService.cpp:433