EMAN::BaldwinWoolfordReconstructor Class Reference

A more mathematically rigorou Fourier reconstructor That did not seem to outperform the FourierReconstructor. More...

#include <reconstructor.h>

Inheritance diagram for EMAN::BaldwinWoolfordReconstructor:

Inheritance graph
[legend]
Collaboration diagram for EMAN::BaldwinWoolfordReconstructor:

Collaboration graph
[legend]

List of all members.

Public Member Functions

 BaldwinWoolfordReconstructor ()
virtual ~BaldwinWoolfordReconstructor ()
 Deconstructor.
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.
virtual TypeDict get_param_types () const
 Get the parameter types of this object.
virtual EMData * finish ()
 Finish reconstruction and return the complete model.
virtual int insert_slice_weights (const Transform &t3d)
virtual int insert_slice (const EMData *const image, const Transform3D &t3d)
 Insert an image slice to the reconstructor.
virtual int insert_slice (const EMData *const image, const Transform &t3d)
 Insert a slice into a 3D volume, in a given orientation.
void insert_pixel (const float &x, const float &y, const float &z, const float dt[2])
void insert_density_at (const float &x, const float &y, const float &z)
virtual void setup ()
 Setup the Fourier reconstructor relies on the correct parameters throws a variety of exceptions if the parameters are unworkable
Exceptions:
InvalidValueException when the x_in parameter is odd or less than zero
InvalidValueException when the y_in parameter is odd or less than zero
InvalidValueException when the x_out parameter is odd or less than zero
InvalidValueException when the y_out parameter is odd or less than zero
InvalidValueException when the z_out parameter is odd or less than zero
InvalidValueException when the pad parameter is less than the greatest dimension of the input images (the greater of x_in and y_in) Note the restraint that the Fourier volumes should be even will be lifted once debugging of the the xform.fourierorigin processor is performed, but however, when this happens a rigorous testing should be performed because no testing has been done in this regard.


Static Public Member Functions

static Reconstructor * NEW ()
 Factory themed method allocating a new FourierReconstructor.

Protected Member Functions

void load_default_settings ()
 Load default settings.

Private Member Functions

 BaldwinWoolfordReconstructor (const BaldwinWoolfordReconstructor &that)
 Disallow copy construction.
BaldwinWoolfordReconstructor & operator= (const BaldwinWoolfordReconstructor &)
 Disallow assignment.

Private Attributes

float * W
float dfreq

Detailed Description

A more mathematically rigorou Fourier reconstructor That did not seem to outperform the FourierReconstructor.

Author:
David Woolford (programming) and Phil Baldwin (mathematics)
Date:
early 2008

Definition at line 563 of file reconstructor.h.

Constructor & Destructor Documentation

EMAN::BaldwinWoolfordReconstructor::BaldwinWoolfordReconstructor (  )  [inline]

Definition at line 566 of file reconstructor.h.

Referenced by NEW().

00566 : W(0) {}

virtual EMAN::BaldwinWoolfordReconstructor::~BaldwinWoolfordReconstructor (  )  [inline, virtual]

Deconstructor.

Definition at line 570 of file reconstructor.h.

References W.

00570                                                         {
00571                         if ( W != 0 )
00572                                 delete [] W;
00573                 }

EMAN::BaldwinWoolfordReconstructor::BaldwinWoolfordReconstructor ( const BaldwinWoolfordReconstructor &  that  )  [private]

Disallow copy construction.

Member Function Documentation

virtual string EMAN::BaldwinWoolfordReconstructor::get_name (  )  const [inline, virtual]

Get the unique name of the reconstructor.

Reimplemented from EMAN::FourierReconstructor.

Definition at line 577 of file reconstructor.h.

00578                 {
00579                         return "baldwinwoolford";
00580                 }

virtual string EMAN::BaldwinWoolfordReconstructor::get_desc (  )  const [inline, virtual]

Get the one line description of the reconstructor.

