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  template<typename TO, typename FROM>
62  static TO clamp(const FROM& in);//deal with integer cast being undefined when converting from outside range
63  public:
64  void openRead(const AString& filename);
65  void writeNew(const AString& filename, const NiftiHeader& header, const int& version = 1, const bool& withRead = false, const bool& swapEndian = false);
66  AString getFilename() const { return m_file.getFilename(); }
67  void overrideDimensions(const std::vector<int64_t>& newDims) { m_dims = newDims; }//HACK: deal with reading/writing CIFTI-1's broken headers
68  void close();
69  const NiftiHeader& getHeader() const { return m_header; }
70  const std::vector<int64_t>& getDimensions() const { return m_dims; }
71  int getNumComponents() const;
72  //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
73  //NOTE: you need to provide storage for all components within the range, if getNumComponents() == 3 and fullDims == 0, you need 3 elements allocated
74  template<typename T>
75  void readData(T* dataOut, const int& fullDims, const std::vector<int64_t>& indexSelect, const bool& tolerateShortRead = false);
76  template<typename T>
77  void writeData(const T* dataIn, const int& fullDims, const std::vector<int64_t>& indexSelect);
78  };
79 
80  template<typename T>
81  void NiftiIO::readData(T* dataOut, const int& fullDims, const std::vector<int64_t>& indexSelect, const bool& tolerateShortRead)
82  {
83  if (fullDims < 0) throw CiftiException("NiftiIO: fulldims must not be negative");
84  if (fullDims > (int)m_dims.size()) throw CiftiException("NiftiIO: fulldims must not be greater than number of dimensions");
85  if ((size_t)fullDims + indexSelect.size() != m_dims.size())
86  {//could be >=, but should catch more stupid mistakes as ==
87  throw CiftiException("NiftiIO: fulldims plus length of indexSelect must equal number of dimensions");
88  }
89  int64_t numElems = getNumComponents();//for now, calculate read size on the fly, as the read call will be the slowest part
90  int curDim;
91  for (curDim = 0; curDim < fullDims; ++curDim)
92  {
93  numElems *= m_dims[curDim];
94  }
95  int64_t numDimSkip = numElems, numSkip = 0;
96  for (; curDim < (int)m_dims.size(); ++curDim)
97  {
98  if (indexSelect[curDim - fullDims] < 0) throw CiftiException("NiftiIO: indices must not be negative");
99  if (indexSelect[curDim - fullDims] >= m_dims[curDim]) throw CiftiException("NiftiIO: index exceeds nifti dimension length");
100  numSkip += indexSelect[curDim - fullDims] * numDimSkip;
101  numDimSkip *= m_dims[curDim];
102  }
103  CiftiMutexLocker locked(&m_mutex);//protect starting with resizing until we are done converting, because we use an internal variable for scratch space
104  //we can't guarantee that the output memory is enough to use as scratch space, as we might be doing a narrowing conversion
105  //we are doing FILE ACCESS, so cpu performance isn't really something to worry about
106  m_scratch.resize(numElems * numBytesPerElem());
107  m_file.seek(numSkip * numBytesPerElem() + m_header.getDataOffset());
108  int64_t numRead = 0;
109  m_file.read(m_scratch.data(), m_scratch.size(), &numRead);
110  if ((numRead != (int64_t)m_scratch.size() && !tolerateShortRead) || numRead < 0)//for now, assume read giving -1 is always a problem
111  {
112  throw CiftiException("error while reading from file '" + m_file.getFilename() + "'");
113  }
114  switch (m_header.getDataType())
115  {
116  case NIFTI_TYPE_UINT8:
117  case NIFTI_TYPE_RGB24://handled by components
118  convertRead(dataOut, (uint8_t*)m_scratch.data(), numElems);
119  break;
120  case NIFTI_TYPE_INT8:
121  convertRead(dataOut, (int8_t*)m_scratch.data(), numElems);
122  break;
123  case NIFTI_TYPE_UINT16:
124  convertRead(dataOut, (uint16_t*)m_scratch.data(), numElems);
125  break;
126  case NIFTI_TYPE_INT16:
127  convertRead(dataOut, (int16_t*)m_scratch.data(), numElems);
128  break;
129  case NIFTI_TYPE_UINT32:
130  convertRead(dataOut, (uint32_t*)m_scratch.data(), numElems);
131  break;
132  case NIFTI_TYPE_INT32:
133  convertRead(dataOut, (int32_t*)m_scratch.data(), numElems);
134  break;
135  case NIFTI_TYPE_UINT64:
136  convertRead(dataOut, (uint64_t*)m_scratch.data(), numElems);
137  break;
138  case NIFTI_TYPE_INT64:
139  convertRead(dataOut, (int64_t*)m_scratch.data(), numElems);
140  break;
141  case NIFTI_TYPE_FLOAT32:
142  case NIFTI_TYPE_COMPLEX64://components
143  convertRead(dataOut, (float*)m_scratch.data(), numElems);
144  break;
145  case NIFTI_TYPE_FLOAT64:
147  convertRead(dataOut, (double*)m_scratch.data(), numElems);
148  break;
149  case NIFTI_TYPE_FLOAT128:
151  convertRead(dataOut, (long double*)m_scratch.data(), numElems);
152  break;
153  default:
154  throw CiftiException("internal error, tell the developers what you just tried to do");
155  }
156  }
157 
158  template<typename T>
159  void NiftiIO::writeData(const T* dataIn, const int& fullDims, const std::vector<int64_t>& indexSelect)
160  {
161  if (fullDims < 0) throw CiftiException("NiftiIO: fulldims must not be negative");
162  if (fullDims > (int)m_dims.size()) throw CiftiException("NiftiIO: fulldims must not be greater than number of dimensions");
163  if ((size_t)fullDims + indexSelect.size() != m_dims.size())
164  {//could be >=, but should catch more stupid mistakes as ==
165  throw CiftiException("NiftiIO: fulldims plus length of indexSelect must equal number of dimensions");
166  }
167  int64_t numElems = getNumComponents();//for now, calculate read size on the fly, as the read call will be the slowest part
168  int curDim;
169  for (curDim = 0; curDim < fullDims; ++curDim)
170  {
171  numElems *= m_dims[curDim];
172  }
173  int64_t numDimSkip = numElems, numSkip = 0;
174  for (; curDim < (int)m_dims.size(); ++curDim)
175  {
176  if (indexSelect[curDim - fullDims] < 0) throw CiftiException("NiftiIO: indices must not be negative");
177  if (indexSelect[curDim - fullDims] >= m_dims[curDim]) throw CiftiException("NiftiIO: index exceeds nifti dimension length");
178  numSkip += indexSelect[curDim - fullDims] * numDimSkip;
179  numDimSkip *= m_dims[curDim];
180  }
181  CiftiMutexLocker locked(&m_mutex);//protect starting with resizing until we are done writing, because we use an internal variable for scratch space
182  //we are doing FILE ACCESS, so cpu performance isn't really something to worry about
183  m_scratch.resize(numElems * numBytesPerElem());
184  m_file.seek(numSkip * numBytesPerElem() + m_header.getDataOffset());
185  switch (m_header.getDataType())
186  {
187  case NIFTI_TYPE_UINT8:
188  case NIFTI_TYPE_RGB24://handled by components
189  convertWrite((uint8_t*)m_scratch.data(), dataIn, numElems);
190  break;
191  case NIFTI_TYPE_INT8:
192  convertWrite((int8_t*)m_scratch.data(), dataIn, numElems);
193  break;
194  case NIFTI_TYPE_UINT16:
195  convertWrite((uint16_t*)m_scratch.data(), dataIn, numElems);
196  break;
197  case NIFTI_TYPE_INT16:
198  convertWrite((int16_t*)m_scratch.data(), dataIn, numElems);
199  break;
200  case NIFTI_TYPE_UINT32:
201  convertWrite((uint32_t*)m_scratch.data(), dataIn, numElems);
202  break;
203  case NIFTI_TYPE_INT32:
204  convertWrite((int32_t*)m_scratch.data(), dataIn, numElems);
205  break;
206  case NIFTI_TYPE_UINT64:
207  convertWrite((uint64_t*)m_scratch.data(), dataIn, numElems);
208  break;
209  case NIFTI_TYPE_INT64:
210  convertWrite((int64_t*)m_scratch.data(), dataIn, numElems);
211  break;
212  case NIFTI_TYPE_FLOAT32:
213  case NIFTI_TYPE_COMPLEX64://components
214  convertWrite((float*)m_scratch.data(), dataIn, numElems);
215  break;
216  case NIFTI_TYPE_FLOAT64:
218  convertWrite((double*)m_scratch.data(), dataIn, numElems);
219  break;
220  case NIFTI_TYPE_FLOAT128:
222  convertWrite((long double*)m_scratch.data(), dataIn, numElems);
223  break;
224  default:
225  throw CiftiException("internal error, tell the developers what you just tried to do");
226  }
227  m_file.write(m_scratch.data(), m_scratch.size());
228  }
229 
230  template<typename TO, typename FROM>
231  void NiftiIO::convertRead(TO* out, FROM* in, const int64_t& count)
232  {
233  if (m_header.isSwapped())
234  {
235  ByteSwapping::swapArray(in, count);
236  }
237  double mult, offset;
238  bool doScale = m_header.getDataScaling(mult, offset);
239  if (std::numeric_limits<TO>::is_integer)//do round to nearest when integer output type
240  {
241  if (doScale)
242  {
243  for (int64_t i = 0; i < count; ++i)
244  {
245  out[i] = clamp<TO, long double>(floor(0.5l + offset + mult * (long double)in[i]));//we don't always need that much precision, but it will still be faster than hard drives
246  }
247  } else {
248  for (int64_t i = 0; i < count; ++i)
249  {
250  out[i] = clamp<TO, double>(floor(0.5 + in[i]));
251  }
252  }
253  } else {
254  if (doScale)
255  {
256  for (int64_t i = 0; i < count; ++i)
257  {
258  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
259  }
260  } else {
261  for (int64_t i = 0; i < count; ++i)
262  {
263  out[i] = (TO)in[i];//explicit cast to make sure the compiler doesn't squawk
264  }
265  }
266  }
267  }
268 
269  template<typename TO, typename FROM>
270  void NiftiIO::convertWrite(TO* out, const FROM* in, const int64_t& count)
271  {
272  double mult, offset;
273  bool doScale = m_header.getDataScaling(mult, offset);
274  if (std::numeric_limits<TO>::is_integer)//do round to nearest when integer output type
275  {//TODO: what about NaN?
276  if (doScale)
277  {
278  for (int64_t i = 0; i < count; ++i)
279  {
280  out[i] = clamp<TO, long double>(floor(0.5l + ((long double)in[i] - offset) / mult));//we don't always need that much precision, but it will still be faster than hard drives
281  }
282  } else {
283  for (int64_t i = 0; i < count; ++i)
284  {
285  out[i] = clamp<TO, double>(floor(0.5 + in[i]));
286  }
287  }
288  } else {
289  if (doScale)
290  {
291  for (int64_t i = 0; i < count; ++i)
292  {
293  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
294  }
295  } else {
296  for (int64_t i = 0; i < count; ++i)
297  {
298  out[i] = (TO)in[i];//explicit cast to make sure the compiler doesn't squawk
299  }
300  }
301  }
302  if (m_header.isSwapped()) ByteSwapping::swapArray(out, count);
303  }
304 
305  template<typename TO, typename FROM>
306  TO NiftiIO::clamp(const FROM& in)
307  {
308  typedef std::numeric_limits<TO> mylimits;
309  if (mylimits::has_infinity && std::isinf(in)) return (TO)in;//in case we use this on float types at some point
310  if (mylimits::max() < in) return mylimits::max();
311  if (mylimits::is_integer)//c++11 can use lowest() instead of this mess
312  {
313  if (mylimits::min() > in) return mylimits::min();
314  } else {
315  if (-mylimits::max() > in) return -mylimits::max();
316  }
317  return (TO)in;
318  }
319 }
320 
321 #endif //__NIFTI_IO_H__
Definition: BinaryFile.h:41
Definition: CiftiException.h:40
Definition: NiftiIO.h:50
const int32_t NIFTI_TYPE_UINT64
Definition: nifti1.h:568
const int32_t NIFTI_TYPE_COMPLEX64
Definition: nifti1.h:554
const int32_t NIFTI_TYPE_FLOAT64
Definition: nifti1.h:556
const int32_t NIFTI_TYPE_INT8
Definition: nifti1.h:560
const int32_t NIFTI_TYPE_UINT8
Definition: nifti1.h:546
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_INT64
Definition: nifti1.h:566
const int32_t NIFTI_TYPE_INT32
Definition: nifti1.h:550
const int32_t NIFTI_TYPE_INT16
Definition: nifti1.h:548
const int32_t NIFTI_TYPE_UINT16
Definition: nifti1.h:562
const int32_t NIFTI_TYPE_COMPLEX128
Definition: nifti1.h:572
const int32_t NIFTI_TYPE_FLOAT128
Definition: nifti1.h:570
const int32_t NIFTI_TYPE_COMPLEX256
Definition: nifti1.h:574
const int32_t NIFTI_TYPE_UINT32
Definition: nifti1.h:564
namespace for all CiftiLib functionality
Definition: CiftiBrainModelsMap.h:42
Definition: NiftiHeader.h:49