#include <reconstructor.h>
Inheritance diagram for EMAN::WienerFourierReconstructor:
Public Member Functions | |
WienerFourierReconstructor () | |
virtual | ~WienerFourierReconstructor () |
virtual void | setup () |
Initialize the reconstructor. | |
virtual int | insert_slice (const EMData *const slice, const Transform &euler, const float weight=1.0) |
Insert an image slice to the reconstructor. | |
virtual EMData * | finish (bool doift=true) |
Finish reconstruction and return the complete model. | |
virtual string | get_name () const |
Get the unique name of this class (especially for factory based instantiation access). | |
virtual string | get_desc () const |
Get a clear, concise description of this class. | |
virtual TypeDict | get_param_types () const |
Static Public Member Functions | |
static Reconstructor * | NEW () |
Static Public Attributes | |
static const string | NAME = "wiener_fourier" |
Private Member Functions | |
WienerFourierReconstructor (const WienerFourierReconstructor &that) | |
WienerFourierReconstructor & | operator= (const WienerFourierReconstructor &) |
void | load_default_settings () |
Definition at line 645 of file reconstructor.h.
|
Definition at line 648 of file reconstructor.h. References load_default_settings(). Referenced by NEW(). 00648 { load_default_settings(); };
|
|
Definition at line 649 of file reconstructor.h.
|
|
|
|
Finish reconstruction and return the complete model.
Reimplemented from EMAN::Reconstructor. Definition at line 1432 of file reconstructor.cpp. References EMAN::EMData::get_data(), EMAN::ReconstructorVolumeData::image, norm(), EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::nz, rdata, EMAN::ReconstructorVolumeData::tmp_data, and EMAN::EMData::update(). 01433 { 01434 float *norm = tmp_data->get_data(); 01435 float *rdata = image->get_data(); 01436 01437 size_t size = nx * ny * nz; 01438 for (size_t i = 0; i < size; i += 2) { 01439 float d = norm[i]; 01440 if (d == 0) { 01441 rdata[i] = 0; 01442 rdata[i + 1] = 0; 01443 } 01444 else { 01445 float w = 1 + 1 / d; 01446 rdata[i] /= d * w; 01447 rdata[i + 1] /= d * w; 01448 } 01449 } 01450 01451 if( tmp_data ) { 01452 delete tmp_data; 01453 tmp_data = 0; 01454 } 01455 image->update(); 01456 return image; 01457 }
|
|
Get a clear, concise description of this class.
Implements EMAN::FactoryBase. Definition at line 670 of file reconstructor.h. 00671 { 00672 return "Experimental - Direct Fourier reconstruction taking into account the Wiener filtration of the individual images."; 00673 }
|
|
Get the unique name of this class (especially for factory based instantiation access).
Implements EMAN::FactoryBase. Definition at line 665 of file reconstructor.h. References NAME. 00666 { 00667 return NAME; 00668 }
|
|
Implements EMAN::FactoryBase. Definition at line 680 of file reconstructor.h. References EMAN::EMObject::BOOL, EMAN::EMObject::FLOAT, EMAN::EMObject::FLOATARRAY, EMAN::EMObject::INT, and EMAN::TypeDict::put(). 00681 { 00682 TypeDict d; 00683 // FIXME:: double check all of thes are need, expecially dlog and weight 00684 d.put("size", EMObject::INT); 00685 d.put("mode", EMObject::INT); 00686 d.put("weight", EMObject::FLOAT); 00687 d.put("use_weights", EMObject::BOOL); 00688 d.put("dlog", EMObject::BOOL); 00689 d.put("padratio", EMObject::FLOAT); 00690 d.put("snr", EMObject::FLOATARRAY); 00691 return d; 00692 }
|
|
Insert an image slice to the reconstructor. To insert multiple image slices, call this function multiple times.
Reimplemented from EMAN::Reconstructor. Definition at line 1460 of file reconstructor.cpp. References abs, EMAN::Util::agauss(), EMAN::Ctf::CTFOS, dt, EMAN::EMData::get_data(), EMAN::Interp::get_gimx(), EMAN::Interp::hyperg(), EMAN::Util::hypot3sq(), EMAN::EMConsts::I2G, EMAN::EMConsts::I3G, EMAN::EMConsts::I4G, EMAN::EMConsts::I5G, EMAN::ReconstructorVolumeData::image, EMAN::EMData::is_complex(), LOGERR, norm(), EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::nz, EMAN::FactoryBase::params, rdata, EMAN::Util::round(), sqrt(), EMAN::Util::square(), square, EMAN::ReconstructorVolumeData::tmp_data, EMAN::EMData::update(), x, and y. 01461 { 01462 if (!slice) { 01463 LOGERR("try to insert NULL slice"); 01464 return 1; 01465 } 01466 01467 int mode = params["mode"]; 01468 float padratio = params["padratio"]; 01469 vector < float >snr = params["snr"]; 01470 01471 if (!slice->is_complex()) { 01472 LOGERR("Only complex slice can be inserted."); 01473 return 1; 01474 } 01475 float *gimx = 0; 01476 if (mode == 5) { 01477 gimx = Interp::get_gimx(); 01478 } 01479 01480 int nxy = nx * ny; 01481 int off[8] = {0,0,0,0,0,0,0,0}; 01482 if (mode == 2) { 01483 off[0] = 0; 01484 off[1] = 2; 01485 off[2] = nx; 01486 off[3] = nx + 2; 01487 off[4] = nxy; 01488 off[5] = nxy + 2; 01489 off[6] = nxy + nx; 01490 off[7] = nxy + nx + 2; 01491 } 01492 01493 float *norm = tmp_data->get_data(); 01494 float *dat = slice->get_data(); 01495 float *rdata = image->get_data(); 01496 01497 int rl = Util::square(ny / 2 - 1); 01498 float dt[2]; 01499 float g[8]; 01500 01501 for (int y = 0; y < ny; y++) { 01502 for (int x = 0; x < nx / 2; x++) { 01503 if ((x * x + Util::square(y - ny / 2)) >= rl) 01504 { 01505 continue; 01506 } 01507 01508 #ifdef _WIN32 01509 int r = Util::round((float)_hypot(x, (float) y - ny / 2) * Ctf::CTFOS / padratio); 01510 #else 01511 int r = Util::round((float)hypot(x, (float) y - ny / 2) * Ctf::CTFOS / padratio); 01512 #endif 01513 if (r >= Ctf::CTFOS * ny / 2) { 01514 r = Ctf::CTFOS * ny / 2 - 1; 01515 } 01516 01517 float weight = snr[r]; 01518 01519 float xx = (x * euler[0][0] + (y - ny / 2) * euler[0][1]); 01520 float yy = (x * euler[1][0] + (y - ny / 2) * euler[1][1]); 01521 float zz = (x * euler[2][0] + (y - ny / 2) * euler[2][1]); 01522 float cc = 1; 01523 01524 if (xx < 0) { 01525 xx = -xx; 01526 yy = -yy; 01527 zz = -zz; 01528 cc = -1.0; 01529 } 01530 01531 yy += ny / 2; 01532 zz += nz / 2; 01533 01534 dt[0] = dat[x * 2 + y * nx] * (1 + 1.0f / weight); 01535 dt[1] = cc * dat[x * 2 + 1 + y * nx] * (1 + 1.0f / weight); 01536 01537 int x0 = 0; 01538 int y0 = 0; 01539 int z0 = 0; 01540 int i = 0; 01541 int l = 0; 01542 float dx = 0; 01543 float dy = 0; 01544 float dz = 0; 01545 01546 int mx0 = 0; 01547 int my0 = 0; 01548 int mz0 = 0; 01549 01550 size_t idx; 01551 switch (mode) { 01552 case 1: 01553 x0 = 2 * (int) floor(xx + 0.5f); 01554 y0 = (int) floor(yy + 0.5f); 01555 z0 = (int) floor(zz + 0.5f); 01556 01557 rdata[x0 + y0 * nx + z0 * nxy] += weight * dt[0]; 01558 rdata[x0 + y0 * nx + z0 * nxy + 1] += weight * dt[1]; 01559 norm[x0 + y0 * nx + z0 * nxy] += weight; 01560 break; 01561 01562 case 2: 01563 x0 = (int) floor(xx); 01564 y0 = (int) floor(yy); 01565 z0 = (int) floor(zz); 01566 01567 dx = xx - x0; 01568 dy = yy - y0; 01569 dz = zz - z0; 01570 01571 weight /= (float)pow((float)(EMConsts::I2G * M_PI), 1.5f); 01572 01573 if (x0 > nx - 2 || y0 > ny - 1 || z0 > nz - 1) { 01574 break; 01575 } 01576 01577 i = (int) (x0 * 2 + y0 * nx + z0 * nxy); 01578 01579 01580 g[0] = Util::agauss(1, dx, dy, dz, EMConsts::I2G); 01581 g[1] = Util::agauss(1, 1 - dx, dy, dz, EMConsts::I2G); 01582 g[2] = Util::agauss(1, dx, 1 - dy, dz, EMConsts::I2G); 01583 g[3] = Util::agauss(1, 1 - dx, 1 - dy, dz, EMConsts::I2G); 01584 g[4] = Util::agauss(1, dx, dy, 1 - dz, EMConsts::I2G); 01585 g[5] = Util::agauss(1, 1 - dx, dy, 1 - dz, EMConsts::I2G); 01586 g[6] = Util::agauss(1, dx, 1 - dy, 1 - dz, EMConsts::I2G); 01587 g[7] = Util::agauss(1, 1 - dx, 1 - dy, 1 - dz, EMConsts::I2G); 01588 01589 for (int j = 0; j < 8; j++) { 01590 int k = i + off[j]; 01591 rdata[k] += weight * dt[0] * g[j]; 01592 rdata[k + 1] += weight * dt[1] * g[j]; 01593 norm[k] += weight * g[j]; 01594 } 01595 01596 break; 01597 case 3: 01598 x0 = 2 * (int) floor(xx + 0.5f); 01599 y0 = (int) floor(yy + 0.5f); 01600 z0 = (int) floor(zz + 0.5f); 01601 01602 weight /= (float)pow((float)(EMConsts::I3G * M_PI), 1.5f); 01603 01604 if (x0 >= nx - 4 || y0 > ny - 3 || z0 > nz - 3 || y0 < 2 || z0 < 2) { 01605 break; 01606 } 01607 01608 l = x0 - 2; 01609 if (x0 == 0) { 01610 l = x0; 01611 } 01612 01613 for (int k = z0 - 1; k <= z0 + 1; k++) { 01614 for (int j = y0 - 1; j <= y0 + 1; j++) { 01615 for (int i = l; i <= x0 + 2; i += 2) { 01616 float r = Util::hypot3sq((float) i / 2 - xx, j - yy, k - zz); 01617 float gg = exp(-r / EMConsts::I3G); 01618 01619 idx = i + j * nx + k * nxy; 01620 rdata[idx] += weight * gg * dt[0]; 01621 rdata[idx + 1] += weight * gg * dt[1]; 01622 norm[idx] += weight * gg; 01623 } 01624 } 01625 } 01626 break; 01627 01628 case 4: 01629 x0 = 2 * (int) floor(xx); 01630 y0 = (int) floor(yy); 01631 z0 = (int) floor(zz); 01632 01633 weight /= (float)pow((float)(EMConsts::I4G * M_PI), 1.5f); 01634 01635 if (x0 >= nx - 4 || y0 > ny - 3 || z0 > nz - 3 || y0 < 2 || z0 < 2) { 01636 break; 01637 } 01638 01639 l = x0 - 2; 01640 if (x0 == 0) { 01641 l = x0; 01642 } 01643 01644 for (int k = z0 - 1; k <= z0 + 2; ++k) { 01645 for (int j = y0 - 1; j <= y0 + 2; ++j) { 01646 for (int i = l; i <= x0 + 4; i += 2) { 01647 float r = Util::hypot3sq((float) i / 2 - xx, j - yy, k - zz); 01648 float gg = exp(-r / EMConsts::I4G); 01649 01650 idx = i + j * nx + k * nxy; 01651 rdata[idx] += weight * gg * dt[0]; 01652 rdata[idx + 1] += weight * gg * dt[1]; 01653 norm[idx] += weight * gg; 01654 } 01655 } 01656 } 01657 break; 01658 01659 case 5: 01660 x0 = (int) floor(xx + .5); 01661 y0 = (int) floor(yy + .5); 01662 z0 = (int) floor(zz + .5); 01663 01664 weight /= (float)pow((float)(EMConsts::I5G * M_PI), 1.5f); 01665 01666 mx0 = -(int) floor((xx - x0) * 39.0f + 0.5) - 78; 01667 my0 = -(int) floor((yy - y0) * 39.0f + 0.5) - 78; 01668 mz0 = -(int) floor((zz - z0) * 39.0f + 0.5) - 78; 01669 x0 *= 2; 01670 01671 if (x0 >= nx - 4 || y0 > ny - 3 || z0 > nz - 3 || y0 < 2 || z0 < 2) 01672 break; 01673 01674 if (x0 == 0) { 01675 l = 0; 01676 mx0 += 78; 01677 } 01678 else if (x0 == 2) { 01679 l = 0; 01680 mx0 += 39; 01681 } 01682 else 01683 l = x0 - 4; 01684 for (int k = z0 - 2, mmz = mz0; k <= z0 + 2; k++, mmz += 39) { 01685 for (int j = y0 - 2, mmy = my0; j <= y0 + 2; j++, mmy += 39) { 01686 for (int i = l, mmx = mx0; i <= x0 + 4; i += 2, mmx += 39) { 01687 size_t ii = i + j * nx + k * nxy; 01688 float gg = weight * gimx[abs(mmx) + abs(mmy) * 100 + abs(mmz) * 10000]; 01689 01690 rdata[ii] += gg * dt[0]; 01691 rdata[ii + 1] += gg * dt[1]; 01692 norm[ii] += gg; 01693 } 01694 } 01695 } 01696 01697 if (x0 <= 2) { 01698 xx = -xx; 01699 yy = -(yy - ny / 2) + ny / 2; 01700 zz = -(zz - nz / 2) + nz / 2; 01701 x0 = (int) floor(xx + 0.5f); 01702 y0 = (int) floor(yy + 0.5f); 01703 z0 = (int) floor(zz + 0.5f); 01704 int mx0 = -(int) floor((xx - x0) * 39.0f + .5); 01705 x0 *= 2; 01706 01707 if (y0 > ny - 3 || z0 > nz - 3 || y0 < 2 || z0 < 2) 01708 break; 01709 01710 for (int k = z0 - 2, mmz = mz0; k <= z0 + 2; k++, mmz += 39) { 01711 for (int j = y0 - 2, mmy = my0; j <= y0 + 2; j++, mmy += 39) { 01712 for (int i = 0, mmx = mx0; i <= x0 + 4; i += 2, mmx += 39) { 01713 size_t ii = i + j * nx + k * nxy; 01714 float gg = 01715 weight * gimx[abs(mmx) + abs(mmy) * 100 + abs(mmz) * 10000]; 01716 01717 rdata[ii] += gg * dt[0]; 01718 rdata[ii + 1] -= gg * dt[1]; 01719 norm[ii] += gg; 01720 } 01721 } 01722 } 01723 } 01724 break; 01725 // mode 6 is now mode 5 without the fast interpolation 01726 case 6: 01727 x0 = 2 * (int) floor(xx + .5); 01728 y0 = (int) floor(yy + .5); 01729 z0 = (int) floor(zz + .5); 01730 01731 weight /= (float)pow((float)(EMConsts::I5G * M_PI), 1.5f); 01732 01733 if (x0 >= nx - 4 || y0 > ny - 3 || z0 > nz - 3 || y0 < 2 || z0 < 2) 01734 break; 01735 01736 if (x0 <= 2) 01737 l = 0; 01738 else 01739 l = x0 - 4; 01740 for (int k = z0 - 2; k <= z0 + 2; ++k) { 01741 for (int j = y0 - 2; j <= y0 + 2; ++j) { 01742 for (int i = l; i <= x0 + 4; i += 2) { 01743 size_t ii = i + j * nx + k * nxy; 01744 float r = Util::hypot3sq((float) i / 2 - xx, (float) j - yy, 01745 (float) k - zz); 01746 float gg = weight * exp(-r / EMConsts::I5G); 01747 01748 rdata[ii] += gg * dt[0]; 01749 rdata[ii + 1] += gg * dt[1]; 01750 norm[ii] += gg; 01751 } 01752 } 01753 } 01754 01755 if (x0 <= 2) { 01756 xx = -xx; 01757 yy = -(yy - ny / 2) + ny / 2; 01758 zz = -(zz - nz / 2) + nz / 2; 01759 x0 = 2 * (int) floor(xx + 0.5f); 01760 y0 = (int) floor(yy + 0.5f); 01761 z0 = (int) floor(zz + 0.5f); 01762 01763 if (y0 > ny - 3 || z0 > nz - 3 || y0 < 2 || z0 < 2) 01764 break; 01765 01766 for (int k = z0 - 2; k <= z0 + 2; ++k) { 01767 for (int j = y0 - 2; j <= y0 + 2; ++j) { 01768 for (int i = 0; i <= x0 + 4; i += 2) { 01769 size_t ii = i + j * nx + k * nxy; 01770 float r = Util::hypot3sq((float) i / 2 - xx, (float) j - yy, 01771 (float) k - zz); 01772 float gg = weight * exp(-r / EMConsts::I5G); 01773 01774 rdata[ii] += gg * dt[0]; 01775 rdata[ii + 1] -= gg * dt[1]; // note the -, complex conj. 01776 norm[ii] += gg; 01777 } 01778 } 01779 } 01780 } 01781 break; 01782 01783 case 7: 01784 x0 = 2 * (int) floor(xx + .5); 01785 y0 = (int) floor(yy + .5); 01786 z0 = (int) floor(zz + .5); 01787 01788 if (x0 >= nx - 4 || y0 > ny - 3 || z0 > nz - 3 || y0 < 2 || z0 < 2) 01789 break; 01790 01791 if (x0 <= 2) 01792 l = 0; 01793 else 01794 l = x0 - 4; 01795 for (int k = z0 - 2; k <= z0 + 2; k++) { 01796 for (int j = y0 - 2; j <= y0 + 2; j++) { 01797 for (int i = l; i <= x0 + 4; i += 2) { 01798 size_t ii = i + j * nx + k * nxy; 01799 float r = (float)sqrt(Util::hypot3sq((float) i / 2 - xx, 01800 (float) j - yy, 01801 (float) k - zz)); 01802 float gg = weight * Interp::hyperg(r); 01803 01804 rdata[ii] += gg * dt[0]; 01805 rdata[ii + 1] += gg * dt[1]; 01806 norm[ii] += gg; 01807 } 01808 } 01809 } 01810 01811 if (x0 <= 2) { 01812 xx = -xx; 01813 yy = -(yy - ny / 2) + ny / 2; 01814 zz = -(zz - nz / 2) + nz / 2; 01815 x0 = 2 * (int) floor(xx + .5); 01816 y0 = (int) floor(yy + .5); 01817 z0 = (int) floor(zz + .5); 01818 01819 if (y0 > ny - 3 || z0 > nz - 3 || y0 < 2 || z0 < 2) 01820 break; 01821 01822 for (int k = z0 - 2; k <= z0 + 2; ++k) { 01823 for (int j = y0 - 2; j <= y0 + 2; ++j) { 01824 for (int i = 0; i <= x0 + 4; i += 2) { 01825 size_t ii = i + j * nx + k * nxy; 01826 float r = sqrt(Util::hypot3sq((float) i / 2 - xx, (float) j - yy, 01827 (float) k - zz)); 01828 float gg = weight * Interp::hyperg(r); 01829 01830 rdata[ii] += gg * dt[0]; 01831 rdata[ii + 1] -= gg * dt[1]; 01832 norm[ii] += gg; 01833 } 01834 } 01835 } 01836 } 01837 break; 01838 } 01839 01840 } 01841 } 01842 01843 image->update(); 01844 tmp_data->update(); 01845 // slice->update(); 01846 01847 return 0; 01848 }
|
|
Definition at line 702 of file reconstructor.h. References EMAN::FactoryBase::params. Referenced by WienerFourierReconstructor(). 00703 { 00704 params["size"] = 0; 00705 params["mode"] = 2; 00706 params["weight"] = 1.0; 00707 params["use_weights"] = true; 00708 params["dlog"] = false; 00709 params["padratio"] = 1.0; 00710 }
|
|
Definition at line 675 of file reconstructor.h. References WienerFourierReconstructor(). 00676 { 00677 return new WienerFourierReconstructor(); 00678 }
|
|
|
|
|
Definition at line 694 of file reconstructor.h. Referenced by get_name(). |