Reimplemented from EMAN::FourierReconstructor.

Definition at line 584 of file reconstructor.h.

00585                 {
00586                         return "Reconstruction via direct Fourier inversion using gridding and delta function based weights";
00587                 }

static Reconstructor* EMAN::BaldwinWoolfordReconstructor::NEW (  )  [inline, static]

Factory themed method allocating a new FourierReconstructor.

Returns:
a Reconstructor pointer to a newly allocated FourierReconstructor

Reimplemented from EMAN::FourierReconstructor.

Definition at line 592 of file reconstructor.h.

References BaldwinWoolfordReconstructor().

00593                 {
00594                         return new BaldwinWoolfordReconstructor();
00595                 }

virtual TypeDict EMAN::BaldwinWoolfordReconstructor::get_param_types (  )  const [inline, virtual]

Get the parameter types of this object.

Returns:
a TypeDict detailing all of the acceptable (and necessary) parameters

Reimplemented from EMAN::FourierReconstructor.

Definition at line 597 of file reconstructor.h.

References EMAN::EMObject::BOOL, EMAN::EMObject::FLOAT, EMAN::EMObject::INT, EMAN::TypeDict::put(), and EMAN::EMObject::STRING.

00598                 {
00599                         TypeDict d;
00600                         d.put("mode", EMObject::INT, "Optional. Fourier pixel insertion mode [1-8] - mode 2 is default.");
00601                         d.put("x_in", EMObject::INT, "Necessary. The x dimension of the input images.");
00602                         d.put("y_in", EMObject::INT, "Necessary. The y dimension of the input images.");
00603                         d.put("zsample", EMObject::INT, "Optional. The z dimension (Fourier sampling) of the reconstructed volume, very useful for tomographic reconstruction. Works for general volumes.");
00604                         d.put("ysample", EMObject::INT, "Optional. The y dimension (Fourier sampling) of the reconstructed volume, works for general volumes. Not commonly specified.");
00605                         d.put("xsample", EMObject::INT, "Optional. The x dimension (Fourier sampling) of the reconstructed volume, works for general volumes. Not commonly specified.");
00606                         d.put("sym", EMObject::STRING, "Optional. The symmetry of the reconstructed volume, c?, d?, oct, tet, icos, h?. Default is c1");
00607                         d.put("maskwidth", EMObject::INT, "The width of the Fourier space kernel used to interpolate data to grid points" );
00608                         d.put("postmultiply", EMObject::BOOL, "A flag that controls whether or not the reconstructed volume is post multiplied in real space by the IFT of the weighting function. Default is on");
00609                         // Currently redundant
00610                         d.put("3damp", EMObject::BOOL, "this doesn't work, fixme dsaw");
00611                         d.put("hard", EMObject::FLOAT, "Optional. The quality metric threshold. Default is 0 (off).");
00612                         d.put("quiet", EMObject::BOOL, "Optional. Toggles writing useful information to standard out. Default is false.");
00613                         return d;
00614                 }

EMData * BaldwinWoolfordReconstructor::finish (  )  [virtual]

Finish reconstruction and return the complete model.

Returns:
The result 3D model.

Reimplemented from EMAN::FourierReconstructor.

Definition at line 933 of file reconstructor.cpp.

References EMAN::EMData::depad(), EMAN::EMData::do_ift_inplace(), EMAN::EMData::get_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::ReconstructorVolumeData::image, EMAN::FactoryBase::params, EMAN::EMData::process_inplace(), EMAN::Dict::set_default(), EMAN::ReconstructorVolumeData::tmp_data, EMAN::EMData::update(), weight, and EMAN::EMData::write_image().

