ProgCtfGroup Class Reference

CTFGroup program class. More...

#include <ctf_group.h>

Inheritance diagram for ProgCtfGroup:

Collaboration diagram for ProgCtfGroup:

Public Member Functions
void	readParams ()
void	defineParams ()
void	show ()
void	usage ()
bool	isIsotropic (CTFDescription &ctf)
void	produceSideInfo ()
void	autoRun ()
void	manualRun ()
void	simpleRun ()
void	writeOutputToDisc ()
void	run ()
Public Member Functions inherited from XmippProgram
const char *	getParam (const char *param, int arg=0)
const char *	getParam (const char *param, const char *subparam, int arg=0)
int	getIntParam (const char *param, int arg=0)
int	getIntParam (const char *param, const char *subparam, int arg=0)
double	getDoubleParam (const char *param, int arg=0)
double	getDoubleParam (const char *param, const char *subparam, int arg=0)
float	getFloatParam (const char *param, int arg=0)
float	getFloatParam (const char *param, const char *subparam, int arg=0)
void	getListParam (const char *param, StringVector &list)
int	getCountParam (const char *param)
bool	checkParam (const char *param)
bool	existsParam (const char *param)
void	addParamsLine (const String &line)
void	addParamsLine (const char *line)
ParamDef *	getParamDef (const char *param) const
virtual void	quit (int exit_code=0) const
virtual int	tryRun ()
void	initProgress (size_t total, size_t stepBin=60)
void	setProgress (size_t value=0)
void	endProgress ()
void	processDefaultComment (const char *param, const char *left)
void	setDefaultComment (const char *param, const char *comment)
virtual void	initComments ()
void	setProgramName (const char *name)
void	addUsageLine (const char *line, bool verbatim=false)
void	clearUsage ()
void	addExampleLine (const char *example, bool verbatim=true)
void	addSeeAlsoLine (const char *seeAlso)
void	addKeywords (const char *keywords)
const char *	name () const
virtual void	usage (int verb=0) const
virtual void	usage (const String &param, int verb=2)
int	version () const
virtual void	show () const
virtual void	read (int argc, const char **argv, bool reportErrors=true)
virtual void	read (int argc, char **argv, bool reportErrors=true)
void	read (const String &argumentsLine)
	XmippProgram ()
	XmippProgram (int argc, const char **argv)
virtual	~XmippProgram ()

Public Attributes
FileName	fn_ctfdat
	Filenames. More...
FileName	fn_root
FileName	fn_split
FileName	format
double	max_error
	Maximum allowed error. More...
double	resol_error
	Resolution for maximum allowed error (in dig. freq.) More...
int	iresol_error
	Resolution for maximum allowed error (in pixels) More...
double	pixel_size
bool	phase_flipped
	Flag for phase-flipped data. More...
size_t	dim
size_t	xpaddim
size_t	ypaddim
size_t	paddim
size_t	ctfxpaddim
double	pad
bool	do_auto
bool	do_discard_anisotropy
bool	do_wiener
	Flag for calculating Wiener filters. More...
double	wiener_constant
	Wiener filter constant. More...
bool	replaceSampling
	Replace ctf.param file sampling rate by this. More...
double	samplingRate
	New ctf sampling rate. More...
int	simpleBins
	Simple algorithm. More...
bool	do1Dctf
double	memory
	Available memory, set mmap to on if more memory is needed. More...
bool	mmapOn
	do not use ram but map data to a file More...
MultidimArray< double >	Mwien
	Matrix with denominator term of Wiener filter. More...
MultidimArray< double >	diff
	auxiliary matrices to speed up process More...
MultidimArray< int >	dd
MultidimArray< double >	mics_ctf2d
std::vector< double >	mics_defocus
MetaDataDb	ImagesMD
MetaDataDb	sortedCtfMD
std::vector< std::vector< int > >	pointer_group2mic
std::vector< MDLabel >	groupbyLabels
Public Attributes inherited from XmippProgram
bool	doRun
bool	runWithoutArgs
int	verbose
	Verbosity level. More...
int	debug

Additional Inherited Members
Protected Member Functions inherited from XmippProgram
void	defineCommons ()
Protected Attributes inherited from XmippProgram
int	errorCode
ProgramDef *	progDef
	Program definition and arguments parser. More...
std::map< String, CommentList >	defaultComments
int	argc
	Original command line arguments. More...
const char **	argv

Detailed Description

CTFGroup program class.

Definition at line 40 of file ctf_group.h.

Member Function Documentation

autoRun()

void ProgCtfGroup::autoRun

(

)

Produce the CTF groups automatically

Definition at line 420 of file ctf_group.cpp.

