EMAN::EMAN2Ctf Class Reference

EMAN2Ctf is the default CTF model used in EMAN2. More...

#include <ctf.h>

Inheritance diagram for EMAN::EMAN2Ctf:

Collaboration diagram for EMAN::EMAN2Ctf:

Public Member Functions
vector< float >	get_snr ()
void	set_snr (const vector< float > &vf)
vector< float >	get_background ()
void	set_background (const vector< float > &vf)
	EMAN2Ctf ()
	EMAN2Ctf (const vector< float > &vf)
	~EMAN2Ctf ()
vector< float >	compute_1d (int size, float ds, CtfType type, XYData *struct_factor=0)
vector< float >	compute_1d_fromimage (int size, float ds, EMData *image)
void	compute_2d_real (EMData *image, CtfType type, XYData *struct_factor=0)
void	compute_2d_complex (EMData *image, CtfType type, XYData *struct_factor=0)
int	from_string (const string &ctf)
string	to_string () const
void	from_dict (const Dict &dict)
Dict	to_dict () const
void	from_vector (const vector< float > &vctf)
vector< float >	to_vector () const
void	copy_from (const Ctf *new_ctf)
bool	equal (const Ctf *ctf1) const
float	zero (int n) const
float	stos2 (float s, float dZ)
float	get_phase () const
void	set_phase (float phase)
Public Member Functions inherited from EMAN::Ctf
virtual	~Ctf ()

Public Attributes
float	dfdiff
float	dfang
float	ampcont
float	dsbg
vector< float >	background
vector< float >	snr
Public Attributes inherited from EMAN::Ctf
float	defocus
float	bfactor
float	voltage
float	cs
float	apix

Private Member Functions
float	lambda () const
float	df (float ang) const
float	calc_noise (float s) const

Additional Inherited Members
Public Types inherited from EMAN::Ctf
enum	CtfType {   CTF_AMP , CTF_SIGN , CTF_BACKGROUND , CTF_SNR ,   CTF_SNR_SMOOTH , CTF_WIENER_FILTER , CTF_TOTAL , CTF_FITREF ,   CTF_NOISERATIO , CTF_INTEN , CTF_POWEVAL , CTF_ALIFILT ,   CTF_ABS }

Detailed Description

EMAN2Ctf is the default CTF model used in EMAN2.

Definition at line 236 of file ctf.h.

Constructor & Destructor Documentation

EMAN2Ctf() [1/2]

EMAN2Ctf::EMAN2Ctf

(

)

Definition at line 488 of file ctf.cpp.

{
        defocus = 0;
        dfdiff = 0;
        dfang = 0;
        bfactor = 0;
        ampcont = 0;
        voltage = 0;
        cs = 0;
        apix = 1.0;
        dsbg=-1;
        background.clear();
        snr.clear();
}

References ampcont, EMAN::Ctf::apix, background, EMAN::Ctf::bfactor, EMAN::Ctf::cs, EMAN::Ctf::defocus, dfang, dfdiff, dsbg, snr, and EMAN::Ctf::voltage.

EMAN2Ctf() [2/2]

EMAN::EMAN2Ctf::EMAN2Ctf ( const vector< float > &  vf )

inline

Definition at line 258 of file ctf.h.

{from_vector(vf);}      //for unpickling

References from_vector().

~EMAN2Ctf()

EMAN2Ctf::~EMAN2Ctf

(

)

Definition at line 504 of file ctf.cpp.

{
}

Member Function Documentation

calc_noise()

float EMAN::EMAN2Ctf::calc_noise ( float  s ) const

inlineprivate

Definition at line 305 of file ctf.h.

                {
                        int si=(int)(s/dsbg);
                        if (si>(int)background.size()||si<0) return background.back();
                        return background[si];
                }

References background, and dsbg.

Referenced by compute_1d(), and compute_2d_complex().

compute_1d()

vector< float > EMAN2Ctf::compute_1d ( int  size, float  ds, CtfType  type, XYData *  struct_factor = 0  )

virtual

Implements EMAN::Ctf.

Definition at line 731 of file ctf.cpp.

