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);
61 VerboseTrackingFlagOff();
63 //--------------------------------------------------------------------
66 //--------------------------------------------------------------------
67 template <class TImageType>
69 clitk::ExtractMediastinalVesselsFilter<TImageType>::
70 GenerateOutputInformation() {
73 m_Input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
75 // ------------------------------------------------------------------
76 // Crop the initial image superiorly and inferiorly.
77 // TODO : add options for x cm above/below
80 // ------------------------------------------------------------------
81 // Binarize the image. Need two thresholds, one high to select
82 // structures (CCL) that are almost not connected (after erosion),
83 // and one low thresholds to select the real contours. Will be
84 // reconstructed later.
85 StartNewStep("Binarize with high threshold = "+toString(GetThresholdHigh()));
86 typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> BinarizeFilterType;
87 typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
88 binarizeFilter->SetInput(m_Input);
89 binarizeFilter->SetLowerThreshold(GetThresholdHigh());
90 binarizeFilter->SetInsideValue(GetForegroundValue());
91 binarizeFilter->SetOutsideValue(GetBackgroundValue());
92 binarizeFilter->Update();
93 m_Mask = binarizeFilter->GetOutput();
94 StopCurrentStep<MaskImageType>(m_Mask);
96 // ------------------------------------------------------------------
97 StartNewStep("Binarize with low threshold = "+toString(GetThresholdLow()));
98 binarizeFilter = BinarizeFilterType::New();
99 binarizeFilter->SetInput(m_Input);
100 binarizeFilter->SetLowerThreshold(GetThresholdLow());
101 binarizeFilter->SetInsideValue(GetForegroundValue());
102 binarizeFilter->SetOutsideValue(GetBackgroundValue());
103 binarizeFilter->Update();
104 MaskImagePointer m_Mask2 = binarizeFilter->GetOutput();
105 StopCurrentStep<MaskImageType>(m_Mask2);
107 // ------------------------------------------------------------------
109 StartNewStep("Detect objects : erosion, then slice by slice reconstruction");
110 std::vector<MaskSlicePointer> slices_mask;
111 clitk::ExtractSlices<MaskImageType>(m_Mask, 2, slices_mask);
112 std::vector<MaskSlicePointer> slices_mask2;
113 clitk::ExtractSlices<MaskImageType>(m_Mask2, 2, slices_mask2);
114 int radius = GetErosionRadius();
116 // List of working slices (debug only)
117 std::vector<MaskSlicePointer> debug_eroded;
118 // std::vector<MaskSlicePointer> debug_labeled;
120 // ------------------------------------------------------------------
121 // Loop Slice by Slice in order to break connectivity between
122 // CCL. Erode and reconstruct all labels at the same time without
124 for(uint i=0; i<slices_mask.size(); i++) {
127 typedef itk::BinaryBallStructuringElement<MaskSliceType::PixelType,2> KernelType;
128 KernelType structuringElement;
129 structuringElement.SetRadius(radius);
130 structuringElement.CreateStructuringElement();
133 // Erosion -> we break the connectivity between structure
134 MaskSliceType::SizeType r;
135 r[0] = r[1] = radius;
136 MaskSlicePointer eroded = clitk::Opening<MaskSliceType>(slices_mask[i],
138 GetBackgroundValue(),
139 GetForegroundValue());
141 //typedef itk::BinaryErodeImageFilter<MaskSliceType, MaskSliceType, KernelType> ErodeFilterType;
142 typedef itk::BinaryMorphologicalOpeningImageFilter<MaskSliceType, MaskSliceType, KernelType> ErodeFilterType;
143 typename ErodeFilterType::Pointer eroder = ErodeFilterType::New();
144 eroder->SetInput(slices_mask[i]);
145 eroder->SetBackgroundValue(GetBackgroundValue());
146 eroder->SetForegroundValue(GetForegroundValue());
147 //eroder->SetBoundaryToForeground(true); // ?? for BinaryErodeImageFilter
148 eroder->SetKernel(structuringElement);
150 MaskSlicePointer eroded = eroder->GetOutput();
153 // Keep slice for debug
154 debug_eroded.push_back(eroded);
157 MaskSlicePointer labeled =
158 clitk::Labelize<MaskSliceType>(eroded, GetBackgroundValue(), true, 1); // Fully connected !
159 // debug_labeled.push_back(labeled);
161 // Make Reconstruction filter : dilation all labels at the same
162 // time, prevent to merge them.
163 typedef clitk::ReconstructWithConditionalGrayscaleDilateImageFilter<MaskSliceType> ReconstructFilterType;
164 typename ReconstructFilterType::Pointer reconstructor = ReconstructFilterType::New();
165 reconstructor->SetInput(labeled);
166 reconstructor->SetIterationNumber(radius+GetDilatationRadius());
167 reconstructor->Update();
168 MaskSlicePointer s = reconstructor->GetOutput();
170 // Remove Initial BG of the second tresholded image
171 s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, slices_mask2[i],
172 GetBackgroundValue(), GetBackgroundValue(), true);
173 m_slice_recon.push_back(s);
177 // Build 3D images from the slice by slice processing
178 MaskImageType::Pointer eroded = clitk::JoinSlices<MaskImageType>(debug_eroded, m_Mask, 2);
179 writeImage<MaskImageType>(eroded, "erode.mhd");
180 //MaskImageType::Pointer l = clitk::JoinSlices<MaskImageType>(debug_labeled, m_Mask, 2);
181 MaskImageType::Pointer r = clitk::JoinSlices<MaskImageType>(m_slice_recon, m_Mask, 2);
182 writeImage<MaskImageType>(r, "recon1.mhd");
184 // ------------------------------------------------------------------
185 // Track the SoughtVessel from the given first point
186 // superiorly. This is done by TrackBifurcationFromPoint
187 MaskImagePointType SoughtVesselSeedPoint;
188 GetAFDB()->GetPoint3D(m_SoughtVesselSeedName, SoughtVesselSeedPoint);
189 MaskImagePointType MaxSlicePoint;
190 if (GetAFDB()->TagExist(m_SoughtVesselSeedName+"Max")) {
191 GetAFDB()->GetPoint3D(m_SoughtVesselSeedName+"Max", MaxSlicePoint);
194 MaxSlicePoint = SoughtVesselSeedPoint;
195 MaxSlicePoint[2] += 1000;
198 // Find the label with the maximum value to set the result
199 typedef itk::MinimumMaximumImageCalculator<MaskImageType> MinMaxFilterType;
200 MinMaxFilterType::Pointer ff = MinMaxFilterType::New();
202 ff->ComputeMaximum();
203 LabelType newLabel = ff->GetMaximum()+1;
205 // the following bifurcations point will the centroids of the
206 // components obtain when (hopefully!) the SoughtVessel
207 // split into CommonArtery and SubclavianArtery.
208 std::vector<MaskImagePointType> bifurcations;
209 // TrackBifurcationFromPoint(r, m_slice_recon, SoughtVesselSeedPoint,
210 // MaxSlicePoint, newLabel, bifurcations);
212 TrackVesselsFromPoint(r, m_slice_recon, SoughtVesselSeedPoint,
213 MaxSlicePoint, newLabel);
215 // Build the final 3D image from the previous slice by slice processing
216 m_SoughtVessel = clitk::JoinSlices<MaskImageType>(m_slice_recon, m_Mask, 2);
217 // writeImage<MaskImageType>(m_SoughtVessel, "recon2.mhd");
219 // Set binary image, (remove other labels).
221 clitk::Binarize<MaskImageType>(m_SoughtVessel, newLabel, newLabel,
222 GetBackgroundValue(), GetForegroundValue());
224 // writeImage<MaskImageType>(m_SoughtVessel, "afterbinarize.mhd");
225 m_SoughtVessel = clitk::AutoCrop<MaskImageType>(m_SoughtVessel, GetBackgroundValue());
227 // writeImage<MaskImageType>(m_SoughtVessel, "afterautocrop.mhd");
229 // Clean the image : Opening (not in Z direction)
230 typename MaskImageType::SizeType rad;
231 rad[0] = rad[1] = GetFinalOpeningRadius();
233 m_SoughtVessel = clitk::Opening<MaskImageType>(m_SoughtVessel, rad,
234 GetBackgroundValue(), GetForegroundValue());
236 // writeImage<MaskImageType>(m_SoughtVessel, "afteropen.mhd");
238 // Clean the image : keep main CCL slice by slice
239 m_SoughtVessel = clitk::SliceBySliceKeepMainCCL<MaskImageType>(m_SoughtVessel,
240 GetBackgroundValue(),
241 GetForegroundValue());
243 //--------------------------------------------------------------------
246 //--------------------------------------------------------------------
247 template <class TImageType>
249 clitk::ExtractMediastinalVesselsFilter<TImageType>::
250 GenerateInputRequestedRegion() {
251 //DD("GenerateInputRequestedRegion (nothing?)");
253 //--------------------------------------------------------------------
256 //--------------------------------------------------------------------
257 template <class TImageType>
259 clitk::ExtractMediastinalVesselsFilter<TImageType>::
261 // Save in the AFDB (not write on the disk here)
262 GetAFDB()->SetImageFilename(GetSoughtVesselName(), GetOutputFilename());
264 // Final Step -> graft output
265 this->GraftNthOutput(0, m_SoughtVessel);
267 //--------------------------------------------------------------------
270 //--------------------------------------------------------------------
271 template <class TImageType>
273 clitk::ExtractMediastinalVesselsFilter<TImageType>::
275 StartNewStep("Crop the input image: SI,AP limits with carina and crop with mediastinum");
277 Need : Trachea, Carina (roi not point),
279 // Get Trachea and Carina
280 MaskImagePointer Trachea = GetAFDB()->template GetImage <MaskImageType>("Trachea");
282 // Compute Carina position
284 MaskImagePointer Carina = GetAFDB()->template GetImage<MaskImageType>("Carina");
285 std::vector<MaskImagePointType> centroids;
286 clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
287 m_CarinaZ = centroids[1][2];
288 // add one slice to include carina
289 m_CarinaZ += Carina->GetSpacing()[2];
290 // We dont need Carina structure from now
292 GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]);
294 // Crop Inf, remove below Carina
296 clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 2, m_CarinaZ, false, GetBackgroundValue());
298 // // Get seed, crop around
299 // MaskImagePointType SoughtVesselSeedPoint;
300 // GetAFDB()->GetPoint3D(m_SoughtVesselSeedName, SoughtVesselSeedPoint);
303 double m_CarinaY = centroids[1][1];
304 m_Input = clitk::CropImageRemoveGreaterThan<ImageType>(m_Input, 1,
305 m_CarinaY+GetMaxDistancePostToCarina(),
306 false, GetBackgroundValue());
308 m_Input = clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 1,
309 m_CarinaY-GetMaxDistanceAntToCarina(),
310 false, GetBackgroundValue());
312 double m_CarinaX = centroids[1][0];
313 m_Input = clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 0,
314 m_CarinaX-GetMaxDistanceRightToCarina(),
315 false, GetBackgroundValue());
317 m_Input = clitk::CropImageRemoveGreaterThan<ImageType>(m_Input, 0,
318 m_CarinaX+GetMaxDistanceLeftToCarina(),
319 false, GetBackgroundValue());
322 // AutoCrop with Mediastinum, generaly only allow to remove few slices (superiorly)
323 m_Mediastinum = GetAFDB()->template GetImage<MaskImageType>("Mediastinum");
324 // Resize like input (sup to carina)
325 m_Mediastinum = clitk::ResizeImageLike<MaskImageType>(m_Mediastinum, m_Input, GetBackgroundValue());
327 m_Mediastinum = clitk::AutoCrop<MaskImageType>(m_Mediastinum, GetBackgroundValue());
329 m_Input = clitk::ResizeImageLike<ImageType>(m_Input, m_Mediastinum, GetBackgroundValue());
332 // writeImage<ImageType>(m_Input, "crop.mhd");
334 StopCurrentStep<ImageType>(m_Input);
336 //--------------------------------------------------------------------
339 //--------------------------------------------------------------------
340 template <class TImageType>
342 clitk::ExtractMediastinalVesselsFilter<TImageType>::
343 TrackBifurcationFromPoint(MaskImagePointer & recon,
344 std::vector<MaskSlicePointer> & slices_recon,
345 MaskImagePointType point3D,
346 MaskImagePointType pointMaxSlice,
348 std::vector<MaskImagePointType> & bifurcations) {
349 StartNewStep("Track the SoughtVessel from the seed point");
351 // Find first slice index
352 MaskImageIndexType index;
353 recon->TransformPhysicalPointToIndex(point3D, index);
354 int numberOfBifurcation = 0;
355 typedef typename MaskSliceType::PointType SlicePointType;
356 SlicePointType previousCenter;
359 MaskImageIndexType indexMaxSlice;
360 recon->TransformPhysicalPointToIndex(pointMaxSlice, indexMaxSlice);
361 uint maxSlice = indexMaxSlice[2];
363 // Get current label at the point3D of interest
364 uint currentSlice=index[2];
366 LabelType previous_slice_label=recon->GetPixel(index);
369 if (GetVerboseTrackingFlag()) {
370 std::cout << "TrackBifurcationFromPoint " << std::endl;
371 std::cout << "\t point3D = " << point3D << std::endl;
372 std::cout << "\t pointMaxSlice = " << pointMaxSlice << std::endl;
373 std::cout << "\t newLabel = " << newLabel << std::endl;
376 // DD(slices_recon.size());
378 if (GetVerboseTrackingFlag()) {
379 std::cout << "currentSlice = " << currentSlice << std::endl;
381 // Consider current reconstructed slice
382 MaskSlicePointer s = slices_recon[currentSlice];
383 MaskSlicePointer previous;
384 if (currentSlice == index[2]) previous = s;
386 previous = slices_recon[currentSlice-1];
389 // Get centroids of the labels in the current slice
390 static const unsigned int Dim = MaskSliceType::ImageDimension;
391 typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
392 typedef itk::LabelMap< LabelObjectType > LabelMapType;
393 typedef itk::LabelImageToLabelMapFilter<MaskSliceType, LabelMapType> ImageToMapFilterType;
394 typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
395 typedef itk::ShapeLabelMapFilter<LabelMapType, MaskSliceType> ShapeFilterType;
396 typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
397 imageToLabelFilter->SetBackgroundValue(GetBackgroundValue());
398 imageToLabelFilter->SetInput(s);
399 statFilter->SetInput(imageToLabelFilter->GetOutput());
400 // statFilter->SetComputeFeretDiameter( true );
401 statFilter->ComputePerimeterOn(); // To be able to get proper Roundness value
402 statFilter->Update();
403 typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
405 // Look what centroid inside the previous largest one
406 std::vector<SlicePointType> centroids;
407 std::vector<LabelType> centroids_label;
408 std::vector<double> labels_size;
409 for(uint c=0; c<labelMap->GetNumberOfLabelObjects(); c++) {
410 int label = labelMap->GetLabels()[c];
412 SlicePointType center = labelMap->GetLabelObject(label)->GetCentroid();
414 // Get label into previous slice
415 typename MaskSliceType::IndexType centerIndex;
416 previous->TransformPhysicalPointToIndex(center, centerIndex);
417 LabelType labelInPreviousSlice = previous->GetPixel(centerIndex);
418 // if this current centroid was in the current label, add it to bifurcations
419 if (labelInPreviousSlice == previous_slice_label) {
420 centroids.push_back(center);
421 centroids_label.push_back(label);
422 labels_size.push_back(labelMap->GetLabelObject(label)->GetPhysicalSize());
423 //DD(labels_size.back());
424 //DD(labelMap->GetLabelObject(label)->GetRoundness());
425 // previousCenter = centroids.back();
429 // -------------------------
430 // If no centroid were found
431 if (centroids.size() == 0) {
432 if (GetVerboseTrackingFlag()) {
433 std::cout << "no centroid" << std::endl;
435 // Last attempt to find -> check if previous centroid is inside a CCL
436 // if in s -> get value, getcentroid add.
438 //DD("Last change to find");
439 //DD(previous_slice_label);
441 typename MaskSliceType::IndexType previousCenterIndex;
442 s->TransformPhysicalPointToIndex(previousCenter, previousCenterIndex);
443 //DD(previousCenter);
444 LabelType labelInSlice = s->GetPixel(previousCenterIndex);
446 if (labelInSlice != GetBackgroundValue()) {
447 centroids.push_back(labelMap->GetLabelObject(labelInSlice)->GetCentroid());
448 centroids_label.push_back(labelInSlice);
449 labels_size.push_back(labelMap->GetLabelObject(labelInSlice)->GetPhysicalSize());
453 // Some centroid were found
454 // If several centroids, we found a bifurcation
455 if (centroids.size() > 1) {
456 if (GetVerboseTrackingFlag()) {
457 std::cout << "Centroid" << centroids.size() << std::endl;
459 // int n = centroids.size();
461 std::vector<ImagePointType> d;
462 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(centroids, currentSlice, m_Mask, d);
466 // try one or all centroids
468 for(uint a<=0; a<nb++; a++) {
470 // Create the list of candidates
472 if (a==nb) { for(uint x=0; x<nb; x++) c.push_back(x); }
478 for(uint x=0; x<c.size(); c++) { size += labels_size[c[x]]; }
483 // new potential bifurcation found
484 numberOfBifurcation++;
485 // If the number of bifurcation is greater than the required one, we stop
486 if (numberOfBifurcation > GetMaxNumberOfFoundBifurcation()) {
488 //DD("max bif reach");
489 for(uint c=0; c<centroids.size(); c++) {
491 clitk::PointsUtils<MaskImageType>::Convert2DTo3D(centroids[c], m_Mask, currentSlice, bif);
492 bifurcations.push_back(bif);
495 // Else we continue along the main (largest) connected component
497 int indexOfLargest = 0;
498 for(uint b=0; b<centroids.size(); b++) {
499 if (labels_size[b] > labels_size[indexOfLargest]) {
503 //DD(indexOfLargest);
504 //DD(labels_size[indexOfLargest]);
505 SlicePointType c = centroids[indexOfLargest];
506 LabelType l = centroids_label[indexOfLargest];
510 centroids.push_back(c);
511 centroids_label.clear();
512 centroids_label.push_back(l);
517 /* ==> here all centroid are considered as ok.*/
519 // REMOVE IF CENTROID=1, REPLACE BY >0
521 // if only one centroids, we change the current image with the current label
522 if (centroids.size() == 1) {
523 if (GetVerboseTrackingFlag()) {
524 std::cout << "Centroid" << centroids.size() << std::endl;
527 //DD(centroids_label[0]);
528 s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, s, centroids_label[0], newLabel, true);
529 slices_recon[currentSlice] = s;
530 previous_slice_label = newLabel;
531 // It can happend that several CCL share this same label. To
532 // prevent this case, we only consider the one that contains
534 MaskSlicePointer temp = clitk::Binarize<MaskSliceType>(s, newLabel, newLabel, GetBackgroundValue(), GetForegroundValue());
535 // writeImage<MaskSliceType>(temp, "bin-"+toString(currentSlice)+".mhd");
536 temp = clitk::Labelize<MaskSliceType>(temp, GetBackgroundValue(), true, 1);
537 //writeImage<MaskSliceType>(temp, "label-"+toString(currentSlice)+".mhd");
538 typename MaskSliceType::IndexType centroids_index;
539 temp->TransformPhysicalPointToIndex(centroids[0], centroids_index);
540 typename MaskSliceType::PixelType v = temp->GetPixel(centroids_index);
542 // It can happend that the centroid is inside the BG, so we keep
543 // the largest CCL (the first);
544 if (v == GetBackgroundValue()) {
548 v = 1; // largest one
552 temp = clitk::Binarize<MaskSliceType>(temp, v, v, GetBackgroundValue(), newLabel);
553 //writeImage<MaskSliceType>(temp, "relabel-"+toString(currentSlice)+".mhd");
555 slices_recon[currentSlice] = s;
557 // I need to recompute the centroid if we have removed some
558 // connected component.
559 clitk::ComputeCentroids<MaskSliceType>(s, GetBackgroundValue(), centroids);
560 previousCenter = centroids[1];
563 if (GetVerboseTrackingFlag()) {
564 std::cout << "End iteration c=" << centroids.size() << std::endl;
567 if (centroids.size() == 0) {
572 if (currentSlice == slices_recon.size()-1) {
573 // DD("end of slices");
577 if (currentSlice == maxSlice) {
578 // DD("end max slice");
589 //--------------------------------------------------------------------
592 //--------------------------------------------------------------------
593 template <class TImageType>
595 clitk::ExtractMediastinalVesselsFilter<TImageType>::
596 TrackVesselsFromPoint(MaskImagePointer & recon,
597 std::vector<MaskSlicePointer> & slices,
598 MaskImagePointType seedPoint,
599 MaskImagePointType pointMaxSlice,
600 LabelType newLabel) {
601 StartNewStep("Track the SoughtVessel from the seed point");
603 // Find first slice index
604 MaskImageIndexType seedIndex;
605 recon->TransformPhysicalPointToIndex(seedPoint, seedIndex);
606 // int numberOfBifurcation = 0;
607 typedef typename MaskSliceType::PointType SlicePointType;
608 SlicePointType previousCentroid;
609 previousCentroid[0] = seedPoint[0];
610 previousCentroid[1] = seedPoint[1];
613 MaskImageIndexType indexMaxSlice;
614 recon->TransformPhysicalPointToIndex(pointMaxSlice, indexMaxSlice);
615 uint maxSlice = std::min((uint)indexMaxSlice[2], (uint)slices.size());
618 // Get current label at the seedPoint of interest
619 uint currentSlice=seedIndex[2];
621 // LabelType previous_slice_label=recon->GetPixel(seedIndex);
623 // Currrent label map variable
624 typedef itk::ShapeLabelObject< LabelType, 2> LabelObjectType;
625 typedef itk::LabelMap< LabelObjectType > LabelMapType;
626 typename LabelMapType::Pointer labelMap;
627 std::vector<typename LabelObjectType::Pointer> shapeObjectsList;
628 std::vector<MaskSlicePointer> shapeObjectsSliceList; // keep slice, useful for 'union'
629 typename LabelObjectType::Pointer previousShapeObject;
631 if (GetVerboseTrackingFlag()) {
632 std::cout << "TrackBifurcationFromPoint " << std::endl;
633 std::cout << "\t seedPoint = " << seedPoint << std::endl;
634 std::cout << "\t pointMaxSlice = " << pointMaxSlice << std::endl;
635 std::cout << "\t newLabel = " << newLabel << std::endl;
640 //std::cout << std::endl;
643 clitk::PointsUtils<MaskImageType>::Convert2DTo3D(previousCentroid, m_Mask, currentSlice-1, c);
646 if (GetVerboseTrackingFlag()) {
647 std::cout << "Loop slice = " << currentSlice << " c = " << c << std::endl;
650 // Consider current reconstructed slice
651 MaskSlicePointer s = slices[currentSlice];
652 MaskSlicePointer previous;
653 shapeObjectsList.clear();
654 shapeObjectsSliceList.clear();
656 // Get shape of all labels in the current slice (it is already labelized)
658 // Normal -> same CCL with different label
659 // PB -> sometimes same label in different CCL ! car growing
660 //ADD les deux l+s ? mais avec max entre chaque ccl number (bof)
662 for each label in s -> map avec label in l; if different -> change
664 MaskSlicePointer ll = clitk::Labelize<MaskSliceType>(s, GetBackgroundValue(), true, 1);
665 // writeImage<MaskSliceType>(s, "slice-"+toString(currentSlice)+".mhd");
666 //writeImage<MaskSliceType>(ll, "slice-label-"+toString(currentSlice)+".mhd");
667 typedef itk::ImageRegionIteratorWithIndex<MaskSliceType> IteratorType;
668 IteratorType its(s, s->GetLargestPossibleRegion());
669 IteratorType itl(ll, ll->GetLargestPossibleRegion());
670 std::map<LabelType, LabelType> labelInL;
671 std::map<LabelType, std::map<LabelType, LabelType> > labelToChange;
674 int currentLabel = newLabel+10;
675 while (!its.IsAtEnd()) {
676 LabelType labs = its.Get();
677 if (labs != GetBackgroundValue()) {
678 LabelType labl = itl.Get();
679 if (labelInL.find(labs) == labelInL.end()) { // Not found in map, first time
681 labelInL[labs] = labl;
686 if (labelInL[labs] != labl) { // I found a labs with a different labl. Store it.
687 if (labelToChange[labs].find(labl) == labelToChange[labs].end()) { // if not already found
691 //DD(labelInL[labs]);
693 labelToChange[labs][labl] = currentLabel;
705 while (!its.IsAtEnd()) {
706 LabelType labs = its.Get();
707 if (labs != GetBackgroundValue()) { // if not BG
708 LabelType labl = itl.Get();
709 if (labelToChange[labs].find(labl) != labelToChange[labs].end()) { // if some labs can change their label
710 its.Set(labelToChange[labs][labl]); // change with the label for <labs-labl>
717 // writeImage<MaskSliceType>(s, "slice-final"+toString(currentSlice)+".mhd");
720 labelMap = clitk::ComputeLabelMap<MaskSliceType, LabelType>(s, GetBackgroundValue(), true);
721 // DD(labelMap->GetNumberOfLabelObjects());
723 // If this is the first slice, get the object that contains the seed
724 if (currentSlice == seedIndex[2]) {
725 // DD("First slice");
726 LabelType l = recon->GetPixel(seedIndex);
728 shapeObjectsList.push_back(labelMap->GetLabelObject(l));
729 shapeObjectsSliceList.push_back(s);
731 previousCentroid = shapeObjectsList[0]->GetCentroid();
732 previousShapeObject = shapeObjectsList[0];
735 previous = slices[currentSlice-1];
736 // Loop on labels to check if centroid is on the previous contour
737 for(uint c=0; c<labelMap->GetNumberOfLabelObjects(); c++) {
738 // Get the current label number
739 int label = labelMap->GetLabels()[c];
742 SlicePointType centroid = labelMap->GetLabelObject(label)->GetCentroid();
743 // Convert centroid into index in previous slice (same coordinate)
744 typename MaskSliceType::IndexType centroidIndex;
745 previous->TransformPhysicalPointToIndex(centroid, centroidIndex);
746 LabelType labelInPreviousSlice = previous->GetPixel(centroidIndex);
747 // if this current centroid was in the current label, we keep it
748 //DD(labelInPreviousSlice);
749 if (labelInPreviousSlice == newLabel) {
750 shapeObjectsList.push_back(labelMap->GetLabelObject(label));
751 shapeObjectsSliceList.push_back(s);
757 // Potentially the previous centroid could be inside another
758 // labels, we consider i
759 typename MaskSliceType::IndexType previousCentroidIndex;
760 s->TransformPhysicalPointToIndex(previousCentroid, previousCentroidIndex);
761 LabelType l = s->GetPixel(previousCentroidIndex);
763 if (l != 0) { // if is not the background label
765 for(uint c=0; c<shapeObjectsList.size(); c++) {
766 if (shapeObjectsList[c]->GetLabel() == l) {
772 shapeObjectsList.push_back(labelMap->GetLabelObject(l));
773 shapeObjectsSliceList.push_back(s);
776 // DD("already inside");
780 // for(uint c=0; c<shapeObjectsList.size(); c++) {
781 // std::cout << c << " " << shapeObjectsList[c]->GetLabel() << " "
782 // << shapeObjectsList[c]->GetCentroid() << std::endl;
786 // If several candidates, add one more with the union of all candidates
787 MaskSlicePointer temp;
788 if (shapeObjectsList.size() > 1) {
791 temp = clitk::Clone<MaskSliceType>(s);
792 // change label to a single label
793 LabelType l = newLabel+1;
794 for(uint c=0; c<shapeObjectsList.size(); c++) {
795 LabelType ll = shapeObjectsList[c]->GetLabel();
796 temp = clitk::SetBackground<MaskSliceType, MaskSliceType>(temp, temp, ll, l, true);
798 // Compute Label object properties
799 labelMap = clitk::ComputeLabelMap<MaskSliceType, LabelType>(temp, GetBackgroundValue(), true);
800 shapeObjectsList.push_back(labelMap->GetLabelObject(l));
801 shapeObjectsSliceList.push_back(temp);
805 for(uint c=0; c<shapeObjectsList.size(); c++) {
806 std::cout << c << " " << shapeObjectsList[c]->GetLabel() << " "
807 << shapeObjectsList[c]->GetCentroid() << std::endl;
812 for(uint c=0; c<shapeObjectsList.size(); c++) {
814 clitk::PointsUtils<MaskImageType>::Convert2DTo3D(shapeObjectsList[c]->GetCentroid(), m_Mask, currentSlice, cc);
815 // std::cout << c << " " << shapeObjectsList[c]->GetLabel() << " "
816 // // << shapeObjectsList[c]->GetCentroid() << " "
818 // << shapeObjectsList[c]->GetPhysicalSize() << " "
819 // << shapeObjectsList[c]->GetRoundness() << std::endl;
823 if (shapeObjectsList.size() == 0) {
827 // Heuristic to select the good one. For each candidate, we consider the size
828 std::vector<double> sizes;
829 std::vector<double> roundness;
830 std::vector<size_t> index_sizes;
831 std::vector<size_t> index_roundness;
832 double previousSize = previousShapeObject->GetPhysicalSize();
834 for(uint c=0; c<shapeObjectsList.size(); c++) {
835 double s = shapeObjectsList[c]->GetPhysicalSize();
836 sizes.push_back(fabs(previousSize-s)/previousSize);
837 roundness.push_back(fabs(1.0-shapeObjectsList[c]->GetRoundness()));
838 index_sizes.push_back(c);
839 index_roundness.push_back(c);
841 //DDV(sizes, sizes.size());
842 //DDV(roundness, roundness.size());
843 // DDV(index_sizes, index_sizes.size());
844 // DDV(index_roundness, index_roundness.size());
845 sort(index_sizes.begin(), index_sizes.end(), index_cmp<std::vector<double>&>(sizes));
846 sort(index_roundness.begin(), index_roundness.end(), index_cmp<std::vector<double>&>(roundness));
847 //DDV(index_sizes, index_sizes.size());
848 // DDV(index_roundness, index_roundness.size());
850 // TEMPORARY GET THE FIRST
851 int best = index_sizes[0];
852 // if (currentSlice == seedIndex[2]) { // first contour => idiot, first = single contour
853 // best = index_roundness[0]; // best is more round
855 LabelType label = shapeObjectsList[best]->GetLabel();
857 s = shapeObjectsSliceList[best];
858 s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, s, label, newLabel, true);
862 // It can happend that several CCL share this same label. To
863 // prevent this case, we only consider the one that contains
866 // TODO -> PREVIOUS CENTROID ???
868 MaskSlicePointer temp = clitk::Binarize<MaskSliceType>(s, newLabel, newLabel, GetBackgroundValue(), GetForegroundValue());
869 temp = clitk::Labelize<MaskSliceType>(temp, GetBackgroundValue(), true, 1);
870 typename MaskSliceType::IndexType centroids_index;
871 temp->TransformPhysicalPointToIndex(shapeObjectsList[best]->GetCentroid(), centroids_index);
872 typename MaskSliceType::PixelType v = temp->GetPixel(centroids_index);
873 if (v == GetBackgroundValue()) {
877 v = 1; // largest one
881 temp = clitk::Binarize<MaskSliceType>(temp, v, v, GetBackgroundValue(), newLabel);
882 //writeImage<MaskSliceType>(temp, "relabel-"+toString(currentSlice)+".mhd");
887 slices[currentSlice] = s;
888 previousCentroid = shapeObjectsList[best]->GetCentroid();
889 previousShapeObject = shapeObjectsList[best];
894 if (currentSlice == maxSlice) {
895 // DD("end max slice");
902 // -------------------------
903 // If no centroid were found
904 if (shapeObjectsList.size() == 0) {
905 // Last attempt to find -> check if previous centroid is inside a CCL
906 // if in s -> get value, getcentroid add.
908 DD("Last change to find");
909 DD(previous_slice_label);
911 typename MaskSliceType::IndexType previousCentroidIndex;
912 s->TransformPhysicalPointToIndex(previousCentroid, previousCentroidIndex);
913 DD(previousCentroid);
914 LabelType labelInSlice = s->GetPixel(previousCentroidIndex);
916 if (labelInSlice != GetBackgroundValue()) {
917 centroids.push_back(labelMap->GetLabelObject(labelInSlice)->GetCentroid());
918 centroids_label.push_back(labelInSlice);
919 labels_size.push_back(labelMap->GetLabelObject(labelInSlice)->GetPhysicalSize());
923 // Some centroid were found
924 // If several centroids, we found a bifurcation
925 if (centroids.size() > 1) {
926 // int n = centroids.size();
928 std::vector<ImagePointType> d;
929 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(centroids, currentSlice, m_Mask, d);
932 // new potential bifurcation found
933 numberOfBifurcation++;
934 // If the number of bifurcation is greater than the required one, we stop
935 if (numberOfBifurcation > GetMaxNumberOfFoundBifurcation()) {
938 for(uint c=0; c<centroids.size(); c++) {
940 clitk::PointsUtils<MaskImageType>::Convert2DTo3D(centroids[c], m_Mask, currentSlice, bif);
941 bifurcations.push_back(bif);
944 // Else we continue along the main (largest) connected component
946 int indexOfLargest = 0;
947 for(uint b=0; b<centroids.size(); b++) {
948 if (labels_size[b] > labels_size[indexOfLargest]) {
953 DD(labels_size[indexOfLargest]);
954 SlicePointType c = centroids[indexOfLargest];
955 LabelType l = centroids_label[indexOfLargest];
959 centroids.push_back(c);
960 centroids_label.clear();
961 centroids_label.push_back(l);
966 /* ==> here all centroid are considered as ok.*/
969 // REMOVE IF CENTROID=1, REPLACE BY >0
971 // if only one centroids, we change the current image with the current label
972 if (centroids.size() == 1) {
973 DD(centroids_label[0]);
974 s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, s, centroids_label[0], newLabel, true);
975 slices[currentSlice] = s;
976 previous_slice_label = newLabel;
977 // It can happend that several CCL share this same label. To
978 // prevent this case, we only consider the one that contains
980 MaskSlicePointer temp = clitk::Binarize<MaskSliceType>(s, newLabel, newLabel, GetBackgroundValue(), GetForegroundValue());
981 // writeImage<MaskSliceType>(temp, "bin-"+toString(currentSlice)+".mhd");
982 temp = clitk::Labelize<MaskSliceType>(temp, GetBackgroundValue(), true, 1);
983 //writeImage<MaskSliceType>(temp, "label-"+toString(currentSlice)+".mhd");
984 typename MaskSliceType::IndexType centroids_index;
985 temp->TransformPhysicalPointToIndex(centroids[0], centroids_index);
986 typename MaskSliceType::PixelType v = temp->GetPixel(centroids_index);
988 // It can happend that the centroid is inside the BG, so we keep
989 // the largest CCL (the first);
990 if (v == GetBackgroundValue()) {
994 v = 1; // largest one
998 temp = clitk::Binarize<MaskSliceType>(temp, v, v, GetBackgroundValue(), newLabel);
999 //writeImage<MaskSliceType>(temp, "relabel-"+toString(currentSlice)+".mhd");
1001 slices[currentSlice] = s;
1003 // I need to recompute the centroid if we have removed some
1004 // connected component.
1005 clitk::ComputeCentroids<MaskSliceType>(s, GetBackgroundValue(), centroids);
1006 previousCentroid = centroids[1];
1009 if (centroids.size() == 0) {
1014 if (currentSlice == slices.size()-1) {
1015 DD("end of slices");
1019 if (currentSlice == maxSlice) {
1020 DD("end max slice");
1032 //--------------------------------------------------------------------
1035 #endif //#define CLITKEXTRACTMEDIASTINALVESSELSFILTER_TXX