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 //--------------------------------------------------------------------
33 template <class TImageType>
34 clitk::ExtractMediastinalVesselsFilter<TImageType>::
35 ExtractMediastinalVesselsFilter():
37 clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
38 itk::ImageToImageFilter<TImageType, MaskImageType>()
40 this->SetNumberOfRequiredInputs(1);
41 SetBackgroundValue(0);
42 SetForegroundValue(1);
43 SetThresholdHigh(140);
46 SetDilatationRadius(9);
47 SetMaxDistancePostToCarina(10);
48 SetMaxDistanceAntToCarina(40);
49 SetMaxDistanceLeftToCarina(35);
50 SetMaxDistanceRightToCarina(35);
51 SetSoughtVesselSeedName("NoSeedNameGiven");
53 //--------------------------------------------------------------------
56 //--------------------------------------------------------------------
57 template <class TImageType>
59 clitk::ExtractMediastinalVesselsFilter<TImageType>::
60 GenerateOutputInformation() {
63 m_Input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
65 // ------------------------------------------------------------------
66 // Crop the initial image superiorly and inferiorly.
67 // TODO : add options for x cm above/below
70 // ------------------------------------------------------------------
71 // Binarize the image. Need two thresholds, one high to select
72 // structures (CCL) that are almost not connected (after erosion),
73 // and one low thresholds to select the real contours. Will be
74 // reconstructed later.
75 StartNewStep("Binarize with high threshold = "+toString(GetThresholdHigh()));
76 typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> BinarizeFilterType;
77 typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
78 binarizeFilter->SetInput(m_Input);
79 binarizeFilter->SetLowerThreshold(GetThresholdHigh());
80 binarizeFilter->SetInsideValue(GetForegroundValue());
81 binarizeFilter->SetOutsideValue(GetBackgroundValue());
82 binarizeFilter->Update();
83 m_Mask = binarizeFilter->GetOutput();
84 StopCurrentStep<MaskImageType>(m_Mask);
86 // ------------------------------------------------------------------
87 StartNewStep("Binarize with low threshold = "+toString(GetThresholdLow()));
88 binarizeFilter = BinarizeFilterType::New();
89 binarizeFilter->SetInput(m_Input);
90 binarizeFilter->SetLowerThreshold(GetThresholdLow());
91 binarizeFilter->SetInsideValue(GetForegroundValue());
92 binarizeFilter->SetOutsideValue(GetBackgroundValue());
93 binarizeFilter->Update();
94 MaskImagePointer m_Mask2 = binarizeFilter->GetOutput();
95 StopCurrentStep<MaskImageType>(m_Mask2);
97 // ------------------------------------------------------------------
99 StartNewStep("Detect objects : erosion, then slice by slice reconstruction");
100 std::vector<MaskSlicePointer> slices_mask;
101 clitk::ExtractSlices<MaskImageType>(m_Mask, 2, slices_mask);
102 std::vector<MaskSlicePointer> slices_mask2;
103 clitk::ExtractSlices<MaskImageType>(m_Mask2, 2, slices_mask2);
104 int radius = GetErosionRadius();
106 // List of working slices (debug only)
107 std::vector<MaskSlicePointer> debug_eroded;
108 // std::vector<MaskSlicePointer> debug_labeled;
110 // ------------------------------------------------------------------
111 // Loop Slice by Slice in order to break connectivity between
112 // CCL. Erode and reconstruct all labels at the same time without
114 for(uint i=0; i<slices_mask.size(); i++) {
116 typedef itk::BinaryBallStructuringElement<MaskSliceType::PixelType,2> KernelType;
117 KernelType structuringElement;
118 structuringElement.SetRadius(radius);
119 structuringElement.CreateStructuringElement();
121 // Erosion -> we break the connectivity between structure
122 typedef itk::BinaryErodeImageFilter<MaskSliceType, MaskSliceType, KernelType> ErodeFilterType;
123 typename ErodeFilterType::Pointer eroder = ErodeFilterType::New();
124 eroder->SetInput(slices_mask[i]);
125 eroder->SetBackgroundValue(GetBackgroundValue());
126 eroder->SetForegroundValue(GetForegroundValue());
127 eroder->SetBoundaryToForeground(true); // ??
128 eroder->SetKernel(structuringElement);
130 MaskSlicePointer eroded = eroder->GetOutput();
131 debug_eroded.push_back(eroded);
134 MaskSlicePointer labeled =
135 clitk::Labelize<MaskSliceType>(eroded, GetBackgroundValue(), true, 1); // Fully connected !
136 // debug_labeled.push_back(labeled);
138 // Make Reconstruction filter : dilation all labels at the same
139 // time, prevent to merge them.
140 typedef clitk::ReconstructWithConditionalGrayscaleDilateImageFilter<MaskSliceType> ReconstructFilterType;
141 typename ReconstructFilterType::Pointer reconstructor = ReconstructFilterType::New();
142 reconstructor->SetInput(labeled);
143 reconstructor->SetIterationNumber(radius+GetDilatationRadius());
144 reconstructor->Update();
145 MaskSlicePointer s = reconstructor->GetOutput();
147 // Remove Initial BG of the second tresholded image
148 s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, slices_mask2[i],
149 GetBackgroundValue(), GetBackgroundValue(), true);
150 m_slice_recon.push_back(s);
154 // Build 3D images from the slice by slice processing
155 MaskImageType::Pointer eroded = clitk::JoinSlices<MaskImageType>(debug_eroded, m_Mask, 2);
156 writeImage<MaskImageType>(eroded, "erode.mhd");
157 //MaskImageType::Pointer l = clitk::JoinSlices<MaskImageType>(debug_labeled, m_Mask, 2);
158 MaskImageType::Pointer r = clitk::JoinSlices<MaskImageType>(m_slice_recon, m_Mask, 2);
159 writeImage<MaskImageType>(r, "recon1.mhd");
161 // ------------------------------------------------------------------
162 // Track the SoughtVessel from the given first point
163 // superiorly. This is done by TrackBifurcationFromPoint
164 MaskImagePointType SoughtVesselSeedPoint;
165 GetAFDB()->GetPoint3D(m_SoughtVesselSeedName, SoughtVesselSeedPoint);
167 // Find the label with the maximum value to set the result
168 typedef itk::MinimumMaximumImageCalculator<MaskImageType> MinMaxFilterType;
169 MinMaxFilterType::Pointer ff = MinMaxFilterType::New();
171 ff->ComputeMaximum();
172 LabelType newLabel = ff->GetMaximum()+1;
174 // the following bifurcations point will the centroids of the
175 // components obtain when (hopefully!) the SoughtVessel
176 // split into CommonArtery and SubclavianArtery.
177 std::vector<MaskImagePointType> bifurcations;
178 TrackBifurcationFromPoint(r, m_slice_recon, SoughtVesselSeedPoint, newLabel, bifurcations);
180 // Build the final 3D image from the previous slice by slice processing
181 m_SoughtVessel = clitk::JoinSlices<MaskImageType>(m_slice_recon, m_Mask, 2);
182 writeImage<MaskImageType>(m_SoughtVessel, "recon2.mhd");
184 // Set binary image, (remove other labels).
185 // TODO: keep labeled image to track SubclavianArtery and CommonArtery
187 clitk::Binarize<MaskImageType>(m_SoughtVessel, newLabel, newLabel,
188 GetBackgroundValue(), GetForegroundValue());
190 writeImage<MaskImageType>(m_SoughtVessel, "afterbinarize.mhd");
192 m_SoughtVessel = clitk::AutoCrop<MaskImageType>(m_SoughtVessel, GetBackgroundValue());
194 writeImage<MaskImageType>(m_SoughtVessel, "afterautocrop.mhd");
196 // Clean the image : Opening (not in Z direction)
197 typename MaskImageType::SizeType rad;
200 m_SoughtVessel = clitk::Opening<MaskImageType>(m_SoughtVessel, rad,
201 GetBackgroundValue(), GetForegroundValue());
203 writeImage<MaskImageType>(m_SoughtVessel, "afteropen.mhd");
205 // Clean the image : keep main CCL slice by slice
206 m_SoughtVessel = clitk::SliceBySliceKeepMainCCL<MaskImageType>(m_SoughtVessel,
207 GetBackgroundValue(),
208 GetForegroundValue());
210 //--------------------------------------------------------------------
213 //--------------------------------------------------------------------
214 template <class TImageType>
216 clitk::ExtractMediastinalVesselsFilter<TImageType>::
217 GenerateInputRequestedRegion() {
218 //DD("GenerateInputRequestedRegion (nothing?)");
220 //--------------------------------------------------------------------
223 //--------------------------------------------------------------------
224 template <class TImageType>
226 clitk::ExtractMediastinalVesselsFilter<TImageType>::
228 // Save in the AFDB (not write on the disk here)
229 GetAFDB()->SetImageFilename(GetSoughtVesselName(), GetOutputFilename());
231 // Final Step -> graft output
232 this->GraftNthOutput(0, m_SoughtVessel);
234 //--------------------------------------------------------------------
237 //--------------------------------------------------------------------
238 template <class TImageType>
240 clitk::ExtractMediastinalVesselsFilter<TImageType>::
242 StartNewStep("Crop the input image: SI,AP limits with carina and crop with mediastinum");
244 Need : Trachea, Carina (roi not point),
246 // Get Trachea and Carina
247 MaskImagePointer Trachea = GetAFDB()->template GetImage <MaskImageType>("Trachea");
249 // Compute Carina position
251 MaskImagePointer Carina = GetAFDB()->template GetImage<MaskImageType>("Carina");
252 std::vector<MaskImagePointType> centroids;
253 clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
254 m_CarinaZ = centroids[1][2];
255 // add one slice to include carina
256 m_CarinaZ += Carina->GetSpacing()[2];
257 // We dont need Carina structure from now
259 GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]);
261 // Crop Inf, remove below Carina
263 clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 2, m_CarinaZ, false, GetBackgroundValue());
266 double m_CarinaY = centroids[1][1];
267 m_Input = clitk::CropImageRemoveGreaterThan<ImageType>(m_Input, 1,
268 m_CarinaY+GetMaxDistancePostToCarina(),
269 false, GetBackgroundValue());
271 m_Input = clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 1,
272 m_CarinaY-GetMaxDistanceAntToCarina(),
273 false, GetBackgroundValue());
276 double m_CarinaX = centroids[1][0];
277 m_Input = clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 0,
278 m_CarinaX-GetMaxDistanceRightToCarina(),
279 false, GetBackgroundValue());
281 m_Input = clitk::CropImageRemoveGreaterThan<ImageType>(m_Input, 0,
282 m_CarinaX+GetMaxDistanceLeftToCarina(),
283 false, GetBackgroundValue());
286 // AutoCrop with Mediastinum, generaly only allow to remove few slices (superiorly)
287 m_Mediastinum = GetAFDB()->template GetImage<MaskImageType>("Mediastinum");
288 // Resize like input (sup to carina)
289 m_Mediastinum = clitk::ResizeImageLike<MaskImageType>(m_Mediastinum, m_Input, GetBackgroundValue());
291 m_Mediastinum = clitk::AutoCrop<MaskImageType>(m_Mediastinum, GetBackgroundValue());
293 m_Input = clitk::ResizeImageLike<ImageType>(m_Input, m_Mediastinum, GetBackgroundValue());
297 StopCurrentStep<ImageType>(m_Input);
299 //--------------------------------------------------------------------
302 //--------------------------------------------------------------------
303 template <class TImageType>
305 clitk::ExtractMediastinalVesselsFilter<TImageType>::
306 TrackBifurcationFromPoint(MaskImagePointer & recon,
307 std::vector<MaskSlicePointer> & slices_recon,
308 MaskImagePointType point3D,
310 std::vector<MaskImagePointType> & bifurcations) {
311 StartNewStep("Track the SoughtVessel from the seed point");
313 // Find first slice index
314 MaskImageIndexType index;
315 recon->TransformPhysicalPointToIndex(point3D, index);
316 int numberOfBifurcation = 0;
317 typedef typename MaskSliceType::PointType SlicePointType;
318 SlicePointType previousCenter;
320 // Get current label at the point3D of interest
321 uint currentSlice=index[2];
323 LabelType previous_slice_label=recon->GetPixel(index);
324 // DD(slices_recon.size());
327 // Consider current reconstructed slice
328 MaskSlicePointer s = slices_recon[currentSlice];
329 MaskSlicePointer previous;
330 if (currentSlice == index[2]) previous = s;
332 previous = slices_recon[currentSlice-1];
335 // Get centroids of the labels in the current slice
336 static const unsigned int Dim = MaskSliceType::ImageDimension;
337 typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
338 typedef itk::LabelMap< LabelObjectType > LabelMapType;
339 typedef itk::LabelImageToLabelMapFilter<MaskSliceType, LabelMapType> ImageToMapFilterType;
340 typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
341 typedef itk::ShapeLabelMapFilter<LabelMapType, MaskSliceType> ShapeFilterType;
342 typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
343 imageToLabelFilter->SetBackgroundValue(GetBackgroundValue());
344 imageToLabelFilter->SetInput(s);
345 statFilter->SetInput(imageToLabelFilter->GetOutput());
346 statFilter->Update();
347 typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
349 // Look what centroid inside the previous largest one
350 std::vector<SlicePointType> centroids;
351 std::vector<LabelType> centroids_label;
352 std::vector<double> labels_size;
353 for(uint c=0; c<labelMap->GetNumberOfLabelObjects(); c++) {
354 int label = labelMap->GetLabels()[c];
356 SlicePointType center = labelMap->GetLabelObject(label)->GetCentroid();
358 // Get label into previous slice
359 typename MaskSliceType::IndexType centerIndex;
360 previous->TransformPhysicalPointToIndex(center, centerIndex);
361 LabelType labelInPreviousSlice = previous->GetPixel(centerIndex);
362 // if this current centroid was in the current label, add it to bifurcations
363 if (labelInPreviousSlice == previous_slice_label) {
364 centroids.push_back(center);
365 centroids_label.push_back(label);
366 labels_size.push_back(labelMap->GetLabelObject(label)->GetPhysicalSize());
370 if (centroids.size() == 0) {
371 // Last attempt to find -> check if previous centroid is inside a CCL
372 // if in s -> get value, getcentroid add.
374 DD("Last change to find");
375 typename MaskSliceType::IndexType previousCenterIndex;
376 s->TransformPhysicalPointToIndex(previousCenter, previousCenterIndex);
378 LabelType labelInSlice = s->GetPixel(previousCenterIndex);
380 if (labelInSlice != GetBackgroundValue()) {
381 centroids.push_back(labelMap->GetLabelObject(labelInSlice)->GetCentroid());
382 centroids_label.push_back(labelInSlice);
383 labels_size.push_back(labelMap->GetLabelObject(labelInSlice)->GetPhysicalSize());
388 // DD(centroids.size());
390 // If several centroids, we found a bifurcation
391 if (centroids.size() > 1) {
392 numberOfBifurcation++;
393 // If the number of bifurcation is greater than the required one, we stop
394 if (numberOfBifurcation > GetMaxNumberOfFoundBifurcation()) {
397 for(uint c=0; c<centroids.size(); c++) {
399 clitk::PointsUtils<MaskImageType>::Convert2DTo3D(centroids[c], m_Mask, currentSlice, bif);
400 bifurcations.push_back(bif);
403 // Else we continue along the main (largest) connected component
405 int indexOfLargest = 0;
406 for(uint b=0; b<centroids.size(); b++) {
407 if (labels_size[b] > labels_size[indexOfLargest]) {
411 SlicePointType c = centroids[indexOfLargest];
412 LabelType l = centroids_label[indexOfLargest];
414 centroids.push_back(c);
415 centroids_label.push_back(l);
419 // if only one centroids, we change the current image with the current label
420 if (centroids.size() == 1) {
421 s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, s, centroids_label[0], newLabel, true);
422 slices_recon[currentSlice] = s;
423 previous_slice_label = newLabel;
424 // It can happend that several CCL share this same label. To
425 // prevent this case, we only consider the one that contains
427 MaskSlicePointer temp = clitk::Binarize<MaskSliceType>(s, newLabel, newLabel, GetBackgroundValue(), GetForegroundValue());
428 // writeImage<MaskSliceType>(temp, "bin-"+toString(currentSlice)+".mhd");
429 temp = clitk::Labelize<MaskSliceType>(temp, GetBackgroundValue(), true, 1);
430 //writeImage<MaskSliceType>(temp, "label-"+toString(currentSlice)+".mhd");
431 typename MaskSliceType::IndexType centroids_index;
432 temp->TransformPhysicalPointToIndex(centroids[0], centroids_index);
433 typename MaskSliceType::PixelType v = temp->GetPixel(centroids_index);
435 // It can happend that the centroid is inside the BG, so we keep
436 // the largest CCL (the first);
437 if (v == GetBackgroundValue()) {
441 v = 1; // largest one
445 temp = clitk::Binarize<MaskSliceType>(temp, v, v, GetBackgroundValue(), newLabel);
446 //writeImage<MaskSliceType>(temp, "relabel-"+toString(currentSlice)+".mhd");
448 slices_recon[currentSlice] = s;
450 // I need to recompute the centroid if we have removed some
451 // connected component.
452 clitk::ComputeCentroids<MaskSliceType>(s, GetBackgroundValue(), centroids);
453 previousCenter = centroids[1];
456 if (centroids.size() == 0) {
461 if (currentSlice == slices_recon.size()-1) {
473 //--------------------------------------------------------------------
476 #endif //#define CLITKEXTRACTMEDIASTINALVESSELSFILTER_TXX