{
    double diff=0.;
    long int ctfMdSize=sortedCtfMD.size();
    int c = XMIPP_MAX(1, ctfMdSize / 60);
    int counter=0;
    bool newgroup;

    //size_t id;
    size_t orderOut, orderIn;
    int defocusGroup;
    int groupNumber = 1;
    std::vector<size_t> vectorID;
    sortedCtfMD.findObjects(vectorID);
    sortedCtfMD.setValueCol(MDL_DEFGROUP,-1);
    //iterate
    std::vector<size_t>::iterator itOut;
    std::vector<size_t>::iterator itIn;
    size_t begin=vectorID.at(0);
    sortedCtfMD.setValue(MDL_DEFGROUP,groupNumber,begin);

    if (verbose!=0)
    {
        std::cout << "\nCompute differences between CTFs" <<std::endl;
        init_progress_bar(ctfMdSize);
    }

    for ( itOut=vectorID.begin()+1 ; itOut < vectorID.end(); itOut++ )
    {
        counter++;
        newgroup = true;
        sortedCtfMD.getValue(MDL_DEFGROUP,defocusGroup,*itOut);
        if(defocusGroup!=-1)
            continue;
        sortedCtfMD.getValue(MDL_ORDER,orderOut,*itOut);//index in mics_ctf2d array
        for ( itIn=vectorID.begin() ; itIn < itOut; itIn++ )
        {
            sortedCtfMD.getValue(MDL_ORDER,orderIn,*itIn);//index in mics_ctf2d array
            for (int iresol=0; iresol<=iresol_error; iresol++)
            {
                //NZYX_ELEM(mics_ctf2d, orderIn, 1, 1, iresol);
                diff = fabs( NZYX_ELEM(mics_ctf2d, orderIn,  0, 0, iresol) -
                             NZYX_ELEM(mics_ctf2d, orderOut, 0, 0, iresol) );
                if (diff > max_error)
                {
                    break;
                }
            }
            if (diff < max_error)
            {
                newgroup=false;
                break;
            }
        }
        if(newgroup)
        {
            groupNumber++;

        }
        sortedCtfMD.setValue(MDL_DEFGROUP,groupNumber,*itOut);
        if (counter % c == 0 && verbose!=0)
            progress_bar(counter);
    }

    progress_bar(ctfMdSize);
}

defineParams()

void ProgCtfGroup::defineParams ( )

virtual

Define parameters.

Reimplemented from XmippProgram.

Definition at line 116 of file ctf_group.cpp.

{
    addUsageLine("Generate CTF (or defocus) groups from a single CTFdat file.");
    addUsageLine("+The automated mode groups all images with absolute difference in");
    addUsageLine("+CTF-values up to a given resolution below a given threshold.");
    addUsageLine("+For example the example bellow groups all images together with");
    addUsageLine("+absolute CTF differences smaller than 0.5 up to 15 Angstroms resolution).");
    addUsageLine("+A complementary manual mode allows to combine different groups or to split)");
    addUsageLine("+groups up even further.");
    addExampleLine("Example of use: Sample using automated mode (resolution = 15 Ang.)",false);
    addExampleLine("   xmipp_ctf_group --ctfdat all_images_new.ctfdat -o CtfGroupsNew/ctfAuto   --wiener --wc -1 --pad 2 --phase_flipped --error 0.5 --resol 15 ");
    addExampleLine("Example of use: Sample using manual mode (after manual editing of ctf_group_split.doc)",false);
    addExampleLine("   xmipp_ctf_group --ctfdat all_images_new.ctfdat -o CtfGroupsNew/ctfManual --wiener --wc -1 --pad 2 --phase_flipped --split CtfGroupsNew/ctf_group_split.doc");


    //    addParamsLine("   -i <sel_file>              : Input selfile");
    addParamsLine("   --ctfdat <ctfdat_file>     : Input CTFdat file for all data");
    addParamsLine("   [-o <oext=\"ctf:stk\">]        : Output root name, you may force format ctf:mrc");
    addParamsLine("   [--pad <float=1>]          : Padding factor ");
    addParamsLine("   [--phase_flipped]          : Output filters for phase-flipped data");
    addParamsLine("   [--discard_anisotropy]     : Exclude anisotropic CTFs from groups");
    addParamsLine("   [--wiener]                 : Also calculate Wiener filters");
    addParamsLine("   [--sampling_rate <s>]       : This sampling rate overwrites the one in the ctf.param files");
    addParamsLine("   [--memory <double=1.>]     : Available memory in Gb");
    addParamsLine("   [--do1Dctf]                : Compute Groups using 1D CTF, select this option is you have many \
                  non astismatic CTFs");
    addParamsLine("   [--wc <float=-1>]          : Wiener-filter constant (if < 0: use FREALIGN default)");
    addParamsLine(" == MODE 1: AUTOMATED: == ");
    addParamsLine("   [--error <float=0.5> ]     : Maximum allowed error");
    addParamsLine("   [--resol <float=-1> ]      : Resol. (in Ang) for error calculation. Default (-1) = Nyquist");
    addParamsLine("   [--simple <bins=-1>]       : Simple algorithm based on bins of size (max-min)/bins");
    addParamsLine(" == MODE 2: MANUAL: == ");
    addParamsLine("   [--split <docfile> ]       : 1-column docfile with defocus values where to split the data ");

}

