EMAN2
Public Member Functions | Static Public Member Functions | Static Public Attributes | Private Member Functions | Private Attributes
EMAN::nn4_ctfwReconstructor Class Reference

nn4_ctfw Direct Fourier Weighted Inversion Reconstructor More...

#include <reconstructor.h>

Inheritance diagram for EMAN::nn4_ctfwReconstructor:
Inheritance graph
[legend]
Collaboration diagram for EMAN::nn4_ctfwReconstructor:
Collaboration graph
[legend]

List of all members.

Public Member Functions

 nn4_ctfwReconstructor ()
 nn4_ctfwReconstructor (const string &symmetry, int size, int npad, float snr, int sign)
virtual ~nn4_ctfwReconstructor ()
virtual void setup ()
 Initialize the reconstructor.
virtual int insert_slice (const EMData *const slice, const Transform &euler, const float weight)
 Insert a slice into a 3D volume, in a given orientation.
virtual EMDatafinish (bool doift=true)
 Finish reconstruction and return the complete model.
virtual string get_name () const
 Get the unique name of this class (especially for factory based instantiation access)
virtual string get_desc () const
 Get a clear, concise description of this class.
TypeDict get_param_types () const
void setup (const string &symmetry, int size, int npad, float snr, int sign)
int insert_padfft_slice_weighted (EMData *padfft, EMData *sigmasq2, const Transform &trans, const float weight)

Static Public Member Functions

static ReconstructorNEW ()

Static Public Attributes

static const string NAME = "nn4_ctfw"

Private Member Functions

void buildFFTVolume ()
void buildNormVolume ()

Private Attributes

EMDatam_volume
EMDatam_wptr
EMDatam_refvol
int m_vnx
int m_vny
int m_vnz
int m_vnzp
int m_vnyp
int m_vnxp
int m_vnxc
int m_vnyc
int m_vnzc
int m_npad
int m_sign
int m_varsnr
int m_weighting
float m_wghta
float m_wghtb
float m_snr
string m_symmetry
int m_nsym

Detailed Description

nn4_ctfw Direct Fourier Weighted Inversion Reconstructor

Definition at line 1313 of file reconstructor.h.


Constructor & Destructor Documentation

nn4_ctfwReconstructor::nn4_ctfwReconstructor ( )

Definition at line 3451 of file reconstructor.cpp.

References m_volume, and m_wptr.

Referenced by NEW().

{
        m_volume  = NULL;
        m_wptr    = NULL;
}
nn4_ctfwReconstructor::nn4_ctfwReconstructor ( const string &  symmetry,
int  size,
int  npad,
float  snr,
int  sign 
)

Definition at line 3457 of file reconstructor.cpp.

References setup().

{
        setup( symmetry, size, npad, snr, sign );
}
nn4_ctfwReconstructor::~nn4_ctfwReconstructor ( ) [virtual]

Definition at line 3462 of file reconstructor.cpp.

{
        //if( m_delete_volume ) checked_delete(m_volume);

        //if( m_delete_weight ) checked_delete( m_wptr );

        //checked_delete( m_result );
}

Member Function Documentation

void nn4_ctfwReconstructor::buildFFTVolume ( ) [private]
void nn4_ctfwReconstructor::buildNormVolume ( ) [private]
EMData * nn4_ctfwReconstructor::finish ( bool  doift = true) [virtual]

Finish reconstruction and return the complete model.

Parameters:
doiftA flag indicating whether the returned object should be guaranteed to be in real-space (true) or should be left in whatever space the reconstructor generated
Returns:
The result 3D model.

Reimplemented from EMAN::Reconstructor.

Definition at line 3599 of file reconstructor.cpp.

References circumf(), EMAN::EMData::depad(), EMAN::EMData::do_ift_inplace(), EMAN::EMData::get_attr(), m_refvol, m_snr, m_vnxc, m_vnyc, m_vnyp, m_vnzc, m_vnzp, m_volume, m_wghta, m_wptr, max, max3d(), min, EMAN::EMData::set_array_offsets(), sqrt(), and EMAN::EMData::symplane0_ctf().

