EMAN2
Classes | Public Member Functions | Static Public Member Functions | Static Public Attributes | Private Member Functions | List of all 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]

Classes

struct  IPCube
 

Public Member Functions

EMData * project3d (EMData *image) const
 Project an 3D image into a 2D image. More...
 
EMData * backproject3d (EMData *image) const
 Back-project a 2D image into a 3D image. More...
 
string get_name () const
 Get the projector's name. More...
 
string get_desc () const
 
TypeDict get_param_types () const
 Get processor parameter information in a dictionary. More...
 
- Public Member Functions inherited from EMAN::Projector
virtual ~Projector ()
 
virtual Dict get_params () const
 Get the projector parameters in a key/value dictionary. More...
 
void set_params (const Dict &new_params)
 Set the projector parameters using a key/value dictionary. More...
 

Static Public Member Functions

static Projector * NEW ()
 

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
 

Additional Inherited Members

- Protected Attributes inherited from EMAN::Projector
Dict params
 

Detailed Description

Pawel Penczek's optimized projection routine.

Definition at line 245 of file projector.h.

Member Function Documentation

◆ backproject3d()

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

2014{
2015
2016 float *images;
2017
2018 if (!imagestack) {
2019 return 0;
2020 }
2021 int ri;
2022 int nx = imagestack->get_xsize();
2023 int ny = imagestack->get_ysize();
2024// int nslices = imagestack->get_zsize();
2025 int dim = Util::get_min(nx,ny);
2026 images = imagestack->get_data();
2027
2028 Vec3i origin(0,0,0);
2029 // If a sensible origin isn't passed in, choose the middle of
2030 // the cube.
2031 if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
2032 else {origin[0] = nx/2;}
2033 if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
2034 else {origin[1] = ny/2;}
2035 if (params.has_key("origin_z")) {origin[1] = params["origin_z"];}
2036 else {origin[2] = dim/2;}
2037
2038 if (params.has_key("radius")) {ri = params["radius"];}
2039 else {ri = dim/2 - 1;}
2040
2041 // Determine the number of rows (x-lines) within the radius
2042 int nn = -1;
2043 prepcubes(nx, ny, dim, ri, origin, nn);
2044 // nn is now the number of rows-1 within the radius
2045 // so we can create and fill the ipcubes
2046 IPCube* ipcube = new IPCube[nn+1];
2047 prepcubes(nx, ny, dim, ri, origin, nn, ipcube);
2048
2049 int nangles = 0;
2050 vector<float> anglelist;
2051 // Do we have a list of angles?
2052 if (params.has_key("anglelist")) {
2053 anglelist = params["anglelist"];
2054 nangles = anglelist.size() / 3;
2055 } else {
2056 Transform* t3d = params["transform"];
2057 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
2058 // This part was modified by David Woolford -
2059 // Before this the code worked only for SPIDER and EMAN angles,
2060 // but the framework of the Transform3D allows for a generic implementation
2061 // as specified here.
2062 Dict p = t3d->get_rotation("spider");
2063 if(t3d) {delete t3d; t3d=0;}
2064
2065 string angletype = "SPIDER";
2066 float phi = p["phi"];
2067 float theta = p["theta"];
2068 float psi = p["psi"];
2069 anglelist.push_back(phi);
2070 anglelist.push_back(theta);
2071 anglelist.push_back(psi);
2072 nangles = 1;
2073 }
2074
2075 // End David Woolford modifications
2076
2077 // initialize return object
2078 EMData* ret = new EMData();
2079 ret->set_size(nx, ny, dim);
2080 ret->to_zero();
2081
2082 // loop over sets of angles
2083 for (int ia = 0; ia < nangles; ia++) {
2084 int indx = 3*ia;
2085 Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]);
2086 Transform rotation(d);
2087 float dm1 = rotation.at(0,0);
2088 float dm4 = rotation.at(1,0);
2089
2090 if (2*(ri+1)+1 > dim) {
2091 // Must check x and y boundaries
2092 LOGERR("backproject3d, pawel, 2*(ri+1)+1 > dim\n");
2093 return 0;
2094 } else {
2095 // No need to check x and y boundaries
2096 for (int i = 0 ; i <= nn; i++) {
2097 int iox = (int)ipcube[i].loc[0]+origin[0];
2098 int ioy = (int)ipcube[i].loc[1]+origin[1];
2099 int ioz = (int)ipcube[i].loc[2]+origin[2];
2100
2101 Vec3f vb = rotation*ipcube[i].loc + origin;
2102 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
2103 float xbb = (j-ipcube[i].start)*dm1 + vb[0];
2104 int iqx = (int)floor(xbb);
2105
2106 float ybb = (j-ipcube[i].start)*dm4 + vb[1];
2107 int iqy = (int)floor(ybb);
2108
2109 float dipx = xbb - iqx;
2110 float dipy = ybb - iqy;
2111
2112 (*ret)(iox,ioy,ioz) += images(iqx,iqy,ia)
2113 + dipy*(images(iqx,iqy+1,ia)-images(iqx,iqy,ia))
2114 + dipx*(images(iqx+1,iqy,ia)-images(iqx,iqy,ia)
2115 + dipy*(images(iqx+1,iqy+1,ia)-images(iqx+1,iqy,ia)
2116 - images(iqx,iqy+1,ia)+images(iqx,iqy,ia)));
2117 iox++;
2118 } // end for j
2119 } // end for i
2120 } // end if
2121 } // end for ia
2122
2123 ret->update();
2124 EMDeleteArray(ipcube);
2125 return ret;
2126}
EMAN::Dict
Dict is a dictionary to store <string, EMObject> pair.
Definition: emobject.h:385
EMAN::Dict::has_key
bool has_key(const string &key) const
Ask the Dictionary if it as a particular key.
Definition: emobject.h:511
EMAN::EMData
EMData stores an image's data and defines core image processing routines.
Definition: emdata.h:82
EMAN::PawelProjector::prepcubes
void prepcubes(int nx, int ny, int nz, int ri, Vec3i origin, int &nn, IPCube *ipcube=NULL) const
Definition: projector.cpp:559
EMAN::Transform
A Transform object is a somewhat specialized object designed specifically for EMAN2/Sparx storage of ...
Definition: transform.h:75
EMAN::Transform::get_rotation
Dict get_rotation(const string &euler_type="eman") const
Get a rotation in any Euler format.
Definition: transform.cpp:829
EMAN::Util::get_min
static int get_min(int f1, int f2)
Get the minimum of 2 numbers.
Definition: util.h:922
EMDeleteArray
void EMDeleteArray(T &x)
Definition: emutil.h:62
NullPointerException
#define NullPointerException(desc)
Definition: exception.h:241
LOGERR
#define LOGERR
Definition: log.h:51
anglelist
#define anglelist(i, j)
Definition: projector.cpp:1607
images
#define images(i, j, k)
Definition: projector.cpp:1897

