[SciPy-User] Scaling up root-finding

[Chris Weisiger]

...yes, sorry, I should have thought that through. For some reason I
thought the quadratic formula wasn't amenable to vectorization, but on
reflection that's not the case. The +- portion is a bit tricky, but given
that we know the resulting value is positive and one of its roots is where
x < 0, I think it's safe to say that we should always use +, not -.

Thanks for the gut-check. Sometimes you just need someone to point out the
obvious!

I'm still interested in knowing if there's some efficient way to handle
arbitrary polynomials, though.

-Chris


On Wed, Sep 4, 2013 at 1:45 PM, Nathaniel Smith <njs at pobox.com> wrote:

> If you know the polynomial coefficients a, b, c and the reported photons,
> then why not just use the textbook quadratic formula to get a (vectorised!)
> closed form for the incident photons?
>
> This only works for a few polynomial degrees of course, and there's a
> small subtlety in selecting the correct root, but it should be more
> efficient than your current approach...
>
> -n
> On 4 Sep 2013 21:39, "Chris Weisiger" <cweisiger at msg.ucsf.edu> wrote:
>
>> We have a somewhat elderly microscopy camera whose response curve is
>> best-characterized as a degree-2 polynomial (starts out steep, then
>> flattens out; we hit the camera's maximum bit depth before the function
>> stops being monotonic though). That is:
>>
>> reported photons = a + b * (incident photons) + c * (incident photons)^2
>>
>> Of course, each pixel is slightly different, giving 512*512 different
>> polynomials. I want to linearize the camera's response in post-production,
>> by inverting the (monotonic) function so that it maps reported photons to
>> incident photons.
>>
>> The obvious method would be to use something from scipy.optimize's
>> root-finding functions. I'm not remotely a numerical analysis expert, so I
>> have no feel for which function would be best for this problem.
>> Performance-wise, a quick test indicates that it would take about 3.5s on
>> my laptop to process a single frame by manually iterating over every pixel
>> in the frame and calling scipy.optimize.newton(that pixel's response
>> function, that pixel's starting value). This isn't great, but it's well
>> within the realm of feasibility.
>>
>> Is there some more efficient method to do this? In the ideal world we'd
>> be able to handle arbitrary-degree polynomials, but I'd be willing to
>> accept being restricted to degree-2. Of course, degree-1 is trivial to
>> implement.
>>
>> Thanks for any advice you care to share.
>>
>> -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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20130904/8cc3b9fa/attachment.html>