EMAN2
projector.cpp
Go to the documentation of this file.
1/*
2 * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
3 * Copyright (c) 2000-2006 Baylor College of Medicine
4 *
5 * This software is issued under a joint BSD/GNU license. You may use the
6 * source code in this file under either license. However, note that the
7 * complete EMAN2 and SPARX software packages have some GPL dependencies,
8 * so you are responsible for compliance with the licenses of these packages
9 * if you opt to use BSD licensing. The warranty disclaimer below holds
10 * in either instance.
11 *
12 * This complete copyright notice must be included in any revised version of the
13 * source code. Additional authorship citations may be added, but existing
14 * author citations must be preserved.
15 *
16 * This program is free software; you can redistribute it and/or modify
17 * it under the terms of the GNU General Public License as published by
18 * the Free Software Foundation; either version 2 of the License, or
19 * (at your option) any later version.
20 *
21 * This program is distributed in the hope that it will be useful,
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 * GNU General Public License for more details.
25 *
26 * You should have received a copy of the GNU General Public License
27 * along with this program; if not, write to the Free Software
28 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29 *
30 * */
31
32#include "projector.h"
33#include "emdata.h"
34#include "interp.h"
35#include "emutil.h"
36#include "plugins/projector_template.h"
37
38#ifdef WIN32
39 #define M_PI 3.14159265358979323846f
40#endif //WIN32
41
42#ifdef EMAN2_USING_CUDA
43#include "cuda/cuda_util.h"
44#include "cuda/cuda_projector.h"
45#endif
46
47using namespace std;
48using namespace EMAN;
49
50const string GaussFFTProjector::NAME = "gauss_fft";
51const string FourierGriddingProjector::NAME = "fourier_gridding";
52const string PawelProjector::NAME = "pawel";
53const string StandardProjector::NAME = "standard";
54const string MaxValProjector::NAME = "maxval";
55const string ChaoProjector::NAME = "chao";
56
57template <> Factory < Projector >::Factory()
58{
59 force_add<GaussFFTProjector>();
60 force_add<PawelProjector>();
61 force_add<StandardProjector>();
62 force_add<MaxValProjector>();
63 force_add<FourierGriddingProjector>();
64 force_add<ChaoProjector>();
65
66// force_add<XYZProjector>();
67}
68
69EMData *GaussFFTProjector::project3d(EMData * image) const
70{
71 Transform* t3d = params["transform"];
72 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
73 if ( image->get_ndim() != 3 ) throw ImageDimensionException("Error, the projection volume must be 3D");
74
75 EMData *f = image;
76 if (!image->is_complex()) {
77 image->process_inplace("xform.phaseorigin.tocorner");
78 f = image->do_fft();
79 f->process_inplace("xform.fourierorigin.tocenter");
80 image->process_inplace("xform.phaseorigin.tocenter");
81 }
82
83 int f_nx = f->get_xsize();
84 int f_ny = f->get_ysize();
85 int f_nz = f->get_zsize();
86
87 if (!f->is_complex() || f_nz != f_ny || f_nx != f_ny + 2) {
88 LOGERR("Cannot project this image");
89 return 0;
90 }
91
92 f->ap2ri();
93
94 EMData *tmp = new EMData();
95 tmp->set_size(f_nx, f_ny, 1);
96 tmp->set_complex(true);
97 tmp->set_ri(true);
98
99 float *data = tmp->get_data();
100
101 Transform r = t3d->get_rotation_transform();
102 r.invert();
103 float scale = t3d->get_scale();
104
105 int mode = params["mode"];
106 float gauss_width = 1;
107// if ( mode == 0 ) mode = 2;
108 if (mode == 2 ) {
109 gauss_width = EMConsts::I2G;
110 }
111 else if (mode == 3) {
112 gauss_width = EMConsts::I3G;
113 }
114 else if (mode == 4) {
115 gauss_width = EMConsts::I4G;
116 }
117 else if (mode == 6 || mode == 7) {
118 gauss_width = EMConsts::I5G;
119 }
120
121 for (int y = 0; y < f_ny; y++) {
122 for (int x = 0; x < f_nx / 2; x++) {
123 int ii = x * 2 + y * f_nx;
124 if (hypot(x, y - f_ny / 2) >= f_ny / 2 - 1) {
125 data[ii] = 0;
126 data[ii + 1] = 0;
127 }
128 else {
129 Vec3f coord(x,(y - f_ny / 2),0);
130 Vec3f soln = r*coord;
131 float xx = soln[0];
132 float yy = soln[1];
133 float zz = soln[2];
134
135 int cc = 1;
136 if (xx < 0) {
137 xx = -xx;
138 yy = -yy;
139 zz = -zz;
140 cc = -1;
141 }
142
143 if (scale != 1.0) {
144 xx *= scale;
145 yy *= scale;
146 zz *= scale;
147 }
148
149 yy += f_ny / 2;
150 zz += f_nz / 2;
151
152 if (yy < 0 || xx < 0 || zz < 0) {
153 data[ii] = 0;
154 data[ii+1] = 0;
155 continue;
156 }
157
158 if (interp_ft_3d(mode, f, xx, yy, zz, data + ii, gauss_width)) {
159 data[ii + 1] *= cc;
160 } else {
161 data[ii] = 0;
162 data[ii+1] = 0;
163 }
164 }
165 }
166 }
167
168// f->update();
169 tmp->update();
170
171 int returnfft=params["returnfft"];
172 EMData *ret;
173 if (returnfft==1){
174 //apperently there is something wrong with translating images in fourier space...
175// printf("return fft!\n");
176 Vec3f trans=t3d->get_trans();
177 ret=tmp->copy();
178 ret->process_inplace("xform.fourierorigin.tocorner");
179 ret->process_inplace("xform", Dict("tx", (float)trans[0], "ty", (float)trans[1]));
180
181// if (t3d->get_mirror() ) ret->process_inplace("xform.flip",Dict("axis","x"));
182
183 }
184 else{
185
186 tmp->process_inplace("xform.fourierorigin.tocorner");
187 ret = tmp->do_ift();
188 ret->process_inplace("xform.phaseorigin.tocenter");
189
190 ret->translate(t3d->get_trans());
191
192 if (t3d->get_mirror() ) ret->process_inplace("xform.flip",Dict("axis","x"));
193
194 Dict filter_d;
195 filter_d["gauss_width"] = gauss_width;
196 // filter_d["ring_width"] = ret->get_xsize() / 2;
197 ret->process_inplace("math.gausskernelfix", filter_d);
198 }
199
200 if( tmp )
201 {
202 delete tmp;
203 tmp = 0;
204 }
205
206 ret->set_attr("xform.projection",t3d);
207 ret->update();
208
209 if(t3d) {delete t3d; t3d=0;}
210
211 return ret;
212}
213
214
215
216bool GaussFFTProjector::interp_ft_3d(int mode, EMData * image, float x, float y,
217 float z, float *data, float gw) const
218{
219 float *rdata = image->get_data();
220 int nx = image->get_xsize();
221 int ny = image->get_ysize();
222 int nz = image->get_zsize();
223
224 if ( mode == 0 ) {
225 int x0 = (int) floor(x);
226 int y0 = (int) floor(y);
227 int z0 = (int) floor(z);
228 float d0,d1,nrm;
229 d0=d1=nrm=0;
230
231 for (int ix=x0; ix<=x0+1; ix++){
232 for (int iy=y0; iy<=y0+1; iy++){
233 for (int iz=z0; iz<=z0+1; iz++){
234 float w=(1.0-fabs(ix-x))*(1.0-fabs(iy-y))*(1.0-fabs(iz-z));
235 d0 += w * rdata[2*ix + iy * nx + iz * nx * ny];
236 d1 += w * rdata[2*ix + iy * nx + iz * nx * ny + 1];
237 nrm+=w;
238 }}}
239// printf("%f\t%f\t%f\n", d0, d1, nrm);
240 data[0]=d0/nrm;
241 data[1]=d1/nrm;
242 return true;
243 }
244 if (mode == 1) {
245 int x0 = 2 * (int) floor(x + 0.5f);
246 int y0 = (int) floor(y + 0.5f);
247 int z0 = (int) floor(z + 0.5f);
248
249 data[0] = rdata[x0 + y0 * nx + z0 * nx * ny];
250 data[1] = rdata[x0 + y0 * nx + z0 * nx * ny + 1];
251 return true;
252 }
253 else if (mode == 2) {
254 int x0 = (int) floor(x);
255 int y0 = (int) floor(y);
256 int z0 = (int) floor(z);
257
258 float dx = x - x0;
259 float dy = y - y0;
260 float dz = z - z0;
261
262 if (x0 <= nx/2- 2 && y0 <= ny - 1 && z0 <= nz - 1) {
263 int i = (int) (x0 * 2 + y0 * nx + z0 * nx * ny);
264
265 float n = Util::agauss(1, dx, dy, dz, gw) +
266 Util::agauss(1, 1 - dx, dy, dz, gw) +
267 Util::agauss(1, dx, 1 - dy, dz, gw) +
268 Util::agauss(1, 1 - dx, 1 - dy, dz, gw) +
269 Util::agauss(1, dx, dy, 1 - dz, gw) +
270 Util::agauss(1, 1 - dx, dy, 1 - dz, gw) +
271 Util::agauss(1, dx, 1 - dy, 1 - dz, gw) + Util::agauss(1, 1 - dx, 1 - dy, 1 - dz, gw);
272
273 data[0] = Util::agauss(rdata[i], dx, dy, dz, gw) +
274 Util::agauss(rdata[i + 2], 1 - dx, dy, dz, gw) +
275 Util::agauss(rdata[i + nx], dx, 1 - dy, dz, gw) +
276 Util::agauss(rdata[i + nx + 2], 1 - dx, 1 - dy, dz, gw) +
277 Util::agauss(rdata[i + nx * ny], dx, dy, 1 - dz, gw) +
278 Util::agauss(rdata[i + 2 + nx * ny], 1 - dx, dy, 1 - dz, gw) +
279 Util::agauss(rdata[i + nx + nx * ny], dx, 1 - dy, 1 - dz, gw) +
280 Util::agauss(rdata[i + 2 + nx + nx * ny], 1 - dx, 1 - dy, 1 - dz, gw) / n;
281
282 i++;
283
284 data[1] = Util::agauss(rdata[i], dx, dy, dz, gw) +
285 Util::agauss(rdata[i + 2], 1 - dx, dy, dz, gw) +
286 Util::agauss(rdata[i + nx], dx, 1 - dy, dz, gw) +
287 Util::agauss(rdata[i + nx + 2], 1 - dx, 1 - dy, dz, gw) +
288 Util::agauss(rdata[i + nx * ny], dx, dy, 1 - dz, gw) +
289 Util::agauss(rdata[i + 2 + nx * ny], 1 - dx, dy, 1 - dz, gw) +
290 Util::agauss(rdata[i + nx + nx * ny], dx, 1 - dy, 1 - dz, gw) +
291 Util::agauss(rdata[i + 2 + nx + nx * ny], 1 - dx, 1 - dy, 1 - dz, gw) / n;
292 return true;
293 }
294 return false;
295 }
296 else if (mode == 3) {
297 int x0 = 2 * (int) floor(x + .5);
298 int y0 = (int) floor(y + .5);
299 int z0 = (int) floor(z + .5);
300
301 float *supp = image->setup4slice();
302
303 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
304 float n = 0;
305
306 if (x0 == 0) {
307 x0 += 4;
308 size_t idx;
309 float r, g;
310 for (int k = z0 - 1; k <= z0 + 1; k++) {
311 for (int j = y0 - 1; j <= y0 + 1; j++) {
312 for (int i = x0 - 2; i <= x0 + 2; i += 2) {
313 r = Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z);
314 g = exp(-r / gw);
315 n += g;
316 idx = i + j * 12 + (size_t)k * 12 * ny;
317 data[0] += supp[idx] * g;
318 data[1] += supp[idx + 1] * g;
319 }
320 }
321 }
322 }
323 else {
324 size_t idx;
325 float r, g;
326 for (int k = z0 - 1; k <= z0 + 1; k++) {
327 for (int j = y0 - 1; j <= y0 + 1; j++) {
328 for (int i = x0 - 2; i <= x0 + 2; i += 2) {
329 r = Util::hypot3sq(i / 2.0f - x, j - y, k - z);
330 g = exp(-r / gw);
331 n += g;
332 idx = i + j * nx + (size_t)k * nx * ny;
333 data[0] += rdata[idx] * g;
334 data[1] += rdata[idx + 1] * g;
335 }
336 }
337 }
338 }
339
340 data[0] /= n;
341 data[1] /= n;
342 return true;
343 }
344 return false;
345 }
346 else if (mode == 4) {
347 int x0 = 2 * (int) floor(x);
348 int y0 = (int) floor(y);
349 int z0 = (int) floor(z);
350
351 float *supp = image->setup4slice();
352
353 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
354 float n = 0;
355
356 if (x0 == 0) {
357 x0 += 4;
358 size_t idx;
359 float r, g;
360 for (int k = z0 - 1; k <= z0 + 2; k++) {
361 for (int j = y0 - 1; j <= y0 + 2; j++) {
362 for (int i = x0 - 2; i <= x0 + 4; i += 2) {
363 r = Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z);
364 g = exp(-r / gw);
365 n += g;
366 idx = i + j * 12 + (size_t)k * 12 * ny;
367 data[0] += supp[idx] * g;
368 data[1] += supp[idx + 1] * g;
369 }
370 }
371 }
372 }
373 else {
374 float r, g;
375 size_t idx;
376 for (int k = z0 - 1; k <= z0 + 2; k++) {
377 for (int j = y0 - 1; j <= y0 + 2; j++) {
378 for (int i = x0 - 2; i <= x0 + 4; i += 2) {
379 r = Util::hypot3sq(i / 2.0f - x, j - y, k - z);
380 g = exp(-r / gw);
381 n += g;
382 idx = i + j * nx + (size_t)k * nx * ny;
383 data[0] += rdata[idx] * g;
384 data[1] += rdata[idx + 1] * g;
385 }
386 }
387 }
388 }
389 data[0] /= n;
390 data[1] /= n;
391 return true;
392 }
393 return false;
394 }
395 else if (mode == 5) {
396 int x0 = 2 * (int) floor(x + .5);
397 int y0 = (int) floor(y + .5);
398 int z0 = (int) floor(z + .5);
399
400 float *supp = image->setup4slice();
401 float *gimx = Interp::get_gimx();
402
403 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
404 int mx0 = -(int) floor((x - x0 / 2) * 39.0 + .5) - 78;
405 int my0 = -(int) floor((y - y0) * 39.0 + .5) - 78;
406 int mz0 = -(int) floor((z - z0) * 39.0 + .5) - 78;
407
408 float n = 0;
409 int mmz = mz0;
410 int mmy = my0;
411 int mmx = mx0;
412
413 if (x0 < 4) {
414 x0 += 4;
415
416 size_t idx;
417 float g;
418 for (int k = z0 - 2; k <= z0 + 2; k++, mmz += 39) {
419 for (int j = y0 - 2; j <= y0 + 2; j++, mmy += 39) {
420 for (int i = x0 - 4; i <= x0 + 4; i += 2, mmx += 39) {
421 g = gimx[abs(mmx) + abs(mmy) * 100 + abs(mmz) * 10000];
422 n += g;
423 idx = i + j * 12 + (size_t)k * 12 * ny;
424 data[0] += supp[idx] * g;
425 data[1] += supp[idx + 1] * g;
426 }
427 }
428 }
429 }
430 else {
431 size_t ii;
432 float g;
433 for (int k = z0 - 2; k <= z0 + 2; k++, mmz += 39) {
434 for (int j = y0 - 2; j <= y0 + 2; j++, mmy += 39) {
435 for (int i = x0 - 4; i <= x0 + 4; i += 2, mmx += 39) {
436 ii = i + j * nx + (size_t)k * nx * ny;
437 g = gimx[abs(mmx) + abs(mmy) * 100 + abs(mmz) * 10000];
438 n += g;
439 data[0] += rdata[ii] * g;
440 data[1] += rdata[ii + 1] * g;
441 }
442 }
443 }
444 }
445
446 data[0] /= n;
447 data[1] /= n;
448 return true;
449 }
450 return false;
451 }
452 else if (mode == 6) {
453
454 int x0 = 2 * (int) floor(x + .5);
455 int y0 = (int) floor(y + .5);
456 int z0 = (int) floor(z + .5);
457
458 float *supp = image->setup4slice();
459
460 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
461 float n = 0;
462
463 if (x0 < 4) {
464 x0 += 4;
465 float r, g;
466 size_t idx;
467 for (int k = z0 - 2; k <= z0 + 2; k++) {
468 for (int j = y0 - 2; j <= y0 + 2; j++) {
469 for (int i = x0 - 4; i <= x0 + 4; i += 2) {
470 r = Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z);
471 g = exp(-r / gw);
472 n += g;
473 idx = i + j * 12 + (size_t)k * 12 * ny;
474 data[0] += supp[idx] * g;
475 data[1] += supp[idx + 1] * g;
476 }
477 }
478 }
479 }
480 else {
481 float r, g;
482 size_t idx;
483 for (int k = z0 - 2; k <= z0 + 2; k++) {
484 for (int j = y0 - 2; j <= y0 + 2; j++) {
485 for (int i = x0 - 4; i <= x0 + 4; i += 2) {
486 r = Util::hypot3sq(i / 2.0f - x, j - y, k - z);
487 g = exp(-r / gw);
488 n += g;
489 idx = i + j * nx + (size_t)k * nx * ny;
490 data[0] += rdata[idx] * g;
491 data[1] += rdata[idx + 1] * g;
492 }
493
494 }
495 }
496 }
497
498 data[0] /= n;
499 data[1] /= n;
500
501 return true;
502 }
503 return false;
504 }
505 else if (mode == 7) {
506 int x0 = 2 * (int) floor(x + .5);
507 int y0 = (int) floor(y + .5);
508 int z0 = (int) floor(z + .5);
509
510 float *supp = image->setup4slice();
511
512 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
513 float n = 0;
514 if (x0 < 4) {
515 x0 += 4;
516 float r, g;
517 size_t idx;
518 for (int k = z0 - 2; k <= z0 + 2; k++) {
519 for (int j = y0 - 2; j <= y0 + 2; j++) {
520 for (int i = x0 - 4; i <= x0 + 4; i += 2) {
521 r = sqrt(Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z));
522 g = Interp::hyperg(r);
523 n += g;
524 idx = i + j * 12 + (size_t)k * 12 * ny;
525 data[0] += supp[idx] * g;
526 data[1] += supp[idx + 1] * g;
527 }
528 }
529 }
530 }
531 else {
532 float r, g;
533 size_t idx;
534 for (int k = z0 - 2; k <= z0 + 2; k++) {
535 for (int j = y0 - 2; j <= y0 + 2; j++) {
536 for (int i = x0 - 4; i <= x0 + 4; i += 2) {
537 r = sqrt(Util::hypot3sq(i / 2.0f - x, j - y, k - z));
538 g = Interp::hyperg(r);
539 n += g;
540 idx = i + j * nx + (size_t)k * nx * ny;
541 data[0] += rdata[idx] * g;
542 data[1] += rdata[idx + 1] * g;
543 }
544
545 }
546 }
547 }
548 data[0] /= n;
549 data[1] /= n;
550 return true;
551 }
552 return false;
553
554 }
555// throw InvalidParameterException("Error, unkown mode");
556 return false;
557}
558
559void PawelProjector::prepcubes(int nx, int ny, int nz, int ri, Vec3i origin,
560 int& nn, IPCube* ipcube) const {
561 const float r = float(ri)*float(ri);
562 const int ldpx = origin[0];
563 const int ldpy = origin[1];
564 const int ldpz = origin[2];
565 //cout<<" ldpx "<<ldpx<<" i2 "<<ldpy<<" i1 "<<ldpz<<endl;
566 float t;
567 nn = -1;
568 for (int i1 = 0; i1 < nz; i1++) {
569 t = float(i1 - ldpz);
570 const float xx = t*t;
571 for (int i2 = 0; i2 < ny; i2++) {
572 t = float(i2 - ldpy);
573 const float yy = t*t + xx;
574 bool first = true;
575 for (int i3 = 0; i3 < nx; i3++) {
576 t = float(i3 - ldpx);
577 const float rc = t*t + yy;
578 if (first) {
579 // first pixel on this line
580 if (rc > r) continue;
581 first = false;
582 nn++;
583 if (ipcube != NULL) {
584 ipcube[nn].start = i3;
585 ipcube[nn].end = i3;
586 ipcube[nn].loc[0] = i3 - ldpx;
587 ipcube[nn].loc[1] = i2 - ldpy;
588 ipcube[nn].loc[2] = i1 - ldpz;
589 //cout<<" start "<<i3<<" i2 "<<i2<<" i1 "<<i1<<endl;
590 }
591 } else {
592 // second or later pixel on this line
593 if (ipcube != NULL) {
594 if (rc <= r) ipcube[nn].end = i3;
595 }
596 }
597 }
598 }
599 }
600}
601
602EMData *PawelProjector::project3d(EMData * image) const
603{
604 if (!image) return 0;
605 int ri;
606 int nx = image->get_xsize();
607 int ny = image->get_ysize();
608 int nz = image->get_zsize();
609 int dim = Util::get_max(nx,ny,nz);
610 int dimc = dim/2;
611 if (nz == 1) {
612 LOGERR("The PawelProjector needs a volume!");
613 return 0;
614 }
615
616 Vec3i origin(0,0,0);
617
618 if (params.has_key("radius")) {ri = params["radius"];}
619 else {ri = dim/2 - 1;}
620
621 Transform* rotation = params["transform"];
622 //int nangles = 0;
623 //vector<float> anglelist;
624 // Do we have a list of angles?
625 /*
626 if (params.has_key("anglelist")) {
627 anglelist = params["anglelist"];
628 nangles = anglelist.size() / 3;
629 } else {*/
630
631 //if ( rotation == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
632 /*
633 Dict p = t3d->get_rotation("spider");
634
635 string angletype = "SPIDER";
636 float phi = p["phi"];
637 float theta = p["theta"];
638 float psi = p["psi"];
639 anglelist.push_back(phi);
640 anglelist.push_back(theta);
641 anglelist.push_back(psi);
642 */
643 int nangles = 1;
644 //}
645
646 //for (int i = 0 ; i <= nn; i++) cout<<" IPCUBE "<<ipcube[i].start<<" "<<ipcube[i].end<<" "<<ipcube[i].loc[0]<<" "<<ipcube[i].loc[1]<<" "<<ipcube[i].loc[2]<<endl;
647 // loop over sets of angles
648 //for (int ia = 0; ia < nangles; ia++) {
649 EMData* ret = new EMData();
650 int ia = 0;
651 //int indx = 3*ia;
652 //Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]);
653 //Transform rotation(d);
654 if (2*(ri+1)+1 > dim) {
655 // initialize return object
656 ret->set_size(dim, dim, nangles);
657 ret->to_zero();
658 origin[0] = dimc;
659 origin[1] = dimc;
660 origin[2] = dimc;
661 Vec3i loc(-dimc,0,0);
662 Vec3i vorg(nx/2,ny/2,nz/2);
663 // This code is for arbitrary dimensions, so must check x and y boundaries
664 for (int l = 0 ; l < dim; l++) { // Z
665 loc[2] = l - dimc;
666 for (int k = 0 ; k < dim; k++) { // Y
667 loc[1] = k - dimc;
668 Vec3f vb = loc*(*rotation) + vorg;
669 for (int j = 0; j < dim; j++) { //X
670 // check for pixels out-of-bounds
671 //cout<<j<<" j "<<k<<" k "<<" iox "<<vb[0]<<" ioy "<<vb[1]<<" ioz "<<vb[2]<<endl;
672 int iox = int(vb[0]);
673 if ((iox >= 0) && (iox < nx-1)) {
674 int ioy = int(vb[1]);
675 if ((ioy >= 0) && (ioy < ny-1)) {
676 int ioz = int(vb[2]);
677 if ((ioz >= 0) && (ioz < nz-1)) {
678 // real work for pixels in bounds
679 //cout<<j<<" j "<<k<<" k "<<" iox "<<iox<<" ioy "<<ioy<<" ioz "<<ioz<<endl;
680 //cout<<" TAKE"<<endl;
681 float dx = vb[0] - iox;
682 float dy = vb[1] - ioy;
683 float dz = vb[2] - ioz;
684 float a1 = (*image)(iox,ioy,ioz);
685 float a2 = (*image)(iox+1,ioy,ioz) - a1;
686 float a3 = (*image)(iox,ioy+1,ioz) - a1;
687 float a4 = (*image)(iox,ioy,ioz+1) - a1;
688 float a5 = -a2 -(*image)(iox,ioy+1,ioz)
689 + (*image)(iox+1,ioy+1,ioz);
690 float a61 = -(*image)(iox,ioy,ioz+1)
691 + (*image)(iox+1,ioy,ioz+1);
692 float a6 = -a2 + a61;
693 float a7 = -a3 - (*image)(iox,ioy,ioz+1)
694 + (*image)(iox,ioy+1,ioz+1);
695 float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
696 + (*image)(iox+1,ioy+1,ioz+1);
697 (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
698 + (a7 + a8*dx)*dy)
699 + a3*dy + dx*(a2 + a5*dy);
700 } //else {cout<<" iox "<<iox<<" ioy "<<int(vb[1])<<" ioz "<<int(vb[2])<<" VIOLAATED Z "<<endl; }
701 }//else {cout<<" iox "<<iox<<" ioy "<<int(vb[1])<<" ioz "<<int(vb[2])<<" VIOLATED Y "<<endl; }
702 }//else {cout<<" iox "<<iox<<" ioy "<<int(vb[1])<<" ioz "<<int(vb[2])<<" VIOLATED X "<<endl; }
703 vb += rotation->get_matrix3_row(0);
704 }
705 }
706 }
707
708 } else {
709 // If a sensible origin isn't passed in, choose the middle of
710 // the cube.
711 if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
712 else {origin[0] = nx/2;}
713 if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
714 else {origin[1] = ny/2;}
715 if (params.has_key("origin_z")) {origin[2] = params["origin_z"];}
716 else {origin[2] = nz/2;}
717 // Determine the number of rows (x-lines) within the radius
718 int nn = -1;
719 prepcubes(nx, ny, nz, ri, origin, nn);
720 // nn is now the number of rows-1 within the radius
721 // so we can create and fill the ipcubes
722 IPCube* ipcube = new IPCube[nn+1];
723 prepcubes(nx, ny, nz, ri, origin, nn, ipcube);
724 // initialize return object
725 ret->set_size(nx, ny, nangles);
726 ret->to_zero();
727 // No need to check x and y boundaries
728 for (int i = 0 ; i <= nn; i++) {
729 int k = ipcube[i].loc[1] + origin[1];
730 Vec3f vb = ipcube[i].loc*(*rotation) + origin;
731 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
732 int iox = int(vb[0]);
733 int ioy = int(vb[1]);
734 int ioz = int(vb[2]);
735 float dx = vb[0] - iox;
736 float dy = vb[1] - ioy;
737 float dz = vb[2] - ioz;
738 float a1 = (*image)(iox,ioy,ioz);
739 float a2 = (*image)(iox+1,ioy,ioz) - a1;
740 float a3 = (*image)(iox,ioy+1,ioz) - a1;
741 float a4 = (*image)(iox,ioy,ioz+1) - a1;
742 float a5 = -a2 -(*image)(iox,ioy+1,ioz)
743 + (*image)(iox+1,ioy+1,ioz);
744 float a61 = -(*image)(iox,ioy,ioz+1)
745 + (*image)(iox+1,ioy,ioz+1);
746 float a6 = -a2 + a61;
747 float a7 = -a3 - (*image)(iox,ioy,ioz+1)
748 + (*image)(iox,ioy+1,ioz+1);
749 float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
750 + (*image)(iox+1,ioy+1,ioz+1);
751 (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
752 + (a7 + a8*dx)*dy)
753 + a3*dy + dx*(a2 + a5*dy);
754 vb += rotation->get_matrix3_row(0);
755 }
756 }
757 EMDeleteArray(ipcube);
758 }
759 //}
760 ret->update();
761 if(rotation) {delete rotation; rotation=0;}
762 return ret;
763}
764
765EMData *StandardProjector::project3d(EMData * image) const
766{
767 Transform* t3d = params["transform"];
768 if ( t3d == NULL ) throw NullPointerException("The transform object containing the angles(required for projection), was not specified");
769// Dict p = t3d->get_rotation();
770 if ( image->get_ndim() == 3 )
771 {
772
773#ifdef EMAN2_USING_CUDA
774 if(EMData::usecuda == 1) {
775 if(!image->isrodataongpu()) image->copy_to_cudaro();
776 //cout << "CUDA PROJ" << endl;
777 Transform* t3d = params["transform"];
778 if ( t3d == NULL ) throw NullPointerException("The transform object containing the angles(required for projection), was not specified");
779 float * m = new float[12];
780 t3d->copy_matrix_into_array(m);
781 image->bindcudaarrayA(true);
782 //EMData* e = new EMData(0,0,image->get_xsize(),image->get_ysize(),1);
783 EMData *e = new EMData();
784 e->set_size_cuda(image->get_xsize(), image->get_ysize(), 1);
785 e->rw_alloc();
786 standard_project(m,e->getcudarwdata(), e->get_xsize(), e->get_ysize());
787 image->unbindcudaarryA();
788 delete [] m;
789
790 e->update();
791 e->set_attr("xform.projection",t3d);
792 e->set_attr("apix_x",(float)image->get_attr("apix_x"));
793 e->set_attr("apix_y",(float)image->get_attr("apix_y"));
794 e->set_attr("apix_z",(float)image->get_attr("apix_z"));
795 //e_>copy_from_device();
796 if(t3d) {delete t3d; t3d=0;}
797 return e;
798 }
799#endif
800 int nx = image->get_xsize();
801 int ny = image->get_ysize();
802 int nz = image->get_zsize();
803
804// Transform3D r(Transform3D::EMAN, az, alt, phi);
805 Transform r = t3d->inverse(); // The inverse is taken here because we are rotating the coordinate system, not the image
806 int xy = nx * ny;
807
808 EMData *proj = new EMData();
809 proj->set_size(nx, ny, 1);
810
811 Vec3i offset(nx/2,ny/2,nz/2);
812
813 float *sdata = image->get_data();
814 float *ddata = proj->get_data();
815 for (int k = -nz / 2; k < nz - nz / 2; k++) {
816 int l = 0;
817 for (int j = -ny / 2; j < ny - ny / 2; j++) {
818 ddata[l]=0;
819 for (int i = -nx / 2; i < nx - nx / 2; i++,l++) {
820
821 Vec3f coord(i,j,k);
822 Vec3f soln = r*coord;
823 soln += offset;
824
828// printf(" ");
829
830 float x2 = soln[0];
831 float y2 = soln[1];
832 float z2 = soln[2];
833
834 float x = (float)Util::fast_floor(x2);
835 float y = (float)Util::fast_floor(y2);
836 float z = (float)Util::fast_floor(z2);
837
838 float t = x2 - x;
839 float u = y2 - y;
840 float v = z2 - z;
841
842 size_t ii = (size_t) ((size_t)x + (size_t)y * nx + (size_t)z * xy);
843//
844 if (x2 < 0 || y2 < 0 || z2 < 0 ) continue;
845 if (x2 > (nx-1) || y2 > (ny-1) || z2 > (nz-1) ) continue;
846
847 if (x2 < (nx - 1) && y2 < (ny - 1) && z2 < (nz - 1)) {
848 ddata[l] +=
849 Util::trilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx],
850 sdata[ii + nx + 1], sdata[ii + xy], sdata[ii + xy + 1], sdata[ii + xy + nx],
851 sdata[ii + xy + nx + 1], t, u, v);
852 }
853 else if ( x2 == (nx - 1) && y2 == (ny - 1) && z2 == (nz - 1) ) {
854 ddata[l] += sdata[ii];
855 }
856 else if ( x2 == (nx - 1) && y2 == (ny - 1) ) {
857 ddata[l] += Util::linear_interpolate(sdata[ii], sdata[ii + xy],v);
858 }
859 else if ( x2 == (nx - 1) && z2 == (nz - 1) ) {
860 ddata[l] += Util::linear_interpolate(sdata[ii], sdata[ii + nx],u);
861 }
862 else if ( y2 == (ny - 1) && z2 == (nz - 1) ) {
863 ddata[l] += Util::linear_interpolate(sdata[ii], sdata[ii + 1],t);
864 }
865 else if ( x2 == (nx - 1) ) {
866 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + nx], sdata[ii + xy], sdata[ii + xy + nx],u,v);
867 }
868 else if ( y2 == (ny - 1) ) {
869 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + xy], sdata[ii + xy + 1],t,v);
870 }
871 else if ( z2 == (nz - 1) ) {
872 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx], sdata[ii + nx + 1],t,u);
873 }
874 }
875 }
876 }
877 proj->update();
878 proj->set_attr("xform.projection",t3d);
879 proj->set_attr("apix_x",(float)image->get_attr("apix_x"));
880 proj->set_attr("apix_y",(float)image->get_attr("apix_y"));
881 proj->set_attr("apix_z",(float)image->get_attr("apix_z"));
882
883 if(t3d) {delete t3d; t3d=0;}
884 return proj;
885 }
886 else if ( image->get_ndim() == 2 ) {
887
888 Transform r = t3d->inverse(); // The inverse is taken here because we are rotating the coordinate system, not the image
889
890 int nx = image->get_xsize();
891 int ny = image->get_ysize();
892
893 EMData *proj = new EMData();
894 proj->set_size(nx, 1, 1);
895 proj->to_zero();
896
897 float *sdata = image->get_data();
898 float *ddata = proj->get_data();
899
900 Vec2f offset(nx/2,ny/2);
901 for (int j = -ny / 2; j < ny - ny / 2; j++) { // j represents a column of pixels in the direction of the angle
902 int l = 0;
903 for (int i = -nx / 2; i < nx - nx / 2; i++,l++) {
904
905 Vec2f coord(i,j);
906 Vec2f soln = r*coord;
907 soln += offset;
908
909 float x2 = soln[0];
910 float y2 = soln[1];
911
912 float x = (float)Util::fast_floor(x2);
913 float y = (float)Util::fast_floor(y2);
914
915 int ii = (int) (x + y * nx);
916 float u = x2 - x;
917 float v = y2 - y;
918
919 if (x2 < 0 || y2 < 0 ) continue;
920 if (x2 > (nx-1) || y2 > (ny-1) ) continue;
921
922 if ( x2 < (nx - 1) && y2 < (ny - 1) ) {
923 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx],sdata[ii + nx + 1], u, v);
924 }
925 else if (x2 == (nx-1) && y2 == (ny-1) ) {
926 ddata[l] += sdata[ii];
927 }
928 else if (x2 == (nx-1) ) {
929 ddata[l] += Util::linear_interpolate(sdata[ii],sdata[ii + nx], v);
930 }
931 else if (y2 == (ny-1) ) {
932 ddata[l] += Util::linear_interpolate(sdata[ii],sdata[ii + 1], u);
933 }
934 }
935 }
936 proj->set_attr("xform.projection",t3d);
937 proj->update();
938 if(t3d) {delete t3d; t3d=0;}
939 return proj;
940 }
941 else throw ImageDimensionException("Standard projection works only for 2D and 3D images");
942}
943
944EMData *MaxValProjector::project3d(EMData * image) const
945{
946 Transform* t3d = params["transform"];
947 if ( t3d == NULL ) throw NullPointerException("The transform object containing the angles(required for projection), was not specified");
948// Dict p = t3d->get_rotation();
949 if ( image->get_ndim() == 3 )
950 {
951
952 int nx = image->get_xsize();
953 int ny = image->get_ysize();
954 int nz = image->get_zsize();
955
956// Transform3D r(Transform3D::EMAN, az, alt, phi);
957 Transform r = t3d->inverse(); // The inverse is taken here because we are rotating the coordinate system, not the image
958 int xy = nx * ny;
959
960 EMData *proj = new EMData();
961 proj->set_size(nx, ny, 1);
962
963 Vec3i offset(nx/2,ny/2,nz/2);
964
965 float *sdata = image->get_data();
966 float *ddata = proj->get_data();
967 for (int k = -nz / 2; k < nz - nz / 2; k++) {
968 int l = 0;
969 for (int j = -ny / 2; j < ny - ny / 2; j++) {
970 ddata[l]=0;
971 for (int i = -nx / 2; i < nx - nx / 2; i++,l++) {
972
973 Vec3f coord(i,j,k);
974 Vec3f soln = r*coord;
975 soln += offset;
976
980// printf(" ");
981
982 float x2 = soln[0];
983 float y2 = soln[1];
984 float z2 = soln[2];
985
986 float x = (float)Util::fast_floor(x2);
987 float y = (float)Util::fast_floor(y2);
988 float z = (float)Util::fast_floor(z2);
989
990 float t = x2 - x;
991 float u = y2 - y;
992 float v = z2 - z;
993
994 size_t ii = (size_t) ((size_t)x + (size_t)y * nx + (size_t)z * xy);
995//
996 if (x2 < 0 || y2 < 0 || z2 < 0 ) continue;
997 if (x2 > (nx-1) || y2 > (ny-1) || z2 > (nz-1) ) continue;
998
999 if (x2 < (nx - 1) && y2 < (ny - 1) && z2 < (nz - 1)) {
1000 ddata[l] = Util::get_max(ddata[l],
1001 Util::trilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx],
1002 sdata[ii + nx + 1], sdata[ii + xy], sdata[ii + xy + 1], sdata[ii + xy + nx],
1003 sdata[ii + xy + nx + 1], t, u, v));
1004 }
1005 else if ( x2 == (nx - 1) && y2 == (ny - 1) && z2 == (nz - 1) ) {
1006 ddata[l] = Util::get_max(ddata[l],sdata[ii]);
1007 }
1008 else if ( x2 == (nx - 1) && y2 == (ny - 1) ) {
1009 ddata[l] = Util::get_max(ddata[l],Util::linear_interpolate(sdata[ii], sdata[ii + xy],v));
1010 }
1011 else if ( x2 == (nx - 1) && z2 == (nz - 1) ) {
1012 ddata[l] = Util::get_max(ddata[l],Util::linear_interpolate(sdata[ii], sdata[ii + nx],u));
1013 }
1014 else if ( y2 == (ny - 1) && z2 == (nz - 1) ) {
1015 ddata[l] = Util::get_max(ddata[l],Util::linear_interpolate(sdata[ii], sdata[ii + 1],t));
1016 }
1017 else if ( x2 == (nx - 1) ) {
1018 ddata[l] = Util::get_max(ddata[l],Util::bilinear_interpolate(sdata[ii], sdata[ii + nx], sdata[ii + xy], sdata[ii + xy + nx],u,v));
1019 }
1020 else if ( y2 == (ny - 1) ) {
1021 ddata[l] = Util::get_max(ddata[l],Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + xy], sdata[ii + xy + 1],t,v));
1022 }
1023 else if ( z2 == (nz - 1) ) {
1024 ddata[l] = Util::get_max(ddata[l],Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx], sdata[ii + nx + 1],t,u));
1025 }
1026 }
1027 }
1028 }
1029 proj->update();
1030 proj->set_attr("xform.projection",t3d);
1031 proj->set_attr("apix_x",(float)image->get_attr("apix_x"));
1032 proj->set_attr("apix_y",(float)image->get_attr("apix_y"));
1033 proj->set_attr("apix_z",(float)image->get_attr("apix_z"));
1034
1035 if(t3d) {delete t3d; t3d=0;}
1036 return proj;
1037 }
1038 else if ( image->get_ndim() == 2 ) {
1039
1040 Transform r = t3d->inverse(); // The inverse is taken here because we are rotating the coordinate system, not the image
1041
1042 int nx = image->get_xsize();
1043 int ny = image->get_ysize();
1044
1045 EMData *proj = new EMData();
1046 proj->set_size(nx, 1, 1);
1047 proj->to_zero();
1048
1049 float *sdata = image->get_data();
1050 float *ddata = proj->get_data();
1051
1052 Vec2f offset(nx/2,ny/2);
1053 for (int j = -ny / 2; j < ny - ny / 2; j++) { // j represents a column of pixels in the direction of the angle
1054 int l = 0;
1055 for (int i = -nx / 2; i < nx - nx / 2; i++,l++) {
1056
1057 Vec2f coord(i,j);
1058 Vec2f soln = r*coord;
1059 soln += offset;
1060
1061 float x2 = soln[0];
1062 float y2 = soln[1];
1063
1064 float x = (float)Util::fast_floor(x2);
1065 float y = (float)Util::fast_floor(y2);
1066
1067 int ii = (int) (x + y * nx);
1068 float u = x2 - x;
1069 float v = y2 - y;
1070
1071 if (x2 < 0 || y2 < 0 ) continue;
1072 if (x2 > (nx-1) || y2 > (ny-1) ) continue;
1073
1074 if ( x2 < (nx - 1) && y2 < (ny - 1) ) {
1075 ddata[l] = Util::get_max(ddata[l],Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx],sdata[ii + nx + 1], u, v));
1076 }
1077 else if (x2 == (nx-1) && y2 == (ny-1) ) {
1078 ddata[l] = Util::get_max(ddata[l],sdata[ii]);
1079 }
1080 else if (x2 == (nx-1) ) {
1081 ddata[l] = Util::get_max(ddata[l],Util::linear_interpolate(sdata[ii],sdata[ii + nx], v));
1082 }
1083 else if (y2 == (ny-1) ) {
1084 ddata[l] = Util::get_max(ddata[l],Util::linear_interpolate(sdata[ii],sdata[ii + 1], u));
1085 }
1086 }
1087 }
1088 proj->set_attr("xform.projection",t3d);
1089 proj->update();
1090 if(t3d) {delete t3d; t3d=0;}
1091 return proj;
1092 }
1093 else throw ImageDimensionException("Standard projection works only for 2D and 3D images");
1094}
1095
1096
1097// EMData *FourierGriddingProjector::project3d(EMData * image) const
1098// {
1099// if (!image) {
1100// return 0;
1101// }
1102// if (3 != image->get_ndim())
1103// throw ImageDimensionException(
1104// "FourierGriddingProjector needs a 3-D volume");
1105// if (image->is_complex())
1106// throw ImageFormatException(
1107// "FourierGriddingProjector requires a real volume");
1108// const int npad = params.has_key("npad") ? int(params["npad"]) : 2;
1109// const int nx = image->get_xsize();
1110// const int ny = image->get_ysize();
1111// const int nz = image->get_zsize();
1112// if (nx != ny || nx != nz)
1113// throw ImageDimensionException(
1114// "FourierGriddingProjector requires nx==ny==nz");
1115// const int m = Util::get_min(nx,ny,nz);
1116// const int n = m*npad;
1117//
1118// int K = params["kb_K"];
1119// if ( K == 0 ) K = 6;
1120// float alpha = params["kb_alpha"];
1121// if ( alpha == 0 ) alpha = 1.25;
1122// Util::KaiserBessel kb(alpha, K, (float)(m/2), K/(2.0f*n), n);
1123//
1124// // divide out gridding weights
1125// EMData* tmpImage = image->copy();
1126// tmpImage->divkbsinh(kb);
1127// // pad and center volume, then FFT and multiply by (-1)**(i+j+k)
1128// //EMData* imgft = tmpImage->pad_fft(npad);
1129// //imgft->center_padded();
1130// EMData* imgft = tmpImage->norm_pad(false, npad);
1131// imgft->do_fft_inplace();
1132// imgft->center_origin_fft();
1133// imgft->fft_shuffle();
1134// delete tmpImage;
1135//
1136// // Do we have a list of angles?
1137// int nangles = 0;
1138// vector<float> anglelist;
1139// // Do we have a list of angles?
1140// if (params.has_key("anglelist")) {
1141// anglelist = params["anglelist"];
1142// nangles = anglelist.size() / 3;
1143// } else {
1144// // This part was modified by David Woolford -
1145// // Before this the code worked only for SPIDER and EMAN angles,
1146// // but the framework of the Transform3D allows for a generic implementation
1147// // as specified here.
1148// Transform3D* t3d = params["t3d"];
1149// if ( t3d == NULL ) throw NullPointerException("The transform3d object (required for projection), was not specified");
1150// Dict p = t3d->get_rotation(Transform3D::SPIDER);
1151//
1152// string angletype = "SPIDER";
1153// float phi = p["phi"];
1154// float theta = p["theta"];
1155// float psi = p["psi"];
1156// anglelist.push_back(phi);
1157// anglelist.push_back(theta);
1158// anglelist.push_back(psi);
1159// nangles = 1;
1160// }
1161//
1162// // End David Woolford modifications
1163//
1164// // initialize return object
1165// EMData* ret = new EMData();
1166// ret->set_size(nx, ny, nangles);
1167// ret->to_zero();
1168// // loop over sets of angles
1169// for (int ia = 0; ia < nangles; ia++) {
1170// int indx = 3*ia;
1171// Transform3D tf(Transform3D::SPIDER, anglelist[indx],anglelist[indx+1],anglelist[indx+2]);
1172// EMData* proj = imgft->extractplane(tf, kb);
1173// if (proj->is_shuffled()) proj->fft_shuffle();
1174// proj->center_origin_fft();
1175// proj->do_ift_inplace();
1176// EMData* winproj = proj->window_center(m);
1177// delete proj;
1178// for (int iy=0; iy < ny; iy++)
1179// for (int ix=0; ix < nx; ix++)
1180// (*ret)(ix,iy,ia) = (*winproj)(ix,iy);
1181// delete winproj;
1182// }
1183// delete imgft;
1184// ret->update();
1185//
1186// return ret;
1187// }
1188
1189
1190EMData *FourierGriddingProjector::project3d(EMData * image) const
1191{
1192 if (!image) {
1193 return 0;
1194 }
1195 if (3 != image->get_ndim())
1196 throw ImageDimensionException(
1197 "FourierGriddingProjector needs a 3-D volume");
1198 if (image->is_complex())
1199 throw ImageFormatException(
1200 "FourierGriddingProjector requires a real volume");
1201 const int npad = params.has_key("npad") ? int(params["npad"]) : 2;
1202 const int nx = image->get_xsize();
1203 const int ny = image->get_ysize();
1204 const int nz = image->get_zsize();
1205 if (nx != ny || nx != nz)
1206 throw ImageDimensionException(
1207 "FourierGriddingProjector requires nx==ny==nz");
1208 const int m = Util::get_min(nx,ny,nz);
1209 const int n = m*npad;
1210
1211 int K = params["kb_K"];
1212 if ( K == 0 ) K = 6;
1213 float alpha = params["kb_alpha"];
1214 if ( alpha == 0 ) alpha = 1.25;
1215 Util::KaiserBessel kb(alpha, K, (float)(m/2), K/(2.0f*n), n);
1216
1217 // divide out gridding weights
1218 EMData* tmpImage = image->copy();
1219 tmpImage->divkbsinh(kb);
1220 // pad and center volume, then FFT and multiply by (-1)**(i+j+k)
1221 //EMData* imgft = tmpImage->pad_fft(npad);
1222 //imgft->center_padded();
1223 EMData* imgft = tmpImage->norm_pad(false, npad);
1224 imgft->do_fft_inplace();
1225 imgft->center_origin_fft();
1226 imgft->fft_shuffle();
1227 delete tmpImage;
1228
1229 // Do we have a list of angles?
1230 int nangles = 0;
1231 vector<float> anglelist;
1232 // Do we have a list of angles?
1233 if (params.has_key("anglelist")) {
1234 anglelist = params["anglelist"];
1235 nangles = anglelist.size() / 3;
1236 } else {
1237 // This part was modified by David Woolford -
1238 // Before this the code worked only for SPIDER and EMAN angles,
1239 // but the framework of the Transform3D allows for a generic implementation
1240 // as specified here.
1241 Transform* t3d = params["transform"];
1242 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
1243 Dict p = t3d->get_rotation("spider");
1244
1245 string angletype = "SPIDER";
1246 float phi = p["phi"];
1247 float theta = p["theta"];
1248 float psi = p["psi"];
1249 anglelist.push_back(phi);
1250 anglelist.push_back(theta);
1251 anglelist.push_back(psi);
1252 nangles = 1;
1253 if(t3d) {delete t3d; t3d=0;}
1254 }
1255
1256 // End David Woolford modifications
1257
1258 // initialize return object
1259 EMData* ret = new EMData();
1260 ret->set_size(nx, ny, nangles);
1261 ret->to_zero();
1262 // loop over sets of angles
1263 for (int ia = 0; ia < nangles; ia++) {
1264 int indx = 3*ia;
1265 Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]);
1266 Transform tf(d);
1267 EMData* proj = imgft->extract_plane(tf, kb);
1268 if (proj->is_shuffled()) proj->fft_shuffle();
1269 proj->center_origin_fft();
1270 proj->do_ift_inplace();
1271 EMData* winproj = proj->window_center(m);
1272 delete proj;
1273 for (int iy=0; iy < ny; iy++)
1274 for (int ix=0; ix < nx; ix++)
1275 (*ret)(ix,iy,ia) = (*winproj)(ix,iy);
1276 delete winproj;
1277 }
1278 delete imgft;
1279
1280 if (!params.has_key("anglelist")) {
1281 Transform* t3d = params["transform"];
1282 ret->set_attr("xform.projection",t3d);
1283 if(t3d) {delete t3d; t3d=0;}
1284 }
1285 ret->update();
1286 return ret;
1287}
1288
1289// BEGIN Chao projectors and backprojector addition (04/25/06)
1290int ChaoProjector::getnnz(Vec3i volsize, int ri, Vec3i origin, int *nrays, int *nnz) const
1291/*
1292 purpose: count the number of voxels within a sphere centered
1293 at origin and with a radius ri.
1294
1295 input:
1296 volsize contains the size information (nx,ny,nz) about the volume
1297 ri radius of the object embedded in the cube.
1298 origin coordinates for the center of the volume
1299
1300 output:
1301 nnz total number of voxels within the sphere (of radius ri)
1302 nrays number of rays in z-direction.
1303*/
1304{
1305 int ix, iy, iz, rs, r2, xs, ys, zs, xx, yy, zz;
1306 int ftm=0, status = 0;
1307
1308 r2 = ri*ri;
1309 *nnz = 0;
1310 *nrays = 0;
1311 int nx = (int)volsize[0];
1312 int ny = (int)volsize[1];
1313 int nz = (int)volsize[2];
1314
1315 int xcent = (int)origin[0];
1316 int ycent = (int)origin[1];
1317 int zcent = (int)origin[2];
1318
1319 // need to add some error checking
1320 for (ix = 1; ix <=nx; ix++) {
1321 xs = ix-xcent;
1322 xx = xs*xs;
1323 for (iy = 1; iy <= ny; iy++) {
1324 ys = iy-ycent;
1325 yy = ys*ys;
1326 ftm = 1;
1327 for (iz = 1; iz <= nz; iz++) {
1328 zs = iz-zcent;
1329 zz = zs*zs;
1330 rs = xx + yy + zz;
1331 if (rs <= r2) {
1332 (*nnz)++;
1333 if (ftm) {
1334 (*nrays)++;
1335 ftm = 0;
1336 }
1337 }
1338 }
1339 } // end for iy
1340 } // end for ix
1341 return status;
1342}
1343
1344#define cube(i,j,k) cube[ ((k-1)*ny + j-1)*nx + i-1 ]
1345#define sphere(i) sphere[(i)-1]
1346#define cord(i,j) cord[((j)-1)*3 + (i) -1]
1347#define ptrs(i) ptrs[(i)-1]
1348#define dm(i) dm[(i)-1]
1349
1350int ChaoProjector:: cb2sph(float *cube, Vec3i volsize, int ri, Vec3i origin,
1351 int nnz0, int *ptrs, int *cord, float *sphere) const
1352{
1353 int xs, ys, zs, xx, yy, zz, rs, r2;
1354 int ix, iy, iz, jnz, nnz, nrays;
1355 int ftm = 0, status = 0;
1356
1357 int xcent = (int)origin[0];
1358 int ycent = (int)origin[1];
1359 int zcent = (int)origin[2];
1360
1361 int nx = (int)volsize[0];
1362 int ny = (int)volsize[1];
1363 int nz = (int)volsize[2];
1364
1365 r2 = ri*ri;
1366 nnz = 0;
1367 nrays = 0;
1368 ptrs(1) = 1;
1369
1370 for (ix = 1; ix <= nx; ix++) {
1371 xs = ix-xcent;
1372 xx = xs*xs;
1373 for ( iy = 1; iy <= ny; iy++ ) {
1374 ys = iy-ycent;
1375 yy = ys*ys;
1376 jnz = 0;
1377
1378 ftm = 1;
1379 // not the most efficient implementation
1380 for (iz = 1; iz <= nz; iz++) {
1381 zs = iz-zcent;
1382 zz = zs*zs;
1383 rs = xx + yy + zz;
1384 if (rs <= r2) {
1385 jnz++;
1386 nnz++;
1387 sphere(nnz) = cube(iz, iy, ix);
1388
1389 // record the coordinates of the first nonzero ===
1390 if (ftm) {
1391 nrays++;
1392 cord(1,nrays) = iz;
1393 cord(2,nrays) = iy;
1394 cord(3,nrays) = ix;
1395 ftm = 0;
1396 }
1397 }
1398 } // end for (iz..)
1399 if (jnz > 0) {
1400 ptrs(nrays+1) = ptrs(nrays) + jnz;
1401 } // endif (jnz)
1402 } // end for iy
1403 } // end for ix
1404 if (nnz != nnz0) status = -1;
1405 return status;
1406}
1407
1408// decompress sphere into a cube
1409int ChaoProjector::sph2cb(float *sphere, Vec3i volsize, int nrays, int ri,
1410 int nnz0, int *ptrs, int *cord, float *cube) const
1411{
1412 int status=0;
1413 int r2, i, j, ix, iy, iz, nnz;
1414
1415 int nx = (int)volsize[0];
1416 int ny = (int)volsize[1];
1417 // int nz = (int)volsize[2];
1418
1419 r2 = ri*ri;
1420 nnz = 0;
1421 ptrs(1) = 1;
1422
1423 // no need to initialize
1424 // for (i = 0; i<nx*ny*nz; i++) cube[i]=0.0;
1425
1426 nnz = 0;
1427 for (j = 1; j <= nrays; j++) {
1428 iz = cord(1,j);
1429 iy = cord(2,j);
1430 ix = cord(3,j);
1431 for (i = ptrs(j); i<=ptrs(j+1)-1; i++, iz++) {
1432 nnz++;
1433 cube(iz,iy,ix) = sphere(nnz);
1434 }
1435 }
1436 if (nnz != nnz0) status = -1;
1437 return status;
1438}
1439
1440#define x(i) x[(i)-1]
1441#define y(i,j) y[(j-1)*nx + i - 1]
1442
1443// project from 3D to 2D (single image)
1444int ChaoProjector::fwdpj3(Vec3i volsize, int nrays, int , float *dm,
1445 Vec3i origin, int ri, int *ptrs, int *cord,
1446 float *x, float *y) const
1447{
1448 /*
1449 purpose: y <--- proj(x)
1450 input : volsize the size (nx,ny,nz) of the volume
1451 nrays number of rays within the compact spherical
1452 representation
1453 nnz number of voxels within the sphere
1454 dm an array of size 9 storing transformation
1455 associated with the projection direction
1456 origin coordinates of the center of the volume
1457 ri radius of the sphere
1458 ptrs the beginning address of each ray
1459 cord the coordinates of the first point in each ray
1460 x 3d input volume
1461 y 2d output image
1462 */
1463
1464 int iqx, iqy, i, j, xc, yc, zc;
1465 float ct, dipx, dipy, dipx1m, dipy1m, xb, yb, dm1, dm4;
1466 int status = 0;
1467
1468 int xcent = origin[0];
1469 int ycent = origin[1];
1470 int zcent = origin[2];
1471
1472 int nx = volsize[0];
1473
1474 dm1 = dm(1);
1475 dm4 = dm(4);
1476
1477 if ( nx > 2*ri ) {
1478 for (i = 1; i <= nrays; i++) {
1479 zc = cord(1,i)-zcent;
1480 yc = cord(2,i)-ycent;
1481 xc = cord(3,i)-xcent;
1482
1483 xb = zc*dm(1)+yc*dm(2)+xc*dm(3) + xcent;
1484 yb = zc*dm(4)+yc*dm(5)+xc*dm(6) + ycent;
1485
1486 for (j = ptrs(i); j< ptrs(i+1); j++) {
1487 iqx = ifix(xb);
1488 iqy = ifix(yb);
1489
1490 ct = x(j);
1491 dipx = xb - (float)(iqx);
1492 dipy = (yb - (float)(iqy)) * ct;
1493
1494 dipy1m = ct - dipy;
1495 dipx1m = 1.0f - dipx;
1496
1497 y(iqx ,iqy) = y(iqx ,iqy) + dipx1m*dipy1m;
1498 y(iqx+1,iqy) = y(iqx+1,iqy) + dipx*dipy1m;
1499 y(iqx+1,iqy+1) = y(iqx+1,iqy+1) + dipx*dipy;
1500 y(iqx ,iqy+1) = y(iqx ,iqy+1) + dipx1m*dipy;
1501
1502 xb += dm1;
1503 yb += dm4;
1504 }
1505 }
1506 }
1507 else {
1508 fprintf(stderr, " nx must be greater than 2*ri\n");
1509 exit(1);
1510 }
1511 return status;
1512}
1513#undef x
1514#undef y
1515
1516#define y(i) y[(i)-1]
1517#define x(i,j) x[((j)-1)*nx + (i) - 1]
1518
1519// backproject from 2D to 3D for a single image
1520int ChaoProjector::bckpj3(Vec3i volsize, int nrays, int , float *dm,
1521 Vec3i origin, int ri, int *ptrs, int *cord,
1522 float *x, float *y) const
1523{
1524 int i, j, iqx,iqy, xc, yc, zc;
1525 float xb, yb, dx, dy, dx1m, dy1m, dxdy;
1526 int status = 0;
1527
1528 int xcent = origin[0];
1529 int ycent = origin[1];
1530 int zcent = origin[2];
1531
1532 int nx = volsize[0];
1533
1534 if ( nx > 2*ri) {
1535 for (i = 1; i <= nrays; i++) {
1536 zc = cord(1,i) - zcent;
1537 yc = cord(2,i) - ycent;
1538 xc = cord(3,i) - xcent;
1539
1540 xb = zc*dm(1)+yc*dm(2)+xc*dm(3) + xcent;
1541 yb = zc*dm(4)+yc*dm(5)+xc*dm(6) + ycent;
1542
1543 for (j = ptrs(i); j <ptrs(i+1); j++) {
1544 iqx = ifix((float)(xb));
1545 iqy = ifix((float)(yb));
1546
1547 dx = xb - (float)(iqx);
1548 dy = yb - (float)(iqy);
1549 dx1m = 1.0f - dx;
1550 dy1m = 1.0f - dy;
1551 dxdy = dx*dy;
1552/*
1553c y(j) = y(j) + dx1m*dy1m*x(iqx , iqy)
1554c & + dx1m*dy *x(iqx , iqy+1)
1555c & + dx *dy1m*x(iqx+1, iqy)
1556c & + dx *dy *x(iqx+1, iqy+1)
1557c
1558c --- faster version of the above commented out
1559c code (derived by summing the following table
1560c of coefficients along the colunms) ---
1561c
1562c 1 dx dy dxdy
1563c ------ -------- -------- -------
1564c x(i,j) -x(i,j) -x(i,j) x(i,j)
1565c x(i,j+1) -x(i,j+1)
1566c x(i+1,j) -x(i+1,j)
1567c x(i+1,j+1)
1568c
1569*/
1570 y(j) += x(iqx,iqy)
1571 + dx*(-x(iqx,iqy)+x(iqx+1,iqy))
1572 + dy*(-x(iqx,iqy)+x(iqx,iqy+1))
1573 + dxdy*( x(iqx,iqy) - x(iqx,iqy+1)
1574 -x(iqx+1,iqy) + x(iqx+1,iqy+1) );
1575
1576 xb += dm(1);
1577 yb += dm(4);
1578 } // end for j
1579 } // end for i
1580 }
1581 else {
1582 fprintf(stderr, "bckpj3: nx must be greater than 2*ri\n");
1583 }
1584
1585 return status;
1586}
1587
1588#undef x
1589#undef y
1590#undef dm
1591
1592// funny F90 style strange rounding function
1593int ChaoProjector::ifix(float a) const
1594{
1595 int ia;
1596
1597 if (a>=0) {
1598 ia = (int)floor(a);
1599 }
1600 else {
1601 ia = (int)ceil(a);
1602 }
1603 return ia;
1604}
1605
1606#define dm(i,j) dm[((j)-1)*9 + (i) -1]
1607#define anglelist(i,j) anglelist[((j)-1)*3 + (i) - 1]
1608
1609// SPIDER stype transformation
1610void ChaoProjector::setdm(vector<float> anglelist, string const , float *dm) const
1611{ // convert Euler angles to transformations, dm is an 9 by nangles array
1612
1613 float psi, theta, phi;
1614 double cthe, sthe, cpsi, spsi, cphi, sphi;
1615 int j;
1616
1617 int nangles = anglelist.size() / 3;
1618
1619 // now convert all angles
1620 for (j = 1; j <= nangles; j++) {
1621 phi = static_cast<float>(anglelist(1,j)*dgr_to_rad);
1622 theta = static_cast<float>(anglelist(2,j)*dgr_to_rad);
1623 psi = static_cast<float>(anglelist(3,j)*dgr_to_rad);
1624
1625 // cout << phi << " " << theta << " " << psi << endl;
1626 cthe = cos(theta);
1627 sthe = sin(theta);
1628 cpsi = cos(psi);
1629 spsi = sin(psi);
1630 cphi = cos(phi);
1631 sphi = sin(phi);
1632
1633 dm(1,j)=static_cast<float>(cphi*cthe*cpsi-sphi*spsi);
1634 dm(2,j)=static_cast<float>(sphi*cthe*cpsi+cphi*spsi);
1635 dm(3,j)=static_cast<float>(-sthe*cpsi);
1636 dm(4,j)=static_cast<float>(-cphi*cthe*spsi-sphi*cpsi);
1637 dm(5,j)=static_cast<float>(-sphi*cthe*spsi+cphi*cpsi);
1638 dm(6,j)=static_cast<float>(sthe*spsi);
1639 dm(7,j)=static_cast<float>(sthe*cphi);
1640 dm(8,j)=static_cast<float>(sthe*sphi);
1641 dm(9,j)=static_cast<float>(cthe);
1642 }
1643}
1644#undef anglelist
1645
1646#define images(i,j,k) images[ ((k-1)*nyvol + j-1)*nxvol + i-1 ]
1647
1648EMData *ChaoProjector::project3d(EMData * vol) const
1649{
1650
1651 int nrays, nnz, status, j;
1652 float *dm;
1653 int *ptrs, *cord;
1654 float *sphere, *images;
1655
1656 int nxvol = vol->get_xsize();
1657 int nyvol = vol->get_ysize();
1658 int nzvol = vol->get_zsize();
1659 Vec3i volsize(nxvol,nyvol,nzvol);
1660
1661 int dim = Util::get_min(nxvol,nyvol,nzvol);
1662 if (nzvol == 1) {
1663 LOGERR("The ChaoProjector needs a volume!");
1664 return 0;
1665 }
1666 Vec3i origin(0,0,0);
1667 // If a sensible origin isn't passed in, choose the middle of
1668 // the cube.
1669 if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
1670 else {origin[0] = nxvol/2+1;}
1671 if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
1672 else {origin[1] = nyvol/2+1;}
1673 if (params.has_key("origin_z")) {origin[2] = params["origin_z"];}
1674 else {origin[2] = nzvol/2+1;}
1675
1676 int ri;
1677 if (params.has_key("radius")) {ri = params["radius"];}
1678 else {ri = dim/2 - 1;}
1679
1680 // retrieve the voxel values
1681 float *cube = vol->get_data();
1682
1683 // count the number of voxels within a sphere centered at icent,
1684 // with radius ri
1685 status = getnnz(volsize, ri, origin, &nrays, &nnz);
1686 // need to check status...
1687
1688 // convert from cube to sphere
1689 sphere = new float[nnz];
1690 ptrs = new int[nrays+1];
1691 cord = new int[3*nrays];
1692 if (sphere == NULL || ptrs == NULL || cord == NULL) {
1693 fprintf(stderr,"ChaoProjector::project3d, failed to allocate!\n");
1694 exit(1);
1695 }
1696 for (int i = 0; i<nnz; i++) sphere[i] = 0.0;
1697 for (int i = 0; i<nrays+1; i++) ptrs[i] = 0;
1698 for (int i = 0; i<3*nrays; i++) cord[i] = 0;
1699
1700 status = cb2sph(cube, volsize, ri, origin, nnz, ptrs, cord, sphere);
1701 // check status
1702
1703 int nangles = 0;
1704 vector<float> anglelist;
1705 string angletype = "SPIDER";
1706 // Do we have a list of angles?
1707 if (params.has_key("anglelist")) {
1708 anglelist = params["anglelist"];
1709 nangles = anglelist.size() / 3;
1710 } else {
1711 Transform* t3d = params["transform"];
1712 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
1713 // This part was modified by David Woolford -
1714 // Before this the code worked only for SPIDER and EMAN angles,
1715 // but the framework of the Transform3D allows for a generic implementation
1716 // as specified here.
1717 Dict p = t3d->get_rotation("spider");
1718 if(t3d) {delete t3d; t3d=0;}
1719
1720 float phi = p["phi"];
1721 float theta = p["theta"];
1722 float psi = p["psi"];
1723 anglelist.push_back(phi);
1724 anglelist.push_back(theta);
1725 anglelist.push_back(psi);
1726 nangles = 1;
1727 }
1728 // End David Woolford modifications
1729
1730 dm = new float[nangles*9];
1731 setdm(anglelist, angletype, dm);
1732
1733 // return images
1734 EMData *ret = new EMData();
1735 ret->set_size(nxvol, nyvol, nangles);
1736 ret->set_complex(false);
1737 ret->set_ri(true);
1738
1739 images = ret->get_data();
1740
1741 for (j = 1; j <= nangles; j++) {
1742 status = fwdpj3(volsize, nrays, nnz , &dm(1,j), origin, ri,
1743 ptrs , cord, sphere, &images(1,1,j));
1744 // check status?
1745 }
1746
1747 // deallocate all temporary work space
1748 EMDeleteArray(dm);
1749 EMDeleteArray(ptrs);
1750 EMDeleteArray(cord);
1751 EMDeleteArray(sphere);
1752
1753 if (!params.has_key("anglelist")) {
1754 Transform* t3d = params["transform"];
1755 ret->set_attr("xform.projection",t3d);
1756 if(t3d) {delete t3d; t3d=0;}
1757 }
1758 ret->update();
1759 return ret;
1760}
1761
1762
1763#undef images
1764
1765#define images(i,j,k) images[ ((k)-1)*nximg*nyimg + ((j)-1)*nximg + (i)-1 ]
1766// backproject from 2D to 3D (multiple images)
1767EMData *ChaoProjector::backproject3d(EMData * imagestack) const
1768{
1769 int nrays, nnz, status, j;
1770 float *dm;
1771 int *ptrs, *cord;
1772 float *sphere, *images, *cube;
1773
1774 int nximg = imagestack->get_xsize();
1775 int nyimg = imagestack->get_ysize();
1776 int nslices = imagestack->get_zsize();
1777
1778 int dim = Util::get_min(nximg,nyimg);
1779 Vec3i volsize(nximg,nyimg,dim);
1780
1781 Vec3i origin(0,0,0);
1782 // If a sensible origin isn't passed in, choose the middle of
1783 // the cube.
1784 if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
1785 else {origin[0] = nximg/2+1;}
1786 if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
1787 else {origin[1] = nyimg/2+1;}
1788 if (params.has_key("origin_z")) {origin[1] = params["origin_z"];}
1789 else {origin[2] = dim/2+1;}
1790
1791 int ri;
1792 if (params.has_key("radius")) {ri = params["radius"];}
1793 else {ri = dim/2 - 1;}
1794
1795 // retrieve the voxel values
1796 images = imagestack->get_data();
1797
1798 // count the number of voxels within a sphere centered at icent,
1799 // with radius ri
1800 status = getnnz(volsize, ri, origin, &nrays, &nnz);
1801 // need to check status...
1802
1803 // convert from cube to sphere
1804 sphere = new float[nnz];
1805 ptrs = new int[nrays+1];
1806 cord = new int[3*nrays];
1807 if (sphere == NULL || ptrs == NULL || cord == NULL) {
1808 fprintf(stderr,"ChaoProjector::backproject3d, failed to allocate!\n");
1809 exit(1);
1810 }
1811 for (int i = 0; i<nnz; i++) sphere[i] = 0.0;
1812 for (int i = 0; i<nrays+1; i++) ptrs[i] = 0;
1813 for (int i = 0; i<3*nrays; i++) cord[i] = 0;
1814
1815 int nangles = 0;
1816 vector<float> anglelist;
1817 string angletype = "SPIDER";
1818 // Do we have a list of angles?
1819 if (params.has_key("anglelist")) {
1820 anglelist = params["anglelist"];
1821 nangles = anglelist.size() / 3;
1822 } else {
1823 Transform* t3d = params["transform"];
1824 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
1825 // This part was modified by David Woolford -
1826 // Before this the code worked only for SPIDER and EMAN angles,
1827 // but the framework of the Transform3D allows for a generic implementation
1828 // as specified here.
1829 // This was broken by david. we need here a loop over all projections and put all angles on stack PAP 06/28/09
1830 Dict p = t3d->get_rotation("spider");
1831 if(t3d) {delete t3d; t3d=0;}
1832
1833 float phi = p["phi"];
1834 float theta = p["theta"];
1835 float psi = p["psi"];
1836 anglelist.push_back(phi);
1837 anglelist.push_back(theta);
1838 anglelist.push_back(psi);
1839 nangles = 1;
1840 }
1841
1842 // End David Woolford modifications
1843
1844 if (nslices != nangles) {
1845 LOGERR("the number of images does not match the number of angles");
1846 return 0;
1847 }
1848
1849 dm = new float[nangles*9];
1850 setdm(anglelist, angletype, dm);
1851
1852 // return volume
1853 EMData *ret = new EMData();
1854 ret->set_size(nximg, nyimg, dim);
1855 ret->set_complex(false);
1856 ret->set_ri(true);
1857 ret->to_zero();
1858
1859 cube = ret->get_data();
1860 // cb2sph should be replaced by something that touches only ptrs and cord
1861 status = cb2sph(cube, volsize, ri, origin, nnz, ptrs, cord, sphere);
1862 // check status
1863
1864 for (j = 1; j <= nangles; j++) {
1865 status = bckpj3(volsize, nrays, nnz, &dm(1,j), origin, ri,
1866 ptrs , cord , &images(1,1,j), sphere);
1867 // check status?
1868 }
1869
1870 status = sph2cb(sphere, volsize, nrays, ri, nnz, ptrs, cord, cube);
1871 // check status?
1872
1873 // deallocate all temporary work space
1874 EMDeleteArray(dm);
1875 EMDeleteArray(ptrs);
1876 EMDeleteArray(cord);
1877 EMDeleteArray(sphere);
1878
1879 ret->update();
1880 return ret;
1881}
1882
1883#undef images
1884#undef cube
1885#undef sphere
1886#undef cord
1887#undef ptrs
1888#undef dm
1889
1890EMData *GaussFFTProjector::backproject3d(EMData * ) const
1891{
1892 // no implementation yet
1893 EMData *ret = new EMData();
1894 return ret;
1895}
1896
1897#define images(i,j,k) images[ (k)*nx*ny + ((j)-1)*nx + (i)-1 ]
1898
1899// EMData *PawelProjector::backproject3d(EMData * imagestack) const
1900// {
1901//
1902// float *images;
1903//
1904// if (!imagestack) {
1905// return 0;
1906// }
1907// int ri;
1908// int nx = imagestack->get_xsize();
1909// int ny = imagestack->get_ysize();
1910// // int nslices = imagestack->get_zsize();
1911// int dim = Util::get_min(nx,ny);
1912// images = imagestack->get_data();
1913//
1914// Vec3i origin(0,0,0);
1915// // If a sensible origin isn't passed in, choose the middle of
1916// // the cube.
1917// if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
1918// else {origin[0] = nx/2;}
1919// if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
1920// else {origin[1] = ny/2;}
1921// if (params.has_key("origin_z")) {origin[1] = params["origin_z"];}
1922// else {origin[2] = dim/2;}
1923//
1924// if (params.has_key("radius")) {ri = params["radius"];}
1925// else {ri = dim/2 - 1;}
1926//
1927// // Determine the number of rows (x-lines) within the radius
1928// int nn = -1;
1929// prepcubes(nx, ny, dim, ri, origin, nn);
1930// // nn is now the number of rows-1 within the radius
1931// // so we can create and fill the ipcubes
1932// IPCube* ipcube = new IPCube[nn+1];
1933// prepcubes(nx, ny, dim, ri, origin, nn, ipcube);
1934//
1935// int nangles = 0;
1936// vector<float> anglelist;
1937// // Do we have a list of angles?
1938// if (params.has_key("anglelist")) {
1939// anglelist = params["anglelist"];
1940// nangles = anglelist.size() / 3;
1941// } else {
1942// Transform3D* t3d = params["t3d"];
1943// if ( t3d == NULL ) throw NullPointerException("The transform3d object (required for projection), was not specified");
1944// // This part was modified by David Woolford -
1945// // Before this the code worked only for SPIDER and EMAN angles,
1946// // but the framework of the Transform3D allows for a generic implementation
1947// // as specified here.
1948// Dict p = t3d->get_rotation(Transform3D::SPIDER);
1949//
1950// string angletype = "SPIDER";
1951// float phi = p["phi"];
1952// float theta = p["theta"];
1953// float psi = p["psi"];
1954// anglelist.push_back(phi);
1955// anglelist.push_back(theta);
1956// anglelist.push_back(psi);
1957// nangles = 1;
1958// }
1959//
1960// // End David Woolford modifications
1961//
1962// // initialize return object
1963// EMData* ret = new EMData();
1964// ret->set_size(nx, ny, dim);
1965// ret->to_zero();
1966//
1967// // loop over sets of angles
1968// for (int ia = 0; ia < nangles; ia++) {
1969// int indx = 3*ia;
1970// Transform3D rotation(Transform3D::SPIDER, float(anglelist[indx]),
1971// float(anglelist[indx+1]),
1972// float(anglelist[indx+2]));
1973// float dm1 = rotation.at(0,0);
1974// float dm4 = rotation.at(1,0);
1975//
1976// if (2*(ri+1)+1 > dim) {
1977// // Must check x and y boundaries
1978// LOGERR("backproject3d, pawel, 2*(ri+1)+1 > dim\n");
1979// return 0;
1980// } else {
1981// // No need to check x and y boundaries
1982// for (int i = 0 ; i <= nn; i++) {
1983// int iox = (int)ipcube[i].loc[0]+origin[0];
1984// int ioy = (int)ipcube[i].loc[1]+origin[1];
1985// int ioz = (int)ipcube[i].loc[2]+origin[2];
1986//
1987// Vec3f vb = rotation*ipcube[i].loc + origin;
1988// for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
1989// float xbb = (j-ipcube[i].start)*dm1 + vb[0];
1990// int iqx = (int)floor(xbb);
1991//
1992// float ybb = (j-ipcube[i].start)*dm4 + vb[1];
1993// int iqy = (int)floor(ybb);
1994//
1995// float dipx = xbb - iqx;
1996// float dipy = ybb - iqy;
1997//
1998// (*ret)(iox,ioy,ioz) += images(iqx,iqy,ia)
1999// + dipy*(images(iqx,iqy+1,ia)-images(iqx,iqy,ia))
2000// + dipx*(images(iqx+1,iqy,ia)-images(iqx,iqy,ia)
2001// + dipy*(images(iqx+1,iqy+1,ia)-images(iqx+1,iqy,ia)
2002// - images(iqx,iqy+1,ia)+images(iqx,iqy,ia)));
2003// iox++;
2004// } // end for j
2005// } // end for i
2006// } // end if
2007// } // end for ia
2008//
2009// ret->update();
2010// EMDeleteArray(ipcube);
2011// return ret;
2012// }
2013EMData *PawelProjector::backproject3d(EMData * imagestack) const
2014{
2015
2016 float *images;
2017
2018 if (!imagestack) {
2019 return 0;
2020 }
2021 int ri;
2022 int nx = imagestack->get_xsize();
2023 int ny = imagestack->get_ysize();
2024// int nslices = imagestack->get_zsize();
2025 int dim = Util::get_min(nx,ny);
2026 images = imagestack->get_data();
2027
2028 Vec3i origin(0,0,0);
2029 // If a sensible origin isn't passed in, choose the middle of
2030 // the cube.
2031 if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
2032 else {origin[0] = nx/2;}
2033 if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
2034 else {origin[1] = ny/2;}
2035 if (params.has_key("origin_z")) {origin[1] = params["origin_z"];}
2036 else {origin[2] = dim/2;}
2037
2038 if (params.has_key("radius")) {ri = params["radius"];}
2039 else {ri = dim/2 - 1;}
2040
2041 // Determine the number of rows (x-lines) within the radius
2042 int nn = -1;
2043 prepcubes(nx, ny, dim, ri, origin, nn);
2044 // nn is now the number of rows-1 within the radius
2045 // so we can create and fill the ipcubes
2046 IPCube* ipcube = new IPCube[nn+1];
2047 prepcubes(nx, ny, dim, ri, origin, nn, ipcube);
2048
2049 int nangles = 0;
2050 vector<float> anglelist;
2051 // Do we have a list of angles?
2052 if (params.has_key("anglelist")) {
2053 anglelist = params["anglelist"];
2054 nangles = anglelist.size() / 3;
2055 } else {
2056 Transform* t3d = params["transform"];
2057 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
2058 // This part was modified by David Woolford -
2059 // Before this the code worked only for SPIDER and EMAN angles,
2060 // but the framework of the Transform3D allows for a generic implementation
2061 // as specified here.
2062 Dict p = t3d->get_rotation("spider");
2063 if(t3d) {delete t3d; t3d=0;}
2064
2065 string angletype = "SPIDER";
2066 float phi = p["phi"];
2067 float theta = p["theta"];
2068 float psi = p["psi"];
2069 anglelist.push_back(phi);
2070 anglelist.push_back(theta);
2071 anglelist.push_back(psi);
2072 nangles = 1;
2073 }
2074
2075 // End David Woolford modifications
2076
2077 // initialize return object
2078 EMData* ret = new EMData();
2079 ret->set_size(nx, ny, dim);
2080 ret->to_zero();
2081
2082 // loop over sets of angles
2083 for (int ia = 0; ia < nangles; ia++) {
2084 int indx = 3*ia;
2085 Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]);
2086 Transform rotation(d);
2087 float dm1 = rotation.at(0,0);
2088 float dm4 = rotation.at(1,0);
2089
2090 if (2*(ri+1)+1 > dim) {
2091 // Must check x and y boundaries
2092 LOGERR("backproject3d, pawel, 2*(ri+1)+1 > dim\n");
2093 return 0;
2094 } else {
2095 // No need to check x and y boundaries
2096 for (int i = 0 ; i <= nn; i++) {
2097 int iox = (int)ipcube[i].loc[0]+origin[0];
2098 int ioy = (int)ipcube[i].loc[1]+origin[1];
2099 int ioz = (int)ipcube[i].loc[2]+origin[2];
2100
2101 Vec3f vb = rotation*ipcube[i].loc + origin;
2102 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
2103 float xbb = (j-ipcube[i].start)*dm1 + vb[0];
2104 int iqx = (int)floor(xbb);
2105
2106 float ybb = (j-ipcube[i].start)*dm4 + vb[1];
2107 int iqy = (int)floor(ybb);
2108
2109 float dipx = xbb - iqx;
2110 float dipy = ybb - iqy;
2111
2112 (*ret)(iox,ioy,ioz) += images(iqx,iqy,ia)
2113 + dipy*(images(iqx,iqy+1,ia)-images(iqx,iqy,ia))
2114 + dipx*(images(iqx+1,iqy,ia)-images(iqx,iqy,ia)
2115 + dipy*(images(iqx+1,iqy+1,ia)-images(iqx+1,iqy,ia)
2116 - images(iqx,iqy+1,ia)+images(iqx,iqy,ia)));
2117 iox++;
2118 } // end for j
2119 } // end for i
2120 } // end if
2121 } // end for ia
2122
2123 ret->update();
2124 EMDeleteArray(ipcube);
2125 return ret;
2126}
2127#undef images
2128
2129EMData *StandardProjector::backproject3d(EMData * ) const
2130{
2131 // no implementation yet
2132 EMData *ret = new EMData();
2133 return ret;
2134}
2135
2136EMData *MaxValProjector::backproject3d(EMData * ) const
2137{
2138 // no implementation yet
2139 EMData *ret = new EMData();
2140 return ret;
2141}
2142
2143
2144EMData *FourierGriddingProjector::backproject3d(EMData * ) const
2145{
2146 // no implementation yet
2147 EMData *ret = new EMData();
2148 return ret;
2149}
2150
2151// End Chao's projector addition 4/25/06
2152
2153void EMAN::dump_projectors()
2154{
2155 dump_factory < Projector > ();
2156}
2157
2158map<string, vector<string> > EMAN::dump_projectors_list()
2159{
2160 return dump_factory_list < Projector > ();
2161}
rdata
#define rdata(i)
Definition: analyzer.cpp:592
EMAN::ChaoProjector::cb2sph
int cb2sph(float *cube, Vec3i volsize, int ri, Vec3i origin, int nnz0, int *ptrs, int *cord, float *sphere) const
Definition: projector.cpp:1350
EMAN::ChaoProjector::backproject3d
EMData * backproject3d(EMData *imagestack) const
Back-project a 2D image into a 3D image.
Definition: projector.cpp:1767
EMAN::ChaoProjector::setdm
void setdm(vector< float > anglelist, string const angletype, float *dm) const
Definition: projector.cpp:1610
EMAN::ChaoProjector::sph2cb
int sph2cb(float *sphere, Vec3i volsize, int nray, int ri, int nnz0, int *ptrs, int *cord, float *cube) const
Definition: projector.cpp:1409
EMAN::ChaoProjector::ifix
int ifix(float a) const
Definition: projector.cpp:1593
EMAN::ChaoProjector::project3d
EMData * project3d(EMData *vol) const
Project an 3D image into a 2D image.
Definition: projector.cpp:1648
EMAN::ChaoProjector::fwdpj3
int fwdpj3(Vec3i volsize, int nray, int nnz, float *dm, Vec3i origin, int ri, int *ptrs, int *cord, float *x, float *y) const
Definition: projector.cpp:1444
EMAN::ChaoProjector::bckpj3
int bckpj3(Vec3i volsize, int nray, int nnz, float *dm, Vec3i origin, int ri, int *ptrs, int *cord, float *x, float *y) const
Definition: projector.cpp:1520
EMAN::ChaoProjector::getnnz
int getnnz(Vec3i volsize, int ri, Vec3i origin, int *nray, int *nnz) const
Definition: projector.cpp:1290
EMAN::Dict
Dict is a dictionary to store <string, EMObject> pair.
Definition: emobject.h:385
EMAN::Dict::has_key
bool has_key(const string &key) const
Ask the Dictionary if it as a particular key.
Definition: emobject.h:511
EMAN::EMConsts::I5G
static const float I5G
Definition: emobject.h:75
EMAN::EMConsts::I4G
static const float I4G
Definition: emobject.h:74
EMAN::EMConsts::I3G
static const float I3G
Definition: emobject.h:73
EMAN::EMConsts::I2G
static const float I2G
Definition: emobject.h:72
EMAN::EMData
EMData stores an image's data and defines core image processing routines.
Definition: emdata.h:82
EMAN::EMData::translate
void translate(float dx, float dy, float dz)
Translate this image.
Definition: emdata.cpp:904
EMAN::EMData::window_center
EMData * window_center(int l)
Window the center of an image.
Definition: emdata.cpp:776
EMAN::EMData::setup4slice
float * setup4slice(bool redo=true)
Set up for fftslice operations.
Definition: emdata.cpp:822
EMAN::Factory
Factory is used to store objects to create new instances.
Definition: emobject.h:725
EMAN::FourierGriddingProjector::project3d
EMData * project3d(EMData *image) const
Project an 3D image into a 2D image.
Definition: projector.cpp:1190
EMAN::FourierGriddingProjector::backproject3d
EMData * backproject3d(EMData *image) const
Back-project a 2D image into a 3D image.
Definition: projector.cpp:2144
EMAN::GaussFFTProjector::interp_ft_3d
bool interp_ft_3d(int mode, EMData *image, float x, float y, float z, float *data, float gauss_width) const
Definition: projector.cpp:216
EMAN::GaussFFTProjector::backproject3d
EMData * backproject3d(EMData *image) const
Back-project a 2D image into a 3D image.
Definition: projector.cpp:1890
EMAN::Interp::hyperg
static float hyperg(float v)
Definition: interp.h:51
EMAN::Interp::get_gimx
static float * get_gimx()
Definition: interp.h:62
EMAN::MaxValProjector::project3d
EMData * project3d(EMData *image) const
Project an 3D image into a 2D image.
Definition: projector.cpp:944
EMAN::MaxValProjector::backproject3d
EMData * backproject3d(EMData *image) const
Back-project a 2D image into a 3D image.
Definition: projector.cpp:2136
EMAN::PawelProjector::backproject3d
EMData * backproject3d(EMData *image) const
Back-project a 2D image into a 3D image.
Definition: projector.cpp:2013
EMAN::PawelProjector::project3d
EMData * project3d(EMData *image) const
Project an 3D image into a 2D image.
Definition: projector.cpp:602
EMAN::PawelProjector::prepcubes
void prepcubes(int nx, int ny, int nz, int ri, Vec3i origin, int &nn, IPCube *ipcube=NULL) const
Definition: projector.cpp:559
EMAN::StandardProjector::project3d
EMData * project3d(EMData *image) const
Project an 3D image into a 2D image.
Definition: projector.cpp:765
EMAN::StandardProjector::backproject3d
EMData * backproject3d(EMData *image) const
Back-project a 2D image into a 3D image.
Definition: projector.cpp:2129
EMAN::Transform
A Transform object is a somewhat specialized object designed specifically for EMAN2/Sparx storage of ...
Definition: transform.h:75
EMAN::Transform::at
float at(int r, int c) const
Get the value stored in the internal transformation matrix at at coordinate (r,c)
Definition: transform.h:396
EMAN::Transform::get_mirror
bool get_mirror() const
Query whether x_mirroring is occurring.
Definition: transform.cpp:1250
EMAN::Transform::copy_matrix_into_array
void copy_matrix_into_array(float *const) const
Definition: transform.cpp:158
EMAN::Transform::get_matrix3_row
Vec3f get_matrix3_row(int i) const
Get a matrix row as a Vec3f required for back compatibility with Tranform3D - see PawelProjector.
Definition: transform.h:465
EMAN::Transform::get_rotation_transform
Transform get_rotation_transform() const
Get the rotation part of the tranformation matrix as a Transform object.
Definition: transform.cpp:688
EMAN::Transform::get_rotation
Dict get_rotation(const string &euler_type="eman") const
Get a rotation in any Euler format.
Definition: transform.cpp:829
EMAN::Transform::get_trans
Vec3f get_trans() const
Get the post trans as a vec3f.
Definition: transform.cpp:1046
EMAN::Transform::get_scale
float get_scale() const
Get the scale that was applied.
Definition: transform.cpp:1145
EMAN::Transform::inverse
Transform inverse() const
Get the inverse of this transformation matrix.
Definition: transform.cpp:1327
EMAN::Transform::invert
void invert()
Get the inverse of this transformation matrix.
Definition: transform.cpp:1293
EMAN::Util::hypot3sq
static int hypot3sq(int x, int y, int z)
Euclidean distance function squared in 3D: f(x,y,z) = (x*x + y*y + z*z);.
Definition: util.h:805
EMAN::Util::agauss
static float agauss(float a, float dx, float dy, float dz, float d)
Calculate Gaussian value.
Definition: util.h:912
EMAN::Util::linear_interpolate
static float linear_interpolate(float p1, float p2, float t)
Calculate linear interpolation.
Definition: util.h:518
EMAN::Util::fast_floor
static int fast_floor(float x)
A fast way to calculate a floor, which is largest integral value not greater than argument.
Definition: util.h:874
EMAN::Util::trilinear_interpolate
static float trilinear_interpolate(float p1, float p2, float p3, float p4, float p5, float p6, float p7, float p8, float t, float u, float v)
Calculate trilinear interpolation.
Definition: util.h:619
EMAN::Util::get_min
static int get_min(int f1, int f2)
Get the minimum of 2 numbers.
Definition: util.h:922
EMAN::Util::bilinear_interpolate
static float bilinear_interpolate(float p1, float p2, float p3, float p4, float t, float u)
Calculate bilinear interpolation.
Definition: util.h:543
EMAN::Util::get_max
static float get_max(float f1, float f2)
Get the maximum of 2 numbers.
Definition: util.h:998
EMAN::Vec2
The Vec2 is precisely the same as Vec3 except it works exclusively in 2D Note there are convenient ty...
Definition: vec3.h:710
cuda_projector.h
standard_project
void standard_project(const float *const matrix, float *data, const int nx, const int ny)
cuda_util.h
sqrt
EMData * sqrt() const
return square root of current image
Definition: emdata_core.cpp:1084
EMDeleteArray
void EMDeleteArray(T &x)
Definition: emutil.h:62
ImageFormatException
#define ImageFormatException(desc)
Definition: exception.h:147
ImageDimensionException
#define ImageDimensionException(desc)
Definition: exception.h:166
NullPointerException
#define NullPointerException(desc)
Definition: exception.h:241
LOGERR
#define LOGERR
Definition: log.h:51
EMAN
E2Exception class.
Definition: aligner.h:40
EMAN::dump_projectors_list
map< string, vector< string > > dump_projectors_list()
Definition: projector.cpp:2158
EMAN::dump_projectors
void dump_projectors()
Definition: projector.cpp:2153
y
#define y(i, j)
Definition: projector.cpp:1516
cube
#define cube(i, j, k)
Definition: projector.cpp:1344
ptrs
#define ptrs(i)
Definition: projector.cpp:1347
sphere
#define sphere(i)
Definition: projector.cpp:1345
cord
#define cord(i, j)
Definition: projector.cpp:1346
x
#define x(i)
Definition: projector.cpp:1517
anglelist
#define anglelist(i, j)
Definition: projector.cpp:1607
images
#define images(i, j, k)
Definition: projector.cpp:1897
dm
#define dm(i)
Definition: projector.cpp:1606
projector_template.h
Generated by doxygen 1.9.2