References anglelist, EMAN::Transform::at(), EMDeleteArray(), EMAN::PawelProjector::IPCube::end, EMAN::Util::get_min(), EMAN::Transform::get_rotation(), EMAN::Dict::has_key(), images, EMAN::PawelProjector::IPCube::loc, LOGERR, NullPointerException, EMAN::Projector::params, prepcubes(), and EMAN::PawelProjector::IPCube::start.

◆ get_desc()

string EMAN::PawelProjector::get_desc ( ) const
inlinevirtual

Implements EMAN::Projector.

Definition at line 256 of file projector.h.

257 {
258 return "Pawel Penczek's optimized real-space projection generation. Minor interpolation artifacts.";
259 }

◆ get_name()

string EMAN::PawelProjector::get_name ( ) const
inlinevirtual

Get the projector's name.

Each projector is indentified by unique name.

Returns
The projector's name.

Implements EMAN::Projector.

Definition at line 251 of file projector.h.

252 {
253 return NAME;
254 }
EMAN::PawelProjector::NAME
static const string NAME
Definition: projector.h:281

References NAME.

◆ get_param_types()

TypeDict EMAN::PawelProjector::get_param_types ( ) const
inlinevirtual

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 266 of file projector.h.

267 {
268 TypeDict d;
269 d.put("transform", EMObject::TRANSFORM);
270 d.put("origin_x", EMObject::INT);
271 d.put("origin_y", EMObject::INT);
272 d.put("origin_z", EMObject::INT);
273 d.put("radius", EMObject::INT);
274 d.put("anglelist", EMObject::FLOATARRAY);
275 d.put("angletype", EMObject::STRING);
276 d.put("theta", EMObject::FLOAT);
277 d.put("psi", EMObject::FLOAT);
278 return d;
279 }

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

◆ NEW()

static Projector * EMAN::PawelProjector::NEW ( )
inlinestatic

Definition at line 261 of file projector.h.

262 {
263 return new PawelProjector();
264 }

◆ prepcubes()

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

Definition at line 559 of file projector.cpp.

