Rebinning to polar coordinates

[arsbbr]

Hi,

I'm trying to rebin some uniform gridded data to polar coordinates. The original data (Z) is noisy but on a fine grid.

### snip

# define grid in cartesian coordinates
x = arange(-500, 500)
y = arange(-500, 500)
X, Y = meshgrid(x, y)
Z = X**2 + Y**2
Z += uniform(-5000, 5000, size=Z.shape))

# transform to polar coordinates

def cart2pol(x, y):
    theta = arctan2(y, x)
    rho = sqrt(x**2 + y**2)
    return (theta, rho)  

THETA_xy, RHO_xy = cart2pol(X, Y)

# define new polar grid
theta = arange(-pi, pi, 0.1)
rho = arange(0.1 ,500, 2.0)    

THETA, RHO = meshgrid(theta, rho)

# Tried to use griddata, but it gives of course very jumpy results, because it
# does not take the average over the pixels in a patch (dTHETA, dRHO) but an interpolation at the
# exact new point.

interpolate.griddata((THETA_xy.ravel(), RHO_xy.ravel()),
   Z.ravel(), (THETA, RHO), method='linear')
       
### snip

Is there a method that rebins the data to a new grid, possibly taking the average
of the Z-values inside a patch (dTHETA, dRHO)?

Thanks for any hints!

[Jerome Kieffer]

On Tue, 28 Aug 2012 03:33:11 -0700 (PDT)
arsbbr <[hidden email]> wrote:

>
> Hi,
>
> I'm trying to rebin some uniform gridded data to polar coordinates. The
> original data (Z) is noisy but on a fine grid.
>
> ### snip
>
> # define grid in cartesian coordinates
> x = arange(-500, 500)
> y = arange(-500, 500)
> X, Y = meshgrid(x, y)
> Z = X**2 + Y**2
> Z += uniform(-5000, 5000, size=Z.shape))
>
> # transform to polar coordinates
>
> def cart2pol(x, y):
>     theta = arctan2(y, x)
>     rho = sqrt(x**2 + y**2)
>     return (theta, rho)  
>
> THETA_xy, RHO_xy = cart2pol(X, Y)
>
> # define new polar grid
> theta = arange(-pi, pi, 0.1)
> rho = arange(0.1 ,500, 2.0)    
>
> THETA, RHO = meshgrid(theta, rho)
>
> # Tried to use griddata, but it gives of course very jumpy results, because
> it
> # does not take the average over the pixels in a patch (dTHETA, dRHO) but an
> interpolation at the
> # exact new point.
>
> interpolate.griddata((THETA_xy.ravel(), RHO_xy.ravel()),
>    Z.ravel(), (THETA, RHO), method='linear')
>        
> ### snip
>
> Is there a method that rebins the data to a new grid, possibly taking the
> average
> of the Z-values inside a patch (dTHETA, dRHO)?

Try to use 2D-weighted histogram (divided by unweighted histograms).
This unfortunately does dot split pixels over various bins so it works badly when the number of output bins is large.

