EMAN2
Public Member Functions | Static Public Member Functions | Static Public Attributes | Private Member Functions
EMAN::ChaoProjector Class Reference

Fast real space projection using Bi-Linear interpolation. More...

#include <projector.h>

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

List of all members.

Public Member Functions

EMDataproject3d (EMData *vol) const
 Project an 3D image into a 2D image.
EMDatabackproject3d (EMData *imagestack) 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 = "chao"

Private Member Functions

int getnnz (Vec3i volsize, int ri, Vec3i origin, int *nray, int *nnz) const
int cb2sph (float *cube, Vec3i volsize, int ri, Vec3i origin, int nnz0, int *ptrs, int *cord, float *sphere) const
int sph2cb (float *sphere, Vec3i volsize, int nray, int ri, int nnz0, int *ptrs, int *cord, float *cube) const
int fwdpj3 (Vec3i volsize, int nray, int nnz, float *dm, Vec3i origin, int ri, int *ptrs, int *cord, float *x, float *y) const
int bckpj3 (Vec3i volsize, int nray, int nnz, float *dm, Vec3i origin, int ri, int *ptrs, int *cord, float *x, float *y) const
int ifix (float a) const
void setdm (vector< float > anglelist, string const angletype, float *dm) const

Detailed Description

Fast real space projection using Bi-Linear interpolation.

(C. Yang)

Definition at line 374 of file projector.h.


Member Function Documentation

EMData * ChaoProjector::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 1741 of file projector.cpp.

References anglelist, bckpj3(), cb2sph(), cord, cube, dm, EMDeleteArray(), EMAN::EMData::get_data(), EMAN::Transform::get_rotation(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), getnnz(), images, LOGERR, nnz, nrays, NullPointerException, phi, ptrs, EMAN::EMData::set_complex(), EMAN::EMData::set_ri(), EMAN::EMData::set_size(), sph2cb(), sphere, status, theta, EMAN::EMData::to_zero(), and EMAN::EMData::update().

{
        int nrays, nnz, status, j;
        float *dm;
        int   *ptrs, *cord;
        float *sphere, *images, *cube;

        int nximg   = imagestack->get_xsize();
        int nyimg   = imagestack->get_ysize();
        int nslices = imagestack->get_zsize();

        int dim = Util::get_min(nximg,nyimg);
        Vec3i volsize(nximg,nyimg,dim);

        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] = nximg/2+1;}
        if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
        else {origin[1] = nyimg/2+1;}
        if (params.has_key("origin_z")) {origin[1] = params["origin_z"];}
        else {origin[2] = dim/2+1;}

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

        // retrieve the voxel values
        images = imagestack->get_data();

        // count the number of voxels within a sphere centered at icent,
        // with radius ri
        status = getnnz(volsize, ri, origin, &nrays, &nnz);
        // need to check status...

        // convert from cube to sphere
        sphere = new float[nnz];
        ptrs   = new int[nrays+1];
        cord   = new int[3*nrays];
        if (sphere == NULL || ptrs == NULL || cord == NULL) {
                fprintf(stderr,"ChaoProjector::backproject3d, failed to allocate!\n");
                exit(1);
        }
        for (int i = 0; i<nnz; i++) sphere[i] = 0.0;
        for (int i = 0; i<nrays+1; i++) ptrs[i] = 0;
        for (int i = 0; i<3*nrays; i++) cord[i] = 0;

        int nangles = 0;
        vector<float> anglelist;
        string angletype = "SPIDER";
        // 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.
                //  This was broken by david.  we need here a loop over all projections and put all angles on stack  PAP 06/28/09
                Dict p = t3d->get_rotation("spider");
                if(t3d) {delete t3d; t3d=0;}

                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

        if (nslices != nangles) {
                LOGERR("the number of images does not match the number of angles");
                return 0;
        }

        dm = new float[nangles*9];
        setdm(anglelist, angletype, dm);

        // return volume
        EMData *ret = new EMData();
        ret->set_size(nximg, nyimg, dim);
        ret->set_complex(false);
        ret->set_ri(true);
        ret->to_zero();

        cube = ret->get_data();
        // cb2sph should be replaced by something that touches only ptrs and cord
        status = cb2sph(cube, volsize, ri, origin, nnz, ptrs, cord, sphere);
        // check status

        for (j = 1; j <= nangles; j++) {
                status = bckpj3(volsize, nrays, nnz, &dm(1,j), origin, ri,
                         ptrs   , cord , &images(1,1,j), sphere);
        // check status?
        }

        status = sph2cb(sphere, volsize, nrays, ri, nnz, ptrs, cord, cube);
        // check status?

        // deallocate all temporary work space
        EMDeleteArray(dm);
        EMDeleteArray(ptrs);
        EMDeleteArray(cord);
        EMDeleteArray(sphere);

        ret->update();
        return ret;
}
int ChaoProjector::bckpj3 ( Vec3i  volsize,
int  nray,
int  nnz,
float *  dm,
Vec3i  origin,
int  ri,
int *  ptrs,
int *  cord,
float *  x,
float *  y 
) const [private]

