- // Try to extract bifurcation in the trachea (bronchi)
- // STILL EXPERIMENTAL
- 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);
- }