{
        Assert(size > 0);
 
//      float tmp_f1 =  sqrt((float) 2) * size / 2;
//      int np = (int) ceil(tmp_f1) + 2;
        int np=size/2;
        vector < float >r;
 
        r.resize(np);
 
        float s = 0;
        float g1=M_PI/2.0*cs*1.0e7*pow(lambda(),3.0f);  // s^4 coefficient for gamma, cached in a variable for simplicity (maybe speed? depends on the compiler)
        float g2=M_PI*lambda()*defocus*10000.0;                 // s^2 coefficient for gamma 
//      float acac=acos(ampcont/100.0);                                 // instead of ac*cos(g)+sqrt(1-ac^2)*sin(g), we can use cos(g-acos(ac)) and save a trig op
        float acac=M_PI/2.0-get_phase();
        switch (type) {
        case CTF_AMP:
                for (int i = 0; i < np; i++) {
                        float gam=-g1*s*s*s*s+g2*s*s;
                        r[i] = cos(gam-acac)*exp(-(bfactor/4.0f * s * s));
                        s += ds;
                }
                break;
 
        case CTF_ABS:
                for (int i = 0; i < np; i++) {
                        float gam=-g1*s*s*s*s+g2*s*s;
                        r[i] = fabs(cos(gam-acac)*exp(-(bfactor/4.0f * s * s)));
                        s += ds;
                }
                break;
                
        case CTF_INTEN:
                for (int i = 0; i < np; i++) {
                        float gam=-g1*s*s*s*s+g2*s*s;
                        r[i] = cos(gam-acac)*exp(-(bfactor/4.0f * s*s));
                        r[i]*=r[i];
                        s += ds;
                }
                break;
                
        case CTF_SIGN:
                for (int i = 0; i < np; i++) {
                        float gam=-g1*s*s*s*s+g2*s*s;
                        r[i] = cos(gam-acac)<0?-1.0:1.0;
                        s += ds;
                }
                break;
 
        case CTF_BACKGROUND:
                for (int i = 0; i < np; i++) {
                        float f = s/dsbg;
                        int j = (int)floor(f);
                        f-=j;
                        if (j>(int)background.size()-2) r[i]=background.back();
                        else r[i]=background[j]*(1.0f-f)+background[j+1]*f;
                        s+=ds;
                }
                break;
 
        case CTF_SNR:
                if (snr.size()<1) {
                        throw InvalidValueException("CTF_SNR"," ERROR: No SNR info in CTF header.");
                        break;
                }
                for (int i = 0; i < np; i++) {
                        float f = s/dsbg;
                        int j = (int)floor(f);
                        f-=j;
                        if (j>(int)snr.size()-2) r[i]=snr.back();
                        else r[i]=snr[j]*(1.0f-f)+snr[j+1]*f;
//                      printf("%d\t%f\n",j,snr[j]);
                        s+=ds;
                }
                break;
        case CTF_SNR_SMOOTH:
                // This apparently complicated routine tries to make a nice smooth and accurate SNR curve. It does this
                // by fitting local regions of the SNR vs the theoretical SNR (theoretical CTF^2/measured background),
                // then taking the slope of the result times the theoretical SNR to produce a local SNR estimate
 
                { // <- is to permit new temporary value allocation
                        vector < float >tsnr;   // theoretical SNR
                        tsnr.resize(np);
                        vector < float >dsnr;   // data SNR
                        dsnr.resize(np);
                        
                        float s0=s;
                        
                        for (int i = 0; i < np; i++) {
                                float gam=-g1*s*s*s*s+g2*s*s;
                                tsnr[i] = cos(gam-acac)*exp(-(bfactor/4.0f * s * s));           // ctf amp
 
                                // background value
                                float f = s/dsbg;
                                int j = (int)floor(f);
                                f-=j;
                                float bg;
                                if (j>(int)background.size()-2) bg=background.back();
                                else bg=background[j]*(1.0f-f)+background[j+1]*f;
                                if (bg <=0) bg=.001f;
 
                                tsnr[i] = tsnr[i]*tsnr[i]/bg;           // This is now a SNR curve
                                if (sf && s) {
                                        tsnr[i] *= sf->get_yatx(s);
                                }
 
                                
                                // This is the SNR computed from the data without fitting
                                if (j>(int)snr.size()-2) dsnr[i]=snr.back();
                                else dsnr[i]=snr[j]*(1.0f-f)+snr[j+1]*f;
                                
                                s+=ds;
                        }
 
                        int npsm=np/25;                 // 1/2 number of points to smooth over, 25 is arbitrary
                        if (npsm<2) npsm=2;
                        
                        s=s0;
                        for (int i = 1; i < np; i++) {
                                // simple linear regression embedded here
                                double sum = 0;
                                double sum_x = 0;
                                double sum_y = 0;
                                double sum_xx = 0;
                                double sum_xy = 0;
 
                                for (int k=max_int(i-npsm,1); k<=min_int(i+npsm,np-1); k++) {
                                        double y = dsnr[k];
                                        double x = tsnr[k];
 
                                        sum_x += x;
                                        sum_y += y;
                                        sum_xx += x * x;
                                        sum_xy += x * y;
                                        sum++;
                                }
 
                                double div = sum * sum_xx - sum_x * sum_x;
//                              if (div == 0) {
//                                      div = 0.0000001f;
//                              }
 
        //                      *intercept = (float) ((sum_xx * sum_y - sum_x * sum_xy) / div);
        //                      *slope = (float) ((sum * sum_xy - sum_x * sum_y) / div);
 
                                if (div!=0.0) r[i]=(float) ((sum * sum_xy - sum_x * sum_y) / div)*tsnr[i];
                                else r[i]=0.0;
                                if (r[i]<0) r[i]=0;
                                
                                s+=ds;
                        }
                        r[0]=0;
                }
                break;
 
        case CTF_WIENER_FILTER:
//              if (!sf) {
//                      LOGERR("CTF computation error, no SF found\n");
//                      return r;
//              }
 
                for (int i = 0; i < np; i++) {
                        float f = s/dsbg;
                        int j = (int)floor(f);
                        float bg;
                        f-=j;
                        if (j>(int)snr.size()-2) {
/*                              r[i]=snr.back();
                                bg=background.back();*/
                                r[i]=0;
                        }
                        else {
                                r[i]=snr[j]*(1.0f-f)+snr[j+1]*f;
                                bg=background[j]*(1.0f-f)+background[j+1]*f;
                        }
                        if (r[i]<0) r[i]=0;
                        r[i]=r[i]/(r[i]+1.0f);          // Full Wiener filter assuming perfect image with noise
//                      r[i]=sqrt(r[i]/bg)/(r[i]+1.0);  // Wiener filter with 1/CTF term (sort of) to restore image then filter
                        s+=ds;
                }
                r[0]=0;
                break;
 
        case CTF_TOTAL:
 
                for (int i = 0; i < np; i++) {
                        float gam=-g1*s*s*s*s+g2*s*s;
                        if (sf) {
                                r[i] = cos(gam-acac)*exp(-(bfactor/4.0f * s*s));
                                r[i] = r[i] * r[i] * sf->get_yatx(s) + calc_noise(s);
                        }
                        else {
                                r[i] = cos(gam-acac)*exp(-(bfactor/4.0f * s*s));
                                r[i] = r[i] * r[i] + calc_noise(s);
                        }
                        s += ds;
                }
                break;
        case CTF_NOISERATIO:
                for (int i = 0; i < np; i++) {
                        float gam=-g1*s*s*s*s+g2*s*s;
                        if (sf) {
                                r[i] = cos(gam-acac)*exp(-(bfactor/4.0f * s*s));
                                r[i] = r[i] * r[i] * sf->get_yatx(s);
                                r[i] = r[i]/(r[i]+calc_noise(s));
                        }
                        else {
                                r[i] = cos(gam-acac)*exp(-(bfactor/4.0f * s*s));
                                r[i] = r[i] * r[i];
                                r[i] = r[i]/(r[i]+calc_noise(s));
                        }
                        
                        s+=ds;
                }
                        
                break;
        default:
                break;
        }
 
        return r;
}

