<commit_before>/** \file Nifti-inl.h
 \brief Definitions of templated IO code for Nifti
 - Written by <NAME>, IoP KCL
 - Based on nifti1_io.h (Thanks to <NAME> et al)
 - This code is released to the public domain. Do with it what you will.
 - This file should NOT be included directly in projects. It is just to 
 - separate the interface from the definition a bit better.
 */

#ifndef NIFTI_NIFTI_INL
#define NIFTI_NIFTI_INL

/**
  *   Internal function to convert the internal NIfTI data to the desired dataype.
  *
  *   Converts a sequence of bytes from a NIfTI image to the templated datatype.
  *   If the file contains complex data but a real datatype is requested then the magnitude is returned.
  *   If the file contains real data but a complex datatype is requested then the imaginary part will be 0.
  *
  *   The number of elements is implicit in the size of the input byte array as the NIfTI stores the element size.
  *   The output array can be larger than required, in which case offset can be non-zero and data will be written
  *   from that point onwards.
  *   
  *   @param T Desired datatype. Valid (scalar) conversions must exist.
  *   @param bytes  - Raw bytes read from NIfTI file.
  *   @param data   - Storage for converted data. Must contain enough space.
  *   @param offset - Location in data to start writing.
  */
template<typename T> void Nifti::convertFromBytes(const std::vector<char> &bytes,
                                                  const typename std::vector<T>::iterator begin,
												  const typename std::vector<T>::iterator end) {
	size_t nEl = bytes.size() / m_typeinfo.size;
	assert(nEl == std::distance(begin, end));
	T sc_sl = static_cast<T>(scaling_slope);
	T sc_in = static_cast<T>(scaling_inter);
	auto el = begin;
	for (size_t i = 0; i < nEl; i++) {
		switch (m_typeinfo.type) {
			case DataType::INT8:       *el = sc_in + sc_sl * static_cast<T>(reinterpret_cast<const int8_t *>(bytes.data())[i]); break;
			case DataType::INT16:      *el = sc_in + sc_sl * static_cast<T>(reinterpret_cast<const int16_t *>(bytes.data())[i]); break;
			case DataType::INT32:      *el = sc_in + sc_sl * static_cast<T>(reinterpret_cast<const int32_t *>(bytes.data())[i]); break;
			case DataType::INT64:      *el = sc_in + sc_sl * static_cast<T>(reinterpret_cast<const int64_t *>(bytes.data())[i]); break;
			case DataType::UINT8:      *el = sc_in + sc_sl * static_cast<T>(reinterpret_cast<const uint8_t *>(bytes.data())[i]); break;
			case DataType::UINT16:     *el = sc_in + sc_sl * static_cast<T>(reinterpret_cast<const uint16_t *>(bytes.data())[i]); break;
			case DataType::UINT32:     *el = sc_in + sc_sl * static_cast<T>(reinterpret_cast<const uint32_t *>(bytes.data())[i]); break;
			case DataType::UINT64:     *el = sc_in + sc_sl * static_cast<T>(reinterpret_cast<const uint64_t *>(bytes.data())[i]); break;
			case DataType::FLOAT32:    *el = sc_in + sc_sl * static_cast<T>(reinterpret_cast<const float *>(bytes.data())[i]); break;
			case DataType::FLOAT64:    *el = sc_in + sc_sl * static_cast<T>(reinterpret_cast<const double *>(bytes.data())[i]); break;
			case DataType::FLOAT128:   *el = sc_in + sc_sl * static_cast<T>(reinterpret_cast<const long double *>(bytes.data())[i]); break;
			case DataType::COMPLEX64:  *el = sc_in + sc_sl * static_cast<T>(abs(reinterpret_cast<const std::complex<float> *>(bytes.data())[i])); break;
			case DataType::COMPLEX128: *el = sc_in + sc_sl * static_cast<T>(abs(reinterpret_cast<const std::complex<double> *>(bytes.data())[i])); break;
			case DataType::COMPLEX256: *el = sc_in + sc_sl * static_cast<T>(abs(reinterpret_cast<const std::complex<long double> *>(bytes.data())[i])); break;
			case DataType::RGB24: case DataType::RGBA32:
				throw(std::runtime_error("Unsupported datatype.")); break;
		}
		el++;
	}
}

