multidim_array.h

Go to the documentation of this file.

/***************************************************************************
 *
 * Authors:     Carlos Oscar S. Sorzano (coss@cnb.csic.es)
 *
 * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
 * 02111-1307  USA
 *
 *  All comments concerning this program package may be sent to the
 *  e-mail address 'xmipp@cnb.csic.es'
 ***************************************************************************/

#ifndef CORE_MULTIDIM_ARRAY_H
#define CORE_MULTIDIM_ARRAY_H

#include "xmipp_macros.h" // to define XMIPP_MMAP
#ifdef XMIPP_MMAP
#include <sys/mman.h>
#endif
#include <complex>
#include "matrix1d.h"
#include "xmipp_random_mode.h"
#include "multidim_array_base.h"
#include "xmipp_memory.h"
#include "axis_view.h"

template<typename T>
class Matrix2D;

extern int bestPrecision(float F, int _width);
extern String floatToString(float F, int _width, int _prec);




// Forward declarations ====================================================
template<typename T>
class MultidimArray;

template<typename T>
void coreArrayByScalar(const MultidimArray<T>& op1, const T& op2,
                       MultidimArray<T>& result, char operation);

template<typename T>
void coreScalarByArray(const T& op1, const MultidimArray<T>& op2,
                       MultidimArray<T>& result, char operation);

template<typename T>
void coreArrayByArray(const MultidimArray<T>& op1, const MultidimArray<T>& op2,
                      MultidimArray<T>& result, char operation);
template<typename T>
void selfCoreArrayByArrayMask(const MultidimArray<T>& op1, const MultidimArray<T>& op2,
                          MultidimArray<T>& result, char operation,
                          const MultidimArray<T> *mask);





template<typename T>
class MultidimArray: public MultidimArrayBase
{
public:
    /* The array itself.
       The array is always a 3D array (Z,Y,X). For vectors the size of the array
       is (1,1,X) and for matrices (1,Y,X). The pixel (i,j) (y,x) is at the
       position data[i*Xdim+j] or data[y*Xdim+x]
    */
    T* data;

public:


    MultidimArray()
    {
        coreInit();
    }

    MultidimArray(size_t Ndim, size_t Zdim, size_t Ydim, size_t Xdim, T *data) {
        this->coreInit();
        this->setDimensions(Xdim, Ydim, Zdim, Ndim);
        this->data = data;
        this->nzyxdimAlloc = this->nzyxdim;
        this->destroyData = false;
    }


    MultidimArray(size_t Ndim, int Zdim, int Ydim, int Xdim)
    {
        coreInit();
        coreAllocate(Ndim, Zdim, Ydim, Xdim);
    }

    MultidimArray( int Zdim, int Ydim, int Xdim)
    {
        coreInit();
        coreAllocate(1UL, Zdim, Ydim, Xdim);
    }

    MultidimArray(int Ydim, int Xdim)
    {
        coreInit();
        coreAllocate(1UL, 1, Ydim, Xdim);
    }

    MultidimArray(int Xdim)
    {
        coreInit();
        coreAllocate(1UL, 1, 1, Xdim);
    }

    MultidimArray(const MultidimArray<T>& V)
    {
        coreInit();
        *this = V;
    }

    MultidimArray(MultidimArray<T>&& V) noexcept
    {
        coreInit();
        swap(V);
    }

    MultidimArray(const Matrix1D<T>& V)
    {
        coreInit();
        coreAllocate(1, 1, 1, V.size());
        for (size_t i = 0; i < V.size(); i++)
            DIRECT_A1D_ELEM(*this,i) = VEC_ELEM(V,i);
    }

    MultidimArray(const std::vector<T> &vector)
    {
        coreInit();
        coreAllocate(1, 1, 1, vector.size());
        for (size_t i = 0; i < vector.size(); i++)
            DIRECT_A1D_ELEM(*this,i) = vector[i];
    }


    virtual ~MultidimArray()
    {
        coreDeallocate();
    }

    void clear()
    {
        coreDeallocate();
        coreInit();
    }



    void coreInit() noexcept
    {
        xdim = yxdim = zyxdim = nzyxdim = 0;
        ydim = zdim = ndim = 0;
        zinit = yinit = xinit = 0;
        data = NULL;
        nzyxdimAlloc = 0;
        destroyData = true;
        mmapOn = false;
        mFd = NULL;
    }

    void swap(MultidimArray<T>& other) noexcept
    {
        std::swap(xdim, other.xdim);
        std::swap(ydim, other.ydim);
        std::swap(zdim, other.zdim);
        std::swap(ndim, other.ndim);
        std::swap(yxdim, other.yxdim);
        std::swap(zyxdim, other.zyxdim);
        std::swap(nzyxdim, other.nzyxdim);
        std::swap(xinit, other.xinit);
        std::swap(yinit, other.yinit);
        std::swap(zinit, other.zinit);
        std::swap(data, other.data);
        std::swap(nzyxdimAlloc, other.nzyxdimAlloc);
        std::swap(destroyData, other.destroyData);
        std::swap(mmapOn, other.mmapOn);
        std::swap(mFd, other.mFd);
    }

    void coreAllocate(size_t _ndim, int _zdim, int _ydim, int _xdim)
    {
        if (_ndim <= 0 || _zdim <= 0 || _ydim<=0 || _xdim<=0)
        {
            clear();
            return;
        }
        if(data!=NULL)
            REPORT_ERROR(ERR_MEM_NOTDEALLOC, "do not allocate space for an image if you have not deallocate it first");

        ndim=_ndim;
        zdim=_zdim;
        ydim=_ydim;
        xdim=_xdim;
        yxdim=(size_t)ydim*xdim;
        zyxdim=yxdim*zdim;
        nzyxdim=zyxdim*ndim;

        coreAllocate();
    }

    void coreAllocate()
    {
        if(data!=NULL)
            REPORT_ERROR(ERR_MEM_NOTDEALLOC, "do not allocate space for an image if you have not deallocate it first");

        if (mmapOn)
            mFd = mmapFile(data, nzyxdim);
        else
        {
            try
            {
                data = new T [nzyxdim];
                if (data == NULL)
                {
                    setMmap(true);
                    mFd = mmapFile(data, nzyxdim);
                }
            }
            catch (std::bad_alloc &)
            {
                setMmap(true);
                mFd = mmapFile(data, nzyxdim);
            }
        }
        memset(data,0,nzyxdim*sizeof(T));
        nzyxdimAlloc = nzyxdim;
    }

    void coreAllocateReuse()
    {
        if(data != NULL && nzyxdim <= nzyxdimAlloc)
            return;
        else if (nzyxdim > nzyxdimAlloc)
            coreDeallocate();

        if (mmapOn)
            mFd = mmapFile(data, nzyxdim);
        else
        {
            data = new T [nzyxdim];
            if (data == NULL)
                REPORT_ERROR(ERR_MEM_NOTENOUGH, "Allocate: No space left");
        }
        memset(data,0,nzyxdim*sizeof(T));
        nzyxdimAlloc = nzyxdim;
    }

    /* Create a temporary file to mmap the data.
     */
    FILE* mmapFile(T* &_data, size_t nzyxDim) const;

    void coreDeallocate() noexcept
    {
        if (data != NULL && destroyData)
        {
            if (mmapOn)
            {
#ifdef XMIPP_MMAP
                munmap(data,nzyxdimAlloc*sizeof(T));
                fclose(mFd);
#else

                REPORT_ERROR(ERR_MMAP,"Mapping not supported in Windows");
#endif

            }
            else
                delete[] data;
        }
        data = NULL;
        destroyData = true;
        nzyxdimAlloc = 0;
    }

    void alias(const MultidimArray<T> &m)
    {
        coreDeallocate();
        copyShape(m);
        this->data=m.data;
        this->nzyxdimAlloc = this->nzyxdim;
        this->destroyData = false;
    }

    void aliasRow(const MultidimArray<T> &m, size_t select_row)
    {
        if (select_row >= YSIZE(m))
            REPORT_ERROR(ERR_MULTIDIM_SIZE, "aliasRow: Selected row cannot be higher than Y size.");

        coreDeallocate();
        setDimensions(XSIZE(m),1, 1, 1);
        this->data = m.data + XSIZE(m)*select_row;
        this->nzyxdimAlloc = this->nzyxdim;
        this->destroyData = false;
    }

    void aliasSlice(const MultidimArray<T> &m, size_t select_slice)
    {
        if (select_slice >= ZSIZE(m))
            REPORT_ERROR(ERR_MULTIDIM_SIZE, "aliasSlice: Selected slice cannot be higher than Z size.");

        coreDeallocate();
        setDimensions(XSIZE(m), YSIZE(m), 1, 1);
        this->data = m.data + XSIZE(m)*YSIZE(m)*(select_slice);
        this->nzyxdimAlloc = this->nzyxdim;
        this->destroyData = false;
    }

    void aliasImageInStack(const MultidimArray<T> &m, size_t select_image)
    {
        if (select_image >= NSIZE(m))
            REPORT_ERROR(ERR_MULTIDIM_SIZE, "aliasImageInStack: Selected image cannot be higher than N size.");
        if (ZSIZE(m)!=1)
            REPORT_ERROR(ERR_MULTIDIM_SIZE, "aliasImageInStack: This function is not meant for volumes");

        coreDeallocate();
        setDimensions(XSIZE(m), YSIZE(m), 1, 1);
        this->data = m.data + XSIZE(m)*YSIZE(m)*(select_image);
        this->nzyxdimAlloc = this->nzyxdim;
        this->destroyData = false;
    }






    /* These "using" declarations must be done due to c++ cannot overload methods that have been
     * declared in base class.
     */
    using MultidimArrayBase::resizeNoCopy;
    using MultidimArrayBase::resize;
    using MultidimArrayBase::maxIndex;

    void resize(size_t Ndim, size_t Zdim, size_t Ydim, size_t Xdim, bool copy=true);

    template<typename T1>
    void resize(const MultidimArray<T1> &v, bool copy = true)
    {
        if (NSIZE(*this) != NSIZE(v) || XSIZE(*this) != XSIZE(v) ||
            YSIZE(*this) != YSIZE(v) || ZSIZE(*this) != ZSIZE(v) || data==NULL)
            resize(NSIZE(v), ZSIZE(v), YSIZE(v), XSIZE(v), copy);

        STARTINGX(*this) = STARTINGX(v);
        STARTINGY(*this) = STARTINGY(v);
        STARTINGZ(*this) = STARTINGZ(v);
    }

    template<typename T1>
    void resizeNoCopy(const MultidimArray<T1> &v)
    {
        if (NSIZE(*this) != NSIZE(v) || XSIZE(*this) != XSIZE(v) ||
            YSIZE(*this) != YSIZE(v) || ZSIZE(*this) != ZSIZE(v) || data==NULL)
            resize(NSIZE(v), ZSIZE(v), YSIZE(v), XSIZE(v), false);

        STARTINGX(*this) = STARTINGX(v);
        STARTINGY(*this) = STARTINGY(v);
        STARTINGZ(*this) = STARTINGZ(v);
    }

#define checkDimension(dim) checkDimensionWithDebug(dim,__FILE__,__LINE__)
    void checkDimensionWithDebug(int dim, const char *file, int line) const
    {
        if (getDim() != dim)
        {
            std::cerr<<" Check for dimension: "  << dim <<std::endl;
            std::cerr << "MultidimArray shape: ";
            printShape(std::cerr);
            std::cerr << std::endl;
            std::cerr << "Check called from file "<<file<<" line "<<line<<std::endl;
            exit(1);
        }
    }

    void selfWindow(int n0,int z0, int y0, int x0,
                    int nF,int zF, int yF, int xF,
                    T init_value = 0)
    {
        if (this->ndim >1)
            REPORT_ERROR(ERR_MULTIDIM_DIM,"stack windowing not implemented");
        if (this->zdim >1)
        {//call 3Dwindow
            selfWindow( z0,  y0,  x0,
                        zF,  yF,  xF,
                        init_value);
        }
        else if (this->ydim >1)
        {//call 2Dwindow
            selfWindow( y0,  x0,
                        yF,  xF,
                        init_value);

        }
        else if (this->xdim >1)
        {//call 1Dwindow
            selfWindow( x0, xF, init_value);
        }
    }

    template <class T1>
    void window(MultidimArray<T1> &result, int n0,int z0, int y0, int x0,
                int nF,int zF, int yF, int xF,
                T1 init_value = 0) const
    {
        if (this->ndim >1)
            REPORT_ERROR(ERR_MULTIDIM_DIM,"stack windowing not implemented");
        if (this->zdim >1)
        {//call 3Dwindow
            window(result, z0,  y0,  x0,
                   zF,  yF,  xF,
                   init_value);
        }
        else if (this->ydim >1)
        {//call 2Dwindow
            window(result, y0,  x0,
                   yF,  xF,
                   init_value);
        }
        else if (this->xdim >1)
        {//call 1Dwindow
            window(result, x0, xF, init_value);
        }
    }

    template <class T1>
    void window(MultidimArray<T1> &result, int z0, int y0, int x0, int zF, int yF, int xF,
                T1 init_value = 0) const
    {
        result.resizeNoCopy(zF - z0 + 1, yF - y0 + 1, xF - x0 + 1);
        STARTINGZ(result) = z0;
        STARTINGY(result) = y0;
        STARTINGX(result) = x0;

        FOR_ALL_ELEMENTS_IN_ARRAY3D(result)
        if ((k >= STARTINGZ(*this) && k <= FINISHINGZ(*this)) &&
            (i >= STARTINGY(*this) && i <= FINISHINGY(*this)) &&
            (j >= STARTINGX(*this) && j <= FINISHINGX(*this)))
            A3D_ELEM(result, k, i, j) = (T1) A3D_ELEM(*this, k, i, j);
        else
            A3D_ELEM(result, k, i, j) = init_value;
    }

    void selfWindow(int z0, int y0, int x0, int zF, int yF, int xF,
                    T init_value = 0)
    {
        if (z0 == STARTINGZ(*this) && zF == FINISHINGZ(*this) &&
            y0 == STARTINGY(*this) && yF == FINISHINGY(*this) &&
            x0 == STARTINGX(*this) && xF == FINISHINGX(*this))
            return;

        MultidimArray<T> result;
        window(result,z0,y0,x0,zF,yF,xF,init_value);
        *this=result;
    }

    template <class T1>
    void window(MultidimArray<T1> &result, int y0, int x0, int yF, int xF,
                T1 init_value = 0) const
    {
        result.resizeNoCopy(yF - y0 + 1, xF - x0 + 1);
        STARTINGY(result) = y0;
        STARTINGX(result) = x0;

        FOR_ALL_ELEMENTS_IN_ARRAY2D(result)
        if (j >= STARTINGX(*this) && j <= FINISHINGX(*this) &&
            i >= STARTINGY(*this) && i <= FINISHINGY(*this))
            A2D_ELEM(result, i, j) = A2D_ELEM(*this, i, j);
        else
            A2D_ELEM(result, i, j) = init_value;
    }

    void selfWindow(int y0, int x0, int yF, int xF,
                    T init_value = 0)
    {
        MultidimArray<T> result;
        window(result,y0,x0,yF,xF,init_value);
        *this=result;
    }

    template <class T1>
    void window(MultidimArray<T1> &result, int x0, int xF, T1 init_value = 0) const
    {
        result.resizeNoCopy(xF - x0 + 1);
        STARTINGX(result) = x0;

        FOR_ALL_ELEMENTS_IN_ARRAY1D(result)
        if (i >= STARTINGX(*this) && i <= FINISHINGX(*this))
            A1D_ELEM(result, i) = A1D_ELEM(*this, i);
        else
            A1D_ELEM(result, i) = init_value;
    }

    void selfWindow(int x0, int xF,
                    T init_value = 0)
    {
        MultidimArray<T> result;
        window(result,x0,xF,init_value);
        *this=result;
    }

    void patch(MultidimArray<T> patchArray, int x, int y)
    {
        int n = XSIZE(patchArray)*sizeof(T);

        for (size_t i=0; i < YSIZE(patchArray); ++i)
            memcpy(&dAij(*this, y+i, x), &dAij(patchArray, i, 0), n);
    }




    T& operator()(const Matrix1D< double >& v) const
    {
        switch (VEC_XSIZE(v))
        {
        case 1:
            return A1D_ELEM((*this), ROUND(XX(v)));
        case 2:
            return A2D_ELEM((*this), ROUND(YY(v)), ROUND(XX(v)));
        case 3:
            return A3D_ELEM((*this), ROUND(ZZ(v)), ROUND(YY(v)), ROUND(XX(v)));
        default:
            REPORT_ERROR(ERR_ARG_INCORRECT,"Cannot handle indexes with dimension larger than 3");
        }
    }

    T& operator()(const Matrix1D< int >& v) const
    {
        switch (VEC_XSIZE(v))
        {
        case 1:
            return A1D_ELEM((*this), XX(v));
        case 2:
            return A2D_ELEM((*this), YY(v), XX(v));
        case 3:
            return A3D_ELEM((*this), ZZ(v), YY(v), XX(v));
        default:
            REPORT_ERROR(ERR_ARG_INCORRECT,"Cannot handle indexes with dimension larger than 3");
        }
    }

    inline T& operator()(size_t n, int k, int i, int j) const
    {
        return NZYX_ELEM(*this, n, k, i, j);
    }

    inline T& operator()(int k, int i, int j) const
    {
        return A3D_ELEM(*this, k, i, j);
    }

    inline T& operator()(int i, int j) const
    {
        return A2D_ELEM(*this, i, j);
    }

    inline T& operator()(int i) const
    {
        return A1D_ELEM(*this, i);
    }

    inline T& operator[](size_t i) const
     {
         return data[i];
     }


    void* getArrayPointer() const
    {
        return (void*) data;

    }

    void getImage(size_t n, MultidimArray<T>& M, size_t n2 = 0) const
    {
        if (XSIZE(*this) == 0)
        {
            M.clear();
            return;
        }

        if (n > NSIZE(*this))
            REPORT_ERROR(ERR_INDEX_OUTOFBOUNDS," Multidimarray getImage: n larger than NSIZE");

        if (ZSIZE(*this) != ZSIZE(M) || YSIZE(*this) != YSIZE(M) || XSIZE(*this) != XSIZE(M))
        {
            if (n2 == 0)
                M.resizeNoCopy(*this);
            else
                REPORT_ERROR(ERR_MULTIDIM_SIZE, "MultidimArray::getImage: Target dimensions do not match source dimensions.");
        }

        if (n2 > NSIZE(M))
            REPORT_ERROR(ERR_INDEX_OUTOFBOUNDS," Multidimarray getImage: n larger than MultidimArray target NSIZE");


        FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(M)
        DIRECT_NZYX_ELEM(M, n2, k, i, j) = DIRECT_NZYX_ELEM(*this, n, k, i, j);

        STARTINGX(M) = STARTINGX(*this);
        STARTINGY(M) = STARTINGY(*this);
        STARTINGZ(M) = STARTINGZ(*this);
    }

    template <typename T1>
    void getSlice(int k, MultidimArray<T1>& M, char axis = 'Z', bool reverse = false, size_t n = 0) const
    {
        if (XSIZE(*this) == 0)
        {
            M.clear();
            return;
        }

        T* ptr=NULL;
        switch (axis)
        {
        case 'Z':
            if (k < STARTINGZ(*this) || k > FINISHINGZ(*this))
                REPORT_ERROR(ERR_INDEX_OUTOFBOUNDS,
                             "Slice: Multidim subscript (k) out of range");

            M.resize(1, 1, YSIZE(*this), XSIZE(*this),false);
            k = k - STARTINGZ(*this);

            if (reverse)
            {
                int zEnd = ZSIZE(*this) - 1;
                FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(M)
                DIRECT_A2D_ELEM(M, i, j) = (T1) DIRECT_NZYX_ELEM(*this, n, k,  zEnd-i, j);
            }
            else
            {
                ptr=&(DIRECT_NZYX_ELEM(*this, n, k, 0, 0));
                FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(M)
                DIRECT_MULTIDIM_ELEM(M, n) = (T1) *(ptr++);
            }

            STARTINGX(M) = STARTINGX(*this);
            STARTINGY(M) = STARTINGY(*this);
            break;
        case 'Y':
            if (k < STARTINGY(*this) || k > FINISHINGY(*this))
                REPORT_ERROR(ERR_INDEX_OUTOFBOUNDS,
                             "Slice: Multidim subscript (i) out of range");

            k = k - STARTINGY(*this);
            M.resizeNoCopy(ZSIZE(*this), XSIZE(*this));

            if (reverse)
            {
                int zEnd = ZSIZE(*this) - 1;
                FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(M)
                DIRECT_A2D_ELEM(M, i, j) = (T1) DIRECT_NZYX_ELEM(*this, n, zEnd-i, k, j);
            }
            else
            {
                FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(M)
                DIRECT_A2D_ELEM(M, i, j) = (T1) DIRECT_NZYX_ELEM(*this, n, i, k, j);
            }

            STARTINGX(M) = STARTINGX(*this);
            STARTINGY(M) = STARTINGZ(*this);
            break;
        case 'X':
            if (k < STARTINGX(*this) || k > FINISHINGX(*this))
                REPORT_ERROR(ERR_INDEX_OUTOFBOUNDS,
                             "Slice: Multidim subscript (j) out of range");

            k = k - STARTINGX(*this);
            M.resizeNoCopy(YSIZE(*this), ZSIZE(*this));

            if (reverse)
            {
                int zEnd = ZSIZE(*this) - 1;
                FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(M)
                DIRECT_A2D_ELEM(M, i, j) = (T1) DIRECT_NZYX_ELEM(*this, n, zEnd-j, i, k);
            }
            else
            {
                FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(M)
                DIRECT_A2D_ELEM(M, i, j) = (T1) DIRECT_NZYX_ELEM(*this, n, j, i, k);
            }

            STARTINGX(M) = STARTINGY(*this);
            STARTINGY(M) = STARTINGZ(*this);
            break;
        default:
            REPORT_ERROR(ERR_VALUE_INCORRECT,
                         formatString("Slice: not supported axis %c", axis));
        }
    }

    void getSliceAsMatrix(size_t k, Matrix2D<T> &m) const;


    void getAliasAsRowVector(Matrix1D<T> &m) const
    {
        m.vdim = NZYXSIZE(*this);
        m.destroyData = false;
        m.row = true;
        m.vdata = data;
    }

