Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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.