MRC File format

Collaboration diagram for MRC File format:

Classes
struct	MRCheadold
struct	MRChead

Functions
int	readMRC (size_t select_img, bool isStack=false)
int	writeMRC (size_t select_img, bool isStack=false, int mode=WRITE_OVERWRITE, const String &bitDepth="", CastWriteMode castMode=CW_CAST)
int	ImageBase::readMRC (size_t select_img, bool isStack=false)
int	ImageBase::writeMRC (size_t select_img, bool isStack=false, int mode=WRITE_OVERWRITE, const String &bitDepth="", CastWriteMode castMode=CW_CAST)

Detailed Description

Function Documentation

readMRC() [1/2]

int readMRC	(	size_t	select_img,
		bool	isStack = false
	)

MRC Reader

readMRC() [2/2]

int ImageBase::readMRC ( size_t  select_img, bool  isStack = false  )

protected

MRC Reader

Definition at line 352 of file rwMRC.cpp.

{
#undef DEBUG
    //#define DEBUG
#ifdef DEBUG
    printf("DEBUG readMRC: Reading MRC file\n");
#endif

    if (select_img == ALL_IMAGES) {
        return readMRC(1, ALL_IMAGES, isStack);
    }
    return readMRC(select_img, 1, isStack);
}

writeMRC() [1/2]

int writeMRC	(	size_t	select_img,
		bool	isStack = false,
		int	mode = WRITE_OVERWRITE,
		const String &	bitDepth = "",
		CastWriteMode	castMode = CW_CAST
	)

MRC Writer

writeMRC() [2/2]

int ImageBase::writeMRC ( size_t  select_img, bool  isStack = false, int  mode = WRITE_OVERWRITE, const String &  bitDepth = "", CastWriteMode  castMode = CW_CAST  )

protected

MRC Writer

header->a,b,c info is related to sampling rate, so it is only written when writing header, so it is initialized to number of voxels to avoid a mistaken value. If sampling is provided a, b and c are overwritten bellow

Definition at line 369 of file rwMRC.cpp.

