librelion-dev-common 1.3+dfsg-2 / usr / include / relion-1.3 / src / image.h

/***************************************************************************
 *
 * Author: "Sjors H.W. Scheres"
 * MRC Laboratory of Molecular Biology
 *
 * 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.
 *
 * This complete copyright notice must be included in any revised version of the
 * source code. Additional authorship citations may be added, but existing
 * author citations must be preserved.
 ***************************************************************************/
/***************************************************************************
 *
 * Authors: Sjors H.W. Scheres (scheres@cnb.csic.es)
 *
 * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
 *
 * Part of this module has been developed by Lorenzo Zampighi and Nelson Tang
 * Dept. Physiology of the David Geffen School of Medicine
 * Univ. of California, Los Angeles.
 *
 * 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 IMAGE_H
#define IMAGE_H

#include <typeinfo>
#include <fcntl.h>
#include <sys/mman.h>
#include <fcntl.h>
#include <unistd.h>
#include "src/funcs.h"
#include "src/memory.h"
#include "src/filename.h"
#include "src/multidim_array.h"
#include "src/transformations.h"
#include "src/metadata_table.h"
#include "src/fftw.h"

/// @defgroup Images Images
//@{

/** Data type.
 * This class defines the datatype of the data inside this image.
 */
typedef enum
{
    Unknown_Type = 0,       // Undefined data type
    UChar = 1,              // Unsigned character or byte type
    SChar = 2,              // Signed character (for CCP4)
    UShort = 3,             // Unsigned integer (2-byte)
    Short = 4,              // Signed integer (2-byte)
    UInt = 5,               // Unsigned integer (4-byte)
    Int = 6,                // Signed integer (4-byte)
    Long = 7,               // Signed integer (4 or 8 byte, depending on system)
    Float = 8,              // Floating point (4-byte)
    Double = 9,             // Double precision floating point (8-byte)
    Bool = 10,              // Boolean (1-byte?)
    LastEntry = 15          // This must be the last entry
} DataType;

/** Write mode
 * This class defines the writing behavior.
 */
typedef enum
{
    WRITE_OVERWRITE, //forget about the old file and overwrite it
    WRITE_APPEND,    //append and object at the end of a stack, so far can not append stacks
    WRITE_REPLACE,   //replace a particular object by another
    WRITE_READONLY   //only can read the file
} WriteMode;


/** Open File struct
 * This struct is used to share the File handlers with Image Collection class
 */
struct fImageHandler
{
    FILE*     fimg;       // Image File handler
    FILE*     fhed;       // Image File header handler
    FileName  ext_name;   // Filename extension
    bool     exist;       // Shows if the file exists
};

/** Returns memory size of datatype
 */
unsigned long gettypesize(DataType type);

/** Convert datatype string to datatypr enun */
int datatypeString2Int(std::string s);

/** Swapping trigger.
 * Threshold file z size above which bytes are swapped.
 */
#define SWAPTRIG     65535

/** Template class for images.
 * The image class is the general image handling class.
 */
template<typename T>
class Image
{
public:
    MultidimArray<T>    data;        // The image data array
    MetaDataTable MDMainHeader;      // metadata for the file

private:
    FileName            filename;    // File name
    FILE*                fimg;       // Image File handler
    FILE*                fhed;       // Image File header handler
    bool                stayOpen;    // To maintain the image file open after read/write
    int                 dataflag;    // Flag to force reading of the data
    unsigned long       i;           // Current image number (may be > NSIZE)
    unsigned long       offset;      // Data offset
    int                 swap;        // Perform byte swapping upon reading
    long int            replaceNsize;// Stack size in the replace case
    bool                _exists;     // does target file exists?
    // equal 0 is not exists or not a stack
    bool                mmapOn;      // Mapping when loading from file
    int                 mFd;         // Handle the file in reading method and mmap
    size_t              mappedSize;  // Size of the mapped file

public:
    /** Empty constructor
     *
     * An empty image is created.
     *
     * @code
     * Image<double> I;
     * @endcode
     */
    Image()
    {
        mmapOn = false;
        clear();
        MDMainHeader.addObject();
    }

    /** Constructor with size
     *
     * A blank image (0.0 filled) is created with the given size. Pay attention
     * to the dimension order: Y and then X.
     *
     * @code
     * Image I(64,64);
     * @endcode
     */
    Image(long int Xdim, long int Ydim, long int Zdim=1, long int Ndim=1)
    {
        mmapOn = false;
        clear();
        data.resize(Ndim, Zdim, Ydim, Xdim);
        MDMainHeader.addObject();
    }

