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
mitkNrrdTensorImageReader.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 #include "mitkNrrdTensorImageReader.h"
18 #include <mitkCustomMimeType.h>
19 #include "mitkDiffusionCoreIOMimeTypes.h"
20 
21 #include "itkImageFileReader.h"
22 #include "itkImageRegionIterator.h"
23 #include "itkMetaDataObject.h"
24 #include "itkNrrdImageIO.h"
25 #include <itkNiftiImageIO.h>
26 
27 #include "mitkITKImageImport.h"
28 #include "mitkImageDataItem.h"
29 #include <mitkLocaleSwitch.h>
30 
31 namespace mitk
32 {
33 
34  NrrdTensorImageReader::NrrdTensorImageReader(const NrrdTensorImageReader& other)
35  : mitk::AbstractFileReader(other)
36  {
37  }
38 
39  NrrdTensorImageReader::NrrdTensorImageReader()
40  : mitk::AbstractFileReader( CustomMimeType( mitk::DiffusionCoreIOMimeTypes::DTI_MIMETYPE() ), mitk::DiffusionCoreIOMimeTypes::DTI_MIMETYPE_DESCRIPTION() )
41  {
42  m_ServiceReg = this->RegisterService();
43  }
44 
45  NrrdTensorImageReader::~NrrdTensorImageReader()
46  {
47  }
48 
49  std::vector<itk::SmartPointer<BaseData> > NrrdTensorImageReader::Read()
50  {
51  std::vector<itk::SmartPointer<mitk::BaseData> > result;
52  std::string location = GetInputLocation();
53 
54  if ( location == "")
55  {
56  throw itk::ImageFileReaderException(__FILE__, __LINE__, "Sorry, the filename is empty!");
57  }
58  else
59  {
60  try
61  {
62  mitk::LocaleSwitch localeSwitch("C");
63 
64  try
65  {
66  std::string fname3 = "temp_dti.nii";
67  itksys::SystemTools::CopyAFile(location.c_str(), fname3.c_str());
68 
69  typedef itk::VectorImage<float,3> ImageType;
70  itk::NiftiImageIO::Pointer io = itk::NiftiImageIO::New();
71  typedef itk::ImageFileReader<ImageType> FileReaderType;
72  FileReaderType::Pointer reader = FileReaderType::New();
73  reader->SetImageIO(io);
74  reader->SetFileName(fname3);
75  reader->Update();
76  ImageType::Pointer img = reader->GetOutput();
77 
78  typedef itk::Image<itk::DiffusionTensor3D<float>,3> VecImgType;
79  VecImgType::Pointer vecImg = VecImgType::New();
80  vecImg->SetSpacing( img->GetSpacing() ); // Set the image spacing
81  vecImg->SetOrigin( img->GetOrigin() ); // Set the image origin
82  vecImg->SetDirection( img->GetDirection() ); // Set the image direction
83  vecImg->SetRegions( img->GetLargestPossibleRegion());
84  vecImg->Allocate();
85 
86  itk::ImageRegionIterator<VecImgType> ot (vecImg, vecImg->GetLargestPossibleRegion() );
87  ot.GoToBegin();
88 
89  itk::ImageRegionIterator<ImageType> it (img, img->GetLargestPossibleRegion() );
90  it.GoToBegin();
91 
92  typedef ImageType::PixelType VarPixType;
93  typedef VecImgType::PixelType FixPixType;
94  int numComponents = img->GetNumberOfComponentsPerPixel();
95 
96  if (numComponents==6)
97  {
98  MITK_INFO << "Trying to load dti as 6-comp nifti ...";
99  while (!it.IsAtEnd())
100  {
101  VarPixType vec = it.Get();
102  FixPixType fixVec(vec.GetDataPointer());
103 
104  itk::DiffusionTensor3D<float> tensor;
105  tensor.SetElement(0, vec.GetElement(0));
106  tensor.SetElement(1, vec.GetElement(1));
107  tensor.SetElement(2, vec.GetElement(2));
108  tensor.SetElement(3, vec.GetElement(3));
109  tensor.SetElement(4, vec.GetElement(4));
110  tensor.SetElement(5, vec.GetElement(5));
111 
112  fixVec = tensor;
113 
114  ot.Set(fixVec);
115  ++ot;
116  ++it;
117  }
118  }
119  else if(numComponents==9)
120  {
121  MITK_INFO << "Trying to load dti as 9-comp nifti ...";
122  while (!it.IsAtEnd())
123  {
124  VarPixType vec = it.Get();
125  itk::DiffusionTensor3D<float> tensor;
126  tensor.SetElement(0, vec.GetElement(0));
127  tensor.SetElement(1, vec.GetElement(1));
128  tensor.SetElement(2, vec.GetElement(2));
129  tensor.SetElement(3, vec.GetElement(4));
130  tensor.SetElement(4, vec.GetElement(5));
131  tensor.SetElement(5, vec.GetElement(8));
132 
133  FixPixType fixVec(tensor);
134  ot.Set(fixVec);
135  ++ot;
136  ++it;
137  }
138  }
139  else if (numComponents==1)
140  {
141  MITK_INFO << "Trying to load dti as 4D nifti ...";
142  typedef itk::Image<float,4> ImageType;
143  typedef itk::ImageFileReader<ImageType> FileReaderType;
144  FileReaderType::Pointer reader = FileReaderType::New();
145  reader->SetImageIO(io);
146  reader->SetFileName(fname3);
147  reader->Update();
148  ImageType::Pointer img = reader->GetOutput();
149 
150  itk::Size<4> size = img->GetLargestPossibleRegion().GetSize();
151 
152  while (!ot.IsAtEnd())
153  {
154  itk::DiffusionTensor3D<float> tensor;
155  ImageType::IndexType idx;
156  idx[0] = ot.GetIndex()[0]; idx[1] = ot.GetIndex()[1]; idx[2] = ot.GetIndex()[2];
157 
158  if (size[3]==6)
159  {
160  for (unsigned int te=0; te<size[3]; te++)
161  {
162  idx[3] = te;
163  tensor.SetElement(te, img->GetPixel(idx));
164  }
165  }
166  else if (size[3]==9)
167  {
168  idx[3] = 0;
169  tensor.SetElement(0, img->GetPixel(idx));
170  idx[3] = 1;
171  tensor.SetElement(1, img->GetPixel(idx));
172  idx[3] = 2;
173  tensor.SetElement(2, img->GetPixel(idx));
174  idx[3] = 4;
175  tensor.SetElement(3, img->GetPixel(idx));
176  idx[3] = 5;
177  tensor.SetElement(4, img->GetPixel(idx));
178  idx[3] = 8;
179  tensor.SetElement(5, img->GetPixel(idx));
180  }
181  else
182  throw itk::ImageFileReaderException(__FILE__, __LINE__, "Unknown number of components for DTI file. Should be 6 or 9!");
183 
184  FixPixType fixVec(tensor);
185  ot.Set(fixVec);
186  ++ot;
187  }
188  }
189  OutputType::Pointer resultImage = OutputType::New();
190  resultImage->InitializeByItk( vecImg.GetPointer() );
191  resultImage->SetVolume( vecImg->GetBufferPointer() );
192  result.push_back( resultImage.GetPointer() );
193  }
194  catch(...)
195  {
196  MITK_INFO << "Trying to load dti as nrrd ...";
197 
198  typedef itk::VectorImage<float,3> ImageType;
199  itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New();
200  typedef itk::ImageFileReader<ImageType> FileReaderType;
201  FileReaderType::Pointer reader = FileReaderType::New();
202  reader->SetImageIO(io);
203  reader->SetFileName(location);
204  reader->Update();
205  ImageType::Pointer img = reader->GetOutput();
206 
207  typedef itk::Image<itk::DiffusionTensor3D<float>,3> VecImgType;
208  VecImgType::Pointer vecImg = VecImgType::New();
209  vecImg->SetSpacing( img->GetSpacing() ); // Set the image spacing
210  vecImg->SetOrigin( img->GetOrigin() ); // Set the image origin
211  vecImg->SetDirection( img->GetDirection() ); // Set the image direction
212  vecImg->SetRegions( img->GetLargestPossibleRegion());
213  vecImg->Allocate();
214 
215  itk::ImageRegionIterator<VecImgType> ot (vecImg, vecImg->GetLargestPossibleRegion() );
216  ot.GoToBegin();
217 
218  itk::ImageRegionIterator<ImageType> it (img, img->GetLargestPossibleRegion() );
219  it.GoToBegin();
220 
221  typedef ImageType::PixelType VarPixType;
222  typedef VecImgType::PixelType FixPixType;
223  int numComponents = img->GetNumberOfComponentsPerPixel();
224 
225  itk::MetaDataDictionary imgMetaDictionary = img->GetMetaDataDictionary();
226  std::vector<std::string> imgMetaKeys = imgMetaDictionary.GetKeys();
227  std::vector<std::string>::const_iterator itKey = imgMetaKeys.begin();
228  std::string metaString;
229 
230  bool readFrame = false;
231  double xx, xy, xz, yx, yy, yz, zx, zy, zz;
232  MeasurementFrameType measFrame;
233  measFrame.SetIdentity();
234  MeasurementFrameType measFrameTransp;
235  measFrameTransp.SetIdentity();
236 
237  for (; itKey != imgMetaKeys.end(); itKey ++)
238  {
239  itk::ExposeMetaData<std::string> (imgMetaDictionary, *itKey, metaString);
240  if (itKey->find("measurement frame") != std::string::npos)
241  {
242  sscanf(metaString.c_str(), " ( %lf , %lf , %lf ) ( %lf , %lf , %lf ) ( %lf , %lf , %lf ) \n", &xx, &xy, &xz, &yx, &yy, &yz, &zx, &zy, &zz);
243 
244  if (xx>10e-10 || xy>10e-10 || xz>10e-10 ||
245  yx>10e-10 || yy>10e-10 || yz>10e-10 ||
246  zx>10e-10 || zy>10e-10 || zz>10e-10 )
247  {
248  readFrame = true;
249 
250  measFrame(0,0) = xx;
251  measFrame(0,1) = xy;
252  measFrame(0,2) = xz;
253  measFrame(1,0) = yx;
254  measFrame(1,1) = yy;
255  measFrame(1,2) = yz;
256  measFrame(2,0) = zx;
257  measFrame(2,1) = zy;
258  measFrame(2,2) = zz;
259 
260  measFrameTransp = measFrame.GetTranspose();
261  }
262  }
263  }
264 
265  if (numComponents==6)
266  {
267  while (!it.IsAtEnd())
268  {
269  // T'=RTR'
270  VarPixType vec = it.Get();
271  FixPixType fixVec(vec.GetDataPointer());
272 
273  if(readFrame)
274  {
275  itk::DiffusionTensor3D<float> tensor;
276  tensor.SetElement(0, vec.GetElement(0));
277  tensor.SetElement(1, vec.GetElement(1));
278  tensor.SetElement(2, vec.GetElement(2));
279  tensor.SetElement(3, vec.GetElement(3));
280  tensor.SetElement(4, vec.GetElement(4));
281  tensor.SetElement(5, vec.GetElement(5));
282 
283  tensor = ConvertMatrixTypeToFixedArrayType(tensor.PreMultiply(measFrame));
284  tensor = ConvertMatrixTypeToFixedArrayType(tensor.PostMultiply(measFrameTransp));
285  fixVec = tensor;
286  }
287 
288  ot.Set(fixVec);
289  ++ot;
290  ++it;
291  }
292  }
293  else if(numComponents==9)
294  {
295  while (!it.IsAtEnd())
296  {
297  VarPixType vec = it.Get();
298  itk::DiffusionTensor3D<float> tensor;
299  tensor.SetElement(0, vec.GetElement(0));
300  tensor.SetElement(1, vec.GetElement(1));
301  tensor.SetElement(2, vec.GetElement(2));
302  tensor.SetElement(3, vec.GetElement(4));
303  tensor.SetElement(4, vec.GetElement(5));
304  tensor.SetElement(5, vec.GetElement(8));
305 
306  if(readFrame)
307  {
308  tensor = ConvertMatrixTypeToFixedArrayType(tensor.PreMultiply(measFrame));
309  tensor = ConvertMatrixTypeToFixedArrayType(tensor.PostMultiply(measFrameTransp));
310  }
311 
312  FixPixType fixVec(tensor);
313  ot.Set(fixVec);
314  ++ot;
315  ++it;
316  }
317  }
318  else if (numComponents==1)
319  {
320  typedef itk::Image<float,4> ImageType;
321  itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New();
322  typedef itk::ImageFileReader<ImageType> FileReaderType;
323  FileReaderType::Pointer reader = FileReaderType::New();
324  reader->SetImageIO(io);
325  reader->SetFileName(location);
326  reader->Update();
327  ImageType::Pointer img = reader->GetOutput();
328 
329  itk::Size<4> size = img->GetLargestPossibleRegion().GetSize();
330 
331  while (!ot.IsAtEnd())
332  {
333  itk::DiffusionTensor3D<float> tensor;
334  ImageType::IndexType idx;
335  idx[0] = ot.GetIndex()[0]; idx[1] = ot.GetIndex()[1]; idx[2] = ot.GetIndex()[2];
336 
337  if (size[3]==6)
338  {
339  for (unsigned int te=0; te<size[3]; te++)
340  {
341  idx[3] = te;
342  tensor.SetElement(te, img->GetPixel(idx));
343  }
344  }
345  else if (size[3]==9)
346  {
347  idx[3] = 0;
348  tensor.SetElement(0, img->GetPixel(idx));
349  idx[3] = 1;
350  tensor.SetElement(1, img->GetPixel(idx));
351  idx[3] = 2;
352  tensor.SetElement(2, img->GetPixel(idx));
353  idx[3] = 4;
354  tensor.SetElement(3, img->GetPixel(idx));
355  idx[3] = 5;
356  tensor.SetElement(4, img->GetPixel(idx));
357  idx[3] = 8;
358  tensor.SetElement(5, img->GetPixel(idx));
359  }
360  else
361  throw itk::ImageFileReaderException(__FILE__, __LINE__, "Unknown number of komponents for DTI file. Should be 6 or 9!");
362 
363  if(readFrame)
364  {
365  tensor = ConvertMatrixTypeToFixedArrayType(tensor.PreMultiply(measFrame));
366  tensor = ConvertMatrixTypeToFixedArrayType(tensor.PostMultiply(measFrameTransp));
367  }
368  FixPixType fixVec(tensor);
369  ot.Set(fixVec);
370  ++ot;
371  }
372  }
373  else
374  {
375  throw itk::ImageFileReaderException(__FILE__, __LINE__, "Image has wrong number of pixel components!");
376  }
377 
378  OutputType::Pointer resultImage = OutputType::New();
379  resultImage->InitializeByItk( vecImg.GetPointer() );
380  resultImage->SetVolume( vecImg->GetBufferPointer() );
381  result.push_back( resultImage.GetPointer() );
382  }
383 
384  }
385  catch(std::exception& e)
386  {
387  throw itk::ImageFileReaderException(__FILE__, __LINE__, e.what());
388  }
389  catch(...)
390  {
391  throw itk::ImageFileReaderException(__FILE__, __LINE__, "Sorry, an error occurred while reading the requested DTI file!");
392  }
393  }
394 
395  return result;
396  }
397 
398 itk::DiffusionTensor3D<float> NrrdTensorImageReader
399 ::ConvertMatrixTypeToFixedArrayType(const itk::DiffusionTensor3D<float>::Superclass::MatrixType & matrix)
400 {
401  /* | 0 1 2 |
402  * | X 3 4 |
403  * | X X 5 |
404  */
405  itk::DiffusionTensor3D<float> arr;
406  arr.SetElement(0,matrix(0,0));
407  arr.SetElement(1,matrix(0,1));
408  arr.SetElement(2,matrix(0,2));
409  arr.SetElement(3,matrix(1,3));
410  arr.SetElement(4,matrix(1,4));
411  arr.SetElement(5,matrix(2,5));
412  return arr;
413 }
414 
415 } //namespace MITK
416 
417 mitk::NrrdTensorImageReader* mitk::NrrdTensorImageReader::Clone() const
418 {
419  return new NrrdTensorImageReader(*this);
420 }
mitk::NrrdTensorImageReader::Read
virtual std::vector< itk::SmartPointer< BaseData > > Read() override
Reads a path or stream and creates a list of BaseData objects.
Definition: mitkNrrdTensorImageReader.cpp:49
mitk::Pointer
itk::SmartPointer< Self > Pointer
Definition: mitkRenderingManager.h:389
MITK_INFO
#define MITK_INFO
Definition: mitkLogMacros.h:22
mitk
DataCollection - Class to facilitate loading/accessing structured data.
Definition: GeometryOverview.dox:1
mitk::TensorImage::New
static Pointer New()
mitk::CustomMimeType
The CustomMimeType class represents a custom mime-type which may be registered as a service object...
Definition: mitkCustomMimeType.h:48
mitk::ImageType
map::core::discrete::Elements< 3 >::InternalImageType ImageType
Definition: mitkMultiModalAffine_default.cpp:26
mitk::AbstractFileReader::RegisterService
us::ServiceRegistration< IFileReader > RegisterService(us::ModuleContext *context=us::GetModuleContext())
Definition: mitkAbstractFileReader.cpp:139
mitk::LocaleSwitch
Convenience class to temporarily change the current locale.
Definition: mitkLocaleSwitch.h:50
mitk::NrrdTensorImageReader::MeasurementFrameType
itk::Matrix< float, 3, 3 > MeasurementFrameType
Definition: mitkNrrdTensorImageReader.h:41
mitk::AbstractFileReader
Base class for creating mitk::BaseData objects from files or streams.
Definition: mitkAbstractFileReader.h:46
PixelType
unsigned short PixelType
Definition: TumorProgressionMapping.cpp:55
mitk::AbstractFileReader::GetInputLocation
virtual std::string GetInputLocation() const override
Get the current input location.
Definition: mitkAbstractFileReader.cpp:264
mitk::New
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.