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 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>
67 template<class ArgsInfoType>
69 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
70 SetArgsInfo(ArgsInfoType mArgsInfo)
72 SetVerboseOption_GGO(mArgsInfo);
73 SetVerboseStep_GGO(mArgsInfo);
74 SetWriteStep_GGO(mArgsInfo);
75 SetVerboseWarningOff_GGO(mArgsInfo);
76 SetAFDBFilename_GGO(mArgsInfo);
78 //--------------------------------------------------------------------
81 //--------------------------------------------------------------------
82 template <class ImageType>
84 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
85 GenerateOutputInformation()
87 Superclass::GenerateOutputInformation();
88 //this->GetOutput(0)->SetRequestedRegion(this->GetOutput(0)->GetLargestPossibleRegion());
91 input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
96 } catch (clitk::ExceptionObject e) {
97 // Do nothing if not found, it will be used anyway to write the result
102 //--------------------------------------------------------------------
105 //--------------------------------------------------------------------
106 template <class ImageType>
108 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
111 StartNewStep("Thinning filter (skeleton)");
113 typedef itk::BinaryThinningImageFilter3D<ImageType, ImageType> ThinningFilterType;
114 typename ThinningFilterType::Pointer thinningFilter = ThinningFilterType::New();
115 thinningFilter->SetInput(input); // input = trachea
116 thinningFilter->Update();
117 skeleton = thinningFilter->GetOutput();
118 StopCurrentStep<ImageType>(skeleton);
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();
134 // Initialize neighborhooditerator
135 typedef itk::NeighborhoodIterator<ImageType> NeighborhoodIteratorType;
136 typename NeighborhoodIteratorType::SizeType radius;
138 NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
140 // Find first label number (must be different from BG and FG)
141 typename ImageType::PixelType label = GetForegroundValue()+1;
142 while ((label == GetBackgroundValue()) || (label == GetForegroundValue())) { label++; }
147 Goal : find position of C, RUL, RML, RLL, LUL, LLL bronchus
148 Carina : ok "easy", track, slice by slice until 2 path into different label
155 // Track from the first point
156 StartNewStep("Start tracking");
157 std::vector<BifurcationType> listOfBifurcations;
158 m_SkeletonTree.set_head(index);
159 TrackFromThisIndex(listOfBifurcations, skeleton, index, label, m_SkeletonTree.begin());
162 // Convert index into physical point coordinates
163 for(unsigned int i=0; i<listOfBifurcations.size(); i++) {
164 skeleton->TransformIndexToPhysicalPoint(listOfBifurcations[i].index,
165 listOfBifurcations[i].point);
168 // Search for the first slice that separate the bronchus
169 // (carina). Labelize slice by slice, stop when the two points of
170 // the skeleton ar not in the same connected component
171 StartNewStep("Search for carina position");
172 typedef clitk::ExtractSliceFilter<ImageType> ExtractSliceFilterType;
173 typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
174 extractSliceFilter->SetInput(input);
175 extractSliceFilter->SetDirection(2);
176 extractSliceFilter->Update();
177 typedef typename ExtractSliceFilterType::SliceType SliceType;
178 std::vector<typename SliceType::Pointer> mInputSlices;
179 extractSliceFilter->GetOutputSlices(mInputSlices);
180 DD(mInputSlices.size());
183 int slice_index = listOfBifurcations[0].index[2]; // first slice from carina in skeleton
185 TreeIterator firstIter = m_SkeletonTree.child(listOfBifurcations[0].treeIter, 0);
186 TreeIterator secondIter = m_SkeletonTree.child(listOfBifurcations[0].treeIter, 1);
187 typename SliceType::IndexType in1;
188 typename SliceType::IndexType in2;
190 // Labelize the current slice
191 typename SliceType::Pointer temp = Labelize<SliceType>(mInputSlices[slice_index],
192 GetBackgroundValue(),
194 0); // min component size=0
195 // Check the value of the two skeleton points;
196 in1[0] = (*firstIter)[0];
197 in1[1] = (*firstIter)[1];
198 typename SliceType::PixelType v1 = temp->GetPixel(in1);
199 in2[0] = (*secondIter)[0];
200 in2[1] = (*secondIter)[1];
201 typename SliceType::PixelType v2 = temp->GetPixel(in2);
203 // Check the label value of the two points
206 stop = true; // We found it !
210 if (slice_index == (int)(mInputSlices.size()-1)) {
211 clitkExceptionMacro("Error while searching for carina, the two skeleton points are always in the same CC ... ???");
220 ImageIndexType carina_index; // middle position in X/Y
221 carina_index[0] = lrint(in2[0] + in1[0])/2.0;
222 carina_index[1] = lrint(in2[1] + in1[1])/2.0;
223 carina_index[2] = slice_index;
224 // Get physical coordinates
225 ImagePointType carina_position;
226 skeleton->TransformIndexToPhysicalPoint(carina_index,
229 // Set and save Carina position
230 if (GetVerboseStep()) {
231 std::cout << "\t Found carina at " << carina_position << " mm" << std::endl;
233 GetAFDB()->SetPoint3D("Carina", carina_position);
235 // Write bifurcation (debug)
236 for(uint i=0; i<listOfBifurcations.size(); i++) {
237 GetAFDB()->SetPoint3D("Bif"+toString(i), listOfBifurcations[i].point);
240 // Set the output (skeleton);
241 this->GraftOutput(skeleton); // not SetNthOutput
243 //--------------------------------------------------------------------
246 //--------------------------------------------------------------------
247 template <class ImageType>
249 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
250 TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations,
251 ImagePointer skeleton,
252 ImageIndexType index,
253 ImagePixelType label,
254 TreeIterator currentNode)
256 // Create NeighborhoodIterator
257 typedef itk::NeighborhoodIterator<ImageType> NeighborhoodIteratorType;
258 typename NeighborhoodIteratorType::SizeType radius;
260 NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
263 std::vector<typename NeighborhoodIteratorType::IndexType> listOfTrackedPoint;
266 nit.SetLocation(index);
267 nit.SetCenterPixel(label);
268 listOfTrackedPoint.clear();
269 for(unsigned int i=0; i<nit.Size(); i++) {
270 if (i != nit.GetCenterNeighborhoodIndex ()) { // Do not observe the current point
271 if (nit.GetPixel(i) == GetForegroundValue()) { // if this is foreground, we continue the tracking
272 listOfTrackedPoint.push_back(nit.GetIndex(i));
276 if (listOfTrackedPoint.size() == 1) {
277 // Add this point to the current path
278 currentNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[0]);
279 index = listOfTrackedPoint[0];
280 skeleton->SetPixel(index, label); // change label in skeleton image
283 if (listOfTrackedPoint.size() == 2) {
284 // m_SkeletonTree->Add(listOfTrackedPoint[0], index); // the parent is 'index'
285 // m_SkeletonTree->Add(listOfTrackedPoint[1], index); // the parent is 'index'
286 BifurcationType bif(index, label, label+1, label+2);
287 bif.treeIter = currentNode;
288 listOfBifurcations.push_back(bif);
289 TreeIterator firstNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[0]);
290 TreeIterator secondNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[1]);
291 TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1, firstNode);
292 TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2, secondNode);
295 //DD(listOfTrackedPoint.size());
296 if (listOfTrackedPoint.size() > 2) {
297 //clitkExceptionMacro("error while tracking trachea bifurcation. Too much bifurcation points ... ?");
298 stop = true; // this it the end of the tracking
300 // Else this it the end of the tracking
306 //--------------------------------------------------------------------
308 /*TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations,
309 ImagePointer skeleton,
310 ImageIndexType index,
311 ImagePixelType label)
313 DD("TrackFromThisIndex");
316 // Create NeighborhoodIterator
317 typedef itk::NeighborhoodIterator<ImageType> NeighborhoodIteratorType;
318 typename NeighborhoodIteratorType::SizeType radius;
320 NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
323 std::vector<typename NeighborhoodIteratorType::IndexType> listOfTrackedPoint;
326 nit.SetLocation(index);
327 // DD((int)nit.GetCenterPixel());
328 nit.SetCenterPixel(label);
329 listOfTrackedPoint.clear();
330 for(unsigned int i=0; i<nit.Size(); i++) {
331 if (i != nit.GetCenterNeighborhoodIndex ()) { // Do not observe the current point
332 // DD(nit.GetIndex(i));
333 if (nit.GetPixel(i) == GetForegroundValue()) { // if this is foreground, we continue the tracking
334 // DD(nit.GetIndex(i));
335 listOfTrackedPoint.push_back(nit.GetIndex(i));
339 // DD(listOfTrackedPoint.size());
340 if (listOfTrackedPoint.size() == 1) {
341 index = listOfTrackedPoint[0];
344 if (listOfTrackedPoint.size() == 2) {
345 BifurcationType bif(index, label, label+1, label+2);
346 listOfBifurcations.push_back(bif);
347 TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1);
348 TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2);
351 if (listOfTrackedPoint.size() > 2) {
352 std::cerr << "too much bifurcation points ... ?" << std::endl;
355 // Else this it the end of the tracking
362 //--------------------------------------------------------------------
365 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX