#include "itkConnectedComponentImageFilter.h"
#include "itkRelabelComponentImageFilter.h"
#include "itkOtsuThresholdImageFilter.h"
+#include "itkBinaryThinningImageFilter3D.h"
+#include "itkImageIteratorWithIndex.h"
+#include "itkBinaryMorphologicalOpeningImageFilter.h"
+#include "itkBinaryMorphologicalClosingImageFilter.h"
//--------------------------------------------------------------------
-template <class TInputImageType, class TMaskImageType>
-clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
+template <class ImageType, class MaskImageType>
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
ExtractLungFilter():
clitk::FilterBase(),
- itk::ImageToImageFilter<TInputImageType, TMaskImageType>()
+ itk::ImageToImageFilter<ImageType, MaskImageType>()
{
+ SetNumberOfSteps(10);
+ m_MaxSeedNumber = 500;
+
// Default global options
this->SetNumberOfRequiredInputs(2);
SetPatientMaskBackgroundValue(0);
SetBackgroundValue(0); // Must be zero
SetForegroundValue(1);
+ SetMinimalComponentSize(100);
// Step 1 default values
SetUpperThreshold(-300);
SetUpperThresholdForTrachea(-900);
SetMultiplierForTrachea(5);
SetThresholdStepSizeForTrachea(64);
+ SetNumberOfSlicesToSkipBeforeSearchingSeed(0);
// Step 3 default values
SetNumberOfHistogramBins(500);
p3->SetLastKeep(2);
p3->UseLastKeepOff();
SetLabelizeParameters3(p3);
+
+ // Step 5
+ FinalOpenCloseOff();
+ SetFinalOpenCloseRadius(1);
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
-template <class TInputImageType, class TMaskImageType>
+template <class ImageType, class MaskImageType>
void
-clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
-SetInput(const TInputImageType * image)
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+SetInput(const ImageType * image)
{
- this->SetNthInput(0, const_cast<TInputImageType *>(image));
+ this->SetNthInput(0, const_cast<ImageType *>(image));
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
-template <class TInputImageType, class TMaskImageType>
+template <class ImageType, class MaskImageType>
void
-clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
-SetInputPatientMask(TMaskImageType * image, MaskImagePixelType bg )
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+SetInputPatientMask(MaskImageType * image, MaskImagePixelType bg )
{
- this->SetNthInput(1, const_cast<TMaskImageType *>(image));
+ this->SetNthInput(1, const_cast<MaskImageType *>(image));
SetPatientMaskBackgroundValue(bg);
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
-template <class TInputImageType, class TMaskImageType>
+template <class ImageType, class MaskImageType>
void
-clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
AddSeed(InternalIndexType s)
{
m_Seeds.push_back(s);
//--------------------------------------------------------------------
-template <class TInputImageType, class TMaskImageType>
+template <class ImageType, class MaskImageType>
template<class ArgsInfoType>
void
-clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
SetArgsInfo(ArgsInfoType mArgsInfo)
{
SetVerboseOption_GGO(mArgsInfo);
SetUpperThreshold_GGO(mArgsInfo);
SetLowerThreshold_GGO(mArgsInfo);
+ SetNumberOfSlicesToSkipBeforeSearchingSeed_GGO(mArgsInfo);
SetLabelizeParameters1_GGO(mArgsInfo);
if (!mArgsInfo.remove1_given) {
GetLabelizeParameters1()->AddLabelToRemove(2);
SetThresholdStepSizeForTrachea_GGO(mArgsInfo);
AddSeed_GGO(mArgsInfo);
+ SetMinimalComponentSize_GGO(mArgsInfo);
+
SetNumberOfHistogramBins_GGO(mArgsInfo);
SetLabelizeParameters2_GGO(mArgsInfo);
SetRadiusForTrachea_GGO(mArgsInfo);
SetLabelizeParameters3_GGO(mArgsInfo);
+
+ SetFinalOpenCloseRadius_GGO(mArgsInfo);
+ SetFinalOpenClose_GGO(mArgsInfo);
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
-template <class TInputImageType, class TMaskImageType>
+template <class ImageType, class MaskImageType>
void
-clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
GenerateOutputInformation()
{
- input = dynamic_cast<const TInputImageType*>(itk::ProcessObject::GetInput(0));
Superclass::GenerateOutputInformation();
-// MaskImagePointer output = this->GetOutput(0);
// Get input pointers
- input = dynamic_cast<const TInputImageType*>(itk::ProcessObject::GetInput(0));
- patient = dynamic_cast<const TMaskImageType*>(itk::ProcessObject::GetInput(1));
+ patient = dynamic_cast<const MaskImageType*>(itk::ProcessObject::GetInput(1));
+ input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
// Check image
- if (!HaveSameSizeAndSpacing<TInputImageType, TMaskImageType>(input, patient)) {
- this->SetLastError("* ERROR * the images (input and patient mask) must have the same size & spacing");
- return;
+ if (!HaveSameSizeAndSpacing<ImageType, MaskImageType>(input, patient)) {
+ clitkExceptionMacro("the 'input' and 'patient' masks must have the same size & spacing.");
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
StartNewStep("Set background to initial image");
- working_input = SetBackground<TInputImageType, TMaskImageType>
+ working_input = SetBackground<ImageType, MaskImageType>
(input, patient, GetPatientMaskBackgroundValue(), -1000);
- StopCurrentStep<InputImageType>(working_input);
+ StopCurrentStep<ImageType>(working_input);
//--------------------------------------------------------------------
//--------------------------------------------------------------------
StartNewStep("Remove Air");
+ // Check threshold
+ if (m_UseLowerThreshold) {
+ if (m_LowerThreshold > m_UpperThreshold) {
+ clitkExceptionMacro("lower threshold cannot be greater than upper threshold.");
+ }
+ }
// Threshold to get air
- typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
+ typedef itk::BinaryThresholdImageFilter<ImageType, InternalImageType> InputBinarizeFilterType;
typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
binarizeFilter->SetInput(working_input);
if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold);
// Labelize and keep right labels
working_image = Labelize<InternalImageType>(working_image, GetBackgroundValue(), true, GetMinimalComponentSize());
+
working_image = RemoveLabels<InternalImageType>
(working_image, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove());
+
typename InternalImageType::Pointer air = KeepLabels<InternalImageType>
(working_image,
GetBackgroundValue(),
GetLabelizeParameters1()->GetUseLastKeep());
// Set Air to BG
- working_input = SetBackground<TInputImageType, InternalImageType>
+ working_input = SetBackground<ImageType, InternalImageType>
(working_input, air, this->GetForegroundValue(), this->GetBackgroundValue());
// End
- StopCurrentStep<InputImageType>(working_input);
+ StopCurrentStep<ImageType>(working_input);
//--------------------------------------------------------------------
//--------------------------------------------------------------------
- StartNewStep("Find the trachea");
- //DD(m_Seeds.size());
- if (m_Seeds.size() == 0) { // try to find seed
- // Search seed (parameters = UpperThresholdForTrachea)
- static const unsigned int Dim = InputImageType::ImageDimension;
- typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
- typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
- typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
- sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-5;
- sliceRegionSize[Dim-1]=5;
- sliceRegion.SetSize(sliceRegionSize);
- sliceRegion.SetIndex(sliceRegionIndex);
- //DD(GetUpperThresholdForTrachea());
- //DD(sliceRegion);
- typedef itk::ImageRegionConstIterator<InputImageType> IteratorType;
- IteratorType it(working_input, sliceRegion);
- it.GoToBegin();
- while (!it.IsAtEnd()) {
- if(it.Get() < GetUpperThresholdForTrachea() ) {
- AddSeed(it.GetIndex());
- // DD(it.GetIndex());
- }
- ++it;
- }
- }
-
- //DD(m_Seeds.size());
- if (m_Seeds.size() != 0) {
- // Explosion controlled region growing
- typedef clitk::ExplosionControlledThresholdConnectedImageFilter<InputImageType, InternalImageType> ImageFilterType;
- typename ImageFilterType::Pointer f= ImageFilterType::New();
- f->SetInput(working_input);
- f->SetVerbose(false);
- f->SetLower(-2000);
- f->SetUpper(GetUpperThresholdForTrachea());
- f->SetMinimumLowerThreshold(-2000);
- f->SetMaximumUpperThreshold(0);
- f->SetAdaptLowerBorder(false);
- f->SetAdaptUpperBorder(true);
- f->SetMinimumSize(5000);
- f->SetReplaceValue(1);
- f->SetMultiplier(GetMultiplierForTrachea());
- f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
- f->SetMinimumThresholdStepSize(1);
- for(unsigned int i=0; i<m_Seeds.size();i++) {
- // std::cout<<"Adding seed " <<m_Seeds[i]<<"..."<<std::endl;
- f->AddSeed(m_Seeds[i]);
- }
- f->Update();
- trachea_tmp = f->GetOutput();
- // Set output
- StopCurrentStep<InternalImageType>(trachea_tmp);
-
- }
- else { // Trachea not found
- this->SetWarning("* WARNING * No seed found for trachea.");
- // Set output
- StopCurrentStep();
- }
+ StartNewStep("Search for the trachea");
+ SearchForTrachea();
//--------------------------------------------------------------------
//--------------------------------------------------------------------
StartNewStep("Extract the lung with Otsu filter");
// Automated Otsu thresholding and relabeling
- typedef itk::OtsuThresholdImageFilter<InputImageType,InternalImageType> OtsuThresholdImageFilterType;
+ typedef itk::OtsuThresholdImageFilter<ImageType,InternalImageType> OtsuThresholdImageFilterType;
typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
otsuFilter->SetInput(working_input);
otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins());
working_image = SetBackground<InternalImageType, InternalImageType>
(working_image, trachea_tmp, 1, -1);
- // Dilate the trachea
- static const unsigned int Dim = InputImageType::ImageDimension;
+ // Dilate the trachea
+ static const unsigned int Dim = ImageType::ImageDimension;
typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
KernelType structuringElement;
structuringElement.SetRadius(GetRadiusForTrachea());
// Set output
StopCurrentStep<InternalImageType>(working_image);
}
-
//--------------------------------------------------------------------
//--------------------------------------------------------------------
working_image = cropFilter2->GetOutput();
StopCurrentStep<InternalImageType>(working_image);
+ //--------------------------------------------------------------------
+ //--------------------------------------------------------------------
+ // Final OpenClose
+ if (GetFinalOpenClose()) {
+ StartNewStep("Open/Close");
+
+ // Structuring element
+ typedef itk::BinaryBallStructuringElement<InternalPixelType, ImageDimension> KernelType;
+ KernelType structuringElement;
+ structuringElement.SetRadius(GetFinalOpenCloseRadius());
+ structuringElement.CreateStructuringElement();
+
+ // Open
+ typedef itk::BinaryMorphologicalOpeningImageFilter<InternalImageType, InternalImageType, KernelType> OpenFilterType;
+ typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
+ openFilter->SetInput(working_image);
+ openFilter->SetBackgroundValue(GetBackgroundValue());
+ openFilter->SetForegroundValue(GetForegroundValue());
+ openFilter->SetKernel(structuringElement);
+
+ // Close
+ typedef itk::BinaryMorphologicalClosingImageFilter<InternalImageType, InternalImageType, KernelType> CloseFilterType;
+ typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
+ closeFilter->SetInput(openFilter->GetOutput());
+ closeFilter->SetSafeBorder(true);
+ closeFilter->SetForegroundValue(GetForegroundValue());
+ closeFilter->SetKernel(structuringElement);
+ closeFilter->Update();
+ working_image = closeFilter->GetOutput();
+ }
+
//--------------------------------------------------------------------
//--------------------------------------------------------------------
StartNewStep("Separate Left/Right lungs");
working_image = statisticsImageFilter->GetOutput();
// Decompose the first label
- static const unsigned int Dim = InputImageType::ImageDimension;
+ static const unsigned int Dim = ImageType::ImageDimension;
if (initialNumberOfLabels<2) {
// Structuring element radius
- typename InputImageType::SizeType radius;
+ typename ImageType::SizeType radius;
for (unsigned int i=0;i<Dim;i++) radius[i]=1;
typedef clitk::DecomposeAndReconstructImageFilter<InternalImageType,InternalImageType> DecomposeAndReconstructFilterType;
typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
working_image = decomposeAndReconstructFilter->GetOutput();
}
- // Retain labels (lungs)
+ // Retain labels ('1' is largset lung, so right. '2' is left)
typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
thresholdFilter->SetInput(working_image);
thresholdFilter->Update();
working_image = thresholdFilter->GetOutput();
StopCurrentStep<InternalImageType> (working_image);
+
+ // Final Cast
+ StartNewStep("Cast the lung mask");
+ typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
+ typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
+ caster->SetInput(working_image);
+ caster->Update();
+ output = caster->GetOutput();
+
+ // Update output info
+ this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
+
+
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
-template <class TInputImageType, class TMaskImageType>
+template <class ImageType, class MaskImageType>
void
-clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
-GenerateData() {
- // Final Cast
- typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
- typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
- caster->SetInput(working_image);
- caster->Update();
- // Set output
- //this->SetNthOutput(0, caster->GetOutput()); // -> no because redo filter otherwise
- this->GraftOutput(caster->GetOutput());
- return;
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+GenerateData()
+{
+ // Set the output
+ this->GraftOutput(output); // not SetNthOutput
}
//--------------------------------------------------------------------
+
+//--------------------------------------------------------------------
+template <class ImageType, class MaskImageType>
+bool
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+SearchForTracheaSeed(int skip)
+{
+ if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
+ // Restart the search until a seed is found, skipping more and more slices
+ bool stop = false;
+ while (!stop) {
+ // Search seed (parameters = UpperThresholdForTrachea)
+ static const unsigned int Dim = ImageType::ImageDimension;
+ typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
+ typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
+ typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
+ sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-skip-5;
+ sliceRegionSize[Dim-1]=5;
+ sliceRegion.SetSize(sliceRegionSize);
+ sliceRegion.SetIndex(sliceRegionIndex);
+ typedef itk::ImageRegionConstIterator<ImageType> IteratorType;
+ IteratorType it(working_input, sliceRegion);
+ it.GoToBegin();
+ while (!it.IsAtEnd()) {
+ if(it.Get() < GetUpperThresholdForTrachea() ) {
+ AddSeed(it.GetIndex());
+ // DD(it.GetIndex());
+ }
+ ++it;
+ }
+
+ // if we do not found : restart
+ stop = (m_Seeds.size() != 0);
+ if (!stop) {
+ if (GetVerboseStep()) {
+ std::cout << "\t No seed found this time. I skip some slices and restart." << std::endl;
+ }
+ 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 << "Number of slices to skip to find trachea to high = " << skip << std::endl;
+ stop = true;
+ }
+ skip += 5;
+ }
+ }
+ }
+ return (m_Seeds.size() != 0);
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType, class MaskImageType>
+void
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+TracheaRegionGrowing()
+{
+ // Explosion controlled region growing
+ typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, InternalImageType> ImageFilterType;
+ typename ImageFilterType::Pointer f= ImageFilterType::New();
+ f->SetInput(working_input);
+ f->SetVerbose(false);
+ f->SetLower(-2000);
+ f->SetUpper(GetUpperThresholdForTrachea());
+ f->SetMinimumLowerThreshold(-2000);
+ f->SetMaximumUpperThreshold(0);
+ f->SetAdaptLowerBorder(false);
+ f->SetAdaptUpperBorder(true);
+ f->SetMinimumSize(5000);
+ f->SetReplaceValue(1);
+ f->SetMultiplier(GetMultiplierForTrachea());
+ f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
+ f->SetMinimumThresholdStepSize(1);
+ for(unsigned int i=0; i<m_Seeds.size();i++) {
+ f->AddSeed(m_Seeds[i]);
+ // DD(m_Seeds[i]);
+ }
+ f->Update();
+
+ // take first (main) connected component
+ trachea_tmp = Labelize<InternalImageType>(f->GetOutput(),
+ GetBackgroundValue(),
+ true,
+ GetMinimalComponentSize());
+ trachea_tmp = KeepLabels<InternalImageType>(trachea_tmp,
+ GetBackgroundValue(),
+ GetForegroundValue(),
+ 1, 1, false);
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType, class MaskImageType>
+double
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+ComputeTracheaVolume()
+{
+ typedef itk::ImageRegionConstIterator<InternalImageType> IteratorType;
+ IteratorType iter(trachea_tmp, trachea_tmp->GetLargestPossibleRegion());
+ iter.GoToBegin();
+ double volume = 0.0;
+ while (!iter.IsAtEnd()) {
+ if (iter.Get() == this->GetForegroundValue()) volume++;
+ ++iter;
+ }
+
+ double voxelsize = trachea_tmp->GetSpacing()[0]*trachea_tmp->GetSpacing()[1]*trachea_tmp->GetSpacing()[2];
+ return volume*voxelsize;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType, class MaskImageType>
+void
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+SearchForTrachea()
+{
+ // Search for seed among n slices, skip some slices before starting
+ // if not found -> skip more and restart
+ // when seed found : perform region growing
+ // compute trachea volume
+ // if volume not plausible -> skip more slices and restart
+
+ bool stop = false;
+ double volume = 0.0;
+ int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed();
+ while (!stop) {
+ stop = SearchForTracheaSeed(skip);
+ if (stop) {
+ TracheaRegionGrowing();
+ volume = ComputeTracheaVolume()/1000; // assume mm3, so divide by 1000 to get cc
+ if ((volume > 10) && (volume < 55 )) { // it is ok
+ // Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ]
+ if (GetVerboseStep()) {
+ std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
+ }
+ stop = true;
+ }
+ else {
+ if (GetVerboseStep()) {
+ std::cout << "\t The volume of the trachea (" << volume
+ << " cc) seems not correct. I skip some slices (" << skip << ") and restart to find seeds."
+ << std::endl;
+ }
+ skip += 5;
+ stop = false;
+ // empty the list of seed
+ m_Seeds.clear();
+ }
+ }
+ }
+
+ if (volume != 0.0) {
+ // Set output
+ StopCurrentStep<InternalImageType>(trachea_tmp);
+ }
+ else { // Trachea not found
+ this->SetWarning("* WARNING * No seed found for trachea.");
+ StopCurrentStep();
+ }
+}
+//--------------------------------------------------------------------
+
#endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX
git://git.creatis.insa-lyon.fr / clitk.git / blobdiff