From: Romulo Pinho Date: Fri, 27 May 2011 10:50:23 +0000 (+0200) Subject: Merge branch 'master' of /home/dsarrut/clitk3.server X-Git-Tag: v1.3.0~338 X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=a515cd4eb114c9310570341ca29f35484356a1e8;hp=12be037aac07fd0f9cb046e78e983ed95f89d877;p=clitk.git Merge branch 'master' of /home/dsarrut/clitk3.server --- diff --git a/common/clitkDicomRT_ROI.cxx b/common/clitkDicomRT_ROI.cxx index 25a2f43..88bab72 100644 --- a/common/clitkDicomRT_ROI.cxx +++ b/common/clitkDicomRT_ROI.cxx @@ -36,7 +36,6 @@ clitk::DicomRT_ROI::DicomRT_ROI() mMeshIsUpToDate = false; mBackgroundValue = 0; mForegroundValue = 1; - mZDelta = 0; } //-------------------------------------------------------------------- @@ -163,6 +162,8 @@ void clitk::DicomRT_ROI::Read(std::map & rois, gdcm::Item cons gdcm::Tag tcsq(0x3006,0x0040); if( !nestedds.FindDataElement( tcsq ) ) { + std::cerr << "Warning. Could not read contour for structure <" << mName << ">, number" << mNumber << " ? I ignore it" << std::endl; + return; } const gdcm::DataElement& csq = nestedds.GetDataElement( tcsq ); gdcm::SmartPointer sqi2 = csq.GetValueAsSQ(); @@ -171,28 +172,15 @@ void clitk::DicomRT_ROI::Read(std::map & rois, gdcm::Item cons } unsigned int nitems = sqi2->GetNumberOfItems(); - bool contour_processed=false; - bool delta_computed=false; - double last_z=0; for(unsigned int i = 0; i < nitems; ++i) { - const gdcm::Item & j = sqi2->GetItem(i+1); // Item start at #1 - DicomRT_Contour::Pointer c = DicomRT_Contour::New(); - bool b = c->Read(j); - if (b) { - mListOfContours.push_back(c); - if (contour_processed) { - double delta=c->GetZ() - last_z; - if (delta_computed) - assert(mZDelta == delta); - else - mZDelta = delta; - } else - contour_processed=true; - last_z=c->GetZ(); + const gdcm::Item & j = sqi2->GetItem(i+1); // Item start at #1 + DicomRT_Contour::Pointer c = DicomRT_Contour::New(); + bool b = c->Read(j); + if (b) { + mListOfContours.push_back(c); + } } - } - } #else void clitk::DicomRT_ROI::Read(std::map & rois, gdcm::SQItem * item) @@ -213,25 +201,20 @@ void clitk::DicomRT_ROI::Read(std::map & rois, gdcm::SQItem * // Read contours [Contour Sequence] gdcm::SeqEntry * contours=item->GetSeqEntry(0x3006,0x0040); - bool contour_processed=false; - bool delta_computed=false; - double last_z=0; - for(gdcm::SQItem* j=contours->GetFirstSQItem(); j!=0; j=contours->GetNextSQItem()) { - DicomRT_Contour::Pointer c = DicomRT_Contour::New(); - bool b = c->Read(j); - if (b) { - mListOfContours.push_back(c); - if (contour_processed) { - double delta=c->GetZ() - last_z; - if (delta_computed) - assert(mZDelta == delta); - else - mZDelta = delta; - } else - contour_processed=true; - last_z=c->GetZ(); + if (contours) { + int i=0; + for(gdcm::SQItem* j=contours->GetFirstSQItem(); j!=0; j=contours->GetNextSQItem()) { + DicomRT_Contour::Pointer c = DicomRT_Contour::New(); + bool b = c->Read(j); + if (b) { + mListOfContours.push_back(c); + } + ++i; } } + else { + std::cerr << "Warning. Could not read contour for structure <" << mName << ">, number" << mNumber << " ? I ignore it" << std::endl; + } } #endif //-------------------------------------------------------------------- @@ -277,9 +260,9 @@ void clitk::DicomRT_ROI::ComputeMesh() //-------------------------------------------------------------------- void clitk::DicomRT_ROI::SetFromBinaryImage(vvImage * image, int n, - std::string name, - std::vector color, - std::string filename) + std::string name, + std::vector color, + std::string filename) { // ROI number [Referenced ROI Number] diff --git a/common/clitkDicomRT_ROI.h b/common/clitkDicomRT_ROI.h index a22322e..0fb3076 100644 --- a/common/clitkDicomRT_ROI.h +++ b/common/clitkDicomRT_ROI.h @@ -64,7 +64,7 @@ public: void SetImage(vvImage * im); DicomRT_Contour* GetContour(int n); - double GetContourSpacing() const {return mZDelta;} + // double GetContourSpacing() const {return mZDelta;} protected: void ComputeMesh(); @@ -79,7 +79,7 @@ protected: double mBackgroundValue; double mForegroundValue; ///Spacing between two contours - double mZDelta; + // double mZDelta; private: DicomRT_ROI(); diff --git a/common/clitkDicomRT_ROI_ConvertToImageFilter.cxx b/common/clitkDicomRT_ROI_ConvertToImageFilter.cxx index 07c7a95..e3b3472 100644 --- a/common/clitkDicomRT_ROI_ConvertToImageFilter.cxx +++ b/common/clitkDicomRT_ROI_ConvertToImageFilter.cxx @@ -19,12 +19,18 @@ #include #include + +// clitk #include "clitkDicomRT_ROI_ConvertToImageFilter.h" +#include "clitkImageCommon.h" + +// vtk #include #include #include #include -#include "clitkImageCommon.h" +#include + //-------------------------------------------------------------------- clitk::DicomRT_ROI_ConvertToImageFilter::DicomRT_ROI_ConvertToImageFilter() @@ -121,11 +127,9 @@ void clitk::DicomRT_ROI_ConvertToImageFilter::Update() std::cerr << "Error. Please provide image info (spacing/origin) with SetImageFilename" << std::endl; exit(0); } - // DD("Update"); // Get Mesh vtkPolyData * mesh = mROI->GetMesh(); - DD(mesh->GetNumberOfCells()); // Get bounds double *bounds=mesh->GetBounds(); @@ -176,7 +180,8 @@ void clitk::DicomRT_ROI_ConvertToImageFilter::Update() // Extrude vtkSmartPointer extrude=vtkSmartPointer::New(); extrude->SetInput(mesh); - extrude->SetVector(0, 0, mROI->GetContourSpacing()); + ///We extrude in the -slice_spacing direction to respect the FOCAL convention (NEEDED !) + extrude->SetVector(0, 0, -mSpacing[2]); // Binarization vtkSmartPointer sts=vtkSmartPointer::New(); @@ -192,6 +197,14 @@ void clitk::DicomRT_ROI_ConvertToImageFilter::Update() stencil->SetInput(mBinaryImage); stencil->ReverseStencilOn(); stencil->Update(); + + /* + vtkSmartPointer w = vtkSmartPointer::New(); + w->SetInput(stencil->GetOutput()); + w->SetFileName("binary2.mhd"); + w->Write(); + */ + mBinaryImage->ShallowCopy(stencil->GetOutput()); if (mWriteOutput) { diff --git a/common/clitkDicomRT_StructureSet.cxx b/common/clitkDicomRT_StructureSet.cxx index f0b690f..8543fb7 100644 --- a/common/clitkDicomRT_StructureSet.cxx +++ b/common/clitkDicomRT_StructureSet.cxx @@ -274,7 +274,12 @@ void clitk::DicomRT_StructureSet::Read(const std::string & filename) mStudyTime = reader.GetValEntry(0x008,0x0020)->GetValue(); mStudyDate = reader.GetValEntry(0x008,0x0030)->GetValue(); mLabel = reader.GetValEntry(0x3006,0x002)->GetValue(); - mName = reader.GetValEntry(0x3006,0x004)->GetValue(); + if (!reader.GetValEntry(0x3006,0x004)) { + mName = "Anonymous"; + } + else { + mName = reader.GetValEntry(0x3006,0x004)->GetValue(); + } mTime = reader.GetValEntry(0x3006,0x009)->GetValue(); //---------------------------------- diff --git a/itk/clitkAutoCropFilter.txx b/itk/clitkAutoCropFilter.txx index 0005eb8..71f9ab0 100644 --- a/itk/clitkAutoCropFilter.txx +++ b/itk/clitkAutoCropFilter.txx @@ -38,7 +38,7 @@ namespace clitk { AutoCropFilter():itk::ImageToImageFilter() { this->SetNumberOfRequiredInputs(1); m_BackgroundValue = 0; - UseBorderOn(); + UseBorderOff(); } //-------------------------------------------------------------------- @@ -97,6 +97,8 @@ namespace clitk { autoCropFilter->SetInput(imageToLabelFilter->GetOutput()); // autoCropFilter->ReleaseDataFlagOff(); if (GetUseBorder()) { + DD("UseBorder seems buggy ?"); + exit(0); typename ImageType::SizeType s; for(uint i=0; iSetCropBorder(s); diff --git a/itk/clitkSegmentationUtils.txx b/itk/clitkSegmentationUtils.txx index 36594aa..f7b897e 100644 --- a/itk/clitkSegmentationUtils.txx +++ b/itk/clitkSegmentationUtils.txx @@ -511,7 +511,6 @@ namespace clitk { typename ImageType::PointType dummy; centroids.push_back(dummy); // label 0 -> no centroid, use dummy point for BG //DS FIXME (not useful ! to change ..) - DD(labelMap->GetNumberOfLabelObjects()); for(uint i=0; iGetNumberOfLabelObjects(); i++) { int label = labelMap->GetLabels()[i]; centroids.push_back(labelMap->GetLabelObject(label)->GetCentroid()); @@ -1061,21 +1060,17 @@ namespace clitk { typedef typename itk::Image SliceType; std::vector slices; clitk::ExtractSlices(input, d, slices); - DD(slices.size()); // Labelize and keep the main one std::vector o; for(uint i=0; i(slices[i], BG, false, 1)); o[i] = clitk::KeepLabels(o[i], BG, FG, 1, 1, true); } // Join slices - DD("join"); typename ImageType::Pointer output; output = clitk::JoinSlices(o, input, d); - DD("return"); return output; } //-------------------------------------------------------------------- diff --git a/segmentation/CMakeLists.txt b/segmentation/CMakeLists.txt index ad4e11e..3b67730 100644 --- a/segmentation/CMakeLists.txt +++ b/segmentation/CMakeLists.txt @@ -43,9 +43,13 @@ IF(CLITK_BUILD_SEGMENTATION) ADD_EXECUTABLE(clitkExtractMediastinum clitkExtractMediastinum.cxx ${clitkExtractMediastinum_GGO_C}) TARGET_LINK_LIBRARIES(clitkExtractMediastinum clitkCommon clitkSegmentationGgoLib ${ITK_LIBRARIES}) - # WRAP_GGO(clitkExtractLymphStations_GGO_C clitkExtractLymphStations.ggo) - # ADD_EXECUTABLE(clitkExtractLymphStations clitkExtractLymphStations.cxx clitkFilterWithAnatomicalFeatureDatabaseManagement.cxx ${clitkExtractLymphStations_GGO_C}) - # TARGET_LINK_LIBRARIES(clitkExtractLymphStations clitkSegmentationGgoLib clitkCommon ITKIO ITKStatistics vtkHybrid) + #WRAP_GGO(clitkExtractLymphStations_GGO_C clitkExtractLymphStations.ggo) + #ADD_EXECUTABLE(clitkExtractLymphStations clitkExtractLymphStations.cxx clitkFilterWithAnatomicalFeatureDatabaseManagement.cxx ${clitkExtractLymphStations_GGO_C}) + #TARGET_LINK_LIBRARIES(clitkExtractLymphStations clitkSegmentationGgoLib clitkCommon vtkHybrid) + + WRAP_GGO(clitkExtractMediastinalVessels_GGO_C clitkExtractMediastinalVessels.ggo) + ADD_EXECUTABLE(clitkExtractMediastinalVessels clitkExtractMediastinalVessels.cxx clitkFilterWithAnatomicalFeatureDatabaseManagement.cxx ${clitkExtractMediastinalVessels_GGO_C}) + TARGET_LINK_LIBRARIES(clitkExtractMediastinalVessels clitkSegmentationGgoLib clitkCommon vtkHybrid) WRAP_GGO(clitkMorphoMath_GGO_C clitkMorphoMath.ggo) ADD_EXECUTABLE(clitkMorphoMath clitkMorphoMath.cxx ${clitkMorphoMath_GGO_C}) @@ -83,9 +87,9 @@ IF(CLITK_BUILD_SEGMENTATION) ADD_EXECUTABLE(clitkFillImageRegion clitkFillImageRegion.cxx clitkFillImageRegionGenericFilter.cxx ${clitkFillImageRegion_GGO_C}) TARGET_LINK_LIBRARIES(clitkFillImageRegion clitkCommon ${ITK_LIBRARIES}) - # WRAP_GGO(clitkTestFilter_GGO_C clitkTestFilter.ggo) - # ADD_EXECUTABLE(clitkTestFilter clitkTestFilter.cxx ${clitkTestFilter_GGO_C}) - # TARGET_LINK_LIBRARIES(clitkTestFilter clitkSegmentationGgoLib clitkCommon vtkHybrid ITKIO) + WRAP_GGO(clitkTestFilter_GGO_C clitkTestFilter.ggo) + ADD_EXECUTABLE(clitkTestFilter clitkTestFilter.cxx ${clitkTestFilter_GGO_C}) + TARGET_LINK_LIBRARIES(clitkTestFilter clitkSegmentationGgoLib clitkCommon vtkHybrid) ENDIF(CLITK_BUILD_SEGMENTATION) diff --git a/segmentation/clitkExtractLymphStation_2RL.txx b/segmentation/clitkExtractLymphStation_2RL.txx index 542540f..7e8c834 100644 --- a/segmentation/clitkExtractLymphStation_2RL.txx +++ b/segmentation/clitkExtractLymphStation_2RL.txx @@ -94,7 +94,7 @@ clitk::ExtractLymphStationsFilter:: ExtractStation_2RL_SetDefaultValues() { SetFuzzyThreshold("2RL", "CommonCarotidArtery", 0.7); - SetFuzzyThreshold("2RL", "BrachioCephalicTrunk", 0.7); + SetFuzzyThreshold("2RL", "BrachioCephalicArtery", 0.7); SetFuzzyThreshold("2RL", "BrachioCephalicVein", 0.3); SetFuzzyThreshold("2RL", "Aorta", 0.7); SetFuzzyThreshold("2RL", "SubclavianArteryRight", 0.5); @@ -208,17 +208,17 @@ ExtractStation_2RL_Ant_Limits() // ----------------------------------------------------- // Remove Ant to H line from the Ant most part of the // CommonCarotidArtery until we reach the first slice of - // BrachioCephalicTrunk + // BrachioCephalicArtery StartNewStep("[Station 2RL] Ant limits with CommonCarotidArtery, H line"); - // First, find the first slice of BrachioCephalicTrunk - MaskImagePointer BrachioCephalicTrunk = GetAFDB()->template GetImage("BrachioCephalicTrunk"); - MaskImagePointType p = BrachioCephalicTrunk->GetOrigin(); // initialise to avoid warning - clitk::FindExtremaPointInAGivenDirection(BrachioCephalicTrunk, GetBackgroundValue(), 2, false, p); - double TopOfBrachioCephalicTrunkZ=p[2] + BrachioCephalicTrunk->GetSpacing()[2]; // Add one slice + // First, find the first slice of BrachioCephalicArtery + MaskImagePointer BrachioCephalicArtery = GetAFDB()->template GetImage("BrachioCephalicArtery"); + MaskImagePointType p = BrachioCephalicArtery->GetOrigin(); // initialise to avoid warning + clitk::FindExtremaPointInAGivenDirection(BrachioCephalicArtery, GetBackgroundValue(), 2, false, p); + double TopOfBrachioCephalicArteryZ=p[2] + BrachioCephalicArtery->GetSpacing()[2]; // Add one slice // Remove CommonCarotidArtery below this point - CommonCarotidArtery = clitk::CropImageRemoveLowerThan(CommonCarotidArtery, 2, TopOfBrachioCephalicTrunkZ, true, GetBackgroundValue()); + CommonCarotidArtery = clitk::CropImageRemoveLowerThan(CommonCarotidArtery, 2, TopOfBrachioCephalicArteryZ, true, GetBackgroundValue()); // Find most Ant points std::vector ccaAntPositionA; @@ -242,26 +242,26 @@ ExtractStation_2RL_Ant_Limits() m_ListOfStations["2L"] = m_Working_Support; // ----------------------------------------------------- - // Ant limit with the BrachioCephalicTrunk - StartNewStep("[Station 2RL] Ant limits with BrachioCephalicTrunk line"); + // Ant limit with the BrachioCephalicArtery + StartNewStep("[Station 2RL] Ant limits with BrachioCephalicArtery line"); - // Remove Ant to BrachioCephalicTrunk + // Remove Ant to BrachioCephalicArtery m_Working_Support = - clitk::SliceBySliceRelativePosition(m_Working_Support, BrachioCephalicTrunk, 2, - GetFuzzyThreshold("2RL", "BrachioCephalicTrunk"), "NotAntTo", false, 2, true, false); + clitk::SliceBySliceRelativePosition(m_Working_Support, BrachioCephalicArtery, 2, + GetFuzzyThreshold("2RL", "BrachioCephalicArtery"), "NotAntTo", false, 2, true, false); // End StopCurrentStep(m_Working_Support); m_ListOfStations["2R"] = m_Working_Support; m_ListOfStations["2L"] = m_Working_Support; // ----------------------------------------------------- - // Ant limit with the BrachioCephalicTrunk H line - StartNewStep("[Station 2RL] Ant limits with BrachioCephalicTrunk, Horizontal line"); + // Ant limit with the BrachioCephalicArtery H line + StartNewStep("[Station 2RL] Ant limits with BrachioCephalicArtery, Horizontal line"); // Find most Ant points std::vector bctAntPositionA; std::vector bctAntPositionB; - clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition(BrachioCephalicTrunk, + clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition(BrachioCephalicArtery, GetBackgroundValue(), 2, 1, true, // Ant 0, // Horizontal line @@ -317,7 +317,7 @@ ExtractStation_2RL_Ant_Limits2() /* Here, we consider the vessels form a kind of anterior barrier. We link all vessels centroids and remove what is post to it. - select the list of structure - vessel1 = BrachioCephalicTrunk + vessel1 = BrachioCephalicArtery vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right) vessel3 = CommonCarotidArtery vessel4 = SubclavianArtery @@ -332,7 +332,7 @@ ExtractStation_2RL_Ant_Limits2() */ // Read structures - MaskImagePointer BrachioCephalicTrunk = GetAFDB()->template GetImage("BrachioCephalicTrunk"); + MaskImagePointer BrachioCephalicArtery = GetAFDB()->template GetImage("BrachioCephalicArtery"); MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage("BrachioCephalicVein"); MaskImagePointer CommonCarotidArtery = GetAFDB()->template GetImage("CommonCarotidArtery"); MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage("SubclavianArtery"); @@ -341,8 +341,8 @@ ExtractStation_2RL_Ant_Limits2() MaskImagePointer Trachea = GetAFDB()->template GetImage("Trachea"); // Resize all structures like support - BrachioCephalicTrunk = - clitk::ResizeImageLike(BrachioCephalicTrunk, m_Working_Support, GetBackgroundValue()); + BrachioCephalicArtery = + clitk::ResizeImageLike(BrachioCephalicArtery, m_Working_Support, GetBackgroundValue()); CommonCarotidArtery = clitk::ResizeImageLike(CommonCarotidArtery, m_Working_Support, GetBackgroundValue()); SubclavianArtery = @@ -357,8 +357,8 @@ ExtractStation_2RL_Ant_Limits2() clitk::ResizeImageLike(Trachea, m_Working_Support, GetBackgroundValue()); // Extract slices - std::vector slices_BrachioCephalicTrunk; - clitk::ExtractSlices(BrachioCephalicTrunk, 2, slices_BrachioCephalicTrunk); + std::vector slices_BrachioCephalicArtery; + clitk::ExtractSlices(BrachioCephalicArtery, 2, slices_BrachioCephalicArtery); std::vector slices_BrachioCephalicVein; clitk::ExtractSlices(BrachioCephalicVein, 2, slices_BrachioCephalicVein); std::vector slices_CommonCarotidArtery; @@ -371,11 +371,11 @@ ExtractStation_2RL_Ant_Limits2() clitk::ExtractSlices(Aorta, 2, slices_Aorta); std::vector slices_Trachea; clitk::ExtractSlices(Trachea, 2, slices_Trachea); - uint n= slices_BrachioCephalicTrunk.size(); + uint n= slices_BrachioCephalicArtery.size(); // Get the boundaries of one slice std::vector bounds; - ComputeImageBoundariesCoordinates(slices_BrachioCephalicTrunk[0], bounds); + ComputeImageBoundariesCoordinates(slices_BrachioCephalicArtery[0], bounds); // For all slices, for all structures, find the centroid and build the contour // List of 3D points (for debug) @@ -388,7 +388,7 @@ ExtractStation_2RL_Ant_Limits2() GetBackgroundValue(), true, 1); slices_SubclavianArtery[i] = Labelize(slices_SubclavianArtery[i], GetBackgroundValue(), true, 1); - slices_BrachioCephalicTrunk[i] = Labelize(slices_BrachioCephalicTrunk[i], + slices_BrachioCephalicArtery[i] = Labelize(slices_BrachioCephalicArtery[i], GetBackgroundValue(), true, 1); slices_BrachioCephalicVein[i] = Labelize(slices_BrachioCephalicVein[i], GetBackgroundValue(), true, 1); @@ -407,7 +407,7 @@ ExtractStation_2RL_Ant_Limits2() std::vector centroids6; ComputeCentroids(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1); ComputeCentroids(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2); - ComputeCentroids(slices_BrachioCephalicTrunk[i], GetBackgroundValue(), centroids3); + ComputeCentroids(slices_BrachioCephalicArtery[i], GetBackgroundValue(), centroids3); ComputeCentroids(slices_Thyroid[i], GetBackgroundValue(), centroids4); ComputeCentroids(slices_Aorta[i], GetBackgroundValue(), centroids5); ComputeCentroids(slices_BrachioCephalicVein[i], GetBackgroundValue(), centroids6); @@ -420,7 +420,7 @@ ExtractStation_2RL_Ant_Limits2() } // BrachioCephalicVein -> when SubclavianArtery has 2 CCL - // (BrachioCephalicTrunk is divided) -> forget BrachioCephalicVein + // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein if ((centroids3.size() ==1) && (centroids2.size() > 2)) { centroids6.clear(); } diff --git a/segmentation/clitkExtractMediastinalVessels.ggo b/segmentation/clitkExtractMediastinalVessels.ggo index 718c24e..c87d960 100644 --- a/segmentation/clitkExtractMediastinalVessels.ggo +++ b/segmentation/clitkExtractMediastinalVessels.ggo @@ -16,6 +16,16 @@ section "I/O" option "afdb" a "Input Anatomical Feature DB" string yes option "input" i "Input CT filename" string yes -option "output" o "Output folder" string yes +option "output" o "Output filename" string yes +option "seed" s "Seed point name in AFDB" string yes +option "name" n "Name of the vessels to sought" string yes -option "threshold" t "Initial threshold" double default="140" no +option "thresholdHigh" t "Initial high threshold" double default="140" no +option "thresholdLow" l "Initial low threshold" double default="55" no +option "erode" e "Erosion radius in pixel" int default="2" no +option "dilate" d "Dilatation radius in pixel" int default="9" no +option "maxPost" - "Max distance post to Carina" double default="10" no +option "maxAnt" - "Max distance ant to Carina" double default="40" no +option "maxLeft" - "Max distance left to Carina" double default="35" no +option "maxRight" - "Max distance right to Carina" double default="35" no +option "bif" - "Max number of bifurcation during traking" int default="0" no diff --git a/segmentation/clitkExtractMediastinalVesselsFilter.h b/segmentation/clitkExtractMediastinalVesselsFilter.h index 665f79c..59ff6f7 100644 --- a/segmentation/clitkExtractMediastinalVesselsFilter.h +++ b/segmentation/clitkExtractMediastinalVesselsFilter.h @@ -88,11 +88,42 @@ namespace clitk { itkGetConstMacro(TemporaryForegroundValue, MaskImagePixelType); itkSetMacro(TemporaryForegroundValue, MaskImagePixelType); - itkGetConstMacro(OutputFolder, std::string); - itkSetMacro(OutputFolder, std::string); + itkGetConstMacro(ThresholdHigh, ImagePixelType); + itkSetMacro(ThresholdHigh, ImagePixelType); - itkGetConstMacro(Threshold, ImagePixelType); - itkSetMacro(Threshold, ImagePixelType); + itkGetConstMacro(ThresholdLow, ImagePixelType); + itkSetMacro(ThresholdLow, ImagePixelType); + + itkGetConstMacro(ErosionRadius, int); + itkSetMacro(ErosionRadius, int); + + itkGetConstMacro(DilatationRadius, int); + itkSetMacro(DilatationRadius, int); + + itkGetConstMacro(MaxDistancePostToCarina, double); + itkSetMacro(MaxDistancePostToCarina, double); + itkGetConstMacro(MaxDistanceAntToCarina, double); + itkSetMacro(MaxDistanceAntToCarina, double); + itkGetConstMacro(MaxDistanceLeftToCarina, double); + itkSetMacro(MaxDistanceLeftToCarina, double); + itkGetConstMacro(MaxDistanceRightToCarina, double); + itkSetMacro(MaxDistanceRightToCarina, double); + + itkSetMacro(DebugFlag, bool); + itkGetConstMacro(DebugFlag, bool); + itkBooleanMacro(DebugFlag); + + itkSetMacro(SoughtVesselSeedName, std::string); + itkGetConstMacro(SoughtVesselSeedName, std::string); + + itkSetMacro(SoughtVesselName, std::string); + itkGetConstMacro(SoughtVesselName, std::string); + + itkSetMacro(OutputFilename, std::string); + itkGetConstMacro(OutputFilename, std::string); + + itkSetMacro(MaxNumberOfFoundBifurcation, int); + itkGetConstMacro(MaxNumberOfFoundBifurcation, int); protected: ExtractMediastinalVesselsFilter(); @@ -102,7 +133,7 @@ namespace clitk { virtual void GenerateInputRequestedRegion(); virtual void GenerateData(); - std::string m_OutputFolder; + bool m_DebugFlag; ImagePointer m_Input; MaskImagePointer m_Working_Support; MaskImagePointer m_Mediastinum; @@ -110,12 +141,26 @@ namespace clitk { MaskImagePixelType m_BackgroundValue; MaskImagePixelType m_ForegroundValue; MaskImagePixelType m_TemporaryForegroundValue; - ImagePixelType m_Threshold; + ImagePixelType m_ThresholdHigh; + ImagePixelType m_ThresholdLow; + int m_ErosionRadius; + int m_DilatationRadius; + double m_MaxDistancePostToCarina; + double m_MaxDistanceAntToCarina; + double m_MaxDistanceLeftToCarina; + double m_MaxDistanceRightToCarina; + int m_MaxNumberOfFoundBifurcation; std::vector m_slice_recon; + std::vector m_slice_recon2; + + // Resulting structures + MaskImageType::Pointer m_SoughtVessel; + std::string m_SoughtVesselSeedName; + std::string m_SoughtVesselName; + std::string m_OutputFilename; - void CropSupInf(); - //void SearchBrachioCephalicArtery(int & BCA_first_slice, LabelType & BCA_first_label); + void CropInputImage(); void TrackBifurcationFromPoint(MaskImagePointer & recon, std::vector & slices_recon, MaskImagePointType BCA_p, diff --git a/segmentation/clitkExtractMediastinalVesselsFilter.txx b/segmentation/clitkExtractMediastinalVesselsFilter.txx index 02576f4..2198529 100644 --- a/segmentation/clitkExtractMediastinalVesselsFilter.txx +++ b/segmentation/clitkExtractMediastinalVesselsFilter.txx @@ -22,13 +22,11 @@ // clitk #include "clitkCommon.h" #include "clitkExtractMediastinalVesselsFilter.h" -#include "clitkAutoCropFilter.h" #include "clitkSegmentationUtils.h" -#include "clitkMorphoMathFilter.h" +#include "clitkReconstructWithConditionalGrayscaleDilateImageFilter.h" // itk #include -#include #include //-------------------------------------------------------------------- @@ -42,8 +40,15 @@ ExtractMediastinalVesselsFilter(): this->SetNumberOfRequiredInputs(1); SetBackgroundValue(0); SetForegroundValue(1); - SetThreshold(140); - SetTemporaryForegroundValue(1); + SetThresholdHigh(140); + SetThresholdLow(55); + SetErosionRadius(2); + SetDilatationRadius(9); + SetMaxDistancePostToCarina(10); + SetMaxDistanceAntToCarina(40); + SetMaxDistanceLeftToCarina(35); + SetMaxDistanceRightToCarina(35); + SetSoughtVesselSeedName("NoSeedNameGiven"); } //-------------------------------------------------------------------- @@ -58,145 +63,149 @@ GenerateOutputInformation() { m_Input = dynamic_cast(itk::ProcessObject::GetInput(0)); // ------------------------------------------------------------------ - // Sup-Inf crop -> Carina - CropSupInf(); + // Crop the initial image superiorly and inferiorly. + // TODO : add options for x cm above/below + CropInputImage(); // ------------------------------------------------------------------ - // Binarize - StartNewStep("Binarize with treshold = "+toString(GetThreshold())); + // Binarize the image. Need two thresholds, one high to select + // structures (CCL) that are almost not connected (after erosion), + // and one low thresholds to select the real contours. Will be + // reconstructed later. + StartNewStep("Binarize with high threshold = "+toString(GetThresholdHigh())); typedef itk::BinaryThresholdImageFilter BinarizeFilterType; typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New(); binarizeFilter->SetInput(m_Input); - binarizeFilter->SetLowerThreshold(GetThreshold()); - binarizeFilter->SetInsideValue(GetTemporaryForegroundValue()); + binarizeFilter->SetLowerThreshold(GetThresholdHigh()); + binarizeFilter->SetInsideValue(GetForegroundValue()); binarizeFilter->SetOutsideValue(GetBackgroundValue()); binarizeFilter->Update(); m_Mask = binarizeFilter->GetOutput(); - clitk::writeImage(m_Mask, "m.mhd"); StopCurrentStep(m_Mask); - // Keep main CCL ? - m_Mask = clitk::Labelize(m_Mask, GetBackgroundValue(), false, 10); - m_Mask = KeepLabels(m_Mask, GetBackgroundValue(), GetTemporaryForegroundValue(), 1, 1, true); - clitk::writeImage(m_Mask, "m2.mhd"); - + // ------------------------------------------------------------------ + StartNewStep("Binarize with low threshold = "+toString(GetThresholdLow())); + binarizeFilter = BinarizeFilterType::New(); + binarizeFilter->SetInput(m_Input); + binarizeFilter->SetLowerThreshold(GetThresholdLow()); + binarizeFilter->SetInsideValue(GetForegroundValue()); + binarizeFilter->SetOutsideValue(GetBackgroundValue()); + binarizeFilter->Update(); + MaskImagePointer m_Mask2 = binarizeFilter->GetOutput(); + StopCurrentStep(m_Mask2); + // ------------------------------------------------------------------ // Extract slices - StartNewStep("Detect vessels (slice by slice reconstruction)"); + StartNewStep("Detect objects : erosion, then slice by slice reconstruction"); std::vector slices_mask; clitk::ExtractSlices(m_Mask, 2, slices_mask); - - DD(slices_mask.size()); - + std::vector slices_mask2; + clitk::ExtractSlices(m_Mask2, 2, slices_mask2); + int radius = GetErosionRadius(); + + // List of working slices (debug only) std::vector debug_eroded; - std::vector debug_labeled; - std::vector debug_slabeled; + // std::vector debug_labeled; - int radius = 3; - DD(radius); // TO PUT IN OPTION - // ------------------------------------------------------------------ - // Loop Slice by Slice -> erode find CCL and reconstruct - clitk::MorphoMathFilter::Pointer f= clitk::MorphoMathFilter::New(); + // Loop Slice by Slice in order to break connectivity between + // CCL. Erode and reconstruct all labels at the same time without + // merging them. for(uint i=0; iSetInput(slices_mask[i]); - f->SetBackgroundValue(GetBackgroundValue()); - f->SetForegroundValue(GetTemporaryForegroundValue()); - f->SetRadius(radius); - f->SetOperationType(0); // Erode - f->VerboseFlagOff(); - f->Update(); - MaskSlicePointer eroded = f->GetOutput(); - // writeImage(eroded, "er-"+toString(i)+".mhd"); + // Erosion kernel + typedef itk::BinaryBallStructuringElement KernelType; + KernelType structuringElement; + structuringElement.SetRadius(radius); + structuringElement.CreateStructuringElement(); + + // Erosion -> we break the connectivity between structure + typedef itk::BinaryErodeImageFilter ErodeFilterType; + typename ErodeFilterType::Pointer eroder = ErodeFilterType::New(); + eroder->SetInput(slices_mask[i]); + eroder->SetBackgroundValue(GetBackgroundValue()); + eroder->SetForegroundValue(GetForegroundValue()); + eroder->SetBoundaryToForeground(true); // ?? + eroder->SetKernel(structuringElement); + eroder->Update(); + MaskSlicePointer eroded = eroder->GetOutput(); debug_eroded.push_back(eroded); - // CCL - int nb; + // Labelize (CCL) MaskSlicePointer labeled = - clitk::LabelizeAndCountNumberOfObjects(eroded, GetBackgroundValue(), false, 1, nb); - - // Relabel, large CCL with large label number - for(int n=nb; n>0; n--) { - // DD(n); - int li = n; - int lo = 2*(nb+1)-li; - labeled = clitk::SetBackground(labeled, labeled, li, lo, true); - } - debug_labeled.push_back(labeled); - - // Create kernel for GrayscaleDilateImageFilter - typedef itk::BinaryBallStructuringElement KernelType; - KernelType k; - k.SetRadius(radius+1); - k.CreateStructuringElement(); - - // Keep the MAX -> we prefer the opposite su change the label - typedef itk::GrayscaleDilateImageFilter FilterType; - FilterType::Pointer m = FilterType::New(); - m->SetKernel(k); - m->SetInput(labeled); - // DD(m->GetAlgorithm()); - // m->SetAlgorithm(3); - m->Update(); - MaskSlicePointer s = m->GetOutput(); - - - // Remove Initial BG - s = clitk::SetBackground(s, slices_mask[i], - GetBackgroundValue(), GetBackgroundValue(), true); - + clitk::Labelize(eroded, GetBackgroundValue(), true, 1); // Fully connected ! + // debug_labeled.push_back(labeled); + + // Make Reconstruction filter : dilation all labels at the same + // time, prevent to merge them. + typedef clitk::ReconstructWithConditionalGrayscaleDilateImageFilter ReconstructFilterType; + typename ReconstructFilterType::Pointer reconstructor = ReconstructFilterType::New(); + reconstructor->SetInput(labeled); + reconstructor->SetIterationNumber(radius+GetDilatationRadius()); + reconstructor->Update(); + MaskSlicePointer s = reconstructor->GetOutput(); + + // Remove Initial BG of the second tresholded image + s = clitk::SetBackground(s, slices_mask2[i], + GetBackgroundValue(), GetBackgroundValue(), true); m_slice_recon.push_back(s); + } // end loop - DD("end loop"); + // Build 3D images from the slice by slice processing MaskImageType::Pointer eroded = clitk::JoinSlices(debug_eroded, m_Mask, 2); - clitk::writeImage(eroded, "eroded.mhd"); - - DD("l"); - MaskImageType::Pointer l = clitk::JoinSlices(debug_labeled, m_Mask, 2); - clitk::writeImage(l, "labeled.mhd"); - - DD("r"); + writeImage(eroded, "erode.mhd"); + //MaskImageType::Pointer l = clitk::JoinSlices(debug_labeled, m_Mask, 2); MaskImageType::Pointer r = clitk::JoinSlices(m_slice_recon, m_Mask, 2); - clitk::writeImage(r, "recon.mhd"); + writeImage(r, "recon1.mhd"); // ------------------------------------------------------------------ - // Loop Slice by Slice -> BCA not found yet - /* MaskImagePointType BCA_p; - GetAFDB()->GetPoint3D("BrachioCephalicArteryFirstInferiorPoint", BCA_p); - DD(BCA_p); - MaskImagePointType bif1; - MaskImagePointType bif2; - TrackBifurcationFromPoint(r, BCA_p, bif1, bif2); - DD(bif1); - DD(bif2); - */ - // Find max label + // Track the SoughtVessel from the given first point + // superiorly. This is done by TrackBifurcationFromPoint + MaskImagePointType SoughtVesselSeedPoint; + GetAFDB()->GetPoint3D(m_SoughtVesselSeedName, SoughtVesselSeedPoint); + + // Find the label with the maximum value to set the result typedef itk::MinimumMaximumImageCalculator MinMaxFilterType; MinMaxFilterType::Pointer ff = MinMaxFilterType::New(); ff->SetImage(r); ff->ComputeMaximum(); LabelType newLabel = ff->GetMaximum()+1; - DD(newLabel); - // Get all centroids of the first slice - std::vector centroids2D; - clitk::ComputeCentroids(m_slice_recon[0], GetBackgroundValue(), centroids2D); - DD(centroids2D.size()); + // the following bifurcations point will the centroids of the + // components obtain when (hopefully!) the SoughtVessel + // split into CommonArtery and SubclavianArtery. std::vector bifurcations; - clitk::PointsUtils::Convert2DListTo3DList(centroids2D, 0, r, bifurcations); - DD(bifurcations.size()); - for(uint i=1; i(m_slice_recon, m_Mask, 2); - clitk::writeImage(rr, "recon"+toString(i)+".mhd"); - } + TrackBifurcationFromPoint(r, m_slice_recon, SoughtVesselSeedPoint, newLabel, bifurcations); + + // Build the final 3D image from the previous slice by slice processing + m_SoughtVessel = clitk::JoinSlices(m_slice_recon, m_Mask, 2); + writeImage(m_SoughtVessel, "recon2.mhd"); + + // Set binary image, (remove other labels). + // TODO: keep labeled image to track SubclavianArtery and CommonArtery + m_SoughtVessel = + clitk::Binarize(m_SoughtVessel, newLabel, newLabel, + GetBackgroundValue(), GetForegroundValue()); + + writeImage(m_SoughtVessel, "afterbinarize.mhd"); + + m_SoughtVessel = clitk::AutoCrop(m_SoughtVessel, GetBackgroundValue()); + writeImage(m_SoughtVessel, "afterautocrop.mhd"); + + // Clean the image : Opening (not in Z direction) + typename MaskImageType::SizeType rad; + rad[0] = rad[1] = 2; + rad[2] = 0; + m_SoughtVessel = clitk::Opening(m_SoughtVessel, rad, + GetBackgroundValue(), GetForegroundValue()); + + writeImage(m_SoughtVessel, "afteropen.mhd"); + + // Clean the image : keep main CCL slice by slice + m_SoughtVessel = clitk::SliceBySliceKeepMainCCL(m_SoughtVessel, + GetBackgroundValue(), + GetForegroundValue()); } //-------------------------------------------------------------------- @@ -216,11 +225,11 @@ template void clitk::ExtractMediastinalVesselsFilter:: GenerateData() { - DD("GenerateData"); - // Final Step -> graft output (if SetNthOutput => redo) - MaskImagePointer BrachioCephalicArtery = - GetAFDB()->template GetImage("BrachioCephalicArtery"); - this->GraftNthOutput(0, BrachioCephalicArtery); + // Save in the AFDB (not write on the disk here) + GetAFDB()->SetImageFilename(GetSoughtVesselName(), GetOutputFilename()); + WriteAFDB(); + // Final Step -> graft output + this->GraftNthOutput(0, m_SoughtVessel); } //-------------------------------------------------------------------- @@ -229,40 +238,52 @@ GenerateData() { template void clitk::ExtractMediastinalVesselsFilter:: -CropSupInf() { - StartNewStep("Inf/Sup limits (carina) and crop with mediastinum"); +CropInputImage() { + StartNewStep("Crop the input image: SI,AP limits with carina and crop with mediastinum"); + /* + Need : Trachea, Carina (roi not point), + */ // Get Trachea and Carina MaskImagePointer Trachea = GetAFDB()->template GetImage ("Trachea"); - // Get or compute Carina + // Compute Carina position double m_CarinaZ; - // Get Carina Z position MaskImagePointer Carina = GetAFDB()->template GetImage("Carina"); - std::vector centroids; clitk::ComputeCentroids(Carina, GetBackgroundValue(), centroids); m_CarinaZ = centroids[1][2]; - // DD(m_CarinaZ); - // add one slice to include carina ? + // add one slice to include carina m_CarinaZ += Carina->GetSpacing()[2]; // We dont need Carina structure from now Carina->Delete(); - GetAFDB()->SetDouble("CarinaZ", m_CarinaZ); + GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]); - // Crop - m_Input = clitk::CropImageRemoveLowerThan(m_Input, 2, - m_CarinaZ, false, GetBackgroundValue()); + // Crop Inf, remove below Carina + m_Input = + clitk::CropImageRemoveLowerThan(m_Input, 2, m_CarinaZ, false, GetBackgroundValue()); - // Crop not post to centroid + // Crop post double m_CarinaY = centroids[1][1]; - DD(m_CarinaY); - m_Input = clitk::CropImageRemoveGreaterThan(m_Input, 1, // OLD ABOVE - m_CarinaY, false, GetBackgroundValue()); - // Crop not ant to centroid + m_Input = clitk::CropImageRemoveGreaterThan(m_Input, 1, + m_CarinaY+GetMaxDistancePostToCarina(), + false, GetBackgroundValue()); + // Crop ant m_Input = clitk::CropImageRemoveLowerThan(m_Input, 1, - m_CarinaY-80, false, GetBackgroundValue()); - - // AutoCrop with Mediastinum + m_CarinaY-GetMaxDistanceAntToCarina(), + false, GetBackgroundValue()); + + // Crop Right + double m_CarinaX = centroids[1][0]; + m_Input = clitk::CropImageRemoveLowerThan(m_Input, 0, + m_CarinaX-GetMaxDistanceRightToCarina(), + false, GetBackgroundValue()); + // Crop Left + m_Input = clitk::CropImageRemoveGreaterThan(m_Input, 0, + m_CarinaX+GetMaxDistanceLeftToCarina(), + false, GetBackgroundValue()); + + /* + // AutoCrop with Mediastinum, generaly only allow to remove few slices (superiorly) m_Mediastinum = GetAFDB()->template GetImage("Mediastinum"); // Resize like input (sup to carina) m_Mediastinum = clitk::ResizeImageLike(m_Mediastinum, m_Input, GetBackgroundValue()); @@ -270,6 +291,7 @@ CropSupInf() { m_Mediastinum = clitk::AutoCrop(m_Mediastinum, GetBackgroundValue()); // Resize input m_Input = clitk::ResizeImageLike(m_Input, m_Mediastinum, GetBackgroundValue()); + */ // End StopCurrentStep(m_Input); @@ -281,39 +303,36 @@ CropSupInf() { template void clitk::ExtractMediastinalVesselsFilter:: -//SearchBrachioCephalicArtery(int & BCA_first_slice, LabelType & BCA_first_label) { TrackBifurcationFromPoint(MaskImagePointer & recon, std::vector & slices_recon, MaskImagePointType point3D, LabelType newLabel, std::vector & bifurcations) { - StartNewStep("Search for BCA first slice and label"); - DD(newLabel); + StartNewStep("Track the SoughtVessel from the seed point"); - // Extract slices - // std::vector slices_recon; - //clitk::ExtractSlices(recon, 2, slices_recon); - - // Find first slice + // Find first slice index MaskImageIndexType index; recon->TransformPhysicalPointToIndex(point3D, index); - DD(point3D); - DD(index); + int numberOfBifurcation = 0; + typedef typename MaskSliceType::PointType SlicePointType; + SlicePointType previousCenter; - uint i=index[2]; + // Get current label at the point3D of interest + uint currentSlice=index[2]; bool found = false; - LabelType previous_largest_label=recon->GetPixel(index); - DD(previous_largest_label); + LabelType previous_slice_label=recon->GetPixel(index); + // DD(slices_recon.size()); do { - DD(i); + // DD(currentSlice); // Consider current reconstructed slice - MaskSlicePointer s = slices_recon[i]; + MaskSlicePointer s = slices_recon[currentSlice]; MaskSlicePointer previous; - if (i==index[2]) previous = s; - else previous = slices_recon[i-1]; + if (currentSlice == index[2]) previous = s; + else { + previous = slices_recon[currentSlice-1]; + } // Get centroids of the labels in the current slice - typedef typename MaskSliceType::PointType SlicePointType; static const unsigned int Dim = MaskSliceType::ImageDimension; typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType; typedef itk::LabelMap< LabelObjectType > LabelMapType; @@ -330,73 +349,128 @@ TrackBifurcationFromPoint(MaskImagePointer & recon, // Look what centroid inside the previous largest one std::vector centroids; std::vector centroids_label; + std::vector labels_size; for(uint c=0; cGetNumberOfLabelObjects(); c++) { int label = labelMap->GetLabels()[c]; - DD(label); + // DD(label); SlicePointType center = labelMap->GetLabelObject(label)->GetCentroid(); - SlicePointType center_previous = center; - center_previous[2] -= m_Input->GetSpacing()[2]; + // DD(center); // Get label into previous slice - typename MaskSliceType::IndexType index; - previous->TransformPhysicalPointToIndex(center_previous, index); - LabelType l = previous->GetPixel(index); - DD(l); - if (l == previous_largest_label) { + typename MaskSliceType::IndexType centerIndex; + previous->TransformPhysicalPointToIndex(center, centerIndex); + LabelType labelInPreviousSlice = previous->GetPixel(centerIndex); + // if this current centroid was in the current label, add it to bifurcations + if (labelInPreviousSlice == previous_slice_label) { centroids.push_back(center); centroids_label.push_back(label); + labels_size.push_back(labelMap->GetLabelObject(label)->GetPhysicalSize()); + } + } + + if (centroids.size() == 0) { + // Last attempt to find -> check if previous centroid is inside a CCL + // if in s -> get value, getcentroid add. + DD(currentSlice); + DD("Last change to find"); + typename MaskSliceType::IndexType previousCenterIndex; + s->TransformPhysicalPointToIndex(previousCenter, previousCenterIndex); + DD(previousCenter); + LabelType labelInSlice = s->GetPixel(previousCenterIndex); + DD(labelInSlice); + if (labelInSlice != GetBackgroundValue()) { + centroids.push_back(labelMap->GetLabelObject(labelInSlice)->GetCentroid()); + centroids_label.push_back(labelInSlice); + labels_size.push_back(labelMap->GetLabelObject(labelInSlice)->GetPhysicalSize()); } } - DD(centroids.size()); + + + // DD(centroids.size()); // If several centroids, we found a bifurcation if (centroids.size() > 1) { - found = true; - for(uint c=0; c::Convert2DTo3D(centroids[c], m_Mask, i, bif); - bifurcations.push_back(bif); - s = clitk::SetBackground(s, s, centroids_label[c], newLabel+c+1, true); - // slices_recon[i] = s; // (useful ?) + numberOfBifurcation++; + // If the number of bifurcation is greater than the required one, we stop + if (numberOfBifurcation > GetMaxNumberOfFoundBifurcation()) { + found = true; + DD("max bif reach"); + for(uint c=0; c::Convert2DTo3D(centroids[c], m_Mask, currentSlice, bif); + bifurcations.push_back(bif); + } + } + // Else we continue along the main (largest) connected component + else { + int indexOfLargest = 0; + for(uint b=0; b labels_size[indexOfLargest]) { + indexOfLargest = b; + } + } + SlicePointType c = centroids[indexOfLargest]; + LabelType l = centroids_label[indexOfLargest]; + centroids.clear(); + centroids.push_back(c); + centroids_label.push_back(l); } - DD("FOUND"); } // if only one centroids, we change the current image with the current label if (centroids.size() == 1) { s = clitk::SetBackground(s, s, centroids_label[0], newLabel, true); - previous_largest_label = newLabel; - /*typedef itk::BinaryThresholdImageFilter BinarizeFilterType; - typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New(); - binarizeFilter->SetInput(s); - binarizeFilter->SetLowerThreshold(centroids_label[0]); - binarizeFilter->SetUpperThreshold(centroids_label[0]+1); - binarizeFilter->SetInsideValue(previous_largest_label); - binarizeFilter->SetOutsideValue(GetBackgroundValue()); - binarizeFilter->Update(); - s = binarizeFilter->GetOutput();*/ - slices_recon[i] = s; // (not useful ?) + slices_recon[currentSlice] = s; + previous_slice_label = newLabel; + // It can happend that several CCL share this same label. To + // prevent this case, we only consider the one that contains + // the centroid. + MaskSlicePointer temp = clitk::Binarize(s, newLabel, newLabel, GetBackgroundValue(), GetForegroundValue()); + // writeImage(temp, "bin-"+toString(currentSlice)+".mhd"); + temp = clitk::Labelize(temp, GetBackgroundValue(), true, 1); + //writeImage(temp, "label-"+toString(currentSlice)+".mhd"); + typename MaskSliceType::IndexType centroids_index; + temp->TransformPhysicalPointToIndex(centroids[0], centroids_index); + typename MaskSliceType::PixelType v = temp->GetPixel(centroids_index); + + // It can happend that the centroid is inside the BG, so we keep + // the largest CCL (the first); + if (v == GetBackgroundValue()) { + DD(currentSlice); + DD("inside BG"); + DD(centroids[0]); + v = 1; // largest one + } + + //DD(v); + temp = clitk::Binarize(temp, v, v, GetBackgroundValue(), newLabel); + //writeImage(temp, "relabel-"+toString(currentSlice)+".mhd"); + s = temp; + slices_recon[currentSlice] = s; + + // I need to recompute the centroid if we have removed some + // connected component. + clitk::ComputeCentroids(s, GetBackgroundValue(), centroids); + previousCenter = centroids[1]; } if (centroids.size() == 0) { - DD("no centroid, I stop"); + DD("ZERO"); found = true; } - if (i == slices_recon.size()-1) found = true; + if (currentSlice == slices_recon.size()-1) { + DD("end of slices"); + found = true; + } // iterate - ++i; + ++currentSlice; } while (!found); - - //MaskImageType::Pointer rr = clitk::JoinSlices(slices_recon, m_Mask, 2); - //clitk::writeImage(rr, "recon2.mhd"); - // End StopCurrentStep(); } //-------------------------------------------------------------------- - #endif //#define CLITKEXTRACTMEDIASTINALVESSELSFILTER_TXX diff --git a/segmentation/clitkExtractMediastinalVesselsGenericFilter.txx b/segmentation/clitkExtractMediastinalVesselsGenericFilter.txx index b85da27..6344344 100644 --- a/segmentation/clitkExtractMediastinalVesselsGenericFilter.txx +++ b/segmentation/clitkExtractMediastinalVesselsGenericFilter.txx @@ -51,7 +51,6 @@ void clitk::ExtractMediastinalVesselsGenericFilter::SetArgsInfo(co SetIOVerbose(mArgsInfo.verbose_flag); if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes(); if (mArgsInfo.input_given) AddInputFilename(mArgsInfo.input_arg); - if (mArgsInfo.output_given) AddOutputFilename(mArgsInfo.output_arg); } //-------------------------------------------------------------------- @@ -68,9 +67,24 @@ SetOptionsFromArgsInfoToFilter(FilterType * f) f->SetWriteStepFlag(mArgsInfo.writeStep_flag); f->SetVerboseMemoryFlag(mArgsInfo.verboseMemory_flag); f->SetAFDBFilename(mArgsInfo.afdb_arg); - f->SetOutputFolder(mArgsInfo.output_arg); - f->SetThreshold(mArgsInfo.threshold_arg); + f->SetThresholdHigh(mArgsInfo.thresholdHigh_arg); + f->SetThresholdLow(mArgsInfo.thresholdLow_arg); + f->SetErosionRadius(mArgsInfo.erode_arg); + f->SetDilatationRadius(mArgsInfo.dilate_arg); + + f->SetMaxDistancePostToCarina(mArgsInfo.maxPost_arg); + f->SetMaxDistanceAntToCarina(mArgsInfo.maxAnt_arg); + f->SetMaxDistanceLeftToCarina(mArgsInfo.maxLeft_arg); + f->SetMaxDistanceRightToCarina(mArgsInfo.maxRight_arg); + + f->SetSoughtVesselSeedName(mArgsInfo.seed_arg); + f->SetSoughtVesselName(mArgsInfo.name_arg); + f->SetMaxNumberOfFoundBifurcation(mArgsInfo.bif_arg); + + // Output filename + this->AddOutputFilename(mArgsInfo.output_arg); + f->SetOutputFilename(mArgsInfo.output_arg); } //-------------------------------------------------------------------- diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index cc56a6d..2ccbb6e 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -115,9 +115,9 @@ IF (CLITK_BUILD_TOOLS) ADD_EXECUTABLE(clitkFooImage clitkFooImage.cxx ${clitkFooImage_GGO_C}) TARGET_LINK_LIBRARIES(clitkFooImage clitkCommon ${ITK_LIBRARIES} ) - WRAP_GGO(clitkMedianImageFilter_GGO_C clitkMedianImageFilter.ggo) - ADD_EXECUTABLE(clitkMedianImageFilter clitkMedianImageFilter.cxx ${clitkMedianImageFilter_GGO_C}) - TARGET_LINK_LIBRARIES(clitkMedianImageFilter clitkCommon ${ITK_LIBRARIES}) + #WRAP_GGO(clitkMedianImageFilter_GGO_C clitkMedianImageFilter.ggo) + #ADD_EXECUTABLE(clitkMedianImageFilter clitkMedianImageFilter.cxx ${clitkMedianImageFilter_GGO_C}) + #TARGET_LINK_LIBRARIES(clitkMedianImageFilter clitkMedianImageFilterLib clitkCommon ${ITK_LIBRARIES}) ADD_EXECUTABLE(clitkResampleImage clitkResampleImage.cxx) TARGET_LINK_LIBRARIES(clitkResampleImage clitkResampleImageLib clitkCommon ${ITK_LIBRARIES}) @@ -126,9 +126,9 @@ IF (CLITK_BUILD_TOOLS) ADD_EXECUTABLE(clitkMinMaxMask clitkMinMaxMask.cxx ${clitkMinMaxMask_GGO_C}) TARGET_LINK_LIBRARIES(clitkMinMaxMask clitkCommon ${ITK_LIBRARIES} ) - # WRAP_GGO(clitkAutoCrop_GGO_C clitkAutoCrop.ggo) - # ADD_EXECUTABLE(clitkAutoCrop clitkAutoCrop.cxx ${clitkAutoCrop_GGO_C}) - # TARGET_LINK_LIBRARIES(clitkAutoCrop clitkCommon ${ITK_LIBRARIES} ) + WRAP_GGO(clitkAutoCrop_GGO_C clitkAutoCrop.ggo) + ADD_EXECUTABLE(clitkAutoCrop clitkAutoCrop.cxx ${clitkAutoCrop_GGO_C}) + TARGET_LINK_LIBRARIES(clitkAutoCrop clitkCommon ${ITK_LIBRARIES} ) WRAP_GGO(clitkDicomRTStruct2BinaryImage_GGO_C clitkDicomRTStruct2BinaryImage.ggo) ADD_EXECUTABLE(clitkDicomRTStruct2BinaryImage clitkDicomRTStruct2BinaryImage.cxx ${clitkDicomRTStruct2BinaryImage_GGO_C}) diff --git a/tools/clitkDicomRTStruct2BinaryImage.cxx b/tools/clitkDicomRTStruct2BinaryImage.cxx index d54c47e..97eff97 100644 --- a/tools/clitkDicomRTStruct2BinaryImage.cxx +++ b/tools/clitkDicomRTStruct2BinaryImage.cxx @@ -30,7 +30,9 @@ int main(int argc, char * argv[]) { // Read and display information clitk::DicomRT_StructureSet::Pointer s = clitk::DicomRT_StructureSet::New(); s->Read(args_info.input_arg); - // s->Print(std::cout); + if (args_info.verboseFile_flag) { + s->Print(std::cout); + } // New filter to convert to binary image clitk::DicomRT_ROI_ConvertToImageFilter filter; diff --git a/tools/clitkDicomRTStruct2BinaryImage.ggo b/tools/clitkDicomRTStruct2BinaryImage.ggo index 5fdd741..d403f06 100644 --- a/tools/clitkDicomRTStruct2BinaryImage.ggo +++ b/tools/clitkDicomRTStruct2BinaryImage.ggo @@ -4,6 +4,7 @@ version "Convert DICOM RT Structure Set (contours) to binary image" option "config" - "Config file" string no option "verbose" v "Verbose" flag off +option "verboseFile" - "Verbose file content" flag off option "input" i "Input Dicom file" string yes option "image" j "Used to read image info (spacing, origin)" string yes option "roi" r "ROI to binarize (if -1 = all roi)" int no default="-1" diff --git a/vv/vvMesh.cxx b/vv/vvMesh.cxx index 72ed60d..5b11d48 100644 --- a/vv/vvMesh.cxx +++ b/vv/vvMesh.cxx @@ -141,6 +141,7 @@ void vvMesh::ComputeMasks(vtkImageData* sample,bool extrude) double * samp_origin=sample->GetOrigin(); double * spacing=sample->GetSpacing(); binary_image->SetSpacing(spacing); + /// Put the origin on a voxel to avoid small skips binary_image->SetOrigin(floor((bounds[0]-samp_origin[0])/spacing[0]-2)*spacing[0]+samp_origin[0], floor((bounds[2]-samp_origin[1])/spacing[1]-2)*spacing[1]+samp_origin[1], @@ -173,10 +174,13 @@ void vvMesh::ComputeMasks(vtkImageData* sample,bool extrude) stencil->SetInput(binary_image); stencil->Update(); this->AddMask(stencil->GetOutput()); - //vtkSmartPointer w = vtkSmartPointer::New(); - //w->SetInput(stencil->GetOutput()); - //w->SetFileName("binary.mhd"); - //w->Write(); + + /* + vtkSmartPointer w = vtkSmartPointer::New(); + w->SetInput(stencil->GetOutput()); + w->SetFileName("binary.mhd"); + w->Write(); + */ } }