--- /dev/null
+/*=========================================================================
+ 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<args_info_clitkExtractAirwayTreeInfo> 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
+//--------------------------------------------------------------------
--- /dev/null
+#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
+
--- /dev/null
+/*=========================================================================
+ 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 IndexType, class PixelType>
+ 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 TImageType>
+ class ITK_EXPORT ExtractAirwayTreeInfoFilter:
+ public clitk::FilterBase,
+ public itk::ImageToImageFilter<TImageType, TImageType>
+ {
+
+ public:
+ /** Standard class typedefs. */
+ typedef itk::ImageToImageFilter<TImageType, TImageType> Superclass;
+ typedef ExtractAirwayTreeInfoFilter Self;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> 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<InternalPixelType, ImageType::ImageDimension> InternalImageType;
+ typedef typename InternalImageType::Pointer InternalImagePointer;
+ typedef typename InternalImageType::IndexType InternalIndexType;
+ typedef LabelizeParameters<InternalPixelType> LabelParamType;
+
+ /** Connect inputs */
+ void SetInput(const ImageType * image);
+
+ // Set all options at a time
+ template<class ArgsInfoType>
+ 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<ImageIndexType,ImagePixelType> BifurcationType;
+ void TrackFromThisIndex(std::vector<BifurcationType> & 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
--- /dev/null
+/*=========================================================================
+ 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 <class ImageType>
+clitk::ExtractAirwayTreeInfoFilter<ImageType>::
+ExtractAirwayTreeInfoFilter():
+ clitk::FilterBase(),
+ itk::ImageToImageFilter<ImageType, ImageType>()
+{
+ // Default global options
+ this->SetNumberOfRequiredInputs(1);
+ SetBackgroundValue(0);
+ SetForegroundValue(1);
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::ExtractAirwayTreeInfoFilter<ImageType>::
+SetInput(const ImageType * image)
+{
+ this->SetNthInput(0, const_cast<ImageType *>(image));
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+template<class ArgsInfoType>
+void
+clitk::ExtractAirwayTreeInfoFilter<ImageType>::
+SetArgsInfo(ArgsInfoType mArgsInfo)
+{
+ SetVerboseOption_GGO(mArgsInfo);
+ SetVerboseStep_GGO(mArgsInfo);
+ SetWriteStep_GGO(mArgsInfo);
+ SetVerboseWarningOff_GGO(mArgsInfo);
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::ExtractAirwayTreeInfoFilter<ImageType>::
+GenerateOutputInformation()
+{
+ Superclass::GenerateOutputInformation();
+ //this->GetOutput(0)->SetRequestedRegion(this->GetOutput(0)->GetLargestPossibleRegion());
+
+ // Get input pointer
+ input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
+
+ StartNewStepOrStop("Find bronchial bifurcations");
+ // Step 1 : extract skeleton
+ // Define the thinning filter
+ typedef itk::BinaryThinningImageFilter3D<ImageType, ImageType> ThinningFilterType;
+ typename ThinningFilterType::Pointer thinningFilter = ThinningFilterType::New();
+ thinningFilter->SetInput(input);
+ thinningFilter->Update();
+ skeleton = thinningFilter->GetOutput();
+ writeImage<ImageType>(skeleton, "skeleton.mhd");
+
+ // Step 2 : tracking
+ DD("tracking");
+
+ // Step 2.1 : find first point for tracking
+ typedef itk::ImageRegionConstIteratorWithIndex<ImageType> 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<ImageType> 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<BifurcationType> listOfBifurcations;
+ TrackFromThisIndex(listOfBifurcations, skeleton, index, label);
+ DD("end track");
+ DD(listOfBifurcations.size());
+ writeImage<ImageType>(skeleton, "skeleton2.mhd");
+
+ for(unsigned int i=0; i<listOfBifurcations.size(); i++) {
+ typename ImageType::PointType p;
+ skeleton->TransformIndexToPhysicalPoint(listOfBifurcations[i].index, p);
+ DD(p);
+ }
+
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::ExtractAirwayTreeInfoFilter<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
+ this->GraftOutput(skeleton); // not SetNthOutput
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::ExtractAirwayTreeInfoFilter<ImageType>::
+TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations,
+ ImagePointer skeleton,
+ ImageIndexType index,
+ ImagePixelType label)
+{
+ DD("TrackFromThisIndex");
+ DD(index);
+ DD((int)label);
+ // Create NeighborhoodIterator
+ typedef itk::NeighborhoodIterator<ImageType> 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;
+ }
+ }
+}
+//--------------------------------------------------------------------
+
+
+#endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX
--- /dev/null
+/*=========================================================================
+ 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 ArgsInfoType>
+ class ITK_EXPORT ExtractAirwayTreeInfoGenericFilter:
+ public ImageToImageGenericFilter<ExtractAirwayTreeInfoGenericFilter<ArgsInfoType> >
+ {
+
+ public:
+ //--------------------------------------------------------------------
+ ExtractAirwayTreeInfoGenericFilter();
+
+ //--------------------------------------------------------------------
+ typedef ExtractAirwayTreeInfoGenericFilter Self;
+ typedef ImageToImageGenericFilter<ExtractAirwayTreeInfoGenericFilter<ArgsInfoType> > Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ //--------------------------------------------------------------------
+ itkNewMacro(Self);
+ itkTypeMacro(ExtractAirwayTreeInfoGenericFilter, LightObject);
+
+ //--------------------------------------------------------------------
+ void SetArgsInfo(const ArgsInfoType & a);
+
+ //--------------------------------------------------------------------
+ // Main function called each time the filter is updated
+ template<class ImageType>
+ void UpdateWithInputImageType();
+
+ protected:
+ template<unsigned int Dim> void InitializeImageType();
+ ArgsInfoType mArgsInfo;
+
+ }; // end class
+ //--------------------------------------------------------------------
+
+} // end namespace clitk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "clitkExtractAirwayTreeInfoGenericFilter.txx"
+#endif
+
+#endif // #define CLITKEXTRACTAIRWAYTREEINFOSGENERICFILTER_H
--- /dev/null
+/*=========================================================================
+ 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<class ArgsInfoType>
+clitk::ExtractAirwayTreeInfoGenericFilter<ArgsInfoType>::ExtractAirwayTreeInfoGenericFilter():
+ ImageToImageGenericFilter<Self>("ExtractAirwayTreeInfo")
+{
+ this->SetFilterBase(NULL);
+ // Default values
+ cmdline_parser_clitkExtractAirwayTreeInfo_init(&mArgsInfo);
+ InitializeImageType<3>();
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<class ArgsInfoType>
+template<unsigned int Dim>
+void clitk::ExtractAirwayTreeInfoGenericFilter<ArgsInfoType>::InitializeImageType()
+{
+ ADD_IMAGE_TYPE(Dim, uchar);
+ // ADD_IMAGE_TYPE(Dim, int);
+ // ADD_IMAGE_TYPE(Dim, float);
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<class ArgsInfoType>
+void clitk::ExtractAirwayTreeInfoGenericFilter<ArgsInfoType>::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<class ArgsInfoType>
+template<class ImageType>
+void clitk::ExtractAirwayTreeInfoGenericFilter<ArgsInfoType>::UpdateWithInputImageType()
+{
+ // Reading input
+ typename ImageType::Pointer input = this->template GetInput<ImageType>(0);
+
+ // Create filter
+ typedef clitk::ExtractAirwayTreeInfoFilter<ImageType> 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<ImageType>(output);
+}
+//--------------------------------------------------------------------
+
+#endif //#define CLITKEXTRACTAIRWAYTREEINFOSGENERICFILTER_TXX