Xmipp

spARMA (Spectral Estimation by means of ARMA model)

Purpose

This program is no longer available as an isolated program but it can be accessed through ctf_estimate_from_micrograph

This program takes a Xmipp image obtained from a electron microscope, computes an spectral estimation of its CTF by means of an 2D ARMA (Auto Regressive Moving Average) model and produces an output (complex) image (which is called a Fourier Xmipp image) with the filter of the ARMA model. This filter is actually an estimation of the CTF of the micrograph. The ARMA filter can then be applied to a random image of zero mean and unitary variance to get simulated noise according to the model, or more frequently the output Xmipp Fourier Image will be used to determine the CTF parameters using Adjust_CTF

Usage

$ spARMA ...

Parameters

  • -i [image] This parameter is the input image from which the ARMA model will be generated
  • -o [Fourier Xmipp image] This is the filter of the ARMA model that you have to apply to a random image to get a simulated microscope noise. The filter generated has the same dimensions of the input image
  • -N_AR [N_AR] The order of the AR part of the model in de rows direction
  • -M_AR [M_AR=N_AR] The order of the AR part of the model in de columns direction
  • -N_MA [n2] The order of the MA part of the model in de rows direction
  • -M_MA [M_MA=N_MA>] The order of the MA part of the model in de columns direction
  • -tol This is a rather technical parameter, related with the algorithm used to find the AR coeficients of the model. It is the allowed tolerance in the eigenvalues of the single value decomposition (SVD) of the linear system matrix used in the calculus. A value of 0 means that all the eigenvalues will be accepted, and it is the value to use in normal running. If problems are experienced in the solutions, the value of this parameter should be increased to reject those eigenvalues lower than the value chosen. This improves the stability of the algorithm an allows it to find the best solution compatible with a rather unstable linear system.

The last four parameters are better understood seeing this picture, where a pixel is represented by a box in the grid:

ARMA-support.jpg

The support region for the AR part of the model are the first a third quadrants (cyan ones), N_AR pixels in the rows direction and M_AR pixels in the columns direction. And the suuport region for the MA part of the model are the second and fourth quadrants (red ones). The black pixel is the actual pixel considered. Summing up, the info of the surrounding pixels is considered to generate the actual value of the black pixel. In this case, the values are : N_AR = 6, M_AR = 4, N_MA = 4, M_MA = 3.

The optimal values for N_AR , M_AR, N_MA, and M_MA were found to be 20 - 30 for N_AR , M_AR, and 15 - 25 for N_MA, M_MA, a quite wide range not substantially modified by changes in size of the images or defocus of the CTF. The recommended values are 24 for AR order and 20 for MA order, as these ones were not statistically distinguishable of the best results in all the test.

Examples and notes

Let´s see a practical example. Suppose we have an image of the microscope noise, like this one called noise.xmp:

noise.jpg

And we want to perform an spectral estimation. Well, we execute the spARMA program:

$ xmipp_spARMA  -i noise.xmp -o ARMAfilter.fft -N_AR 40 -M_AR 40 -N_MA 20 -M_MA 20

When the program ends, ARMAfilter.fft contains the filter of the ARMA model to apply to a random image to get the simulated noise. Actually, this filter is the transfer function (Fourier Transform) of the filter that the ARMA model represents.

You can do two things with this file:

I. Compare it with other forms of spectral estimation

Like periodogram averaging (this essentially is an averaging of the Fourier Transforms of the image divided by the number of pixels). FourierFilter program is used to generate the amplitude of the ARMA filter:

$ fourierfilter -i input.xmp -o output.xmp -fourier_mask ARMAfilter.fft -save_amplitude ARMAfilter_amplitude.xmp

The amplitude file can be seen like any other xmipp image, as it's not a complex one, but a normal one. When you save the amplitude, what is really saved is the logarithm10 of the squared amplitude plus 1. However, the ARMAfilter.fft contains values so that their square is the estimation of the spectral density function.

Here you can see the comparison between images, which represent the amplitude (log of absolute value squared) of the Periodogram averaging (left) and ARMA estimation (right) of the image.

periodogram.jpg AR-40-MA-20.jpg

II. Generate a simulated noise image

This is done again with FourierFilter. You have to generate a random image with zero mean and unitary variance in every pixel, of the same dimensions of the original image:

random.jpg

And then apply the ARMA filter:

$ fourierfilter -i random.xmp -o simulated-noise.xmp -fourier_mask ARMAfilter.fft

Here you have the simulated noise image, whose spectral density is similar to the original one from the microscope. This can be useful to simulate microscope noise in artificially generated images

simulated_noise.jpg

  • The values of the simulated noise image and the original one doesn´t have to be similar. What is similar is the spectral density.

  • If you want a purely autoregresive (AR) model, you must include the MA parameters, but set them to zero
    $ spARMA -i noise.xmp -o ARMAfilter.fft -N_AR 40 -M_AR 40 -N_MA 0 -M_MA 0
    The same applies if you want a purely moving-average MA filter, set the AR parameters to zero.

  • Memory usage of the program is not very high. Usually 13-20 Mb of memory space are required for normal execution. With huge images this requirement could go up to 60 Mb. The program makes intensive use of the processor. The most time consuming step is the generation of the ARMA filter in the Fourier space, which is O(n*n*(nAR+nMA )), where n is the image size, nAR is the number of AR coefficients (roughly N_AR*M_AR) and nMA is the number of MA coefficients (roughly N_AR*M_AR).

  • This program determines the ARMA model following the method depicted in:
       R. L. Kashyap, "Characterization and Estimation of Two-Dimensional
       ARMA Models, " IEEE Transactions on Information Theory, Vol. IT-30, No.
       5, pp. 736-745, September 1984.

  • 2D AR models are explained in:
       Dudgeon, D. E., Russell M. Mersereau. "Multidimensional Digital Signal Processing", Section 6.5.3, pp. 325.
       Prentice Hall Signal Processing Series. Ed. Alan Oppenheim. 
       
-- AlfredoSolano - 17 Jan 2007