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

[Sebastian Haase]

On 8/2/07, David Huard <david.huard at gmail.com> wrote:
> Hi Sebastian,
>
> I know next to nothing about image transforms but it appears to me you could
> solve your problem by the following.
>
> a: image array
>
> fa = a.flatten()
>
> sa = sort(fa)
> c = a1.cumsum()/sa[-1]  # This is the cumulative distribution function [0,1]
> c = c*(max-min)+min     # Scale the CDF to whatever scale you're using (eg.
> [0,256])
>
> # Apply the scaling to the original image
> y = interp(fa, aa, c).
>
> # Reshape to the original form.
> new_a = y.reshape (a.shape)

Hi David,
thanks for the reply.
It certainly seems simple enough this way.

a1 = sa   I assume.
Is interp part of scipy.interpolate ? I don't know it.
What is aa ?

The only concern I have is that this might eat lots of memory and time
for large images. Sometimes I have 3D data with multiple 100MB in
size.
For example,  when the image's dtype is float or float32, I was hoping
that a "coarser" histogram would suffice.
Maybe, since you know "interp", you might know of another function,
that can lookup values from a function defined by a (sparse) sequence
of value pairs.
How about looking up multiple values at one (I think this is called
vectorizing !?) ?

Thanks again,
Sebastian






>
> 2007/7/31, Sebastian Haase <haase at msg.ucsf.edu>:
> >
> > 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
> > _______________________________________________
> > SciPy-user mailing list
> > SciPy-user at scipy.org
> > http://projects.scipy.org/mailman/listinfo/scipy-user
> >
>
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>
>