w4m - 1 year ago 62

Python Question

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 Source

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