grids.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'
 ***************************************************************************/
/* ------------------------------------------------------------------------- */
/* GRIDS                                                                     */
/* ------------------------------------------------------------------------- */

#ifndef _CORE_GRIDS_HH
#define _CORE_GRIDS_HH

#include <vector>

#include <core/xmipp_image.h>
#include <core/geometry.h>
#include <core/args.h>
#include "core/matrix2d.h"

/* Forward declarations ---------------------------------------------------- */
template <class T>
class GridVolumeT;
template <class T>
GridVolumeT<T> operator -(T f, const GridVolumeT<T> &GV);
template <class T>
GridVolumeT<T> operator /(T f, const GridVolumeT<T> &GV);
template <class T>
std::ostream& operator << (std::ostream &o, const GridVolumeT<T> &GV);
class Basis;

/*****************************************************************************/
/* Simple Grids                                                              */
/*****************************************************************************/
class SimpleGrid
{
    /* Structure --------------------------------------------------------------- */
public:
    Matrix2D<double> basis;                // These are the 3 unitary vectors
    // which define the grid. First
    // column is vector X, second
    // column is vector Y, and the third
    // is Z
    Matrix2D<double> inv_basis;            // Inverse matrix of the basis
    // It is used to save computation
    // time
    Matrix1D<double> lowest;
    Matrix1D<double> highest;
    double relative_size;
    Matrix1D<double> origin;
    double R2;

    /* Prototypes --------------------------------------------------------------- */
public:
    SimpleGrid();

    friend std::ostream& operator <<(std::ostream& o, const SimpleGrid &grid);

    void assign(const SimpleGrid &SG);

    void set_X(const Matrix1D<double> &v)
    {
        basis.setCol(0, v);
    }

    void set_Y(const Matrix1D<double> &v)
    {
        basis.setCol(1, v);
    }

    void set_Z(const Matrix1D<double> &v)
    {
        basis.setCol(2, v);
    }

    void get_X(Matrix1D<double> &v) const
    {
        basis.getCol(0, v);
    }

    void get_Y(Matrix1D<double> &v) const
    {
        basis.getCol(1, v);
    }

    void get_Z(Matrix1D<double> &v) const
    {
        basis.getCol(2, v);
    }

    void getSize(int &Zdim, int &Ydim, int &Xdim) const
    {
        Zdim = (int)(ZZ(highest) - ZZ(lowest)) + 1;
        Ydim = (int)(YY(highest) - YY(lowest)) + 1;
        Xdim = (int)(XX(highest) - XX(lowest)) + 1;
    }

    int get_number_of_samples() const;

    void set_interest_radius(double _R)
    {
        if (_R == -1)
            R2 = -1;
        else
            R2 = _R * _R;
    }

    double get_interest_radius() const
    {
        if (R2 == -1)
            return R2;
        else
            return sqrt(R2);
    }

    void set_relative_size(double size)
    {
        relative_size = size;
    }

    double get_relative_size() const
    {
        return relative_size;
    }

    void prepare_grid();

    void grid2universe(const Matrix1D<double> &gv, Matrix1D<double> &uv) const
    {
        SPEED_UP_temps012;
        uv.resize(3);
        M3x3_BY_V3x1(uv, basis, gv);
        V3_BY_CT(uv, uv, relative_size);
        V3_PLUS_V3(uv, uv, origin);
    }

    void universe2grid(const Matrix1D<double> &uv, Matrix1D<double> &gv) const
    {
        SPEED_UP_temps012;
        gv.resize(3);
        V3_MINUS_V3(gv, uv, origin);
        V3_BY_CT(gv, gv, 1 / relative_size);
        M3x3_BY_V3x1(gv, inv_basis, gv);
    }

    bool is_interesting(const Matrix1D<double> &uv) const
    {
        if (R2 == -1)
            return true;
        else
            return XX(uv)*XX(uv) + YY(uv)*YY(uv) + ZZ(uv)*ZZ(uv) < R2;
    }

    void Gproject_to_plane(const Matrix1D<double> &gr,
                           double rot, double tilt, double psi, Matrix1D<double> &result) const
    {
        grid2universe(gr, result);
        Uproject_to_plane(result, rot, tilt, psi, result);
    }

    void Gproject_to_plane(const Matrix1D<double> &gr,
                           const Matrix2D<double> &euler, Matrix1D<double> &result) const
    {
        grid2universe(gr, result);
        Uproject_to_plane(result, euler, result);
    }

    void Gdir_project_to_plane(const Matrix1D<double> &gr,
                               double rot, double tilt, double psi, Matrix1D<double> &result) const
    {
        grid2universe(gr, result);
        V3_MINUS_V3(result, result, origin);
        Uproject_to_plane(result, rot, tilt, psi, result);
    }

    void Gdir_project_to_plane(const Matrix1D<double> &gr,
                               const Matrix2D<double> &euler, Matrix1D<double> &result) const
    {
        grid2universe(gr, result);
        V3_MINUS_V3(result, result, origin);
        Uproject_to_plane(result, euler, result);
    }
};

/*****************************************************************************/
/* Complex Grids                                                             */
/*****************************************************************************/
class Grid
{
    /* Structure --------------------------------------------------------------- */
    std::vector<SimpleGrid>   LG;                 // List of grids
public:
    /* Protoypes --------------------------------------------------------------- */
    void add_grid(const SimpleGrid &SG)
    {
        LG.push_back(SG);
    }

    const SimpleGrid & operator()(size_t n) const
    {
        if (n>LG.size())
            REPORT_ERROR(ERR_VALUE_INCORRECT, "The Grid hasn't got so many Simple Grids");
        return LG[n];
    }

    SimpleGrid& operator()(size_t n)
    {
        if (n>LG.size())
            REPORT_ERROR(ERR_VALUE_INCORRECT, "The Grid hasn't got so many Simple Grids");
        return LG[n];
    }

    void get_SimpleGrid(int n, SimpleGrid& G) const
    {
        G = (*this)(n);
    }

    size_t GridsNo() const
    {
        return LG.size();
    }

    friend std::ostream& operator << (std::ostream& o, const Grid &grid)
    {
        o << "Complex Grid -------------------------------------\n";
        for (size_t i = 0; i < grid.GridsNo(); i++)
            o << grid(i);
        return o;
    }

    void assign(const Grid &G)
    {
        *this = G;
    }

    void clear()
    {
        LG.clear();
    }

    void voxel_corners(Matrix1D<double> &Gcorner1, Matrix1D<double> &Gcorner2,
                       const Matrix2D<double> *V = nullptr) const;
};

/*****************************************************************************/
/* Some useful Grids                                                         */
/*****************************************************************************/
#define  CC 0
#define FCC 1
#define BCC 2

SimpleGrid Create_CC_grid(double relative_size,
                          const Matrix1D<double> &corner1,
                          const Matrix1D<double> &corner2,
                          const Matrix1D<double> &origin);

Grid Create_CC_grid(double relative_size,
                    const Matrix1D<double> &corner1,
                    const Matrix1D<double> &corner2);

Grid Create_CC_grid(double relative_size,
                    int Zdim, int Ydim, int Xdim);

Grid Create_BCC_grid(double relative_size,
                     const Matrix1D<double> &corner1, const Matrix1D<double> &corner2);

Grid Create_FCC_grid(double relative_size,
                     const Matrix1D<double> &corner1, const Matrix1D<double> &corner2);

SimpleGrid Create_grid_within_sphere(double relative_size,
                                     const Matrix1D<double> &origin,
                                     const Matrix1D<double> &X, const Matrix1D<double> &Y,
                                     const Matrix1D<double> &Z, double R2);

Grid Create_CC_grid(double relative_size, double R);

Grid Create_BCC_grid(double relative_size, double R);

Grid Create_FCC_grid(double relative_size, double R);

/*****************************************************************************/
/* Grid Volumes                                                              */
/*****************************************************************************/
template <class T>
class GridVolumeT
{
    // Structure ---------------------------------------------------------------
    std::vector<Image<T> * > LV;               // List of volumes
    Grid                 G;                 // Grid associated to this volume

public:
    // Protoypes ---------------------------------------------------------------
    GridVolumeT()
    {
        LV.clear();
        (Grid) G;
    }

    GridVolumeT(const GridVolumeT &_RV)
    {
        *this = _RV;
    }

    GridVolumeT(const Grid &_grid)
    {
        adapt_to_grid(_grid);
    }

