#include <reconstructor.h>
Inheritance diagram for EMAN::WienerFourierReconstructor:
Public Member Functions | |
WienerFourierReconstructor () | |
Default constructor calls load_default_settings(). | |
virtual | ~WienerFourierReconstructor () |
Deconstructor calls free_memory(). | |
virtual int | insert_slice (const EMData *const slice, const Transform &euler, const float weight=1.0) |
Insert a slice into a 3D volume, in a given orientation. | |
virtual int | determine_slice_agreement (EMData *slice, const Transform &euler, const float weight=1.0, bool sub=true) |
Compares a slice to the current reconstruction volume and computes a normalization factor and quality. | |
virtual EMData * | finish (bool doift=true) |
Get the reconstructed volume Normally will return the volume in real-space with the requested size. | |
virtual string | get_name () const |
Get the unique name of the reconstructor. | |
virtual string | get_desc () const |
Get the one line description of the reconstructor. | |
Static Public Member Functions | |
Reconstructor * | NEW () |
Factory incorporation uses the pointer of this function. | |
Static Public Attributes | |
const string | NAME = "wiener_fourier" |
Protected Member Functions | |
virtual void | do_insert_slice_work (const EMData *const input_slice, const Transform &euler, const float weight) |
A function to perform the nuts and bolts of inserting an image slice. | |
virtual void | do_compare_slice_work (EMData *input_slice, const Transform &euler, float weight) |
A function to perform the nuts and bolts of comparing an image slice. | |
virtual bool | pixel_at (const float &xx, const float &yy, const float &zz, float *dt) |
This is a mode-2 pixel extractor. | |
Private Member Functions | |
WienerFourierReconstructor (const WienerFourierReconstructor &that) | |
Disallow copy construction. | |
WienerFourierReconstructor & | operator= (const WienerFourierReconstructor &) |
Disallow assignment. |
It will perform a reconstruction with a nonisotropic Wiener filter applied to the final reconstruction, and will 'undo' the Wiener filter on the individual class-averages if ctf_wiener_filtered is set. This represents something which was not possible to accomplish in EMAN1, and should produce superior results, with proper anisotropic filtering. Still, the filtration makes the assumption that the original SNR estimates were accurate, and that the data will average completely coherently, which is not truly the case. This may produce models which are somewhat underfiltered in the Wiener sense, but since B-factor corrections are not applied in the ctf.auto averager, this effect is likely already more than compensated for.
Definition at line 547 of file reconstructor.h.
|
Default constructor calls load_default_settings().
Definition at line 553 of file reconstructor.h. 00553 {};
|
|
Deconstructor calls free_memory().
Definition at line 558 of file reconstructor.h. 00558 { }
|
|
Disallow copy construction.
|
|
Compares a slice to the current reconstruction volume and computes a normalization factor and quality. Normalization and quality are returned via attributes set in the passed slice. You may freely mix calls to determine_slice_agreement with calls to insert_slice, but note that determine_slice_agreement can only use information from slices that have already been inserted. Attributes set in the slice are: reconstruct_norm the relative normalization factor which should be applied before inserting the slice reconstruct_qual a scaled quality factor (larger better) for this slice as compared to the existing reconstruction reconstruct_absqual the absolute (not scaled based on weight) quality factor comparing this slice to the existing reconstruction reconstruct_weight the summed weights from all voxels intersecting with the inserted slice, larger -> more overlap with other slices
Reimplemented from EMAN::FourierReconstructor. Definition at line 1064 of file reconstructor.cpp. References EMAN::EMData::copy(), do_compare_slice_work(), do_insert_slice_work(), EMAN::EMData::get_attr(), EMAN::EMData::get_attr_default(), NullPointerException, EMAN::FourierReconstructor::preprocess_slice(), EMAN::EMData::set_attr(), EMAN::Transform::set_mirror(), EMAN::Transform::set_scale(), EMAN::Transform::set_trans(), and weight. 01065 { 01066 // Are these exceptions really necessary? (d.woolford) 01067 if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL"); 01068 01069 Transform * rotation; 01070 rotation = new Transform(arg); // assignment operator 01071 01072 EMData *slice; 01073 if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy(); 01074 else slice = preprocess_slice( input_slice, *rotation); 01075 01076 01077 // We must use only the rotational component of the transform, scaling, translation and mirroring 01078 // are not implemented in Fourier space, but are in preprocess_slice 01079 rotation->set_scale(1.0); 01080 rotation->set_mirror(false); 01081 rotation->set_trans(0,0,0); 01082 01083 // tmp_data->write_image("dbug.hdf",0); 01084 01085 // Remove the current slice first (not threadsafe, but otherwise performance would be awful) 01086 if (sub) do_insert_slice_work(slice, *rotation, -weight); 01087 01088 // Compare 01089 do_compare_slice_work(slice, *rotation,weight); 01090 01091 input_slice->set_attr("reconstruct_norm",slice->get_attr("reconstruct_norm")); 01092 input_slice->set_attr("reconstruct_absqual",slice->get_attr("reconstruct_absqual")); 01093 // input_slice->set_attr("reconstruct_qual",slice->get_attr("reconstruct_qual")); 01094 input_slice->set_attr("reconstruct_weight",slice->get_attr("reconstruct_weight")); 01095 01096 // Now put the slice back 01097 if (sub) do_insert_slice_work(slice, *rotation, weight); 01098 01099 01100 delete rotation; 01101 delete slice; 01102 01103 // image->update(); 01104 return 0; 01105 01106 }
|
|
A function to perform the nuts and bolts of comparing an image slice.
Reimplemented from EMAN::FourierReconstructor. Definition at line 1108 of file reconstructor.cpp. References dot(), dt, EMAN::EMData::get_data(), EMAN::Symmetry3D::get_symmetries(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), pixel_at(), power(), EMAN::EMData::set_attr(), sqrt(), EMAN::Vec3f, weight, x, and y. Referenced by determine_slice_agreement(). 01109 { 01110 01111 float dt[3]; // This stores the complex and weight from the volume 01112 float dt2[2]; // This stores the local image complex 01113 float *dat = input_slice->get_data(); 01114 vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]); 01115 01116 float inx=(float)(input_slice->get_xsize()); // x/y dimensions of the input image 01117 float iny=(float)(input_slice->get_ysize()); 01118 01119 double dot=0; // summed pixel*weight dot product 01120 double vweight=0; // sum of weights 01121 double power=0; // sum of inten*weight from volume 01122 double power2=0; // sum of inten*weight from image 01123 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) { 01124 Transform t3d = arg*(*it); 01125 for (int y = -iny/2; y < iny/2; y++) { 01126 for (int x = 0; x <= inx/2; x++) { 01127 if (x==0 && y==0) continue; // We don't want to use the Fourier origin 01128 01129 float rx = (float) x/(inx-2); // coords relative to Nyquist=.5 01130 float ry = (float) y/iny; 01131 01132 // if ((rx * rx + Util::square(ry - max_input_dim / 2)) > rl) 01133 // continue; 01134 01135 Vec3f coord(rx,ry,0); 01136 coord = coord*t3d; // transpose multiplication 01137 float xx = coord[0]; // transformed coordinates in terms of Nyquist 01138 float yy = coord[1]; 01139 float zz = coord[2]; 01140 01141 01142 if (fabs(xx)>0.5 || fabs(yy)>=0.5 || fabs(zz)>=0.5) continue; 01143 01144 // Map back to actual pixel coordinates in output volume 01145 xx=xx*(nx-2); 01146 yy=yy*ny; 01147 zz=zz*nz; 01148 01149 01150 int idx = (int)(x * 2 + inx*(y<0?iny+y:y)); 01151 dt2[0] = dat[idx]; 01152 dt2[1] = dat[idx+1]; 01153 01154 // value returned indirectly in dt 01155 if (!pixel_at(xx,yy,zz,dt) || dt[2]<=0) continue; 01156 01157 // printf("%f\t%f\t%f\t%f\t%f\n",dt[0],dt[1],dt[2],dt2[0],dt2[1]); 01158 dot+=(dt[0]*dt2[0]+dt[1]*dt2[1])*dt[2]; 01159 vweight+=dt[2]; 01160 power+=(dt[0]*dt[0]+dt[1]*dt[1])*dt[2]; 01161 power2+=(dt2[0]*dt2[0]+dt2[1]*dt2[1])*dt[2]; 01162 } 01163 } 01164 } 01165 01166 dot/=sqrt(power*power2); // normalize the dot product 01167 // input_slice->set_attr("reconstruct_norm",(float)(power2<=0?1.0:sqrt(power/power2)/(inx*iny))); 01168 input_slice->set_attr("reconstruct_norm",(float)(power2<=0?1.0:sqrt(power/power2))); 01169 input_slice->set_attr("reconstruct_absqual",(float)dot); 01170 float rw=weight<=0?1.0f:1.0f/weight; 01171 input_slice->set_attr("reconstruct_qual",(float)(dot*rw/((rw-1.0)*dot+1.0))); // here weight is a proxy for SNR 01172 input_slice->set_attr("reconstruct_weight",(float)vweight/(float)(subnx*subny*subnz)); 01173 // printf("** %g\t%g\t%g\t%g ##\n",dot,vweight,power,power2); 01174 //printf("** %f %f %f ##\n",(float)(power2<=0?1.0:sqrt(power/power2)/(inx*iny)),(float)dot,(float)(dot*weight/((weight-1.0)*dot+1.0))); 01175 }
|
|
A function to perform the nuts and bolts of inserting an image slice.
Reimplemented from EMAN::FourierReconstructor. Definition at line 1009 of file reconstructor.cpp. References EMAN::EMData::get_attr(), EMAN::EMData::get_attr_default(), EMAN::EMData::get_complex_at(), EMAN::Symmetry3D::get_symmetries(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::FourierPixelInserter3D::insert_pixel(), EMAN::Util::linear_interpolate(), sub(), EMAN::Vec3f, weight, x, and y. Referenced by determine_slice_agreement(), and insert_slice(). 01010 { 01011 01012 vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]); 01013 01014 float inx=(float)(input_slice->get_xsize()); // x/y dimensions of the input image 01015 float iny=(float)(input_slice->get_ysize()); 01016 01017 int undo_wiener=(int)input_slice->get_attr_default("ctf_wiener_filtered",0); // indicates whether we need to undo a wiener filter before insertion 01018 // if (undo_wiener) throw UnexpectedBehaviorException("wiener_fourier does not yet accept already Wiener filtered class-averages. Suggest using ctf.auto averager for now."); 01019 01020 vector<float> snr=input_slice->get_attr("ctf_snr_total"); 01021 float sub=1.0; 01022 if (inweight<0) sub=-1.0; 01023 float weight; 01024 01025 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) { 01026 Transform t3d = arg*(*it); 01027 for (int y = -iny/2; y < iny/2; y++) { 01028 for (int x = 0; x <= inx/2; x++) { 01029 01030 float rx = (float) x/(inx-2.0f); // coords relative to Nyquist=.5 01031 float ry = (float) y/iny; 01032 01033 // This deals with the SNR weight 01034 float rn = (float)hypot(rx,ry); 01035 if (rn>=.5) continue; // no SNR in the corners, and we're going to mask them later anyway 01036 rn*=snr.size()*2.0f; 01037 int rni=(int)floor(rn); 01038 if ((unsigned int)rni>=snr.size()-1) weight=snr[snr.size()-1]*sub; 01039 else { 01040 rn-=rni; 01041 weight=Util::linear_interpolate(snr[rni],snr[rni+1],rn); 01042 } 01043 // if (weight>500.0) printf("%f %d %d %f %f %d %f\n",weight,x,y,rx,ry,rni); 01044 01045 Vec3f coord(rx,ry,0); 01046 coord = coord*t3d; // transpose multiplication 01047 float xx = coord[0]; // transformed coordinates in terms of Nyquist 01048 float yy = coord[1]; 01049 float zz = coord[2]; 01050 01051 // Map back to real pixel coordinates in output volume 01052 xx=xx*(nx-2); 01053 yy=yy*ny; 01054 zz=zz*nz; 01055 01056 // printf("%f\n",weight); 01057 if (undo_wiener) inserter->insert_pixel(xx,yy,zz,(input_slice->get_complex_at(x,y))*((weight+1.0f)/weight),weight*sub); 01058 else inserter->insert_pixel(xx,yy,zz,input_slice->get_complex_at(x,y),weight*sub); 01059 } 01060 } 01061 } 01062 }
|
|
Get the reconstructed volume Normally will return the volume in real-space with the requested size. The calling application is responsible for removing any padding.
Reimplemented from EMAN::FourierReconstructor. Definition at line 1251 of file reconstructor.cpp. References EMAN::EMData::depad(), EMAN::EMData::do_ift_inplace(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::Dict::has_key(), EMAN::ReconstructorVolumeData::normalize_threed(), EMAN::EMData::process_inplace(), EMAN::Dict::set_default(), EMAN::EMData::update(), and EMAN::EMData::write_image(). 01252 { 01253 01254 bool sqrtnorm=params.set_default("sqrtnorm",false); 01255 normalize_threed(sqrtnorm,true); // true is the wiener filter 01256 01257 if (doift) { 01258 image->do_ift_inplace(); 01259 image->depad(); 01260 image->process_inplace("xform.phaseorigin.tocenter"); 01261 } 01262 01263 image->update(); 01264 01265 if (params.has_key("savenorm") && strlen((const char *)params["savenorm"])>0) { 01266 if (tmp_data->get_ysize()%2==0 && tmp_data->get_zsize()%2==0) tmp_data->process_inplace("xform.fourierorigin.tocenter"); 01267 tmp_data->write_image((const char *)params["savenorm"]); 01268 } 01269 01270 delete tmp_data; 01271 tmp_data=0; 01272 EMData *ret=image; 01273 image=0; 01274 01275 return ret; 01276 }
|
|
Get the one line description of the reconstructor.
Reimplemented from EMAN::FourierReconstructor. Definition at line 608 of file reconstructor.h. 00609 { 00610 return "Reconstruction via direct Fourier methods using one of a variety of different kernels, most of which are Gaussian based. This version also incorporates a nonisotropic Wiener filter based on SNR estimates stored in the class-average headers by the ctf.auto averager."; 00611 }
|
|
Get the unique name of the reconstructor.
Reimplemented from EMAN::FourierReconstructor. Definition at line 601 of file reconstructor.h. 00602 {
00603 return NAME;
00604 }
|
|
Insert a slice into a 3D volume, in a given orientation.
Reimplemented from EMAN::FourierReconstructor. Definition at line 973 of file reconstructor.cpp. References EMAN::EMData::copy(), do_insert_slice_work(), EMAN::EMData::get_attr_default(), EMAN::EMData::has_attr(), NotExistingObjectException, NullPointerException, EMAN::FourierReconstructor::preprocess_slice(), EMAN::Transform::set_mirror(), EMAN::Transform::set_scale(), EMAN::Transform::set_trans(), and weight. 00974 { 00975 // Are these exceptions really necessary? (d.woolford) 00976 if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL"); 00977 00978 Transform * rotation; 00979 /* if ( input_slice->has_attr("xform.projection") ) { 00980 rotation = (Transform*) (input_slice->get_attr("xform.projection")); // assignment operator 00981 } else {*/ 00982 rotation = new Transform(arg); // assignment operator 00983 // } 00984 00985 if (!input_slice->has_attr("ctf_snr_total")) 00986 throw NotExistingObjectException("ctf_snr_total","No SNR information present in class-average. Must use the ctf.auto or ctfw.auto averager."); 00987 00988 EMData *slice; 00989 if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy(); 00990 else slice = preprocess_slice( input_slice, *rotation); 00991 00992 00993 // We must use only the rotational component of the transform, scaling, translation and mirroring 00994 // are not implemented in Fourier space, but are in preprocess_slice 00995 rotation->set_scale(1.0); 00996 rotation->set_mirror(false); 00997 rotation->set_trans(0,0,0); 00998 00999 // Finally to the pixel wise slice insertion 01000 do_insert_slice_work(slice, *rotation, weight); 01001 01002 delete rotation; rotation=0; 01003 delete slice; 01004 01005 // image->update(); 01006 return 0; 01007 }
|
|
Factory incorporation uses the pointer of this function.
Reimplemented from EMAN::FourierReconstructor. Definition at line 616 of file reconstructor.h. 00617 { 00618 return new WienerFourierReconstructor(); 00619 }
|
|
Disallow assignment.
|
|
This is a mode-2 pixel extractor.
Reimplemented from EMAN::FourierReconstructor. Definition at line 1177 of file reconstructor.cpp. References dt, EMAN::Util::fast_exp(), EMAN::EMData::get_complex_index(), EMAN::EMData::get_complex_index_fast(), EMAN::EMData::get_data(), EMAN::Util::hypot3sq(), norm(), and rdata. Referenced by do_compare_slice_work(). 01178 { 01179 int x0 = (int) floor(xx); 01180 int y0 = (int) floor(yy); 01181 int z0 = (int) floor(zz); 01182 01183 float *rdata=image->get_data(); 01184 float *norm=tmp_data->get_data(); 01185 float normsum=0,normsum2=0; 01186 01187 dt[0]=dt[1]=dt[2]=0.0; 01188 01189 if (nx==subnx) { // normal full reconstruction 01190 if (x0<-nx2-1 || y0<-ny2-1 || z0<-nz2-1 || x0>nx2 || y0>ny2 || z0>nz2 ) return false; 01191 01192 // no error checking on add_complex_fast, so we need to be careful here 01193 int x1=x0+1; 01194 int y1=y0+1; 01195 int z1=z0+1; 01196 if (x0<-nx2) x0=-nx2; 01197 if (x1>nx2) x1=nx2; 01198 if (y0<-ny2) y0=-ny2; 01199 if (y1>ny2) y1=ny2; 01200 if (z0<-nz2) z0=-nz2; 01201 if (z1>nz2) z1=nz2; 01202 01203 size_t idx=0; 01204 float r, gg; 01205 for (int k = z0 ; k <= z1; k++) { 01206 for (int j = y0 ; j <= y1; j++) { 01207 for (int i = x0; i <= x1; i ++) { 01208 r = Util::hypot3sq((float) i - xx, j - yy, k - zz); 01209 idx=image->get_complex_index_fast(i,j,k); 01210 gg = Util::fast_exp(-r / EMConsts::I2G); 01211 01212 dt[0]+=gg*rdata[idx]; 01213 dt[1]+=(i<0?-1.0f:1.0f)*gg*rdata[idx+1]; 01214 dt[2]+=norm[idx/2]*gg; 01215 normsum2+=gg; 01216 normsum+=gg*norm[idx/2]; 01217 } 01218 } 01219 } 01220 if (normsum==0) return false; 01221 dt[0]/=normsum; 01222 dt[1]/=normsum; 01223 dt[2]/=normsum2; 01224 // printf("%1.2f,%1.2f,%1.2f\t%1.3f\t%1.3f\t%1.3f\t%1.3f\t%1.3f\n",xx,yy,zz,dt[0],dt[1],dt[2],rdata[idx],rdata[idx+1]); 01225 return true; 01226 } 01227 else { // for subvolumes, not optimized yet 01228 size_t idx; 01229 float r, gg; 01230 for (int k = z0 ; k <= z0 + 1; k++) { 01231 for (int j = y0 ; j <= y0 + 1; j++) { 01232 for (int i = x0; i <= x0 + 1; i ++) { 01233 r = Util::hypot3sq((float) i - xx, j - yy, k - zz); 01234 idx=image->get_complex_index(i,j,k,subx0,suby0,subz0,nx,ny,nz); 01235 gg = Util::fast_exp(-r / EMConsts::I2G)*norm[idx/2]; 01236 01237 dt[0]+=gg*rdata[idx]; 01238 dt[1]+=(i<0?-1.0f:1.0f)*gg*rdata[idx+1]; 01239 dt[2]+=norm[idx/2]; 01240 normsum+=gg; 01241 } 01242 } 01243 } 01244 01245 if (normsum==0) return false; 01246 return true; 01247 } 01248 }
|
|
Reimplemented from EMAN::FourierReconstructor. Definition at line 79 of file reconstructor.cpp. |