isIsotropic()

bool ProgCtfGroup::isIsotropic

(

CTFDescription &

ctf

)

Check anisotropy of a single CTF

Definition at line 392 of file ctf_group.cpp.

{
    double cosp, sinp, ctfp, diff;
    Matrix1D<double> freq(2);

    cosp = COSD(ctf.azimuthal_angle);
    sinp = SIND(ctf.azimuthal_angle);

    for (double digres = 0; digres < resol_error; digres+= 0.001)
    {
        XX(freq) = cosp * digres;
        YY(freq) = sinp * digres;
        digfreq2contfreq(freq, freq, pixel_size);
        ctf.precomputeValues(XX(freq), YY(freq));
        ctfp = ctf.getValueAt();
        ctf.precomputeValues(YY(freq), XX(freq));
        diff = ABS(ctfp - ctf.getValueAt());
        if (diff > max_error)
        {
            std::cout<<" Anisotropy!"<<digres<<" "<<max_error<<" "<<diff<<" "<<ctfp
            <<" "<<ctf.getValueAt()<<std::endl;
            return false;
        }
    }
    return true;
}

manualRun()

void ProgCtfGroup::manualRun

(

)

Split based onto defocus values given in a docfile

Definition at line 487 of file ctf_group.cpp.

{
    MetaDataDb DF;
    int groupNumber = 1;
    DF.read(fn_split);
    int counter=0;
    DF.setValueCol(MDL_DEFGROUP,-2);
    sortedCtfMD.setValueCol(MDL_DEFGROUP,-1);
    //#define DEBUG
#ifdef DEBUG

    sortedCtfMD.write("sortedCtfMD1.xmd");
#endif
#undef DEBUG

    MetaDataDb unionMD;
    DF.unionAll(sortedCtfMD);
    unionMD.sort(DF,MDL_CTF_DEFOCUSA,false);
    int n = unionMD.size();
    init_progress_bar(n);
    int c = XMIPP_MAX(1, n / 60);
    int defGroup;

    for (size_t objId : unionMD.ids())
    {
        unionMD.getValue(MDL_DEFGROUP,defGroup,objId);
        if(defGroup==-2)
            groupNumber++;
        else
            unionMD.setValue(MDL_DEFGROUP,groupNumber,objId);
        if (counter % c == 0)
            progress_bar(counter);
    }
    progress_bar(n);
    //#define DEBUG
#ifdef DEBUG

    unionMD.write("unionMD.xmd");
#endif
#undef DEBUG

    sortedCtfMD.importObjects(unionMD, MDValueNE(MDL_DEFGROUP, -2));
    //#define DEBUG
#ifdef DEBUG

    sortedCtfMD.write("sortedCtfMD3.xmd");
#endif
#undef DEBUG

}

produceSideInfo()

void ProgCtfGroup::produceSideInfo

(

)

Produce side information.

Definition at line 153 of file ctf_group.cpp.