template<typename T> void Nifti::convertFromBytes(const std::vector<char> &bytes,
                                                  const typename std::vector<std::complex<T>>::iterator begin,
												  const typename std::vector<std::complex<T>>::iterator end) {
	size_t nEl = bytes.size() / m_typeinfo.size;
	assert(nEl == std::distance(begin, end));
	T sc_sl = static_cast<T>(scaling_slope);
	T sc_in = static_cast<T>(scaling_inter);
	auto el = begin;
	for (size_t i = 0; i < nEl; i++) {
		switch (m_typeinfo.type) {
			case DataType::INT8:       *el = sc_in + sc_sl * std::complex<T>(static_cast<T>(reinterpret_cast<const int8_t *>(bytes.data())[i]), 0.); break;
			case DataType::INT16:      *el = sc_in + sc_sl * std::complex<T>(static_cast<T>(reinterpret_cast<const int16_t *>(bytes.data())[i]), 0.); break;
			case DataType::INT32:      *el = sc_in + sc_sl * std::complex<T>(static_cast<T>(reinterpret_cast<const int32_t *>(bytes.data())[i]), 0.); break;
			case DataType::INT64:      *el = sc_in + sc_sl * std::complex<T>(static_cast<T>(reinterpret_cast<const int64_t *>(bytes.data())[i]), 0.); break;
			case DataType::UINT8:      *el = sc_in + sc_sl * std::complex<T>(static_cast<T>(reinterpret_cast<const uint8_t *>(bytes.data())[i]), 0.); break;
			case DataType::UINT16:     *el = sc_in + sc_sl * std::complex<T>(static_cast<T>(reinterpret_cast<const uint16_t *>(bytes.data())[i]), 0.); break;
			case DataType::UINT32:     *el = sc_in + sc_sl * std::complex<T>(static_cast<T>(reinterpret_cast<const uint32_t *>(bytes.data())[i]), 0.); break;
			case DataType::UINT64:     *el = sc_in + sc_sl * std::complex<T>(static_cast<T>(reinterpret_cast<const uint64_t *>(bytes.data())[i]), 0.); break;
			case DataType::FLOAT32:    *el = sc_in + sc_sl * std::complex<T>(static_cast<T>(reinterpret_cast<const float *>(bytes.data())[i]), 0.); break;
			case DataType::FLOAT64:    *el = sc_in + sc_sl * std::complex<T>(static_cast<T>(reinterpret_cast<const double *>(bytes.data())[i]), 0.); break;
			case DataType::FLOAT128:   *el = sc_in + sc_sl * std::complex<T>(static_cast<T>(reinterpret_cast<const long double *>(bytes.data())[i]), 0.); break;
			case DataType::COMPLEX64:  *el = sc_in + sc_sl * static_cast<std::complex<T> >(reinterpret_cast<const std::complex<float> *>(bytes.data())[i]); break;
			case DataType::COMPLEX128: *el = sc_in + sc_sl * static_cast<std::complex<T> >(reinterpret_cast<const std::complex<double> *>(bytes.data())[i]); break;
			case DataType::COMPLEX256: *el = sc_in + sc_sl * static_cast<std::complex<T> >(reinterpret_cast<const std::complex<long double> *>(bytes.data())[i]); break;
			case DataType::RGB24: case DataType::RGBA32:
				throw(std::runtime_error("Unsupported datatype.")); break;
		}
		el++;
	}
}

/**
  *   Internal function to convert from arbitrary data to the internal NIfTI dataype.
  *
  *   Converts a from the templated datatype to a sequence of bytes.
  *   If complex data is passed in but the file datatype is not complex, then the absolute magnitude is taken.
  *   If non-complex data is passed in but the file datatype is complex, then the imaginary part will be 0.
  *   
  *   The number of elements is implicit in the size of the input byte array as the NIfTI stores the element size.
  *   The output array can be larger than required, in which case offset can be non-zero and data will be written
  *   from that point onwards.
  *   
  *   @param T Desired datatype. Valid (scalar) conversions must exist.
  *   @param bytes  - Raw bytes read from NIfTI file.
  *   @param data   - Storage for converted data. Must contain enough space.
  *   @param offset - Location in data to start writing.
  */
