How to compute residuals from OpenCV fitLine?

Question

I'm fitting a line to 3D points with OpenCV fitLine. What's the best way to calculate residuals of the resulting fit? Or, since I want residuals in addition to the fit, is there a better method than fitLine?

The following works, but there must be a better (faster) way.

# fit points
u, v, w, x, y, z = cv2.fitLine(points, cv2.DIST_L2, 0, 1, 0.01)
v = np.array(u[0], v[0], w[0])
p = np.array(x[0], y[0], z[0])

# rotate fit to z axis
k = np.cross(v, [0, 0, 1])
mag = np.linalg.norm(k)
R, _ = cv2.Rodrigues(k * np.arcsin(mag) / mag)

# rotate points and calculate distance to z-axis
rot_points = np.dot(R, (points-p).T)
err = rot_points[0]**2 + rot_points[1]**2

I'm assuming that fitLine computes the residuals err while estimating the line, so it is a waste to have to recompute them myself. Basically, knowing that I want the line and the residuals, is there a better alternative than fitLine, which only returns the line?

Answer 1

I am not aware of any direct method to get the sum of residuals from cv2.fitLine itself and I would be focusing solely on speeding up the existing code. Now, upon benchmarking with a relatively high number of points, it shows up that the most of the runtime is spent at the last two lines, where we get rot_points and err. Also, it seems that we are not exactly using the last row of rot_points for calculating err, so hopefully we can shave off some percentage of runtime by slicing into the first two rows only.

Let's dive into researching efficient ways to obtain rot_points and err.

1) rot_points = np.dot(R, (points-p).T)

This step involves broadcasting in points-p, which is reduced by matrix-multiplication later on. Now, broadcasting involves heavy memory usage, which in this case can be skipped if we split up the matrix-multiplication of R with points and p separately. Also, let's bring in the first two rows of slicing as discussed earlier. Thus, we could get the first two rows of rot_points like so -

rot_points_row2 = np.dot(R[:2], (points.T)) - np.dot(R[:2],p[:,None]) 

2) err = rot_points[0]**2 + rot_points[1]**2

This second step could be sped up with np.einsum for efficient squaring and sum-reduction, like so -

err = np.einsum('ij,ij->j',rot_points_row2,rot_points_row2)

For a relatively small number of points as 2000, the step to calculate mag : mag = np.linalg.norm(k) might become significant too in terms of runtime. So, to speed that up, one could alternatively use np.einsum again, like so -

mag = np.sqrt(np.einsum('i,i->',k,k))

---

Runtime test

Let's use a random array of 2000 points in 3D space as the input points and look at the associated runtime numbers with the original approach and the proposed ones for the last two lines.

In [44]: # Setup input points
    ...: N = 2000
    ...: points = np.random.rand(N,3)
    ...: 
    ...: u, v, w, x, y, z = cv2.fitLine(points, cv2.DIST_L2, 0, 1, 0.01)
    ...: v = np.array([u[0], v[0], w[0]])
    ...: p = np.array([x[0], y[0], z[0]])
    ...: 
    ...: # rotate fit to z axis
    ...: k = np.cross(v, [0, 0, 1])
    ...: mag = np.linalg.norm(k)
    ...: R, _ = cv2.Rodrigues(k * np.arcsin(mag) / mag)
    ...: 
    ...: # rotate points and calculate distance to z-axis
    ...: rot_points = np.dot(R, (points-p).T)
    ...: err = rot_points[0]**2 + rot_points[1]**2
    ...: 

Let's run our proposed methods and verify their outputs against the original outputs -

In [45]: rot_points_row2 = np.dot(R[:2], (points.T)) - np.dot(R[:2],p[:,None]) 
    ...: err2 = np.einsum('ij,ij->j',rot_points_row2,rot_points_row2)
    ...: 

In [46]: np.allclose(rot_points[:2],rot_points_row2)
Out[46]: True

In [47]: np.allclose(err,err2)
Out[47]: True

Finally and most importantly, let's time these sections of code -

In [48]: %timeit np.dot(R, (points-p).T) # Original code
10000 loops, best of 3: 79.5 µs per loop

In [49]: %timeit np.dot(R[:2], (points.T)) - np.dot(R[:2],p[:,None]) # Proposed
10000 loops, best of 3: 49.7 µs per loop

In [50]: %timeit rot_points[0]**2 + rot_points[1]**2  # Original code
100000 loops, best of 3: 12.6 µs per loop

In [51]: %timeit np.einsum('ij,ij->j',rot_points_row2,rot_points_row2) # Proposed
100000 loops, best of 3: 11.7 µs per loop

When the number of points is increased to a huge number, the runtimes look more promising. With N = 5000000 points, we get -

In [59]: %timeit np.dot(R, (points-p).T) # Original code
1 loops, best of 3: 410 ms per loop

In [60]: %timeit np.dot(R[:2], (points.T)) - np.dot(R[:2],p[:,None]) # Proposed
1 loops, best of 3: 254 ms per loop

In [61]: %timeit rot_points[0]**2 + rot_points[1]**2  # Original code
10 loops, best of 3: 144 ms per loop

In [62]: %timeit np.einsum('ij,ij->j',rot_points_row2,rot_points_row2) # Proposed
10 loops, best of 3: 77.5 ms per loop

Answer 2

Just wanted to mention that I compared PCA with fitLine for speed. fitLine was faster than PCA for all the cases where the number of points was greater than ~1000.

PCA directly gives you the eigenvectors, so you don't need the Rodrigues (but this time should be negligible). So, may be, the key to optimization lies in the rest of the code, unless there's a faster way to fit the model, other than fitLine and PCA.

I don't remember the PCA related math very well, so I could be wrong in the paragraph that follows.

Eigenvalues give us the variance along each dimension of the new eigenspace. If we consider a simple 2D case, I think you can take the smaller eigenvalue and multiply it with the number of points N (or is it N-1 ?) in the dataset to obtain the squared-residual-sum. Similarly we can extend it to the 3D case. As PCA gives us the eigenvalues, it takes simple scalar multiplications and additions to get the squared-residual-sum.

I'm adding the code (c++) for your reference.

RNG rng;
float a = -0.1, b = 0.1;
int rows = 3, cols = 1024*2;

Mat data = Mat::zeros(rows, cols, CV_32F);

for (int i = 0; i < cols; i++)
{
    Vec3f& v = data.at<Vec3f>(i);
    v = Vec3f(i+rng.uniform(a, b), i+rng.uniform(a, b), i+rng.uniform(a, b));
}

Mat datat = data.t();

Vec6f line;
fitLine(datat, line, CV_DIST_L2, 0, 1, 0.01);

PCA pca(datat, Mat(), CV_PCA_DATA_AS_ROW);

cout << "fitLine:\n" << line << endl;
cout << "\nPCA:\n" << pca.eigenvalues << endl;
cout << pca.eigenvectors << endl;