CiftiLib
A C++ library for CIFTI-2 and CIFTI-1 files
NiftiIO.h
1 #ifndef __NIFTI_IO_H__
2 #define __NIFTI_IO_H__
3 
4 /*LICENSE_START*/
5 /*
6  * Copyright (c) 2014, Washington University School of Medicine
7  * All rights reserved.
8  *
9  * Redistribution and use in source and binary forms, with or without modification,
10  * are permitted provided that the following conditions are met:
11  *
12  * 1. Redistributions of source code must retain the above copyright notice,
13  * this list of conditions and the following disclaimer.
14  *
15  * 2. Redistributions in binary form must reproduce the above copyright notice,
16  * this list of conditions and the following disclaimer in the documentation
17  * and/or other materials provided with the distribution.
18  *
19  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
21  * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
23  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
24  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
25  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
26  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
27  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
28  * EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29  */
30 
31 #include "Common/AString.h"
32 
33 #include "Common/ByteSwapping.h"
34 #include "Common/BinaryFile.h"
35 #include "Common/CiftiException.h"
36 #include "Common/CiftiMutex.h"
37 #include "Nifti/NiftiHeader.h"
38 
39 //include MultiDimIterator from a private include directory, in case people want to use it with NiftiIO
40 #include "Common/MultiDimIterator.h"
41 
42 #include <cmath>
43 #include <limits>
44 #include <vector>
45 
46 namespace cifti
47 {
48 
49  class NiftiIO
50  {
51  BinaryFile m_file;
52  NiftiHeader m_header;
53  std::vector<int64_t> m_dims;
54  std::vector<char> m_scratch;//scratch memory for byteswapping, type conversion, etc
55  CiftiMutex m_mutex;
56  int numBytesPerElem();//for resizing scratch
57  template<typename TO, typename FROM>
58  void convertRead(TO* out, FROM* in, const int64_t& count);//for reading from file
59  template<typename TO, typename FROM>
60  void convertWrite(TO* out, const FROM* in, const int64_t& count);//for writing to file
61  public:
62  void openRead(const AString& filename);
63  void writeNew(const AString& filename, const NiftiHeader& header, const int& version = 1, const bool& withRead = false, const bool& swapEndian = false);
64  AString getFilename() const { return m_file.getFilename(); }
65  void overrideDimensions(const std::vector<int64_t>& newDims) { m_dims = newDims; }//HACK: deal with reading/writing CIFTI-1's broken headers
66  void close();
67  const NiftiHeader& getHeader() const { return m_header; }
68  const std::vector<int64_t>& getDimensions() const { return m_dims; }
69  int getNumComponents() const;
70  //to read/write 1 frame of a standard volume file, call with fullDims = 3, indexSelect containing indexes for any of dims 4-7 that exist
71  //NOTE: you need to provide storage for all components within the range, if getNumComponents() == 3 and fullDims == 0, you need 3 elements allocated
72  template<typename T>
73  void readData(T* dataOut, const int& fullDims, const std::vector<int64_t>& indexSelect, const bool& tolerateShortRead = false);
74  template<typename T>
75  void writeData(const T* dataIn, const int& fullDims, const std::vector<int64_t>& indexSelect);
76  };
77 
78  template<typename T>
79  void NiftiIO::readData(T* dataOut, const int& fullDims, const std::vector<int64_t>& indexSelect, const bool& tolerateShortRead)
80  {
81  if (fullDims < 0) throw CiftiException("NiftiIO: fulldims must not be negative");
82  if (fullDims > (int)m_dims.size()) throw CiftiException("NiftiIO: fulldims must not be greater than number of dimensions");
83  if ((size_t)fullDims + indexSelect.size() != m_dims.size())
84  {//could be >=, but should catch more stupid mistakes as ==
85  throw CiftiException("NiftiIO: fulldims plus length of indexSelect must equal number of dimensions");
86  }
87  int64_t numElems = getNumComponents();//for now, calculate read size on the fly, as the read call will be the slowest part
88  int curDim;
89  for (curDim = 0; curDim < fullDims; ++curDim)
90  {
91  numElems *= m_dims[curDim];
92  }
93  int64_t numDimSkip = numElems, numSkip = 0;
94  for (; curDim < (int)m_dims.size(); ++curDim)
95  {
96  if (indexSelect[curDim - fullDims] < 0) throw CiftiException("NiftiIO: indices must not be negative");
97  if (indexSelect[curDim - fullDims] >= m_dims[curDim]) throw CiftiException("NiftiIO: index exceeds nifti dimension length");
98  numSkip += indexSelect[curDim - fullDims] * numDimSkip;
99  numDimSkip *= m_dims[curDim];
100  }
101  CiftiMutexLocker locked(&m_mutex);//protect starting with resizing until we are done converting, because we use an internal variable for scratch space
102  //we can't guarantee that the output memory is enough to use as scratch space, as we might be doing a narrowing conversion
103  //we are doing FILE ACCESS, so cpu performance isn't really something to worry about
104  m_scratch.resize(numElems * numBytesPerElem());
105  m_file.seek(numSkip * numBytesPerElem() + m_header.getDataOffset());
106  int64_t numRead = 0;
107  m_file.read(m_scratch.data(), m_scratch.size(), &numRead);
108  if ((numRead != (int64_t)m_scratch.size() && !tolerateShortRead) || numRead < 0)//for now, assume read giving -1 is always a problem
109  {
110  throw CiftiException("error while reading from file '" + m_file.getFilename() + "'");
111  }
112  switch (m_header.getDataType())
113  {
114  case NIFTI_TYPE_UINT8:
115  case NIFTI_TYPE_RGB24://handled by components
116  convertRead(dataOut, (uint8_t*)m_scratch.data(), numElems);
117  break;
118  case NIFTI_TYPE_INT8:
119  convertRead(dataOut, (int8_t*)m_scratch.data(), numElems);
120  break;
121  case NIFTI_TYPE_UINT16:
122  convertRead(dataOut, (uint16_t*)m_scratch.data(), numElems);
123  break;
124  case NIFTI_TYPE_INT16:
125  convertRead(dataOut, (int16_t*)m_scratch.data(), numElems);
126  break;
127  case NIFTI_TYPE_UINT32:
128  convertRead(dataOut, (uint32_t*)m_scratch.data(), numElems);
129  break;
130  case NIFTI_TYPE_INT32:
131  convertRead(dataOut, (int32_t*)m_scratch.data(), numElems);
132  break;
133  case NIFTI_TYPE_UINT64:
134  convertRead(dataOut, (uint64_t*)m_scratch.data(), numElems);
135  break;
136  case NIFTI_TYPE_INT64:
137  convertRead(dataOut, (int64_t*)m_scratch.data(), numElems);
138  break;
139  case NIFTI_TYPE_FLOAT32:
140  case NIFTI_TYPE_COMPLEX64://components
141  convertRead(dataOut, (float*)m_scratch.data(), numElems);
142  break;
143  case NIFTI_TYPE_FLOAT64:
145  convertRead(dataOut, (double*)m_scratch.data(), numElems);
146  break;
147  case NIFTI_TYPE_FLOAT128:
149  convertRead(dataOut, (long double*)m_scratch.data(), numElems);
150  break;
151  default:
152  throw CiftiException("internal error, tell the developers what you just tried to do");
153  }
154  }
155 
156  template<typename T>
157  void NiftiIO::writeData(const T* dataIn, const int& fullDims, const std::vector<int64_t>& indexSelect)
158  {
159  if (fullDims < 0) throw CiftiException("NiftiIO: fulldims must not be negative");
160  if (fullDims > (int)m_dims.size()) throw CiftiException("NiftiIO: fulldims must not be greater than number of dimensions");
161  if ((size_t)fullDims + indexSelect.size() != m_dims.size())
162  {//could be >=, but should catch more stupid mistakes as ==
163  throw CiftiException("NiftiIO: fulldims plus length of indexSelect must equal number of dimensions");
164  }
165  int64_t numElems = getNumComponents();//for now, calculate read size on the fly, as the read call will be the slowest part
166  int curDim;
167  for (curDim = 0; curDim < fullDims; ++curDim)
168  {
169  numElems *= m_dims[curDim];
170  }
171  int64_t numDimSkip = numElems, numSkip = 0;
172  for (; curDim < (int)m_dims.size(); ++curDim)
173  {
174  if (indexSelect[curDim - fullDims] < 0) throw CiftiException("NiftiIO: indices must not be negative");
175  if (indexSelect[curDim - fullDims] >= m_dims[curDim]) throw CiftiException("NiftiIO: index exceeds nifti dimension length");
176  numSkip += indexSelect[curDim - fullDims] * numDimSkip;
177  numDimSkip *= m_dims[curDim];
178  }
179  CiftiMutexLocker locked(&m_mutex);//protect starting with resizing until we are done writing, because we use an internal variable for scratch space
180  //we are doing FILE ACCESS, so cpu performance isn't really something to worry about
181  m_scratch.resize(numElems * numBytesPerElem());
182  m_file.seek(numSkip * numBytesPerElem() + m_header.getDataOffset());
183  switch (m_header.getDataType())
184  {
185  case NIFTI_TYPE_UINT8:
186  case NIFTI_TYPE_RGB24://handled by components
187  convertWrite((uint8_t*)m_scratch.data(), dataIn, numElems);
188  break;
189  case NIFTI_TYPE_INT8:
190  convertWrite((int8_t*)m_scratch.data(), dataIn, numElems);
191  break;
192  case NIFTI_TYPE_UINT16:
193  convertWrite((uint16_t*)m_scratch.data(), dataIn, numElems);
194  break;
195  case NIFTI_TYPE_INT16:
196  convertWrite((int16_t*)m_scratch.data(), dataIn, numElems);
197  break;
198  case NIFTI_TYPE_UINT32:
199  convertWrite((uint32_t*)m_scratch.data(), dataIn, numElems);
200  break;
201  case NIFTI_TYPE_INT32:
202  convertWrite((int32_t*)m_scratch.data(), dataIn, numElems);
203  break;
204  case NIFTI_TYPE_UINT64:
205  convertWrite((uint64_t*)m_scratch.data(), dataIn, numElems);
206  break;
207  case NIFTI_TYPE_INT64:
208  convertWrite((int64_t*)m_scratch.data(), dataIn, numElems);
209  break;
210  case NIFTI_TYPE_FLOAT32:
211  case NIFTI_TYPE_COMPLEX64://components
212  convertWrite((float*)m_scratch.data(), dataIn, numElems);
213  break;
214  case NIFTI_TYPE_FLOAT64:
216  convertWrite((double*)m_scratch.data(), dataIn, numElems);
217  break;
218  case NIFTI_TYPE_FLOAT128:
220  convertWrite((long double*)m_scratch.data(), dataIn, numElems);
221  break;
222  default:
223  throw CiftiException("internal error, tell the developers what you just tried to do");
224  }
225  m_file.write(m_scratch.data(), m_scratch.size());
226  }
227 
228  template<typename TO, typename FROM>
229  void NiftiIO::convertRead(TO* out, FROM* in, const int64_t& count)
230  {
231  if (m_header.isSwapped())
232  {
233  ByteSwapping::swapArray(in, count);
234  }
235  double mult, offset;
236  bool doScale = m_header.getDataScaling(mult, offset);
237  if (std::numeric_limits<TO>::is_integer)//do round to nearest when integer output type
238  {
239  if (doScale)
240  {
241  for (int64_t i = 0; i < count; ++i)
242  {
243  out[i] = (TO)floor(0.5 + offset + mult * (long double)in[i]);//we don't always need that much precision, but it will still be faster than hard drives
244  }
245  } else {
246  for (int64_t i = 0; i < count; ++i)
247  {
248  out[i] = (TO)floor(0.5 + in[i]);
249  }
250  }
251  } else {
252  if (doScale)
253  {
254  for (int64_t i = 0; i < count; ++i)
255  {
256  out[i] = (TO)(offset + mult * (long double)in[i]);//we don't always need that much precision, but it will still be faster than hard drives
257  }
258  } else {
259  for (int64_t i = 0; i < count; ++i)
260  {
261  out[i] = (TO)in[i];//explicit cast to make sure the compiler doesn't squawk
262  }
263  }
264  }
265  }
266 
267  template<typename TO, typename FROM>
268  void NiftiIO::convertWrite(TO* out, const FROM* in, const int64_t& count)
269  {
270  double mult, offset;
271  bool doScale = m_header.getDataScaling(mult, offset);
272  if (std::numeric_limits<TO>::is_integer)//do round to nearest when integer output type
273  {
274  if (doScale)
275  {
276  for (int64_t i = 0; i < count; ++i)
277  {
278  out[i] = (TO)floor(0.5 + ((long double)in[i] - offset) / mult);//we don't always need that much precision, but it will still be faster than hard drives
279  }
280  } else {
281  for (int64_t i = 0; i < count; ++i)
282  {
283  out[i] = (TO)floor(0.5 + in[i]);
284  }
285  }
286  } else {
287  if (doScale)
288  {
289  for (int64_t i = 0; i < count; ++i)
290  {
291  out[i] = (TO)(((long double)in[i] - offset) / mult);//we don't always need that much precision, but it will still be faster than hard drives
292  }
293  } else {
294  for (int64_t i = 0; i < count; ++i)
295  {
296  out[i] = (TO)in[i];//explicit cast to make sure the compiler doesn't squawk
297  }
298  }
299  }
300  if (m_header.isSwapped()) ByteSwapping::swapArray(out, count);
301  }
302 
303 }
304 
305 #endif //__NIFTI_IO_H__
Definition: CiftiException.h:40
namespace for all CiftiLib functionality
Definition: CiftiBrainModelsMap.h:41
const int32_t NIFTI_TYPE_COMPLEX256
Definition: nifti1.h:574
const int32_t NIFTI_TYPE_INT8
Definition: nifti1.h:560
const int32_t NIFTI_TYPE_COMPLEX128
Definition: nifti1.h:572
const int32_t NIFTI_TYPE_UINT64
Definition: nifti1.h:568
const int32_t NIFTI_TYPE_UINT8
Definition: nifti1.h:546
const int32_t NIFTI_TYPE_FLOAT128
Definition: nifti1.h:570
Definition: NiftiIO.h:49
const int32_t NIFTI_TYPE_INT32
Definition: nifti1.h:550
const int32_t NIFTI_TYPE_UINT32
Definition: nifti1.h:564
const int32_t NIFTI_TYPE_INT64
Definition: nifti1.h:566
const int32_t NIFTI_TYPE_FLOAT32
Definition: nifti1.h:552
const int32_t NIFTI_TYPE_RGB24
Definition: nifti1.h:558
const int32_t NIFTI_TYPE_COMPLEX64
Definition: nifti1.h:554
Definition: NiftiHeader.h:48
const int32_t NIFTI_TYPE_INT16
Definition: nifti1.h:548
Definition: BinaryFile.h:40
const int32_t NIFTI_TYPE_UINT16
Definition: nifti1.h:562
const int32_t NIFTI_TYPE_FLOAT64
Definition: nifti1.h:556