gridding.rst

gridding.rst (text/plain), 15.55 KB, created by Grant Goodyear (RETIRED) on 2006-01-26 15:52:39 UTC

========
Gridding
========

.. contents::

Gridding in EMAN2
=================

Gridding in EMAN2 currently shows up in 2-d rotations and 3-d projections.
At the moment, all of the gridding routines assume that a Kaiser-Bessel
window function is used.  Because gridding involves a lot of steps, 
high-order functionality exists as python code in Golem, while the
core procedures needed for gridding are written in C++.

Gridding projections using Golem
--------------------------------

Projections may be created (and written to disk) using code like the
following::

    e = getImage("volname.spi")
    anglelist = voea(delta_theta)
    create_write_gridprojections(e, anglelist, kb_K, kb_alpha)

where kb_K and kb_alpha are Kaiser-Bessel window parameters.  This 
snippet of code will generate a series of files named "proj****.tst"
(by default) that contain the projections of ``e`` along the Euler
angles generated by ``voea``.  The ``create_write_gridprojections``
code may be found in ``Golem/projection.py``.

Gridding rotations using Golem
------------------------------

Two-dimensional images may be rotated using gridding (and also subsequently
translated) via the ``grotshift2d`` function::

    e = getImage("image.spi")
    m = e.get_xsize()
    n = m*npad
    kb = Util.KaiserBessel(alpha, K, m/2, K/(2.*n), n)
    rotshift = grotshift2d(e, angle, kb, npad, dx, dy)

where ``kb`` is a KaiserBessel window function (which only needs to be created
once, not for each image to be rotated, if all of the images are the same
size), ``angle`` is the rotation angle (in radians), and ``dx`` & ``dy`` are
the translation amounts (in pixels).  The ``grotshift2d`` function may
be found in ``Golem/fundamentals.py``.  The translation itself does not
involve gridding; instead it breaks the translation into integer- and
fractional-pixel parts, uses a Fourier phase multiplication to handle
the fractional translation, and performs a cyclic integer shift in real
space for the integer translation.

Gridding C++ code in EMAN2
--------------------------

Projection
..........

Gridding projections of a volume are accomplished using the 
``FourierGriddingProjector`` in ``libEM/projector.{h,cpp}``.  Projectors 
in EMAN2 are a pain, since they are objects, but the important part 
of the code is the ``FourierGriddingProjector::project3d`` method
in ``projector.cpp``.  That code calls ``EMData::extractplane``
(which is in ``libEM/emdata.{h,cpp}``) to do the Fourier convolutions.  

Rotation
........

Gridding rotations are accomplished using ``EMData::fouriergridrot2d``,
which may be found in ``libEM/emdata.{h,cpp}``.  The ``fouriergridrot2d``
code itself is a simple rotation code, while the actual gridding
convolution is found in ``EMData::extractpoint`` (also in 
``libEM/emdata.{h,cpp}``).


Kaiser-Bessel window
....................

A Kaiser-Bessel window function is actually a class in EMAN2.  Making it
a class comes in handy, since one can then have separate K-B objects that
have been defined with different values for ``alpha`` or ``K``, for example.
The code is the ``KaiserBessel`` class in ``libEM/util.{h,cpp}``.  See
the example in `Gridding Rotations using Golem`_ to see how a Kaiser-Bessel
window object may be created in Python.


70S
===

Our intitial testing was performed with the 70S ribosome using the 
phenix gui with some home-written "tasks" that access either EMAN2 or
spider routines.

Initial data
------------

The 70S data we use, `model.tcp`_, is a 75×75×75 volume that I obtained
from Zhong.  Zhong tells me that the pixel size of this volume is 4.78Å.
An advantage of this volume is that 2-d projections of this volume are
very distinctive.  A drawback is the fact that the dynamic range of this
volume is very small, ranging only from 0.0 to 3.0.

.. _model.tcp: data/70S/model.tcp

Noisy volume
------------

We also created a `noisy version`_ of the model using Golem::

    e = getImage('model.tcp')
    m = e.get_xsize()
    e += model_gauss_noise(10.5, m, m, m)
    dropImage(e, 'noise.tcp')

(We also created a version with a sigma of `20.5`_ pixels.)

.. _noisy version: data/70S/noise.tcp
.. _20.5: data/70S/noise20.tcp

Of course, it isn't clear that there is really a good reason to do this,
since what one expects to have is noisy projections, not a noisy volume.

Phenix testing strategies
-------------------------