    GridVolumeT& operator = (const GridVolumeT& RV)
    {
        if (this != &RV)
        {
            clear();
            G = RV.G;
            for (size_t i = 0; i < RV.VolumesNo(); i++)
            {
                auto  *V = new Image<T>;
                *V = RV(i);
                LV.push_back(V);
            }
        }
        return *this;
    }

    void assign(const GridVolumeT& RV)
    {
        *this = RV;
    }

    ~GridVolumeT()
    {
        clear();
    }

    void adapt_to_grid(const Grid &_grid)
    {
        // Clear old list of volumes
        for (size_t i = 0; i < VolumesNo(); i++)
            if (LV[i]!=nullptr)
                delete LV[i];
        LV.clear();

        // Keep the incoming Grid at the same time the old grid is forgotten
        G = _grid;

        // Generate a volume for each subgrid
        int Zdim;
        int Ydim;
        int Xdim;
        Image<T> *                 Vol_aux;
        for (size_t i = 0; i < G.GridsNo(); i++)
        {
            SimpleGrid & grid = G(i);
            grid.getSize(Zdim, Ydim, Xdim);
            Vol_aux = new Image<T>;
            (*Vol_aux)().resize(Zdim, Ydim, Xdim);  // Using this function
            // after empty creation the volume
            // is zero-valued.
            STARTINGX((*Vol_aux)()) = (int) XX(grid.lowest); // This values are already
            STARTINGY((*Vol_aux)()) = (int) YY(grid.lowest); // integer although they
            STARTINGZ((*Vol_aux)()) = (int) ZZ(grid.lowest); // are stored as float
            LV.push_back(Vol_aux);
        }
    }

    void resize(const Matrix1D<double> &corner1,
                const Matrix1D<double> &corner2)
    {
        Image<T> *         Vol_aux;
        std::vector<Image<T> * > LV_aux;

        for (size_t n = 0; n < G.GridsNo(); n++)
        {
            SimpleGrid &grid = G(n);

            // Resize grid
            grid.universe2grid(corner1, grid.lowest);
            grid.lowest.selfFLOOR();
            grid.universe2grid(corner2, grid.highest);
            grid.highest.selfCEIL();

            // Resize auxiliary volume
            int Zdim;
            int Ydim;
            int Xdim;
            grid.getSize(Zdim, Ydim, Xdim);
            Vol_aux = new Image<T>;
            (*Vol_aux)().resize(Zdim, Ydim, Xdim);
            STARTINGX((*Vol_aux)()) = (int) XX(grid.lowest); // This values are already
            STARTINGY((*Vol_aux)()) = (int) YY(grid.lowest); // integer although they
            STARTINGZ((*Vol_aux)()) = (int) ZZ(grid.lowest); // are stored as float

            // Copy values in common
            Image<T> * origin = LV[n];
            SPEED_UP_tempsInt;
            FOR_ALL_ELEMENTS_IN_COMMON_IN_ARRAY3D(VOLMATRIX(*Vol_aux), VOLMATRIX(*origin))
            {
                VOLVOXEL(*Vol_aux, k, i, j) = VOLVOXEL(*origin, k, i, j);
            }

            // Extract old volume and push new one
            delete LV[n];
            LV_aux.push_back(Vol_aux);
        }
        LV = LV_aux;
    }

    template <class T1>
    void resize(const GridVolumeT<T1> &GV)
    {
        clear();
        for (size_t n = 0; n < GV.VolumesNo(); n++)
        {
            SimpleGrid grid;
            grid = GV.grid(n);
            G.add_grid(grid);

            Image<T> *Vol_aux;
            Vol_aux = new Image<T>;
            (*Vol_aux)().resize(GV(n)());
            LV.push_back(Vol_aux);
        }
    }

    void initZeros()
    {
        for (size_t i = 0; i < VolumesNo(); i++)
            (*this)(i)().initZeros();
    }

    void clear()
    {
        for (size_t i = 0; i < VolumesNo(); i++)
            delete LV[i];
        LV.clear();
        G.clear();
    }

    Image<T> & operator()(size_t n)
    {
        if (n>LV.size())
            REPORT_ERROR(ERR_VALUE_INCORRECT, "The Grid Volume hasn't got so many Simple Volumes");
        return *(LV[n]);
    }

