5

I have a numpy array with rectangular data on a grid and would like to fit a 2-D spline to it to reproduce the large-scale variations while getting rid of all/most of the noise. The data also has some regions marked as invalid with values of NaN.

I tried using the scipy.interpolate.RectBivariateSpline function, but the gaps screw up the results. So I tried using the LSQBivariateSpline function from the same package, hoping that when I set the weights of all NaN pixels to 0 it would simply ignore them. However, that's when I ran into the following error that I don't know how to avoid:

My code is:

# some preparation, loading data and stuff
# all my data is stored in 'data'

# Create the knots (10 knots in each direction, making 100 total
xcoord = numpy.linspace(5, data.shape[0]-5, 10)
ycoord = numpy.linspace(5, data.shape[1]-5, 10)

# Create all weights, and set them to 0 when the data is NaN
weights = numpy.ones(data.shape)
weights[numpy.isnan(data)] = 1e-15  # weights must be >0

# LSQBivariateSpline needs x and y coordinates as 1-D arrays
x, y = numpy.indices(data.shape)
spline_fit = scipy.interpolate.LSQBivariateSpline(x.ravel(), y.ravel(), data.ravel(), 
                                               xcoord, ycoord, 
                                               w=weights.ravel(),
                                               bbox=[None, None, None, None], 
                                               kx=2, ky=2)

output of the code is the following error message:

The coefficients of the spline returned have been computed as the
minimal norm least-squares solution of a (numerically) rank deficient
system (deficiency=25). If deficiency is large, the results may be
inaccurate. Deficiency may strongly depend on the value of eps.
done!
Traceback (most recent call last):
  File "./fitpg.py", line 513, in <module>
    fit_pupilghost(prebinned, (center_y, center_x), (r_inner, r_outer), dr)
  File "./fitpg.py", line 298, in fit_pupilghost
    fitout = pupil2d(radius[:,y], angle[:,y])
  File "/usr/local/lib64/python2.7/site-packages/scipy/interpolate/fitpack2.py", line 545, in __call__
    raise ValueError("Error code returned by bispev: %s" % ier)
ValueError: Error code returned by bispev: 10

The input matrix ('data') that I enter is roughly 1000 x 1000 px, that should be more than enough to constrain the splines at the 100 knots. Increasing the number of knots to 100 in each direction amek the code run waaaay slower, but except the deficiency number nothing changes. I also tried to increase and decrease the eps-value with values between 1-e30 to 0.9 (default is 1e-16

I also tried to google around for the error code, but couldn't come up with a good hit, so that didn't help either.

Any idea what might be wrong here? Or is there a workaround / better way to solve this problem?

Any help will be greatly appreciated.

Thanks

rkotulla
  • 51
  • 1
  • 3
  • Apart from BivariateSplines, see scipy.ndimage.morphology.binary_fill_holes under [finding-blank-regions-in-image](http://stackoverflow.com/questions/3310681/finding-blank-regions-in-image) . – denis Feb 08 '13 at 10:49

1 Answers1

7

The spline fitting code does not handle NaNs in any special way. Since any number that comes into touch with a NaN also becomes NaN, this means their presence poisons the whole calculation, and so you don't get any result.

What you can do, is to replace the NaN values with some (arbitrary) finite value, in addition to setting a zero weight, for instance

weights = numpy.ones(data.shape)
mask = numpy.isnan(data)
weights[mask] = 0
data[mask] = 0   # arbitrary

Since the weight is small, the value chosen shouldn't matter. You can also try setting the corresponding weights to a small value such as your 1e-15, if it for some reason doesn't like zero weights.

pv.
  • 33,875
  • 8
  • 55
  • 49