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