From b676e0ed11db2ee1f2c16e76b35bf9fdf8be90a1 Mon Sep 17 00:00:00 2001 From: dsarrut Date: Wed, 8 Sep 2010 13:05:43 +0000 Subject: [PATCH] extract trachea skeleton --- segmentation/clitkExtractAirwayTreeInfo.cxx | 44 ++++ segmentation/clitkExtractAirwayTreeInfo.ggo | 18 ++ .../clitkExtractAirwayTreeInfoFilter.h | 152 ++++++++++++ .../clitkExtractAirwayTreeInfoFilter.txx | 227 ++++++++++++++++++ .../clitkExtractAirwayTreeInfoGenericFilter.h | 70 ++++++ ...litkExtractAirwayTreeInfoGenericFilter.txx | 100 ++++++++ 6 files changed, 611 insertions(+) create mode 100644 segmentation/clitkExtractAirwayTreeInfo.cxx create mode 100644 segmentation/clitkExtractAirwayTreeInfo.ggo create mode 100644 segmentation/clitkExtractAirwayTreeInfoFilter.h create mode 100644 segmentation/clitkExtractAirwayTreeInfoFilter.txx create mode 100644 segmentation/clitkExtractAirwayTreeInfoGenericFilter.h create mode 100644 segmentation/clitkExtractAirwayTreeInfoGenericFilter.txx diff --git a/segmentation/clitkExtractAirwayTreeInfo.cxx b/segmentation/clitkExtractAirwayTreeInfo.cxx new file mode 100644 index 0000000..e5a60c4 --- /dev/null +++ b/segmentation/clitkExtractAirwayTreeInfo.cxx @@ -0,0 +1,44 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html +======================================================================-====*/ + +// clitk +#include "clitkExtractAirwayTreeInfo_ggo.h" +#include "clitkExtractAirwayTreeInfoGenericFilter.h" + +//-------------------------------------------------------------------- +int main(int argc, char * argv[]) +{ + + // Init command line + GGO(clitkExtractAirwayTreeInfo, args_info); + CLITK_INIT; + + // Filter + typedef clitk::ExtractAirwayTreeInfoGenericFilter FilterType; + FilterType::Pointer filter = FilterType::New(); + + filter->SetArgsInfo(args_info); + filter->Update(); + + if (filter->HasError()) { + std::cout << filter->GetLastError() << std::endl; + } + + return EXIT_SUCCESS; +} // This is the end, my friend +//-------------------------------------------------------------------- diff --git a/segmentation/clitkExtractAirwayTreeInfo.ggo b/segmentation/clitkExtractAirwayTreeInfo.ggo new file mode 100644 index 0000000..75fe817 --- /dev/null +++ b/segmentation/clitkExtractAirwayTreeInfo.ggo @@ -0,0 +1,18 @@ +#File clitkExtractAirwayTreeInfos.ggo +package "clitkExtractAirwayTreeInfos" +version "1.0" +purpose "From trachea binary image, extract skeleton and carena point." + +option "config" - "Config file" string no +option "imagetypes" - "Display allowed image types" flag off +option "verbose" v "Verbose" flag off +option "verboseStep" - "Verbose each step" flag off +option "writeStep" w "Write image at each step" flag off +option "verboseOption" - "Display options values" flag off +option "verboseWarningOff" - "Do not display warning" flag off + +section "I/O" + +option "input" i "Input trachea mask image filename" string yes +option "output" o "Output skeleton mask image filename" string yes + diff --git a/segmentation/clitkExtractAirwayTreeInfoFilter.h b/segmentation/clitkExtractAirwayTreeInfoFilter.h new file mode 100644 index 0000000..f537311 --- /dev/null +++ b/segmentation/clitkExtractAirwayTreeInfoFilter.h @@ -0,0 +1,152 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html + ======================================================================-====*/ + +#ifndef CLITKEXTRACTAIRWAYTREEINFOSFILTER_H +#define CLITKEXTRACTAIRWAYTREEINFOSFILTER_H + +// clitk +#include "clitkFilterBase.h" +#include "clitkDecomposeAndReconstructImageFilter.h" +#include "clitkExplosionControlledThresholdConnectedImageFilter.h" +#include "clitkSegmentationUtils.h" + +// itk +#include "itkStatisticsImageFilter.h" + +namespace clitk { + + //-------------------------------------------------------------------- + /* + From a trachea binary image, compute the skeleton and track the + path to find the carena. + */ + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + class Bifurcation + { + public: + Bifurcation(IndexType _index, PixelType _l, PixelType _l1, PixelType _l2) { + index = _index; + _l = l; + _l1 = l1; + _l2 = l2; + } + IndexType index; + PixelType l; + PixelType l1; + PixelType l2; + }; + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + class ITK_EXPORT ExtractAirwayTreeInfoFilter: + public clitk::FilterBase, + public itk::ImageToImageFilter + { + + public: + /** Standard class typedefs. */ + typedef itk::ImageToImageFilter Superclass; + typedef ExtractAirwayTreeInfoFilter Self; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(ExtractAirwayTreeInfoFilter, ImageToImageFilter); + FILTERBASE_INIT; + + /** Some convenient typedefs */ + typedef TImageType ImageType; + typedef typename ImageType::ConstPointer ImageConstPointer; + typedef typename ImageType::Pointer ImagePointer; + typedef typename ImageType::RegionType ImageRegionType; + typedef typename ImageType::PixelType ImagePixelType; + typedef typename ImageType::SizeType ImageSizeType; + typedef typename ImageType::IndexType ImageIndexType; + typedef typename ImageType::PointType ImagePointType; + + itkStaticConstMacro(ImageDimension, unsigned int, ImageType::ImageDimension); + typedef int InternalPixelType; + typedef itk::Image InternalImageType; + typedef typename InternalImageType::Pointer InternalImagePointer; + typedef typename InternalImageType::IndexType InternalIndexType; + typedef LabelizeParameters LabelParamType; + + /** Connect inputs */ + void SetInput(const ImageType * image); + + // Set all options at a time + template + void SetArgsInfo(ArgsInfoType arg); + + // Background / Foreground + itkGetConstMacro(BackgroundValue, ImagePixelType); + itkGetConstMacro(ForegroundValue, ImagePixelType); + + // Get results + itkGetConstMacro(FirstTracheaPoint, ImagePointType); + + protected: + ExtractAirwayTreeInfoFilter(); + virtual ~ExtractAirwayTreeInfoFilter() {} + + // Main members + ImageConstPointer input; + ImagePointer skeleton; + ImagePointer working_input; + + // Global options + itkSetMacro(BackgroundValue, ImagePixelType); + itkSetMacro(ForegroundValue, ImagePixelType); + ImagePixelType m_BackgroundValue; + ImagePixelType m_ForegroundValue; + + // Results + ImagePointType m_FirstTracheaPoint; + + virtual void GenerateOutputInformation(); + virtual void GenerateData(); + + typedef Bifurcation BifurcationType; + void TrackFromThisIndex(std::vector & listOfBifurcations, + ImagePointer skeleton, + ImageIndexType index, + ImagePixelType label); + private: + ExtractAirwayTreeInfoFilter(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + }; // end class + //-------------------------------------------------------------------- + +} // end namespace clitk +//-------------------------------------------------------------------- + +#ifndef ITK_MANUAL_INSTANTIATION +#include "clitkExtractAirwayTreeInfoFilter.txx" +#endif + +#endif diff --git a/segmentation/clitkExtractAirwayTreeInfoFilter.txx b/segmentation/clitkExtractAirwayTreeInfoFilter.txx new file mode 100644 index 0000000..ba99b3d --- /dev/null +++ b/segmentation/clitkExtractAirwayTreeInfoFilter.txx @@ -0,0 +1,227 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html + ======================================================================-====*/ + +#ifndef CLITKEXTRACTAIRWAYTREEINFOSFILTER_TXX +#define CLITKEXTRACTAIRWAYTREEINFOSFILTER_TXX + +// clitk +#include "clitkImageCommon.h" +#include "clitkSetBackgroundImageFilter.h" +#include "clitkSegmentationUtils.h" +#include "clitkAutoCropFilter.h" + +// itk +#include "itkBinaryThresholdImageFilter.h" +#include "itkConnectedComponentImageFilter.h" +#include "itkRelabelComponentImageFilter.h" +#include "itkOtsuThresholdImageFilter.h" +#include "itkBinaryThinningImageFilter3D.h" +#include "itkImageIteratorWithIndex.h" + + +//-------------------------------------------------------------------- +template +clitk::ExtractAirwayTreeInfoFilter:: +ExtractAirwayTreeInfoFilter(): + clitk::FilterBase(), + itk::ImageToImageFilter() +{ + // Default global options + this->SetNumberOfRequiredInputs(1); + SetBackgroundValue(0); + SetForegroundValue(1); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractAirwayTreeInfoFilter:: +SetInput(const ImageType * image) +{ + this->SetNthInput(0, const_cast(image)); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +template +void +clitk::ExtractAirwayTreeInfoFilter:: +SetArgsInfo(ArgsInfoType mArgsInfo) +{ + SetVerboseOption_GGO(mArgsInfo); + SetVerboseStep_GGO(mArgsInfo); + SetWriteStep_GGO(mArgsInfo); + SetVerboseWarningOff_GGO(mArgsInfo); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractAirwayTreeInfoFilter:: +GenerateOutputInformation() +{ + Superclass::GenerateOutputInformation(); + //this->GetOutput(0)->SetRequestedRegion(this->GetOutput(0)->GetLargestPossibleRegion()); + + // Get input pointer + input = dynamic_cast(itk::ProcessObject::GetInput(0)); + + StartNewStepOrStop("Find bronchial bifurcations"); + // Step 1 : extract skeleton + // Define the thinning filter + typedef itk::BinaryThinningImageFilter3D ThinningFilterType; + typename ThinningFilterType::Pointer thinningFilter = ThinningFilterType::New(); + thinningFilter->SetInput(input); + thinningFilter->Update(); + skeleton = thinningFilter->GetOutput(); + writeImage(skeleton, "skeleton.mhd"); + + // Step 2 : tracking + DD("tracking"); + + // Step 2.1 : find first point for tracking + typedef itk::ImageRegionConstIteratorWithIndex 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 ImageType::IndexType index = it.GetIndex(); + skeleton->TransformIndexToPhysicalPoint(index, m_FirstTracheaPoint); + DD(index); + DD(m_FirstTracheaPoint); + + // Step 2.2 : initialize neighborhooditerator + typedef itk::NeighborhoodIterator 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 ImageType::PixelType label = GetForegroundValue()+1; + while ((label == GetBackgroundValue()) || (label == GetForegroundValue())) { label++; } + DD(label); + + // Track from the first point + std::vector listOfBifurcations; + TrackFromThisIndex(listOfBifurcations, skeleton, index, label); + DD("end track"); + DD(listOfBifurcations.size()); + writeImage(skeleton, "skeleton2.mhd"); + + for(unsigned int i=0; iTransformIndexToPhysicalPoint(listOfBifurcations[i].index, p); + DD(p); + } + +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractAirwayTreeInfoFilter:: +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 + this->GraftOutput(skeleton); // not SetNthOutput +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractAirwayTreeInfoFilter:: +TrackFromThisIndex(std::vector & listOfBifurcations, + ImagePointer skeleton, + ImageIndexType index, + ImagePixelType label) +{ + DD("TrackFromThisIndex"); + DD(index); + DD((int)label); + // Create NeighborhoodIterator + typedef itk::NeighborhoodIterator NeighborhoodIteratorType; + typename NeighborhoodIteratorType::SizeType radius; + radius.Fill(1); + NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion()); + + // Track + std::vector listOfTrackedPoint; + bool stop = false; + while (!stop) { + nit.SetLocation(index); + // DD((int)nit.GetCenterPixel()); + nit.SetCenterPixel(label); + listOfTrackedPoint.clear(); + for(unsigned int i=0; i 2) { + std::cerr << "too much bifurcation points ... ?" << std::endl; + exit(0); + } + // Else this it the end of the tracking + } + stop = true; + } + } +} +//-------------------------------------------------------------------- + + +#endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX diff --git a/segmentation/clitkExtractAirwayTreeInfoGenericFilter.h b/segmentation/clitkExtractAirwayTreeInfoGenericFilter.h new file mode 100644 index 0000000..a2f51aa --- /dev/null +++ b/segmentation/clitkExtractAirwayTreeInfoGenericFilter.h @@ -0,0 +1,70 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html +======================================================================-====*/ + +#ifndef CLITKEXTRACTAIRWAYTREEINFOSGENERICFILTER_H +#define CLITKEXTRACTAIRWAYTREEINFOSGENERICFILTER_H + +#include "clitkIO.h" +#include "clitkImageToImageGenericFilter.h" +#include "clitkExtractAirwayTreeInfoFilter.h" + +//-------------------------------------------------------------------- +namespace clitk +{ + + template + class ITK_EXPORT ExtractAirwayTreeInfoGenericFilter: + public ImageToImageGenericFilter > + { + + public: + //-------------------------------------------------------------------- + ExtractAirwayTreeInfoGenericFilter(); + + //-------------------------------------------------------------------- + typedef ExtractAirwayTreeInfoGenericFilter Self; + typedef ImageToImageGenericFilter > Superclass; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer ConstPointer; + + //-------------------------------------------------------------------- + itkNewMacro(Self); + itkTypeMacro(ExtractAirwayTreeInfoGenericFilter, LightObject); + + //-------------------------------------------------------------------- + void SetArgsInfo(const ArgsInfoType & a); + + //-------------------------------------------------------------------- + // Main function called each time the filter is updated + template + void UpdateWithInputImageType(); + + protected: + template void InitializeImageType(); + ArgsInfoType mArgsInfo; + + }; // end class + //-------------------------------------------------------------------- + +} // end namespace clitk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "clitkExtractAirwayTreeInfoGenericFilter.txx" +#endif + +#endif // #define CLITKEXTRACTAIRWAYTREEINFOSGENERICFILTER_H diff --git a/segmentation/clitkExtractAirwayTreeInfoGenericFilter.txx b/segmentation/clitkExtractAirwayTreeInfoGenericFilter.txx new file mode 100644 index 0000000..47511e4 --- /dev/null +++ b/segmentation/clitkExtractAirwayTreeInfoGenericFilter.txx @@ -0,0 +1,100 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html + ======================================================================-====*/ + +#ifndef CLITKEXTRACTAIRWAYTREEINFOSGENERICFILTER_TXX +#define CLITKEXTRACTAIRWAYTREEINFOSGENERICFILTER_TXX + +#include "clitkImageCommon.h" + +//-------------------------------------------------------------------- +template +clitk::ExtractAirwayTreeInfoGenericFilter::ExtractAirwayTreeInfoGenericFilter(): + ImageToImageGenericFilter("ExtractAirwayTreeInfo") +{ + this->SetFilterBase(NULL); + // Default values + cmdline_parser_clitkExtractAirwayTreeInfo_init(&mArgsInfo); + InitializeImageType<3>(); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +template +void clitk::ExtractAirwayTreeInfoGenericFilter::InitializeImageType() +{ + ADD_IMAGE_TYPE(Dim, uchar); + // ADD_IMAGE_TYPE(Dim, int); + // ADD_IMAGE_TYPE(Dim, float); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void clitk::ExtractAirwayTreeInfoGenericFilter::SetArgsInfo(const ArgsInfoType & a) +{ + mArgsInfo=a; + SetIOVerbose(mArgsInfo.verbose_flag); + if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes(); + if (mArgsInfo.input_given) AddInputFilename(mArgsInfo.input_arg); + if (mArgsInfo.output_given) AddOutputFilename(mArgsInfo.output_arg); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +// Update with the number of dimensions and the pixeltype +//-------------------------------------------------------------------- +template +template +void clitk::ExtractAirwayTreeInfoGenericFilter::UpdateWithInputImageType() +{ + // Reading input + typename ImageType::Pointer input = this->template GetInput(0); + + // Create filter + typedef clitk::ExtractAirwayTreeInfoFilter FilterType; + typename FilterType::Pointer filter = FilterType::New(); + + // Set the filter (needed for example for threaded monitoring) + this->SetFilterBase(filter); + + // Set global Options + filter->SetStopOnError(this->GetStopOnError()); + filter->SetArgsInfo(mArgsInfo); + filter->SetInput(input); + + // Go ! + filter->Update(); + + // Check if error + if (filter->HasError()) { + SetLastError(filter->GetLastError()); + // No output + return; + } + + // Write/Save results + typename ImageType::Pointer output = filter->GetOutput(); + this->template SetNextOutput(output); +} +//-------------------------------------------------------------------- + +#endif //#define CLITKEXTRACTAIRWAYTREEINFOSGENERICFILTER_TXX -- 2.47.1