template<typename T> void Nifti::convertToBytes(const typename std::vector<T>::iterator begin,
                                                const typename std::vector<T>::iterator end,
												std::vector<char> &bytes) {
	size_t nEl = std::distance(begin, end);
	assert(nEl == bytes.size() / m_typeinfo.size);
	T sc_sl = static_cast<T>(scaling_slope);
	T sc_in = static_cast<T>(scaling_inter);
	auto el = begin;
	for (size_t i = 0; i < nEl; i++) {
		switch (m_typeinfo.type) {
			case DataType::INT8:              reinterpret_cast<char *>(bytes.data())[i] = static_cast<int8_t>(*el / sc_sl - sc_in); break;
			case DataType::INT16:            reinterpret_cast<short *>(bytes.data())[i] = static_cast<int16_t>(*el / sc_sl - sc_in); break;
			case DataType::INT32:              reinterpret_cast<int *>(bytes.data())[i] = static_cast<int32_t>(*el / sc_sl - sc_in); break;
			case DataType::INT64:             reinterpret_cast<long *>(bytes.data())[i] = static_cast<int64_t>(*el / sc_sl - sc_in); break;
			case DataType::UINT8:    reinterpret_cast<unsigned char *>(bytes.data())[i] = static_cast<uint8_t>(*el / sc_sl - sc_in); break;
			case DataType::UINT16:  reinterpret_cast<unsigned short *>(bytes.data())[i] = static_cast<uint16_t>(*el / sc_sl - sc_in); break;
			case DataType::UINT32:    reinterpret_cast<unsigned int *>(bytes.data())[i] = static_cast<uint32_t>(*el / sc_sl - sc_in); break;
			case DataType::UINT64:   reinterpret_cast<unsigned long *>(bytes.data())[i] = static_cast<uint64_t>(*el / sc_sl - sc_in); break;
			case DataType::FLOAT32:          reinterpret_cast<float *>(bytes.data())[i] = static_cast<float>(*el / sc_sl - sc_in); break;
			case DataType::FLOAT64:         reinterpret_cast<double *>(bytes.data())[i] = static_cast<double>(*el / sc_sl - sc_in); break;
			case DataType::FLOAT128:   reinterpret_cast<long double *>(bytes.data())[i] = static_cast<long double>(*el / sc_sl - sc_in); break;
			case DataType::COMPLEX64:  reinterpret_cast<std::complex<float> *>(bytes.data())[i] = std::complex<float>(static_cast<float>(*el / sc_sl - sc_in)); break;
			case DataType::COMPLEX128: reinterpret_cast<std::complex<double> *>(bytes.data())[i] = std::complex<double>(static_cast<double>(*el / sc_sl - sc_in)); break;
			case DataType::COMPLEX256: reinterpret_cast<std::complex<long double> *>(bytes.data())[i] = std::complex<long double>(static_cast<long double>(*el / sc_sl - sc_in)); break;
			case DataType::RGB24: case DataType::RGBA32:
				throw(std::runtime_error("Unsupported datatype.")); break;
		}
		el++;
	}
}

