Correct Amplitude3D parameters. More...

#include <ctf_correct_wiener3d.h>

Inheritance diagram for ProgCtfCorrectAmplitude3D:

Collaboration diagram for ProgCtfCorrectAmplitude3D:

Public Member Functions
void	readParams ()

void	show ()

void	defineParams ()

void	produceSideInfo ()

void	generateCTF1D (const FileName &fnCTF, size_t nr_steps, MultidimArray< double > &CTF1D)

void	generateWienerFilters ()

void	generateVolumes ()

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	fnIn
	Metadata with volume, ctf and number of images in that volume. More...

FileName	fnRoot
	Rootname for output files. More...

double	wienerConstant
	Wiener filter constant. More...

double	minFreq
	Low resolution cutoff to apply Wiener filter. More...

bool	isFlipped
	Flag for phase flipped images. More...

size_t	Zdim
	Dimensions of the volumes. More...

size_t	Ydim

size_t	Xdim

FourierFilter	ctf
	Side Info: CTF. More...

MetaDataVec	ctfdat
	Side Info: ctfdat. More...

std::vector< MultidimArray< double > >	Vctfs1D
	The 3D CTFs and Wiener filters. More...

std::vector< MultidimArray< double > >	Vwien1D

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

Correct Amplitude3D parameters.

Definition at line 34 of file ctf_correct_wiener3d.h.

Member Function Documentation

◆ defineParams()

void ProgCtfCorrectAmplitude3D::defineParams ( )

virtual

Define Parameters

Reimplemented from XmippProgram.

Definition at line 56 of file ctf_correct_wiener3d.cpp.

 {
     addUsageLine("Wiener filtering of volumes");
     addUsageLine("+The program combines a set of volumes, each one with its one CTF and produces a deconvolved Wiener volume");
 
     addParamsLine("  -i <metadataFile>           : Metadata with the volumes, ctfs, number of images in that group");
     addParamsLine("                              :+The metadata labels are _image, _CTFModel, _class_count");
     addParamsLine("  --oroot <file>              : Output rootname ");
     addParamsLine("                              :+oroot+_deconvolved.vol contains the combination of all volumes");
     addParamsLine("                              :+oroot+_ctffiltered_group01.vol contains each volume obtained after filtering the deconvolved one");
     addParamsLine("                              :+With verbose==2: oroot+_wien01.txt contains the Wiener filters in Fourier");
     addParamsLine("  [--minFreq <Ang=-1>]        : Apply Wiener filter only beyond this resolution (in Angstrom)");
     addParamsLine("  [--phase_flipped]           : Use this if the maps were reconstructed from phase corrected images ");
     addParamsLine("  [--wienerConstant <K=0.05>] : Wiener constant (to be multiplied by the total number of images) ");
     addExampleLine("xmipp_ctf_correct_wiener3d -i ctf_correct3d.xmd --oroot volumeCorrected");
     addExampleLine("In the following link you can find an example of input file:",false);
     addExampleLine(" ",false);
     addExampleLine("http://sourceforge.net/p/testxmipp/code/ci/master/tree/input/ctf_correct3d.xmd?format=raw",false);
 }

◆ generateCTF1D()

void ProgCtfCorrectAmplitude3D::generateCTF1D	(	const FileName &	fnCTF,
		size_t	nr_steps,
		MultidimArray< double > &	CTF1D
	)

Generate 1D CTFs.

Definition at line 94 of file ctf_correct_wiener3d.cpp.

 {
     // Read the CTF
     ctf.FilterBand = CTF;
     ctf.ctf.enable_CTFnoise = false;
     ctf.ctf.read(fnCTF);
     ctf.ctf.produceSideInfo();
 
     double maxres = ( 0.5 * sqrt(3.) ) / ctf.ctf.Tm;
     double stepsize = maxres / (double)nr_steps;
     CTF1D.resizeNoCopy(nr_steps);
     double freq = 0.;
 
     freq=0.;
     for (size_t step=0; step < nr_steps; step++)
     {
         if ( (minFreq < 0) || (1./freq < minFreq) )
         {
             ctf.ctf.precomputeValues(freq, 0.0);
             A1D_ELEM(CTF1D,step)=ctf.ctf.getValueAt();
             if (isFlipped)
                 A1D_ELEM(CTF1D,step)=fabs(CTF1D(step));
         }
         else
             A1D_ELEM(CTF1D,step)=1.;
         freq+=stepsize;
     }
 }

