X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=segmentation%2FclitkExtractLungFilter.txx;h=a0acdd5ad2ab48aa3669d58f8d9d22171b74d329;hb=d7f456d86ca398a89ccf9de43ab68a2b50b8ca1f;hp=b7b28665baff5006bc2de75cfc158fd8be9821fe;hpb=0793c4dd3d451bb9bdc00e3172fa792d1ddaa8dc;p=clitk.git diff --git a/segmentation/clitkExtractLungFilter.txx b/segmentation/clitkExtractLungFilter.txx index b7b2866..a0acdd5 100644 --- a/segmentation/clitkExtractLungFilter.txx +++ b/segmentation/clitkExtractLungFilter.txx @@ -47,6 +47,8 @@ #include "itkOrientImageFilter.h" #include "itkSpatialOrientationAdapter.h" #include "itkImageDuplicator.h" +#include "itkRelabelComponentImageFilter.h" + #include //-------------------------------------------------------------------- @@ -67,6 +69,7 @@ ExtractLungFilter(): SetForegroundValue(1); SetMinimalComponentSize(100); VerboseRegionGrowingFlagOff(); + RemoveSmallLabelBeforeSeparationFlagOn(); // Step 1 default values SetUpperThreshold(-300); @@ -259,6 +262,9 @@ GenerateOutputInformation() StartNewStep("Search for the trachea"); SearchForTrachea(); PrintMemory(GetVerboseMemoryFlag(), "After SearchForTrachea"); + if (m_Seeds.empty()) { + clitkExceptionMacro("No seeds for trachea... Aborting."); + } //-------------------------------------------------------------------- //-------------------------------------------------------------------- @@ -433,7 +439,7 @@ GenerateOutputInformation() //-------------------------------------------------------------------- //-------------------------------------------------------------------- StartNewStep("Separate Left/Right lungs"); - PrintMemory(GetVerboseMemoryFlag(), "Before Separate"); + PrintMemory(GetVerboseMemoryFlag(), "Before Separate"); // Initial label working_mask = Labelize(working_mask, GetBackgroundValue(), @@ -448,9 +454,33 @@ GenerateOutputInformation() statisticsImageFilter->SetInput(working_mask); statisticsImageFilter->Update(); unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum(); - working_mask = statisticsImageFilter->GetOutput(); - + working_mask = statisticsImageFilter->GetOutput(); PrintMemory(GetVerboseMemoryFlag(), "After count label"); + + // If already 2 labels, but a too big differences, remove the + // smalest one (sometimes appends with the stomach + if (initialNumberOfLabels >1) { + if (GetRemoveSmallLabelBeforeSeparationFlag()) { + DD(GetRemoveSmallLabelBeforeSeparationFlag()); + typedef itk::RelabelComponentImageFilter RelabelFilterType; + typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New(); + relabelFilter->SetInput(working_mask); + relabelFilter->SetMinimumObjectSize(10); + relabelFilter->Update(); + const std::vector & a = relabelFilter->GetSizeOfObjectsInPhysicalUnits(); + std::vector remove_label; + for(unsigned int i=1; i(working_mask, GetBackgroundValue(), remove_label); + statisticsImageFilter->SetInput(working_mask); + statisticsImageFilter->Update(); + initialNumberOfLabels = statisticsImageFilter->GetMaximum(); + } + } // Decompose the first label if (initialNumberOfLabels<2) { @@ -580,9 +610,6 @@ SearchForTracheaSeed(int skip) bool is_orientation_superior(itk::SpatialOrientation::ValidCoordinateOrientationFlags orientation) { itk::SpatialOrientation::CoordinateTerms sup = itk::SpatialOrientation::ITK_COORDINATE_Superior; - std::cout << "orientation: " << std::hex << orientation << "; superior: " << std::hex << sup << std::endl; - std::cout << std::dec; - bool primary = (orientation & 0x0000ff) == sup; bool secondary = ((orientation & 0x00ff00) >> 8) == sup; bool tertiary = ((orientation & 0xff0000) >> 16) == sup; @@ -596,7 +623,7 @@ clitk::ExtractLungFilter:: SearchForTracheaSeed2(int numberOfSlices) { if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user) - if (GetVerboseFlag()) + if (GetVerboseRegionGrowingFlag()) std::cout << "SearchForTracheaSeed2(" << numberOfSlices << ", " << GetMaxElongation() << ")" << std::endl; typedef unsigned char MaskPixelType; @@ -631,6 +658,9 @@ SearchForTracheaSeed2(int numberOfSlices) opening->Update(); typename SlicerFilterType::Pointer slicer = SlicerFilterType::New(); +#if ITK_VERSION_MAJOR >= 4 + slicer->SetDirectionCollapseToIdentity(); +#endif slicer->SetInput(opening->GetOutput()); // label result @@ -702,18 +732,21 @@ SearchForTracheaSeed2(int numberOfSlices) writer->SetFileName(file_name.str().c_str()); writer->Update(); } - - typename LabelImageType::LabelObjectContainerType shapes_map = label_map->GetLabelObjectContainer(); - typename LabelImageType::LabelObjectContainerType::const_iterator s; + 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; - for (s = shapes_map.begin(); s != shapes_map.end(); s++) { - shape = s->second; - if (shape->GetSize() > 150 && shape->GetSize() <= max_size) { + 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++; @@ -746,7 +779,7 @@ SearchForTracheaSeed2(int numberOfSlices) p2[2] = prev_e_centre[2]; double mag = (p2 - p1).GetNorm(); - if (GetVerboseFlag()) { + if (GetVerboseRegionGrowingFlag()) { cout.precision(3); cout << index[2] << ": "; cout << "region(" << max_e_centre[0] << ", " << max_e_centre[1] << ", " << max_e_centre[2] << "); "; @@ -766,6 +799,11 @@ SearchForTracheaSeed2(int numberOfSlices) prev_e_centre= max_e_centre; } + else { + if (GetVerboseRegionGrowingFlag()) { + cout << "No shapes found at slice " << index[2] << std::endl; + } + } } size_t longest = 0; @@ -778,9 +816,11 @@ SearchForTracheaSeed2(int numberOfSlices) } } - if (GetVerboseFlag()) - std::cout << "seed at: " << trachea_centre << std::endl; - m_Seeds.push_back(trachea_centre); + if (longest > 0) { + if (GetVerboseRegionGrowingFlag()) + std::cout << "seed at: " << trachea_centre << std::endl; + m_Seeds.push_back(trachea_centre); + } } return (m_Seeds.size() != 0); @@ -802,7 +842,7 @@ TracheaRegionGrowing() f->SetUpper(GetUpperThresholdForTrachea()); f->SetMinimumLowerThreshold(-2000); // f->SetMaximumUpperThreshold(0); // MAYBE TO CHANGE ??? - f->SetMaximumUpperThreshold(-700); // MAYBE TO CHANGE ??? + f->SetMaximumUpperThreshold(-300); // MAYBE TO CHANGE ??? f->SetAdaptLowerBorder(false); f->SetAdaptUpperBorder(true); f->SetMinimumSize(5000); @@ -892,16 +932,17 @@ SearchForTrachea() std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl; } } - else { - 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(); + 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