We used the phenix_ gui to perform much of our testing.  We had two
nearly-identical phenix strategies, `grid_recon_fsc`_ and `interp_recon_fsc`_
(see the figure below), where each strategy computes a set of projection Euler
angles using EMAN2 (``voea``), generates a spider-acceptable doc file
containing those angles (``anglist2doc``), generates 2-d projections along
those angles in EMAN2 using either gridding (``gridproj``) or interpolation
(``interpproj``), performs a 3-d reconstruction using spider either with
gridding (``spidbprg``) or interpolation (``spidbp3d``), and finally
computes and writes out the FSC of the original model and the reconstructed
volume (``fsc``).

.. _phenix: phenix.html
.. _grid_recon_fsc: data/70S/grid_recon_fsc...wxGUI.StrategyClass.StrategyClass
.. _interp_recon_fsc: data/70S/interp_recon_fsc...wxGUI.StrategyClass.StrategyClass


.. figure:: data/70S/phenix_grid.png
   :align: center

   Phenix strategy for a gridding-based reconstruction.  The 
   interpolation-based strategy looks the same, except that the
   projection task is ``interpproj`` instead of ``gridproj``, and
   the reconstruction task is ``spidbp3d`` instead of ``spidbprg``.


Testing
-------

Random rotation/translation tests
.................................

For a 75×75×75 volume it is easy enough to compute a rough estimate of
the number of projections needed to reconstruct the original volume::

    [4πr³/3]/[πr²] = 4r/3 ≃ 49 projections

so for a bit of overkill we used 195 projections in our reconstructions
(corresponding to δθ=10˚ in voea).  In the gridding reconstructions
during the projection stage we randomly rotate (also using gridding)
and translate the images, and then translate and rotate them back, in order to
simulate the steps involved in an alignment.  

=== ======= ============================ ==================-=====
G/I Noise σ fsc file                     reconstructed volume
=== ======= ============================ =======================
G   10.5    `fsc_grid_noise_195.out0`_   `grid_noise_195_0.tst`_
I   10.5    `fsc_interp_noise_195.out0`_ `interp_noise_195_0.tst`_
I   --      `fsc_interp_195.out0`_       `interp_195_0.tst`_
G   --      `fsc_grid_195.out0`_         `grid_195_0.tst`_
G   20.5    `fsc_grid_noise20_195.out0`_ `grid_noise20_195_0.tst`_
=== ======= ============================ =======================

.. _fsc_grid_noise_195.out0: data/70S/fsc_grid_noise_195.out0
.. _fsc_interp_noise_195.out0: data/70S/fsc_interp_noise_195.out0
.. _fsc_interp_195.out0: data/70S/fsc_interp_195.out0
.. _fsc_grid_195.out0: data/70S/fsc_grid_195.out0
.. _fsc_grid_noise20_195.out0: data/70S/fsc_grid_noise20_195.out0
.. _grid_noise_195.tst: data/70S/grid_noise_195.tst
.. _interp_noise_195.tst: data/70S/interp_noise_195.tst
.. _interp_195.tst: data/70S/interp_195.tst
.. _grid_195.tst: data/70S/grid_195.tst
.. _grid_noise20_195.tst: data/70S/grid_noise20_195.tst

Single-angle 2-d tests
......................

After running the reconstruction tests, we dropped back to double-check
2-d rotation performance using simple python scripts `rot.py`_ (interpolation)
and `gridrot.py`_ (gridding) with input images `proj0001.tcp`_ and
`noise0001.tcp`_.  The latter, noisy image was created by::

    e = getImage('proj0001.tcp')
    m = e.get_xsize()
    e += model_gauss_noise(10.5, m, m)
    dropImage(e, 'noise0001.tcp')

.. _rot.py: data/70S/rot.py
.. _gridrot.py: data/70S/gridrot.py
.. _proj0001.tcp: data/70S/proj0001.tcp
.. _noise0001.tcp: data/70S/noise0001.tcp

===  ===  =======  ============== ============== ==================
∠    G/I  noise σ  fsc file       rot file       rotback file
===  ===  =======  ============== ============== ==================
π/6  G    --       `g30fsc.out`_  `g30rot.tcp`_  `g30rotback.tcp`_
π/6  I    --       `i30fsc.out`_  `i30rot.tcp`_  `i30rotback.tcp`_
π/6  G    10.5     `gn30fsc.out`_ `gn30rot.tcp`_ `gn30rotback.tcp`_
π/6  I    10.5     `in30fsc.out`_ `in30rot.tcp`_ `in30rotback.tcp`_
π/4  G    --       `g45fsc.out`_  `g45rot.tcp`_  `g45rotback.tcp`_
π/4  I    --       `i45fsc.out`_  `i45rot.tcp`_  `i45rotback.tcp`_
π/4  G    10.5     `gn45fsc.out`_ `gn45rot.tcp`_ `gn45rotback.tcp`_
π/4  I    10.5     `in45fsc.out`_ `in45rot.tcp`_ `in45rotback.tcp`_
===  ===  =======  ============== ============== ==================