    template <typename T1>
    void setSlice(int k, const MultidimArray<T1>& v, size_t n=0)
    {
        if (xdim == 0)
            return;

        if (k < STARTINGZ(*this) || k > FINISHINGZ(*this))
            REPORT_ERROR(ERR_INDEX_OUTOFBOUNDS,formatString(
                             "setSlice: MultidimArray subscript (k=%d) out of range [%d, %d]",
                             k,STARTINGZ(*this),FINISHINGZ(*this)));

        if (v.rowNumber() != YSIZE(*this) || v.colNumber() != XSIZE(*this))
            REPORT_ERROR(ERR_MULTIDIM_DIM,
                         "setSlice: MultidimArray dimensions different from the matrix ones");

        k-=STARTINGZ(*this);
        T *ptr=&(DIRECT_NZYX_ELEM(*this, n, k, 0, 0));
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(v)
        *(ptr++) = (T) DIRECT_MULTIDIM_ELEM(v, n);
    }

    template <typename T1>
    void reslice(MultidimArray<T1>& out, AxisView face, bool flip = false, size_t n = 0) const
    {
        ArrayDim aDim, aDimOut;
        getDimensions(aDim);

        char axis='Z';
        bool reverse=false;

        aDimOut = aDim;

        if (face == VIEW_Y_NEG || face == VIEW_Y_POS)
        {
            axis = 'Y';
            aDimOut.ydim = aDim.zdim;
            aDimOut.zdim = aDim.ydim;
            reverse = (face == VIEW_Y_NEG);
        }
        else if (face == VIEW_X_NEG || face == VIEW_X_POS)
        {
            axis = 'X';
            aDimOut.xdim = aDim.zdim;
            aDimOut.zdim = aDim.xdim;
            reverse = (face == VIEW_X_NEG);
        }

        flip = flip^reverse;

        out.resize(aDimOut, false);

        MultidimArray<T1> imTemp;

        int index;

        for (size_t k = 0; k < aDimOut.zdim; k++)
        {
            imTemp.aliasSlice(out, k);
            index = k + (aDimOut.zdim - 1 - 2*k) * (int)flip;
            this->getSlice(index, imTemp, axis, !reverse);
        }
    }

    void reslice(AxisView face, bool flip = false, size_t n = 0)
    {
        MultidimArray<T> mTemp;
        reslice(mTemp, face, flip, n);
        *this = mTemp;
    }

    void getCol(size_t j, MultidimArray<T>& v) const
    {
        if (xdim == 0 || ydim == 0)
        {
            v.clear();
            return;
        }

        if (j >= xdim)
            REPORT_ERROR(ERR_INDEX_OUTOFBOUNDS,"getCol: Matrix subscript (j) greater than matrix dimension");

        v.resizeNoCopy(ydim);
        for (size_t i = 0; i < ydim; i++)
            DIRECT_A1D_ELEM(v,i) = DIRECT_A2D_ELEM(*this,i, j);
    }

    void setCol(size_t j, const MultidimArray<T>& v)
    {
        if (xdim == 0 || ydim == 0)
            REPORT_ERROR(ERR_MULTIDIM_EMPTY, "setCol: Target matrix is empty");

        if (j>= xdim)
            REPORT_ERROR(ERR_INDEX_OUTOFBOUNDS, "setCol: Matrix subscript (j) out of range");

        if (v.xdim != ydim)
            REPORT_ERROR(ERR_MULTIDIM_SIZE,
                         "setCol: Vector dimension different from matrix one");

        for (size_t i = 0; i < ydim; i++)
            DIRECT_A2D_ELEM(*this,i, j) = DIRECT_A1D_ELEM(v,i);
    }

    void getRow(size_t i, MultidimArray<T>& v) const
    {
        if (xdim == 0 || ydim == 0)
        {
            v.clear();
            return;
        }

        if (i >= ydim)
            REPORT_ERROR(ERR_INDEX_OUTOFBOUNDS, "getRow: Matrix subscript (i) greater than matrix dimension");

        v.resizeNoCopy(xdim);
        memcpy(&A1D_ELEM(v,0),&A2D_ELEM(*this,i,0),xdim*sizeof(T));
    }

    void setRow(int i, const MultidimArray<T>& v)
    {
        if (xdim == 0 || ydim == 0)
            REPORT_ERROR(ERR_MULTIDIM_EMPTY, "setRow: Target matrix is empty");

        if (i < 0 || i >= ydim)
            REPORT_ERROR(ERR_INDEX_OUTOFBOUNDS, "setRow: Matrix subscript (i) out of range");

        if (v.xdim != xdim)
            REPORT_ERROR(ERR_MULTIDIM_SIZE,
                         "setRow: Vector dimension different from matrix one");

        memcpy(&A2D_ELEM(*this,i,0),&A1D_ELEM(v,0),xdim*sizeof(T));
    }

    void getReal(MultidimArray< double > & realImg) const
    {
        REPORT_ERROR(ERR_TYPE_INCORRECT, "MultidimArray: Non complex datatype.");
    }

    void getImag(MultidimArray< double > & imagImg) const
    {
        REPORT_ERROR(ERR_TYPE_INCORRECT, "MultidimArray: Non complex datatype.");
    }


    void toPhysical(int k_log, int i_log, int j_log,
                    int& k_phys, int& i_phys, int& j_phys) const
    {
        k_phys = k_log - STARTINGZ(*this);
        i_phys = i_log - STARTINGY(*this);
        j_phys = j_log - STARTINGX(*this);
    }

    void toLogical(int k_phys, int i_phys, int j_phys,
                   int& k_log, int& i_log, int& j_log) const
    {
        k_log = k_phys + STARTINGZ(*this);
        i_log = i_phys + STARTINGY(*this);
        j_log = j_phys + STARTINGX(*this);
    }

    void toPhysical(int i_log, int j_log, int& i_phys, int& j_phys) const
    {
        i_phys = i_log - STARTINGY(*this);
        j_phys = j_log - STARTINGX(*this);
    }

    void toLogical(int i_phys, int j_phys, int &i_log, int& j_log) const
    {
        i_log = i_phys + STARTINGY(*this);
        j_log = j_phys + STARTINGX(*this);
    }

    void toPhysical(int i_log, int& i_phys) const
    {
        i_phys = i_log - STARTINGX(*this);
    }

    void toLogical(int i_phys, int& i_log) const
    {
        i_log = i_phys + STARTINGX(*this);
    }

    T interpolatedElement3D(double x, double y, double z, T outside_value = (T) 0) const
    {
        int x0 = FLOOR(x);
        double fx = x - x0;
        int x1 = x0 + 1;

        int y0 = FLOOR(y);
        double fy = y - y0;
        int y1 = y0 + 1;

        int z0 = FLOOR(z);
        double fz = z - z0;
        int z1 = z0 + 1;

        double doutside_value=outside_value;
        double d000 = (OUTSIDE3D(z0, y0, x0)) ? doutside_value : A3D_ELEM(*this, z0, y0, x0);
        double d001 = (OUTSIDE3D(z0, y0, x1)) ? doutside_value : A3D_ELEM(*this, z0, y0, x1);
        double d010 = (OUTSIDE3D(z0, y1, x0)) ? doutside_value : A3D_ELEM(*this, z0, y1, x0);
        double d011 = (OUTSIDE3D(z0, y1, x1)) ? doutside_value : A3D_ELEM(*this, z0, y1, x1);
        double d100 = (OUTSIDE3D(z1, y0, x0)) ? doutside_value : A3D_ELEM(*this, z1, y0, x0);
        double d101 = (OUTSIDE3D(z1, y0, x1)) ? doutside_value : A3D_ELEM(*this, z1, y0, x1);
        double d110 = (OUTSIDE3D(z1, y1, x0)) ? doutside_value : A3D_ELEM(*this, z1, y1, x0);
        double d111 = (OUTSIDE3D(z1, y1, x1)) ? doutside_value : A3D_ELEM(*this, z1, y1, x1);

        double dx00 = LIN_INTERP(fx, d000, d001);
        double dx01 = LIN_INTERP(fx, d100, d101);
        double dx10 = LIN_INTERP(fx, d010, d011);
        double dx11 = LIN_INTERP(fx, d110, d111);
        double dxy0 = LIN_INTERP(fy, dx00, dx10);
        double dxy1 = LIN_INTERP(fy, dx01, dx11);

        return (T) LIN_INTERP(fz, dxy0, dxy1);
    }

    T interpolatedElement2D(double x, double y, T outside_value = (T) 0) const
    {
        int x0 = floor(x);
        double fx = x - x0;
        int x1 = x0 + 1;
        int y0 = floor(y);
        double fy = y - y0;
        int y1 = y0 + 1;

        int i0=STARTINGY(*this);
        int j0=STARTINGX(*this);
        int iF=FINISHINGY(*this);
        int jF=FINISHINGX(*this);

#define ASSIGNVAL2D(d,i,j) \
     if ((j) < j0 || (j) > jF || (i) < i0 || (i) > iF) \
      d=outside_value;\
        else \
         d=A2D_ELEM(*this, i, j);

        double d00, d10, d11, d01;
        ASSIGNVAL2D(d00,y0,x0);
        ASSIGNVAL2D(d01,y0,x1);
        ASSIGNVAL2D(d10,y1,x0);
        ASSIGNVAL2D(d11,y1,x1);

        double d0 = LIN_INTERP(fx, d00, d01);
        double d1 = LIN_INTERP(fx, d10, d11);
        return (T) LIN_INTERP(fy, d0, d1);
    }

    T interpolatedElement2DOutsideZero(double x, double y) const
    {
        int x0 = floor(x);
        double fx = x - x0;
        int x1 = x0 + 1;
        int y0 = floor(y);
        double fy = y - y0;
        int y1 = y0 + 1;

        int i0=STARTINGY(*this);
        int j0=STARTINGX(*this);
        int iF=FINISHINGY(*this);
        int jF=FINISHINGX(*this);

        int b;
#define ASSIGNVAL2DNODIV(d,i,j) \
        b = ((j) < j0 || (j) > jF || (i) < i0 || (i) > iF); \
        b ? d=0 : d=A2D_ELEM(*this, i, j);

        double d00, d10, d11, d01;
        ASSIGNVAL2DNODIV(d00,y0,x0);
        ASSIGNVAL2DNODIV(d01,y0,x1);
        ASSIGNVAL2DNODIV(d10,y1,x0);
        ASSIGNVAL2DNODIV(d11,y1,x1);

        double d0 = LIN_INTERP(fx, d00, d01);
        double d1 = LIN_INTERP(fx, d10, d11);
        return (T) LIN_INTERP(fy, d0, d1);
    }

    T interpolatedElement1D(double x, T outside_value = (T) 0) const
    {
        int x0 = floor(x);
        double fx = x - x0;
        int x1 = x0 + 1;

        int j0=STARTINGX(*this);
        int jF=FINISHINGX(*this);

#define ASSIGNVAL1D(d,j) \
     if ((j) < j0 || (j) > jF) \
      d=outside_value;\
        else \
         d=A1D_ELEM(*this, j);

        double d0, d1;
        ASSIGNVAL1D(d0,x0);
        ASSIGNVAL1D(d1,x1);

        return (T) LIN_INTERP(fx, d0, d1);
    }

