#include "clitkSetBackgroundImageFilter.h"
#include "clitkSegmentationUtils.h"
#include "clitkAutoCropFilter.h"
+#include "clitkCropLikeImageFilter.h"
+#include "clitkFillMaskFilter.h"
// itk
#include "itkBinaryThresholdImageFilter.h"
#include "itkOtsuThresholdImageFilter.h"
#include "itkBinaryThinningImageFilter3D.h"
#include "itkImageIteratorWithIndex.h"
-
+#include "itkBinaryMorphologicalOpeningImageFilter.h"
+#include "itkBinaryMorphologicalClosingImageFilter.h"
//--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
+template <class ImageType>
+clitk::ExtractLungFilter<ImageType>::
ExtractLungFilter():
clitk::FilterBase(),
+ clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
itk::ImageToImageFilter<ImageType, MaskImageType>()
{
SetNumberOfSteps(10);
m_MaxSeedNumber = 500;
// Default global options
- this->SetNumberOfRequiredInputs(2);
+ this->SetNumberOfRequiredInputs(1);
SetPatientMaskBackgroundValue(0);
SetBackgroundValue(0); // Must be zero
SetForegroundValue(1);
p3->UseLastKeepOff();
SetLabelizeParameters3(p3);
- // Step 5 : find bronchial bifurcations
- FindBronchialBifurcationsOn();
+ // Step 5
+ OpenCloseOff();
+ SetOpenCloseRadius(1);
+
+ // Step 6
+ FillHolesOn();
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
+template <class ImageType>
void
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
+clitk::ExtractLungFilter<ImageType>::
SetInput(const ImageType * image)
{
this->SetNthInput(0, const_cast<ImageType *>(image));
//--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
-void
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
-SetInputPatientMask(MaskImageType * image, MaskImagePixelType bg )
-{
- this->SetNthInput(1, const_cast<MaskImageType *>(image));
- SetPatientMaskBackgroundValue(bg);
-}
-//--------------------------------------------------------------------
-
-
-//--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
+template <class ImageType>
void
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
+clitk::ExtractLungFilter<ImageType>::
AddSeed(InternalIndexType s)
{
m_Seeds.push_back(s);
//--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
+template <class ImageType>
template<class ArgsInfoType>
void
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
+clitk::ExtractLungFilter<ImageType>::
SetArgsInfo(ArgsInfoType mArgsInfo)
{
SetVerboseOption_GGO(mArgsInfo);
SetWriteStep_GGO(mArgsInfo);
SetVerboseWarningOff_GGO(mArgsInfo);
+ SetAFDBFilename_GGO(mArgsInfo);
+ SetOutputLungFilename_GGO(mArgsInfo);
+ SetOutputTracheaFilename_GGO(mArgsInfo);
+
SetUpperThreshold_GGO(mArgsInfo);
SetLowerThreshold_GGO(mArgsInfo);
SetNumberOfSlicesToSkipBeforeSearchingSeed_GGO(mArgsInfo);
SetRadiusForTrachea_GGO(mArgsInfo);
SetLabelizeParameters3_GGO(mArgsInfo);
+
+ SetOpenCloseRadius_GGO(mArgsInfo);
+ SetOpenClose_GGO(mArgsInfo);
+
+ SetFillHoles_GGO(mArgsInfo);
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
+template <class ImageType>
void
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
+clitk::ExtractLungFilter<ImageType>::
GenerateOutputInformation()
{
Superclass::GenerateOutputInformation();
- //this->GetOutput(0)->SetRequestedRegion(this->GetOutput(0)->GetLargestPossibleRegion());
+
+ // Read DB
+ LoadAFDB();
// Get input pointers
- patient = dynamic_cast<const MaskImageType*>(itk::ProcessObject::GetInput(1));
input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
+ patient = GetAFDB()->template GetImage <MaskImageType>("patient");
- // Check image
- if (!HaveSameSizeAndSpacing<ImageType, MaskImageType>(input, patient)) {
- this->SetLastError("* ERROR * the images (input and patient mask) must have the same size & spacing");
- return;
- }
-
//--------------------------------------------------------------------
//--------------------------------------------------------------------
- StartNewStepOrStop("Set background to initial image");
+ // Crop input like patient image (must have the same spacing)
+ StartNewStep("Crop input image to 'patient' extends");
+ typedef clitk::CropLikeImageFilter<ImageType> CropImageFilter;
+ typename CropImageFilter::Pointer cropFilter = CropImageFilter::New();
+ cropFilter->SetInput(input);
+ cropFilter->SetCropLikeImage(patient);
+ cropFilter->Update();
+ working_input = cropFilter->GetOutput();
+ DD(working_input->GetLargestPossibleRegion());
+ StopCurrentStep<ImageType>(working_input);
+
+ //--------------------------------------------------------------------
+ //--------------------------------------------------------------------
+ StartNewStep("Set background to initial image");
working_input = SetBackground<ImageType, MaskImageType>
- (input, patient, GetPatientMaskBackgroundValue(), -1000);
+ (working_input, patient, GetPatientMaskBackgroundValue(), -1000);
StopCurrentStep<ImageType>(working_input);
//--------------------------------------------------------------------
//--------------------------------------------------------------------
- StartNewStepOrStop("Remove Air");
+ StartNewStep("Remove Air");
// Check threshold
if (m_UseLowerThreshold) {
if (m_LowerThreshold > m_UpperThreshold) {
- this->SetLastError("ERROR: lower threshold cannot be greater than upper threshold.");
- return;
+ clitkExceptionMacro("lower threshold cannot be greater than upper threshold.");
}
}
// Threshold to get air
//--------------------------------------------------------------------
//--------------------------------------------------------------------
- StartNewStepOrStop("Search for the trachea");
+ StartNewStep("Search for the trachea");
SearchForTrachea();
//--------------------------------------------------------------------
//--------------------------------------------------------------------
- StartNewStepOrStop("Extract the lung with Otsu filter");
+ StartNewStep("Extract the lung with Otsu filter");
// Automated Otsu thresholding and relabeling
typedef itk::OtsuThresholdImageFilter<ImageType,InternalImageType> OtsuThresholdImageFilterType;
typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
//--------------------------------------------------------------------
//--------------------------------------------------------------------
- StartNewStepOrStop("Select labels");
+ StartNewStep("Select labels");
// Keep right labels
working_image = LabelizeAndSelectLabels<InternalImageType>
(working_image,
//--------------------------------------------------------------------
//--------------------------------------------------------------------
if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
- StartNewStepOrStop("Remove the trachea");
+ StartNewStep("Remove the trachea");
// Set the trachea
working_image = SetBackground<InternalImageType, InternalImageType>
(working_image, trachea_tmp, 1, -1);
//--------------------------------------------------------------------
//--------------------------------------------------------------------
- typedef clitk::AutoCropFilter<InternalImageType> CropFilterType;
- typename CropFilterType::Pointer cropFilter = CropFilterType::New();
+ typedef clitk::AutoCropFilter<InternalImageType> AutoCropFilterType;
+ typename AutoCropFilterType::Pointer autocropFilter = AutoCropFilterType::New();
if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
- StartNewStepOrStop("Croping trachea");
- cropFilter->SetInput(trachea_tmp);
- cropFilter->Update(); // Needed
+ StartNewStep("Cropping trachea");
+ autocropFilter->SetInput(trachea_tmp);
+ autocropFilter->Update(); // Needed
typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
- caster->SetInput(cropFilter->GetOutput());
+ caster->SetInput(autocropFilter->GetOutput());
caster->Update();
trachea = caster->GetOutput();
StopCurrentStep<MaskImageType>(trachea);
//--------------------------------------------------------------------
//--------------------------------------------------------------------
- StartNewStepOrStop("Croping lung");
- typename CropFilterType::Pointer cropFilter2 = CropFilterType::New(); // Needed to reset pipeline
- cropFilter2->SetInput(working_image);
- cropFilter2->Update();
- working_image = cropFilter2->GetOutput();
+ StartNewStep("Cropping lung");
+ typename AutoCropFilterType::Pointer autocropFilter2 = AutoCropFilterType::New(); // Needed to reset pipeline
+ autocropFilter2->SetInput(working_image);
+ autocropFilter2->Update();
+ working_image = autocropFilter2->GetOutput();
+ DD(working_image->GetLargestPossibleRegion());
StopCurrentStep<InternalImageType>(working_image);
//--------------------------------------------------------------------
//--------------------------------------------------------------------
- StartNewStepOrStop("Separate Left/Right lungs");
+ // Final OpenClose
+ if (GetOpenClose()) {
+ StartNewStep("Open/Close");
+
+ // Structuring element
+ typedef itk::BinaryBallStructuringElement<InternalPixelType, ImageDimension> KernelType;
+ KernelType structuringElement;
+ structuringElement.SetRadius(GetOpenCloseRadius());
+ 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();
+ }
+
+ //--------------------------------------------------------------------
+ //--------------------------------------------------------------------
+ // Fill Lungs
+ if (GetFillHoles()) {
+ StartNewStep("Fill Holes");
+ typedef clitk::FillMaskFilter<InternalImageType> FillMaskFilterType;
+ typename FillMaskFilterType::Pointer fillMaskFilter = FillMaskFilterType::New();
+ fillMaskFilter->SetInput(working_image);
+ fillMaskFilter->AddDirection(2);
+ //fillMaskFilter->AddDirection(1);
+ fillMaskFilter->Update();
+ working_image = fillMaskFilter->GetOutput();
+ StopCurrentStep<InternalImageType>(working_image);
+ }
+
+ //--------------------------------------------------------------------
+ //--------------------------------------------------------------------
+ StartNewStep("Separate Left/Right lungs");
// Initial label
working_image = Labelize<InternalImageType>(working_image,
GetBackgroundValue(),
// Decompose the first label
static const unsigned int Dim = ImageType::ImageDimension;
if (initialNumberOfLabels<2) {
+ DD(initialNumberOfLabels);
// Structuring element 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();
decomposeAndReconstructFilter->SetInput(working_image);
- decomposeAndReconstructFilter->SetVerbose(false);
+ decomposeAndReconstructFilter->SetVerbose(true);
decomposeAndReconstructFilter->SetRadius(radius);
decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
StopCurrentStep<InternalImageType> (working_image);
// Final Cast
- StartNewStepOrStop("Final cast");
+ StartNewStep("Cast the lung mask");
typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
caster->SetInput(working_image);
// Update output info
this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
-
- // Try to extract bifurcation in the trachea (bronchi)
- // STILL EXPERIMENTAL
- if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
-
- if (GetFindBronchialBifurcations()) {
- StartNewStepOrStop("Find bronchial bifurcations");
- // Step 1 : extract skeleton
- // Define the thinning filter
- typedef itk::BinaryThinningImageFilter3D<MaskImageType, MaskImageType> ThinningFilterType;
- typename ThinningFilterType::Pointer thinningFilter = ThinningFilterType::New();
- thinningFilter->SetInput(trachea);
- thinningFilter->Update();
- typename MaskImageType::Pointer skeleton = thinningFilter->GetOutput();
- writeImage<MaskImageType>(skeleton, "skeleton.mhd");
-
- // Step 2 : tracking
- DD("tracking");
-
- // Step 2.1 : find first point for tracking
- typedef itk::ImageRegionConstIteratorWithIndex<MaskImageType> IteratorType;
- IteratorType it(skeleton, skeleton->GetLargestPossibleRegion());
- it.GoToReverseBegin();
- while ((!it.IsAtEnd()) && (it.Get() == GetBackgroundValue())) {
- --it;
- }
- if (it.IsAtEnd()) {
- this->SetLastError("ERROR: first point in the skeleton not found ! Abort");
- return;
- }
- DD(skeleton->GetLargestPossibleRegion().GetIndex());
- typename MaskImageType::IndexType index = it.GetIndex();
- DD(index);
-
- // Step 2.2 : initialize neighborhooditerator
- typedef itk::NeighborhoodIterator<MaskImageType> NeighborhoodIteratorType;
- typename NeighborhoodIteratorType::SizeType radius;
- radius.Fill(1);
- NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
- DD(nit.GetSize());
- DD(nit.Size());
-
- // Find first label number (must be different from BG and FG)
- typename MaskImageType::PixelType label = GetForegroundValue()+1;
- while ((label == GetBackgroundValue()) || (label == GetForegroundValue())) { label++; }
- DD(label);
-
- // Track from the first point
- std::vector<BifurcationType> listOfBifurcations;
- TrackFromThisIndex(listOfBifurcations, skeleton, index, label);
- DD("end track");
- DD(listOfBifurcations.size());
- writeImage<MaskImageType>(skeleton, "skeleton2.mhd");
-
- for(unsigned int i=0; i<listOfBifurcations.size(); i++) {
- typename MaskImageType::PointType p;
- skeleton->TransformIndexToPhysicalPoint(listOfBifurcations[i].index, p);
- DD(p);
- }
-
- }
- }
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
+template <class ImageType>
void
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
+clitk::ExtractLungFilter<ImageType>::
GenerateData()
{
- // Do not put some "startnewstep" here, because the object if
- // modified and the filter's pipeline it do two times. But it is
- // required to quit if MustStop was set before.
- if (GetMustStop()) return;
-
- // If everything goes well, set the output
+ // Set the output
this->GraftOutput(output); // not SetNthOutput
+ // Store image filenames into AFDB
+ GetAFDB()->SetImageFilename("lungs", this->GetOutputLungFilename());
+ GetAFDB()->SetImageFilename("trachea", this->GetOutputTracheaFilename());
+ WriteAFDB();
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
-void
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
-TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations,
- MaskImagePointer skeleton,
- MaskImageIndexType index,
- MaskImagePixelType label)
-{
- DD("TrackFromThisIndex");
- DD(index);
- DD((int)label);
- // Create NeighborhoodIterator
- typedef itk::NeighborhoodIterator<MaskImageType> NeighborhoodIteratorType;
- typename NeighborhoodIteratorType::SizeType radius;
- radius.Fill(1);
- NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
-
- // Track
- std::vector<typename NeighborhoodIteratorType::IndexType> listOfTrackedPoint;
- bool stop = false;
- while (!stop) {
- nit.SetLocation(index);
- // DD((int)nit.GetCenterPixel());
- nit.SetCenterPixel(label);
- listOfTrackedPoint.clear();
- for(unsigned int i=0; i<nit.Size(); i++) {
- if (i != nit.GetCenterNeighborhoodIndex ()) { // Do not observe the current point
- // DD(nit.GetIndex(i));
- if (nit.GetPixel(i) == GetForegroundValue()) { // if this is foreground, we continue the tracking
- // DD(nit.GetIndex(i));
- listOfTrackedPoint.push_back(nit.GetIndex(i));
- }
- }
- }
- // DD(listOfTrackedPoint.size());
- if (listOfTrackedPoint.size() == 1) {
- index = listOfTrackedPoint[0];
- }
- else {
- if (listOfTrackedPoint.size() == 2) {
- BifurcationType bif(index, label, label+1, label+2);
- listOfBifurcations.push_back(bif);
- TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1);
- TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2);
- }
- else {
- if (listOfTrackedPoint.size() > 2) {
- std::cerr << "too much bifurcation points ... ?" << std::endl;
- exit(0);
- }
- // Else this it the end of the tracking
- }
- stop = true;
- }
- }
-}
-//--------------------------------------------------------------------
-
-
-//--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
+template <class ImageType>
bool
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
+clitk::ExtractLungFilter<ImageType>::
SearchForTracheaSeed(int skip)
{
if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
//--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
+template <class ImageType>
void
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
+clitk::ExtractLungFilter<ImageType>::
TracheaRegionGrowing()
{
// Explosion controlled region growing
f->SetLower(-2000);
f->SetUpper(GetUpperThresholdForTrachea());
f->SetMinimumLowerThreshold(-2000);
- f->SetMaximumUpperThreshold(0);
+ // f->SetMaximumUpperThreshold(0); // MAYBE TO CHANGE ???
+ f->SetMaximumUpperThreshold(-800); // MAYBE TO CHANGE ???
f->SetAdaptLowerBorder(false);
f->SetAdaptUpperBorder(true);
f->SetMinimumSize(5000);
f->SetMultiplier(GetMultiplierForTrachea());
f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
f->SetMinimumThresholdStepSize(1);
+ f->VerboseOn();
for(unsigned int i=0; i<m_Seeds.size();i++) {
f->AddSeed(m_Seeds[i]);
// DD(m_Seeds[i]);
}
f->Update();
+ writeImage<InternalImageType>(f->GetOutput(), "trg.mhd");
+
// take first (main) connected component
- writeImage<InternalImageType>(f->GetOutput(), "t1.mhd");
trachea_tmp = Labelize<InternalImageType>(f->GetOutput(),
GetBackgroundValue(),
true,
GetMinimalComponentSize());
- writeImage<InternalImageType>(trachea_tmp, "t2.mhd");
trachea_tmp = KeepLabels<InternalImageType>(trachea_tmp,
GetBackgroundValue(),
GetForegroundValue(),
1, 1, false);
- writeImage<InternalImageType>(trachea_tmp, "t3.mhd");
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
+template <class ImageType>
double
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
+clitk::ExtractLungFilter<ImageType>::
ComputeTracheaVolume()
{
typedef itk::ImageRegionConstIterator<InternalImageType> IteratorType;
//--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
+template <class ImageType>
void
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
+clitk::ExtractLungFilter<ImageType>::
SearchForTrachea()
{
// Search for seed among n slices, skip some slices before starting
stop = SearchForTracheaSeed(skip);
if (stop) {
TracheaRegionGrowing();
- volume = ComputeTracheaVolume(); // assume mm3
- if ((volume > 10000) && (volume < 55000 )) { // it is ok
+ 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;
git://git.creatis.insa-lyon.fr / clitk.git / blobdiff