References Assert, background, EMAN::Ctf::bfactor, calc_noise(), EMAN::Ctf::cs, EMAN::Ctf::CTF_ABS, EMAN::Ctf::CTF_AMP, EMAN::Ctf::CTF_BACKGROUND, EMAN::Ctf::CTF_INTEN, EMAN::Ctf::CTF_NOISERATIO, EMAN::Ctf::CTF_SIGN, EMAN::Ctf::CTF_SNR, EMAN::Ctf::CTF_SNR_SMOOTH, EMAN::Ctf::CTF_TOTAL, EMAN::Ctf::CTF_WIENER_FILTER, EMAN::Ctf::defocus, div(), dsbg, get_phase(), EMAN::XYData::get_yatx(), InvalidValueException, lambda(), max_int(), min_int(), snr, x, and y.

Referenced by EMAN::CtfSimProcessor::process(), and EMAN::CTFCorrProcessor::process_inplace().

compute_1d_fromimage()

vector< float > EMAN2Ctf::compute_1d_fromimage ( int  size, float  ds, EMData *  image  )

virtual

Implements EMAN::Ctf.

Definition at line 685 of file ctf.cpp.

{
        size/=2;        // to be consistent with the next routine
        vector <float> ret(size);
        vector <float> norm(size);
        
        if (!image->is_complex() || !image->is_ri()) ImageFormatException("EMAN2Ctf::compute_1d_fromimage requires a complex, ri image");
        int isinten=image->get_attr_default("is_intensity",0);
        if (isinten) ImageFormatException("EMAN2Ctf::compute_1d_fromimage does not take intensity images");
        
        int nx=image->get_attr("nx");
        int ny=image->get_attr("ny");
        float *data=image->get_data();
        
        int x,y,i;
        for (y=i=0; y<ny; y++) {
                for (x=0; x<nx; x+=2,i+=2) {
                        if (x==0 && y>ny/2) continue;
                        float r=(float)(Util::hypot_fast(x/2,y<ny/2?y:ny-y));           // origin at 0,0; periodic
                        float a=(float)atan2((float)(y<ny/2?y:ny-y),(float)(x/2));              // angle to point (radians)
                        r=stos2(r*ds,-dfdiff/2*cos(2.0*a-2.0*M_PI/180.0*dfang))/ds;             // angle-based defocus -> convert s to remove astigmatism, dfang in degrees
                        int f=int(r);   // safe truncation, so floor isn't needed
                        r-=float(f);    // r is now the fractional spacing between bins
                        
                        float v=Util::square_sum(data[i],data[i+1]);                            //intensity
                        if (f>=0 && f<size) {
                                ret[f]+=v*(1.0f-r);
                                norm[f]+=(1.0f-r);
                                if (f<size-1) {
                                        ret[f+1]+=v*r;
                                        norm[f+1]+=r;
                                }
                        }
                }
        }
 
        for (i=0; i<size; i++) ret[i]/=norm[i]?norm[i]:1.0f;    // Normalize
 
        
        return ret;
}

References dfang, dfdiff, EMAN::Util::hypot_fast(), ImageFormatException, EMAN::Util::square_sum(), stos2(), x, and y.

compute_2d_complex()

void EMAN2Ctf::compute_2d_complex ( EMData *  image, CtfType  type, XYData *  struct_factor = 0  )

virtual

Implements EMAN::Ctf.

Definition at line 964 of file ctf.cpp.