Definition at line 1494 of file projector.cpp.

References cord, dm, ifix(), nrays, nx, ptrs, status, x, and y.

{
    int       i, j, iqx,iqy, xc, yc, zc;
    float     xb, yb, dx, dy, dx1m, dy1m, dxdy;
    int       status = 0;

    int xcent = origin[0];
    int ycent = origin[1];
    int zcent = origin[2];

    int nx = volsize[0];

    if ( nx > 2*ri) {
        for (i = 1; i <= nrays; i++) {
            zc = cord(1,i) - zcent;
            yc = cord(2,i) - ycent;
            xc = cord(3,i) - xcent;

            xb = zc*dm(1)+yc*dm(2)+xc*dm(3) + xcent;
            yb = zc*dm(4)+yc*dm(5)+xc*dm(6) + ycent;

            for (j = ptrs(i); j <ptrs(i+1); j++) {
                iqx = ifix((float)(xb));
                iqy = ifix((float)(yb));

                dx = xb - (float)(iqx);
                dy = yb - (float)(iqy);
                dx1m = 1.0f - dx;
                dy1m = 1.0f - dy;
                dxdy = dx*dy;
/*
c               y(j) = y(j) + dx1m*dy1m*x(iqx  , iqy)
c     &                     + dx1m*dy  *x(iqx  , iqy+1)
c     &                     + dx  *dy1m*x(iqx+1, iqy)
c     &                     + dx  *dy  *x(iqx+1, iqy+1)
c
c              --- faster version of the above commented out
c                  code (derived by summing the following table
c                  of coefficients along  the colunms) ---
c
c                        1         dx        dy      dxdy
c                     ------   --------  --------  -------
c                      x(i,j)   -x(i,j)   -x(i,j)    x(i,j)
c                                        x(i,j+1) -x(i,j+1)
c                              x(i+1,j)           -x(i+1,j)
c                                                x(i+1,j+1)
c
*/
               y(j) += x(iqx,iqy)
                    +  dx*(-x(iqx,iqy)+x(iqx+1,iqy))
                    +  dy*(-x(iqx,iqy)+x(iqx,iqy+1))
                    +  dxdy*( x(iqx,iqy) - x(iqx,iqy+1)
                             -x(iqx+1,iqy) + x(iqx+1,iqy+1) );

               xb += dm(1);
               yb += dm(4);
            } // end for j
        } // end for i
     }
    else {
        fprintf(stderr, "bckpj3: nx must be greater than 2*ri\n");
    }

    return status;
}
int ChaoProjector::cb2sph ( float *  cube,
Vec3i  volsize,
int  ri,
Vec3i  origin,
int  nnz0,
int *  ptrs,
int *  cord,
float *  sphere 
) const [private]

Definition at line 1324 of file projector.cpp.

References cord, cube, nnz, nrays, nx, ny, ptrs, sphere, and status.

