transform_window.cpp

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'
 ***************************************************************************/

#include "core/xmipp_metadata_program.h"
#include "data/unitCell.h"

class ProgWindow: public XmippMetadataProgram {
public:
    typedef enum {
        CORNERMODE, SIZEMODE, CROPMODE, UNITCELLMODE
    } WindowMode;

    bool physical_coords;
    int sizeX, sizeY, sizeZ;
    int cropX, cropY, cropZ;
    int x0, y0, z0;
    int xF, yF, zF;
    String sym;
    //unit cell option parameters
    double rmin, rmax, expand, offset, sampling, x_origin, y_origin, z_origin;
    double padValue;
    String padType;
    WindowMode windowMode;
    FileName tempFn;
    bool renameTempFn;

    void defineParams() {
        each_image_produces_an_output = true;
        allow_apply_geo = true;
        XmippMetadataProgram::defineParams();
        addUsageLine("This program takes a region from the input images");
        addUsageLine(" ");
        addUsageLine(
                "But... what is windowing? Let's suppose we have an image 32x32,");
        addUsageLine(
                "and we window it to 16x16, then the central part which contains");
        addUsageLine(
                "this 16x16 square is returned, but if we window to 64x64 then the");
        addUsageLine(
                "image is kept in the center of the final image and it is padded");
        addUsageLine(
                "with 0's until the new size is reached (the padding value may be");
        addUsageLine(
                "modified using --pad_value, --corner_pad_value or --average_pad_value).");
        addUsageLine(" ");
        addUsageLine(
                "With this idea in mind of padding or cutting, if you define two");
        addUsageLine(
                "logical corners in the input image (the most negative -the top");
        addUsageLine(
                "left in an image- and the most positive -bottom right in an image-)");
        addUsageLine(
                "you can store the new image or volume generated by windowing with");
        addUsageLine("a square box.");
        addParamsLine(
                "  --corners <...>                        : Windows corners (by default indexes are logical)");
        addParamsLine(
                "                                         : 2D: <x0> <y0> <xF> <yF>");
        addParamsLine(
                "                                         : 3D: <x0> <y0> <z0> <xF> <yF> <zF>");
        addParamsLine(
                "  or --size <sizeX> <sizeY=0> <sizeZ=0>  : Output dimension. The volume is windowed");
        addParamsLine(
                "                                         : (expanded or cutted) in all directions such");
        addParamsLine(
                "                                         : that the origin of the volume remains the");
        addParamsLine(
                "                                         : same. If the Y and Z dimensions are not");
        addParamsLine(
                "                                         : specified they are assumed to be the same as");
        addParamsLine(
                "                                         : the X dimension.");
        addParamsLine(
                "  or --crop <sizeX> <sizeY=0> <sizeZ=0>  : Crop this amount of pixels in each direction.");
        addParamsLine(
                "                                         : Half of the pixels will be cropped from the left");
        addParamsLine(
                "                                         : and the other half from the right");
        addParamsLine(
                "                                         : if only one is given, the other two");
        addParamsLine(
                "                                         : are supposed to be the same");
        addParamsLine(
                "  or --unitcell <sym> <rmin=0> <rmax=0> <expandFactor=0> <offset=0> <sampling=1.> <x_origin=-1.> <y_origin=-1.> <z_origin=-1.> : Extract a unit cell from volume");
        addParamsLine(
                "                                         : sym = particle symmetry");
        addParamsLine(
                "                                         : expandFactor = expand unitcell by this factor");
        addParamsLine(
                "                                         : rmax, rmin = cut a shell with these radii (pixels)");
        addParamsLine(
                "                                         : offset= for CN symmetry rotate unit cell by this angle");
        addParamsLine(
                "                                         : sampling=if the output is a mrc file it will use this value to fill the header");
        addParamsLine(
                "                                         : x_origin = origin x coordinate introduced by the user with the input volume (in pixels)");
        addParamsLine(
                "                                         : y_origin = origin y coordinate introduced by the user with the input volume (in pixels)");
        addParamsLine(
                "                                         : z_origin = origin z coordinate introduced by the user with the input volume (in pixels)");
        addParamsLine(
                "  [--physical]                           : use physical instead of logical coordinates");
        addParamsLine("    requires --corners;");
        addParamsLine(
                "  [--pad <padtype=value>]                : value used for padding");
        addParamsLine("   where <padtype>");
        addParamsLine(
                "         value <v=0>                     : use this value for padding");
        addParamsLine(
                "         corner                          : use the top-left corner for padding");
        addParamsLine(
                "         avg                             : use the image average for padding");
        addExampleLine(
                "Window a single image to 16x16, overwriting input image",
                false);
        addExampleLine("xmipp_transform_window -i g0ta0001.xmp --size 16");
        addExampleLine("Window a single volume to 32x64x64", false);
        addExampleLine("xmipp_transform_window -i g0ta.vol --size 64 64 32");
        addExampleLine("The same using logical indexes", false);
        addExampleLine(
                "xmipp_transform_window -i g0ta.vol --corners -32 -32 -16 31 31 15");
        addExampleLine(
                "Note that r0 and rF are not symmetric because the volume is of an",
                false);
        addExampleLine(
                "even size, if we wanted to get a 33x65x65 volume, the right indexes",
                false);
        addExampleLine("would be", false);
        addExampleLine(
                "xmipp_transform_window -i g0ta.vol --corners -32 -32 -16 32 32 16");
        addExampleLine("Reduce the volume by 10 pixels on each direction",
                false);
        addExampleLine("xmipp_transform_window -i g0ta.vol --crop 10");
        addExampleLine(
                "Enlarge the volume by 10 pixels on each direction (negative crop)",
                false);
        addExampleLine("xmipp_transform_window -i g0ta.vol --crop -10");
        addExampleLine("scipion xmipp_transform_window -i postprocess.mrc:mrc  -o postprocess_win.mrc  --unitcell i2 0 300 20 0 0.34");
        addKeywords("window, crop, resize, corner, padding");
    }