◆ generateVolumes()

void ProgCtfCorrectAmplitude3D::generateVolumes ( )

Generate deconvolved volume

Definition at line 209 of file ctf_correct_wiener3d.cpp.

 {
     Image<double> V;
     FileName fnVol, fnCTF;
     MultidimArray<std::complex<double> > fft,fft_out;
     Matrix1D<int>    idx(3);
     Matrix1D<double> freq(3);
 
     int ii = 0;
     double Xdim2=Xdim*Xdim;
     double Ydim2=Ydim*Ydim;
     double Zdim2=Zdim*Zdim;
     for (size_t objId : ctfdat.ids())
     {
         ctfdat.getValue(MDL_IMAGE,fnVol,objId);
         ctfdat.getValue(MDL_CTF_MODEL,fnCTF,objId);
         V.read(fnVol);
         FourierTransform(V(),fft);
         if (ii == 0)
             fft_out.initZeros(fft);
         MultidimArray<double>& Vwien1D_ii=Vwien1D[ii];
         FOR_ALL_ELEMENTS_IN_ARRAY3D(fft)
         {
             XX(idx) = j;
             YY(idx) = i;
             ZZ(idx) = k;
             FFT_idx2digfreq(fft, idx, freq);
             auto ires= (int)round(OVERSAMPLE*sqrt(XX(freq)*XX(freq)*Xdim2+
                                                  YY(freq)*YY(freq)*Ydim2+
                                                  ZZ(freq)*ZZ(freq)*Zdim2));
             A3D_ELEM(fft_out,k,i,j)+=A1D_ELEM(Vwien1D_ii,ires)*A3D_ELEM(fft,k,i,j);
         }
         ii++;
     }
 
     // Inverse Fourier Transform
     InverseFourierTransform(fft_out,V());
     fnVol=fnRoot+"_deconvolved.vol";
     V.write(fnVol);
 
     // Calculate CTF-affected volumes
     ii = 0;
     for (size_t objId : ctfdat.ids())
     {
         ctfdat.getValue(MDL_IMAGE,fnVol,objId);
         ctfdat.getValue(MDL_CTF_MODEL,fnCTF,objId);
         MultidimArray<double>& Vwien1D_ii=Vwien1D[ii];
         FOR_ALL_ELEMENTS_IN_ARRAY3D(fft)
         {
             XX(idx) = j;
             YY(idx) = i;
             ZZ(idx) = k;
             FFT_idx2digfreq(fft, idx, freq);
             auto ires= (int)round(OVERSAMPLE*sqrt(XX(freq)*XX(freq)*Xdim2+
                                                  YY(freq)*YY(freq)*Ydim2+
                                                  ZZ(freq)*ZZ(freq)*Zdim2));
             A3D_ELEM(fft,k,i,j)+=A1D_ELEM(Vwien1D_ii,ires)*A3D_ELEM(fft_out,k,i,j);
         }
         ii++;
 
         fnVol=fnRoot+"_ctffiltered_group";
         fnVol.compose(fnVol,ii,"vol");
         InverseFourierTransform(fft,V());
         V.write(fnVol);
     }
 }

◆ generateWienerFilters()

void ProgCtfCorrectAmplitude3D::generateWienerFilters ( )

Generate Wiener filters