    void get_volume(size_t n, Image<T> &V)
    {
        V = (*this)(n);
    }

    const Image<T> & operator()(size_t n) const
    {
        if (n>LV.size())
            REPORT_ERROR(ERR_VALUE_INCORRECT, "The Grid Volume hasn't got so many Simple Volumes");
        return *(LV[n]);
    }

    const SimpleGrid & grid(size_t n) const
    {
        return G(n);
    }

    void get_SimpleGrid(size_t n, SimpleGrid &G)
    {
        G = grid(n);
    }

    const Grid & grid() const
    {
        return G;
    }

    void get_Grid(Grid &G)
    {
        G = grid();
    }

    size_t VolumesNo() const
    {
        return LV.size();
    }

#define GRIDVOLUME_BY_SCALAR(op) \
    GridVolumeT<T> result; \
    result.G = G; \
    result.LV.reserve(VolumesNo()); \
    for (size_t i=0; i<VolumesNo(); i++) \
        array_by_scalar((*this)(i)(),f,result(i)(),op); \
    return result;

    GridVolumeT<T> operator + (T f) const
    {
        GRIDVOLUME_BY_SCALAR('+');
    }

    GridVolumeT<T> operator - (T f) const
    {
        GRIDVOLUME_BY_SCALAR('-');
    }

    GridVolumeT<T> operator *(T f) const
    {
        GRIDVOLUME_BY_SCALAR('*');
    }

    GridVolumeT<T> operator / (T f) const
    {
        GRIDVOLUME_BY_SCALAR('/');
    }

    GridVolumeT<T> friend operator + (T f, const GridVolumeT<T> &GV)
    {
        return GV + f;
    }

    GridVolumeT<T> friend operator *(T f, const GridVolumeT<T> &GV)
    {
        return GV*f;
    }

#define GRIDVOL_BY_GRIDVOL(op) \
    GridVolumeT<T> result; \
    Image<T> * Vol_aux; \
    \
    if (VolumesNo()!=GV.VolumesNo()) \
        REPORT_ERROR(ERR_GRID_SIZE,(std::string)"GridVolume::"+op+": Different number of subvolumes");\
    \
    result.G = G;\
    result.LV.reserve(VolumesNo());\
    \
    for (size_t i=0; i<VolumesNo(); i++) { \
        try { \
            Vol_aux = new Image<T>; \
            arrayByArray((*this)(i)(),GV(i)(),(*Vol_aux)(),op); \
            result.LV.push_back(Vol_aux); \
        } catch (XmippError XE) {\
            std::cout << XE.what(); \
            REPORT_ERROR(ERR_GRID_SIZE,(std::string)"GridVolume::"+op+": Different shape of volume " +\
                         integerToString(i)); \
        } \
    } \
    \
    return result;

    GridVolumeT<T> operator + (const GridVolumeT<T> &GV)
    {
        GRIDVOL_BY_GRIDVOL('+');
    }

    GridVolumeT<T> operator - (const GridVolumeT<T> &GV)
    {
        GRIDVOL_BY_GRIDVOL('-');
    }

    GridVolumeT<T> operator *(const GridVolumeT<T> &GV)
    {
        GRIDVOL_BY_GRIDVOL('*');
    }

    GridVolumeT<T> operator / (const GridVolumeT<T> &GV)
    {
        GRIDVOL_BY_GRIDVOL('/');
    }

#define GRIDVOL_BY_GRIDVOLASSIG(op) \
    if (VolumesNo()!=GV.VolumesNo()) \
        REPORT_ERROR(ERR_GRID_SIZE,(std::string)"GridVolume::"+op+"=: Different number of subvolumes");\
    \
    for (size_t i=0; i<VolumesNo(); i++) { \
        try { \
            arrayByArray((*this)(i)(),GV(i)(),(*this)(i)(),op); \
        } catch (XmippError XE) {\
            std::cout << XE.what(); \
            REPORT_ERROR(ERR_GRID_SIZE,(std::string)"GridVolume::"+op+"=: Different shape of volume " +\
                         integerToString(i)); \
        } \
    }

    void operator += (const GridVolumeT<T> &GV)
    {
        GRIDVOL_BY_GRIDVOLASSIG('+');
    }

