[SciPy-user] Numarray image processing

[Perry Greenfield]

On May 10, 2005, at 11:29 AM, Jeremy Sanders wrote:

> Hi -
>
> I'm trying to think of an efficient way to do the following in 
> numarray. I have an image, and I want to caclulate the mean of the 
> pixels at a particular radius in the image (from a fixed point), and 
> subtract this mean away from those pixels at that radius.
>
> Can anyone provide insight into how this can be done without looping 
> over the pixels?
>
> I've managed to generate an image containing the radius of each pixel, 
> but this isn't much use.
>
> Thanks
>
> Jeremy

How you want to do this depends on whether you want to do this for only 
one radius or many. As it turns out, I just did an example for a 
tutorial that computes the azimuthally-averaged radial. So here it is 
in the off-chance that it is useful. Note that it is not as fast as C 
code would be primarily because of the sort needed. It's a bit tricky.

 >>> # calculate average radial profile
 >>> y, x = indices((512,512)) # first determine radii of all pixels
 >>> r = sqrt((x-257.)**2+(y-258)**2
 >>> ind = argsort(r.flat) # get sorted indices
                           # needed to arange image values too
 >>> sr = r.flat[ind] # sorted radii
 >>> sim = im.flat[ind] # image values sorted by radii
 >>> ri = sr.astype(Int16) # integer part of radii (bin size = 1)

 >>> deltar = ri[1:] - ri[-1] # assume all radii represented
                              # (more work if not)
 >>> rind = where(deltar)[0] # location of changed radius
 >>> nr = rind[1:] - rind[:-1] # number in radius bin
 >>> sim = sim*1. # turn into double to avoid integer overflow
                  # (and single precision rounding)
 >>> csim = cumsum(sim) # cumulative sum to figure out sums for each 
radii bin
 >>> tbin = csim[rind[1:]] - csim[rind[:-1]] # sum for image values in 
radius bins
 >>> radialprofile = tbin/nr # the answer (note missing 0 radius value)


This is only part of what you need. The profile here can be used to 
generate a profile image. How you do that depends on how much 
interpolation you want. If none, this would approximate by taking the 
nearest radius value

 >>> radiusimage = radialprofile[(r+.5).astype(Int16)]

If you only need it for one radius then all you need do is construct a 
mask from the radius image and total all values in the masked image 
divided by the total of the mask.

Is this what you were looking for?

Perry Greenfield