EMAN2
project3d_Cart.cpp
Go to the documentation of this file.
00001 #include "mpi.h"
00002 #include "stdlib.h"
00003 #include "emdata.h"
00004 
00005 #include "utilcomm_Cart.h"
00006 #include "project3d_Cart.h"
00007 #include "project3d.h"
00008 
00009 using namespace EMAN;
00010 
00011 #define cord(i,j)   cord[((j)-1)*3 + (i) -1]
00012 #define ptrs(i)     ptrs[(i)-1]
00013 #define dm(i)       dm[(i)-1]
00014 
00015 int getcb2sph(Vec3i volsize, int ri, Vec3i origin, int nnz0, int *ptrs, int *cord) 
00016 /* This function is similar to 'cb2sph' except we don't actually deal with the cube and sphere data. We only use this function to get vectors ptrs and coord.
00017 
00018         Input:
00019          volsize - dimensions of 3D cube
00020               ri - radius of sphere
00021           origin - origin of sphere
00022             nnz0 - number of non-zeros
00023 
00024         Output:
00025               ptrs - vector containing pointers to the beginning of each ray
00026               cord - vector containing the x,y,z coordinates of the beginning of each ray
00027 */
00028 
00029 {
00030     int    xs, ys, zs, xx, yy, zz, rs, r2;
00031     int    ix, iy, iz, jnz, nnz, nrays;
00032     int    ftm = 0, status = 0;  
00033 
00034     int xcent = (int)origin[0];
00035     int ycent = (int)origin[1];
00036     int zcent = (int)origin[2];
00037 
00038     int nx = (int)volsize[0];
00039     int ny = (int)volsize[1];
00040     int nz = (int)volsize[2];
00041 
00042     r2      = ri*ri;
00043     nnz     = 0;
00044     nrays    = 0;
00045     ptrs(1) = 1;
00046 
00047     for (ix = 1; ix <= nx; ix++) {
00048        xs  = ix-xcent;
00049        xx  = xs*xs;
00050        for ( iy = 1; iy <= ny; iy++ ) {
00051            ys = iy-ycent;
00052            yy = ys*ys;
00053            jnz = 0;
00054 
00055            ftm = 1;
00056            // not the most efficient implementation
00057            for (iz = 1; iz <= nz; iz++) {
00058                zs = iz-zcent;
00059                zz = zs*zs;
00060                rs = xx + yy + zz;
00061                if (rs <= r2) {
00062                   jnz++;
00063                   nnz++;
00064      //             sphere(nnz) = cube(iz, iy, ix); 
00065 
00066                   //  record the coordinates of the first nonzero ===
00067                   if (ftm) {
00068                      nrays++;
00069                      cord(1,nrays) = iz; 
00070                      cord(2,nrays) = iy; 
00071                      cord(3,nrays) = ix;
00072                      ftm = 0;
00073                   }
00074                }
00075             } // end for (iz..)
00076             if (jnz > 0) {
00077                 ptrs(nrays+1) = ptrs(nrays) + jnz;
00078             }  // endif (jnz)
00079        } // end for iy
00080     } // end for ix
00081     if (nnz != nnz0) status = -1;
00082     return status;
00083 }
00084 
00085 int sphpart(MPI_Comm comm_2d, int nrays, int *ptrs, int *nnzbase, int *ptrstart)
00086 /* This function divides the volume among the column processors by computing nraysloc for each processor and the starting pointer for the first ray.
00087 
00088         Input:
00089             comm_2d - Cartesian communicator
00090               nrays - total number of rays
00091                ptrs - vector containing pointers
00092             nnzbase - ideal volume partition of nnz
00093         Output:
00094            ptrstart - vector containing all the starting pointers for each column processor group
00095            nraysloc - actual number of local rays
00096 */
00097 {
00098         int ROW = 0, COL = 1;
00099         int dims[2], periods[2], mycoords[2];
00100         int nraysloc = 0;
00101         
00102         // Get information associated with comm_2d
00103         MPI_Cart_get(comm_2d, 2, dims, periods, mycoords);
00104 
00105         int count = 1;
00106         if (mycoords[COL] == 0){ // First column starts out with the first ray
00107           nraysloc++;
00108         }
00109         ptrstart[0] = 0;
00110 
00111         for(int i=1; i<nrays ; i++){
00112           if (ptrs[i]-1 <= nnzbase[count] && ptrs[i+1]-1 >= nnzbase[count]){ 
00113                 //nnzbase is between or equal to ptrs
00114 
00115             if (nnzbase[count]-(ptrs[i]-1)>= ptrs[i+1]-1-nnzbase[count]){
00116               if(mycoords[COL] == count-1){  // ray belongs to count-1
00117                 nraysloc++;
00118               }
00119               ptrstart[count] = i+1;
00120               count++;
00121 
00122             } else { //nnzbase[count]-(ptrs[i]-1)>= ptrs[i+1]-1-nnzbase[count]
00123                 if(mycoords[COL] == count){// ray belongs to count and it's a new starting ray
00124                   nraysloc++;
00125                 }
00126                 ptrstart[count] = i;
00127                 count++;
00128             }
00129 
00130         }
00131         else {  //ptrs[i]-1 > nnzbase[count] so ray belongs to count-1
00132           if(mycoords[COL] == count-1){
00133             nraysloc++;
00134           }
00135 
00136         }
00137         } // end for
00138         ptrstart[dims[COL]] = nrays;
00139         return nraysloc;
00140 
00141 }
00142 
00143 
00144 #define x(i)        x[(i)-1]
00145 #define y(i,j)      y[((j)-1)*nx + (i) - 1]
00146 
00147 int fwdpj3_Cart(Vec3i volsize, int nraysloc, int nnzloc, float *dm, Vec3i origin, int ri, int *ptrs, int *cord, int myptrstart, float *x, float *y)
00148 /* This function will project from 3D to 2D for a single image for one portion of the 3D volume.  This function is to be run on every processor. y = P*x
00149 
00150         Input:
00151           volsize - size (nx,ny,nz) of 3D cube volume
00152           nraysloc - local number of rays within the compact spherical representation
00153           nnzloc - local number of voxels
00154           dm - an array of size 9 storing transformation 
00155           origin, ri - origin and radius of sphere
00156           ptrs, cord - pointers and coordinates for first ray
00157           myptrstart - determines starting rays
00158           x - portion of the 3D input volume
00159         Output:
00160           y - 2D output image
00161 */
00162 {
00163     int    iqx, iqy, i, j, xc, yc, zc, jj;
00164     float  ct, dipx, dipy, dipx1m, dipy1m, xb, yb, dm1, dm4;
00165     int    status = 0;
00166     
00167     // Phi: adding the shift parameters that get passed in as the last two entries of dm
00168     float sx, sy;
00169 
00170     sx = dm(7);
00171     sy = dm(8);
00172 
00173     int xcent = origin[0];
00174     int ycent = origin[1];
00175     int zcent = origin[2];
00176 
00177     int nx = volsize[0];
00178     int ny = volsize[1];
00179 
00180     dm1 = dm(1);
00181     dm4 = dm(4);
00182  
00183     if ( nx > 2*ri ) {
00184         for (i = 1; i <= nraysloc; i++) {
00185 
00186             zc = cord(1,myptrstart+i)-zcent;
00187             yc = cord(2,myptrstart+i)-ycent;
00188             xc = cord(3,myptrstart+i)-xcent;
00189             xb = zc* dm(1) +yc* dm(2) +xc* dm(3) + xcent + sx;
00190             yb = zc* dm(4) +yc* dm(5) +xc* dm(6) + ycent + sy;
00191 
00192             for (j = ptrs(myptrstart+i); j< ptrs(myptrstart+i+1); j++) {
00193                jj = j-ptrs[myptrstart]+1;
00194                iqx = ifix(xb);
00195                iqy = ifix(yb);
00196 
00197                 // jj is the index on the local volume
00198                ct   = x(jj);
00199 
00200                // dipx =  xb - (float)(iqx);
00201                // dipy = (yb - (float)(iqy)) * ct;
00202                    dipx =  xb - iqx;
00203                    dipy = (yb - iqy) * ct;
00204 
00205                dipy1m = ct - dipy;
00206                dipx1m = 1.0 - dipx;
00207 
00208                         if (iqx <= nx && iqy <= ny && iqx >= 1 && iqy >= 1) 
00209                // y(iqx  ,iqy)   = y(iqx  ,iqy)   + dipx1m*dipy1m;
00210                y(iqx  ,iqy)   +=  dipx1m*dipy1m;
00211                         if (iqx + 1 <= nx && iqy <= ny && iqx >= 0 && iqy >= 1) 
00212                // y(iqx+1,iqy)   = y(iqx+1,iqy)   + dipx*dipy1m; 
00213                y(iqx+1,iqy)   +=  dipx*dipy1m; 
00214                         if (iqx + 1 <= nx && iqy + 1 <= ny && iqx >= 0 && iqy >= 0) 
00215                // y(iqx+1,iqy+1) = y(iqx+1,iqy+1) + dipx*dipy;         
00216                y(iqx+1,iqy+1) +=  dipx*dipy;         
00217                         if (iqx <= nx && iqy + 1 <= ny && iqx >= 1 && iqy >= 0) 
00218                // y(iqx  ,iqy+1) = y(iqx  ,iqy+1) + dipx1m*dipy;
00219                y(iqx  ,iqy+1) +=  dipx1m*dipy;
00220                xb += dm1;
00221                yb += dm4;
00222            }
00223         }
00224     }
00225     else {
00226         fprintf(stderr, " nx must be greater than 2*ri\n");
00227         exit(1);
00228     }
00229     return status;
00230 }
00231 #undef x
00232 #undef y
00233 
00234 #define y(i)        y[(i)-1]
00235 #define x(i,j)      x[((j)-1)*nx + (i) - 1]
00236 
00237 int bckpj3_Cart(Vec3i volsize, int nraysloc, int nnzloc, float *dm, Vec3i origin, int ri, int *ptrs, int *cord, int myptrstart, float *x, float *y)
00238 /* This function will backproject from 2D to 3D for a single image for one portion of the 3D volume.  This function is to be run on every processor. y = P'x
00239 
00240         Input:
00241           volsize - size of 3D cube volume
00242           nraysloc - local number of rays
00243           nnzloc - local number of voxels
00244           dm - projection matrix
00245           origin, ri - origin and radius of sphere
00246           ptrs, cord - pointers and coordinates for first ray
00247           myptrstart - determines starting rays
00248           x - 2D image
00249         Output:
00250           y - portion of the 3D volume
00251 */
00252 {
00253     int       i, j, iqx,iqy, xc, yc, zc, jj;
00254     float     xb, yb, dx, dy, dx1m, dy1m, dxdy;
00255     int       status = 0; 
00256 
00257     int xcent = origin[0];
00258     int ycent = origin[1];
00259     int zcent = origin[2];
00260 
00261     int nx = volsize[0];
00262     int ny = volsize[1];
00263 
00264     // Phi: adding the shift parameters that get passed in as the last two entries of dm
00265     float sx, sy;
00266 
00267     sx = dm(7);
00268     sy = dm(8);
00269 
00270     if ( nx > 2*ri) {
00271         for (i = 1; i <= nraysloc; i++) {
00272             zc = cord(1,myptrstart+i) - zcent;
00273             yc = cord(2,myptrstart+i) - ycent;
00274             xc = cord(3,myptrstart+i) - xcent;
00275 
00276             xb = zc*dm(1)+yc*dm(2)+xc*dm(3) + xcent + sx;
00277             yb = zc*dm(4)+yc*dm(5)+xc*dm(6) + ycent + sy;
00278 
00279             for (j = ptrs(myptrstart+i); j <ptrs(myptrstart+i+1); j++) {
00280                 jj = j-ptrs[myptrstart]+1;
00281 
00282                 // jj is the index on the local volume
00283                 iqx = ifix(xb);
00284                 iqy = ifix(yb);
00285 
00286                 dx = xb - iqx;
00287                 dy = yb - iqy;
00288                 dx1m = 1.0 - dx;
00289                 dy1m = 1.0 - dy;
00290                 dxdy = dx*dy;
00291 /*
00292 c               y(j) = y(j) + dx1m*dy1m*x(iqx  , iqy)
00293 c     &                     + dx1m*dy  *x(iqx  , iqy+1)
00294 c     &                     + dx  *dy1m*x(iqx+1, iqy)
00295 c     &                     + dx  *dy  *x(iqx+1, iqy+1)  
00296 c
00297 c              --- faster version of the above commented out
00298 c                  code (derived by summing the following table 
00299 c                  of coefficients along  the colunms) ---
00300 c
00301 c                        1         dx        dy      dxdy
00302 c                     ------   --------  --------  -------
00303 c                      x(i,j)   -x(i,j)   -x(i,j)    x(i,j)  
00304 c                                        x(i,j+1) -x(i,j+1)
00305 c                              x(i+1,j)           -x(i+1,j)
00306 c                                                x(i+1,j+1) 
00307 c
00308 */
00309                 // Phi: add index checking, now that shifts are being used
00310                 if ( iqx <= nx && iqy <= ny && iqx >= 1 && iqy >= 1 ) {
00311                     y(jj) += x(iqx,iqy);
00312                     if ( iqx + 1 <= nx && iqx + 1 >= 1 ) {
00313                         y(jj) += dx*(-x(iqx,iqy)+x(iqx+1,iqy));
00314                     }
00315                     if ( iqy + 1 <= ny && iqy + 1 >= 1 ) {
00316                         y(jj) += dy*(-x(iqx,iqy)+x(iqx,iqy+1));
00317                     }
00318                     if ( iqx + 1 <= nx && iqy + 1 <= ny && iqx + 1 >= 1 && iqy + 1 >= 1 ) {
00319                         y(jj) += dxdy*( x(iqx,iqy) - x(iqx,iqy+1) -x(iqx+1,iqy) + x(iqx+1,iqy+1) );
00320                     }
00321                 }
00322 
00323 //                y(j) += x(iqx,iqy)
00324 //                     +  dx*(-x(iqx,iqy)+x(iqx+1,iqy))
00325 //                     +  dy*(-x(iqx,iqy)+x(iqx,iqy+1))
00326 //                     +  dxdy*( x(iqx,iqy) - x(iqx,iqy+1) 
00327 //                              -x(iqx+1,iqy) + x(iqx+1,iqy+1) );
00328 
00329                xb += dm(1);
00330                yb += dm(4);
00331             } // end for j
00332         } // end for i
00333      }
00334     else {
00335         fprintf(stderr, "bckpj3: nx must be greater than 2*ri\n");
00336     }
00337 
00338     return status;
00339 }
00340 
00341 #undef x
00342 #undef y
00343