X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=segmentation%2FclitkExtractMediastinalVesselsFilter.txx;h=178194871e3b536e5b8755ec71167248dee8d40a;hb=15e781be4b87a75be8fe01e56da5eb4435c97786;hp=02576f498289c8ef0b3e04e9dadb6fd94fee64ca;hpb=22f9a65075ca86a805d06f3a1a65151150d570d3;p=clitk.git diff --git a/segmentation/clitkExtractMediastinalVesselsFilter.txx b/segmentation/clitkExtractMediastinalVesselsFilter.txx index 02576f4..1781948 100644 --- a/segmentation/clitkExtractMediastinalVesselsFilter.txx +++ b/segmentation/clitkExtractMediastinalVesselsFilter.txx @@ -22,15 +22,21 @@ // clitk #include "clitkCommon.h" #include "clitkExtractMediastinalVesselsFilter.h" -#include "clitkAutoCropFilter.h" #include "clitkSegmentationUtils.h" -#include "clitkMorphoMathFilter.h" +#include "clitkReconstructWithConditionalGrayscaleDilateImageFilter.h" // itk #include -#include #include +template struct index_cmp { + index_cmp(const T varr) : arr(varr) {} + bool operator()(const size_t a, const size_t b) const + { return arr[a] < arr[b]; } + const T arr; +}; + + //-------------------------------------------------------------------- template clitk::ExtractMediastinalVesselsFilter:: @@ -42,8 +48,16 @@ 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"); + SetFinalOpeningRadius(1); } //-------------------------------------------------------------------- @@ -58,145 +72,172 @@ 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"); - debug_eroded.push_back(eroded); - - // CCL - int nb; - 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(); + /*// Erosion kernel + typedef itk::BinaryBallStructuringElement KernelType; + KernelType structuringElement; + structuringElement.SetRadius(radius); + structuringElement.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); + // Erosion -> we break the connectivity between structure + MaskSliceType::SizeType r; + r[0] = r[1] = radius; + MaskSlicePointer eroded = clitk::Opening(slices_mask[i], + r, + GetBackgroundValue(), + GetForegroundValue()); + /* + //typedef itk::BinaryErodeImageFilter ErodeFilterType; + typedef itk::BinaryMorphologicalOpeningImageFilter ErodeFilterType; + typename ErodeFilterType::Pointer eroder = ErodeFilterType::New(); + eroder->SetInput(slices_mask[i]); + eroder->SetBackgroundValue(GetBackgroundValue()); + eroder->SetForegroundValue(GetForegroundValue()); + //eroder->SetBoundaryToForeground(true); // ?? for BinaryErodeImageFilter + eroder->SetKernel(structuringElement); + eroder->Update(); + MaskSlicePointer eroded = eroder->GetOutput(); + */ + + // Keep slice for debug + debug_eroded.push_back(eroded); + // Labelize (CCL) + MaskSlicePointer labeled = + 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); + MaskImagePointType MaxSlicePoint; + if (GetAFDB()->TagExist(m_SoughtVesselSeedName+"Max")) { + GetAFDB()->GetPoint3D(m_SoughtVesselSeedName+"Max", MaxSlicePoint); + } + else { + MaxSlicePoint = SoughtVesselSeedPoint; + MaxSlicePoint[2] += 1000; + } + + // 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, + // MaxSlicePoint, newLabel, bifurcations); + + TrackVesselsFromPoint(r, m_slice_recon, SoughtVesselSeedPoint, + MaxSlicePoint, newLabel); + + // 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). + 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] = GetFinalOpeningRadius(); + 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 +257,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 +270,51 @@ 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,7 +322,9 @@ CropSupInf() { m_Mediastinum = clitk::AutoCrop(m_Mediastinum, GetBackgroundValue()); // Resize input m_Input = clitk::ResizeImageLike(m_Input, m_Mediastinum, GetBackgroundValue()); + */ + // writeImage(m_Input, "crop.mhd"); // End StopCurrentStep(m_Input); } @@ -281,39 +335,42 @@ CropSupInf() { template void clitk::ExtractMediastinalVesselsFilter:: -//SearchBrachioCephalicArtery(int & BCA_first_slice, LabelType & BCA_first_label) { TrackBifurcationFromPoint(MaskImagePointer & recon, std::vector & slices_recon, MaskImagePointType point3D, + MaskImagePointType pointMaxSlice, 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; + + // Max slice + MaskImageIndexType indexMaxSlice; + recon->TransformPhysicalPointToIndex(pointMaxSlice, indexMaxSlice); + uint maxSlice = indexMaxSlice[2]; - 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; @@ -323,74 +380,610 @@ TrackBifurcationFromPoint(MaskImagePointer & recon, typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New(); imageToLabelFilter->SetBackgroundValue(GetBackgroundValue()); imageToLabelFilter->SetInput(s); - statFilter->SetInput(imageToLabelFilter->GetOutput()); + statFilter->SetInput(imageToLabelFilter->GetOutput()); + // statFilter->SetComputeFeretDiameter( true ); + statFilter->ComputePerimeterOn(); // To be able to get proper Roundness value statFilter->Update(); typename LabelMapType::Pointer labelMap = statFilter->GetOutput(); // 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()); + //DD(labels_size.back()); + //DD(labelMap->GetLabelObject(label)->GetRoundness()); + // previousCenter = centroids.back(); } } - DD(centroids.size()); - + + // ------------------------- + // If no centroid were found + 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"); + //DD(previous_slice_label); + //DD(newLabel); + 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()); + } + } + + // Some centroid were found // 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 ?) + // int n = centroids.size(); + // Debug point + std::vector d; + clitk::PointsUtils::Convert2DListTo3DList(centroids, currentSlice, m_Mask, d); + // DDV(d, d.size()); + + /* + // try one or all centroids + std::vector< + for(uint a<=0; a c; + if (a==nb) { for(uint x=0; x 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; + } + } + //DD(indexOfLargest); + //DD(labels_size[indexOfLargest]); + SlicePointType c = centroids[indexOfLargest]; + LabelType l = centroids_label[indexOfLargest]; + //DD(l); + //DD(c); + centroids.clear(); + centroids.push_back(c); + centroids_label.clear(); + centroids_label.push_back(l); } - DD("FOUND"); } + + + /* ==> here all centroid are considered as ok.*/ + + // REMOVE IF CENTROID=1, REPLACE BY >0 // if only one centroids, we change the current image with the current label if (centroids.size() == 1) { + //DD(centroids_label[0]); 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; + } + + if (currentSlice == maxSlice) { + // DD("end max slice"); + found = true; + } // iterate - ++i; + ++currentSlice; + } while (!found); + + // End + StopCurrentStep(); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractMediastinalVesselsFilter:: +TrackVesselsFromPoint(MaskImagePointer & recon, + std::vector & slices, + MaskImagePointType seedPoint, + MaskImagePointType pointMaxSlice, + LabelType newLabel) { + StartNewStep("Track the SoughtVessel from the seed point"); + + // Find first slice index + MaskImageIndexType seedIndex; + recon->TransformPhysicalPointToIndex(seedPoint, seedIndex); + // int numberOfBifurcation = 0; + typedef typename MaskSliceType::PointType SlicePointType; + SlicePointType previousCentroid; + previousCentroid[0] = seedPoint[0]; + previousCentroid[1] = seedPoint[1]; + + // Max slice + MaskImageIndexType indexMaxSlice; + recon->TransformPhysicalPointToIndex(pointMaxSlice, indexMaxSlice); + uint maxSlice = std::min((uint)indexMaxSlice[2], (uint)slices.size()); + // DD(maxSlice); + + // Get current label at the seedPoint of interest + uint currentSlice=seedIndex[2]; + bool found = false; + // LabelType previous_slice_label=recon->GetPixel(seedIndex); + + // Currrent label map variable + typedef itk::ShapeLabelObject< LabelType, 2> LabelObjectType; + typedef itk::LabelMap< LabelObjectType > LabelMapType; + typename LabelMapType::Pointer labelMap; + std::vector shapeObjectsList; + std::vector shapeObjectsSliceList; // keep slice, useful for 'union' + typename LabelObjectType::Pointer previousShapeObject; + + do { + // Debug + //std::cout << std::endl; + //DD(currentSlice); + ImagePointType c; + clitk::PointsUtils::Convert2DTo3D(previousCentroid, m_Mask, currentSlice-1, c); + //DD(c); + + // Consider current reconstructed slice + MaskSlicePointer s = slices[currentSlice]; + MaskSlicePointer previous; + shapeObjectsList.clear(); + shapeObjectsSliceList.clear(); + + // Get shape of all labels in the current slice (it is already labelized) + + // Normal -> same CCL with different label + // PB -> sometimes same label in different CCL ! car growing + //ADD les deux l+s ? mais avec max entre chaque ccl number (bof) + /* + for each label in s -> map avec label in l; if different -> change + */ + MaskSlicePointer ll = clitk::Labelize(s, GetBackgroundValue(), true, 1); + // writeImage(s, "slice-"+toString(currentSlice)+".mhd"); + //writeImage(ll, "slice-label-"+toString(currentSlice)+".mhd"); + typedef itk::ImageRegionIteratorWithIndex IteratorType; + IteratorType its(s, s->GetLargestPossibleRegion()); + IteratorType itl(ll, ll->GetLargestPossibleRegion()); + std::map labelInL; + std::map > labelToChange; + its.GoToBegin(); + itl.GoToBegin(); + int currentLabel = newLabel+10; + while (!its.IsAtEnd()) { + LabelType labs = its.Get(); + if (labs != GetBackgroundValue()) { + LabelType labl = itl.Get(); + if (labelInL.find(labs) == labelInL.end()) { // Not found in map, first time + // DD("first"); + labelInL[labs] = labl; + //DD(labs); + //DD(labl); + } + else { + if (labelInL[labs] != labl) { // I found a labs with a different labl. Store it. + if (labelToChange[labs].find(labl) == labelToChange[labs].end()) { // if not already found + //DD("found"); + //DD(labs); + //DD(labl); + //DD(labelInL[labs]); + //DD(currentLabel); + labelToChange[labs][labl] = currentLabel; + ++currentLabel; + } + } + } + } + ++its; + ++itl; + } + + its.GoToBegin(); + itl.GoToBegin(); + while (!its.IsAtEnd()) { + LabelType labs = its.Get(); + if (labs != GetBackgroundValue()) { // if not BG + LabelType labl = itl.Get(); + if (labelToChange[labs].find(labl) != labelToChange[labs].end()) { // if some labs can change their label + its.Set(labelToChange[labs][labl]); // change with the label for + } + } + ++its; + ++itl; + } // end while + + // writeImage(s, "slice-final"+toString(currentSlice)+".mhd"); + + + labelMap = clitk::ComputeLabelMap(s, GetBackgroundValue(), true); + // DD(labelMap->GetNumberOfLabelObjects()); + + // If this is the first slice, get the object that contains the seed + if (currentSlice == seedIndex[2]) { + // DD("First slice"); + LabelType l = recon->GetPixel(seedIndex); + // DD(l); + shapeObjectsList.push_back(labelMap->GetLabelObject(l)); + shapeObjectsSliceList.push_back(s); + previous = s; + previousCentroid = shapeObjectsList[0]->GetCentroid(); + previousShapeObject = shapeObjectsList[0]; + } + else { + previous = slices[currentSlice-1]; + // Loop on labels to check if centroid is on the previous contour + for(uint c=0; cGetNumberOfLabelObjects(); c++) { + // Get the current label number + int label = labelMap->GetLabels()[c]; + //DD(label); + // Get the centroids + SlicePointType centroid = labelMap->GetLabelObject(label)->GetCentroid(); + // Convert centroid into index in previous slice (same coordinate) + typename MaskSliceType::IndexType centroidIndex; + previous->TransformPhysicalPointToIndex(centroid, centroidIndex); + LabelType labelInPreviousSlice = previous->GetPixel(centroidIndex); + // if this current centroid was in the current label, we keep it + //DD(labelInPreviousSlice); + if (labelInPreviousSlice == newLabel) { + shapeObjectsList.push_back(labelMap->GetLabelObject(label)); + shapeObjectsSliceList.push_back(s); + } + } + } + + + // Potentially the previous centroid could be inside another + // labels, we consider i + typename MaskSliceType::IndexType previousCentroidIndex; + s->TransformPhysicalPointToIndex(previousCentroid, previousCentroidIndex); + LabelType l = s->GetPixel(previousCentroidIndex); + //DD(l); + if (l != 0) { // if is not the background label + int index = -1; + for(uint c=0; cGetLabel() == l) { + index = c; + } + } + if (index == -1) { + //DD("not inside"); + shapeObjectsList.push_back(labelMap->GetLabelObject(l)); + shapeObjectsSliceList.push_back(s); + } + else { + // DD("already inside"); + } + } + + // for(uint c=0; cGetLabel() << " " + // << shapeObjectsList[c]->GetCentroid() << std::endl; + // } + + + // If several candidates, add one more with the union of all candidates + MaskSlicePointer temp; + if (shapeObjectsList.size() > 1) { + //DD("add union"); + // Copy the slice + temp = clitk::Clone(s); + // change label to a single label + LabelType l = newLabel+1; + for(uint c=0; cGetLabel(); + temp = clitk::SetBackground(temp, temp, ll, l, true); + } + // Compute Label object properties + labelMap = clitk::ComputeLabelMap(temp, GetBackgroundValue(), true); + shapeObjectsList.push_back(labelMap->GetLabelObject(l)); + shapeObjectsSliceList.push_back(temp); + } + + /* + for(uint c=0; cGetLabel() << " " + << shapeObjectsList[c]->GetCentroid() << std::endl; + } + */ + + + for(uint c=0; c::Convert2DTo3D(shapeObjectsList[c]->GetCentroid(), m_Mask, currentSlice, cc); + // std::cout << c << " " << shapeObjectsList[c]->GetLabel() << " " + // // << shapeObjectsList[c]->GetCentroid() << " " + // << cc << " " + // << shapeObjectsList[c]->GetPhysicalSize() << " " + // << shapeObjectsList[c]->GetRoundness() << std::endl; + } + + + if (shapeObjectsList.size() == 0) { + found = true; + } + else { + // Heuristic to select the good one. For each candidate, we consider the size + std::vector sizes; + std::vector roundness; + std::vector index_sizes; + std::vector index_roundness; + double previousSize = previousShapeObject->GetPhysicalSize(); + //DD(previousSize); + for(uint c=0; cGetPhysicalSize(); + sizes.push_back(fabs(previousSize-s)/previousSize); + roundness.push_back(fabs(1.0-shapeObjectsList[c]->GetRoundness())); + index_sizes.push_back(c); + index_roundness.push_back(c); + } + //DDV(sizes, sizes.size()); + //DDV(roundness, roundness.size()); + // DDV(index_sizes, index_sizes.size()); + // DDV(index_roundness, index_roundness.size()); + sort(index_sizes.begin(), index_sizes.end(), index_cmp&>(sizes)); + sort(index_roundness.begin(), index_roundness.end(), index_cmp&>(roundness)); + //DDV(index_sizes, index_sizes.size()); + // DDV(index_roundness, index_roundness.size()); + + // TEMPORARY GET THE FIRST + int best = index_sizes[0]; + // if (currentSlice == seedIndex[2]) { // first contour => idiot, first = single contour + // best = index_roundness[0]; // best is more round + // } + LabelType label = shapeObjectsList[best]->GetLabel(); + // DD(label); + s = shapeObjectsSliceList[best]; + s = clitk::SetBackground(s, s, label, newLabel, true); + + // HERE + + // It can happend that several CCL share this same label. To + // prevent this case, we only consider the one that contains + // the centroid. + + // TODO -> PREVIOUS CENTROID ??? + + MaskSlicePointer temp = clitk::Binarize(s, newLabel, newLabel, GetBackgroundValue(), GetForegroundValue()); + temp = clitk::Labelize(temp, GetBackgroundValue(), true, 1); + typename MaskSliceType::IndexType centroids_index; + temp->TransformPhysicalPointToIndex(shapeObjectsList[best]->GetCentroid(), centroids_index); + typename MaskSliceType::PixelType v = temp->GetPixel(centroids_index); + 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; + + + // end + slices[currentSlice] = s; + previousCentroid = shapeObjectsList[best]->GetCentroid(); + previousShapeObject = shapeObjectsList[best]; + } + + ++currentSlice; + + if (currentSlice == maxSlice) { + // DD("end max slice"); + found = true; + } + } while (!found); + + /* + // ------------------------- + // If no centroid were found + if (shapeObjectsList.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"); + DD(previous_slice_label); + DD(newLabel); + typename MaskSliceType::IndexType previousCentroidIndex; + s->TransformPhysicalPointToIndex(previousCentroid, previousCentroidIndex); + DD(previousCentroid); + LabelType labelInSlice = s->GetPixel(previousCentroidIndex); + 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()); + } + } + + // Some centroid were found + // If several centroids, we found a bifurcation + if (centroids.size() > 1) { + // int n = centroids.size(); + // Debug point + std::vector d; + clitk::PointsUtils::Convert2DListTo3DList(centroids, currentSlice, m_Mask, d); + DDV(d, d.size()); + + // new potential bifurcation found + 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; + } + } + DD(indexOfLargest); + DD(labels_size[indexOfLargest]); + SlicePointType c = centroids[indexOfLargest]; + LabelType l = centroids_label[indexOfLargest]; + DD(l); + DD(c); + centroids.clear(); + centroids.push_back(c); + centroids_label.clear(); + centroids_label.push_back(l); + } + } + */ + + /* ==> here all centroid are considered as ok.*/ + + /* + // REMOVE IF CENTROID=1, REPLACE BY >0 + + // if only one centroids, we change the current image with the current label + if (centroids.size() == 1) { + DD(centroids_label[0]); + s = clitk::SetBackground(s, s, centroids_label[0], newLabel, true); + slices[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[currentSlice] = s; + + // I need to recompute the centroid if we have removed some + // connected component. + clitk::ComputeCentroids(s, GetBackgroundValue(), centroids); + previousCentroid = centroids[1]; + } + if (centroids.size() == 0) { + DD("ZERO"); + found = true; + } - //MaskImageType::Pointer rr = clitk::JoinSlices(slices_recon, m_Mask, 2); - //clitk::writeImage(rr, "recon2.mhd"); + if (currentSlice == slices.size()-1) { + DD("end of slices"); + found = true; + } + + if (currentSlice == maxSlice) { + DD("end max slice"); + found = true; + } + + // iterate + ++currentSlice; + } while (!found); + */ // End StopCurrentStep(); @@ -398,5 +991,4 @@ TrackBifurcationFromPoint(MaskImagePointer & recon, //-------------------------------------------------------------------- - #endif //#define CLITKEXTRACTMEDIASTINALVESSELSFILTER_TXX