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

David Huard david.huard@gmail....
Thu Aug 2 10:07:28 CDT 2007

```oups...

y = numpy.interp(flattened_a, bins, cdf).

where

hist, bins = numpy.histogram(flattened_a, nbins, normed=True)
cdf = hist.cumsum()

2007/8/2, David Huard <david.huard@gmail.com>:
>
> Sebastian,
>
> 2007/8/2, Sebastian Haase <haase@msg.ucsf.edu>:
> >
> >
> > a1 = sa   I assume.
>
> Yes. Missed that one.
>
> What is aa ?\
>
>
>
> Sorry, that's what happen with crappy notation. Here is a clearer version.
>
>
> flattened_a = a.flatten()
>
> sorted_a = numpy.sort(flattened_a)
> cdf = sorted_a.cumsum()/sorted_a[-1]  # This is the cumulative
> distribution function [0,1]
> cdf = cdf*(max-min)+min     # Scale the CDF to whatever scale you're using
> (eg. [0,256])
>
> # Apply the scaling to the original image
> y = numpy.interp(flattened_a, sorted_a, cdf).
>
> # Reshape to the original form.
> new_a = y.reshape (a.shape)
>
>
> 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.
>
>
>
> Do you want the final image to have the same resolution ? If you want to
> optimize the equalization process,
> you could do
> y = numpy.interp(flattened_a, bins, cdf).
> hist, bins = numpy.histogram(flattened_a, nbins, normed=True)
> cdf = hist.cumsum()
>
> There were posts a while back with C implementation of histogram that is
> quite fast even for large matrices.
>
> 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 !?) ?
>
>
> Not sure what you mean, here. My 2c: get it to work then worry about
> optimization...
>
> Cheers,
>
> David
>
> Thanks again,
> > Sebastian
> >
> >
> >
> >
> >
> >
> > >
> > > 2007/7/31, Sebastian Haase <haase@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@scipy.org
> > > > http://projects.scipy.org/mailman/listinfo/scipy-user
> > > >
> > >
> > >
> > > _______________________________________________
> > > SciPy-user mailing list
> > > SciPy-user@scipy.org
> > > http://projects.scipy.org/mailman/listinfo/scipy-user
> > >
> > >
> > _______________________________________________
> > SciPy-user mailing list
> > SciPy-user@scipy.org
> > http://projects.scipy.org/mailman/listinfo/scipy-user
> >
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://projects.scipy.org/pipermail/scipy-user/attachments/20070802/977a7ff4/attachment-0001.html
```