PyTom: Align projections and reconstruct a tomogram
Prior to a 3D reconstruction projections must be aligned to a common coordinate system because despite the best efforts to keep the feature of interest in the center during acquisition movement is inevitable. Moreover, the electron optical lense system causes image rotations and even subtle focus-dependent magnification changes, which need to be compensated for.
Currently, PyTom only supports alignment by fiducials (gold beads). The positions
of markers in the projection images need to be provided in an EM-file. This
file can be obtained by interactive marker localization in the
EM packages. There is also a function in AV3 to convert IMOD
marker files to an EM file (av3_wimp2em.m).
The PyTom alignment function simply minimizes the residual error of the
coordinates projected from a 3D model and the oberserved coordinates in the
projections. In more detail, 3D coordinates of the marker model, the projection
shifts, and the image rotation is first computed using an analytical approach,
in which a constant rotation (=tilt axis) and magnification is assumed
function determines a tilt axis (to x-axis) between 0 and 180 degrees; by definition
for a tilt axis x there is a second solution, x+180 deg, that would satisfy the
marker equations equally well. You need prior knowledge to decide, which tilt
axis is correct - assuming the wrong tilt axis will yield reconstructions
with the inverted handedness. In the alignment it can be specified if the solution
180<x<360 deg is supposed to be correct rather than the default solution
To determine the 3D coordinate system of the markers the coordinates of one
marker point need to be defined. Thus, a reference marker needs to be chosen
irefmark). To this marker reference coordinates (
can be assigned or default values can be chosen. The default value is
y_iref are the coordinates of the marker in a reference projection
After approximate determination of the image rotations, translations, and 3D
coordinates of the markers an optimization algorithm
TiltAlignment.alignFromFiducials) determines the tilt-angle specific image
rotation, tilt-angle specific magnification, and again translations and 3D coordinates.
For the reference projection
iref the magnification is defined as 1.
The projections are aligned by inverse application of the determined translations, rotations, and magnifications to each projection. The transformation is carried out as a single interpolation in real space (default interpolation method: 3rd order spline).
If a tomogram were reconstructed from the unprocessed projections the low-frequency information would be artificially enhanced as illustrated in Fig.1. To obtain a truthful 3D image the information below the so-called Crowther frequency need to be attenuated. This is accomplished by weighting the projections in Fourier space proportional to its frequency perpendicular to the tilt axis. This can be approximated as a 'ramp' function (analytical weighting). A more accurate weighting scheme would be to explicitly compute the overlapping information in Fourier space (exact weighting, not yet implemented in PyTom). In many cases it makes sense to furthermore apply a low-pass filter to the data to eliminate noise, which is dominant in the higher frequencies.
Fig. 1.Sampling of Fourier space by projections in 2D. The Fourier transform of each projection samples a slice of Fourier space (black lines); its normal is determined by the tilt angle. The signal ‘leak’ of each sample point into the neighboring area is reciprocally proportional to the diameter D of the object. For reconstruction, the Fourier coefficients on the 3D grid (intersected lines) need to be approximated from the sampling points of the projections by an appropriate algorithm. In principle, the 3D signal can be reconstructed without any gaps to a resolution k_C if the entire tilt range of +/-90 deg was accessible. Restricting the tilting to +/-90 deg gives rise to a ‘missing wedge’.
In PyTom reconstructions are performed using backprojection. The algorithm can use different interpolation kernels. By default we use a 3rd order spline.
In our tutorial data-set you will find a directory
reconstructTomo where you can perform reconstruction of a full tomogram by weighted backprojection.
Projections will be aligned and weighted during this step so that the data you start with are unaligned and unweighted projections. Essentially raw projections only sorted according to the tilt angle.
Running this script will generate a full
tomogram.em file, either 2x binned (4x downscaled) (
512,512,128) or in original size (
reconstructTomo.sh is essentially the only script you have to run to obtain the aligned, weighted projections and the 2x binned tomogram.
reconstructTomogram.py --tiltSeriesPath ../projections/ \ --tiltSeriesPrefix tomo01_sorted \ --firstIndex 1 \ --lastIndex 41 \ --referenceIndex 21 \ --markerFile ../projections/markfile_temp.mark \ --referenceMarkerIndex 1 \ --projectionTargets ./alignedProjections/tomo01 \ --projectionBinning 4 \ --lowpassFilter 0.5 \ --tomogramFile tomogram.em \ --fileType em \ --tomogramSizeX 512 \ --tomogramSizeY 512 \ --tomogramSizeZ 128 \ --reconstructionCenterX 0 \ --reconstructionCenterY 0 \ --reconstructionCenterZ 0
reconstructTomogram.py, simply run
EMformat, .e.g., generated using TOM)
MRC- no tomogram written if not specified)
Now PyTom also includes the support of a new reconstruction algorithm called Iterative Nonuniform fast Fourier transform based Reconstruction method (INFR). For more details about this algorithm, please check: Iterative reconstruction of cryo-electron tomograms using nonuniform fast Fourier transforms, Y. Chen et al., JSB 2014.
The way to use it is simple:
pytom PathToPytom/bin/reconstructTomogram.py --tiltSeriesPath MyTiltSeriesDir \ --tiltSeriesPrefix NameOfTiltSeriesFiles --firstIndex 1 --lastIndex MyLastIndex \ --referenceIndex MyReferenceProjectionIndex --markerFile MyMarkerFile \ --referenceMarkerIndex MyReferenceMarkerIndex --projectionTargets \ MyFileNameOfAlignedProjections --projectionBinning MyBinning --weightingType 0
Sometimes you might want to do simple things on your reconstructed tomogram like reducing its size or low-pass filtering it without reconstructing it again. Here is an example how to reduce the size in z:
from pytom_volume import read, subvolume
v = read('tomogram.em')
#this function will cut out a subvolume the size 512,512,256 out of the original volume starting at position 0,0,128
s = subvolume(v,0,0,128,512,512,256)
Make sure you store the position at which the subvolume was cut. This will coordinate is required for reconstructing subtomograms at full resolution (see Reconstruct Subtomograms Tutorial).
If you want to lowpass filter the tomogram / subvolume to a specified resolution you can do that interactively as well:
$ipytom volume = read('subvolume.em')
from pytom.basic.filter import lowpassFilter
filteredVolume = lowpassFilter(volume,100,10) #filter the volume to the 100th pixel in Fourier space with 10 pixel smoothing