[SciPy-user] histogram equalization using scipy.interpolate.splrep

[Sebastian Haase]

Hi,
I'm trying to implement an efficiant numpy / scipy based
implementation of an image analysis procedure called "histogram
equalization".
http://en.wikipedia.org/wiki/Histogram_equalization

(I don't want to use the Python Imaging Library !)

For that I calculate the intesity histogram of my 2D or 3D images
(Let's call this 1D array `h`). Then I calculate `hh=
h.astype(N.float).cumsum()`

Now I need to create a new 2D / 3D image replacing each pixel value
with the corresponding "look-up value" in `hh`.
The means essentially doing something like

`a2 = hh[a]`  -- if `a` was the original image.
Of course this syntax does not work, because
a) the array index lookup probably is not possible using 2D/3D indices
like this - right !?
and
b) my histogram is just sampling / approximating all the actually
occuring values in `a`. In other words:  the histogram is based on
creating bins, and nearby pixel values in `a` are then counted into
same bin.  Consequently there is no simple  "array index lookup" to
get the needed value out of the `h` array. Instead one needs to
interpolate.

This is were I thought the scipy.interpolate package might come in handy ;-)
In fact, if `a` was a 1D "image" I think this would work:
 >>> rep = scipy.interpolate.splrep(x,y)   # x,y being the
horizontal,count axis in my histogram
 >>> aa = scipy.interpolate.splev(a,rep)

But for a.ndim being 2 I get:
Traceback (most recent call last):
  File "<input>", line 1, in ?
  File "C:\PrWinN\scipy\interpolate\fitpack.py", line 443, in splev
    y,ier=_fitpack._spl_(x,der,t,c,k)
ValueError: object too deep for desired array

Using N.vectorize(lambda x: scipy.interpolate.splev(x,rep))
works but is slow.

I there a fast vectorized version of some interpolation already in
scipy.interpolate ?
Am I missing something ?

Thanks,
Sebastian Haase