NiftiImageData.h

Go to the documentation of this file.

/*
SyneRBI Synergistic Image Reconstruction Framework (SIRF)
Copyright 2017 - 2020 University College London
 
This is software developed for the Collaborative Computational
Project in Synergistic Reconstruction for Biomedical Imaging (formerly CCP PETMR)
(http://www.ccpsynerbi.ac.uk/).
 
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
 
*/
 
#pragma once
 
#include <nifti1_io.h>
#include <vector>
#include <string>
#include <memory>
#include <iostream>
#include <sstream>
#include "sirf/common/ANumRef.h"
#include "sirf/common/ImageData.h"
#include <_reg_tools.h>
 
namespace sirf {
 
template<class dataType>
class NiftiImageData : public ImageData
{
public:
 
    typedef ImageData::Iterator BaseIter;
    typedef ImageData::Iterator_const BaseIter_const;
    class Iterator : public BaseIter {
    public:
        Iterator(dataType *iter) : _iter(iter)
        {}
        Iterator& operator=(const Iterator& iter)
        {
            _iter = iter._iter;
            _ref.copy(iter._ref);
            return *this;
        }
        virtual Iterator& operator++()
        {
            ++_iter;
            return *this;
        }
        virtual bool operator==(const BaseIter& an_iter) const
        {
            const Iterator& iter = (const Iterator&)an_iter;
            return _iter == iter._iter;
        }
        virtual bool operator!=(const BaseIter& an_iter) const
        {
            const Iterator& iter = (const Iterator&)an_iter;
            return _iter != iter._iter;
        }
        virtual FloatRef& operator*()
        {
            dataType& v = *_iter;
            _ref.set_ptr(&v);
            return _ref;
        }
    private:
        dataType *_iter;
        FloatRef _ref;
    };
    class Iterator_const : public BaseIter_const {
    public:
        Iterator_const(const dataType *iter) : _iter(iter)
        {}
        Iterator_const& operator=(const Iterator_const& iter)
        {
            _iter = iter._iter;
            _ref.copy(iter._ref);
            return *this;
        }
        virtual Iterator_const& operator++()
        {
            ++_iter;
            return *this;
        }
        virtual bool operator==(const BaseIter_const& an_iter) const
        {
            const Iterator_const& iter = (const Iterator_const&)an_iter;
            return _iter == iter._iter;
        }
        virtual bool operator!=(const BaseIter_const& an_iter) const
        {
            const Iterator_const& iter = (const Iterator_const&)an_iter;
            return _iter != iter._iter;
        }
        virtual const FloatRef& operator*() const
        {
            const dataType& v = *_iter;
            _ref.set_ptr((void*)&v);
            return _ref;
        }
    private:
        const dataType *_iter;
        mutable FloatRef _ref;
    };
 
    NiftiImageData() {}
 
    virtual ~NiftiImageData() {}
 
    NiftiImageData(const NiftiImageData& to_copy);
 
    NiftiImageData& operator=(const NiftiImageData& to_copy);
 
    NiftiImageData(const ImageData& to_copy);
 
    NiftiImageData& operator=(const ImageData& to_copy);
 
    NiftiImageData(const std::string &filename);
 
    NiftiImageData(const nifti_image &image_nifti);
 
    template<class inputType>
    NiftiImageData(const inputType * const data, const VoxelisedGeometricalInfo3D &geom, const bool is_tensor = false)
    {
        this->_nifti_image = create_from_geom_info(geom, is_tensor);
 
        // Set the datatype
        if      (typeid(inputType) == typeid(bool))               this->set_up_data(DT_BINARY);
        else if (typeid(inputType) == typeid(signed char))        this->set_up_data(DT_INT8);
        else if (typeid(inputType) == typeid(signed short))       this->set_up_data(DT_INT16);
        else if (typeid(inputType) == typeid(signed int))         this->set_up_data(DT_INT32);
        else if (typeid(inputType) == typeid(float))              this->set_up_data(DT_FLOAT32);
        else if (typeid(inputType) == typeid(double))             this->set_up_data(DT_FLOAT64);
        else if (typeid(inputType) == typeid(unsigned char))      this->set_up_data(DT_UINT8);
        else if (typeid(inputType) == typeid(unsigned short))     this->set_up_data(DT_UINT16);
        else if (typeid(inputType) == typeid(unsigned int))       this->set_up_data(DT_UINT32);
        else if (typeid(inputType) == typeid(signed long long))   this->set_up_data(DT_INT64);
        else if (typeid(inputType) == typeid(unsigned long long)) this->set_up_data(DT_UINT64);
        else if (typeid(inputType) == typeid(long double))        this->set_up_data(DT_FLOAT128);
        else {
            std::stringstream ss;
            ss << "NiftiImageData constructor from raw array: ";
            ss << typeid (inputType).name();
            ss << " (bytes per voxel: ";
            ss << sizeof(inputType) << ").";
            throw std::runtime_error(ss.str());
        }
 
        for (unsigned i=0; i<_nifti_image->nvox; ++i)
            this->_data[i] = dataType(data[i]);
    }
 