    /** Clear.
     * Initialize everything to 0
     */
    void clear()
    {
        if (mmapOn)
        {
            munmap(data.data-offset,mappedSize);
            close(mFd);
            data.data = NULL;
        }
        else
            data.clear();

        dataflag = -1;
        i = 0;
        filename = "";
        offset = 0;
        swap = 0;
        clearHeader();
        replaceNsize=0;
        mmapOn = false;
    }

    /** Clear the header of the image
     */
    void clearHeader()
    {
        MDMainHeader.clear();
    }

    /** Destructor.
     */
    ~Image()
    {
        clear();
    }


    /** Specific read functions for different file formats
      */
#include "src/rwSPIDER.h"
#include "src/rwMRC.h"
#include "src/rwIMAGIC.h"

    /** Is this file an image
     *
     *  Check whether a real-space image can be read
     *
     */
    bool isImage(const FileName &name)
    {
        return !read(name, false);
    }

    /** Rename the image
      */
    void rename (const FileName &name)
    {
        filename = name;
    }

    /** General read function
     * you can read a single image from a single image file
     * or a single image file from an stack, in the second case
     * the select slide may come in the image name or in the select_img parameter
     * file name takes precedence over select_img
     * If -1 is given the whole object is read
     *
     */
    int read(const FileName &name, bool readdata=true, long int select_img=-1, bool mapData = false, bool is_2D = false)
    {

        int err = 0;
        fImageHandler* hFile = openFile(name);
        err = _read(name, hFile, readdata, select_img, mapData, is_2D);
        closeFile(hFile);

        // Negative errors are bad
        return err;
    }

    /** Read function from a file that has already been opened
     *
     */
    int readFromOpenFile(const FileName &name, fImageHandler* hFile, long int select_img, bool is_2D = false)
    {
    	int err = 0;
    	err = _read(name, hFile, true, select_img, false, is_2D);
    	// Reposition file pointer for a next read
    	rewind(fimg);
    	return err;
    }

    /** General write function
     * select_img= which slice should I replace
     * overwrite = 0, append slice
     * overwrite = 1 overwrite slice
     */
    void write(FileName name="",
               long int select_img=-1,
               bool isStack=false,
               int mode=WRITE_OVERWRITE)
    {

        const FileName &fname = (name == "") ? filename : name;
        fImageHandler* hFile = openFile(fname, mode);
        _write(fname, hFile, select_img, isStack, mode);
        closeFile(hFile);

    }

    /** Cast a page of data from type dataType to type Tdest
     *    input pointer  char *
     */
    void castPage2T(char * page, T * ptrDest, DataType datatype, size_t pageSize )
    {
        switch (datatype)
        {
        case Unknown_Type:
            REPORT_ERROR("ERROR: datatype is Unknown_Type");
        case UChar:
            {
                if (typeid(T) == typeid(unsigned char))
                    memcpy(ptrDest, page, pageSize*sizeof(T));
                else
                {
                    unsigned char * ptr = (unsigned char *) page;
                    for(int i=0; i<pageSize; i++)
                        ptrDest[i]=(T) ptr[i];
                }
                break;
            }
        case SChar:
                {
                    if (typeid(T) == typeid(signed char))
                {
                    memcpy(ptrDest, page, pageSize*sizeof(T));
                    }
                    else
                    {
                        signed char * ptr = (signed char *) page;
                        for(int i=0; i<pageSize; i++)
                            ptrDest[i]=(T) ptr[i];
                    }
                break;
            }
        case UShort:
                {
                    if (typeid(T) == typeid(unsigned short))
                {
                    memcpy(ptrDest, page, pageSize*sizeof(T));
                    }
                    else
                    {
                        unsigned short * ptr = (unsigned short *) page;
                        for(int i=0; i<pageSize; i++)
                            ptrDest[i]=(T) ptr[i];
                    }
                break;
            }
        case Short:
                {
                    if (typeid(T) == typeid(short))
                {
                    memcpy(ptrDest, page, pageSize*sizeof(T));
                    }
                    else
                    {
                        short * ptr = (short *) page;
                        for(int i=0; i<pageSize; i++)
                            ptrDest[i]=(T) ptr[i];
                    }
                break;
            }
        case UInt:
                {
                    if (typeid(T) == typeid(unsigned int))
                {
                    memcpy(ptrDest, page, pageSize*sizeof(T));
                    }
                    else
                    {
                        unsigned int * ptr = (unsigned int *) page;
                        for(int i=0; i<pageSize; i++)
                            ptrDest[i]=(T) ptr[i];
                    }
                break;
            }
        case Int:
                {
                    if (typeid(T) == typeid(int))
                {
                    memcpy(ptrDest, page, pageSize*sizeof(T));
                    }
                    else
                    {
                        int * ptr = (int *) page;
                        for(int i=0; i<pageSize; i++)
                            ptrDest[i]=(T) ptr[i];
                    }
                break;
            }
        case Long:
                {
                    if (typeid(T) == typeid(long))
                {
                    memcpy(ptrDest, page, pageSize*sizeof(T));
                    }
                    else
                    {
                        long * ptr = (long *) page;
                        for(int i=0; i<pageSize; i++)
                            ptrDest[i]=(T) ptr[i];
                    }
                break;
            }
        case Float:
                {
                    if (typeid(T) == typeid(float))
                {
                    memcpy(ptrDest, page, pageSize*sizeof(T));
                    }
                    else
                    {
                        float * ptr = (float *) page;
                        for(int i=0; i<pageSize; i++)
                            ptrDest[i]=(T) ptr[i];
                    }
                break;
            }
        case Double:
                {
                    if (typeid(T) == typeid(double))
                {
                    memcpy(ptrDest, page, pageSize*sizeof(T));
                    }
                    else
                    {
                        double * ptr = (double *) page;
                        for(int i=0; i<pageSize; i++)
                            ptrDest[i]=(T) ptr[i];
                    }
                break;
            }
        default:
                {
                    std::cerr<<"Datatype= "<<datatype<<std::endl;
                    REPORT_ERROR(" ERROR: cannot cast datatype to T");
                    break;
                }
            }

    }

