00001 #include "mpi.h"
00002 #include "stdlib.h"
00003 #include "emdata.h"
00004
00005 #include "utilcomm_Cart.h"
00006 #include "project3d_Cart.h"
00007 #include "project3d.h"
00008
00009 using namespace EMAN;
00010
00011 #define cord(i,j) cord[((j)-1)*3 + (i) -1]
00012 #define ptrs(i) ptrs[(i)-1]
00013 #define dm(i) dm[(i)-1]
00014
00015 int getcb2sph(Vec3i volsize, int ri, Vec3i origin, int nnz0, int *ptrs, int *cord)
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
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
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
00065
00066
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 }
00076 if (jnz > 0) {
00077 ptrs(nrays+1) = ptrs(nrays) + jnz;
00078 }
00079 }
00080 }
00081 if (nnz != nnz0) status = -1;
00082 return status;
00083 }
00084
00085 int sphpart(MPI_Comm comm_2d, int nrays, int *ptrs, int *nnzbase, int *ptrstart)
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097 {
00098 int ROW = 0, COL = 1;
00099 int dims[2], periods[2], mycoords[2];
00100 int nraysloc = 0;
00101
00102
00103 MPI_Cart_get(comm_2d, 2, dims, periods, mycoords);
00104
00105 int count = 1;
00106 if (mycoords[COL] == 0){
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
00114
00115 if (nnzbase[count]-(ptrs[i]-1)>= ptrs[i+1]-1-nnzbase[count]){
00116 if(mycoords[COL] == count-1){
00117 nraysloc++;
00118 }
00119 ptrstart[count] = i+1;
00120 count++;
00121
00122 } else {
00123 if(mycoords[COL] == count){
00124 nraysloc++;
00125 }
00126 ptrstart[count] = i;
00127 count++;
00128 }
00129
00130 }
00131 else {
00132 if(mycoords[COL] == count-1){
00133 nraysloc++;
00134 }
00135
00136 }
00137 }
00138 ptrstart[dims[COL]] = nrays;
00139 return nraysloc;
00140
00141 }
00142
00143
00144 #define x(i) x[(i)-1]
00145 #define y(i,j) y[((j)-1)*nx + (i) - 1]
00146
00147 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)
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
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
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
00198 ct = x(jj);
00199
00200
00201
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
00210 y(iqx ,iqy) += dipx1m*dipy1m;
00211 if (iqx + 1 <= nx && iqy <= ny && iqx >= 0 && iqy >= 1)
00212
00213 y(iqx+1,iqy) += dipx*dipy1m;
00214 if (iqx + 1 <= nx && iqy + 1 <= ny && iqx >= 0 && iqy >= 0)
00215
00216 y(iqx+1,iqy+1) += dipx*dipy;
00217 if (iqx <= nx && iqy + 1 <= ny && iqx >= 1 && iqy >= 0)
00218
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 }
00231 #undef x
00232 #undef y
00233
00234 #define y(i) y[(i)-1]
00235 #define x(i,j) x[((j)-1)*nx + (i) - 1]
00236
00237 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)
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
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
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
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
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
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
00324
00325
00326
00327
00328
00329 xb += dm(1);
00330 yb += dm(4);
00331 }
00332 }
00333 }
00334 else {
00335 fprintf(stderr, "bckpj3: nx must be greater than 2*ri\n");
00336 }
00337
00338 return status;
00339 }
00340
00341 #undef x
00342 #undef y
00343