    static std::shared_ptr<nifti_image> create_from_geom_info(const VoxelisedGeometricalInfo3D &geom, const bool is_tensor=false, const NREG_TRANS_TYPE tensor_type=NREG_TRANS_TYPE::DEF_FIELD);
 
    static void construct_NiftiImageData_from_complex_im_real_component(std::shared_ptr<NiftiImageData> &out_sptr, const std::shared_ptr<const ImageData> in_sptr);
 
    static void construct_NiftiImageData_from_complex_im_imag_component(std::shared_ptr<NiftiImageData> &out_sptr, const std::shared_ptr<const ImageData> in_sptr);
 
    static void construct_NiftiImageData_from_complex_im(std::shared_ptr<NiftiImageData> &out_real_sptr, std::shared_ptr<NiftiImageData> &out_imag_sptr, const std::shared_ptr<const ImageData> in_sptr);
 
    bool operator==(const NiftiImageData &other) const;
 
    bool operator!=(const NiftiImageData &other) const;
 
    NiftiImageData& operator+=(const NiftiImageData &rhs);
 
    friend NiftiImageData operator+(NiftiImageData lhs, const NiftiImageData& rhs)
    {
        lhs += rhs;
        return lhs;
    }
 
    NiftiImageData& operator-=(const NiftiImageData &rhs);
 
    friend NiftiImageData operator-(NiftiImageData lhs, const NiftiImageData& rhs)
    {
        lhs -= rhs;
        return lhs;
    }
 
    NiftiImageData& operator+=(const float);
 
    friend NiftiImageData operator+(NiftiImageData lhs, const float val)
    {
        lhs += val;
        return lhs;
    }
 
    NiftiImageData& operator-=(const float);
 
    friend NiftiImageData operator-(NiftiImageData lhs, const float val)
    {
        lhs -= val;
        return lhs;
    }
 
    NiftiImageData& operator*=(const float);
    NiftiImageData& operator*=(const NiftiImageData &rhs);
 
    friend NiftiImageData operator*(NiftiImageData lhs, const float val)
    {
        lhs *= val;
        return lhs;
    }
 
    friend NiftiImageData operator*(NiftiImageData lhs, const NiftiImageData& rhs)
    {
        lhs *= rhs;
        return lhs;
    }
 
    NiftiImageData& operator/=(const float);
    // /// Division operator
    NiftiImageData& operator/=(const NiftiImageData &rhs);
    
    // /// Division operator
    friend NiftiImageData operator/(NiftiImageData lhs, const float val)
    {
        lhs /= val;
        return lhs;
    }
 
    friend NiftiImageData operator/(NiftiImageData lhs, const NiftiImageData& rhs)
    {
        lhs /= rhs;
        return lhs;
    }
 
    float operator()(const int index) const;
 
    float &operator()(const int index);
 
    float operator()(const int index[7]) const;
 
    float &operator()(const int index[7]);
 
    float operator()(const int x, const int y, const int z, const int t=0, const int u=0, const int v=0, const int w=0) const;
 
    float &operator()(const int x, const int y, const int z, const int t=0, const int u=0, const int v=0, const int w=0);
 
    bool is_initialised() const { return (_nifti_image && _data && _nifti_image->datatype == DT_FLOAT32 ? true : false); }
 
    std::shared_ptr<const nifti_image> get_raw_nifti_sptr() const;
 
    std::shared_ptr<nifti_image> get_raw_nifti_sptr();
 