{
        m_volume->set_array_offsets(0, 1, 1);
        m_wptr->set_array_offsets(0, 1, 1);
        cout <<  "  will set refvol  "  <<endl;
        //m_refvol->set_array_offsets(0, 1, 1);
        m_volume->symplane0_ctf(m_wptr);

        int box = 7;
        int vol = box*box*box;
        int kc = (box-1)/2;
        vector< float > pow_a( 3*kc+1, 1.0 );
        for( unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
        pow_a[3*kc]=0.0;


        float max = max3d( kc, pow_a );
        float alpha = ( 1.0f - 1.0f/(float)vol ) / max;
        float osnr = 1.0f/m_snr;


    vector<float> sigma2(m_vnyc+1, 0.0f);
    vector<float> count(m_vnyc+1, 0.0f);

        int ix,iy,iz;
        // compute sigma2
        for (iz = 1; iz <= m_vnzp; iz++) {
                int   izp = (iz<=m_vnzc) ? iz - 1 : iz-m_vnzp-1;
                float argz = float(izp*izp);
                for (iy = 1; iy <= m_vnyp; iy++) {
                        int   iyp = (iy<=m_vnyc) ? iy - 1 : iy-m_vnyp-1;
            float argy = argz + float(iyp*iyp);
                        for (ix = 0; ix <= m_vnxc; ix++) {
                            if(ix>0 || (izp>=0 && (iyp>=0 || izp!=0))) {  //Skip Friedel related values
                    float r = std::sqrt(argy + float(ix*ix));
                    int  ir = int(r);
                    if (ir <= m_vnyc) {
                        float frac = r - float(ir);
                        float qres = 1.0f - frac;
                        float temp = (*m_wptr)(ix,iy,iz);
                        //cout<<" WEIGHTS "<<jx<<"  "<<jy<<"  "<<ir<<"  "<<temp<<"  "<<frac<<endl;
                        //cout<<" WEIGHTS "<<ix<<"  "<<iy-1<<"  "<<iz-1<<"  "<<temp<<"  "<<endl;
                        sigma2[ir]   += temp*qres;
                        sigma2[ir+1] += temp*frac;
                        count[ir]    += qres;
                        count[ir+1]  += frac;
                    }
                }
            }
        }
    }
    for (ix = 0; ix <= m_vnyc+1; ix++) {
        if( sigma2[ix] > 0.0f )  sigma2[ix] = count[ix]/sigma2[ix];
        cout<<"  1/sigma2  "<< ix <<"   "<<sigma2[ix]<<endl;
    }
    // now counter will serve to keep fsc-derived stuff

    for (ix = 0; ix <= m_vnyc+1; ix++)
#ifdef _WIN32
        count[ix] = _cpp_max(0.0f, _cpp_min( 0.999f, (*m_refvol)(ix) ) );
#else
            count[ix] = std::max(0.0f, std::min( 0.999f, (*m_refvol)(ix) ) );
#endif  //_WIN32
    for (ix = 0; ix <= m_vnyc+1; ix++)  count[ix] = count[ix]/(1.0f - count[ix]) * sigma2[ix];
    for (ix = 0; ix <= m_vnyc+1; ix++)  {
        if ( count[ix] >0.0f) count[ix] = 1.0f/count[ix];  //fudge?
    }
for (ix = 0; ix <= m_vnyc+1; ix++)  cout<<"  tau2  "<< ix <<"   "<<count[ix]<<endl;
        // normalize
        for (iz = 1; iz <= m_vnzp; iz++) {
                int   izp = (iz<=m_vnzc) ? iz - 1 : iz-m_vnzp-1;
                float argz = float(izp*izp);
                for (iy = 1; iy <= m_vnyp; iy++) {
                        int   iyp = (iy<=m_vnyc) ? iy - 1 : iy-m_vnyp-1;
            float argy = argz + float(iyp*iyp);
                        for (ix = 0; ix <= m_vnxc; ix++) {
                    float r = std::sqrt(argy + float(ix*ix));
                    int  ir = int(r);
                    if (ir <= m_vnyc) {
                        float frac = r - float(ir);
                        float qres = 1.0f - frac;
                        osnr = qres*count[ir] + frac*count[ir+1];
                        if(osnr == 0.0f)  osnr = 1.0f/(0.001*(*m_wptr)(ix,iy,iz));
                        //cout<<"  "<<iz<<"   "<<iy<<"   "<<"   "<<ix<<"   "<<(*m_wptr)(ix,iy,iz)<<"   "<<osnr<<"      "<<(*m_volume)(2*ix,iy,iz)<<"      "<<(*m_volume)(2*ix+1,iy,iz)<<endl;
                                            float tmp=((*m_wptr)(ix,iy,iz)+osnr);
                                            //if( m_varsnr )  tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+freq*osnr)*m_sign;
                                            //else {
                                            //cout<<"  "<<iz<<"  "<<iy<<"  "<<"  "<<ix<<"  "<<iz<<"  "<<"  "<<(*m_wptr)(ix,iy,iz)<<"  "<<osnr<<"  "<<endl;
                                            if(tmp>0.0f) {
                                                tmp = (-2*((ix+iy+iz)%2)+1)/tmp;
                                            /*
                        int ir = int(freq);
                        float df = freq - float(ir);
                        float add = (1.0f - df)*(*m_refvol)(ir) + df*(*m_refvol)(ir+1);
                        cout<<"  "<<iz<<"  "<<iy<<"  "<<"  "<<ix<<"  "<<ir<<"  "<<"  "<<(*m_wptr)(ix,iy,iz)<<"  "<<add<<"  "<<endl;
                                            tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+add)*m_sign;
                                            */
                                        //}

                        /*if( false) {//m_weighting == ESTIMATE ) {
                        cout <<  "  ESTIMATE  "  <<endl;
                                int cx = ix;
                                int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
                                int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
                                float sum = 0.0;
                                for( int ii = -kc; ii <= kc; ++ii ) {
                                        int nbrcx = cx + ii;
                                        if( nbrcx >= m_vnxc ) continue;
                                        for( int jj= -kc; jj <= kc; ++jj ) {
                                                int nbrcy = cy + jj;
                                                if( nbrcy <= -m_vnyc || nbrcy >= m_vnyc ) continue;
                                                for( int kk = -kc; kk <= kc; ++kk ) {
                                                        int nbrcz = cz + jj;
                                                        if( nbrcz <= -m_vnyc || nbrcz >= m_vnyc ) continue;
                                                        if( nbrcx < 0 ) {
                                                                nbrcx = -nbrcx;
                                                                nbrcy = -nbrcy;
                                                                nbrcz = -nbrcz;
                                                        }

                                                        int nbrix = nbrcx;
                                                        int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 + m_vnyp;
                                                        int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 + m_vnzp;
                                                        if( (*m_wptr)( nbrix, nbriy, nbriz ) == 0.0 ) {
                                                                int c = 3*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
                                                                sum = sum + pow_a[c];
                                                                  // if(ix%20==0 && iy%20==0 && iz%20==0)
                                                                 //   std::cout << boost::format( "%4d %4d %4d %4d %10.3f" ) % nbrix % nbriy % nbriz % c % sum << std::endl;
                                                        }
                                                }
                                        }
                                }
                                float wght = 1.0f / ( 1.0f - alpha * sum );
/
                        if(ix%10==0 && iy%10==0)
                        {
                            std::cout << boost::format( "%4d %4d %4d " ) % ix % iy %iz;
                            std::cout << boost::format( "%10.3f %10.3f %10.3f " )  % tmp % wght % sum;
                            std::  << boost::format( "%10.3f %10.3e " ) % pow_b[r] % alpha;
                            std::cout << std::endl;
                        }
 /
                                tmp = tmp * wght;
                                }*/
                                (*m_volume)(2*ix,iy,iz)   *= tmp;
                                (*m_volume)(2*ix+1,iy,iz) *= tmp;
                                } else {
                                (*m_volume)(2*ix,iy,iz)   = 0.0f;
                                (*m_volume)(2*ix+1,iy,iz) = 0.0f;
                                }
                                }
                        }
                }
        }

        // back fft
        m_volume->do_ift_inplace();
        int npad = m_volume->get_attr("npad");
        m_volume->depad();
        circumf( m_volume, npad );
        m_volume->set_array_offsets( 0, 0, 0 );

        return 0;
}
virtual string EMAN::nn4_ctfwReconstructor::get_desc ( ) const [inline, virtual]

Get a clear, concise description of this class.

Returns:
a clear, concise description of this class

Implements EMAN::FactoryBase.

Definition at line 1342 of file reconstructor.h.

                {
                        return "Direct Fourier inversion reconstruction routine";
                }
virtual string EMAN::nn4_ctfwReconstructor::get_name ( ) const [inline, virtual]

Get the unique name of this class (especially for factory based instantiation access)

Returns:
the unique name of this class

Implements EMAN::FactoryBase.

Definition at line 1337 of file reconstructor.h.

References NAME.

                {
                        return NAME;
                }
TypeDict EMAN::nn4_ctfwReconstructor::get_param_types ( ) const [inline, virtual]
Returns:
a TypeDict defining and describing the feasible parameters of this class

Implements EMAN::FactoryBase.

Definition at line 1353 of file reconstructor.h.

References EMAN::EMObject::EMDATA, EMAN::EMObject::FLOAT, EMAN::EMObject::INT, EMAN::TypeDict::put(), and EMAN::EMObject::STRING.

                {
                        TypeDict d;
                        d.put("size",           EMObject::INT);
                        d.put("npad",           EMObject::INT);
                        d.put("sign",           EMObject::INT);
                        d.put("symmetry",       EMObject::STRING);
                        d.put("snr",            EMObject::FLOAT);
                        d.put("fftvol",         EMObject::EMDATA);
                        d.put("weight",         EMObject::EMDATA);
                        d.put("refvol",         EMObject::EMDATA);
                        d.put("weighting",  EMObject::INT);
                        d.put("varsnr",     EMObject::INT);
                        return d;
                }
int nn4_ctfwReconstructor::insert_padfft_slice_weighted ( EMData padfft,
EMData sigmasq2,
const Transform trans,
const float  weight 
)

Definition at line 3586 of file reconstructor.cpp.

References Assert, EMAN::Transform::get_sym_proj(), m_symmetry, m_volume, m_wptr, and EMAN::EMData::nn_ctfw().

Referenced by insert_slice().

{
        Assert( padfft != NULL );
        //float tmp = padfft->get_attr_default("ctf_applied", 0);
        //int   ctf_applied = (int) tmp;

        vector<Transform> tsym = t.get_sym_proj(m_symmetry);
        for (unsigned int isym=0; isym < tsym.size(); isym++)  m_volume->nn_ctfw(m_wptr, padfft, sigmasq2, tsym[isym], weight);

        return 0;
}
int nn4_ctfwReconstructor::insert_slice ( const EMData *const  slice,
const Transform euler,
const float  weight 
) [virtual]

Insert a slice into a 3D volume, in a given orientation.

Returns:
0 if successful, 1 otherwise
Parameters:
slicethe image slice to be inserted into the 3D volume
eulerEuler angle of this image slice.
weightA weighting factor for this slice, generally the number of particles in a class-average. May be ignored by some reconstructors
Returns:
0 if OK. 1 if error.
Exceptions:
NullPointerExceptionif the input EMData pointer is null
ImageFormatExceptionif the image is complex as opposed to real

for (int ix = 0; ix <= 20; ix++) cout <<" "<<(*sigmasq2)(ix,1) <<" "<<(*padfft)(ix,1);

Reimplemented from EMAN::Reconstructor.

Definition at line 3555 of file reconstructor.cpp.

References checked_delete(), EMAN::EMData::get_attr(), EMAN::EMData::get_attr_default(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), insert_padfft_slice_weighted(), LOGERR, m_npad, m_vnx, and EMAN::padfft_slice().

{
        // sanity checks
        if (!slice) {
                LOGERR("try to insert NULL slice");
                return 1;
        }
        if(weight >0.0f) {
                int padffted= slice->get_attr_default("padffted", 0);
                if( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=m_vnx)  ) {
                        // FIXME: Why doesn't this throw an exception?
                        LOGERR("Tried to insert a slice that is the wrong size.");
                        return 1;
                }

                EMData* padfft = NULL;
                EMData* sigmasq2;
                sigmasq2 = slice->get_attr("sigmasq2");

                if( padffted != 0 ) padfft = new EMData(*slice);
                else                padfft = padfft_slice( slice, t, m_npad );

//cout <<endl;
                insert_padfft_slice_weighted( padfft, sigmasq2, t, weight );

                checked_delete( padfft );
        }
        return 0;
}
static Reconstructor* EMAN::nn4_ctfwReconstructor::NEW ( ) [inline, static]

