EMAN::CtfAverager Class Reference

CtfAverager is the base Averager class for CTF correction or SNR weighting. More...

#include <averager.h>

Inheritance diagram for EMAN::CtfAverager:

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

Collaboration graph
[legend]
List of all members.

Public Member Functions

 CtfAverager ()
void add_image (EMData *image)
 To add an image to the Averager.
EMDatafinish ()
 Finish up the averaging and return the result.
vector< float > get_snr () const

Protected Attributes

XYDatasf
EMDatacurves
bool need_snr
const char * outfile

Private Attributes

vector< float > snr
EMDataimage0_fft
EMDataimage0_copy
vector< vector< float > > ctf
vector< vector< float > > ctfn
float * snri
float * snrn
float * tdr
float * tdi
float * tn
float filter
int nimg
int nx
int ny
int nz

Detailed Description

CtfAverager is the base Averager class for CTF correction or SNR weighting.

Definition at line 288 of file averager.h.


Constructor & Destructor Documentation

CtfAverager::CtfAverager (  ) 

Definition at line 685 of file averager.cpp.

00685                          :
00686         sf(0), curves(0), need_snr(false), outfile(0),
00687         image0_fft(0), image0_copy(0), snri(0), snrn(0),
00688         tdr(0), tdi(0), tn(0),
00689         filter(0), nimg(0), nx(0), ny(0), nz(0)
00690 {
00691 
00692 }


Member Function Documentation

void CtfAverager::add_image ( EMData image  )  [virtual]

To add an image to the Averager.

This image will be averaged in this function.

Parameters:
image The image to be averaged.

Implements EMAN::Averager.

Definition at line 694 of file averager.cpp.

References EMAN::EMData::ap2ri(), EMAN::Ctf::apix, EMAN::Ctf::compute_1d(), EMAN::EMData::copy_head(), ctf, EMAN::Ctf::CTF_AMP, EMAN::Ctf::CTF_SNR, ctfn, EMAN::Ctf::CTFOS, curves, EMAN::EMData::do_fft(), EMAN::Util::fast_floor(), filter, EMAN::EMData::get_data(), EMAN::Averager::get_name(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), image0_copy, image0_fft, EMAN::EMUtil::is_same_size(), LOGERR, LOGWARN, nimg, EMAN::Averager::params, EMAN::Averager::result, EMAN::EMData::set_size(), sf, snri, snrn, sqrt(), tdi, tdr, tn, EMAN::EMData::to_zero(), EMAN::EMData::update(), x, and y.

