Collectives™ on Stack Overflow

Find centralized, trusted content and collaborate around the technologies you use most.

Learn more about Collectives

Teams

Q&A for work

Connect and share knowledge within a single location that is structured and easy to search.

Learn more about Teams

I want to plot some points with the rbf function like here to get the density distribution of the points:

if i run the following code, it works fine:

from scipy.interpolate.rbf import Rbf  # radial basis functions
import cv2
import matplotlib.pyplot as plt
import numpy as np
# import data
x = [1, 1, 2 ,3, 2, 7, 8, 6, 6, 7, 6.5, 7.5, 9, 8, 9, 8.5]
y = [0, 2, 5, 6, 1, 2, 9, 2, 3, 3, 2.5, 2, 8, 8, 9, 8.5]
d = np.ones(len(x))
print(d)
ti = np.linspace(-1,10)
xx, yy = np.meshgrid(ti, ti)
rbf = Rbf(x, y, d, function='gaussian')
jet = cm = plt.get_cmap('jet')
zz = rbf(xx, yy)
plt.pcolor(xx, yy, zz, cmap=jet)
plt.colorbar()
# Plotting the original points.
plot3 = plt.plot(x, y, 'ko', markersize=5)  # the original points.
plt.show()

Now I want to change my input data:

# import data
input = "testProbe.jpg"
image = cv2.imread(input)   # load the image
gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)  # convert it to grayscale
# threshold the image to reveal light regions in the gray image
thresh = cv2.threshold(gray, 145, 200, cv2.THRESH_BINARY)[1]
x, y = np.where(thresh > 0)
d = np.ones(len(x))

And I get the following error message:

Traceback (most recent call last):
  File "C:/Users/.../Python/pythonprojects/03a_test_rbfWithScipy.py", line 32, in <module>
    rbf = Rbf(x, y, z, function='gaussian')
  File "C:\Users\...\Python\Anaconda2\lib\site-packages\scipy\interpolate\rbf.py", line 220, in __init__
    self.nodes = linalg.solve(self.A, self.di)
  File "C:\Users\...\Python\Anaconda2\lib\site-packages\scipy\interpolate\rbf.py", line 226, in A
    r = self._call_norm(self.xi, self.xi)
  File "C:\Users\...\Python\Anaconda2\lib\site-packages\scipy\interpolate\rbf.py", line 236, in _call_norm
    return self.norm(x1, x2)
  File "C:\Users\...\Python\Anaconda2\lib\site-packages\scipy\interpolate\rbf.py", line 118, in _euclidean_norm
    return np.sqrt(((x1 - x2)**2).sum(axis=0))
MemoryError

The error message appears really fast and when I look at the task manager, the PC is not running at full capacity. He issues this message immediately. Something can't be right then, can't it?

I tried it with other threshold values thresh = cv2.threshold(gray, 250, 255, cv2.THRESH_BINARY)[1] to get less values and it works although I get the following error message:

C:\Users\...\Python\Anaconda2\lib\site-packages\scipy\interpolate\rbf.py:220: LinAlgWarning: 
scipy.linalg.solve Ill-conditioned matrix detected. Result is not guaranteed to be accurate. 
Reciprocal condition number2.246772e-22 self.nodes = linalg.solve(self.A, self.di)

Any ideas?

If I only cut out a small part of the picture, then it works too. Maybe it makes some preliminary calculation for the memory and gives me the memory error

Have you tried increasing the threshold value such that, the rbf function uses few points? How big is the image in pixels? – jadelord Aug 3, 2018 at 5:44 @jadelord I tried it with other threshold values and it works although I get the following error message: C:\Users\...\Python\Anaconda2\lib\site-packages\scipy\interpolate\rbf.py:220: LinAlgWarning: scipy.linalg.solve Ill-conditioned matrix detected. Result is not guaranteed to be accurate. Reciprocal condition number2.246772e-22 self.nodes = linalg.solve(self.A, self.di) – Dark Aug 3, 2018 at 7:20 I got something similar (for my problem), but it was a division by zero. May be the threshold level was too extreme. I have two alternative solutions. I will answer below. – jadelord Aug 6, 2018 at 16:32

It appears to me Rbf interpolation is quite computationally intensive and leads to O(N^2) operations (correct me if I am wrong). So in order to avoid memory error you can do one of the following, instead of

zz = rbf(xx, yy)

1. Iterate with nditer

Slow, but works for small arrays:

for iz, ix, iy in np.nditer(
    [zz, xx, yy],
    flags=['external_loop', 'buffered'], 
    op_flags=['readwrite']
    iz[...] = rbf(ix, iy)

2. Use dask.array

Faster option and uses threads

import dask.array as da
n1 = xx.shape[1]
ix = da.from_array(xx, chunks=(1, n1))
iy = da.from_array(yy, chunks=(1, n1))
iz = da.map_blocks(rbf, ix, iy)
zz = iz.compute()

Hope this works for you.

P.S: As suggested on another question and scipy docs on Rbf, Rbf also allows you to replace the norm calculation function with a callable (see LowLevelCallable and the Pythran blog). It is a bit involved, but it might also solve your problem, while increasing the performance as well.

P.P.S: I managed to speedup by a factor of 50 by making this modification:

import numpy as np
def euclidean_norm_numpy(x1, x2):
    return np.linalg.norm(x1 - x2, axis=0)
rbf = Rbf(x, y, d, function='gaussian', norm=euclidean_norm_numpy)
                thank you! I tried your second option and it's 10s faster! But if I want to run my big image the memory error appears again.
– Dark
                Aug 7, 2018 at 8:41
                I am glad it works at least that much. Play with the chunks parameter until it can fit in your memory. I am also contemplating to use other interpolation schemes. Rbf is extremely slow for my application.
– jadelord
                Aug 8, 2018 at 11:31
                can you please explain what that code does? So I can better understand and play with the parameters.
– Dark
                Aug 9, 2018 at 6:39
                @Dark: See (dask.readthedocs.io/en/latest/array.html) Dask is a multithreading/multiprocessing framework. It breaks up a large array which does not fit in the memory into chunks.
– jadelord
                Aug 9, 2018 at 11:01
        

Thanks for contributing an answer to Stack Overflow!

  • Please be sure to answer the question. Provide details and share your research!

But avoid

  • Asking for help, clarification, or responding to other answers.
  • Making statements based on opinion; back them up with references or personal experience.

To learn more, see our tips on writing great answers.