{
    int    xs, ys, zs, xx, yy, zz, rs, r2;
    int    ix, iy, iz, jnz, nnz, nrays;
    int    ftm = 0, status = 0;

    int xcent = (int)origin[0];
    int ycent = (int)origin[1];
    int zcent = (int)origin[2];

    int nx = (int)volsize[0];
    int ny = (int)volsize[1];
    int nz = (int)volsize[2];

    r2      = ri*ri;
    nnz     = 0;
    nrays    = 0;
    ptrs(1) = 1;

    for (ix = 1; ix <= nx; ix++) {
       xs  = ix-xcent;
       xx  = xs*xs;
       for ( iy = 1; iy <= ny; iy++ ) {
           ys = iy-ycent;
           yy = ys*ys;
           jnz = 0;

           ftm = 1;
           // not the most efficient implementation
           for (iz = 1; iz <= nz; iz++) {
               zs = iz-zcent;
               zz = zs*zs;
               rs = xx + yy + zz;
               if (rs <= r2) {
                  jnz++;
                  nnz++;
                  sphere(nnz) = cube(iz, iy, ix);

                  //  record the coordinates of the first nonzero ===
                  if (ftm) {
                     nrays++;
                     cord(1,nrays) = iz;
                     cord(2,nrays) = iy;
                     cord(3,nrays) = ix;
                     ftm = 0;
                  }
               }
            } // end for (iz..)
            if (jnz > 0) {
                ptrs(nrays+1) = ptrs(nrays) + jnz;
            }  // endif (jnz)
       } // end for iy
    } // end for ix
    if (nnz != nnz0) status = -1;
    return status;
}
int ChaoProjector::fwdpj3 ( Vec3i  volsize,
int  nray,
int  nnz,
float *  dm,
Vec3i  origin,
int  ri,
int *  ptrs,
int *  cord,
float *  x,
float *  y 
) const [private]

Definition at line 1418 of file projector.cpp.

References cord, dm, ifix(), nrays, nx, ptrs, status, x, and y.

{
    /*
        purpose:  y <--- proj(x)
        input  :  volsize  the size (nx,ny,nz) of the volume
                  nrays    number of rays within the compact spherical
                           representation
                  nnz      number of voxels within the sphere
                  dm       an array of size 9 storing transformation
                           associated with the projection direction
                  origin   coordinates of the center of the volume
                  ri       radius of the sphere
                  ptrs     the beginning address of each ray
                  cord     the coordinates of the first point in each ray
                  x        3d input volume
                  y        2d output image
    */

    int    iqx, iqy, i, j, xc, yc, zc;
    float  ct, dipx, dipy, dipx1m, dipy1m, xb, yb, dm1, dm4;
    int    status = 0;

    int xcent = origin[0];
    int ycent = origin[1];
    int zcent = origin[2];

    int nx = volsize[0];

    dm1 = dm(1);
    dm4 = dm(4);

    if ( nx > 2*ri ) {
        for (i = 1; i <= nrays; i++) {
            zc = cord(1,i)-zcent;
            yc = cord(2,i)-ycent;
            xc = cord(3,i)-xcent;

            xb = zc*dm(1)+yc*dm(2)+xc*dm(3) + xcent;
            yb = zc*dm(4)+yc*dm(5)+xc*dm(6) + ycent;

            for (j = ptrs(i); j< ptrs(i+1); j++) {
               iqx = ifix(xb);
               iqy = ifix(yb);

               ct   = x(j);
               dipx =  xb - (float)(iqx);
               dipy = (yb - (float)(iqy)) * ct;

               dipy1m = ct - dipy;
               dipx1m = 1.0f - dipx;

               y(iqx  ,iqy)   = y(iqx  ,iqy)   + dipx1m*dipy1m;
               y(iqx+1,iqy)   = y(iqx+1,iqy)   + dipx*dipy1m;
               y(iqx+1,iqy+1) = y(iqx+1,iqy+1) + dipx*dipy;
               y(iqx  ,iqy+1) = y(iqx  ,iqy+1) + dipx1m*dipy;

               xb += dm1;
               yb += dm4;
           }
        }
    }
    else {
        fprintf(stderr, " nx must be greater than 2*ri\n");
        exit(1);
    }
    return status;
}
string EMAN::ChaoProjector::get_desc ( ) const [inline, virtual]

Implements EMAN::Projector.

Definition at line 385 of file projector.h.

                {
                        return "Fast real space projection generation with Bi-Linear interpolation.";
                }
string EMAN::ChaoProjector::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 380 of file projector.h.

References NAME.

                {
                        return NAME;
                }
TypeDict EMAN::ChaoProjector::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 395 of file projector.h.

