PyTom: Align projections and reconstruct a tomogram


Pytom currently supports the alignment of projections by fiducials and reconstruction by weighted backprojection. Both, alignment by gold beads and reonstruction by weighted backprojection are still the most common approaches in the field. Nevertheless, we are working on the implementation of further alignment and reconstruction methods.

Theoretical Background

Alignment of tilt series

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 TOM/AV3 or 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 (TiltAlignment.computeCoarseAlignment). This 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 0<x<180 (handflip=True).

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 (r) can be assigned or default values can be chosen. The default value is r=(x_iref-center,y_iref-center,0). x_iref and y_iref are the coordinates of the marker in a reference projection iref.

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).

Weighting of projections

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’.

3D Reconstruction of tomogram

In PyTom reconstructions are performed using backprojection. The algorithm can use different interpolation kernels. By default we use a 3rd order spline.

Reconstruct a tomogram using

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 (2048,2048,512).
The is essentially the only script you have to run to obtain the aligned, weighted projections and the 2x binned tomogram.	--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 
In order to get acquainted with the parameters of, simply run --help
In detail the provided parameters for reconstruction are: All further steps shown in the tutorial are based on particles detected within this tomogram.

Reconstruct a tomogram using INFR algorithm

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/reconstruction/
It requires the following parameters: Note that INFR requires a lot of memory. Currently it is limited to 1k x 1k projections on most computers (i.e., depending on the amount of memory you have).

To create the aligned projections you can also use PyTom. Simply use the script for WBP reconstruction, but ommit weighting and reconstruction (--weightingType 0):
    pytom PathToPytom/bin/ --tiltSeriesPath MyTiltSeriesDir \
        --tiltSeriesPrefix NameOfTiltSeriesFiles --firstIndex 1 --lastIndex MyLastIndex \
        --referenceIndex MyReferenceProjectionIndex --markerFile MyMarkerFile \
        --referenceMarkerIndex MyReferenceMarkerIndex --projectionTargets \
        MyFileNameOfAlignedProjections --projectionBinning MyBinning --weightingType 0

Basic interactive manipulation of volumes

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)[0] #filter the volume to the 100th pixel in Fourier space with 10 pixel smoothing