w4m w4m - 4 months ago 14
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
t2 = iradius
// 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

Answer

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

iradius = (icol_centered**2 + irow_centered**2)**0.5
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]])
>>> iradius
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
    return reduce(np.add.outer, vecs)