560 {
561 const float r = float(ri)*float(ri);
562 const int ldpx = origin[0];
563 const int ldpy = origin[1];
564 const int ldpz = origin[2];
565 //cout<<" ldpx "<<ldpx<<" i2 "<<ldpy<<" i1 "<<ldpz<<endl;
566 float t;
567 nn = -1;
568 for (int i1 = 0; i1 < nz; i1++) {
569 t = float(i1 - ldpz);
570 const float xx = t*t;
571 for (int i2 = 0; i2 < ny; i2++) {
572 t = float(i2 - ldpy);
573 const float yy = t*t + xx;
574 bool first = true;
575 for (int i3 = 0; i3 < nx; i3++) {
576 t = float(i3 - ldpx);
577 const float rc = t*t + yy;
578 if (first) {
579 // first pixel on this line
580 if (rc > r) continue;
581 first = false;
582 nn++;
583 if (ipcube != NULL) {
584 ipcube[nn].start = i3;
585 ipcube[nn].end = i3;
586 ipcube[nn].loc[0] = i3 - ldpx;
587 ipcube[nn].loc[1] = i2 - ldpy;
588 ipcube[nn].loc[2] = i1 - ldpz;
589 //cout<<" start "<<i3<<" i2 "<<i2<<" i1 "<<i1<<endl;
590 }
591 } else {
592 // second or later pixel on this line
593 if (ipcube != NULL) {
594 if (rc <= r) ipcube[nn].end = i3;
595 }
596 }
597 }
598 }
599 }
600}

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

Referenced by backproject3d(), and project3d().

◆ project3d()

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

