- if (GetFindBronchialBifurcations()) {
- StartNewStep("Find bronchial bifurcations");
- // Step 1 : extract skeleton
- typedef itk::BinaryThinningImageFilter3D<MaskImageType, MaskImageType> ThinningFilterType;
- typename ThinningFilterType::Pointer thinningFilter = ThinningFilterType::New();
- thinningFilter->SetInput(trachea);
- thinningFilter->Update();
- typename MaskImageType::Pointer skeleton = thinningFilter->GetOutput();
-
- // 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()) {
- clitkExceptionMacro("first point in the trachea skeleton not found.");
- }
- typename MaskImageType::IndexType index = it.GetIndex();
-
- // Step 2.2 : initialize neighborhooditerator
- typedef itk::NeighborhoodIterator<MaskImageType> NeighborhoodIteratorType;
- typename NeighborhoodIteratorType::SizeType radius;
- radius.Fill(1);
- NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
-
- // Find first label number (must be different from BG and FG)
- typename MaskImageType::PixelType label = GetForegroundValue()+1;
- while ((label == GetBackgroundValue()) || (label == GetForegroundValue())) { label++; }
-
- // Track from the first point
- std::vector<BifurcationType> listOfBifurcations;
- m_SkeletonTree.set_head(index);
- TrackFromThisIndex(listOfBifurcations, skeleton, index, label, m_SkeletonTree.begin());
- DD("end track");
- DD(listOfBifurcations.size());
- DD(m_SkeletonTree.size());
-
- for(unsigned int i=0; i<listOfBifurcations.size(); i++) {
- skeleton->TransformIndexToPhysicalPoint(listOfBifurcations[i].index,
- listOfBifurcations[i].point);
- }
-
- // Search for the first slice that separate the bronchus (carena)
- typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
- typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
- extractSliceFilter->SetInput(trachea);
- extractSliceFilter->SetDirection(2);
- extractSliceFilter->Update();
- typedef typename ExtractSliceFilterType::SliceType SliceType;
- std::vector<typename SliceType::Pointer> mInputSlices;
- extractSliceFilter->GetOutputSlices(mInputSlices);
-
- bool stop = false;
- DD(listOfBifurcations[0].index);
- DD(listOfBifurcations[1].index);
- int slice_index = listOfBifurcations[0].index[2]; // first slice from carena in skeleton
- int i=0;
- TreeIterator firstIter = m_SkeletonTree.child(listOfBifurcations[1].treeIter, 0);
- TreeIterator secondIter = m_SkeletonTree.child(listOfBifurcations[1].treeIter, 1);
- typename SliceType::IndexType in1;
- typename SliceType::IndexType in2;
- while (!stop) {
- DD(slice_index);
-
- // Labelize the current slice
- typename SliceType::Pointer temp = Labelize<SliceType>(mInputSlices[slice_index],
- GetBackgroundValue(),
- true,
- GetMinimalComponentSize());
- // Check the value of the two skeleton points;
- in1[0] = (*firstIter)[0];
- in1[1] = (*firstIter)[1];
- typename SliceType::PixelType v1 = temp->GetPixel(in1);
- DD(in1);
- DD((int)v1);
- in2[0] = (*secondIter)[0];
- in2[1] = (*secondIter)[1];
- typename SliceType::PixelType v2 = temp->GetPixel(in2);
- DD(in2);
- DD((int)v2);
-
- // TODO IF NOT FOUND ????
-
- if (v1 != v2) {
- stop = true;
- }
- else {
- i++;
- --slice_index;
- ++firstIter;
- ++secondIter;
- }
- }
- MaskImageIndexType carena_index;
- carena_index[0] = lrint(in2[0] + in1[0])/2.0;
- carena_index[1] = lrint(in2[1] + in1[1])/2.0;
- carena_index[2] = slice_index;
- MaskImagePointType carena_position;
- DD(carena_index);
- skeleton->TransformIndexToPhysicalPoint(carena_index,
- carena_position);
- DD(carena_position);
-
- // Set and save Carina position
- StartNewStep("Save carina position");
- // Try to load the DB
- try {
- LoadAFDB();
- } catch (clitk::ExceptionObject e) {
- // Do nothing if not found, it will be used anyway to write the result
- }
- GetAFDB()->SetPoint3D("Carena", carena_position);
- }
- }