Definition at line 125 of file ctf_correct_wiener3d.cpp.

 {
     MultidimArray<double> CTF1D, sumterm;
     std::ofstream  fh;
     double res;
     double tot_nr_imgs = 0;
     // Oversample the 1D CTF and Wiener filter vectors OVERSAMPLE times
     // Use 0.55*sqrt(3) to make sure all pixels fit in...
     auto nr_steps= (size_t)ceil(OVERSAMPLE * 0.55 * sqrt((double)(Zdim*Zdim + Ydim*Ydim + Xdim*Xdim)));
 
     Vctfs1D.clear();
     Vwien1D.clear();
     int ii = 0;
     FileName fnVol, fnCTF;
     for (size_t objId : ctfdat.ids())
     {
         // Calculate 1D CTF
         ctfdat.getValue(MDL_IMAGE,fnVol,objId);
         ctfdat.getValue(MDL_CTF_MODEL,fnCTF,objId);
         generateCTF1D(fnCTF,nr_steps,CTF1D);
         Vctfs1D.push_back(CTF1D);
 
         // Get the number of images contributing to this group
         size_t NumberOfImages;
         ctfdat.getValue(MDL_CLASS_COUNT,NumberOfImages,objId);
         tot_nr_imgs += NumberOfImages;
 
         // Calculate denominator of the Wiener filter
         if (ii==0)
             sumterm.resizeNoCopy(CTF1D);
         FOR_ALL_ELEMENTS_IN_ARRAY1D(sumterm)
         {
             double CTF1Di=A1D_ELEM(CTF1D,i);
             A1D_ELEM(sumterm,i) += NumberOfImages * CTF1Di * CTF1Di;
         }
         ii++;
     }
 
     // Add (normalized) Wiener filter constant
     sumterm += tot_nr_imgs*wienerConstant;
 
     double maxwien=0, minwien=1e38;
     ii=0;
     for (size_t objId : ctfdat.ids())
     {
         size_t NumberOfImages;
         ctfdat.getValue(MDL_CLASS_COUNT,NumberOfImages,objId);
         CTF1D=Vctfs1D[ii];
         CTF1D*=NumberOfImages;
         CTF1D/=sumterm;
         Vwien1D.push_back(CTF1D);
 
         // Write CTF and Wiener filter curves to disc
         if (verbose>2)
         {
             FileName fn_tmp;
             fn_tmp = fnRoot + "_wien";
             fn_tmp.compose(fn_tmp, ii+1, "txt");
             fh.open(fn_tmp.c_str(), std::ios::out);
             if (!fh)
                 REPORT_ERROR(ERR_IO_NOWRITE, fn_tmp);
             for (size_t step = 0; step < nr_steps; step++)
             {
                 res = (step * sqrt(3.) ) /
                       (OVERSAMPLE * sqrt( (double) (Zdim*Zdim + Ydim*Ydim + Xdim*Xdim) ) );
                 if (res<=0.5)
                     fh << res << " " << Vwien1D[ii](step) << " " << Vctfs1D[ii](step) << "\n";
             }
             fh.close();
         }
         ii++;
     }
 
     // Some output to screen
     std::cerr <<" ---------------------------------------------------"<<std::endl;
     std::cerr <<" + Number of defocus groups      = "<<ctfdat.size()<<std::endl;
     std::cerr <<" + Total number of images        = "<<tot_nr_imgs<<std::endl;
     std::cerr <<" + Normalized Wiener constant    = "<<tot_nr_imgs*wienerConstant<<std::endl;
     std::cerr <<" + Minimum Wiener filter value   = "<<minwien<<std::endl;
     std::cerr <<" + Maximum Wiener filter value   = "<<maxwien<<std::endl;
     std::cerr <<" ---------------------------------------------------"<<std::endl;
 }

◆ produceSideInfo()

void ProgCtfCorrectAmplitude3D::produceSideInfo ( )

Produce side information. The CTFdat, nr_imgs docfile and selection file with envelopes are read.

Definition at line 77 of file ctf_correct_wiener3d.cpp.

 {
     // Read the CTFdat
     ctfdat.read(fnIn);
 
     // Get dimensions of the volumes
     size_t id = ctfdat.firstRowId();
     FileName fnVol, fnCTF;
     ctfdat.getValue(MDL_IMAGE,fnVol, id);
     ctfdat.getValue(MDL_CTF_MODEL,fnCTF, id);
     Image<double> V;
     V.read(fnVol, HEADER);
     size_t Ndim;
     V.getDimensions(Xdim,Ydim,Zdim,Ndim);
 }

◆ readParams()

void ProgCtfCorrectAmplitude3D::readParams ( )

virtual

Read parameters

Reimplemented from XmippProgram.

