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 "clitkAutoCropFilter.h"
26 #include "clitkSegmentationUtils.h"
27 #include "clitkMorphoMathFilter.h"
30 #include <itkBinaryThresholdImageFilter.h>
31 #include <itkGrayscaleDilateImageFilter.h>
32 #include <itkMinimumMaximumImageCalculator.h>
34 //--------------------------------------------------------------------
35 template <class TImageType>
36 clitk::ExtractMediastinalVesselsFilter<TImageType>::
37 ExtractMediastinalVesselsFilter():
39 clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
40 itk::ImageToImageFilter<TImageType, MaskImageType>()
42 this->SetNumberOfRequiredInputs(1);
43 SetBackgroundValue(0);
44 SetForegroundValue(1);
46 SetTemporaryForegroundValue(1);
48 //--------------------------------------------------------------------
51 //--------------------------------------------------------------------
52 template <class TImageType>
54 clitk::ExtractMediastinalVesselsFilter<TImageType>::
55 GenerateOutputInformation() {
58 m_Input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
60 // ------------------------------------------------------------------
61 // Sup-Inf crop -> Carina
64 // ------------------------------------------------------------------
66 StartNewStep("Binarize with treshold = "+toString(GetThreshold()));
67 typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> BinarizeFilterType;
68 typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
69 binarizeFilter->SetInput(m_Input);
70 binarizeFilter->SetLowerThreshold(GetThreshold());
71 binarizeFilter->SetInsideValue(GetTemporaryForegroundValue());
72 binarizeFilter->SetOutsideValue(GetBackgroundValue());
73 binarizeFilter->Update();
74 m_Mask = binarizeFilter->GetOutput();
75 clitk::writeImage<MaskImageType>(m_Mask, "m.mhd");
76 StopCurrentStep<MaskImageType>(m_Mask);
79 m_Mask = clitk::Labelize<MaskImageType>(m_Mask, GetBackgroundValue(), false, 10);
80 m_Mask = KeepLabels<MaskImageType>(m_Mask, GetBackgroundValue(), GetTemporaryForegroundValue(), 1, 1, true);
81 clitk::writeImage<MaskImageType>(m_Mask, "m2.mhd");
83 // ------------------------------------------------------------------
85 StartNewStep("Detect vessels (slice by slice reconstruction)");
86 std::vector<MaskSlicePointer> slices_mask;
87 clitk::ExtractSlices<MaskImageType>(m_Mask, 2, slices_mask);
89 DD(slices_mask.size());
91 std::vector<MaskSlicePointer> debug_eroded;
92 std::vector<MaskSlicePointer> debug_labeled;
93 std::vector<MaskSlicePointer> debug_slabeled;
96 DD(radius); // TO PUT IN OPTION
98 // ------------------------------------------------------------------
99 // Loop Slice by Slice -> erode find CCL and reconstruct
100 clitk::MorphoMathFilter<MaskSliceType>::Pointer f= clitk::MorphoMathFilter<MaskSliceType>::New();
101 for(uint i=0; i<slices_mask.size(); i++) {
103 f->SetInput(slices_mask[i]);
104 f->SetBackgroundValue(GetBackgroundValue());
105 f->SetForegroundValue(GetTemporaryForegroundValue());
106 f->SetRadius(radius);
107 f->SetOperationType(0); // Erode
110 MaskSlicePointer eroded = f->GetOutput();
111 // writeImage<MaskSliceType>(eroded, "er-"+toString(i)+".mhd");
112 debug_eroded.push_back(eroded);
116 MaskSlicePointer labeled =
117 clitk::LabelizeAndCountNumberOfObjects<MaskSliceType>(eroded, GetBackgroundValue(), false, 1, nb);
119 // Relabel, large CCL with large label number
120 for(int n=nb; n>0; n--) {
123 int lo = 2*(nb+1)-li;
124 labeled = clitk::SetBackground<MaskSliceType, MaskSliceType>(labeled, labeled, li, lo, true);
126 debug_labeled.push_back(labeled);
128 // Create kernel for GrayscaleDilateImageFilter
129 typedef itk::BinaryBallStructuringElement<MaskSliceType::PixelType,MaskSliceType::ImageDimension > KernelType;
131 k.SetRadius(radius+1);
132 k.CreateStructuringElement();
134 // Keep the MAX -> we prefer the opposite su change the label
135 typedef itk::GrayscaleDilateImageFilter<MaskSliceType, MaskSliceType, KernelType> FilterType;
136 FilterType::Pointer m = FilterType::New();
138 m->SetInput(labeled);
139 // DD(m->GetAlgorithm());
140 // m->SetAlgorithm(3);
142 MaskSlicePointer s = m->GetOutput();
146 s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, slices_mask[i],
147 GetBackgroundValue(), GetBackgroundValue(), true);
149 m_slice_recon.push_back(s);
153 MaskImageType::Pointer eroded = clitk::JoinSlices<MaskImageType>(debug_eroded, m_Mask, 2);
154 clitk::writeImage<MaskImageType>(eroded, "eroded.mhd");
157 MaskImageType::Pointer l = clitk::JoinSlices<MaskImageType>(debug_labeled, m_Mask, 2);
158 clitk::writeImage<MaskImageType>(l, "labeled.mhd");
161 MaskImageType::Pointer r = clitk::JoinSlices<MaskImageType>(m_slice_recon, m_Mask, 2);
162 clitk::writeImage<MaskImageType>(r, "recon.mhd");
164 // ------------------------------------------------------------------
165 // Loop Slice by Slice -> BCA not found yet
166 /* MaskImagePointType BCA_p;
167 GetAFDB()->GetPoint3D("BrachioCephalicArteryFirstInferiorPoint", BCA_p);
169 MaskImagePointType bif1;
170 MaskImagePointType bif2;
171 TrackBifurcationFromPoint(r, BCA_p, bif1, bif2);
176 typedef itk::MinimumMaximumImageCalculator<MaskImageType> MinMaxFilterType;
177 MinMaxFilterType::Pointer ff = MinMaxFilterType::New();
179 ff->ComputeMaximum();
180 LabelType newLabel = ff->GetMaximum()+1;
183 // Get all centroids of the first slice
184 std::vector<MaskSlicePointType> centroids2D;
185 clitk::ComputeCentroids<MaskSliceType>(m_slice_recon[0], GetBackgroundValue(), centroids2D);
186 DD(centroids2D.size());
187 std::vector<MaskImagePointType> bifurcations;
188 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(centroids2D, 0, r, bifurcations);
189 DD(bifurcations.size());
190 for(uint i=1; i<bifurcations.size()+1; i++) {
192 DD(bifurcations.size());
193 TrackBifurcationFromPoint(r, m_slice_recon, bifurcations[i], newLabel+i, bifurcations);
195 DD(bifurcations.size());
196 MaskImageType::Pointer rr = clitk::JoinSlices<MaskImageType>(m_slice_recon, m_Mask, 2);
197 clitk::writeImage<MaskImageType>(rr, "recon"+toString(i)+".mhd");
201 //--------------------------------------------------------------------
204 //--------------------------------------------------------------------
205 template <class TImageType>
207 clitk::ExtractMediastinalVesselsFilter<TImageType>::
208 GenerateInputRequestedRegion() {
209 //DD("GenerateInputRequestedRegion (nothing?)");
211 //--------------------------------------------------------------------
214 //--------------------------------------------------------------------
215 template <class TImageType>
217 clitk::ExtractMediastinalVesselsFilter<TImageType>::
220 // Final Step -> graft output (if SetNthOutput => redo)
221 MaskImagePointer BrachioCephalicArtery =
222 GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
223 this->GraftNthOutput(0, BrachioCephalicArtery);
225 //--------------------------------------------------------------------
228 //--------------------------------------------------------------------
229 template <class TImageType>
231 clitk::ExtractMediastinalVesselsFilter<TImageType>::
233 StartNewStep("Inf/Sup limits (carina) and crop with mediastinum");
234 // Get Trachea and Carina
235 MaskImagePointer Trachea = GetAFDB()->template GetImage <MaskImageType>("Trachea");
237 // Get or compute Carina
239 // Get Carina Z position
240 MaskImagePointer Carina = GetAFDB()->template GetImage<MaskImageType>("Carina");
242 std::vector<MaskImagePointType> centroids;
243 clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
244 m_CarinaZ = centroids[1][2];
246 // add one slice to include carina ?
247 m_CarinaZ += Carina->GetSpacing()[2];
248 // We dont need Carina structure from now
250 GetAFDB()->SetDouble("CarinaZ", m_CarinaZ);
253 m_Input = clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 2,
254 m_CarinaZ, false, GetBackgroundValue());
256 // Crop not post to centroid
257 double m_CarinaY = centroids[1][1];
259 m_Input = clitk::CropImageRemoveGreaterThan<ImageType>(m_Input, 1, // OLD ABOVE
260 m_CarinaY, false, GetBackgroundValue());
261 // Crop not ant to centroid
262 m_Input = clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 1,
263 m_CarinaY-80, false, GetBackgroundValue());
265 // AutoCrop with Mediastinum
266 m_Mediastinum = GetAFDB()->template GetImage<MaskImageType>("Mediastinum");
267 // Resize like input (sup to carina)
268 m_Mediastinum = clitk::ResizeImageLike<MaskImageType>(m_Mediastinum, m_Input, GetBackgroundValue());
270 m_Mediastinum = clitk::AutoCrop<MaskImageType>(m_Mediastinum, GetBackgroundValue());
272 m_Input = clitk::ResizeImageLike<ImageType>(m_Input, m_Mediastinum, GetBackgroundValue());
275 StopCurrentStep<ImageType>(m_Input);
277 //--------------------------------------------------------------------
280 //--------------------------------------------------------------------
281 template <class TImageType>
283 clitk::ExtractMediastinalVesselsFilter<TImageType>::
284 //SearchBrachioCephalicArtery(int & BCA_first_slice, LabelType & BCA_first_label) {
285 TrackBifurcationFromPoint(MaskImagePointer & recon,
286 std::vector<MaskSlicePointer> & slices_recon,
287 MaskImagePointType point3D,
289 std::vector<MaskImagePointType> & bifurcations) {
290 StartNewStep("Search for BCA first slice and label");
294 // std::vector<MaskSlicePointer> slices_recon;
295 //clitk::ExtractSlices<MaskImageType>(recon, 2, slices_recon);
298 MaskImageIndexType index;
299 recon->TransformPhysicalPointToIndex(point3D, index);
305 LabelType previous_largest_label=recon->GetPixel(index);
306 DD(previous_largest_label);
309 // Consider current reconstructed slice
310 MaskSlicePointer s = slices_recon[i];
311 MaskSlicePointer previous;
312 if (i==index[2]) previous = s;
313 else previous = slices_recon[i-1];
315 // Get centroids of the labels in the current slice
316 typedef typename MaskSliceType::PointType SlicePointType;
317 static const unsigned int Dim = MaskSliceType::ImageDimension;
318 typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
319 typedef itk::LabelMap< LabelObjectType > LabelMapType;
320 typedef itk::LabelImageToLabelMapFilter<MaskSliceType, LabelMapType> ImageToMapFilterType;
321 typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
322 typedef itk::ShapeLabelMapFilter<LabelMapType, MaskSliceType> ShapeFilterType;
323 typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
324 imageToLabelFilter->SetBackgroundValue(GetBackgroundValue());
325 imageToLabelFilter->SetInput(s);
326 statFilter->SetInput(imageToLabelFilter->GetOutput());
327 statFilter->Update();
328 typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
330 // Look what centroid inside the previous largest one
331 std::vector<SlicePointType> centroids;
332 std::vector<LabelType> centroids_label;
333 for(uint c=0; c<labelMap->GetNumberOfLabelObjects(); c++) {
334 int label = labelMap->GetLabels()[c];
336 SlicePointType center = labelMap->GetLabelObject(label)->GetCentroid();
337 SlicePointType center_previous = center;
338 center_previous[2] -= m_Input->GetSpacing()[2];
339 // Get label into previous slice
340 typename MaskSliceType::IndexType index;
341 previous->TransformPhysicalPointToIndex(center_previous, index);
342 LabelType l = previous->GetPixel(index);
344 if (l == previous_largest_label) {
345 centroids.push_back(center);
346 centroids_label.push_back(label);
349 DD(centroids.size());
351 // If several centroids, we found a bifurcation
352 if (centroids.size() > 1) {
354 for(uint c=0; c<centroids.size(); c++) {
356 clitk::PointsUtils<MaskImageType>::Convert2DTo3D(centroids[c], m_Mask, i, bif);
357 bifurcations.push_back(bif);
358 s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, s, centroids_label[c], newLabel+c+1, true);
359 // slices_recon[i] = s; // (useful ?)
364 // if only one centroids, we change the current image with the current label
365 if (centroids.size() == 1) {
366 s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, s, centroids_label[0], newLabel, true);
367 previous_largest_label = newLabel;
368 /*typedef itk::BinaryThresholdImageFilter<MaskSliceType, MaskSliceType> BinarizeFilterType;
369 typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
370 binarizeFilter->SetInput(s);
371 binarizeFilter->SetLowerThreshold(centroids_label[0]);
372 binarizeFilter->SetUpperThreshold(centroids_label[0]+1);
373 binarizeFilter->SetInsideValue(previous_largest_label);
374 binarizeFilter->SetOutsideValue(GetBackgroundValue());
375 binarizeFilter->Update();
376 s = binarizeFilter->GetOutput();*/
377 slices_recon[i] = s; // (not useful ?)
380 if (centroids.size() == 0) {
381 DD("no centroid, I stop");
385 if (i == slices_recon.size()-1) found = true;
392 //MaskImageType::Pointer rr = clitk::JoinSlices<MaskImageType>(slices_recon, m_Mask, 2);
393 //clitk::writeImage<MaskImageType>(rr, "recon2.mhd");
398 //--------------------------------------------------------------------
402 #endif //#define CLITKEXTRACTMEDIASTINALVESSELSFILTER_TXX