    /** Cast page from T to datatype
     *  input pointer char *
     */
    void castPage2Datatype(T * srcPtr, char * page, DataType datatype, size_t pageSize )
    {
        switch (datatype)
        {
        case Float:
            {
                if (typeid(T) == typeid(float))
                {
                    memcpy(page, srcPtr, pageSize*sizeof(T));
                }
                else
                {
                    float * ptr = (float *) page;
                    for(int i=0; i<pageSize; i++)
                        ptr[i] = (float)srcPtr[i];
                }
                break;
            }
        case Double:
                {
                    if (typeid(T) == typeid(double))
                {
                    memcpy(page, srcPtr, pageSize*sizeof(T));
                    }
                    else
                    {
                        double * ptr = (double *) page;
                        for(int i=0; i<pageSize; i++)
                            ptr[i] = (double)srcPtr[i];
                    }
                break;
            }
       case UShort:
                {
                    if (typeid(T) == typeid(unsigned short))
                {
                    memcpy(page, srcPtr, pageSize*sizeof(T));
                    }
                    else
                    {
                        unsigned short * ptr = (unsigned short *) page;
                        for(int i=0; i<pageSize; i++)
                            ptr[i] = (unsigned short)srcPtr[i];
                    }
                break;
            }
        case UChar:
                {
                    if (typeid(T) == typeid(unsigned char))
                {
                    memcpy(page, srcPtr, pageSize*sizeof(T));
                    }
                    else
                    {
                        unsigned char * ptr = (unsigned char *) page;
                        for(int i=0; i<pageSize; i++)
                            ptr[i] = (unsigned char)srcPtr[i];
                    }
                break;
            }
        default:
                {
                    std::cerr<<"outputDatatype= "<<datatype<<std::endl;
                    REPORT_ERROR(" ERROR: cannot cast T to outputDatatype");
                    break;
                }
            }
    }

    /** Check file Datatype is same as T type to use mmap.
     */
    bool checkMmapT(DataType datatype)
    {

        switch (datatype)
        {
        case Unknown_Type:
            REPORT_ERROR("ERROR: datatype is Unknown_Type");
        case UChar:
            {
                if (typeid(T) == typeid(unsigned char))
                    return 1;
                else
                    return 0;
                break;
            }
        case SChar:
            {
                if (typeid(T) == typeid(signed char))
                    return 1;
                else
                    return 0;
                break;
            }
        case UShort:
            {
                if (typeid(T) == typeid(unsigned short))
                    return 1;
                else
                    return 0;
                break;
            }
        case Short:
            {
                if (typeid(T) == typeid(short))
                    return 1;
                else
                    return 0;
                break;
            }
        case UInt:
            {
                if (typeid(T) == typeid(unsigned int))
                    return 1;
                else
                    return 0;
                break;
            }
        case Int:
            {
                if (typeid(T) == typeid(int))
                    return 1;
                else
                    return 0;
                break;
            }
        case Long:
            {
                if (typeid(T) == typeid(long))
                    return 1;
                else
                    return 0;
                break;
            }
        case Float:
            {
                if (typeid(T) == typeid(float))
                    return 1;
                else
                    return 0;
                break;
            }
        case Double:
            {
                if (typeid(T) == typeid(double))
                    return 1;
                else
                    return 0;
                break;
            }
        default:
            {
                std::cerr<<"Datatype= "<<datatype<<std::endl;
                REPORT_ERROR(" ERROR: cannot cast datatype to T");
                break;
            }
        }
        //               int * iTemp = (int*) map;
        //                ptrDest = reinterpret_cast<T*> (iTemp);
    }