{
    FileName  fnt_ctf,aux;
    CTFDescription ctf;
    MultidimArray<double> Mctf;
    MultidimArray<std::complex<double> >  ctfmask;
    size_t ydim, zdim, ndim;
    double avgdef;

    ImagesMD.read(fn_ctfdat);
    getImageSize(ImagesMD,dim,ydim,zdim,ndim);
    //set output format
    ImagesMD.getValue(MDL_IMAGE,aux,ImagesMD.firstRowId());
    if(format=="")
        format=aux.getFileFormat();

    if ( dim != ydim )
        REPORT_ERROR(ERR_MULTIDIM_SIZE,"Only squared images are allowed!");

    if (simpleBins>0)
        return;

    paddim=xpaddim = ROUND(pad*dim);
    ypaddim=1;
    ctfxpaddim =  (size_t)(sqrt(2.) *  xpaddim + 1);
    Mctf.resize(ypaddim,ctfxpaddim);

    if (do_wiener)
    {
        Mwien.resize(paddim,paddim);
        Mwien.initZeros();
    }

    diff.resize(paddim,paddim);
    dd.resize(paddim,paddim);
    {
        double d;
        int ii,jj;
        FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(diff)
        {
            if (i <= paddim/2 )
                ii=i;
            else
                ii=(paddim-i);

            if (j <= paddim/2)
                jj=j;
            else
                jj=(paddim-j);

            d      = sqrt(ii*ii+jj*jj);
            auto idd = (int) d ;
            dAij(diff, i,j)=d-idd;
            dAij(dd, i,j)=idd;
        }
    }


    MetaDataDb ctfMD;
    groupCTFMetaData(ImagesMD, ctfMD, groupbyLabels);

    int nCTFs = ctfMD.size();
    //how much memory do I need to store them
    double _sizeGb = (double) ypaddim * xpaddim * sizeof(double) * nCTFs /1073741824.;
    mmapOn = _sizeGb > memory;
    mics_ctf2d.setMmap(mmapOn);
    mics_ctf2d.resize(nCTFs,1,ypaddim,ctfxpaddim);

    int c = XMIPP_MAX(1, nCTFs / 60);
    init_progress_bar(nCTFs);

    //use this sampling instead of the one in the CTFparam file
    if(replaceSampling)
    {
        ctfMD.setValueCol(MDL_CTF_SAMPLING_RATE, samplingRate);
    }
    ctf.readFromMetadataRow(ctfMD,ctfMD.firstRowId());

    //do not read directly Tm from metadata because it may be not there
    pixel_size = ctf.Tm;

    if (do_auto)
    {
        if (resol_error < 0)
        {
            // Set to Nyquist
            resol_error = 2. * pixel_size;
        }
        // Set resolution limits in dig freq:
        resol_error = pixel_size / resol_error;
        resol_error = XMIPP_MIN(0.5, resol_error);
        // and in pixels:

        iresol_error = ROUND(resol_error * paddim);
    }

    //fill multiarray with ctfs
    size_t count;
    size_t counter=0;
    if (verbose!=0)
        std::cout << "\nFill multiarray with ctfs" <<std::endl;
    for (size_t objId : ctfMD.ids())
    {
        ctf.readFromMetadataRow(ctfMD, objId);
        ctf.enable_CTF = true;
        ctf.enable_CTFnoise = false;
        ctf.produceSideInfo();
        if (pixel_size != ctf.Tm)
            REPORT_ERROR(ERR_VALUE_INCORRECT,
                         "Cannot mix CTFs with different sampling rates!");
        ctf.Tm /= sqrt(2.);
        if (!do_discard_anisotropy || isIsotropic(ctf))
        {
            avgdef = (ctf.DeltafU + ctf.DeltafV)/2.;
            ctf.DeltafU = avgdef;
            ctf.DeltafV = avgdef;
            ctf.generateCTF(ypaddim, ctfxpaddim, ctfmask);
            FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(ctfmask)
            {
                if (phase_flipped)
                    dAij(Mctf, i, j) = fabs(dAij(ctfmask, i, j).real());
                else
                    dAij(Mctf, i, j) =      dAij(ctfmask, i, j).real();
            }

            //#define DEBUG
#ifdef  DEBUG
            {
                MetaDataVec md1;
                size_t id;
                static int counter=0;

                FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(Mctf)
                {
                    id=md1.addObject();
                    md1.setValue(MDL_ORDER,(size_t)(id-1),id);
                    md1.setValue(MDL_RESOLUTION_FRC,dAij(Mctf, 0, j),id );
                }
                std::stringstream ss;
                ss << counter;//add number to the stream
                md1.write(std::string("NEW/CTF_") +  ss.str());
                counter++;
            }
            ctf.write("new.ctfparam");
#endif
#undef DEBUG
            // Fill vectors
            ctfMD.setValue(MDL_ORDER,counter, objId);
            ctfMD.setValue(MDL_CTF_DEFOCUSA,avgdef, objId);
            mics_ctf2d.setSlice(0,Mctf,counter++);
        }
        else
        {
            std::cout<<" Discard CTF "<<fnt_ctf<<" because of too large anisotropy"<<std::endl;
            ctfMD.removeObject(objId);
        }
        if (counter % c == 0 && verbose!=0)
            progress_bar(counter);
    }
    // Precalculate denominator term of the Wiener filter
    if (do_wiener)
    {

        if (verbose!=0)
            std::cout << "\nPrecalculate denominator term of the Wiener filter" <<std::endl;

        double sumimg = 0.;
        double result;
        for (size_t objId : ctfMD.ids())
        {
            ctfMD.getValue(MDL_COUNT,count, objId);
            ctfMD.getValue(MDL_ORDER,counter,objId);
            auto dCount = (double)count;
            sumimg += dCount;
            FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(Mwien)
            {
                //change DIRECT_N__X_ELEM by DIRECT_N_YX_ELEM if you want to process a 2D ctf
                result =         dAij(diff, i,j)  * DIRECT_N__X_ELEM(mics_ctf2d, counter, dAij(dd,i,j)+1) +
                                 (1.-dAij(diff, i,j)) * DIRECT_N__X_ELEM(mics_ctf2d, counter, dAij(dd,i,j));
                dAij(Mwien,i,j) += dCount * result *result;

            }
            //#define DEBUG
#ifdef DEBUG

            {
                std::cerr << "no_micro_per_Group order " << count
                << " " << counter
                << " " << NZYX_ELEM(mics_ctf2d, counter, 0, 0, 141)
                << " " << dAij(Mwien,100,100)
                << std::endl;
                Image<double> save;
                save()=Mwien;
                save.write("vienertempnew.spi");
                std::cout << "Press any key\n";
                char c;
                std::cin >> c;
            }
#endif
#undef DEBUG

        }
        // Divide by sumimg (Wiener filter is for summing images, not averaging!)
        //#define DEBUG
#ifdef DEBUG

        {
            Image<double> save;
            save()=Mwien;
            save.write("vienerB.spi");
        }
#endif
#undef DEBUG
        Mwien /= sumimg;
        // Add Wiener constant
        if (wiener_constant < 0.)
        {
            // Use Grigorieff's default for Wiener filter constant: 10% of average over all Mwien terms
            // Grigorieff JSB 157(1) (2006), pp 117-125
            wiener_constant = 0.1 * Mwien.computeAvg();
        }
        Mwien += wiener_constant;
        //#define DEBUG
#ifdef DEBUG

        {
            Image<double> save;
            save()=Mwien;
            save.write("vienerA.spi");
        }
#endif
#undef DEBUG

    }
    // Sort by average defocus
    sortedCtfMD.sort(ctfMD,MDL_CTF_DEFOCUSA,false);

}

