EMAN2
Public Member Functions | Static Public Member Functions | Static Public Attributes | Private Member Functions | Private Attributes
EMAN::GaussFFTProjector Class Reference

Gaussian FFT 3D projection. More...

#include <projector.h>

Inheritance diagram for EMAN::GaussFFTProjector:
Inheritance graph
[legend]
Collaboration diagram for EMAN::GaussFFTProjector:
Collaboration graph
[legend]

List of all members.

Public Member Functions

 GaussFFTProjector ()
EMDataproject3d (EMData *image) const
 Project an 3D image into a 2D image.
EMDatabackproject3d (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 ProjectorNEW ()

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

Detailed Description

Gaussian FFT 3D projection.

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.


Constructor & Destructor Documentation

EMAN::GaussFFTProjector::GaussFFTProjector ( ) [inline]

Definition at line 153 of file projector.h.

Referenced by NEW().

                                   :alt(0), az(0), phi(0)
                {
                }

Member Function Documentation

EMData * GaussFFTProjector::backproject3d ( EMData image) const [virtual]

Back-project a 2D image into a 3D image.

Returns:
A 3D image from the backprojection.

Implements EMAN::Projector.

Definition at line 1864 of file projector.cpp.

{
    // no implementation yet
    EMData *ret = new EMData();
    return ret;
}
string EMAN::GaussFFTProjector::get_desc ( ) const [inline, virtual]

Implements EMAN::Projector.

Definition at line 175 of file projector.h.

                {
                        return "Projections using a Gaussian kernel in Fourier space. Produces artifacts, not recommended.";
                }
string EMAN::GaussFFTProjector::get_name ( ) const [inline, virtual]

Get the projector's name.

Each projector is indentified by unique name.

Returns:
The projector's name.

Implements EMAN::Projector.

Definition at line 170 of file projector.h.

References NAME.

                {
                        return NAME;
                }
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.

Returns:
A dictionary containing the parameter info.

Reimplemented from EMAN::Projector.

Definition at line 185 of file projector.h.

References EMAN::EMObject::INT, EMAN::TypeDict::put(), and EMAN::EMObject::TRANSFORM.

                {
                        TypeDict d;
                        d.put("transform", EMObject::TRANSFORM);
                        d.put("mode", EMObject::INT);
                        return d;
                }
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::EMData::get_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), nx, ny, rdata, EMAN::EMData::setup4slice(), and sqrt().