    virtual void write(const std::string &filename, const int datatype) const;
 
    virtual void write(const std::string &filename) const { this->write(filename,-1); }
 
    float get_max() const;
 
    float get_min() const;
 
    float get_mean() const;
 
    float get_variance() const;
 
    float get_standard_deviation() const;
 
    float get_sum() const;
 
    unsigned get_nan_count() const;
 
    void fill(const float v);
 
    void fill(const dataType *v);
 
    void fill(const NiftiImageData &im);
 
    float get_norm(const NiftiImageData&) const;
 
    const int* get_dimensions() const;
 
    size_t get_num_voxels() const;
 
    void print_header() const;
 
    static void print_headers(const std::vector<const NiftiImageData*> &ims);
 
    void crop(const int min_index[7], const int max_index[7]);
 
    void pad(const int min_index[7], const int max_index[7], const dataType val = 0);
 
    int get_1D_index(const int idx[7]) const;
 
    int get_original_datatype() const { return _original_datatype; }
 
    static bool are_equal_to_given_accuracy(const NiftiImageData &im1, const NiftiImageData &im2, const float required_accuracy_compared_to_max);
 
    bool is_in_bounds(const int index[7]) const;
 
    bool is_in_bounds(const int index) const;
 
    bool is_same_size(const NiftiImageData &im) const;
 
    static bool do_nifti_image_metadata_match(const NiftiImageData &im1, const NiftiImageData &im2, bool verbose);
 
    static void dump_headers(const std::vector<const NiftiImageData*> &ims);
 
    template<typename T>
    static void dump_nifti_element(const std::vector<const NiftiImageData*> &ims, const std::string &name, const T &call_back);
 
    template<typename T>
    static void dump_nifti_element(const std::vector<const NiftiImageData*> &ims, const std::string &name, const T &call_back, const unsigned num_elems);
 
    static std::string get_headers(const std::vector<const NiftiImageData<dataType>*> &ims);
    template<typename T>
    static std::string get_nifti_element(const std::vector<const NiftiImageData*> &ims, const std::string &name, const T &call_back);
    template<typename T>
    static std::string get_nifti_element(const std::vector<const NiftiImageData*> &ims, const std::string &name, const T &call_back, const unsigned num_elems);
 
    void set_voxel_spacing(const float factors[3], const int interpolation_order);
 
    void kernel_convolution(const float sigma, NREG_CONV_KERNEL_TYPE conv_type = GAUSSIAN_KERNEL);
 
    bool get_contains_nans() const { return (this->get_nan_count() > 0); }
 
    void flip_along_axis(const unsigned axis);
 
    void mirror_along_axis(const unsigned axis);
 
    dataType get_inner_product(const NiftiImageData &other) const;
 
    void normalise_zero_and_one();
 
    void standardise();
 
protected:
 
    enum NiftiImageDataType { _general, _3D, _3DTensor, _3DDisp, _3DDef};
 
    enum MathsType { ADD, sub, mul, div};
 
    std::shared_ptr<nifti_image>  _nifti_image;
 
    float *_data = NULL;
 
    int _original_datatype = -1;
 
    void check_dimensions(const enum NiftiImageDataType image_type = _general);
 
    void set_up_data(const int original_datatype);
 
    void maths(const NiftiImageData& c, const MathsType type);
 
    void maths(const float val, const MathsType type);
 
    static void open_nifti_image(std::shared_ptr<nifti_image> &image, const std::string &filename);
 
    static void copy_nifti_image(std::shared_ptr<nifti_image> &output_image_sptr, const std::shared_ptr<nifti_image> &image_to_copy_sptr);
 
private:
 
    void change_datatype(const int datatype);
 
