[SciPy-User] Help optimizing an algorithm

[jkhilmer at chemistry.montana.edu]

Chris,

Can't you extend Zach's suggestion to pixel categories?  Rather than a 2^9
x 2^9 matrix, have a vector of len=3 to describe the different possible
types of pixel response curves.

Jonathan



On Wed, Jan 30, 2013 at 11:09 AM, Chris Weisiger <cweisiger at 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.
>
> -Chris
>
>
> On Wed, Jan 30, 2013 at 9:57 AM, Zachary Pincus <zachary.pincus at 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 at gmail.com> wrote:
>> > On Wed, Jan 30, 2013 at 12:29 PM, Chris Weisiger <
>> cweisiger at 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 at scipy.org
>> > > http://mail.scipy.org/mailman/listinfo/scipy-user
>> > >
>> > _______________________________________________
>> > SciPy-User mailing list
>> > SciPy-User at scipy.org
>> > http://mail.scipy.org/mailman/listinfo/scipy-user
>> >
>> > _______________________________________________
>> > SciPy-User mailing list
>> > SciPy-User at scipy.org
>> > http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20130130/cd1fc8e3/attachment.html>