readParams()

void ProgCtfGroup::readParams ( )

virtual

Read parameters from command line.

Reimplemented from XmippProgram.

Definition at line 34 of file ctf_group.cpp.

{
    //    fn_sel        = getParam("-i");
    fn_ctfdat     = getParam("--ctfdat");
    fn_root       = getParam("-o");

    format = fn_root.getFileFormat();
    fn_root = fn_root.removeFileFormat();

    phase_flipped = checkParam("--phase_flipped");
    do_discard_anisotropy = checkParam("--discard_anisotropy");
    do_auto       = !checkParam("--split");
    pad           = XMIPP_MAX(1., getDoubleParam("--pad"));
    if (do_auto)
    {
        max_error     = getDoubleParam("--error");
        resol_error   = getDoubleParam("--resol");
    }
    else
    {
        fn_split = getParam("--split");
    }
    do_wiener        = checkParam("--wiener");
    replaceSampling  = checkParam("--sampling_rate");
    if(replaceSampling)
        samplingRate       =  getDoubleParam("--sampling_rate");
    memory           = getDoubleParam("--memory");
    do1Dctf          = checkParam("--do1Dctf");
    wiener_constant  = getDoubleParam("--wc");
    simpleBins       = getIntParam("--simple");
}

run()

void ProgCtfGroup::run ( )

virtual

Do the job

Reimplemented from XmippProgram.

Definition at line 758 of file ctf_group.cpp.

{
    produceSideInfo();
    if (simpleBins>0)
        simpleRun();
    else
    {
        if (do_auto)
            autoRun();
        else
            manualRun();
        writeOutputToDisc();
    }
}

show()

void ProgCtfGroup::show

(

)

Show.

Definition at line 67 of file ctf_group.cpp.

{
    std::cout << "  Input ctfdat file       : "<< fn_ctfdat << std::endl;
    std::cout << "  Output rootname         : "<< fn_root << std::endl;
    if (pad > 1.)
    {
        std::cout << "  Padding factor          : "<< pad << std::endl;
    }
    if (do_discard_anisotropy)
    {
        std::cout << " -> Exclude anisotropic CTFs from the groups"<<std::endl;
    }
    if (do_auto)
    {
        std::cout << " -> Using automated mode for making groups"<<std::endl;
        std::cout << " -> With a maximum allowed error of "<<max_error
        <<" at "<<resol_error<<" dig. freq."<<std::endl;
    }
    else
    {
        std::cout << " -> Group based on defocus values in "<<fn_split<<std::endl;
    }
    if (phase_flipped)
    {
        std::cout << " -> Assume that data are PHASE FLIPPED"<<std::endl;
    }
    else
    {
        std::cout << " -> Assume that data are NOT PHASE FLIPPED"<<std::endl;
    }
    if (do_wiener)
    {
        std::cout << " -> Also calculate Wiener filters, with constant= "<<wiener_constant<<std::endl;
    }
    std::cout << "  Available memory (Gb)        : "<< memory << std::endl;
    if (do1Dctf)
    {
        std::cout << " -> compute CTF groups using 1D CTFs= "<<std::endl;
    }
    else
    {
        std::cout << " -> compute CTF groups using 2D CTFs= "<<std::endl;
    }
    if (simpleBins>0)
        std::cout << "Simple algorithm: " << simpleBins << std::endl;
    std::cout << "----------------------------------------------------------"<<std::endl;
}

