+//--------------------------------------------------------------------
+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");
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType, class MaskImageType>
+double
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+ComputeTracheaVolume()
+{
+ typedef itk::ImageRegionConstIterator<InternalImageType> IteratorType;
+ IteratorType iter(trachea_tmp, trachea_tmp->GetLargestPossibleRegion());
+ iter.GoToBegin();
+ double volume = 0.0;
+ while (!iter.IsAtEnd()) {
+ if (iter.Get() == this->GetForegroundValue()) volume++;
+ ++iter;
+ }
+
+ double voxelsize = trachea_tmp->GetSpacing()[0]*trachea_tmp->GetSpacing()[1]*trachea_tmp->GetSpacing()[2];
+ return volume*voxelsize;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType, class MaskImageType>
+void
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+SearchForTrachea()
+{
+ // Search for seed among n slices, skip some slices before starting
+ // if not found -> skip more and restart
+
+ // when seed found : perform region growing
+ // compute trachea volume
+ // if volume not plausible -> skip more slices and restart
+
+ bool stop = false;
+ double volume = 0.0;
+ int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed();
+ while (!stop) {
+ stop = SearchForTracheaSeed(skip);
+ if (stop) {
+ TracheaRegionGrowing();
+ volume = ComputeTracheaVolume(); // assume mm3
+ if ((volume > 10000) && (volume < 55000 )) { // it is ok
+ // Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ]
+ if (GetVerboseStep()) {
+ std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
+ }
+ stop = true;
+ }
+ else {
+ if (GetVerboseStep()) {
+ std::cout << "\t The volume of the trachea (" << volume
+ << " cc) seems not correct. I skip some slices (" << skip << ") and restart to find seeds."
+ << std::endl;
+ }
+ skip += 5;
+ stop = false;
+ // empty the list of seed
+ m_Seeds.clear();
+ }
+ }
+ }
+
+ if (volume != 0.0) {
+ // Set output
+ StopCurrentStep<InternalImageType>(trachea_tmp);
+ }
+ else { // Trachea not found
+ this->SetWarning("* WARNING * No seed found for trachea.");
+ StopCurrentStep();
+ }
+}
+//--------------------------------------------------------------------
+