00934 {
00935         tmp_data->write_image("density.mrc");
00936         image->process_inplace("xform.fourierorigin.tocorner");
00937         image->do_ift_inplace();
00938         image->depad();
00939         image->process_inplace("xform.phaseorigin.tocenter");
00940 
00941         if ( (bool) params.set_default("postmultiply", false) )
00942         {
00943                 cout << "POST MULTIPLYING" << endl;
00944         // now undo the Fourier convolution with real space division
00945                 float* d = image->get_data();
00946                 float N = (float) image->get_xsize()/2.0f;
00947                 N *= N;
00948                 size_t rnx = image->get_xsize();
00949                 size_t rny = image->get_ysize();
00950                 size_t rnxy = rnx*rny;
00951                 int cx = image->get_xsize()/2;
00952                 int cy = image->get_ysize()/2;
00953                 int cz = image->get_zsize()/2;
00954                 size_t idx;
00955                 for (int k = 0; k < image->get_zsize(); ++k ){
00956                         for (int j = 0; j < image->get_ysize(); ++j ) {
00957                                 for (int i =0; i < image->get_xsize(); ++i ) {
00958                                         float xd = (float)(i-cx); xd *= xd;
00959                                         float yd = (float)(j-cy); yd *= yd;
00960                                         float zd = (float)(k-cz); zd *= zd;
00961                                         float weight = exp((xd+yd+zd)/N);
00962                                         idx = k*rnxy + j*rnx + i;
00963                                         d[idx] *=  weight;
00964                                 }
00965                         }
00966                 }
00967         }
00968         image->update();
00969         return  image;
00970 }

int BaldwinWoolfordReconstructor::insert_slice_weights ( const Transform &  t3d  )  [virtual]

Reimplemented from EMAN::Reconstructor.

Definition at line 974 of file reconstructor.cpp.

References EMAN::Symmetry3D::get_symmetries(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::ReconstructorVolumeData::image, insert_density_at(), EMAN::EMData::is_fftodd(), EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::ny, EMAN::FactoryBase::params, EMAN::ReconstructorVolumeData::tmp_data, x, and y.

00975 {
00976         bool fftodd = image->is_fftodd();
00977         int rnx = nx-2*!fftodd;
00978 
00979         float y_scale = 1.0, x_scale = 1.0;
00980 
00981         if ( ny != rnx  )
00982         {
00983                 if ( rnx > ny ) y_scale = (float) rnx / (float) ny;
00984                 else x_scale = (float) ny / (float) rnx;
00985         }
00986 
00987         int tnx = tmp_data->get_xsize();
00988         int tny = tmp_data->get_ysize();
00989         int tnz = tmp_data->get_zsize();
00990 
00991         vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
00992         for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
00993                 Transform n3d = t3d*(*it);
00994 
00995                 for (int y = 0; y < tny; y++) {
00996                         for (int x = 0; x < tnx; x++) {
00997 
00998                                 float rx = (float) x;
00999                                 float ry = (float) y;
01000 
01001                                 if ( ny != rnx )
01002                                 {
01003                                         if ( rnx > ny ) ry *= y_scale;
01004                                         else rx *= x_scale;
01005                                 }
01006 //                              float xx = rx * n3d[0][0] + (ry - tny/2) * n3d[1][0];
01007 //                              float yy = rx * n3d[0][1] + (ry - tny/2) * n3d[1][1];
01008 //                              float zz = rx * n3d[0][2] + (ry - tny/2) * n3d[1][2];
01009 
01010                                 Vec3f coord(rx,(ry - tny/2),0);
01011                                 coord = coord*n3d; // transpose multiplication
01012                                 float xx = coord[0];
01013                                 float yy = coord[1];
01014                                 float zz = coord[2];
01015 
01016                                 if (xx < 0 ){
01017                                         xx = -xx;
01018                                         yy = -yy;
01019                                         zz = -zz;
01020                                 }
01021 
01022                                 yy += tny/2;
01023                                 zz += tnz/2;
01024                                 insert_density_at(xx,yy,zz);
01025                         }
01026                 }
01027         }
01028 
01029         return 0;
01030 }

virtual int EMAN::BaldwinWoolfordReconstructor::insert_slice ( const EMData *const   slice,
const Transform3D &  euler 
) [inline, virtual]

