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
46namespace 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;
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;
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;
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;
142 case NIFTI_TYPE_COMPLEX64://components
143 convertRead(dataOut, (float*)m_scratch.data(), numElems);
144 break;
147 convertRead(dataOut, (double*)m_scratch.data(), numElems);
148 break;
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;
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;
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;
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;
213 case NIFTI_TYPE_COMPLEX64://components
214 convertWrite((float*)m_scratch.data(), dataIn, numElems);
215 break;
218 convertWrite((double*)m_scratch.data(), dataIn, numElems);
219 break;
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