-/*=========================================================================
- 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));
-
- StartNewStep("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()) {
- throw clitk::ExceptionObject("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);
- }
-
- // The first bifurcation is supposed to be the carina
- skeleton->TransformIndexToPhysicalPoint(listOfBifurcations[0].index, m_CarinaPoint);
- DD(listOfBifurcations[0].index);
- DD(m_CarinaPoint);
-}
-//--------------------------------------------------------------------
-
-
-//--------------------------------------------------------------------
-template <class ImageType>
-void
-clitk::ExtractAirwayTreeInfoFilter<ImageType>::
-GenerateData()
-{
- // 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
git://git.creatis.insa-lyon.fr / clitk.git / commitdiff