#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 | |
static Projector * | NEW () |
Static Public Attributes | |
static 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.
EMAN::GaussFFTProjector::GaussFFTProjector | ( | ) | [inline] |
Back-project a 2D image into a 3D image.
Implements EMAN::Projector.
Definition at line 2004 of file projector.cpp.
string EMAN::GaussFFTProjector::get_desc | ( | ) | const [inline, virtual] |
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 }
string EMAN::GaussFFTProjector::get_name | ( | ) | const [inline, virtual] |
Get the projector's name.
Each projector is indentified by unique name.
Implements EMAN::Projector.
Definition at line 170 of file projector.h.
References NAME.
00171 { 00172 return NAME; 00173 }
TypeDict EMAN::GaussFFTProjector::get_param_types | ( | ) | const [inline, virtual] |
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::EMObject::INT, EMAN::TypeDict::put(), and EMAN::EMObject::TRANSFORM.
00186 { 00187 TypeDict d; 00188 d.put("transform", EMObject::TRANSFORM); 00189 d.put("mode", EMObject::INT); 00190 return d; 00191 }
bool GaussFFTProjector::interp_ft_3d | ( | int | mode, | |
EMData * | image, | |||
float | x, | |||
float | y, | |||
float | z, | |||
float * | data, | |||
float | gauss_width | |||
) | const [private] |
Definition at line 208 of file projector.cpp.
References abs, EMAN::Util::agauss(), 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(), and sqrt().
Referenced by project3d().
00210 { 00211 float *rdata = image->get_data(); 00212 int nx = image->get_xsize(); 00213 int ny = image->get_ysize(); 00214 int nz = image->get_zsize(); 00215 00216 if ( mode == 0 ) mode = 2; 00217 00218 if (mode == 1) { 00219 int x0 = 2 * (int) floor(x + 0.5f); 00220 int y0 = (int) floor(y + 0.5f); 00221 int z0 = (int) floor(z + 0.5f); 00222 00223 data[0] = rdata[x0 + y0 * nx + z0 * nx * ny]; 00224 data[1] = rdata[x0 + y0 * nx + z0 * nx * ny + 1]; 00225 return true; 00226 } 00227 else if (mode == 2) { 00228 int x0 = (int) floor(x); 00229 int y0 = (int) floor(y); 00230 int z0 = (int) floor(z); 00231 00232 float dx = x - x0; 00233 float dy = y - y0; 00234 float dz = z - z0; 00235 00236 if (x0 <= nx/2- 2 && y0 <= ny - 1 && z0 <= nz - 1) { 00237 int i = (int) (x0 * 2 + y0 * nx + z0 * nx * ny); 00238 00239 float n = Util::agauss(1, dx, dy, dz, gw) + 00240 Util::agauss(1, 1 - dx, dy, dz, gw) + 00241 Util::agauss(1, dx, 1 - dy, dz, gw) + 00242 Util::agauss(1, 1 - dx, 1 - dy, dz, gw) + 00243 Util::agauss(1, dx, dy, 1 - dz, gw) + 00244 Util::agauss(1, 1 - dx, dy, 1 - dz, gw) + 00245 Util::agauss(1, dx, 1 - dy, 1 - dz, gw) + Util::agauss(1, 1 - dx, 1 - dy, 1 - dz, gw); 00246 00247 data[0] = Util::agauss(rdata[i], dx, dy, dz, gw) + 00248 Util::agauss(rdata[i + 2], 1 - dx, dy, dz, gw) + 00249 Util::agauss(rdata[i + nx], dx, 1 - dy, dz, gw) + 00250 Util::agauss(rdata[i + nx + 2], 1 - dx, 1 - dy, dz, gw) + 00251 Util::agauss(rdata[i + nx * ny], dx, dy, 1 - dz, gw) + 00252 Util::agauss(rdata[i + 2 + nx * ny], 1 - dx, dy, 1 - dz, gw) + 00253 Util::agauss(rdata[i + nx + nx * ny], dx, 1 - dy, 1 - dz, gw) + 00254 Util::agauss(rdata[i + 2 + nx + nx * ny], 1 - dx, 1 - dy, 1 - dz, gw) / n; 00255 00256 i++; 00257 00258 data[1] = Util::agauss(rdata[i], dx, dy, dz, gw) + 00259 Util::agauss(rdata[i + 2], 1 - dx, dy, dz, gw) + 00260 Util::agauss(rdata[i + nx], dx, 1 - dy, dz, gw) + 00261 Util::agauss(rdata[i + nx + 2], 1 - dx, 1 - dy, dz, gw) + 00262 Util::agauss(rdata[i + nx * ny], dx, dy, 1 - dz, gw) + 00263 Util::agauss(rdata[i + 2 + nx * ny], 1 - dx, dy, 1 - dz, gw) + 00264 Util::agauss(rdata[i + nx + nx * ny], dx, 1 - dy, 1 - dz, gw) + 00265 Util::agauss(rdata[i + 2 + nx + nx * ny], 1 - dx, 1 - dy, 1 - dz, gw) / n; 00266 return true; 00267 } 00268 return false; 00269 } 00270 else if (mode == 3) { 00271 int x0 = 2 * (int) floor(x + .5); 00272 int y0 = (int) floor(y + .5); 00273 int z0 = (int) floor(z + .5); 00274 00275 float *supp = image->setup4slice(); 00276 00277 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) { 00278 float n = 0; 00279 00280 if (x0 == 0) { 00281 x0 += 4; 00282 size_t idx; 00283 float r, g; 00284 for (int k = z0 - 1; k <= z0 + 1; k++) { 00285 for (int j = y0 - 1; j <= y0 + 1; j++) { 00286 for (int i = x0 - 2; i <= x0 + 2; i += 2) { 00287 r = Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z); 00288 g = exp(-r / gw); 00289 n += g; 00290 idx = i + j * 12 + (size_t)k * 12 * ny; 00291 data[0] += supp[idx] * g; 00292 data[1] += supp[idx + 1] * g; 00293 } 00294 } 00295 } 00296 } 00297 else { 00298 size_t idx; 00299 float r, g; 00300 for (int k = z0 - 1; k <= z0 + 1; k++) { 00301 for (int j = y0 - 1; j <= y0 + 1; j++) { 00302 for (int i = x0 - 2; i <= x0 + 2; i += 2) { 00303 r = Util::hypot3sq(i / 2.0f - x, j - y, k - z); 00304 g = exp(-r / gw); 00305 n += g; 00306 idx = i + j * nx + (size_t)k * nx * ny; 00307 data[0] += rdata[idx] * g; 00308 data[1] += rdata[idx + 1] * g; 00309 } 00310 } 00311 } 00312 } 00313 00314 data[0] /= n; 00315 data[1] /= n; 00316 return true; 00317 } 00318 return false; 00319 } 00320 else if (mode == 4) { 00321 int x0 = 2 * (int) floor(x); 00322 int y0 = (int) floor(y); 00323 int z0 = (int) floor(z); 00324 00325 float *supp = image->setup4slice(); 00326 00327 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) { 00328 float n = 0; 00329 00330 if (x0 == 0) { 00331 x0 += 4; 00332 size_t idx; 00333 float r, g; 00334 for (int k = z0 - 1; k <= z0 + 2; k++) { 00335 for (int j = y0 - 1; j <= y0 + 2; j++) { 00336 for (int i = x0 - 2; i <= x0 + 4; i += 2) { 00337 r = Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z); 00338 g = exp(-r / gw); 00339 n += g; 00340 idx = i + j * 12 + (size_t)k * 12 * ny; 00341 data[0] += supp[idx] * g; 00342 data[1] += supp[idx + 1] * g; 00343 } 00344 } 00345 } 00346 } 00347 else { 00348 float r, g; 00349 size_t idx; 00350 for (int k = z0 - 1; k <= z0 + 2; k++) { 00351 for (int j = y0 - 1; j <= y0 + 2; j++) { 00352 for (int i = x0 - 2; i <= x0 + 4; i += 2) { 00353 r = Util::hypot3sq(i / 2.0f - x, j - y, k - z); 00354 g = exp(-r / gw); 00355 n += g; 00356 idx = i + j * nx + (size_t)k * nx * ny; 00357 data[0] += rdata[idx] * g; 00358 data[1] += rdata[idx + 1] * g; 00359 } 00360 } 00361 } 00362 } 00363 data[0] /= n; 00364 data[1] /= n; 00365 return true; 00366 } 00367 return false; 00368 } 00369 else if (mode == 5) { 00370 int x0 = 2 * (int) floor(x + .5); 00371 int y0 = (int) floor(y + .5); 00372 int z0 = (int) floor(z + .5); 00373 00374 float *supp = image->setup4slice(); 00375 float *gimx = Interp::get_gimx(); 00376 00377 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) { 00378 int mx0 = -(int) floor((x - x0 / 2) * 39.0 + .5) - 78; 00379 int my0 = -(int) floor((y - y0) * 39.0 + .5) - 78; 00380 int mz0 = -(int) floor((z - z0) * 39.0 + .5) - 78; 00381 00382 float n = 0; 00383 int mmz = mz0; 00384 int mmy = my0; 00385 int mmx = mx0; 00386 00387 if (x0 < 4) { 00388 x0 += 4; 00389 00390 size_t idx; 00391 float g; 00392 for (int k = z0 - 2; k <= z0 + 2; k++, mmz += 39) { 00393 for (int j = y0 - 2; j <= y0 + 2; j++, mmy += 39) { 00394 for (int i = x0 - 4; i <= x0 + 4; i += 2, mmx += 39) { 00395 g = gimx[abs(mmx) + abs(mmy) * 100 + abs(mmz) * 10000]; 00396 n += g; 00397 idx = i + j * 12 + (size_t)k * 12 * ny; 00398 data[0] += supp[idx] * g; 00399 data[1] += supp[idx + 1] * g; 00400 } 00401 } 00402 } 00403 } 00404 else { 00405 size_t ii; 00406 float g; 00407 for (int k = z0 - 2; k <= z0 + 2; k++, mmz += 39) { 00408 for (int j = y0 - 2; j <= y0 + 2; j++, mmy += 39) { 00409 for (int i = x0 - 4; i <= x0 + 4; i += 2, mmx += 39) { 00410 ii = i + j * nx + (size_t)k * nx * ny; 00411 g = gimx[abs(mmx) + abs(mmy) * 100 + abs(mmz) * 10000]; 00412 n += g; 00413 data[0] += rdata[ii] * g; 00414 data[1] += rdata[ii + 1] * g; 00415 } 00416 } 00417 } 00418 } 00419 00420 data[0] /= n; 00421 data[1] /= n; 00422 return true; 00423 } 00424 return false; 00425 } 00426 else if (mode == 6) { 00427 00428 int x0 = 2 * (int) floor(x + .5); 00429 int y0 = (int) floor(y + .5); 00430 int z0 = (int) floor(z + .5); 00431 00432 float *supp = image->setup4slice(); 00433 00434 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) { 00435 float n = 0; 00436 00437 if (x0 < 4) { 00438 x0 += 4; 00439 float r, g; 00440 size_t idx; 00441 for (int k = z0 - 2; k <= z0 + 2; k++) { 00442 for (int j = y0 - 2; j <= y0 + 2; j++) { 00443 for (int i = x0 - 4; i <= x0 + 4; i += 2) { 00444 r = Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z); 00445 g = exp(-r / gw); 00446 n += g; 00447 idx = i + j * 12 + (size_t)k * 12 * ny; 00448 data[0] += supp[idx] * g; 00449 data[1] += supp[idx + 1] * g; 00450 } 00451 } 00452 } 00453 } 00454 else { 00455 float r, g; 00456 size_t idx; 00457 for (int k = z0 - 2; k <= z0 + 2; k++) { 00458 for (int j = y0 - 2; j <= y0 + 2; j++) { 00459 for (int i = x0 - 4; i <= x0 + 4; i += 2) { 00460 r = Util::hypot3sq(i / 2.0f - x, j - y, k - z); 00461 g = exp(-r / gw); 00462 n += g; 00463 idx = i + j * nx + (size_t)k * nx * ny; 00464 data[0] += rdata[idx] * g; 00465 data[1] += rdata[idx + 1] * g; 00466 } 00467 00468 } 00469 } 00470 } 00471 00472 data[0] /= n; 00473 data[1] /= n; 00474 00475 return true; 00476 } 00477 return false; 00478 } 00479 else if (mode == 7) { 00480 int x0 = 2 * (int) floor(x + .5); 00481 int y0 = (int) floor(y + .5); 00482 int z0 = (int) floor(z + .5); 00483 00484 float *supp = image->setup4slice(); 00485 00486 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) { 00487 float n = 0; 00488 if (x0 < 4) { 00489 x0 += 4; 00490 float r, g; 00491 size_t idx; 00492 for (int k = z0 - 2; k <= z0 + 2; k++) { 00493 for (int j = y0 - 2; j <= y0 + 2; j++) { 00494 for (int i = x0 - 4; i <= x0 + 4; i += 2) { 00495 r = sqrt(Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z)); 00496 g = Interp::hyperg(r); 00497 n += g; 00498 idx = i + j * 12 + (size_t)k * 12 * ny; 00499 data[0] += supp[idx] * g; 00500 data[1] += supp[idx + 1] * g; 00501 } 00502 } 00503 } 00504 } 00505 else { 00506 float r, g; 00507 size_t idx; 00508 for (int k = z0 - 2; k <= z0 + 2; k++) { 00509 for (int j = y0 - 2; j <= y0 + 2; j++) { 00510 for (int i = x0 - 4; i <= x0 + 4; i += 2) { 00511 r = sqrt(Util::hypot3sq(i / 2.0f - x, j - y, k - z)); 00512 g = Interp::hyperg(r); 00513 n += g; 00514 idx = i + j * nx + (size_t)k * nx * ny; 00515 data[0] += rdata[idx] * g; 00516 data[1] += rdata[idx + 1] * g; 00517 } 00518 00519 } 00520 } 00521 } 00522 data[0] /= n; 00523 data[1] /= n; 00524 return true; 00525 } 00526 return false; 00527 00528 } 00529 // throw InvalidParameterException("Error, unkown mode"); 00530 return false; 00531 }
static Projector* EMAN::GaussFFTProjector::NEW | ( | ) | [inline, static] |
Definition at line 180 of file projector.h.
References GaussFFTProjector().
00181 { 00182 return new GaussFFTProjector(); 00183 }
Project an 3D image into a 2D image.
Implements EMAN::Projector.
Definition at line 73 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(), EMAN::EMConsts::I2G, EMAN::EMConsts::I3G, EMAN::EMConsts::I4G, EMAN::EMConsts::I5G, ImageDimensionException, interp_ft_3d(), EMAN::Transform::invert(), EMAN::EMData::is_complex(), LOGERR, NullPointerException, EMAN::Projector::params, 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(), x, and y.
00074 { 00075 Transform* t3d = params["transform"]; 00076 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified"); 00077 if ( image->get_ndim() != 3 ) throw ImageDimensionException("Error, the projection volume must be 3D"); 00078 00079 EMData *f = image; 00080 if (!image->is_complex()) { 00081 image->process_inplace("xform.phaseorigin.tocorner"); 00082 f = image->do_fft(); 00083 f->process_inplace("xform.fourierorigin.tocenter"); 00084 image->process_inplace("xform.phaseorigin.tocenter"); 00085 } 00086 00087 int f_nx = f->get_xsize(); 00088 int f_ny = f->get_ysize(); 00089 int f_nz = f->get_zsize(); 00090 00091 if (!f->is_complex() || f_nz != f_ny || f_nx != f_ny + 2) { 00092 LOGERR("Cannot project this image"); 00093 return 0; 00094 } 00095 00096 f->ap2ri(); 00097 00098 EMData *tmp = new EMData(); 00099 tmp->set_size(f_nx, f_ny, 1); 00100 tmp->set_complex(true); 00101 tmp->set_ri(true); 00102 00103 float *data = tmp->get_data(); 00104 00105 Transform r = t3d->get_rotation_transform(); 00106 r.invert(); 00107 float scale = t3d->get_scale(); 00108 00109 int mode = params["mode"]; 00110 float gauss_width = 1; 00111 if ( mode == 0 ) mode = 2; 00112 if (mode == 2 ) { 00113 gauss_width = EMConsts::I2G; 00114 } 00115 else if (mode == 3) { 00116 gauss_width = EMConsts::I3G; 00117 } 00118 else if (mode == 4) { 00119 gauss_width = EMConsts::I4G; 00120 } 00121 else if (mode == 6 || mode == 7) { 00122 gauss_width = EMConsts::I5G; 00123 } 00124 00125 for (int y = 0; y < f_ny; y++) { 00126 for (int x = 0; x < f_nx / 2; x++) { 00127 int ii = x * 2 + y * f_nx; 00128 #ifdef _WIN32 00129 if (_hypot(x, y - f_ny / 2) >= f_ny / 2 - 1) { 00130 #else 00131 if (hypot(x, y - f_ny / 2) >= f_ny / 2 - 1) { 00132 #endif //_WIN32 00133 data[ii] = 0; 00134 data[ii + 1] = 0; 00135 } 00136 else { 00137 Vec3f coord(x,(y - f_ny / 2),0); 00138 Vec3f soln = r*coord; 00139 float xx = soln[0]; 00140 float yy = soln[1]; 00141 float zz = soln[2]; 00142 00143 int cc = 1; 00144 if (xx < 0) { 00145 xx = -xx; 00146 yy = -yy; 00147 zz = -zz; 00148 cc = -1; 00149 } 00150 00151 if (scale != 1.0) { 00152 xx *= scale; 00153 yy *= scale; 00154 zz *= scale; 00155 } 00156 00157 yy += f_ny / 2; 00158 zz += f_nz / 2; 00159 00160 if (yy < 0 || xx < 0 || zz < 0) { 00161 data[ii] = 0; 00162 data[ii+1] = 0; 00163 continue; 00164 } 00165 00166 if (interp_ft_3d(mode, f, xx, yy, zz, data + ii, gauss_width)) { 00167 data[ii + 1] *= cc; 00168 } else { 00169 data[ii] = 0; 00170 data[ii+1] = 0; 00171 } 00172 } 00173 } 00174 } 00175 00176 f->update(); 00177 tmp->update(); 00178 00179 tmp->process_inplace("xform.fourierorigin.tocorner"); 00180 EMData *ret = tmp->do_ift(); 00181 ret->process_inplace("xform.phaseorigin.tocenter"); 00182 00183 ret->translate(t3d->get_trans()); 00184 00185 if (t3d->get_mirror() ) ret->process_inplace("xform.flip",Dict("axis","x")); 00186 00187 Dict filter_d; 00188 filter_d["gauss_width"] = gauss_width; 00189 filter_d["ring_width"] = ret->get_xsize() / 2; 00190 ret->process_inplace("math.gausskernelfix", filter_d); 00191 00192 if( tmp ) 00193 { 00194 delete tmp; 00195 tmp = 0; 00196 } 00197 00198 ret->set_attr("xform.projection",t3d); 00199 ret->update(); 00200 00201 if(t3d) {delete t3d; t3d=0;} 00202 00203 return ret; 00204 }
void EMAN::GaussFFTProjector::set_params | ( | const Dict & | new_params | ) | [inline] |
Set the projector parameters using a key/value dictionary.
Reimplemented from EMAN::Projector.
Definition at line 162 of file projector.h.
References alt, az, EMAN::Projector::params, phi, and EMAN::Projector::set_params().
00163 { 00164 Projector::set_params(new_params); 00165 alt = params["alt"]; 00166 az = params["az"]; 00167 phi = params["phi"]; 00168 }
float EMAN::GaussFFTProjector::alt [private] |
float EMAN::GaussFFTProjector::az [private] |
const string GaussFFTProjector::NAME = "gauss_fft" [static] |
float EMAN::GaussFFTProjector::phi [private] |