Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

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 333 of file averager.h.


Constructor & Destructor Documentation

CtfAverager::CtfAverager  ) 
 

Definition at line 935 of file averager.cpp.

References nx, and ny.

00935                          :
00936         sf(0), curves(0), need_snr(false), outfile(0),
00937         image0_fft(0), image0_copy(0), snri(0), snrn(0),
00938         tdr(0), tdi(0), tn(0),
00939         filter(0), nimg(0), nx(0), ny(0), nz(0)
00940 {
00941 
00942 }


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 944 of file averager.cpp.

References EMAN::EMData::ap2ri(), EMAN::Ctf::apix, EMAN::Ctf::compute_1d(), EMAN::EMData::copy_head(), ctf, ctfn, 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, nx, ny, EMAN::EMData::set_size(), snri, snrn, sqrt(), tdi, tdr, tn, EMAN::EMData::to_zero(), EMAN::EMData::update(), x, and y.

00945 {
00946         if (!image) {
00947                 return;
00948         }
00949 
00950         if (nimg >= 1 && !EMUtil::is_same_size(image, result)) {
00951                 LOGERR("%sAverager can only process same-size Image",
00952                                                          get_name().c_str());
00953                 return;
00954         }
00955 
00956         if (image->get_zsize() != 1) {
00957                 LOGERR("%sAverager: Only 2D images are currently supported",
00958                                                          get_name().c_str());
00959         }
00960 
00961         string alg_name = get_name();
00962 
00963         if (alg_name == "CtfCW" || alg_name == "CtfCWauto") {
00964                 if (image->get_ctf() != 0 && !image->has_ctff()) {
00965                         LOGERR("%sAverager: Attempted CTF Correction with no ctf parameters",
00966                                                                  get_name().c_str());
00967                 }
00968         }
00969         else {
00970                 if (image->get_ctf() != 0) {
00971                         LOGERR("%sAverager: Attempted CTF Correction with no ctf parameters",
00972                                                                  get_name().c_str());
00973                 }
00974         }
00975 
00976         nimg++;
00977 
00978 
00979         if (nimg == 1) {
00980                 image0_fft = image->do_fft();
00981 
00982                 nx = image0_fft->get_xsize();
00983                 ny = image0_fft->get_ysize();
00984                 nz = image0_fft->get_zsize();
00985 
00986                 result = new EMData();
00987                 result->set_size(nx - 2, ny, nz);
00988 
00989 
00990                 if (alg_name == "Weighting" && curves) {
00991                         if (!sf) {
00992                                 LOGWARN("CTF curve in file will contain relative, not absolute SNR!");
00993                         }
00994                         curves->set_size(Ctf::CTFOS * ny / 2, 3, 1);
00995                         curves->to_zero();
00996                 }
00997 
00998 
00999                 if (alg_name == "CtfC") {
01000                         filter = params["filter"];
01001                         if (filter == 0) {
01002                                 filter = 22.0f;
01003                         }
01004                         float apix_y = image->get_attr_dict().get("apix_y");
01005                         float ds = 1.0f / (apix_y * ny * Ctf::CTFOS);
01006                         filter = 1.0f / (filter * ds);
01007                 }
01008 
01009                 if (alg_name == "CtfCWauto") {
01010                         int nxy2 = nx * ny/2;
01011 
01012                         snri = new float[ny / 2];
01013                         snrn = new float[ny / 2];
01014                         tdr = new float[nxy2];
01015                         tdi = new float[nxy2];
01016                         tn = new float[nxy2];
01017 
01018                         for (int i = 0; i < ny / 2; i++) {
01019                                 snri[i] = 0;
01020                                 snrn[i] = 0;
01021                         }
01022 
01023                         for (int i = 0; i < nxy2; i++) {
01024                                 tdr[i] = 1;
01025                                 tdi[i] = 1;
01026                                 tn[i] = 1;
01027                         }
01028                 }
01029 
01030                 image0_copy = image0_fft->copy_head();
01031                 image0_copy->ap2ri();
01032                 image0_copy->to_zero();
01033         }
01034 
01035         Ctf::CtfType curve_type = Ctf::CTF_AMP;
01036         if (alg_name == "CtfCWauto") {
01037                 curve_type = Ctf::CTF_AMP;
01038         }
01039 
01040         float *src = image->get_data();
01041         image->ap2ri();
01042         Ctf *image_ctf = image->get_ctf();
01043         int ny2 = image->get_ysize();
01044 
01045         vector<float> ctf1 = image_ctf->compute_1d(ny2,1.0f/(image_ctf->apix*image->get_ysize()), curve_type);
01046 
01047         if (ctf1.size() == 0) {
01048                 LOGERR("Unexpected CTF correction problem");
01049         }
01050 
01051         ctf.push_back(ctf1);
01052 
01053         vector<float> ctfn1;
01054         if (sf) {
01055                 ctfn1 = image_ctf->compute_1d(ny2,1.0f/(image_ctf->apix*image->get_ysize()), Ctf::CTF_SNR, sf);
01056         }
01057         else {
01058                 ctfn1 = image_ctf->compute_1d(ny2,1.0f/(image_ctf->apix*image->get_ysize()), Ctf::CTF_SNR);
01059         }
01060 
01061         ctfn.push_back(ctfn1);
01062 
01063         if (alg_name == "CtfCWauto") {
01064                 int j = 0;
01065                 for (int y = 0; y < ny; y++) {
01066                         for (int x = 0; x < nx / 2; x++, j += 2) {
01067 #ifdef  _WIN32
01068                                 float r = (float)_hypot((float)x, (float)(y - ny / 2.0f));
01069 #else
01070                                 float r = (float)hypot((float)x, (float)(y - ny / 2.0f));
01071 #endif  //_WIN32
01072                                 int l = static_cast < int >(Util::fast_floor(r));
01073 
01074                                 if (l >= 0 && l < ny / 2) {
01075                                         int k = y*nx/2 + x;
01076                                         tdr[k] *= src[j];
01077                                         tdi[k] *= src[j + 1];
01078 #ifdef  _WIN32
01079                                         tn[k] *= (float)_hypot(src[j], src[j + 1]);
01080 #else
01081                                         tn[k] *= (float)hypot(src[j], src[j + 1]);
01082 #endif  //_WIN32
01083                                 }
01084                         }
01085                 }
01086         }
01087 
01088 
01089         float *tmp_data = image0_copy->get_data();
01090 
01091         int j = 0;
01092         for (int y = 0; y < ny; y++) {
01093                 for (int x = 0; x < nx / 2; x++, j += 2) {
01094                         float r = Ctf::CTFOS * sqrt(x * x + (y - ny / 2.0f) * (y - ny / 2.0f));
01095                         int l = static_cast < int >(Util::fast_floor(r));
01096                         r -= l;
01097 
01098                         float f = 0;
01099                         if (l <= Ctf::CTFOS * ny / 2 - 2) {
01100                                 f = (ctf1[l] * (1 - r) + ctf1[l + 1] * r);
01101                         }
01102                         tmp_data[j] += src[j] * f;
01103                         tmp_data[j + 1] += src[j + 1] * f;
01104                 }
01105         }
01106 
01107         EMData *image_fft = image->do_fft();
01108         image_fft->update();
01109         if(image_ctf) {delete image_ctf; image_ctf=0;}
01110 }

