+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+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(output); // not SetNthOutput
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType, class MaskImageType>
+void
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations,
+ MaskImagePointer skeleton,
+ MaskImageIndexType index,
+ MaskImagePixelType label)
+{
+ DD("TrackFromThisIndex");
+ DD(index);
+ DD((int)label);
+ // Create NeighborhoodIterator
+ typedef itk::NeighborhoodIterator<MaskImageType> 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;
+ }
+ }
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType, class MaskImageType>
+bool
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+SearchForTracheaSeed(int skip)
+{
+ if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
+ // Restart the search until a seed is found, skipping more and more slices
+ bool stop = false;
+ while (!stop) {
+ // Search seed (parameters = UpperThresholdForTrachea)
+ static const unsigned int Dim = ImageType::ImageDimension;
+ typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
+ typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
+ typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
+ sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-skip-5;
+ sliceRegionSize[Dim-1]=5;
+ sliceRegion.SetSize(sliceRegionSize);
+ sliceRegion.SetIndex(sliceRegionIndex);
+ typedef itk::ImageRegionConstIterator<ImageType> IteratorType;
+ IteratorType it(working_input, sliceRegion);
+ it.GoToBegin();
+ while (!it.IsAtEnd()) {
+ if(it.Get() < GetUpperThresholdForTrachea() ) {
+ AddSeed(it.GetIndex());
+ // DD(it.GetIndex());
+ }
+ ++it;
+ }
+
+ // if we do not found : restart
+ stop = (m_Seeds.size() != 0);
+ if (!stop) {
+ if (GetVerboseStep()) {
+ std::cout << "\t No seed found this time. I skip some slices and restart." << std::endl;
+ }
+ if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
+ // we want to skip more than a half of the image, it is probably a bug
+ std::cerr << "Number of slices to skip to find trachea to high = " << skip << std::endl;
+ stop = true;
+ }
+ skip += 5;
+ }
+ }
+ }
+ return (m_Seeds.size() != 0);
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType, class MaskImageType>
+void
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+TracheaRegionGrowing()
+{
+ // Explosion controlled region growing
+ typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, InternalImageType> ImageFilterType;
+ typename ImageFilterType::Pointer f= ImageFilterType::New();
+ f->SetInput(working_input);
+ f->SetVerbose(false);
+ f->SetLower(-2000);
+ f->SetUpper(GetUpperThresholdForTrachea());
+ f->SetMinimumLowerThreshold(-2000);
+ f->SetMaximumUpperThreshold(0);
+ f->SetAdaptLowerBorder(false);
+ f->SetAdaptUpperBorder(true);
+ f->SetMinimumSize(5000);
+ f->SetReplaceValue(1);
+ f->SetMultiplier(GetMultiplierForTrachea());
+ f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
+ f->SetMinimumThresholdStepSize(1);
+ for(unsigned int i=0; i<m_Seeds.size();i++) {
+ f->AddSeed(m_Seeds[i]);
+ // DD(m_Seeds[i]);
+ }
+ f->Update();
+
+ // take first (main) connected component
+ writeImage<InternalImageType>(f->GetOutput(), "t1.mhd");
+ trachea_tmp = Labelize<InternalImageType>(f->GetOutput(),
+ GetBackgroundValue(),
+ true,
+ GetMinimalComponentSize());
+ writeImage<InternalImageType>(trachea_tmp, "t2.mhd");
+ trachea_tmp = KeepLabels<InternalImageType>(trachea_tmp,
+ GetBackgroundValue(),
+ GetForegroundValue(),
+ 1, 1, false);
+ writeImage<InternalImageType>(trachea_tmp, "t3.mhd");