#include <averager.h>
Inheritance diagram for EMAN::CtfAverager:
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 |
Definition at line 374 of file averager.h.
|
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 }
|
|
To add an image to the Averager. This image will be averaged in this function.
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, 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. 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 }
|
|
Finish up the averaging and return the result.
Implements EMAN::Averager. Definition at line 1169 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. 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 }
|
|
Definition at line 381 of file averager.h. 00382 {
00383 return snr;
00384 }
|
|
Definition at line 396 of file averager.h. Referenced by add_image(), and finish(). |
|
Definition at line 397 of file averager.h. Referenced by add_image(), and finish(). |
|
Definition at line 388 of file averager.h. Referenced by add_image(), and finish(). |
|
Definition at line 405 of file averager.h. Referenced by add_image(), and finish(). |
|
Definition at line 394 of file averager.h. Referenced by add_image(), and finish(). |
|
Definition at line 393 of file averager.h. Referenced by add_image(). |
|
Definition at line 389 of file averager.h. |
|
Definition at line 406 of file averager.h. Referenced by add_image(), and finish(). |
|
Definition at line 407 of file averager.h. |
|
Definition at line 408 of file averager.h. |
|
Definition at line 409 of file averager.h. |
|
Definition at line 390 of file averager.h. Referenced by finish(). |
|
Definition at line 387 of file averager.h. |
|
Definition at line 392 of file averager.h. Referenced by finish(). |
|
Definition at line 399 of file averager.h. Referenced by add_image(), and finish(). |
|
Definition at line 400 of file averager.h. Referenced by add_image(), and finish(). |
|
Definition at line 402 of file averager.h. Referenced by add_image(), and finish(). |
|
Definition at line 401 of file averager.h. Referenced by add_image(), and finish(). |
|
Definition at line 403 of file averager.h. Referenced by add_image(), and finish(). |