EMData * CtfAverager::finish  )  [virtual]
 

Finish up the averaging and return the result.

Returns:
The averaged image.

Implements EMAN::Averager.

Definition at line 1112 of file averager.cpp.

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

01113 {
01114         int j = 0;
01115         for (int y = 0; y < ny; y++) {
01116                 for (int x = 0; x < nx / 2; x++, j += 2) {
01117 #ifdef  _WIN32
01118                         float r = (float) _hypot(x, y - ny / 2.0f);
01119 #else
01120                         float r = (float) hypot(x, y - ny / 2.0f);
01121 #endif
01122                         int l = static_cast < int >(Util::fast_floor(r));
01123                         if (l >= 0 && l < ny / 2) {
01124                                 int k = y*nx/2 + x;
01125                                 snri[l] += (tdr[k] + tdi[k]/tn[k]);
01126                                 snrn[l] += 1;
01127                         }
01128                 }
01129         }
01130 
01131         for (int i = 0; i < ny / 2; i++) {
01132                 snri[i] *= nimg / snrn[i];
01133         }
01134 
01135         if(strcmp(outfile, "") != 0) {
01136                 Util::save_data(0, 1, snri, ny / 2, outfile);
01137         }
01138 
01139 
01140         float *cd = 0;
01141         if (curves) {
01142                 cd = curves->get_data();
01143         }
01144 
01145         for (int i = 0; i < Ctf::CTFOS * ny / 2; i++) {
01146                 float ctf0 = 0;
01147                 for (int j = 0; j < nimg; j++) {
01148                         ctf0 += ctfn[j][i];
01149                         if (ctf[j][i] == 0) {
01150                                 ctf[j][i] = 1.0e-12f;
01151                         }
01152 
01153                         if (curves) {
01154                                 cd[i] += ctf[j][i] * ctfn[j][i];
01155                                 cd[i + Ctf::CTFOS * ny / 2] += ctfn[j][i];
01156                                 cd[i + 2 * Ctf::CTFOS * ny / 2] += ctfn[j][i];
01157                         }
01158                 }
01159 
01160                 string alg_name = get_name();
01161 
01162                 if (alg_name == "CtfCW" && need_snr) {
01163                         snr[i] = ctf0;
01164                 }
01165 
01166                 float ctf1 = ctf0;
01167                 if (alg_name == "CtfCWauto") {
01168                         ctf1 = snri[i / Ctf::CTFOS];
01169                 }
01170 
01171                 if (ctf1 <= 0.0001f) {
01172                         ctf1 = 0.1f;
01173                 }
01174 
01175                 if (alg_name == "CtfC") {
01176                         for (int j = 0; j < nimg; j++) {
01177                                 ctf[j][i] = exp(-i * i / (filter * filter)) * ctfn[j][i] / (fabs(ctf[j][i]) * ctf1);
01178                         }
01179                 }
01180                 else if (alg_name == "Weighting") {
01181                         for (int j = 0; j < nimg; j++) {
01182                                 ctf[j][i] = ctfn[j][i] / ctf1;
01183                         }
01184                 }
01185                 else if (alg_name == "CtfCW") {
01186                         for (int j = 0; j < nimg; j++) {
01187                                 ctf[j][i] = (ctf1 / (ctf1 + 1)) * ctfn[j][i] / (ctf[j][i] * ctf1);
01188                         }
01189                 }
01190                 else if (alg_name == "CtfCWauto") {
01191                         for (int j = 0; j < nimg; j++) {
01192                                 ctf[j][i] = ctf1 * ctfn[j][i] / (fabs(ctf[j][i]) * ctf0);
01193                         }
01194                 }
01195         }
01196 
01197 
01198         if (curves) {
01199                 for (int i = 0; i < Ctf::CTFOS * ny / 2; i++) {
01200                         cd[i] /= cd[i + Ctf::CTFOS * ny / 2];
01201                 }
01202                 curves->update();
01203         }
01204 
01205         image0_copy->update();
01206 
01207         float *result_data = result->get_data();
01208         EMData *tmp_ift = image0_copy->do_ift();
01209         float *tmp_ift_data = tmp_ift->get_data();
01210         memcpy(result_data, tmp_ift_data, (nx - 2) * ny * sizeof(float));
01211 
01212         tmp_ift->update();
01213         result->update();
01214         result->set_attr("ptcl_repr",nimg);
01215 
01216         if( image0_copy )
01217         {
01218                 delete image0_copy;
01219                 image0_copy = 0;
01220         }
01221 
01222         if (snri) {
01223                 delete[]snri;
01224                 snri = 0;
01225         }
01226 
01227         if (snrn) {
01228                 delete[]snrn;
01229                 snrn = 0;
01230         }
01231 
01232         if( snri )
01233         {
01234                 delete [] snri;
01235                 snri = 0;
01236         }
01237         if( snrn )
01238         {
01239                 delete [] snrn;
01240                 snrn = 0;
01241         }
01242         if( tdr )
01243         {
01244                 delete [] tdr;
01245                 tdr = 0;
01246         }
01247         if( tdi )
01248         {
01249                 delete [] tdi;
01250                 tdi = 0;
01251         }
01252         if( tn )
01253         {
01254                 delete [] tn;
01255                 tn = 0;
01256         }
01257 
01258         return result;
01259 }

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

