# [SciPy-User] Help optimizing an algorithm

Zachary Pincus zachary.pincus@yale....
Wed Jan 30 11:57:59 CST 2013

```If g is different for each pixel, the look-up table approach will probably still work. You'll need a 3D look-up table mapping the function for each pixel:

2x2 case:
table = numpy.array([[[0,1,200], [0,1,2]], [[1,2,3], [0,100,200]]])
input = numpy.array([[0,2],[1,1]])

I can't figure out exactly how to do this right now and my flight's boarding, but perhaps some fancy-indexing wizards can help. Also, scipy.ndimage.interpolate could be used in various contexts to do look-up/function interpolation in this context, all in parallel.

Zach

On Jan 30, 2013, at 10:47 AM, Chris Weisiger wrote:

> Right, I should have clarified that g is different for each pixel. It looks like scipy.interpolate.interp1d ought to do exactly what I want, though I'll have to handle the bounds conditions (where the input data is outside the range of the interpolation function that interp1d generates) myself. Thanks for the help!
>
> -Chris
>
>
> On Wed, Jan 30, 2013 at 9:38 AM, <josef.pktd@gmail.com> wrote:
> On Wed, Jan 30, 2013 at 12:29 PM, Chris Weisiger <cweisiger@msg.ucsf.edu> wrote:
> > We have a camera at our lab that has a nonlinear (but monotonic) response to
> > light. I'm attempting to linearize the data output by the camera. I'm doing
> > this by sampling the response curve of the camera, generating a linear fit
> > of the sample, and mapping new data to the linear fit by way of the sample.
> > In other words, we have the following functions:
> >
> > f(x): the response curve of the camera (maps photon intensity to reported
> > counts by the camera)
> > g(x): an approximation of f(x), composed of line segments
> > h(x): a linear fit of g(x)
> >
> > We get a new pixel value Y in -- this is counts reported by the camera. We
> > invert g() to get the approximate photon intensity for that many counts. And
> > then we plug that photon intensity into the linear fit.
> >
> > Right now I believe I have a working algorithm, but it's very slow (which in
> > turn makes testing for validity slow), largely because inverting g()
> > involves iterating over each datapoint in the approximation to find the two
> > that bracket Y so that I can linearly interpolate between them. Having to
> > iterate over every pixel in the image in Python isn't doing me any favors
> > either; we typically deal with 528x512 images so that's 270k iterations per
> > image.
> >
> > If anyone has any suggestions for optimizations I could make, I'd love to
> > hear them. My current algorithm can be seen here:
> > http://pastebin.com/mwaxWHGy
>
>
> np.searchsorted or scipy.interp1d
>
> If g is the same for all pixels, then there is no loop necessary and
> can be done fully vectorized
>
> Josef
>
> >
> > -Chris
> >
> > _______________________________________________
> > SciPy-User mailing list
> > SciPy-User@scipy.org
> > http://mail.scipy.org/mailman/listinfo/scipy-user
> >
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user

```