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 1964 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().

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

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

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 723 of file projector.cpp.

References anglelist, EMDeleteArray(), EMAN::PawelProjector::IPCube::end, EMAN::Transform::get_matrix3_row(), EMAN::Util::get_min(), 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().

00724 {
00725         if (!image) {
00726                 return 0;
00727         }
00728         int ri;
00729         int nx = image->get_xsize();
00730         int ny = image->get_ysize();
00731         int nz = image->get_zsize();
00732         int dim = Util::get_min(nx,ny,nz);
00733         if (nz == 1) {
00734                 LOGERR("The PawelProjector needs a volume!");
00735                 return 0;
00736         }
00737         Vec3i origin(0,0,0);
00738         // If a sensible origin isn't passed in, choose the middle of
00739         // the cube.
00740         if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
00741         else {origin[0] = nx/2;}
00742         if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
00743         else {origin[1] = ny/2;}
00744         if (params.has_key("origin_z")) {origin[2] = params["origin_z"];}
00745         else {origin[2] = nz/2;}
00746 
00747         if (params.has_key("radius")) {ri = params["radius"];}
00748         else {ri = dim/2 - 1;}
00749 
00750         // Determine the number of rows (x-lines) within the radius
00751         int nn = -1;
00752         prepcubes(nx, ny, nz, ri, origin, nn);
00753         // nn is now the number of rows-1 within the radius
00754         // so we can create and fill the ipcubes
00755         IPCube* ipcube = new IPCube[nn+1];
00756         prepcubes(nx, ny, nz, ri, origin, nn, ipcube);
00757 
00758         Transform* rotation = params["transform"];
00759         int nangles = 0;
00760         vector<float> anglelist;
00761         // Do we have a list of angles?
00762         /*
00763         if (params.has_key("anglelist")) {
00764                 anglelist = params["anglelist"];
00765                 nangles = anglelist.size() / 3;
00766         } else {*/
00767 
00768                 if ( rotation == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
00769                 /*
00770                 Dict p = t3d->get_rotation("spider");
00771 
00772                 string angletype = "SPIDER";
00773                 float phi = p["phi"];
00774                 float theta = p["theta"];
00775                 float psi = p["psi"];
00776                 anglelist.push_back(phi);
00777                 anglelist.push_back(theta);
00778                 anglelist.push_back(psi);
00779                 */
00780                 nangles = 1;
00781         //}
00782 
00783         // initialize return object
00784         EMData* ret = new EMData();
00785         ret->set_size(nx, ny, nangles);
00786         ret->to_zero();
00787 
00788         // loop over sets of angles
00789         for (int ia = 0; ia < nangles; ia++) {
00790                 //int indx = 3*ia;
00791                 //Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]);
00792                 //Transform rotation(d);
00793                 if (2*(ri+1)+1 > dim) {
00794                         // Must check x and y boundaries
00795                         for (int i = 0 ; i <= nn; i++) {
00796                                 int k = ipcube[i].loc[1] + origin[1];
00797                                 Vec3f vb = ipcube[i].loc*(*rotation) + origin;
00798                                 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
00799                                         // check for pixels out-of-bounds
00800                                         int iox = int(vb[0]);
00801                                         if ((iox >= 0) && (iox < nx-1)) {
00802                                                 int ioy = int(vb[1]);
00803                                                 if ((ioy >= 0) && (ioy < ny-1)) {
00804                                                         int ioz = int(vb[2]);
00805                                                         if ((ioz >= 0) && (ioz < nz-1)) {
00806                                                                 // real work for pixels in bounds
00807                                                                 float dx = vb[0] - iox;
00808                                                                 float dy = vb[1] - ioy;
00809                                                                 float dz = vb[2] - ioz;
00810                                                                 float a1 = (*image)(iox,ioy,ioz);
00811                                                                 float a2 = (*image)(iox+1,ioy,ioz) - a1;
00812                                                                 float a3 = (*image)(iox,ioy+1,ioz) - a1;
00813                                                                 float a4 = (*image)(iox,ioy,ioz+1) - a1;
00814                                                                 float a5 = -a2 -(*image)(iox,ioy+1,ioz)
00815                                                                                 + (*image)(iox+1,ioy+1,ioz);
00816                                                                 float a61 = -(*image)(iox,ioy,ioz+1)
00817                                                                                 + (*image)(iox+1,ioy,ioz+1);
00818                                                                 float a6 = -a2 + a61;
00819                                                                 float a7 = -a3 - (*image)(iox,ioy,ioz+1)
00820                                                                                 + (*image)(iox,ioy+1,ioz+1);
00821                                                                 float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
00822                                                                                 + (*image)(iox+1,ioy+1,ioz+1);
00823                                                                 (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
00824                                                                                 + (a7 + a8*dx)*dy)
00825                                                                                 + a3*dy + dx*(a2 + a5*dy);
00826                                                         }
00827                                                 }
00828                                         }
00829                                         vb += rotation->get_matrix3_row(0);
00830                                 }
00831                         }
00832 
00833                 } else {
00834                         // No need to check x and y boundaries
00835                         for (int i = 0 ; i <= nn; i++) {
00836                                 int k = ipcube[i].loc[1] + origin[1];
00837                                 Vec3f vb = ipcube[i].loc*(*rotation) + origin;
00838                                 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
00839                                         int iox = int(vb[0]);
00840                                         int ioy = int(vb[1]);
00841                                         int ioz = int(vb[2]);
00842                                         float dx = vb[0] - iox;
00843                                         float dy = vb[1] - ioy;
00844                                         float dz = vb[2] - ioz;
00845                                         float a1 = (*image)(iox,ioy,ioz);
00846                                         float a2 = (*image)(iox+1,ioy,ioz) - a1;
00847                                         float a3 = (*image)(iox,ioy+1,ioz) - a1;
00848                                         float a4 = (*image)(iox,ioy,ioz+1) - a1;
00849                                         float a5 = -a2 -(*image)(iox,ioy+1,ioz)
00850                                                         + (*image)(iox+1,ioy+1,ioz);
00851                                         float a61 = -(*image)(iox,ioy,ioz+1)
00852                                                         + (*image)(iox+1,ioy,ioz+1);
00853                                         float a6 = -a2 + a61;
00854                                         float a7 = -a3 - (*image)(iox,ioy,ioz+1)
00855                                                         + (*image)(iox,ioy+1,ioz+1);
00856                                         float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
00857                                                         + (*image)(iox+1,ioy+1,ioz+1);
00858                                         (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
00859                                                         + (a7 + a8*dx)*dy)
00860                                                         + a3*dy + dx*(a2 + a5*dy);
00861                                         vb += rotation->get_matrix3_row(0);
00862                                 }
00863                         }
00864                 }
00865         }
00866         ret->update();
00867         EMDeleteArray(ipcube);
00868         if(rotation) {delete rotation; rotation=0;}
00869         return ret;
00870 }


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 Mon Mar 7 18:03:34 2011 for EMAN2 by  doxygen 1.4.7