{
        if (!image) {
                LOGERR("image is null. cannot computer 2D complex CTF");
                return;
        }
 
        if (image->is_complex() == false) {
                LOGERR("compute_2d_complex can only work on complex images");
                return;
        }
 
        int nx = image->get_xsize();
        int ny = image->get_ysize();
 
        if ((ny%2==1 && nx!=ny+1) || (ny%2==0 && nx != ny + 2)) {
                printf("nx=%d  ny=%d\n",nx,ny);
                LOGERR("compute_2d_complex only works on (nx, nx-2) images");
                return;
        }
 
        float ds = 1.0f / (apix * ny);
        image->to_one();
 
        float *d = image->get_data();
        float g1=M_PI/2.0*cs*1.0e7*pow(lambda(),3.0f);  // s^4 coefficient for gamma, cached in a variable for simplicity (maybe speed? depends on the compiler)
        float g2=M_PI*lambda()*10000.0;                                 // s^2 coefficient for gamma 
//      float acac=acos(ampcont/100.0);                                 // instead of ac*cos(g)+sqrt(1-ac^2)*sin(g), we can use cos(g-acos(ac)) and save a trig op
        float acac=M_PI/2.0-get_phase();
 
        if (type == CTF_BACKGROUND) {
                for (int y = -ny/2; y < ny/2; y++) {
                        int y2=(y+ny)%ny;
                        int ynx = y2 * nx;
 
                        for (int x = 0; x < nx / 2; x++) {
                                float s = (float) Util::hypot_fast(x, y ) * ds;
                                d[x * 2 + ynx] = calc_noise(s);
                                d[x * 2 + ynx + 1] = 0;                 // The phase is somewhat arbitrary
                        }
                }
        }
        else if (type == CTF_AMP) {
                for (int y = -ny/2; y < ny/2; y++) {
                        int y2=(y+ny)%ny;
                        int ynx = y2 * nx;
 
                        for (int x = 0; x < nx / 2; x++) {
                                float s = (float)Util::hypot_fast(x,y ) * ds;
                                float gam;
                                if (dfdiff==0) gam=-g1*s*s*s*s+g2*defocus*s*s;
                                else gam=-g1*s*s*s*s+g2*df(atan2((float)y,(float)x))*s*s;
                                float v = cos(gam-acac)*exp(-(bfactor/4.0f * s*s));
                                d[x * 2 + ynx] = v;
                                d[x * 2 + ynx + 1] = 0;
                        }
                }
        }
        else if (type == CTF_ABS) {
                for (int y = -ny/2; y < ny/2; y++) {
                        int y2=(y+ny)%ny;
                        int ynx = y2 * nx;
 
                        for (int x = 0; x < nx / 2; x++) {
                                float s = (float)Util::hypot_fast(x,y ) * ds;
                                float gam;
                                if (dfdiff==0) gam=-g1*s*s*s*s+g2*defocus*s*s;
                                else gam=-g1*s*s*s*s+g2*df(atan2((float)y,(float)x))*s*s;
                                float v = cos(gam-acac)*exp(-(bfactor/4.0f * s*s));
                                d[x * 2 + ynx] = fabs(v);
                                d[x * 2 + ynx + 1] = 0;
                        }
                }
        }
        else if (type == CTF_INTEN) {
                for (int y = -ny/2; y < ny/2; y++) {
                        int y2=(y+ny)%ny;
                        int ynx = y2 * nx;
 
                        for (int x = 0; x < nx / 2; x++) {
                                float s = (float)Util::hypot_fast(x,y ) * ds;
                                float gam;
                                if (dfdiff==0) gam=-g1*s*s*s*s+g2*defocus*s*s;
                                else gam=-g1*s*s*s*s+g2*df(atan2((float)y,(float)x))*s*s;
                                float v = cos(gam-acac)*exp(-(bfactor/4.0f * s*s));
                                d[x * 2 + ynx] = v*v;
                                d[x * 2 + ynx + 1] = 0;
                        }
                }
        }
        else if (type == CTF_POWEVAL) {
                for (int y = -ny/2; y < ny/2; y++) {
                        int y2=(y+ny)%ny;
                        int ynx = y2 * nx;
 
                        for (int x = 0; x < nx / 2; x++) {
                                float s = (float)Util::hypot_fast(x,y ) * ds;
                                float gam;
                                // mf is used to gradually "turn on" the curve as we approach 20 A
                                float mf=1.0;
                                if (s<.04) mf=0.0;
                                else if (s<.05) mf=1.0-exp(-pow((s-.04f)*300.0f,2.0f));
                                if (dfdiff==0) gam=-g1*s*s*s*s+g2*defocus*s*s;
                                else gam=-g1*s*s*s*s+g2*df(atan2((float)y,(float)x))*s*s;
                                float v = cos(gam-acac);
                                if (v>0.9) v=exp(-(50.0f/4.0f * s*s));
                                else v=0;
                                d[x * 2 + ynx] = mf*v*v;
                                d[x * 2 + ynx + 1] = 0;
                        }
                }
        }
        else if (type == CTF_SIGN) {
                for (int y = -ny/2; y < ny/2; y++) {
                        int y2=(y+ny)%ny;
                        int ynx = y2 * nx;
                        for (int x = 0; x < nx / 2; x++) {
                                float s = (float)Util::hypot_fast(x,y ) * ds;
                                float gam;
                                if (dfdiff==0) gam=-g1*s*s*s*s+g2*defocus*s*s;
                                else gam=-g1*s*s*s*s+g2*df(atan2((float)y,(float)x))*s*s;
                                float v = cos(gam-acac);
                                d[x * 2 + ynx] = v<0?-1.0:1.0;
                                d[x * 2 + ynx + 1] = 0;
                        }
                }
        }
        else if (type == CTF_FITREF) {
                for (int y = -ny/2; y < ny/2; y++) {
                        int y2=(y+ny)%ny;
                        int ynx = y2 * nx;
 
                        for (int x = 0; x < nx / 2; x++) {
                                float s = (float)Util::hypot_fast(x,y ) * ds;
                                // We exclude very low frequencies from consideration due to the strong structure factor
                                if (s<.04) {
                                        d[x * 2 + ynx] = 0;
                                        d[x * 2 + ynx + 1] = 0;
                                        continue;
                                }
                                // Rather than suddenly "turning on" the CTF, we do it gradually
                                float mf=1.0;
                                if (s<.05) {
                                        mf=1.0-exp(-pow((s-.04f)*300.0f,2.0f));
                                }
                                float gam;
                                if (dfdiff==0) gam=-g1*s*s*s*s+g2*defocus*s*s;
                                else gam=-g1*s*s*s*s+g2*df(atan2((float)y,(float)x))*s*s;
                                float v = cos(gam-acac)*exp(-(bfactor/4.0f * s*s));
                                d[x * 2 + ynx] = mf*v*v;
                                d[x * 2 + ynx + 1] = 0;
                        }
                }
        }
        else if (type == CTF_SNR) {
 
                for (int y = -ny/2; y < ny/2; y++) {
                        int y2=(y+ny)%ny;
                        int ynx = y2 * nx;
 
                        for (int x = 0; x < nx / 2; x++) {
 
                                float s = (float)Util::hypot_fast(x,y ) * ds;
                                float f = s/dsbg;
                                int j = (int)floor(f);
                                f-=j;
                                if (j>(int)snr.size()-2) d[x*2+ynx]=snr.back();
                                else d[x*2+ynx]=snr[j]*(1.0f-f)+snr[j+1]*f;
                                d[x * 2 + ynx + 1] = 0;
                        }
                }
                d[0]=0;
        }
        else if (type == CTF_SNR_SMOOTH) {
                for (int y = -ny/2; y < ny/2; y++) {
                        int y2=(y+ny)%ny;
                        int ynx = y2 * nx;
 
                        for (int x = 0; x < nx / 2; x++) {
 
                                float s = (float)Util::hypot_fast(x,y ) * ds;
                                float f = s/dsbg;
                                int j = (int)floor(f);
                                f-=j;
                                if (j>(int)snr.size()-2) d[x*2+ynx]=snr.back();
                                else d[x*2+ynx]=snr[j]*(1.0f-f)+snr[j+1]*f;
                                d[x * 2 + ynx + 1] = 0;
                        }
                }
                d[0]=0;
        }
        else if (type == CTF_WIENER_FILTER) {
                if (dsbg==0) printf("Warning, DSBG set to 0\n");
                for (int y = -ny/2; y < ny/2; y++) {
                        int y2=(y+ny)%ny;
                        int ynx = y2 * nx;
 
                        for (int x = 0; x < nx / 2; x++) {
 
                                float s = (float)Util::hypot_fast(x,y ) * ds;
                                float f = s/dsbg;
                                int j = (int)floor(f);
                                float bg,snrf;
                                f-=j;
                                if (j>(int)snr.size()-2) {
/*                                      bg=background.back();
                                        d[x*2+ynx]=snr.back()/(snr.back()+1.0);*/
                                        d[x*2+ynx]=0;
                                }
                                else {
                                        bg=background[j]*(1.0f-f)+background[j+1]*f;
                                        snrf=snr[j]*(1.0f-f)+snr[j+1]*f;
 
//                                      printf("%d\t%f\n",x,sqrt(snrf/bg)/(snrf+1.0));
                                        if (snrf<0) snrf=0.0;
//                                      d[x*2+ynx]=sqrt(snrf/bg)/(snrf+1.0);    // Note that this is a Wiener filter with a 1/CTF term to compensate for the filtration already applied, but needs to be multiplied by the structure factor
                                        d[x*2+ynx]=snrf/(snrf+1);       // This is just the simple Wiener filter
 
                                }
                                d[x * 2 + ynx + 1] = 0;
                        }
                }
                d[0]=0;
        }
        else if (type == CTF_TOTAL) {
 
                for (int y = -ny/2; y < ny/2; y++) {
                        int y2=(y+ny)%ny;
                        int ynx = y2 * nx;
 
                        for (int x = 0; x < nx / 2; x++) {
                                float s = (float)Util::hypot_fast(x,y ) * ds;
                                float gam;
                                if (dfdiff==0) gam=-g1*s*s*s*s+g2*defocus*s*s;
                                else gam=-g1*s*s*s*s+g2*df(atan2((float)y,(float)x))*s*s;
                                float v = cos(gam-acac)*exp(-(bfactor/4.0f * s*s));
                                d[x * 2 + ynx] = v*v+calc_noise(s);
                                d[x * 2 + ynx + 1] = 0;
//                              if ((x==101||x==100) && y2==79) printf("%f %f %f %f\n",s,gam,v,d[x * 2 + ynx]);
                        }
                }
        }
        else if (type == CTF_ALIFILT) {
                double vt=0;
                for (int y = -ny/2; y < ny/2; y++) {
                        int y2=(y+ny)%ny;
                        int ynx = y2 * nx;
 
                        for (int x = 0; x < nx / 2; x++) {
                                float s = (float)Util::hypot_fast(x,y ) * ds;
                                float gam;
                                if (dfdiff==0) gam=-g1*s*s*s*s+g2*defocus*s*s;
                                else gam=-g1*s*s*s*s+g2*df(atan2((float)y,(float)x))*s*s;
                                // this make values near the zero crossings negative to cancel out 0,0 correlation
//                              float v = (pow(cos(gam-acac),2.0)-0.5)*exp(-(bfactor/4.0f * s*s));
                                // Basically just a CTF weight
                                float v = (pow((float)cos(gam-acac),2.0f))*exp(-(bfactor/4.0f * s*s));
                                //v*=v;
                                //v=fabs(v);
                                d[x * 2 + ynx] = v;
                                d[x * 2 + ynx + 1] = 0;
                                vt+=v;
                        }
                }
                vt/=nx*ny/2;
//              for (size_t i=0; i<image->get_xsize()*image->get_ysize(); i+=2) d[i]=(d[i]-vt)<0?-1.0:1.0;
//              for (size_t i=0; i<image->get_xsize()*image->get_ysize(); i+=2) d[i]=(d[i]-vt);
        }
        else printf("Unknown CTF image mode\n");
 
        image->update();
}

