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


Constructor & Destructor Documentation

CtfAverager::CtfAverager (  ) 

Definition at line 992 of file averager.cpp.

00992                          :
00993         sf(0), curves(0), need_snr(false), outfile(0),
00994         image0_fft(0), image0_copy(0), snri(0), snrn(0),
00995         tdr(0), tdi(0), tn(0),
00996         filter(0), nimg(0), nx(0), ny(0), nz(0)
00997 {
00998 
00999 }


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 1001 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.

01002 {
01003         if (!image) {
01004                 return;
01005         }
01006 
01007         if (nimg >= 1 && !EMUtil::is_same_size(image, result)) {
01008                 LOGERR("%sAverager can only process same-size Image",
01009                                                          get_name().c_str());
01010                 return;
01011         }
01012 
01013         if (image->get_zsize() != 1) {
01014                 LOGERR("%sAverager: Only 2D images are currently supported",
01015                                                          get_name().c_str());
01016         }
01017 
01018         string alg_name = get_name();
01019 
01020         if (alg_name == "CtfCW" || alg_name == "CtfCWauto") {
01021                 if (image->get_ctf() != 0 && !image->has_ctff()) {
01022                         LOGERR("%sAverager: Attempted CTF Correction with no ctf parameters",
01023                                                                  get_name().c_str());
01024                 }
01025         }
01026         else {
01027                 if (image->get_ctf() != 0) {
01028                         LOGERR("%sAverager: Attempted CTF Correction with no ctf parameters",
01029                                                                  get_name().c_str());
01030                 }
01031         }
01032 
01033         nimg++;
01034 
01035 
01036         if (nimg == 1) {
01037                 image0_fft = image->do_fft();
01038 
01039                 nx = image0_fft->get_xsize();
01040                 ny = image0_fft->get_ysize();
01041                 nz = image0_fft->get_zsize();
01042 
01043                 result = image->copy_head();
01044                 result->set_size(nx - 2, ny, nz);
01045 
01046 
01047                 if (alg_name == "Weighting" && curves) {
01048                         if (!sf) {
01049                                 LOGWARN("CTF curve in file will contain relative, not absolute SNR!");
01050                         }
01051                         curves->set_size(Ctf::CTFOS * ny / 2, 3, 1);
01052                         curves->to_zero();
01053                 }
01054 
01055 
01056                 if (alg_name == "CtfC") {
01057                         filter = params["filter"];
01058                         if (filter == 0) {
01059                                 filter = 22.0f;
01060                         }
01061                         float apix_y = image->get_attr_dict().get("apix_y");
01062                         float ds = 1.0f / (apix_y * ny * Ctf::CTFOS);
01063                         filter = 1.0f / (filter * ds);
01064                 }
01065 
01066                 if (alg_name == "CtfCWauto") {
01067                         int nxy2 = nx * ny/2;
01068 
01069                         snri = new float[ny / 2];
01070                         snrn = new float[ny / 2];
01071                         tdr = new float[nxy2];
01072                         tdi = new float[nxy2];
01073                         tn = new float[nxy2];
01074 
01075                         for (int i = 0; i < ny / 2; i++) {
01076                                 snri[i] = 0;
01077                                 snrn[i] = 0;
01078                         }
01079 
01080                         for (int i = 0; i < nxy2; i++) {
01081                                 tdr[i] = 1;
01082                                 tdi[i] = 1;
01083                                 tn[i] = 1;
01084                         }
01085                 }
01086 
01087                 image0_copy = image0_fft->copy_head();
01088                 image0_copy->ap2ri();
01089                 image0_copy->to_zero();
01090         }
01091 
01092         Ctf::CtfType curve_type = Ctf::CTF_AMP;
01093         if (alg_name == "CtfCWauto") {
01094                 curve_type = Ctf::CTF_AMP;
01095         }
01096 
01097         float *src = image->get_data();
01098         image->ap2ri();
01099         Ctf *image_ctf = image->get_ctf();
01100         int ny2 = image->get_ysize();
01101 
01102         vector<float> ctf1 = image_ctf->compute_1d(ny2,1.0f/(image_ctf->apix*image->get_ysize()), curve_type);
01103 
01104         if (ctf1.size() == 0) {
01105                 LOGERR("Unexpected CTF correction problem");
01106         }
01107 
01108         ctf.push_back(ctf1);
01109 
01110         vector<float> ctfn1;
01111         if (sf) {
01112                 ctfn1 = image_ctf->compute_1d(ny2,1.0f/(image_ctf->apix*image->get_ysize()), Ctf::CTF_SNR, sf);
01113         }
01114         else {
01115                 ctfn1 = image_ctf->compute_1d(ny2,1.0f/(image_ctf->apix*image->get_ysize()), Ctf::CTF_SNR);
01116         }
01117 
01118         ctfn.push_back(ctfn1);
01119 
01120         if (alg_name == "CtfCWauto") {
01121                 int j = 0;
01122                 for (int y = 0; y < ny; y++) {
01123                         for (int x = 0; x < nx / 2; x++, j += 2) {
01124 #ifdef  _WIN32
01125                                 float r = (float)_hypot((float)x, (float)(y - ny / 2.0f));
01126 #else
01127                                 float r = (float)hypot((float)x, (float)(y - ny / 2.0f));
01128 #endif  //_WIN32
01129                                 int l = static_cast < int >(Util::fast_floor(r));
01130 
01131                                 if (l >= 0 && l < ny / 2) {
01132                                         int k = y*nx/2 + x;
01133                                         tdr[k] *= src[j];
01134                                         tdi[k] *= src[j + 1];
01135 #ifdef  _WIN32
01136                                         tn[k] *= (float)_hypot(src[j], src[j + 1]);
01137 #else
01138                                         tn[k] *= (float)hypot(src[j], src[j + 1]);
01139 #endif  //_WIN32
01140                                 }
01141                         }
01142                 }
01143         }
01144 
01145 
01146         float *tmp_data = image0_copy->get_data();
01147 
01148         int j = 0;
01149         for (int y = 0; y < ny; y++) {
01150                 for (int x = 0; x < nx / 2; x++, j += 2) {
01151                         float r = Ctf::CTFOS * sqrt(x * x + (y - ny / 2.0f) * (y - ny / 2.0f));
01152                         int l = static_cast < int >(Util::fast_floor(r));
01153                         r -= l;
01154 
01155                         float f = 0;
01156                         if (l <= Ctf::CTFOS * ny / 2 - 2) {
01157                                 f = (ctf1[l] * (1 - r) + ctf1[l + 1] * r);
01158                         }
01159                         tmp_data[j] += src[j] * f;
01160                         tmp_data[j + 1] += src[j + 1] * f;
01161                 }
01162         }
01163 
01164         EMData *image_fft = image->do_fft();
01165         image_fft->update();
01166         if(image_ctf) {delete image_ctf; image_ctf=0;}
01167 }