simpleRun()

void ProgCtfGroup::simpleRun

(

)

Split based on bins

Definition at line 721 of file ctf_group.cpp.

{
    // Maximum and minimum defocus
    double minDefocus=ImagesMD.getColumnMin(MDL_CTF_DEFOCUSU);
    double maxDefocus=ImagesMD.getColumnMax(MDL_CTF_DEFOCUSU);
    double iStepDefocus=1.0/((maxDefocus-minDefocus)/simpleBins);

    double defocus;
    for (size_t objId : ImagesMD.ids())
    {
        ImagesMD.getValue(MDL_CTF_DEFOCUSU,defocus, objId);
        int defGroup=floor((defocus-minDefocus)*iStepDefocus)+1;
        ImagesMD.setValue(MDL_DEFGROUP,defGroup, objId);
    }

    unlink( (fn_root+"_images.sel").c_str());
    FileName imagesInDefoculGroup;
    MetaDataVec auxMetaData;
    auxMetaData.clear();
    auxMetaData.setComment("images (particles) per defocus group, block name is defocusgroup No");
    FileName fnGroup;
    if (verbose>0)
    {
        std::cerr << "Writing defocus groups ...\n";
        init_progress_bar(simpleBins);
    }
    for(int i=1;i<=simpleBins+1; i++)
    {
        auxMetaData.importObjects(ImagesMD,MDValueEQ(MDL_DEFGROUP,i));
        fnGroup.assign( formatString("ctfGroup%06d@%s_images.sel", i, fn_root.c_str()) );
        auxMetaData.write( fnGroup, i > 1 ? MD_APPEND : MD_OVERWRITE);
        if (verbose>0)
            progress_bar(i);
    }
    progress_bar(simpleBins);
}

usage()

void ProgCtfGroup::usage

(

)

Usage.

writeOutputToDisc()

void ProgCtfGroup::writeOutputToDisc

(

)

Write output

Definition at line 538 of file ctf_group.cpp.

