1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
9 This software is distributed WITHOUT ANY WARRANTY; without even
10 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 PURPOSE. See the copyright notices for more information.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ======================================================================-====*/
19 #ifndef CLITKEXTRACTLYMPHSTATIONSFILTER_TXX
20 #define CLITKEXTRACTLYMPHSTATIONSFILTER_TXX
23 #include "clitkCommon.h"
24 #include "clitkExtractMediastinalVesselsFilter.h"
25 #include "clitkSegmentationUtils.h"
26 #include "clitkReconstructWithConditionalGrayscaleDilateImageFilter.h"
29 #include <itkBinaryThresholdImageFilter.h>
30 #include <itkMinimumMaximumImageCalculator.h>
32 template<class T> struct index_cmp {
33 index_cmp(const T varr) : arr(varr) {}
34 bool operator()(const size_t a, const size_t b) const
35 { return arr[a] < arr[b]; }
40 //--------------------------------------------------------------------
41 template <class TImageType>
42 clitk::ExtractMediastinalVesselsFilter<TImageType>::
43 ExtractMediastinalVesselsFilter():
45 clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
46 itk::ImageToImageFilter<TImageType, MaskImageType>()
48 this->SetNumberOfRequiredInputs(1);
49 SetBackgroundValue(0);
50 SetForegroundValue(1);
51 SetThresholdHigh(140);
54 SetDilatationRadius(9);
55 SetMaxDistancePostToCarina(10);
56 SetMaxDistanceAntToCarina(40);
57 SetMaxDistanceLeftToCarina(35);
58 SetMaxDistanceRightToCarina(35);
59 SetSoughtVesselSeedName("NoSeedNameGiven");
60 SetFinalOpeningRadius(1);
62 //--------------------------------------------------------------------
65 //--------------------------------------------------------------------
66 template <class TImageType>
68 clitk::ExtractMediastinalVesselsFilter<TImageType>::
69 GenerateOutputInformation() {
72 m_Input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
74 // ------------------------------------------------------------------
75 // Crop the initial image superiorly and inferiorly.
76 // TODO : add options for x cm above/below
79 // ------------------------------------------------------------------
80 // Binarize the image. Need two thresholds, one high to select
81 // structures (CCL) that are almost not connected (after erosion),
82 // and one low thresholds to select the real contours. Will be
83 // reconstructed later.
84 StartNewStep("Binarize with high threshold = "+toString(GetThresholdHigh()));
85 typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> BinarizeFilterType;
86 typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
87 binarizeFilter->SetInput(m_Input);
88 binarizeFilter->SetLowerThreshold(GetThresholdHigh());
89 binarizeFilter->SetInsideValue(GetForegroundValue());
90 binarizeFilter->SetOutsideValue(GetBackgroundValue());
91 binarizeFilter->Update();
92 m_Mask = binarizeFilter->GetOutput();
93 StopCurrentStep<MaskImageType>(m_Mask);
95 // ------------------------------------------------------------------
96 StartNewStep("Binarize with low threshold = "+toString(GetThresholdLow()));
97 binarizeFilter = BinarizeFilterType::New();
98 binarizeFilter->SetInput(m_Input);
99 binarizeFilter->SetLowerThreshold(GetThresholdLow());
100 binarizeFilter->SetInsideValue(GetForegroundValue());
101 binarizeFilter->SetOutsideValue(GetBackgroundValue());
102 binarizeFilter->Update();
103 MaskImagePointer m_Mask2 = binarizeFilter->GetOutput();
104 StopCurrentStep<MaskImageType>(m_Mask2);
106 // ------------------------------------------------------------------
108 StartNewStep("Detect objects : erosion, then slice by slice reconstruction");
109 std::vector<MaskSlicePointer> slices_mask;
110 clitk::ExtractSlices<MaskImageType>(m_Mask, 2, slices_mask);
111 std::vector<MaskSlicePointer> slices_mask2;
112 clitk::ExtractSlices<MaskImageType>(m_Mask2, 2, slices_mask2);
113 int radius = GetErosionRadius();
115 // List of working slices (debug only)
116 std::vector<MaskSlicePointer> debug_eroded;
117 // std::vector<MaskSlicePointer> debug_labeled;
119 // ------------------------------------------------------------------
120 // Loop Slice by Slice in order to break connectivity between
121 // CCL. Erode and reconstruct all labels at the same time without
123 for(uint i=0; i<slices_mask.size(); i++) {
126 typedef itk::BinaryBallStructuringElement<MaskSliceType::PixelType,2> KernelType;
127 KernelType structuringElement;
128 structuringElement.SetRadius(radius);
129 structuringElement.CreateStructuringElement();
132 // Erosion -> we break the connectivity between structure
133 MaskSliceType::SizeType r;
134 r[0] = r[1] = radius;
135 MaskSlicePointer eroded = clitk::Opening<MaskSliceType>(slices_mask[i],
137 GetBackgroundValue(),
138 GetForegroundValue());
140 //typedef itk::BinaryErodeImageFilter<MaskSliceType, MaskSliceType, KernelType> ErodeFilterType;
141 typedef itk::BinaryMorphologicalOpeningImageFilter<MaskSliceType, MaskSliceType, KernelType> ErodeFilterType;
142 typename ErodeFilterType::Pointer eroder = ErodeFilterType::New();
143 eroder->SetInput(slices_mask[i]);
144 eroder->SetBackgroundValue(GetBackgroundValue());
145 eroder->SetForegroundValue(GetForegroundValue());
146 //eroder->SetBoundaryToForeground(true); // ?? for BinaryErodeImageFilter
147 eroder->SetKernel(structuringElement);
149 MaskSlicePointer eroded = eroder->GetOutput();
152 // Keep slice for debug
153 debug_eroded.push_back(eroded);
156 MaskSlicePointer labeled =
157 clitk::Labelize<MaskSliceType>(eroded, GetBackgroundValue(), true, 1); // Fully connected !
158 // debug_labeled.push_back(labeled);
160 // Make Reconstruction filter : dilation all labels at the same
161 // time, prevent to merge them.
162 typedef clitk::ReconstructWithConditionalGrayscaleDilateImageFilter<MaskSliceType> ReconstructFilterType;
163 typename ReconstructFilterType::Pointer reconstructor = ReconstructFilterType::New();
164 reconstructor->SetInput(labeled);
165 reconstructor->SetIterationNumber(radius+GetDilatationRadius());
166 reconstructor->Update();
167 MaskSlicePointer s = reconstructor->GetOutput();
169 // Remove Initial BG of the second tresholded image
170 s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, slices_mask2[i],
171 GetBackgroundValue(), GetBackgroundValue(), true);
172 m_slice_recon.push_back(s);
176 // Build 3D images from the slice by slice processing
177 MaskImageType::Pointer eroded = clitk::JoinSlices<MaskImageType>(debug_eroded, m_Mask, 2);
178 writeImage<MaskImageType>(eroded, "erode.mhd");
179 //MaskImageType::Pointer l = clitk::JoinSlices<MaskImageType>(debug_labeled, m_Mask, 2);
180 MaskImageType::Pointer r = clitk::JoinSlices<MaskImageType>(m_slice_recon, m_Mask, 2);
181 writeImage<MaskImageType>(r, "recon1.mhd");
183 // ------------------------------------------------------------------
184 // Track the SoughtVessel from the given first point
185 // superiorly. This is done by TrackBifurcationFromPoint
186 MaskImagePointType SoughtVesselSeedPoint;
187 GetAFDB()->GetPoint3D(m_SoughtVesselSeedName, SoughtVesselSeedPoint);
188 MaskImagePointType MaxSlicePoint;
189 if (GetAFDB()->TagExist(m_SoughtVesselSeedName+"Max")) {
190 GetAFDB()->GetPoint3D(m_SoughtVesselSeedName+"Max", MaxSlicePoint);
193 MaxSlicePoint = SoughtVesselSeedPoint;
194 MaxSlicePoint[2] += 1000;
197 // Find the label with the maximum value to set the result
198 typedef itk::MinimumMaximumImageCalculator<MaskImageType> MinMaxFilterType;
199 MinMaxFilterType::Pointer ff = MinMaxFilterType::New();
201 ff->ComputeMaximum();
202 LabelType newLabel = ff->GetMaximum()+1;
204 // the following bifurcations point will the centroids of the
205 // components obtain when (hopefully!) the SoughtVessel
206 // split into CommonArtery and SubclavianArtery.
207 std::vector<MaskImagePointType> bifurcations;
208 // TrackBifurcationFromPoint(r, m_slice_recon, SoughtVesselSeedPoint,
209 // MaxSlicePoint, newLabel, bifurcations);
211 TrackVesselsFromPoint(r, m_slice_recon, SoughtVesselSeedPoint,
212 MaxSlicePoint, newLabel);
214 // Build the final 3D image from the previous slice by slice processing
215 m_SoughtVessel = clitk::JoinSlices<MaskImageType>(m_slice_recon, m_Mask, 2);
216 writeImage<MaskImageType>(m_SoughtVessel, "recon2.mhd");
218 // Set binary image, (remove other labels).
220 clitk::Binarize<MaskImageType>(m_SoughtVessel, newLabel, newLabel,
221 GetBackgroundValue(), GetForegroundValue());
223 // writeImage<MaskImageType>(m_SoughtVessel, "afterbinarize.mhd");
224 m_SoughtVessel = clitk::AutoCrop<MaskImageType>(m_SoughtVessel, GetBackgroundValue());
226 // writeImage<MaskImageType>(m_SoughtVessel, "afterautocrop.mhd");
228 // Clean the image : Opening (not in Z direction)
229 typename MaskImageType::SizeType rad;
230 rad[0] = rad[1] = GetFinalOpeningRadius();
232 m_SoughtVessel = clitk::Opening<MaskImageType>(m_SoughtVessel, rad,
233 GetBackgroundValue(), GetForegroundValue());
235 // writeImage<MaskImageType>(m_SoughtVessel, "afteropen.mhd");
237 // Clean the image : keep main CCL slice by slice
238 m_SoughtVessel = clitk::SliceBySliceKeepMainCCL<MaskImageType>(m_SoughtVessel,
239 GetBackgroundValue(),
240 GetForegroundValue());
242 //--------------------------------------------------------------------
245 //--------------------------------------------------------------------
246 template <class TImageType>
248 clitk::ExtractMediastinalVesselsFilter<TImageType>::
249 GenerateInputRequestedRegion() {
250 //DD("GenerateInputRequestedRegion (nothing?)");
252 //--------------------------------------------------------------------
255 //--------------------------------------------------------------------
256 template <class TImageType>
258 clitk::ExtractMediastinalVesselsFilter<TImageType>::
260 // Save in the AFDB (not write on the disk here)
261 GetAFDB()->SetImageFilename(GetSoughtVesselName(), GetOutputFilename());
263 // Final Step -> graft output
264 this->GraftNthOutput(0, m_SoughtVessel);
266 //--------------------------------------------------------------------
269 //--------------------------------------------------------------------
270 template <class TImageType>
272 clitk::ExtractMediastinalVesselsFilter<TImageType>::
274 StartNewStep("Crop the input image: SI,AP limits with carina and crop with mediastinum");
276 Need : Trachea, Carina (roi not point),
278 // Get Trachea and Carina
279 MaskImagePointer Trachea = GetAFDB()->template GetImage <MaskImageType>("Trachea");
281 // Compute Carina position
283 MaskImagePointer Carina = GetAFDB()->template GetImage<MaskImageType>("Carina");
284 std::vector<MaskImagePointType> centroids;
285 clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
286 m_CarinaZ = centroids[1][2];
287 // add one slice to include carina
288 m_CarinaZ += Carina->GetSpacing()[2];
289 // We dont need Carina structure from now
291 GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]);
293 // Crop Inf, remove below Carina
295 clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 2, m_CarinaZ, false, GetBackgroundValue());
298 double m_CarinaY = centroids[1][1];
299 m_Input = clitk::CropImageRemoveGreaterThan<ImageType>(m_Input, 1,
300 m_CarinaY+GetMaxDistancePostToCarina(),
301 false, GetBackgroundValue());
303 m_Input = clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 1,
304 m_CarinaY-GetMaxDistanceAntToCarina(),
305 false, GetBackgroundValue());
307 double m_CarinaX = centroids[1][0];
308 m_Input = clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 0,
309 m_CarinaX-GetMaxDistanceRightToCarina(),
310 false, GetBackgroundValue());
312 m_Input = clitk::CropImageRemoveGreaterThan<ImageType>(m_Input, 0,
313 m_CarinaX+GetMaxDistanceLeftToCarina(),
314 false, GetBackgroundValue());
317 // AutoCrop with Mediastinum, generaly only allow to remove few slices (superiorly)
318 m_Mediastinum = GetAFDB()->template GetImage<MaskImageType>("Mediastinum");
319 // Resize like input (sup to carina)
320 m_Mediastinum = clitk::ResizeImageLike<MaskImageType>(m_Mediastinum, m_Input, GetBackgroundValue());
322 m_Mediastinum = clitk::AutoCrop<MaskImageType>(m_Mediastinum, GetBackgroundValue());
324 m_Input = clitk::ResizeImageLike<ImageType>(m_Input, m_Mediastinum, GetBackgroundValue());
327 // writeImage<ImageType>(m_Input, "crop.mhd");
329 StopCurrentStep<ImageType>(m_Input);
331 //--------------------------------------------------------------------
334 //--------------------------------------------------------------------
335 template <class TImageType>
337 clitk::ExtractMediastinalVesselsFilter<TImageType>::
338 TrackBifurcationFromPoint(MaskImagePointer & recon,
339 std::vector<MaskSlicePointer> & slices_recon,
340 MaskImagePointType point3D,
341 MaskImagePointType pointMaxSlice,
343 std::vector<MaskImagePointType> & bifurcations) {
344 StartNewStep("Track the SoughtVessel from the seed point");
346 // Find first slice index
347 MaskImageIndexType index;
348 recon->TransformPhysicalPointToIndex(point3D, index);
349 int numberOfBifurcation = 0;
350 typedef typename MaskSliceType::PointType SlicePointType;
351 SlicePointType previousCenter;
354 MaskImageIndexType indexMaxSlice;
355 recon->TransformPhysicalPointToIndex(pointMaxSlice, indexMaxSlice);
356 uint maxSlice = indexMaxSlice[2];
358 // Get current label at the point3D of interest
359 uint currentSlice=index[2];
361 LabelType previous_slice_label=recon->GetPixel(index);
362 // DD(slices_recon.size());
365 // Consider current reconstructed slice
366 MaskSlicePointer s = slices_recon[currentSlice];
367 MaskSlicePointer previous;
368 if (currentSlice == index[2]) previous = s;
370 previous = slices_recon[currentSlice-1];
373 // Get centroids of the labels in the current slice
374 static const unsigned int Dim = MaskSliceType::ImageDimension;
375 typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
376 typedef itk::LabelMap< LabelObjectType > LabelMapType;
377 typedef itk::LabelImageToLabelMapFilter<MaskSliceType, LabelMapType> ImageToMapFilterType;
378 typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
379 typedef itk::ShapeLabelMapFilter<LabelMapType, MaskSliceType> ShapeFilterType;
380 typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
381 imageToLabelFilter->SetBackgroundValue(GetBackgroundValue());
382 imageToLabelFilter->SetInput(s);
383 statFilter->SetInput(imageToLabelFilter->GetOutput());
384 // statFilter->SetComputeFeretDiameter( true );
385 statFilter->ComputePerimeterOn(); // To be able to get proper Roundness value
386 statFilter->Update();
387 typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
389 // Look what centroid inside the previous largest one
390 std::vector<SlicePointType> centroids;
391 std::vector<LabelType> centroids_label;
392 std::vector<double> labels_size;
393 for(uint c=0; c<labelMap->GetNumberOfLabelObjects(); c++) {
394 int label = labelMap->GetLabels()[c];
396 SlicePointType center = labelMap->GetLabelObject(label)->GetCentroid();
398 // Get label into previous slice
399 typename MaskSliceType::IndexType centerIndex;
400 previous->TransformPhysicalPointToIndex(center, centerIndex);
401 LabelType labelInPreviousSlice = previous->GetPixel(centerIndex);
402 // if this current centroid was in the current label, add it to bifurcations
403 if (labelInPreviousSlice == previous_slice_label) {
404 centroids.push_back(center);
405 centroids_label.push_back(label);
406 labels_size.push_back(labelMap->GetLabelObject(label)->GetPhysicalSize());
407 //DD(labels_size.back());
408 //DD(labelMap->GetLabelObject(label)->GetRoundness());
409 // previousCenter = centroids.back();
413 // -------------------------
414 // If no centroid were found
415 if (centroids.size() == 0) {
416 // Last attempt to find -> check if previous centroid is inside a CCL
417 // if in s -> get value, getcentroid add.
419 //DD("Last change to find");
420 //DD(previous_slice_label);
422 typename MaskSliceType::IndexType previousCenterIndex;
423 s->TransformPhysicalPointToIndex(previousCenter, previousCenterIndex);
424 //DD(previousCenter);
425 LabelType labelInSlice = s->GetPixel(previousCenterIndex);
427 if (labelInSlice != GetBackgroundValue()) {
428 centroids.push_back(labelMap->GetLabelObject(labelInSlice)->GetCentroid());
429 centroids_label.push_back(labelInSlice);
430 labels_size.push_back(labelMap->GetLabelObject(labelInSlice)->GetPhysicalSize());
434 // Some centroid were found
435 // If several centroids, we found a bifurcation
436 if (centroids.size() > 1) {
437 // int n = centroids.size();
439 std::vector<ImagePointType> d;
440 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(centroids, currentSlice, m_Mask, d);
444 // try one or all centroids
446 for(uint a<=0; a<nb++; a++) {
448 // Create the list of candidates
450 if (a==nb) { for(uint x=0; x<nb; x++) c.push_back(x); }
456 for(uint x=0; x<c.size(); c++) { size += labels_size[c[x]]; }
461 // new potential bifurcation found
462 numberOfBifurcation++;
463 // If the number of bifurcation is greater than the required one, we stop
464 if (numberOfBifurcation > GetMaxNumberOfFoundBifurcation()) {
466 //DD("max bif reach");
467 for(uint c=0; c<centroids.size(); c++) {
469 clitk::PointsUtils<MaskImageType>::Convert2DTo3D(centroids[c], m_Mask, currentSlice, bif);
470 bifurcations.push_back(bif);
473 // Else we continue along the main (largest) connected component
475 int indexOfLargest = 0;
476 for(uint b=0; b<centroids.size(); b++) {
477 if (labels_size[b] > labels_size[indexOfLargest]) {
481 //DD(indexOfLargest);
482 //DD(labels_size[indexOfLargest]);
483 SlicePointType c = centroids[indexOfLargest];
484 LabelType l = centroids_label[indexOfLargest];
488 centroids.push_back(c);
489 centroids_label.clear();
490 centroids_label.push_back(l);
495 /* ==> here all centroid are considered as ok.*/
497 // REMOVE IF CENTROID=1, REPLACE BY >0
499 // if only one centroids, we change the current image with the current label
500 if (centroids.size() == 1) {
501 //DD(centroids_label[0]);
502 s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, s, centroids_label[0], newLabel, true);
503 slices_recon[currentSlice] = s;
504 previous_slice_label = newLabel;
505 // It can happend that several CCL share this same label. To
506 // prevent this case, we only consider the one that contains
508 MaskSlicePointer temp = clitk::Binarize<MaskSliceType>(s, newLabel, newLabel, GetBackgroundValue(), GetForegroundValue());
509 // writeImage<MaskSliceType>(temp, "bin-"+toString(currentSlice)+".mhd");
510 temp = clitk::Labelize<MaskSliceType>(temp, GetBackgroundValue(), true, 1);
511 //writeImage<MaskSliceType>(temp, "label-"+toString(currentSlice)+".mhd");
512 typename MaskSliceType::IndexType centroids_index;
513 temp->TransformPhysicalPointToIndex(centroids[0], centroids_index);
514 typename MaskSliceType::PixelType v = temp->GetPixel(centroids_index);
516 // It can happend that the centroid is inside the BG, so we keep
517 // the largest CCL (the first);
518 if (v == GetBackgroundValue()) {
522 v = 1; // largest one
526 temp = clitk::Binarize<MaskSliceType>(temp, v, v, GetBackgroundValue(), newLabel);
527 //writeImage<MaskSliceType>(temp, "relabel-"+toString(currentSlice)+".mhd");
529 slices_recon[currentSlice] = s;
531 // I need to recompute the centroid if we have removed some
532 // connected component.
533 clitk::ComputeCentroids<MaskSliceType>(s, GetBackgroundValue(), centroids);
534 previousCenter = centroids[1];
537 if (centroids.size() == 0) {
542 if (currentSlice == slices_recon.size()-1) {
543 // DD("end of slices");
547 if (currentSlice == maxSlice) {
548 // DD("end max slice");
559 //--------------------------------------------------------------------
562 //--------------------------------------------------------------------
563 template <class TImageType>
565 clitk::ExtractMediastinalVesselsFilter<TImageType>::
566 TrackVesselsFromPoint(MaskImagePointer & recon,
567 std::vector<MaskSlicePointer> & slices,
568 MaskImagePointType seedPoint,
569 MaskImagePointType pointMaxSlice,
570 LabelType newLabel) {
571 StartNewStep("Track the SoughtVessel from the seed point");
573 // Find first slice index
574 MaskImageIndexType seedIndex;
575 recon->TransformPhysicalPointToIndex(seedPoint, seedIndex);
576 // int numberOfBifurcation = 0;
577 typedef typename MaskSliceType::PointType SlicePointType;
578 SlicePointType previousCentroid;
579 previousCentroid[0] = seedPoint[0];
580 previousCentroid[1] = seedPoint[1];
583 MaskImageIndexType indexMaxSlice;
584 recon->TransformPhysicalPointToIndex(pointMaxSlice, indexMaxSlice);
585 uint maxSlice = std::min((uint)indexMaxSlice[2], (uint)slices.size());
588 // Get current label at the seedPoint of interest
589 uint currentSlice=seedIndex[2];
591 // LabelType previous_slice_label=recon->GetPixel(seedIndex);
593 // Currrent label map variable
594 typedef itk::ShapeLabelObject< LabelType, 2> LabelObjectType;
595 typedef itk::LabelMap< LabelObjectType > LabelMapType;
596 typename LabelMapType::Pointer labelMap;
597 std::vector<typename LabelObjectType::Pointer> shapeObjectsList;
598 std::vector<MaskSlicePointer> shapeObjectsSliceList; // keep slice, useful for 'union'
599 typename LabelObjectType::Pointer previousShapeObject;
603 std::cout << std::endl;
606 clitk::PointsUtils<MaskImageType>::Convert2DTo3D(previousCentroid, m_Mask, currentSlice-1, c);
609 // Consider current reconstructed slice
610 MaskSlicePointer s = slices[currentSlice];
611 MaskSlicePointer previous;
612 shapeObjectsList.clear();
613 shapeObjectsSliceList.clear();
615 // Get shape of all labels in the current slice (it is already labelized)
617 // Normal -> same CCL with different label
618 // PB -> sometimes same label in different CCL ! car growing
619 //ADD les deux l+s ? mais avec max entre chaque ccl number (bof)
621 for each label in s -> map avec label in l; if different -> change
623 MaskSlicePointer ll = clitk::Labelize<MaskSliceType>(s, GetBackgroundValue(), true, 1);
624 // writeImage<MaskSliceType>(s, "slice-"+toString(currentSlice)+".mhd");
625 //writeImage<MaskSliceType>(ll, "slice-label-"+toString(currentSlice)+".mhd");
626 typedef itk::ImageRegionIteratorWithIndex<MaskSliceType> IteratorType;
627 IteratorType its(s, s->GetLargestPossibleRegion());
628 IteratorType itl(ll, ll->GetLargestPossibleRegion());
629 std::map<LabelType, LabelType> labelInL;
630 std::map<LabelType, std::map<LabelType, LabelType> > labelToChange;
633 int currentLabel = newLabel+10;
634 while (!its.IsAtEnd()) {
635 LabelType labs = its.Get();
636 if (labs != GetBackgroundValue()) {
637 LabelType labl = itl.Get();
638 if (labelInL.find(labs) == labelInL.end()) { // Not found in map, first time
640 labelInL[labs] = labl;
645 if (labelInL[labs] != labl) { // I found a labs with a different labl. Store it.
646 if (labelToChange[labs].find(labl) == labelToChange[labs].end()) { // if not already found
650 //DD(labelInL[labs]);
652 labelToChange[labs][labl] = currentLabel;
664 while (!its.IsAtEnd()) {
665 LabelType labs = its.Get();
666 if (labs != GetBackgroundValue()) { // if not BG
667 LabelType labl = itl.Get();
668 if (labelToChange[labs].find(labl) != labelToChange[labs].end()) { // if some labs can change their label
669 its.Set(labelToChange[labs][labl]); // change with the label for <labs-labl>
676 // writeImage<MaskSliceType>(s, "slice-final"+toString(currentSlice)+".mhd");
679 labelMap = clitk::ComputeLabelMap<MaskSliceType, LabelType>(s, GetBackgroundValue(), true);
680 // DD(labelMap->GetNumberOfLabelObjects());
682 // If this is the first slice, get the object that contains the seed
683 if (currentSlice == seedIndex[2]) {
684 // DD("First slice");
685 LabelType l = recon->GetPixel(seedIndex);
687 shapeObjectsList.push_back(labelMap->GetLabelObject(l));
688 shapeObjectsSliceList.push_back(s);
690 previousCentroid = shapeObjectsList[0]->GetCentroid();
691 previousShapeObject = shapeObjectsList[0];
694 previous = slices[currentSlice-1];
695 // Loop on labels to check if centroid is on the previous contour
696 for(uint c=0; c<labelMap->GetNumberOfLabelObjects(); c++) {
697 // Get the current label number
698 int label = labelMap->GetLabels()[c];
701 SlicePointType centroid = labelMap->GetLabelObject(label)->GetCentroid();
702 // Convert centroid into index in previous slice (same coordinate)
703 typename MaskSliceType::IndexType centroidIndex;
704 previous->TransformPhysicalPointToIndex(centroid, centroidIndex);
705 LabelType labelInPreviousSlice = previous->GetPixel(centroidIndex);
706 // if this current centroid was in the current label, we keep it
707 //DD(labelInPreviousSlice);
708 if (labelInPreviousSlice == newLabel) {
709 shapeObjectsList.push_back(labelMap->GetLabelObject(label));
710 shapeObjectsSliceList.push_back(s);
716 // Potentially the previous centroid could be inside another
717 // labels, we consider i
718 typename MaskSliceType::IndexType previousCentroidIndex;
719 s->TransformPhysicalPointToIndex(previousCentroid, previousCentroidIndex);
720 LabelType l = s->GetPixel(previousCentroidIndex);
722 if (l != 0) { // if is not the background label
724 for(uint c=0; c<shapeObjectsList.size(); c++) {
725 if (shapeObjectsList[c]->GetLabel() == l) {
731 shapeObjectsList.push_back(labelMap->GetLabelObject(l));
732 shapeObjectsSliceList.push_back(s);
735 // DD("already inside");
739 // for(uint c=0; c<shapeObjectsList.size(); c++) {
740 // std::cout << c << " " << shapeObjectsList[c]->GetLabel() << " "
741 // << shapeObjectsList[c]->GetCentroid() << std::endl;
745 // If several candidates, add one more with the union of all candidates
746 MaskSlicePointer temp;
747 if (shapeObjectsList.size() > 1) {
750 temp = clitk::Clone<MaskSliceType>(s);
751 // change label to a single label
752 LabelType l = newLabel+1;
753 for(uint c=0; c<shapeObjectsList.size(); c++) {
754 LabelType ll = shapeObjectsList[c]->GetLabel();
755 temp = clitk::SetBackground<MaskSliceType, MaskSliceType>(temp, temp, ll, l, true);
757 // Compute Label object properties
758 labelMap = clitk::ComputeLabelMap<MaskSliceType, LabelType>(temp, GetBackgroundValue(), true);
759 shapeObjectsList.push_back(labelMap->GetLabelObject(l));
760 shapeObjectsSliceList.push_back(temp);
764 for(uint c=0; c<shapeObjectsList.size(); c++) {
765 std::cout << c << " " << shapeObjectsList[c]->GetLabel() << " "
766 << shapeObjectsList[c]->GetCentroid() << std::endl;
771 for(uint c=0; c<shapeObjectsList.size(); c++) {
773 clitk::PointsUtils<MaskImageType>::Convert2DTo3D(shapeObjectsList[c]->GetCentroid(), m_Mask, currentSlice, cc);
774 // std::cout << c << " " << shapeObjectsList[c]->GetLabel() << " "
775 // // << shapeObjectsList[c]->GetCentroid() << " "
777 // << shapeObjectsList[c]->GetPhysicalSize() << " "
778 // << shapeObjectsList[c]->GetRoundness() << std::endl;
782 if (shapeObjectsList.size() == 0) {
786 // Heuristic to select the good one. For each candidate, we consider the size
787 std::vector<double> sizes;
788 std::vector<double> roundness;
789 std::vector<size_t> index_sizes;
790 std::vector<size_t> index_roundness;
791 double previousSize = previousShapeObject->GetPhysicalSize();
793 for(uint c=0; c<shapeObjectsList.size(); c++) {
794 double s = shapeObjectsList[c]->GetPhysicalSize();
795 sizes.push_back(fabs(previousSize-s)/previousSize);
796 roundness.push_back(fabs(1.0-shapeObjectsList[c]->GetRoundness()));
797 index_sizes.push_back(c);
798 index_roundness.push_back(c);
800 //DDV(sizes, sizes.size());
801 //DDV(roundness, roundness.size());
802 // DDV(index_sizes, index_sizes.size());
803 // DDV(index_roundness, index_roundness.size());
804 sort(index_sizes.begin(), index_sizes.end(), index_cmp<std::vector<double>&>(sizes));
805 sort(index_roundness.begin(), index_roundness.end(), index_cmp<std::vector<double>&>(roundness));
806 //DDV(index_sizes, index_sizes.size());
807 // DDV(index_roundness, index_roundness.size());
809 // TEMPORARY GET THE FIRST
810 int best = index_sizes[0];
811 // if (currentSlice == seedIndex[2]) { // first contour => idiot, first = single contour
812 // best = index_roundness[0]; // best is more round
814 LabelType label = shapeObjectsList[best]->GetLabel();
816 s = shapeObjectsSliceList[best];
817 s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, s, label, newLabel, true);
821 // It can happend that several CCL share this same label. To
822 // prevent this case, we only consider the one that contains
825 // TODO -> PREVIOUS CENTROID ???
827 MaskSlicePointer temp = clitk::Binarize<MaskSliceType>(s, newLabel, newLabel, GetBackgroundValue(), GetForegroundValue());
828 temp = clitk::Labelize<MaskSliceType>(temp, GetBackgroundValue(), true, 1);
829 typename MaskSliceType::IndexType centroids_index;
830 temp->TransformPhysicalPointToIndex(shapeObjectsList[best]->GetCentroid(), centroids_index);
831 typename MaskSliceType::PixelType v = temp->GetPixel(centroids_index);
832 if (v == GetBackgroundValue()) {
836 v = 1; // largest one
840 temp = clitk::Binarize<MaskSliceType>(temp, v, v, GetBackgroundValue(), newLabel);
841 //writeImage<MaskSliceType>(temp, "relabel-"+toString(currentSlice)+".mhd");
846 slices[currentSlice] = s;
847 previousCentroid = shapeObjectsList[best]->GetCentroid();
848 previousShapeObject = shapeObjectsList[best];
853 if (currentSlice == maxSlice) {
854 // DD("end max slice");
861 // -------------------------
862 // If no centroid were found
863 if (shapeObjectsList.size() == 0) {
864 // Last attempt to find -> check if previous centroid is inside a CCL
865 // if in s -> get value, getcentroid add.
867 DD("Last change to find");
868 DD(previous_slice_label);
870 typename MaskSliceType::IndexType previousCentroidIndex;
871 s->TransformPhysicalPointToIndex(previousCentroid, previousCentroidIndex);
872 DD(previousCentroid);
873 LabelType labelInSlice = s->GetPixel(previousCentroidIndex);
875 if (labelInSlice != GetBackgroundValue()) {
876 centroids.push_back(labelMap->GetLabelObject(labelInSlice)->GetCentroid());
877 centroids_label.push_back(labelInSlice);
878 labels_size.push_back(labelMap->GetLabelObject(labelInSlice)->GetPhysicalSize());
882 // Some centroid were found
883 // If several centroids, we found a bifurcation
884 if (centroids.size() > 1) {
885 // int n = centroids.size();
887 std::vector<ImagePointType> d;
888 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(centroids, currentSlice, m_Mask, d);
891 // new potential bifurcation found
892 numberOfBifurcation++;
893 // If the number of bifurcation is greater than the required one, we stop
894 if (numberOfBifurcation > GetMaxNumberOfFoundBifurcation()) {
897 for(uint c=0; c<centroids.size(); c++) {
899 clitk::PointsUtils<MaskImageType>::Convert2DTo3D(centroids[c], m_Mask, currentSlice, bif);
900 bifurcations.push_back(bif);
903 // Else we continue along the main (largest) connected component
905 int indexOfLargest = 0;
906 for(uint b=0; b<centroids.size(); b++) {
907 if (labels_size[b] > labels_size[indexOfLargest]) {
912 DD(labels_size[indexOfLargest]);
913 SlicePointType c = centroids[indexOfLargest];
914 LabelType l = centroids_label[indexOfLargest];
918 centroids.push_back(c);
919 centroids_label.clear();
920 centroids_label.push_back(l);
925 /* ==> here all centroid are considered as ok.*/
928 // REMOVE IF CENTROID=1, REPLACE BY >0
930 // if only one centroids, we change the current image with the current label
931 if (centroids.size() == 1) {
932 DD(centroids_label[0]);
933 s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, s, centroids_label[0], newLabel, true);
934 slices[currentSlice] = s;
935 previous_slice_label = newLabel;
936 // It can happend that several CCL share this same label. To
937 // prevent this case, we only consider the one that contains
939 MaskSlicePointer temp = clitk::Binarize<MaskSliceType>(s, newLabel, newLabel, GetBackgroundValue(), GetForegroundValue());
940 // writeImage<MaskSliceType>(temp, "bin-"+toString(currentSlice)+".mhd");
941 temp = clitk::Labelize<MaskSliceType>(temp, GetBackgroundValue(), true, 1);
942 //writeImage<MaskSliceType>(temp, "label-"+toString(currentSlice)+".mhd");
943 typename MaskSliceType::IndexType centroids_index;
944 temp->TransformPhysicalPointToIndex(centroids[0], centroids_index);
945 typename MaskSliceType::PixelType v = temp->GetPixel(centroids_index);
947 // It can happend that the centroid is inside the BG, so we keep
948 // the largest CCL (the first);
949 if (v == GetBackgroundValue()) {
953 v = 1; // largest one
957 temp = clitk::Binarize<MaskSliceType>(temp, v, v, GetBackgroundValue(), newLabel);
958 //writeImage<MaskSliceType>(temp, "relabel-"+toString(currentSlice)+".mhd");
960 slices[currentSlice] = s;
962 // I need to recompute the centroid if we have removed some
963 // connected component.
964 clitk::ComputeCentroids<MaskSliceType>(s, GetBackgroundValue(), centroids);
965 previousCentroid = centroids[1];
968 if (centroids.size() == 0) {
973 if (currentSlice == slices.size()-1) {
978 if (currentSlice == maxSlice) {
991 //--------------------------------------------------------------------
994 #endif //#define CLITKEXTRACTMEDIASTINALVESSELSFILTER_TXX