Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

project3d_Cart.cpp File Reference

#include "mpi.h"
#include "stdlib.h"
#include "emdata.h"
#include "utilcomm_Cart.h"
#include "project3d_Cart.h"
#include "project3d.h"

Include dependency graph for project3d_Cart.cpp:

Include dependency graph

Go to the source code of this file.

Defines

#define cord(i, j)   cord[((j)-1)*3 + (i) -1]
#define ptrs(i)   ptrs[(i)-1]
#define dm(i)   dm[(i)-1]
#define x(i)   x[(i)-1]
#define y(i, j)   y[((j)-1)*nx + (i) - 1]
#define y(i)   y[(i)-1]
#define x(i, j)   x[((j)-1)*nx + (i) - 1]

Functions

int getcb2sph (Vec3i volsize, int ri, Vec3i origin, int nnz0, int *ptrs, int *cord)
int sphpart (MPI_Comm comm_2d, int nrays, int *ptrs, int *nnzbase, int *ptrstart)
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)
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)


Define Documentation

#define cord i,
 )     cord[((j)-1)*3 + (i) -1]
 

Definition at line 11 of file project3d_Cart.cpp.

#define dm  )     dm[(i)-1]
 

Definition at line 13 of file project3d_Cart.cpp.

#define ptrs  )     ptrs[(i)-1]
 

Definition at line 12 of file project3d_Cart.cpp.

#define x i,
 )     x[((j)-1)*nx + (i) - 1]
 

Definition at line 235 of file project3d_Cart.cpp.

#define x  )     x[(i)-1]
 

Definition at line 235 of file project3d_Cart.cpp.

#define y  )     y[(i)-1]
 

Definition at line 234 of file project3d_Cart.cpp.

#define y i,
 )     y[((j)-1)*nx + (i) - 1]
 

Definition at line 234 of file project3d_Cart.cpp.


Function Documentation

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
 

Definition at line 237 of file project3d_Cart.cpp.

00240              :
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 }

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
 

Definition at line 147 of file project3d_Cart.cpp.

00150              :
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 }
#undef x

int getcb2sph Vec3i  volsize,
int  ri,
Vec3i  origin,
int  nnz0,
int *  ptrs,
int *  cord
 

Definition at line 15 of file project3d_Cart.cpp.

00018              :
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 }

int sphpart MPI_Comm  comm_2d,
int  nrays,
int *  ptrs,
int *  nnzbase,
int *  ptrstart
 

Definition at line 85 of file project3d_Cart.cpp.

00088              :
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 }


Generated on Tue Jun 11 13:47:12 2013 for EMAN2 by  doxygen 1.3.9.1