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 using namespace std;
00047 using namespace EMAN;
00048
00049 const string GaussFFTProjector::NAME = "gauss_fft";
00050 const string FourierGriddingProjector::NAME = "fourier_gridding";
00051 const string PawelProjector::NAME = "pawel";
00052 const string StandardProjector::NAME = "standard";
00053 const string ChaoProjector::NAME = "chao";
00054
00055 #ifdef EMAN2_USING_CUDA
00056 const string CudaStandardProjector::NAME = "cuda_standard";
00057 #endif // EMAN2_USING_CUDA
00058
00059 template <> Factory < Projector >::Factory()
00060 {
00061 force_add<GaussFFTProjector>();
00062 force_add<PawelProjector>();
00063 force_add<StandardProjector>();
00064 force_add<FourierGriddingProjector>();
00065 force_add<ChaoProjector>();
00066 #ifdef EMAN2_USING_CUDA
00067 force_add<CudaStandardProjector>();
00068 #endif // EMAN2_USING_CUDA
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 + 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 + 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 + 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 + 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 + 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 + 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 + 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 + 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 + 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 + 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, int ny, int nz, int ri, Vec3i origin,
00534 int& nn, IPCube* ipcube) const {
00535 const float r = float(ri*ri);
00536 const int ldpx = origin[0];
00537 const int ldpy = origin[1];
00538 const int ldpz = origin[2];
00539 float t;
00540 nn = -1;
00541 for (int i1 = 0; i1 < nz; i1++) {
00542 t = float(i1 - ldpz);
00543 const float xx = t*t;
00544 for (int i2 = 0; i2 < ny; i2++) {
00545 t = float(i2 - ldpy);
00546 const float yy = t*t + xx;
00547 bool first = true;
00548 for (int i3 = 0; i3 < nz; i3++) {
00549 t = float(i3 - ldpx);
00550 const float rc = t*t + yy;
00551 if (first) {
00552
00553 if (rc > r) continue;
00554 first = false;
00555 nn++;
00556 if (ipcube != NULL) {
00557 ipcube[nn].start = i3;
00558 ipcube[nn].end = i3;
00559 ipcube[nn].loc[0] = i3 - ldpx;
00560 ipcube[nn].loc[1] = i2 - ldpy;
00561 ipcube[nn].loc[2] = i1 - ldpz;
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
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725 EMData *PawelProjector::project3d(EMData * image) const
00726 {
00727 if (!image) {
00728 return 0;
00729 }
00730 int ri;
00731 int nx = image->get_xsize();
00732 int ny = image->get_ysize();
00733 int nz = image->get_zsize();
00734 int dim = Util::get_min(nx,ny,nz);
00735 if (nz == 1) {
00736 LOGERR("The PawelProjector needs a volume!");
00737 return 0;
00738 }
00739 Vec3i origin(0,0,0);
00740
00741
00742 if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
00743 else {origin[0] = nx/2;}
00744 if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
00745 else {origin[1] = ny/2;}
00746 if (params.has_key("origin_z")) {origin[2] = params["origin_z"];}
00747 else {origin[2] = nz/2;}
00748
00749 if (params.has_key("radius")) {ri = params["radius"];}
00750 else {ri = dim/2 - 1;}
00751
00752
00753 int nn = -1;
00754 prepcubes(nx, ny, nz, ri, origin, nn);
00755
00756
00757 IPCube* ipcube = new IPCube[nn+1];
00758 prepcubes(nx, ny, nz, ri, origin, nn, ipcube);
00759
00760 Transform* rotation = params["transform"];
00761 int nangles = 0;
00762 vector<float> anglelist;
00763
00764
00765
00766
00767
00768
00769
00770 if ( rotation == NULL ) throw NullPointerException("The transform3d object (required for projection), was not specified");
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782 nangles = 1;
00783
00784
00785
00786 EMData* ret = new EMData();
00787 ret->set_size(nx, ny, nangles);
00788 ret->to_zero();
00789
00790
00791 for (int ia = 0; ia < nangles; ia++) {
00792
00793
00794
00795 if (2*(ri+1)+1 > dim) {
00796
00797 for (int i = 0 ; i <= nn; i++) {
00798 int k = ipcube[i].loc[1] + origin[1];
00799 Vec3f vb = ipcube[i].loc*(*rotation) + origin;
00800 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
00801
00802 int iox = int(vb[0]);
00803 if ((iox >= 0) && (iox < nx-1)) {
00804 int ioy = int(vb[1]);
00805 if ((ioy >= 0) && (ioy < ny-1)) {
00806 int ioz = int(vb[2]);
00807 if ((ioz >= 0) && (ioz < nz-1)) {
00808
00809 float dx = vb[0] - iox;
00810 float dy = vb[1] - ioy;
00811 float dz = vb[2] - ioz;
00812 float a1 = (*image)(iox,ioy,ioz);
00813 float a2 = (*image)(iox+1,ioy,ioz) - a1;
00814 float a3 = (*image)(iox,ioy+1,ioz) - a1;
00815 float a4 = (*image)(iox,ioy,ioz+1) - a1;
00816 float a5 = -a2 -(*image)(iox,ioy+1,ioz)
00817 + (*image)(iox+1,ioy+1,ioz);
00818 float a61 = -(*image)(iox,ioy,ioz+1)
00819 + (*image)(iox+1,ioy,ioz+1);
00820 float a6 = -a2 + a61;
00821 float a7 = -a3 - (*image)(iox,ioy,ioz+1)
00822 + (*image)(iox,ioy+1,ioz+1);
00823 float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
00824 + (*image)(iox+1,ioy+1,ioz+1);
00825 (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
00826 + (a7 + a8*dx)*dy)
00827 + a3*dy + dx*(a2 + a5*dy);
00828 }
00829 }
00830 }
00831 vb += rotation->get_matrix3_row(0);
00832 }
00833 }
00834
00835 } else {
00836
00837 for (int i = 0 ; i <= nn; i++) {
00838 int k = ipcube[i].loc[1] + origin[1];
00839 Vec3f vb = ipcube[i].loc*(*rotation) + origin;
00840 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
00841 int iox = int(vb[0]);
00842 int ioy = int(vb[1]);
00843 int ioz = int(vb[2]);
00844 float dx = vb[0] - iox;
00845 float dy = vb[1] - ioy;
00846 float dz = vb[2] - ioz;
00847 float a1 = (*image)(iox,ioy,ioz);
00848 float a2 = (*image)(iox+1,ioy,ioz) - a1;
00849 float a3 = (*image)(iox,ioy+1,ioz) - a1;
00850 float a4 = (*image)(iox,ioy,ioz+1) - a1;
00851 float a5 = -a2 -(*image)(iox,ioy+1,ioz)
00852 + (*image)(iox+1,ioy+1,ioz);
00853 float a61 = -(*image)(iox,ioy,ioz+1)
00854 + (*image)(iox+1,ioy,ioz+1);
00855 float a6 = -a2 + a61;
00856 float a7 = -a3 - (*image)(iox,ioy,ioz+1)
00857 + (*image)(iox,ioy+1,ioz+1);
00858 float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
00859 + (*image)(iox+1,ioy+1,ioz+1);
00860 (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
00861 + (a7 + a8*dx)*dy)
00862 + a3*dy + dx*(a2 + a5*dy);
00863 vb += rotation->get_matrix3_row(0);
00864 }
00865 }
00866 }
00867 }
00868 ret->update();
00869 EMDeleteArray(ipcube);
00870 if(rotation) {delete rotation; rotation=0;}
00871 return ret;
00872 }
00873
00874 EMData *StandardProjector::project3d(EMData * image) const
00875 {
00876 Transform* t3d = params["transform"];
00877 if ( t3d == NULL ) throw NullPointerException("The transform object containing the angles(required for projection), was not specified");
00878
00879 if ( image->get_ndim() == 3 )
00880 {
00881
00882
00883
00884
00885 int nx = image->get_xsize();
00886 int ny = image->get_ysize();
00887 int nz = image->get_zsize();
00888
00889
00890 Transform r = t3d->inverse();
00891 int xy = nx * ny;
00892
00893 EMData *proj = new EMData();
00894 proj->set_size(nx, ny, 1);
00895
00896 Vec3i offset(nx/2,ny/2,nz/2);
00897
00898 float *sdata = image->get_data();
00899 float *ddata = proj->get_data();
00900 for (int k = -nz / 2; k < nz - nz / 2; k++) {
00901 int l = 0;
00902 for (int j = -ny / 2; j < ny - ny / 2; j++) {
00903 ddata[l]=0;
00904 for (int i = -nx / 2; i < nx - nx / 2; i++,l++) {
00905
00906 Vec3f coord(i,j,k);
00907 Vec3f soln = r*coord;
00908 soln += offset;
00909
00913
00914
00915 float x2 = soln[0];
00916 float y2 = soln[1];
00917 float z2 = soln[2];
00918
00919 float x = (float)Util::fast_floor(x2);
00920 float y = (float)Util::fast_floor(y2);
00921 float z = (float)Util::fast_floor(z2);
00922
00923 float t = x2 - x;
00924 float u = y2 - y;
00925 float v = z2 - z;
00926
00927 size_t ii = (size_t) ((size_t)x + (size_t)y * nx + (size_t)z * xy);
00928
00929 if (x2 < 0 || y2 < 0 || z2 < 0 ) continue;
00930 if (x2 > (nx-1) || y2 > (ny-1) || z2 > (nz-1) ) continue;
00931
00932 if (x2 < (nx - 1) && y2 < (ny - 1) && z2 < (nz - 1)) {
00933 ddata[l] +=
00934 Util::trilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx],
00935 sdata[ii + nx + 1], sdata[ii + xy], sdata[ii + xy + 1], sdata[ii + xy + nx],
00936 sdata[ii + xy + nx + 1], t, u, v);
00937 }
00938 else if ( x2 == (nx - 1) && y2 == (ny - 1) && z2 == (nz - 1) ) {
00939 ddata[l] += sdata[ii];
00940 }
00941 else if ( x2 == (nx - 1) && y2 == (ny - 1) ) {
00942 ddata[l] += Util::linear_interpolate(sdata[ii], sdata[ii + xy],v);
00943 }
00944 else if ( x2 == (nx - 1) && z2 == (nz - 1) ) {
00945 ddata[l] += Util::linear_interpolate(sdata[ii], sdata[ii + nx],u);
00946 }
00947 else if ( y2 == (ny - 1) && z2 == (nz - 1) ) {
00948 ddata[l] += Util::linear_interpolate(sdata[ii], sdata[ii + 1],t);
00949 }
00950 else if ( x2 == (nx - 1) ) {
00951 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + nx], sdata[ii + xy], sdata[ii + xy + nx],u,v);
00952 }
00953 else if ( y2 == (ny - 1) ) {
00954 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + xy], sdata[ii + xy + 1],t,v);
00955 }
00956 else if ( z2 == (nz - 1) ) {
00957 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx], sdata[ii + nx + 1],t,u);
00958 }
00959 }
00960 }
00961 }
00962 proj->update();
00963 proj->set_attr("xform.projection",t3d);
00964 proj->set_attr("apix_x",(float)image->get_attr("apix_x"));
00965 proj->set_attr("apix_y",(float)image->get_attr("apix_y"));
00966 proj->set_attr("apix_z",(float)image->get_attr("apix_z"));
00967
00968 if(t3d) {delete t3d; t3d=0;}
00969 return proj;
00970 }
00971 else if ( image->get_ndim() == 2 ) {
00972
00973 Transform r = t3d->inverse();
00974
00975 int nx = image->get_xsize();
00976 int ny = image->get_ysize();
00977
00978 EMData *proj = new EMData();
00979 proj->set_size(nx, 1, 1);
00980 proj->to_zero();
00981
00982 float *sdata = image->get_data();
00983 float *ddata = proj->get_data();
00984
00985 Vec2f offset(nx/2,ny/2);
00986 for (int j = -ny / 2; j < ny - ny / 2; j++) {
00987 int l = 0;
00988 for (int i = -nx / 2; i < nx - nx / 2; i++,l++) {
00989
00990 Vec2f coord(i,j);
00991 Vec2f soln = r*coord;
00992 soln += offset;
00993
00994 float x2 = soln[0];
00995 float y2 = soln[1];
00996
00997 float x = (float)Util::fast_floor(x2);
00998 float y = (float)Util::fast_floor(y2);
00999
01000 int ii = (int) (x + y * nx);
01001 float u = x2 - x;
01002 float v = y2 - y;
01003
01004 if (x2 < 0 || y2 < 0 ) continue;
01005 if (x2 > (nx-1) || y2 > (ny-1) ) continue;
01006
01007 if ( x2 < (nx - 1) && y2 < (ny - 1) ) {
01008 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx],sdata[ii + nx + 1], u, v);
01009 }
01010 else if (x2 == (nx-1) && y2 == (ny-1) ) {
01011 ddata[l] += sdata[ii];
01012 }
01013 else if (x2 == (nx-1) ) {
01014 ddata[l] += Util::linear_interpolate(sdata[ii],sdata[ii + nx], v);
01015 }
01016 else if (y2 == (ny-1) ) {
01017 ddata[l] += Util::linear_interpolate(sdata[ii],sdata[ii + 1], u);
01018 }
01019 }
01020 }
01021 proj->set_attr("xform.projection",t3d);
01022 proj->update();
01023 if(t3d) {delete t3d; t3d=0;}
01024 return proj;
01025 }
01026 else throw ImageDimensionException("Standard projection works only for 2D and 3D images");
01027 }
01028
01029 #ifdef EMAN2_USING_CUDA
01030 #include "cuda/cuda_projector.h"
01031 #include "cuda/cuda_util.h"
01032 EMData *CudaStandardProjector::project3d(EMData * image) const
01033 {
01034
01035 device_init();
01036 Transform* t3d = params["transform"];
01037 if ( t3d == NULL ) throw NullPointerException("The transform object containing the angles(required for projection), was not specified");
01038 float * m = new float[12];
01039 t3d->copy_matrix_into_array(m);
01040 image->bind_cuda_texture();
01041 EMData* e = new EMData(image->get_xsize(),image->get_ysize());
01042 EMDataForCuda tmp = e->get_data_struct_for_cuda();
01043 standard_project(m,&tmp);
01044 delete [] m;
01045 e->set_attr("xform.projection",t3d);
01046 if(t3d) {delete t3d; t3d=0;}
01047 e->gpu_update();
01048 return e;
01049 }
01050 #endif // EMAN2_USING_CUDA
01051
01052
01053
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 EMData *FourierGriddingProjector::project3d(EMData * image) const
01145 {
01146 if (!image) {
01147 return 0;
01148 }
01149 if (3 != image->get_ndim())
01150 throw ImageDimensionException(
01151 "FourierGriddingProjector needs a 3-D volume");
01152 if (image->is_complex())
01153 throw ImageFormatException(
01154 "FourierGriddingProjector requires a real volume");
01155 const int npad = params.has_key("npad") ? int(params["npad"]) : 2;
01156 const int nx = image->get_xsize();
01157 const int ny = image->get_ysize();
01158 const int nz = image->get_zsize();
01159 if (nx != ny || nx != nz)
01160 throw ImageDimensionException(
01161 "FourierGriddingProjector requires nx==ny==nz");
01162 const int m = Util::get_min(nx,ny,nz);
01163 const int n = m*npad;
01164
01165 int K = params["kb_K"];
01166 if ( K == 0 ) K = 6;
01167 float alpha = params["kb_alpha"];
01168 if ( alpha == 0 ) alpha = 1.25;
01169 Util::KaiserBessel kb(alpha, K, (float)(m/2), K/(2.0f*n), n);
01170
01171
01172 EMData* tmpImage = image->copy();
01173 tmpImage->divkbsinh(kb);
01174
01175
01176
01177 EMData* imgft = tmpImage->norm_pad(false, npad);
01178 imgft->do_fft_inplace();
01179 imgft->center_origin_fft();
01180 imgft->fft_shuffle();
01181 delete tmpImage;
01182
01183
01184 int nangles = 0;
01185 vector<float> anglelist;
01186
01187 if (params.has_key("anglelist")) {
01188 anglelist = params["anglelist"];
01189 nangles = anglelist.size() / 3;
01190 } else {
01191
01192
01193
01194
01195 Transform* t3d = params["transform"];
01196 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
01197 Dict p = t3d->get_rotation("spider");
01198
01199 string angletype = "SPIDER";
01200 float phi = p["phi"];
01201 float theta = p["theta"];
01202 float psi = p["psi"];
01203 anglelist.push_back(phi);
01204 anglelist.push_back(theta);
01205 anglelist.push_back(psi);
01206 nangles = 1;
01207 if(t3d) {delete t3d; t3d=0;}
01208 }
01209
01210
01211
01212
01213 EMData* ret = new EMData();
01214 ret->set_size(nx, ny, nangles);
01215 ret->to_zero();
01216
01217 for (int ia = 0; ia < nangles; ia++) {
01218 int indx = 3*ia;
01219 Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]);
01220 Transform tf(d);
01221 EMData* proj = imgft->extract_plane(tf, kb);
01222 if (proj->is_shuffled()) proj->fft_shuffle();
01223 proj->center_origin_fft();
01224 proj->do_ift_inplace();
01225 EMData* winproj = proj->window_center(m);
01226 delete proj;
01227 for (int iy=0; iy < ny; iy++)
01228 for (int ix=0; ix < nx; ix++)
01229 (*ret)(ix,iy,ia) = (*winproj)(ix,iy);
01230 delete winproj;
01231 }
01232 delete imgft;
01233
01234 if (!params.has_key("anglelist")) {
01235 Transform* t3d = params["transform"];
01236 ret->set_attr("xform.projection",t3d);
01237 if(t3d) {delete t3d; t3d=0;}
01238 }
01239 ret->update();
01240 return ret;
01241 }
01242
01243
01244 int ChaoProjector::getnnz(Vec3i volsize, int ri, Vec3i origin, int *nrays, int *nnz) const
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258 {
01259 int ix, iy, iz, rs, r2, xs, ys, zs, xx, yy, zz;
01260 int ftm=0, status = 0;
01261
01262 r2 = ri*ri;
01263 *nnz = 0;
01264 *nrays = 0;
01265 int nx = (int)volsize[0];
01266 int ny = (int)volsize[1];
01267 int nz = (int)volsize[2];
01268
01269 int xcent = (int)origin[0];
01270 int ycent = (int)origin[1];
01271 int zcent = (int)origin[2];
01272
01273
01274 for (ix = 1; ix <=nx; ix++) {
01275 xs = ix-xcent;
01276 xx = xs*xs;
01277 for (iy = 1; iy <= ny; iy++) {
01278 ys = iy-ycent;
01279 yy = ys*ys;
01280 ftm = 1;
01281 for (iz = 1; iz <= nz; iz++) {
01282 zs = iz-zcent;
01283 zz = zs*zs;
01284 rs = xx + yy + zz;
01285 if (rs <= r2) {
01286 (*nnz)++;
01287 if (ftm) {
01288 (*nrays)++;
01289 ftm = 0;
01290 }
01291 }
01292 }
01293 }
01294 }
01295 return status;
01296 }
01297
01298 #define cube(i,j,k) cube[ ((k-1)*ny + j-1)*nx + i-1 ]
01299 #define sphere(i) sphere[(i)-1]
01300 #define cord(i,j) cord[((j)-1)*3 + (i) -1]
01301 #define ptrs(i) ptrs[(i)-1]
01302 #define dm(i) dm[(i)-1]
01303
01304 int ChaoProjector:: cb2sph(float *cube, Vec3i volsize, int ri, Vec3i origin,
01305 int nnz0, int *ptrs, int *cord, float *sphere) const
01306 {
01307 int xs, ys, zs, xx, yy, zz, rs, r2;
01308 int ix, iy, iz, jnz, nnz, nrays;
01309 int ftm = 0, status = 0;
01310
01311 int xcent = (int)origin[0];
01312 int ycent = (int)origin[1];
01313 int zcent = (int)origin[2];
01314
01315 int nx = (int)volsize[0];
01316 int ny = (int)volsize[1];
01317 int nz = (int)volsize[2];
01318
01319 r2 = ri*ri;
01320 nnz = 0;
01321 nrays = 0;
01322 ptrs(1) = 1;
01323
01324 for (ix = 1; ix <= nx; ix++) {
01325 xs = ix-xcent;
01326 xx = xs*xs;
01327 for ( iy = 1; iy <= ny; iy++ ) {
01328 ys = iy-ycent;
01329 yy = ys*ys;
01330 jnz = 0;
01331
01332 ftm = 1;
01333
01334 for (iz = 1; iz <= nz; iz++) {
01335 zs = iz-zcent;
01336 zz = zs*zs;
01337 rs = xx + yy + zz;
01338 if (rs <= r2) {
01339 jnz++;
01340 nnz++;
01341 sphere(nnz) = cube(iz, iy, ix);
01342
01343
01344 if (ftm) {
01345 nrays++;
01346 cord(1,nrays) = iz;
01347 cord(2,nrays) = iy;
01348 cord(3,nrays) = ix;
01349 ftm = 0;
01350 }
01351 }
01352 }
01353 if (jnz > 0) {
01354 ptrs(nrays+1) = ptrs(nrays) + jnz;
01355 }
01356 }
01357 }
01358 if (nnz != nnz0) status = -1;
01359 return status;
01360 }
01361
01362
01363 int ChaoProjector::sph2cb(float *sphere, Vec3i volsize, int nrays, int ri,
01364 int nnz0, int *ptrs, int *cord, float *cube) const
01365 {
01366 int status=0;
01367 int r2, i, j, ix, iy, iz, nnz;
01368
01369 int nx = (int)volsize[0];
01370 int ny = (int)volsize[1];
01371
01372
01373 r2 = ri*ri;
01374 nnz = 0;
01375 ptrs(1) = 1;
01376
01377
01378
01379
01380 nnz = 0;
01381 for (j = 1; j <= nrays; j++) {
01382 iz = cord(1,j);
01383 iy = cord(2,j);
01384 ix = cord(3,j);
01385 for (i = ptrs(j); i<=ptrs(j+1)-1; i++, iz++) {
01386 nnz++;
01387 cube(iz,iy,ix) = sphere(nnz);
01388 }
01389 }
01390 if (nnz != nnz0) status = -1;
01391 return status;
01392 }
01393
01394 #define x(i) x[(i)-1]
01395 #define y(i,j) y[(j-1)*nx + i - 1]
01396
01397
01398 int ChaoProjector::fwdpj3(Vec3i volsize, int nrays, int , float *dm,
01399 Vec3i origin, int ri, int *ptrs, int *cord,
01400 float *x, float *y) const
01401 {
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418 int iqx, iqy, i, j, xc, yc, zc;
01419 float ct, dipx, dipy, dipx1m, dipy1m, xb, yb, dm1, dm4;
01420 int status = 0;
01421
01422 int xcent = origin[0];
01423 int ycent = origin[1];
01424 int zcent = origin[2];
01425
01426 int nx = volsize[0];
01427
01428 dm1 = dm(1);
01429 dm4 = dm(4);
01430
01431 if ( nx > 2*ri ) {
01432 for (i = 1; i <= nrays; i++) {
01433 zc = cord(1,i)-zcent;
01434 yc = cord(2,i)-ycent;
01435 xc = cord(3,i)-xcent;
01436
01437 xb = zc*dm(1)+yc*dm(2)+xc*dm(3) + xcent;
01438 yb = zc*dm(4)+yc*dm(5)+xc*dm(6) + ycent;
01439
01440 for (j = ptrs(i); j< ptrs(i+1); j++) {
01441 iqx = ifix(xb);
01442 iqy = ifix(yb);
01443
01444 ct = x(j);
01445 dipx = xb - (float)(iqx);
01446 dipy = (yb - (float)(iqy)) * ct;
01447
01448 dipy1m = ct - dipy;
01449 dipx1m = 1.0f - dipx;
01450
01451 y(iqx ,iqy) = y(iqx ,iqy) + dipx1m*dipy1m;
01452 y(iqx+1,iqy) = y(iqx+1,iqy) + dipx*dipy1m;
01453 y(iqx+1,iqy+1) = y(iqx+1,iqy+1) + dipx*dipy;
01454 y(iqx ,iqy+1) = y(iqx ,iqy+1) + dipx1m*dipy;
01455
01456 xb += dm1;
01457 yb += dm4;
01458 }
01459 }
01460 }
01461 else {
01462 fprintf(stderr, " nx must be greater than 2*ri\n");
01463 exit(1);
01464 }
01465 return status;
01466 }
01467 #undef x
01468 #undef y
01469
01470 #define y(i) y[(i)-1]
01471 #define x(i,j) x[((j)-1)*nx + (i) - 1]
01472
01473
01474 int ChaoProjector::bckpj3(Vec3i volsize, int nrays, int , float *dm,
01475 Vec3i origin, int ri, int *ptrs, int *cord,
01476 float *x, float *y) const
01477 {
01478 int i, j, iqx,iqy, xc, yc, zc;
01479 float xb, yb, dx, dy, dx1m, dy1m, dxdy;
01480 int status = 0;
01481
01482 int xcent = origin[0];
01483 int ycent = origin[1];
01484 int zcent = origin[2];
01485
01486 int nx = volsize[0];
01487
01488 if ( nx > 2*ri) {
01489 for (i = 1; i <= nrays; i++) {
01490 zc = cord(1,i) - zcent;
01491 yc = cord(2,i) - ycent;
01492 xc = cord(3,i) - xcent;
01493
01494 xb = zc*dm(1)+yc*dm(2)+xc*dm(3) + xcent;
01495 yb = zc*dm(4)+yc*dm(5)+xc*dm(6) + ycent;
01496
01497 for (j = ptrs(i); j <ptrs(i+1); j++) {
01498 iqx = ifix((float)(xb));
01499 iqy = ifix((float)(yb));
01500
01501 dx = xb - (float)(iqx);
01502 dy = yb - (float)(iqy);
01503 dx1m = 1.0f - dx;
01504 dy1m = 1.0f - dy;
01505 dxdy = dx*dy;
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522
01523
01524 y(j) += x(iqx,iqy)
01525 + dx*(-x(iqx,iqy)+x(iqx+1,iqy))
01526 + dy*(-x(iqx,iqy)+x(iqx,iqy+1))
01527 + dxdy*( x(iqx,iqy) - x(iqx,iqy+1)
01528 -x(iqx+1,iqy) + x(iqx+1,iqy+1) );
01529
01530 xb += dm(1);
01531 yb += dm(4);
01532 }
01533 }
01534 }
01535 else {
01536 fprintf(stderr, "bckpj3: nx must be greater than 2*ri\n");
01537 }
01538
01539 return status;
01540 }
01541
01542 #undef x
01543 #undef y
01544 #undef dm
01545
01546
01547 int ChaoProjector::ifix(float a) const
01548 {
01549 int ia;
01550
01551 if (a>=0) {
01552 ia = (int)floor(a);
01553 }
01554 else {
01555 ia = (int)ceil(a);
01556 }
01557 return ia;
01558 }
01559
01560 #define dm(i,j) dm[((j)-1)*9 + (i) -1]
01561 #define anglelist(i,j) anglelist[((j)-1)*3 + (i) - 1]
01562
01563
01564 void ChaoProjector::setdm(vector<float> anglelist, string const , float *dm) const
01565 {
01566
01567 float psi, theta, phi;
01568 double cthe, sthe, cpsi, spsi, cphi, sphi;
01569 int j;
01570
01571 int nangles = anglelist.size() / 3;
01572
01573
01574 for (j = 1; j <= nangles; j++) {
01575 phi = static_cast<float>(anglelist(1,j)*dgr_to_rad);
01576 theta = static_cast<float>(anglelist(2,j)*dgr_to_rad);
01577 psi = static_cast<float>(anglelist(3,j)*dgr_to_rad);
01578
01579
01580 cthe = cos(theta);
01581 sthe = sin(theta);
01582 cpsi = cos(psi);
01583 spsi = sin(psi);
01584 cphi = cos(phi);
01585 sphi = sin(phi);
01586
01587 dm(1,j)=static_cast<float>(cphi*cthe*cpsi-sphi*spsi);
01588 dm(2,j)=static_cast<float>(sphi*cthe*cpsi+cphi*spsi);
01589 dm(3,j)=static_cast<float>(-sthe*cpsi);
01590 dm(4,j)=static_cast<float>(-cphi*cthe*spsi-sphi*cpsi);
01591 dm(5,j)=static_cast<float>(-sphi*cthe*spsi+cphi*cpsi);
01592 dm(6,j)=static_cast<float>(sthe*spsi);
01593 dm(7,j)=static_cast<float>(sthe*cphi);
01594 dm(8,j)=static_cast<float>(sthe*sphi);
01595 dm(9,j)=static_cast<float>(cthe);
01596 }
01597 }
01598 #undef anglelist
01599
01600 #define images(i,j,k) images[ ((k-1)*nyvol + j-1)*nxvol + i-1 ]
01601
01602 EMData *ChaoProjector::project3d(EMData * vol) const
01603 {
01604
01605 int nrays, nnz, status, j;
01606 float *dm;
01607 int *ptrs, *cord;
01608 float *sphere, *images;
01609
01610 int nxvol = vol->get_xsize();
01611 int nyvol = vol->get_ysize();
01612 int nzvol = vol->get_zsize();
01613 Vec3i volsize(nxvol,nyvol,nzvol);
01614
01615 int dim = Util::get_min(nxvol,nyvol,nzvol);
01616 if (nzvol == 1) {
01617 LOGERR("The ChaoProjector needs a volume!");
01618 return 0;
01619 }
01620 Vec3i origin(0,0,0);
01621
01622
01623 if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
01624 else {origin[0] = nxvol/2+1;}
01625 if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
01626 else {origin[1] = nyvol/2+1;}
01627 if (params.has_key("origin_z")) {origin[2] = params["origin_z"];}
01628 else {origin[2] = nzvol/2+1;}
01629
01630 int ri;
01631 if (params.has_key("radius")) {ri = params["radius"];}
01632 else {ri = dim/2 - 1;}
01633
01634
01635 float *cube = vol->get_data();
01636
01637
01638
01639 status = getnnz(volsize, ri, origin, &nrays, &nnz);
01640
01641
01642
01643 sphere = new float[nnz];
01644 ptrs = new int[nrays+1];
01645 cord = new int[3*nrays];
01646 if (sphere == NULL || ptrs == NULL || cord == NULL) {
01647 fprintf(stderr,"ChaoProjector::project3d, failed to allocate!\n");
01648 exit(1);
01649 }
01650 for (int i = 0; i<nnz; i++) sphere[i] = 0.0;
01651 for (int i = 0; i<nrays+1; i++) ptrs[i] = 0;
01652 for (int i = 0; i<3*nrays; i++) cord[i] = 0;
01653
01654 status = cb2sph(cube, volsize, ri, origin, nnz, ptrs, cord, sphere);
01655
01656
01657 int nangles = 0;
01658 vector<float> anglelist;
01659 string angletype = "SPIDER";
01660
01661 if (params.has_key("anglelist")) {
01662 anglelist = params["anglelist"];
01663 nangles = anglelist.size() / 3;
01664 } else {
01665 Transform* t3d = params["transform"];
01666 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
01667
01668
01669
01670
01671 Dict p = t3d->get_rotation("spider");
01672 if(t3d) {delete t3d; t3d=0;}
01673
01674 float phi = p["phi"];
01675 float theta = p["theta"];
01676 float psi = p["psi"];
01677 anglelist.push_back(phi);
01678 anglelist.push_back(theta);
01679 anglelist.push_back(psi);
01680 nangles = 1;
01681 }
01682
01683
01684 dm = new float[nangles*9];
01685 setdm(anglelist, angletype, dm);
01686
01687
01688 EMData *ret = new EMData();
01689 ret->set_size(nxvol, nyvol, nangles);
01690 ret->set_complex(false);
01691 ret->set_ri(true);
01692
01693 images = ret->get_data();
01694
01695 for (j = 1; j <= nangles; j++) {
01696 status = fwdpj3(volsize, nrays, nnz , &dm(1,j), origin, ri,
01697 ptrs , cord, sphere, &images(1,1,j));
01698
01699 }
01700
01701
01702 EMDeleteArray(dm);
01703 EMDeleteArray(ptrs);
01704 EMDeleteArray(cord);
01705 EMDeleteArray(sphere);
01706
01707 if (!params.has_key("anglelist")) {
01708 Transform* t3d = params["transform"];
01709 ret->set_attr("xform.projection",t3d);
01710 if(t3d) {delete t3d; t3d=0;}
01711 }
01712 ret->update();
01713 return ret;
01714 }
01715
01716
01717 #undef images
01718
01719 #define images(i,j,k) images[ ((k)-1)*nximg*nyimg + ((j)-1)*nximg + (i)-1 ]
01720
01721 EMData *ChaoProjector::backproject3d(EMData * imagestack) const
01722 {
01723 int nrays, nnz, status, j;
01724 float *dm;
01725 int *ptrs, *cord;
01726 float *sphere, *images, *cube;
01727
01728 int nximg = imagestack->get_xsize();
01729 int nyimg = imagestack->get_ysize();
01730 int nslices = imagestack->get_zsize();
01731
01732 int dim = Util::get_min(nximg,nyimg);
01733 Vec3i volsize(nximg,nyimg,dim);
01734
01735 Vec3i origin(0,0,0);
01736
01737
01738 if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
01739 else {origin[0] = nximg/2+1;}
01740 if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
01741 else {origin[1] = nyimg/2+1;}
01742 if (params.has_key("origin_z")) {origin[1] = params["origin_z"];}
01743 else {origin[2] = dim/2+1;}
01744
01745 int ri;
01746 if (params.has_key("radius")) {ri = params["radius"];}
01747 else {ri = dim/2 - 1;}
01748
01749
01750 images = imagestack->get_data();
01751
01752
01753
01754 status = getnnz(volsize, ri, origin, &nrays, &nnz);
01755
01756
01757
01758 sphere = new float[nnz];
01759 ptrs = new int[nrays+1];
01760 cord = new int[3*nrays];
01761 if (sphere == NULL || ptrs == NULL || cord == NULL) {
01762 fprintf(stderr,"ChaoProjector::backproject3d, failed to allocate!\n");
01763 exit(1);
01764 }
01765 for (int i = 0; i<nnz; i++) sphere[i] = 0.0;
01766 for (int i = 0; i<nrays+1; i++) ptrs[i] = 0;
01767 for (int i = 0; i<3*nrays; i++) cord[i] = 0;
01768
01769 int nangles = 0;
01770 vector<float> anglelist;
01771 string angletype = "SPIDER";
01772
01773 if (params.has_key("anglelist")) {
01774 anglelist = params["anglelist"];
01775 nangles = anglelist.size() / 3;
01776 } else {
01777 Transform* t3d = params["transform"];
01778 if ( t3d == NULL ) throw NullPointerException("The transform3d object (required for projection), was not specified");
01779
01780
01781
01782
01783
01784 Dict p = t3d->get_rotation("spider");
01785 if(t3d) {delete t3d; t3d=0;}
01786
01787 float phi = p["phi"];
01788 float theta = p["theta"];
01789 float psi = p["psi"];
01790 anglelist.push_back(phi);
01791 anglelist.push_back(theta);
01792 anglelist.push_back(psi);
01793 nangles = 1;
01794 }
01795
01796
01797
01798 if (nslices != nangles) {
01799 LOGERR("the number of images does not match the number of angles");
01800 return 0;
01801 }
01802
01803 dm = new float[nangles*9];
01804 setdm(anglelist, angletype, dm);
01805
01806
01807 EMData *ret = new EMData();
01808 ret->set_size(nximg, nyimg, dim);
01809 ret->set_complex(false);
01810 ret->set_ri(true);
01811 ret->to_zero();
01812
01813 cube = ret->get_data();
01814
01815 status = cb2sph(cube, volsize, ri, origin, nnz, ptrs, cord, sphere);
01816
01817
01818 for (j = 1; j <= nangles; j++) {
01819 status = bckpj3(volsize, nrays, nnz, &dm(1,j), origin, ri,
01820 ptrs , cord , &images(1,1,j), sphere);
01821
01822 }
01823
01824 status = sph2cb(sphere, volsize, nrays, ri, nnz, ptrs, cord, cube);
01825
01826
01827
01828 EMDeleteArray(dm);
01829 EMDeleteArray(ptrs);
01830 EMDeleteArray(cord);
01831 EMDeleteArray(sphere);
01832
01833 ret->update();
01834 return ret;
01835 }
01836
01837 #undef images
01838 #undef cube
01839 #undef sphere
01840 #undef cord
01841 #undef ptrs
01842 #undef dm
01843
01844 EMData *GaussFFTProjector::backproject3d(EMData * ) const
01845 {
01846
01847 EMData *ret = new EMData();
01848 return ret;
01849 }
01850
01851 #define images(i,j,k) images[ (k)*nx*ny + ((j)-1)*nx + (i)-1 ]
01852
01853
01854
01855
01856
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 EMData *PawelProjector::backproject3d(EMData * imagestack) const
01968 {
01969
01970 float *images;
01971
01972 if (!imagestack) {
01973 return 0;
01974 }
01975 int ri;
01976 int nx = imagestack->get_xsize();
01977 int ny = imagestack->get_ysize();
01978
01979 int dim = Util::get_min(nx,ny);
01980 images = imagestack->get_data();
01981
01982 Vec3i origin(0,0,0);
01983
01984
01985 if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
01986 else {origin[0] = nx/2;}
01987 if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
01988 else {origin[1] = ny/2;}
01989 if (params.has_key("origin_z")) {origin[1] = params["origin_z"];}
01990 else {origin[2] = dim/2;}
01991
01992 if (params.has_key("radius")) {ri = params["radius"];}
01993 else {ri = dim/2 - 1;}
01994
01995
01996 int nn = -1;
01997 prepcubes(nx, ny, dim, ri, origin, nn);
01998
01999
02000 IPCube* ipcube = new IPCube[nn+1];
02001 prepcubes(nx, ny, dim, ri, origin, nn, ipcube);
02002
02003 int nangles = 0;
02004 vector<float> anglelist;
02005
02006 if (params.has_key("anglelist")) {
02007 anglelist = params["anglelist"];
02008 nangles = anglelist.size() / 3;
02009 } else {
02010 Transform* t3d = params["transform"];
02011 if ( t3d == NULL ) throw NullPointerException("The transform3d object (required for projection), was not specified");
02012
02013
02014
02015
02016 Dict p = t3d->get_rotation("spider");
02017 if(t3d) {delete t3d; t3d=0;}
02018
02019 string angletype = "SPIDER";
02020 float phi = p["phi"];
02021 float theta = p["theta"];
02022 float psi = p["psi"];
02023 anglelist.push_back(phi);
02024 anglelist.push_back(theta);
02025 anglelist.push_back(psi);
02026 nangles = 1;
02027 }
02028
02029
02030
02031
02032 EMData* ret = new EMData();
02033 ret->set_size(nx, ny, dim);
02034 ret->to_zero();
02035
02036
02037 for (int ia = 0; ia < nangles; ia++) {
02038 int indx = 3*ia;
02039 Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]);
02040 Transform rotation(d);
02041 float dm1 = rotation.at(0,0);
02042 float dm4 = rotation.at(1,0);
02043
02044 if (2*(ri+1)+1 > dim) {
02045
02046 LOGERR("backproject3d, pawel, 2*(ri+1)+1 > dim\n");
02047 return 0;
02048 } else {
02049
02050 for (int i = 0 ; i <= nn; i++) {
02051 int iox = (int)ipcube[i].loc[0]+origin[0];
02052 int ioy = (int)ipcube[i].loc[1]+origin[1];
02053 int ioz = (int)ipcube[i].loc[2]+origin[2];
02054
02055 Vec3f vb = rotation*ipcube[i].loc + origin;
02056 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
02057 float xbb = (j-ipcube[i].start)*dm1 + vb[0];
02058 int iqx = (int)floor(xbb);
02059
02060 float ybb = (j-ipcube[i].start)*dm4 + vb[1];
02061 int iqy = (int)floor(ybb);
02062
02063 float dipx = xbb - iqx;
02064 float dipy = ybb - iqy;
02065
02066 (*ret)(iox,ioy,ioz) += images(iqx,iqy,ia)
02067 + dipy*(images(iqx,iqy+1,ia)-images(iqx,iqy,ia))
02068 + dipx*(images(iqx+1,iqy,ia)-images(iqx,iqy,ia)
02069 + dipy*(images(iqx+1,iqy+1,ia)-images(iqx+1,iqy,ia)
02070 - images(iqx,iqy+1,ia)+images(iqx,iqy,ia)));
02071 iox++;
02072 }
02073 }
02074 }
02075 }
02076
02077 ret->update();
02078 EMDeleteArray(ipcube);
02079 return ret;
02080 }
02081 #undef images
02082
02083 EMData *StandardProjector::backproject3d(EMData * ) const
02084 {
02085
02086 EMData *ret = new EMData();
02087 return ret;
02088 }
02089
02090 #ifdef EMAN2_USING_CUDA
02091 EMData *CudaStandardProjector::backproject3d(EMData * ) const
02092 {
02093
02094 EMData *ret = new EMData();
02095 return ret;
02096 }
02097 #endif //EMAN2_USING_CUDA
02098
02099 EMData *FourierGriddingProjector::backproject3d(EMData * ) const
02100 {
02101
02102 EMData *ret = new EMData();
02103 return ret;
02104 }
02105
02106
02107
02108 void EMAN::dump_projectors()
02109 {
02110 dump_factory < Projector > ();
02111 }
02112
02113 map<string, vector<string> > EMAN::dump_projectors_list()
02114 {
02115 return dump_factory_list < Projector > ();
02116 }
02117
02118