#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 1849 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 206 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().
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 }
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 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(), 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.
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 }
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] |