Insert an image slice to the reconstructor.

To insert multiple image slices, call this function multiple times.

Parameters:
slice Image slice.
euler Euler angle of this image slice.
Returns:
0 if OK. 1 if error.

Reimplemented from EMAN::FourierReconstructor.

Definition at line 622 of file reconstructor.h.

References UnexpectedBehaviorException.

00622                                                                                             {
00623                         throw UnexpectedBehaviorException("Transform3D interface is redundant");
00624                 }

int BaldwinWoolfordReconstructor::insert_slice ( const EMData *const   slice,
const Transform &  t3d 
) [virtual]

Insert a slice into a 3D volume, in a given orientation.

Returns:
0 if successful, 1 otherwise
Parameters:
slice the image slice to be inserted into the 3D volume
t3d the Transform that stores the image Euler angles
Exceptions:
NullPointerException if the input EMData pointer is null
ImageFormatException if the image is complex as opposed to real

Reimplemented from EMAN::FourierReconstructor.

Definition at line 1126 of file reconstructor.cpp.

References EMAN::EMData::do_fft_inplace(), dt, EMAN::EMData::get_attr(), EMAN::EMData::get_data(), EMAN::Transform::get_mirror(), EMAN::Transform::get_scale(), EMAN::Symmetry3D::get_symmetries(), EMAN::Transform::get_trans_2d(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::EMData::has_attr(), EMAN::ReconstructorVolumeData::image, insert_pixel(), EMAN::EMData::is_fftodd(), EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::ny, EMAN::FactoryBase::params, EMAN::EMData::process(), EMAN::EMData::process_inplace(), EMAN::Transform::set_mirror(), EMAN::Transform::set_rotation(), EMAN::Transform::set_scale(), EMAN::Transform::set_trans(), EMAN::ReconstructorVolumeData::tmp_data, x, and y.

01127 {
01128         Transform * rotation;
01129         if ( input_slice->has_attr("xform.projection") ) {
01130                 rotation = (Transform*) (input_slice->get_attr("xform.projection")); // assignment operator
01131         } else {
01132                 rotation = new Transform(t); // assignment operator
01133         }
01134         Transform tmp(*rotation);
01135         tmp.set_rotation(Dict("type","eman")); // resets the rotation to 0 implicitly
01136 
01137         Vec2f trans = tmp.get_trans_2d();
01138         float scale = tmp.get_scale();
01139         bool mirror = tmp.get_mirror();
01140         EMData* slice = 0;
01141         if (trans[0] != 0 || trans[1] != 0 || scale != 1.0 ) {
01142                 slice = input_slice->process("math.transform",Dict("transform",&tmp));
01143         } else if ( mirror == true ) {
01144                 slice = input_slice->process("xform.flip",Dict("axis","x"));
01145         }
01146         if ( slice == 0 ) {
01147                 slice = input_slice->process("xform.phaseorigin.tocorner");
01148         } else {
01149                 slice->process_inplace("xform.phaseorigin.tocorner");
01150         }
01151 
01152         slice->do_fft_inplace();
01153         slice->process_inplace("xform.fourierorigin.tocenter");
01154         float *dat = slice->get_data();
01155         float dt[2];
01156 
01157         bool fftodd = image->is_fftodd();
01158         int rnx = nx-2*!fftodd;
01159 
01160         float y_scale = 1.0, x_scale = 1.0;
01161 
01162         if ( ny != rnx  )
01163         {
01164                 if ( rnx > ny ) y_scale = (float) rnx / (float) ny;
01165                 else x_scale = (float) ny / (float) rnx;
01166         }
01167 
01168         int tnx = tmp_data->get_xsize();
01169         int tny = tmp_data->get_ysize();
01170         int tnz = tmp_data->get_zsize();
01171 
01172         vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
01173 //      float weight = params.set_default("weight",1.0f);
01174 
01175         rotation->set_scale(1.0); rotation->set_mirror(false); rotation->set_trans(0,0,0);
01176         for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
01177                 Transform t3d = (*rotation)*(*it);
01178 
01179                 for (int y = 0; y < tny; y++) {
01180                         for (int x = 0; x < tnx; x++) {
01181                                 float rx = (float) x;
01182                                 float ry = (float) y;
01183 
01184                                 if ( ny != rnx )
01185                                 {
01186                                         if ( rnx > ny ) ry *= y_scale;
01187                                         else rx *= x_scale;
01188                                 }
01189 
01190 //                              float xx = rx * n3d[0][0] + (ry - tny/2) * n3d[1][0];
01191 //                              float yy = rx * n3d[0][1] + (ry - tny/2) * n3d[1][1];
01192 //                              float zz = rx * n3d[0][2] + (ry - tny/2) * n3d[1][2];
01193 
01194                                 Vec3f coord(rx,(ry - tny/2),0);
01195                                 coord = coord*t3d; // transpose multiplication
01196                                 float xx = coord[0];
01197                                 float yy = coord[1];
01198                                 float zz = coord[2];
01199 
01200 
01201                                 float cc = 1;
01202                                 if (xx < 0 ){
01203                                         xx = -xx;
01204                                         yy = -yy;
01205                                         zz = -zz;
01206                                         cc = -1;
01207                                 }
01208 
01209                                 yy += tny/2;
01210                                 zz += tnz/2;
01211 
01212                                 int idx = x * 2 + y * (slice->get_xsize());
01213                                 dt[0] = dat[idx];
01214                                 dt[1] = cc * dat[idx+1];
01215 
01216                                 insert_pixel(xx,yy,zz,dt);
01217                         }
01218                 }
01219         }
01220 
01221         if(rotation) {delete rotation; rotation=0;}
01222         delete slice;
01223 
01224         return 0;
01225 }

