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.
EMData * finish ()
 Finish up the averaging and return the result.
vector< float > get_snr () const

Protected Attributes

XYData * sf
EMData * curves
bool need_snr
const char * outfile

Private Attributes

vector< float > snr
EMData * image0_fft
EMData * image0_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 281 of file averager.h.

Constructor & Destructor Documentation

CtfAverager::CtfAverager (  ) 

Definition at line 675 of file averager.cpp.

00675                          :
00676         sf(0), curves(0), need_snr(false), outfile(0),
00677         image0_fft(0), image0_copy(0), snri(0), snrn(0),
00678         tdr(0), tdi(0), tn(0),
00679         filter(0), nimg(0), nx(0), ny(0), nz(0)
00680 {
00681 
00682 }

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 684 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::EMObject::f, EMAN::Util::fast_floor(), filter, EMAN::Dict::get(), EMAN::EMData::get_attr_dict(), EMAN::EMData::get_ctf(), EMAN::EMData::get_data(), EMAN::Averager::get_name(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::EMData::has_ctff(), image0_copy, image0_fft, EMAN::EMUtil::is_same_size(), LOGERR, LOGWARN, nimg, EMAN::Averager::params, EMAN::EMData::set_size(), sf, snri, snrn, sqrt(), tdi, tdr, tn, EMAN::EMData::to_zero(), EMAN::EMData::update(), x, and y.

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

EMData * CtfAverager::finish (  )  [virtual]

Finish up the averaging and return the result.

Returns:
The averaged image.

Implements EMAN::Averager.

Definition at line 852 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::Util::save_data(), snr, snri, snrn, tdi, tdr, tn, EMAN::EMData::update(), x, and y.

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

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

Definition at line 288 of file averager.h.

References snr.

00289                 {
00290                         return snr;
00291                 }

Member Data Documentation

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

Definition at line 294 of file averager.h.

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

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

Definition at line 295 of file averager.h.

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

bool EMAN::CtfAverager::need_snr [protected]

Definition at line 296 of file averager.h.

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

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

Definition at line 297 of file averager.h.

Referenced by finish().

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

Definition at line 299 of file averager.h.

Referenced by finish(), and get_snr().

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

Definition at line 300 of file averager.h.

Referenced by add_image().

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

Definition at line 301 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 303 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 304 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 306 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 307 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 308 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 309 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 310 of file averager.h.

Referenced by add_image(), and finish().

float EMAN::CtfAverager::filter [private]

Definition at line 312 of file averager.h.

Referenced by add_image(), and finish().

int EMAN::CtfAverager::nimg [private]

Definition at line 313 of file averager.h.

Referenced by add_image(), and finish().

int EMAN::CtfAverager::nx [private]

Definition at line 314 of file averager.h.

int EMAN::CtfAverager::ny [private]

Definition at line 315 of file averager.h.

int EMAN::CtfAverager::nz [private]

Definition at line 316 of file averager.h.

The documentation for this class was generated from the following files:
Generated on Sat Nov 21 02:20:16 2009 for EMAN2 by  doxygen 1.5.6