EMAN2
Classes | Public Member Functions | Static Public Member Functions | Static Public Attributes | Private Member Functions
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.

Classes

struct  IPCube

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

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

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

{

        float *images;

        if (!imagestack) {
                return 0;
        }
        int ri;
        int nx      = imagestack->get_xsize();
        int ny      = imagestack->get_ysize();
//     int nslices = imagestack->get_zsize();
        int dim = Util::get_min(nx,ny);
        images  = imagestack->get_data();

        Vec3i origin(0,0,0);
    // If a sensible origin isn't passed in, choose the middle of
    // the cube.
        if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
        else {origin[0] = nx/2;}
        if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
        else {origin[1] = ny/2;}
        if (params.has_key("origin_z")) {origin[1] = params["origin_z"];}
        else {origin[2] = dim/2;}

        if (params.has_key("radius")) {ri = params["radius"];}
        else {ri = dim/2 - 1;}

    // Determine the number of rows (x-lines) within the radius
        int nn = -1;
        prepcubes(nx, ny, dim, ri, origin, nn);
    // nn is now the number of rows-1 within the radius
    // so we can create and fill the ipcubes
        IPCube* ipcube = new IPCube[nn+1];
        prepcubes(nx, ny, dim, ri, origin, nn, ipcube);

        int nangles = 0;
        vector<float> anglelist;
        // Do we have a list of angles?
        if (params.has_key("anglelist")) {
                anglelist = params["anglelist"];
                nangles = anglelist.size() / 3;
        } else {
                Transform* t3d = params["transform"];
                if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
                // This part was modified by David Woolford -
                // Before this the code worked only for SPIDER and EMAN angles,
                // but the framework of the Transform3D allows for a generic implementation
                // as specified here.
                Dict p = t3d->get_rotation("spider");
                if(t3d) {delete t3d; t3d=0;}

                string angletype = "SPIDER";
                float phi = p["phi"];
                float theta = p["theta"];
                float psi = p["psi"];
                anglelist.push_back(phi);
                anglelist.push_back(theta);
                anglelist.push_back(psi);
                nangles = 1;
        }

        // End David Woolford modifications

    // initialize return object
        EMData* ret = new EMData();
        ret->set_size(nx, ny, dim);
        ret->to_zero();

    // loop over sets of angles
        for (int ia = 0; ia < nangles; ia++) {
                int indx = 3*ia;
                Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]);
                Transform rotation(d);
                float dm1 = rotation.at(0,0);
                float dm4 = rotation.at(1,0);

                if (2*(ri+1)+1 > dim) {
          // Must check x and y boundaries
                        LOGERR("backproject3d, pawel, 2*(ri+1)+1 > dim\n");
                        return 0;
                } else {
          // No need to check x and y boundaries
                        for (int i = 0 ; i <= nn; i++) {
                                int iox = (int)ipcube[i].loc[0]+origin[0];
                                int ioy = (int)ipcube[i].loc[1]+origin[1];
                                int ioz = (int)ipcube[i].loc[2]+origin[2];

                                Vec3f vb = rotation*ipcube[i].loc + origin;
                                for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
                                        float xbb = (j-ipcube[i].start)*dm1 + vb[0];
                                        int   iqx = (int)floor(xbb);

                                        float ybb = (j-ipcube[i].start)*dm4 + vb[1];
                                        int   iqy = (int)floor(ybb);

                                        float dipx = xbb - iqx;
                                        float dipy = ybb - iqy;

                                        (*ret)(iox,ioy,ioz) += images(iqx,iqy,ia)
                                                        + dipy*(images(iqx,iqy+1,ia)-images(iqx,iqy,ia))
                                                        + dipx*(images(iqx+1,iqy,ia)-images(iqx,iqy,ia)
                                                                        + dipy*(images(iqx+1,iqy+1,ia)-images(iqx+1,iqy,ia)
                                                                                        - images(iqx,iqy+1,ia)+images(iqx,iqy,ia)));
                                        iox++;
                                } // end for j
                        } // end for i
                } // end if
        } // end for ia

        ret->update();
        EMDeleteArray(ipcube);
        return ret;
}
string EMAN::PawelProjector::get_desc ( ) const [inline, virtual]

Implements EMAN::Projector.

Definition at line 259 of file projector.h.

                {
                        return "Pawel Penczek's optimized real-space projection generation. Minor interpolation artifacts.";
                }
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.

                {
                        return NAME;
                }
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.

                {
                        TypeDict d;
                        d.put("transform", EMObject::TRANSFORM);
                        d.put("origin_x",  EMObject::INT);
                        d.put("origin_y",  EMObject::INT);
                        d.put("origin_z",  EMObject::INT);
                        d.put("radius",    EMObject::INT);
                        d.put("anglelist", EMObject::FLOATARRAY);
                        d.put("angletype", EMObject::STRING);
                        d.put("theta",     EMObject::FLOAT);
                        d.put("psi",       EMObject::FLOAT);
                        return d;
                }