{
        float *rdata = image->get_data();
        int nx = image->get_xsize();
        int ny = image->get_ysize();
        int nz = image->get_zsize();

        if ( mode == 0 ) mode = 2;

        if (mode == 1) {
                int x0 = 2 * (int) floor(x + 0.5f);
                int y0 = (int) floor(y + 0.5f);
                int z0 = (int) floor(z + 0.5f);

                data[0] = rdata[x0 + y0 * nx + z0 * nx * ny];
                data[1] = rdata[x0 + y0 * nx + z0 * nx * ny + 1];
                return true;
        }
        else if (mode == 2) {
                int x0 = (int) floor(x);
                int y0 = (int) floor(y);
                int z0 = (int) floor(z);

                float dx = x - x0;
                float dy = y - y0;
                float dz = z - z0;

                if (x0 <= nx/2- 2 && y0 <= ny - 1 && z0 <= nz - 1) {
                        int i = (int) (x0 * 2 + y0 * nx + z0 * nx * ny);

                        float n = Util::agauss(1, dx, dy, dz, gw) +
                                Util::agauss(1, 1 - dx, dy, dz, gw) +
                                Util::agauss(1, dx, 1 - dy, dz, gw) +
                                Util::agauss(1, 1 - dx, 1 - dy, dz, gw) +
                                Util::agauss(1, dx, dy, 1 - dz, gw) +
                                Util::agauss(1, 1 - dx, dy, 1 - dz, gw) +
                                Util::agauss(1, dx, 1 - dy, 1 - dz, gw) + Util::agauss(1, 1 - dx, 1 - dy, 1 - dz, gw);

                        data[0] = Util::agauss(rdata[i], dx, dy, dz, gw) +
                                Util::agauss(rdata[i + 2], 1 - dx, dy, dz, gw) +
                                Util::agauss(rdata[i + nx], dx, 1 - dy, dz, gw) +
                                Util::agauss(rdata[i + nx + 2], 1 - dx, 1 - dy, dz, gw) +
                                Util::agauss(rdata[i + nx * ny], dx, dy, 1 - dz, gw) +
                                Util::agauss(rdata[i + 2 + nx * ny], 1 - dx, dy, 1 - dz, gw) +
                                Util::agauss(rdata[i + nx + nx * ny], dx, 1 - dy, 1 - dz, gw) +
                                Util::agauss(rdata[i + 2 + nx + nx * ny], 1 - dx, 1 - dy, 1 - dz, gw) / n;

                        i++;

                        data[1] = Util::agauss(rdata[i], dx, dy, dz, gw) +
                                Util::agauss(rdata[i + 2], 1 - dx, dy, dz, gw) +
                                Util::agauss(rdata[i + nx], dx, 1 - dy, dz, gw) +
                                Util::agauss(rdata[i + nx + 2], 1 - dx, 1 - dy, dz, gw) +
                                Util::agauss(rdata[i + nx * ny], dx, dy, 1 - dz, gw) +
                                Util::agauss(rdata[i + 2 + nx * ny], 1 - dx, dy, 1 - dz, gw) +
                                Util::agauss(rdata[i + nx + nx * ny], dx, 1 - dy, 1 - dz, gw) +
                                Util::agauss(rdata[i + 2 + nx + nx * ny], 1 - dx, 1 - dy, 1 - dz, gw) / n;
                        return true;
                }
                return false;
        }
        else if (mode == 3) {
                int x0 = 2 * (int) floor(x + .5);
                int y0 = (int) floor(y + .5);
                int z0 = (int) floor(z + .5);

                float *supp = image->setup4slice();

                if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
                        float n = 0;

                        if (x0 == 0) {
                                x0 += 4;
                                size_t idx;
                                float r, g;
                                for (int k = z0 - 1; k <= z0 + 1; k++) {
                                        for (int j = y0 - 1; j <= y0 + 1; j++) {
                                                for (int i = x0 - 2; i <= x0 + 2; i += 2) {
                                                        r = Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z);
                                                        g = exp(-r / gw);
                                                        n += g;
                                                        idx = i + j * 12 + (size_t)k * 12 * ny;
                                                        data[0] += supp[idx] * g;
                                                        data[1] += supp[idx + 1] * g;
                                                }
                                        }
                                }
                        }
                        else {
                                size_t idx;
                                float r, g;
                                for (int k = z0 - 1; k <= z0 + 1; k++) {
                                        for (int j = y0 - 1; j <= y0 + 1; j++) {
                                                for (int i = x0 - 2; i <= x0 + 2; i += 2) {
                                                        r = Util::hypot3sq(i / 2.0f - x, j - y, k - z);
                                                        g = exp(-r / gw);
                                                        n += g;
                                                        idx = i + j * nx + (size_t)k * nx * ny;
                                                        data[0] += rdata[idx] * g;
                                                        data[1] += rdata[idx + 1] * g;
                                                }
                                        }
                                }
                        }

                        data[0] /= n;
                        data[1] /= n;
                        return true;
                }
                return false;
        }
        else if (mode == 4) {
                int x0 = 2 * (int) floor(x);
                int y0 = (int) floor(y);
                int z0 = (int) floor(z);

                float *supp = image->setup4slice();

                if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
                        float n = 0;

                        if (x0 == 0) {
                                x0 += 4;
                                size_t idx;
                                float r, g;
                                for (int k = z0 - 1; k <= z0 + 2; k++) {
                                        for (int j = y0 - 1; j <= y0 + 2; j++) {
                                                for (int i = x0 - 2; i <= x0 + 4; i += 2) {
                                                        r = Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z);
                                                        g = exp(-r / gw);
                                                        n += g;
                                                        idx = i + j * 12 + (size_t)k * 12 * ny;
                                                        data[0] += supp[idx] * g;
                                                        data[1] += supp[idx + 1] * g;
                                                }
                                        }
                                }
                        }
                        else {
                                float r, g;
                                size_t idx;
                                for (int k = z0 - 1; k <= z0 + 2; k++) {
                                        for (int j = y0 - 1; j <= y0 + 2; j++) {
                                                for (int i = x0 - 2; i <= x0 + 4; i += 2) {
                                                        r = Util::hypot3sq(i / 2.0f - x, j - y, k - z);
                                                        g = exp(-r / gw);
                                                        n += g;
                                                        idx = i + j * nx + (size_t)k * nx * ny;
                                                        data[0] += rdata[idx] * g;
                                                        data[1] += rdata[idx + 1] * g;
                                                }
                                        }
                                }
                        }
                        data[0] /= n;
                        data[1] /= n;
                        return true;
                }
                return false;
        }
        else if (mode == 5) {
                int x0 = 2 * (int) floor(x + .5);
                int y0 = (int) floor(y + .5);
                int z0 = (int) floor(z + .5);

                float *supp = image->setup4slice();
                float *gimx = Interp::get_gimx();

                if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
                        int mx0 = -(int) floor((x - x0 / 2) * 39.0 + .5) - 78;
                        int my0 = -(int) floor((y - y0) * 39.0 + .5) - 78;
                        int mz0 = -(int) floor((z - z0) * 39.0 + .5) - 78;

                        float n = 0;
                        int mmz = mz0;
                        int mmy = my0;
                        int mmx = mx0;

                        if (x0 < 4) {
                                x0 += 4;

                                size_t idx;
                                float g;
                                for (int k = z0 - 2; k <= z0 + 2; k++, mmz += 39) {
                                        for (int j = y0 - 2; j <= y0 + 2; j++, mmy += 39) {
                                                for (int i = x0 - 4; i <= x0 + 4; i += 2, mmx += 39) {
                                                        g = gimx[abs(mmx) + abs(mmy) * 100 + abs(mmz) * 10000];
                                                        n += g;
                                                        idx = i + j * 12 + (size_t)k * 12 * ny;
                                                        data[0] += supp[idx] * g;
                                                        data[1] += supp[idx + 1] * g;
                                                }
                                        }
                                }
                        }
                        else {
                                size_t ii;
                                float g;
                                for (int k = z0 - 2; k <= z0 + 2; k++, mmz += 39) {
                                        for (int j = y0 - 2; j <= y0 + 2; j++, mmy += 39) {
                                                for (int i = x0 - 4; i <= x0 + 4; i += 2, mmx += 39) {
                                                        ii = i + j * nx + (size_t)k * nx * ny;
                                                        g = gimx[abs(mmx) + abs(mmy) * 100 + abs(mmz) * 10000];
                                                        n += g;
                                                        data[0] += rdata[ii] * g;
                                                        data[1] += rdata[ii + 1] * g;
                                                }
                                        }
                                }
                        }

                        data[0] /= n;
                        data[1] /= n;
                        return true;
                }
                return false;
        }
        else if (mode == 6) {

                int x0 = 2 * (int) floor(x + .5);
                int y0 = (int) floor(y + .5);
                int z0 = (int) floor(z + .5);

                float *supp = image->setup4slice();

                if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
                        float n = 0;

                        if (x0 < 4) {
                                x0 += 4;
                                float r, g;
                                size_t idx;
                                for (int k = z0 - 2; k <= z0 + 2; k++) {
                                        for (int j = y0 - 2; j <= y0 + 2; j++) {
                                                for (int i = x0 - 4; i <= x0 + 4; i += 2) {
                                                        r = Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z);
                                                        g = exp(-r / gw);
                                                        n += g;
                                                        idx = i + j * 12 + (size_t)k * 12 * ny;
                                                        data[0] += supp[idx] * g;
                                                        data[1] += supp[idx + 1] * g;
                                                }
                                        }
                                }
                        }
                        else {
                                float r, g;
                                size_t idx;
                                for (int k = z0 - 2; k <= z0 + 2; k++) {
                                        for (int j = y0 - 2; j <= y0 + 2; j++) {
                                                for (int i = x0 - 4; i <= x0 + 4; i += 2) {
                                                        r = Util::hypot3sq(i / 2.0f - x, j - y, k - z);
                                                        g = exp(-r / gw);
                                                        n += g;
                                                        idx = i + j * nx + (size_t)k * nx * ny;
                                                        data[0] += rdata[idx] * g;
                                                        data[1] += rdata[idx + 1] * g;
                                                }

                                        }
                                }
                        }

                        data[0] /= n;
                        data[1] /= n;

                        return true;
                }
                return false;
        }
        else if (mode == 7) {
                int x0 = 2 * (int) floor(x + .5);
                int y0 = (int) floor(y + .5);
                int z0 = (int) floor(z + .5);

                float *supp = image->setup4slice();

                if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
                        float n = 0;
                        if (x0 < 4) {
                                x0 += 4;
                                float r, g;
                                size_t idx;
                                for (int k = z0 - 2; k <= z0 + 2; k++) {
                                        for (int j = y0 - 2; j <= y0 + 2; j++) {
                                                for (int i = x0 - 4; i <= x0 + 4; i += 2) {
                                                        r = sqrt(Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z));
                                                        g = Interp::hyperg(r);
                                                        n += g;
                                                        idx = i + j * 12 + (size_t)k * 12 * ny;
                                                        data[0] += supp[idx] * g;
                                                        data[1] += supp[idx + 1] * g;
                                                }
                                        }
                                }
                        }
                        else {
                                float r, g;
                                size_t idx;
                                for (int k = z0 - 2; k <= z0 + 2; k++) {
                                        for (int j = y0 - 2; j <= y0 + 2; j++) {
                                                for (int i = x0 - 4; i <= x0 + 4; i += 2) {
                                                        r = sqrt(Util::hypot3sq(i / 2.0f - x, j - y, k - z));
                                                        g = Interp::hyperg(r);
                                                        n += g;
                                                        idx = i + j * nx + (size_t)k * nx * ny;
                                                        data[0] += rdata[idx] * g;
                                                        data[1] += rdata[idx + 1] * g;
                                                }

                                        }
                                }
                        }
                        data[0] /= n;
                        data[1] /= n;
                        return true;
                }
                return false;

        }
//      throw InvalidParameterException("Error, unkown mode");
        return false;
}
static Projector* EMAN::GaussFFTProjector::NEW ( ) [inline, static]

Definition at line 180 of file projector.h.

References GaussFFTProjector().

                {
                        return new GaussFFTProjector();
                }
EMData * GaussFFTProjector::project3d ( EMData image) const [virtual]

Project an 3D image into a 2D image.

Returns:
A 2D image from the projection.

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(), ImageDimensionException, 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(), x, and y.

{
        Transform* t3d = params["transform"];
        if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
        if ( image->get_ndim() != 3 ) throw ImageDimensionException("Error, the projection volume must be 3D");

        EMData *f = image;
        if (!image->is_complex()) {
                image->process_inplace("xform.phaseorigin.tocorner");
                f = image->do_fft();
                f->process_inplace("xform.fourierorigin.tocenter");
                image->process_inplace("xform.phaseorigin.tocenter");
        }

        int f_nx = f->get_xsize();
        int f_ny = f->get_ysize();
        int f_nz = f->get_zsize();

        if (!f->is_complex() || f_nz != f_ny || f_nx != f_ny + 2) {
                LOGERR("Cannot project this image");
                return 0;
        }

        f->ap2ri();

        EMData *tmp = new EMData();
        tmp->set_size(f_nx, f_ny, 1);
        tmp->set_complex(true);
        tmp->set_ri(true);

        float *data = tmp->get_data();

        Transform r = t3d->get_rotation_transform();
        r.invert();
        float scale = t3d->get_scale();

        int mode = params["mode"];
        float gauss_width = 1;
        if ( mode == 0 ) mode = 2;
        if (mode == 2 ) {
                gauss_width = EMConsts::I2G;
        }
        else if (mode == 3) {
                gauss_width = EMConsts::I3G;
        }
        else if (mode == 4) {
                gauss_width = EMConsts::I4G;
        }
        else if (mode == 6 || mode == 7) {
                gauss_width = EMConsts::I5G;
        }

        for (int y = 0; y < f_ny; y++) {
                for (int x = 0; x < f_nx / 2; x++) {
                        int ii = x * 2 + y * f_nx;
#ifdef  _WIN32
                        if (_hypot(x, y - f_ny / 2) >= f_ny / 2 - 1) {
#else
                        if (hypot(x, y - f_ny / 2) >= f_ny / 2 - 1) {
#endif  //_WIN32
                                data[ii] = 0;
                                data[ii + 1] = 0;
                        }
                        else {
                                Vec3f coord(x,(y - f_ny / 2),0);
                                Vec3f soln = r*coord;
                                float xx = soln[0];
                                float yy = soln[1];
                                float zz = soln[2];

                                int cc = 1;
                                if (xx < 0) {
                                        xx = -xx;
                                        yy = -yy;
                                        zz = -zz;
                                        cc = -1;
                                }

                                if (scale != 1.0) {
                                        xx *= scale;
                                        yy *= scale;
                                        zz *= scale;
                                }

                                yy += f_ny / 2;
                                zz += f_nz / 2;

                                if (yy < 0 || xx < 0 || zz < 0) {
                                        data[ii] = 0;
                                        data[ii+1] = 0;
                                        continue;
                                }

                                if (interp_ft_3d(mode, f, xx, yy, zz, data + ii, gauss_width)) {
                                        data[ii + 1] *= cc;
                                } else {
                                        data[ii] = 0;
                                        data[ii+1] = 0;
                                }
                        }
                }
        }

        f->update();
        tmp->update();

        tmp->process_inplace("xform.fourierorigin.tocorner");
        EMData *ret = tmp->do_ift();
        ret->process_inplace("xform.phaseorigin.tocenter");

        ret->translate(t3d->get_trans());

        if (t3d->get_mirror() ) ret->process_inplace("xform.flip",Dict("axis","x"));

        Dict filter_d;
        filter_d["gauss_width"] = gauss_width;
        filter_d["ring_width"] = ret->get_xsize() / 2;
        ret->process_inplace("math.gausskernelfix", filter_d);

        if( tmp )
        {
                delete tmp;
                tmp = 0;
        }

        ret->set_attr("xform.projection",t3d);
        ret->update();

        if(t3d) {delete t3d; t3d=0;}

        return ret;
}
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, and phi.

                {
                        Projector::set_params(new_params);
                        alt = params["alt"];
                        az = params["az"];
                        phi = params["phi"];
                }

Member Data Documentation

Definition at line 196 of file projector.h.

Referenced by set_params().

float EMAN::GaussFFTProjector::az [private]

Definition at line 196 of file projector.h.

Referenced by set_params().

const string GaussFFTProjector::NAME = "gauss_fft" [static]

Definition at line 193 of file projector.h.

Referenced by get_name().

Definition at line 196 of file projector.h.

Referenced by set_params().


The documentation for this class was generated from the following files: