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