00695 {
00696         if (!image) {
00697                 return;
00698         }
00699 
00700         if (nimg >= 1 && !EMUtil::is_same_size(image, result)) {
00701                 LOGERR("%sAverager can only process same-size Image",
00702                                                          get_name().c_str());
00703                 return;
00704         }
00705 
00706         if (image->get_zsize() != 1) {
00707                 LOGERR("%sAverager: Only 2D images are currently supported",
00708                                                          get_name().c_str());
00709         }
00710 
00711         string alg_name = get_name();
00712 
00713         if (alg_name == "CtfCW" || alg_name == "CtfCWauto") {
00714                 if (image->get_ctf() != 0 && !image->has_ctff()) {
00715                         LOGERR("%sAverager: Attempted CTF Correction with no ctf parameters",
00716                                                                  get_name().c_str());
00717                 }
00718         }
00719         else {
00720                 if (image->get_ctf() != 0) {
00721                         LOGERR("%sAverager: Attempted CTF Correction with no ctf parameters",
00722                                                                  get_name().c_str());
00723                 }
00724         }
00725 
00726         nimg++;
00727 
00728 
00729         if (nimg == 1) {
00730                 image0_fft = image->do_fft();
00731 
00732                 nx = image0_fft->get_xsize();
00733                 ny = image0_fft->get_ysize();
00734                 nz = image0_fft->get_zsize();
00735 
00736                 result = new EMData();
00737                 result->set_size(nx - 2, ny, nz);
00738 
00739 
00740                 if (alg_name == "Weighting" && curves) {
00741                         if (!sf) {
00742                                 LOGWARN("CTF curve in file will contain relative, not absolute SNR!");
00743                         }
00744                         curves->set_size(Ctf::CTFOS * ny / 2, 3, 1);
00745                         curves->to_zero();
00746                 }
00747 
00748 
00749                 if (alg_name == "CtfC") {
00750                         filter = params["filter"];
00751                         if (filter == 0) {
00752                                 filter = 22.0f;
00753                         }
00754                         float apix_y = image->get_attr_dict().get("apix_y");
00755                         float ds = 1.0f / (apix_y * ny * Ctf::CTFOS);
00756                         filter = 1.0f / (filter * ds);
00757                 }
00758 
00759                 if (alg_name == "CtfCWauto") {
00760                         int nxy2 = nx * ny/2;
00761 
00762                         snri = new float[ny / 2];
00763                         snrn = new float[ny / 2];
00764                         tdr = new float[nxy2];
00765                         tdi = new float[nxy2];
00766                         tn = new float[nxy2];
00767 
00768                         for (int i = 0; i < ny / 2; i++) {
00769                                 snri[i] = 0;
00770                                 snrn[i] = 0;
00771                         }
00772 
00773                         for (int i = 0; i < nxy2; i++) {
00774                                 tdr[i] = 1;
00775                                 tdi[i] = 1;
00776                                 tn[i] = 1;
00777                         }
00778                 }
00779 
00780                 image0_copy = image0_fft->copy_head();
00781                 image0_copy->ap2ri();
00782                 image0_copy->to_zero();
00783         }
00784 
00785         Ctf::CtfType curve_type = Ctf::CTF_AMP;
00786         if (alg_name == "CtfCWauto") {
00787                 curve_type = Ctf::CTF_AMP;
00788         }
00789 
00790         float *src = image->get_data();
00791         image->ap2ri();
00792         Ctf *image_ctf = image->get_ctf();
00793         int ny2 = image->get_ysize();
00794 
00795         vector<float> ctf1 = image_ctf->compute_1d(ny2,1.0f/(image_ctf->apix*image->get_ysize()), curve_type);
00796 
00797         if (ctf1.size() == 0) {
00798                 LOGERR("Unexpected CTF correction problem");
00799         }
00800 
00801         ctf.push_back(ctf1);
00802 
00803         vector<float> ctfn1;
00804         if (sf) {
00805                 ctfn1 = image_ctf->compute_1d(ny2,1.0f/(image_ctf->apix*image->get_ysize()), Ctf::CTF_SNR, sf);
00806         }
00807         else {
00808                 ctfn1 = image_ctf->compute_1d(ny2,1.0f/(image_ctf->apix*image->get_ysize()), Ctf::CTF_SNR);
00809         }
00810 
00811         ctfn.push_back(ctfn1);
00812 
00813         if (alg_name == "CtfCWauto") {
00814                 int j = 0;
00815                 for (int y = 0; y < ny; y++) {
00816                         for (int x = 0; x < nx / 2; x++, j += 2) {
00817 #ifdef  _WIN32
00818                                 float r = (float)_hypot((float)x, (float)(y - ny / 2.0f));
00819 #else
00820                                 float r = (float)hypot((float)x, (float)(y - ny / 2.0f));
00821 #endif  //_WIN32
00822                                 int l = static_cast < int >(Util::fast_floor(r));
00823 
00824                                 if (l >= 0 && l < ny / 2) {
00825                                         int k = y*nx/2 + x;
00826                                         tdr[k] *= src[j];
00827                                         tdi[k] *= src[j + 1];
00828 #ifdef  _WIN32
00829                                         tn[k] *= (float)_hypot(src[j], src[j + 1]);
00830 #else
00831                                         tn[k] *= (float)hypot(src[j], src[j + 1]);
00832 #endif  //_WIN32
00833                                 }
00834                         }
00835                 }
00836         }
00837 
00838 
00839         float *tmp_data = image0_copy->get_data();
00840 
00841         int j = 0;
00842         for (int y = 0; y < ny; y++) {
00843                 for (int x = 0; x < nx / 2; x++, j += 2) {
00844                         float r = Ctf::CTFOS * sqrt(x * x + (y - ny / 2.0f) * (y - ny / 2.0f));
00845                         int l = static_cast < int >(Util::fast_floor(r));
00846                         r -= l;
00847 
00848                         float f = 0;
00849                         if (l <= Ctf::CTFOS * ny / 2 - 2) {
00850                                 f = (ctf1[l] * (1 - r) + ctf1[l + 1] * r);
00851                         }
00852                         tmp_data[j] += src[j] * f;
00853                         tmp_data[j + 1] += src[j + 1] * f;
00854                 }
00855         }
00856 
00857         EMData *image_fft = image->do_fft();
00858         image_fft->update();
00859         if(image_ctf) {delete image_ctf; image_ctf=0;}
00860 }