    /** Write an entire page as datatype
     *
     * A page of datasize_n elements T is cast to datatype and written to fimg
     * The memory for the casted page is allocated and freed internally.
     */
    void writePageAsDatatype(FILE * fimg, DataType datatype, size_t datasize_n )
    {
        size_t datasize = datasize_n * gettypesize(datatype);
        char * fdata = (char *) askMemory(datasize);
        castPage2Datatype(MULTIDIM_ARRAY(data), fdata, datatype, datasize_n);
        fwrite( fdata, datasize, 1, fimg );
        freeMemory(fdata, datasize);
    }

    /** Swap an entire page
      * input pointer char *
      */
    void swapPage(char * page, size_t pageNrElements, DataType datatype)
    {
        unsigned long datatypesize = gettypesize(datatype);
#ifdef DEBUG

        std::cerr<<"DEBUG swapPage: Swapping image data with swap= "
        << swap<<" datatypesize= "<<datatypesize
        << " pageNrElements " << pageNrElements
        << " datatype " << datatype
        <<std::endl;
        ;
#endif

        // Swap bytes if required
        if ( swap == 1 )
        {
            for ( unsigned long i=0; i<pageNrElements; i+=datatypesize )
                swapbytes(page+i, datatypesize);
        }
        else if ( swap > 1 )
        {
            for ( unsigned long i=0; i<pageNrElements; i+=swap )
                swapbytes(page+i, swap);
        }
    }

    /** Read the raw data
      */
    int readData(FILE* fimg, long int select_img, DataType datatype, unsigned long pad)
    {
        //#define DEBUG
#ifdef DEBUG
        std::cerr<<"entering readdata"<<std::endl;
        std::cerr<<" readData flag= "<<dataflag<<std::endl;
#endif

        if ( dataflag < 1 )
            return 0;

        size_t myoffset, readsize, readsize_n, pagemax = 1073741824; //1Gb
        size_t datatypesize=gettypesize(datatype);
        size_t pagesize  =ZYXSIZE(data)*datatypesize;
        size_t haveread_n=0;

        //Multidimarray mmapOn is priority over image mmapOn
        if(data.mmapOn)
            mmapOn = false;

        // Flag to know that data is not going to be mapped although mmapOn is true
        if (mmapOn && !checkMmapT(datatype))
        {
            std::cout << "WARNING: Image Class. File datatype and image declaration not compatible with mmap. Loading into memory." <<std::endl;
            mmapOn = false;
            mFd = -1;
        }

        if (mmapOn)
        {
            if ( NSIZE(data) > 1 )
            {
                REPORT_ERROR("Image Class::ReadData: mmap with multiple \
                             images file not compatible. Try selecting a unique image.");
            }

            fclose(fimg);

            //if ( ( mFd = open(filename.c_str(), O_RDWR, S_IREAD | S_IWRITE) ) == -1 )
            if ( ( mFd = open(filename.c_str(), O_RDWR, S_IRUSR | S_IWUSR) ) == -1 )
                REPORT_ERROR("Image Class::ReadData: Error opening the image file.");

            char * map;
            mappedSize = pagesize+offset;

            if ( (map = (char*) mmap(0,mappedSize, PROT_READ | PROT_WRITE, MAP_SHARED, mFd, 0)) == (void*) -1 )
                REPORT_ERROR("Image Class::ReadData: mmap of image file failed.");
            data.data = reinterpret_cast<T*> (map+offset);
        }
        else
        {
            // Reset select to get the correct offset
            if ( select_img < 0 )
                select_img = 0;

            char* page = NULL;

            // Allocate memory for image data (Assume xdim, ydim, zdim and ndim are already set
            //if memory already allocated use it (no resize allowed)
            data.coreAllocateReuse();
            myoffset = offset + select_img*(pagesize + pad);
            //#define DEBUG

#ifdef DEBUG

            data.printShape();
            printf("DEBUG: Page size: %ld offset= %d \n", pagesize, offset);
            printf("DEBUG: Swap = %d  Pad = %ld  Offset = %ld\n", swap, pad, offset);
            printf("DEBUG: myoffset = %d select_img= %d \n", myoffset, select_img);
#endif

            if (pagesize > pagemax)
                page = (char *) askMemory(pagemax*sizeof(char));
            else
                page = (char *) askMemory(pagesize*sizeof(char));

            int error_fseek = fseek( fimg, myoffset, SEEK_SET );
            if (error_fseek != 0)
            	return -1;

            for ( size_t myn=0; myn<NSIZE(data); myn++ )
            {
                for (size_t myj=0; myj<pagesize; myj+=pagemax )//pagesize size of object
                {
                    // Read next page. Divide pages larger than pagemax
                    readsize = pagesize - myj;
                    if ( readsize > pagemax )
                        readsize = pagemax;
                    readsize_n = readsize/datatypesize;

                    //Read page from disc
                    size_t result = fread( page, readsize, 1, fimg );
                    if (result != 1)
                    	return -2;

                    //swap per page
                    if (swap)
                        swapPage(page, readsize, datatype);
                    // cast to T per page
                    castPage2T(page, MULTIDIM_ARRAY(data) + haveread_n, datatype, readsize_n);
                    haveread_n += readsize_n;
                }
                if ( pad > 0 )
                {
                    //fread( padpage, pad, 1, fimg);
                	error_fseek = fseek( fimg, pad, SEEK_CUR );
                    if (error_fseek != 0)
                    	return -1;
                }
            }
            //if ( pad > 0 )
            //    freeMemory(padpage, pad*sizeof(char));
            if ( page > 0 )
                freeMemory(page, pagesize*sizeof(char));

#ifdef DEBUG

            printf("DEBUG img_read_data: Finished reading and converting data\n");
#endif

        }
        return 0;
    }

