EMAN2
utilcomm2d.cpp
Go to the documentation of this file.
00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <mpi.h>
00004 #include <math.h>
00005 #include <string.h>
00006 
00007 int setpart_gc1(MPI_Comm comm_2d, int nangs, int *psize, int *nbase)
00008 /* This function provides the partition for the set of 2D images along the first column of processors
00009 
00010         Input:
00011             comm_2d - Cartesian communicator
00012               nangs - number of 2D images
00013         Output:
00014               psize - vector containing the partition distribution
00015               nbase - vector containing the base of the partitions
00016            nangsloc - number of local angles
00017 */
00018 {
00019         int ROW = 0, COL = 1;
00020         int dims[2], periods[2], mycoords[2];
00021         int nangsloc, nrem;
00022 
00023         // Get information associated with comm_2d
00024         MPI_Cart_get(comm_2d, 2, dims, periods, mycoords);
00025  
00026         nangsloc = nangs/dims[ROW];
00027         nrem = nangs - dims[ROW]*nangsloc;
00028         if (mycoords[ROW] < nrem) nangsloc++;
00029 
00030         for (int i = 0; i < dims[ROW]; i++) {
00031                 psize[i] = nangs/dims[ROW];
00032                 if (i < nrem) psize[i]++;
00033         }
00034  
00035         nbase[0] = 0; 
00036         for (int i = 1; i < dims[ROW]; i++) {
00037                 nbase[i] = nbase[i-1] + psize[i-1];
00038         }
00039    
00040         return nangsloc;
00041 }
00042 
00043 int getnnz(int *volsize, int ri, int *origin, int *nrays, int *nnz) 
00044 /*
00045    purpose: count the number of voxels within a sphere centered
00046             at origin and with a radius ri.
00047 
00048    input:
00049      volsize vector containing size of the volume/cube
00050      ri      radius of the object embedded in the cube.
00051      origin  coordinates for the center of the volume
00052    
00053      output:
00054      nnz    total number of voxels within the sphere (of radius ri)
00055      nrays  number of rays in z-direction. 
00056 */
00057 {
00058     int  ix, iy, iz, rs, r2, xs, ys, zs, xx, yy, zz;
00059     int  ftm=0, status = 0;
00060 
00061     r2    = ri*ri;
00062     *nnz  = 0;
00063     *nrays = 0;
00064     int nx = volsize[0];
00065     int ny = volsize[1];
00066     int nz = volsize[2];
00067     // int nx = (int)volsize[0];
00068     // int ny = (int)volsize[1];
00069     // int nz = (int)volsize[2];
00070 
00071     int xcent = origin[0]; 
00072     int ycent = origin[1]; 
00073     int zcent = origin[2]; 
00074     // int xcent = (int)origin[0]; 
00075     // int ycent = (int)origin[1]; 
00076     // int zcent = (int)origin[2]; 
00077 
00078     // need to add some error checking 
00079     for (ix = 1; ix <=nx; ix++) { 
00080         xs  = ix-xcent;
00081         xx  = xs*xs;
00082         for (iy = 1; iy <= ny; iy++) {
00083             ys = iy-ycent;
00084             yy = ys*ys; 
00085             ftm = 1;
00086             for (iz = 1; iz <= nz; iz++) {
00087                 zs = iz-zcent;
00088                 zz = zs*zs;
00089                 rs = xx + yy + zz;
00090                 if (rs <= r2) {
00091                     (*nnz)++;
00092                     if (ftm) {
00093                        (*nrays)++;
00094                        ftm = 0;
00095                     }
00096                 }
00097             }
00098         } // end for iy
00099     } // end for ix
00100     return status;
00101 }
00102 
00103 
00104 int setpart_gr1(MPI_Comm comm_2d, int nnz, int *nnzpart, int *nnzbase)
00105 /* This function provides the ideal partition for nnz in the 3D volume along the first row of processors.
00106 
00107         Input:
00108             comm_2d - Cartesian communicator
00109                 nnz - total number of non-zeros in the volume
00110         Output:
00111             nnzpart - vector containing the partition distribution
00112             nnzbase - vector containing the base of the partitions
00113              nnzloc - initial local nnz
00114 */
00115 {
00116         int ROW = 0, COL = 1;
00117         int dims[2], periods[2], mycoords[2];
00118         int nnzloc, nrem;
00119         
00120         // Get information associated with comm_2d
00121         MPI_Cart_get(comm_2d, 2, dims, periods, mycoords);
00122  
00123         nnzloc = nnz/dims[COL];
00124         nrem = nnz - dims[COL]*nnzloc;
00125         if (mycoords[COL] < nrem) nnzloc++;
00126 
00127         for (int i = 0; i < dims[COL]; i++) {
00128                 nnzpart[i] = nnz/dims[COL];
00129                 if (i < nrem) nnzpart[i]++;
00130         }
00131  
00132         nnzbase[0] = 0; 
00133         for (int i = 1; i < dims[COL]; i++) {
00134                 nnzbase[i] = nnzbase[i-1] + nnzpart[i-1];
00135         }
00136         nnzbase[dims[COL]] = nnz;
00137         return nnzloc;
00138 }
00139 
00140 #define cord(i,j)   cord[((j)-1)*3 + (i) -1]
00141 #define ptrs(i)     ptrs[(i)-1]
00142 #define dm(i)       dm[(i)-1]
00143 
00144 int getcb2sph(int *volsize, int ri, int *origin, int nnz0, int *ptrs, int *cord) 
00145 /* 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.
00146 
00147         Input:
00148          volsize - dimensions of 3D cube
00149               ri - radius of sphere
00150           origin - origin of sphere
00151             nnz0 - number of non-zeros
00152 
00153         Output:
00154               ptrs - vector containing pointers to the beginning of each ray
00155               cord - vector containing the x,y,z coordinates of the beginning of each ray
00156 */
00157 
00158 {
00159     int    xs, ys, zs, xx, yy, zz, rs, r2;
00160     int    ix, iy, iz, jnz, nnz, nrays;
00161     int    ftm = 0, status = 0;  
00162 
00163     int xcent = (int)origin[0];
00164     int ycent = (int)origin[1];
00165     int zcent = (int)origin[2];
00166 
00167     int nx = (int)volsize[0];
00168     int ny = (int)volsize[1];
00169     int nz = (int)volsize[2];
00170 
00171     r2      = ri*ri;
00172     nnz     = 0;
00173     nrays    = 0;
00174     ptrs(1) = 1;
00175 
00176     for (ix = 1; ix <= nx; ix++) {
00177        xs  = ix-xcent;
00178        xx  = xs*xs;
00179        for ( iy = 1; iy <= ny; iy++ ) {
00180            ys = iy-ycent;
00181            yy = ys*ys;
00182            jnz = 0;
00183 
00184            ftm = 1;
00185            // not the most efficient implementation
00186            for (iz = 1; iz <= nz; iz++) {
00187                zs = iz-zcent;
00188                zz = zs*zs;
00189                rs = xx + yy + zz;
00190                if (rs <= r2) {
00191                   jnz++;
00192                   nnz++;
00193      //             sphere(nnz) = cube(iz, iy, ix); 
00194 
00195                   //  record the coordinates of the first nonzero ===
00196                   if (ftm) {
00197                      nrays++;
00198                      cord(1,nrays) = iz; 
00199                      cord(2,nrays) = iy; 
00200                      cord(3,nrays) = ix;
00201                      ftm = 0;
00202                   }
00203                }
00204             } // end for (iz..)
00205             if (jnz > 0) {
00206                 ptrs(nrays+1) = ptrs(nrays) + jnz;
00207             }  // endif (jnz)
00208        } // end for iy
00209     } // end for ix
00210     if (nnz != nnz0) status = -1;
00211     return status;
00212 }
00213 
00214 int sphpart(MPI_Comm comm_2d, int nrays, int *ptrs, int *nnzbase, int *ptrstart)
00215 /* This function divides the volume among the column processors by computing nraysloc for each processor and the starting pointer for the first ray.
00216 
00217         Input:
00218             comm_2d - Cartesian communicator
00219               nrays - total number of rays
00220                ptrs - vector containing pointers
00221             nnzbase - ideal volume partition of nnz
00222         Output:
00223            ptrstart - vector containing all the starting pointers for each column processor group
00224            nraysloc - actual number of local rays
00225 */
00226 {
00227         int ROW = 0, COL = 1;
00228         int dims[2], periods[2], mycoords[2];
00229         int nraysloc = 0;
00230         
00231         // Get information associated with comm_2d
00232         MPI_Cart_get(comm_2d, 2, dims, periods, mycoords);
00233 
00234         int count = 1;
00235         if (mycoords[COL] == 0){ // First column starts out with the first ray
00236           nraysloc++;
00237         }
00238         ptrstart[0] = 0;
00239 
00240         for(int i=1; i<nrays ; i++){
00241           if (ptrs[i]-1 <= nnzbase[count] && ptrs[i+1]-1 >= nnzbase[count]){ 
00242                 //nnzbase is between or equal to ptrs
00243 
00244             if (nnzbase[count]-(ptrs[i]-1)>= ptrs[i+1]-1-nnzbase[count]){
00245               if(mycoords[COL] == count-1){  // ray belongs to count-1
00246                 nraysloc++;
00247               }
00248               ptrstart[count] = i+1;
00249               count++;
00250 
00251             } else { //nnzbase[count]-(ptrs[i]-1)>= ptrs[i+1]-1-nnzbase[count]
00252                 if(mycoords[COL] == count){// ray belongs to count and it's a new starting ray
00253                   nraysloc++;
00254                 }
00255                 ptrstart[count] = i;
00256                 count++;
00257             }
00258 
00259         }
00260         else {  //ptrs[i]-1 > nnzbase[count] so ray belongs to count-1
00261           if(mycoords[COL] == count-1){
00262             nraysloc++;
00263           }
00264 
00265         }
00266         } // end for
00267         ptrstart[dims[COL]] = nrays;
00268         return nraysloc;
00269 
00270 }
00271 
00272 int make_proj_mat(float phi, float theta, float psi, float * dm)
00273 /* This function computes the projection matrix provided the Euler angles.
00274         Input:
00275           phi, theta, psi - Euler angles
00276         Output:
00277           dm - projection matrix
00278 */
00279 {
00280 //     float cphi=cos(phi);
00281 //     float sphi=sin(phi);
00282 //     float cthe=cos(theta);
00283 //     float sthe=sin(theta);
00284 //     float cpsi=cos(psi);
00285 //     float spsi=sin(psi);
00286 
00287     double cphi, sphi, cthe, sthe, cpsi, spsi;
00288     double dphi = phi;
00289     double dthe = theta;
00290     double dpsi = psi;
00291     sincos(dphi, &sphi, &cphi);
00292     sincos(dthe, &sthe, &cthe);
00293     sincos(dpsi, &spsi, &cpsi);
00294 
00295     dm[0]=cphi*cthe*cpsi-sphi*spsi;
00296     dm[1]=sphi*cthe*cpsi+cphi*spsi;
00297     dm[2]=-sthe*cpsi;
00298     dm[3]=-cphi*cthe*spsi-sphi*cpsi;
00299     dm[4]=-sphi*cthe*spsi+cphi*cpsi;
00300     dm[5]=sthe*spsi;
00301 
00302     return 0;
00303 }
00304 
00305 int ifix(float a)
00306 /* This function will floor non-negative numbers and ceil negative numbers.
00307 */
00308 {
00309     int ia = 0;
00310     if (a >= 0.0) {
00311        ia = (int)floor(a);
00312     }
00313     else {
00314        ia = (int)ceil(a);
00315     }
00316     return ia;
00317 }
00318 
00319 #define x(i)        x[(i)-1]
00320 #define y(i,j)      y[((j)-1)*nx + (i) - 1]
00321 
00322 int fwdpj3_Cart(int *volsize, int nraysloc, int nnzloc, float *dm, int *origin, int ri, int *ptrs, int *cord, int myptrstart, float *x, float *y)
00323 /* 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
00324 
00325         Input:
00326           volsize - size (nx,ny,nz) of 3D cube volume
00327           nraysloc - local number of rays within the compact spherical representation
00328           nnzloc - local number of voxels
00329           dm - an array of size 9 storing transformation 
00330           origin, ri - origin and radius of sphere
00331           ptrs, cord - pointers and coordinates for first ray
00332           myptrstart - determines starting rays
00333           x - portion of the 3D input volume
00334         Output:
00335           y - 2D output image
00336 */
00337 {
00338     int    iqx, iqy, i, j, xc, yc, zc, jj;
00339     float  ct, dipx, dipy, dipx1m, dipy1m, xb, yb, dm1, dm4;
00340     int    status = 0;
00341     
00342     // Phi: adding the shift parameters that get passed in as the last two entries of dm
00343     float sx, sy;
00344 
00345     sx = dm(7);
00346     sy = dm(8);
00347 
00348     int xcent = origin[0];
00349     int ycent = origin[1];
00350     int zcent = origin[2];
00351 
00352     int nx = volsize[0];
00353     int ny = volsize[1];
00354 
00355     dm1 = dm(1);
00356     dm4 = dm(4);
00357  
00358     if ( nx > 2*ri ) {
00359         for (i = 1; i <= nraysloc; i++) {
00360 
00361             zc = cord(1,myptrstart+i)-zcent;
00362             yc = cord(2,myptrstart+i)-ycent;
00363             xc = cord(3,myptrstart+i)-xcent;
00364             xb = zc* dm(1) +yc* dm(2) +xc* dm(3) + xcent + sx;
00365             yb = zc* dm(4) +yc* dm(5) +xc* dm(6) + ycent + sy;
00366 
00367             for (j = ptrs(myptrstart+i); j< ptrs(myptrstart+i+1); j++) {
00368                jj = j-ptrs[myptrstart]+1;
00369                iqx = ifix(xb);
00370                iqy = ifix(yb);
00371 
00372                 // jj is the index on the local volume
00373                ct   = x(jj);
00374 
00375                // dipx =  xb - (float)(iqx);
00376                // dipy = (yb - (float)(iqy)) * ct;
00377                    dipx =  xb - iqx;
00378                    dipy = (yb - iqy) * ct;
00379 
00380                dipy1m = ct - dipy;
00381                dipx1m = 1.0 - dipx;
00382 
00383                         if (iqx <= nx && iqy <= ny && iqx >= 1 && iqy >= 1) 
00384                // y(iqx  ,iqy)   = y(iqx  ,iqy)   + dipx1m*dipy1m;
00385                y(iqx  ,iqy)   +=  dipx1m*dipy1m;
00386                         if (iqx + 1 <= nx && iqy <= ny && iqx >= 0 && iqy >= 1) 
00387                // y(iqx+1,iqy)   = y(iqx+1,iqy)   + dipx*dipy1m; 
00388                y(iqx+1,iqy)   +=  dipx*dipy1m; 
00389                         if (iqx + 1 <= nx && iqy + 1 <= ny && iqx >= 0 && iqy >= 0) 
00390                // y(iqx+1,iqy+1) = y(iqx+1,iqy+1) + dipx*dipy;         
00391                y(iqx+1,iqy+1) +=  dipx*dipy;         
00392                         if (iqx <= nx && iqy + 1 <= ny && iqx >= 1 && iqy >= 0) 
00393                // y(iqx  ,iqy+1) = y(iqx  ,iqy+1) + dipx1m*dipy;
00394                y(iqx  ,iqy+1) +=  dipx1m*dipy;
00395                xb += dm1;
00396                yb += dm4;
00397            }
00398         }
00399     }
00400     else {
00401         fprintf(stderr, " nx must be greater than 2*ri\n");
00402         exit(1);
00403     }
00404     return status;
00405 }
00406 #undef x
00407 #undef y
00408 
00409 #define y(i)        y[(i)-1]
00410 #define x(i,j)      x[((j)-1)*nx + (i) - 1]
00411 
00412 int bckpj3_Cart(int *volsize, int nraysloc, int nnzloc, float *dm, int *origin, int ri, int *ptrs, int *cord, int myptrstart, float *x, float *y)
00413 /* 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
00414 
00415         Input:
00416           volsize - size of 3D cube volume
00417           nraysloc - local number of rays
00418           nnzloc - local number of voxels
00419           dm - projection matrix
00420           origin, ri - origin and radius of sphere
00421           ptrs, cord - pointers and coordinates for first ray
00422           myptrstart - determines starting rays
00423           x - 2D image
00424         Output:
00425           y - portion of the 3D volume
00426 */
00427 {
00428     int       i, j, iqx,iqy, xc, yc, zc, jj;
00429     float     xb, yb, dx, dy, dx1m, dy1m, dxdy;
00430     int       status = 0; 
00431 
00432     int xcent = origin[0];
00433     int ycent = origin[1];
00434     int zcent = origin[2];
00435 
00436     int nx = volsize[0];
00437     int ny = volsize[1];
00438 
00439     // Phi: adding the shift parameters that get passed in as the last two entries of dm
00440     float sx, sy;
00441 
00442     sx = dm(7);
00443     sy = dm(8);
00444 
00445     if ( nx > 2*ri) {
00446         for (i = 1; i <= nraysloc; i++) {
00447             zc = cord(1,myptrstart+i) - zcent;
00448             yc = cord(2,myptrstart+i) - ycent;
00449             xc = cord(3,myptrstart+i) - xcent;
00450 
00451             xb = zc*dm(1)+yc*dm(2)+xc*dm(3) + xcent + sx;
00452             yb = zc*dm(4)+yc*dm(5)+xc*dm(6) + ycent + sy;
00453 
00454             for (j = ptrs(myptrstart+i); j <ptrs(myptrstart+i+1); j++) {
00455                 jj = j-ptrs[myptrstart]+1;
00456 
00457                 // jj is the index on the local volume
00458                 iqx = ifix(xb);
00459                 iqy = ifix(yb);
00460 
00461                 dx = xb - iqx;
00462                 dy = yb - iqy;
00463                 dx1m = 1.0 - dx;
00464                 dy1m = 1.0 - dy;
00465                 dxdy = dx*dy;
00466 /*
00467 c               y(j) = y(j) + dx1m*dy1m*x(iqx  , iqy)
00468 c     &                     + dx1m*dy  *x(iqx  , iqy+1)
00469 c     &                     + dx  *dy1m*x(iqx+1, iqy)
00470 c     &                     + dx  *dy  *x(iqx+1, iqy+1)  
00471 c
00472 c              --- faster version of the above commented out
00473 c                  code (derived by summing the following table 
00474 c                  of coefficients along  the colunms) ---
00475 c
00476 c                        1         dx        dy      dxdy
00477 c                     ------   --------  --------  -------
00478 c                      x(i,j)   -x(i,j)   -x(i,j)    x(i,j)  
00479 c                                        x(i,j+1) -x(i,j+1)
00480 c                              x(i+1,j)           -x(i+1,j)
00481 c                                                x(i+1,j+1) 
00482 c
00483 */
00484                 // Phi: add index checking, now that shifts are being used
00485                 if ( iqx <= nx && iqy <= ny && iqx >= 1 && iqy >= 1 ) {
00486                     y(jj) += x(iqx,iqy);
00487                     if ( iqx + 1 <= nx && iqx + 1 >= 1 ) {
00488                         y(jj) += dx*(-x(iqx,iqy)+x(iqx+1,iqy));
00489                     }
00490                     if ( iqy + 1 <= ny && iqy + 1 >= 1 ) {
00491                         y(jj) += dy*(-x(iqx,iqy)+x(iqx,iqy+1));
00492                     }
00493                     if ( iqx + 1 <= nx && iqy + 1 <= ny && iqx + 1 >= 1 && iqy + 1 >= 1 ) {
00494                         y(jj) += dxdy*( x(iqx,iqy) - x(iqx,iqy+1) -x(iqx+1,iqy) + x(iqx+1,iqy+1) );
00495                     }
00496                 }
00497 
00498 //                y(j) += x(iqx,iqy)
00499 //                     +  dx*(-x(iqx,iqy)+x(iqx+1,iqy))
00500 //                     +  dy*(-x(iqx,iqy)+x(iqx,iqy+1))
00501 //                     +  dxdy*( x(iqx,iqy) - x(iqx,iqy+1) 
00502 //                              -x(iqx+1,iqy) + x(iqx+1,iqy+1) );
00503 
00504                xb += dm(1);
00505                yb += dm(4);
00506             } // end for j
00507         } // end for i
00508      }
00509     else {
00510         fprintf(stderr, "bckpj3: nx must be greater than 2*ri\n");
00511     }
00512 
00513     return status;
00514 }
00515 
00516 #undef x
00517 #undef y
00518 
00519 
00520 
00521 
00522 #define x(i)        x[(i)-1]
00523 #define y(i,j)      y[((j)-1)*nx + (i) - 1]
00524 
00525 // project from 3D to 2D (single image)
00526 int fwdpj3(int *volsize, int nrays, int   nnz, float *dm, int *origin, int ri, int *ptrs, int *cord, float *x, float  *y)
00527 {
00528     /*
00529         purpose:  y <--- proj(x)
00530         input  :  volsize  the size (nx,ny,nz) of the volume
00531                   nrays    number of rays within the compact spherical 
00532                            representation
00533                   nnz      number of voxels within the sphere
00534                   dm       an array of size 9 storing transformation 
00535                            associated with the projection direction
00536                   origin   coordinates of the center of the volume
00537                   ri       radius of the sphere
00538                   ptrs     the beginning address of each ray
00539                   cord     the coordinates of the first point in each ray
00540                   x        3d input volume
00541                   y        2d output image 
00542     */
00543 
00544     int    iqx, iqy, i, j, xc, yc, zc;
00545     float  ct, dipx, dipy, dipx1m, dipy1m, xb, yb, dm1, dm4;
00546     int    status = 0;
00547     
00548     // Phi: adding the shift parameters that get passed in as the last two entries of dm
00549     float sx, sy;
00550 
00551     sx = dm(7);
00552     sy = dm(8);
00553 
00554     int xcent = origin[0];
00555     int ycent = origin[1];
00556     int zcent = origin[2];
00557 
00558     int nx = volsize[0];
00559     int ny = volsize[1];
00560 
00561     dm1 = dm(1);
00562     dm4 = dm(4);
00563  
00564     if ( nx > 2*ri ) {
00565         for (i = 1; i <= nrays; i++) {
00566 
00567             zc = cord(1,i)-zcent;
00568             yc = cord(2,i)-ycent;
00569             xc = cord(3,i)-xcent;
00570             xb = zc* dm(1) +yc* dm(2) +xc* dm(3) + xcent + sx;
00571             yb = zc* dm(4) +yc* dm(5) +xc* dm(6) + ycent + sy;
00572 
00573             for (j = ptrs(i); j< ptrs(i+1); j++) {
00574                iqx = ifix(xb);
00575                iqy = ifix(yb);
00576 
00577                ct   = x(j);
00578 
00579                // dipx =  xb - (float)(iqx);
00580                // dipy = (yb - (float)(iqy)) * ct;
00581                    dipx =  xb - iqx;
00582                    dipy = (yb - iqy) * ct;
00583 
00584                dipy1m = ct - dipy;
00585                dipx1m = 1.0 - dipx;
00586 
00587                         if (iqx <= nx && iqy <= ny && iqx >= 1 && iqy >= 1) 
00588                // y(iqx  ,iqy)   = y(iqx  ,iqy)   + dipx1m*dipy1m;
00589                y(iqx  ,iqy)   +=  dipx1m*dipy1m;
00590                         if (iqx + 1 <= nx && iqy <= ny && iqx >= 0 && iqy >= 1) 
00591                // y(iqx+1,iqy)   = y(iqx+1,iqy)   + dipx*dipy1m; 
00592                y(iqx+1,iqy)   +=  dipx*dipy1m; 
00593                         if (iqx + 1 <= nx && iqy + 1 <= ny && iqx >= 0 && iqy >= 0) 
00594                // y(iqx+1,iqy+1) = y(iqx+1,iqy+1) + dipx*dipy;         
00595                y(iqx+1,iqy+1) +=  dipx*dipy;         
00596                         if (iqx <= nx && iqy + 1 <= ny && iqx >= 1 && iqy >= 0) 
00597                // y(iqx  ,iqy+1) = y(iqx  ,iqy+1) + dipx1m*dipy;
00598                y(iqx  ,iqy+1) +=  dipx1m*dipy;
00599                xb += dm1;
00600                yb += dm4;
00601            }
00602         }
00603     }
00604     else {
00605         fprintf(stderr, " nx must be greater than 2*ri\n");
00606         exit(1);
00607     }
00608     return status;
00609 }
00610 #undef x
00611 #undef y
00612 
00613 
00614 #define y(i)        y[(i)-1]
00615 #define x(i,j)      x[((j)-1)*nx + (i) - 1]
00616 
00617 // backproject from 2D to 3D for a single image
00618 int bckpj3(int *volsize, int nrays, int   nnz, float *dm, int *origin, int ri, int *ptrs, int *cord, float *x, float *y)
00619 {
00620     int       i, j, iqx,iqy, xc, yc, zc;
00621     float     xb, yb, dx, dy, dx1m, dy1m, dxdy;
00622     int       status = 0; 
00623 
00624     int xcent = origin[0];
00625     int ycent = origin[1];
00626     int zcent = origin[2];
00627 
00628     int nx = volsize[0];
00629     int ny = volsize[1];
00630 
00631     // Phi: adding the shift parameters that get passed in as the last two entries of dm
00632     float sx, sy;
00633 
00634     sx = dm(7);
00635     sy = dm(8);
00636 
00637 
00638     if ( nx > 2*ri) {
00639         for (i = 1; i <= nrays; i++) {
00640             zc = cord(1,i) - zcent;
00641             yc = cord(2,i) - ycent;
00642             xc = cord(3,i) - xcent;
00643 
00644             xb = zc*dm(1)+yc*dm(2)+xc*dm(3) + xcent + sx;
00645             yb = zc*dm(4)+yc*dm(5)+xc*dm(6) + ycent + sy;
00646 
00647             for (j = ptrs(i); j <ptrs(i+1); j++) {
00648                 iqx = ifix(xb);
00649                 iqy = ifix(yb);
00650 
00651                 dx = xb - iqx;
00652                 dy = yb - iqy;
00653                 dx1m = 1.0 - dx;
00654                 dy1m = 1.0 - dy;
00655                 dxdy = dx*dy;
00656 /*
00657 c               y(j) = y(j) + dx1m*dy1m*x(iqx  , iqy)
00658 c     &                     + dx1m*dy  *x(iqx  , iqy+1)
00659 c     &                     + dx  *dy1m*x(iqx+1, iqy)
00660 c     &                     + dx  *dy  *x(iqx+1, iqy+1)  
00661 c
00662 c              --- faster version of the above commented out
00663 c                  code (derived by summing the following table 
00664 c                  of coefficients along  the colunms) ---
00665 c
00666 c                        1         dx        dy      dxdy
00667 c                     ------   --------  --------  -------
00668 c                      x(i,j)   -x(i,j)   -x(i,j)    x(i,j)  
00669 c                                        x(i,j+1) -x(i,j+1)
00670 c                              x(i+1,j)           -x(i+1,j)
00671 c                                                x(i+1,j+1) 
00672 c
00673 */
00674                 // Phi: add index checking, now that shifts are being used
00675                 if ( iqx <= nx && iqy <= ny && iqx >= 1 && iqy >= 1 ) {
00676                     y(j) += x(iqx,iqy);
00677                     if ( iqx + 1 <= nx && iqx + 1 >= 1 ) {
00678                         y(j) += dx*(-x(iqx,iqy)+x(iqx+1,iqy));
00679                     }
00680                     if ( iqy + 1 <= ny && iqy + 1 >= 1 ) {
00681                         y(j) += dy*(-x(iqx,iqy)+x(iqx,iqy+1));
00682                     }
00683                     if ( iqx + 1 <= nx && iqy + 1 <= ny && iqx + 1 >= 1 && iqy + 1 >= 1 ) {
00684                         y(j) += dxdy*( x(iqx,iqy) - x(iqx,iqy+1) -x(iqx+1,iqy) + x(iqx+1,iqy+1) );
00685                     }
00686                 }
00687 
00688 //                y(j) += x(iqx,iqy)
00689 //                     +  dx*(-x(iqx,iqy)+x(iqx+1,iqy))
00690 //                     +  dy*(-x(iqx,iqy)+x(iqx,iqy+1))
00691 //                     +  dxdy*( x(iqx,iqy) - x(iqx,iqy+1) 
00692 //                              -x(iqx+1,iqy) + x(iqx+1,iqy+1) );
00693 
00694                xb += dm(1);
00695                yb += dm(4);
00696             } // end for j
00697         } // end for i
00698      }
00699     else {
00700         fprintf(stderr, "bckpj3: nx must be greater than 2*ri\n");
00701     }
00702 
00703     return status;
00704 }
00705 
00706 #undef x
00707 #undef y
00708 
00709 #undef dm
00710