EMAN2
Functions
project3d.h File Reference
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

int fwdpj3 (Vec3i volsize, int nrays, int nnz, float *dm, Vec3i origin, int ri, int *ptrs, int *cord, float *x, float *y)
int bckpj3 (Vec3i volsize, int nrays, int nnz, float *dm, Vec3i origin, int ri, int *ptrs, int *cord, float *x, float *y)
int sph2cb (float *sphere, Vec3i volsize, int nrays, int ri, int nnz0, int *ptrs, int *cord, float *cube)
int cb2sph (float *cube, Vec3i volsize, int ri, Vec3i origin, int nnz0, int *ptrs, int *cord, float *sphere)
int getnnz (Vec3i volsize, int ri, Vec3i origin, int *nrays, int *nnz)
int ifix (float a)
int make_proj_mat (float phi, float theta, float psi, float *dm)

Function Documentation

int bckpj3 ( Vec3i  volsize,
int  nrays,
int  nnz,
float *  dm,
Vec3i  origin,
int  ri,
int *  ptrs,
int *  cord,
float *  x,
float *  y 
)

Definition at line 199 of file project3d.cpp.

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

Referenced by EMAN::ChaoProjector::backproject3d(), fgcalc(), main(), and recons3d_sirt_mpi().

{
    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];
    int ny = volsize[1];

    // Phi: adding the shift parameters that get passed in as the last two entries of dm
    float sx, sy;

    sx = dm(7);
    sy = dm(8);


    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 + sx;
            yb = zc*dm(4)+yc*dm(5)+xc*dm(6) + ycent + sy;

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

                dx = xb - iqx;
                dy = yb - iqy;
                dx1m = 1.0 - dx;
                dy1m = 1.0 - 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
*/
                // Phi: add index checking, now that shifts are being used
                if ( iqx <= nx && iqy <= ny && iqx >= 1 && iqy >= 1 ) {
                    y(j) += x(iqx,iqy);
                    if ( iqx + 1 <= nx && iqx + 1 >= 1 ) {
                        y(j) += dx*(-x(iqx,iqy)+x(iqx+1,iqy));
                    }
                    if ( iqy + 1 <= ny && iqy + 1 >= 1 ) {
                        y(j) += dy*(-x(iqx,iqy)+x(iqx,iqy+1));
                    }
                    if ( iqx + 1 <= nx && iqy + 1 <= ny && iqx + 1 >= 1 && iqy + 1 >= 1 ) {
                        y(j) += dxdy*( x(iqx,iqy) - x(iqx,iqy+1) -x(iqx+1,iqy) + x(iqx+1,iqy+1) );
                    }
                }

//                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 cb2sph ( float *  cube,
Vec3i  volsize,
int  ri,
Vec3i  origin,
int  nnz0,
int *  ptrs,
int *  cord,
float *  sphere 
)

Definition at line 12 of file project3d.cpp.

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

Referenced by ali3d_d(), EMAN::ChaoProjector::backproject3d(), EMAN::ChaoProjector::project3d(), recons3d_sirt_mpi(), and unified().