void BaldwinWoolfordReconstructor::insert_pixel ( const float &  x,
const float &  y,
const float &  z,
const float  dt[2] 
)

Definition at line 1227 of file reconstructor.cpp.

References dfreq, dist(), EMAN::Util::fast_floor(), EMAN::EMData::get_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::ReconstructorVolumeData::image, InvalidValueException, OutofRangeException, EMAN::FactoryBase::params, EMAN::Dict::set_default(), EMAN::ReconstructorVolumeData::tmp_data, W, and weight.

Referenced by insert_slice().

01228 {
01229         int xl = Util::fast_floor(x);
01230         int yl = Util::fast_floor(y);
01231         int zl = Util::fast_floor(z);
01232 
01233         // w is the windowing width
01234         int w = params.set_default("maskwidth",2);
01235         float wsquared = (float) w*w;
01236         float dw = 1.0f/w;
01237         dw *= dw;
01238 
01239         int wmox = w-1;
01240         int wmoy = w-1;
01241         int wmoz = w-1;
01242 
01243         // If any coordinate is incedental with a vertex, then
01244         // make sure there is symmetry in density accruing.
01245         // i.e. the window width must be equal in both directions
01246         if ( ((float) xl) == x ) wmox = w;
01247         if ( ((float) yl) == y ) wmoy = w;
01248         if ( ((float) zl) == z ) wmoz = w;
01249 
01250         float* we = tmp_data->get_data();
01251         int tnx = tmp_data->get_xsize();
01252         int tny = tmp_data->get_ysize();
01253         int tnz = tmp_data->get_zsize();
01254         int tnxy = tnx*tny;
01255 
01256         int rnx = 2*tnx;
01257         int rnxy = 2*tnxy;
01258 
01259         int mode = params.set_default("mode",1);
01260 
01261         float* d = image->get_data();
01262         for(int k = zl-wmoz; k <= zl+w; ++k ) {
01263                 for(int j = yl-wmoy; j <= yl+w; ++j) {
01264                         for( int i = xl-wmox; i <= xl+w; ++i) {
01265                                 float fac = 1.0;
01266                                 int ic = i, jc = j, kc = k;
01267 
01268                                 // Fourier space is periodic, which is enforced
01269                                 // by the next 6 if statements. These if statements
01270                                 // assume that the Fourier DC component is at
01271                                 // (0,ny/2,nz/2).
01272                                 float negfac=1.0;
01273                                 if ( i <= 0 ) {
01274                                         if ( x != 0 && i == 0 ) fac = 1.0;
01275                                         else if ( x == 0 && i < 0) continue;
01276                                         if (i < 0 ) {
01277                                                 continue;
01278                                                 ic = -i;
01279                                                 jc = tny-jc;
01280                                                 kc = tnz-kc;
01281                                                 negfac=-1.0;
01282                                         }
01283                                 }
01284                                 if ( ic >= tnx ) ic = 2*tnx-ic-1;
01285                                 if ( jc < 0 ) jc = tny+jc;
01286                                 if ( jc >= tny ) jc = jc-tny;
01287                                 if ( kc < 0 ) kc = tnz+kc;
01288                                 if ( kc >= tnz ) kc = kc-tnz;
01289 //                              if ( ic >= tnx ) continue;
01290 //                              if ( jc < 0 ) continue;
01291 //                              if ( jc >= tny ) continue;
01292 //                              if ( kc < 0 ) continue;
01293 //                              if ( kc >= tnz ) continue;
01294 
01295                                 float zd = (z-(float)k);
01296                                 float yd = (y-(float)j);
01297                                 float xd = (x-(float)i);
01298                                 zd *= zd; yd *= yd; xd *= xd;
01299                                 float distsquared = xd+yd+zd;
01300 //                              float f = fac*exp(-dw*(distsquared));
01301                                 float f = fac*exp(-2.467f*(distsquared));
01302                                 float weight = f/we[kc*tnxy+jc*tnx+ic];
01303                                 // debug - this error should never occur
01304                                 if ( (kc*rnxy+jc*rnx+2*ic+1) >= rnxy*tnz ) throw OutofRangeException(0,rnxy*tnz,kc*rnxy+jc*rnx+2*ic+1, "in pixel insertion" );
01305                                 size_t k = kc*rnxy+jc*rnx+2*ic;
01306 
01307                                 float factor, dist,residual;
01308                                 int sizeW,sizeWmid,idx;
01309                                 switch (mode) {
01310                                         case 0:
01311                                                 d[k] += weight*f*dt[0];
01312                                                 d[k+1] += negfac*weight*f*dt[1];
01313                                                 cout << "hello" << endl;
01314                                         break;
01315 
01316                                         case 1:
01317                                                 // We enforce a spherical kernel
01318                                                 if ( distsquared > wsquared ) continue;
01319 
01320                                                 sizeW = (int)(1+2*w/dfreq);
01321                                                 sizeWmid = sizeW/2;
01322 
01323                                                 dist = sqrtf(distsquared);
01324                                                 idx = (int)(sizeWmid + dist/dfreq);
01325                                                 if (idx >= sizeW) throw InvalidValueException(idx, "idx was greater than or equal to sizeW");
01326                                                 residual = dist/dfreq - (int)(dist/dfreq);
01327                                                 if ( fabs(residual) > 1) throw InvalidValueException(residual, "Residual was too big");
01328 
01329                                                 factor = (W[idx]*(1.0f-residual)+W[idx+1]*residual)*weight;
01330 
01331                                                 d[k] += dt[0]*factor;
01332                                                 d[k+1] += dt[1]*factor;
01333                                         break;
01334 
01335                                         default:
01336                                                 throw InvalidValueException(mode, "The mode was unsupported in BaldwinWoolfordReconstructor::insert_pixel");
01337                                         break;
01338                                 }
01339                         }
01340                 }
01341         }
01342 }

