evaluating B-Splines made with scipy.signal.cspline2d?
Hi group! in brief: I'm looking for an efficient way to evaluate B-splines generated by scipy.signal.cspline2d() at arbitrary float coordinates to interpolate image data. I am doing image processing and need a method to interpolate an image at arbitrary float coordinates. Scipy.signal kindly provides a fast and efficient routine to calculate spline coefficients, but the module seems not to adress using this spline for interpolation. I have written some code to calculate the interpolated value at arbitrary positions, but this code is cumbersome and slow - I have to calculate the value of the basis function 8 times (I suspect that takes a good while using signal.bspline(), haven't done any measurements), plus a fair bit of matrix manipulation to pick the relevant window of spline coefficients to multiply with the basis function values - my code (without any frills or checks, just to convey the basic idea) looks like this: # assume the spline coefficients are in a 2D array 'cf' # and (x,y) is the position to interpolate at def cf_matrix ( cf , x , y ) : return cf [ x - 1 : x + 3 , y - 1 : y + 3 ] def base_matrix ( x , y ) : x0 = x - floor ( x ) y0 = y - floor ( y ) rng = arange ( 1 , -3 , -1 ) xv = rng + x0 yv = rng + y0 xb = signal.bspline ( xv , 3 ) yb = signal.bspline ( yv , 3 ) xx , yy = meshgrid ( xb , yb ) return xx * yy def interpolate ( cf , x , y ) : m = cf_matrix ( cf , x , y ) b = base_matrix ( x , y ) return sum ( m * b ) obviously this is fine to use for a few hundred points or so, but for real images it's not practical. In Scipy.interpolate, there is scipy.interpolate.bisplev() to evaluate the splines generated by scipy.interpolate.bisplrep(), but these splines are of the more general form with arbitraryly spaced knot points, therefore the relevant data structures contain the knot vectors, and I assume that all sorts of performance gains that could be derived from the fact that the splines from the signal package are using equally-spaced samples are inapplicable in the general case. I'd appreciate suggestions on how to proceed or hints at what I'm missing.
Hi Kay, This doesn't answer your specific question, but look at scipy.ndimage.map_coordinates() for general-purpose spline interpolation of regularly-spaced (e.g. image) data. If you want to repeatedly interpolate the same data, you can get the spline coefficients with: scipy.ndimage.spline_filter() and pass them to map_coordinates() with the "prefilter=False" option. It is curious that while scipy.signal has cspline1d() and cspline1d_eval(), there is no cspline2d_eval() function... hopefully someone else can weight in on what's going on here. Zach On Nov 9, 2011, at 4:16 AM, Kay F. Jahnke wrote:
Hi group!
in brief: I'm looking for an efficient way to evaluate B-splines generated by scipy.signal.cspline2d() at arbitrary float coordinates to interpolate image data.
I am doing image processing and need a method to interpolate an image at arbitrary float coordinates. Scipy.signal kindly provides a fast and efficient routine to calculate spline coefficients, but the module seems not to adress using this spline for interpolation.
I have written some code to calculate the interpolated value at arbitrary positions, but this code is cumbersome and slow - I have to calculate the value of the basis function 8 times (I suspect that takes a good while using signal.bspline(), haven't done any measurements), plus a fair bit of matrix manipulation to pick the relevant window of spline coefficients to multiply with the basis function values - my code (without any frills or checks, just to convey the basic idea) looks like this:
# assume the spline coefficients are in a 2D array 'cf' # and (x,y) is the position to interpolate at
def cf_matrix ( cf , x , y ) : return cf [ x - 1 : x + 3 , y - 1 : y + 3 ]
def base_matrix ( x , y ) : x0 = x - floor ( x ) y0 = y - floor ( y ) rng = arange ( 1 , -3 , -1 ) xv = rng + x0 yv = rng + y0 xb = signal.bspline ( xv , 3 ) yb = signal.bspline ( yv , 3 ) xx , yy = meshgrid ( xb , yb ) return xx * yy
def interpolate ( cf , x , y ) : m = cf_matrix ( cf , x , y ) b = base_matrix ( x , y ) return sum ( m * b )
obviously this is fine to use for a few hundred points or so, but for real images it's not practical.
In Scipy.interpolate, there is scipy.interpolate.bisplev() to evaluate the splines generated by scipy.interpolate.bisplrep(), but these splines are of the more general form with arbitraryly spaced knot points, therefore the relevant data structures contain the knot vectors, and I assume that all sorts of performance gains that could be derived from the fact that the splines from the signal package are using equally-spaced samples are inapplicable in the general case.
I'd appreciate suggestions on how to proceed or hints at what I'm missing.
_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
thank you very much for your helpful hint. I composed a longer reply, but the web interface insists I have lines longer than 80 chars (which I haven't) and bounces my reply. Kay
I did have lines over 80 chars after all... Zachary Pincus <zachary.pincus <at> yale.edu> writes:
This doesn't answer your specific question, but look at scipy.ndimage.map_coordinates() for general-purpose spline interpolation of regularly-spaced (e.g. image) data. If you want to repeatedly interpolate the same data, you can get the spline coefficients with: scipy.ndimage.spline_filter() and pass them to map_coordinates() with the "prefilter=False" option.
Thank you very much for your helpful hint. I tried out the code you suggested and it seems to do just what I want. It seems to me that I can also process the coefficient matrix I get from signal.cspline2d() with ndimage.map_coordinates() as well as being able to have ndimage.spline_filter() generate the coefficients, though the coefficient matrices the two routines provide aren't identical. The resulting program is fast; both the calculation of the coefficient matrix for a 512X512 image and the interpolation of an equally-sized output took under a second on my system.
It is curious that while scipy.signal has cspline1d() and cspline1d_eval(), there is no cspline2d_eval() function... hopefully someone else can weight in on what's going on here.
curious indeed, but ndimage seems to solve the problem. Still signal.cspline2d_eval() should be put on the wish-list :-)
On Wed, Nov 9, 2011 at 10:19 AM, Kay F. Jahnke <_kfj@yahoo.com> wrote:
I did have lines over 80 chars after all...
Zachary Pincus <zachary.pincus <at> yale.edu> writes:
This doesn't answer your specific question, but look at scipy.ndimage.map_coordinates() for general-purpose spline interpolation of regularly-spaced (e.g. image) data. If you want to repeatedly interpolate the same data, you can get the spline coefficients with: scipy.ndimage.spline_filter() and pass them to map_coordinates() with the "prefilter=False" option.
Thank you very much for your helpful hint. I tried out the code you suggested and it seems to do just what I want. It seems to me that I can also process the coefficient matrix I get from signal.cspline2d() with ndimage.map_coordinates() as well as being able to have ndimage.spline_filter() generate the coefficients, though the coefficient matrices the two routines provide aren't identical.
I believe they both use the same algorithm, i.e., uniform b-splines and prefiltering to get the coefficients, but the boundary conditions may be different. I think the one in signal uses reflection at the ends, which is the common case. Of course, one of the routines may also be buggy ;) <snip> Chuck
participants (3)
-
Charles R Harris -
Kay F. Jahnke -
Zachary Pincus