EMData * CtfAverager::finish (  )  [virtual]

Finish up the averaging and return the result.

Returns:
The averaged image.

Implements EMAN::Averager.

Definition at line 862 of file averager.cpp.

References ctf, ctfn, EMAN::Ctf::CTFOS, curves, EMAN::EMData::do_ift(), EMAN::Util::fast_floor(), filter, EMAN::EMData::get_data(), EMAN::Averager::get_name(), image0_copy, need_snr, nimg, outfile, EMAN::Averager::result, EMAN::Util::save_data(), snr, snri, snrn, tdi, tdr, tn, EMAN::EMData::update(), x, and y.

00863 {
00864         int j = 0;
00865         for (int y = 0; y < ny; y++) {
00866                 for (int x = 0; x < nx / 2; x++, j += 2) {
00867 #ifdef  _WIN32
00868                         float r = (float) _hypot(x, y - ny / 2.0f);
00869 #else
00870                         float r = (float) hypot(x, y - ny / 2.0f);
00871 #endif
00872                         int l = static_cast < int >(Util::fast_floor(r));
00873                         if (l >= 0 && l < ny / 2) {
00874                                 int k = y*nx/2 + x;
00875                                 snri[l] += (tdr[k] + tdi[k]/tn[k]);
00876                                 snrn[l] += 1;
00877                         }
00878                 }
00879         }
00880 
00881         for (int i = 0; i < ny / 2; i++) {
00882                 snri[i] *= nimg / snrn[i];
00883         }
00884 
00885         if(strcmp(outfile, "") != 0) {
00886                 Util::save_data(0, 1, snri, ny / 2, outfile);
00887         }
00888 
00889 
00890         float *cd = 0;
00891         if (curves) {
00892                 cd = curves->get_data();
00893         }
00894 
00895         for (int i = 0; i < Ctf::CTFOS * ny / 2; i++) {
00896                 float ctf0 = 0;
00897                 for (int j = 0; j < nimg; j++) {
00898                         ctf0 += ctfn[j][i];
00899                         if (ctf[j][i] == 0) {
00900                                 ctf[j][i] = 1.0e-12f;
00901                         }
00902 
00903                         if (curves) {
00904                                 cd[i] += ctf[j][i] * ctfn[j][i];
00905                                 cd[i + Ctf::CTFOS * ny / 2] += ctfn[j][i];
00906                                 cd[i + 2 * Ctf::CTFOS * ny / 2] += ctfn[j][i];
00907                         }
00908                 }
00909 
00910                 string alg_name = get_name();
00911 
00912                 if (alg_name == "CtfCW" && need_snr) {
00913                         snr[i] = ctf0;
00914                 }
00915 
00916                 float ctf1 = ctf0;
00917                 if (alg_name == "CtfCWauto") {
00918                         ctf1 = snri[i / Ctf::CTFOS];
00919                 }
00920 
00921                 if (ctf1 <= 0.0001f) {
00922                         ctf1 = 0.1f;
00923                 }
00924 
00925                 if (alg_name == "CtfC") {
00926                         for (int j = 0; j < nimg; j++) {
00927                                 ctf[j][i] = exp(-i * i / (filter * filter)) * ctfn[j][i] / (fabs(ctf[j][i]) * ctf1);
00928                         }
00929                 }
00930                 else if (alg_name == "Weighting") {
00931                         for (int j = 0; j < nimg; j++) {
00932                                 ctf[j][i] = ctfn[j][i] / ctf1;
00933                         }
00934                 }
00935                 else if (alg_name == "CtfCW") {
00936                         for (int j = 0; j < nimg; j++) {
00937                                 ctf[j][i] = (ctf1 / (ctf1 + 1)) * ctfn[j][i] / (ctf[j][i] * ctf1);
00938                         }
00939                 }
00940                 else if (alg_name == "CtfCWauto") {
00941                         for (int j = 0; j < nimg; j++) {
00942                                 ctf[j][i] = ctf1 * ctfn[j][i] / (fabs(ctf[j][i]) * ctf0);
00943                         }
00944                 }
00945         }
00946 
00947 
00948         if (curves) {
00949                 for (int i = 0; i < Ctf::CTFOS * ny / 2; i++) {
00950                         cd[i] /= cd[i + Ctf::CTFOS * ny / 2];
00951                 }
00952                 curves->update();
00953         }
00954 
00955         image0_copy->update();
00956 
00957         float *result_data = result->get_data();
00958         EMData *tmp_ift = image0_copy->do_ift();
00959         float *tmp_ift_data = tmp_ift->get_data();
00960         memcpy(result_data, tmp_ift_data, (nx - 2) * ny * sizeof(float));
00961 
00962         tmp_ift->update();
00963         result->update();
00964 
00965         if( image0_copy )
00966         {
00967                 delete image0_copy;
00968                 image0_copy = 0;
00969         }
00970 
00971         if (snri) {
00972                 delete[]snri;
00973                 snri = 0;
00974         }
00975 
00976         if (snrn) {
00977                 delete[]snrn;
00978                 snrn = 0;
00979         }
00980 
00981         if( snri )
00982         {
00983                 delete [] snri;
00984                 snri = 0;
00985         }
00986         if( snrn )
00987         {
00988                 delete [] snrn;
00989                 snrn = 0;
00990         }
00991         if( tdr )
00992         {
00993                 delete [] tdr;
00994                 tdr = 0;
00995         }
00996         if( tdi )
00997         {
00998                 delete [] tdi;
00999                 tdi = 0;
01000         }
01001         if( tn )
01002         {
01003                 delete [] tn;
01004                 tn = 0;
01005         }
01006 
01007         return result;
01008 }

