w4m - 1 year ago 74
Python Question

# Is there any python (numpy) synonym to i-variables (e.g., irow, iradius, itheta, etc.) in DM scripting?

I work with electron microscopy image processing mostly with digital microscopy (DM) scripting, and recently start to learn Python because of its wider versatility, rich open libraries, and cross-platform ability.

Does anybody know if there are any similar tools in Python (numpy) to index 2D (image), or 3D (spectrum image) arrays similar to DM's i-variables?

The i-variables are briefly introduced on page 11 of this tutorial about DM-scripting:
http://portal.tugraz.at/portal/page/portal/Files/felmi/images/DM-Script/DM-basic-scripting_bs.pdf

They are easy way to index any image-like 2D or 3D object, which is very convenient for image processing, e.g., generate mask functions

For example, the following DM-script

``````image t1 := RealImage ("test1", 4, 5, 5)
image t2 := RealImage ("test2", 4, 5, 5)
image t3 := RealImage ("test3", 4, 5, 5)
t1 = irow
// the value in each pixel equals to the row index
// the value in each pixel equals to the radius
// (i.e., distance to the center pixel)
t3 = itheta
// the value in each pixel quals to the angle (radian)
// to the center pixel (i.e., angle in polar representation)
t1.showimage(); t2.showimage(); t3.showimage()
``````

result in the following images (expressed here in spreadsheet, or say matrix form):

``````t1 =
0   0   0   0   0
1   1   1   1   1
2   2   2   2   2
3   3   3   3   3
4   4   4   4   4

t2=
3.5355339   2.9154758   2.5495098   2.5495098   2.9154758
2.9154758   2.1213202   1.5811388   1.5811388   2.1213202
2.5495098   1.5811388   0.70710677  0.70710677  1.5811388
2.5495098   1.5811388   0.70710677  0.70710677  1.5811388
2.9154758   2.1213202   1.5811388   1.5811388   2.1213202

t3=
-2.3561945  -2.1112158  -1.7681919  -1.3734008  -1.0303768
-2.6011732  -2.3561945  -1.8925469  -1.2490457  -0.78539819
-2.9441972  -2.8198421  -2.3561945  -0.78539819 -0.32175055
2.9441972   2.8198421   2.3561945   0.78539819  0.32175055
2.6011732   2.3561945   1.8925469   1.2490457   0.78539819
``````

The equivalent way to do this in NumPy is to use the `numpy.indices` function.

So to do the same in NumPy as you did in DM (remember coordinates are always (y,x) in numpy, and that irow, say, is index of y-coordinates):

``````from __future__ import division
import numpy as np

test1 = np.random.random((5,5))
irow, icol = np.indices(test1.shape)

# to use iradius and itheta we need to get icol, irow centered
irow_centered = irow - test1.shape[0] / 2.0
icol_centered = icol - test1.shape[1] / 2.0

itheta = np.arctan2(irow_centered, icol_centered)
``````

Then

``````>>> irow
array([[0, 0, 0, 0, 0],
[1, 1, 1, 1, 1],
[2, 2, 2, 2, 2],
[3, 3, 3, 3, 3],
[4, 4, 4, 4, 4]])
array([[ 3.53553391,  2.91547595,  2.54950976,  2.54950976,  2.91547595],
[ 2.91547595,  2.12132034,  1.58113883,  1.58113883,  2.12132034],
[ 2.54950976,  1.58113883,  0.70710678,  0.70710678,  1.58113883],
[ 2.54950976,  1.58113883,  0.70710678,  0.70710678,  1.58113883],
[ 2.91547595,  2.12132034,  1.58113883,  1.58113883,  2.12132034]])
>>> itheta
array([[-2.35619449, -2.11121583, -1.76819189, -1.37340077, -1.03037683],
[-2.60117315, -2.35619449, -1.89254688, -1.24904577, -0.78539816],
[-2.94419709, -2.8198421 , -2.35619449, -0.78539816, -0.32175055],
[ 2.94419709,  2.8198421 ,  2.35619449,  0.78539816,  0.32175055],
[ 2.60117315,  2.35619449,  1.89254688,  1.24904577,  0.78539816]])
``````

It's also possible to do this using the `mgrid` and `ogrid` functions. `mgrid` will return the coordinates in a fully populated array (the same shape as the source image, identical to numpy.indices) while `ogrid` returns row and column vectors of the right shape that often get automatically broadcast correctly.

If you're using this to create masks for Fourier-transformed images there's also the np.fft.fftfreq function which will give you the frequencies of the pixels in the Fourier transform. I use the following to get the frequency squared at each pixel for a given image shape:

``````def get_freq_squared(shape, packed_complex=True):
"""Return the frequency squared for an n-d array.
The returned image will match shape, if packed_complex is true the last
dimension will be treated as 0,f,...,nf rather than 0,f,..., nf,...,-f
"""
vecs = [np.fft.fftfreq(s, 1.0 / s) ** 2 for s in shape]

if packed_complex:
s = shape[-1]
olds = (s-1)*2
vecs[-1] = rfftfreq(olds, 1.0/olds)**2