00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include "projector.h"
00037 #include "emdata.h"
00038 #include "interp.h"
00039 #include "emutil.h"
00040 #include "plugins/projector_template.h"
00041
00042 #ifdef WIN32
00043 #define M_PI 3.14159265358979323846f
00044 #endif //WIN32
00045
00046 #ifdef EMAN2_USING_CUDA
00047 #include "cuda/cuda_util.h"
00048 #include "cuda/cuda_projector.h"
00049 #endif
00050
00051 using namespace std;
00052 using namespace EMAN;
00053
00054 const string GaussFFTProjector::NAME = "gauss_fft";
00055 const string FourierGriddingProjector::NAME = "fourier_gridding";
00056 const string PawelProjector::NAME = "pawel";
00057 const string StandardProjector::NAME = "standard";
00058 const string ChaoProjector::NAME = "chao";
00059
00060 template <> Factory < Projector >::Factory()
00061 {
00062 force_add<GaussFFTProjector>();
00063 force_add<PawelProjector>();
00064 force_add<StandardProjector>();
00065 force_add<FourierGriddingProjector>();
00066 force_add<ChaoProjector>();
00067
00068
00069 }
00070
00071 EMData *GaussFFTProjector::project3d(EMData * image) const
00072 {
00073 Transform* t3d = params["transform"];
00074 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
00075 if ( image->get_ndim() != 3 ) throw ImageDimensionException("Error, the projection volume must be 3D");
00076
00077 EMData *f = image;
00078 if (!image->is_complex()) {
00079 image->process_inplace("xform.phaseorigin.tocorner");
00080 f = image->do_fft();
00081 f->process_inplace("xform.fourierorigin.tocenter");
00082 image->process_inplace("xform.phaseorigin.tocenter");
00083 }
00084
00085 int f_nx = f->get_xsize();
00086 int f_ny = f->get_ysize();
00087 int f_nz = f->get_zsize();
00088
00089 if (!f->is_complex() || f_nz != f_ny || f_nx != f_ny + 2) {
00090 LOGERR("Cannot project this image");
00091 return 0;
00092 }
00093
00094 f->ap2ri();
00095
00096 EMData *tmp = new EMData();
00097 tmp->set_size(f_nx, f_ny, 1);
00098 tmp->set_complex(true);
00099 tmp->set_ri(true);
00100
00101 float *data = tmp->get_data();
00102
00103 Transform r = t3d->get_rotation_transform();
00104 r.invert();
00105 float scale = t3d->get_scale();
00106
00107 int mode = params["mode"];
00108 float gauss_width = 1;
00109 if ( mode == 0 ) mode = 2;
00110 if (mode == 2 ) {
00111 gauss_width = EMConsts::I2G;
00112 }
00113 else if (mode == 3) {
00114 gauss_width = EMConsts::I3G;
00115 }
00116 else if (mode == 4) {
00117 gauss_width = EMConsts::I4G;
00118 }
00119 else if (mode == 6 || mode == 7) {
00120 gauss_width = EMConsts::I5G;
00121 }
00122
00123 for (int y = 0; y < f_ny; y++) {
00124 for (int x = 0; x < f_nx / 2; x++) {
00125 int ii = x * 2 + y * f_nx;
00126 #ifdef _WIN32
00127 if (_hypot(x, y - f_ny / 2) >= f_ny / 2 - 1) {
00128 #else
00129 if (hypot(x, y - f_ny / 2) >= f_ny / 2 - 1) {
00130 #endif
00131 data[ii] = 0;
00132 data[ii + 1] = 0;
00133 }
00134 else {
00135 Vec3f coord(x,(y - f_ny / 2),0);
00136 Vec3f soln = r*coord;
00137 float xx = soln[0];
00138 float yy = soln[1];
00139 float zz = soln[2];
00140
00141 int cc = 1;
00142 if (xx < 0) {
00143 xx = -xx;
00144 yy = -yy;
00145 zz = -zz;
00146 cc = -1;
00147 }
00148
00149 if (scale != 1.0) {
00150 xx *= scale;
00151 yy *= scale;
00152 zz *= scale;
00153 }
00154
00155 yy += f_ny / 2;
00156 zz += f_nz / 2;
00157
00158 if (yy < 0 || xx < 0 || zz < 0) {
00159 data[ii] = 0;
00160 data[ii+1] = 0;
00161 continue;
00162 }
00163
00164 if (interp_ft_3d(mode, f, xx, yy, zz, data + ii, gauss_width)) {
00165 data[ii + 1] *= cc;
00166 } else {
00167 data[ii] = 0;
00168 data[ii+1] = 0;
00169 }
00170 }
00171 }
00172 }
00173
00174 f->update();
00175 tmp->update();
00176
00177 tmp->process_inplace("xform.fourierorigin.tocorner");
00178 EMData *ret = tmp->do_ift();
00179 ret->process_inplace("xform.phaseorigin.tocenter");
00180
00181 ret->translate(t3d->get_trans());
00182
00183 if (t3d->get_mirror() ) ret->process_inplace("xform.flip",Dict("axis","x"));
00184
00185 Dict filter_d;
00186 filter_d["gauss_width"] = gauss_width;
00187 filter_d["ring_width"] = ret->get_xsize() / 2;
00188 ret->process_inplace("math.gausskernelfix", filter_d);
00189
00190 if( tmp )
00191 {
00192 delete tmp;
00193 tmp = 0;
00194 }
00195
00196 ret->set_attr("xform.projection",t3d);
00197 ret->update();
00198
00199 if(t3d) {delete t3d; t3d=0;}
00200
00201 return ret;
00202 }
00203
00204
00205
00206 bool GaussFFTProjector::interp_ft_3d(int mode, EMData * image, float x, float y,
00207 float z, float *data, float gw) const
00208 {
00209 float *rdata = image->get_data();
00210 int nx = image->get_xsize();
00211 int ny = image->get_ysize();
00212 int nz = image->get_zsize();
00213
00214 if ( mode == 0 ) mode = 2;
00215
00216 if (mode == 1) {
00217 int x0 = 2 * (int) floor(x + 0.5f);
00218 int y0 = (int) floor(y + 0.5f);
00219 int z0 = (int) floor(z + 0.5f);
00220
00221 data[0] = rdata[x0 + y0 * nx + z0 * nx * ny];
00222 data[1] = rdata[x0 + y0 * nx + z0 * nx * ny + 1];
00223 return true;
00224 }
00225 else if (mode == 2) {
00226 int x0 = (int) floor(x);
00227 int y0 = (int) floor(y);
00228 int z0 = (int) floor(z);
00229
00230 float dx = x - x0;
00231 float dy = y - y0;
00232 float dz = z - z0;
00233
00234 if (x0 <= nx/2- 2 && y0 <= ny - 1 && z0 <= nz - 1) {
00235 int i = (int) (x0 * 2 + y0 * nx + z0 * nx * ny);
00236
00237 float n = Util::agauss(1, dx, dy, dz, gw) +
00238 Util::agauss(1, 1 - dx, dy, dz, gw) +
00239 Util::agauss(1, dx, 1 - dy, dz, gw) +
00240 Util::agauss(1, 1 - dx, 1 - dy, dz, gw) +
00241 Util::agauss(1, dx, dy, 1 - dz, gw) +
00242 Util::agauss(1, 1 - dx, dy, 1 - dz, gw) +
00243 Util::agauss(1, dx, 1 - dy, 1 - dz, gw) + Util::agauss(1, 1 - dx, 1 - dy, 1 - dz, gw);
00244
00245 data[0] = Util::agauss(rdata[i], dx, dy, dz, gw) +
00246 Util::agauss(rdata[i + 2], 1 - dx, dy, dz, gw) +
00247 Util::agauss(rdata[i + nx], dx, 1 - dy, dz, gw) +
00248 Util::agauss(rdata[i + nx + 2], 1 - dx, 1 - dy, dz, gw) +
00249 Util::agauss(rdata[i + nx * ny], dx, dy, 1 - dz, gw) +
00250 Util::agauss(rdata[i + 2 + nx * ny], 1 - dx, dy, 1 - dz, gw) +
00251 Util::agauss(rdata[i + nx + nx * ny], dx, 1 - dy, 1 - dz, gw) +
00252 Util::agauss(rdata[i + 2 + nx + nx * ny], 1 - dx, 1 - dy, 1 - dz, gw) / n;
00253
00254 i++;
00255
00256 data[1] = Util::agauss(rdata[i], dx, dy, dz, gw) +
00257 Util::agauss(rdata[i + 2], 1 - dx, dy, dz, gw) +
00258 Util::agauss(rdata[i + nx], dx, 1 - dy, dz, gw) +
00259 Util::agauss(rdata[i + nx + 2], 1 - dx, 1 - dy, dz, gw) +
00260 Util::agauss(rdata[i + nx * ny], dx, dy, 1 - dz, gw) +
00261 Util::agauss(rdata[i + 2 + nx * ny], 1 - dx, dy, 1 - dz, gw) +
00262 Util::agauss(rdata[i + nx + nx * ny], dx, 1 - dy, 1 - dz, gw) +
00263 Util::agauss(rdata[i + 2 + nx + nx * ny], 1 - dx, 1 - dy, 1 - dz, gw) / n;
00264 return true;
00265 }
00266 return false;
00267 }
00268 else if (mode == 3) {
00269 int x0 = 2 * (int) floor(x + .5);
00270 int y0 = (int) floor(y + .5);
00271 int z0 = (int) floor(z + .5);
00272
00273 float *supp = image->setup4slice();
00274
00275 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
00276 float n = 0;
00277
00278 if (x0 == 0) {
00279 x0 += 4;
00280 size_t idx;
00281 float r, g;
00282 for (int k = z0 - 1; k <= z0 + 1; k++) {
00283 for (int j = y0 - 1; j <= y0 + 1; j++) {
00284 for (int i = x0 - 2; i <= x0 + 2; i += 2) {
00285 r = Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z);
00286 g = exp(-r / gw);
00287 n += g;
00288 idx = i + j * 12 + (size_t)k * 12 * ny;
00289 data[0] += supp[idx] * g;
00290 data[1] += supp[idx + 1] * g;
00291 }
00292 }
00293 }
00294 }
00295 else {
00296 size_t idx;
00297 float r, g;
00298 for (int k = z0 - 1; k <= z0 + 1; k++) {
00299 for (int j = y0 - 1; j <= y0 + 1; j++) {
00300 for (int i = x0 - 2; i <= x0 + 2; i += 2) {
00301 r = Util::hypot3sq(i / 2.0f - x, j - y, k - z);
00302 g = exp(-r / gw);
00303 n += g;
00304 idx = i + j * nx + (size_t)k * nx * ny;
00305 data[0] += rdata[idx] * g;
00306 data[1] += rdata[idx + 1] * g;
00307 }
00308 }
00309 }
00310 }
00311
00312 data[0] /= n;
00313 data[1] /= n;
00314 return true;
00315 }
00316 return false;
00317 }
00318 else if (mode == 4) {
00319 int x0 = 2 * (int) floor(x);
00320 int y0 = (int) floor(y);
00321 int z0 = (int) floor(z);
00322
00323 float *supp = image->setup4slice();
00324
00325 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
00326 float n = 0;
00327
00328 if (x0 == 0) {
00329 x0 += 4;
00330 size_t idx;
00331 float r, g;
00332 for (int k = z0 - 1; k <= z0 + 2; k++) {
00333 for (int j = y0 - 1; j <= y0 + 2; j++) {
00334 for (int i = x0 - 2; i <= x0 + 4; i += 2) {
00335 r = Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z);
00336 g = exp(-r / gw);
00337 n += g;
00338 idx = i + j * 12 + (size_t)k * 12 * ny;
00339 data[0] += supp[idx] * g;
00340 data[1] += supp[idx + 1] * g;
00341 }
00342 }
00343 }
00344 }
00345 else {
00346 float r, g;
00347 size_t idx;
00348 for (int k = z0 - 1; k <= z0 + 2; k++) {
00349 for (int j = y0 - 1; j <= y0 + 2; j++) {
00350 for (int i = x0 - 2; i <= x0 + 4; i += 2) {
00351 r = Util::hypot3sq(i / 2.0f - x, j - y, k - z);
00352 g = exp(-r / gw);
00353 n += g;
00354 idx = i + j * nx + (size_t)k * nx * ny;
00355 data[0] += rdata[idx] * g;
00356 data[1] += rdata[idx + 1] * g;
00357 }
00358 }
00359 }
00360 }
00361 data[0] /= n;
00362 data[1] /= n;
00363 return true;
00364 }
00365 return false;
00366 }
00367 else if (mode == 5) {
00368 int x0 = 2 * (int) floor(x + .5);
00369 int y0 = (int) floor(y + .5);
00370 int z0 = (int) floor(z + .5);
00371
00372 float *supp = image->setup4slice();
00373 float *gimx = Interp::get_gimx();
00374
00375 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
00376 int mx0 = -(int) floor((x - x0 / 2) * 39.0 + .5) - 78;
00377 int my0 = -(int) floor((y - y0) * 39.0 + .5) - 78;
00378 int mz0 = -(int) floor((z - z0) * 39.0 + .5) - 78;
00379
00380 float n = 0;
00381 int mmz = mz0;
00382 int mmy = my0;
00383 int mmx = mx0;
00384
00385 if (x0 < 4) {
00386 x0 += 4;
00387
00388 size_t idx;
00389 float g;
00390 for (int k = z0 - 2; k <= z0 + 2; k++, mmz += 39) {
00391 for (int j = y0 - 2; j <= y0 + 2; j++, mmy += 39) {
00392 for (int i = x0 - 4; i <= x0 + 4; i += 2, mmx += 39) {
00393 g = gimx[abs(mmx) + abs(mmy) * 100 + abs(mmz) * 10000];
00394 n += g;
00395 idx = i + j * 12 + (size_t)k * 12 * ny;
00396 data[0] += supp[idx] * g;
00397 data[1] += supp[idx + 1] * g;
00398 }
00399 }
00400 }
00401 }
00402 else {
00403 size_t ii;
00404 float g;
00405 for (int k = z0 - 2; k <= z0 + 2; k++, mmz += 39) {
00406 for (int j = y0 - 2; j <= y0 + 2; j++, mmy += 39) {
00407 for (int i = x0 - 4; i <= x0 + 4; i += 2, mmx += 39) {
00408 ii = i + j * nx + (size_t)k * nx * ny;
00409 g = gimx[abs(mmx) + abs(mmy) * 100 + abs(mmz) * 10000];
00410 n += g;
00411 data[0] += rdata[ii] * g;
00412 data[1] += rdata[ii + 1] * g;
00413 }
00414 }
00415 }
00416 }
00417
00418 data[0] /= n;
00419 data[1] /= n;
00420 return true;
00421 }
00422 return false;
00423 }
00424 else if (mode == 6) {
00425
00426 int x0 = 2 * (int) floor(x + .5);
00427 int y0 = (int) floor(y + .5);
00428 int z0 = (int) floor(z + .5);
00429
00430 float *supp = image->setup4slice();
00431
00432 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
00433 float n = 0;
00434
00435 if (x0 < 4) {
00436 x0 += 4;
00437 float r, g;
00438 size_t idx;
00439 for (int k = z0 - 2; k <= z0 + 2; k++) {
00440 for (int j = y0 - 2; j <= y0 + 2; j++) {
00441 for (int i = x0 - 4; i <= x0 + 4; i += 2) {
00442 r = Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z);
00443 g = exp(-r / gw);
00444 n += g;
00445 idx = i + j * 12 + (size_t)k * 12 * ny;
00446 data[0] += supp[idx] * g;
00447 data[1] += supp[idx + 1] * g;
00448 }
00449 }
00450 }
00451 }
00452 else {
00453 float r, g;
00454 size_t idx;
00455 for (int k = z0 - 2; k <= z0 + 2; k++) {
00456 for (int j = y0 - 2; j <= y0 + 2; j++) {
00457 for (int i = x0 - 4; i <= x0 + 4; i += 2) {
00458 r = Util::hypot3sq(i / 2.0f - x, j - y, k - z);
00459 g = exp(-r / gw);
00460 n += g;
00461 idx = i + j * nx + (size_t)k * nx * ny;
00462 data[0] += rdata[idx] * g;
00463 data[1] += rdata[idx + 1] * g;
00464 }
00465
00466 }
00467 }
00468 }
00469
00470 data[0] /= n;
00471 data[1] /= n;
00472
00473 return true;
00474 }
00475 return false;
00476 }
00477 else if (mode == 7) {
00478 int x0 = 2 * (int) floor(x + .5);
00479 int y0 = (int) floor(y + .5);
00480 int z0 = (int) floor(z + .5);
00481
00482 float *supp = image->setup4slice();
00483
00484 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
00485 float n = 0;
00486 if (x0 < 4) {
00487 x0 += 4;
00488 float r, g;
00489 size_t idx;
00490 for (int k = z0 - 2; k <= z0 + 2; k++) {
00491 for (int j = y0 - 2; j <= y0 + 2; j++) {
00492 for (int i = x0 - 4; i <= x0 + 4; i += 2) {
00493 r = sqrt(Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z));
00494 g = Interp::hyperg(r);
00495 n += g;
00496 idx = i + j * 12 + (size_t)k * 12 * ny;
00497 data[0] += supp[idx] * g;
00498 data[1] += supp[idx + 1] * g;
00499 }
00500 }
00501 }
00502 }
00503 else {
00504 float r, g;
00505 size_t idx;
00506 for (int k = z0 - 2; k <= z0 + 2; k++) {
00507 for (int j = y0 - 2; j <= y0 + 2; j++) {
00508 for (int i = x0 - 4; i <= x0 + 4; i += 2) {
00509 r = sqrt(Util::hypot3sq(i / 2.0f - x, j - y, k - z));
00510 g = Interp::hyperg(r);
00511 n += g;
00512 idx = i + j * nx + (size_t)k * nx * ny;
00513 data[0] += rdata[idx] * g;
00514 data[1] += rdata[idx + 1] * g;
00515 }
00516
00517 }
00518 }
00519 }
00520 data[0] /= n;
00521 data[1] /= n;
00522 return true;
00523 }
00524 return false;
00525
00526 }
00527
00528 return false;
00529 }
00530
00531 void PawelProjector::prepcubes(int nx, int ny, int nz, int ri, Vec3i origin,
00532 int& nn, IPCube* ipcube) const {
00533 const float r = float(ri)*float(ri);
00534 const int ldpx = origin[0];
00535 const int ldpy = origin[1];
00536 const int ldpz = origin[2];
00537
00538 float t;
00539 nn = -1;
00540 for (int i1 = 0; i1 < nz; i1++) {
00541 t = float(i1 - ldpz);
00542 const float xx = t*t;
00543 for (int i2 = 0; i2 < ny; i2++) {
00544 t = float(i2 - ldpy);
00545 const float yy = t*t + xx;
00546 bool first = true;
00547 for (int i3 = 0; i3 < nx; i3++) {
00548 t = float(i3 - ldpx);
00549 const float rc = t*t + yy;
00550 if (first) {
00551
00552 if (rc > r) continue;
00553 first = false;
00554 nn++;
00555 if (ipcube != NULL) {
00556 ipcube[nn].start = i3;
00557 ipcube[nn].end = i3;
00558 ipcube[nn].loc[0] = i3 - ldpx;
00559 ipcube[nn].loc[1] = i2 - ldpy;
00560 ipcube[nn].loc[2] = i1 - ldpz;
00561
00562 }
00563 } else {
00564
00565 if (ipcube != NULL) {
00566 if (rc <= r) ipcube[nn].end = i3;
00567 }
00568 }
00569 }
00570 }
00571 }
00572 }
00573
00574
00575 EMData *PawelProjector::project3d(EMData * image) const
00576 {
00577 if (!image) {
00578 return 0;
00579 }
00580 int ri;
00581 int nx = image->get_xsize();
00582 int ny = image->get_ysize();
00583 int nz = image->get_zsize();
00584 int dim = Util::get_max(nx,ny,nz);
00585 if (nz == 1) {
00586 LOGERR("The PawelProjector needs a volume!");
00587 return 0;
00588 }
00589
00590 Vec3i origin(0,0,0);
00591
00592
00593 if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
00594 else {origin[0] = nx/2;}
00595 if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
00596 else {origin[1] = ny/2;}
00597 if (params.has_key("origin_z")) {origin[2] = params["origin_z"];}
00598 else {origin[2] = nz/2;}
00599
00600 if (params.has_key("radius")) {ri = params["radius"];}
00601 else {ri = dim/2 - 1;}
00602
00603
00604 int nn = -1;
00605 prepcubes(nx, ny, nz, ri, origin, nn);
00606
00607
00608 IPCube* ipcube = new IPCube[nn+1];
00609 prepcubes(nx, ny, nz, ri, origin, nn, ipcube);
00610
00611 Transform* rotation = params["transform"];
00612 int nangles = 0;
00613 vector<float> anglelist;
00614
00615
00616
00617
00618
00619
00620
00621 if ( rotation == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633 nangles = 1;
00634
00635
00636
00637 EMData* ret = new EMData();
00638 ret->set_size(nx, ny, nangles);
00639 ret->to_zero();
00640
00641
00642 for (int ia = 0; ia < nangles; ia++) {
00643
00644
00645
00646 if (2*(ri+1)+1 > dim) {
00647
00648 for (int i = 0 ; i <= nn; i++) {
00649 int k = ipcube[i].loc[1] + origin[1];
00650 Vec3f vb = ipcube[i].loc*(*rotation) + origin;
00651 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
00652
00653 int iox = int(vb[0]);
00654 if ((iox >= 0) && (iox < nx-1)) {
00655 int ioy = int(vb[1]);
00656 if ((ioy >= 0) && (ioy < ny-1)) {
00657 int ioz = int(vb[2]);
00658 if ((ioz >= 0) && (ioz < nz-1)) {
00659
00660
00661
00662 float dx = vb[0] - iox;
00663 float dy = vb[1] - ioy;
00664 float dz = vb[2] - ioz;
00665 float a1 = (*image)(iox,ioy,ioz);
00666 float a2 = (*image)(iox+1,ioy,ioz) - a1;
00667 float a3 = (*image)(iox,ioy+1,ioz) - a1;
00668 float a4 = (*image)(iox,ioy,ioz+1) - a1;
00669 float a5 = -a2 -(*image)(iox,ioy+1,ioz)
00670 + (*image)(iox+1,ioy+1,ioz);
00671 float a61 = -(*image)(iox,ioy,ioz+1)
00672 + (*image)(iox+1,ioy,ioz+1);
00673 float a6 = -a2 + a61;
00674 float a7 = -a3 - (*image)(iox,ioy,ioz+1)
00675 + (*image)(iox,ioy+1,ioz+1);
00676 float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
00677 + (*image)(iox+1,ioy+1,ioz+1);
00678 (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
00679 + (a7 + a8*dx)*dy)
00680 + a3*dy + dx*(a2 + a5*dy);
00681 }
00682 }
00683 }
00684 vb += rotation->get_matrix3_row(0);
00685 }
00686 }
00687
00688 } else {
00689
00690 for (int i = 0 ; i <= nn; i++) {
00691 int k = ipcube[i].loc[1] + origin[1];
00692 Vec3f vb = ipcube[i].loc*(*rotation) + origin;
00693 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
00694 int iox = int(vb[0]);
00695 int ioy = int(vb[1]);
00696 int ioz = int(vb[2]);
00697 float dx = vb[0] - iox;
00698 float dy = vb[1] - ioy;
00699 float dz = vb[2] - ioz;
00700 float a1 = (*image)(iox,ioy,ioz);
00701 float a2 = (*image)(iox+1,ioy,ioz) - a1;
00702 float a3 = (*image)(iox,ioy+1,ioz) - a1;
00703 float a4 = (*image)(iox,ioy,ioz+1) - a1;
00704 float a5 = -a2 -(*image)(iox,ioy+1,ioz)
00705 + (*image)(iox+1,ioy+1,ioz);
00706 float a61 = -(*image)(iox,ioy,ioz+1)
00707 + (*image)(iox+1,ioy,ioz+1);
00708 float a6 = -a2 + a61;
00709 float a7 = -a3 - (*image)(iox,ioy,ioz+1)
00710 + (*image)(iox,ioy+1,ioz+1);
00711 float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
00712 + (*image)(iox+1,ioy+1,ioz+1);
00713 (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
00714 + (a7 + a8*dx)*dy)
00715 + a3*dy + dx*(a2 + a5*dy);
00716 vb += rotation->get_matrix3_row(0);
00717 }
00718 }
00719 }
00720 }
00721 ret->update();
00722 EMDeleteArray(ipcube);
00723 if(rotation) {delete rotation; rotation=0;}
00724 return ret;
00725 }
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877 EMData *StandardProjector::project3d(EMData * image) const
00878 {
00879 Transform* t3d = params["transform"];
00880 if ( t3d == NULL ) throw NullPointerException("The transform object containing the angles(required for projection), was not specified");
00881
00882 if ( image->get_ndim() == 3 )
00883 {
00884
00885 #ifdef EMAN2_USING_CUDA
00886 if(EMData::usecuda == 1) {
00887 if(!image->isrodataongpu()) image->copy_to_cudaro();
00888
00889 Transform* t3d = params["transform"];
00890 if ( t3d == NULL ) throw NullPointerException("The transform object containing the angles(required for projection), was not specified");
00891 float * m = new float[12];
00892 t3d->copy_matrix_into_array(m);
00893 image->bindcudaarrayA(true);
00894
00895 EMData *e = new EMData();
00896 e->set_size_cuda(image->get_xsize(), image->get_ysize(), 1);
00897 e->rw_alloc();
00898 standard_project(m,e->getcudarwdata(), e->get_xsize(), e->get_ysize());
00899 image->unbindcudaarryA();
00900 delete [] m;
00901
00902 e->update();
00903 e->set_attr("xform.projection",t3d);
00904 e->set_attr("apix_x",(float)image->get_attr("apix_x"));
00905 e->set_attr("apix_y",(float)image->get_attr("apix_y"));
00906 e->set_attr("apix_z",(float)image->get_attr("apix_z"));
00907
00908 if(t3d) {delete t3d; t3d=0;}
00909 return e;
00910 }
00911 #endif
00912 int nx = image->get_xsize();
00913 int ny = image->get_ysize();
00914 int nz = image->get_zsize();
00915
00916
00917 Transform r = t3d->inverse();
00918 int xy = nx * ny;
00919
00920 EMData *proj = new EMData();
00921 proj->set_size(nx, ny, 1);
00922
00923 Vec3i offset(nx/2,ny/2,nz/2);
00924
00925 float *sdata = image->get_data();
00926 float *ddata = proj->get_data();
00927 for (int k = -nz / 2; k < nz - nz / 2; k++) {
00928 int l = 0;
00929 for (int j = -ny / 2; j < ny - ny / 2; j++) {
00930 ddata[l]=0;
00931 for (int i = -nx / 2; i < nx - nx / 2; i++,l++) {
00932
00933 Vec3f coord(i,j,k);
00934 Vec3f soln = r*coord;
00935 soln += offset;
00936
00940
00941
00942 float x2 = soln[0];
00943 float y2 = soln[1];
00944 float z2 = soln[2];
00945
00946 float x = (float)Util::fast_floor(x2);
00947 float y = (float)Util::fast_floor(y2);
00948 float z = (float)Util::fast_floor(z2);
00949
00950 float t = x2 - x;
00951 float u = y2 - y;
00952 float v = z2 - z;
00953
00954 size_t ii = (size_t) ((size_t)x + (size_t)y * nx + (size_t)z * xy);
00955
00956 if (x2 < 0 || y2 < 0 || z2 < 0 ) continue;
00957 if (x2 > (nx-1) || y2 > (ny-1) || z2 > (nz-1) ) continue;
00958
00959 if (x2 < (nx - 1) && y2 < (ny - 1) && z2 < (nz - 1)) {
00960 ddata[l] +=
00961 Util::trilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx],
00962 sdata[ii + nx + 1], sdata[ii + xy], sdata[ii + xy + 1], sdata[ii + xy + nx],
00963 sdata[ii + xy + nx + 1], t, u, v);
00964 }
00965 else if ( x2 == (nx - 1) && y2 == (ny - 1) && z2 == (nz - 1) ) {
00966 ddata[l] += sdata[ii];
00967 }
00968 else if ( x2 == (nx - 1) && y2 == (ny - 1) ) {
00969 ddata[l] += Util::linear_interpolate(sdata[ii], sdata[ii + xy],v);
00970 }
00971 else if ( x2 == (nx - 1) && z2 == (nz - 1) ) {
00972 ddata[l] += Util::linear_interpolate(sdata[ii], sdata[ii + nx],u);
00973 }
00974 else if ( y2 == (ny - 1) && z2 == (nz - 1) ) {
00975 ddata[l] += Util::linear_interpolate(sdata[ii], sdata[ii + 1],t);
00976 }
00977 else if ( x2 == (nx - 1) ) {
00978 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + nx], sdata[ii + xy], sdata[ii + xy + nx],u,v);
00979 }
00980 else if ( y2 == (ny - 1) ) {
00981 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + xy], sdata[ii + xy + 1],t,v);
00982 }
00983 else if ( z2 == (nz - 1) ) {
00984 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx], sdata[ii + nx + 1],t,u);
00985 }
00986 }
00987 }
00988 }
00989 proj->update();
00990 proj->set_attr("xform.projection",t3d);
00991 proj->set_attr("apix_x",(float)image->get_attr("apix_x"));
00992 proj->set_attr("apix_y",(float)image->get_attr("apix_y"));
00993 proj->set_attr("apix_z",(float)image->get_attr("apix_z"));
00994
00995 if(t3d) {delete t3d; t3d=0;}
00996 return proj;
00997 }
00998 else if ( image->get_ndim() == 2 ) {
00999
01000 Transform r = t3d->inverse();
01001
01002 int nx = image->get_xsize();
01003 int ny = image->get_ysize();
01004
01005 EMData *proj = new EMData();
01006 proj->set_size(nx, 1, 1);
01007 proj->to_zero();
01008
01009 float *sdata = image->get_data();
01010 float *ddata = proj->get_data();
01011
01012 Vec2f offset(nx/2,ny/2);
01013 for (int j = -ny / 2; j < ny - ny / 2; j++) {
01014 int l = 0;
01015 for (int i = -nx / 2; i < nx - nx / 2; i++,l++) {
01016
01017 Vec2f coord(i,j);
01018 Vec2f soln = r*coord;
01019 soln += offset;
01020
01021 float x2 = soln[0];
01022 float y2 = soln[1];
01023
01024 float x = (float)Util::fast_floor(x2);
01025 float y = (float)Util::fast_floor(y2);
01026
01027 int ii = (int) (x + y * nx);
01028 float u = x2 - x;
01029 float v = y2 - y;
01030
01031 if (x2 < 0 || y2 < 0 ) continue;
01032 if (x2 > (nx-1) || y2 > (ny-1) ) continue;
01033
01034 if ( x2 < (nx - 1) && y2 < (ny - 1) ) {
01035 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx],sdata[ii + nx + 1], u, v);
01036 }
01037 else if (x2 == (nx-1) && y2 == (ny-1) ) {
01038 ddata[l] += sdata[ii];
01039 }
01040 else if (x2 == (nx-1) ) {
01041 ddata[l] += Util::linear_interpolate(sdata[ii],sdata[ii + nx], v);
01042 }
01043 else if (y2 == (ny-1) ) {
01044 ddata[l] += Util::linear_interpolate(sdata[ii],sdata[ii + 1], u);
01045 }
01046 }
01047 }
01048 proj->set_attr("xform.projection",t3d);
01049 proj->update();
01050 if(t3d) {delete t3d; t3d=0;}
01051 return proj;
01052 }
01053 else throw ImageDimensionException("Standard projection works only for 2D and 3D images");
01054 }
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149 EMData *FourierGriddingProjector::project3d(EMData * image) const
01150 {
01151 if (!image) {
01152 return 0;
01153 }
01154 if (3 != image->get_ndim())
01155 throw ImageDimensionException(
01156 "FourierGriddingProjector needs a 3-D volume");
01157 if (image->is_complex())
01158 throw ImageFormatException(
01159 "FourierGriddingProjector requires a real volume");
01160 const int npad = params.has_key("npad") ? int(params["npad"]) : 2;
01161 const int nx = image->get_xsize();
01162 const int ny = image->get_ysize();
01163 const int nz = image->get_zsize();
01164 if (nx != ny || nx != nz)
01165 throw ImageDimensionException(
01166 "FourierGriddingProjector requires nx==ny==nz");
01167 const int m = Util::get_min(nx,ny,nz);
01168 const int n = m*npad;
01169
01170 int K = params["kb_K"];
01171 if ( K == 0 ) K = 6;
01172 float alpha = params["kb_alpha"];
01173 if ( alpha == 0 ) alpha = 1.25;
01174 Util::KaiserBessel kb(alpha, K, (float)(m/2), K/(2.0f*n), n);
01175
01176
01177 EMData* tmpImage = image->copy();
01178 tmpImage->divkbsinh(kb);
01179
01180
01181
01182 EMData* imgft = tmpImage->norm_pad(false, npad);
01183 imgft->do_fft_inplace();
01184 imgft->center_origin_fft();
01185 imgft->fft_shuffle();
01186 delete tmpImage;
01187
01188
01189 int nangles = 0;
01190 vector<float> anglelist;
01191
01192 if (params.has_key("anglelist")) {
01193 anglelist = params["anglelist"];
01194 nangles = anglelist.size() / 3;
01195 } else {
01196
01197
01198
01199
01200 Transform* t3d = params["transform"];
01201 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
01202 Dict p = t3d->get_rotation("spider");
01203
01204 string angletype = "SPIDER";
01205 float phi = p["phi"];
01206 float theta = p["theta"];
01207 float psi = p["psi"];
01208 anglelist.push_back(phi);
01209 anglelist.push_back(theta);
01210 anglelist.push_back(psi);
01211 nangles = 1;
01212 if(t3d) {delete t3d; t3d=0;}
01213 }
01214
01215
01216
01217
01218 EMData* ret = new EMData();
01219 ret->set_size(nx, ny, nangles);
01220 ret->to_zero();
01221
01222 for (int ia = 0; ia < nangles; ia++) {
01223 int indx = 3*ia;
01224 Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]);
01225 Transform tf(d);
01226 EMData* proj = imgft->extract_plane(tf, kb);
01227 if (proj->is_shuffled()) proj->fft_shuffle();
01228 proj->center_origin_fft();
01229 proj->do_ift_inplace();
01230 EMData* winproj = proj->window_center(m);
01231 delete proj;
01232 for (int iy=0; iy < ny; iy++)
01233 for (int ix=0; ix < nx; ix++)
01234 (*ret)(ix,iy,ia) = (*winproj)(ix,iy);
01235 delete winproj;
01236 }
01237 delete imgft;
01238
01239 if (!params.has_key("anglelist")) {
01240 Transform* t3d = params["transform"];
01241 ret->set_attr("xform.projection",t3d);
01242 if(t3d) {delete t3d; t3d=0;}
01243 }
01244 ret->update();
01245 return ret;
01246 }
01247
01248
01249 int ChaoProjector::getnnz(Vec3i volsize, int ri, Vec3i origin, int *nrays, int *nnz) const
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263 {
01264 int ix, iy, iz, rs, r2, xs, ys, zs, xx, yy, zz;
01265 int ftm=0, status = 0;
01266
01267 r2 = ri*ri;
01268 *nnz = 0;
01269 *nrays = 0;
01270 int nx = (int)volsize[0];
01271 int ny = (int)volsize[1];
01272 int nz = (int)volsize[2];
01273
01274 int xcent = (int)origin[0];
01275 int ycent = (int)origin[1];
01276 int zcent = (int)origin[2];
01277
01278
01279 for (ix = 1; ix <=nx; ix++) {
01280 xs = ix-xcent;
01281 xx = xs*xs;
01282 for (iy = 1; iy <= ny; iy++) {
01283 ys = iy-ycent;
01284 yy = ys*ys;
01285 ftm = 1;
01286 for (iz = 1; iz <= nz; iz++) {
01287 zs = iz-zcent;
01288 zz = zs*zs;
01289 rs = xx + yy + zz;
01290 if (rs <= r2) {
01291 (*nnz)++;
01292 if (ftm) {
01293 (*nrays)++;
01294 ftm = 0;
01295 }
01296 }
01297 }
01298 }
01299 }
01300 return status;
01301 }
01302
01303 #define cube(i,j,k) cube[ ((k-1)*ny + j-1)*nx + i-1 ]
01304 #define sphere(i) sphere[(i)-1]
01305 #define cord(i,j) cord[((j)-1)*3 + (i) -1]
01306 #define ptrs(i) ptrs[(i)-1]
01307 #define dm(i) dm[(i)-1]
01308
01309 int ChaoProjector:: cb2sph(float *cube, Vec3i volsize, int ri, Vec3i origin,
01310 int nnz0, int *ptrs, int *cord, float *sphere) const
01311 {
01312 int xs, ys, zs, xx, yy, zz, rs, r2;
01313 int ix, iy, iz, jnz, nnz, nrays;
01314 int ftm = 0, status = 0;
01315
01316 int xcent = (int)origin[0];
01317 int ycent = (int)origin[1];
01318 int zcent = (int)origin[2];
01319
01320 int nx = (int)volsize[0];
01321 int ny = (int)volsize[1];
01322 int nz = (int)volsize[2];
01323
01324 r2 = ri*ri;
01325 nnz = 0;
01326 nrays = 0;
01327 ptrs(1) = 1;
01328
01329 for (ix = 1; ix <= nx; ix++) {
01330 xs = ix-xcent;
01331 xx = xs*xs;
01332 for ( iy = 1; iy <= ny; iy++ ) {
01333 ys = iy-ycent;
01334 yy = ys*ys;
01335 jnz = 0;
01336
01337 ftm = 1;
01338
01339 for (iz = 1; iz <= nz; iz++) {
01340 zs = iz-zcent;
01341 zz = zs*zs;
01342 rs = xx + yy + zz;
01343 if (rs <= r2) {
01344 jnz++;
01345 nnz++;
01346 sphere(nnz) = cube(iz, iy, ix);
01347
01348
01349 if (ftm) {
01350 nrays++;
01351 cord(1,nrays) = iz;
01352 cord(2,nrays) = iy;
01353 cord(3,nrays) = ix;
01354 ftm = 0;
01355 }
01356 }
01357 }
01358 if (jnz > 0) {
01359 ptrs(nrays+1) = ptrs(nrays) + jnz;
01360 }
01361 }
01362 }
01363 if (nnz != nnz0) status = -1;
01364 return status;
01365 }
01366
01367
01368 int ChaoProjector::sph2cb(float *sphere, Vec3i volsize, int nrays, int ri,
01369 int nnz0, int *ptrs, int *cord, float *cube) const
01370 {
01371 int status=0;
01372 int r2, i, j, ix, iy, iz, nnz;
01373
01374 int nx = (int)volsize[0];
01375 int ny = (int)volsize[1];
01376
01377
01378 r2 = ri*ri;
01379 nnz = 0;
01380 ptrs(1) = 1;
01381
01382
01383
01384
01385 nnz = 0;
01386 for (j = 1; j <= nrays; j++) {
01387 iz = cord(1,j);
01388 iy = cord(2,j);
01389 ix = cord(3,j);
01390 for (i = ptrs(j); i<=ptrs(j+1)-1; i++, iz++) {
01391 nnz++;
01392 cube(iz,iy,ix) = sphere(nnz);
01393 }
01394 }
01395 if (nnz != nnz0) status = -1;
01396 return status;
01397 }
01398
01399 #define x(i) x[(i)-1]
01400 #define y(i,j) y[(j-1)*nx + i - 1]
01401
01402
01403 int ChaoProjector::fwdpj3(Vec3i volsize, int nrays, int , float *dm,
01404 Vec3i origin, int ri, int *ptrs, int *cord,
01405 float *x, float *y) const
01406 {
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423 int iqx, iqy, i, j, xc, yc, zc;
01424 float ct, dipx, dipy, dipx1m, dipy1m, xb, yb, dm1, dm4;
01425 int status = 0;
01426
01427 int xcent = origin[0];
01428 int ycent = origin[1];
01429 int zcent = origin[2];
01430
01431 int nx = volsize[0];
01432
01433 dm1 = dm(1);
01434 dm4 = dm(4);
01435
01436 if ( nx > 2*ri ) {
01437 for (i = 1; i <= nrays; i++) {
01438 zc = cord(1,i)-zcent;
01439 yc = cord(2,i)-ycent;
01440 xc = cord(3,i)-xcent;
01441
01442 xb = zc*dm(1)+yc*dm(2)+xc*dm(3) + xcent;
01443 yb = zc*dm(4)+yc*dm(5)+xc*dm(6) + ycent;
01444
01445 for (j = ptrs(i); j< ptrs(i+1); j++) {
01446 iqx = ifix(xb);
01447 iqy = ifix(yb);
01448
01449 ct = x(j);
01450 dipx = xb - (float)(iqx);
01451 dipy = (yb - (float)(iqy)) * ct;
01452
01453 dipy1m = ct - dipy;
01454 dipx1m = 1.0f - dipx;
01455
01456 y(iqx ,iqy) = y(iqx ,iqy) + dipx1m*dipy1m;
01457 y(iqx+1,iqy) = y(iqx+1,iqy) + dipx*dipy1m;
01458 y(iqx+1,iqy+1) = y(iqx+1,iqy+1) + dipx*dipy;
01459 y(iqx ,iqy+1) = y(iqx ,iqy+1) + dipx1m*dipy;
01460
01461 xb += dm1;
01462 yb += dm4;
01463 }
01464 }
01465 }
01466 else {
01467 fprintf(stderr, " nx must be greater than 2*ri\n");
01468 exit(1);
01469 }
01470 return status;
01471 }
01472 #undef x
01473 #undef y
01474
01475 #define y(i) y[(i)-1]
01476 #define x(i,j) x[((j)-1)*nx + (i) - 1]
01477
01478
01479 int ChaoProjector::bckpj3(Vec3i volsize, int nrays, int , float *dm,
01480 Vec3i origin, int ri, int *ptrs, int *cord,
01481 float *x, float *y) const
01482 {
01483 int i, j, iqx,iqy, xc, yc, zc;
01484 float xb, yb, dx, dy, dx1m, dy1m, dxdy;
01485 int status = 0;
01486
01487 int xcent = origin[0];
01488 int ycent = origin[1];
01489 int zcent = origin[2];
01490
01491 int nx = volsize[0];
01492
01493 if ( nx > 2*ri) {
01494 for (i = 1; i <= nrays; i++) {
01495 zc = cord(1,i) - zcent;
01496 yc = cord(2,i) - ycent;
01497 xc = cord(3,i) - xcent;
01498
01499 xb = zc*dm(1)+yc*dm(2)+xc*dm(3) + xcent;
01500 yb = zc*dm(4)+yc*dm(5)+xc*dm(6) + ycent;
01501
01502 for (j = ptrs(i); j <ptrs(i+1); j++) {
01503 iqx = ifix((float)(xb));
01504 iqy = ifix((float)(yb));
01505
01506 dx = xb - (float)(iqx);
01507 dy = yb - (float)(iqy);
01508 dx1m = 1.0f - dx;
01509 dy1m = 1.0f - dy;
01510 dxdy = dx*dy;
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529 y(j) += x(iqx,iqy)
01530 + dx*(-x(iqx,iqy)+x(iqx+1,iqy))
01531 + dy*(-x(iqx,iqy)+x(iqx,iqy+1))
01532 + dxdy*( x(iqx,iqy) - x(iqx,iqy+1)
01533 -x(iqx+1,iqy) + x(iqx+1,iqy+1) );
01534
01535 xb += dm(1);
01536 yb += dm(4);
01537 }
01538 }
01539 }
01540 else {
01541 fprintf(stderr, "bckpj3: nx must be greater than 2*ri\n");
01542 }
01543
01544 return status;
01545 }
01546
01547 #undef x
01548 #undef y
01549 #undef dm
01550
01551
01552 int ChaoProjector::ifix(float a) const
01553 {
01554 int ia;
01555
01556 if (a>=0) {
01557 ia = (int)floor(a);
01558 }
01559 else {
01560 ia = (int)ceil(a);
01561 }
01562 return ia;
01563 }
01564
01565 #define dm(i,j) dm[((j)-1)*9 + (i) -1]
01566 #define anglelist(i,j) anglelist[((j)-1)*3 + (i) - 1]
01567
01568
01569 void ChaoProjector::setdm(vector<float> anglelist, string const , float *dm) const
01570 {
01571
01572 float psi, theta, phi;
01573 double cthe, sthe, cpsi, spsi, cphi, sphi;
01574 int j;
01575
01576 int nangles = anglelist.size() / 3;
01577
01578
01579 for (j = 1; j <= nangles; j++) {
01580 phi = static_cast<float>(anglelist(1,j)*dgr_to_rad);
01581 theta = static_cast<float>(anglelist(2,j)*dgr_to_rad);
01582 psi = static_cast<float>(anglelist(3,j)*dgr_to_rad);
01583
01584
01585 cthe = cos(theta);
01586 sthe = sin(theta);
01587 cpsi = cos(psi);
01588 spsi = sin(psi);
01589 cphi = cos(phi);
01590 sphi = sin(phi);
01591
01592 dm(1,j)=static_cast<float>(cphi*cthe*cpsi-sphi*spsi);
01593 dm(2,j)=static_cast<float>(sphi*cthe*cpsi+cphi*spsi);
01594 dm(3,j)=static_cast<float>(-sthe*cpsi);
01595 dm(4,j)=static_cast<float>(-cphi*cthe*spsi-sphi*cpsi);
01596 dm(5,j)=static_cast<float>(-sphi*cthe*spsi+cphi*cpsi);
01597 dm(6,j)=static_cast<float>(sthe*spsi);
01598 dm(7,j)=static_cast<float>(sthe*cphi);
01599 dm(8,j)=static_cast<float>(sthe*sphi);
01600 dm(9,j)=static_cast<float>(cthe);
01601 }
01602 }
01603 #undef anglelist
01604
01605 #define images(i,j,k) images[ ((k-1)*nyvol + j-1)*nxvol + i-1 ]
01606
01607 EMData *ChaoProjector::project3d(EMData * vol) const
01608 {
01609
01610 int nrays, nnz, status, j;
01611 float *dm;
01612 int *ptrs, *cord;
01613 float *sphere, *images;
01614
01615 int nxvol = vol->get_xsize();
01616 int nyvol = vol->get_ysize();
01617 int nzvol = vol->get_zsize();
01618 Vec3i volsize(nxvol,nyvol,nzvol);
01619
01620 int dim = Util::get_min(nxvol,nyvol,nzvol);
01621 if (nzvol == 1) {
01622 LOGERR("The ChaoProjector needs a volume!");
01623 return 0;
01624 }
01625 Vec3i origin(0,0,0);
01626
01627
01628 if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
01629 else {origin[0] = nxvol/2+1;}
01630 if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
01631 else {origin[1] = nyvol/2+1;}
01632 if (params.has_key("origin_z")) {origin[2] = params["origin_z"];}
01633 else {origin[2] = nzvol/2+1;}
01634
01635 int ri;
01636 if (params.has_key("radius")) {ri = params["radius"];}
01637 else {ri = dim/2 - 1;}
01638
01639
01640 float *cube = vol->get_data();
01641
01642
01643
01644 status = getnnz(volsize, ri, origin, &nrays, &nnz);
01645
01646
01647
01648 sphere = new float[nnz];
01649 ptrs = new int[nrays+1];
01650 cord = new int[3*nrays];
01651 if (sphere == NULL || ptrs == NULL || cord == NULL) {
01652 fprintf(stderr,"ChaoProjector::project3d, failed to allocate!\n");
01653 exit(1);
01654 }
01655 for (int i = 0; i<nnz; i++) sphere[i] = 0.0;
01656 for (int i = 0; i<nrays+1; i++) ptrs[i] = 0;
01657 for (int i = 0; i<3*nrays; i++) cord[i] = 0;
01658
01659 status = cb2sph(cube, volsize, ri, origin, nnz, ptrs, cord, sphere);
01660
01661
01662 int nangles = 0;
01663 vector<float> anglelist;
01664 string angletype = "SPIDER";
01665
01666 if (params.has_key("anglelist")) {
01667 anglelist = params["anglelist"];
01668 nangles = anglelist.size() / 3;
01669 } else {
01670 Transform* t3d = params["transform"];
01671 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
01672
01673
01674
01675
01676 Dict p = t3d->get_rotation("spider");
01677 if(t3d) {delete t3d; t3d=0;}
01678
01679 float phi = p["phi"];
01680 float theta = p["theta"];
01681 float psi = p["psi"];
01682 anglelist.push_back(phi);
01683 anglelist.push_back(theta);
01684 anglelist.push_back(psi);
01685 nangles = 1;
01686 }
01687
01688
01689 dm = new float[nangles*9];
01690 setdm(anglelist, angletype, dm);
01691
01692
01693 EMData *ret = new EMData();
01694 ret->set_size(nxvol, nyvol, nangles);
01695 ret->set_complex(false);
01696 ret->set_ri(true);
01697
01698 images = ret->get_data();
01699
01700 for (j = 1; j <= nangles; j++) {
01701 status = fwdpj3(volsize, nrays, nnz , &dm(1,j), origin, ri,
01702 ptrs , cord, sphere, &images(1,1,j));
01703
01704 }
01705
01706
01707 EMDeleteArray(dm);
01708 EMDeleteArray(ptrs);
01709 EMDeleteArray(cord);
01710 EMDeleteArray(sphere);
01711
01712 if (!params.has_key("anglelist")) {
01713 Transform* t3d = params["transform"];
01714 ret->set_attr("xform.projection",t3d);
01715 if(t3d) {delete t3d; t3d=0;}
01716 }
01717 ret->update();
01718 return ret;
01719 }
01720
01721
01722 #undef images
01723
01724 #define images(i,j,k) images[ ((k)-1)*nximg*nyimg + ((j)-1)*nximg + (i)-1 ]
01725
01726 EMData *ChaoProjector::backproject3d(EMData * imagestack) const
01727 {
01728 int nrays, nnz, status, j;
01729 float *dm;
01730 int *ptrs, *cord;
01731 float *sphere, *images, *cube;
01732
01733 int nximg = imagestack->get_xsize();
01734 int nyimg = imagestack->get_ysize();
01735 int nslices = imagestack->get_zsize();
01736
01737 int dim = Util::get_min(nximg,nyimg);
01738 Vec3i volsize(nximg,nyimg,dim);
01739
01740 Vec3i origin(0,0,0);
01741
01742
01743 if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
01744 else {origin[0] = nximg/2+1;}
01745 if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
01746 else {origin[1] = nyimg/2+1;}
01747 if (params.has_key("origin_z")) {origin[1] = params["origin_z"];}
01748 else {origin[2] = dim/2+1;}
01749
01750 int ri;
01751 if (params.has_key("radius")) {ri = params["radius"];}
01752 else {ri = dim/2 - 1;}
01753
01754
01755 images = imagestack->get_data();
01756
01757
01758
01759 status = getnnz(volsize, ri, origin, &nrays, &nnz);
01760
01761
01762
01763 sphere = new float[nnz];
01764 ptrs = new int[nrays+1];
01765 cord = new int[3*nrays];
01766 if (sphere == NULL || ptrs == NULL || cord == NULL) {
01767 fprintf(stderr,"ChaoProjector::backproject3d, failed to allocate!\n");
01768 exit(1);
01769 }
01770 for (int i = 0; i<nnz; i++) sphere[i] = 0.0;
01771 for (int i = 0; i<nrays+1; i++) ptrs[i] = 0;
01772 for (int i = 0; i<3*nrays; i++) cord[i] = 0;
01773
01774 int nangles = 0;
01775 vector<float> anglelist;
01776 string angletype = "SPIDER";
01777
01778 if (params.has_key("anglelist")) {
01779 anglelist = params["anglelist"];
01780 nangles = anglelist.size() / 3;
01781 } else {
01782 Transform* t3d = params["transform"];
01783 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
01784
01785
01786
01787
01788
01789 Dict p = t3d->get_rotation("spider");
01790 if(t3d) {delete t3d; t3d=0;}
01791
01792 float phi = p["phi"];
01793 float theta = p["theta"];
01794 float psi = p["psi"];
01795 anglelist.push_back(phi);
01796 anglelist.push_back(theta);
01797 anglelist.push_back(psi);
01798 nangles = 1;
01799 }
01800
01801
01802
01803 if (nslices != nangles) {
01804 LOGERR("the number of images does not match the number of angles");
01805 return 0;
01806 }
01807
01808 dm = new float[nangles*9];
01809 setdm(anglelist, angletype, dm);
01810
01811
01812 EMData *ret = new EMData();
01813 ret->set_size(nximg, nyimg, dim);
01814 ret->set_complex(false);
01815 ret->set_ri(true);
01816 ret->to_zero();
01817
01818 cube = ret->get_data();
01819
01820 status = cb2sph(cube, volsize, ri, origin, nnz, ptrs, cord, sphere);
01821
01822
01823 for (j = 1; j <= nangles; j++) {
01824 status = bckpj3(volsize, nrays, nnz, &dm(1,j), origin, ri,
01825 ptrs , cord , &images(1,1,j), sphere);
01826
01827 }
01828
01829 status = sph2cb(sphere, volsize, nrays, ri, nnz, ptrs, cord, cube);
01830
01831
01832
01833 EMDeleteArray(dm);
01834 EMDeleteArray(ptrs);
01835 EMDeleteArray(cord);
01836 EMDeleteArray(sphere);
01837
01838 ret->update();
01839 return ret;
01840 }
01841
01842 #undef images
01843 #undef cube
01844 #undef sphere
01845 #undef cord
01846 #undef ptrs
01847 #undef dm
01848
01849 EMData *GaussFFTProjector::backproject3d(EMData * ) const
01850 {
01851
01852 EMData *ret = new EMData();
01853 return ret;
01854 }
01855
01856 #define images(i,j,k) images[ (k)*nx*ny + ((j)-1)*nx + (i)-1 ]
01857
01858
01859
01860
01861
01862
01863
01864
01865
01866
01867
01868
01869
01870
01871
01872
01873
01874
01875
01876
01877
01878
01879
01880
01881
01882
01883
01884
01885
01886
01887
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897
01898
01899
01900
01901
01902
01903
01904
01905
01906
01907
01908
01909
01910
01911
01912
01913
01914
01915
01916
01917
01918
01919
01920
01921
01922
01923
01924
01925
01926
01927
01928
01929
01930
01931
01932
01933
01934
01935
01936
01937
01938
01939
01940
01941
01942
01943
01944
01945
01946
01947
01948
01949
01950
01951
01952
01953
01954
01955
01956
01957
01958
01959
01960
01961
01962
01963
01964
01965
01966
01967
01968
01969
01970
01971
01972 EMData *PawelProjector::backproject3d(EMData * imagestack) const
01973 {
01974
01975 float *images;
01976
01977 if (!imagestack) {
01978 return 0;
01979 }
01980 int ri;
01981 int nx = imagestack->get_xsize();
01982 int ny = imagestack->get_ysize();
01983
01984 int dim = Util::get_min(nx,ny);
01985 images = imagestack->get_data();
01986
01987 Vec3i origin(0,0,0);
01988
01989
01990 if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
01991 else {origin[0] = nx/2;}
01992 if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
01993 else {origin[1] = ny/2;}
01994 if (params.has_key("origin_z")) {origin[1] = params["origin_z"];}
01995 else {origin[2] = dim/2;}
01996
01997 if (params.has_key("radius")) {ri = params["radius"];}
01998 else {ri = dim/2 - 1;}
01999
02000
02001 int nn = -1;
02002 prepcubes(nx, ny, dim, ri, origin, nn);
02003
02004
02005 IPCube* ipcube = new IPCube[nn+1];
02006 prepcubes(nx, ny, dim, ri, origin, nn, ipcube);
02007
02008 int nangles = 0;
02009 vector<float> anglelist;
02010
02011 if (params.has_key("anglelist")) {
02012 anglelist = params["anglelist"];
02013 nangles = anglelist.size() / 3;
02014 } else {
02015 Transform* t3d = params["transform"];
02016 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
02017
02018
02019
02020
02021 Dict p = t3d->get_rotation("spider");
02022 if(t3d) {delete t3d; t3d=0;}
02023
02024 string angletype = "SPIDER";
02025 float phi = p["phi"];
02026 float theta = p["theta"];
02027 float psi = p["psi"];
02028 anglelist.push_back(phi);
02029 anglelist.push_back(theta);
02030 anglelist.push_back(psi);
02031 nangles = 1;
02032 }
02033
02034
02035
02036
02037 EMData* ret = new EMData();
02038 ret->set_size(nx, ny, dim);
02039 ret->to_zero();
02040
02041
02042 for (int ia = 0; ia < nangles; ia++) {
02043 int indx = 3*ia;
02044 Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]);
02045 Transform rotation(d);
02046 float dm1 = rotation.at(0,0);
02047 float dm4 = rotation.at(1,0);
02048
02049 if (2*(ri+1)+1 > dim) {
02050
02051 LOGERR("backproject3d, pawel, 2*(ri+1)+1 > dim\n");
02052 return 0;
02053 } else {
02054
02055 for (int i = 0 ; i <= nn; i++) {
02056 int iox = (int)ipcube[i].loc[0]+origin[0];
02057 int ioy = (int)ipcube[i].loc[1]+origin[1];
02058 int ioz = (int)ipcube[i].loc[2]+origin[2];
02059
02060 Vec3f vb = rotation*ipcube[i].loc + origin;
02061 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
02062 float xbb = (j-ipcube[i].start)*dm1 + vb[0];
02063 int iqx = (int)floor(xbb);
02064
02065 float ybb = (j-ipcube[i].start)*dm4 + vb[1];
02066 int iqy = (int)floor(ybb);
02067
02068 float dipx = xbb - iqx;
02069 float dipy = ybb - iqy;
02070
02071 (*ret)(iox,ioy,ioz) += images(iqx,iqy,ia)
02072 + dipy*(images(iqx,iqy+1,ia)-images(iqx,iqy,ia))
02073 + dipx*(images(iqx+1,iqy,ia)-images(iqx,iqy,ia)
02074 + dipy*(images(iqx+1,iqy+1,ia)-images(iqx+1,iqy,ia)
02075 - images(iqx,iqy+1,ia)+images(iqx,iqy,ia)));
02076 iox++;
02077 }
02078 }
02079 }
02080 }
02081
02082 ret->update();
02083 EMDeleteArray(ipcube);
02084 return ret;
02085 }
02086 #undef images
02087
02088 EMData *StandardProjector::backproject3d(EMData * ) const
02089 {
02090
02091 EMData *ret = new EMData();
02092 return ret;
02093 }
02094
02095 EMData *FourierGriddingProjector::backproject3d(EMData * ) const
02096 {
02097
02098 EMData *ret = new EMData();
02099 return ret;
02100 }
02101
02102
02103
02104 void EMAN::dump_projectors()
02105 {
02106 dump_factory < Projector > ();
02107 }
02108
02109 map<string, vector<string> > EMAN::dump_projectors_list()
02110 {
02111 return dump_factory_list < Projector > ();
02112 }
02113
02114