EMAN::PawelProjector Class Reference

Pawel Penczek's optimized projection routine. More...

#include <projector.h>

Inheritance diagram for EMAN::PawelProjector:

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

Collaboration graph
[legend]
List of all members.

Public Member Functions

EMDataproject3d (EMData *image) const
 Project an 3D image into a 2D image.
EMDatabackproject3d (EMData *image) const
 Back-project a 2D image into a 3D image.
string get_name () const
 Get the projector's name.
string get_desc () const
TypeDict get_param_types () const
 Get processor parameter information in a dictionary.

Static Public Member Functions

static ProjectorNEW ()

Static Public Attributes

static const string NAME = "pawel"

Private Member Functions

void prepcubes (int nx, int ny, int nz, int ri, Vec3i origin, int &nn, IPCube *ipcube=NULL) const

Classes

struct  IPCube

Detailed Description

Pawel Penczek's optimized projection routine.

Definition at line 248 of file projector.h.


Member Function Documentation

EMData * PawelProjector::backproject3d ( EMData image  )  const [virtual]

Back-project a 2D image into a 3D image.

Returns:
A 3D image from the backprojection.

Implements EMAN::Projector.

Definition at line 1972 of file projector.cpp.

References anglelist, EMAN::Transform::at(), EMDeleteArray(), EMAN::PawelProjector::IPCube::end, EMAN::EMData::get_data(), EMAN::Util::get_min(), EMAN::Transform::get_rotation(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::Dict::has_key(), images, EMAN::PawelProjector::IPCube::loc, LOGERR, nn(), NullPointerException, nx, ny, EMAN::Projector::params, phi, prepcubes(), EMAN::EMData::set_size(), EMAN::Dict::size(), EMAN::PawelProjector::IPCube::start, theta, EMAN::EMData::to_zero(), and EMAN::EMData::update().

01973 {
01974 
01975         float *images;
01976 
01977         if (!imagestack) {
01978                 return 0;
01979         }
01980         int ri;
01981         int nx      = imagestack->get_xsize();
01982         int ny      = imagestack->get_ysize();
01983 //     int nslices = imagestack->get_zsize();
01984         int dim = Util::get_min(nx,ny);
01985         images  = imagestack->get_data();
01986 
01987         Vec3i origin(0,0,0);
01988     // If a sensible origin isn't passed in, choose the middle of
01989     // the cube.
01990         if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
01991         else {origin[0] = nx/2;}
01992         if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
01993         else {origin[1] = ny/2;}
01994         if (params.has_key("origin_z")) {origin[1] = params["origin_z"];}
01995         else {origin[2] = dim/2;}
01996 
01997         if (params.has_key("radius")) {ri = params["radius"];}
01998         else {ri = dim/2 - 1;}
01999 
02000     // Determine the number of rows (x-lines) within the radius
02001         int nn = -1;
02002         prepcubes(nx, ny, dim, ri, origin, nn);
02003     // nn is now the number of rows-1 within the radius
02004     // so we can create and fill the ipcubes
02005         IPCube* ipcube = new IPCube[nn+1];
02006         prepcubes(nx, ny, dim, ri, origin, nn, ipcube);
02007 
02008         int nangles = 0;
02009         vector<float> anglelist;
02010         // Do we have a list of angles?
02011         if (params.has_key("anglelist")) {
02012                 anglelist = params["anglelist"];
02013                 nangles = anglelist.size() / 3;
02014         } else {
02015                 Transform* t3d = params["transform"];
02016                 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
02017                 // This part was modified by David Woolford -
02018                 // Before this the code worked only for SPIDER and EMAN angles,
02019                 // but the framework of the Transform3D allows for a generic implementation
02020                 // as specified here.
02021                 Dict p = t3d->get_rotation("spider");
02022                 if(t3d) {delete t3d; t3d=0;}
02023 
02024                 string angletype = "SPIDER";
02025                 float phi = p["phi"];
02026                 float theta = p["theta"];
02027                 float psi = p["psi"];
02028                 anglelist.push_back(phi);
02029                 anglelist.push_back(theta);
02030                 anglelist.push_back(psi);
02031                 nangles = 1;
02032         }
02033 
02034         // End David Woolford modifications
02035 
02036     // initialize return object
02037         EMData* ret = new EMData();
02038         ret->set_size(nx, ny, dim);
02039         ret->to_zero();
02040 
02041     // loop over sets of angles
02042         for (int ia = 0; ia < nangles; ia++) {
02043                 int indx = 3*ia;
02044                 Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]);
02045                 Transform rotation(d);
02046                 float dm1 = rotation.at(0,0);
02047                 float dm4 = rotation.at(1,0);
02048 
02049                 if (2*(ri+1)+1 > dim) {
02050           // Must check x and y boundaries
02051                         LOGERR("backproject3d, pawel, 2*(ri+1)+1 > dim\n");
02052                         return 0;
02053                 } else {
02054           // No need to check x and y boundaries
02055                         for (int i = 0 ; i <= nn; i++) {
02056                                 int iox = (int)ipcube[i].loc[0]+origin[0];
02057                                 int ioy = (int)ipcube[i].loc[1]+origin[1];
02058                                 int ioz = (int)ipcube[i].loc[2]+origin[2];
02059 
02060                                 Vec3f vb = rotation*ipcube[i].loc + origin;
02061                                 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
02062                                         float xbb = (j-ipcube[i].start)*dm1 + vb[0];
02063                                         int   iqx = (int)floor(xbb);
02064 
02065                                         float ybb = (j-ipcube[i].start)*dm4 + vb[1];
02066                                         int   iqy = (int)floor(ybb);
02067 
02068                                         float dipx = xbb - iqx;
02069                                         float dipy = ybb - iqy;
02070 
02071                                         (*ret)(iox,ioy,ioz) += images(iqx,iqy,ia)
02072                                                         + dipy*(images(iqx,iqy+1,ia)-images(iqx,iqy,ia))
02073                                                         + dipx*(images(iqx+1,iqy,ia)-images(iqx,iqy,ia)
02074                                                                         + dipy*(images(iqx+1,iqy+1,ia)-images(iqx+1,iqy,ia)
02075                                                                                         - images(iqx,iqy+1,ia)+images(iqx,iqy,ia)));
02076                                         iox++;
02077                                 } // end for j
02078                         } // end for i
02079                 } // end if
02080         } // end for ia
02081 
02082         ret->update();
02083         EMDeleteArray(ipcube);
02084         return ret;
02085 }

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