w=histogram2d(THETA_xy.ravel(),RHO_xy.ravel(),(int(2*pi/0.1),500//2),weights=Z.ravel())
u=histogram2d(THETA_xy.ravel(),RHO_xy.ravel(),(int(2*pi/0.1),500//2))
imshow (w[0]/u[0])

If you want to split pixels, I am working on more sophisticated algorithm to do that in
https://github.com/kif/pyFAI  Tell me if you are interested.

Cheers,
--
Jérôme Kieffer
On-Line Data analysis / Software Group
ISDD / ESRF
tel [phone]
_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user

[arsbbr]

Dear Jerome,

thank you for your tip! Unfortunately I'm getting big gaps for my real data in the small rho region. The optimal algorithm should interpolate in the small rho region and bin in the high rho region.
I guess your split pixel solution would be best. But this a whole API with a lot of corrections and methods for real experimental data. Is it possible to use your sophisticated algorithm only for the transformation from Cartesian to a evenly stepped polar grid? Speed is not the issue with my data, so a numpy solution would be fine.

Thanks a lot!

Jerome Kieffer wrote

On Tue, 28 Aug 2012 03:33:11 -0700 (PDT)
arsbbr <[e-mail]> wrote:

>
> Hi,
>
> I'm trying to rebin some uniform gridded data to polar coordinates. The
> original data (Z) is noisy but on a fine grid.
>
> ### snip
>
> # define grid in cartesian coordinates
> x = arange(-500, 500)
> y = arange(-500, 500)
> X, Y = meshgrid(x, y)
> Z = X**2 + Y**2
> Z += uniform(-5000, 5000, size=Z.shape))
>
> # transform to polar coordinates
>
> def cart2pol(x, y):
>     theta = arctan2(y, x)
>     rho = sqrt(x**2 + y**2)
>     return (theta, rho)  
>
> THETA_xy, RHO_xy = cart2pol(X, Y)
>
> # define new polar grid
> theta = arange(-pi, pi, 0.1)
> rho = arange(0.1 ,500, 2.0)    
>
> THETA, RHO = meshgrid(theta, rho)
>
> # Tried to use griddata, but it gives of course very jumpy results, because
> it
> # does not take the average over the pixels in a patch (dTHETA, dRHO) but an
> interpolation at the
> # exact new point.
>
> interpolate.griddata((THETA_xy.ravel(), RHO_xy.ravel()),
>    Z.ravel(), (THETA, RHO), method='linear')
>        
> ### snip
>
> Is there a method that rebins the data to a new grid, possibly taking the
> average
> of the Z-values inside a patch (dTHETA, dRHO)?

Try to use 2D-weighted histogram (divided by unweighted histograms).
This unfortunately does dot split pixels over various bins so it works badly when the number of output bins is large.

w=histogram2d(THETA_xy.ravel(),RHO_xy.ravel(),(int(2*pi/0.1),500//2),weights=Z.ravel())
u=histogram2d(THETA_xy.ravel(),RHO_xy.ravel(),(int(2*pi/0.1),500//2))
imshow (w[0]/u[0])

If you want to split pixels, I am working on more sophisticated algorithm to do that in
https://github.com/kif/pyFAI  Tell me if you are interested.

Cheers,
--
Jérôme Kieffer
On-Line Data analysis / Software Group
ISDD / ESRF
tel [phone]
_______________________________________________
SciPy-User mailing list
[e-mail]
http://mail.scipy.org/mailman/listinfo/scipy-user

[Jerome Kieffer]

On Fri, 31 Aug 2012 00:36:21 -0700 (PDT)
arsbbr <[hidden email]> wrote:

> thank you for your tip! Unfortunately I'm getting big gaps for my real data
> in the small rho region. The optimal algorithm should interpolate in the
> small rho region and bin in the high rho region.
> I guess your split pixel solution would be best. But this a whole API with a
> lot of corrections and methods for real experimental data. Is it possible to
> use your sophisticated algorithm only for the transformation from Cartesian
> to a evenly stepped polar grid? Speed is not the issue with my data, so a
> numpy solution would be fine.

I wish this would be possible using only numpy but I failed.
I had to write it in a low-level (Cython) way to get decent speed because the algorithm is iterating over all input pixels.

What you are looking for is in:
https://github.com/kif/pyFAI/blob/master/src/splitBBox.pyx

Use the function histoBBox2d(weights,
                        pos0, delta_pos0,
                        pos1, delta_pos1,
                bins=(100, 36))

If you refer to the manual of numpy.histogram2d:
-> bins has the same meaning, i.e. the shape of the output image
-> pos0 and pos1 are the same a x and y; they should contain the coordinated in (rho,theta) of the center of input pixels.
-> delta_pos0 and delta_pos1 are the distance in (rho,theta) from pixel center to the pixel corner. This is of course new compared to numpy
-> weights is as in numpy the intensity of input pixels.
All array need to have the same size.

The output is:
-> your transformed image (maybe it's transposed), i.e. weighted_histogram/unweighted_histogram
-> axis in rho (or theta) 1D
-> axis in theta (or rho) 1D
-> weighted histogram (2D)
-> unweighted histogram (2D)

splitBBox should be stand alone: just use this setup.py to compile the module (needs Cython)
########################################
try:
    from setuptools import setup
except ImportError:
    from distutils.core import setup

from distutils.core import  Extension
from Cython.Distutils import build_ext

# for numpy
from numpy.distutils.misc_util import get_numpy_include_dirs


split_ext = Extension("splitBBox",
                    include_dirs=get_numpy_include_dirs(),
                    sources=['splitBBox.pyx'],

setup(name='histogram',
      version="0.3.0",
      author="Jerome Kieffer",
      author_email="[hidden email]",
      description='test for azim int',
      ext_modules=[split_ext],
      cmdclass={'build_ext': build_ext},
      )
#####################################################
Nota:
I am working on X-ray diffraction so
Your "rho" could be named "2theta", tth or q.
Your "theta" is usually named "chi" in my comments.
Different sciences; different conventions :)

I hope this helps, cheers,

--
Jérôme Kieffer
On-Line Data analysis / Software Group
ISDD / ESRF
tel [phone]
_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user