static Projector* EMAN::PawelProjector::NEW ( ) [inline, static]

Definition at line 264 of file projector.h.

                {
                        return new PawelProjector();
                }
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(), nx, ny, EMAN::PawelProjector::IPCube::start, and t.

                                                                      {
        const float r = float(ri)*float(ri);
        const int ldpx = origin[0];
        const int ldpy = origin[1];
        const int ldpz = origin[2];
        //cout<<"  ldpx  "<<ldpx<<"  i2  "<<ldpy<<"  i1 "<<ldpz<<endl;
        float t;
        nn = -1;
        for (int i1 = 0; i1 < nz; i1++) {
                t = float(i1 - ldpz);
                const float xx = t*t;
                for (int i2 = 0; i2 < ny; i2++) {
                        t = float(i2 - ldpy);
                        const float yy = t*t + xx;
                        bool first = true;
                        for (int i3 = 0; i3 < nx; i3++) {
                                t = float(i3 - ldpx);
                                const float rc = t*t + yy;
                                if (first) {
                                        // first pixel on this line
                                        if (rc > r) continue;
                                        first = false;
                                        nn++;
                                        if (ipcube != NULL) {
                                                ipcube[nn].start = i3;
                                                ipcube[nn].end   = i3;
                                                ipcube[nn].loc[0] = i3 - ldpx;
                                                ipcube[nn].loc[1] = i2 - ldpy;
                                                ipcube[nn].loc[2] = i1 - ldpz;
                                                //cout<<"  start  "<<i3<<"  i2  "<<i2<<"  i1 "<<i1<<endl;
                                        }
                                } else {
                                        // second or later pixel on this line
                                        if (ipcube != NULL) {
                                                if (rc <= r) ipcube[nn].end = i3;
                                        }
                                }
                        }
                }
        }
}
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 576 of file projector.cpp.