{
    //(1) compute no of micrographs, no of images , minimum defocus ,maximum defocus, average defocus per ctf group
    MetaDataDb ctfInfo,ctfImagesGroup,auxMetaData;

    const AggregateOperation MyaggregateOperations[] =
        {
            AGGR_COUNT, AGGR_SUM, AGGR_MIN,         AGGR_MAX,          AGGR_AVG
        };
    std::vector<AggregateOperation> aggregateOperations(MyaggregateOperations,MyaggregateOperations+5);

    const MDLabel MyoperateLabels[]       =
        {
            MDL_COUNT,MDL_COUNT, MDL_CTF_DEFOCUSA, MDL_CTF_DEFOCUSA, MDL_CTF_DEFOCUSA
        };
    std::vector<MDLabel> operateLabels(MyoperateLabels,MyoperateLabels+5);

    const MDLabel MyresultLabels[]        =
        {
            MDL_DEFGROUP,MDL_COUNT, MDL_SUM,  MDL_MIN,          MDL_MAX,          MDL_AVG
        };
    std::vector<MDLabel> resultLabels(MyresultLabels,MyresultLabels+6);

    ctfInfo.aggregate(sortedCtfMD,aggregateOperations,operateLabels,resultLabels);
    ctfInfo.setComment("N. of micrographs, N. of particles, min defocus, max defocus and avg defocus");
    ctfInfo.write("groups@"+fn_root+"Info.xmd");
    size_t numberDefGroups=ctfInfo.size();
    MetaDataVec MD;
    size_t idctf = MD.addObject();
    MD.setValue(MDL_COUNT,numberDefGroups,idctf);
    MD.setColumnFormat(false);
    MD.write("numberGroups@"+fn_root+"Info.xmd",MD_APPEND);

    //(2)save auxiliary file for defocus split

    double maxDef,minDef;
    auto it = ctfInfo.ids().begin();
    size_t id1,id2,id;
    auxMetaData.clear();
    auxMetaData.setComment(formatString("Defocus values to split into %lu ctf groups", ctfInfo.size()));

    id1 = *it;
    ++it;
    const auto totalSize = ctfInfo.ids().end();
    while (it != totalSize)
    {
        id2 = *it;
        ctfInfo.getValue(MDL_MIN,minDef,id1);
        ctfInfo.getValue(MDL_MAX,maxDef,id2);
        id1=id2;

        id=auxMetaData.addObject();
        auxMetaData.setValue(MDL_CTF_DEFOCUSA,(minDef+maxDef)/2.,id);
        ++it;
    }
    auxMetaData.write(fn_root+"_split.doc");


    //(3) make block-sel per image group
    ImagesMD.read(fn_ctfdat);
    if (ImagesMD.containsLabel(MDL_DEFGROUP))
        ImagesMD.removeLabel(MDL_DEFGROUP);
    if (ImagesMD.containsLabel(MDL_CTF_MODEL))
        ctfImagesGroup.join1(ImagesMD, sortedCtfMD, MDL_CTF_MODEL, INNER );
    else
    {
        ctfImagesGroup.join1(ImagesMD,
                             sortedCtfMD,
                             groupbyLabels,LEFT);
        //ctfImagesGroup.joinNatural(ImagesMD, sortedCtfMD);
//#define DEBUG
#ifdef DEBUG
    ImagesMD.write("/tmp/ImagesMD.sqlite");
    sortedCtfMD.write("/tmp/sortedCtfMD.sqlite");
    ctfImagesGroup.write("/tmp/ctfImagesGroup.sqlite");
    ImagesMD.write("/tmp/ImagesMD.xmd");
    sortedCtfMD.write("/tmp/sortedCtfMD.xmd");
    ctfImagesGroup.write("/tmp/ctfImagesGroup.xmd");

#endif
    }

    //
    unlink( (fn_root+"_images.sel").c_str());
    FileName imagesInDefoculGroup;
    auxMetaData.clear();
    auxMetaData.setComment("images (particles) per defocus group, block name is defocusgroup No");
    int ctfInfoSize;
    ctfInfoSize = (int) ctfInfo.size();
    for(int i=1;i<= ctfInfoSize; i++)
    {
        auxMetaData.importObjects(ctfImagesGroup,MDValueEQ(MDL_DEFGROUP,i));
        imagesInDefoculGroup.assign( formatString("ctfGroup%06d@%s_images.sel", i, fn_root.c_str()) );
        auxMetaData.write( imagesInDefoculGroup, i > 1 ? MD_APPEND : MD_OVERWRITE);
    }

    //(4)create average ctf
    int olddefGroup=-1;
    int defGroup=-1;
    size_t order, count;
    double sumimg=0.;

    MultidimArray<double> ctf2D(paddim,paddim);
    Image<double> Ictf2D;
    Ictf2D.data.alias(ctf2D);

    FileName outFileNameCTF,outFileNameWIEN,outFileName;
    outFileNameCTF = fn_root + "_ctf."+format;
    outFileNameWIEN = fn_root + "_wien."+format;
    if (verbose!=0)
    {
        std::cout << "Saving CTF Images" <<std::endl;
    }
    for (size_t objId : sortedCtfMD.ids())
    {
        sortedCtfMD.getValueOrDefault(MDL_DEFGROUP,defGroup,objId,-1);
        if (olddefGroup<0)
            olddefGroup=defGroup;
        sortedCtfMD.getValue(MDL_ORDER,order,objId);
        sortedCtfMD.getValue(MDL_COUNT,count,objId);
        auto dCount = (double)count;

        if (defGroup != olddefGroup)
        {
            if (sumimg!=0)
            {
                ctf2D /= sumimg;
                outFileName.compose(olddefGroup,outFileNameCTF);
                //save CTF
                Ictf2D.write(outFileName);
                //save winer filter
                if (do_wiener)
                {
                    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(ctf2D)
                    {
                        dAij(ctf2D,i,j) /= dAij(Mwien,i,j);
                    }
                    outFileName.compose(olddefGroup,outFileNameWIEN);
                    Ictf2D.write(outFileName);
                }

                ctf2D.initZeros();
                olddefGroup=defGroup;
                sumimg=0.;
            }
        }

        sumimg += dCount;
        FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(ctf2D)
        {
            double result;
            //interpolate ctf point from table
            //change DIRECT_N__X_ELEM by DIRECT_N_YX_ELEM if you want to process a 2D ctf

            result =     dAij(diff, i,j)      * DIRECT_N__X_ELEM(mics_ctf2d, order, dAij(dd,i,j)+1) +
                         (1.-dAij(diff, i,j)) * DIRECT_N__X_ELEM(mics_ctf2d, order, dAij(dd,i,j));
            dAij(ctf2D,i,j) += dCount * result ;
        }
    }
    //Save last CTF
    if (defGroup == olddefGroup)
    {
        if (sumimg!=0)
        {
            ctf2D /= sumimg;
            outFileName.compose(defGroup,outFileNameCTF);
            //save CTF
            Ictf2D.write(outFileName);
            //save winer filter
            if (do_wiener)
            {
                FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(ctf2D)
                {
                    dAij(ctf2D,i,j) /= dAij(Mwien,i,j);
                }
                outFileName.compose(defGroup,outFileNameWIEN);
                Ictf2D.write(outFileName);
            }
        }
    }

}