vector< float > EMAN::CtfAverager::get_snr (  )  const [inline]

Definition at line 295 of file averager.h.

References snr.

00296                 {
00297                         return snr;
00298                 }


Member Data Documentation

vector<vector<float> > EMAN::CtfAverager::ctf [private]

Definition at line 310 of file averager.h.

Referenced by add_image(), and finish().

vector<vector<float> > EMAN::CtfAverager::ctfn [private]

Definition at line 311 of file averager.h.

Referenced by add_image(), and finish().

EMData* EMAN::CtfAverager::curves [protected]

Definition at line 302 of file averager.h.

Referenced by add_image(), finish(), and EMAN::WeightingAverager::set_params().

float EMAN::CtfAverager::filter [private]

Definition at line 319 of file averager.h.

Referenced by add_image(), and finish().

EMData* EMAN::CtfAverager::image0_copy [private]

Definition at line 308 of file averager.h.

Referenced by add_image(), and finish().

EMData* EMAN::CtfAverager::image0_fft [private]

Definition at line 307 of file averager.h.

Referenced by add_image().

bool EMAN::CtfAverager::need_snr [protected]

Definition at line 303 of file averager.h.

Referenced by finish(), and EMAN::CtfCWAverager::set_params().

int EMAN::CtfAverager::nimg [private]

Definition at line 320 of file averager.h.

Referenced by add_image(), and finish().

int EMAN::CtfAverager::nx [private]

Definition at line 321 of file averager.h.

int EMAN::CtfAverager::ny [private]

Definition at line 322 of file averager.h.

int EMAN::CtfAverager::nz [private]

Definition at line 323 of file averager.h.

const char* EMAN::CtfAverager::outfile [protected]

Definition at line 304 of file averager.h.

Referenced by finish().

XYData* EMAN::CtfAverager::sf [protected]

Definition at line 301 of file averager.h.

Referenced by add_image(), and EMAN::WeightingAverager::set_params().

vector< float > EMAN::CtfAverager::snr [mutable, private]

Definition at line 306 of file averager.h.

Referenced by finish(), and get_snr().

float* EMAN::CtfAverager::snri [private]

Definition at line 313 of file averager.h.

Referenced by add_image(), and finish().

float* EMAN::CtfAverager::snrn [private]

Definition at line 314 of file averager.h.

Referenced by add_image(), and finish().

float* EMAN::CtfAverager::tdi [private]

Definition at line 316 of file averager.h.

Referenced by add_image(), and finish().

float* EMAN::CtfAverager::tdr [private]

Definition at line 315 of file averager.h.

Referenced by add_image(), and finish().

float* EMAN::CtfAverager::tn [private]

Definition at line 317 of file averager.h.

Referenced by add_image(), and finish().


The documentation for this class was generated from the following files:
Generated on Tue May 25 17:15:48 2010 for EMAN2 by  doxygen 1.4.7