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);
// -----------------------------------------------------
// 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<MaskImageType>("BrachioCephalicTrunk");
- MaskImagePointType p = BrachioCephalicTrunk->GetOrigin(); // initialise to avoid warning
- clitk::FindExtremaPointInAGivenDirection<MaskImageType>(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<MaskImageType>("BrachioCephalicArtery");
+ MaskImagePointType p = BrachioCephalicArtery->GetOrigin(); // initialise to avoid warning
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(BrachioCephalicArtery, GetBackgroundValue(), 2, false, p);
+ double TopOfBrachioCephalicArteryZ=p[2] + BrachioCephalicArtery->GetSpacing()[2]; // Add one slice
// Remove CommonCarotidArtery below this point
- CommonCarotidArtery = clitk::CropImageRemoveLowerThan<MaskImageType>(CommonCarotidArtery, 2, TopOfBrachioCephalicTrunkZ, true, GetBackgroundValue());
+ CommonCarotidArtery = clitk::CropImageRemoveLowerThan<MaskImageType>(CommonCarotidArtery, 2, TopOfBrachioCephalicArteryZ, true, GetBackgroundValue());
// Find most Ant points
std::vector<MaskImagePointType> ccaAntPositionA;
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<MaskImageType>(m_Working_Support, BrachioCephalicTrunk, 2,
- GetFuzzyThreshold("2RL", "BrachioCephalicTrunk"), "NotAntTo", false, 2, true, false);
+ clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, BrachioCephalicArtery, 2,
+ GetFuzzyThreshold("2RL", "BrachioCephalicArtery"), "NotAntTo", false, 2, true, false);
// End
StopCurrentStep<MaskImageType>(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<MaskImagePointType> bctAntPositionA;
std::vector<MaskImagePointType> bctAntPositionB;
- clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition<MaskImageType>(BrachioCephalicTrunk,
+ clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition<MaskImageType>(BrachioCephalicArtery,
GetBackgroundValue(), 2,
1, true, // Ant
0, // Horizontal line
/* 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
*/
// Read structures
- MaskImagePointer BrachioCephalicTrunk = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicTrunk");
+ MaskImagePointer BrachioCephalicArtery = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
MaskImagePointer CommonCarotidArtery = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArtery");
MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
// Resize all structures like support
- BrachioCephalicTrunk =
- clitk::ResizeImageLike<MaskImageType>(BrachioCephalicTrunk, m_Working_Support, GetBackgroundValue());
+ BrachioCephalicArtery =
+ clitk::ResizeImageLike<MaskImageType>(BrachioCephalicArtery, m_Working_Support, GetBackgroundValue());
CommonCarotidArtery =
clitk::ResizeImageLike<MaskImageType>(CommonCarotidArtery, m_Working_Support, GetBackgroundValue());
SubclavianArtery =
clitk::ResizeImageLike<MaskImageType>(Trachea, m_Working_Support, GetBackgroundValue());
// Extract slices
- std::vector<MaskSlicePointer> slices_BrachioCephalicTrunk;
- clitk::ExtractSlices<MaskImageType>(BrachioCephalicTrunk, 2, slices_BrachioCephalicTrunk);
+ std::vector<MaskSlicePointer> slices_BrachioCephalicArtery;
+ clitk::ExtractSlices<MaskImageType>(BrachioCephalicArtery, 2, slices_BrachioCephalicArtery);
std::vector<MaskSlicePointer> slices_BrachioCephalicVein;
clitk::ExtractSlices<MaskImageType>(BrachioCephalicVein, 2, slices_BrachioCephalicVein);
std::vector<MaskSlicePointer> slices_CommonCarotidArtery;
clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices_Aorta);
std::vector<MaskSlicePointer> slices_Trachea;
clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices_Trachea);
- uint n= slices_BrachioCephalicTrunk.size();
+ uint n= slices_BrachioCephalicArtery.size();
// Get the boundaries of one slice
std::vector<MaskSlicePointType> bounds;
- ComputeImageBoundariesCoordinates<MaskSliceType>(slices_BrachioCephalicTrunk[0], bounds);
+ ComputeImageBoundariesCoordinates<MaskSliceType>(slices_BrachioCephalicArtery[0], bounds);
// For all slices, for all structures, find the centroid and build the contour
// List of 3D points (for debug)
GetBackgroundValue(), true, 1);
slices_SubclavianArtery[i] = Labelize<MaskSliceType>(slices_SubclavianArtery[i],
GetBackgroundValue(), true, 1);
- slices_BrachioCephalicTrunk[i] = Labelize<MaskSliceType>(slices_BrachioCephalicTrunk[i],
+ slices_BrachioCephalicArtery[i] = Labelize<MaskSliceType>(slices_BrachioCephalicArtery[i],
GetBackgroundValue(), true, 1);
slices_BrachioCephalicVein[i] = Labelize<MaskSliceType>(slices_BrachioCephalicVein[i],
GetBackgroundValue(), true, 1);
std::vector<MaskSlicePointType> centroids6;
ComputeCentroids<MaskSliceType>(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1);
ComputeCentroids<MaskSliceType>(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2);
- ComputeCentroids<MaskSliceType>(slices_BrachioCephalicTrunk[i], GetBackgroundValue(), centroids3);
+ ComputeCentroids<MaskSliceType>(slices_BrachioCephalicArtery[i], GetBackgroundValue(), centroids3);
ComputeCentroids<MaskSliceType>(slices_Thyroid[i], GetBackgroundValue(), centroids4);
ComputeCentroids<MaskSliceType>(slices_Aorta[i], GetBackgroundValue(), centroids5);
ComputeCentroids<MaskSliceType>(slices_BrachioCephalicVein[i], GetBackgroundValue(), centroids6);
}
// 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();
}
// clitk
#include "clitkCommon.h"
#include "clitkExtractMediastinalVesselsFilter.h"
-#include "clitkAutoCropFilter.h"
#include "clitkSegmentationUtils.h"
-#include "clitkMorphoMathFilter.h"
+#include "clitkReconstructWithConditionalGrayscaleDilateImageFilter.h"
// itk
#include <itkBinaryThresholdImageFilter.h>
-#include <itkGrayscaleDilateImageFilter.h>
#include <itkMinimumMaximumImageCalculator.h>
//--------------------------------------------------------------------
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");
}
//--------------------------------------------------------------------
m_Input = dynamic_cast<ImageType*>(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<ImageType, MaskImageType> 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<MaskImageType>(m_Mask, "m.mhd");
StopCurrentStep<MaskImageType>(m_Mask);
- // Keep main CCL ?
- m_Mask = clitk::Labelize<MaskImageType>(m_Mask, GetBackgroundValue(), false, 10);
- m_Mask = KeepLabels<MaskImageType>(m_Mask, GetBackgroundValue(), GetTemporaryForegroundValue(), 1, 1, true);
- clitk::writeImage<MaskImageType>(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<MaskImageType>(m_Mask2);
+
// ------------------------------------------------------------------
// Extract slices
- StartNewStep("Detect vessels (slice by slice reconstruction)");
+ StartNewStep("Detect objects : erosion, then slice by slice reconstruction");
std::vector<MaskSlicePointer> slices_mask;
clitk::ExtractSlices<MaskImageType>(m_Mask, 2, slices_mask);
-
- DD(slices_mask.size());
-
+ std::vector<MaskSlicePointer> slices_mask2;
+ clitk::ExtractSlices<MaskImageType>(m_Mask2, 2, slices_mask2);
+ int radius = GetErosionRadius();
+
+ // List of working slices (debug only)
std::vector<MaskSlicePointer> debug_eroded;
- std::vector<MaskSlicePointer> debug_labeled;
- std::vector<MaskSlicePointer> debug_slabeled;
+ // std::vector<MaskSlicePointer> debug_labeled;
- int radius = 3;
- DD(radius); // TO PUT IN OPTION
-
// ------------------------------------------------------------------
- // Loop Slice by Slice -> erode find CCL and reconstruct
- clitk::MorphoMathFilter<MaskSliceType>::Pointer f= clitk::MorphoMathFilter<MaskSliceType>::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; i<slices_mask.size(); i++) {
- // Erosion
- f->SetInput(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<MaskSliceType>(eroded, "er-"+toString(i)+".mhd");
+ // Erosion kernel
+ typedef itk::BinaryBallStructuringElement<MaskSliceType::PixelType,2> KernelType;
+ KernelType structuringElement;
+ structuringElement.SetRadius(radius);
+ structuringElement.CreateStructuringElement();
+
+ // Erosion -> we break the connectivity between structure
+ typedef itk::BinaryErodeImageFilter<MaskSliceType, MaskSliceType, KernelType> 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<MaskSliceType>(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<MaskSliceType, MaskSliceType>(labeled, labeled, li, lo, true);
- }
- debug_labeled.push_back(labeled);
-
- // Create kernel for GrayscaleDilateImageFilter
- typedef itk::BinaryBallStructuringElement<MaskSliceType::PixelType,MaskSliceType::ImageDimension > KernelType;
- KernelType k;
- k.SetRadius(radius+1);
- k.CreateStructuringElement();
-
- // Keep the MAX -> we prefer the opposite su change the label
- typedef itk::GrayscaleDilateImageFilter<MaskSliceType, MaskSliceType, KernelType> 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<MaskSliceType, MaskSliceType>(s, slices_mask[i],
- GetBackgroundValue(), GetBackgroundValue(), true);
-
+ clitk::Labelize<MaskSliceType>(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<MaskSliceType> 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<MaskSliceType, MaskSliceType>(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<MaskImageType>(debug_eroded, m_Mask, 2);
- clitk::writeImage<MaskImageType>(eroded, "eroded.mhd");
-
- DD("l");
- MaskImageType::Pointer l = clitk::JoinSlices<MaskImageType>(debug_labeled, m_Mask, 2);
- clitk::writeImage<MaskImageType>(l, "labeled.mhd");
-
- DD("r");
+ writeImage<MaskImageType>(eroded, "erode.mhd");
+ //MaskImageType::Pointer l = clitk::JoinSlices<MaskImageType>(debug_labeled, m_Mask, 2);
MaskImageType::Pointer r = clitk::JoinSlices<MaskImageType>(m_slice_recon, m_Mask, 2);
- clitk::writeImage<MaskImageType>(r, "recon.mhd");
+ writeImage<MaskImageType>(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<MaskImageType> 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<MaskSlicePointType> centroids2D;
- clitk::ComputeCentroids<MaskSliceType>(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<MaskImagePointType> bifurcations;
- clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(centroids2D, 0, r, bifurcations);
- DD(bifurcations.size());
- for(uint i=1; i<bifurcations.size()+1; i++) {
- DD(i);
- DD(bifurcations.size());
- TrackBifurcationFromPoint(r, m_slice_recon, bifurcations[i], newLabel+i, bifurcations);
- DD("end track");
- DD(bifurcations.size());
- MaskImageType::Pointer rr = clitk::JoinSlices<MaskImageType>(m_slice_recon, m_Mask, 2);
- clitk::writeImage<MaskImageType>(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<MaskImageType>(m_slice_recon, m_Mask, 2);
+ writeImage<MaskImageType>(m_SoughtVessel, "recon2.mhd");
+
+ // Set binary image, (remove other labels).
+ // TODO: keep labeled image to track SubclavianArtery and CommonArtery
+ m_SoughtVessel =
+ clitk::Binarize<MaskImageType>(m_SoughtVessel, newLabel, newLabel,
+ GetBackgroundValue(), GetForegroundValue());
+
+ writeImage<MaskImageType>(m_SoughtVessel, "afterbinarize.mhd");
+
+ m_SoughtVessel = clitk::AutoCrop<MaskImageType>(m_SoughtVessel, GetBackgroundValue());
+ writeImage<MaskImageType>(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<MaskImageType>(m_SoughtVessel, rad,
+ GetBackgroundValue(), GetForegroundValue());
+
+ writeImage<MaskImageType>(m_SoughtVessel, "afteropen.mhd");
+
+ // Clean the image : keep main CCL slice by slice
+ m_SoughtVessel = clitk::SliceBySliceKeepMainCCL<MaskImageType>(m_SoughtVessel,
+ GetBackgroundValue(),
+ GetForegroundValue());
}
//--------------------------------------------------------------------
void
clitk::ExtractMediastinalVesselsFilter<TImageType>::
GenerateData() {
- DD("GenerateData");
- // Final Step -> graft output (if SetNthOutput => redo)
- MaskImagePointer BrachioCephalicArtery =
- GetAFDB()->template GetImage<MaskImageType>("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);
}
//--------------------------------------------------------------------
template <class TImageType>
void
clitk::ExtractMediastinalVesselsFilter<TImageType>::
-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 <MaskImageType>("Trachea");
- // Get or compute Carina
+ // Compute Carina position
double m_CarinaZ;
- // Get Carina Z position
MaskImagePointer Carina = GetAFDB()->template GetImage<MaskImageType>("Carina");
-
std::vector<MaskImagePointType> centroids;
clitk::ComputeCentroids<MaskImageType>(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<ImageType>(m_Input, 2,
- m_CarinaZ, false, GetBackgroundValue());
+ // Crop Inf, remove below Carina
+ m_Input =
+ clitk::CropImageRemoveLowerThan<ImageType>(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<ImageType>(m_Input, 1, // OLD ABOVE
- m_CarinaY, false, GetBackgroundValue());
- // Crop not ant to centroid
+ m_Input = clitk::CropImageRemoveGreaterThan<ImageType>(m_Input, 1,
+ m_CarinaY+GetMaxDistancePostToCarina(),
+ false, GetBackgroundValue());
+ // Crop ant
m_Input = clitk::CropImageRemoveLowerThan<ImageType>(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<ImageType>(m_Input, 0,
+ m_CarinaX-GetMaxDistanceRightToCarina(),
+ false, GetBackgroundValue());
+ // Crop Left
+ m_Input = clitk::CropImageRemoveGreaterThan<ImageType>(m_Input, 0,
+ m_CarinaX+GetMaxDistanceLeftToCarina(),
+ false, GetBackgroundValue());
+
+ /*
+ // AutoCrop with Mediastinum, generaly only allow to remove few slices (superiorly)
m_Mediastinum = GetAFDB()->template GetImage<MaskImageType>("Mediastinum");
// Resize like input (sup to carina)
m_Mediastinum = clitk::ResizeImageLike<MaskImageType>(m_Mediastinum, m_Input, GetBackgroundValue());
m_Mediastinum = clitk::AutoCrop<MaskImageType>(m_Mediastinum, GetBackgroundValue());
// Resize input
m_Input = clitk::ResizeImageLike<ImageType>(m_Input, m_Mediastinum, GetBackgroundValue());
+ */
// End
StopCurrentStep<ImageType>(m_Input);
template <class TImageType>
void
clitk::ExtractMediastinalVesselsFilter<TImageType>::
-//SearchBrachioCephalicArtery(int & BCA_first_slice, LabelType & BCA_first_label) {
TrackBifurcationFromPoint(MaskImagePointer & recon,
std::vector<MaskSlicePointer> & slices_recon,
MaskImagePointType point3D,
LabelType newLabel,
std::vector<MaskImagePointType> & bifurcations) {
- StartNewStep("Search for BCA first slice and label");
- DD(newLabel);
+ StartNewStep("Track the SoughtVessel from the seed point");
- // Extract slices
- // std::vector<MaskSlicePointer> slices_recon;
- //clitk::ExtractSlices<MaskImageType>(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;
// Look what centroid inside the previous largest one
std::vector<SlicePointType> centroids;
std::vector<LabelType> centroids_label;
+ std::vector<double> labels_size;
for(uint c=0; c<labelMap->GetNumberOfLabelObjects(); 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<centroids.size(); c++) {
- ImagePointType bif;
- clitk::PointsUtils<MaskImageType>::Convert2DTo3D(centroids[c], m_Mask, i, bif);
- bifurcations.push_back(bif);
- s = clitk::SetBackground<MaskSliceType, MaskSliceType>(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<centroids.size(); c++) {
+ ImagePointType bif;
+ clitk::PointsUtils<MaskImageType>::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<centroids.size(); b++) {
+ if (labels_size[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<MaskSliceType, MaskSliceType>(s, s, centroids_label[0], newLabel, true);
- previous_largest_label = newLabel;
- /*typedef itk::BinaryThresholdImageFilter<MaskSliceType, MaskSliceType> 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<MaskSliceType>(s, newLabel, newLabel, GetBackgroundValue(), GetForegroundValue());
+ // writeImage<MaskSliceType>(temp, "bin-"+toString(currentSlice)+".mhd");
+ temp = clitk::Labelize<MaskSliceType>(temp, GetBackgroundValue(), true, 1);
+ //writeImage<MaskSliceType>(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<MaskSliceType>(temp, v, v, GetBackgroundValue(), newLabel);
+ //writeImage<MaskSliceType>(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<MaskSliceType>(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<MaskImageType>(slices_recon, m_Mask, 2);
- //clitk::writeImage<MaskImageType>(rr, "recon2.mhd");
-
// End
StopCurrentStep();
}
//--------------------------------------------------------------------
-
#endif //#define CLITKEXTRACTMEDIASTINALVESSELSFILTER_TXX
git://git.creatis.insa-lyon.fr / clitk.git / commitdiff