    /** Data access
     *
     * This operator can be used to access the data multidimarray.
     * In this way we could resize an image just by
     * resizing its associated matrix or we could add two images by adding their
     * matrices.
     * @code
     * I().resize(128, 128);
     * I2() = I1() + I2();
     * @endcode
     */
    MultidimArray<T>& operator()()
    {
        return data;
    }
    const MultidimArray<T>& operator()() const
    {
        return data;
    }

    /** Pixel access
    *
    * This operator is used to access a pixel within a 2D image. This is a
    * logical access, so you could access to negative positions if the image
    * has been defined so (see the general explanation for the class).
    *
    * @code
    * std::cout << "Grey level of pixel (-3,-3) of the image = " << I(-3, -3)
    * << std::endl;
    *
    * I(-3, -3) = I(-3, -2);
    * @endcode
    */
    T& operator()(int i, int j) const
    {
        return A2D_ELEM(data, i, j);
    }
    /** Set pixel
     * (direct access) needed by swig
     */
    void setPixel(int i, int j, T v)
    {
        IMGPIXEL(*this,i,j)=v;
    }

    /** Get pixel
     * (direct acces) needed by swig
     */
    T getPixel(int i, int j) const
    {
        return IMGPIXEL(*this,i,j);
    }

    /** Voxel access
     *
     * This operator is used to access a voxel within a 3D image. This is a
     * logical access, so you could access to negative positions if the image
     * has been defined so (see the general explanation for the class).
     *
     * @code
     * std::cout << "Grey level of pixel (-3,-3, 1) of the volume = " << I(-3, -3, 1)
     * << std::endl;
     *
     * I(-3, -3, 1) = I(-3, -2, 0);
     * @endcode
     */
    T& operator()(int k, int i, int j) const
    {
        return A3D_ELEM(data, k, i, j);
    }

    /** Get file name
     *
     * @code
     * std::cout << "Image name = " << I.name() << std::endl;
     * @endcode
     */
    const FileName & name() const
    {
        return filename;
    }

    /** Get Image dimensions
     */
    void getDimensions(int &Xdim, int &Ydim, int &Zdim, long int &Ndim) const
    {
        Xdim = XSIZE(data);
        Ydim = YSIZE(data);
        Zdim = ZSIZE(data);
        Ndim = NSIZE(data);
    }

    long unsigned int getSize() const
    {
        return NZYXSIZE(data);
    }

    /* Is there label in the main header */
    bool mainContainsLabel(EMDLabel label) const
    {
        return MDMainHeader.containsLabel(label);
    }

    /** Data type
        *
        * @code
        * std::cout << "datatype= " << dataType() << std::endl;
        * @endcode
        */
    int dataType() const
    {
        int dummy;
        MDMainHeader.getValue(EMDL_IMAGE_DATATYPE, dummy);
        return dummy;
    }

