EMAN2
project3d.cpp
Go to the documentation of this file.
00001 #include "mpi.h"
00002 
00003 #include "emdata.h"
00004 #include "project3d.h"
00005 
00006 #define cube(i,j,k) cube[ (((k)-1)*ny + (j)-1)*nx + (i)-1 ] 
00007 #define sphere(i)   sphere[(i)-1]
00008 #define cord(i,j)   cord[((j)-1)*3 + (i) -1]
00009 #define ptrs(i)     ptrs[(i)-1]
00010 #define dm(i)       dm[(i)-1]
00011 
00012 int cb2sph(float *cube, Vec3i volsize, int    ri, Vec3i origin, 
00013            int    nnz0, int     *ptrs, int *cord, float *sphere) 
00014 {
00015     int    xs, ys, zs, xx, yy, zz, rs, r2;
00016     int    ix, iy, iz, jnz, nnz, nrays;
00017     int    ftm = 0, status = 0;  
00018 
00019     int xcent = (int)origin[0];
00020     int ycent = (int)origin[1];
00021     int zcent = (int)origin[2];
00022 
00023     int nx = (int)volsize[0];
00024     int ny = (int)volsize[1];
00025     int nz = (int)volsize[2];
00026 
00027     r2      = ri*ri;
00028     nnz     = 0;
00029     nrays    = 0;
00030     ptrs(1) = 1;
00031 
00032     for (ix = 1; ix <= nx; ix++) {
00033        xs  = ix-xcent;
00034        xx  = xs*xs;
00035        for ( iy = 1; iy <= ny; iy++ ) {
00036            ys = iy-ycent;
00037            yy = ys*ys;
00038            jnz = 0;
00039 
00040            ftm = 1;
00041            // not the most efficient implementation
00042            for (iz = 1; iz <= nz; iz++) {
00043                zs = iz-zcent;
00044                zz = zs*zs;
00045                rs = xx + yy + zz;
00046                if (rs <= r2) {
00047                   jnz++;
00048                   nnz++;
00049                   sphere(nnz) = cube(iz, iy, ix); 
00050 
00051                   //  record the coordinates of the first nonzero ===
00052                   if (ftm) {
00053                      nrays++;
00054                      cord(1,nrays) = iz; 
00055                      cord(2,nrays) = iy; 
00056                      cord(3,nrays) = ix;
00057                      ftm = 0;
00058                   }
00059                }
00060             } // end for (iz..)
00061             if (jnz > 0) {
00062                 ptrs(nrays+1) = ptrs(nrays) + jnz;
00063             }  // endif (jnz)
00064        } // end for iy
00065     } // end for ix
00066     if (nnz != nnz0) status = -1;
00067     return status;
00068 }
00069 
00070 // decompress sphere into a cube
00071 int sph2cb(float *sphere, Vec3i volsize, int  nrays, int    ri, 
00072            int      nnz0, int     *ptrs, int  *cord, float *cube)
00073 {
00074     int       status=0;
00075     int       r2, i, j, ix, iy, iz,  nnz;
00076 
00077     int nx = (int)volsize[0];
00078     int ny = (int)volsize[1];
00079     // int nz = (int)volsize[2];
00080 
00081     r2      = ri*ri;
00082     nnz     = 0;
00083     ptrs(1) = 1;
00084 
00085     // no need to initialize
00086     // for (i = 0; i<nx*ny*nz; i++) cube[i]=0.0;
00087 
00088     nnz = 0;
00089     for (j = 1; j <= nrays; j++) {
00090        iz = cord(1,j);
00091        iy = cord(2,j);
00092        ix = cord(3,j);
00093        for (i = ptrs(j); i<=ptrs(j+1)-1; i++, iz++) {
00094            nnz++;
00095            cube(iz,iy,ix) = sphere(nnz);
00096        }
00097     }
00098     if (nnz != nnz0) status = -1;
00099     return status;
00100 }
00101 
00102 #define x(i)        x[(i)-1]
00103 #define y(i,j)      y[((j)-1)*nx + (i) - 1]
00104 
00105 // project from 3D to 2D (single image)
00106 int fwdpj3(Vec3i volsize, int nrays, int   nnz, float *dm, 
00107            Vec3i  origin, int    ri, int *ptrs, int *cord, 
00108            float      *x, float  *y)
00109 {
00110     /*
00111         purpose:  y <--- proj(x)
00112         input  :  volsize  the size (nx,ny,nz) of the volume
00113                   nrays    number of rays within the compact spherical 
00114                            representation
00115                   nnz      number of voxels within the sphere
00116                   dm       an array of size 9 storing transformation 
00117                            associated with the projection direction
00118                   origin   coordinates of the center of the volume
00119                   ri       radius of the sphere
00120                   ptrs     the beginning address of each ray
00121                   cord     the coordinates of the first point in each ray
00122                   x        3d input volume
00123                   y        2d output image 
00124     */
00125 
00126     int    iqx, iqy, i, j, xc, yc, zc;
00127     float  ct, dipx, dipy, dipx1m, dipy1m, xb, yb, dm1, dm4;
00128     int    status = 0;
00129     
00130     // Phi: adding the shift parameters that get passed in as the last two entries of dm
00131     float sx, sy;
00132 
00133     sx = dm(7);
00134     sy = dm(8);
00135 
00136     int xcent = origin[0];
00137     int ycent = origin[1];
00138     int zcent = origin[2];
00139 
00140     int nx = volsize[0];
00141     int ny = volsize[1];
00142 
00143     dm1 = dm(1);
00144     dm4 = dm(4);
00145  
00146     if ( nx > 2*ri ) {
00147         for (i = 1; i <= nrays; i++) {
00148 
00149             zc = cord(1,i)-zcent;
00150             yc = cord(2,i)-ycent;
00151             xc = cord(3,i)-xcent;
00152             xb = zc* dm(1) +yc* dm(2) +xc* dm(3) + xcent + sx;
00153             yb = zc* dm(4) +yc* dm(5) +xc* dm(6) + ycent + sy;
00154 
00155             for (j = ptrs(i); j< ptrs(i+1); j++) {
00156                iqx = ifix(xb);
00157                iqy = ifix(yb);
00158 
00159                ct   = x(j);
00160 
00161                // dipx =  xb - (float)(iqx);
00162                // dipy = (yb - (float)(iqy)) * ct;
00163                    dipx =  xb - iqx;
00164                    dipy = (yb - iqy) * ct;
00165 
00166                dipy1m = ct - dipy;
00167                dipx1m = 1.0 - dipx;
00168 
00169                         if (iqx <= nx && iqy <= ny && iqx >= 1 && iqy >= 1) 
00170                // y(iqx  ,iqy)   = y(iqx  ,iqy)   + dipx1m*dipy1m;
00171                y(iqx  ,iqy)   +=  dipx1m*dipy1m;
00172                         if (iqx + 1 <= nx && iqy <= ny && iqx >= 0 && iqy >= 1) 
00173                // y(iqx+1,iqy)   = y(iqx+1,iqy)   + dipx*dipy1m; 
00174                y(iqx+1,iqy)   +=  dipx*dipy1m; 
00175                         if (iqx + 1 <= nx && iqy + 1 <= ny && iqx >= 0 && iqy >= 0) 
00176                // y(iqx+1,iqy+1) = y(iqx+1,iqy+1) + dipx*dipy;         
00177                y(iqx+1,iqy+1) +=  dipx*dipy;         
00178                         if (iqx <= nx && iqy + 1 <= ny && iqx >= 1 && iqy >= 0) 
00179                // y(iqx  ,iqy+1) = y(iqx  ,iqy+1) + dipx1m*dipy;
00180                y(iqx  ,iqy+1) +=  dipx1m*dipy;
00181                xb += dm1;
00182                yb += dm4;
00183            }
00184         }
00185     }
00186     else {
00187         fprintf(stderr, " nx must be greater than 2*ri\n");
00188         exit(1);
00189     }
00190     return status;
00191 }
00192 #undef x
00193 #undef y
00194 
00195 #define y(i)        y[(i)-1]
00196 #define x(i,j)      x[((j)-1)*nx + (i) - 1]
00197 
00198 // backproject from 2D to 3D for a single image
00199 int bckpj3(Vec3i volsize, int nrays, int   nnz, float *dm, 
00200            Vec3i  origin, int    ri, int *ptrs, int *cord, 
00201            float      *x, float *y)
00202 {
00203     int       i, j, iqx,iqy, xc, yc, zc;
00204     float     xb, yb, dx, dy, dx1m, dy1m, dxdy;
00205     int       status = 0; 
00206 
00207     int xcent = origin[0];
00208     int ycent = origin[1];
00209     int zcent = origin[2];
00210 
00211     int nx = volsize[0];
00212     int ny = volsize[1];
00213 
00214     // Phi: adding the shift parameters that get passed in as the last two entries of dm
00215     float sx, sy;
00216 
00217     sx = dm(7);
00218     sy = dm(8);
00219 
00220 
00221     if ( nx > 2*ri) {
00222         for (i = 1; i <= nrays; i++) {
00223             zc = cord(1,i) - zcent;
00224             yc = cord(2,i) - ycent;
00225             xc = cord(3,i) - xcent;
00226 
00227             xb = zc*dm(1)+yc*dm(2)+xc*dm(3) + xcent + sx;
00228             yb = zc*dm(4)+yc*dm(5)+xc*dm(6) + ycent + sy;
00229 
00230             for (j = ptrs(i); j <ptrs(i+1); j++) {
00231                 iqx = ifix(xb);
00232                 iqy = ifix(yb);
00233 
00234                 dx = xb - iqx;
00235                 dy = yb - iqy;
00236                 dx1m = 1.0 - dx;
00237                 dy1m = 1.0 - dy;
00238                 dxdy = dx*dy;
00239 /*
00240 c               y(j) = y(j) + dx1m*dy1m*x(iqx  , iqy)
00241 c     &                     + dx1m*dy  *x(iqx  , iqy+1)
00242 c     &                     + dx  *dy1m*x(iqx+1, iqy)
00243 c     &                     + dx  *dy  *x(iqx+1, iqy+1)  
00244 c
00245 c              --- faster version of the above commented out
00246 c                  code (derived by summing the following table 
00247 c                  of coefficients along  the colunms) ---
00248 c
00249 c                        1         dx        dy      dxdy
00250 c                     ------   --------  --------  -------
00251 c                      x(i,j)   -x(i,j)   -x(i,j)    x(i,j)  
00252 c                                        x(i,j+1) -x(i,j+1)
00253 c                              x(i+1,j)           -x(i+1,j)
00254 c                                                x(i+1,j+1) 
00255 c
00256 */
00257                 // Phi: add index checking, now that shifts are being used
00258                 if ( iqx <= nx && iqy <= ny && iqx >= 1 && iqy >= 1 ) {
00259                     y(j) += x(iqx,iqy);
00260                     if ( iqx + 1 <= nx && iqx + 1 >= 1 ) {
00261                         y(j) += dx*(-x(iqx,iqy)+x(iqx+1,iqy));
00262                     }
00263                     if ( iqy + 1 <= ny && iqy + 1 >= 1 ) {
00264                         y(j) += dy*(-x(iqx,iqy)+x(iqx,iqy+1));
00265                     }
00266                     if ( iqx + 1 <= nx && iqy + 1 <= ny && iqx + 1 >= 1 && iqy + 1 >= 1 ) {
00267                         y(j) += dxdy*( x(iqx,iqy) - x(iqx,iqy+1) -x(iqx+1,iqy) + x(iqx+1,iqy+1) );
00268                     }
00269                 }
00270 
00271 //                y(j) += x(iqx,iqy)
00272 //                     +  dx*(-x(iqx,iqy)+x(iqx+1,iqy))
00273 //                     +  dy*(-x(iqx,iqy)+x(iqx,iqy+1))
00274 //                     +  dxdy*( x(iqx,iqy) - x(iqx,iqy+1) 
00275 //                              -x(iqx+1,iqy) + x(iqx+1,iqy+1) );
00276 
00277                xb += dm(1);
00278                yb += dm(4);
00279             } // end for j
00280         } // end for i
00281      }
00282     else {
00283         fprintf(stderr, "bckpj3: nx must be greater than 2*ri\n");
00284     }
00285 
00286     return status;
00287 }
00288 
00289 #undef x
00290 #undef y
00291 #undef dm
00292 
00293 int getnnz(Vec3i volsize, int ri, Vec3i origin, int *nrays, int *nnz) 
00294 /*
00295    purpose: count the number of voxels within a sphere centered
00296             at origin and with a radius ri.
00297 
00298      input:
00299      volsize contains the size information (nx,ny,nz) about the volume
00300      ri      radius of the object embedded in the cube.
00301      origin  coordinates for the center of the volume
00302    
00303      output:
00304      nnz    total number of voxels within the sphere (of radius ri)
00305      nrays  number of rays in z-direction. 
00306 */
00307 {
00308     int  ix, iy, iz, rs, r2, xs, ys, zs, xx, yy, zz;
00309     int  ftm=0, status = 0;
00310 
00311     r2    = ri*ri;
00312     *nnz  = 0;
00313     *nrays = 0;
00314     int nx = volsize[0];
00315     int ny = volsize[1];
00316     int nz = volsize[2];
00317     // int nx = (int)volsize[0];
00318     // int ny = (int)volsize[1];
00319     // int nz = (int)volsize[2];
00320 
00321     int xcent = origin[0]; 
00322     int ycent = origin[1]; 
00323     int zcent = origin[2]; 
00324     // int xcent = (int)origin[0]; 
00325     // int ycent = (int)origin[1]; 
00326     // int zcent = (int)origin[2]; 
00327 
00328     // need to add some error checking 
00329     for (ix = 1; ix <=nx; ix++) { 
00330         xs  = ix-xcent;
00331         xx  = xs*xs;
00332         for (iy = 1; iy <= ny; iy++) {
00333             ys = iy-ycent;
00334             yy = ys*ys; 
00335             ftm = 1;
00336             for (iz = 1; iz <= nz; iz++) {
00337                 zs = iz-zcent;
00338                 zz = zs*zs;
00339                 rs = xx + yy + zz;
00340                 if (rs <= r2) {
00341                     (*nnz)++;
00342                     if (ftm) {
00343                        (*nrays)++;
00344                        ftm = 0;
00345                     }
00346                 }
00347             }
00348         } // end for iy
00349     } // end for ix
00350     return status;
00351 }
00352 
00353 // This might be a convenient function to use in fgcalc, rather than repeating the same code six times
00354 int make_proj_mat(float phi, float theta, float psi, float * dm)
00355 {
00356 //     float cphi=cos(phi);
00357 //     float sphi=sin(phi);
00358 //     float cthe=cos(theta);
00359 //     float sthe=sin(theta);
00360 //     float cpsi=cos(psi);
00361 //     float spsi=sin(psi);
00362 
00363     double cphi, sphi, cthe, sthe, cpsi, spsi;
00364     double dphi = phi;
00365     double dthe = theta;
00366     double dpsi = psi;
00367     sincos(dphi, &sphi, &cphi);
00368     sincos(dthe, &sthe, &cthe);
00369     sincos(dpsi, &spsi, &cpsi);
00370 
00371     dm[0]=cphi*cthe*cpsi-sphi*spsi;
00372     dm[1]=sphi*cthe*cpsi+cphi*spsi;
00373     dm[2]=-sthe*cpsi;
00374     dm[3]=-cphi*cthe*spsi-sphi*cpsi;
00375     dm[4]=-sphi*cthe*spsi+cphi*cpsi;
00376     dm[5]=sthe*spsi;
00377 
00378     return 0;
00379 }