.. _g30fsc.out: data/70S/g30fsc.out
.. _g30rot.tcp: data/70S/g30rot.tcp
.. _g30rotback.tcp: data/70S/g30rotback.tcp
.. _i30fsc.out: data/70S/i30fsc.out
.. _i30rot.tcp: data/70S/i30rot.tcp
.. _i30rotback.tcp: data/70S/i30rotback.tcp
.. _gn30fsc.out: data/70S/gn30fsc.out
.. _gn30rot.tcp: data/70S/gn30rot.tcp
.. _gn30rotback.tcp: data/70S/gn30rotback.tcp
.. _in30fsc.out: data/70S/in30fsc.out
.. _in30rot.tcp: data/70S/in30rot.tcp
.. _in30rotback.tcp: data/70S/in30rotback.tcp
.. _g45fsc.out: data/70S/g45fsc.out
.. _g45rot.tcp: data/70S/g45rot.tcp
.. _g45rotback.tcp: data/70S/g45rotback.tcp
.. _i45fsc.out: data/70S/i45fsc.out
.. _i45rot.tcp: data/70S/i45rot.tcp
.. _i45rotback.tcp: data/70S/i45rotback.tcp
.. _gn45fsc.out: data/70S/gn45fsc.out
.. _gn45rot.tcp: data/70S/gn45rot.tcp
.. _gn45rotback.tcp: data/70S/gn45rotback.tcp
.. _in45fsc.out: data/70S/in45fsc.out
.. _in45rot.tcp: data/70S/in45rot.tcp
.. _in45rotback.tcp: data/70S/in45rotback.tcp

.. figure:: data/70S/rotfsc.png
   :align: center

   FRC curves comparing the original and rotated (45° forwards and back)
   images.  The no-noise gridding result is essentially 1 for frequencies
   less than 0.45.

GroEL
=====

Initial data
------------

Zhong used a `GroEL pdb`_ file (and spider script `b03.gen`_) to generate 
both `1Å`_ and `0.5Å`_ spider volumes.

.. _GroEL pdb: data/GroEL/GroEL.pdb
.. _b03.gen: data/GroEL/b03.gen
.. _1Å: data/GroEL/GroEL.grl
.. _0.5Å: data/GroEL/GroEL_pt5.grl

Ultimately, the 0.5Å data was low-pass filtered in spider::

    fq np
    GroEL_pt5
    GroEL_filt
    7
    0.23,0.27

and then symmetrized (using the spider ``@symv``
procedure) and decimated ``(2,2,2)`` to generate `GroEL_symdc.grl`_.

.. _GroEL_symdc.grl: data/GroEL/GroEL_symdc.grl

Projections
-----------

Two sets of projections were computed from `GroEL_symdc.grl`_: one using
spider's ``pj 3q`` command (which uses interpolation), and another using
spider's ``pj rg`` command, which uses gridding.  In both cases, a set of 
angles (`angro.grl`_) was required.  This set of angles was generated with
considerable pain.  A symmetric-subunit set (GroEL has 7-fold symmetry) of
projection Euler angles was first constructed using spider's ``vo ea`` command
for the following (with all angles in degrees):

=======   ==    =======     ===========
θ range   δθ    φ range     # of angles
=======   ==    =======     ===========
0-10       1    0-51.42     38
11-75      5    0-51.42     63
76-89      2    0-51.42     167
90         1    0-25.71     24
=======   ==    =======     ===========

for a total of 292 Euler angles.  A single doc file (`angles.grl`_) containing
all of these angles was constructed using ``vo ea`` and ``doc merge``.  
Then ``vo ra`` and ``doc merge`` were used to generate a doc file
(`angro.grl`_) containing Euler angles corresponding to all seven symmetric
subunits (using the angles psi, theta, phi in the doc file `syms.grl`_).

.. _angles.grl: data/GroEL/angles.grl
.. _angro.grl:  data/GroEL/angro.grl
.. _syms.grl:   data/GroEL/syms.grl

Reconstructions
---------------

Both sets of projections were used in reconstructions both with
and without gridding using spider.  (Some of the reconstructions
were run by hand, but many were run using the batch script
`b01.grl`_.)  The table below lists the various reconstructions that
were computed, including the command used to generate the projections,
the command used to generate the reconstructions, whether or not the
reconstruction imposed the 7-fold symmetry of `syms.grl`_, the 
file containing the reconstructed volume, and the file containing
the FSC curve (generated by the spider ``rf 3`` command) comparing the
original and reconstructed volumes.