void BaldwinWoolfordReconstructor::insert_density_at ( const float &  x,
const float &  y,
const float &  z 
)

Definition at line 1032 of file reconstructor.cpp.

References EMAN::Util::fast_floor(), EMAN::EMData::get_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), OutofRangeException, EMAN::FactoryBase::params, EMAN::Dict::set_default(), and EMAN::ReconstructorVolumeData::tmp_data.

Referenced by insert_slice_weights().

01033 {
01034         int xl = Util::fast_floor(x);
01035         int yl = Util::fast_floor(y);
01036         int zl = Util::fast_floor(z);
01037 
01038         // w is the windowing width
01039         int w = params.set_default("maskwidth",2);
01040         float wsquared = (float) w*w;
01041         float dw = 1.0f/w;
01042         dw *= dw;
01043 
01044         // w minus one - this control the number of
01045         // pixels/voxels to the left of the main pixel
01046         // that will have density
01047         int wmox = w-1;
01048         int wmoy = w-1;
01049         int wmoz = w-1;
01050 
01051         // If any coordinate is incedental with a vertex, then
01052         // make sure there is symmetry in density accruing.
01053         // i.e. the window width must be equal in both directions
01054         if ( ((float) xl) == x ) wmox = w;
01055         if ( ((float) yl) == y ) wmoy = w;
01056         if ( ((float) zl) == z ) wmoz = w;
01057 
01058         float* d = tmp_data->get_data();
01059         int tnx = tmp_data->get_xsize();
01060         int tny = tmp_data->get_ysize();
01061         int tnz = tmp_data->get_zsize();
01062         size_t tnxy = tnx*tny;
01063 
01064         int mode = params.set_default("mode",1);
01065 
01066         for(int k = zl-wmoz; k <= zl+w; ++k ) {
01067                 for(int j = yl-wmoy; j <= yl+w; ++j) {
01068                         for( int i = xl-wmox; i <= xl+w; ++i) {
01069                                 float fac = 1.0;
01070                                 int ic = i, jc = j, kc = k;
01071 
01072                                 // Fourier space is periodic, which is enforced
01073                                 // by the next 6 if statements. These if statements
01074                                 // assume that the Fourier DC components is at
01075                                 // (0,ny/2,nz/2).
01076                                 if ( i <= 0 ) {
01077 
01078                                         if ( x != 0 && i == 0 ) fac = 1.0;
01079                                         else if ( x == 0 && i < 0) continue;
01080 //                                      if (i < 0 ) ic = -i;
01081                                         if (i < 0 ) {
01082                                                 continue;
01083                                                 ic = -i;
01084                                                 jc = tny-jc;
01085                                                 kc = tnz-kc;
01086                                         }
01087                                 }
01088                                 if ( ic >= tnx ) ic = 2*tnx-ic-1;
01089                                 if ( jc < 0 ) jc = tny+jc;
01090                                 if ( jc >= tny ) jc = jc-tny;
01091                                 if ( kc < 0 ) kc = tnz+kc;
01092                                 if ( kc >= tnz ) kc = kc-tnz;
01093 //                              if ( ic >= tnx ) continue;
01094 //                              if ( jc < 0 ) continue;
01095 //                              if ( jc >= tny ) continue;
01096 //                              if ( kc < 0 ) continue;
01097 //                              if ( kc >= tnz ) continue;
01098                                 // This shouldn't happen
01099                                 // Debug remove later
01100                                 if ( ic < 0 ) { cout << "wo 1" << endl; }
01101                                 if ( ic >= tnx  ){ cout << "wo 2" << endl; }
01102                                 if ( jc < 0 ) { cout << "wo 3" << endl; }
01103                                 if ( jc >= tny ) { cout << "wo 4" << endl; }
01104                                 if ( kc < 0 ) { cout << "wo 5" << endl; }
01105                                 if ( kc >= tnz ) { cout << "wo 6" << endl; }
01106 
01107 
01108                                 float zd = (z-(float)k);
01109                                 float yd = (y-(float)j);
01110                                 float xd = (x-(float)i);
01111                                 zd *= zd; yd *= yd; xd *= xd;
01112                                 float distsquared = xd+yd+zd;
01113                                 // We enforce a spherical kernel
01114                                 if ( mode == 1 && distsquared > wsquared ) continue;
01115 
01116 //                              float f = fac*exp(-dw*(distsquared));
01117                                 float f = fac*exp(-2.467f*(distsquared));
01118                                 // Debug - this error should never occur.
01119                                 if ( (kc*tnxy+jc*tnx+ic) >= tnxy*tnz ) throw OutofRangeException(0,tnxy*tnz,kc*tnxy+jc*tnx+ic, "in density insertion" );
01120                                 d[kc*tnxy+jc*tnx+ic] += f;
01121                         }
01122                 }
01123         }
01124 }

void BaldwinWoolfordReconstructor::setup (  )  [virtual]

Setup the Fourier reconstructor relies on the correct parameters throws a variety of exceptions if the parameters are unworkable

Exceptions:
InvalidValueException when the x_in parameter is odd or less than zero
InvalidValueException when the y_in parameter is odd or less than zero
InvalidValueException when the x_out parameter is odd or less than zero
InvalidValueException when the y_out parameter is odd or less than zero
InvalidValueException when the z_out parameter is odd or less than zero
InvalidValueException when the pad parameter is less than the greatest dimension of the input images (the greater of x_in and y_in) Note the restraint that the Fourier volumes should be even will be lifted once debugging of the the xform.fourierorigin processor is performed, but however, when this happens a rigorous testing should be performed because no testing has been done in this regard.

Reimplemented from EMAN::FourierReconstructor.

Definition at line 919 of file reconstructor.cpp.

References dfreq, EMAN::Util::getBaldwinGridWeights(), EMAN::FourierReconstructor::max_input_dim, EMAN::FactoryBase::params, EMAN::Dict::set_default(), EMAN::FourierReconstructor::setup(), and W.

00920 {
00921         //This is a bit of a hack - but for now it suffices
00922         params.set_default("mode",1);
00923         FourierReconstructor::setup();
00924         // Set up the Baldwin Kernel
00925         int P = (int)((1.0+0.25)*max_input_dim+1);
00926         float r = (float)(max_input_dim+1)/(float)P;
00927         dfreq = 0.2f;
00928         if (W != 0) delete [] W;
00929         int maskwidth = params.set_default("maskwidth",2);
00930         W = Util::getBaldwinGridWeights(maskwidth, (float)P, r,dfreq,0.5f,0.2f);
00931 }

void EMAN::BaldwinWoolfordReconstructor::load_default_settings (  )  [inline, protected]

Load default settings.

Reimplemented from EMAN::FourierReconstructor.

Definition at line 635 of file reconstructor.h.

References EMAN::FactoryBase::params.

00636                 {
00637                         params["mode"] = 1;
00638                         params["x_in"] = 0;
00639                         params["y_in"] = 0;
00640                         params["zsample"] = 0;
00641                         params["ysample"] = 0;
00642                         params["xsample"] = 0;
00643                         params["sym"] = "c1";
00644                         params["maskwidth"] = 3;
00645 
00646                         // Currently redundant
00647                         params["3damp"] = false;
00648                         params["hard"] = 0.05;
00649                         params["quiet"] = false;
00650                 }

BaldwinWoolfordReconstructor& EMAN::BaldwinWoolfordReconstructor::operator= ( const BaldwinWoolfordReconstructor &   )  [private]

Disallow assignment.

Member Data Documentation

float* EMAN::BaldwinWoolfordReconstructor::W [private]

Definition at line 660 of file reconstructor.h.

Referenced by insert_pixel(), setup(), and ~BaldwinWoolfordReconstructor().

float EMAN::BaldwinWoolfordReconstructor::dfreq [private]

Definition at line 661 of file reconstructor.h.

Referenced by insert_pixel(), and setup().

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