    template<typename newType>
    static void change_datatype(NiftiImageData<dataType> &im)
    {
        if (im.get_raw_nifti_sptr()->datatype == DT_BINARY)   return change_datatype<newType,bool>              (im);
        if (im.get_raw_nifti_sptr()->datatype == DT_INT8)     return change_datatype<newType,signed char>       (im);
        if (im.get_raw_nifti_sptr()->datatype == DT_INT16)    return change_datatype<newType,signed short>      (im);
        if (im.get_raw_nifti_sptr()->datatype == DT_INT32)    return change_datatype<newType,signed int>        (im);
        if (im.get_raw_nifti_sptr()->datatype == DT_FLOAT32)  return change_datatype<newType,float>             (im);
        if (im.get_raw_nifti_sptr()->datatype == DT_FLOAT64)  return change_datatype<newType,double>            (im);
        if (im.get_raw_nifti_sptr()->datatype == DT_UINT8)    return change_datatype<newType,unsigned char>     (im);
        if (im.get_raw_nifti_sptr()->datatype == DT_UINT16)   return change_datatype<newType,unsigned short>    (im);
        if (im.get_raw_nifti_sptr()->datatype == DT_UINT32)   return change_datatype<newType,unsigned int>      (im);
        if (im.get_raw_nifti_sptr()->datatype == DT_INT64)    return change_datatype<newType,signed long long>  (im);
        if (im.get_raw_nifti_sptr()->datatype == DT_UINT64)   return change_datatype<newType,unsigned long long>(im);
        if (im.get_raw_nifti_sptr()->datatype == DT_FLOAT128) return change_datatype<newType,long double>       (im);
 
        std::stringstream ss;
        ss << "change_datatype not implemented for your data type: ";
        ss << nifti_datatype_string(im.get_raw_nifti_sptr()->datatype);
        ss << " (bytes per voxel: ";
        ss << im.get_raw_nifti_sptr()->nbyper << ").";
        throw std::runtime_error(ss.str());
    }
 
    template<typename newType, typename oldType>
    static void change_datatype(NiftiImageData &image)
    {
        // If the two types are equal, nothing to be done.
        if (typeid (newType) == typeid(oldType))
            return;
 
        nifti_image *im = image.get_raw_nifti_sptr().get();
 
        // Copy the original array
        oldType *originalArray = static_cast<oldType*>(malloc(im->nvox*im->nbyper));
        memcpy(originalArray, im->data, im->nvox*im->nbyper);
        // Reset image array
        free(im->data);
 
        // Set the datatype
        if      (typeid(newType) == typeid(bool))               im->datatype = DT_BINARY;
        else if (typeid(newType) == typeid(signed char))        im->datatype = DT_INT8;
        else if (typeid(newType) == typeid(signed short))       im->datatype = DT_INT16;
        else if (typeid(newType) == typeid(signed int))         im->datatype = DT_INT32;
        else if (typeid(newType) == typeid(float))              im->datatype = DT_FLOAT32;
        else if (typeid(newType) == typeid(double))             im->datatype = DT_FLOAT64;
        else if (typeid(newType) == typeid(unsigned char))      im->datatype = DT_UINT8;
        else if (typeid(newType) == typeid(unsigned short))     im->datatype = DT_UINT16;
        else if (typeid(newType) == typeid(unsigned int))       im->datatype = DT_UINT32;
        else if (typeid(newType) == typeid(signed long long))   im->datatype = DT_INT64;
        else if (typeid(newType) == typeid(unsigned long long)) im->datatype = DT_UINT64;
        else if (typeid(newType) == typeid(long double))        im->datatype = DT_FLOAT128;
        else {
            std::stringstream ss;
            ss << "change_datatype not implemented for your new data type: ";
            ss << typeid (newType).name();
            ss << " (bytes per voxel: ";
            ss << sizeof(newType) << ").";
            throw std::runtime_error(ss.str());
        }
 
        // Set the nbyper and swap size from the datatype
        nifti_datatype_sizes(im->datatype, &im->nbyper, &im->swapsize);
 
        // Copy data
        im->data = static_cast<void*>(calloc(im->nvox,sizeof(newType)));
        newType *dataPtr = static_cast<newType*>(im->data);
        for (size_t i = 0; i < im->nvox; i++)
           dataPtr[i] = newType(originalArray[i]);
 
        free(originalArray);
        return;
    }
 
