RNifti
Fast R and C++ Access to NIfTI Images
NiftiImage.h
1 #ifndef _NIFTI_IMAGE_H_
2 #define _NIFTI_IMAGE_H_
3 
4 
5 #ifdef USING_R
6 
7 #include <Rcpp.h>
8 
9 // Defined since R 3.1.0, according to Tomas Kalibera, but there's no reason to break compatibility with 3.0.x
10 #ifndef MAYBE_SHARED
11 #define MAYBE_SHARED(x) (NAMED(x) > 1)
12 #endif
13 
14 #else
15 
16 #define R_NegInf -INFINITY
17 
18 #include <stdint.h>
19 #include <cstddef>
20 #include <cmath>
21 #include <string>
22 #include <sstream>
23 #include <vector>
24 #include <stdexcept>
25 #include <algorithm>
26 #include <map>
27 #include <locale>
28 #include <limits>
29 
30 #endif
31 
32 
33 #include "niftilib/nifti1_io.h"
34 
43 namespace RNifti {
44 
45 namespace internal {
46 
47 struct vec3
48 {
49  float v[3];
50 
51  vec3 operator-() const
52  {
53  vec3 r;
54  r.v[0] = -v[0];
55  r.v[1] = -v[1];
56  r.v[2] = -v[2];
57  return r;
58  }
59 };
60 
61 inline mat33 topLeftCorner (const mat44 &matrix)
62 {
63  mat33 newMatrix;
64  for (int i=0; i<3; i++)
65  {
66  for (int j=0; j<3; j++)
67  newMatrix.m[i][j] = matrix.m[i][j];
68  }
69  return newMatrix;
70 }
71 
72 inline vec3 matrixVectorProduct (const mat33 &matrix, const vec3 &vector)
73 {
74  vec3 newVector;
75  for (int i=0; i<3; i++)
76  {
77  newVector.v[i] = 0.0;
78  for (int j=0; j<3; j++)
79  newVector.v[i] += matrix.m[i][j] * vector.v[j];
80  }
81  return newVector;
82 }
83 
84 } // internal namespace
85 
93 {
94 public:
99  struct Block
100  {
101  const NiftiImage &image;
102  const int dimension;
103  const int index;
112  Block (const NiftiImage &image, const int dimension, const int index)
114  {
115  if (dimension != image->ndim)
116  throw std::runtime_error("Blocks must be along the last dimension in the image");
117  }
118 
127  Block & operator= (const NiftiImage &source)
128  {
129  if (source->datatype != image->datatype)
130  throw std::runtime_error("New data does not have the same datatype as the target block");
131  if (source->scl_slope != image->scl_slope || source->scl_inter != image->scl_inter)
132  throw std::runtime_error("New data does not have the same scale parameters as the target block");
133 
134  size_t blockSize = 1;
135  for (int i=1; i<dimension; i++)
136  blockSize *= image->dim[i];
137 
138  if (blockSize != source->nvox)
139  throw std::runtime_error("New data does not have the same size as the target block");
140 
141  blockSize *= image->nbyper;
142  memcpy(static_cast<char*>(image->data) + blockSize*index, source->data, blockSize);
143  return *this;
144  }
145 
153  template <typename TargetType>
154  std::vector<TargetType> getData (const bool useSlope = true) const;
155  };
156 
157 #ifdef USING_R
158 
164  static short sexpTypeToNiftiType (const int sexpType)
165  {
166  if (sexpType == INTSXP || sexpType == LGLSXP)
167  return DT_INT32;
168  else if (sexpType == REALSXP)
169  return DT_FLOAT64;
170  else
171  throw std::runtime_error("Array elements must be numeric");
172  }
173 #endif
174 
180  static mat33 xformToRotation (const mat44 matrix)
181  {
182  float qb, qc, qd, qfac;
183  nifti_mat44_to_quatern(matrix, &qb, &qc, &qd, NULL, NULL, NULL, NULL, NULL, NULL, &qfac);
184  mat44 rotationMatrix = nifti_quatern_to_mat44(qb, qc, qd, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, qfac);
185  return internal::topLeftCorner(rotationMatrix);
186  }
187 
193  static std::string xformToString (const mat44 matrix)
194  {
195  int icode, jcode, kcode;
196  nifti_mat44_to_orientation(matrix, &icode, &jcode, &kcode);
197 
198  int codes[3] = { icode, jcode, kcode };
199  std::string result("---");
200  for (int i=0; i<3; i++)
201  {
202  switch (codes[i])
203  {
204  case NIFTI_L2R: result[i] = 'R'; break;
205  case NIFTI_R2L: result[i] = 'L'; break;
206  case NIFTI_P2A: result[i] = 'A'; break;
207  case NIFTI_A2P: result[i] = 'P'; break;
208  case NIFTI_I2S: result[i] = 'S'; break;
209  case NIFTI_S2I: result[i] = 'I'; break;
210  }
211  }
212  return result;
213  }
214 
221  static int fileVersion (const std::string &path)
222  {
223  nifti_1_header *header = nifti_read_header(path.c_str(), NULL, false);
224  if (header == NULL)
225  return -1;
226  else
227  {
228  int version = NIFTI_VERSION(*header);
229  if (version == 0)
230  {
231  // NIfTI-2 has a 540-byte header - check for this or its byte-swapped equivalent
232  if (header->sizeof_hdr == 540 || header->sizeof_hdr == 469893120)
233  {
234  // The magic number has moved in NIfTI-2, so find it by byte offset
235  const char *magic = (char *) header + 4;
236  if (strncmp(magic,"ni2",3) == 0 || strncmp(magic,"n+2",3) == 0)
237  version = 2;
238  }
239  else if (!nifti_hdr_looks_good(header))
240  {
241  // Not plausible as ANALYZE, so return -1
242  version = -1;
243  }
244  }
245  free(header);
246  return version;
247  }
248  }
249 
250 
251 protected:
252  nifti_image *image;
253  int *refCount;
261  void acquire (nifti_image * const image);
262 
267  void acquire (const NiftiImage &source)
268  {
269  refCount = source.refCount;
270  acquire(source.image);
271  }
272 
277  void release ();
278 
283  void copy (const nifti_image *source);
284 
289  void copy (const NiftiImage &source);
290 
295  void copy (const Block &source);
296 
297 
298 #ifdef USING_R
299 
305  void initFromNiftiS4 (const Rcpp::RObject &object, const bool copyData = true);
306 
312  void initFromMriImage (const Rcpp::RObject &object, const bool copyData = true);
313 
318  void initFromList (const Rcpp::RObject &object);
319 
325  void initFromArray (const Rcpp::RObject &object, const bool copyData = true);
326 
327 #endif
328 
333  void updatePixdim (const std::vector<float> &pixdim);
334 
339  void setPixunits (const std::vector<std::string> &pixunits);
340 
341 public:
346  : image(NULL), refCount(NULL) {}
347 
354  NiftiImage (const NiftiImage &source, const bool copy = true)
355  : image(NULL), refCount(NULL)
356  {
357  if (copy)
358  this->copy(source);
359  else
360  acquire(source);
361 #ifndef NDEBUG
362  Rc_printf("Creating NiftiImage with pointer %p (from NiftiImage)\n", this->image);
363 #endif
364  }
365 
370  NiftiImage (const Block &source)
371  : image(NULL), refCount(NULL)
372  {
373  this->copy(source);
374 #ifndef NDEBUG
375  Rc_printf("Creating NiftiImage with pointer %p (from Block)\n", this->image);
376 #endif
377  }
378 
385  NiftiImage (nifti_image * const image, const bool copy = false)
386  : image(NULL), refCount(NULL)
387  {
388  if (copy)
389  this->copy(image);
390  else
391  acquire(image);
392 #ifndef NDEBUG
393  Rc_printf("Creating NiftiImage with pointer %p (from pointer)\n", this->image);
394 #endif
395  }
396 
403  NiftiImage (const std::string &path, const bool readData = true)
404  : image(NULL), refCount(NULL)
405  {
406  acquire(nifti_image_read(path.c_str(), readData));
407  if (image == NULL)
408  throw std::runtime_error("Failed to read image from path " + path);
409 #ifndef NDEBUG
410  Rc_printf("Creating NiftiImage with pointer %p (from string)\n", this->image);
411 #endif
412  }
413 
420  NiftiImage (const std::string &path, const std::vector<int> &volumes);
421 
422 #ifdef USING_R
423 
433  NiftiImage (const SEXP object, const bool readData = true, const bool readOnly = false);
434 #endif
435 
440  virtual ~NiftiImage () { release(); }
441 
445  operator const nifti_image* () const { return image; }
446 
450  operator nifti_image* () { return image; }
451 
455  const nifti_image * operator-> () const { return image; }
456 
460  nifti_image * operator-> () { return image; }
461 
467  {
468  copy(source);
469 #ifndef NDEBUG
470  Rc_printf("Creating NiftiImage with pointer %p (from NiftiImage)\n", this->image);
471 #endif
472  return *this;
473  }
474 
480  NiftiImage & operator= (const Block &source)
481  {
482  copy(source);
483 #ifndef NDEBUG
484  Rc_printf("Creating NiftiImage with pointer %p (from Block)\n", this->image);
485 #endif
486  return *this;
487  }
488 
496  NiftiImage & setPersistence (const bool persistent) { return *this; }
497 
502  bool isNull () const { return (image == NULL); }
503 
508  bool isShared () const { return (refCount != NULL && *refCount > 1); }
509 
516  bool isPersistent () const { return false; }
517 
523  bool isDataScaled () const { return (image != NULL && image->scl_slope != 0.0 && (image->scl_slope != 1.0 || image->scl_inter != 0.0)); }
524 
529  int nDims () const
530  {
531  if (image == NULL)
532  return 0;
533  else
534  return image->ndim;
535  }
536 
541  std::vector<int> dim () const
542  {
543  if (image == NULL)
544  return std::vector<int>();
545  else
546  return std::vector<int>(image->dim+1, image->dim+image->ndim+1);
547  }
548 
553  std::vector<float> pixdim () const
554  {
555  if (image == NULL)
556  return std::vector<float>();
557  else
558  return std::vector<float>(image->pixdim+1, image->pixdim+image->ndim+1);
559  }
560 
568  {
569  int ndim = image->ndim;
570  while (image->dim[ndim] < 2)
571  ndim--;
572  image->dim[0] = image->ndim = ndim;
573 
574  return *this;
575  }
576 
585  template <typename TargetType>
586  std::vector<TargetType> getData (const bool useSlope = true) const;
587 
595  NiftiImage & changeDatatype (const short datatype, const bool useSlope = false);
596 
604  NiftiImage & changeDatatype (const std::string &datatype, const bool useSlope = false);
605 
614  template <typename SourceType>
615  NiftiImage & replaceData (const std::vector<SourceType> &data, const short datatype = DT_NONE);
616 
622  {
623  nifti_image_unload(image);
624  return *this;
625  }
626 
633  NiftiImage & rescale (const std::vector<float> &scales);
634 
643  NiftiImage & reorient (const int i, const int j, const int k);
644 
653  NiftiImage & reorient (const std::string &orientation);
654 
655 #ifdef USING_R
656 
661  NiftiImage & update (const Rcpp::RObject &object);
662 #endif
663 
669  mat44 xform (const bool preferQuaternion = true) const;
670 
675  int nBlocks () const
676  {
677  if (image == NULL)
678  return 0;
679  else
680  return image->dim[image->ndim];
681  }
682 
690  const Block block (const int i) const { return Block(*this, nDims(), i); }
691 
699  Block block (const int i) { return Block(*this, nDims(), i); }
700 
706  const Block slice (const int i) const { return Block(*this, 3, i); }
707 
713  Block slice (const int i) { return Block(*this, 3, i); }
714 
720  const Block volume (const int i) const { return Block(*this, 4, i); }
721 
727  Block volume (const int i) { return Block(*this, 4, i); }
728 
734  void toFile (const std::string fileName, const short datatype = DT_NONE) const;
735 
741  void toFile (const std::string fileName, const std::string &datatype) const;
742 
743 #ifdef USING_R
744 
749  Rcpp::RObject toArray () const;
750 
756  Rcpp::RObject toPointer (const std::string label) const;
757 
764  Rcpp::RObject toArrayOrPointer (const bool internal, const std::string label) const;
765 
770  Rcpp::RObject headerToList () const;
771 
772 #endif
773 
774 };
775 
776 // Include helper functions
777 #include "lib/NiftiImage_internal.h"
778 
779 inline void NiftiImage::acquire (nifti_image * const image)
780 {
781  // If we're taking ownership of a new image, release the old one
782  if (this->image != NULL && this->image != image)
783  release();
784 
785  // Set the internal pointer and create or update the reference counter
786  this->image = image;
787  if (image != NULL)
788  {
789  if (this->refCount == NULL)
790  this->refCount = new int(1);
791  else
792  (*this->refCount)++;
793 
794 #ifndef NDEBUG
795  Rc_printf("Acquiring pointer %p (reference count is %d)\n", this->image, *this->refCount);
796 #endif
797  }
798 }
799 
800 inline void NiftiImage::release ()
801 {
802  if (this->image != NULL)
803  {
804  if (this->refCount != NULL)
805  {
806  (*this->refCount)--;
807 #ifndef NDEBUG
808  Rc_printf("Releasing pointer %p (reference count is %d)\n", this->image, *this->refCount);
809 #endif
810  if (*this->refCount < 1)
811  {
812  nifti_image_free(this->image);
813  this->image = NULL;
814  delete this->refCount;
815  this->refCount = NULL;
816  }
817  }
818  else
819  Rc_printf("Releasing untracked object %p", this->image);
820  }
821 }
822 
823 inline void NiftiImage::copy (const nifti_image *source)
824 {
825  if (source == NULL)
826  acquire(NULL);
827  else
828  {
829  acquire(nifti_copy_nim_info(source));
830  if (source->data != NULL)
831  {
832  size_t dataSize = nifti_get_volsize(source);
833  image->data = calloc(1, dataSize);
834  memcpy(image->data, source->data, dataSize);
835  }
836  }
837 }
838 
839 inline void NiftiImage::copy (const NiftiImage &source)
840 {
841  const nifti_image *sourceStruct = source;
842  copy(sourceStruct);
843 }
844 
845 inline void NiftiImage::copy (const Block &source)
846 {
847  const nifti_image *sourceStruct = source.image;
848  if (sourceStruct == NULL)
849  acquire(NULL);
850  else
851  {
852  acquire(nifti_copy_nim_info(sourceStruct));
853  image->dim[0] = source.image->dim[0] - 1;
854  image->dim[source.dimension] = 1;
855  image->pixdim[source.dimension] = 1.0;
856  nifti_update_dims_from_array(image);
857 
858  if (sourceStruct->data != NULL)
859  {
860  size_t blockSize = nifti_get_volsize(image);
861  image->data = calloc(1, blockSize);
862  memcpy(image->data, static_cast<char*>(source.image->data) + blockSize*source.index, blockSize);
863  }
864  }
865 }
866 
867 #ifdef USING_R
868 
869 // Convert an S4 "nifti" object, as defined in the oro.nifti package, to a "nifti_image" struct
870 inline void NiftiImage::initFromNiftiS4 (const Rcpp::RObject &object, const bool copyData)
871 {
872  nifti_1_header header;
873  header.sizeof_hdr = 348;
874 
875  const std::vector<short> dims = object.slot("dim_");
876  for (int i=0; i<8; i++)
877  header.dim[i] = dims[i];
878 
879  header.intent_p1 = object.slot("intent_p1");
880  header.intent_p2 = object.slot("intent_p2");
881  header.intent_p3 = object.slot("intent_p3");
882  header.intent_code = object.slot("intent_code");
883 
884  header.datatype = object.slot("datatype");
885  header.bitpix = object.slot("bitpix");
886 
887  header.slice_start = object.slot("slice_start");
888  header.slice_end = object.slot("slice_end");
889  header.slice_code = Rcpp::as<int>(object.slot("slice_code"));
890  header.slice_duration = object.slot("slice_duration");
891 
892  const std::vector<float> pixdims = object.slot("pixdim");
893  for (int i=0; i<8; i++)
894  header.pixdim[i] = pixdims[i];
895  header.xyzt_units = Rcpp::as<int>(object.slot("xyzt_units"));
896 
897  header.vox_offset = object.slot("vox_offset");
898 
899  header.scl_slope = object.slot("scl_slope");
900  header.scl_inter = object.slot("scl_inter");
901  header.toffset = object.slot("toffset");
902 
903  header.cal_max = object.slot("cal_max");
904  header.cal_min = object.slot("cal_min");
905  header.glmax = header.glmin = 0;
906 
907  strncpy(header.descrip, Rcpp::as<std::string>(object.slot("descrip")).c_str(), 79);
908  header.descrip[79] = '\0';
909  strncpy(header.aux_file, Rcpp::as<std::string>(object.slot("aux_file")).c_str(), 23);
910  header.aux_file[23] = '\0';
911  strncpy(header.intent_name, Rcpp::as<std::string>(object.slot("intent_name")).c_str(), 15);
912  header.intent_name[15] = '\0';
913  strncpy(header.magic, Rcpp::as<std::string>(object.slot("magic")).c_str(), 3);
914  header.magic[3] = '\0';
915 
916  header.qform_code = object.slot("qform_code");
917  header.sform_code = object.slot("sform_code");
918 
919  header.quatern_b = object.slot("quatern_b");
920  header.quatern_c = object.slot("quatern_c");
921  header.quatern_d = object.slot("quatern_d");
922  header.qoffset_x = object.slot("qoffset_x");
923  header.qoffset_y = object.slot("qoffset_y");
924  header.qoffset_z = object.slot("qoffset_z");
925 
926  const std::vector<float> srow_x = object.slot("srow_x");
927  const std::vector<float> srow_y = object.slot("srow_y");
928  const std::vector<float> srow_z = object.slot("srow_z");
929  for (int i=0; i<4; i++)
930  {
931  header.srow_x[i] = srow_x[i];
932  header.srow_y[i] = srow_y[i];
933  header.srow_z[i] = srow_z[i];
934  }
935 
936  if (header.datatype == DT_UINT8 || header.datatype == DT_INT16 || header.datatype == DT_INT32 || header.datatype == DT_INT8 || header.datatype == DT_UINT16 || header.datatype == DT_UINT32)
937  header.datatype = DT_INT32;
938  else if (header.datatype == DT_FLOAT32 || header.datatype == DT_FLOAT64)
939  header.datatype = DT_FLOAT64; // This assumes that sizeof(double) == 8
940  else
941  throw std::runtime_error("Data type is not supported");
942 
943  acquire(nifti_convert_nhdr2nim(header, NULL));
944 
945  const SEXP data = PROTECT(object.slot(".Data"));
946  if (!copyData || Rf_length(data) <= 1)
947  this->image->data = NULL;
948  else
949  {
950  const size_t dataSize = nifti_get_volsize(this->image);
951  this->image->data = calloc(1, dataSize);
952  if (header.datatype == DT_INT32)
953  {
954  Rcpp::IntegerVector intData(data);
955  std::copy(intData.begin(), intData.end(), static_cast<int32_t*>(this->image->data));
956  }
957  else
958  {
959  Rcpp::DoubleVector doubleData(data);
960  std::copy(doubleData.begin(), doubleData.end(), static_cast<double*>(this->image->data));
961  }
962  }
963  UNPROTECT(1);
964 }
965 
966 inline void NiftiImage::initFromMriImage (const Rcpp::RObject &object, const bool copyData)
967 {
968  Rcpp::Reference mriImage(object);
969  Rcpp::Function getXform = mriImage.field("getXform");
970  Rcpp::NumericMatrix xform = getXform();
971 
972  acquire(NULL);
973 
974  if (Rf_length(mriImage.field("tags")) > 0)
975  initFromList(mriImage.field("tags"));
976 
977  Rcpp::RObject data = mriImage.field("data");
978  if (data.inherits("SparseArray"))
979  {
980  Rcpp::Language call("as.array", data);
981  data = call.eval();
982  }
983 
984  const short datatype = (Rf_isNull(data) ? DT_INT32 : sexpTypeToNiftiType(data.sexp_type()));
985 
986  int dims[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
987  const std::vector<int> dimVector = mriImage.field("imageDims");
988  const int nDims = std::min(7, int(dimVector.size()));
989  dims[0] = nDims;
990  size_t nVoxels = 1;
991  for (int i=0; i<nDims; i++)
992  {
993  dims[i+1] = dimVector[i];
994  nVoxels *= dimVector[i];
995  }
996 
997  if (this->image == NULL)
998  acquire(nifti_make_new_nim(dims, datatype, FALSE));
999  else
1000  {
1001  std::copy(dims, dims+8, this->image->dim);
1002  this->image->datatype = datatype;
1003  nifti_datatype_sizes(image->datatype, &image->nbyper, NULL);
1004  }
1005 
1006  if (copyData && !Rf_isNull(data))
1007  {
1008  // NB: nifti_get_volsize() will not be right here if there were tags
1009  const size_t dataSize = nVoxels * image->nbyper;
1010  this->image->data = calloc(1, dataSize);
1011  if (datatype == DT_INT32)
1012  memcpy(this->image->data, INTEGER(data), dataSize);
1013  else
1014  memcpy(this->image->data, REAL(data), dataSize);
1015  }
1016  else
1017  this->image->data = NULL;
1018 
1019  const std::vector<float> pixdimVector = mriImage.field("voxelDims");
1020  const int pixdimLength = pixdimVector.size();
1021  for (int i=0; i<std::min(pixdimLength,nDims); i++)
1022  this->image->pixdim[i+1] = std::abs(pixdimVector[i]);
1023 
1024  const std::vector<std::string> pixunitsVector = mriImage.field("voxelDimUnits");
1025  setPixunits(pixunitsVector);
1026 
1027  if (xform.rows() != 4 || xform.cols() != 4)
1028  this->image->qform_code = this->image->sform_code = 0;
1029  else
1030  {
1031  mat44 matrix;
1032  for (int i=0; i<4; i++)
1033  {
1034  for (int j=0; j<4; j++)
1035  matrix.m[i][j] = static_cast<float>(xform(i,j));
1036  }
1037 
1038  this->image->qto_xyz = matrix;
1039  this->image->qto_ijk = nifti_mat44_inverse(image->qto_xyz);
1040  nifti_mat44_to_quatern(image->qto_xyz, &image->quatern_b, &image->quatern_c, &image->quatern_d, &image->qoffset_x, &image->qoffset_y, &image->qoffset_z, NULL, NULL, NULL, &image->qfac);
1041 
1042  this->image->sto_xyz = matrix;
1043  this->image->sto_ijk = nifti_mat44_inverse(image->sto_xyz);
1044 
1045  this->image->qform_code = this->image->sform_code = 2;
1046  }
1047 }
1048 
1049 inline void NiftiImage::initFromList (const Rcpp::RObject &object)
1050 {
1051  Rcpp::List list(object);
1052  nifti_1_header *header = nifti_make_new_header(NULL, DT_FLOAT64);
1053 
1054  internal::updateHeader(header, list);
1055 
1056  acquire(nifti_convert_nhdr2nim(*header, NULL));
1057  this->image->data = NULL;
1058  free(header);
1059 }
1060 
1061 inline void NiftiImage::initFromArray (const Rcpp::RObject &object, const bool copyData)
1062 {
1063  int dims[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
1064  const std::vector<int> dimVector = object.attr("dim");
1065 
1066  const int nDims = std::min(7, int(dimVector.size()));
1067  dims[0] = nDims;
1068  for (int i=0; i<nDims; i++)
1069  dims[i+1] = dimVector[i];
1070 
1071  const short datatype = sexpTypeToNiftiType(object.sexp_type());
1072  acquire(nifti_make_new_nim(dims, datatype, int(copyData)));
1073 
1074  if (copyData)
1075  {
1076  const size_t dataSize = nifti_get_volsize(image);
1077  if (datatype == DT_INT32)
1078  memcpy(this->image->data, INTEGER(object), dataSize);
1079  else
1080  memcpy(this->image->data, REAL(object), dataSize);
1081  }
1082  else
1083  this->image->data = NULL;
1084 
1085  if (object.hasAttribute("pixdim"))
1086  {
1087  const std::vector<float> pixdimVector = object.attr("pixdim");
1088  const int pixdimLength = pixdimVector.size();
1089  for (int i=0; i<std::min(pixdimLength,nDims); i++)
1090  this->image->pixdim[i+1] = pixdimVector[i];
1091  }
1092 
1093  if (object.hasAttribute("pixunits"))
1094  {
1095  const std::vector<std::string> pixunitsVector = object.attr("pixunits");
1096  setPixunits(pixunitsVector);
1097  }
1098 }
1099 
1100 inline NiftiImage::NiftiImage (const SEXP object, const bool readData, const bool readOnly)
1101  : image(NULL), refCount(NULL)
1102 {
1103  Rcpp::RObject imageObject(object);
1104  bool resolved = false;
1105 
1106  if (imageObject.hasAttribute(".nifti_image_ptr"))
1107  {
1108  Rcpp::XPtr<NiftiImage> imagePtr(SEXP(imageObject.attr(".nifti_image_ptr")));
1109  NiftiImage *ptr = imagePtr;
1110  if (ptr != NULL)
1111  {
1112  if (MAYBE_SHARED(object) && !readOnly)
1113  copy(*ptr);
1114  else
1115  acquire(*ptr);
1116  resolved = true;
1117 
1118  if (imageObject.hasAttribute("dim"))
1119  update(imageObject);
1120  }
1121  else if (Rf_isString(object))
1122  throw std::runtime_error("Internal image is not valid");
1123  else
1124  Rf_warning("Ignoring invalid internal pointer");
1125  }
1126 
1127  if (!resolved)
1128  {
1129  if (Rf_isNull(object))
1130  acquire(NULL);
1131  else if (Rf_isString(object))
1132  {
1133  const std::string path = Rcpp::as<std::string>(object);
1134  acquire(nifti_image_read(path.c_str(), readData));
1135  if (this->image == NULL)
1136  throw std::runtime_error("Failed to read image from path " + path);
1137  }
1138  else if (imageObject.inherits("nifti"))
1139  initFromNiftiS4(imageObject, readData);
1140  else if (imageObject.inherits("anlz"))
1141  throw std::runtime_error("Cannot currently convert objects of class \"anlz\"");
1142  else if (imageObject.inherits("MriImage"))
1143  initFromMriImage(imageObject, readData);
1144  else if (Rf_isVectorList(object))
1145  initFromList(imageObject);
1146  else if (imageObject.hasAttribute("dim"))
1147  initFromArray(imageObject, readData);
1148  else if (imageObject.hasAttribute("class"))
1149  throw std::runtime_error("Cannot convert object of class \"" + Rcpp::as<std::string>(imageObject.attr("class")) + "\" to a nifti_image");
1150  else
1151  throw std::runtime_error("Cannot convert unclassed non-array object");
1152  }
1153 
1154  if (this->image != NULL)
1155  nifti_update_dims_from_array(this->image);
1156 
1157 #ifndef NDEBUG
1158  Rc_printf("Creating NiftiImage with pointer %p (from SEXP)\n", this->image);
1159 #endif
1160 }
1161 
1162 #endif // USING_R
1163 
1164 inline NiftiImage::NiftiImage (const std::string &path, const std::vector<int> &volumes)
1165  : image(NULL), refCount(NULL)
1166 {
1167  if (volumes.empty())
1168  throw std::runtime_error("The vector of volumes is empty");
1169 
1170  nifti_brick_list brickList;
1171  acquire(nifti_image_read_bricks(path.c_str(), volumes.size(), &volumes[0], &brickList));
1172  if (image == NULL)
1173  throw std::runtime_error("Failed to read image from path " + path);
1174 
1175  size_t brickSize = image->nbyper * image->nx * image->ny * image->nz;
1176  image->data = calloc(1, nifti_get_volsize(image));
1177  for (int i=0; i<brickList.nbricks; i++)
1178  memcpy((char *) image->data + i * brickSize, brickList.bricks[i], brickSize);
1179  nifti_free_NBL(&brickList);
1180 
1181 #ifndef NDEBUG
1182  Rc_printf("Creating NiftiImage with pointer %p (from string and volume vector)\n", this->image);
1183 #endif
1184 }
1185 
1186 inline void NiftiImage::updatePixdim (const std::vector<float> &pixdim)
1187 {
1188  const int nDims = image->dim[0];
1189  const std::vector<float> origPixdim(image->pixdim+1, image->pixdim+4);
1190 
1191  for (int i=1; i<8; i++)
1192  image->pixdim[i] = 0.0;
1193 
1194  const int pixdimLength = pixdim.size();
1195  for (int i=0; i<std::min(pixdimLength,nDims); i++)
1196  image->pixdim[i+1] = pixdim[i];
1197 
1198  if (!std::equal(origPixdim.begin(), origPixdim.begin() + std::min(3,nDims), pixdim.begin()))
1199  {
1200  mat33 scaleMatrix;
1201  for (int i=0; i<3; i++)
1202  {
1203  for (int j=0; j<3; j++)
1204  {
1205  if (i != j)
1206  scaleMatrix.m[i][j] = 0.0;
1207  else if (i >= nDims)
1208  scaleMatrix.m[i][j] = 1.0;
1209  else
1210  scaleMatrix.m[i][j] = pixdim[i] / origPixdim[i];
1211  }
1212  }
1213 
1214  if (image->qform_code > 0)
1215  {
1216  mat33 prod = nifti_mat33_mul(scaleMatrix, internal::topLeftCorner(image->qto_xyz));
1217  for (int i=0; i<3; i++)
1218  {
1219  for (int j=0; j<3; j++)
1220  image->qto_xyz.m[i][j] = prod.m[i][j];
1221  }
1222  image->qto_ijk = nifti_mat44_inverse(image->qto_xyz);
1223  nifti_mat44_to_quatern(image->qto_xyz, &image->quatern_b, &image->quatern_c, &image->quatern_d, &image->qoffset_x, &image->qoffset_y, &image->qoffset_z, NULL, NULL, NULL, &image->qfac);
1224  }
1225 
1226  if (image->sform_code > 0)
1227  {
1228  mat33 prod = nifti_mat33_mul(scaleMatrix, internal::topLeftCorner(image->sto_xyz));
1229  for (int i=0; i<3; i++)
1230  {
1231  for (int j=0; j<3; j++)
1232  image->sto_xyz.m[i][j] = prod.m[i][j];
1233  }
1234  image->sto_ijk = nifti_mat44_inverse(image->sto_xyz);
1235  }
1236  }
1237 }
1238 
1239 inline void NiftiImage::setPixunits (const std::vector<std::string> &pixunits)
1240 {
1241  for (size_t i=0; i<pixunits.size(); i++)
1242  {
1243  if (pixunits[i] == "m")
1244  image->xyz_units = NIFTI_UNITS_METER;
1245  else if (pixunits[i] == "mm")
1246  image->xyz_units = NIFTI_UNITS_MM;
1247  else if (pixunits[i] == "um")
1248  image->xyz_units = NIFTI_UNITS_MICRON;
1249  else if (pixunits[i] == "s")
1250  image->time_units = NIFTI_UNITS_SEC;
1251  else if (pixunits[i] == "ms")
1252  image->time_units = NIFTI_UNITS_MSEC;
1253  else if (pixunits[i] == "us")
1254  image->time_units = NIFTI_UNITS_USEC;
1255  else if (pixunits[i] == "Hz")
1256  image->time_units = NIFTI_UNITS_HZ;
1257  else if (pixunits[i] == "ppm")
1258  image->time_units = NIFTI_UNITS_PPM;
1259  else if (pixunits[i] == "rad/s")
1260  image->time_units = NIFTI_UNITS_RADS;
1261  }
1262 }
1263 
1264 inline NiftiImage & NiftiImage::rescale (const std::vector<float> &scales)
1265 {
1266  std::vector<float> pixdim(image->pixdim+1, image->pixdim+4);
1267 
1268  for (int i=0; i<std::min(3, int(scales.size())); i++)
1269  {
1270  if (scales[i] != 1.0)
1271  {
1272  pixdim[i] /= scales[i];
1273  image->dim[i+1] = static_cast<int>(std::floor(image->dim[i+1] * scales[i]));
1274  }
1275  }
1276 
1278  nifti_update_dims_from_array(image);
1279 
1280  // Data vector is now the wrong size, so drop it
1281  nifti_image_unload(image);
1282 
1283  image->scl_slope = 0.0;
1284  image->scl_inter = 0.0;
1285 
1286  return *this;
1287 }
1288 
1289 inline NiftiImage & NiftiImage::reorient (const int icode, const int jcode, const int kcode)
1290 {
1291  if (this->isNull())
1292  return *this;
1293  if (image->qform_code == 0 && image->sform_code == 0)
1294  {
1295  Rf_warning("Image qform and sform codes are both zero, so it cannot be reoriented");
1296  return *this;
1297  }
1298 
1299  int used[6] = { 0, 0, 0, 0, 0, 0 };
1300  used[icode-1] = 1;
1301  used[jcode-1] = 1;
1302  used[kcode-1] = 1;
1303  if (used[0]+used[1] != 1 || used[2]+used[3] != 1 || used[4]+used[5] != 1)
1304  throw std::runtime_error("Each canonical axis should be used exactly once");
1305 
1306  const int codes[3] = { icode, jcode, kcode };
1307  const mat44 native = this->xform();
1308 
1309  // Calculate the origin, which requires inverting the current xform
1310  // Here we use a simplified formula that exploits blockwise inversion and the nature of xforms
1311  internal::vec3 origin;
1312  for (int i=0; i<3; i++)
1313  origin.v[i] = native.m[i][3];
1314  origin = -internal::matrixVectorProduct(nifti_mat33_inverse(internal::topLeftCorner(native)), origin);
1315 
1316  // Create a target xform (rotation matrix only)
1317  mat33 target;
1318  for (int j=0; j<3; j++)
1319  {
1320  for (int i=0; i<3; i++)
1321  target.m[i][j] = 0.0;
1322 
1323  switch (codes[j])
1324  {
1325  case NIFTI_L2R: target.m[0][j] = 1.0; break;
1326  case NIFTI_R2L: target.m[0][j] = -1.0; break;
1327  case NIFTI_P2A: target.m[1][j] = 1.0; break;
1328  case NIFTI_A2P: target.m[1][j] = -1.0; break;
1329  case NIFTI_I2S: target.m[2][j] = 1.0; break;
1330  case NIFTI_S2I: target.m[2][j] = -1.0; break;
1331  }
1332  }
1333 
1334  // Extract (inverse of) canonical axis matrix from native xform
1335  int nicode, njcode, nkcode;
1336  nifti_mat44_to_orientation(native, &nicode, &njcode, &nkcode);
1337  int ncodes[3] = { nicode, njcode, nkcode };
1338  mat33 nativeAxesTransposed;
1339  for (int i=0; i<3; i++)
1340  {
1341  for (int j=0; j<3; j++)
1342  nativeAxesTransposed.m[i][j] = 0.0;
1343 
1344  switch (ncodes[i])
1345  {
1346  case NIFTI_L2R: nativeAxesTransposed.m[i][0] = 1.0; break;
1347  case NIFTI_R2L: nativeAxesTransposed.m[i][0] = -1.0; break;
1348  case NIFTI_P2A: nativeAxesTransposed.m[i][1] = 1.0; break;
1349  case NIFTI_A2P: nativeAxesTransposed.m[i][1] = -1.0; break;
1350  case NIFTI_I2S: nativeAxesTransposed.m[i][2] = 1.0; break;
1351  case NIFTI_S2I: nativeAxesTransposed.m[i][2] = -1.0; break;
1352  }
1353  }
1354 
1355  // Check for no-op case
1356  if (icode == nicode && jcode == njcode && kcode == nkcode)
1357  return *this;
1358 
1359  // The transform is t(approx_old_xform) %*% target_xform
1360  // The new xform is old_xform %*% transform
1361  // NB: "transform" is really 4x4, but the last row is simple and the last column is filled below
1362  mat33 transform = nifti_mat33_mul(nativeAxesTransposed, target);
1363  mat44 result;
1364  for (int i=0; i<4; i++)
1365  {
1366  for (int j=0; j<3; j++)
1367  result.m[i][j] = native.m[i][0] * transform.m[0][j] + native.m[i][1] * transform.m[1][j] + native.m[i][2] * transform.m[2][j];
1368 
1369  result.m[3][i] = i == 3 ? 1.0 : 0.0;
1370  }
1371 
1372  // Extract the mapping between dimensions and the signs
1373  // These vectors are all indexed in the target space, except "revsigns"
1374  int locs[3], signs[3], newdim[3], revsigns[3];
1375  float newpixdim[3];
1376  double maxes[3] = { R_NegInf, R_NegInf, R_NegInf };
1377  internal::vec3 offset;
1378  for (int j=0; j<3; j++)
1379  {
1380  // Find the largest absolute value in each column, which gives the old dimension corresponding to each new dimension
1381  for (int i=0; i<3; i++)
1382  {
1383  const double value = static_cast<double>(transform.m[i][j]);
1384  if (fabs(value) > maxes[j])
1385  {
1386  maxes[j] = fabs(value);
1387  signs[j] = value > 0.0 ? 1 : -1;
1388  locs[j] = i;
1389  }
1390  }
1391 
1392  // Obtain the sign for the reverse mapping
1393  revsigns[locs[j]] = signs[j];
1394 
1395  // Permute dim and pixdim
1396  newdim[j] = image->dim[locs[j]+1];
1397  newpixdim[j] = image->pixdim[locs[j]+1];
1398 
1399  // Flip and/or permute the origin
1400  if (signs[j] < 0)
1401  offset.v[j] = image->dim[locs[j]+1] - origin.v[locs[j]] - 1.0;
1402  else
1403  offset.v[j] = origin.v[locs[j]];
1404  }
1405 
1406  // Convert the origin back to an xform offset and insert it
1407  offset = -internal::matrixVectorProduct(internal::topLeftCorner(result), offset);
1408  for (int i=0; i<3; i++)
1409  result.m[i][3] = offset.v[i];
1410 
1411  // Update the xforms with nonzero codes
1412  if (image->qform_code > 0)
1413  {
1414  image->qto_xyz = result;
1415  image->qto_ijk = nifti_mat44_inverse(image->qto_xyz);
1416  nifti_mat44_to_quatern(image->qto_xyz, &image->quatern_b, &image->quatern_c, &image->quatern_d, &image->qoffset_x, &image->qoffset_y, &image->qoffset_z, NULL, NULL, NULL, &image->qfac);
1417  }
1418  if (image->sform_code > 0)
1419  {
1420  image->sto_xyz = result;
1421  image->sto_ijk = nifti_mat44_inverse(image->sto_xyz);
1422  }
1423 
1424  // Calculate strides: the step in target space associated with each dimension in source space
1425  ptrdiff_t strides[3];
1426  strides[locs[0]] = 1;
1427  strides[locs[1]] = strides[locs[0]] * image->dim[locs[0]+1];
1428  strides[locs[2]] = strides[locs[1]] * image->dim[locs[1]+1];
1429 
1430  // Permute the data (if present)
1431  if (image->data != NULL)
1432  {
1433  size_t volSize = size_t(image->nx * image->ny * image->nz);
1434  size_t nVolumes = std::max(size_t(1), image->nvox / volSize);
1435 
1436  const std::vector<double> oldData = this->getData<double>();
1437  std::vector<double> newData(image->nvox);
1438 
1439  // Where the sign is negative we need to start at the end of the dimension
1440  size_t volStart = 0;
1441  for (int i=0; i<3; i++)
1442  {
1443  if (revsigns[i] < 0)
1444  volStart += (image->dim[i+1] - 1) * strides[i];
1445  }
1446 
1447  // Iterate over the data and place it into a new vector
1448  std::vector<double>::const_iterator it = oldData.begin();
1449  for (size_t v=0; v<nVolumes; v++)
1450  {
1451  for (int k=0; k<image->nz; k++)
1452  {
1453  ptrdiff_t offset = k * strides[2] * revsigns[2];
1454  for (int j=0; j<image->ny; j++)
1455  {
1456  for (int i=0; i<image->nx; i++)
1457  {
1458  newData[volStart + offset] = *it++;
1459  offset += strides[0] * revsigns[0];
1460  }
1461  offset += strides[1] * revsigns[1] - image->nx * strides[0] * revsigns[0];
1462  }
1463  }
1464  volStart += volSize;
1465  }
1466 
1467  // Replace the existing data in the image
1468  this->replaceData(newData);
1469  }
1470 
1471  // Copy new dims and pixdims in
1472  // NB: Old dims are used above, so this must happen last
1473  std::copy(newdim, newdim+3, image->dim+1);
1474  std::copy(newpixdim, newpixdim+3, image->pixdim+1);
1475  nifti_update_dims_from_array(image);
1476 
1477  return *this;
1478 }
1479 
1480 inline NiftiImage & NiftiImage::reorient (const std::string &orientation)
1481 {
1482  if (orientation.length() != 3)
1483  throw std::runtime_error("Orientation string should have exactly three characters");
1484 
1485  int codes[3];
1486  for (int i=0; i<3; i++)
1487  {
1488  switch(orientation[i])
1489  {
1490  case 'r': case 'R': codes[i] = NIFTI_L2R; break;
1491  case 'l': case 'L': codes[i] = NIFTI_R2L; break;
1492  case 'a': case 'A': codes[i] = NIFTI_P2A; break;
1493  case 'p': case 'P': codes[i] = NIFTI_A2P; break;
1494  case 's': case 'S': codes[i] = NIFTI_I2S; break;
1495  case 'i': case 'I': codes[i] = NIFTI_S2I; break;
1496 
1497  default:
1498  throw std::runtime_error("Orientation string is invalid");
1499  }
1500  }
1501 
1502  return reorient(codes[0], codes[1], codes[2]);
1503 }
1504 
1505 #ifdef USING_R
1506 
1507 inline NiftiImage & NiftiImage::update (const Rcpp::RObject &object)
1508 {
1509  if (Rf_isVectorList(object))
1510  {
1511  Rcpp::List list(object);
1512  nifti_1_header *header = NULL;
1513  if (this->isNull())
1514  {
1515  header = nifti_make_new_header(NULL, DT_FLOAT64);
1516  internal::updateHeader(header, list, true);
1517  }
1518  else
1519  {
1520  header = (nifti_1_header *) calloc(1, sizeof(nifti_1_header));
1521  *header = nifti_convert_nim2nhdr(image);
1522  internal::updateHeader(header, list, true);
1523  }
1524 
1525  if (header != NULL)
1526  {
1527  // Retain the data pointer, but otherwise overwrite the stored object with one created from the header
1528  // The file names can't be preserved through the round-trip, so free them
1529  void *dataPtr = image->data;
1530  nifti_image *tempImage = nifti_convert_nhdr2nim(*header, NULL);
1531 
1532  if (image->fname != NULL)
1533  free(image->fname);
1534  if (image->iname != NULL)
1535  free(image->iname);
1536 
1537  memcpy(image, tempImage, sizeof(nifti_image));
1538  image->num_ext = 0;
1539  image->ext_list = NULL;
1540  image->data = dataPtr;
1541 
1542  nifti_image_free(tempImage);
1543  free(header);
1544  }
1545  }
1546  else if (object.hasAttribute("dim"))
1547  {
1548  for (int i=0; i<8; i++)
1549  image->dim[i] = 0;
1550  const std::vector<int> dimVector = object.attr("dim");
1551 
1552  const int nDims = std::min(7, int(dimVector.size()));
1553  image->dim[0] = nDims;
1554  for (int i=0; i<nDims; i++)
1555  image->dim[i+1] = dimVector[i];
1556 
1557  if (object.hasAttribute("pixdim"))
1558  {
1559  const std::vector<float> pixdimVector = object.attr("pixdim");
1560  updatePixdim(pixdimVector);
1561  }
1562 
1563  if (object.hasAttribute("pixunits"))
1564  {
1565  const std::vector<std::string> pixunitsVector = object.attr("pixunits");
1566  setPixunits(pixunitsVector);
1567  }
1568 
1569  // This NIfTI-1 library function clobbers dim[0] if the last dimension is unitary; we undo that here
1570  nifti_update_dims_from_array(image);
1571  image->dim[0] = image->ndim = nDims;
1572 
1573  image->datatype = NiftiImage::sexpTypeToNiftiType(object.sexp_type());
1574  nifti_datatype_sizes(image->datatype, &image->nbyper, NULL);
1575 
1576  nifti_image_unload(image);
1577 
1578  const size_t dataSize = nifti_get_volsize(image);
1579  image->data = calloc(1, dataSize);
1580  if (image->datatype == DT_INT32)
1581  memcpy(image->data, INTEGER(object), dataSize);
1582  else
1583  memcpy(image->data, REAL(object), dataSize);
1584 
1585  image->scl_slope = 0.0;
1586  image->scl_inter = 0.0;
1587  }
1588 
1589  return *this;
1590 }
1591 
1592 #endif // USING_R
1593 
1594 inline mat44 NiftiImage::xform (const bool preferQuaternion) const
1595 {
1596  if (image == NULL)
1597  {
1598  mat44 matrix;
1599  for (int i=0; i<4; i++)
1600  {
1601  for (int j=0; j<4; j++)
1602  matrix.m[i][j] = 0.0;
1603  }
1604  return matrix;
1605  }
1606  else if (image->qform_code <= 0 && image->sform_code <= 0)
1607  {
1608  // No qform or sform so use pixdim (NB: other software may assume differently)
1609  mat44 matrix;
1610  for (int i=0; i<4; i++)
1611  {
1612  for (int j=0; j<4; j++)
1613  {
1614  if (i != j)
1615  matrix.m[i][j] = 0.0;
1616  else if (i == 3)
1617  matrix.m[3][3] = 1.0;
1618  else
1619  matrix.m[i][j] = (image->pixdim[i+1]==0.0 ? 1.0 : image->pixdim[i+1]);
1620  }
1621  }
1622  return matrix;
1623  }
1624  else if ((preferQuaternion && image->qform_code > 0) || image->sform_code <= 0)
1625  return image->qto_xyz;
1626  else
1627  return image->sto_xyz;
1628 }
1629 
1630 template <typename TargetType>
1631 inline std::vector<TargetType> NiftiImage::Block::getData (const bool useSlope) const
1632 {
1633  if (image.isNull())
1634  return std::vector<TargetType>();
1635 
1636  size_t blockSize = 1;
1637  for (int i=1; i<dimension; i++)
1638  blockSize *= image->dim[i];
1639 
1640  std::vector<TargetType> data(blockSize);
1641  internal::DataConverter<TargetType> *converter = NULL;
1642  if (useSlope && image.isDataScaled())
1643  converter = new internal::DataConverter<TargetType>(image->scl_slope, image->scl_inter);
1644 
1645  internal::convertData<TargetType>(image->data, image->datatype, blockSize, data.begin(), blockSize*index, converter);
1646 
1647  delete converter;
1648  return data;
1649 }
1650 
1651 template <typename TargetType>
1652 inline std::vector<TargetType> NiftiImage::getData (const bool useSlope) const
1653 {
1654  if (this->isNull())
1655  return std::vector<TargetType>();
1656 
1657  std::vector<TargetType> data(image->nvox);
1658  internal::DataConverter<TargetType> *converter = NULL;
1659  if (useSlope && this->isDataScaled())
1660  converter = new internal::DataConverter<TargetType>(image->scl_slope, image->scl_inter);
1661 
1662  internal::convertData<TargetType>(image->data, image->datatype, image->nvox, data.begin(), 0, converter);
1663 
1664  delete converter;
1665  return data;
1666 }
1667 
1668 inline NiftiImage & NiftiImage::changeDatatype (const short datatype, const bool useSlope)
1669 {
1670  if (this->isNull() || image->datatype == datatype)
1671  return *this;
1672 
1673  if (useSlope && this->isDataScaled())
1674  throw std::runtime_error("Resetting the slope and intercept for an image with them already set is not supported");
1675 
1676  if (image->data != NULL)
1677  {
1678  int bytesPerPixel;
1679  nifti_datatype_sizes(datatype, &bytesPerPixel, NULL);
1680  void *data = calloc(image->nvox, bytesPerPixel);
1681 
1682  switch (datatype)
1683  {
1684  case DT_UINT8:
1685  {
1686  internal::DataConverter<uint8_t> converter(useSlope ? internal::DataConverter<uint8_t>::IndexMode : internal::DataConverter<uint8_t>::CastMode);
1687  internal::convertData(image->data, image->datatype, image->nvox, static_cast<uint8_t *>(data), 0, &converter);
1688  image->scl_slope = static_cast<float>(converter.getSlope());
1689  image->scl_inter = static_cast<float>(converter.getIntercept());
1690  break;
1691  }
1692 
1693  case DT_INT16:
1694  {
1695  internal::DataConverter<int16_t> converter(useSlope ? internal::DataConverter<int16_t>::IndexMode : internal::DataConverter<int16_t>::CastMode);
1696  internal::convertData(image->data, image->datatype, image->nvox, static_cast<int16_t *>(data), 0, &converter);
1697  image->scl_slope = static_cast<float>(converter.getSlope());
1698  image->scl_inter = static_cast<float>(converter.getIntercept());
1699  break;
1700  }
1701 
1702  case DT_INT32:
1703  {
1704  internal::DataConverter<int32_t> converter(useSlope ? internal::DataConverter<int32_t>::IndexMode : internal::DataConverter<int32_t>::CastMode);
1705  internal::convertData(image->data, image->datatype, image->nvox, static_cast<int32_t *>(data), 0, &converter);
1706  image->scl_slope = static_cast<float>(converter.getSlope());
1707  image->scl_inter = static_cast<float>(converter.getIntercept());
1708  break;
1709  }
1710 
1711  case DT_FLOAT32:
1712  {
1713  internal::DataConverter<float> converter(useSlope ? internal::DataConverter<float>::IndexMode : internal::DataConverter<float>::CastMode);
1714  internal::convertData(image->data, image->datatype, image->nvox, static_cast<float *>(data), 0, &converter);
1715  image->scl_slope = static_cast<float>(converter.getSlope());
1716  image->scl_inter = static_cast<float>(converter.getIntercept());
1717  break;
1718  }
1719 
1720  case DT_FLOAT64:
1721  {
1722  internal::DataConverter<double> converter(useSlope ? internal::DataConverter<double>::IndexMode : internal::DataConverter<double>::CastMode);
1723  internal::convertData(image->data, image->datatype, image->nvox, static_cast<double *>(data), 0, &converter);
1724  image->scl_slope = static_cast<float>(converter.getSlope());
1725  image->scl_inter = static_cast<float>(converter.getIntercept());
1726  break;
1727  }
1728 
1729  case DT_INT8:
1730  {
1731  internal::DataConverter<int8_t> converter(useSlope ? internal::DataConverter<int8_t>::IndexMode : internal::DataConverter<int8_t>::CastMode);
1732  internal::convertData(image->data, image->datatype, image->nvox, static_cast<int8_t *>(data), 0, &converter);
1733  image->scl_slope = static_cast<float>(converter.getSlope());
1734  image->scl_inter = static_cast<float>(converter.getIntercept());
1735  break;
1736  }
1737 
1738  case DT_UINT16:
1739  {
1740  internal::DataConverter<uint16_t> converter(useSlope ? internal::DataConverter<uint16_t>::IndexMode : internal::DataConverter<uint16_t>::CastMode);
1741  internal::convertData(image->data, image->datatype, image->nvox, static_cast<uint16_t *>(data), 0, &converter);
1742  image->scl_slope = static_cast<float>(converter.getSlope());
1743  image->scl_inter = static_cast<float>(converter.getIntercept());
1744  break;
1745  }
1746 
1747  case DT_UINT32:
1748  {
1749  internal::DataConverter<uint32_t> converter(useSlope ? internal::DataConverter<uint32_t>::IndexMode : internal::DataConverter<uint32_t>::CastMode);
1750  internal::convertData(image->data, image->datatype, image->nvox, static_cast<uint32_t *>(data), 0, &converter);
1751  image->scl_slope = static_cast<float>(converter.getSlope());
1752  image->scl_inter = static_cast<float>(converter.getIntercept());
1753  break;
1754  }
1755 
1756  case DT_INT64:
1757  {
1758  internal::DataConverter<int64_t> converter(useSlope ? internal::DataConverter<int64_t>::IndexMode : internal::DataConverter<int64_t>::CastMode);
1759  internal::convertData(image->data, image->datatype, image->nvox, static_cast<int64_t *>(data), 0, &converter);
1760  image->scl_slope = static_cast<float>(converter.getSlope());
1761  image->scl_inter = static_cast<float>(converter.getIntercept());
1762  break;
1763  }
1764 
1765  case DT_UINT64:
1766  {
1767  internal::DataConverter<uint64_t> converter(useSlope ? internal::DataConverter<uint64_t>::IndexMode : internal::DataConverter<uint64_t>::CastMode);
1768  internal::convertData(image->data, image->datatype, image->nvox, static_cast<uint64_t *>(data), 0, &converter);
1769  image->scl_slope = static_cast<float>(converter.getSlope());
1770  image->scl_inter = static_cast<float>(converter.getIntercept());
1771  break;
1772  }
1773 
1774  default:
1775  throw std::runtime_error("Unsupported data type (" + std::string(nifti_datatype_string(datatype)) + ")");
1776  }
1777 
1778  nifti_image_unload(image);
1779  image->data = data;
1780  }
1781 
1782  image->datatype = datatype;
1783  nifti_datatype_sizes(datatype, &image->nbyper, &image->swapsize);
1784 
1785  return *this;
1786 }
1787 
1788 inline NiftiImage & NiftiImage::changeDatatype (const std::string &datatype, const bool useSlope)
1789 {
1790  return changeDatatype(internal::stringToDatatype(datatype), useSlope);
1791 }
1792 
1793 template <typename SourceType>
1794 inline NiftiImage & NiftiImage::replaceData (const std::vector<SourceType> &data, const short datatype)
1795 {
1796  if (this->isNull())
1797  return *this;
1798  else if (data.size() != image->nvox)
1799  throw std::runtime_error("New data length does not match the number of voxels in the image");
1800 
1801  if (datatype != DT_NONE)
1802  {
1803  nifti_image_unload(image);
1804  image->datatype = datatype;
1805  nifti_datatype_sizes(datatype, &image->nbyper, &image->swapsize);
1806  }
1807 
1808  if (image->data == NULL)
1809  image->data = calloc(image->nvox, image->nbyper);
1810  internal::replaceData<SourceType>(data.begin(), data.end(), image->data, image->datatype);
1811 
1812  image->scl_slope = 0.0;
1813  image->scl_inter = 0.0;
1814  image->cal_min = static_cast<float>(*std::min_element(data.begin(), data.end()));
1815  image->cal_max = static_cast<float>(*std::max_element(data.begin(), data.end()));
1816 
1817  return *this;
1818 }
1819 
1820 inline void NiftiImage::toFile (const std::string fileName, const short datatype) const
1821 {
1822  const bool changingDatatype = (datatype != DT_NONE && !this->isNull() && datatype != image->datatype);
1823 
1824  // Copy the source image only if the datatype will be changed
1825  NiftiImage imageToWrite(*this, changingDatatype);
1826 
1827  if (changingDatatype)
1828  imageToWrite.changeDatatype(datatype, true);
1829 
1830  const int status = nifti_set_filenames(imageToWrite, fileName.c_str(), false, true);
1831  if (status != 0)
1832  throw std::runtime_error("Failed to set filenames for NIfTI object");
1833  nifti_image_write(imageToWrite);
1834 }
1835 
1836 inline void NiftiImage::toFile (const std::string fileName, const std::string &datatype) const
1837 {
1838  toFile(fileName, internal::stringToDatatype(datatype));
1839 }
1840 
1841 #ifdef USING_R
1842 
1843 inline Rcpp::RObject NiftiImage::toArray () const
1844 {
1845  Rcpp::RObject array;
1846 
1847  if (this->isNull())
1848  return array;
1849  else if (this->isDataScaled())
1850  {
1851  array = internal::imageDataToArray<REALSXP>(image);
1852  std::transform(REAL(array), REAL(array)+Rf_length(array), REAL(array), internal::DataConverter<double>(image->scl_slope,image->scl_inter));
1853  }
1854  else
1855  {
1856  switch (image->datatype)
1857  {
1858  case DT_UINT8:
1859  case DT_INT16:
1860  case DT_INT32:
1861  case DT_INT8:
1862  case DT_UINT16:
1863  case DT_UINT32:
1864  case DT_INT64:
1865  case DT_UINT64:
1866  array = internal::imageDataToArray<INTSXP>(image);
1867  break;
1868 
1869  case DT_FLOAT32:
1870  case DT_FLOAT64:
1871  array = internal::imageDataToArray<REALSXP>(image);
1872  break;
1873 
1874  default:
1875  throw std::runtime_error("Unsupported data type (" + std::string(nifti_datatype_string(image->datatype)) + ")");
1876  }
1877  }
1878 
1879  internal::addAttributes(array, *this);
1880  array.attr("class") = Rcpp::CharacterVector::create("niftiImage", "array");
1881  return array;
1882 }
1883 
1884 inline Rcpp::RObject NiftiImage::toPointer (const std::string label) const
1885 {
1886  if (this->isNull())
1887  return Rcpp::RObject();
1888  else
1889  {
1890  Rcpp::RObject string = Rcpp::wrap(label);
1891  internal::addAttributes(string, *this, false);
1892  string.attr("class") = Rcpp::CharacterVector::create("internalImage", "niftiImage");
1893  return string;
1894  }
1895 }
1896 
1897 inline Rcpp::RObject NiftiImage::toArrayOrPointer (const bool internal, const std::string label) const
1898 {
1899  return (internal ? toPointer(label) : toArray());
1900 }
1901 
1902 inline Rcpp::RObject NiftiImage::headerToList () const
1903 {
1904  if (this->image == NULL)
1905  return Rcpp::RObject();
1906 
1907  nifti_1_header header = nifti_convert_nim2nhdr(this->image);
1908  Rcpp::List result;
1909 
1910  result["sizeof_hdr"] = header.sizeof_hdr;
1911 
1912  result["dim_info"] = int(header.dim_info);
1913  result["dim"] = std::vector<short>(header.dim, header.dim+8);
1914 
1915  result["intent_p1"] = header.intent_p1;
1916  result["intent_p2"] = header.intent_p2;
1917  result["intent_p3"] = header.intent_p3;
1918  result["intent_code"] = header.intent_code;
1919 
1920  result["datatype"] = header.datatype;
1921  result["bitpix"] = header.bitpix;
1922 
1923  result["slice_start"] = header.slice_start;
1924  result["pixdim"] = std::vector<float>(header.pixdim, header.pixdim+8);
1925  result["vox_offset"] = header.vox_offset;
1926  result["scl_slope"] = header.scl_slope;
1927  result["scl_inter"] = header.scl_inter;
1928  result["slice_end"] = header.slice_end;
1929  result["slice_code"] = int(header.slice_code);
1930  result["xyzt_units"] = int(header.xyzt_units);
1931  result["cal_max"] = header.cal_max;
1932  result["cal_min"] = header.cal_min;
1933  result["slice_duration"] = header.slice_duration;
1934  result["toffset"] = header.toffset;
1935  result["descrip"] = std::string(header.descrip, 80);
1936  result["aux_file"] = std::string(header.aux_file, 24);
1937 
1938  result["qform_code"] = header.qform_code;
1939  result["sform_code"] = header.sform_code;
1940  result["quatern_b"] = header.quatern_b;
1941  result["quatern_c"] = header.quatern_c;
1942  result["quatern_d"] = header.quatern_d;
1943  result["qoffset_x"] = header.qoffset_x;
1944  result["qoffset_y"] = header.qoffset_y;
1945  result["qoffset_z"] = header.qoffset_z;
1946  result["srow_x"] = std::vector<float>(header.srow_x, header.srow_x+4);
1947  result["srow_y"] = std::vector<float>(header.srow_y, header.srow_y+4);
1948  result["srow_z"] = std::vector<float>(header.srow_z, header.srow_z+4);
1949 
1950  result["intent_name"] = std::string(header.intent_name, 16);
1951  result["magic"] = std::string(header.magic, 4);
1952 
1953  internal::addAttributes(result, *this, false, false);
1954  result.attr("class") = Rcpp::CharacterVector::create("niftiHeader");
1955 
1956  return result;
1957 }
1958 
1959 #endif // USING_R
1960 
1961 } // namespace
1962 
1963 #endif
NiftiImage & reorient(const int i, const int j, const int k)
Reorient the image by permuting dimensions and potentially reversing some.
Definition: NiftiImage.h:1289
NiftiImage()
Default constructor.
Definition: NiftiImage.h:345
const int index
The location along dimension.
Definition: NiftiImage.h:103
static std::string xformToString(const mat44 matrix)
Convert a 4x4 xform matrix to a string describing its canonical axes.
Definition: NiftiImage.h:193
std::vector< int > dim() const
Return the dimensions of the image.
Definition: NiftiImage.h:541
const nifti_image * operator->() const
Allows a NiftiImage object to be treated as a pointer to a const nifti_image.
Definition: NiftiImage.h:455
int * refCount
A reference counter, shared with other objects wrapping the same pointer.
Definition: NiftiImage.h:253
std::vector< float > pixdim() const
Return the dimensions of the pixels or voxels in the image.
Definition: NiftiImage.h:553
NiftiImage(const std::string &path, const bool readData=true)
Initialise using a path string.
Definition: NiftiImage.h:403
NiftiImage & rescale(const std::vector< float > &scales)
Rescale the image, changing its image dimensions and pixel dimensions.
Definition: NiftiImage.h:1264
NiftiImage(nifti_image *const image, const bool copy=false)
Initialise using an existing nifti_image pointer.
Definition: NiftiImage.h:385
Definition: NiftiImage.h:43
Block & operator=(const NiftiImage &source)
Copy assignment operator, which allows a block in one image to be replaced with the contents of anoth...
Definition: NiftiImage.h:127
bool isNull() const
Determine whether or not the wrapped pointer is NULL.
Definition: NiftiImage.h:502
Inner class referring to a subset of an image.
Definition: NiftiImage.h:99
Block slice(const int i)
Extract a slice block from a 3D image.
Definition: NiftiImage.h:713
bool isShared() const
Determine whether the wrapped pointer is shared with another NiftiImage.
Definition: NiftiImage.h:508
NiftiImage & replaceData(const std::vector< SourceType > &data, const short datatype=DT_NONE)
Replace the pixel data in the image with the contents of a vector.
Definition: NiftiImage.h:1794
void acquire(const NiftiImage &source)
Acquire the same pointer as another NiftiImage, incrementing the shared reference count...
Definition: NiftiImage.h:267
void acquire(nifti_image *const image)
Acquire the specified pointer to a nifti_image struct, taking (possibly shared) responsibility for fr...
Definition: NiftiImage.h:779
bool isDataScaled() const
Determine whether nontrivial scale and slope parameters are set.
Definition: NiftiImage.h:523
const Block slice(const int i) const
Extract a slice block from a 3D image.
Definition: NiftiImage.h:706
const Block volume(const int i) const
Extract a volume block from a 4D image.
Definition: NiftiImage.h:720
bool isPersistent() const
Determine whether or not the image is marked as persistent.
Definition: NiftiImage.h:516
std::vector< TargetType > getData(const bool useSlope=true) const
Extract a vector of data from the image, casting it to any required element type. ...
Definition: NiftiImage.h:1652
std::vector< TargetType > getData(const bool useSlope=true) const
Extract a vector of data from a block, casting it to any required element type.
Definition: NiftiImage.h:1631
void setPixunits(const std::vector< std::string > &pixunits)
Modify the pixel dimension units.
Definition: NiftiImage.h:1239
NiftiImage(const Block &source)
Initialise from a block, copying in the data.
Definition: NiftiImage.h:370
virtual ~NiftiImage()
Destructor which decrements the reference counter, and releases the wrapped pointer if the counter dr...
Definition: NiftiImage.h:440
Thin wrapper around a C-style nifti_image struct that allows C++-style destruction.
Definition: NiftiImage.h:92
NiftiImage(const NiftiImage &source, const bool copy=true)
Copy constructor.
Definition: NiftiImage.h:354
int nBlocks() const
Return the number of blocks in the image.
Definition: NiftiImage.h:675
void release()
Release the currently wrapped pointer, if it is not NULL, decrementing the reference count and releas...
Definition: NiftiImage.h:800
const int dimension
The dimension along which the block applies (which should be the last)
Definition: NiftiImage.h:102
const NiftiImage & image
The parent image.
Definition: NiftiImage.h:101
NiftiImage & operator=(const NiftiImage &source)
Copy assignment operator, which copies from its argument.
Definition: NiftiImage.h:466
static int fileVersion(const std::string &path)
Get the NIfTI format version used by the file at the specified path.
Definition: NiftiImage.h:221
static mat33 xformToRotation(const mat44 matrix)
Extract the pure rotation part of a 4x4 xform matrix.
Definition: NiftiImage.h:180
void toFile(const std::string fileName, const short datatype=DT_NONE) const
Write the image to a NIfTI-1 file.
Definition: NiftiImage.h:1820
void updatePixdim(const std::vector< float > &pixdim)
Modify the pixel dimensions, and potentially the xform matrices to match.
Definition: NiftiImage.h:1186
const Block block(const int i) const
Extract a block from the image.
Definition: NiftiImage.h:690
NiftiImage & setPersistence(const bool persistent)
Mark the image as persistent, so that it can be passed back to R.
Definition: NiftiImage.h:496
Definition: NiftiImage.h:47
NiftiImage & changeDatatype(const short datatype, const bool useSlope=false)
Change the datatype of the image, casting the pixel data if present.
Definition: NiftiImage.h:1668
void copy(const nifti_image *source)
Copy the contents of a nifti_image to create a new image, acquiring the new pointer.
Definition: NiftiImage.h:823
mat44 xform(const bool preferQuaternion=true) const
Obtain an xform matrix, indicating the orientation of the image.
Definition: NiftiImage.h:1594
nifti_image * image
The wrapped nifti_image pointer.
Definition: NiftiImage.h:252
int nDims() const
Return the number of dimensions in the image.
Definition: NiftiImage.h:529
Block volume(const int i)
Extract a volume block from a 4D image.
Definition: NiftiImage.h:727
NiftiImage & drop()
Drop unitary dimensions.
Definition: NiftiImage.h:567
Block(const NiftiImage &image, const int dimension, const int index)
Standard constructor for this class.
Definition: NiftiImage.h:112
NiftiImage & dropData()
Drop the data from the image, retaining only the metadata.
Definition: NiftiImage.h:621
Block block(const int i)
Extract a block from the image.
Definition: NiftiImage.h:699