template<typename T> void Nifti::convertToBytes(const typename std::vector<std::complex<T>>::iterator begin,
												const typename std::vector<std::complex<T>>::iterator end,
												std::vector<char> &bytes) {
	size_t nEl = std::distance(begin, end);
	assert(nEl == bytes.size() / m_typeinfo.size);
	T sc_sl = static_cast<T>(scaling_slope);
	T sc_in = static_cast<T>(scaling_inter);
	auto el = begin;
	for (int i = 0; i < nEl; i++) {
		switch (m_typeinfo.type) {
			case DataType::INT8:              reinterpret_cast<char *>(bytes.data())[i] = static_cast<int8_t>(abs(*el / sc_sl - sc_in)); break;
			case DataType::INT16:            reinterpret_cast<short *>(bytes.data())[i] = static_cast<int16_t>(abs(*el / sc_sl - sc_in)); break;
			case DataType::INT32:              reinterpret_cast<int *>(bytes.data())[i] = static_cast<int32_t>(abs(*el / sc_sl - sc_in)); break;
			case DataType::INT64:             reinterpret_cast<long *>(bytes.data())[i] = static_cast<int64_t>(abs(*el / sc_sl - sc_in)); break;
			case DataType::UINT8:    reinterpret_cast<unsigned char *>(bytes.data())[i] = static_cast<uint8_t>(abs(*el / sc_sl - sc_in)); break;
			case DataType::UINT16:  reinterpret_cast<unsigned short *>(bytes.data())[i] = static_cast<uint16_t>(abs(*el / sc_sl - sc_in)); break;
			case DataType::UINT32:    reinterpret_cast<unsigned int *>(bytes.data())[i] = static_cast<uint32_t>(abs(*el / sc_sl - sc_in)); break;
			case DataType::UINT64:   reinterpret_cast<unsigned long *>(bytes.data())[i] = static_cast<uint64_t>(abs(*el / sc_sl - sc_in)); break;
			case DataType::FLOAT32:          reinterpret_cast<float *>(bytes.data())[i] = static_cast<float>(abs(*el / sc_sl - sc_in)); break;
			case DataType::FLOAT64:         reinterpret_cast<double *>(bytes.data())[i] = static_cast<double>(abs(*el / sc_sl - sc_in)); break;
			case DataType::FLOAT128:   reinterpret_cast<long double *>(bytes.data())[i] = static_cast<long double>(abs(*el / sc_sl - sc_in)); break;
			case DataType::COMPLEX64:  reinterpret_cast<std::complex<float> *>(bytes.data())[i] = static_cast<std::complex<float> >(*el / sc_sl - sc_in); break;
			case DataType::COMPLEX128: reinterpret_cast<std::complex<double> *>(bytes.data())[i] = static_cast<std::complex<double> >(*el / sc_sl - sc_in); break;
			case DataType::COMPLEX256: reinterpret_cast<std::complex<long double> *>(bytes.data())[i] = static_cast<std::complex<long double> >(*el / sc_sl - sc_in); break;
			case DataType::RGB24: case DataType::RGBA32:
				throw(std::runtime_error("Unsupported datatype.")); break;
		}
		el++;
	}
}

/*
 *   Core IO routine. Depending on the mode of the file, reads/writes a
 *   contiguous subregion of the region and converts to/from the desired datatype.
 *
 *   start and size can be short of the full number of dimensions, e.g. if you
 *   only want a part of the first volume in a multi-volume image, then they can
 *   just be 3D. The only limitation is that they must not have more dimensions
 *   than the image.
 *
 *   @param start The voxel indices of the first desired voxel.
 *   @param size  The size of the desired subregion.
 *   @param data  Storage for the data to read/write. Must be sufficiently large.
 */
