Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
mitkDiffusionImageNiftiReaderService.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 __mitkDiffusionImageNiftiReaderService_cpp
18 #define __mitkDiffusionImageNiftiReaderService_cpp
19 
21 
22 #include <iostream>
23 #include <fstream>
24 
25 
26 // Diffusion properties
28 #include <mitkBValueMapProperty.h>
30 #include <mitkProperties.h>
31 
32 // ITK includes
33 #include <itkImageRegionIterator.h>
34 #include <itkMetaDataObject.h>
35 #include "itksys/SystemTools.hxx"
36 #include "itkImageFileReader.h"
37 #include "itkMetaDataObject.h"
38 #include "itkNiftiImageIO.h"
39 
40 #include "mitkCustomMimeType.h"
42 
43 #include <mitkITKImageImport.h>
44 #include <mitkImageWriteAccessor.h>
45 #include <mitkImageDataItem.h>
46 #include <mitkLocaleSwitch.h>
47 #include "mitkIOUtil.h"
48 
49 
50 namespace mitk
51 {
52 
55  : AbstractFileReader(other)
56 {
57 }
58 
59 
60 DiffusionImageNiftiReaderService* DiffusionImageNiftiReaderService::Clone() const
61 {
62  return new DiffusionImageNiftiReaderService(*this);
63 }
64 
65 
68 {}
69 
72  : mitk::AbstractFileReader( CustomMimeType( mitk::DiffusionCoreIOMimeTypes::DWI_NIFTI_MIMETYPE() ), mitk::DiffusionCoreIOMimeTypes::DWI_NIFTI_MIMETYPE_DESCRIPTION() )
73 {
74  m_ServiceReg = this->RegisterService();
75 }
76 
77 std::vector<itk::SmartPointer<mitk::BaseData> >
80 {
81  std::vector<itk::SmartPointer<mitk::BaseData> > result;
82 
83  // Since everything is completely read in GenerateOutputInformation() it is stored
84  // in a cache variable. A timestamp is associated.
85  // If the timestamp of the cache variable is newer than the MTime, we only need to
86  // assign the cache variable to the DataObject.
87  // Otherwise, the tree must be read again from the file and OuputInformation must
88  // be updated!
89 
90  if(m_OutputCache.IsNull()) InternalRead();
91 
92  result.push_back(m_OutputCache.GetPointer());
93  return result;
94 }
95 
96 
98 {
99  OutputType::Pointer outputForCache = OutputType::New();
100  if ( this->GetInputLocation() == "")
101  {
102  throw itk::ImageFileReaderException(__FILE__, __LINE__, "Sorry, the filename to be read is empty!");
103  }
104  else
105  {
106  try
107  {
108  mitk::LocaleSwitch localeSwitch("C");
109 
110  MITK_INFO << "DiffusionImageNiftiReaderService: reading image information";
111  VectorImageType::Pointer itkVectorImage;
112 
113  std::string ext = this->GetMimeType()->GetExtension( this->GetInputLocation() );
114  ext = itksys::SystemTools::LowerCase( ext );
115 
116  if(ext == ".fsl" || ext == ".fslgz")
117  {
118  // create temporary file with correct ending for nifti-io
119  std::string fname3 = "temp_dwi";
120  fname3 += ext == ".fsl" ? ".nii" : ".nii.gz";
121  itksys::SystemTools::CopyAFile(this->GetInputLocation().c_str(), fname3.c_str());
122 
123  // create reader and read file
124  typedef itk::Image<DiffusionPixelType,4> ImageType4D;
126  typedef itk::ImageFileReader<ImageType4D> FileReaderType;
128  reader->SetFileName(fname3);
129  reader->SetImageIO(io2);
130  reader->Update();
131  ImageType4D::Pointer img4 = reader->GetOutput();
132 
133  // delete temporary file
134  itksys::SystemTools::RemoveFile(fname3.c_str());
135 
136  // convert 4D file to vector image
137  itkVectorImage = VectorImageType::New();
138 
139  VectorImageType::SpacingType spacing;
140  ImageType4D::SpacingType spacing4 = img4->GetSpacing();
141  for(int i=0; i<3; i++)
142  spacing[i] = spacing4[i];
143  itkVectorImage->SetSpacing( spacing ); // Set the image spacing
144 
146  ImageType4D::PointType origin4 = img4->GetOrigin();
147  for(int i=0; i<3; i++)
148  origin[i] = origin4[i];
149  itkVectorImage->SetOrigin( origin ); // Set the image origin
150 
151  VectorImageType::DirectionType direction;
152  ImageType4D::DirectionType direction4 = img4->GetDirection();
153  for(int i=0; i<3; i++)
154  for(int j=0; j<3; j++)
155  direction[i][j] = direction4[i][j];
156  itkVectorImage->SetDirection( direction ); // Set the image direction
157 
158  VectorImageType::RegionType region;
159  ImageType4D::RegionType region4 = img4->GetLargestPossibleRegion();
160 
161  VectorImageType::RegionType::SizeType size;
162  ImageType4D::RegionType::SizeType size4 = region4.GetSize();
163 
164  for(int i=0; i<3; i++)
165  size[i] = size4[i];
166 
167  VectorImageType::RegionType::IndexType index;
168  ImageType4D::RegionType::IndexType index4 = region4.GetIndex();
169  for(int i=0; i<3; i++)
170  index[i] = index4[i];
171 
172  region.SetSize(size);
173  region.SetIndex(index);
174  itkVectorImage->SetRegions( region );
175 
176  itkVectorImage->SetVectorLength(size4[3]);
177  itkVectorImage->Allocate();
178 
179  itk::ImageRegionIterator<VectorImageType> it ( itkVectorImage, itkVectorImage->GetLargestPossibleRegion() );
180  typedef VectorImageType::PixelType VecPixType;
181  for (it.GoToBegin(); !it.IsAtEnd(); ++it)
182  {
183  VecPixType vec = it.Get();
184  VectorImageType::IndexType currentIndex = it.GetIndex();
185  for(int i=0; i<3; i++)
186  index4[i] = currentIndex[i];
187  for(unsigned int ind=0; ind<vec.Size(); ind++)
188  {
189  index4[3] = ind;
190  vec[ind] = img4->GetPixel(index4);
191  }
192  it.Set(vec);
193  }
194  }
195  else if(ext == ".nii" || ext == ".nii.gz")
196  {
197  // create reader and read file
198  typedef itk::Image<DiffusionPixelType,4> ImageType4D;
200  typedef itk::ImageFileReader<ImageType4D> FileReaderType;
202  reader->SetFileName( this->GetInputLocation() );
203  reader->SetImageIO(io2);
204  reader->Update();
205  ImageType4D::Pointer img4 = reader->GetOutput();
206 
207 
208  // convert 4D file to vector image
209  itkVectorImage = VectorImageType::New();
210 
211  VectorImageType::SpacingType spacing;
212  ImageType4D::SpacingType spacing4 = img4->GetSpacing();
213  for(int i=0; i<3; i++)
214  spacing[i] = spacing4[i];
215  itkVectorImage->SetSpacing( spacing ); // Set the image spacing
216 
218  ImageType4D::PointType origin4 = img4->GetOrigin();
219  for(int i=0; i<3; i++)
220  origin[i] = origin4[i];
221  itkVectorImage->SetOrigin( origin ); // Set the image origin
222 
223  VectorImageType::DirectionType direction;
224  ImageType4D::DirectionType direction4 = img4->GetDirection();
225  for(int i=0; i<3; i++)
226  for(int j=0; j<3; j++)
227  direction[i][j] = direction4[i][j];
228  itkVectorImage->SetDirection( direction ); // Set the image direction
229 
230  VectorImageType::RegionType region;
231  ImageType4D::RegionType region4 = img4->GetLargestPossibleRegion();
232 
233  VectorImageType::RegionType::SizeType size;
234  ImageType4D::RegionType::SizeType size4 = region4.GetSize();
235 
236  for(int i=0; i<3; i++)
237  size[i] = size4[i];
238 
239  VectorImageType::RegionType::IndexType index;
240  ImageType4D::RegionType::IndexType index4 = region4.GetIndex();
241  for(int i=0; i<3; i++)
242  index[i] = index4[i];
243 
244  region.SetSize(size);
245  region.SetIndex(index);
246  itkVectorImage->SetRegions( region );
247 
248  itkVectorImage->SetVectorLength(size4[3]);
249  itkVectorImage->Allocate();
250 
251  itk::ImageRegionIterator<VectorImageType> it ( itkVectorImage, itkVectorImage->GetLargestPossibleRegion() );
252  typedef VectorImageType::PixelType VecPixType;
253  for (it.GoToBegin(); !it.IsAtEnd(); ++it)
254  {
255  VecPixType vec = it.Get();
256  VectorImageType::IndexType currentIndex = it.GetIndex();
257  for(int i=0; i<3; i++)
258  index4[i] = currentIndex[i];
259  for(unsigned int ind=0; ind<vec.Size(); ind++)
260  {
261  index4[3] = ind;
262  vec[ind] = img4->GetPixel(index4);
263  }
264  it.Set(vec);
265  }
266  }
267 
268  // Diffusion Image information START
271  MeasurementFrameType MeasurementFrame;
272  float BValue = -1;
273  // Diffusion Image information END
274 
275  if(ext == ".fsl" || ext == ".fslgz" || ext == ".nii" || ext == ".nii.gz")
276  {
277  // some parsing depending on the extension
278  bool useFSLstyle( true );
279  std::string bvecsExtension("");
280  std::string bvalsExtension("");
281 
282  std::string base = itksys::SystemTools::GetFilenamePath( this->GetInputLocation() ) + "/"
284 
285  // check for possible file names
286  {
287  if( useFSLstyle && itksys::SystemTools::FileExists( std::string( base + ".bvec").c_str() )
288  && itksys::SystemTools::FileExists( std::string( base + ".bval").c_str() )
289  )
290  {
291  useFSLstyle = false;
292  bvecsExtension = ".bvec";
293  bvalsExtension = ".bval";
294  }
295 
296  if( useFSLstyle && itksys::SystemTools::FileExists( std::string( base + ".bvecs").c_str() )
297  && itksys::SystemTools::FileExists( std::string( base + ".bvals").c_str() )
298  )
299  {
300  useFSLstyle = false;
301  bvecsExtension = ".bvecs";
302  bvalsExtension = ".bvals";
303  }
304  }
305 
306  std::string line;
307  std::vector<float> bvec_entries;
308  std::string fname = this->GetInputLocation();
309  if( useFSLstyle )
310  {
311  fname += ".bvecs";
312  }
313  else
314  {
315  fname = std::string( base + bvecsExtension);
316  }
317  std::ifstream myfile (fname.c_str());
318  if (myfile.is_open())
319  {
320  while ( myfile.good() )
321  {
322  std::string line2;
323  getline (myfile,line2);
324 
325  std::stringstream iss;
326  iss << line2;
327 
328  while(getline(iss,line, ' '))
329  {
330  // remove any potenial control sequences that might be introduced by lines ending in a single space
331  line.erase(std::remove_if(line.begin(), line.end(),
332  [](char c) { return std::isspace(c) || std::iscntrl(c); } ), line.end());
333 
334  if (line.length() > 0 ) // otherwise string contained only control characters before, empty string are converted to 0 by atof resulting in a broken b-value list
335  {
336  bvec_entries.push_back(atof(line.c_str()));
337  }
338  }
339  }
340  myfile.close();
341  }
342  else
343  {
344  MITK_INFO << "Unable to open bvecs file. Expected name: " << fname;
345  }
346 
347  std::vector<float> bval_entries;
348  std::string fname2 = this->GetInputLocation();
349  if( useFSLstyle )
350  {
351  fname2 += ".bvals";
352  }
353  else
354  {
355  fname2 = std::string( base + bvalsExtension);
356  }
357  std::ifstream myfile2 (fname2.c_str());
358  if (myfile2.is_open())
359  {
360  while ( myfile2.good() )
361  {
362  getline (myfile2,line, ' ');
363  // remove any potenial control sequences that might be introduced by lines ending in a single space
364  line.erase(std::remove_if(line.begin(), line.end(),
365  [](char c) { return std::isspace(c) || std::iscntrl(c); } ), line.end());
366 
367  if (line.length() > 0 ) // otherwise string contained only control characters before, empty string are converted to 0 by atof resulting in a broken b-value list
368  {
369  bval_entries.push_back(atof(line.c_str()));
370  }
371 
372 
373  }
374  myfile2.close();
375  }
376  else
377  {
378  MITK_INFO << "Unable to open bvals file. Expected name: " << fname2;
379  }
380 
381 
382  BValue = -1;
383  unsigned int numb = bval_entries.size();
384  for(unsigned int i=0; i<numb; i++)
385  {
386 
387  // Take the first entry in bvals as the reference b-value
388  if(BValue == -1 && bval_entries.at(i) != 0)
389  {
390  BValue = bval_entries.at(i);
391  }
392 
393  float b_val = bval_entries.at(i);
394 
395 
396  vnl_vector_fixed< double, 3 > vec;
397  vec[0] = bvec_entries.at(i);
398  vec[1] = bvec_entries.at(i+numb);
399  vec[2] = bvec_entries.at(i+2*numb);
400 
401  // Adjust the vector length to encode gradient strength
402  float factor = b_val/BValue;
403  if(vec.magnitude() > 0)
404  {
405  vec[0] = sqrt(factor)*vec[0];
406  vec[1] = sqrt(factor)*vec[1];
407  vec[2] = sqrt(factor)*vec[2];
408  }
409 
410  DiffusionVectors->InsertElement(i,vec);
411  }
412 
413  for(int i=0; i<3; i++)
414  for(int j=0; j<3; j++)
415  MeasurementFrame[i][j] = i==j ? 1 : 0;
416  }
417 
418  outputForCache = mitk::GrabItkImageMemory( itkVectorImage);
419 
420  // create BValueMap
422  outputForCache->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( DiffusionVectors ) );
423  outputForCache->SetProperty( mitk::DiffusionPropertyHelper::ORIGINALGRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( OriginalDiffusionVectors ) );
424  outputForCache->SetProperty( mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str(), mitk::MeasurementFrameProperty::New( MeasurementFrame ) );
425  outputForCache->SetProperty( mitk::DiffusionPropertyHelper::BVALUEMAPPROPERTYNAME.c_str(), mitk::BValueMapProperty::New( BValueMap ) );
426  outputForCache->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( BValue ) );
427 
428 
429  // Since we have already read the tree, we can store it in a cache variable
430  // so that it can be assigned to the DataObject in GenerateData();
431  m_OutputCache = outputForCache;
432  m_CacheTime.Modified();
433  }
434  catch(std::exception& e)
435  {
436  MITK_INFO << "Std::Exception while reading file!!";
437  MITK_INFO << e.what();
438  throw itk::ImageFileReaderException(__FILE__, __LINE__, e.what());
439  }
440  catch(...)
441  {
442  MITK_INFO << "Exception while reading file!!";
443  throw itk::ImageFileReaderException(__FILE__, __LINE__, "Sorry, an error occurred while reading the requested vessel tree file!");
444  }
445  }
446 }
447 
448 } //namespace MITK
449 
450 #endif
mitk::Point3D PointType
static const std::string REFERENCEBVALUEPROPERTYNAME
static char * line
Definition: svm.cpp:2884
itk::SmartPointer< Self > Pointer
std::map< unsigned int, std::vector< unsigned int > > BValueMap
The BValueMap contains seperated IndicesVectors for each b value (index for GradientDirectionContaine...
#define MITK_INFO
Definition: mitkLogMacros.h:22
std::string GetFilenameWithoutExtension(const std::string &path) const
Provides the filename minus the extension.
static BValueMap CreateBValueMap(const GradientDirectionsContainerType *gdc, float referenceBValue)
DataCollection - Class to facilitate loading/accessing structured data.
Image::Pointer GrabItkImageMemory(itk::SmartPointer< ItkOutputImageType > &itkimage, mitk::Image *mitkImage=nullptr, const BaseGeometry *geometry=nullptr, bool update=true)
Grabs the memory of an itk::Image (with a specific type) and puts it into an mitk::Image.The memory is managed by the mitk::Image after calling this function. The itk::Image remains valid until the mitk::Image decides to free the memory.
static const std::string MEASUREMENTFRAMEPROPERTYNAME
mitk::DiffusionPropertyHelper::MeasurementFrameType MeasurementFrameType
The CustomMimeType class represents a custom mime-type which may be registered as a service object...
virtual std::vector< itk::SmartPointer< BaseData > > Read() override
Reads a path or stream and creates a list of BaseData objects.
us::ServiceRegistration< IFileReader > RegisterService(us::ModuleContext *context=us::GetModuleContext())
Convenience class to temporarily change the current locale.
const CustomMimeType * GetMimeType() const
static Pointer New()
static Pointer New()
Base class for creating mitk::BaseData objects from files or streams.
std::string GetExtension(const std::string &path) const
Provides the first matching extension.
unsigned short PixelType
virtual std::string GetInputLocation() const override
Get the current input location.
static const std::string GRADIENTCONTAINERPROPERTYNAME
static const std::string BVALUEMAPPROPERTYNAME
static const std::string ORIGINALGRADIENTCONTAINERPROPERTYNAME
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.