Definition at line 1347 of file reconstructor.h.

References nn4_ctfwReconstructor().

                {
                        return new nn4_ctfwReconstructor();
                }
void nn4_ctfwReconstructor::setup ( ) [virtual]

Initialize the reconstructor.

Implements EMAN::Reconstructor.

Definition at line 3471 of file reconstructor.cpp.

References EMAN::Dict::has_key(), m_varsnr, EMAN::FactoryBase::params, and sign.

Referenced by nn4_ctfwReconstructor().

{
        if( ! params.has_key("size") ) throw std::logic_error("Error: image size is not given");

        int size = params["size"];
        int npad = params.has_key("npad") ? int(params["npad"]) : 2;
        // int sign = params.has_key("sign") ? int(params["sign"]) : 1;
        int sign = 1;
        string symmetry = params.has_key("symmetry")? params["symmetry"].to_str() : "c1";

        float snr = params["snr"];

        m_varsnr = params.has_key("varsnr") ? int(params["varsnr"]) : 0;
        setup( symmetry, size, npad, snr, sign );

}
void nn4_ctfwReconstructor::setup ( const string &  symmetry,
int  size,
int  npad,
float  snr,
int  sign 
)

Definition at line 3488 of file reconstructor.cpp.

References buildFFTVolume(), buildNormVolume(), ESTIMATE, EMAN::Transform::get_nsym(), EMAN::Dict::has_key(), m_npad, m_nsym, m_refvol, m_sign, m_snr, m_symmetry, m_vnx, m_vnxc, m_vnxp, m_vny, m_vnyc, m_vnyp, m_vnz, m_vnzc, m_vnzp, m_weighting, m_wghta, m_wghtb, NONE, EMAN::FactoryBase::params, and sign.

