#include "mpi.h"
#include "emdata.h"
#include "project3d.h"
Include dependency graph for project3d.cpp:
Go to the source code of this file.
Defines | |
#define | cube(i, j, k) cube[ (((k)-1)*ny + (j)-1)*nx + (i)-1 ] |
#define | sphere(i) sphere[(i)-1] |
#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 | cb2sph (float *cube, Vec3i volsize, int ri, Vec3i origin, int nnz0, int *ptrs, int *cord, float *sphere) |
int | sph2cb (float *sphere, Vec3i volsize, int nrays, int ri, int nnz0, int *ptrs, int *cord, float *cube) |
int | fwdpj3 (Vec3i volsize, int nrays, int nnz, float *dm, Vec3i origin, int ri, int *ptrs, int *cord, float *x, float *y) |
int | bckpj3 (Vec3i volsize, int nrays, int nnz, float *dm, Vec3i origin, int ri, int *ptrs, int *cord, float *x, float *y) |
int | getnnz (Vec3i volsize, int ri, Vec3i origin, int *nrays, int *nnz) |
int | make_proj_mat (float phi, float theta, float psi, float *dm) |
#define cord | ( | i, | |||
j | ) | cord[((j)-1)*3 + (i) -1] |
Definition at line 8 of file project3d.cpp.
Definition at line 6 of file project3d.cpp.
#define dm | ( | i | ) | dm[(i)-1] |
Definition at line 10 of file project3d.cpp.
#define ptrs | ( | i | ) | ptrs[(i)-1] |
Definition at line 9 of file project3d.cpp.
#define sphere | ( | i | ) | sphere[(i)-1] |
Definition at line 7 of file project3d.cpp.
#define x | ( | i, | |||
j | ) | x[((j)-1)*nx + (i) - 1] |
Definition at line 196 of file project3d.cpp.
#define x | ( | i | ) | x[(i)-1] |
Definition at line 196 of file project3d.cpp.
#define y | ( | i | ) | y[(i)-1] |
Definition at line 195 of file project3d.cpp.
#define y | ( | i, | |||
j | ) | y[((j)-1)*nx + (i) - 1] |
Definition at line 195 of file project3d.cpp.
int bckpj3 | ( | Vec3i | volsize, | |
int | nrays, | |||
int | nnz, | |||
float * | dm, | |||
Vec3i | origin, | |||
int | ri, | |||
int * | ptrs, | |||
int * | cord, | |||
float * | x, | |||
float * | y | |||
) |
Definition at line 199 of file project3d.cpp.
Referenced by fgcalc(), main(), and recons3d_sirt_mpi().
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 }
int cb2sph | ( | float * | cube, | |
Vec3i | volsize, | |||
int | ri, | |||
Vec3i | origin, | |||
int | nnz0, | |||
int * | ptrs, | |||
int * | cord, | |||
float * | sphere | |||
) |
Definition at line 12 of file project3d.cpp.
References cord, cube, nnz, nrays, nx, ny, ptrs, sphere, and status.
Referenced by ali3d_d(), recons3d_sirt_mpi(), and unified().
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 }
int fwdpj3 | ( | Vec3i | volsize, | |
int | nrays, | |||
int | nnz, | |||
float * | dm, | |||
Vec3i | origin, | |||
int | ri, | |||
int * | ptrs, | |||
int * | cord, | |||
float * | x, | |||
float * | y | |||
) |
Definition at line 106 of file project3d.cpp.
Referenced by fcalc(), fgcalc(), main(), and recons3d_sirt_mpi().
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 }
Definition at line 293 of file project3d.cpp.
Referenced by ali3d_d(), main(), recons3d_CGLS_mpi_Cart(), recons3d_HyBR_mpi_Cart(), recons3d_sirt_mpi(), recons3d_sirt_mpi_Cart(), and unified().
00295 : 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 }
int make_proj_mat | ( | float | phi, | |
float | theta, | |||
float | psi, | |||
float * | dm | |||
) |
Definition at line 354 of file project3d.cpp.
Referenced by LBD_Cart(), main(), recons3d_CGLS_mpi_Cart(), recons3d_sirt_mpi(), and recons3d_sirt_mpi_Cart().
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 }
int sph2cb | ( | float * | sphere, | |
Vec3i | volsize, | |||
int | nrays, | |||
int | ri, | |||
int | nnz0, | |||
int * | ptrs, | |||
int * | cord, | |||
float * | cube | |||
) |
Definition at line 71 of file project3d.cpp.
References cord, cube, nnz, nx, ny, ptrs, sphere, and status.
Referenced by recons3d_CGLS_mpi_Cart(), recons3d_HyBR_mpi_Cart(), recons3d_sirt_mpi(), recons3d_sirt_mpi_Cart(), and unified().
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 }