    // ------------------------------------------------------------------------------ //
    // Pure virtual methods from ImageData
    // ------------------------------------------------------------------------------ //
public:
    std::unique_ptr<NiftiImageData> clone() const
    {
    return std::unique_ptr<NiftiImageData>(this->clone_impl());
    }
    virtual Iterator& begin()
    {
        _begin.reset(new Iterator(_data));
        return *_begin;
    }
    virtual Iterator_const& begin() const
    {
        _begin_const.reset(new Iterator_const(_data));
        return *_begin_const;
    }
    virtual Iterator& end()
    {
        _end.reset(new Iterator(_data+_nifti_image->nvox));
        return *_end;
    }
    virtual Iterator_const& end() const
    {
        _end_const.reset(new Iterator_const(_data+_nifti_image->nvox));
        return *_end_const;
    }
    /*
    unsigned int items() const { return 1; }
    virtual void dot      (const DataContainer& a_x, void* ptr) const;
    virtual void axpby    (const void* ptr_a, const DataContainer& a_x, const void* ptr_b, const DataContainer& a_y);
    virtual void xapyb    (const DataContainer& a_x, const void* ptr_a, const DataContainer& a_y, const void* ptr_b);
    virtual void xapyb    (const DataContainer& a_x, const DataContainer& a_a, const DataContainer& a_y, const DataContainer& a_b);
    virtual float norm() const;
    virtual void scale(float s);
    virtual void multiply (const DataContainer& a_x, const DataContainer& a_y);
    virtual void divide   (const DataContainer& a_x, const DataContainer& a_y);
    virtual void maximum(const DataContainer& x, const DataContainer& y);
    virtual void minimum(const DataContainer& x, const DataContainer& y);
*/
protected:
    virtual NiftiImageData* clone_impl() const
    {
    return new NiftiImageData(*this);
    }
    virtual ObjectHandle<DataContainer>* new_data_container_handle() const
    {
        return new ObjectHandle<DataContainer>
            (std::shared_ptr<DataContainer>(new NiftiImageData));
    }
 
public:
    unsigned int items() const { return 1; }
    virtual void sum      (void* ptr) const;
    virtual void max      (void* ptr) const;
    virtual void dot      (const DataContainer& a_x, void* ptr) const;
    virtual void axpby    (const void* ptr_a, const DataContainer& a_x, const void* ptr_b, const DataContainer& a_y);
    virtual void xapyb    (const DataContainer& a_x, const void* ptr_a, const DataContainer& a_y, const void* ptr_b);
    virtual void xapyb    (const DataContainer& a_x, const DataContainer& a_a, const DataContainer& a_y, const DataContainer& a_b);
    virtual void xapyb(
        const DataContainer& a_x, const void* ptr_a,
        const DataContainer& a_y, const DataContainer& a_b);
    virtual float norm() const;
    virtual void multiply (const DataContainer& a_x, const DataContainer& a_y);
    virtual void divide   (const DataContainer& a_x, const DataContainer& a_y);
    virtual void maximum(const DataContainer& x, const DataContainer& y);
    virtual void minimum(const DataContainer& x, const DataContainer& y);
    virtual void power(const DataContainer& x, const DataContainer& y);
    virtual void multiply(const DataContainer& a_x, const void* a_y);
    virtual void add(const DataContainer& a_x, const void* a_y);
    virtual void maximum(const DataContainer& x, const void* a_y);
    virtual void minimum(const DataContainer& x, const void* a_y);
    virtual void power(const DataContainer& x, const void* a_y);
    virtual void exp(const DataContainer& x);
    virtual void log(const DataContainer& x);
    virtual void sqrt(const DataContainer& x);
    virtual void sign(const DataContainer& x);
    virtual void abs(const DataContainer& x);
 
    virtual Dimensions dimensions() const
    {
        Dimensions dim;
        int *d = _nifti_image->dim;
        dim["x"] = d[1];
        dim["y"] = d[2];
        dim["z"] = d[3];
        dim["t"] = d[4];
        dim["u"] = d[5];
        dim["v"] = d[6];
        dim["w"] = d[7];
        return dim;
    }
    void unary_op(const DataContainer& a_x, dataType(*f)(dataType));
    void semibinary_op(const DataContainer& a_x, const void* a_y, dataType(*f)(dataType, dataType));
    void binary_op(const DataContainer& a_x, const DataContainer& a_y, dataType(*f)(dataType, dataType));
    virtual void set_up_geom_info();
protected:
    mutable std::shared_ptr<Iterator> _begin;
    mutable std::shared_ptr<Iterator> _end;
    mutable std::shared_ptr<Iterator_const> _begin_const;
    mutable std::shared_ptr<Iterator_const> _end_const;
};
}