Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
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.