#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:
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 cord | ( | i, | |||
j | ) | cord[((j)-1)*3 + (i) -1] |
Definition at line 11 of file project3d_Cart.cpp.
#define dm | ( | i | ) | dm[(i)-1] |
Definition at line 13 of file project3d_Cart.cpp.
#define ptrs | ( | i | ) | ptrs[(i)-1] |
Definition at line 12 of file project3d_Cart.cpp.
#define x | ( | i, | |||
j | ) | x[((j)-1)*nx + (i) - 1] |
Definition at line 235 of file project3d_Cart.cpp.
#define x | ( | i | ) | x[(i)-1] |
Definition at line 235 of file project3d_Cart.cpp.
#define y | ( | i | ) | y[(i)-1] |
Definition at line 234 of file project3d_Cart.cpp.
#define y | ( | i, | |||
j | ) | y[((j)-1)*nx + (i) - 1] |
Definition at line 234 of file project3d_Cart.cpp.
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.
Referenced by LBD_Cart(), main(), recons3d_CGLS_mpi_Cart(), and recons3d_sirt_mpi_Cart().
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.
Referenced by LBD_Cart(), main(), recons3d_CGLS_mpi_Cart(), and recons3d_sirt_mpi_Cart().
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
Definition at line 15 of file project3d_Cart.cpp.
Referenced by main(), recons3d_CGLS_mpi_Cart(), recons3d_HyBR_mpi_Cart(), and recons3d_sirt_mpi_Cart().
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.
Referenced by main(), recons3d_CGLS_mpi_Cart(), recons3d_HyBR_mpi_Cart(), and recons3d_sirt_mpi_Cart().
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 }