[SciPy-User] Help optimizing an algorithm
Wed Jan 30 11:48:16 CST 2013
> 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
Don't have the time to fully spell this out now (in an airport), but the general gist of what you should do is:
(1) Make a look-up table mapping input pixel values (as indices) to the desired linearized values. So in a simple example with three different possible incoming pixel values 0, 1, and 2, which should be mapped to values 0, 10, and 200, you would have an array like so:
table = numpy.array([0, 10, 200]).
To make the input/output relationship between the pixel values and index positions clear, you then have:
table = 0
table = 10
table = 200
which is exactly the functional relationship you want.
(2) Use fancy-indexing to replace input pixel values with look-up table values. (Read up on the numpy fancy-indexing tutorials that google can point you to.)
table = numpy.array([0, 10, 200])
input = numpy.array([[0, 1, 2, 0], [0, 0, 0, 1], [1, 1, 1, 1]])
output = table[input]
array([[ 0, 10, 200, 0],
[ 0, 0, 0, 10],
[ 10, 10, 10, 10]])
Basically, what this is saying is to treat "input" as an array of indices into the table array, and generate an output array that contains the value of the table at each index.
It's plenty fast, too (using an example 12-bit camera that produces 2048x2048-pixel images):
table = numpy.empty(dtype=numpy.float, shape=(4096,))
input = numpy.empty(dtype=numpy.uint16, shape=(2048,2048))
timeit table[input] # note timeit is from ipython, which you should be using if you aren't.
10 loops, best of 3: 43.3 ms per loop
Note that numpy.take does basically the same thing, but there are less special cases, so it's faster:
In : timeit numpy.take(table, input)
10 loops, best of 3: 27.4 ms per loop
Obviously this only works if your camera produces integer pixel values that can be used as indices into an array. But I don't know of any cameras that don't do this.
More information about the SciPy-User