Line intersection

[ErichZimmer]

Hi Ron,

Thank you for making this wonderful package as I learned a bit about PTV and image segmentation and triangulation from it. Earlier today, I was adding some more epipolar functions to my camera calibration module that will soon be part of OpenPIV Python's tomographic PIV module. One of these functions find the intersection of multiple lines using the linear least squares method. I attempted to get a functioning algorithm working, but had some technical difficulties. I noticed that you had a near identical function located here and wondered if you had any luck. Here is my function for comparison.

def _multi_line_intersect(
lines: list
) -> np.ndarray:
    P = np.zeros(3, dtype="float64")
    M = np.zeros((3,3), dtype="float64")

    n_cams = len(lines)
    
    for i in range(n_cams):
        temp = np.zeros((3,3))
        
        t = lines[i][0]
        r = lines[i][1]
                        
        temp[0, 0] = 1 - t[0] * t[0]
        temp[0, 1] = -t[0] * t[1]
        temp[0, 2] = -t[0] * t[2]
        
        temp[1, 0] = -t[1] * t[0]
        temp[1, 1] = 1 - t[1] * t[1]
        temp[1, 2] = -t[1] * t[2]
        
        temp[2, 0] = -t[2] * t[0]
        temp[2, 1] = -t[2] * t[1]
        temp[2, 2] = 1 - t[2] * t[2]
        
        M += temp
        P += np.dot(temp, r)
    
    object_points = np.linalg.lstsq(M, P, rcond=None)[0]

    return object_points  

For some reason, this function gives totally irrelevant results, even though it should be producing the approximated solution where all lines intersect in a least squares sense. See http://en.wikipedia.org/wiki/Line-line_intersection.

[ErichZimmer]

I found a solution after I told ChatGPT to rephrase the Wikipedia article in layman's terms ChatGPT sure is great for explaining topics that one is not very familiar with, especially if you have enough knowledge to correct its' mistakes. Anyhw, here is the modified function, which you can use if you like.

def _multi_line_intersect_chatgpt(
    lines: list
) -> np.ndarray:
    n_cams = len(lines)
    n_points = lines[0][1].shape[1]
    
    # cycle through each point to add to the least squares problem
    object_points = np.zeros([3, n_points], dtype="float64")
    
    # TODO: Optimize this loop as python loops would be slow if we had
    # 100,000 particles to iterate
    for particle in range(n_points):
        M = [] # A
        P = [] # b
        for cam in range(n_cams):
            t = lines[cam][0]
            r = lines[cam][1][:, particle]

            # This was very important, didn't work with out it
            r /= np.linalg.norm(r)

            # A more elegant solution to the aforementioned function
            temp = np.eye(3,3)
            temp -= np.outer(r, r) # use r instead of t, oops

            # Fill M (A) and P (b). The array shape is very important here and
            # was overlooked on my previous attempt. I still don't know what
            # went wrong, but this works so far...
            M.append(temp)
            P.append(temp.dot(t))

        M = np.asarray(M).reshape(3*n_cams, 3)
        P = np.asarray(P).reshape(3*n_cams)
        
        object_points[:, particle] = np.linalg.lstsq(M, P, rcond=None)[0]

    return object_points

[ErichZimmer]

Here is an updated function fixing my original implementation. Note that the system of equations can be solved with very simple algorithms such as LU decomposition or inverse matrix algorithms. This comes in handy since singular value decompositions with column pivoting householder reflections are hard(ish) to code in c++ without external libraries (this is mainly for a tomo-PIV module for OpenPIV-c--qt). Again, this code snip-bit is free for use and modification :)

def _multi_line_intersect(
    lines: list
) -> np.ndarray:
    # make sure that the direction vector array is 2D
    assert(len(lines[0][1].shape) == 2)
    
    n_cams = len(lines)
    n_points = lines[0][1].shape[1]
    
    A = np.zeros([3,3], dtype="float64")
    b = np.zeros(3, dtype="float64")
    
    I = np.eye(3)
    
    object_points = np.zeros([3, n_points], dtype="float64")
    
    # TODO: Optimize this loop as python loops would be slow if we had
    # 100,000 particles to iterate
    for particle in range(n_points):
        A[:, :] = 0.0
        b[:] = 0.0
        
        for cam in range(n_cams):
            t = lines[cam][0]
            r = lines[cam][1][:, particle]

            # This was very important, didn't work without it since my direction vectors are not normalized
            r /= np.linalg.norm(r)
            
            temp = I - np.outer(r, r)
        
            # Update A matrix and b vector
            A += temp
            b += temp.dot(t)
                
        object_points[:, particle] = np.linalg.lstsq(A, b, rcond=None)[0]

    return object_points