    void readParams() {
        XmippMetadataProgram::readParams();
        padType = getParam("--pad");
        if (padType == "value")
            padValue = getDoubleParam("--pad", "value");
        if (checkParam("--corners")) {
            int nparams = getCountParam("--corners");
            if (nparams == 4 || nparams == 6) {
                x0 = getIntParam("--corners", 0);
                y0 = getIntParam("--corners", 1);
                if (nparams == 4) {
                    xF = getIntParam("--corners", 2);
                    yF = getIntParam("--corners", 3);
                    zF = z0 = 0;
                } else {
                    z0 = getIntParam("--corners", 2);
                    xF = getIntParam("--corners", 3);
                    yF = getIntParam("--corners", 4);
                    zF = getIntParam("--corners", 5);
                }
            } else
                REPORT_ERROR(ERR_ARG_INCORRECT,
                        "Incorrect number of arguments after --corners");
            physical_coords = checkParam("--physical");
            windowMode = CORNERMODE;
        } else if (checkParam("--size")) {
            sizeX = getIntParam("--size", 0);
            sizeY = getIntParam("--size", 1);
            sizeZ = getIntParam("--size", 2);
            if (sizeY == 0)
                sizeY = sizeX;
            if (sizeZ == 0)
                sizeZ = sizeX;
            x0 = FIRST_XMIPP_INDEX(sizeX);
            y0 = FIRST_XMIPP_INDEX(sizeY);
            z0 = FIRST_XMIPP_INDEX(sizeZ);
            xF = LAST_XMIPP_INDEX(sizeX);
            yF = LAST_XMIPP_INDEX(sizeY);
            zF = LAST_XMIPP_INDEX(sizeZ);
            windowMode = SIZEMODE;
            physical_coords = false;
        }
        //<sym> <rmin=0> <rmax=0> <expandFactor=0> <offset=0> <x_offset=-1.> <y_offset=-1.> <z_offset=-1.>: Extract a unit cell from volume");

        else if (checkParam("--crop")) {
            cropX = getIntParam("--crop", 0);
            cropY = getIntParam("--crop", 1);
            cropZ = getIntParam("--crop", 2);
            if (cropY == 0)
                cropY = cropX;
            if (cropZ == 0)
                cropZ = cropX;
            windowMode = CROPMODE;
            physical_coords = false;
        } else if (checkParam("--unitcell")) {
            int nparams = getCountParam("--unitcell");
            if (nparams != 9) REPORT_ERROR(ERR_ARG_INCORRECT, "Incorrect number of arguments after --unitcell");
            sym = getParam("--unitcell", 0);
            rmin = getDoubleParam("--unitcell", 1);
            rmax = getDoubleParam("--unitcell", 2);
            expand = getDoubleParam("--unitcell", 3);
            offset = DEG2RAD(getDoubleParam("--unitcell", 4));
            sampling = getDoubleParam("--unitcell", 5);
            x_origin = getDoubleParam("--unitcell", 6);
            y_origin = getDoubleParam("--unitcell", 7);
            z_origin = getDoubleParam("--unitcell", 8);

            windowMode = UNITCELLMODE;
        }
    }