References EMAN::EMObject::FLOAT, EMAN::EMObject::FLOATARRAY, EMAN::EMObject::INT, EMAN::TypeDict::put(), 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("anglelist", EMObject::FLOATARRAY);
                        d.put("radius",    EMObject::FLOAT);
                        return d;
                }
int ChaoProjector::getnnz ( Vec3i  volsize,
int  ri,
Vec3i  origin,
int *  nray,
int *  nnz 
) const [private]

Definition at line 1264 of file projector.cpp.

References nnz, nrays, nx, ny, and status.

{
        int  ix, iy, iz, rs, r2, xs, ys, zs, xx, yy, zz;
        int  ftm=0, status = 0;

        r2    = ri*ri;
        *nnz  = 0;
        *nrays = 0;
        int nx = (int)volsize[0];
        int ny = (int)volsize[1];
        int nz = (int)volsize[2];

        int xcent = (int)origin[0];
        int ycent = (int)origin[1];
        int zcent = (int)origin[2];

        // need to add some error checking
        for (ix = 1; ix <=nx; ix++) {
            xs  = ix-xcent;
            xx  = xs*xs;
            for (iy = 1; iy <= ny; iy++) {
                ys = iy-ycent;
                yy = ys*ys;
                ftm = 1;
                for (iz = 1; iz <= nz; iz++) {
                    zs = iz-zcent;
                    zz = zs*zs;
                    rs = xx + yy + zz;
                    if (rs <= r2) {
                        (*nnz)++;
                        if (ftm) {
                           (*nrays)++;
                           ftm = 0;
                        }
                    }
                }
            } // end for iy
        } // end for ix
        return status;
}
int ChaoProjector::ifix ( float  a) const [private]

Definition at line 1567 of file projector.cpp.

{
    int ia;

    if (a>=0) {
       ia = (int)floor(a);
    }
    else {
       ia = (int)ceil(a);
    }
    return ia;
}
static Projector* EMAN::ChaoProjector::NEW ( ) [inline, static]

Definition at line 390 of file projector.h.

                {
                        return new ChaoProjector();
                }
EMData * ChaoProjector::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 1622 of file projector.cpp.

References anglelist, cb2sph(), cord, cube, dm, EMDeleteArray(), fwdpj3(), EMAN::EMData::get_data(), EMAN::Transform::get_rotation(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), getnnz(), images, LOGERR, nnz, nrays, NullPointerException, phi, ptrs, EMAN::EMData::set_attr(), EMAN::EMData::set_complex(), EMAN::EMData::set_ri(), EMAN::EMData::set_size(), sphere, status, theta, and EMAN::EMData::update().

{

        int nrays, nnz, status, j;
        float *dm;
        int   *ptrs, *cord;
        float *sphere, *images;

        int nxvol = vol->get_xsize();
        int nyvol = vol->get_ysize();
        int nzvol = vol->get_zsize();
        Vec3i volsize(nxvol,nyvol,nzvol);

        int dim = Util::get_min(nxvol,nyvol,nzvol);
        if (nzvol == 1) {
                LOGERR("The ChaoProjector needs a volume!");
                return 0;
        }
        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] = nxvol/2+1;}
        if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
        else {origin[1] = nyvol/2+1;}
        if (params.has_key("origin_z")) {origin[2] = params["origin_z"];}
        else {origin[2] = nzvol/2+1;}

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

        // retrieve the voxel values
        float *cube = vol->get_data();

        // count the number of voxels within a sphere centered at icent,
        // with radius ri
        status = getnnz(volsize, ri, origin, &nrays, &nnz);
        // need to check status...

        // convert from cube to sphere
        sphere = new float[nnz];
        ptrs   = new int[nrays+1];
        cord   = new int[3*nrays];
        if (sphere == NULL || ptrs == NULL || cord == NULL) {
                fprintf(stderr,"ChaoProjector::project3d, failed to allocate!\n");
                exit(1);
        }
        for (int i = 0; i<nnz; i++) sphere[i] = 0.0;
        for (int i = 0; i<nrays+1; i++) ptrs[i] = 0;
        for (int i = 0; i<3*nrays; i++) cord[i] = 0;

        status = cb2sph(cube, volsize, ri, origin, nnz, ptrs, cord, sphere);
        // check status

        int nangles = 0;
        vector<float> anglelist;
        string angletype = "SPIDER";
        // 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;}

                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

        dm = new float[nangles*9];
        setdm(anglelist, angletype, dm);

                // return images
        EMData *ret = new EMData();
        ret->set_size(nxvol, nyvol, nangles);
        ret->set_complex(false);
        ret->set_ri(true);

        images = ret->get_data();

        for (j = 1; j <= nangles; j++) {
                status = fwdpj3(volsize, nrays, nnz   , &dm(1,j), origin, ri,
                                                ptrs   ,  cord, sphere, &images(1,1,j));
        // check status?
        }

        // deallocate all temporary work space
        EMDeleteArray(dm);
        EMDeleteArray(ptrs);
        EMDeleteArray(cord);
        EMDeleteArray(sphere);

        if (!params.has_key("anglelist")) {
                Transform* t3d = params["transform"];
                ret->set_attr("xform.projection",t3d);
                if(t3d) {delete t3d; t3d=0;}
        }
        ret->update();
        return ret;
}
void ChaoProjector::setdm ( vector< float >  anglelist,
string const  angletype,
float *  dm 
) const [private]