{
        m_weighting = ESTIMATE;
        if( params.has_key("weighting") ) {
                if( int( params["weighting"])==0 ) m_weighting = NONE;
        }



        m_wghta = 0.2f;
        m_wghtb = 0.004f;

        m_symmetry = symmetry;
        m_npad = npad;
        m_sign = sign;
        m_nsym = Transform::get_nsym(m_symmetry);

        m_snr = snr;

        m_vnx = size;
        m_vny = size;
        m_vnz = size;

        m_vnxp = size*npad;
        m_vnyp = size*npad;
        m_vnzp = size*npad;

        m_vnxc = m_vnxp/2;
        m_vnyc = m_vnyp/2;
        m_vnzc = m_vnzp/2;

        buildFFTVolume();
        buildNormVolume();
        m_refvol = params["refvol"];
}

Member Data Documentation

Definition at line 1384 of file reconstructor.h.

Referenced by buildFFTVolume(), insert_slice(), and setup().

Definition at line 1391 of file reconstructor.h.

Referenced by setup().

Definition at line 1380 of file reconstructor.h.

Referenced by finish(), and setup().

Definition at line 1385 of file reconstructor.h.

Referenced by setup().

Definition at line 1389 of file reconstructor.h.

Referenced by finish(), and setup().

Definition at line 1390 of file reconstructor.h.

