Medical Imaging Interaction Toolkit  2018.4.99-c4b6bb11
Medical Imaging Interaction Toolkit
itkLineHistogramBasedMassImageFilter.cpp
Go to the documentation of this file.
1 /*============================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center (DKFZ)
6 All rights reserved.
7 
8 Use of this source code is governed by a 3-clause BSD license that can be
9 found in the LICENSE file.
10 
11 ============================================================================*/
12 
13 #ifndef ITKLINEHISTOGRAMBASEDMASSIMAGEFILTER_CPP
14 #define ITKLINEHISTOGRAMBASEDMASSIMAGEFILTER_CPP
15 
17 #include <itkImageRegionIteratorWithIndex.h>
18 #include <itkBinaryContourImageFilter.h>
19 #include <mitkImageCast.h>
20 #include <mitkIOUtil.h>
21 
22 
23 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
26 {
27 
28  if(m_ImageMask.IsNull())
29  {
30  itkExceptionMacro("No binary mask image provided.")
31  }
32 
33  if(m_BinaryContour.IsNull()){
34 
35 
36  typename TMaskImageType::RegionType region = m_ImageMask->GetLargestPossibleRegion();
37  unsigned int xdim = region.GetSize(0);
38  unsigned int ydim = region.GetSize(1);
39  unsigned int zdim = region.GetSize(2);
40 
41  // ensure border pixels are zero so that an closed contour is created
42  typename TMaskImageType::Pointer mask_copy = TMaskImageType::New();
43  mask_copy->SetSpacing(m_ImageMask->GetSpacing());
44  mask_copy->SetDirection(m_ImageMask->GetDirection());
45  mask_copy->SetOrigin(m_ImageMask->GetOrigin());
46  mask_copy->SetRegions(region);
47  mask_copy->Allocate();
48  mask_copy->FillBuffer(0);
49 
50  itk::ImageRegionIteratorWithIndex<TMaskImageType> oit(mask_copy, mask_copy->GetLargestPossibleRegion());
51  itk::ImageRegionConstIteratorWithIndex<TMaskImageType> mit(m_ImageMask, m_ImageMask->GetLargestPossibleRegion());
52  while(!mit.IsAtEnd())
53  {
54  if(mit.Value() != 0 //is under the mask
55  && mit.GetIndex()[0] != 0 //is not at the border
56  && mit.GetIndex()[1] != 0
57  && mit.GetIndex()[2] != 0
58  && mit.GetIndex()[0] != xdim-1
59  && mit.GetIndex()[1] != ydim-1
60  && mit.GetIndex()[2] != zdim-1)
61  {
62  oit.Set(1);
63  }
64  ++mit;
65  ++oit;
66  }
67 
68  typedef itk::BinaryContourImageFilter<TMaskImageType,TMaskImageType> BinaryContourImagefilterType;
69  typename BinaryContourImagefilterType::Pointer filter = BinaryContourImagefilterType::New();
70  filter->SetInput(mask_copy);
71  filter->SetBackgroundValue(0);
72  filter->SetForegroundValue(1);
73  filter->SetFullyConnected(true);
74  filter->Update();
75 
76  m_BinaryContour = filter->GetOutput();
77  mitk::Image::Pointer outimg;
78  mitk::CastToMitkImage(m_BinaryContour,outimg);
79  }
80 
81  if(m_BinaryContour.IsNull())
82  {
83  itkExceptionMacro("No binary contour image provided.")
84  }
85 
86  m_CenterOfMask = GetCenterOfMass(m_ImageMask);
87 
88  if(m_CenterOfMask.is_zero())
89  {
90  itkExceptionMacro("Center of mass is corrupt.")
91  }
92 
93 
94 }
95 
96 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
98 ::GetCenterOfMass(const TMaskImageType * maskImage)
99 {
101  mitk::CastToMitkImage(maskImage, img);
102 
103  typedef itk::ImageRegionConstIterator< TMaskImageType > IteratorType;
104  IteratorType iter( maskImage, maskImage->GetLargestPossibleRegion() );
105  iter.GoToBegin();
106 
107  vnl_vector_fixed<double,3> mean_index;
108  mean_index.fill(0);
109  unsigned int count = 0;
110  while ( !iter.IsAtEnd() )
111  {
112  if ( iter.Get() != 0 )
113  {
114  mitk::Point3D current_index_pos;
115  img->GetGeometry()->IndexToWorld(iter.GetIndex(),current_index_pos);
116  mean_index += current_index_pos.GetVnlVector();
117  count++;
118  }
119  ++iter;
120  }
121 
122  mean_index /= count;
123  return mean_index;
124 }
125 
126 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
128 ::ThreadedGenerateData(const typename Superclass::OutputImageRegionType &outputRegionForThread, ThreadIdType /*threadId*/)
129 {
130 
131  TOutputImageType * outimage = this->GetOutput();
132 
133 
134  itk::ImageRegionConstIteratorWithIndex<TMaskImageType> cit(m_BinaryContour, outputRegionForThread);
135  itk::ImageRegionConstIteratorWithIndex<TInputImageType> iit(this->GetInput(), outputRegionForThread);
136  itk::ImageRegionIteratorWithIndex<TOutputImageType > oit(outimage,outputRegionForThread);
137  typedef typename itk::ImageRegionIteratorWithIndex<TInputImageType >::IndexType IndexType;
138 
139  std::vector<IndexType> target_world_indices;
140 
141  while(!cit.IsAtEnd())
142  {
143  if(cit.Value() != 0 )
144  target_world_indices.push_back(cit.GetIndex());
145  ++cit;
146  }
147 
149  mitk::CastToMitkImage(this->GetInput(), image);
150  mitk::BaseGeometry * transform = image->GetGeometry();
151 
152  while(!target_world_indices.empty())
153  {
154  vnl_vector<double> u_v(3,1), x_v(3,1), p_v(3,1);
155  double line_histo_area = 0;
156 
157  // w->i skull point
158  mitk::Point3D skull_point, tmp_point;
159  image->GetGeometry()->IndexToWorld(target_world_indices.back(),skull_point);
160 
161  // project centerpoint to slice
162  // x = p + s*u
163  u_v = skull_point.GetVnlVector() - m_CenterOfMask;
164  u_v.normalize();
165 
166  p_v = m_CenterOfMask;
167 
168  //step width
169  double s_s = 0.1;
170 
171  // i->w center point
172  mitk::Point3D center_point;
173  center_point[0] = m_CenterOfMask[0];
174  center_point[1] = m_CenterOfMask[1];
175  center_point[2] = m_CenterOfMask[2];
176  IndexType tmp_index;
177  transform->WorldToIndex(center_point,tmp_index);
178  oit.SetIndex(tmp_index);
179 
180  std::vector<IndexType> under_line_indices;
181 
182  while(true)
183  {
184 
185  // store current_index
186  under_line_indices.push_back(tmp_index);
187 
188  // set histo value
189  iit.SetIndex(tmp_index);
190  line_histo_area += iit.Value();
191 
192  // break in end reached
193  if(tmp_index == target_world_indices.back())
194  {
195  break;
196  }
197 
198  // get next voxel index under the line
199  while(tmp_index == oit.GetIndex()) // if new voxel index reached brake
200  {
201  x_v = p_v + s_s * u_v; // next
202  tmp_point.GetVnlVector().set(x_v.data_block());
203  transform->WorldToIndex(tmp_point, tmp_index);
204  s_s+=0.1;
205  }
206 
207  oit.SetIndex(tmp_index);
208  }
209 
210 
211  while (!under_line_indices.empty()) {
212  IndexType current_index = under_line_indices.back();
213  under_line_indices.pop_back();
214 
215  oit.SetIndex(current_index);
216  oit.Set(line_histo_area / (s_s * u_v).magnitude() );
217  }
218 
219  target_world_indices.pop_back();
220 
221  }
222 
223 }
224 
225 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
227 {
228  this->m_ImageMask = maskimage;
229 }
230 
231 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
233 {
234  this->m_BinaryContour = binarycontour;
235 }
236 
237 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
239 {
240  this->SetNumberOfIndexedOutputs(1);
241  this->SetNumberOfIndexedInputs(1);
242 
243  m_CenterOfMask.fill(0);
244 }
245 
246 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
248 {
249 
250 }
251 
252 
253 
254 #endif
mitk::Image::Pointer image
void WorldToIndex(const mitk::Point3D &pt_mm, mitk::Point3D &pt_units) const
Convert world coordinates (in mm) of a point to (continuous!) index coordinates.
void CastToMitkImage(const itk::SmartPointer< ItkOutputImageType > &itkimage, itk::SmartPointer< mitk::Image > &mitkoutputimage)
Cast an itk::Image (with a specific type) to an mitk::Image.
Definition: mitkImageCast.h:74
BaseGeometry Describes the geometry of a data object.