template<typename T> void Nifti::readWriteVoxels(const Eigen::Ref<ArrayXs> &start, const Eigen::Ref<ArrayXs> &size, std::vector<T> &data) {
	if (start.rows() != size.rows()) throw(std::out_of_range("Start and size must have same dimension in image: " + imagePath()));<commit_msg>Replaced conversion for/switch loops with switch/for loops and macros for readability. Ideally these would templated too, but this means defining 8 new template functions to cope with all the real/complex conversion combinations. Macros seem a better option now, especially as they are undef'd in the functions.<commit_after>/** \file Nifti-inl.h
 \brief Definitions of templated IO code for Nifti
 - Written by <NAME>, IoP KCL
 - Based on nifti1_io.h (Thanks to <NAME> et al)
 - This code is released to the public domain. Do with it what you will.
 - This file should NOT be included directly in projects. It is just to 
 - separate the interface from the definition a bit better.
 */

#ifndef NIFTI_NIFTI_INL
#define NIFTI_NIFTI_INL

/**
  *   Internal function to convert the internal NIfTI data to the desired dataype.
  *
  *   Converts a sequence of bytes from a NIfTI image to the templated datatype.
  *   If the file contains complex data but a real datatype is requested then the magnitude is returned.
  *   If the file contains real data but a complex datatype is requested then the imaginary part will be 0.
  *
  *   The number of elements is implicit in the size of the input byte array as the NIfTI stores the element size.
  *   The output array can be larger than required, in which case offset can be non-zero and data will be written
  *   from that point onwards.
  *   
  *   @param T Desired datatype. Valid (scalar) conversions must exist.
  *   @param bytes  - Raw bytes read from NIfTI file.
  *   @param data   - Storage for converted data. Must contain enough space.
  *   @param offset - Location in data to start writing.
  */
template<typename T> void Nifti::convertFromBytes(const std::vector<char> &bytes,
                                                  const typename std::vector<T>::iterator begin,
												  const typename std::vector<T>::iterator end) {
	size_t nEl = bytes.size() / m_typeinfo.size;
	assert(nEl == std::distance(begin, end));
	T sc_sl = static_cast<T>(scaling_slope);
	T sc_in = static_cast<T>(scaling_inter);
	auto el = begin;
	#define DECL_PTR( BYTE_TYPE ) auto p = reinterpret_cast<const BYTE_TYPE *>(bytes.data())
	#define REAL_LOOP( BYTE_TYPE ) DECL_PTR(BYTE_TYPE); while (el != end) { *el = sc_in + sc_sl * static_cast<T>(*p); el++; p++; }
	#define COMP_LOOP( BYTE_TYPE ) DECL_PTR(std::complex<BYTE_TYPE>); while (el != end) { *el = sc_in + sc_sl * static_cast<T>(abs(*p)); el++; p++; }
	switch (m_typeinfo.type) {
		case DataType::INT8:       { REAL_LOOP(int8_t); }; break;
		case DataType::INT16:      { REAL_LOOP(int16_t); }; break;
		case DataType::INT32:      { REAL_LOOP(int32_t); }; break;
		case DataType::INT64:      { REAL_LOOP(int64_t); }; break;
		case DataType::UINT8:      { REAL_LOOP(uint8_t); }; break;
		case DataType::UINT16:     { REAL_LOOP(uint16_t); }; break;
		case DataType::UINT32:     { REAL_LOOP(uint32_t); }; break;
		case DataType::UINT64:     { REAL_LOOP(uint64_t); }; break;
		case DataType::FLOAT32:    { REAL_LOOP(float); }; break;
		case DataType::FLOAT64:    { REAL_LOOP(double); }; break;
		case DataType::FLOAT128:   { REAL_LOOP(long double); }; break;
		case DataType::COMPLEX64:  { COMP_LOOP(float); }; break;
		case DataType::COMPLEX128: { COMP_LOOP(double); }; break;
		case DataType::COMPLEX256: { COMP_LOOP(long double); }; break;
		case DataType::RGB24: case DataType::RGBA32:
			throw(std::runtime_error("Unsupported datatype.")); break;
	}
	#undef COMP_LOOP
	#undef REAL_LOOP
	#undef DECL_PTR
}

template<typename T> void Nifti::convertFromBytes(const std::vector<char> &bytes,
                                                  const typename std::vector<std::complex<T>>::iterator begin,
												  const typename std::vector<std::complex<T>>::iterator end) {
	size_t nEl = bytes.size() / m_typeinfo.size;
	assert(nEl == std::distance(begin, end));
	T sc_sl = static_cast<T>(scaling_slope);
	T sc_in = static_cast<T>(scaling_inter);
	auto el = begin;
	#define DECL_PTR( BYTE_TYPE ) auto p = reinterpret_cast<const BYTE_TYPE *>(bytes.data())
	#define REAL_LOOP( BYTE_TYPE ) DECL_PTR(BYTE_TYPE); while (el != end) { *el = std::complex<T>(sc_in + sc_sl * static_cast<T>(*p), 0.); el++; p++; }
	#define COMP_LOOP( BYTE_TYPE ) DECL_PTR(std::complex<BYTE_TYPE>); while (el != end) { *el = sc_in + sc_sl * static_cast<std::complex<T>>(*p); el++; p++; }
	switch (m_typeinfo.type) {
		case DataType::INT8:       { REAL_LOOP(int8_t); }; break;
		case DataType::INT16:      { REAL_LOOP(int16_t); }; break;
		case DataType::INT32:      { REAL_LOOP(int32_t); }; break;
		case DataType::INT64:      { REAL_LOOP(int64_t); }; break;
		case DataType::UINT8:      { REAL_LOOP(uint8_t); }; break;
		case DataType::UINT16:     { REAL_LOOP(uint16_t); }; break;
		case DataType::UINT32:     { REAL_LOOP(uint32_t); }; break;
		case DataType::UINT64:     { REAL_LOOP(uint64_t); }; break;
		case DataType::FLOAT32:    { REAL_LOOP(float); }; break;
		case DataType::FLOAT64:    { REAL_LOOP(double); }; break;
		case DataType::FLOAT128:   { REAL_LOOP(long double); }; break;
		case DataType::COMPLEX64:  { COMP_LOOP(float); }; break;
		case DataType::COMPLEX128: { COMP_LOOP(double); }; break;
		case DataType::COMPLEX256: { COMP_LOOP(long double); }; break;
		case DataType::RGB24: case DataType::RGBA32:
			throw(std::runtime_error("Unsupported datatype.")); break;
	}
	#undef COMP_LOOP
	#undef REAL_LOOP
	#undef DECL_PTR
}

/**
  *   Internal function to convert from arbitrary data to the internal NIfTI dataype.
  *
  *   Converts a from the templated datatype to a sequence of bytes.
  *   If complex data is passed in but the file datatype is not complex, then the absolute magnitude is taken.
  *   If non-complex data is passed in but the file datatype is complex, then the imaginary part will be 0.
  *   
  *   The number of elements is implicit in the size of the input byte array as the NIfTI stores the element size.
  *   The output array can be larger than required, in which case offset can be non-zero and data will be written
  *   from that point onwards.
  *   
  *   @param T Desired datatype. Valid (scalar) conversions must exist.
  *   @param bytes  - Raw bytes read from NIfTI file.
  *   @param data   - Storage for converted data. Must contain enough space.
  *   @param offset - Location in data to start writing.
  */
template<typename T> void Nifti::convertToBytes(const typename std::vector<T>::iterator begin,
                                                const typename std::vector<T>::iterator end,
												std::vector<char> &bytes) {
	size_t nEl = std::distance(begin, end);
	assert(nEl == bytes.size() / m_typeinfo.size);
	T sc_sl = static_cast<T>(scaling_slope);
	T sc_in = static_cast<T>(scaling_inter);
	auto el = begin;
	#define DECL_PTR( BYTE_TYPE ) auto p = reinterpret_cast<BYTE_TYPE *>(bytes.data())
	#define REAL_LOOP( BYTE_TYPE ) DECL_PTR( BYTE_TYPE ); while (el != end) { *p = static_cast<BYTE_TYPE>(*el / sc_sl - sc_in); el++; p++; }
	#define COMP_LOOP( BYTE_TYPE ) DECL_PTR( std::complex<BYTE_TYPE> ); while (el != end) { *p = std::complex<BYTE_TYPE>(*el / sc_sl - sc_in, 0.); el++; p++; }
	switch (m_typeinfo.type) {
		case DataType::INT8:       { REAL_LOOP(int8_t); }; break;
		case DataType::INT16:      { REAL_LOOP(int16_t); }; break;
		case DataType::INT32:      { REAL_LOOP(int32_t); }; break;
		case DataType::INT64:      { REAL_LOOP(int64_t); }; break;
		case DataType::UINT8:      { REAL_LOOP(uint8_t); }; break;
		case DataType::UINT16:     { REAL_LOOP(uint16_t); }; break;
		case DataType::UINT32:     { REAL_LOOP(uint32_t); }; break;
		case DataType::UINT64:     { REAL_LOOP(uint64_t); }; break;
		case DataType::FLOAT32:    { REAL_LOOP(float); }; break;
		case DataType::FLOAT64:    { REAL_LOOP(double); }; break;
		case DataType::FLOAT128:   { REAL_LOOP(long double); }; break;
		case DataType::COMPLEX64:  { COMP_LOOP(float); }; break;
		case DataType::COMPLEX128: { COMP_LOOP(double); }; break;
		case DataType::COMPLEX256: { COMP_LOOP(long double); }; break;
		case DataType::RGB24: case DataType::RGBA32:
			throw(std::runtime_error("Unsupported datatype.")); break;
	}
	#undef COMP_LOOP
	#undef REAL_LOOP
	#undef DECL_PTR
}

template<typename T> void Nifti::convertToBytes(const typename std::vector<std::complex<T>>::iterator begin,
												const typename std::vector<std::complex<T>>::iterator end,
												std::vector<char> &bytes) {
	size_t nEl = std::distance(begin, end);
	assert(nEl == bytes.size() / m_typeinfo.size);
	T sc_sl = static_cast<T>(scaling_slope);
	T sc_in = static_cast<T>(scaling_inter);
	auto el = begin;
	#define DECL_PTR( BYTE_TYPE ) auto p = reinterpret_cast<BYTE_TYPE *>(bytes.data())
	#define REAL_LOOP( BYTE_TYPE ) DECL_PTR( BYTE_TYPE ); while (el != end) { *p = static_cast<BYTE_TYPE>(abs(*el) / sc_sl - sc_in); el++; p++; }
	#define COMP_LOOP( BYTE_TYPE ) DECL_PTR( std::complex<BYTE_TYPE> ); while (el != end) { *p = static_cast<std::complex<BYTE_TYPE>>(*el / sc_sl - sc_in); el++; p++; }
	switch (m_typeinfo.type) {
		case DataType::INT8:       { REAL_LOOP(int8_t); }; break;
		case DataType::INT16:      { REAL_LOOP(int16_t); }; break;
		case DataType::INT32:      { REAL_LOOP(int32_t); }; break;
		case DataType::INT64:      { REAL_LOOP(int64_t); }; break;
		case DataType::UINT8:      { REAL_LOOP(uint8_t); }; break;
		case DataType::UINT16:     { REAL_LOOP(uint16_t); }; break;
		case DataType::UINT32:     { REAL_LOOP(uint32_t); }; break;
		case DataType::UINT64:     { REAL_LOOP(uint64_t); }; break;
		case DataType::FLOAT32:    { REAL_LOOP(float); }; break;
		case DataType::FLOAT64:    { REAL_LOOP(double); }; break;
		case DataType::FLOAT128:   { REAL_LOOP(long double); }; break;
		case DataType::COMPLEX64:  { COMP_LOOP(float); }; break;
		case DataType::COMPLEX128: { COMP_LOOP(double); }; break;
		case DataType::COMPLEX256: { COMP_LOOP(long double); }; break;
		case DataType::RGB24: case DataType::RGBA32:
			throw(std::runtime_error("Unsupported datatype.")); break;
	}
	#undef COMP_LOOP
	#undef REAL_LOOP
	#undef DECL_PTR
}

/*
 *   Core IO routine. Depending on the mode of the file, reads/writes a
 *   contiguous subregion of the region and converts to/from the desired datatype.
 *
 *   start and size can be short of the full number of dimensions, e.g. if you
 *   only want a part of the first volume in a multi-volume image, then they can
 *   just be 3D. The only limitation is that they must not have more dimensions
 *   than the image.
 *
 *   @param start The voxel indices of the first desired voxel.
 *   @param size  The size of the desired subregion.
 *   @param data  Storage for the data to read/write. Must be sufficiently large.
 */
template<typename T> void Nifti::readWriteVoxels(const Eigen::Ref<ArrayXs> &start, const Eigen::Ref<ArrayXs> &size, std::vector<T> &data) {
	if (start.rows() != size.rows()) throw(std::out_of_range("Start and size must have same dimension in image: " + imagePath()));
{'fraction_non_alphanumeric': 0.1435522119551885, 'fraction_numerical': 0.016436339847902527, 'mean_word_length': 4.148179330667228, 'pattern_counts': {'":': 0, '<': 238, '<?xml version=': 0, '>': 234, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1165215', 'n_tokens_mistral': 8971, 'n_tokens_neox': 8329, 'n_words': 2234}