EMData * CtfAverager::finish (  )  [virtual]

Finish up the averaging and return the result.

Returns:
The averaged image.

Implements EMAN::Averager.

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

01170 {
01171         int j = 0;
01172         for (int y = 0; y < ny; y++) {
01173                 for (int x = 0; x < nx / 2; x++, j += 2) {
01174 #ifdef  _WIN32
01175                         float r = (float) _hypot(x, y - ny / 2.0f);
01176 #else
01177                         float r = (float) hypot(x, y - ny / 2.0f);
01178 #endif
01179                         int l = static_cast < int >(Util::fast_floor(r));
01180                         if (l >= 0 && l < ny / 2) {
01181                                 int k = y*nx/2 + x;
01182                                 snri[l] += (tdr[k] + tdi[k]/tn[k]);
01183                                 snrn[l] += 1;
01184                         }
01185                 }
01186         }
01187 
01188         for (int i = 0; i < ny / 2; i++) {
01189                 snri[i] *= nimg / snrn[i];
01190         }
01191 
01192         if(strcmp(outfile, "") != 0) {
01193                 Util::save_data(0, 1, snri, ny / 2, outfile);
01194         }
01195 
01196 
01197         float *cd = 0;
01198         if (curves) {
01199                 cd = curves->get_data();
01200         }
01201 
01202         for (int i = 0; i < Ctf::CTFOS * ny / 2; i++) {
01203                 float ctf0 = 0;
01204                 for (int j = 0; j < nimg; j++) {
01205                         ctf0 += ctfn[j][i];
01206                         if (ctf[j][i] == 0) {
01207                                 ctf[j][i] = 1.0e-12f;
01208                         }
01209 
01210                         if (curves) {
01211                                 cd[i] += ctf[j][i] * ctfn[j][i];
01212                                 cd[i + Ctf::CTFOS * ny / 2] += ctfn[j][i];
01213                                 cd[i + 2 * Ctf::CTFOS * ny / 2] += ctfn[j][i];
01214                         }
01215                 }
01216 
01217                 string alg_name = get_name();
01218 
01219                 if (alg_name == "CtfCW" && need_snr) {
01220                         snr[i] = ctf0;
01221                 }
01222 
01223                 float ctf1 = ctf0;
01224                 if (alg_name == "CtfCWauto") {
01225                         ctf1 = snri[i / Ctf::CTFOS];
01226                 }
01227 
01228                 if (ctf1 <= 0.0001f) {
01229                         ctf1 = 0.1f;
01230                 }
01231 
01232                 if (alg_name == "CtfC") {
01233                         for (int j = 0; j < nimg; j++) {
01234                                 ctf[j][i] = exp(-i * i / (filter * filter)) * ctfn[j][i] / (fabs(ctf[j][i]) * ctf1);
01235                         }
01236                 }
01237                 else if (alg_name == "Weighting") {
01238                         for (int j = 0; j < nimg; j++) {
01239                                 ctf[j][i] = ctfn[j][i] / ctf1;
01240                         }
01241                 }
01242                 else if (alg_name == "CtfCW") {
01243                         for (int j = 0; j < nimg; j++) {
01244                                 ctf[j][i] = (ctf1 / (ctf1 + 1)) * ctfn[j][i] / (ctf[j][i] * ctf1);
01245                         }
01246                 }
01247                 else if (alg_name == "CtfCWauto") {
01248                         for (int j = 0; j < nimg; j++) {
01249                                 ctf[j][i] = ctf1 * ctfn[j][i] / (fabs(ctf[j][i]) * ctf0);
01250                         }
01251                 }
01252         }
01253 
01254 
01255         if (curves) {
01256                 for (int i = 0; i < Ctf::CTFOS * ny / 2; i++) {
01257                         cd[i] /= cd[i + Ctf::CTFOS * ny / 2];
01258                 }
01259                 curves->update();
01260         }
01261 
01262         image0_copy->update();
01263 
01264         float *result_data = result->get_data();
01265         EMData *tmp_ift = image0_copy->do_ift();
01266         float *tmp_ift_data = tmp_ift->get_data();
01267         memcpy(result_data, tmp_ift_data, (nx - 2) * ny * sizeof(float));
01268 
01269         tmp_ift->update();
01270         result->update();
01271         result->set_attr("ptcl_repr",nimg);
01272 
01273         if( image0_copy )
01274         {
01275                 delete image0_copy;
01276                 image0_copy = 0;
01277         }
01278 
01279         if (snri) {
01280                 delete[]snri;
01281                 snri = 0;
01282         }
01283 
01284         if (snrn) {
01285                 delete[]snrn;
01286                 snrn = 0;
01287         }
01288 
01289         if( snri )
01290         {
01291                 delete [] snri;
01292                 snri = 0;
01293         }
01294         if( snrn )
01295         {
01296                 delete [] snrn;
01297                 snrn = 0;
01298         }
01299         if( tdr )
01300         {
01301                 delete [] tdr;
01302                 tdr = 0;
01303         }
01304         if( tdi )
01305         {
01306                 delete [] tdi;
01307                 tdi = 0;
01308         }
01309         if( tn )
01310         {
01311                 delete [] tn;
01312                 tn = 0;
01313         }
01314 
01315         return result;
01316 }

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