    void operator -= (const GridVolumeT<T> &GV)
    {
        GRIDVOL_BY_GRIDVOLASSIG('-');
    }

    void operator *= (const GridVolumeT<T> &GV)
    {
        GRIDVOL_BY_GRIDVOLASSIG('*');
    }

    void operator /= (const GridVolumeT<T> &GV)
    {
        GRIDVOL_BY_GRIDVOLASSIG('/');
    }

    void write(const FileName &fn) const
    {
        Image<T>    V;
        float temp_float;
        size_t floatsize;
        const std::type_info &typeinfoT = typeid(T); // We need to know what kind
        // of variable is T
        const std::type_info &typeinfoD = typeid(double);
        const std::type_info &typeinfoI = typeid(int);

        floatsize = (size_t) sizeof(float);

        if (VolumesNo() == 0)
            return;

        // Create the writing volume ............................................
        size_t Zdim = 0;
        size_t Ydim = 0;
        size_t Xdim = 0;
        for (size_t v = 0; v < VolumesNo(); v++)
        {
            const Image<T> & this_vol = (*this)(v);
            Zdim += ZSIZE(this_vol());
            Ydim = XMIPP_MAX(Ydim, YSIZE(this_vol()));
            Xdim = XMIPP_MAX(Xdim, XSIZE(this_vol()));
        }

        // Check if there is enough space for the control slice
        if (Xdim*Ydim < 25)
            Ydim = (int) CEIL(25.0f / Xdim);

        // A slice is added for control information for each subvolume
        VOLMATRIX(V).initZeros(Zdim + VolumesNo(), Ydim, Xdim);

        // Write Grid volume ....................................................
#define PACK_DOUBLE(v) \
{jj=pos%Xdim; ii=pos/Xdim; pos++; VOLVOXEL(V,sli,ii,jj)=(T)(v);}
#define PACK_INT(v) \
    {jj=pos%Xdim; ii=pos/Xdim; pos++; \
        temp_float = (float) (v); \
        memcpy( &(VOLVOXEL(V,sli,ii,jj)) , &temp_float, floatsize); \
    }

        int sli = 0;
        for (size_t v = 0; v < VolumesNo(); v++)
        {
            int pos;                // Position inside the control slice
            int ii;                 // Position inside the control slice
            int jj;                 // Position inside the control slice
            size_t k;               // Auxiliar counters
            size_t i;               // Auxiliar counters
            size_t j;               // Auxiliar counters

            // Choose grid and volume
            const SimpleGrid & this_grid = grid(v);
            const Image<T> & this_vol = (*this)(v);

            // Store Grid data ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
            pos = 0;

            if (typeinfoT == typeinfoD)
            {
                for (i = 0; i < 3; i++)
                    for (j = 0; j < 3; j++)
                        PACK_DOUBLE((this_grid.basis)(i, j));
                for (i = 0; i < 3; i++)
                    PACK_DOUBLE((this_grid.lowest)(i));
                for (i = 0; i < 3; i++)
                    PACK_DOUBLE((this_grid.highest)(i));
                PACK_DOUBLE(this_grid.relative_size);
                for (i = 0; i < 3; i++)
                    PACK_DOUBLE((this_grid.origin)(i));
                PACK_DOUBLE(this_grid.R2);

                // Store volume control ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
                PACK_DOUBLE(ZSIZE(this_vol()));
                PACK_DOUBLE(YSIZE(this_vol()));
                PACK_DOUBLE(XSIZE(this_vol()));
                PACK_DOUBLE(STARTINGZ(this_vol()));
                PACK_DOUBLE(STARTINGY(this_vol()));
                PACK_DOUBLE(STARTINGX(this_vol()));
            }
            else if (typeinfoT == typeinfoI)
            {
                // We use a trick to save the grid information in the volume
                // If the following if is true the trick can not be used
                if (sizeof(float) != sizeof(int))
                    REPORT_ERROR(ERR_TYPE_INCORRECT,
                                 "GridVolume is integer and (sizeof(float)!= sizeof(int)");

                for (i = 0; i < 3; i++)
                    for (j = 0; j < 3; j++)
                        PACK_INT((this_grid.basis)(i, j));
                for (i = 0; i < 3; i++)
                    PACK_INT((this_grid.lowest)(i));
                for (i = 0; i < 3; i++)
                    PACK_INT((this_grid.highest)(i));
                PACK_INT(this_grid.relative_size);
                for (i = 0; i < 3; i++)
                    PACK_INT((this_grid.origin)(i));
                PACK_INT(this_grid.R2);

                // Store volume control ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
                PACK_INT(ZSIZE(this_vol()));
                PACK_INT(YSIZE(this_vol()));
                PACK_INT(XSIZE(this_vol()));
                PACK_INT(STARTINGZ(this_vol()));
                PACK_INT(STARTINGY(this_vol()));
                PACK_INT(STARTINGX(this_vol()));
            }
            else
                REPORT_ERROR(ERR_TYPE_INCORRECT, "GridVolume must be double or int\n");

            sli++;

            // Write the whole volume ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
            for (k = 0; k < ZSIZE(VOLMATRIX(this_vol)); k++)
            {
                for (i = 0; i < YSIZE(VOLMATRIX(this_vol)); i++)
                    for (j = 0; j < XSIZE(VOLMATRIX(this_vol)); j++)
                        DIRECT_VOLVOXEL(V, sli, i, j) = DIRECT_VOLVOXEL(this_vol, k, i, j);
                sli++;
            }
        }
#undef PACK_DOUBLE
#undef PACK_INT

        // Effectively write the volume .........................................
        V.write(fn);
    }

