RNifti Fast R and C++ Access to NIfTI Images

NiftiImage.h
1 #ifndef _NIFTI_IMAGE_H_
2 #define _NIFTI_IMAGE_H_
3
4 #include <Rcpp.h>
5
6 #include "niftilib/nifti1_io.h"
7
16 namespace RNifti {
17
22 class NiftiImage
23 {
24 public:
29 struct Block
30 {
31 const NiftiImage &image;
32 const int dimension;
33 const int index;
42 Block (const NiftiImage &image, const int dimension, const int index)
43 : image(image), dimension(dimension), index(index)
44 {
45 if (dimension != image->ndim)
46 throw std::runtime_error("Blocks must be along the last dimension in the image");
47 }
48
57 Block & operator= (const NiftiImage &source)
58 {
59 if (source->datatype != image->datatype)
60 throw std::runtime_error("New data does not have the same datatype as the target block");
61
62 size_t blockSize = 1;
63 for (int i=1; i<dimension; i++)
64 blockSize *= image->dim[i];
65
66 if (blockSize != source->nvox)
67 throw std::runtime_error("New data does not have the same size as the target block");
68
69 blockSize *= image->nbyper;
70 memcpy(static_cast<char*>(image->data) + blockSize*index, source->data, blockSize);
71 return *this;
72 }
73
77 template <typename TargetType>
78 std::vector<TargetType> getData () const;
79 };
80
87 static short sexpTypeToNiftiType (const int sexpType)
88 {
89 if (sexpType == INTSXP || sexpType == LGLSXP)
90 return DT_INT32;
91 else if (sexpType == REALSXP)
92 return DT_FLOAT64;
93 else
94 throw std::runtime_error("Array elements must be numeric");
95 }
96
97 protected:
98 nifti_image *image;
99 bool persistent;
105 void copy (const nifti_image *source);
106
111 void copy (const NiftiImage &source);
112
117 void copy (const Block &source);
118
124 void initFromNiftiS4 (const Rcpp::RObject &object, const bool copyData = true);
125
131 void initFromMriImage (const Rcpp::RObject &object, const bool copyData = true);
132
137 void initFromList (const Rcpp::RObject &object);
138
144 void initFromArray (const Rcpp::RObject &object, const bool copyData = true);
145
150 void updatePixdim (const std::vector<float> &pixdim);
151
156 void setPixunits (const std::vector<std::string> &pixunits);
157
158 public:
162 NiftiImage ()
163 : image(NULL), persistent(false) {}
164
169 NiftiImage (const NiftiImage &source)
170 : image(NULL), persistent(false)
171 {
172 this->copy(source);
173 #ifndef NDEBUG
174 Rprintf("Creating NiftiImage with pointer %p (from NiftiImage)\n", this->image);
175 #endif
176 }
177
182 NiftiImage (const Block &source)
183 : image(NULL), persistent(false)
184 {
185 this->copy(source);
186 #ifndef NDEBUG
187 Rprintf("Creating NiftiImage with pointer %p (from Block)\n", this->image);
188 #endif
189 }
190
197 NiftiImage (nifti_image * const image, const bool copy = false)
198 : image(NULL), persistent(false)
199 {
200 if (copy)
201 this->copy(image);
202 else
203 this->image = image;
204 #ifndef NDEBUG
205 Rprintf("Creating NiftiImage with pointer %p (from pointer)\n", this->image);
206 #endif
207 }
208
215 NiftiImage (const std::string &path, const bool readData = true)
216 : persistent(false)
217 {
218 this->image = nifti_image_read(path.c_str(), readData);
219 if (this->image == NULL)
220 throw std::runtime_error("Failed to read image from path " + path);
221 #ifndef NDEBUG
222 Rprintf("Creating NiftiImage with pointer %p (from string)\n", this->image);
223 #endif
224 }
225
231 NiftiImage (const SEXP object, const bool readData = true);
232
236 virtual ~NiftiImage ()
237 {
238 if (!persistent)
239 {
240 #ifndef NDEBUG
241 Rprintf("Freeing NiftiImage with pointer %p\n", this->image);
242 #endif
243 nifti_image_free(image);
244 }
245 }
246
250 operator const nifti_image* () const { return image; }
251
255 operator nifti_image* () { return image; }
256
260 const nifti_image * operator-> () const { return image; }
261
265 nifti_image * operator-> () { return image; }
266
271 NiftiImage & operator= (const NiftiImage &source)
272 {
273 copy(source);
274 #ifndef NDEBUG
275 Rprintf("Creating NiftiImage with pointer %p (from NiftiImage)\n", this->image);
276 #endif
277 return *this;
278 }
279
285 NiftiImage & operator= (const Block &source)
286 {
287 copy(source);
288 #ifndef NDEBUG
289 Rprintf("Creating NiftiImage with pointer %p (from Block)\n", this->image);
290 #endif
291 return *this;
292 }
293
298 void setPersistence (const bool persistent)
299 {
300 this->persistent = persistent;
301 #ifndef NDEBUG
302 if (persistent)
303 Rprintf("Setting NiftiImage with pointer %p to be persistent\n", this->image);
304 #endif
305 }
306
310 bool isNull () const { return (image == NULL); }
311
315 bool isPersistent () const { return persistent; }
316
320 int nDims () const
321 {
322 if (image == NULL)
323 return 0;
324 else
325 return image->ndim;
326 }
327
332 std::vector<int> dim () const
333 {
334 if (image == NULL)
335 return std::vector<int>();
336 else
337 return std::vector<int>(image->dim+1, image->dim+image->ndim+1);
338 }
339
344 std::vector<float> pixdim () const
345 {
346 if (image == NULL)
347 return std::vector<float>();
348 else
349 return std::vector<float>(image->pixdim+1, image->pixdim+image->ndim+1);
350 }
351
358 NiftiImage & drop ()
359 {
360 int ndim = image->ndim;
361 while (image->dim[ndim] < 2)
362 ndim--;
363 image->dim[0] = image->ndim = ndim;
364
365 return *this;
366 }
367
371 template <typename TargetType>
372 std::vector<TargetType> getData () const;
373
378 void changeDatatype (const short datatype);
379
384 void changeDatatype (const std::string &datatype);
385
393 template <typename SourceType>
394 void replaceData (const std::vector<SourceType> &data, const short datatype = DT_NONE);
395
399 void dropData () { nifti_image_unload(image); }
400
406 void rescale (const std::vector<float> &scales);
407
412 void update (const SEXP array);
413
419 mat44 xform (const bool preferQuaternion = true) const;
420
424 int nBlocks () const
425 {
426 if (image == NULL)
427 return 0;
428 else
429 return image->dim[image->ndim];
430 }
431
439 const Block block (const int i) const { return Block(*this, nDims(), i); }
440
448 Block block (const int i) { return Block(*this, nDims(), i); }
449
455 const Block slice (const int i) const { return Block(*this, 3, i); }
456
462 Block slice (const int i) { return Block(*this, 3, i); }
463
469 const Block volume (const int i) const { return Block(*this, 4, i); }
470
476 Block volume (const int i) { return Block(*this, 4, i); }
477
483 void toFile (const std::string fileName, const short datatype = DT_NONE) const;
484
490 void toFile (const std::string fileName, const std::string &datatype) const;
491
496 Rcpp::RObject toArray () const;
497
503 Rcpp::RObject toPointer (const std::string label) const;
504
511 Rcpp::RObject toArrayOrPointer (const bool internal, const std::string label) const;
512
517 Rcpp::RObject headerToList () const;
518 };
519
520 // Include helper functions
521 #include "lib/NiftiImage_internal.h"
522
523 inline void NiftiImage::copy (const nifti_image *source)
524 {
525 if (image != NULL)
526 nifti_image_free(image);
527
528 if (source == NULL)
529 image = NULL;
530 else
531 {
532 image = nifti_copy_nim_info(source);
533 if (source->data != NULL)
534 {
535 size_t dataSize = nifti_get_volsize(source);
536 image->data = calloc(1, dataSize);
537 memcpy(image->data, source->data, dataSize);
538 }
539 }
540 }
541
542 inline void NiftiImage::copy (const NiftiImage &source)
543 {
544 const nifti_image *sourceStruct = source;
545 copy(sourceStruct);
546 }
547
548 inline void NiftiImage::copy (const Block &source)
549 {
550 if (image != NULL)
551 nifti_image_free(image);
552
553 const nifti_image *sourceStruct = source.image;
554 if (sourceStruct == NULL)
555 image = NULL;
556 else
557 {
558 image = nifti_copy_nim_info(sourceStruct);
559 image->dim[0] = source.image->dim[0] - 1;
560 image->dim[source.dimension] = 1;
561 image->pixdim[source.dimension] = 1.0;
562 nifti_update_dims_from_array(image);
563
564 if (sourceStruct->data != NULL)
565 {
566 size_t blockSize = nifti_get_volsize(image);
567 image->data = calloc(1, blockSize);
568 memcpy(image->data, static_cast<char*>(source.image->data) + blockSize*source.index, blockSize);
569 }
570 }
571 }
572
573 // Convert an S4 "nifti" object, as defined in the oro.nifti package, to a "nifti_image" struct
574 inline void NiftiImage::initFromNiftiS4 (const Rcpp::RObject &object, const bool copyData)
575 {
576 nifti_1_header header;
577 header.sizeof_hdr = 348;
578
579 const std::vector<short> dims = object.slot("dim_");
580 for (int i=0; i<8; i++)
581 header.dim[i] = dims[i];
582
583 header.intent_p1 = object.slot("intent_p1");
584 header.intent_p2 = object.slot("intent_p2");
585 header.intent_p3 = object.slot("intent_p3");
586 header.intent_code = object.slot("intent_code");
587
588 header.datatype = object.slot("datatype");
589 header.bitpix = object.slot("bitpix");
590
591 header.slice_start = object.slot("slice_start");
592 header.slice_end = object.slot("slice_end");
593 header.slice_code = Rcpp::as<int>(object.slot("slice_code"));
594 header.slice_duration = object.slot("slice_duration");
595
596 const std::vector<float> pixdims = object.slot("pixdim");
597 for (int i=0; i<8; i++)
598 header.pixdim[i] = pixdims[i];
599 header.xyzt_units = Rcpp::as<int>(object.slot("xyzt_units"));
600
601 header.vox_offset = object.slot("vox_offset");
602
603 header.scl_slope = object.slot("scl_slope");
604 header.scl_inter = object.slot("scl_inter");
605 header.toffset = object.slot("toffset");
606
607 header.cal_max = object.slot("cal_max");
608 header.cal_min = object.slot("cal_min");
609 header.glmax = header.glmin = 0;
610
611 strncpy(header.descrip, Rcpp::as<std::string>(object.slot("descrip")).c_str(), 79);
612 header.descrip[79] = '\0';
613 strncpy(header.aux_file, Rcpp::as<std::string>(object.slot("aux_file")).c_str(), 23);
614 header.aux_file[23] = '\0';
615 strncpy(header.intent_name, Rcpp::as<std::string>(object.slot("intent_name")).c_str(), 15);
616 header.intent_name[15] = '\0';
617 strncpy(header.magic, Rcpp::as<std::string>(object.slot("magic")).c_str(), 3);
618 header.magic[3] = '\0';
619
620 header.qform_code = object.slot("qform_code");
621 header.sform_code = object.slot("sform_code");
622
623 header.quatern_b = object.slot("quatern_b");
624 header.quatern_c = object.slot("quatern_c");
625 header.quatern_d = object.slot("quatern_d");
626 header.qoffset_x = object.slot("qoffset_x");
627 header.qoffset_y = object.slot("qoffset_y");
628 header.qoffset_z = object.slot("qoffset_z");
629
630 const std::vector<float> srow_x = object.slot("srow_x");
631 const std::vector<float> srow_y = object.slot("srow_y");
632 const std::vector<float> srow_z = object.slot("srow_z");
633 for (int i=0; i<4; i++)
634 {
635 header.srow_x[i] = srow_x[i];
636 header.srow_y[i] = srow_y[i];
637 header.srow_z[i] = srow_z[i];
638 }
639
640 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)
641 header.datatype = DT_INT32;
642 else if (header.datatype == DT_FLOAT32 || header.datatype == DT_FLOAT64)
643 header.datatype = DT_FLOAT64; // This assumes that sizeof(double) == 8
644 else
645 throw std::runtime_error("Data type is not supported");
646
647 this->image = nifti_convert_nhdr2nim(header, NULL);
648
649 const SEXP data = PROTECT(object.slot(".Data"));
650 if (!copyData || Rf_length(data) <= 1)
651 this->image->data = NULL;
652 else
653 {
654 const size_t dataSize = nifti_get_volsize(this->image);
655 this->image->data = calloc(1, dataSize);
656 if (header.datatype == DT_INT32)
657 {
658 Rcpp::IntegerVector intData(data);
659 std::copy(intData.begin(), intData.end(), static_cast<int32_t*>(this->image->data));
660 }
661 else
662 {
663 Rcpp::DoubleVector doubleData(data);
664 std::copy(doubleData.begin(), doubleData.end(), static_cast<double*>(this->image->data));
665 }
666 }
667 UNPROTECT(1);
668 }
669
670 inline void NiftiImage::initFromMriImage (const Rcpp::RObject &object, const bool copyData)
671 {
672 Rcpp::Reference mriImage(object);
673 Rcpp::Function getXform = mriImage.field("getXform");
674 Rcpp::NumericMatrix xform = getXform();
675
676 this->image = NULL;
677
678 if (Rf_length(mriImage.field("tags")) > 0)
679 initFromList(mriImage.field("tags"));
680
681 Rcpp::RObject data = mriImage.field("data");
682 if (data.inherits("SparseArray"))
683 {
684 Rcpp::Language call("as.array", data);
685 data = call.eval();
686 }
687
688 const short datatype = (Rf_isNull(data) ? DT_INT32 : sexpTypeToNiftiType(data.sexp_type()));
689
690 int dims[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
691 const std::vector<int> dimVector = mriImage.field("imageDims");
692 const int nDims = std::min(7, int(dimVector.size()));
693 dims[0] = nDims;
694 size_t nVoxels = 1;
695 for (int i=0; i<nDims; i++)
696 {
697 dims[i+1] = dimVector[i];
698 nVoxels *= dimVector[i];
699 }
700
701 if (this->image == NULL)
702 this->image = nifti_make_new_nim(dims, datatype, FALSE);
703 else
704 {
705 std::copy(dims, dims+8, this->image->dim);
706 this->image->datatype = datatype;
707 nifti_datatype_sizes(image->datatype, &image->nbyper, NULL);
708 }
709
710 if (copyData && !Rf_isNull(data))
711 {
712 // NB: nifti_get_volsize() will not be right here if there were tags
713 const size_t dataSize = nVoxels * image->nbyper;
714 this->image->data = calloc(1, dataSize);
715 if (datatype == DT_INT32)
716 memcpy(this->image->data, INTEGER(data), dataSize);
717 else
718 memcpy(this->image->data, REAL(data), dataSize);
719 }
720 else
721 this->image->data = NULL;
722
723 const std::vector<float> pixdimVector = mriImage.field("voxelDims");
724 const int pixdimLength = pixdimVector.size();
725 for (int i=0; i<std::min(pixdimLength,nDims); i++)
726 this->image->pixdim[i+1] = std::abs(pixdimVector[i]);
727
728 const std::vector<std::string> pixunitsVector = mriImage.field("voxelDimUnits");
729 setPixunits(pixunitsVector);
730
731 if (xform.rows() != 4 || xform.cols() != 4)
732 this->image->qform_code = this->image->sform_code = 0;
733 else
734 {
735 mat44 matrix;
736 for (int i=0; i<4; i++)
737 {
738 for (int j=0; j<4; j++)
739 matrix.m[i][j] = static_cast<float>(xform(i,j));
740 }
741
742 this->image->qto_xyz = matrix;
743 this->image->qto_ijk = nifti_mat44_inverse(image->qto_xyz);
744 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);
745
746 this->image->sto_xyz = matrix;
747 this->image->sto_ijk = nifti_mat44_inverse(image->sto_xyz);
748
749 this->image->qform_code = this->image->sform_code = 2;
750 }
751 }
752
753 inline void NiftiImage::initFromList (const Rcpp::RObject &object)
754 {
755 Rcpp::List list(object);
756 nifti_1_header *header = nifti_make_new_header(NULL, DT_FLOAT64);
757
758 const Rcpp::CharacterVector _names = list.names();
759 std::set<std::string> names;
760 for (Rcpp::CharacterVector::const_iterator it=_names.begin(); it!=_names.end(); it++)
761 names.insert(Rcpp::as<std::string>(*it));
762
763 internal::copyIfPresent(list, names, "sizeof_hdr", header->sizeof_hdr);
764
765 internal::copyIfPresent(list, names, "dim_info", header->dim_info);
766 if (names.count("dim") == 1)
767 {
768 std::vector<short> dim = list["dim"];
769 for (size_t i=0; i<std::min(dim.size(),size_t(8)); i++)
770 header->dim[i] = dim[i];
771 }
772
773 internal::copyIfPresent(list, names, "intent_p1", header->intent_p1);
774 internal::copyIfPresent(list, names, "intent_p2", header->intent_p2);
775 internal::copyIfPresent(list, names, "intent_p3", header->intent_p3);
776 internal::copyIfPresent(list, names, "intent_code", header->intent_code);
777
778 internal::copyIfPresent(list, names, "datatype", header->datatype);
779 internal::copyIfPresent(list, names, "bitpix", header->bitpix);
780
781 internal::copyIfPresent(list, names, "slice_start", header->slice_start);
782 if (names.count("pixdim") == 1)
783 {
784 std::vector<float> pixdim = list["pixdim"];
785 for (size_t i=0; i<std::min(pixdim.size(),size_t(8)); i++)
786 header->pixdim[i] = pixdim[i];
787 }
788 internal::copyIfPresent(list, names, "vox_offset", header->vox_offset);
789 internal::copyIfPresent(list, names, "scl_slope", header->scl_slope);
790 internal::copyIfPresent(list, names, "scl_inter", header->scl_inter);
791 internal::copyIfPresent(list, names, "slice_end", header->slice_end);
792 internal::copyIfPresent(list, names, "slice_code", header->slice_code);
793 internal::copyIfPresent(list, names, "xyzt_units", header->xyzt_units);
794 internal::copyIfPresent(list, names, "cal_max", header->cal_max);
795 internal::copyIfPresent(list, names, "cal_min", header->cal_min);
796 internal::copyIfPresent(list, names, "slice_duration", header->slice_duration);
797 internal::copyIfPresent(list, names, "toffset", header->toffset);
798
799 if (names.count("descrip") == 1)
800 strcpy(header->descrip, Rcpp::as<std::string>(list["descrip"]).substr(0,79).c_str());
801 if (names.count("aux_file") == 1)
802 strcpy(header->aux_file, Rcpp::as<std::string>(list["aux_file"]).substr(0,23).c_str());
803
804 internal::copyIfPresent(list, names, "qform_code", header->qform_code);
805 internal::copyIfPresent(list, names, "sform_code", header->sform_code);
806 internal::copyIfPresent(list, names, "quatern_b", header->quatern_b);
807 internal::copyIfPresent(list, names, "quatern_c", header->quatern_c);
808 internal::copyIfPresent(list, names, "quatern_d", header->quatern_d);
809 internal::copyIfPresent(list, names, "qoffset_x", header->qoffset_x);
810 internal::copyIfPresent(list, names, "qoffset_y", header->qoffset_y);
811 internal::copyIfPresent(list, names, "qoffset_z", header->qoffset_z);
812
813 if (names.count("srow_x") == 1)
814 {
815 std::vector<float> srow_x = list["srow_x"];
816 for (size_t i=0; i<std::min(srow_x.size(),size_t(4)); i++)
817 header->srow_x[i] = srow_x[i];
818 }
819 if (names.count("srow_y") == 1)
820 {
821 std::vector<float> srow_y = list["srow_y"];
822 for (size_t i=0; i<std::min(srow_y.size(),size_t(4)); i++)
823 header->srow_y[i] = srow_y[i];
824 }
825 if (names.count("srow_z") == 1)
826 {
827 std::vector<float> srow_z = list["srow_z"];
828 for (size_t i=0; i<std::min(srow_z.size(),size_t(4)); i++)
829 header->srow_z[i] = srow_z[i];
830 }
831
832 if (names.count("intent_name") == 1)
833 strcpy(header->intent_name, Rcpp::as<std::string>(list["intent_name"]).substr(0,15).c_str());
834 if (names.count("magic") == 1)
835 strcpy(header->magic, Rcpp::as<std::string>(list["magic"]).substr(0,3).c_str());
836
837 this->image = nifti_convert_nhdr2nim(*header, NULL);
838 this->image->data = NULL;
839 free(header);
840 }
841
842 inline void NiftiImage::initFromArray (const Rcpp::RObject &object, const bool copyData)
843 {
844 int dims[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
845 const std::vector<int> dimVector = object.attr("dim");
846
847 const int nDims = std::min(7, int(dimVector.size()));
848 dims[0] = nDims;
849 for (int i=0; i<nDims; i++)
850 dims[i+1] = dimVector[i];
851
852 const short datatype = sexpTypeToNiftiType(object.sexp_type());
853 this->image = nifti_make_new_nim(dims, datatype, int(copyData));
854
855 if (copyData)
856 {
857 const size_t dataSize = nifti_get_volsize(image);
858 if (datatype == DT_INT32)
859 memcpy(this->image->data, INTEGER(object), dataSize);
860 else
861 memcpy(this->image->data, REAL(object), dataSize);
862 }
863 else
864 this->image->data = NULL;
865
866 if (object.hasAttribute("pixdim"))
867 {
868 const std::vector<float> pixdimVector = object.attr("pixdim");
869 const int pixdimLength = pixdimVector.size();
870 for (int i=0; i<std::min(pixdimLength,nDims); i++)
871 this->image->pixdim[i+1] = pixdimVector[i];
872 }
873
874 if (object.hasAttribute("pixunits"))
875 {
876 const std::vector<std::string> pixunitsVector = object.attr("pixunits");
877 setPixunits(pixunitsVector);
878 }
879 }
880
881 inline NiftiImage::NiftiImage (const SEXP object, const bool readData)
882 : persistent(false)
883 {
884 Rcpp::RObject imageObject(object);
885 bool resolved = false;
886
887 if (imageObject.hasAttribute(".nifti_image_ptr"))
888 {
889 Rcpp::XPtr<NiftiImage> imagePtr(SEXP(imageObject.attr(".nifti_image_ptr")));
890 if (imagePtr.get() != NULL)
891 {
892 this->image = imagePtr->image;
893 this->persistent = true;
894 resolved = true;
895
896 if (imageObject.hasAttribute("dim"))
897 update(object);
898 }
899 else if (Rf_isString(object))
900 throw std::runtime_error("Internal image is not valid");
901 else
902 Rf_warning("Ignoring invalid internal pointer");
903 }
904
905 if (!resolved)
906 {
907 if (Rf_isNull(object))
908 this->image = NULL;
909 else if (Rf_isString(object))
910 {
911 const std::string path = Rcpp::as<std::string>(object);
912 this->image = nifti_image_read(path.c_str(), readData);
913 if (this->image == NULL)
914 throw std::runtime_error("Failed to read image from path " + path);
915 }
916 else if (imageObject.inherits("nifti"))
917 initFromNiftiS4(imageObject, readData);
918 else if (imageObject.inherits("anlz"))
919 throw std::runtime_error("Cannot currently convert objects of class \"anlz\"");
920 else if (imageObject.inherits("MriImage"))
921 initFromMriImage(imageObject, readData);
922 else if (Rf_isVectorList(object))
923 initFromList(imageObject);
924 else if (imageObject.hasAttribute("dim"))
925 initFromArray(imageObject, readData);
926 else
927 throw std::runtime_error("Cannot convert object of class \"" + Rcpp::as<std::string>(imageObject.attr("class")) + "\" to a nifti_image");
928 }
929
930 if (this->image != NULL)
931 nifti_update_dims_from_array(this->image);
932
933 #ifndef NDEBUG
934 Rprintf("Creating NiftiImage with pointer %p (from SEXP)\n", this->image);
935 #endif
936 }
937
938 inline void NiftiImage::updatePixdim (const std::vector<float> &pixdim)
939 {
940 const int nDims = image->dim[0];
941 const std::vector<float> origPixdim(image->pixdim+1, image->pixdim+4);
942
943 for (int i=1; i<8; i++)
944 image->pixdim[i] = 0.0;
945
946 const int pixdimLength = pixdim.size();
947 for (int i=0; i<std::min(pixdimLength,nDims); i++)
948 image->pixdim[i+1] = pixdim[i];
949
950 if (!std::equal(origPixdim.begin(), origPixdim.begin() + std::min(3,nDims), pixdim.begin()))
951 {
952 mat33 scaleMatrix;
953 for (int i=0; i<3; i++)
954 {
955 for (int j=0; j<3; j++)
956 {
957 if (i != j)
958 scaleMatrix.m[i][j] = 0.0;
959 else if (i >= nDims)
960 scaleMatrix.m[i][j] = 1.0;
961 else
962 scaleMatrix.m[i][j] = pixdim[i] / origPixdim[i];
963 }
964 }
965
966 if (image->qform_code > 0)
967 {
968 mat33 prod = nifti_mat33_mul(scaleMatrix, internal::topLeftCorner(image->qto_xyz));
969 for (int i=0; i<3; i++)
970 {
971 for (int j=0; j<3; j++)
972 image->qto_xyz.m[i][j] = prod.m[i][j];
973 }
974 image->qto_ijk = nifti_mat44_inverse(image->qto_xyz);
975 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);
976 }
977
978 if (image->sform_code > 0)
979 {
980 mat33 prod = nifti_mat33_mul(scaleMatrix, internal::topLeftCorner(image->sto_xyz));
981 for (int i=0; i<3; i++)
982 {
983 for (int j=0; j<3; j++)
984 image->sto_xyz.m[i][j] = prod.m[i][j];
985 }
986 image->sto_ijk = nifti_mat44_inverse(image->sto_xyz);
987 }
988 }
989 }
990
991 inline void NiftiImage::setPixunits (const std::vector<std::string> &pixunits)
992 {
993 for (size_t i=0; i<pixunits.size(); i++)
994 {
995 if (pixunits[i] == "m")
996 image->xyz_units = NIFTI_UNITS_METER;
997 else if (pixunits[i] == "mm")
998 image->xyz_units = NIFTI_UNITS_MM;
999 else if (pixunits[i] == "um")
1000 image->xyz_units = NIFTI_UNITS_MICRON;
1001 else if (pixunits[i] == "s")
1002 image->time_units = NIFTI_UNITS_SEC;
1003 else if (pixunits[i] == "ms")
1004 image->time_units = NIFTI_UNITS_MSEC;
1005 else if (pixunits[i] == "us")
1006 image->time_units = NIFTI_UNITS_USEC;
1007 else if (pixunits[i] == "Hz")
1008 image->time_units = NIFTI_UNITS_HZ;
1009 else if (pixunits[i] == "ppm")
1010 image->time_units = NIFTI_UNITS_PPM;
1011 else if (pixunits[i] == "rad/s")
1012 image->time_units = NIFTI_UNITS_RADS;
1013 }
1014 }
1015
1016 inline void NiftiImage::rescale (const std::vector<float> &scales)
1017 {
1018 std::vector<float> pixdim(image->pixdim+1, image->pixdim+4);
1019
1020 for (int i=0; i<std::min(3, int(scales.size())); i++)
1021 {
1022 if (scales[i] != 1.0)
1023 {
1024 pixdim[i] /= scales[i];
1025 image->dim[i+1] = static_cast<int>(std::floor(image->dim[i+1] * scales[i]));
1026 }
1027 }
1028
1029 updatePixdim(pixdim);
1030 nifti_update_dims_from_array(image);
1031
1032 // Data vector is now the wrong size, so drop it
1033 free(image->data);
1034 }
1035
1036 inline void NiftiImage::update (const SEXP array)
1037 {
1038 Rcpp::RObject object(array);
1039 if (!object.hasAttribute("dim"))
1040 return;
1041
1042 for (int i=0; i<8; i++)
1043 image->dim[i] = 0;
1044 const std::vector<int> dimVector = object.attr("dim");
1045
1046 const int nDims = std::min(7, int(dimVector.size()));
1047 image->dim[0] = nDims;
1048 for (int i=0; i<nDims; i++)
1049 image->dim[i+1] = dimVector[i];
1050
1051 if (object.hasAttribute("pixdim"))
1052 {
1053 const std::vector<float> pixdimVector = object.attr("pixdim");
1054 updatePixdim(pixdimVector);
1055 }
1056
1057 if (object.hasAttribute("pixunits"))
1058 {
1059 const std::vector<std::string> pixunitsVector = object.attr("pixunits");
1060 setPixunits(pixunitsVector);
1061 }
1062
1063 // This NIfTI-1 library function clobbers dim[0] if the last dimension is unitary; we undo that here
1064 nifti_update_dims_from_array(image);
1065 image->dim[0] = image->ndim = nDims;
1066
1067 image->datatype = NiftiImage::sexpTypeToNiftiType(object.sexp_type());
1068 nifti_datatype_sizes(image->datatype, &image->nbyper, NULL);
1069
1070 free(image->data);
1071
1072 const size_t dataSize = nifti_get_volsize(image);
1073 image->data = calloc(1, dataSize);
1074 if (image->datatype == DT_INT32)
1075 memcpy(image->data, INTEGER(object), dataSize);
1076 else
1077 memcpy(image->data, REAL(object), dataSize);
1078 }
1079
1080 inline mat44 NiftiImage::xform (const bool preferQuaternion) const
1081 {
1082 if (image == NULL)
1083 {
1084 mat44 matrix;
1085 for (int i=0; i<4; i++)
1086 {
1087 for (int j=0; j<4; j++)
1088 matrix.m[i][j] = 0.0;
1089 }
1090 return matrix;
1091 }
1092 else if (image->qform_code <= 0 && image->sform_code <= 0)
1093 {
1094 // No qform or sform so use pixdim (NB: other software may assume differently)
1095 mat44 matrix;
1096 for (int i=0; i<4; i++)
1097 {
1098 for (int j=0; j<4; j++)
1099 {
1100 if (i != j)
1101 matrix.m[i][j] = 0.0;
1102 else if (i == 3)
1103 matrix.m[3][3] = 1.0;
1104 else
1105 matrix.m[i][j] = (image->pixdim[i+1]==0.0 ? 1.0 : image->pixdim[i+1]);
1106 }
1107 }
1108 return matrix;
1109 }
1110 else if ((preferQuaternion && image->qform_code > 0) || image->sform_code <= 0)
1111 return image->qto_xyz;
1112 else
1113 return image->sto_xyz;
1114 }
1115
1116 template <typename TargetType>
1117 inline std::vector<TargetType> NiftiImage::Block::getData () const
1118 {
1119 if (image.isNull())
1120 return std::vector<TargetType>();
1121
1122 size_t blockSize = 1;
1123 for (int i=1; i<dimension; i++)
1124 blockSize *= image->dim[i];
1125
1126 std::vector<TargetType> data(blockSize);
1127 internal::convertData<TargetType>(image->data, image->datatype, blockSize, data.begin(), blockSize*index);
1128 return data;
1129 }
1130
1131 template <typename TargetType>
1132 inline std::vector<TargetType> NiftiImage::getData () const
1133 {
1134 if (this->isNull())
1135 return std::vector<TargetType>();
1136
1137 std::vector<TargetType> data(image->nvox);
1138 internal::convertData<TargetType>(image->data, image->datatype, image->nvox, data.begin());
1139 return data;
1140 }
1141
1142 inline void NiftiImage::changeDatatype (const short datatype)
1143 {
1144 if (this->isNull() || image->datatype == datatype)
1145 return;
1146
1147 if (image->data != NULL)
1148 {
1149 int bytesPerPixel;
1150 nifti_datatype_sizes(datatype, &bytesPerPixel, NULL);
1151 void *data = calloc(image->nvox, bytesPerPixel);
1152
1153 switch (datatype)
1154 {
1155 case DT_UINT8:
1156 internal::convertData<uint8_t>(image->data, image->datatype, image->nvox, static_cast<uint8_t *>(data));
1157 break;
1158
1159 case DT_INT16:
1160 internal::convertData<int16_t>(image->data, image->datatype, image->nvox, static_cast<int16_t *>(data));
1161 break;
1162
1163 case DT_INT32:
1164 internal::convertData<int32_t>(image->data, image->datatype, image->nvox, static_cast<int32_t *>(data));
1165 break;
1166
1167 case DT_FLOAT32:
1168 internal::convertData<float>(image->data, image->datatype, image->nvox, static_cast<float *>(data));
1169 break;
1170
1171 case DT_FLOAT64:
1172 internal::convertData<double>(image->data, image->datatype, image->nvox, static_cast<double *>(data));
1173 break;
1174
1175 case DT_INT8:
1176 internal::convertData<int8_t>(image->data, image->datatype, image->nvox, static_cast<int8_t *>(data));
1177 break;
1178
1179 case DT_UINT16:
1180 internal::convertData<uint16_t>(image->data, image->datatype, image->nvox, static_cast<uint16_t *>(data));
1181 break;
1182
1183 case DT_UINT32:
1184 internal::convertData<uint32_t>(image->data, image->datatype, image->nvox, static_cast<uint32_t *>(data));
1185 break;
1186
1187 case DT_INT64:
1188 internal::convertData<int64_t>(image->data, image->datatype, image->nvox, static_cast<int64_t *>(data));
1189 break;
1190
1191 case DT_UINT64:
1192 internal::convertData<uint64_t>(image->data, image->datatype, image->nvox, static_cast<uint64_t *>(data));
1193 break;
1194
1195 default:
1196 throw std::runtime_error("Unsupported data type (" + std::string(nifti_datatype_string(datatype)) + ")");
1197 }
1198
1199 nifti_image_unload(image);
1200 image->data = data;
1201 }
1202
1203 image->datatype = datatype;
1204 nifti_datatype_sizes(datatype, &image->nbyper, &image->swapsize);
1205 }
1206
1207 inline void NiftiImage::changeDatatype (const std::string &datatype)
1208 {
1209 changeDatatype(internal::stringToDatatype(datatype));
1210 }
1211
1212 template <typename SourceType>
1213 inline void NiftiImage::replaceData (const std::vector<SourceType> &data, const short datatype)
1214 {
1215 if (this->isNull())
1216 return;
1217 else if (data.size() != image->nvox)
1218 throw std::runtime_error("New data length does not match the number of voxels in the image");
1219
1220 if (datatype != DT_NONE)
1221 {
1222 nifti_image_unload(image);
1223 image->datatype = datatype;
1224 nifti_datatype_sizes(datatype, &image->nbyper, &image->swapsize);
1225 }
1226
1227 if (image->data == NULL)
1228 image->data = calloc(image->nvox, image->nbyper);
1229 internal::replaceData<SourceType>(data.begin(), data.end(), image->data, image->datatype);
1230
1231 image->scl_slope = 0.0;
1232 image->scl_inter = 0.0;
1233 image->cal_min = static_cast<float>(*std::min_element(data.begin(), data.end()));
1234 image->cal_max = static_cast<float>(*std::max_element(data.begin(), data.end()));
1235 }
1236
1237 inline void NiftiImage::toFile (const std::string fileName, const short datatype) const
1238 {
1239 // Copy the source image only if the datatype will be changed
1240 NiftiImage imageToWrite(image, datatype != DT_NONE);
1241
1242 if (datatype == DT_NONE)
1243 imageToWrite.setPersistence(true);
1244 else
1245 imageToWrite.changeDatatype(datatype);
1246
1247 const int status = nifti_set_filenames(imageToWrite, fileName.c_str(), false, true);
1248 if (status != 0)
1249 throw std::runtime_error("Failed to set filenames for NIfTI object");
1250 nifti_image_write(imageToWrite);
1251 }
1252
1253 inline void NiftiImage::toFile (const std::string fileName, const std::string &datatype) const
1254 {
1255 toFile(fileName, internal::stringToDatatype(datatype));
1256 }
1257
1258 inline Rcpp::RObject NiftiImage::toArray () const
1259 {
1260 Rcpp::RObject array;
1261
1262 if (this->isNull())
1263 return array;
1264
1265 switch (image->datatype)
1266 {
1267 case DT_UINT8:
1268 case DT_INT16:
1269 case DT_INT32:
1270 case DT_INT8:
1271 case DT_UINT16:
1272 case DT_UINT32:
1273 case DT_INT64:
1274 case DT_UINT64:
1275 array = internal::imageDataToArray<INTSXP>(image);
1276 break;
1277
1278 case DT_FLOAT32:
1279 case DT_FLOAT64:
1280 array = internal::imageDataToArray<REALSXP>(image);
1281 break;
1282
1283 default:
1284 throw std::runtime_error("Unsupported data type (" + std::string(nifti_datatype_string(image->datatype)) + ")");
1285 }
1286
1287 internal::addAttributes(array, image);
1288 array.attr("class") = Rcpp::CharacterVector::create("niftiImage", "array");
1289 return array;
1290 }
1291
1292 inline Rcpp::RObject NiftiImage::toPointer (const std::string label) const
1293 {
1294 if (this->isNull())
1295 return Rcpp::RObject();
1296 else
1297 {
1298 Rcpp::RObject string = Rcpp::wrap(label);
1299 internal::addAttributes(string, image, false);
1300 string.attr("class") = Rcpp::CharacterVector::create("internalImage", "niftiImage");
1301 return string;
1302 }
1303 }
1304
1305 inline Rcpp::RObject NiftiImage::toArrayOrPointer (const bool internal, const std::string label) const
1306 {
1307 return (internal ? toPointer(label) : toArray());
1308 }
1309
1310 inline Rcpp::RObject NiftiImage::headerToList () const
1311 {
1312 if (this->image == NULL)
1313 return Rcpp::RObject();
1314
1315 nifti_1_header header = nifti_convert_nim2nhdr(this->image);
1316 Rcpp::List result;
1317
1318 result["sizeof_hdr"] = header.sizeof_hdr;
1319
1320 result["dim_info"] = int(header.dim_info);
1321 result["dim"] = std::vector<short>(header.dim, header.dim+8);
1322
1323 result["intent_p1"] = header.intent_p1;
1324 result["intent_p2"] = header.intent_p2;
1325 result["intent_p3"] = header.intent_p3;
1326 result["intent_code"] = header.intent_code;
1327
1328 result["datatype"] = header.datatype;
1329 result["bitpix"] = header.bitpix;
1330
1331 result["slice_start"] = header.slice_start;
1332 result["pixdim"] = std::vector<float>(header.pixdim, header.pixdim+8);
1333 result["vox_offset"] = header.vox_offset;
1334 result["scl_slope"] = header.scl_slope;
1335 result["scl_inter"] = header.scl_inter;
1336 result["slice_end"] = header.slice_end;
1337 result["slice_code"] = int(header.slice_code);
1338 result["xyzt_units"] = int(header.xyzt_units);
1339 result["cal_max"] = header.cal_max;
1340 result["cal_min"] = header.cal_min;
1341 result["slice_duration"] = header.slice_duration;
1342 result["toffset"] = header.toffset;
1343 result["descrip"] = std::string(header.descrip, 80);
1344 result["aux_file"] = std::string(header.aux_file, 24);
1345
1346 result["qform_code"] = header.qform_code;
1347 result["sform_code"] = header.sform_code;
1348 result["quatern_b"] = header.quatern_b;
1349 result["quatern_c"] = header.quatern_c;
1350 result["quatern_d"] = header.quatern_d;
1351 result["qoffset_x"] = header.qoffset_x;
1352 result["qoffset_y"] = header.qoffset_y;
1353 result["qoffset_z"] = header.qoffset_z;
1354 result["srow_x"] = std::vector<float>(header.srow_x, header.srow_x+4);
1355 result["srow_y"] = std::vector<float>(header.srow_y, header.srow_y+4);
1356 result["srow_z"] = std::vector<float>(header.srow_z, header.srow_z+4);
1357
1358 result["intent_name"] = std::string(header.intent_name, 16);
1359 result["magic"] = std::string(header.magic, 4);
1360
1361 result.attr("class") = Rcpp::CharacterVector::create("niftiHeader");
1362
1363 return result;
1364 }
1365
1366 } // namespace
1367
1368 #endif
RNifti::NiftiImage::NiftiImage
NiftiImage()
Default constructor.
Definition: NiftiImage.h:162
RNifti::NiftiImage::dropData
void dropData()
Drop the data from the image, retaining only the metadata.
Definition: NiftiImage.h:399
RNifti::NiftiImage::Block::index
const int index
The location along dimension.
Definition: NiftiImage.h:33
RNifti::NiftiImage::dim
std::vector< int > dim() const
Return the dimensions of the image.
Definition: NiftiImage.h:332
RNifti::NiftiImage::operator->
const nifti_image * operator->() const
Allows a NiftiImage object to be treated as a pointer to a const nifti_image.
Definition: NiftiImage.h:260
RNifti::NiftiImage::headerToList
Rcpp::RObject headerToList() const
Create an R list containing raw image metadata.
Definition: NiftiImage.h:1310
RNifti::NiftiImage::persistent
bool persistent
Marker of persistence, which determines whether the nifti_image should be freed on destruction...
Definition: NiftiImage.h:99
RNifti::NiftiImage::pixdim
std::vector< float > pixdim() const
Return the dimensions of the pixels or voxels in the image.
Definition: NiftiImage.h:344
RNifti::NiftiImage::NiftiImage
NiftiImage(const std::string &path, const bool readData=true)
Initialise using a path string.
Definition: NiftiImage.h:215
RNifti::NiftiImage::NiftiImage
NiftiImage(nifti_image *const image, const bool copy=false)
Initialise using an existing nifti_image pointer.
Definition: NiftiImage.h:197
RNifti::NiftiImage::rescale
void rescale(const std::vector< float > &scales)
Rescale the image, changing its image dimensions and pixel dimensions.
Definition: NiftiImage.h:1016
RNifti
Definition: NiftiImage.h:16
RNifti::NiftiImage::initFromMriImage
void initFromMriImage(const Rcpp::RObject &object, const bool copyData=true)
Initialise the object from a reference object of class "MriImage".
Definition: NiftiImage.h:670
RNifti::NiftiImage::Block::operator=
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:57
RNifti::NiftiImage::isNull
bool isNull() const
Determine whether or not the internal pointer is NULL.
Definition: NiftiImage.h:310
RNifti::NiftiImage::Block
Inner class referring to a subset of an image.
Definition: NiftiImage.h:29
RNifti::NiftiImage::slice
Block slice(const int i)
Extract a slice block from a 3D image.
Definition: NiftiImage.h:462
RNifti::NiftiImage::update
void update(const SEXP array)
Update the image from an R array.
Definition: NiftiImage.h:1036
RNifti::NiftiImage::toPointer
Rcpp::RObject toPointer(const std::string label) const
Create an internal image to pass back to R.
Definition: NiftiImage.h:1292
RNifti::NiftiImage::getData
std::vector< TargetType > getData() const
Extract a vector of data from the image, casting it to any required element type. ...
Definition: NiftiImage.h:1132
RNifti::NiftiImage::initFromArray
void initFromArray(const Rcpp::RObject &object, const bool copyData=true)
Initialise the object from an R array.
Definition: NiftiImage.h:842
RNifti::NiftiImage::initFromNiftiS4
void initFromNiftiS4(const Rcpp::RObject &object, const bool copyData=true)
Initialise the object from an S4 object of class "nifti".
Definition: NiftiImage.h:574
RNifti::NiftiImage::slice
const Block slice(const int i) const
Extract a slice block from a 3D image.
Definition: NiftiImage.h:455
RNifti::NiftiImage::volume
const Block volume(const int i) const
Extract a volume block from a 4D image.
Definition: NiftiImage.h:469
RNifti::NiftiImage::isPersistent
bool isPersistent() const
Determine whether or not the image is marked as persistent.
Definition: NiftiImage.h:315
RNifti::NiftiImage::replaceData
void 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:1213
RNifti::NiftiImage::setPixunits
void setPixunits(const std::vector< std::string > &pixunits)
Modify the pixel dimension units.
Definition: NiftiImage.h:991
RNifti::NiftiImage::NiftiImage
NiftiImage(const Block &source)
Initialise from a block, copying in the data.
Definition: NiftiImage.h:182
RNifti::NiftiImage::~NiftiImage
virtual ~NiftiImage()
Destructor which frees the wrapped pointer, unless the object is marked as persistent.
Definition: NiftiImage.h:236
RNifti::NiftiImage
Thin wrapper around a C-style nifti_image struct that allows C++-style destruction.
Definition: NiftiImage.h:22
RNifti::NiftiImage::setPersistence
void setPersistence(const bool persistent)
Marked the image as persistent, so that it can be passed back to R.
Definition: NiftiImage.h:298
RNifti::NiftiImage::NiftiImage
NiftiImage(const NiftiImage &source)
Copy constructor.
Definition: NiftiImage.h:169
RNifti::NiftiImage::nBlocks
int nBlocks() const
Return the number of blocks in the image.
Definition: NiftiImage.h:424
RNifti::NiftiImage::changeDatatype
void changeDatatype(const short datatype)
Change the datatype of the image, casting the pixel data if present.
Definition: NiftiImage.h:1142
RNifti::NiftiImage::Block::dimension
const int dimension
The dimension along which the block applies (which should be the last)
Definition: NiftiImage.h:32
RNifti::NiftiImage::Block::image
const NiftiImage & image
The parent image.
Definition: NiftiImage.h:31
RNifti::NiftiImage::initFromList
void initFromList(const Rcpp::RObject &object)
Initialise the object from an R list with named elements, which can only contain metadata.
Definition: NiftiImage.h:753
RNifti::NiftiImage::toFile
void toFile(const std::string fileName, const short datatype=DT_NONE) const
Write the image to a NIfTI-1 file.
Definition: NiftiImage.h:1237
RNifti::NiftiImage::updatePixdim
void updatePixdim(const std::vector< float > &pixdim)
Modify the pixel dimensions, and potentially the xform matrices to match.
Definition: NiftiImage.h:938
RNifti::NiftiImage::block
const Block block(const int i) const
Extract a block from the image.
Definition: NiftiImage.h:439
RNifti::NiftiImage::copy
void copy(const nifti_image *source)
Copy the contents of a nifti_image to create a new image.
Definition: NiftiImage.h:523
RNifti::NiftiImage::xform
mat44 xform(const bool preferQuaternion=true) const
Obtain an xform matrix, indicating the orientation of the image.
Definition: NiftiImage.h:1080
RNifti::NiftiImage::image
nifti_image * image
The wrapped nifti_image pointer.
Definition: NiftiImage.h:98
RNifti::NiftiImage::Block::getData
std::vector< TargetType > getData() const
Extract a vector of data from a block, casting it to any required element type.
Definition: NiftiImage.h:1117
RNifti::NiftiImage::toArrayOrPointer
Rcpp::RObject toArrayOrPointer(const bool internal, const std::string label) const
A conditional method that calls either toArray or toPointer.
Definition: NiftiImage.h:1305
RNifti::NiftiImage::sexpTypeToNiftiType
static short sexpTypeToNiftiType(const int sexpType)
Convert between R SEXP object type and nifti_image datatype codes.
Definition: NiftiImage.h:87
RNifti::NiftiImage::nDims
int nDims() const
Return the number of dimensions in the image.
Definition: NiftiImage.h:320
RNifti::NiftiImage::volume
Block volume(const int i)
Extract a volume block from a 4D image.
Definition: NiftiImage.h:476
RNifti::NiftiImage::drop
NiftiImage & drop()
Drop unitary dimensions.
Definition: NiftiImage.h:358
RNifti::NiftiImage::toArray
Rcpp::RObject toArray() const
Create an R array from the image.
Definition: NiftiImage.h:1258
RNifti::NiftiImage::Block::Block
Block(const NiftiImage &image, const int dimension, const int index)
Standard constructor for this class.
Definition: NiftiImage.h:42
RNifti::NiftiImage::block
Block block(const int i)
Extract a block from the image.
Definition: NiftiImage.h:448