    /** Sampling RateX
    *
    * @code
    * std::cout << "sampling= " << samplingRateX() << std::endl;
    * @endcode
    */
    double samplingRateX(const long int n = 0) const
    {
        double dummy = 1.;
        MDMainHeader.getValue(EMDL_IMAGE_SAMPLINGRATE_X, dummy);
        return dummy;
    }

    /** Sampling RateY
    *
    * @code
    * std::cout << "sampling= " << samplingRateY() << std::endl;
    * @endcode
    */
    double samplingRateY(const long int n = 0) const
    {
        double dummy = 1.;
        MDMainHeader.getValue(EMDL_IMAGE_SAMPLINGRATE_Y, dummy);
        return dummy;
    }

    /** Set file name
     */
    void setName(const FileName &_filename)
    {
        filename = _filename;
    }

    /* Set image statistics in the main header
     *
     */
    void setStatisticsInHeader()
    {
    	double avg,stddev,minval,maxval;
    	data.computeStats(avg, stddev, minval, maxval);
    	MDMainHeader.setValue(EMDL_IMAGE_STATS_AVG, avg);
    	MDMainHeader.setValue(EMDL_IMAGE_STATS_STDDEV, stddev);
    	MDMainHeader.setValue(EMDL_IMAGE_STATS_MIN, minval);
    	MDMainHeader.setValue(EMDL_IMAGE_STATS_MAX, maxval);
    }

    void setSamplingRateInHeader(double rate_x, double rate_y = -1., double rate_z = -1.)
    {
    	MDMainHeader.setValue(EMDL_IMAGE_SAMPLINGRATE_X, rate_x);
    	if (rate_y < 0.)
    		rate_y = rate_x;
    	MDMainHeader.setValue(EMDL_IMAGE_SAMPLINGRATE_Y, rate_y);
    	if (ZSIZE(data)>1)
    	{
    	   	if (rate_z < 0.)
    	    	rate_z = rate_x;
    	   	MDMainHeader.setValue(EMDL_IMAGE_SAMPLINGRATE_Z, rate_z);
    	}
    }

    /** Show image properties
      */
    friend std::ostream& operator<<(std::ostream& o, const Image<T>& I)
    {
        o << "Image type   : ";
            o << "Real-space image" << std::endl;

        o << "Reversed     : ";
        if (I.swap)
            o << "TRUE"  << std::endl;
        else
            o << "FALSE" << std::endl;

        o << "Data type    : ";
        switch (I.dataType())
        {
        case Unknown_Type:
            o << "Undefined data type";
            break;
        case UChar:
            o << "Unsigned character or byte type";
            break;
        case SChar:
            o << "Signed character (for CCP4)";
            break;
        case UShort:
            o << "Unsigned integer (2-byte)";
            break;
        case Short:
            o << "Signed integer (2-byte)";
            break;
        case UInt:
            o << "Unsigned integer (4-byte)";
            break;
        case Int:
            o << "Signed integer (4-byte)";
            break;
        case Long:
            o << "Signed integer (4 or 8 byte, depending on system)";
            break;
        case Float:
            o << "Floating point (4-byte)";
            break;
        case Double:
            o << "Double precision floating point (8-byte)";
            break;
        case Bool:
            o << "Boolean (1-byte?)";
            break;
        }
        o << std::endl;

        o << "dimensions   : " << NSIZE(I()) << " x " << ZSIZE(I()) << " x " << YSIZE(I()) << " x " << XSIZE(I());
        o << "  (noObjects x slices x rows x columns)" << std::endl;
        return o;
    }

    /** Sum this object with other file and keep in this object
      */
    void sumWithFile(const FileName &fn)
    {
        Image<T> aux;
        aux.read(fn);
        (*this)()+=aux();
    }