    T interpolatedElementBSpline3D(double x, double y, double z,
                                   int SplineDegree = 3) const;

    T interpolatedElementBSpline2D(double x, double y, int SplineDegree = 3) const;

    T interpolatedElementBSpline2D_Degree3(double x, double y) const;

    T interpolatedElementBSpline1D(double x, int SplineDegree = 3) const;



    void printStats(std::ostream& out = std::cout) const
    {
        T minval, maxval;
        double avgval, devval;

        computeStats(avgval, devval, minval, maxval);

        out.setf(std::ios::showpoint);
        int old_prec = out.precision(7);

        out << " min= ";
        out.width(9);
        out << minval;
        out << " max= ";
        out.width(9);
        out << maxval;
        out << " avg= ";
        out.width(9);
        out << avgval;
        out << " dev= ";
        out.width(9);
        out << devval;

        out.precision(old_prec);
    }

    T computeMax() const
    {
        if (NZYXSIZE(*this) <= 0)
            return static_cast< T >(0);

        T maxval = data[0];

        T* ptr=NULL;
        size_t n;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        if (*ptr > maxval)
            maxval = *ptr;

        return maxval;
    }

    void maxIndex(int& jmax) const
    {
        size_t zeroLong=0;
        int zeroInt=0;
        maxIndex(zeroLong,zeroInt,zeroInt,jmax);
    }

    T computeMin() const
    {
        if (NZYXSIZE(*this) <= 0)
            return static_cast< T >(0);

        T minval = data[0];

        T* ptr=NULL;
        size_t n;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        if (*ptr < minval)
            minval = *ptr;

        return minval;
    }

    void minIndex(int &lmin, int& kmin, int& imin, int& jmin) const
    {
        if (XSIZE(*this) == 0)
        {
            lmin = kmin = imin = jmin = -1;
            return;
        }

        kmin = STARTINGZ(*this);
        imin = STARTINGY(*this);
        jmin = STARTINGX(*this);
        lmin = 0;
        size_t n=0;
        T minval = DIRECT_MULTIDIM_ELEM(*this, n);

        FOR_ALL_NZYX_ELEMENTS_IN_MULTIDIMARRAY(*this)
        {
            T val=DIRECT_MULTIDIM_ELEM(*this,n);
            if (val > minval)
            {
                minval = val;
                lmin = l;
                kmin = k;
                imin = i;
                jmin = j;
            }
            ++n;
        }
    }

    void minIndex(int& kmin, int& imin, int& jmin) const
    {
        int zeroInt=0;
        minIndex(zeroInt,kmin,imin,jmin);
    }

    void minIndex(int& imin, int& jmin) const
    {
        int zeroInt=0;
        minIndex(zeroInt,zeroInt,imin,jmin);
    }

    void minIndex(int& jmin) const
    {
        int zeroInt=0;
        minIndex(zeroInt,zeroInt,zeroInt,jmin);
    }

    void maxIndex(size_t &lmax, int& kmax, int& imax, int& jmax) const
    {
        if (XSIZE(*this) == 0)
        {
            lmax = kmax = imax = jmax = -1;
            return;
        }

        kmax = STARTINGZ(*this);
        imax = STARTINGY(*this);
        jmax = STARTINGX(*this);
        lmax = 0;
        size_t n=0;
        T maxval = DIRECT_MULTIDIM_ELEM(*this, n);

        FOR_ALL_NZYX_ELEMENTS_IN_MULTIDIMARRAY(*this)
        {
            T val=DIRECT_MULTIDIM_ELEM(*this, n);
            if (val > maxval)
            {
                maxval = val;
                lmax = l;
                kmax = k;
                imax = i;
                jmax = j;
            }
            ++n;
        }
    }

    void computeDoubleMinMax(double& minval, double& maxval) const
    {
        if (NZYXSIZE(*this) <= 0)
            return;

        T* ptr=NULL;
        size_t n;
        T Tmin=DIRECT_MULTIDIM_ELEM(*this,0);
        T Tmax=Tmin;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        {
            T val=*ptr;
            if (val < Tmin)
                Tmin = val;
            else if (val > Tmax)
                Tmax = val;
        }
        minval = static_cast< double >(Tmin);
        maxval = static_cast< double >(Tmax);
    }

    void computeDoubleMinMaxRange(double& minval, double& maxval,size_t offset, size_t size) const
    {
        if (NZYXSIZE(*this) <= 0)
            return;

        minval = maxval = static_cast< double >(data[offset]);

        T* ptr=NULL;
        T val;
        size_t n;

        for (n=offset,ptr=data+offset; n<size; ++n, ++ptr)
        {
            val = *ptr;
            if (val < minval)
                minval = static_cast< double >(val);
            else if (val > maxval)
                maxval = static_cast< double >(val);
        }
    }

    double computeAvg() const
    {
        if (NZYXSIZE(*this) <= 0)
            return 0;

        double sum = 0;

        T* ptr=NULL;
        size_t n;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        sum += static_cast< double >(*ptr);

        return sum / NZYXSIZE(*this);
    }

    double computeStddev() const
    {
        if (NZYXSIZE(*this) <= 1)
            return 0;

        double avg = 0, stddev = 0;

        T* ptr=NULL;
        size_t n;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        {
            double val=static_cast< double >(*ptr);
            avg += val;
            stddev += val * val;
        }

        avg /= NZYXSIZE(*this);
        stddev = stddev / NZYXSIZE(*this) - avg * avg;
        stddev *= NZYXSIZE(*this) / (NZYXSIZE(*this) - 1);

        // Foreseeing numerical instabilities
        stddev = sqrt(static_cast<double>((ABS(stddev))));

        return stddev;
    }

    void computeStats(double& avg, double& stddev, T& minval, T& maxval) const
    {
        if (NZYXSIZE(*this) <= 0)
            return;

        avg = 0;
        stddev = 0;

        minval = maxval = data[0];

        T* ptr=NULL;
        size_t n;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        {
            T Tval=*ptr;
            double val=Tval;
            avg += val;
            stddev += val * val;

            if (Tval > maxval)
                maxval = Tval;
            else if (Tval < minval)
                minval = Tval;
        }

        avg /= NZYXSIZE(*this);

        if (NZYXSIZE(*this) > 1)
        {
            stddev = stddev / NZYXSIZE(*this) - avg * avg;
            stddev *= NZYXSIZE(*this) / (NZYXSIZE(*this) - 1);

            // Foreseeing numerical instabilities
            stddev = sqrt(static_cast< double >(ABS(stddev)));
        }
        else
            stddev = 0;
    }

    template<typename U>
    void computeAvgStdev(U& avg, U& stddev) const
    {
        static_assert(
                std::is_same<double, U>::value || std::is_same<float, U>::value,
                "U must be a floating presiont type");
        if (NZYXSIZE(*this) <= 0)
            return;

        avg = 0;
        stddev = 0;
        const size_t nMax = nzyxdim;
        for (size_t n = 0; n < nMax; ++n)
        {
            U v = (U)data[n];
            avg += v;
            stddev += v * v;
        }

        avg /= NZYXSIZE(*this);

        if (NZYXSIZE(*this) > 1)
        {
            stddev = stddev / NZYXSIZE(*this) - avg * avg;
            stddev *= NZYXSIZE(*this) / (NZYXSIZE(*this) - 1);

            // Foreseeing numerical instabilities
            stddev = sqrt(fabs(stddev));
        }
        else
            stddev = 0;
    }

    void computeAvgStdev_within_binary_mask(const MultidimArray< int >& mask,
                                            double& avg, double& stddev) const
    {
        double sum1 = 0;
        double sum2 = 0;
        double N = 0;

        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(*this)
        {
            if (DIRECT_MULTIDIM_ELEM(mask, n) != 0)
            {
                ++N;
                double aux=DIRECT_MULTIDIM_ELEM(*this, n);
                sum1 += aux;
                sum2 += aux*aux;
            }
        }

        // average and standard deviation
        avg  = sum1 / N;
        if (N > 1)
            stddev = sqrt(fabs(sum2 / N - avg * avg) * N / (N - 1));
        else
            stddev = 0;
    }

    void computeMedian_within_binary_mask(const MultidimArray< int >& mask, double& median) const;

    void computeStats(double& avg,
                      double& stddev,
                      T& min_val,
                      T& max_val,
                      Matrix1D< int >& corner1,
                      Matrix1D< int >& corner2,
                      size_t n = 0)
    {
        (*this).checkDimension(2);
        min_val = max_val = NZYX_ELEM((*this), n, 0, YY(corner1), XX(corner1));

        Matrix1D< double > r(3);
        double N = 0, sum = 0, sum2 = 0;

        FOR_ALL_ELEMENTS_IN_ARRAY2D_BETWEEN(corner1, corner2)
        {
            sum += (*this)(r);
            sum2 += (*this)(r) * (*this)(r);
            N++;

            if ((*this)(r) < min_val)
                min_val = (*this)(r);
            else if ((*this)(r) > max_val)
                max_val = (*this)(r);
        }

        if (N != 0)
        {
            avg = sum / N;
            stddev = sqrt(sum2 / N - avg * avg);
        }
        else
        {
            avg = stddev = 0;
        }
    }

    double computeMedian() const
    {
        if (XSIZE(*this) == 0)
            return 0;

        if (XSIZE(*this) == 1)
            return DIRECT_MULTIDIM_ELEM(*this,0);

        // Initialise data
        MultidimArray<T> temp;
        this->sort(temp);

        // Get median
        if (NZYXSIZE(*this)%2==0)
            return 0.5*(DIRECT_MULTIDIM_ELEM(temp,NZYXSIZE(*this)/2-1)+
                        DIRECT_MULTIDIM_ELEM(temp,NZYXSIZE(*this)/2  ));
        else
            return DIRECT_MULTIDIM_ELEM(temp,NZYXSIZE(*this)/2);
    }

    void rangeAdjust(T minF, T maxF)
    {
        if (NZYXSIZE(*this) <= 0)
            return;

        double min0=0., max0=0.;
        computeDoubleMinMax(min0, max0);

        // If max0==min0, it means that the vector is a constant one, so the
        // only possible transformation is to a fixed minF
        double slope;
        if (max0 != min0)
            slope = static_cast< double >(maxF - minF) /
                    static_cast< double >(max0 - min0);
        else
            slope = 0;

        T* ptr=NULL;
        size_t n;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        *ptr = minF + static_cast< T >(slope *
                                       static_cast< double >(*ptr - min0));
    }

    // This function must be explicitly implemented outside
    void rangeAdjust(T minF, T maxF, MultidimArray<int> &mask)
    {
        if (MULTIDIM_SIZE(*this) <= 0)
            return;

        double min0=0., max0=0.;
        bool first=true;
        T* ptr=NULL;
        size_t n;
        int * ptrMask=MULTIDIM_ARRAY(mask);
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        {
            if (*ptrMask)
            {
                T val= *ptr;
                if (first)
                {
                    min0=max0=(double)val;
                    first=false;
                }
                else
                {
                    min0=XMIPP_MIN(min0,val);
                    max0=XMIPP_MAX(max0,val);
                }
            }
            ptrMask++;
        }

        // If max0==min0, it means that the vector is a constant one, so the
        // only possible transformation is to a fixed minF
        double slope;
        if (max0 != min0)
            slope = static_cast< double >(maxF - minF) /
                    static_cast< double >(max0 - min0);
        else
            slope = 0;

        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        *ptr = minF + static_cast< T >(slope *
                                       static_cast< double >(*ptr - min0));
    }