References EMDeleteArray(), EMAN::PawelProjector::IPCube::end, EMAN::Transform::get_matrix3_row(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::PawelProjector::IPCube::loc, LOGERR, nn(), nx, ny, EMAN::EMData::set_size(), EMAN::EMData::to_zero(), and EMAN::EMData::update().

{
        if (!image)  return 0;
        int ri;
        int nx = image->get_xsize();
        int ny = image->get_ysize();
        int nz = image->get_zsize();
        int dim = Util::get_max(nx,ny,nz);
        int dimc = dim/2;
        if (nz == 1) {
                LOGERR("The PawelProjector needs a volume!");
                return 0;
        }

        Vec3i origin(0,0,0);

        if (params.has_key("radius")) {ri = params["radius"];}
        else {ri = dim/2 - 1;}

        Transform* rotation = params["transform"];
        //int nangles = 0;
        //vector<float> anglelist;
        // Do we have a list of angles?
        /*
        if (params.has_key("anglelist")) {
                anglelist = params["anglelist"];
                nangles = anglelist.size() / 3;
        } else {*/

                //if ( rotation == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
                /*
                Dict p = t3d->get_rotation("spider");

                string angletype = "SPIDER";
                float phi = p["phi"];
                float theta = p["theta"];
                float psi = p["psi"];
                anglelist.push_back(phi);
                anglelist.push_back(theta);
                anglelist.push_back(psi);
                */
                int nangles = 1;
        //}

        //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;
        // loop over sets of angles
        //for (int ia = 0; ia < nangles; ia++) {
        EMData* ret = new EMData();
        int ia = 0;
                //int indx = 3*ia;
                //Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]);
                //Transform rotation(d);
                if (2*(ri+1)+1 > dim) {
                        // initialize return object
                        ret->set_size(dim, dim, nangles);
                        ret->to_zero();
                        origin[0] = dimc;
                        origin[1] = dimc;
                        origin[2] = dimc;
                        Vec3i loc(-dimc,0,0);
                        Vec3i vorg(nx/2,ny/2,nz/2);
                        // This code is for arbitrary dimensions, so must check x and y boundaries
                        for (int l = 0 ; l < dim; l++) {  // Z
                                loc[2] = l - dimc;
                                for (int k = 0 ; k < dim; k++) {  // Y
                                        loc[1] = k - dimc;
                                        Vec3f vb = loc*(*rotation) + vorg;
                                        for (int j = 0; j < dim; j++) {  //X
                                                // check for pixels out-of-bounds
                                                //cout<<j<<"  j  "<<k<<"  k  "<<"  iox  "<<vb[0]<<"  ioy  "<<vb[1]<<"  ioz "<<vb[2]<<endl;
                                                int iox = int(vb[0]);
                                                if ((iox >= 0) && (iox < nx-1)) {
                                                        int ioy = int(vb[1]);
                                                        if ((ioy >= 0) && (ioy < ny-1)) {
                                                                int ioz = int(vb[2]);
                                                                if ((ioz >= 0) && (ioz < nz-1)) {
                                                                        // real work for pixels in bounds
                                                                //cout<<j<<"  j  "<<k<<"  k  "<<"  iox  "<<iox<<"  ioy  "<<ioy<<"  ioz "<<ioz<<endl;
                                                                //cout<<"  TAKE"<<endl;
                                                                        float dx = vb[0] - iox;
                                                                        float dy = vb[1] - ioy;
                                                                        float dz = vb[2] - ioz;
                                                                        float a1 = (*image)(iox,ioy,ioz);
                                                                        float a2 = (*image)(iox+1,ioy,ioz) - a1;
                                                                        float a3 = (*image)(iox,ioy+1,ioz) - a1;
                                                                        float a4 = (*image)(iox,ioy,ioz+1) - a1;
                                                                        float a5 = -a2 -(*image)(iox,ioy+1,ioz)
                                                                                                + (*image)(iox+1,ioy+1,ioz);
                                                                        float a61 = -(*image)(iox,ioy,ioz+1)
                                                                                                + (*image)(iox+1,ioy,ioz+1);
                                                                        float a6 = -a2 + a61;
                                                                        float a7 = -a3 - (*image)(iox,ioy,ioz+1)
                                                                                                + (*image)(iox,ioy+1,ioz+1);
                                                                        float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
                                                                                                + (*image)(iox+1,ioy+1,ioz+1);
                                                                        (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
                                                                                                + (a7 + a8*dx)*dy)
                                                                                                + a3*dy + dx*(a2 + a5*dy);
                                                                } //else {cout<<"  iox  "<<iox<<"  ioy  "<<int(vb[1])<<"  ioz "<<int(vb[2])<<"  VIOLAATED Z "<<endl; }
                                                        }//else {cout<<"  iox  "<<iox<<"  ioy  "<<int(vb[1])<<"  ioz "<<int(vb[2])<<"  VIOLATED Y "<<endl; }
                                                }//else {cout<<"  iox  "<<iox<<"  ioy  "<<int(vb[1])<<"  ioz "<<int(vb[2])<<"  VIOLATED X "<<endl; }
                                                vb += rotation->get_matrix3_row(0);
                                        }
                                }
                        }

                } else {
                        // If a sensible origin isn't passed in, choose the middle of
                        // the cube.
                        if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
                        else {origin[0] = nx/2;}
                        if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
                        else {origin[1] = ny/2;}
                        if (params.has_key("origin_z")) {origin[2] = params["origin_z"];}
                        else {origin[2] = nz/2;}
                        // Determine the number of rows (x-lines) within the radius
                        int nn = -1;
                        prepcubes(nx, ny, nz, ri, origin, nn);
                        // nn is now the number of rows-1 within the radius
                        // so we can create and fill the ipcubes
                        IPCube* ipcube = new IPCube[nn+1];
                        prepcubes(nx, ny, nz, ri, origin, nn, ipcube);
                        // initialize return object
                        ret->set_size(nx, ny, nangles);
                        ret->to_zero();
                        // No need to check x and y boundaries
                        for (int i = 0 ; i <= nn; i++) {
                                int k = ipcube[i].loc[1] + origin[1];
                                Vec3f vb = ipcube[i].loc*(*rotation) + origin;
                                for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
                                        int iox = int(vb[0]);
                                        int ioy = int(vb[1]);
                                        int ioz = int(vb[2]);
                                        float dx = vb[0] - iox;
                                        float dy = vb[1] - ioy;
                                        float dz = vb[2] - ioz;
                                        float a1 = (*image)(iox,ioy,ioz);
                                        float a2 = (*image)(iox+1,ioy,ioz) - a1;
                                        float a3 = (*image)(iox,ioy+1,ioz) - a1;
                                        float a4 = (*image)(iox,ioy,ioz+1) - a1;
                                        float a5 = -a2 -(*image)(iox,ioy+1,ioz)
                                                                + (*image)(iox+1,ioy+1,ioz);
                                        float a61 = -(*image)(iox,ioy,ioz+1)
                                                                + (*image)(iox+1,ioy,ioz+1);
                                        float a6 = -a2 + a61;
                                        float a7 = -a3 - (*image)(iox,ioy,ioz+1)
                                                                + (*image)(iox,ioy+1,ioz+1);
                                        float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
                                                                + (*image)(iox+1,ioy+1,ioz+1);
                                        (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
                                                                + (a7 + a8*dx)*dy)
                                                                + a3*dy + dx*(a2 + a5*dy);
                                        vb += rotation->get_matrix3_row(0);
                                }
                        }
                        EMDeleteArray(ipcube);
                }
        //}
        ret->update();
        if(rotation) {delete rotation; rotation=0;}
        return ret;
}

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: