X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=segmentation%2FclitkExtractLungFilter.txx;h=3de027cc3b8210d0dbbac0e5d2a24d1a4aa02cfa;hb=44ebb81a500e42cf19964d209e14a59812a0dd4d;hp=92b2d4f75cf2e8daa0c6633e64c23c7812886ffd;hpb=14f79c6048942f31c8ec28831e8422f99debf54e;p=clitk.git diff --git a/segmentation/clitkExtractLungFilter.txx b/segmentation/clitkExtractLungFilter.txx index 92b2d4f..3de027c 100644 --- a/segmentation/clitkExtractLungFilter.txx +++ b/segmentation/clitkExtractLungFilter.txx @@ -24,7 +24,6 @@ #include "clitkSetBackgroundImageFilter.h" #include "clitkSegmentationUtils.h" #include "clitkAutoCropFilter.h" -#include "clitkExtractSliceFilter.h" // itk #include "itkBinaryThresholdImageFilter.h" @@ -33,13 +32,14 @@ #include "itkOtsuThresholdImageFilter.h" #include "itkBinaryThinningImageFilter3D.h" #include "itkImageIteratorWithIndex.h" +#include "itkBinaryMorphologicalOpeningImageFilter.h" +#include "itkBinaryMorphologicalClosingImageFilter.h" //-------------------------------------------------------------------- template clitk::ExtractLungFilter:: ExtractLungFilter(): clitk::FilterBase(), - clitk::FilterWithAnatomicalFeatureDatabaseManagement(), itk::ImageToImageFilter() { SetNumberOfSteps(10); @@ -84,8 +84,9 @@ ExtractLungFilter(): p3->UseLastKeepOff(); SetLabelizeParameters3(p3); - // Step 5 : find bronchial bifurcations - FindBronchialBifurcationsOn(); + // Step 5 + FinalOpenCloseOff(); + SetFinalOpenCloseRadius(1); } //-------------------------------------------------------------------- @@ -158,7 +159,8 @@ SetArgsInfo(ArgsInfoType mArgsInfo) SetRadiusForTrachea_GGO(mArgsInfo); SetLabelizeParameters3_GGO(mArgsInfo); - SetAFDBFilename_GGO(mArgsInfo); + SetFinalOpenCloseRadius_GGO(mArgsInfo); + SetFinalOpenClose_GGO(mArgsInfo); } //-------------------------------------------------------------------- @@ -170,7 +172,6 @@ clitk::ExtractLungFilter:: GenerateOutputInformation() { Superclass::GenerateOutputInformation(); - //this->GetOutput(0)->SetRequestedRegion(this->GetOutput(0)->GetLargestPossibleRegion()); // Get input pointers patient = dynamic_cast(itk::ProcessObject::GetInput(1)); @@ -335,6 +336,37 @@ GenerateOutputInformation() working_image = cropFilter2->GetOutput(); StopCurrentStep(working_image); + //-------------------------------------------------------------------- + //-------------------------------------------------------------------- + // Final OpenClose + if (GetFinalOpenClose()) { + StartNewStep("Open/Close"); + + // Structuring element + typedef itk::BinaryBallStructuringElement KernelType; + KernelType structuringElement; + structuringElement.SetRadius(GetFinalOpenCloseRadius()); + structuringElement.CreateStructuringElement(); + + // Open + typedef itk::BinaryMorphologicalOpeningImageFilter OpenFilterType; + typename OpenFilterType::Pointer openFilter = OpenFilterType::New(); + openFilter->SetInput(working_image); + openFilter->SetBackgroundValue(GetBackgroundValue()); + openFilter->SetForegroundValue(GetForegroundValue()); + openFilter->SetKernel(structuringElement); + + // Close + typedef itk::BinaryMorphologicalClosingImageFilter CloseFilterType; + typename CloseFilterType::Pointer closeFilter = CloseFilterType::New(); + closeFilter->SetInput(openFilter->GetOutput()); + closeFilter->SetSafeBorder(true); + closeFilter->SetForegroundValue(GetForegroundValue()); + closeFilter->SetKernel(structuringElement); + closeFilter->Update(); + working_image = closeFilter->GetOutput(); + } + //-------------------------------------------------------------------- //-------------------------------------------------------------------- StartNewStep("Separate Left/Right lungs"); @@ -395,125 +427,7 @@ GenerateOutputInformation() // Update output info this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion()); - // Try to extract bifurcation in the trachea (bronchi) - if (m_Seeds.size() != 0) { // if ==0 ->no trachea found - - if (GetFindBronchialBifurcations()) { - StartNewStep("Find bronchial bifurcations"); - // Step 1 : extract skeleton - typedef itk::BinaryThinningImageFilter3D 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 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 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 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; iTransformIndexToPhysicalPoint(listOfBifurcations[i].index, - listOfBifurcations[i].point); - } - // Search for the first slice that separate the bronchus (carena) - typedef clitk::ExtractSliceFilter ExtractSliceFilterType; - typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New(); - extractSliceFilter->SetInput(trachea); - extractSliceFilter->SetDirection(2); - extractSliceFilter->Update(); - typedef typename ExtractSliceFilterType::SliceType SliceType; - std::vector 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(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); - } - } } //-------------------------------------------------------------------- @@ -524,72 +438,12 @@ void clitk::ExtractLungFilter:: GenerateData() { - // If everything goes well, set the output + // Set the output this->GraftOutput(output); // not SetNthOutput } //-------------------------------------------------------------------- -//-------------------------------------------------------------------- -template -void -clitk::ExtractLungFilter:: -TrackFromThisIndex(std::vector & listOfBifurcations, - MaskImagePointer skeleton, - MaskImageIndexType index, - MaskImagePixelType label, - TreeIterator currentNode) -{ - // Create NeighborhoodIterator - typedef itk::NeighborhoodIterator NeighborhoodIteratorType; - typename NeighborhoodIteratorType::SizeType radius; - radius.Fill(1); - NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion()); - - // Track - std::vector listOfTrackedPoint; - bool stop = false; - while (!stop) { - nit.SetLocation(index); - nit.SetCenterPixel(label); - listOfTrackedPoint.clear(); - for(unsigned int i=0; iAdd(listOfTrackedPoint[0], index); // the parent is 'index' - // m_SkeletonTree->Add(listOfTrackedPoint[1], index); // the parent is 'index' - BifurcationType bif(index, label, label+1, label+2); - bif.treeIter = currentNode; - listOfBifurcations.push_back(bif); - TreeIterator firstNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[0]); - TreeIterator secondNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[1]); - TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1, firstNode); - TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2, secondNode); - } - else { - if (listOfTrackedPoint.size() > 2) { - clitkExceptionMacro("error while tracking trachea bifurcation. Too much bifurcation points ... ?"); - } - // Else this it the end of the tracking - } - stop = true; - } - } -} -//-------------------------------------------------------------------- - - //-------------------------------------------------------------------- template bool