603{
604 if (!image) return 0;
605 int ri;
606 int nx = image->get_xsize();
607 int ny = image->get_ysize();
608 int nz = image->get_zsize();
609 int dim = Util::get_max(nx,ny,nz);
610 int dimc = dim/2;
611 if (nz == 1) {
612 LOGERR("The PawelProjector needs a volume!");
613 return 0;
614 }
615
616 Vec3i origin(0,0,0);
617
618 if (params.has_key("radius")) {ri = params["radius"];}
619 else {ri = dim/2 - 1;}
620
621 Transform* rotation = params["transform"];
622 //int nangles = 0;
623 //vector<float> anglelist;
624 // Do we have a list of angles?
625 /*
626 if (params.has_key("anglelist")) {
627 anglelist = params["anglelist"];
628 nangles = anglelist.size() / 3;
629 } else {*/
630
631 //if ( rotation == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
632 /*
633 Dict p = t3d->get_rotation("spider");
634
635 string angletype = "SPIDER";
636 float phi = p["phi"];
637 float theta = p["theta"];
638 float psi = p["psi"];
639 anglelist.push_back(phi);
640 anglelist.push_back(theta);
641 anglelist.push_back(psi);
642 */
643 int nangles = 1;
644 //}
645
646 //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;
647 // loop over sets of angles
648 //for (int ia = 0; ia < nangles; ia++) {
649 EMData* ret = new EMData();
650 int ia = 0;
651 //int indx = 3*ia;
652 //Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]);
653 //Transform rotation(d);
654 if (2*(ri+1)+1 > dim) {
655 // initialize return object
656 ret->set_size(dim, dim, nangles);
657 ret->to_zero();
658 origin[0] = dimc;
659 origin[1] = dimc;
660 origin[2] = dimc;
661 Vec3i loc(-dimc,0,0);
662 Vec3i vorg(nx/2,ny/2,nz/2);
663 // This code is for arbitrary dimensions, so must check x and y boundaries
664 for (int l = 0 ; l < dim; l++) { // Z
665 loc[2] = l - dimc;
666 for (int k = 0 ; k < dim; k++) { // Y
667 loc[1] = k - dimc;
668 Vec3f vb = loc*(*rotation) + vorg;
669 for (int j = 0; j < dim; j++) { //X
670 // check for pixels out-of-bounds
671 //cout<<j<<" j "<<k<<" k "<<" iox "<<vb[0]<<" ioy "<<vb[1]<<" ioz "<<vb[2]<<endl;
672 int iox = int(vb[0]);
673 if ((iox >= 0) && (iox < nx-1)) {
674 int ioy = int(vb[1]);
675 if ((ioy >= 0) && (ioy < ny-1)) {
676 int ioz = int(vb[2]);
677 if ((ioz >= 0) && (ioz < nz-1)) {
678 // real work for pixels in bounds
679 //cout<<j<<" j "<<k<<" k "<<" iox "<<iox<<" ioy "<<ioy<<" ioz "<<ioz<<endl;
680 //cout<<" TAKE"<<endl;
681 float dx = vb[0] - iox;
682 float dy = vb[1] - ioy;
683 float dz = vb[2] - ioz;
684 float a1 = (*image)(iox,ioy,ioz);
685 float a2 = (*image)(iox+1,ioy,ioz) - a1;
686 float a3 = (*image)(iox,ioy+1,ioz) - a1;
687 float a4 = (*image)(iox,ioy,ioz+1) - a1;
688 float a5 = -a2 -(*image)(iox,ioy+1,ioz)
689 + (*image)(iox+1,ioy+1,ioz);
690 float a61 = -(*image)(iox,ioy,ioz+1)
691 + (*image)(iox+1,ioy,ioz+1);
692 float a6 = -a2 + a61;
693 float a7 = -a3 - (*image)(iox,ioy,ioz+1)
694 + (*image)(iox,ioy+1,ioz+1);
695 float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
696 + (*image)(iox+1,ioy+1,ioz+1);
697 (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
698 + (a7 + a8*dx)*dy)
699 + a3*dy + dx*(a2 + a5*dy);
700 } //else {cout<<" iox "<<iox<<" ioy "<<int(vb[1])<<" ioz "<<int(vb[2])<<" VIOLAATED Z "<<endl; }
701 }//else {cout<<" iox "<<iox<<" ioy "<<int(vb[1])<<" ioz "<<int(vb[2])<<" VIOLATED Y "<<endl; }
702 }//else {cout<<" iox "<<iox<<" ioy "<<int(vb[1])<<" ioz "<<int(vb[2])<<" VIOLATED X "<<endl; }
703 vb += rotation->get_matrix3_row(0);
704 }
705 }
706 }
707
708 } else {
709 // If a sensible origin isn't passed in, choose the middle of
710 // the cube.
711 if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
712 else {origin[0] = nx/2;}
713 if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
714 else {origin[1] = ny/2;}
715 if (params.has_key("origin_z")) {origin[2] = params["origin_z"];}
716 else {origin[2] = nz/2;}
717 // Determine the number of rows (x-lines) within the radius
718 int nn = -1;
719 prepcubes(nx, ny, nz, ri, origin, nn);
720 // nn is now the number of rows-1 within the radius
721 // so we can create and fill the ipcubes
722 IPCube* ipcube = new IPCube[nn+1];
723 prepcubes(nx, ny, nz, ri, origin, nn, ipcube);
724 // initialize return object
725 ret->set_size(nx, ny, nangles);
726 ret->to_zero();
727 // No need to check x and y boundaries
728 for (int i = 0 ; i <= nn; i++) {
729 int k = ipcube[i].loc[1] + origin[1];
730 Vec3f vb = ipcube[i].loc*(*rotation) + origin;
731 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
732 int iox = int(vb[0]);
733 int ioy = int(vb[1]);
734 int ioz = int(vb[2]);
735 float dx = vb[0] - iox;
736 float dy = vb[1] - ioy;
737 float dz = vb[2] - ioz;
738 float a1 = (*image)(iox,ioy,ioz);
739 float a2 = (*image)(iox+1,ioy,ioz) - a1;
740 float a3 = (*image)(iox,ioy+1,ioz) - a1;
741 float a4 = (*image)(iox,ioy,ioz+1) - a1;
742 float a5 = -a2 -(*image)(iox,ioy+1,ioz)
743 + (*image)(iox+1,ioy+1,ioz);
744 float a61 = -(*image)(iox,ioy,ioz+1)
745 + (*image)(iox+1,ioy,ioz+1);
746 float a6 = -a2 + a61;
747 float a7 = -a3 - (*image)(iox,ioy,ioz+1)
748 + (*image)(iox,ioy+1,ioz+1);
749 float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
750 + (*image)(iox+1,ioy+1,ioz+1);
751 (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
752 + (a7 + a8*dx)*dy)
753 + a3*dy + dx*(a2 + a5*dy);
754 vb += rotation->get_matrix3_row(0);
755 }
756 }
757 EMDeleteArray(ipcube);
758 }
759 //}
760 ret->update();
761 if(rotation) {delete rotation; rotation=0;}
762 return ret;
763}
EMAN::Transform::get_matrix3_row
Vec3f get_matrix3_row(int i) const
Get a matrix row as a Vec3f required for back compatibility with Tranform3D - see PawelProjector.
Definition: transform.h:465
EMAN::Util::get_max
static float get_max(float f1, float f2)
Get the maximum of 2 numbers.
Definition: util.h:998

References EMDeleteArray(), EMAN::PawelProjector::IPCube::end, EMAN::Transform::get_matrix3_row(), EMAN::Util::get_max(), EMAN::Dict::has_key(), EMAN::PawelProjector::IPCube::loc, LOGERR, EMAN::Projector::params, and prepcubes().

Member Data Documentation

◆ NAME

const string PawelProjector::NAME = "pawel"
static

Definition at line 281 of file projector.h.

Referenced by get_name().

The documentation for this class was generated from the following files:
Generated by doxygen 1.9.2