{
    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 fwdpj3 ( Vec3i  volsize,
int  nrays,
int  nnz,
float *  dm,
Vec3i  origin,
int  ri,
int *  ptrs,
int *  cord,
float *  x,
float *  y 
)

Definition at line 106 of file project3d.cpp.

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

Referenced by fcalc(), fgcalc(), main(), EMAN::ChaoProjector::project3d(), and recons3d_sirt_mpi().

{
    /*
        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;
    
    // Phi: adding the shift parameters that get passed in as the last two entries of dm
    float sx, sy;

    sx = dm(7);
    sy = dm(8);

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

    int nx = volsize[0];
    int ny = volsize[1];

    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 + sx;
            yb = zc* dm(4) +yc* dm(5) +xc* dm(6) + ycent + sy;

            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;
                   dipx =  xb - iqx;
                   dipy = (yb - iqy) * ct;

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

                        if (iqx <= nx && iqy <= ny && iqx >= 1 && iqy >= 1) 
               // y(iqx  ,iqy)   = y(iqx  ,iqy)   + dipx1m*dipy1m;
               y(iqx  ,iqy)   +=  dipx1m*dipy1m;
                        if (iqx + 1 <= nx && iqy <= ny && iqx >= 0 && iqy >= 1) 
               // y(iqx+1,iqy)   = y(iqx+1,iqy)   + dipx*dipy1m; 
               y(iqx+1,iqy)   +=  dipx*dipy1m; 
                        if (iqx + 1 <= nx && iqy + 1 <= ny && iqx >= 0 && iqy >= 0) 
               // y(iqx+1,iqy+1) = y(iqx+1,iqy+1) + dipx*dipy;         
               y(iqx+1,iqy+1) +=  dipx*dipy;         
                        if (iqx <= nx && iqy + 1 <= ny && iqx >= 1 && iqy >= 0) 
               // y(iqx  ,iqy+1) = y(iqx  ,iqy+1) + dipx1m*dipy;
               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;
}
int getnnz ( Vec3i  volsize,
int  ri,
Vec3i  origin,
int *  nrays,
int *  nnz 
)

Definition at line 293 of file project3d.cpp.

References nx, ny, and status.

Referenced by ali3d_d(), EMAN::ChaoProjector::backproject3d(), main(), EMAN::ChaoProjector::project3d(), recons3d_CGLS_mpi_Cart(), recons3d_HyBR_mpi_Cart(), recons3d_sirt_mpi(), recons3d_sirt_mpi_Cart(), and unified().

{
    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 = volsize[0];
    int ny = volsize[1];
    int nz = volsize[2];
    // int nx = (int)volsize[0];
    // int ny = (int)volsize[1];
    // int nz = (int)volsize[2];

    int xcent = origin[0]; 
    int ycent = origin[1]; 
    int zcent = origin[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 ifix ( float  a)

Definition at line 41 of file utilnum.cpp.

Referenced by EMAN::ChaoProjector::bckpj3(), bckpj3(), bckpj3_Cart(), EMAN::ChaoProjector::fwdpj3(), fwdpj3(), and fwdpj3_Cart().

{
    int ia = 0;
    if (a >= 0.0) {
       ia = (int)floor(a);
    }
    else {
       ia = (int)ceil(a);
    }
    return ia;
}
int make_proj_mat ( float  phi,
float  theta,
float  psi,
float *  dm 
)

Definition at line 354 of file project3d.cpp.

{
//     float cphi=cos(phi);
//     float sphi=sin(phi);
//     float cthe=cos(theta);
//     float sthe=sin(theta);
//     float cpsi=cos(psi);
//     float spsi=sin(psi);

    double cphi, sphi, cthe, sthe, cpsi, spsi;
    double dphi = phi;
    double dthe = theta;
    double dpsi = psi;
    sincos(dphi, &sphi, &cphi);
    sincos(dthe, &sthe, &cthe);
    sincos(dpsi, &spsi, &cpsi);

    dm[0]=cphi*cthe*cpsi-sphi*spsi;
    dm[1]=sphi*cthe*cpsi+cphi*spsi;
    dm[2]=-sthe*cpsi;
    dm[3]=-cphi*cthe*spsi-sphi*cpsi;
    dm[4]=-sphi*cthe*spsi+cphi*cpsi;
    dm[5]=sthe*spsi;

    return 0;
}
int sph2cb ( float *  sphere,
Vec3i  volsize,
int  nrays,
int  ri,
int  nnz0,
int *  ptrs,
int *  cord,
float *  cube 
)

Definition at line 71 of file project3d.cpp.

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

Referenced by EMAN::ChaoProjector::backproject3d(), recons3d_CGLS_mpi_Cart(), recons3d_HyBR_mpi_Cart(), recons3d_sirt_mpi(), recons3d_sirt_mpi_Cart(), and unified().

{
    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;
}