#include "itkImageIteratorWithIndex.h"
#include "itkBinaryMorphologicalOpeningImageFilter.h"
#include "itkBinaryMorphologicalClosingImageFilter.h"
+#include "itkConstantPadImageFilter.h"
//--------------------------------------------------------------------
template <class ImageType>
// Step 6
FillHolesFlagOn();
+ AutoCropOn();
}
//--------------------------------------------------------------------
StartNewStep("Set background to initial image");
working_input = SetBackground<ImageType, MaskImageType>
(working_input, patient, GetPatientMaskBackgroundValue(), -1000, true);
+
+ // Pad images with air to prevent patient touching the image border
+ static const unsigned int Dim = ImageType::ImageDimension;
+ typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
+ typename PadFilterType::Pointer padFilter = PadFilterType::New();
+ padFilter->SetInput(working_input);
+ padFilter->SetConstant(-1000);
+ typename ImageType::SizeType bounds;
+ for (unsigned i = 0; i < Dim - 1; ++i)
+ bounds[i] = 1;
+ bounds[Dim - 1] = 0;
+ padFilter->SetPadLowerBound(bounds);
+ padFilter->SetPadUpperBound(bounds);
+ padFilter->Update();
+ working_input = padFilter->GetOutput();
+
StopCurrentStep<ImageType>(working_input);
PrintMemory(GetVerboseMemoryFlag(), "After set bg"); // OK, additional mem = 0
//--------------------------------------------------------------------
PrintMemory(GetVerboseMemoryFlag(), "before autocropfilter");
if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
- trachea = clitk::AutoCrop<MaskImageType>(trachea, GetBackgroundValue());
+ if (GetAutoCrop())
+ trachea = clitk::AutoCrop<MaskImageType>(trachea, GetBackgroundValue());
+ else
+ {
+ // Remove Padding region
+ typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
+ typename CropFilterType::Pointer cropFilter = CropFilterType::New();
+ cropFilter->SetInput(trachea);
+ cropFilter->SetLowerBoundaryCropSize(bounds);
+ cropFilter->SetUpperBoundaryCropSize(bounds);
+ cropFilter->Update();
+ trachea = cropFilter->GetOutput();
+ }
StopCurrentStep<MaskImageType>(trachea);
PrintMemory(GetVerboseMemoryFlag(), "after delete trachea");
}
//--------------------------------------------------------------------
StartNewStep("Cropping lung");
PrintMemory(GetVerboseMemoryFlag(), "Before Autocropfilter");
- working_mask = clitk::AutoCrop<MaskImageType>(working_mask, GetBackgroundValue());
+ if (GetAutoCrop())
+ working_mask = clitk::AutoCrop<MaskImageType>(working_mask, GetBackgroundValue());
+ else
+ {
+ // Remove Padding region
+ typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
+ typename CropFilterType::Pointer cropFilter = CropFilterType::New();
+ cropFilter->SetInput(working_mask);
+ cropFilter->SetLowerBoundaryCropSize(bounds);
+ cropFilter->SetUpperBoundaryCropSize(bounds);
+ cropFilter->Update();
+ working_mask = cropFilter->GetOutput();
+ }
StopCurrentStep<MaskImageType>(working_mask);
//--------------------------------------------------------------------
StopCurrentStep<MaskImageType>(working_mask);
}
- //--------------------------------------------------------------------
- //--------------------------------------------------------------------
- StartNewStep("Separate Left/Right lungs");
- PrintMemory(GetVerboseMemoryFlag(), "Before Separate");
- // Initial label
- working_mask = Labelize<MaskImageType>(working_mask,
- GetBackgroundValue(),
- false,
- GetMinimalComponentSize());
-
- PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
-
- // Count the labels
- typedef itk::StatisticsImageFilter<MaskImageType> StatisticsImageFilterType;
- typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
- statisticsImageFilter->SetInput(working_mask);
- statisticsImageFilter->Update();
- unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
- working_mask = statisticsImageFilter->GetOutput();
+ if (GetSeparateLungsFlag()) {
+ //--------------------------------------------------------------------
+ //--------------------------------------------------------------------
+ StartNewStep("Separate Left/Right lungs");
+ PrintMemory(GetVerboseMemoryFlag(), "Before Separate");
+ // Initial label
+ working_mask = Labelize<MaskImageType>(working_mask,
+ GetBackgroundValue(),
+ false,
+ GetMinimalComponentSize());
+
+ PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
+
+ // Count the labels
+ typedef itk::StatisticsImageFilter<MaskImageType> StatisticsImageFilterType;
+ typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
+ statisticsImageFilter->SetInput(working_mask);
+ statisticsImageFilter->Update();
+ unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
+ working_mask = statisticsImageFilter->GetOutput();
+
+ PrintMemory(GetVerboseMemoryFlag(), "After count label");
- PrintMemory(GetVerboseMemoryFlag(), "After count label");
-
- // Decompose the first label
- static const unsigned int Dim = ImageType::ImageDimension;
- if (initialNumberOfLabels<2) {
- // Structuring element radius
- typename ImageType::SizeType radius;
- for (unsigned int i=0;i<Dim;i++) radius[i]=1;
- typedef clitk::DecomposeAndReconstructImageFilter<MaskImageType,MaskImageType> DecomposeAndReconstructFilterType;
- typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
- decomposeAndReconstructFilter->SetInput(working_mask);
- decomposeAndReconstructFilter->SetVerbose(false);
- decomposeAndReconstructFilter->SetRadius(radius);
- decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
- decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
- decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
- decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
- decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
- decomposeAndReconstructFilter->SetFullyConnected(true);
- decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
- decomposeAndReconstructFilter->Update();
- working_mask = decomposeAndReconstructFilter->GetOutput();
+ // Decompose the first label
+ if (initialNumberOfLabels<2) {
+ // Structuring element radius
+ typename ImageType::SizeType radius;
+ for (unsigned int i=0;i<Dim;i++) radius[i]=1;
+ typedef clitk::DecomposeAndReconstructImageFilter<MaskImageType,MaskImageType> DecomposeAndReconstructFilterType;
+ typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
+ decomposeAndReconstructFilter->SetInput(working_mask);
+ decomposeAndReconstructFilter->SetVerbose(false);
+ decomposeAndReconstructFilter->SetRadius(radius);
+ decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
+ decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
+ decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
+ decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
+ decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
+ decomposeAndReconstructFilter->SetFullyConnected(true);
+ decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
+ decomposeAndReconstructFilter->Update();
+ working_mask = decomposeAndReconstructFilter->GetOutput();
+ }
+ PrintMemory(GetVerboseMemoryFlag(), "After decomposeAndReconstructFilter");
}
- PrintMemory(GetVerboseMemoryFlag(), "After decomposeAndReconstructFilter");
// Retain labels ('1' is largset lung, so right. '2' is left)
typedef itk::ThresholdImageFilter<MaskImageType> ThresholdImageFilterType;
//--------------------------------------------------------------------
template <class ImageType>
-void
+void
clitk::ExtractLungFilter<ImageType>::
SearchForTrachea()
{
double volume = 0.0;
int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed();
while (!stop) {
- stop = SearchForTracheaSeed(skip);
- if (stop) {
+ stop = true;
+ if (SearchForTracheaSeed(skip)) {
TracheaRegionGrowing();
volume = ComputeTracheaVolume()/1000; // assume mm3, so divide by 1000 to get cc
if (GetWriteStepFlag()) {
if (GetVerboseStepFlag()) {
std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
}
- stop = true;
}
else {
if (GetVerboseStepFlag()) {
// empty the list of seed
m_Seeds.clear();
}
+ if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
+ // we want to skip more than a half of the image, it is probably a bug
+ std::cerr << "2 : Number of slices to skip to find trachea too high = " << skip << std::endl;
+ stop = true;
+ }
}
else {
stop = true;