How to project pixels on eigenvectors in OpenCV?

Viewed 251

given a contour, I can extract mean, and eigenvectors by performing PCA. Then I want to project all pixels inside a contour on an eigenvector. Below are my code and my images

my input images my input images

Read image, extract contours, and draw the first component

import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt

src = cv.imread(cv.samples.findFile('/Users/bryan/Desktop/lung.png'))
gray = cv.cvtColor(src, cv.COLOR_BGR2GRAY)

def drawAxis(img, p_, q_, colour, scale):
    p = list(p_)
    q = list(q_)
    ## [visualization1]
    angle = atan2(p[1] - q[1], p[0] - q[0])  # angle in radians
    hypotenuse = sqrt((p[1] - q[1]) * (p[1] - q[1]) + (p[0] - q[0]) * (p[0] - q[0]))

    # Here we lengthen the arrow by a factor of scale
    q[0] = p[0] - scale * hypotenuse * cos(angle)
    q[1] = p[1] - scale * hypotenuse * sin(angle)
    cv.line(img, (int(p[0]), int(p[1])), (int(q[0]), int(q[1])), colour, 1, cv.LINE_AA)

    # create the arrow hooks
    p[0] = q[0] + 9 * cos(angle + pi / 4)
    p[1] = q[1] + 9 * sin(angle + pi / 4)
    cv.line(img, (int(p[0]), int(p[1])), (int(q[0]), int(q[1])), colour, 1, cv.LINE_AA)

    p[0] = q[0] + 9 * cos(angle - pi / 4)
    p[1] = q[1] + 9 * sin(angle - pi / 4)
    cv.line(img, (int(p[0]), int(p[1])), (int(q[0]), int(q[1])), colour, 1, cv.LINE_AA)
    
def getOrientation(pts, img, scale_factor=25):
    ## [pca]
    # Construct a buffer used by the pca analysis
    sz = len(pts)
    data_pts = np.empty((sz, 2), dtype=np.float64)
    for i in range(data_pts.shape[0]):
        data_pts[i, 0] = pts[i, 0, 0]
        data_pts[i, 1] = pts[i, 0, 1]

    # Perform PCA analysis
    mean = np.empty(0)
    mean, eigenvectors, eigenvalues = cv.PCACompute2(data_pts, mean)

    # Store the center of the object
    cntr = (int(mean[0, 0]), int(mean[0, 1]))
    ## [pca]

    ## [visualization]
    # Draw the principal components
    cv.circle(img, cntr, 3, (255, 0, 255), 2)

    p1 = (
        cntr[0] + scale_factor * eigenvectors[0, 0],
        cntr[1] + scale_factor * eigenvectors[0, 1])
    p2 = (
        cntr[0] - scale_factor * eigenvectors[1, 0],
        cntr[1] - scale_factor * eigenvectors[1, 1])
    drawAxis(img, cntr, p1, (0, 255, 0), 4)
    
    ## [visualization]

    # doing projections along eigenvectors
    dim1_ = []
    for _ in data_pts:
        p = make_vector_projection(_, np.array(p1))
        dim1_.append(p.astype(int))
    dim1 = np.array(dim1_)

    dim2_ = []
    for _ in data_pts:
        p = make_vector_projection(_, np.array(p2))
        dim2_.append(p.astype(int))
    dim2 = np.array(dim2_)

    return mean, eigenvectors, eigenvalues, p1, p2, dim1, dim2

for i, c in enumerate(contours):
    mean, evecs, evalues, p1, p2, dim1, dim2 = getOrientation(c, src)

    
plt.figure(figsize=(10, 10))
plt.axis('equal')
plt.gca().invert_yaxis()
plt.imshow(src)

enter image description here

Look as I expected but I want to double-check, I compute angles between a random vector in the first dimension and the first eigenvector. I defined

def unit_vector(vector):
    """ Returns the unit vector of the vector.  """
    return vector / np.linalg.norm(vector)


def angle_between(v1, v2):
    """ Returns the angle in radians between vectors 'v1' and 'v2'::

            >>> angle_between((1, 0, 0), (0, 1, 0))
            1.5707963267948966
            >>> angle_between((1, 0, 0), (1, 0, 0))
            0.0
            >>> angle_between((1, 0, 0), (-1, 0, 0))
            3.141592653589793
    """
    v1_u = unit_vector(v1)
    v2_u = unit_vector(v2)
    return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))



for i, c in enumerate(contours):
    mean, evecs, evalues, p1, p2, dim1, dim2 = getOrientation(c, src)
    draw_point = dim1[45].astype(int)  # extract the first dimension
    print(draw_point, evecs[0], angle_between(draw_point, p1))
    print('#'* 10)

I got output like below, angle_between is in radian, close to zero means they are paralleled

[ 97 148] [ 0.14189901 -0.98988114] 0.002321780300502494
##########
[332 163] [-0.22199134 -0.97504864] 0.0006249775357550807
##########

My problem occurs when I want to plot projected points on eigenvectors. Points I plotted don't lie on green line (my eigenvector). My code is

point_colors = [
    (0, 0, 255), # blue
    (0, 255, 0) # green
]

for i, c in enumerate(contours):
    mean, evecs, evalues, p1, p2, dim1, dim2 = getOrientation(c, src)
    draw_point = dim1[45].astype(int)
    print(draw_point, evecs[0], angle_between(draw_point, p1))
    cv.circle(src, (draw_point[1], draw_point[0]), 7, point_colors[i], 2) # plot point
    print('#'* 10)

and output output

My question is, I want to plot all projected pixels on the eigenvector lines, but my calculation doesn't seem correct as points don't lie on the eigenvector lines. Could you help?

0 Answers
Related