References EMAN::Ctf::apix, background, EMAN::Ctf::bfactor, calc_noise(), EMAN::Ctf::cs, EMAN::Ctf::CTF_ABS, EMAN::Ctf::CTF_ALIFILT, EMAN::Ctf::CTF_AMP, EMAN::Ctf::CTF_BACKGROUND, EMAN::Ctf::CTF_FITREF, EMAN::Ctf::CTF_INTEN, EMAN::Ctf::CTF_POWEVAL, EMAN::Ctf::CTF_SIGN, EMAN::Ctf::CTF_SNR, EMAN::Ctf::CTF_SNR_SMOOTH, EMAN::Ctf::CTF_TOTAL, EMAN::Ctf::CTF_WIENER_FILTER, EMAN::Ctf::defocus, df(), dfdiff, dsbg, get_phase(), EMAN::Util::hypot_fast(), lambda(), LOGERR, snr, x, and y.

compute_2d_real()

void EMAN2Ctf::compute_2d_real ( EMData *  image, CtfType  type, XYData *  struct_factor = 0  )

virtual

Implements EMAN::Ctf.

Definition at line 956 of file ctf.cpp.

{
 
 
}

copy_from()

void EMAN2Ctf::copy_from ( const Ctf *  new_ctf )

virtual

Implements EMAN::Ctf.