    void read(const FileName &fn, const std::string &basisName)
    {
        Image<T>       V;
        Image<T>       * sV;
        SimpleGrid     sG;
        size_t         sli = 0;

        float temp_float;
        size_t floatsize;
        const std::type_info &typeinfoT = typeid(T); // We need to know what kind
        // of variable is T
        const std::type_info &typeinfoD = typeid(double);
        const std::type_info &typeinfoI = typeid(int);

        floatsize = (size_t) sizeof(float);
        // We use a trick to save the grid information in the volume
        // If the following if is true the trick can not be used
        if ((typeid(T) == typeid(int)) && (sizeof(float) != sizeof(int)))
            REPORT_ERROR(ERR_TYPE_INCORRECT,"Error: GridVolume is integer and (sizeof(float)!= sizeof(int)");

        // Allocate memory ......................................................
        sG.basis.resize(3, 3);
        sG.lowest.resize(3);
        sG.highest.resize(3);
        sG.origin.resize(3);

        // Read Reconstructing volume from file .................................
        V.read(fn);
        if (basisName=="voxels")
        {
            V().setXmippOrigin();
            sV=new Image<double>;
            *sV=V;
            LV.push_back(sV);
            G=Create_CC_grid(1.0,ZSIZE(V()),YSIZE(V()),XSIZE(V()));
        }
        else
        {
#define UNPACK_DOUBLE(v,cast) \
            {jj=pos%VOLMATRIX(V).xdim; ii=pos/VOLMATRIX(V).xdim; pos++; \
            (v)=(cast)VOLVOXEL(V,sli,ii,jj);}
#define UNPACK_INT(v,cast) \
            {jj=pos%VOLMATRIX(V).xdim; ii=pos/VOLMATRIX(V).xdim; pos++; \
            memcpy( &temp_float, &(VOLVOXEL(V,sli,ii,jj)),floatsize);\
            (v)=(cast)temp_float;}

            while (sli < ZSIZE(V()))
            {
                int pos;                // Position inside the control slice
                int ii;                 // Position inside the control slice
                int jj;                 // Position inside the control slice
                size_t k;               // Auxiliary counters
                size_t i;               // Auxiliary counters
                size_t j;               // Auxiliary counters
                size_t Zdim;
                size_t Ydim;
                size_t Xdim;
                int Zinit;
                int Yinit;
                int Xinit;

                // Read Grid data ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
                pos = 0;
                if (typeinfoT == typeinfoD)
                {
                    for (i = 0; i < 3; i++)
                        for (j = 0; j < 3; j++)
                            UNPACK_DOUBLE((sG.basis)(i, j), double);
                    for (i = 0; i < 3; i++)
                        UNPACK_DOUBLE((sG.lowest)(i), int);
                    for (i = 0; i < 3; i++)
                        UNPACK_DOUBLE((sG.highest)(i), int);
                    UNPACK_DOUBLE(sG.relative_size, double);
                    for (i = 0; i < 3; i++)
                        UNPACK_DOUBLE((sG.origin)(i), double);
                    UNPACK_DOUBLE(sG.R2, double);
                }
                else if (typeinfoT == typeinfoI)
                {
                    // We use a trick to save the grid information in the volume
                    // If the following if is true the trick can not be used
                    if (sizeof(float) != sizeof(int))
                        REPORT_ERROR(ERR_TYPE_INCORRECT,
                                     "GridVolume is integer and (sizeof(float)!= sizeof(int)");

                    for (i = 0; i < 3; i++)
                        for (j = 0; j < 3; j++)
                            UNPACK_INT((sG.basis)(i, j), double);
                    for (i = 0; i < 3; i++)
                        UNPACK_INT((sG.lowest)(i), int);
                    for (i = 0; i < 3; i++)
                        UNPACK_INT((sG.highest)(i), int);
                    UNPACK_INT(sG.relative_size, double);
                    for (i = 0; i < 3; i++)
                        UNPACK_INT((sG.origin)(i), double);
                    UNPACK_INT(sG.R2, double);
                }
                sG.inv_basis = sG.basis.inv();

                // Store Grid in the list of the grid volume
                G.add_grid(sG);

                // Read Volume Control Information ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
                if (typeinfoT == typeinfoD)
                {
                    UNPACK_DOUBLE(Zdim, int);
                    UNPACK_DOUBLE(Ydim, int);
                    UNPACK_DOUBLE(Xdim, int);
                    UNPACK_DOUBLE(Zinit, int);
                    UNPACK_DOUBLE(Yinit, int);
                    UNPACK_DOUBLE(Xinit, int);
                }
                else if (typeinfoT == typeinfoI)
                {
                    UNPACK_INT(Zdim, int);
                    UNPACK_INT(Ydim, int);
                    UNPACK_INT(Xdim, int);
                    UNPACK_INT(Zinit, int);
                    UNPACK_INT(Yinit, int);
                    UNPACK_INT(Xinit, int);
                }

                // Set volume size and origin
                sV = new Image<T>;
                VOLMATRIX(*sV).initZeros(Zdim, Ydim, Xdim);
                STARTINGZ(VOLMATRIX(*sV)) = Zinit;
                STARTINGY(VOLMATRIX(*sV)) = Yinit;
                STARTINGX(VOLMATRIX(*sV)) = Xinit;
#ifdef DEBUG

                std::cout << "The read grid is \n" << sG;
                std::cout << "Volume dimensions: " << Zdim << " x " << Ydim << " x "
                << Xdim << std::endl;
                std::cout << "Volume init: " << Zinit << " x " << Yinit << " x "
                << Xinit << std::endl;
#endif

                sli++;

                // Read volume ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
                for (k = 0; k < ZSIZE(VOLMATRIX(*sV)); k++)
                {
                    for (i = 0; i < YSIZE(VOLMATRIX(*sV)); i++)
                        for (j = 0; j < XSIZE(VOLMATRIX(*sV)); j++)
                        {
#ifdef DEBUG
                            std::cout << "Reading from file position (" << sli << "," << i
                            << "," << j << ") to subvolume position ("
                            << k << "," << i << "," << j << ")\n";
#endif

                            DIRECT_VOLVOXEL(*sV, k, i, j) = DIRECT_VOLVOXEL(V, sli, i, j);
                        }
                    sli++;
                }

                // Store volume in the list
                LV.push_back(sV);
            }
#undef UNPACK_DOUBLE
#undef UNPACK_INT
        }
    }
#undef DEBUG
};

typedef GridVolumeT<double> GridVolume;

// Show a grid volume ------------------------------------------------------
template <class T>
std::ostream& operator << (std::ostream &o, const GridVolumeT<T> &GV)
{
    o << "Grid Volume -----------\n";
    o << GV.G;
    o << "Number of volumes= " << GV.VolumesNo() << std::endl;
    for (int i = 0; i < GV.VolumesNo(); i++)
    {
        o << "Volume " << i << "------------" << std::endl;
        o << GV(i)();
    }
    return o;
}
#endif