    /** Open file function
      * Open the image file and returns its file hander.
      */
    fImageHandler* openFile(const FileName &name, int mode = WRITE_READONLY)
    {
        fImageHandler* hFile = new fImageHandler;
        FileName fileName, headName = "";
        FileName ext_name = name.getFileFormat();

        long int dump;
        name.decompose(dump, fileName);
        // Subtract 1 to have numbering 0...N-1 instead of 1...N
        if (dump > 0)
        	dump--;

        fileName = fileName.removeFileFormat();

        size_t found = fileName.find_first_of("%");
        if (found!=std::string::npos)
          fileName = fileName.substr(0, found) ;

        hFile->exist = exists(fileName);

        std::string wmChar;

        switch (mode)
        {
        case WRITE_READONLY:
            if (!hFile->exist)
                REPORT_ERROR((std::string) "Cannot read file " + fileName + " It does not exist" );
            wmChar = "r";
            break;
        case WRITE_OVERWRITE:
            wmChar = "w";
            break;
        case WRITE_APPEND:
            if (_exists = exists(fileName))
                wmChar = "r+";
            else
                wmChar = "w+";
            break;
        case WRITE_REPLACE:
            wmChar = "r+";
            break;
        }


        if (ext_name.contains("img") || ext_name.contains("hed"))
        {
            fileName = fileName.withoutExtension();
            headName = fileName.addExtension("hed");
            fileName = fileName.addExtension("img");
        }

        // Open image file
        if ( ( hFile->fimg = fopen(fileName.c_str(), wmChar.c_str()) ) == NULL )
            REPORT_ERROR((std::string)"Image::openFile cannot open: " + name);

        if (headName != "")
        {
            if ( ( hFile->fhed = fopen(headName.c_str(), wmChar.c_str()) ) == NULL )
                REPORT_ERROR((std::string)"Image::openFile cannot open: " + headName);
        }
        else
            hFile->fhed = NULL;

        hFile->ext_name =ext_name;

        return hFile;
    }

    /** Close file function.
      * Close the image file according to its name and file handler.
      */
    void closeFile(fImageHandler* hFile = NULL)
    {
        FileName ext_name;
        FILE* fimg, *fhed;

        if (hFile != NULL)
        {
            ext_name = hFile->ext_name;
            fimg = hFile->fimg;
            fhed = hFile->fhed;
        }
        else
        {
            ext_name = filename.getFileFormat();
            fimg = this->fimg;
            fhed = this->fhed;
        }

        if (fclose(fimg) != 0 )
            REPORT_ERROR((std::string)"Can not close image file "+ filename);

        if (fhed != NULL &&  fclose(fhed) != 0 )
            REPORT_ERROR((std::string)"Can not close header file of "
                         + filename);

        delete hFile;
    }

private:
    int _read(const FileName &name, fImageHandler* hFile, bool readdata=true, long int select_img = -1,
              bool mapData = false, bool is_2D = false)
    {
        int err = 0;

        // Check whether to read the data or only the header
        dataflag = ( readdata ) ? 1 : -1;

        // Check whether to map the data or not
        mmapOn = mapData;

        FileName ext_name = hFile->ext_name;
        fimg = hFile->fimg;
        fhed = hFile->fhed;

        long int dump;
        name.decompose(dump, filename);
        // Subtract 1 to have numbering 0...N-1 instead of 1...N
        if (dump > 0)
        	dump--;
        filename = name;

        if (select_img == -1)
            select_img = dump;

#undef DEBUG
        //#define DEBUG
#ifdef DEBUG

        std::cerr << "READ\n" <<
        "name="<<name <<std::endl;
        std::cerr << "ext= "<<ext_name <<std::endl;
        std::cerr << " now reading: "<< filename <<" dataflag= "<<dataflag
        << " select_img "  << select_img << std::endl;
#endif
#undef DEBUG

        //Just clear the header before reading
        MDMainHeader.clear();
        MDMainHeader.addObject();

        if (ext_name.contains("spi") || ext_name.contains("xmp")  ||
            ext_name.contains("stk") || ext_name.contains("vol"))//mrc stack MUST go BEFORE plain MRC
            err = readSPIDER(select_img);
        else if (ext_name.contains("mrcs") || (is_2D && ext_name.contains("mrc")) )//mrc stack MUST go BEFORE plain MRC
            err = readMRC(select_img, true);
        else if (select_img >= 0 && ext_name.contains("mrc"))
        	REPORT_ERROR("Image::read ERROR: stacks of images in MRC-format should have extension .mrcs; .mrc extensions are reserved for 3D maps.");
        else if (ext_name.contains("mrc")) // mrc 3D map
            err = readMRC(select_img, false);
        else if (ext_name.contains("img") || ext_name.contains("hed"))//
            err = readIMAGIC(select_img);//imagic is always an stack
        else
            err = readSPIDER(select_img);

        // Negative errors are bad.
        return err;
    }