{
    MRChead*  header = (MRChead *) askMemory(sizeof(MRChead));

    // Cast T to datatype
    DataType wDType,myTypeID = myT();

    if (bitDepth == "")
    {
        castMode = CW_CAST;
        switch(myTypeID)
        {
        case DT_Double:
        case DT_Float:
        case DT_Int:
        case DT_UInt:
            wDType = DT_Float;
            header->mode = 2;
            break;
        case DT_UShort:
            wDType = DT_UShort;
            header->mode = 6;
            break;
        case DT_Short:
            wDType = DT_Short;
            header->mode = 1;
            break;
        case DT_SChar:
            castMode = CW_CONVERT;
            /* no break */
        case DT_UChar:
            wDType = DT_UChar;
            header->mode = 0;
            break;
        case DT_CFloat:
        case DT_CDouble:
            wDType = DT_CFloat;
            header->mode = 4;
            break;
        case DT_HalfFloat:
            wDType = DT_HalfFloat;
            header->mode = 12;
            break;
        //case DT_UHalfByte:
        default:
            wDType = DT_Unknown;
            (void)wDType; // to suppress dead assignment warning
            REPORT_ERROR(ERR_TYPE_INCORRECT,(std::string)"ERROR: Unsupported data type by MRC format.");
        }
    }
    else //Convert to other data type
    {
        // Default Value
        wDType = (bitDepth == "default") ? DT_Float : datatypeRAW(bitDepth);

        switch (wDType)
        {
        case DT_Double:
        case DT_Int:
        case DT_UInt:
        case DT_Float:
            header->mode = 2;
            break;
        case DT_UChar:
            header->mode = 0;
            break;
        case DT_UShort:
            header->mode = 6;
            break;
        case DT_Short:
            header->mode = 1;
            break;
        case DT_CFloat:
        case DT_CDouble:
            header->mode = 4;
            break;
        case DT_HalfFloat:
            header->mode = 12;
            break;
        default:
            REPORT_ERROR(ERR_TYPE_INCORRECT,"ERROR: incorrect MRC bits depth value.");
        }
    }

    if (mmapOnWrite)
    {
        MDMainHeader.setValue(MDL_DATATYPE,(int) wDType);
        if (!checkMmapT(wDType))
        {
            if (dataMode < DATA && castMode == CW_CAST) // This means ImageGeneric wants to know which DataType must use in mapFile2Write
                return 0;
            else //Mapping is an extra. When not available, go on and do not report an error.
            {
                /* In this case we cannot map the file because required and feasible datatypes are
                 * not compatible. Then we denote to MapFile2Write the same incoming datatype to
                 * keep using this Image object as usual, without mapping on write.
                 */
                mmapOnWrite = false;
                dataMode = DATA;
                MDMainHeader.setValue(MDL_DATATYPE,(int) myTypeID);

                // In case Image size great then, at least, map the multidimarray
                if (mdaBase->nzyxdim*gettypesize(wDType) > tiff_map_min_size)
                    mdaBase->setMmap(true);

                // Allocate memory for image data (Assume xdim, ydim, zdim and ndim are already set
                //if memory already allocated use it (no resize allowed)
                mdaBase->coreAllocateReuse();

                return 0;
            }
        }
        else
            dataMode = DATA;
    }


    /*
         if ( transform != NoTransform )
             img_convert_fourier(p, CentHerm);
     */

    // Map the parameters
    strncpy(header->map, "MAP ", 4);
    // FIXME TO BE DONE WITH rwCCP4!!
    //set_CCP4_machine_stamp(header->machst);
    char* machine_stamp;
    machine_stamp = (char *)(header->machst);
    if(IsLittleEndian())
    {
        machine_stamp[0] = 68;
        machine_stamp[1] = 65;
    }
    else
    {
        machine_stamp[0] = machine_stamp[1] = 17;
    }
    //                    case LittleVAX:
    //                        machine_stamp[0] = 34;
    //                        machine_stamp[1] = 65;
    //                        break;

    size_t Xdim, Ydim, Zdim, Ndim;
    getDimensions(Xdim, Ydim, Zdim, Ndim);

    header->mx = header->nx = Xdim;
    header->my = header->ny = Ydim;
    header->mz = header->nz = Zdim;

    // Obtaining sampling rate for each dimension and calculating cube size
    // By default sampling rate is 1.0 if no real value was found
    double sampling;
    MDMainHeader.getValueOrDefault(MDL_SAMPLINGRATE_X, sampling, 1.0);
    header->a = (float)(Xdim * sampling);
    MDMainHeader.getValueOrDefault(MDL_SAMPLINGRATE_Y, sampling, 1.0);
    header->b = (float)(Ydim * sampling);
    MDMainHeader.getValueOrDefault(MDL_SAMPLINGRATE_Z, sampling, 1.0);
    header->c = (float)(Zdim * sampling);

    if ( transform == CentHerm )
        header->nx = Xdim/2 + 1;        // If a transform, physical storage is nx/2 + 1

    header->alpha = 90.;
    header->beta  = 90.;
    header->gamma = 90.;

    //    header->mx = 0;//(int) (ua/ux + 0.5);
    //    header->my = 0;//(int) (ub/uy + 0.5);
    //    header->mz = 0;//(int) (uc/uz + 0.5);
    header->mapc = 1;
    header->mapr = 2;
    header->maps = 3;
    double aux,aux2;

    //    header->a = 0.;// ua;
    //    header->b = 0.;// ub;
    //    header->c = 0.;// uc;

    if (!MDMainHeader.empty())
    {
#define SET_MAIN_HEADER_VALUE(field, label)  MDMainHeader.getValueOrDefault(label, aux, 0.); header->field = (float)aux
        SET_MAIN_HEADER_VALUE(amin, MDL_MIN);
        SET_MAIN_HEADER_VALUE(amax, MDL_MAX);
        SET_MAIN_HEADER_VALUE(amean, MDL_AVG);
        SET_MAIN_HEADER_VALUE(arms, MDL_STDDEV);

        if ((dataMode == _HEADER_ALL || dataMode == _DATA_ALL))
        {
#define SET_HEADER_SHIFT(field, label)  MD[0]->getValueOrDefault(label, aux, 0.); header->field = -(int) round(aux)
            SET_HEADER_SHIFT(nxStart, MDL_SHIFT_X);
            SET_HEADER_SHIFT(nyStart, MDL_SHIFT_Y);
            SET_HEADER_SHIFT(nzStart, MDL_SHIFT_Z);
#define SET_HEADER_ORIGIN(field, label1, label2)  MD[0]->getValueOrDefault(label1, aux, 0.);MDMainHeader.getValueOrDefault(label2, aux2, 0.);\
              header->field = (float) (aux * aux2)

            SET_HEADER_ORIGIN(xOrigin, MDL_ORIGIN_X, MDL_SAMPLINGRATE_X);
            SET_HEADER_ORIGIN(yOrigin, MDL_ORIGIN_Y, MDL_SAMPLINGRATE_Y);
            SET_HEADER_ORIGIN(zOrigin, MDL_ORIGIN_Z, MDL_SAMPLINGRATE_Z);

#define SET_HEADER_CELL_DIM(field, label1, dimSize)  MDMainHeader.getValueOrDefault(label1, aux, 0.);\
              header->field = (float) (aux * dimSize)

            SET_HEADER_CELL_DIM(a, MDL_SAMPLINGRATE_X, Xdim);
            SET_HEADER_CELL_DIM(b, MDL_SAMPLINGRATE_Y, Ydim);
            SET_HEADER_CELL_DIM(c, MDL_SAMPLINGRATE_Z, Zdim);
        }
        else
        {
            header->nxStart = header->nyStart = header->nzStart = 0;
            header->xOrigin = header->yOrigin = header->zOrigin = 0;
        }
    }

    header->nsymbt = 0;
    header->nlabl = 10; // or zero?
    //strncpy(header->labels, p->label.c_str(), 799);

    offset = MRCSIZE + header->nsymbt;
    size_t datasize, datasize_n;
    datasize_n = Xdim*Ydim*Zdim;
    datasize = datasize_n * gettypesize(wDType);

    //#define DEBUG
#ifdef DEBUG

    printf("DEBUG rwMRC: Offset = %ld,  Datasize_n = %ld\n", offset, datasize_n);
#endif

    size_t imgStart = 0;

    if (Ndim > 1 || filename.contains(":mrcs")) // If format is forced through ":" flag suffix, then ignore the stack behavior in header
        isStack = true;

    bool isVolStk = isStack && Zdim > 1;
    size_t nDimHeader = Ndim;

    if (isStack)
    {
        imgStart = IMG_INDEX(select_img);

        if( mode == WRITE_APPEND )
        {
            imgStart = replaceNsize;
            nDimHeader = replaceNsize + Ndim;
        }
        else if( mode == WRITE_REPLACE && select_img + Ndim - 1 > replaceNsize)
        {
            nDimHeader = select_img + Ndim - 1;
        }
        //        else if (Ndim > replaceNsize)
        //            nDimHeader = Ndim;


        if (isVolStk)
        {
            header->ispg = 401;
            header->mz = Zdim;
            header->nz = Zdim * nDimHeader;
        }
        else
        {
            header->ispg = 0;
            header->nz = nDimHeader;
        }
    }
    else // To set in the header that the file is a volume not a stack
        header->ispg = (Zdim>1)? 1:0;

    //locking
    FileLock flock;
    flock.lock(fimg);

    // Write header when needed
    if(!isStack || replaceNsize < nDimHeader)
    {
        if ( swapWrite )
            swapPage((char *) header, MRCSIZE - 800, DT_Float);
        fwrite( header, MRCSIZE, 1, fimg );
    }
    freeMemory(header, sizeof(MRChead) );

    // Jump to the selected imgStart position
    fseek( fimg,offset + (datasize)*imgStart, SEEK_SET);

    size_t imgEnd = (isStack)? Ndim : 1;

    if (checkMmapT(wDType) && !mmapOnWrite && dataMode >= DATA) {
        writeData(fimg, 0, wDType, datasize_n * imgEnd, castMode);
    } else {
        for ( size_t i = 0; i < imgEnd; i++ )
        {
            // If to also write the image data or jump its size
            if (dataMode >= DATA)
            {
                if (mmapOnWrite && Ndim == 1) // Can map one image at a time only
                {
                    mappedOffset = ftell(fimg);
                    mappedSize = mappedOffset + datasize;
                    fseek(fimg, datasize-1, SEEK_CUR);
                    fputc(0, fimg);
                }
                else
                    writeData(fimg, i*datasize_n, wDType, datasize_n, castMode);
            }
            else
                fseek(fimg, datasize, SEEK_CUR);
        }    
    }

    

    // Unlock the file
    flock.unlock();

    if (mmapOnWrite)
        mmapFile();

    return(0);
}