    //As written this will only work for T=double
    //nevertheless since this is used is better
    //to use T than double or will create problem for int multidim arrays
    void rangeAdjust(const MultidimArray<T> &example,
                     const MultidimArray<int> *mask=NULL)
    {
        if (NZYXSIZE(*this) <= 0)
            return;

        double avgExample=0., stddevExample=0., avgThis, stddevThis;
        if (mask!=NULL)
        {
            example.computeAvgStdev_within_binary_mask(*mask,avgExample,stddevExample);
            computeAvgStdev_within_binary_mask(*mask,avgThis,stddevThis);
        }
        else
        {
            computeAvgStdev(avgThis,stddevThis);
            example.computeAvgStdev(avgExample,stddevExample);
        }

        // y=a+bx
        double b=stddevThis>0? stddevExample/stddevThis:0;
        double a=avgExample-avgThis*b;

        size_t n;
        T *ptr=NULL;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        *ptr = static_cast< T >(a+b * static_cast< double > (*ptr));
    }

    // This function must be explicitly implemented outside.
    template<typename U>
    void statisticsAdjust(U avgF, U stddevF)
    {
        static_assert(
                std::is_same<double, U>::value || std::is_same<float, U>::value,
                "U must be a floating presiont type");
        U avg0 = 0;
        U stddev0 = 0;

        if (NZYXSIZE(*this) == 0)
            return;

        computeAvgStdev(avg0, stddev0);

        U a = (stddev0 != 0) ? (stddevF / stddev0) : 0;
        U b = avgF - a * avg0;

        T* ptr=&DIRECT_MULTIDIM_ELEM(*this,0);
        size_t nmax=(nzyxdim/4)*4;
        for (size_t n=0; n<nmax; n+=4, ptr+=4)
        {
            *(ptr  )= static_cast< T >(a * (*(ptr  )) + b);
            *(ptr+1)= static_cast< T >(a * (*(ptr+1)) + b);
            *(ptr+2)= static_cast< T >(a * (*(ptr+2)) + b);
            *(ptr+3)= static_cast< T >(a * (*(ptr+3)) + b);
        }
        for (size_t n=nmax; n<nzyxdim; ++n, ptr+=1)
            *(ptr  )= static_cast< T >(a * (*(ptr  )) + b);
    }


    inline friend void coreArrayByArray(const MultidimArray<T>& op1,
                                        const MultidimArray<T>& op2, MultidimArray<T>& result,
                                        char operation)
    {
        T* ptrResult=NULL;
        T* ptrOp1=NULL;
        T* ptrOp2=NULL;
        size_t n;
        // Loop unrolling
        const size_t unroll=4;
        size_t nmax=unroll*(op1.nzyxdim/unroll);
        switch (operation)
        {
        case '+':
            for (n=0, ptrResult=result.data, ptrOp1=op1.data,ptrOp2=op2.data;
                 n<nmax; n+=unroll, ptrResult+=unroll, ptrOp1+=unroll, ptrOp2+=unroll)
            {
                *ptrResult = *ptrOp1 + *ptrOp2;
                *(ptrResult+1) = *(ptrOp1+1) + *(ptrOp2+1);
                *(ptrResult+2) = *(ptrOp1+2) + *(ptrOp2+2);
                *(ptrResult+3) = *(ptrOp1+3) + *(ptrOp2+3);
            }
            for (n=nmax, ptrResult=result.data+nmax, ptrOp1=op1.data+nmax, ptrOp2=op2.data+nmax;
                 n<op1.zyxdim; ++n, ++ptrResult, ++ptrOp1, ++ptrOp2)
                *ptrResult = *ptrOp1 + *ptrOp2;
            break;
        case '-':
                for (n=0, ptrResult=result.data, ptrOp1=op1.data,ptrOp2=op2.data;
                     n<nmax; n+=unroll, ptrResult+=unroll, ptrOp1+=unroll, ptrOp2+=unroll)
                {
                    *ptrResult = *ptrOp1 - *ptrOp2;
                    *(ptrResult+1) = *(ptrOp1+1) - *(ptrOp2+1);
                    *(ptrResult+2) = *(ptrOp1+2) - *(ptrOp2+2);
                    *(ptrResult+3) = *(ptrOp1+3) - *(ptrOp2+3);
                }
            for (n=nmax, ptrResult=result.data+nmax, ptrOp1=op1.data+nmax, ptrOp2=op2.data+nmax;
                 n<op1.zyxdim; ++n, ++ptrResult, ++ptrOp1, ++ptrOp2)
                *ptrResult = *ptrOp1 - *ptrOp2;
            break;
        case '*':
                for (n=0, ptrResult=result.data, ptrOp1=op1.data,ptrOp2=op2.data;
                     n<nmax; n+=unroll, ptrResult+=unroll, ptrOp1+=unroll, ptrOp2+=unroll)
                {
                    *ptrResult = *ptrOp1 * *ptrOp2;
                    *(ptrResult+1) = *(ptrOp1+1) * *(ptrOp2+1);
                    *(ptrResult+2) = *(ptrOp1+2) * *(ptrOp2+2);
                    *(ptrResult+3) = *(ptrOp1+3) * *(ptrOp2+3);
                }
            for (n=nmax, ptrResult=result.data+nmax, ptrOp1=op1.data+nmax, ptrOp2=op2.data+nmax;
                 n<op1.zyxdim; ++n, ++ptrResult, ++ptrOp1, ++ptrOp2)
                *ptrResult = *ptrOp1 * *ptrOp2;
            break;
        case '/':
                for (n=0, ptrResult=result.data, ptrOp1=op1.data,ptrOp2=op2.data;
                     n<nmax; n+=unroll, ptrResult+=unroll, ptrOp1+=unroll, ptrOp2+=unroll)
                {
                    *ptrResult = *ptrOp1 / *ptrOp2;
                    *(ptrResult+1) = *(ptrOp1+1) / *(ptrOp2+1);
                    *(ptrResult+2) = *(ptrOp1+2) / *(ptrOp2+2);
                    *(ptrResult+3) = *(ptrOp1+3) / *(ptrOp2+3);
                }
            for (n=nmax, ptrResult=result.data+nmax, ptrOp1=op1.data+nmax, ptrOp2=op2.data+nmax;
                 n<op1.zyxdim; ++n, ++ptrResult, ++ptrOp1, ++ptrOp2)
                *ptrResult = *ptrOp1 / *ptrOp2;
            break;
        }
    }

    inline friend void selfCoreArrayByArrayMask(const MultidimArray<T>& op1,
                                            const MultidimArray<T>& op2, MultidimArray<T>& result,
                                            char operation, const MultidimArray<T>* mask)
    {
        T* ptrResult=NULL;
        T* ptrOp1=NULL;
        T* ptrOp2=NULL;
        T* ptrMask=NULL;
        T zero;
        zero = T(0);
        size_t n;
        switch (operation)
        {
        case '+':
            for (n=0, ptrResult=result.data, ptrOp1=op1.data,ptrOp2=op2.data,
                 ptrMask=mask->data;
                 n<op1.nzyxdim; n++, ptrResult++,
                 ptrOp1++, ptrOp2++, ptrMask++)
            {
                if ((*ptrMask)==zero)
                    *ptrResult += (*ptrOp1);
                else
                    *ptrResult += (*ptrOp2);//(*ptrOp2) * (*ptrMask);
            }
            break;
        case '-':
            for (n=0, ptrResult=result.data, ptrOp1=op1.data,ptrOp2=op2.data,
                 ptrMask=mask->data;
                 n<op1.nzyxdim; n++, ptrResult++,
                 ptrOp1++, ptrOp2++, ptrMask++)
            {
                if ((*ptrMask)==zero)
                    *ptrResult -= (*ptrOp1);
                else
                    *ptrResult -= (*ptrOp2) * (*ptrMask);
            }
            break;
            //NOTE: not clear if this should be the default case for * and / operators
        case '*':
            for (n=0, ptrResult=result.data, ptrOp1=op1.data,ptrOp2=op2.data,
                 ptrMask=mask->data;
                 n<op1.nzyxdim; n++, ptrResult++,
                 ptrOp1++, ptrOp2++, ptrMask++)
            {
                if ((*ptrMask)==zero)
                    *ptrResult *= (*ptrOp1);
                else
                    *ptrResult *= (*ptrOp2) * (*ptrMask);
            }
            break;
        case '/':
            for (n=0, ptrResult=result.data, ptrOp1=op1.data,ptrOp2=op2.data,
                 ptrMask=mask->data;
                 n<op1.nzyxdim; n++, ptrResult++,
                 ptrOp1++, ptrOp2++, ptrMask++)
            {
                if ((*ptrMask)==zero)
                    *ptrResult /= (*ptrOp1);
                else
                    *ptrResult /= (*ptrOp2) * (*ptrMask);
            }
            break;
        }
    }


    inline friend void arrayByArray(const MultidimArray<T>& op1,
                                    const MultidimArray<T>& op2,
                                    MultidimArray<T>& result,
                                    char operation)
    {
        if (!op1.sameShape(op2))
            REPORT_ERROR(ERR_MULTIDIM_SIZE,
                         formatString("Array_by_array: different shapes (%c)", operation));
        if (result.data == NULL || !result.sameShape(op1))
            result.resizeNoCopy(op1);
        coreArrayByArray(op1, op2, result, operation);
    }
    inline friend void selfArrayByArrayMask(const MultidimArray<T>& op1,
                                            const MultidimArray<T>& op2,
                                            MultidimArray<T>& result,
                                            char operation, const MultidimArray<T> *mask=NULL)
    {
        if (!op1.sameShape(op2) || !result.sameShape(op1))
            REPORT_ERROR(ERR_MULTIDIM_SIZE,
                         formatString("Array_by_array: different shapes (%c)", operation));
        selfCoreArrayByArrayMask(op1, op2, result, operation, mask);
    }

    MultidimArray<T> operator+(const MultidimArray<T>& op1) const
    {
        MultidimArray<T> tmp;
        arrayByArray(*this, op1, tmp, '+');
        return tmp;
    }

    MultidimArray<T> operator-(const MultidimArray<T>& op1) const
    {
        MultidimArray<T> tmp;
        arrayByArray(*this, op1, tmp, '-');
        return tmp;
    }

    MultidimArray<T> operator*(const MultidimArray<T>& op1) const
    {
        MultidimArray<T> tmp;
        arrayByArray(*this, op1, tmp, '*');
        return tmp;
    }

    MultidimArray<T> operator/(const MultidimArray<T>& op1) const
    {
        MultidimArray<T> tmp;
        arrayByArray(*this, op1, tmp, '/');
        return tmp;
    }

    void operator+=(const MultidimArray<T>& op1)
    {
        arrayByArray(*this, op1, *this, '+');
    }

    void operator-=(const MultidimArray<T>& op1)
    {
        arrayByArray(*this, op1, *this, '-');
    }

    void operator*=(const MultidimArray<T>& op1)
    {
        arrayByArray(*this, op1, *this, '*');
    }

    void operator/=(const MultidimArray<T>& op1)
    {
        arrayByArray(*this, op1, *this, '/');
    }

    double dotProduct(const MultidimArray<T>& op1)
    {
        if (!sameShape(op1))
            REPORT_ERROR(ERR_MULTIDIM_SIZE,"The two arrays for dot product are not of the same shape");
        double dot=0;
        size_t n;
        T* ptrOp1=NULL;
        T* ptrOp2=NULL;
        const size_t unroll=4;
        size_t nmax=unroll*(MULTIDIM_SIZE(*this)/unroll);
        for (n=0, ptrOp1=data,ptrOp2=op1.data;
             n<nmax; n+=unroll, ptrOp1+=unroll, ptrOp2+=unroll)
        {
            dot += *ptrOp1 * *ptrOp2;
            dot += *(ptrOp1+1) * *(ptrOp2+1);
            dot += *(ptrOp1+2) * *(ptrOp2+2);
            dot += *(ptrOp1+3) * *(ptrOp2+3);
        }
        for (n=nmax, ptrOp1=data+nmax, ptrOp2=op1.data+nmax;
             n<MULTIDIM_SIZE(*this); ++n, ++ptrOp1, ++ptrOp2)
            dot += *ptrOp1 * *ptrOp2;
        return dot;
    }


    inline friend void coreArrayByScalar(const MultidimArray<T>& op1,
                                         const T& op2,
                                         MultidimArray<T>& result,
                                         char operation)
    {
        T* ptrResult=NULL;
        T* ptrOp1=NULL;
        size_t n;
        switch (operation)
        {
        case '+':
            for (n=0, ptrResult=result.data, ptrOp1=op1.data;
                 n<op1.zyxdim; ++n, ++ptrResult, ++ptrOp1)
                *ptrResult = *ptrOp1 + op2;
            break;
        case '-':
                for (n=0, ptrResult=result.data, ptrOp1=op1.data;
                     n<op1.zyxdim; ++n, ++ptrResult, ++ptrOp1)
                    *ptrResult = *ptrOp1 - op2;
            break;
        case '*':
                for (n=0, ptrResult=result.data, ptrOp1=op1.data;
                     n<op1.zyxdim; ++n, ++ptrResult, ++ptrOp1)
                    *ptrResult = *ptrOp1 * op2;
            break;
        case '/':
                for (n=0, ptrResult=result.data, ptrOp1=op1.data;
                     n<op1.zyxdim; ++n, ++ptrResult, ++ptrOp1)
                    *ptrResult = *ptrOp1 / op2;
            break;
        case '=':
                for (n=0, ptrResult=result.data, ptrOp1=op1.data;
                     n<op1.zyxdim; ++n, ++ptrResult, ++ptrOp1)
                    *ptrResult = *ptrOp1 == op2;
            break;
        }
    }

    inline friend void arrayByScalar(const MultidimArray<T>& op1,
                                     T op2,
                                     MultidimArray<T>& result,
                                     char operation)
    {
        if (result.data == NULL || !result.sameShape(op1))
            result.resizeNoCopy(op1);
        coreArrayByScalar(op1, op2, result, operation);
    }

    MultidimArray<T> operator+(T op1) const
    {
        MultidimArray<T> tmp;
        arrayByScalar(*this, op1, tmp, '+');
        return tmp;
    }

    MultidimArray<T> operator-(T op1) const
    {
        MultidimArray<T> tmp;
        arrayByScalar(*this, op1, tmp, '-');
        return tmp;
    }

    MultidimArray<T> operator*(T op1) const
    {
        MultidimArray<T> tmp;
        arrayByScalar(*this, op1, tmp, '*');
        return tmp;
    }

    MultidimArray<T> operator/(T op1) const
    {
        MultidimArray<T> tmp;
        arrayByScalar(*this, op1, tmp, '/');
        return tmp;
    }


    void equal(T op1, MultidimArray<char> &result) const
    {
        result.resizeNoCopy(*this);
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(result)
        DIRECT_MULTIDIM_ELEM(result,n) = DIRECT_MULTIDIM_ELEM(*this,n) == op1;
    }

    void operator+=(const T& op1)
    {
        arrayByScalar(*this, op1, *this, '+');
    }

    void operator-=(const T& op1)
    {
        arrayByScalar(*this, op1, *this, '-');
    }

    void operator*=(const T& op1)
    {
        arrayByScalar(*this, op1, *this, '*');
    }

    void operator/=(const T& op1)
    {
        arrayByScalar(*this, op1, *this, '/');
    }


    inline friend void coreScalarByArray(const T& op1,
                                         const MultidimArray<T>& op2,
                                         MultidimArray<T>& result,
                                         char operation)
    {
        T* ptrResult=NULL;
        T* ptrOp2=NULL;
        size_t n;
        for (n=0, ptrResult=result.data, ptrOp2=op2.data;
             n<op2.zyxdim; ++n, ++ptrResult, ++ptrOp2)
            switch (operation)
            {
            case '+':
                *ptrResult = op1 + *ptrOp2;
                break;
            case '-':
                *ptrResult = op1 - *ptrOp2;
                break;
            case '*':
                *ptrResult = op1 * *ptrOp2;
                break;
            case '/':
                *ptrResult = op1 / *ptrOp2;
                break;
            }
    }

    inline friend void scalarByArray(T op1,
                                     const MultidimArray<T>& op2,
                                     MultidimArray<T>& result,
                                     char operation)
    {
        if (result.data == NULL || !result.sameShape(op2))
            result.resizeNoCopy(op2);
        coreScalarByArray(op1, op2, result, operation);
    }

    friend MultidimArray<T> operator+(T op1, const MultidimArray<T>& op2)
    {
        MultidimArray<T> tmp;
        scalarByArray(op1, op2, tmp, '+');
        return tmp;
    }

    friend MultidimArray<T> operator-(T op1, const MultidimArray<T>& op2)
    {
        MultidimArray<T> tmp;
        scalarByArray(op1, op2, tmp, '-');
        return tmp;
    }

    friend MultidimArray<T> operator*(T op1, const MultidimArray<T>& op2)
    {
        MultidimArray<T> tmp;
        scalarByArray(op1, op2, tmp, '*');
        return tmp;
    }

    friend MultidimArray<T> operator/(T op1, const MultidimArray<T>& op2)
    {
        MultidimArray<T> tmp;
        scalarByArray(op1, op2, tmp, '/');
        return tmp;
    }


    void initConstant(T val)
    {
        T* ptr=NULL;
        size_t n;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        *ptr = val;
    }

    template <typename T1>
    void initZeros(const MultidimArray<T1>& op)
    {
        if (data == NULL || !sameShape(op))
            resizeNoCopy(op);
        memset(data,0,nzyxdim*sizeof(T));
    }

    inline void initZeros()
    {
        memset(data,0,nzyxdim*sizeof(T));
    }

    inline void initZeros(size_t Ndim, size_t Zdim, size_t Ydim, size_t Xdim)
    {
        if (xdim!=Xdim || ydim!=Ydim || zdim!=Zdim || ndim!=Ndim)
            resize(Ndim, Zdim,Ydim,Xdim,false);
        memset(data,0,nzyxdim*sizeof(T));
    }

    void initZeros(int Xdim)
    {
        initZeros(1, 1, 1, Xdim);
    }

    void initZeros(int Ydim, int Xdim)
    {
        initZeros(1, 1, Ydim, Xdim);
    }

    void initZeros(int Zdim, int Ydim, int Xdim)
    {
        initZeros(1, Zdim, Ydim, Xdim);
    }

    void initLinear(T minF, T maxF, int n = 1, const String& mode = "incr")
    {
        double slope;
        int steps;

        if (mode == "incr")
        {
            steps = 1 + (int) FLOOR((double) ABS((maxF - minF)) / ((double) n));
            slope = n * SGN(maxF - minF);
        }
        else if (mode == "steps")
        {
            steps = n;
            slope = (maxF - minF) / (steps - 1);
        }
        else
            REPORT_ERROR(ERR_VALUE_INCORRECT, "Init_linear: Mode not supported (" + mode +
                         ")");

        if (steps == 0)
            clear();
        else
        {
            resizeNoCopy(steps);
            for (int i = 0; i < steps; i++)
                A1D_ELEM(*this, i) = (T)((double) minF + slope * i);
        }
    }

    void initRandom(double op1, double op2, RandomMode mode = RND_UNIFORM);

    void addNoise(double op1,
                  double op2,
                  const String& mode = "uniform",
                  double df = 3.) const;


    T*** adaptForNumericalRecipes3D(size_t n = 0) const
    {
        T*** m = NULL;
        ask_Tvolume(m, 1, ZSIZE(*this), 1, YSIZE(*this), 1, XSIZE(*this));

        FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(*this)
        m[k+1][i+1][j+1] = DIRECT_NZYX_ELEM(*this, n, k, i, j);

        return m;
    }

    void killAdaptationForNumericalRecipes3D(T*** m) const
    {
        free_Tvolume(m, 1, ZSIZE(*this), 1, YSIZE(*this), 1, XSIZE(*this));
    }

    T** adaptForNumericalRecipes2D(size_t n = 0) const
    {
        T** m = NULL;
        ask_Tmatrix(m, 1, YSIZE(*this), 1, XSIZE(*this));

        FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(*this)
        m[i+1][j+1] = DIRECT_NZYX_ELEM(*this, n, 0, i, j);

        return m;
    }

    T* adaptForNumericalRecipes22D() const
    {
        return MULTIDIM_ARRAY(*this) - 1 - XSIZE(*this);
    }

    void loadFromNumericalRecipes2D(T** m, int Ydim, int Xdim)
    {
        resizeNoCopy(Ydim, Xdim);

        for (int i = 1; i <= Ydim; i++)
            for (int j = 1; j <= Xdim; j++)
                DIRECT_A2D_ELEM(*this,i - 1, j - 1) = m[i][j];
    }

    void killAdaptationForNumericalRecipes2D(T** m) const
    {
        free_Tmatrix(m, 1, YSIZE(*this), 1, XSIZE(*this));
    }

    void killAdaptationForNumericalRecipes22D(T** m) const
        {}

    T* adaptForNumericalRecipes1D() const
    {
        return MULTIDIM_ARRAY(*this) - 1;
    }

    void killAdaptationForNumericalRecipes1D(T* m) const
        {}

    void centerOfMass(Matrix1D< double >& center, void * mask=NULL, size_t n = 0)
    {
        center.initZeros(3);
        double mass = 0;
        MultidimArray< int >* imask = (MultidimArray< int >*) mask;

        FOR_ALL_ELEMENTS_IN_ARRAY3D(*this)
        {
            if ((imask == NULL || NZYX_ELEM(*imask, n, k, i, j)) &&
                A3D_ELEM(*this, k, i, j) > 0)
            {
                XX(center) += j * NZYX_ELEM(*this, n, k, i, j);
                YY(center) += i * NZYX_ELEM(*this, n, k, i, j);
                ZZ(center) += k * NZYX_ELEM(*this, n, k, i, j);

                mass += NZYX_ELEM(*this, n, k, i, j);
            }
        }

        if (mass != 0)
            center /= mass;
    }

    void sort(MultidimArray<T> &result) const;

    void indexSort(MultidimArray< int > &indx) const;