Definition at line 653 of file ctf.cpp.

{
        if (new_ctf) {
                const EMAN2Ctf *c = static_cast<const EMAN2Ctf *>(new_ctf);
                defocus = c->defocus;
                dfdiff = c->dfdiff;
                dfang = c->dfang;
                bfactor = c->bfactor;
                ampcont = c->ampcont;
                voltage = c->voltage;
                cs = c->cs;
                apix = c->apix;
                dsbg = c->dsbg;
                background = c->background;
                snr = c->snr;
        }
}

References ampcont, EMAN::Ctf::apix, background, EMAN::Ctf::bfactor, EMAN::Ctf::cs, EMAN::Ctf::defocus, dfang, dfdiff, dsbg, snr, and EMAN::Ctf::voltage.

Referenced by EMAN::CTFCorrProcessor::process_inplace().

df()

float EMAN::EMAN2Ctf::df ( float  ang ) const

inlineprivate

Definition at line 300 of file ctf.h.

                                                 {
                        if (dfdiff==0.0) return defocus;
                        return defocus+dfdiff/2.0f*cos(2.0f*ang-2.0f*M_PI/180.0f*dfang);
                }

References EMAN::Ctf::defocus, dfang, and dfdiff.

Referenced by compute_2d_complex().

equal()

bool EMAN2Ctf::equal ( const Ctf *  ctf1 ) const

virtual

Implements EMAN::Ctf.

Definition at line 1239 of file ctf.cpp.

{
        if (ctf1) {
                const EMAN2Ctf *c = static_cast<const EMAN2Ctf *>(ctf1);
                if (defocus != c->defocus ||
                        dfdiff != c->dfdiff ||
                        dfang != c->dfang ||
                        bfactor != c->bfactor ||
                        ampcont != c->ampcont ||
                        voltage != c->voltage ||
                        cs != c->cs ||
                        apix != c->apix ||
                        dsbg != c->dsbg ||
                        background.size() != c->background.size() ||
                        snr.size() != c->snr.size()
                        ) return false;
 
                for (unsigned int i=0; i<background.size(); i++) {
                        if (background[i]!=c->background[i]) return false;
                }
                for (unsigned int i=0; i<snr.size(); i++) {
                        if (snr[i]!=c->snr[i]) return false;
                }
                return true;
        }
        return false;
}

References ampcont, EMAN::Ctf::apix, background, EMAN::Ctf::bfactor, EMAN::Ctf::cs, EMAN::Ctf::defocus, dfang, dfdiff, dsbg, snr, and EMAN::Ctf::voltage.

from_dict()

void EMAN2Ctf::from_dict ( const Dict &  dict )

virtual

Implements EMAN::Ctf.

Definition at line 567 of file ctf.cpp.

{
//      printf("st\n"); 
        defocus = (float)dict["defocus"];
//      printf("df\n");
        dfdiff = (float)dict["dfdiff"];
//      printf("dfd\n");
        dfang = (float)dict["dfang"];
//      printf("dfa\n");
        bfactor = (float)dict["bfactor"];
//      printf("bf\n");
        ampcont = (float)dict["ampcont"];
//      printf("ac\n");
        voltage = (float)dict["voltage"];
//      printf("vo\n");
        cs = (float)dict["cs"];
//      printf("cs\n");
        apix = (float)dict["apix"];
//      printf("ap\n");
        dsbg = (float)dict["dsbg"];
//      printf("ds\n");
        background = dict["background"];
//      printf("bg\n");
        snr = dict["snr"];
//      printf("done\n");
}

References ampcont, EMAN::Ctf::apix, background, EMAN::Ctf::bfactor, EMAN::Ctf::cs, EMAN::Ctf::defocus, dfang, dfdiff, dsbg, snr, and EMAN::Ctf::voltage.

from_string()

int EMAN2Ctf::from_string ( const string &  ctf )

virtual

Implements EMAN::Ctf.

Definition at line 509 of file ctf.cpp.

{
        Assert(ctf != "");
        char type=' ';
        int pos=-1,i,j;
        int bglen=0,snrlen=0;
        float v;
        const char *s=ctf.c_str();
 
        sscanf(s, "%c%f %f %f %f %f %f %f %f %f %d%n",
                                   &type,&defocus, &dfdiff,&dfang,&bfactor,&ampcont,&voltage, &cs, &apix,&dsbg,&bglen,&pos);
        if (type!='E') throw InvalidValueException(type,"Trying to initialize Ctf object with bad string");
        if (pos==-1) throw InvalidValueException(s," Invalid CTF string");
        
 
        background.resize(bglen);
        for (i=0; i<bglen; i++) {
                if (sscanf(s+pos,",%f%n",&v,&j)<1) return(1);
                background[i]=v;
                pos+=j;
        }
 
        sscanf(s+pos," %d%n",&snrlen,&j);
        pos+=j;
        snr.resize(snrlen);
        for (i=0; i<snrlen; i++) {
                if (sscanf(s+pos,",%f%n",&v,&j)<1) return(1);
                snr[i]=v;
                pos+=j;
        }
 
        return 0;
 
}