    void show() {
        if (verbose == 0)
            return;
        XmippMetadataProgram::show();
        switch (windowMode) {
        case SIZEMODE:
            std::cout << "New size: (XxYxZ)=" << sizeX << "x" << sizeY << "x"
                    << sizeZ << std::endl;
            break;
        case CROPMODE:
            std::cout << "Crop: (XxYxZ)=" << cropX << "x" << cropY << "x"
                    << cropZ << std::endl;
            break;
        case CORNERMODE:
            std::cout << "New window: from (z0,y0,x0)=(" << z0 << "," << y0
                    << "," << x0 << ") to (zF,yF,xF)=(" << zF << "," << yF
                    << "," << xF << ")\n" << "Physical: " << physical_coords
                    << std::endl;
            break;
        case UNITCELLMODE:
            std::cout << "Sym: " << sym << ", rmin: " << rmin << ", rmax: " << rmax
                      << ", Expand Factor: " << expand << ", offset: " << offset
                      << ", sampling: " << sampling << "; Origin in Angstroms: "
                      << " x: " << x_origin * sampling << ", y: " << y_origin * sampling << ", z: " << z_origin * sampling
                      << std::endl;
            break;
        }
    }

    void preProcess() {
        create_empty_stackfile = false;
        renameTempFn = false;

        /* Since dimensions are usually changed in this program, we cannot reuse the same file
         * but we can create a temporary file */
        if (input_is_stack && fn_out.empty()) {
//            tempFn.initUniqueName("temp_XXXXXX");
            fn_out = fn_in + "xtw_tmpfile" + "." + fn_in.getExtension();
            renameTempFn = true;
        }
    }

    void unitcell(ImageGeneric &in3Dmap,
                  ImageGeneric & out3DDmap) {
        //1) init class and compute auxiliary vectors
        UnitCell UC(sym, rmin, rmax, expand, offset, sampling, x_origin, y_origin, z_origin);
        //2) mask image
        UC.maskUnitCell(in3Dmap,out3DDmap);
        //3) shift them
        //4) retain only positive cases
    }
    void processImage(const FileName &fnImg, const FileName &fnImgOut,
            const MDRow &rowIn, MDRow &rowOut) {
        ImageGeneric Iin;
        bool createTempFile = (fnImg == fnImgOut);
        if (apply_geo) {
            Iin.readApplyGeo(fnImg, rowIn);
        } else {
            Iin.read(fnImg, DATA, ALL_IMAGES, true);
        }

        double init_value(padValue);
        if (padType == "avg")
            init_value = Iin().computeAvg();
        else if (padType == "corner")
            init_value = Iin()(0, 0, 0, 0);

        Iin().setXmippOrigin();
        if (ZSIZE(Iin()()) == 1)
            zF = z0 = 0;

        String oext = fnImgOut.getExtension();
        DataType dataType = Iin.getDatatype();
        ImageGeneric result(dataType);
        if (windowMode == UNITCELLMODE) {
            Iin.convert2Datatype(DT_Float);
            unitcell(Iin, result);
        } else if (windowMode == CROPMODE) {
            int xl = cropX / 2;
            int xr = cropX - xl;
            int yl = cropY / 2;
            int yr = cropY - yl;
            int zl = cropZ / 2;
            int zr = cropZ - zl;
            if (ZSIZE(Iin()()) == 1) {
                zl = zr = 0;
            }
            //call to a generic 4D function;
            result.mapFile2Write(
            FINISHINGX(Iin()()) - xr - (STARTINGX(Iin()()) + xl) + 1,
            FINISHINGY(Iin()()) - yr - (STARTINGY(Iin()()) + yl) + 1,
            FINISHINGZ(Iin()()) - zr - (STARTINGZ(Iin()()) + zl) + 1, fnImgOut,
                    createTempFile);
            Iin().window(result(),
            STARTINGZ(Iin()()) + zl, STARTINGY(Iin()()) + yl,
            STARTINGX(Iin()()) + xl,
            FINISHINGZ(Iin()()) - zr, FINISHINGY(Iin()()) - yr,
            FINISHINGX(Iin()()) - xr);
        } else if (!physical_coords) {
            result.mapFile2Write(xF - x0 + 1, yF - y0 + 1, zF - z0 + 1,
                    fnImgOut, createTempFile);
            Iin().window(result(), z0, y0, x0, zF, yF, xF, init_value);
        } else {
            result.mapFile2Write(xF - x0 + 1, yF - y0 + 1, zF - z0 + 1,
                    fnImgOut, createTempFile);
            Iin().window(result(),
            STARTINGZ(Iin()()) + z0,
            STARTINGY(Iin()()) + y0,
            STARTINGX(Iin()()) + x0,
            STARTINGZ(Iin()()) + zF,
            STARTINGY(Iin()()) + yF,
            STARTINGX(Iin()()) + xF, init_value);
        }
        Iin.clear(); // Close the input file
        result.write(fnImgOut);
    }

    void postProcess() {
        /* Since dimensions are usually changed in this program, we cannot reuse the same file
         * but we can create a temporary file. Now it must be deleted */
        if (renameTempFn) {
            if (std::rename(fn_out.c_str(), fn_in.c_str()) != 0)
                REPORT_ERROR(ERR_IO,
                        formatString("Error renaming file '%s' to '%s'.",
                                fn_out.c_str(), fn_in.c_str()));
        }
    }
};