Definition at line 1584 of file projector.cpp.

References anglelist, dgr_to_rad, dm, phi, and theta.

{ // convert Euler angles to transformations, dm is an 9 by nangles array

        float  psi, theta, phi;
        double cthe, sthe, cpsi, spsi, cphi, sphi;
        int    j;

        int nangles = anglelist.size() / 3;

        // now convert all angles
        for (j = 1; j <= nangles; j++) {
                phi   = static_cast<float>(anglelist(1,j)*dgr_to_rad);
                theta = static_cast<float>(anglelist(2,j)*dgr_to_rad);
                psi   = static_cast<float>(anglelist(3,j)*dgr_to_rad);

                //              cout << phi << " " << theta << " " << psi << endl;
                cthe  = cos(theta);
                sthe  = sin(theta);
                cpsi  = cos(psi);
                spsi  = sin(psi);
                cphi  = cos(phi);
                sphi  = sin(phi);

                dm(1,j)=static_cast<float>(cphi*cthe*cpsi-sphi*spsi);
                dm(2,j)=static_cast<float>(sphi*cthe*cpsi+cphi*spsi);
                dm(3,j)=static_cast<float>(-sthe*cpsi);
                dm(4,j)=static_cast<float>(-cphi*cthe*spsi-sphi*cpsi);
                dm(5,j)=static_cast<float>(-sphi*cthe*spsi+cphi*cpsi);
                dm(6,j)=static_cast<float>(sthe*spsi);
                dm(7,j)=static_cast<float>(sthe*cphi);
                dm(8,j)=static_cast<float>(sthe*sphi);
                dm(9,j)=static_cast<float>(cthe);
        }
}
int ChaoProjector::sph2cb ( float *  sphere,
Vec3i  volsize,
int  nray,
int  ri,
int  nnz0,
int *  ptrs,
int *  cord,
float *  cube 
) const [private]

Definition at line 1383 of file projector.cpp.

References cord, cube, nnz, nrays, nx, ny, ptrs, sphere, and status.

{
    int       status=0;
    int       r2, i, j, ix, iy, iz,  nnz;

    int nx = (int)volsize[0];
    int ny = (int)volsize[1];
    // int nz = (int)volsize[2];

    r2      = ri*ri;
    nnz     = 0;
    ptrs(1) = 1;

    // no need to initialize
    // for (i = 0; i<nx*ny*nz; i++) cube[i]=0.0;

    nnz = 0;
    for (j = 1; j <= nrays; j++) {
       iz = cord(1,j);
       iy = cord(2,j);
       ix = cord(3,j);
       for (i = ptrs(j); i<=ptrs(j+1)-1; i++, iz++) {
           nnz++;
           cube(iz,iy,ix) = sphere(nnz);
       }
    }
    if (nnz != nnz0) status = -1;
    return status;
}

Member Data Documentation

const string ChaoProjector::NAME = "chao" [static]

Definition at line 407 of file projector.h.

Referenced by get_name().


The documentation for this class was generated from the following files: