multidim_array_base.h

Go to the documentation of this file.

/***************************************************************************
 *
 * Authors:    David Strelak (davidstrelak@gmail.com)
 *
 * 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 XMIPPCORE_CORE_MULTIDIM_ARRAY_BASE_H_
#define XMIPPCORE_CORE_MULTIDIM_ARRAY_BASE_H_

#include <stddef.h>
#include <iostream>
#include "xmipp_array_dim.h"
#include "xmipp_array_coord.h"

template<typename T>
class Matrix1D;

#define STARTINGX(v) ((v).xinit)

#define FINISHINGX(v) ((v).xinit + (int)(v).xdim - 1)

#define STARTINGY(v) ((v).yinit)

#define FINISHINGY(v) ((v).yinit + (int)(v).ydim - 1)

#define STARTINGZ(v) ((v).zinit)

#define FINISHINGZ(v) ((v).zinit + (int)(v).zdim - 1)

#define INSIDEX(v, x) ((x) >= STARTINGX(v) && (x) <= FINISHINGX(v))

#define INSIDEY(v, y) ((y) >= STARTINGY(v) && (y) <= FINISHINGY(v))

#define INSIDEZ(v, z) ((z) >= STARTINGZ(v) && (z) <= FINISHINGZ(v))

#define INSIDEXY(v, x, y) (INSIDEX(v, x) && INSIDEY(v, y))

#define INSIDEXYZ(v, x, y, z) (INSIDEX(v, x) && INSIDEY(v, y) && INSIDEZ(v,z))

#define XSIZE(v) ((v).xdim)

#define YSIZE(v) ((v).ydim)

#define ZSIZE(v) ((v).zdim)

#define NSIZE(v) ((v).ndim)

#define YXSIZE(v) ((v).yxdim)

#define ZYXSIZE(v) ((v).zyxdim)

#define MULTIDIM_SIZE(v) ((v).nzyxdim)

#define NZYXSIZE(v) ((v).nzyxdim)

#ifndef MULTIDIM_ARRAY
#define MULTIDIM_ARRAY(v) ((v).data)
#endif

#define DIRECT_NZYX_ELEM(v, l, k, i, j) ((v).data[(l)*ZYXSIZE(v)+(k)*YXSIZE(v)+((i)*XSIZE(v))+(j)])

#define DIRECT_ZYX_ELEM(v, k, i, j) ((v).data[(k)*YXSIZE(v)+((i)*XSIZE(v))+(j)])

#define DIRECT_N_YX_ELEM(v, l, i, j) ((v).data[(l)*ZYXSIZE(v)             +((i)*XSIZE(v))+(j)])

#define DIRECT_N__X_ELEM(v, l, j) ((v).data[(l)*ZYXSIZE(v)+(j)])

#define NZYX_ELEM(v, l, k, i, j)  \
    DIRECT_NZYX_ELEM((v), (l), (k) - STARTINGZ(v), (i) - STARTINGY(v), (j) - STARTINGX(v))

#define DIRECT_MULTIDIM_ELEM(v,n) ((v).data[(n)])

#define FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(v) \
    for (size_t n=0; n<NZYXSIZE(v); ++n)

#define FOR_ALL_DIRECT_NZYX_ELEMENTS_IN_MULTIDIMARRAY(V) \
    for (size_t l=0; l<NSIZE(V); ++l) \
        for (size_t k=0; k<ZSIZE(V); ++k) \
            for (size_t i=0; i<YSIZE(V); ++i)      \
                for (size_t j=0; j<XSIZE(V); ++j)

#define FOR_ALL_NZYX_ELEMENTS_IN_MULTIDIMARRAY(V) \
    for (size_t l=0; l<NSIZE(V); ++l) \
        for (int k=STARTINGZ(V); k<=FINISHINGZ(V); ++k) \
            for (int i=STARTINGY(V); i<=FINISHINGY(V); ++i)     \
                for (int j=STARTINGX(V); j<=FINISHINGX(V); ++j)

#define FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(v,n,ptr) \
    for ((n)=0, (ptr)=(v).data; (n)<NZYXSIZE(v); ++(n), ++(ptr))

#define DIRECT_A3D_ELEM(v,k,i,j) ((v).data[YXSIZE(v)*(k)+((i)*XSIZE(v))+(j)])

#define dAkij(V, k, i, j) DIRECT_A3D_ELEM(V, k, i, j)

#define A3D_ELEM(V, k, i, j) \
    DIRECT_A3D_ELEM((V),(k) - STARTINGZ(V), (i) - STARTINGY(V), (j) - STARTINGX(V))

#define FOR_ALL_ELEMENTS_IN_ARRAY3D(V) \
    for (int k=STARTINGZ(V); k<=FINISHINGZ(V); ++k) \
        for (int i=STARTINGY(V); i<=FINISHINGY(V); ++i) \
            for (int j=STARTINGX(V); j<=FINISHINGX(V); ++j)

#define FOR_ALL_ELEMENTS_IN_COMMON_IN_ARRAY3D(V1, V2) \
    ispduptmp0 = XMIPP_MAX(STARTINGZ(V1), STARTINGZ(V2)); \
    ispduptmp1 = XMIPP_MIN(FINISHINGZ(V1),FINISHINGZ(V2)); \
    ispduptmp2 = XMIPP_MAX(STARTINGY(V1), STARTINGY(V2)); \
    ispduptmp3 = XMIPP_MIN(FINISHINGY(V1),FINISHINGY(V2)); \
    ispduptmp4 = XMIPP_MAX(STARTINGX(V1), STARTINGX(V2)); \
    ispduptmp5 = XMIPP_MIN(FINISHINGX(V1),FINISHINGX(V2)); \
    for (int k=ispduptmp0; k<=ispduptmp1; ++k) \
        for (int i=ispduptmp2; i<=ispduptmp3; ++i) \
            for (int j=ispduptmp4; j<=ispduptmp5; ++j)

#define FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(V) \
    for (size_t k=0; k<ZSIZE(V); ++k) \
        for (size_t i=0; i<YSIZE(V); ++i) \
            for (size_t j=0; j<XSIZE(V); ++j)

#define DIRECT_A2D_ELEM(v,i,j) ((v).data[(i)*(v).xdim+(j)])

#define dAij(M, i, j) DIRECT_A2D_ELEM(M, i, j)

#define A2D_ELEM(v, i, j) \
    DIRECT_A2D_ELEM(v, (i) - STARTINGY(v), (j) - STARTINGX(v))

#define SAME_SHAPE2D(v1, v2) \
    (XSIZE(v1) == XSIZE(v2) && \
     YSIZE(v1) == YSIZE(v2) && \
     STARTINGX(v1) == STARTINGX(v2) && \
     STARTINGY(v1) == STARTINGY(v2))

#define SAME_SHAPE3D(v1, v2) \
    (XSIZE(v1) == XSIZE(v2) && \
     YSIZE(v1) == YSIZE(v2) && \
     ZSIZE(v1) == ZSIZE(v2) && \
     STARTINGX(v1) == STARTINGX(v2) && \
     STARTINGY(v1) == STARTINGY(v2) && \
     STARTINGZ(v1) == STARTINGZ(v2))


#define FOR_ALL_ELEMENTS_IN_ARRAY2D(m) \
    for (int i=STARTINGY(m); i<=FINISHINGY(m); ++i) \
        for (int j=STARTINGX(m); j<=FINISHINGX(m); ++j)

#define FOR_ALL_ELEMENTS_IN_COMMON_IN_ARRAY2D(m1, m2) \
    ispduptmp2 = XMIPP_MAX(STARTINGY(m1), STARTINGY(m2)); \
    ispduptmp3 = XMIPP_MIN(FINISHINGY(m1), FINISHINGY(m2)); \
    ispduptmp4 = XMIPP_MAX(STARTINGX(m1), STARTINGX(m2)); \
    ispduptmp5 = XMIPP_MIN(FINISHINGX(m1), FINISHINGX(m2)); \
    for (int i=ispduptmp2; i<=ispduptmp3; ++i) \
        for (int j=ispduptmp4; j<=ispduptmp5; ++j)

#define FOR_ALL_ELEMENTS_IN_ARRAY3D_BETWEEN(corner1, corner2) \
    for (ZZ(r)=ZZ((corner1)); ZZ(r)<=ZZ((corner2)); ++ZZ(r)) \
        for (YY(r)=YY((corner1)); YY(r)<=YY((corner2)); ++YY(r)) \
            for (XX(r)=XX((corner1)); XX(r)<=XX((corner2)); ++XX(r))

#define FOR_ALL_ELEMENTS_IN_ARRAY2D_BETWEEN(corner1, corner2) \
    for (YY(r)=YY((corner1)); YY(r)<=YY((corner2)); ++YY(r)) \
        for (XX(r)=XX((corner1)); XX(r)<=XX((corner2)); ++XX(r))

#define FOR_ALL_ELEMENTS_IN_ARRAY1D_BETWEEN(corner1, corner2) \
    for (XX(r)=(int) XX((corner1)); XX(r)<=(int) XX((corner2)); ++XX(r))

#define FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(m) \
    for (size_t i=0; i<YSIZE(m); ++i) \
        for (size_t j=0; j<XSIZE(m); ++j)

#define DIRECT_A1D_ELEM(v, i) ((v).data[(i)])

#define dAi(v, i) DIRECT_A1D_ELEM(v, i)

#define A1D_ELEM(v, i) DIRECT_A1D_ELEM(v, (i) - ((v).xinit))

#define FOR_ALL_ELEMENTS_IN_ARRAY1D(v) \
    for (int i=STARTINGX(v); i<=FINISHINGX(v); ++i)

#define FOR_ALL_ELEMENTS_IN_COMMON_IN_ARRAY1D(v1, v2) \
    ispduptmp4 = XMIPP_MAX(STARTINGX(v1), STARTINGX(v2)); \
    ispduptmp5 = XMIPP_MIN(FINISHINGX(v1), FINISHINGX(v2)); \
    for (int i=ispduptmp4; i<=ispduptmp5; ++i)

#define FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY1D(v) \
    for (size_t i=0; i<v.xdim; ++i)

#define OUTSIDE(i,j) \
            ((j) < STARTINGX(*this) || (j) > FINISHINGX(*this) || \
             (i) < STARTINGY(*this) || (i) > FINISHINGY(*this))

#define OUTSIDE3D(k, i,j) \
            ((j) < STARTINGX(*this) || (j) > FINISHINGX(*this) || \
             (i) < STARTINGY(*this) || (i) > FINISHINGY(*this) || \
             (k) < STARTINGZ(*this) || (k) > FINISHINGZ(*this))


// Look up table lenght to be used in interpolation.
#define     LOOKUP_TABLE_LEN        6


class MultidimArrayBase
{
public:
    // Destroy data
    bool destroyData;

    // Number of images
    size_t ndim;

    // Number of elements in Z
    size_t zdim;

    // Number of elements in Y
    size_t ydim;

    // Number of elements in X
    size_t xdim;

    // Number of elements in YX
    size_t yxdim;

    // Number of elements in ZYX
    size_t zyxdim;

    // Number of elements in NZYX
    size_t nzyxdim;

    // Z init
    int zinit;

    // Y init
    int yinit;

    // X init
    int xinit;

    //Alloc memory or map to a file
    bool     mmapOn;
    //Mapped file handler
    FILE*      mFd;
    // Number of elements in NZYX in allocated memory
    size_t nzyxdimAlloc;
public:
    virtual ~MultidimArrayBase()
    {}

    // Virtual declarations to be used from MultidimArrayGeneric
    virtual void clear() = 0;
    virtual void selfReverseX() = 0;
    virtual void selfReverseY() = 0;
    virtual void selfReverseZ() = 0;
    virtual double computeAvg() const = 0;
    virtual void computeDoubleMinMaxRange(double& minval, double& maxval,size_t offset, size_t size) const = 0;
    virtual void maxIndex(size_t &lmax, int& kmax, int& imax, int& jmax) const = 0;
    virtual void coreAllocateReuse() = 0;
    virtual void coreDeallocate()= 0;

    /* return the value of the data pointer
     */
    virtual void * getArrayPointer() const = 0;



    void setNdim(int Ndim);

    void setZdim(int Zdim);

    void setYdim(int Ydim);

    void setXdim(int Xdim);

    void setDimensions(int Xdim, int Ydim, int Zdim, size_t Ndim);

    void setDimensions(ArrayDim &newDim);

    void getDimensions(size_t& Xdim, size_t& Ydim, size_t& Zdim, size_t &Ndim) const;
    void getDimensions(ArrayDim &idim) const;
    ArrayDim getDimensions() const;


    void getDimensions(int* size) const;

    size_t getSize() const;

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

    void resize(size_t Zdim, size_t Ydim, size_t Xdim)
    {
        resize(1, Zdim, Ydim, Xdim);
    }

    void resize(size_t Ydim, size_t Xdim)
    {
        resize(1, 1, Ydim, Xdim);
    }

    void resize(size_t Xdim)
    {
        resize(1, 1, 1, Xdim);
    }

    void resize(ArrayDim &adim, bool copy=true);

    void resizeNoCopy(size_t Ndim, size_t Zdim, size_t Ydim, size_t Xdim)
    {
        resize(Ndim, Zdim, Ydim, Xdim, false);
    }

    void resizeNoCopy(size_t Zdim, size_t Ydim, size_t Xdim)
    {
        resize(1, Zdim, Ydim, Xdim, false);
    }

    void resizeNoCopy(size_t Ydim, size_t Xdim)
    {
        resize(1, 1, Ydim, Xdim, false);
    }

    void resizeNoCopy(size_t Xdim)
    {
        resize(1, 1, 1, Xdim, false);
    }

    inline size_t rowNumber() const
    {
        return ydim;
    }

    inline size_t colNumber() const
    {
        return xdim;
    }

    void copyShape(const MultidimArrayBase &m);

    inline bool sameShape(const MultidimArrayBase &op) const
    {
        return (NSIZE(*this) == NSIZE(op) &&
                XSIZE(*this) == XSIZE(op) &&
                YSIZE(*this) == YSIZE(op) &&
                ZSIZE(*this) == ZSIZE(op) &&
                STARTINGX(*this) == STARTINGX(op) &&
                STARTINGY(*this) == STARTINGY(op) &&
                STARTINGZ(*this) == STARTINGZ(op));
    }

    void setXmippOrigin();

    void resetOrigin();

    void moveOriginTo(int k, int i, int j);

    void moveOriginTo(int i, int j);

    inline int startingZ() const
    {
        return zinit;
    }

    inline int finishingZ() const
    {
        return zinit + zdim - 1;
    }

    inline int startingY() const
    {
        return yinit;
    }

    inline int finishingY() const
    {
        return yinit + ydim - 1;
    }

    inline int startingX() const
    {
        return xinit;
    }

    inline int finishingX() const
    {
        return xinit + xdim - 1;
    }

    bool isCorner(const Matrix1D< double >& v) const;

    inline bool outside(int k, int i, int j) const
    {
        return (j < STARTINGX(*this) || j > FINISHINGX(*this) ||
                i < STARTINGY(*this) || i > FINISHINGY(*this) ||
                k < STARTINGZ(*this) || k > FINISHINGZ(*this));
    }

    inline bool outside(int i, int j) const
    {
        return (j < STARTINGX(*this) || j > FINISHINGX(*this) ||
                i < STARTINGY(*this) || i > FINISHINGY(*this));
    }

    inline bool outside(int i) const
    {
        return (i < STARTINGX(*this) || i > FINISHINGX(*this));
    }

    bool outside(const Matrix1D<double> &r) const;



    inline int getDim() const
    {
        if (NZYXSIZE(*this) < 1)
            return 0;
        if (NSIZE(*this) > 1)
            return 4;
        if (ZSIZE(*this) > 1)
            return 3;
        if (YSIZE(*this) > 1)
            return 2;
        return 1;
    }

    void setMmap(bool mmap)
    {
        coreDeallocate();
        mmapOn = mmap;
    }

    void maxIndex(ArrayCoord &pos) const
    {
        maxIndex(pos.n, pos.z, pos.y, pos.x);
    }

    void maxIndex(int& kmax, int& imax, int& jmax) const
    {
        size_t dum;
        maxIndex(dum, kmax, imax, jmax);
    }

    void maxIndex(int& imax, int& jmax) const
    {
        size_t dum;
        int idum;
        maxIndex(dum, idum, imax, jmax);
    }

    void maxIndex(int& jmax) const
    {
        size_t dum;
        int idum;
        maxIndex(dum, idum, idum, jmax);
    }

    void printShape(std::ostream& out = std::cout) const;
};

#endif /* XMIPPCORE_CORE_MULTIDIM_ARRAY_BASE_H_ */