Implements EMAN::Projector.

Definition at line 259 of file projector.h.

00260                 {
00261                         return "Pawel Penczek's optimized real-space projection generation. Minor interpolation artifacts.";
00262                 }

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

Get the projector's name.

Each projector is indentified by unique name.

Returns:
The projector's name.

Implements EMAN::Projector.

Definition at line 254 of file projector.h.

References NAME.

00255                 {
00256                         return NAME;
00257                 }

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

Get processor parameter information in a dictionary.

Each parameter has one record in the dictionary. Each record contains its name, data-type, and description.

Returns:
A dictionary containing the parameter info.

Reimplemented from EMAN::Projector.

Definition at line 269 of file projector.h.

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

00270                 {
00271                         TypeDict d;
00272                         d.put("transform", EMObject::TRANSFORM);
00273                         d.put("origin_x", EMObject::INT);
00274                         d.put("origin_y", EMObject::INT);
00275                         d.put("origin_z", EMObject::INT);
00276                         d.put("radius", EMObject::INT);
00277                         d.put("anglelist", EMObject::FLOATARRAY);
00278                         d.put("angletype", EMObject::STRING);
00279                         d.put("theta", EMObject::FLOAT);
00280                         d.put("psi", EMObject::FLOAT);
00281                         return d;
00282                 }

static Projector* EMAN::PawelProjector::NEW (  )  [inline, static]

Definition at line 264 of file projector.h.

00265                 {
00266                         return new PawelProjector();
00267                 }

void PawelProjector::prepcubes ( int  nx,
int  ny,
int  nz,
int  ri,
Vec3i  origin,
int &  nn,
IPCube ipcube = NULL 
) const [private]

Definition at line 531 of file projector.cpp.

References EMAN::PawelProjector::IPCube::end, EMAN::PawelProjector::IPCube::loc, and t.

Referenced by backproject3d(), and project3d().

00532                                                                       {
00533         const float r = float(ri)*float(ri);
00534         const int ldpx = origin[0];
00535         const int ldpy = origin[1];
00536         const int ldpz = origin[2];
00537         //cout<<"  ldpx  "<<ldpx<<"  i2  "<<ldpy<<"  i1 "<<ldpz<<endl;
00538         float t;
00539         nn = -1;
00540         for (int i1 = 0; i1 < nz; i1++) {
00541                 t = float(i1 - ldpz);
00542                 const float xx = t*t;
00543                 for (int i2 = 0; i2 < ny; i2++) {
00544                         t = float(i2 - ldpy);
00545                         const float yy = t*t + xx;
00546                         bool first = true;
00547                         for (int i3 = 0; i3 < nx; i3++) {
00548                                 t = float(i3 - ldpx);
00549                                 const float rc = t*t + yy;
00550                                 if (first) {
00551                                         // first pixel on this line
00552                                         if (rc > r) continue;
00553                                         first = false;
00554                                         nn++;
00555                                         if (ipcube != NULL) {
00556                                                 ipcube[nn].start = i3;
00557                                                 ipcube[nn].end   = i3;
00558                                                 ipcube[nn].loc[0] = i3 - ldpx;
00559                                                 ipcube[nn].loc[1] = i2 - ldpy;
00560                                                 ipcube[nn].loc[2] = i1 - ldpz;
00561                                                 //cout<<"  start  "<<i3<<"  i2  "<<i2<<"  i1 "<<i1<<endl;
00562                                         }
00563                                 } else {
00564                                         // second or later pixel on this line
00565                                         if (ipcube != NULL) {
00566                                                 if (rc <= r) ipcube[nn].end = i3;
00567                                         }
00568                                 }
00569                         }
00570                 }
00571         }
00572 }

EMData * PawelProjector::project3d ( EMData image  )  const [virtual]

Project an 3D image into a 2D image.

Returns:
A 2D image from the projection.

Implements EMAN::Projector.

Definition at line 575 of file projector.cpp.

References anglelist, EMDeleteArray(), EMAN::PawelProjector::IPCube::end, EMAN::Transform::get_matrix3_row(), EMAN::Util::get_max(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::Dict::has_key(), EMAN::PawelProjector::IPCube::loc, LOGERR, nn(), NullPointerException, nx, ny, EMAN::Projector::params, prepcubes(), EMAN::EMData::set_size(), EMAN::EMData::to_zero(), and EMAN::EMData::update().

00576 {
00577         if (!image) {
00578                 return 0;
00579         }
00580         int ri;
00581         int nx = image->get_xsize();
00582         int ny = image->get_ysize();
00583         int nz = image->get_zsize();
00584         int dim = Util::get_max(nx,ny,nz);
00585         if (nz == 1) {
00586                 LOGERR("The PawelProjector needs a volume!");
00587                 return 0;
00588         }
00589 
00590         Vec3i origin(0,0,0);
00591         // If a sensible origin isn't passed in, choose the middle of
00592         // the cube.
00593         if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
00594         else {origin[0] = nx/2;}
00595         if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
00596         else {origin[1] = ny/2;}
00597         if (params.has_key("origin_z")) {origin[2] = params["origin_z"];}
00598         else {origin[2] = nz/2;}
00599 
00600         if (params.has_key("radius")) {ri = params["radius"];}
00601         else {ri = dim/2 - 1;}
00602 
00603         // Determine the number of rows (x-lines) within the radius
00604         int nn = -1;
00605         prepcubes(nx, ny, nz, ri, origin, nn);
00606         // nn is now the number of rows-1 within the radius
00607         // so we can create and fill the ipcubes
00608         IPCube* ipcube = new IPCube[nn+1];
00609         prepcubes(nx, ny, nz, ri, origin, nn, ipcube);
00610 
00611         Transform* rotation = params["transform"];
00612         int nangles = 0;
00613         vector<float> anglelist;
00614         // Do we have a list of angles?
00615         /*
00616         if (params.has_key("anglelist")) {
00617                 anglelist = params["anglelist"];
00618                 nangles = anglelist.size() / 3;
00619         } else {*/
00620 
00621                 if ( rotation == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
00622                 /*
00623                 Dict p = t3d->get_rotation("spider");
00624 
00625                 string angletype = "SPIDER";
00626                 float phi = p["phi"];
00627                 float theta = p["theta"];
00628                 float psi = p["psi"];
00629                 anglelist.push_back(phi);
00630                 anglelist.push_back(theta);
00631                 anglelist.push_back(psi);
00632                 */
00633                 nangles = 1;
00634         //}
00635 
00636         // initialize return object
00637         EMData* ret = new EMData();
00638         ret->set_size(nx, ny, nangles);
00639         ret->to_zero();
00640         //for (int i = 0 ; i <= nn; i++)  cout<<" IPCUBE "<<ipcube[i].start<<"  "<<ipcube[i].end<<"  "<<ipcube[i].loc[0]<<"  "<<ipcube[i].loc[1]<<"  "<<ipcube[i].loc[2]<<endl;
00641         // loop over sets of angles
00642         for (int ia = 0; ia < nangles; ia++) {
00643                 //int indx = 3*ia;
00644                 //Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]);
00645                 //Transform rotation(d);
00646                 if (2*(ri+1)+1 > dim) {
00647                         // Must check x and y boundaries
00648                         for (int i = 0 ; i <= nn; i++) {
00649                                 int k = ipcube[i].loc[1] + origin[1];
00650                                 Vec3f vb = ipcube[i].loc*(*rotation) + origin;
00651                                 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
00652                                         // check for pixels out-of-bounds
00653                                         int iox = int(vb[0]);
00654                                         if ((iox >= 0) && (iox < nx-1)) {
00655                                                 int ioy = int(vb[1]);
00656                                                 if ((ioy >= 0) && (ioy < ny-1)) {
00657                                                         int ioz = int(vb[2]);
00658                                                         if ((ioz >= 0) && (ioz < nz-1)) {
00659                                                                 // real work for pixels in bounds
00660                                         //cout<<"  iox  "<<iox<<"  ioy  "<<int(vb[1])<<"  ioz "<<int(vb[2])<<endl;
00661                                         //cout<<"  TAKE"<<endl;
00662                                                                 float dx = vb[0] - iox;
00663                                                                 float dy = vb[1] - ioy;
00664                                                                 float dz = vb[2] - ioz;
00665                                                                 float a1 = (*image)(iox,ioy,ioz);
00666                                                                 float a2 = (*image)(iox+1,ioy,ioz) - a1;
00667                                                                 float a3 = (*image)(iox,ioy+1,ioz) - a1;
00668                                                                 float a4 = (*image)(iox,ioy,ioz+1) - a1;
00669                                                                 float a5 = -a2 -(*image)(iox,ioy+1,ioz)
00670                                                                                 + (*image)(iox+1,ioy+1,ioz);
00671                                                                 float a61 = -(*image)(iox,ioy,ioz+1)
00672                                                                                 + (*image)(iox+1,ioy,ioz+1);
00673                                                                 float a6 = -a2 + a61;
00674                                                                 float a7 = -a3 - (*image)(iox,ioy,ioz+1)
00675                                                                                 + (*image)(iox,ioy+1,ioz+1);
00676                                                                 float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
00677                                                                                 + (*image)(iox+1,ioy+1,ioz+1);
00678                                                                 (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
00679                                                                                 + (a7 + a8*dx)*dy)
00680                                                                                 + a3*dy + dx*(a2 + a5*dy);
00681                                                         } //else {cout<<"  iox  "<<iox<<"  ioy  "<<int(vb[1])<<"  ioz "<<int(vb[2])<<"  VIOLAATED Z "<<endl; }
00682                                                 }//else {cout<<"  iox  "<<iox<<"  ioy  "<<int(vb[1])<<"  ioz "<<int(vb[2])<<"  VIOLATED Y "<<endl; }
00683                                         }//else {cout<<"  iox  "<<iox<<"  ioy  "<<int(vb[1])<<"  ioz "<<int(vb[2])<<"  VIOLATED X "<<endl; }
00684                                         vb += rotation->get_matrix3_row(0);
00685                                 }
00686                         }
00687 
00688                 } else {
00689                         // No need to check x and y boundaries
00690                         for (int i = 0 ; i <= nn; i++) {
00691                                 int k = ipcube[i].loc[1] + origin[1];
00692                                 Vec3f vb = ipcube[i].loc*(*rotation) + origin;
00693                                 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
00694                                         int iox = int(vb[0]);
00695                                         int ioy = int(vb[1]);
00696                                         int ioz = int(vb[2]);
00697                                         float dx = vb[0] - iox;
00698                                         float dy = vb[1] - ioy;
00699                                         float dz = vb[2] - ioz;
00700                                         float a1 = (*image)(iox,ioy,ioz);
00701                                         float a2 = (*image)(iox+1,ioy,ioz) - a1;
00702                                         float a3 = (*image)(iox,ioy+1,ioz) - a1;
00703                                         float a4 = (*image)(iox,ioy,ioz+1) - a1;
00704                                         float a5 = -a2 -(*image)(iox,ioy+1,ioz)
00705                                                         + (*image)(iox+1,ioy+1,ioz);
00706                                         float a61 = -(*image)(iox,ioy,ioz+1)
00707                                                         + (*image)(iox+1,ioy,ioz+1);
00708                                         float a6 = -a2 + a61;
00709                                         float a7 = -a3 - (*image)(iox,ioy,ioz+1)
00710                                                         + (*image)(iox,ioy+1,ioz+1);
00711                                         float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
00712                                                         + (*image)(iox+1,ioy+1,ioz+1);
00713                                         (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
00714                                                         + (a7 + a8*dx)*dy)
00715                                                         + a3*dy + dx*(a2 + a5*dy);
00716                                         vb += rotation->get_matrix3_row(0);
00717                                 }
00718                         }
00719                 }
00720         }
00721         ret->update();
00722         EMDeleteArray(ipcube);
00723         if(rotation) {delete rotation; rotation=0;}
00724         return ret;
00725 }


Member Data Documentation

const string PawelProjector::NAME = "pawel" [static]

Definition at line 284 of file projector.h.

Referenced by get_name().


The documentation for this class was generated from the following files:
Generated on Thu May 3 10:11:06 2012 for EMAN2 by  doxygen 1.4.7