#include <projector.h>
Inheritance diagram for EMAN::GaussFFTProjector:


Public Member Functions | |
| GaussFFTProjector () | |
| EMData * | project3d (EMData *image) const |
| Project an 3D image into a 2D image. | |
| EMData * | backproject3d (EMData *image) const |
| Back-project a 2D image into a 3D image. | |
| void | set_params (const Dict &new_params) |
| Set the projector parameters using a key/value dictionary. | |
| string | get_name () const |
| Get the projector's name. | |
| string | get_desc () const |
| TypeDict | get_param_types () const |
| Get processor parameter information in a dictionary. | |
Static Public Member Functions | |
| Projector * | NEW () |
Static Public Attributes | |
| const string | NAME = "gauss_fft" |
Private Member Functions | |
| bool | interp_ft_3d (int mode, EMData *image, float x, float y, float z, float *data, float gauss_width) const |
Private Attributes | |
| float | alt |
| float | az |
| float | phi |
use integer 'mode' to determine the gaussian width and the way to interpolate a point in a 3d complex image. valid mode range: [1,7]. the gauss widths are as follows with mode from 1 to 7:
mode 1: 0; mode 2: 4.0 / (M_PI * M_PI); mode 3: 6.4 / (M_PI * M_PI); mode 4: 8.8 / (M_PI * M_PI); mode 5: 0; mode 6: 10.4 / (M_PI * M_PI); mode 7: 10.4 / (M_PI * M_PI);
Definition at line 150 of file projector.h.
|
|
Definition at line 153 of file projector.h. References phi.
|
|
|
Back-project a 2D image into a 3D image.
Implements EMAN::Projector. Definition at line 1844 of file projector.cpp.
|
|
|
Implements EMAN::Projector. Definition at line 175 of file projector.h. 00176 {
00177 return "Projections using a Gaussian kernel in Fourier space. Produces artifacts, not recommended.";
00178 }
|
|
|
Get the projector's name. Each projector is indentified by unique name.
Implements EMAN::Projector. Definition at line 170 of file projector.h. 00171 {
00172 return NAME;
00173 }
|
|
|
Get processor parameter information in a dictionary. Each parameter has one record in the dictionary. Each record contains its name, data-type, and description.
Reimplemented from EMAN::Projector. Definition at line 185 of file projector.h. References EMAN::TypeDict::put(). 00186 {
00187 TypeDict d;
00188 d.put("transform", EMObject::TRANSFORM);
00189 d.put("mode", EMObject::INT);
00190 return d;
00191 }
|
|
||||||||||||||||||||||||||||||||
|
Definition at line 206 of file projector.cpp. References abs, EMAN::Util::agauss(), data, EMAN::EMData::get_data(), EMAN::Interp::get_gimx(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::Interp::hyperg(), EMAN::Util::hypot3sq(), nx, ny, rdata, EMAN::EMData::setup4slice(), sqrt(), x, and y. Referenced by project3d(). 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 // throw InvalidParameterException("Error, unkown mode");
00528 return false;
00529 }
|
|
|
Definition at line 180 of file projector.h. 00181 {
00182 return new GaussFFTProjector();
00183 }
|
|
|
Project an 3D image into a 2D image.
Implements EMAN::Projector. Definition at line 71 of file projector.cpp. References EMAN::EMData::ap2ri(), data, EMAN::EMData::do_fft(), EMAN::EMData::do_ift(), EMAN::EMData::get_data(), EMAN::Transform::get_mirror(), EMAN::EMData::get_ndim(), EMAN::Transform::get_rotation_transform(), EMAN::Transform::get_scale(), EMAN::Transform::get_trans(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), ImageDimensionException, interp_ft_3d(), EMAN::Transform::invert(), EMAN::EMData::is_complex(), LOGERR, NullPointerException, EMAN::EMData::process_inplace(), EMAN::EMData::set_attr(), EMAN::EMData::set_complex(), EMAN::EMData::set_ri(), EMAN::EMData::set_size(), EMAN::EMData::translate(), EMAN::EMData::update(), EMAN::Vec3f, x, and y. 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 }
|
|
|
Set the projector parameters using a key/value dictionary.
Reimplemented from EMAN::Projector. Definition at line 162 of file projector.h. References phi. 00163 {
00164 Projector::set_params(new_params);
00165 alt = params["alt"];
00166 az = params["az"];
00167 phi = params["phi"];
00168 }
|
|
|
Definition at line 196 of file projector.h. |
|
|
Definition at line 196 of file projector.h. |
|
|
Definition at line 54 of file projector.cpp. |
|
|
Definition at line 196 of file projector.h. |
1.3.9.1