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://www.centreleonberard.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 CLITKEXTRACTAIRWAYSTREEINFOSFILTER_TXX
20 #define CLITKEXTRACTAIRWAYSTREEINFOSFILTER_TXX
23 #include "clitkImageCommon.h"
24 #include "clitkSetBackgroundImageFilter.h"
25 #include "clitkSegmentationUtils.h"
26 #include "clitkAutoCropFilter.h"
27 #include "clitkExtractSliceFilter.h"
30 #include "itkBinaryThresholdImageFilter.h"
31 #include "itkConnectedComponentImageFilter.h"
32 #include "itkRelabelComponentImageFilter.h"
33 #include "itkOtsuThresholdImageFilter.h"
34 #include "itkBinaryThinningImageFilter3D.h"
35 #include "itkImageIteratorWithIndex.h"
38 //--------------------------------------------------------------------
39 template <class ImageType>
40 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
41 ExtractAirwaysTreeInfoFilter():
43 clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
44 itk::ImageToImageFilter<ImageType, ImageType>()
46 // Default global options
47 this->SetNumberOfRequiredInputs(1);
48 SetBackgroundValue(0);
49 SetForegroundValue(1);
51 //--------------------------------------------------------------------
54 //--------------------------------------------------------------------
55 template <class ImageType>
57 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
58 SetInput(const ImageType * image)
60 this->SetNthInput(0, const_cast<ImageType *>(image));
62 //--------------------------------------------------------------------
65 //--------------------------------------------------------------------
66 template <class ImageType>
68 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
69 GenerateOutputInformation()
71 Superclass::GenerateOutputInformation();
72 //this->GetOutput(0)->SetRequestedRegion(this->GetOutput(0)->GetLargestPossibleRegion());
75 input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
80 } catch (clitk::ExceptionObject e) {
81 // Do nothing if not found, it will be used anyway to write the result
86 //--------------------------------------------------------------------
89 //--------------------------------------------------------------------
90 template <class ImageType>
92 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
96 // Crop just abov lungs ?
98 // read skeleton if already exist
99 bool foundSkeleton = true;
101 skeleton = GetAFDB()->template GetImage <ImageType>("skeleton");
103 catch (clitk::ExceptionObject e) {
104 DD("did not found skeleton");
105 foundSkeleton = false;
109 if (!foundSkeleton) {
110 StartNewStep("Thinning filter (skeleton)");
111 typedef itk::BinaryThinningImageFilter3D<ImageType, ImageType> ThinningFilterType;
112 typename ThinningFilterType::Pointer thinningFilter = ThinningFilterType::New();
113 thinningFilter->SetInput(input); // input = trachea
114 thinningFilter->Update();
115 skeleton = thinningFilter->GetOutput();
116 StopCurrentStep<ImageType>(skeleton);
117 writeImage<ImageType>(skeleton, "skeleton.mhd");
120 // Find first point for tracking
121 StartNewStep("Find first point for tracking");
122 typedef itk::ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
123 IteratorType it(skeleton, skeleton->GetLargestPossibleRegion());
124 it.GoToReverseBegin();
125 while ((!it.IsAtEnd()) && (it.Get() == GetBackgroundValue())) {
129 clitkExceptionMacro("first point in the trachea skeleton not found.");
131 typename ImageType::IndexType index = it.GetIndex();
135 // Initialize neighborhooditerator
136 typedef itk::NeighborhoodIterator<ImageType> NeighborhoodIteratorType;
137 typename NeighborhoodIteratorType::SizeType radius;
139 NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
142 // Find a label number different from BG and FG
143 typename ImageType::PixelType label = GetForegroundValue()+1;
144 while ((label == GetBackgroundValue()) || (label == GetForegroundValue())) { label++; }
149 Goal : find position of C, RUL, RML, RLL, LUL, LLL bronchus
150 Carina : ok "easy", track, slice by slice until 2 path into different label
156 // NEW Track from the first point
157 StartNewStep("Start tracking");
160 skeleton->TransformIndexToPhysicalPoint(n.index, n.point); // à mettre dans TrackFromThisIndex2
161 mFullSkeletonTree.set_head(n);
162 StructuralTreeNodeType sn;
164 skeleton->TransformIndexToPhysicalPoint(sn.index, sn.point);
165 mStructuralSkeletonTree.set_head(sn);
166 TrackFromThisIndex2(index, skeleton, label, mFullSkeletonTree.begin(), mStructuralSkeletonTree.begin());
169 // Reput FG instead of label in the skeleton image
170 skeleton = clitk::SetBackground<ImageType, ImageType>(skeleton, skeleton, label, GetForegroundValue(), true);
173 typename StructuralTreeType::iterator sit = mStructuralSkeletonTree.begin();
174 while (sit != mStructuralSkeletonTree.end()) {
179 // compute weight : n longueurs à chaque bifurfaction.
180 // parcours de FullTreeNodeType, à partir de leaf, remonter de proche en proche la distance eucl.
184 typename FullTreeType::iterator fit = mFullSkeletonTree.begin();
185 while (fit != mFullSkeletonTree.end()) {
190 DD("compute weight");
191 typename FullTreeType::post_order_iterator pit = mFullSkeletonTree.begin_post();
192 while (pit != mFullSkeletonTree.end_post()) {
195 if (pit != mFullSkeletonTree.begin()) {
196 typename FullTreeType::iterator parent = mFullSkeletonTree.parent(pit);
197 double d = pit->point.EuclideanDistanceTo(parent->point);
198 // DD(parent->weight);
201 if (pit.number_of_children() > 1) {
202 DD(pit.number_of_children());
206 DD(mFullSkeletonTree.child(pit, 0)->weight);
207 DD(mFullSkeletonTree.child(pit, 1)->weight);
211 double previous = pit->weight;
212 for(uint i=0; i<pit.number_of_children(); i++) {
213 double d = pit->point.EuclideanDistanceTo(mFullSkeletonTree.child(pit, i)->point);
214 // pit->weight = pit->weight + mFullSkeletonTree.child(pit, i)->weight + d;
216 pit->weight = pit->weight + mFullSkeletonTree.child(pit, i)->weight + d;
218 DD(pit.number_of_children());
221 DD(mFullSkeletonTree.child(pit, 0)->weight);
222 DD(mFullSkeletonTree.child(pit, 1)->weight);
223 pit->weight = std::max(pit->weight, previous+mFullSkeletonTree.child(pit, i)->weight + d);
230 fit = mFullSkeletonTree.begin();
231 while (fit != mFullSkeletonTree.end()) {
232 std::cout << "p = " << fit->point
233 << " " << fit->weight
234 << " child= " << fit.number_of_children()
236 for(uint i=0; i<fit.number_of_children(); i++) {
237 std::cout << " " << mFullSkeletonTree.child(fit, i)->weight;
239 std::cout << std::endl;
243 /* Selection criteria ?
248 DD("=========================================");
249 fit = mFullSkeletonTree.begin();
250 while (fit != mFullSkeletonTree.end()) {
251 if (fit.number_of_children() > 1) {
252 std::cout << "ppp = " << fit->point
253 << " " << fit->weight << std::endl;
254 for(uint i=1; i<fit.number_of_children(); i++) {
255 double w1 = mFullSkeletonTree.child(fit, i-1)->weight;
256 double w2 = mFullSkeletonTree.child(fit, i)->weight;
257 // if (w1 <= 0.1) break;
258 // if (w2 <= 0.1) break;
275 // Track from the first point
276 StartNewStep("Start tracking OLD");
277 std::vector<BifurcationType> listOfBifurcations;
278 m_SkeletonTree.set_head(index);
279 TrackFromThisIndex(listOfBifurcations, skeleton, index, label, m_SkeletonTree.begin());
282 // Convert index into physical point coordinates
283 for(unsigned int i=0; i<listOfBifurcations.size(); i++) {
284 skeleton->TransformIndexToPhysicalPoint(listOfBifurcations[i].index,
285 listOfBifurcations[i].point);
288 // Search for the first slice that separate the bronchus
289 // (carina). Labelize slice by slice, stop when the two points of
290 // the skeleton ar not in the same connected component
291 StartNewStep("Search for carina position");
292 typedef clitk::ExtractSliceFilter<ImageType> ExtractSliceFilterType;
293 typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
294 extractSliceFilter->SetInput(input);
295 extractSliceFilter->SetDirection(2);
296 extractSliceFilter->Update();
297 typedef typename ExtractSliceFilterType::SliceType SliceType;
298 std::vector<typename SliceType::Pointer> mInputSlices;
299 extractSliceFilter->GetOutputSlices(mInputSlices);
300 DD(mInputSlices.size());
303 DD("REDO !!!!!!!!!!!!");
305 => chercher la bif qui a les plus important sous-arbres
306 - iterate sur m_SkeletonTree depth first, remonter au parent, add value = distance
308 NONNNNN m_SkeletonTree n'est pas la structure mais l'ensemble avec ts le skeleton
310 TreeIterator itt = m_SkeletonTree.begin();
311 while (itt != m_SkeletonTree.end()) {
316 typename tree<NodeType>::post_order_iterator pit = m_SkeletonTree.begin_post();
317 while (pit != m_SkeletonTree.end_post()) {
318 typename tree<NodeType>::iterator parent = m_SkeletonTree.parent(pit);
319 double d = pit->point.EuclideanDistanceTo(parent);
326 itt = m_SkeletonTree.begin();
327 while (itt != m_SkeletonTree.end()) {
336 int slice_index = listOfBifurcations[0].index[2]; // first slice from carina in skeleton
338 TreeIterator firstIter = m_SkeletonTree.child(listOfBifurcations[0].treeIter, 0);
339 TreeIterator secondIter = m_SkeletonTree.child(listOfBifurcations[0].treeIter, 1);
340 DD(firstIter.number_of_children());
341 DD(secondIter.number_of_children());
342 typename SliceType::IndexType in1;
343 typename SliceType::IndexType in2;
345 // Labelize the current slice
346 typename SliceType::Pointer temp = Labelize<SliceType>(mInputSlices[slice_index],
347 GetBackgroundValue(),
349 0); // min component size=0
352 // Check the value of the two skeleton points;
353 in1[0] = (*firstIter)[0];
354 in1[1] = (*firstIter)[1];
355 typename SliceType::PixelType v1 = temp->GetPixel(in1);
356 in2[0] = (*secondIter)[0];
357 in2[1] = (*secondIter)[1];
358 typename SliceType::PixelType v2 = temp->GetPixel(in2);
360 // Check the label value of the two points
363 stop = true; // We found it !
367 if (slice_index == (int)(mInputSlices.size()-1)) {
368 clitkExceptionMacro("Error while searching for carina, the two skeleton points are always in the same CC ... ???");
377 ImageIndexType carina_index; // middle position in X/Y
378 carina_index[0] = lrint(in2[0] + in1[0])/2.0;
379 carina_index[1] = lrint(in2[1] + in1[1])/2.0;
380 carina_index[2] = slice_index;
381 // Get physical coordinates
382 ImagePointType carina_position;
383 skeleton->TransformIndexToPhysicalPoint(carina_index,
386 // Set and save Carina position
387 if (GetVerboseStepFlag()) {
388 std::cout << "\t Found carina at " << carina_position << " mm" << std::endl;
390 GetAFDB()->SetPoint3D("carina", carina_position);
392 // Write bifurcation (debug)
393 for(uint i=0; i<listOfBifurcations.size(); i++) {
394 GetAFDB()->SetPoint3D("Bif"+toString(i), listOfBifurcations[i].point);
398 // Set the output (skeleton);
399 this->GraftOutput(skeleton); // not SetNthOutput
401 //--------------------------------------------------------------------
404 //--------------------------------------------------------------------
406 Start from the pixel "index" in the image "skeleton" and track the
407 path from neighbor pixels. Put the tracked path in the tree at the
408 currentNode position. Label is used to mark the FG pixel as already
409 visited. Progress recursively when several neighbors are found.
411 template <class ImageType>
413 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
414 TrackFromThisIndex2(ImageIndexType index, ImagePointer skeleton,
415 ImagePixelType label,
416 typename FullTreeType::iterator currentNode,
417 typename StructuralTreeType::iterator currentSNode)
419 // Create NeighborhoodIterator
420 typedef itk::NeighborhoodIterator<ImageType> NeighborhoodIteratorType;
421 typename NeighborhoodIteratorType::SizeType radius;
423 NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
426 std::vector<typename NeighborhoodIteratorType::IndexType> listOfTrackedPoint;
429 nit.SetLocation(index);
430 nit.SetCenterPixel(label);
431 listOfTrackedPoint.clear();
432 for(unsigned int i=0; i<nit.Size(); i++) {
433 if (i != nit.GetCenterNeighborhoodIndex ()) { // Do not observe the current point
434 if (nit.GetPixel(i) == GetForegroundValue()) {
435 // if pixel value is FG, we keep this point
436 listOfTrackedPoint.push_back(nit.GetIndex(i));
440 // If only neighbor is found, we keep the point and continue
441 if (listOfTrackedPoint.size() == 1) {
443 index = n.index = listOfTrackedPoint[0];
444 skeleton->TransformIndexToPhysicalPoint(n.index, n.point);
445 currentNode = mFullSkeletonTree.append_child(currentNode, n);
446 skeleton->SetPixel(n.index, label); // change label in skeleton image
449 if (listOfTrackedPoint.size() >= 2) {
450 for(uint i=0; i<listOfTrackedPoint.size(); i++) {
452 n.index = listOfTrackedPoint[i];
453 skeleton->TransformIndexToPhysicalPoint(n.index, n.point);
454 typename FullTreeType::iterator node = mFullSkeletonTree.append_child(currentNode, n);
455 StructuralTreeNodeType sn;
456 sn.index = listOfTrackedPoint[i];
457 skeleton->TransformIndexToPhysicalPoint(sn.index, sn.point);
458 typename StructuralTreeType::iterator snode = mStructuralSkeletonTree.append_child(currentSNode, sn);
459 TrackFromThisIndex2(listOfTrackedPoint[i], skeleton, label, node, snode);
462 stop = true; // this it the end of the tracking
466 //--------------------------------------------------------------------
468 //--------------------------------------------------------------------
469 template <class ImageType>
471 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
472 TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations,
473 ImagePointer skeleton,
474 ImageIndexType index,
475 ImagePixelType label,
476 TreeIterator currentNode)
478 // Create NeighborhoodIterator
479 typedef itk::NeighborhoodIterator<ImageType> NeighborhoodIteratorType;
480 typename NeighborhoodIteratorType::SizeType radius;
482 NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
485 std::vector<typename NeighborhoodIteratorType::IndexType> listOfTrackedPoint;
488 nit.SetLocation(index);
489 nit.SetCenterPixel(label);
490 listOfTrackedPoint.clear();
491 for(unsigned int i=0; i<nit.Size(); i++) {
492 if (i != nit.GetCenterNeighborhoodIndex ()) { // Do not observe the current point
493 if (nit.GetPixel(i) == GetForegroundValue()) { // if this is foreground, we continue the tracking
494 listOfTrackedPoint.push_back(nit.GetIndex(i));
498 if (listOfTrackedPoint.size() == 1) {
499 // Add this point to the current path
500 currentNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[0]);
501 index = listOfTrackedPoint[0];
502 skeleton->SetPixel(index, label); // change label in skeleton image
505 if (listOfTrackedPoint.size() == 2) {
506 // m_SkeletonTree->Add(listOfTrackedPoint[0], index); // the parent is 'index'
507 // m_SkeletonTree->Add(listOfTrackedPoint[1], index); // the parent is 'index'
508 DD("BifurcationType");
509 DD(listOfTrackedPoint[0]);
510 DD(listOfTrackedPoint[1]);
511 BifurcationType bif(index, label, label+1, label+2);
512 bif.treeIter = currentNode;
513 listOfBifurcations.push_back(bif);
514 TreeIterator firstNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[0]);
515 TreeIterator secondNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[1]);
516 TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1, firstNode);
517 TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2, secondNode);
520 //DD(listOfTrackedPoint.size());
521 if (listOfTrackedPoint.size() > 2) {
522 //clitkExceptionMacro("error while tracking trachea bifurcation. Too much bifurcation points ... ?");
523 stop = true; // this it the end of the tracking
525 // Else this it the end of the tracking
531 //--------------------------------------------------------------------
533 /*TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations,
534 ImagePointer skeleton,
535 ImageIndexType index,
536 ImagePixelType label)
538 DD("TrackFromThisIndex");
541 // Create NeighborhoodIterator
542 typedef itk::NeighborhoodIterator<ImageType> NeighborhoodIteratorType;
543 typename NeighborhoodIteratorType::SizeType radius;
545 NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
548 std::vector<typename NeighborhoodIteratorType::IndexType> listOfTrackedPoint;
551 nit.SetLocation(index);
552 // DD((int)nit.GetCenterPixel());
553 nit.SetCenterPixel(label);
554 listOfTrackedPoint.clear();
555 for(unsigned int i=0; i<nit.Size(); i++) {
556 if (i != nit.GetCenterNeighborhoodIndex ()) { // Do not observe the current point
557 // DD(nit.GetIndex(i));
558 if (nit.GetPixel(i) == GetForegroundValue()) { // if this is foreground, we continue the tracking
559 // DD(nit.GetIndex(i));
560 listOfTrackedPoint.push_back(nit.GetIndex(i));
564 // DD(listOfTrackedPoint.size());
565 if (listOfTrackedPoint.size() == 1) {
566 index = listOfTrackedPoint[0];
569 if (listOfTrackedPoint.size() == 2) {
570 BifurcationType bif(index, label, label+1, label+2);
571 listOfBifurcations.push_back(bif);
572 TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1);
573 TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2);
576 if (listOfTrackedPoint.size() > 2) {
577 std::cerr << "too much bifurcation points ... ?" << std::endl;
580 // Else this it the end of the tracking
587 //--------------------------------------------------------------------
590 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX