Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

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

ProjectorNEW ()

Static Public Attributes

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

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 2127 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, phi, prepcubes(), EMAN::EMData::set_size(), EMAN::Dict::size(), EMAN::PawelProjector::IPCube::start, theta, EMAN::EMData::to_zero(), EMAN::EMData::update(), EMAN::Vec3f, and EMAN::Vec3i.

02128 {
02129 
02130         float *images;
02131 
02132         if (!imagestack) {
02133                 return 0;
02134         }
02135         int ri;
02136         int nx      = imagestack->get_xsize();
02137         int ny      = imagestack->get_ysize();
02138 //     int nslices = imagestack->get_zsize();
02139         int dim = Util::get_min(nx,ny);
02140         images  = imagestack->get_data();
02141 
02142         Vec3i origin(0,0,0);
02143     // If a sensible origin isn't passed in, choose the middle of
02144     // the cube.
02145         if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
02146         else {origin[0] = nx/2;}
02147         if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
02148         else {origin[1] = ny/2;}
02149         if (params.has_key("origin_z")) {origin[1] = params["origin_z"];}
02150         else {origin[2] = dim/2;}
02151 
02152         if (params.has_key("radius")) {ri = params["radius"];}
02153         else {ri = dim/2 - 1;}
02154 
02155     // Determine the number of rows (x-lines) within the radius
02156         int nn = -1;
02157         prepcubes(nx, ny, dim, ri, origin, nn);
02158     // nn is now the number of rows-1 within the radius
02159     // so we can create and fill the ipcubes
02160         IPCube* ipcube = new IPCube[nn+1];
02161         prepcubes(nx, ny, dim, ri, origin, nn, ipcube);
02162 
02163         int nangles = 0;
02164         vector<float> anglelist;
02165         // Do we have a list of angles?
02166         if (params.has_key("anglelist")) {
02167                 anglelist = params["anglelist"];
02168                 nangles = anglelist.size() / 3;
02169         } else {
02170                 Transform* t3d = params["transform"];
02171                 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
02172                 // This part was modified by David Woolford -
02173                 // Before this the code worked only for SPIDER and EMAN angles,
02174                 // but the framework of the Transform3D allows for a generic implementation
02175                 // as specified here.
02176                 Dict p = t3d->get_rotation("spider");
02177                 if(t3d) {delete t3d; t3d=0;}
02178 
02179                 string angletype = "SPIDER";
02180                 float phi = p["phi"];
02181                 float theta = p["theta"];
02182                 float psi = p["psi"];
02183                 anglelist.push_back(phi);
02184                 anglelist.push_back(theta);
02185                 anglelist.push_back(psi);
02186                 nangles = 1;
02187         }
02188 
02189         // End David Woolford modifications
02190 
02191     // initialize return object
02192         EMData* ret = new EMData();
02193         ret->set_size(nx, ny, dim);
02194         ret->to_zero();
02195 
02196     // loop over sets of angles
02197         for (int ia = 0; ia < nangles; ia++) {
02198                 int indx = 3*ia;
02199                 Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]);
02200                 Transform rotation(d);
02201                 float dm1 = rotation.at(0,0);
02202                 float dm4 = rotation.at(1,0);
02203 
02204                 if (2*(ri+1)+1 > dim) {
02205           // Must check x and y boundaries
02206                         LOGERR("backproject3d, pawel, 2*(ri+1)+1 > dim\n");
02207                         return 0;
02208                 } else {
02209           // No need to check x and y boundaries
02210                         for (int i = 0 ; i <= nn; i++) {
02211                                 int iox = (int)ipcube[i].loc[0]+origin[0];
02212                                 int ioy = (int)ipcube[i].loc[1]+origin[1];
02213                                 int ioz = (int)ipcube[i].loc[2]+origin[2];
02214 
02215                                 Vec3f vb = rotation*ipcube[i].loc + origin;
02216                                 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
02217                                         float xbb = (j-ipcube[i].start)*dm1 + vb[0];
02218                                         int   iqx = (int)floor(xbb);
02219 
02220                                         float ybb = (j-ipcube[i].start)*dm4 + vb[1];
02221                                         int   iqy = (int)floor(ybb);
02222 
02223                                         float dipx = xbb - iqx;
02224                                         float dipy = ybb - iqy;
02225 
02226                                         (*ret)(iox,ioy,ioz) += images(iqx,iqy,ia)
02227                                                         + dipy*(images(iqx,iqy+1,ia)-images(iqx,iqy,ia))
02228                                                         + dipx*(images(iqx+1,iqy,ia)-images(iqx,iqy,ia)
02229                                                                         + dipy*(images(iqx+1,iqy+1,ia)-images(iqx+1,iqy,ia)
02230                                                                                         - images(iqx,iqy+1,ia)+images(iqx,iqy,ia)));
02231                                         iox++;
02232                                 } // end for j
02233                         } // end for i
02234                 } // end if
02235         } // end for ia
02236 
02237         ret->update();
02238         EMDeleteArray(ipcube);
02239         return ret;
02240 }

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.

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::TypeDict::put().

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                 }

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

References EMAN::PawelProjector::IPCube::end, EMAN::PawelProjector::IPCube::loc, nn(), EMAN::PawelProjector::IPCube::start, t, and EMAN::Vec3i.

Referenced by backproject3d(), and project3d().

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

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

References 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, prepcubes(), EMAN::EMData::set_size(), EMAN::PawelProjector::IPCube::start, EMAN::EMData::to_zero(), EMAN::EMData::update(), EMAN::Vec3f, and EMAN::Vec3i.

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


Member Data Documentation

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

Definition at line 56 of file projector.cpp.


The documentation for this class was generated from the following files:
Generated on Tue Jun 11 13:49:20 2013 for EMAN2 by  doxygen 1.3.9.1