00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include "projector.h"
00037 #include "emdata.h"
00038 #include "interp.h"
00039 #include "emutil.h"
00040 #include "plugins/projector_template.h"
00041
00042 #ifdef WIN32
00043 #define M_PI 3.14159265358979323846f
00044 #endif //WIN32
00045
00046 #ifdef EMAN2_USING_CUDA
00047 #include "cuda/cuda_util.h"
00048 #include "cuda/cuda_projector.h"
00049 #endif
00050
00051 using namespace std;
00052 using namespace EMAN;
00053
00054 const string GaussFFTProjector::NAME = "gauss_fft";
00055 const string FourierGriddingProjector::NAME = "fourier_gridding";
00056 const string PawelProjector::NAME = "pawel";
00057 const string StandardProjector::NAME = "standard";
00058 const string ChaoProjector::NAME = "chao";
00059
00060 template <> Factory < Projector >::Factory()
00061 {
00062 force_add<GaussFFTProjector>();
00063 force_add<PawelProjector>();
00064 force_add<StandardProjector>();
00065 force_add<FourierGriddingProjector>();
00066 force_add<ChaoProjector>();
00067
00068
00069 }
00070
00071 EMData *GaussFFTProjector::project3d(EMData * image) const
00072 {
00073 Transform* t3d = params["transform"];
00074 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
00075 if ( image->get_ndim() != 3 ) throw ImageDimensionException("Error, the projection volume must be 3D");
00076
00077 EMData *f = image;
00078 if (!image->is_complex()) {
00079 image->process_inplace("xform.phaseorigin.tocorner");
00080 f = image->do_fft();
00081 f->process_inplace("xform.fourierorigin.tocenter");
00082 image->process_inplace("xform.phaseorigin.tocenter");
00083 }
00084
00085 int f_nx = f->get_xsize();
00086 int f_ny = f->get_ysize();
00087 int f_nz = f->get_zsize();
00088
00089 if (!f->is_complex() || f_nz != f_ny || f_nx != f_ny + 2) {
00090 LOGERR("Cannot project this image");
00091 return 0;
00092 }
00093
00094 f->ap2ri();
00095
00096 EMData *tmp = new EMData();
00097 tmp->set_size(f_nx, f_ny, 1);
00098 tmp->set_complex(true);
00099 tmp->set_ri(true);
00100
00101 float *data = tmp->get_data();
00102
00103 Transform r = t3d->get_rotation_transform();
00104 r.invert();
00105 float scale = t3d->get_scale();
00106
00107 int mode = params["mode"];
00108 float gauss_width = 1;
00109 if ( mode == 0 ) mode = 2;
00110 if (mode == 2 ) {
00111 gauss_width = EMConsts::I2G;
00112 }
00113 else if (mode == 3) {
00114 gauss_width = EMConsts::I3G;
00115 }
00116 else if (mode == 4) {
00117 gauss_width = EMConsts::I4G;
00118 }
00119 else if (mode == 6 || mode == 7) {
00120 gauss_width = EMConsts::I5G;
00121 }
00122
00123 for (int y = 0; y < f_ny; y++) {
00124 for (int x = 0; x < f_nx / 2; x++) {
00125 int ii = x * 2 + y * f_nx;
00126 #ifdef _WIN32
00127 if (_hypot(x, y - f_ny / 2) >= f_ny / 2 - 1) {
00128 #else
00129 if (hypot(x, y - f_ny / 2) >= f_ny / 2 - 1) {
00130 #endif //_WIN32
00131 data[ii] = 0;
00132 data[ii + 1] = 0;
00133 }
00134 else {
00135 Vec3f coord(x,(y - f_ny / 2),0);
00136 Vec3f soln = r*coord;
00137 float xx = soln[0];
00138 float yy = soln[1];
00139 float zz = soln[2];
00140
00141 int cc = 1;
00142 if (xx < 0) {
00143 xx = -xx;
00144 yy = -yy;
00145 zz = -zz;
00146 cc = -1;
00147 }
00148
00149 if (scale != 1.0) {
00150 xx *= scale;
00151 yy *= scale;
00152 zz *= scale;
00153 }
00154
00155 yy += f_ny / 2;
00156 zz += f_nz / 2;
00157
00158 if (yy < 0 || xx < 0 || zz < 0) {
00159 data[ii] = 0;
00160 data[ii+1] = 0;
00161 continue;
00162 }
00163
00164 if (interp_ft_3d(mode, f, xx, yy, zz, data + ii, gauss_width)) {
00165 data[ii + 1] *= cc;
00166 } else {
00167 data[ii] = 0;
00168 data[ii+1] = 0;
00169 }
00170 }
00171 }
00172 }
00173
00174 f->update();
00175 tmp->update();
00176
00177 tmp->process_inplace("xform.fourierorigin.tocorner");
00178 EMData *ret = tmp->do_ift();
00179 ret->process_inplace("xform.phaseorigin.tocenter");
00180
00181 ret->translate(t3d->get_trans());
00182
00183 if (t3d->get_mirror() ) ret->process_inplace("xform.flip",Dict("axis","x"));
00184
00185 Dict filter_d;
00186 filter_d["gauss_width"] = gauss_width;
00187 filter_d["ring_width"] = ret->get_xsize() / 2;
00188 ret->process_inplace("math.gausskernelfix", filter_d);
00189
00190 if( tmp )
00191 {
00192 delete tmp;
00193 tmp = 0;
00194 }
00195
00196 ret->set_attr("xform.projection",t3d);
00197 ret->update();
00198
00199 if(t3d) {delete t3d; t3d=0;}
00200
00201 return ret;
00202 }
00203
00204
00205
00206 bool GaussFFTProjector::interp_ft_3d(int mode, EMData * image, float x, float y,
00207 float z, float *data, float gw) const
00208 {
00209 float *rdata = image->get_data();
00210 int nx = image->get_xsize();
00211 int ny = image->get_ysize();
00212 int nz = image->get_zsize();
00213
00214 if ( mode == 0 ) mode = 2;
00215
00216 if (mode == 1) {
00217 int x0 = 2 * (int) floor(x + 0.5f);
00218 int y0 = (int) floor(y + 0.5f);
00219 int z0 = (int) floor(z + 0.5f);
00220
00221 data[0] = rdata[x0 + y0 * nx + z0 * nx * ny];
00222 data[1] = rdata[x0 + y0 * nx + z0 * nx * ny + 1];
00223 return true;
00224 }
00225 else if (mode == 2) {
00226 int x0 = (int) floor(x);
00227 int y0 = (int) floor(y);
00228 int z0 = (int) floor(z);
00229
00230 float dx = x - x0;
00231 float dy = y - y0;
00232 float dz = z - z0;
00233
00234 if (x0 <= nx/2- 2 && y0 <= ny - 1 && z0 <= nz - 1) {
00235 int i = (int) (x0 * 2 + y0 * nx + z0 * nx * ny);
00236
00237 float n = Util::agauss(1, dx, dy, dz, gw) +
00238 Util::agauss(1, 1 - dx, dy, dz, gw) +
00239 Util::agauss(1, dx, 1 - dy, dz, gw) +
00240 Util::agauss(1, 1 - dx, 1 - dy, dz, gw) +
00241 Util::agauss(1, dx, dy, 1 - dz, gw) +
00242 Util::agauss(1, 1 - dx, dy, 1 - dz, gw) +
00243 Util::agauss(1, dx, 1 - dy, 1 - dz, gw) + Util::agauss(1, 1 - dx, 1 - dy, 1 - dz, gw);
00244
00245 data[0] = Util::agauss(rdata[i], dx, dy, dz, gw) +
00246 Util::agauss(rdata[i + 2], 1 - dx, dy, dz, gw) +
00247 Util::agauss(rdata[i + nx], dx, 1 - dy, dz, gw) +
00248 Util::agauss(rdata[i + nx + 2], 1 - dx, 1 - dy, dz, gw) +
00249 Util::agauss(rdata[i + nx * ny], dx, dy, 1 - dz, gw) +
00250 Util::agauss(rdata[i + 2 + nx * ny], 1 - dx, dy, 1 - dz, gw) +
00251 Util::agauss(rdata[i + nx + nx * ny], dx, 1 - dy, 1 - dz, gw) +
00252 Util::agauss(rdata[i + 2 + nx + nx * ny], 1 - dx, 1 - dy, 1 - dz, gw) / n;
00253
00254 i++;
00255
00256 data[1] = Util::agauss(rdata[i], dx, dy, dz, gw) +
00257 Util::agauss(rdata[i + 2], 1 - dx, dy, dz, gw) +
00258 Util::agauss(rdata[i + nx], dx, 1 - dy, dz, gw) +
00259 Util::agauss(rdata[i + nx + 2], 1 - dx, 1 - dy, dz, gw) +
00260 Util::agauss(rdata[i + nx * ny], dx, dy, 1 - dz, gw) +
00261 Util::agauss(rdata[i + 2 + nx * ny], 1 - dx, dy, 1 - dz, gw) +
00262 Util::agauss(rdata[i + nx + nx * ny], dx, 1 - dy, 1 - dz, gw) +
00263 Util::agauss(rdata[i + 2 + nx + nx * ny], 1 - dx, 1 - dy, 1 - dz, gw) / n;
00264 return true;
00265 }
00266 return false;
00267 }
00268 else if (mode == 3) {
00269 int x0 = 2 * (int) floor(x + .5);
00270 int y0 = (int) floor(y + .5);
00271 int z0 = (int) floor(z + .5);
00272
00273 float *supp = image->setup4slice();
00274
00275 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
00276 float n = 0;
00277
00278 if (x0 == 0) {
00279 x0 += 4;
00280 size_t idx;
00281 float r, g;
00282 for (int k = z0 - 1; k <= z0 + 1; k++) {
00283 for (int j = y0 - 1; j <= y0 + 1; j++) {
00284 for (int i = x0 - 2; i <= x0 + 2; i += 2) {
00285 r = Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z);
00286 g = exp(-r / gw);
00287 n += g;
00288 idx = i + j * 12 + (size_t)k * 12 * ny;
00289 data[0] += supp[idx] * g;
00290 data[1] += supp[idx + 1] * g;
00291 }
00292 }
00293 }
00294 }
00295 else {
00296 size_t idx;
00297 float r, g;
00298 for (int k = z0 - 1; k <= z0 + 1; k++) {
00299 for (int j = y0 - 1; j <= y0 + 1; j++) {
00300 for (int i = x0 - 2; i <= x0 + 2; i += 2) {
00301 r = Util::hypot3sq(i / 2.0f - x, j - y, k - z);
00302 g = exp(-r / gw);
00303 n += g;
00304 idx = i + j * nx + (size_t)k * nx * ny;
00305 data[0] += rdata[idx] * g;
00306 data[1] += rdata[idx + 1] * g;
00307 }
00308 }
00309 }
00310 }
00311
00312 data[0] /= n;
00313 data[1] /= n;
00314 return true;
00315 }
00316 return false;
00317 }
00318 else if (mode == 4) {
00319 int x0 = 2 * (int) floor(x);
00320 int y0 = (int) floor(y);
00321 int z0 = (int) floor(z);
00322
00323 float *supp = image->setup4slice();
00324
00325 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
00326 float n = 0;
00327
00328 if (x0 == 0) {
00329 x0 += 4;
00330 size_t idx;
00331 float r, g;
00332 for (int k = z0 - 1; k <= z0 + 2; k++) {
00333 for (int j = y0 - 1; j <= y0 + 2; j++) {
00334 for (int i = x0 - 2; i <= x0 + 4; i += 2) {
00335 r = Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z);
00336 g = exp(-r / gw);
00337 n += g;
00338 idx = i + j * 12 + (size_t)k * 12 * ny;
00339 data[0] += supp[idx] * g;
00340 data[1] += supp[idx + 1] * g;
00341 }
00342 }
00343 }
00344 }
00345 else {
00346 float r, g;
00347 size_t idx;
00348 for (int k = z0 - 1; k <= z0 + 2; k++) {
00349 for (int j = y0 - 1; j <= y0 + 2; j++) {
00350 for (int i = x0 - 2; i <= x0 + 4; i += 2) {
00351 r = Util::hypot3sq(i / 2.0f - x, j - y, k - z);
00352 g = exp(-r / gw);
00353 n += g;
00354 idx = i + j * nx + (size_t)k * nx * ny;
00355 data[0] += rdata[idx] * g;
00356 data[1] += rdata[idx + 1] * g;
00357 }
00358 }
00359 }
00360 }
00361 data[0] /= n;
00362 data[1] /= n;
00363 return true;
00364 }
00365 return false;
00366 }
00367 else if (mode == 5) {
00368 int x0 = 2 * (int) floor(x + .5);
00369 int y0 = (int) floor(y + .5);
00370 int z0 = (int) floor(z + .5);
00371
00372 float *supp = image->setup4slice();
00373 float *gimx = Interp::get_gimx();
00374
00375 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
00376 int mx0 = -(int) floor((x - x0 / 2) * 39.0 + .5) - 78;
00377 int my0 = -(int) floor((y - y0) * 39.0 + .5) - 78;
00378 int mz0 = -(int) floor((z - z0) * 39.0 + .5) - 78;
00379
00380 float n = 0;
00381 int mmz = mz0;
00382 int mmy = my0;
00383 int mmx = mx0;
00384
00385 if (x0 < 4) {
00386 x0 += 4;
00387
00388 size_t idx;
00389 float g;
00390 for (int k = z0 - 2; k <= z0 + 2; k++, mmz += 39) {
00391 for (int j = y0 - 2; j <= y0 + 2; j++, mmy += 39) {
00392 for (int i = x0 - 4; i <= x0 + 4; i += 2, mmx += 39) {
00393 g = gimx[abs(mmx) + abs(mmy) * 100 + abs(mmz) * 10000];
00394 n += g;
00395 idx = i + j * 12 + (size_t)k * 12 * ny;
00396 data[0] += supp[idx] * g;
00397 data[1] += supp[idx + 1] * g;
00398 }
00399 }
00400 }
00401 }
00402 else {
00403 size_t ii;
00404 float g;
00405 for (int k = z0 - 2; k <= z0 + 2; k++, mmz += 39) {
00406 for (int j = y0 - 2; j <= y0 + 2; j++, mmy += 39) {
00407 for (int i = x0 - 4; i <= x0 + 4; i += 2, mmx += 39) {
00408 ii = i + j * nx + (size_t)k * nx * ny;
00409 g = gimx[abs(mmx) + abs(mmy) * 100 + abs(mmz) * 10000];
00410 n += g;
00411 data[0] += rdata[ii] * g;
00412 data[1] += rdata[ii + 1] * g;
00413 }
00414 }
00415 }
00416 }
00417
00418 data[0] /= n;
00419 data[1] /= n;
00420 return true;
00421 }
00422 return false;
00423 }
00424 else if (mode == 6) {
00425
00426 int x0 = 2 * (int) floor(x + .5);
00427 int y0 = (int) floor(y + .5);
00428 int z0 = (int) floor(z + .5);
00429
00430 float *supp = image->setup4slice();
00431
00432 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
00433 float n = 0;
00434
00435 if (x0 < 4) {
00436 x0 += 4;
00437 float r, g;
00438 size_t idx;
00439 for (int k = z0 - 2; k <= z0 + 2; k++) {
00440 for (int j = y0 - 2; j <= y0 + 2; j++) {
00441 for (int i = x0 - 4; i <= x0 + 4; i += 2) {
00442 r = Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z);
00443 g = exp(-r / gw);
00444 n += g;
00445 idx = i + j * 12 + (size_t)k * 12 * ny;
00446 data[0] += supp[idx] * g;
00447 data[1] += supp[idx + 1] * g;
00448 }
00449 }
00450 }
00451 }
00452 else {
00453 float r, g;
00454 size_t idx;
00455 for (int k = z0 - 2; k <= z0 + 2; k++) {
00456 for (int j = y0 - 2; j <= y0 + 2; j++) {
00457 for (int i = x0 - 4; i <= x0 + 4; i += 2) {
00458 r = Util::hypot3sq(i / 2.0f - x, j - y, k - z);
00459 g = exp(-r / gw);
00460 n += g;
00461 idx = i + j * nx + (size_t)k * nx * ny;
00462 data[0] += rdata[idx] * g;
00463 data[1] += rdata[idx + 1] * g;
00464 }
00465
00466 }
00467 }
00468 }
00469
00470 data[0] /= n;
00471 data[1] /= n;
00472
00473 return true;
00474 }
00475 return false;
00476 }
00477 else if (mode == 7) {
00478 int x0 = 2 * (int) floor(x + .5);
00479 int y0 = (int) floor(y + .5);
00480 int z0 = (int) floor(z + .5);
00481
00482 float *supp = image->setup4slice();
00483
00484 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
00485 float n = 0;
00486 if (x0 < 4) {
00487 x0 += 4;
00488 float r, g;
00489 size_t idx;
00490 for (int k = z0 - 2; k <= z0 + 2; k++) {
00491 for (int j = y0 - 2; j <= y0 + 2; j++) {
00492 for (int i = x0 - 4; i <= x0 + 4; i += 2) {
00493 r = sqrt(Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z));
00494 g = Interp::hyperg(r);
00495 n += g;
00496 idx = i + j * 12 + (size_t)k * 12 * ny;
00497 data[0] += supp[idx] * g;
00498 data[1] += supp[idx + 1] * g;
00499 }
00500 }
00501 }
00502 }
00503 else {
00504 float r, g;
00505 size_t idx;
00506 for (int k = z0 - 2; k <= z0 + 2; k++) {
00507 for (int j = y0 - 2; j <= y0 + 2; j++) {
00508 for (int i = x0 - 4; i <= x0 + 4; i += 2) {
00509 r = sqrt(Util::hypot3sq(i / 2.0f - x, j - y, k - z));
00510 g = Interp::hyperg(r);
00511 n += g;
00512 idx = i + j * nx + (size_t)k * nx * ny;
00513 data[0] += rdata[idx] * g;
00514 data[1] += rdata[idx + 1] * g;
00515 }
00516
00517 }
00518 }
00519 }
00520 data[0] /= n;
00521 data[1] /= n;
00522 return true;
00523 }
00524 return false;
00525
00526 }
00527
00528 return false;
00529 }
00530
00531 void PawelProjector::prepcubes(int, int ny, int nz, int ri, Vec3i origin,
00532 int& nn, IPCube* ipcube) const {
00533 const float r = float(ri*ri);
00534 const int ldpx = origin[0];
00535 const int ldpy = origin[1];
00536 const int ldpz = origin[2];
00537 float t;
00538 nn = -1;
00539 for (int i1 = 0; i1 < nz; i1++) {
00540 t = float(i1 - ldpz);
00541 const float xx = t*t;
00542 for (int i2 = 0; i2 < ny; i2++) {
00543 t = float(i2 - ldpy);
00544 const float yy = t*t + xx;
00545 bool first = true;
00546 for (int i3 = 0; i3 < nz; i3++) {
00547 t = float(i3 - ldpx);
00548 const float rc = t*t + yy;
00549 if (first) {
00550
00551 if (rc > r) continue;
00552 first = false;
00553 nn++;
00554 if (ipcube != NULL) {
00555 ipcube[nn].start = i3;
00556 ipcube[nn].end = i3;
00557 ipcube[nn].loc[0] = i3 - ldpx;
00558 ipcube[nn].loc[1] = i2 - ldpy;
00559 ipcube[nn].loc[2] = i1 - ldpz;
00560 }
00561 } else {
00562
00563 if (ipcube != NULL) {
00564 if (rc <= r) ipcube[nn].end = i3;
00565 }
00566 }
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 EMData *PawelProjector::project3d(EMData * image) const
00724 {
00725 if (!image) {
00726 return 0;
00727 }
00728 int ri;
00729 int nx = image->get_xsize();
00730 int ny = image->get_ysize();
00731 int nz = image->get_zsize();
00732 int dim = Util::get_min(nx,ny,nz);
00733 if (nz == 1) {
00734 LOGERR("The PawelProjector needs a volume!");
00735 return 0;
00736 }
00737 Vec3i origin(0,0,0);
00738
00739
00740 if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
00741 else {origin[0] = nx/2;}
00742 if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
00743 else {origin[1] = ny/2;}
00744 if (params.has_key("origin_z")) {origin[2] = params["origin_z"];}
00745 else {origin[2] = nz/2;}
00746
00747 if (params.has_key("radius")) {ri = params["radius"];}
00748 else {ri = dim/2 - 1;}
00749
00750
00751 int nn = -1;
00752 prepcubes(nx, ny, nz, ri, origin, nn);
00753
00754
00755 IPCube* ipcube = new IPCube[nn+1];
00756 prepcubes(nx, ny, nz, ri, origin, nn, ipcube);
00757
00758 Transform* rotation = params["transform"];
00759 int nangles = 0;
00760 vector<float> anglelist;
00761
00762
00763
00764
00765
00766
00767
00768 if ( rotation == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780 nangles = 1;
00781
00782
00783
00784 EMData* ret = new EMData();
00785 ret->set_size(nx, ny, nangles);
00786 ret->to_zero();
00787
00788
00789 for (int ia = 0; ia < nangles; ia++) {
00790
00791
00792
00793 if (2*(ri+1)+1 > dim) {
00794
00795 for (int i = 0 ; i <= nn; i++) {
00796 int k = ipcube[i].loc[1] + origin[1];
00797 Vec3f vb = ipcube[i].loc*(*rotation) + origin;
00798 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
00799
00800 int iox = int(vb[0]);
00801 if ((iox >= 0) && (iox < nx-1)) {
00802 int ioy = int(vb[1]);
00803 if ((ioy >= 0) && (ioy < ny-1)) {
00804 int ioz = int(vb[2]);
00805 if ((ioz >= 0) && (ioz < nz-1)) {
00806
00807 float dx = vb[0] - iox;
00808 float dy = vb[1] - ioy;
00809 float dz = vb[2] - ioz;
00810 float a1 = (*image)(iox,ioy,ioz);
00811 float a2 = (*image)(iox+1,ioy,ioz) - a1;
00812 float a3 = (*image)(iox,ioy+1,ioz) - a1;
00813 float a4 = (*image)(iox,ioy,ioz+1) - a1;
00814 float a5 = -a2 -(*image)(iox,ioy+1,ioz)
00815 + (*image)(iox+1,ioy+1,ioz);
00816 float a61 = -(*image)(iox,ioy,ioz+1)
00817 + (*image)(iox+1,ioy,ioz+1);
00818 float a6 = -a2 + a61;
00819 float a7 = -a3 - (*image)(iox,ioy,ioz+1)
00820 + (*image)(iox,ioy+1,ioz+1);
00821 float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
00822 + (*image)(iox+1,ioy+1,ioz+1);
00823 (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
00824 + (a7 + a8*dx)*dy)
00825 + a3*dy + dx*(a2 + a5*dy);
00826 }
00827 }
00828 }
00829 vb += rotation->get_matrix3_row(0);
00830 }
00831 }
00832
00833 } else {
00834
00835 for (int i = 0 ; i <= nn; i++) {
00836 int k = ipcube[i].loc[1] + origin[1];
00837 Vec3f vb = ipcube[i].loc*(*rotation) + origin;
00838 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
00839 int iox = int(vb[0]);
00840 int ioy = int(vb[1]);
00841 int ioz = int(vb[2]);
00842 float dx = vb[0] - iox;
00843 float dy = vb[1] - ioy;
00844 float dz = vb[2] - ioz;
00845 float a1 = (*image)(iox,ioy,ioz);
00846 float a2 = (*image)(iox+1,ioy,ioz) - a1;
00847 float a3 = (*image)(iox,ioy+1,ioz) - a1;
00848 float a4 = (*image)(iox,ioy,ioz+1) - a1;
00849 float a5 = -a2 -(*image)(iox,ioy+1,ioz)
00850 + (*image)(iox+1,ioy+1,ioz);
00851 float a61 = -(*image)(iox,ioy,ioz+1)
00852 + (*image)(iox+1,ioy,ioz+1);
00853 float a6 = -a2 + a61;
00854 float a7 = -a3 - (*image)(iox,ioy,ioz+1)
00855 + (*image)(iox,ioy+1,ioz+1);
00856 float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
00857 + (*image)(iox+1,ioy+1,ioz+1);
00858 (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
00859 + (a7 + a8*dx)*dy)
00860 + a3*dy + dx*(a2 + a5*dy);
00861 vb += rotation->get_matrix3_row(0);
00862 }
00863 }
00864 }
00865 }
00866 ret->update();
00867 EMDeleteArray(ipcube);
00868 if(rotation) {delete rotation; rotation=0;}
00869 return ret;
00870 }
00871
00872 EMData *StandardProjector::project3d(EMData * image) const
00873 {
00874 Transform* t3d = params["transform"];
00875 if ( t3d == NULL ) throw NullPointerException("The transform object containing the angles(required for projection), was not specified");
00876
00877 if ( image->get_ndim() == 3 )
00878 {
00879
00880 #ifdef EMAN2_USING_CUDA
00881 if(EMData::usecuda == 1) {
00882 if(!image->isrodataongpu()) image->copy_to_cudaro();
00883
00884 Transform* t3d = params["transform"];
00885 if ( t3d == NULL ) throw NullPointerException("The transform object containing the angles(required for projection), was not specified");
00886 float * m = new float[12];
00887 t3d->copy_matrix_into_array(m);
00888 image->bindcudaarrayA(true);
00889
00890 EMData *e = new EMData();
00891 e->set_size_cuda(image->get_xsize(), image->get_ysize(), 1);
00892 e->rw_alloc();
00893 standard_project(m,e->getcudarwdata(), e->get_xsize(), e->get_ysize());
00894 image->unbindcudaarryA();
00895 delete [] m;
00896
00897 e->update();
00898 e->set_attr("xform.projection",t3d);
00899 e->set_attr("apix_x",(float)image->get_attr("apix_x"));
00900 e->set_attr("apix_y",(float)image->get_attr("apix_y"));
00901 e->set_attr("apix_z",(float)image->get_attr("apix_z"));
00902
00903 if(t3d) {delete t3d; t3d=0;}
00904 return e;
00905 }
00906 #endif
00907 int nx = image->get_xsize();
00908 int ny = image->get_ysize();
00909 int nz = image->get_zsize();
00910
00911
00912 Transform r = t3d->inverse();
00913 int xy = nx * ny;
00914
00915 EMData *proj = new EMData();
00916 proj->set_size(nx, ny, 1);
00917
00918 Vec3i offset(nx/2,ny/2,nz/2);
00919
00920 float *sdata = image->get_data();
00921 float *ddata = proj->get_data();
00922 for (int k = -nz / 2; k < nz - nz / 2; k++) {
00923 int l = 0;
00924 for (int j = -ny / 2; j < ny - ny / 2; j++) {
00925 ddata[l]=0;
00926 for (int i = -nx / 2; i < nx - nx / 2; i++,l++) {
00927
00928 Vec3f coord(i,j,k);
00929 Vec3f soln = r*coord;
00930 soln += offset;
00931
00935
00936
00937 float x2 = soln[0];
00938 float y2 = soln[1];
00939 float z2 = soln[2];
00940
00941 float x = (float)Util::fast_floor(x2);
00942 float y = (float)Util::fast_floor(y2);
00943 float z = (float)Util::fast_floor(z2);
00944
00945 float t = x2 - x;
00946 float u = y2 - y;
00947 float v = z2 - z;
00948
00949 size_t ii = (size_t) ((size_t)x + (size_t)y * nx + (size_t)z * xy);
00950
00951 if (x2 < 0 || y2 < 0 || z2 < 0 ) continue;
00952 if (x2 > (nx-1) || y2 > (ny-1) || z2 > (nz-1) ) continue;
00953
00954 if (x2 < (nx - 1) && y2 < (ny - 1) && z2 < (nz - 1)) {
00955 ddata[l] +=
00956 Util::trilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx],
00957 sdata[ii + nx + 1], sdata[ii + xy], sdata[ii + xy + 1], sdata[ii + xy + nx],
00958 sdata[ii + xy + nx + 1], t, u, v);
00959 }
00960 else if ( x2 == (nx - 1) && y2 == (ny - 1) && z2 == (nz - 1) ) {
00961 ddata[l] += sdata[ii];
00962 }
00963 else if ( x2 == (nx - 1) && y2 == (ny - 1) ) {
00964 ddata[l] += Util::linear_interpolate(sdata[ii], sdata[ii + xy],v);
00965 }
00966 else if ( x2 == (nx - 1) && z2 == (nz - 1) ) {
00967 ddata[l] += Util::linear_interpolate(sdata[ii], sdata[ii + nx],u);
00968 }
00969 else if ( y2 == (ny - 1) && z2 == (nz - 1) ) {
00970 ddata[l] += Util::linear_interpolate(sdata[ii], sdata[ii + 1],t);
00971 }
00972 else if ( x2 == (nx - 1) ) {
00973 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + nx], sdata[ii + xy], sdata[ii + xy + nx],u,v);
00974 }
00975 else if ( y2 == (ny - 1) ) {
00976 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + xy], sdata[ii + xy + 1],t,v);
00977 }
00978 else if ( z2 == (nz - 1) ) {
00979 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx], sdata[ii + nx + 1],t,u);
00980 }
00981 }
00982 }
00983 }
00984 proj->update();
00985 proj->set_attr("xform.projection",t3d);
00986 proj->set_attr("apix_x",(float)image->get_attr("apix_x"));
00987 proj->set_attr("apix_y",(float)image->get_attr("apix_y"));
00988 proj->set_attr("apix_z",(float)image->get_attr("apix_z"));
00989
00990 if(t3d) {delete t3d; t3d=0;}
00991 return proj;
00992 }
00993 else if ( image->get_ndim() == 2 ) {
00994
00995 Transform r = t3d->inverse();
00996
00997 int nx = image->get_xsize();
00998 int ny = image->get_ysize();
00999
01000 EMData *proj = new EMData();
01001 proj->set_size(nx, 1, 1);
01002 proj->to_zero();
01003
01004 float *sdata = image->get_data();
01005 float *ddata = proj->get_data();
01006
01007 Vec2f offset(nx/2,ny/2);
01008 for (int j = -ny / 2; j < ny - ny / 2; j++) {
01009 int l = 0;
01010 for (int i = -nx / 2; i < nx - nx / 2; i++,l++) {
01011
01012 Vec2f coord(i,j);
01013 Vec2f soln = r*coord;
01014 soln += offset;
01015
01016 float x2 = soln[0];
01017 float y2 = soln[1];
01018
01019 float x = (float)Util::fast_floor(x2);
01020 float y = (float)Util::fast_floor(y2);
01021
01022 int ii = (int) (x + y * nx);
01023 float u = x2 - x;
01024 float v = y2 - y;
01025
01026 if (x2 < 0 || y2 < 0 ) continue;
01027 if (x2 > (nx-1) || y2 > (ny-1) ) continue;
01028
01029 if ( x2 < (nx - 1) && y2 < (ny - 1) ) {
01030 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx],sdata[ii + nx + 1], u, v);
01031 }
01032 else if (x2 == (nx-1) && y2 == (ny-1) ) {
01033 ddata[l] += sdata[ii];
01034 }
01035 else if (x2 == (nx-1) ) {
01036 ddata[l] += Util::linear_interpolate(sdata[ii],sdata[ii + nx], v);
01037 }
01038 else if (y2 == (ny-1) ) {
01039 ddata[l] += Util::linear_interpolate(sdata[ii],sdata[ii + 1], u);
01040 }
01041 }
01042 }
01043 proj->set_attr("xform.projection",t3d);
01044 proj->update();
01045 if(t3d) {delete t3d; t3d=0;}
01046 return proj;
01047 }
01048 else throw ImageDimensionException("Standard projection works only for 2D and 3D images");
01049 }
01050
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 transform 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 transform 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 EMData *FourierGriddingProjector::backproject3d(EMData * ) const
02091 {
02092
02093 EMData *ret = new EMData();
02094 return ret;
02095 }
02096
02097
02098
02099 void EMAN::dump_projectors()
02100 {
02101 dump_factory < Projector > ();
02102 }
02103
02104 map<string, vector<string> > EMAN::dump_projectors_list()
02105 {
02106 return dump_factory_list < Projector > ();
02107 }
02108
02109