Definition at line 31 of file ctf_correct_wiener3d.cpp.

 {
     fnIn  = getParam("-i");
     fnRoot = getParam("--oroot");
     minFreq  = getDoubleParam("--minFreq");
     wienerConstant = getDoubleParam("--wienerConstant");
     isFlipped = checkParam("--phase_flipped");
 }

◆ run()

void ProgCtfCorrectAmplitude3D::run ( )

virtual

Do the whole thing ...

Reimplemented from XmippProgram.

Definition at line 277 of file ctf_correct_wiener3d.cpp.

 {
     produceSideInfo();
     generateWienerFilters();
     generateVolumes();
 }

◆ show()

void ProgCtfCorrectAmplitude3D::show ( )

Show.

Definition at line 41 of file ctf_correct_wiener3d.cpp.

 {
     if (verbose==0)
         return;
     std::cout
     << "Input file:      " << fnIn << std::endl
     << "Output rootname: " << fnRoot << std::endl
     << "Wiener constant: " << wienerConstant << std::endl
     << "Phase flipped:   " << isFlipped << std::endl
     ;
     if (minFreq>0)
         std::cout << "Apply Wiener filter only beyond:   " << minFreq << std::endl;
 }

Member Data Documentation

◆ ctf

FourierFilter ProgCtfCorrectAmplitude3D::ctf

Side Info: CTF.

Definition at line 55 of file ctf_correct_wiener3d.h.

◆ ctfdat

MetaDataVec ProgCtfCorrectAmplitude3D::ctfdat

Side Info: ctfdat.

Definition at line 58 of file ctf_correct_wiener3d.h.

◆ fnIn

FileName ProgCtfCorrectAmplitude3D::fnIn

Metadata with volume, ctf and number of images in that volume.

Definition at line 38 of file ctf_correct_wiener3d.h.

◆ fnRoot

FileName ProgCtfCorrectAmplitude3D::fnRoot

Rootname for output files.

Definition at line 40 of file ctf_correct_wiener3d.h.

◆ isFlipped

bool ProgCtfCorrectAmplitude3D::isFlipped

Flag for phase flipped images.

Definition at line 49 of file ctf_correct_wiener3d.h.

◆ minFreq

double ProgCtfCorrectAmplitude3D::minFreq

Low resolution cutoff to apply Wiener filter.

Definition at line 46 of file ctf_correct_wiener3d.h.

◆ Vctfs1D

std::vector< MultidimArray<double> > ProgCtfCorrectAmplitude3D::Vctfs1D

The 3D CTFs and Wiener filters.

Definition at line 61 of file ctf_correct_wiener3d.h.

◆ Vwien1D

std::vector< MultidimArray<double> > ProgCtfCorrectAmplitude3D::Vwien1D

Definition at line 61 of file ctf_correct_wiener3d.h.

◆ wienerConstant

double ProgCtfCorrectAmplitude3D::wienerConstant

Wiener filter constant.

Definition at line 43 of file ctf_correct_wiener3d.h.

◆ Xdim

size_t ProgCtfCorrectAmplitude3D::Xdim

Definition at line 52 of file ctf_correct_wiener3d.h.

◆ Ydim

size_t ProgCtfCorrectAmplitude3D::Ydim

Definition at line 52 of file ctf_correct_wiener3d.h.

◆ Zdim

size_t ProgCtfCorrectAmplitude3D::Zdim

Dimensions of the volumes.

Definition at line 52 of file ctf_correct_wiener3d.h.

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

xmipp/libraries/reconstruction/ctf_correct_wiener3d.h
xmipp/libraries/reconstruction/ctf_correct_wiener3d.cpp

Public Member Functions

Public Attributes

Additional Inherited Members

Detailed Description

Member Function Documentation

◆ defineParams()

◆ generateCTF1D()

◆ generateVolumes()

◆ generateWienerFilters()

◆ produceSideInfo()

◆ readParams()

◆ run()

◆ show()

Member Data Documentation

◆ ctf

◆ ctfdat

◆ fnIn

◆ fnRoot

◆ isFlipped

◆ minFreq

◆ Vctfs1D

◆ Vwien1D

◆ wienerConstant

◆ Xdim

◆ Ydim

◆ Zdim