Member Data Documentation

ctfxpaddim

size_t ProgCtfGroup::ctfxpaddim

Definition at line 63 of file ctf_group.h.

dd

MultidimArray<int> ProgCtfGroup::dd

Definition at line 104 of file ctf_group.h.

diff

MultidimArray<double> ProgCtfGroup::diff

auxiliary matrices to speed up process

Definition at line 103 of file ctf_group.h.

dim

size_t ProgCtfGroup::dim

Definition at line 63 of file ctf_group.h.

do1Dctf

bool ProgCtfGroup::do1Dctf

Compute 1D CTF using avg(defocusU + defocusV)/2 as defocus. This approach speed up the computation and is recommended for very large data sets

Definition at line 91 of file ctf_group.h.

do_auto

bool ProgCtfGroup::do_auto

Definition at line 69 of file ctf_group.h.

do_discard_anisotropy

bool ProgCtfGroup::do_discard_anisotropy

Definition at line 72 of file ctf_group.h.

do_wiener

bool ProgCtfGroup::do_wiener

Flag for calculating Wiener filters.

Definition at line 75 of file ctf_group.h.

fn_ctfdat

FileName ProgCtfGroup::fn_ctfdat

Filenames.

Definition at line 45 of file ctf_group.h.

fn_root

FileName ProgCtfGroup::fn_root

Definition at line 45 of file ctf_group.h.

fn_split

FileName ProgCtfGroup::fn_split

Definition at line 45 of file ctf_group.h.

format

FileName ProgCtfGroup::format

Definition at line 45 of file ctf_group.h.

groupbyLabels

std::vector<MDLabel> ProgCtfGroup::groupbyLabels

Definition at line 129 of file ctf_group.h.

ImagesMD

MetaDataDb ProgCtfGroup::ImagesMD

Definition at line 114 of file ctf_group.h.

iresol_error

int ProgCtfGroup::iresol_error

Resolution for maximum allowed error (in pixels)

Definition at line 54 of file ctf_group.h.

max_error

double ProgCtfGroup::max_error

Maximum allowed error.

Definition at line 48 of file ctf_group.h.

memory

double ProgCtfGroup::memory

Available memory, set mmap to on if more memory is needed.

Definition at line 94 of file ctf_group.h.

mics_ctf2d

MultidimArray<double> ProgCtfGroup::mics_ctf2d

Definition at line 108 of file ctf_group.h.

mics_defocus

std::vector<double> ProgCtfGroup::mics_defocus

Definition at line 111 of file ctf_group.h.

mmapOn

bool ProgCtfGroup::mmapOn

do not use ram but map data to a file

Definition at line 97 of file ctf_group.h.

Mwien

MultidimArray<double> ProgCtfGroup::Mwien

Matrix with denominator term of Wiener filter.

Definition at line 100 of file ctf_group.h.

pad

double ProgCtfGroup::pad

Definition at line 66 of file ctf_group.h.

paddim

size_t ProgCtfGroup::paddim

Definition at line 63 of file ctf_group.h.

phase_flipped

bool ProgCtfGroup::phase_flipped

Flag for phase-flipped data.

Definition at line 60 of file ctf_group.h.

pixel_size

double ProgCtfGroup::pixel_size

Definition at line 57 of file ctf_group.h.

pointer_group2mic

std::vector< std::vector <int> > ProgCtfGroup::pointer_group2mic

Definition at line 126 of file ctf_group.h.

replaceSampling

bool ProgCtfGroup::replaceSampling

Replace ctf.param file sampling rate by this.

Definition at line 81 of file ctf_group.h.

resol_error

double ProgCtfGroup::resol_error

Resolution for maximum allowed error (in dig. freq.)

Definition at line 51 of file ctf_group.h.

samplingRate

double ProgCtfGroup::samplingRate

New ctf sampling rate.

Definition at line 84 of file ctf_group.h.

simpleBins

int ProgCtfGroup::simpleBins

Simple algorithm.

Definition at line 87 of file ctf_group.h.

sortedCtfMD

MetaDataDb ProgCtfGroup::sortedCtfMD

Definition at line 117 of file ctf_group.h.

wiener_constant

double ProgCtfGroup::wiener_constant

Wiener filter constant.

Definition at line 78 of file ctf_group.h.

xpaddim

size_t ProgCtfGroup::xpaddim

Definition at line 63 of file ctf_group.h.

ypaddim

size_t ProgCtfGroup::ypaddim

Definition at line 63 of file ctf_group.h.

---

The documentation for this class was generated from the following files:

• xmipp/libraries/reconstruction/ctf_group.h
• xmipp/libraries/reconstruction/ctf_group.cpp