    void cumlativeDensityFunction(MultidimArray<double> &cdf)
    {
        MultidimArray<int> indx;
        indexSort(indx);
        cdf.resizeNoCopy(*this);
        double iSize=1.0/XSIZE(indx);
        FOR_ALL_ELEMENTS_IN_ARRAY1D(indx)
        {
            int ii=A1D_ELEM(indx,i)-1;
            A1D_ELEM(cdf,ii)=(i+1)*iSize;
        }
    }

    void threshold(const String& type,
                   T a,
                   T b=0,
                   MultidimArray<int> * mask = NULL )
    {
        int mode;

        if (type == "abs_above")
            mode = 1;
        else if (type == "abs_below")
            mode = 2;
        else if (type == "above")
            mode = 3;
        else if (type == "below")
            mode = 4;
        else if (type == "range")
            mode = 5;
        else if (type == "soft")
            mode = 6;
        else
            REPORT_ERROR(ERR_VALUE_INCORRECT,
                         formatString("Threshold: mode not supported (%s)", type.c_str() ));

        T* ptr=NULL;
        size_t n;
        T ma=-a;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        {
            if (mask == NULL || DIRECT_MULTIDIM_ELEM(*mask,n) > 0 )
            {
                switch (mode)
                {
                case 1:
                    if (ABS(*ptr) > a)
                        *ptr = SGN(*ptr) * b;
                    break;
                case 2:
                    if (ABS(*ptr) < a)
                        *ptr = SGN(*ptr) * b;
                    break;
                case 3:
                    if (*ptr > a)
                        *ptr = b;
                    break;
                case 4:
                    if (*ptr < a)
                        *ptr = b;
                    break;
                case 5:
                    if (*ptr < a)
                        *ptr = a;
                    else if (*ptr > b)
                        *ptr = b;
                    break;
                case 6:
                    if (*ptr<ma)
                        *ptr+=a;
                    else if (*ptr>a)
                        *ptr-=a;
                    else
                        *ptr=0;
                }
            }
        }
    }

    size_t countThreshold(const String& type,
                          T a,
                          T b,
                          MultidimArray<int> * mask = NULL )
    {
        int mode;

        if (type == "abs_above")
            mode = 1;
        else if (type == "abs_below")
            mode = 2;
        else if (type == "above")
            mode = 3;
        else if (type == "below")
            mode = 4;
        else if (type == "range")
            mode = 5;
        else
            REPORT_ERROR(ERR_VALUE_INCORRECT,
                         formatString("CountThreshold: mode not supported (%s)", type.c_str()));

        size_t ret = 0;

        T* ptr=NULL;
        size_t n;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        if (mask == NULL || DIRECT_MULTIDIM_ELEM(*mask,n) > 0 )
        {
            switch (mode)
            {
            case 1:
                if (ABS(*ptr) > a)
                    ret++;
                break;
            case 2:
                if (ABS(*ptr) < a)
                    ret++;
                break;
            case 3:
                if (*ptr > a)
                    ret++;
                break;
            case 4:
                if (*ptr < a)
                    ret++;
                break;
            case 5:
                if (*ptr >= a && *ptr <= b)
                    ret++;
                break;
            }
        }
        return ret;
    }

    void substitute(T oldv,
                    T newv,
                    double accuracy = XMIPP_EQUAL_ACCURACY,
                    MultidimArray<int> * mask = NULL )
    {
        T* ptr=NULL;
        size_t n;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        if (mask == NULL || DIRECT_MULTIDIM_ELEM(*mask,n) > 0 )
            if (ABS(*ptr - oldv) <= accuracy)
                *ptr = newv;
    }

    void randomSubstitute(T oldv,
                          T avgv,
                          T sigv,
                          double accuracy = XMIPP_EQUAL_ACCURACY,
                          MultidimArray<int> * mask = NULL );

    void binarize(double val = 0,
                  double accuracy = XMIPP_EQUAL_ACCURACY,
                  MultidimArray<int> * mask = NULL )
    {
        T* ptr=NULL;
        size_t n;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        if (mask == NULL || DIRECT_MULTIDIM_ELEM(*mask,n) > 0 )
        {
            if (*ptr <= val + accuracy)
                *ptr = 0;
            else
                *ptr = 1;
        }
    }

    void binarizeRange(double valMin = 0, double valMax = 255,
                       MultidimArray<int> * mask = NULL )
    {
        T* ptr=NULL;
        size_t n;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        if (mask == NULL || DIRECT_MULTIDIM_ELEM(*mask,n) > 0 )
        {
            if ( (*ptr < valMax) &&  (*ptr > valMin) )
                *ptr = 1;
            else
                *ptr = 0;
        }
    }

    void selfROUND()
    {
        T* ptr=NULL;
        size_t n;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        *ptr = ROUND(*ptr);
    }

    void selfCEIL()
    {
        T* ptr=NULL;
        size_t n;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        *ptr = CEIL(*ptr);
    }

    void selfFLOOR()
    {
        T* ptr=NULL;
        size_t n;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        *ptr = FLOOR(*ptr);
    }

    void selfABS()
    {
        T* ptr=NULL;
        size_t n;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        *ptr = ABS(*ptr);
    }

    void selfNormalizeInterval(double minPerc=0.25, double maxPerc=0.75, int Npix=1000);

    friend void MultidimArrayMax(const MultidimArray<T>& v1, const MultidimArray<T>& v2,
                                 MultidimArray<T>& result)
    {
        if (!v1.sameShape(v2))
            REPORT_ERROR(ERR_MULTIDIM_SIZE, "MAX: arrays of different shape");

        result.resizeNoCopy(v1);
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(result)
        DIRECT_MULTIDIM_ELEM(result,n) = XMIPP_MAX(
                                             DIRECT_MULTIDIM_ELEM(v1,n),
                                             DIRECT_MULTIDIM_ELEM(v2,n));
    }

    friend void MultidimArrayMIN(const MultidimArray<T>& v1, const MultidimArray<T>& v2,
                                 MultidimArray<T>& result)
    {
        if (!v1.sameShape(v2))
            REPORT_ERROR(ERR_MULTIDIM_SIZE, "MIN: arrays of different shape");

        result.resizeNoCopy(v1);
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(result)
        DIRECT_MULTIDIM_ELEM(result,n) = XMIPP_MIN(
                                             DIRECT_MULTIDIM_ELEM(v1,n),
                                             DIRECT_MULTIDIM_ELEM(v2,n));
    }

    void selfSQRT()
    {
        T* ptr=NULL;
        size_t n;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        *ptr = static_cast< T >(sqrt(static_cast< double >(*ptr)));
    }

    double sum() const
    {
        double sum = 0;
        T* ptr=NULL;
        size_t n;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        sum += *ptr;
        return sum;
    }

    double sum2() const
    {
        double sum = 0;

        // Unroll the loop
        const size_t unroll=4;
        size_t nmax=(NZYXSIZE(*this)/unroll)*unroll;
        T* ptr = MULTIDIM_ARRAY(*this);
        for (size_t n=0; n<nmax; n+=unroll, ptr+=unroll)
        {
            sum+= (*ptr)*(*ptr);
            T* ptr1=ptr+1;
            sum+= (*ptr1)*(*ptr1);
            T* ptr2=ptr+2;
            sum+= (*ptr2)*(*ptr2);
            T* ptr3=ptr+3;
            sum+= (*ptr3)*(*ptr3);
        }
        // Do the remaining elements
        for (size_t n=nmax; n<NZYXSIZE(*this); ++n, ++ptr)
            sum+=(*ptr)*(*ptr);
        return sum;
    }

    void selfLog10()
    {
        T* ptr=NULL;
        size_t n;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        *ptr = static_cast< T >(log10(static_cast< double >(*ptr)));
    }

    void selfLog()
    {
        T* ptr=NULL;
        size_t n;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(*this,n,ptr)
        *ptr = static_cast< T >(log(static_cast< double >(*ptr)));
    }

    void selfReverseX()
    {
        size_t xsize=XSIZE(*this);
        size_t halfSizeX = (xsize-2)/2;
        size_t xsize_1=xsize-1;
        for (size_t k = 0; k < ZSIZE(*this); k++)
            for (size_t i = 0; i < YSIZE(*this); i++)
                for (size_t j = 0; j <=halfSizeX; j++)
                {
                    T aux;
                    T& d1=DIRECT_ZYX_ELEM(*this, k, i, j);
                    T& d2=DIRECT_ZYX_ELEM(*this, k, i, xsize_1 - j);
                    SWAP(d1,d2,aux);
                }
        STARTINGX(*this) = -FINISHINGX(*this);
    }

    void selfReverseY()
    {
        size_t ysize=YSIZE(*this);
        size_t halfSizeY = (ysize-2)/2;
        size_t ysize_1=ysize-1;
        for (size_t k = 0; k < ZSIZE(*this); k++)
            for (size_t i = 0; i <= halfSizeY; i++)
                for (size_t j = 0; j < XSIZE(*this); j++)
                {
                    T aux;
                    T& d1=DIRECT_ZYX_ELEM(*this, k, i, j);
                    T& d2=DIRECT_ZYX_ELEM(*this, k, ysize_1 - i, j);
                    SWAP(d1,d2,aux);
                }
        STARTINGY(*this) = -FINISHINGY(*this);
    }

    void selfReverseZ()
    {
        size_t zsize=ZSIZE(*this);
        size_t halfSizeZ = (zsize-2)/2;
        size_t zsize_1=zsize-1;
        for (size_t k = 0; k <= halfSizeZ; k++)
            for (size_t i = 0; i <YSIZE(*this); i++)
                for (size_t j = 0; j < XSIZE(*this); j++)
                {
                    T aux;
                    T& d1=DIRECT_ZYX_ELEM(*this, k, i, j);
                    T& d2=DIRECT_ZYX_ELEM(*this, zsize_1- k, i, j);
                    SWAP(d1,d2,aux);
                }
        STARTINGZ(*this) = -FINISHINGZ(*this);
    }

    void profile(int x0, int y0, int xF, int yF, int N,
                 MultidimArray< double >& profile) const
    {
        checkDimension(2);
        profile.initZeros(N);
        double tx_step = (double)(xF - x0) / (N - 1);
        double ty_step = (double)(yF - y0) / (N - 1);
        double tx = x0, ty = y0;

        for (int i = 0; i < N; i++)
        {
            profile(i) = interpolatedElement2D(tx, ty);
            tx += tx_step;
            ty += ty_step;
        }
    }

    void showWithGnuPlot(const String& xlabel, const String& title);

    void edit();

    void write(const FileName& fn) const;


    MultidimArray<T>& operator=(const MultidimArray<T>& op1)
    {
        if (&op1 != this)
        {
            if (data == NULL || !sameShape(op1))
                resizeNoCopy(op1);
            memcpy(data,op1.data,MULTIDIM_SIZE(op1)*sizeof(T));
        }
        return *this;
    }

    MultidimArray<T>& operator=(MultidimArray<T>&& other) noexcept
    {
        if (&other != this)
        {
            coreDeallocate();
            coreInit();
            swap(other);
        }
        return *this;
    }

    MultidimArray<T>& operator=(const Matrix2D<T>& op1);

    MultidimArray<T> operator-() const
    {
        MultidimArray<T> tmp(*this);
        T* ptr;
        size_t n;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(tmp,n,ptr)
        *ptr = -(*ptr);
        return tmp;
    }