    void _write(const FileName &name, fImageHandler* hFile, long int select_img=-1,
                bool isStack=false, int mode=WRITE_OVERWRITE)
    {
        int err = 0;

        FileName ext_name = hFile->ext_name;
        fimg = hFile->fimg;
        fhed = hFile->fhed;
        _exists = hFile->exist;

        filename = name;

        long int aux;
        FileName filNamePlusExt(name);
        name.decompose(aux, filNamePlusExt);
        // Subtract 1 to have numbering 0...N-1 instead of 1...N
        if (aux > 0)
        	aux--;

        if (select_img == -1)
            select_img = aux;

        size_t found = filNamePlusExt.find_first_of("%");

        std::string imParam = "";

        if (found!=std::string::npos)
        {
            imParam =  filNamePlusExt.substr(found+1).c_str();
            filNamePlusExt = filNamePlusExt.substr(0, found) ;
        }

        found = filNamePlusExt.find_first_of(":");
        if ( found!=std::string::npos)
            filNamePlusExt   = filNamePlusExt.substr(0, found);

//#define DEBUG
#ifdef DEBUG

        std::cerr << "write" <<std::endl;
        std::cerr<<"extension for write= "<<ext_name<<std::endl;
        std::cerr<<"filename= "<<filename<<std::endl;
        std::cerr<<"mode= "<<mode<<std::endl;
        std::cerr<<"isStack= "<<isStack<<std::endl;
        std::cerr<<"select_img= "<<select_img<<std::endl;
#endif
#undef DEBUG
        // Check that image is not empty
        if (getSize() < 1)
            REPORT_ERROR("write Image ERROR: image is empty!");

        // CHECK FOR INCONSISTENCIES BETWEEN data.xdim and x, etc???
        int Xdim, Ydim, Zdim;
        long int Ndim;
        this->getDimensions(Xdim,Ydim, Zdim, Ndim);

        Image<T> auxI;
        replaceNsize=0;//reset replaceNsize in case image is reused
        if(select_img == -1 && mode == WRITE_REPLACE)
            REPORT_ERROR("write: Please specify object to be replaced");
        else if(!_exists && mode == WRITE_REPLACE)
        {
            std:: stringstream replace_number;
            replace_number << select_img;
            REPORT_ERROR((std::string)"Cannot replace object number: "
                         + replace_number.str()
                         + " in file " +filename
                         + ". It does not exist");
        }
        else if (_exists && (mode == WRITE_REPLACE || mode == WRITE_APPEND))
        {
            auxI.dataflag = -2;
            auxI.read(filNamePlusExt,false);
            int _Xdim, _Ydim, _Zdim;
            long int _Ndim;
            auxI.getDimensions(_Xdim,_Ydim, _Zdim, _Ndim);
            replaceNsize=_Ndim;
            if(Xdim!=_Xdim ||
               Ydim!=_Ydim ||
               Zdim!=_Zdim
              )
            	REPORT_ERROR("write: target and source objects have different size");
            if(mode==WRITE_REPLACE && select_img>_Ndim)
                REPORT_ERROR("write: cannot replace image stack is not large enough");
            if(auxI.replaceNsize <1 &&
               (mode==WRITE_REPLACE || mode==WRITE_APPEND))
                REPORT_ERROR("write: output file is not an stack");
        }
        else if(!_exists && mode==WRITE_APPEND)
        {
            ;
        }
        else if (mode == WRITE_READONLY)//If new file we are in the WRITE_OVERWRITE mode
        {
            REPORT_ERROR( (std::string) "File " + name
                         + " opened in read-only mode. Cannot write.");
        }

        /*
         * SELECT FORMAT
         */
        if(ext_name.contains("spi") || ext_name.contains("xmp") ||
           ext_name.contains("stk") || ext_name.contains("vol"))
            err = writeSPIDER(select_img,isStack,mode);
        else if (ext_name.contains("mrcs"))
            writeMRC(select_img,true,mode);
        else if (ext_name.contains("mrc"))
            writeMRC(select_img,false,mode);
        else if (ext_name.contains("img") || ext_name.contains("hed"))
            writeIMAGIC(select_img,mode);
        else
            err = writeSPIDER(select_img,isStack,mode);
        if ( err < 0 )
        {
            std::cerr << " Filename = " << filename << " Extension= " << ext_name << std::endl;
            REPORT_ERROR((std::string)"Error writing file "+ filename + " Extension= " + ext_name);
        }

        /* If initially the file did not existed, once the first image is written,
         * then the file exists
         */
        if (!_exists)
            hFile->exist = _exists = true;
    }


};

// Some image-specific operations

// For image normalization
void normalise(Image<double> &I, int bg_radius, double white_dust_stddev, double black_dust_stddev);

void calculateBackgroundAvgStddev(Image<double> &I, double &avg, double &stddev, int bg_radius);

// For dust removal
void removeDust(Image<double> &I, bool is_white, double thresh, double avg, double stddev);

// for contrast inversion
void invert_contrast(Image<double> &I);

// for image re-scaling
void rescale(Image<double> &I, int mysize);

// for image re-windowing
void rewindow(Image<double> &I, int mysize);



/// @defgroup ImageFormats Image Formats
/// @ingroup Images
// Functions belonging to this topic are commented in rw*.h
//@}
#endif