References ampcont, EMAN::Ctf::apix, Assert, background, EMAN::Ctf::bfactor, EMAN::Ctf::cs, EMAN::Ctf::defocus, dfang, dfdiff, dsbg, InvalidValueException, snr, and EMAN::Ctf::voltage.

from_vector()

void EMAN2Ctf::from_vector ( const vector< float > &  vctf )

virtual

Implements EMAN::Ctf.

Definition at line 612 of file ctf.cpp.

{
        int i;
        defocus = vctf[0];
        dfdiff = vctf[1];
        dfang = vctf[2];
        bfactor = vctf[3];
        ampcont = vctf[4];
        voltage = vctf[5];
        cs = vctf[6];
        apix = vctf[7];
        dsbg = vctf[8];
        background.resize((int)vctf[9]);
        for (i=0; i<(int)vctf[9]; i++) background[i]=vctf[i+10];
        snr.resize((int)vctf[i+10]);
        for (int j=0; j<(int)vctf[i+10]; j++) snr[j]=vctf[i+j+11];
}

References ampcont, EMAN::Ctf::apix, background, EMAN::Ctf::bfactor, EMAN::Ctf::cs, EMAN::Ctf::defocus, dfang, dfdiff, dsbg, snr, and EMAN::Ctf::voltage.

Referenced by EMAN2Ctf().

get_background()

vector< float > EMAN::EMAN2Ctf::get_background ( )

inline

Definition at line 253 of file ctf.h.

{ return background;}

References background.

get_phase()

float EMAN2Ctf::get_phase ( ) const

virtual

Implements EMAN::Ctf.

Definition at line 671 of file ctf.cpp.

                                {
        if (ampcont>-100.0 && ampcont<=100.0) return asin(ampcont/100.0);
        if (ampcont>100.0) return M_PI-asin(2.0f-ampcont/100.0);
        return -M_PI-asin(-2.0-ampcont/100.0);
}

References ampcont.

Referenced by compute_1d(), compute_2d_complex(), and zero().

get_snr()

vector< float > EMAN::EMAN2Ctf::get_snr ( )

inline

Definition at line 251 of file ctf.h.

{ return snr;}

References snr.

lambda()

float EMAN::EMAN2Ctf::lambda ( ) const

inlineprivate

Definition at line 293 of file ctf.h.

                {
                        float lambda = 12.2639f / sqrt(voltage * 1000.0f + 0.97845f * voltage * voltage);
                        return lambda;
                }

References lambda(), sqrt(), and EMAN::Ctf::voltage.

Referenced by compute_1d(), compute_2d_complex(), lambda(), stos2(), and zero().

set_background()

void EMAN::EMAN2Ctf::set_background ( const vector< float > &  vf )

inline

Definition at line 254 of file ctf.h.

{background = vf;}

References background.

set_phase()

void EMAN2Ctf::set_phase ( float  phase )

virtual

Implements EMAN::Ctf.

Definition at line 677 of file ctf.cpp.

                                    {
        if (phase>=-M_PI/2.0f && phase<M_PI/2.0f) ampcont=sin(phase)*100.0;
        else if (phase>=M_PI/2.0f && phase<=M_PI) ampcont=(2.0-sin(phase))*100.0;
        else if (phase>-M_PI && phase<-M_PI/2.0f) ampcont=(-2.0-sin(phase))*100.0;
        else printf("Phase = %0.4f deg  out of range\n",phase*180.0/M_PI);
        
}

References ampcont, and phase().

set_snr()

void EMAN::EMAN2Ctf::set_snr ( const vector< float > &  vf )

inline

Definition at line 252 of file ctf.h.

{snr = vf;}

References snr.

stos2()

float EMAN::EMAN2Ctf::stos2 ( float  s, float  dZ  )

inline

Definition at line 280 of file ctf.h.

                                                     {
                        float lmb=lambda();
                        return sqrt((defocus*1.0e4f-sqrt(defocus*defocus*1.0e8f -2.0e11f*cs*s*s*(defocus-dZ)*lmb*lmb+1.0e14f*cs*cs*pow(s*lmb,4.0f)))/(1.0e7f*cs*lmb*lmb));
                }

References EMAN::Ctf::cs, EMAN::Ctf::defocus, lambda(), and sqrt().

Referenced by compute_1d_fromimage().

to_dict()

Dict EMAN2Ctf::to_dict ( ) const

virtual

Implements EMAN::Ctf.

Definition at line 594 of file ctf.cpp.

{
        Dict dict;
        dict["defocus"] = defocus;
        dict["dfdiff"] = dfdiff;
        dict["dfang"] = dfang;
        dict["bfactor"] = bfactor;
        dict["ampcont"] = ampcont;
        dict["voltage"] = voltage;
        dict["cs"] = cs;
        dict["apix"] = apix;
        dict["dsbg"] = dsbg;
        dict["background"] = background;
        dict["snr"] = snr;
 
        return dict;
}

References ampcont, EMAN::Ctf::apix, background, EMAN::Ctf::bfactor, EMAN::Ctf::cs, EMAN::Ctf::defocus, dfang, dfdiff, dsbg, snr, and EMAN::Ctf::voltage.

to_string()

string EMAN2Ctf::to_string ( ) const

virtual

Implements EMAN::Ctf.

Definition at line 544 of file ctf.cpp.

{
        char ctf[256];
        sprintf(ctf, "E%1.10g %1.10g %1.10g %1.10g %1.10g %1.10g %1.10g %1.10g %1.10g %d",
                        defocus, dfdiff, dfang, bfactor, ampcont, voltage, cs, apix, dsbg,(int)background.size());
 
        string ret=ctf;
        for (int i=0; i<(int)background.size(); i++) {
                sprintf(ctf,",%1.3g",background[i]);
                ret+=ctf;
        }
 
        sprintf(ctf, " %d",(int)snr.size());
        ret+=ctf;
        for (int i=0; i<(int)snr.size(); i++) {
                sprintf(ctf,",%1.3g",snr[i]);
                ret+=ctf;
        }
 
 
        return ret;
}

References ampcont, EMAN::Ctf::apix, background, EMAN::Ctf::bfactor, EMAN::Ctf::cs, EMAN::Ctf::defocus, dfang, dfdiff, dsbg, snr, and EMAN::Ctf::voltage.

to_vector()

vector< float > EMAN2Ctf::to_vector ( ) const

virtual

Implements EMAN::Ctf.

Definition at line 630 of file ctf.cpp.

{
        vector<float> vctf;
 
        vctf.push_back(defocus);
        vctf.push_back(dfdiff);
        vctf.push_back(dfang);
        vctf.push_back(bfactor);
        vctf.push_back(ampcont);
        vctf.push_back(voltage);
        vctf.push_back(cs);
        vctf.push_back(apix);
        vctf.push_back(dsbg);
        vctf.push_back((float)background.size());
        for (unsigned int i=0; i<background.size(); i++) vctf.push_back(background[i]);
        vctf.push_back((float)snr.size());
        for (unsigned int j=0; j<snr.size(); j++) vctf.push_back(snr[j]);
 
        return vctf;
}

References ampcont, EMAN::Ctf::apix, background, EMAN::Ctf::bfactor, EMAN::Ctf::cs, EMAN::Ctf::defocus, dfang, dfdiff, dsbg, snr, and EMAN::Ctf::voltage.

zero()

float EMAN2Ctf::zero ( int  n ) const

virtual

Implements EMAN::Ctf.

Definition at line 1267 of file ctf.cpp.

                                {
vector <float>zeroes;
float lam=lambda();
float acac=M_PI/2.0-get_phase();
 
int m=n*2;
if (m<15) m=15;         // A bit arbitrary. I believe this will always get us the zeroes we need in any practical situation
for (int i=-m; i<m; i++) {
        float r1=defocus*defocus*1.0e8 + cs*lam*1.0e7 - 2.0*cs*i*lam*1.0e7 - 2*cs*1.0e7*lam*acac/M_PI;
//      printf("%f\n",r1);
        if (r1<0) continue;
        float r2=defocus*1.0e4+sqrt(r1);
        if (r2>0) zeroes.push_back(sqrt(r2/(cs*lam*lam*1.0e7)));
        float r2a=defocus*1.0e4-sqrt(r1);
        if (r2a>0) zeroes.push_back(sqrt(r2a/(cs*lam*lam*1.0e7)));
}
if (zeroes.size()==0) return 0.0f;
std::sort(zeroes.begin(),zeroes.end());
 
return zeroes[n];
}

References EMAN::Ctf::cs, EMAN::Ctf::defocus, get_phase(), lambda(), and sqrt().

Member Data Documentation

ampcont

float EMAN::EMAN2Ctf::ampcont

Definition at line 243 of file ctf.h.

Referenced by copy_from(), EMAN2Ctf(), equal(), from_dict(), from_string(), from_vector(), get_phase(), EMAN::CtfSimProcessor::process(), EMAN::CTFCorrProcessor::process_inplace(), set_phase(), to_dict(), to_string(), and to_vector().

background

vector<float> EMAN::EMAN2Ctf::background

Definition at line 248 of file ctf.h.

Referenced by calc_noise(), compute_1d(), compute_2d_complex(), copy_from(), EMAN2Ctf(), equal(), from_dict(), from_string(), from_vector(), get_background(), set_background(), to_dict(), to_string(), and to_vector().

dfang

float EMAN::EMAN2Ctf::dfang

Definition at line 241 of file ctf.h.

Referenced by compute_1d_fromimage(), copy_from(), df(), EMAN2Ctf(), equal(), from_dict(), from_string(), from_vector(), to_dict(), to_string(), and to_vector().

dfdiff

float EMAN::EMAN2Ctf::dfdiff

Definition at line 240 of file ctf.h.

Referenced by compute_1d_fromimage(), compute_2d_complex(), copy_from(), df(), EMAN2Ctf(), equal(), from_dict(), from_string(), from_vector(), EMAN::CTFCorrProcessor::process_inplace(), to_dict(), to_string(), and to_vector().

dsbg

float EMAN::EMAN2Ctf::dsbg

Definition at line 247 of file ctf.h.

Referenced by calc_noise(), compute_1d(), compute_2d_complex(), copy_from(), EMAN2Ctf(), equal(), from_dict(), from_string(), from_vector(), EMAN::CtfSimProcessor::process(), EMAN::CTFCorrProcessor::process_inplace(), to_dict(), to_string(), and to_vector().

snr

vector<float> EMAN::EMAN2Ctf::snr

Definition at line 249 of file ctf.h.

Referenced by compute_1d(), compute_2d_complex(), copy_from(), EMAN2Ctf(), equal(), from_dict(), from_string(), from_vector(), get_snr(), set_snr(), to_dict(), to_string(), and to_vector().

---

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

• ctf.h
• ctf.cpp