    friend std::istream& operator>>(std::istream& in, MultidimArray<T>& v)
    {
        T* ptr;
        size_t n;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(v,n,ptr)
        in >> *ptr;
        return in;
    }

    void copy(Matrix2D<T>& op1) const;

    bool equal(const MultidimArray<T>& op,
               double accuracy = XMIPP_EQUAL_ACCURACY) const
    {
        if (!sameShape(op) || data==NULL || op.data == NULL)
            return false;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(*this)
        {
            if (fabs(DIRECT_MULTIDIM_ELEM(*this,n) -
                     DIRECT_MULTIDIM_ELEM(op,n)) > accuracy)
                return false;
        }
        return true;
    }
};




template<typename T1, typename T2>
void typeCast(const MultidimArray<T1>& v1,  MultidimArray<T2>& v2)
{
    if (NZYXSIZE(v1) == 0)
    {
        v2.clear();
        return;
    }

    v2.resizeNoCopy(v1);
    T1* ptr1 = MULTIDIM_ARRAY(v1);

    // Unroll the loop
    const size_t unroll=4;
    size_t nmax=(NZYXSIZE(v1)/unroll)*unroll;
    for (size_t n=0; n<nmax; n+=unroll, ptr1+=unroll)
    {
        DIRECT_MULTIDIM_ELEM(v2,n)   = static_cast< T2 >(*ptr1);
        DIRECT_MULTIDIM_ELEM(v2,n+1) = static_cast< T2 >(*(ptr1+1));
        DIRECT_MULTIDIM_ELEM(v2,n+2) = static_cast< T2 >(*(ptr1+2));
        DIRECT_MULTIDIM_ELEM(v2,n+3) = static_cast< T2 >(*(ptr1+3));
    }
    // Do the remaining elements
    for (size_t n=nmax; n<NZYXSIZE(v1); ++n, ++ptr1)
        DIRECT_MULTIDIM_ELEM(v2,n)   = static_cast< T2 >(*ptr1);
}

template<typename T1>
void typeCastComplex(const MultidimArray<T1>& v1,  MultidimArray<std::complex<double> >& v2)
{
    if (NZYXSIZE(v1) == 0)
    {
        v2.clear();
        return;
    }

    v2.resizeNoCopy(v1);
    T1* ptr1 = MULTIDIM_ARRAY(v1);
    double * ptr2 = (double*) MULTIDIM_ARRAY(v2);

    // Unroll the loop
    const size_t unroll=4;
    size_t nmax=(NZYXSIZE(v1)/unroll)*unroll;
    for (size_t n=0; n<nmax; n+=unroll, ptr1+=unroll)
    {
        *(ptr2++) = static_cast<double>(*ptr1);
        *(ptr2++) = 0;
        *(ptr2++) = static_cast< double >(*(ptr1+1));
        *(ptr2++) = 0;
        *(ptr2++) = static_cast< double >(*(ptr1+2));
        *(ptr2++) = 0;
        *(ptr2++) = static_cast< double >(*(ptr1+3));
        *(ptr2++) = 0;
    }
    // Do the remaining elements
    for (size_t n=nmax; n<NZYXSIZE(v1); ++n, ++ptr1)
    {
        *(ptr2++) = static_cast< double >(*ptr1);
        *(ptr2++) = 0;
    }
}

template<typename T1>
void typeCast(const MultidimArray<T1>& v1,  MultidimArray<T1>& v2)
{
    v2=v1;
}

template<typename T>
void typeCast(const MultidimArray<T>& v1, Matrix1D<T> &v2)
{
    v2.resizeNoCopy(XSIZE(v1));
    memcpy(&VEC_ELEM(v2,0),&DIRECT_A1D_ELEM(v1,0),MULTIDIM_SIZE(v1)*sizeof(T));
    v2.row=false;
}


template<typename T>
bool operator==(const MultidimArray<T>& op1, const MultidimArray<T>& op2)
{
    return op1.equal(op2);
}

template<typename T>
bool operator!=(const MultidimArray<T>& op1, const MultidimArray<T>& op2)
{
    return !(op1==op2);
}

template<typename T>
void cutToCommonSize(MultidimArray<T>& V1, MultidimArray<T>& V2)
{
    int z0 = XMIPP_MAX(STARTINGZ(V1), STARTINGZ(V2));
    int zF = XMIPP_MIN(FINISHINGZ(V1), FINISHINGZ(V2));
    int y0 = XMIPP_MAX(STARTINGY(V1), STARTINGY(V2));
    int yF = XMIPP_MIN(FINISHINGY(V1), FINISHINGY(V2));
    int x0 = XMIPP_MAX(STARTINGX(V1), STARTINGX(V2));
    int xF = XMIPP_MIN(FINISHINGX(V1), FINISHINGX(V2));

    V1.window(z0, y0, x0, zF, yF, xF);
    V2.window(z0, y0, x0, zF, yF, xF);
}


void sincos(const MultidimArray<double> &x, MultidimArray<double> &s, MultidimArray<double> &c);

void planeFit(const MultidimArray<double> &z, const MultidimArray<double> &x, const MultidimArray<double> &y,
              double &p0, double &p1, double &p2);

/*
   mod    Modulus after division.
    mod(x,y) is x - n.*y where n = floor(x./y) if y ~= 0.  If y is not an
    integer and the quotient x./y is within roundoff error of an integer,
    then n is that integer.  The inputs x and y must be real arrays of the
    same size, or real scalars.

    By convention:
       MOD(x, m, 0) is m = x.
       MOD(x,m,x) is m = 0.
       MOD(x,m,y), for x~=y and y~=0, m has the same sign as y.
 */
template <typename T>
void mod(const MultidimArray<T> &x, MultidimArray<T> &m, double y)
{
    m.resizeNoCopy(x);
    double *ptr=NULL;
    double *ptrm=MULTIDIM_ARRAY(m);
    size_t n;
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(x,n,ptr)
    *(ptrm++) = (*ptr) - std::floor((*ptr)/(y))*(y);
}

template <typename T>
std::ostream& operator<< (std::ostream& ostrm, const MultidimArray<T>& v)
{
    if (v.xdim == 0)
        ostrm << "NULL Array\n";
    else
        ostrm << std::endl;

    double max_val = ABS(DIRECT_A3D_ELEM(v , 0, 0, 0));

    T* ptr;
    size_t n;
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(v,n,ptr)
    max_val = XMIPP_MAX(max_val, ABS(*ptr));

    int prec = bestPrecision(max_val, 10);

    if (YSIZE(v)==1 && ZSIZE(v)==1)
    {
        for (int j = STARTINGX(v); j <= FINISHINGX(v); j++)
            ostrm << floatToString((double) A3D_ELEM(v, 0, 0, j), 10, prec)
            << std::endl;
    }
    else
    {
        for (size_t l = 0; l < NSIZE(v); l++)
        {
            if (NSIZE(v)>1)
                ostrm << "Image No. " << l << std::endl;
            for (int k = STARTINGZ(v); k <= FINISHINGZ(v); k++)
            {
                if (ZSIZE(v)>1)
                    ostrm << "Slice No. " << k << std::endl;
                for (int i = STARTINGY(v); i <= FINISHINGY(v); i++)
                {
                    for (int j = STARTINGX(v); j <= FINISHINGX(v); j++)
                    {
                        ostrm << floatToString((double) A3D_ELEM(v, k, i, j), 10, prec) << ' ';
                    }
                    ostrm << std::endl;
                }
            }
        }
    }

    return ostrm;
}

template<typename T>
void window2D(const MultidimArray<T> &Ibig, MultidimArray<T> &Ismall,
        size_t y0, size_t x0, size_t yF, size_t xF);

template <typename T>
double correlationIndex(const MultidimArray< T >& x,
                        const MultidimArray< T >& y,
                        const MultidimArray< int >* mask = NULL,
                        MultidimArray< double >* Contributions = NULL)
{
    SPEED_UP_tempsInt;

    double retval = 0, aux;
    double mean_x, mean_y;
    double stddev_x, stddev_y;

    long N = 0;

    if (mask == NULL)
    {
        x.computeAvgStdev(mean_x, stddev_x);
        y.computeAvgStdev(mean_y, stddev_y);
    }
    else
    {
        x.computeAvgStdev_within_binary_mask(*mask, mean_x,stddev_x);
        y.computeAvgStdev_within_binary_mask(*mask, mean_y,stddev_y);
    }
    if (ABS(stddev_x)<XMIPP_EQUAL_ACCURACY ||
        ABS(stddev_y)<XMIPP_EQUAL_ACCURACY)
        return 0;

    // If contributions are desired. Please, be careful interpreting individual
    // contributions to the covariance! One pixel value afect others.
    if (Contributions != NULL)
    {
        FOR_ALL_ELEMENTS_IN_COMMON_IN_ARRAY3D(x, y)
        {
            if (mask != NULL)
                if (!A3D_ELEM(*mask,k, i, j))
                    continue;

            aux = (A3D_ELEM(x, k, i, j) - mean_x) * (A3D_ELEM(y, k, i, j) -
                    mean_y);
            A3D_ELEM(*Contributions, k, i, j) = aux;
            retval += aux;
            ++N;
        }

        FOR_ALL_ELEMENTS_IN_ARRAY3D(*Contributions)
        A3D_ELEM(*Contributions, k, i, j) /= ((stddev_x * stddev_y) * N);
    }
    else
    {
        if (mask==NULL && x.sameShape(y))
        {
            FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(x)
                retval += DIRECT_MULTIDIM_ELEM(x, n)*DIRECT_MULTIDIM_ELEM(y, n);
            N=MULTIDIM_SIZE(x);
            retval-=N*mean_x*mean_y;
        }
        else
        {
            FOR_ALL_ELEMENTS_IN_COMMON_IN_ARRAY3D(x, y)
            {
                if (mask != NULL)
                    if (!A3D_ELEM(*mask,k, i, j))
                        continue;

                retval += (A3D_ELEM(x, k, i, j) - mean_x) *
                          (A3D_ELEM(y, k, i, j) - mean_y);
                ++N;
            }
        }
    }

    if (N != 0)
        return retval / ((stddev_x * stddev_y) * N);
    else
        return 0;
}

// Specializations cases for complex numbers
template<>
std::ostream& operator<<(std::ostream& ostrm, const MultidimArray< std::complex<double> >& v);
template<>
void MultidimArray< std::complex< double > >::computeDoubleMinMax(double& minval, double& maxval) const;
template<>
void MultidimArray< std::complex< double > >::computeDoubleMinMaxRange(double& minval, double& maxval, size_t pos, size_t size) const;
template<>
void MultidimArray< std::complex< double > >::rangeAdjust(std::complex< double > minF, std::complex< double > maxF);
template<>
double MultidimArray< std::complex< double > >::computeAvg() const;
template<>
void MultidimArray< std::complex< double > >::maxIndex(size_t &lmax, int& kmax, int& imax, int& jmax) const;
template<>
bool operator==(const MultidimArray< std::complex< double > >& op1,
                const MultidimArray< std::complex< double > >& op2);
template<>
double MultidimArray<double>::interpolatedElement2D(double x, double y, double outside_value) const;
template<>
void MultidimArray< std::complex< double > >::getReal(MultidimArray<double> & realImg) const;
template<>
void MultidimArray< std::complex< double > >::getImag(MultidimArray<double> & imagImg) const;

#endif