+ typedef InternalIndexType SeedType;
+ SeedType trachea_centre, shape_centre, max_e_centre, prev_e_centre;
+ typedef std::list<SeedType> PointListType;
+ typedef std::list<PointListType> SequenceListType;
+ PointListType* current_sequence = NULL;
+ SequenceListType sequence_list;
+
+ prev_e_centre.Fill(0);
+ std::ostringstream file_name;
+ index[0] = index[1] = 0;
+ size[0] = size[1] = 512;
+ size[2] = 0;
+ while (nslices--) {
+ index[2] = start;
+ start += increment;
+
+ region.SetIndex(index);
+ region.SetSize(size);
+ slicer->SetExtractionRegion(region);
+ slicer->Update();
+ label_filter->SetInput(slicer->GetOutput());
+ label_filter->Update();
+
+ label_to_map_filter->SetInput(label_filter->GetOutput());
+ label_to_map_filter->Update();
+ typename LabelImageType::Pointer label_map = label_to_map_filter->GetOutput();
+
+ if (GetWriteStepFlag()) {
+ map_to_label_filter->SetInput(label_map);
+ writer->SetInput(map_to_label_filter->GetOutput());
+ file_name.str("");
+ file_name << "labels_";
+ file_name.width(3);
+ file_name.fill('0');
+ file_name << index[2] << ".mhd";
+ writer->SetFileName(file_name.str().c_str());
+ writer->Update();
+ }
+
+ typename ShapeLabelType::Pointer shape, max_e_shape;
+ double max_elongation = GetMaxElongation();
+ double max_size = size[0]*size[1]/128;
+ double max_e = 0;
+ int nshapes = 0;
+ unsigned int nlables = label_map->GetNumberOfLabelObjects();
+ for (unsigned int j = 0; j < nlables; j++) {
+ shape = label_map->GetNthLabelObject(j);
+ if (shape->Size() > 150 && shape->Size() <= max_size) {
+#if ITK_VERSION_MAJOR < 4
+ double e = 1 - 1/shape->GetBinaryElongation();
+#else
+ double e = 1 - 1/shape->GetElongation();
+#endif
+ //double area = 1 - r->Area() ;
+ if (e < max_elongation) {
+ nshapes++;
+ shape_centre[0] = (shape->GetCentroid()[0] - origin[0])/spacing[0];
+ shape_centre[1] = (shape->GetCentroid()[1] - origin[1])/spacing[1];
+ shape_centre[2] = index[2];
+ //double d = 1 - (shape_centre - image_centre).Magnitude()/max_dist;
+ double dx = shape_centre[0] - image_centre[0];
+ double d = 1 - dx*2/size[0];
+ e = e + d;
+ if (e > max_e)
+ {
+ max_e = e;
+ max_e_shape = shape;
+ max_e_centre = shape_centre;
+ }
+ }
+ }
+ }
+
+ if (nshapes > 0)
+ {
+ itk::Point<typename SeedType::IndexValueType, ImageType::ImageDimension> p1, p2;
+ p1[0] = max_e_centre[0];
+ p1[1] = max_e_centre[1];
+ p1[2] = max_e_centre[2];
+
+ p2[0] = prev_e_centre[0];
+ p2[1] = prev_e_centre[1];
+ p2[2] = prev_e_centre[2];
+
+ double mag = (p2 - p1).GetNorm();
+ if (GetVerboseRegionGrowingFlag()) {
+ cout.precision(3);
+ cout << index[2] << ": ";
+ cout << "region(" << max_e_centre[0] << ", " << max_e_centre[1] << ", " << max_e_centre[2] << "); ";
+ cout << "prev_region(" << prev_e_centre[0] << ", " << prev_e_centre[1] << ", " << prev_e_centre[2] << "); ";
+ cout << "mag(" << mag << "); " << endl;
+ }
+
+ if (mag > 5)
+ {
+ PointListType point_list;
+ point_list.push_back(max_e_centre);
+ sequence_list.push_back(point_list);
+ current_sequence = &sequence_list.back();
+ }
+ else if (current_sequence)
+ current_sequence->push_back(max_e_centre);
+
+ prev_e_centre= max_e_centre;
+ }
+ else {
+ if (GetVerboseRegionGrowingFlag()) {
+ cout << "No shapes found at slice " << index[2] << std::endl;
+ }
+ }
+ }
+
+ size_t longest = 0;
+ for (typename SequenceListType::iterator s = sequence_list.begin(); s != sequence_list.end(); s++)
+ {
+ if (s->size() > longest)
+ {
+ longest = s->size();
+ trachea_centre = s->front();
+ }
+ }
+
+ if (longest > 0) {
+ if (GetVerboseRegionGrowingFlag())
+ std::cout << "seed at: " << trachea_centre << std::endl;
+ m_Seeds.push_back(trachea_centre);
+ }
+ }
+
+ return (m_Seeds.size() != 0);
+}
+//--------------------------------------------------------------------
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::ExtractLungFilter<ImageType>::
+TracheaRegionGrowing()
+{
+ // Explosion controlled region growing
+ PrintMemory(GetVerboseMemoryFlag(), "Before ExplosionControlledThresholdConnectedImageFilter");
+ typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, MaskImageType> ImageFilterType;
+ typename ImageFilterType::Pointer f= ImageFilterType::New();
+ f->SetInput(working_input);
+ f->SetLower(-2000);
+ f->SetUpper(GetUpperThresholdForTrachea());
+ f->SetMinimumLowerThreshold(-2000);
+ // f->SetMaximumUpperThreshold(0); // MAYBE TO CHANGE ???
+ f->SetMaximumUpperThreshold(-700); // MAYBE TO CHANGE ???
+ f->SetAdaptLowerBorder(false);
+ f->SetAdaptUpperBorder(true);
+ f->SetMinimumSize(5000);
+ f->SetReplaceValue(1);
+ f->SetMultiplier(GetMultiplierForTrachea());
+ f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
+ f->SetMinimumThresholdStepSize(1);
+ f->SetVerbose(GetVerboseRegionGrowingFlag());
+ for(unsigned int i=0; i<m_Seeds.size();i++) {
+ f->AddSeed(m_Seeds[i]);
+ // DD(m_Seeds[i]);
+ }
+ f->Update();
+ PrintMemory(GetVerboseMemoryFlag(), "After RG update");
+
+ // take first (main) connected component
+ trachea = Labelize<MaskImageType>(f->GetOutput(),
+ GetBackgroundValue(),
+ true,
+ 1000);//GetMinimalComponentSize());
+ PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
+ trachea = KeepLabels<MaskImageType>(trachea,
+ GetBackgroundValue(),
+ GetForegroundValue(),
+ 1, 1, false);
+ PrintMemory(GetVerboseMemoryFlag(), "After KeepLabels");
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+double
+clitk::ExtractLungFilter<ImageType>::
+ComputeTracheaVolume()
+{
+ typedef itk::ImageRegionConstIterator<InternalImageType> IteratorType;
+ IteratorType iter(trachea, trachea->GetLargestPossibleRegion());
+ iter.GoToBegin();
+ double volume = 0.0;
+ while (!iter.IsAtEnd()) {
+ if (iter.Get() == this->GetForegroundValue()) volume++;
+ ++iter;
+ }
+
+ double voxelsize = trachea->GetSpacing()[0]*trachea->GetSpacing()[1]*trachea->GetSpacing()[2];
+ return volume*voxelsize;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::ExtractLungFilter<ImageType>::
+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 has_seed;
+ bool stop = false;
+ double volume = 0.0;
+ int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed();
+ while (!stop) {
+ stop = true;
+ if (GetTracheaSeedAlgorithm() == 0)
+ has_seed = SearchForTracheaSeed(skip);
+ else
+ has_seed = SearchForTracheaSeed2(GetNumSlices());
+
+ if (has_seed) {
+ TracheaRegionGrowing();
+ volume = ComputeTracheaVolume()/1000; // assume mm3, so divide by 1000 to get cc
+ if (GetWriteStepFlag()) {
+ writeImage<MaskImageType>(trachea, "step-trachea-"+toString(skip)+".mhd");
+ }
+ if (GetTracheaVolumeMustBeCheckedFlag()) {
+ if ((volume > 10) && (volume < 65 )) { // depend on image size ...
+ // Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ]
+ if (GetVerboseStepFlag())
+ {
+ std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
+ }
+ }
+ else
+ if (GetTracheaSeedAlgorithm() == 0) {
+ if (GetVerboseStepFlag()) {
+ 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 (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 << "2 : Number of slices to skip to find trachea too high = " << skip << std::endl;
+ stop = true;
+ }
+ }
+ else {
+ stop = true;
+ }
+ }
+ }
+
+ if (volume != 0.0) {
+ // Set output
+ StopCurrentStep<MaskImageType>(trachea);
+ }
+ else { // Trachea not found
+ this->SetWarning("* WARNING * No seed found for trachea.");
+ StopCurrentStep();
+ }
+}
+//--------------------------------------------------------------------
+