.. _b01.grl: data/GroEL/b01.grl

=========== ============== ======== ====================== ====================
Projections Reconstruction Symmetry Volume                 rf 3 file
=========== ============== ======== ====================== ====================
pj 3q       bp 3d          no                              `rf3.grl`_
pj 3q       bp 3f          no       `GroEL_bp3f.grl`_      `rf3_bp3f.grl`_
pj 3q       bp 3f          yes      `GroEL_bp3fsym.grl`_   `rf3_bp3fsym.grl`_
pj 3q       bp rg          no       `GroEL_bprg.grl`_      `rf3_bprg.grl`_
pj 3q       bp 3d diam 256 no       `GroEL_bp3ddiam.grl`_  `rf3_bp3ddiam.grl`_
pj rg       bp 3f          no       `GroEL_gbp3f.grl`_     `rf3_gbp3f.grl`_
pj rg       bp 3f          yes      `GroEL_gbp3fsym.grl`_  `rf3_gbp3fsym.grl`_
pj rg       bp rg          no       `GroEL_gbprg.grl`_     `rf3_gbprg.grl`_
pj rg       bp 3d diam 256 no       `GroEL_gbp3ddiam.grl`_ `rf3_gbp3ddiam.grl`_
=========== ============== ======== ====================== ====================

.. _rf3.grl:   data/GroEL/rf3.grl
.. _GroEL_bp3f.grl: data/GroEL/GroEL_bp3f.grl
.. _rf3_bp3f.grl: data/GroEL/rf3_bp3f.grl
.. _GroEL_bp3fsym.grl: data/GroEL/GroEL_bp3fsym.grl
.. _rf3_bp3fsym.grl: data/GroEL/rf3_bp3fsym.grl
.. _GroEL_bprg.grl: data/GroEL/GroEL_bprg.grl
.. _rf3_bprg.grl: data/GroEL/rf3_bprg.grl
.. _GroEL_bp3ddiam.grl: data/GroEL/GroEL_bp3ddiam.grl
.. _rf3_bp3ddiam.grl: data/GroEL/rf3_bp3ddiam.grl
.. _GroEL_gbp3f.grl: data/GroEL/GroEL_gbp3f.grl
.. _rf3_gbp3f.grl: data/GroEL/rf3_gbp3f.grl
.. _GroEL_gbp3fsym.grl: data/GroEL/GroEL_gbp3fsym.grl
.. _rf3_gbp3fsym.grl: data/GroEL/rf3_gbp3fsym.grl
.. _GroEL_gbprg.grl: data/GroEL/GroEL_gbprg.grl
.. _rf3_gbprg.grl: data/GroEL/rf3_gbprg.grl
.. _GroEL_gbp3ddiam.grl: data/GroEL/GroEL_gbp3ddiam.grl
.. _rf3_gbp3ddiam.grl: data/GroEL/rf3_gbp3ddiam.grl

An example of using ``rf 3`` is ::

    rf 3
    GroEL_symdc
    GroEL_gbprg
    1.
    0.1,1
    c
    90
    3
    rf3_gbprg

Looking at the FSC curves:

.. image:: data/GroEL/rf3comp.png

It is clear that gridding notably improves the results, at least by this
measure, although I'm not sure that we understand the funny dip in the
gridding curves at low frequencies.

Visualization
-------------

Okay, at least in principle gridding works great.  To sell it, though, it 
would be nice if one could actually see manifestly improved results in
the reconstructed volumes.  For some reason the reconstructed volumes
all seem to have different dynamic ranges, so we ran the spider ``fv``
command to determine a consistent threshold for each of the images.  
Through trial and error it was determined that a volume of |80^3| voxels 
sufficed as the input volume for the ``fv`` command:

.. |80^3| replace:: 80\ :sup:`3`

====================== =========
volume file            threshold
====================== =========
`GroEL_symdc.grl`_     0.0338
`GroEL_bp3f.grl`_      0.0366
`GroEL_gbp3f.grl`_     0.0357
`GroEL_gbprg.grl`_     0.0187
`GroEL_bprg.grl`_      0.0192
`GroEL_gbp3ddiam.grl`_ 26.4
`GroEL_bp3ddiam.grl`_  27.1
====================== =========

We used ``chimera`` to look at the various reconstructed volumes, but we 
couldn't see any noticeable improvement.  In fact, we also pulled a single
helix out of the GroEL pdb file, thinking that perhaps GroEL was too big
to digest, and the resulting volume was rather blobby and not something
that obviously looked like a helix, so we've put this project on hold for
now.