Referenced by insert_padfft_slice_weighted(), and setup().

Definition at line 1386 of file reconstructor.h.

Referenced by setup().

Definition at line 1381 of file reconstructor.h.

Referenced by insert_slice(), and setup().

Definition at line 1383 of file reconstructor.h.

Referenced by buildNormVolume(), finish(), and setup().

Definition at line 1382 of file reconstructor.h.

Referenced by buildFFTVolume(), and setup().

Definition at line 1381 of file reconstructor.h.

Referenced by setup().

Definition at line 1383 of file reconstructor.h.

Referenced by finish(), and setup().

Definition at line 1382 of file reconstructor.h.

Referenced by buildFFTVolume(), buildNormVolume(), finish(), and setup().

Definition at line 1381 of file reconstructor.h.

Referenced by setup().

Definition at line 1383 of file reconstructor.h.

Referenced by finish(), and setup().

Definition at line 1382 of file reconstructor.h.

Referenced by buildFFTVolume(), buildNormVolume(), finish(), and setup().

Definition at line 1387 of file reconstructor.h.

Referenced by setup().

Definition at line 1388 of file reconstructor.h.

Referenced by finish(), and setup().

Definition at line 1388 of file reconstructor.h.

Referenced by setup().

const string nn4_ctfwReconstructor::NAME = "nn4_ctfw" [static]

Definition at line 1375 of file reconstructor.h.

Referenced by get_name().


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