00001 #include "mpi.h"
00002
00003 #include "emdata.h"
00004 #include "project3d.h"
00005
00006 #define cube(i,j,k) cube[ (((k)-1)*ny + (j)-1)*nx + (i)-1 ]
00007 #define sphere(i) sphere[(i)-1]
00008 #define cord(i,j) cord[((j)-1)*3 + (i) -1]
00009 #define ptrs(i) ptrs[(i)-1]
00010 #define dm(i) dm[(i)-1]
00011
00012 int cb2sph(float *cube, Vec3i volsize, int ri, Vec3i origin,
00013 int nnz0, int *ptrs, int *cord, float *sphere)
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
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
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 }
00061 if (jnz > 0) {
00062 ptrs(nrays+1) = ptrs(nrays) + jnz;
00063 }
00064 }
00065 }
00066 if (nnz != nnz0) status = -1;
00067 return status;
00068 }
00069
00070
00071 int sph2cb(float *sphere, Vec3i volsize, int nrays, int ri,
00072 int nnz0, int *ptrs, int *cord, float *cube)
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
00080
00081 r2 = ri*ri;
00082 nnz = 0;
00083 ptrs(1) = 1;
00084
00085
00086
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 }
00101
00102 #define x(i) x[(i)-1]
00103 #define y(i,j) y[((j)-1)*nx + (i) - 1]
00104
00105
00106 int fwdpj3(Vec3i volsize, int nrays, int nnz, float *dm,
00107 Vec3i origin, int ri, int *ptrs, int *cord,
00108 float *x, float *y)
00109 {
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
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
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
00162
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
00171 y(iqx ,iqy) += dipx1m*dipy1m;
00172 if (iqx + 1 <= nx && iqy <= ny && iqx >= 0 && iqy >= 1)
00173
00174 y(iqx+1,iqy) += dipx*dipy1m;
00175 if (iqx + 1 <= nx && iqy + 1 <= ny && iqx >= 0 && iqy >= 0)
00176
00177 y(iqx+1,iqy+1) += dipx*dipy;
00178 if (iqx <= nx && iqy + 1 <= ny && iqx >= 1 && iqy >= 0)
00179
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 }
00192 #undef x
00193 #undef y
00194
00195 #define y(i) y[(i)-1]
00196 #define x(i,j) x[((j)-1)*nx + (i) - 1]
00197
00198
00199 int bckpj3(Vec3i volsize, int nrays, int nnz, float *dm,
00200 Vec3i origin, int ri, int *ptrs, int *cord,
00201 float *x, float *y)
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
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
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
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
00272
00273
00274
00275
00276
00277 xb += dm(1);
00278 yb += dm(4);
00279 }
00280 }
00281 }
00282 else {
00283 fprintf(stderr, "bckpj3: nx must be greater than 2*ri\n");
00284 }
00285
00286 return status;
00287 }
00288
00289 #undef x
00290 #undef y
00291 #undef dm
00292
00293 int getnnz(Vec3i volsize, int ri, Vec3i origin, int *nrays, int *nnz)
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
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
00318
00319
00320
00321 int xcent = origin[0];
00322 int ycent = origin[1];
00323 int zcent = origin[2];
00324
00325
00326
00327
00328
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 }
00349 }
00350 return status;
00351 }
00352
00353
00354 int make_proj_mat(float phi, float theta, float psi, float * dm)
00355 {
00356
00357
00358
00359
00360
00361
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 }