Definition at line 340 of file averager.h.

00341                 {
00342                         return snr;
00343                 }


Member Data Documentation

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

Definition at line 355 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 356 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 347 of file averager.h.

Referenced by add_image(), and finish().

float EMAN::CtfAverager::filter [private]
 

Definition at line 364 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 353 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 352 of file averager.h.

Referenced by add_image().

bool EMAN::CtfAverager::need_snr [protected]
 

Definition at line 348 of file averager.h.

int EMAN::CtfAverager::nimg [private]
 

Definition at line 365 of file averager.h.

Referenced by add_image(), and finish().

int EMAN::CtfAverager::nx [private]
 

Definition at line 366 of file averager.h.

int EMAN::CtfAverager::ny [private]
 

Definition at line 367 of file averager.h.

int EMAN::CtfAverager::nz [private]
 

Definition at line 368 of file averager.h.

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

Definition at line 349 of file averager.h.

Referenced by finish().

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

Definition at line 346 of file averager.h.

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

Definition at line 351 of file averager.h.

Referenced by finish().

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

Definition at line 358 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 359 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 361 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 360 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 362 of file averager.h.

Referenced by add_image(), and finish().


The documentation for this class was generated from the following files:
Generated on Tue Jul 12 13:51:13 2011 for EMAN2 by  doxygen 1.3.9.1