Definition at line 381 of file averager.h.

References snr.

00382                 {
00383                         return snr;
00384                 }


Member Data Documentation

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

Definition at line 396 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 397 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 388 of file averager.h.

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

float EMAN::CtfAverager::filter [private]

Definition at line 405 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 394 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 393 of file averager.h.

Referenced by add_image().

bool EMAN::CtfAverager::need_snr [protected]

Definition at line 389 of file averager.h.

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

int EMAN::CtfAverager::nimg [private]

Definition at line 406 of file averager.h.

Referenced by add_image(), and finish().

int EMAN::CtfAverager::nx [private]

Definition at line 407 of file averager.h.

int EMAN::CtfAverager::ny [private]

Definition at line 408 of file averager.h.

int EMAN::CtfAverager::nz [private]

Definition at line 409 of file averager.h.

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

Definition at line 390 of file averager.h.

Referenced by finish().

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

Definition at line 387 of file averager.h.

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

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

Definition at line 392 of file averager.h.

Referenced by finish(), and get_snr().

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

Definition at line 399 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 400 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 402 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 401 of file averager.h.

Referenced by add_image(), and finish().

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

Definition at line 403 of file averager.h.

Referenced by add_image(), and finish().


The documentation for this class was generated from the following files:
Generated on Fri Aug 10 16:32:09 2012 for EMAN2 by  doxygen 1.4.7