# [SciPy-User] Help optimizing an algorithm

josef.pktd@gmai... josef.pktd@gmai...
Wed Jan 30 12:27:37 CST 2013

```On Wed, Jan 30, 2013 at 1:09 PM, Chris Weisiger <cweisiger@msg.ucsf.edu> wrote:
> Thanks, but considering that in practice we can have values up to 2^16 from
> this camera, this seems likely to create an approximately 2^16 x 2^9 x 2^9
> array, because the camera's outputs are unsigned 16-bit integers. With two
> bytes per element, that'd mean 32.768 GB of RAM for the lookup table. Of
> course it'd be much more feasible if the input range can be constrained to
> not be the full range of data.
>
> Another thing I should note is that the nonlinearity in the camera response
> is fairly localized -- there's one region about 20 counts wide, and another
> about 500 counts wide, and the rest of the response is basically linear. So
> outside of those two regions, we can just linearly interpolate and be
> confident we're getting the right value.

somwhere in between

If the break points of the line segments are at the same location (for
the relevant pixels), then the call to np.searchsorted could still be
vectorized, and then used as index into the relevant parts of g.

However, sometimes using large intermediaries can easily be beaten by
a python loop in speed.

Josef

>
> -Chris
>
>
> On Wed, Jan 30, 2013 at 9:57 AM, Zachary Pincus <zachary.pincus@yale.edu>
> wrote:
>>
>> 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
>>
>> _______________________________________________
>> 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
>
```