0
\$\begingroup\$

I have this code that generates a random line and a random circle.

import numpy as np
import matplotlib.pyplot as plt
def random_2D_position(c_bx, c_by, r, region_side):
    while True:
        v_pt = region_side * np.random.random((2, ))
        if np.linalg.norm(v_pt - np.array([c_bx, c_by])) > r:
            return v_pt
    return
def plot_region(c_bx, c_by, r, v_l_1, v_l_2, region_side):
    plt.plot(v_l_1[0], v_l_1[1], 's')
    plt.plot(v_l_2[0], v_l_2[1], 'o')
    circle = plt.Circle((c_bx, c_by), r, fill=False)
    plt.gca().add_artist(circle)
    plt.plot([v_l_1[0], v_l_2[0]],[v_l_1[1], v_l_2[1]], '-')
    plt.xlim((0, region_side))
    plt.ylim((0, region_side))
    return
region_side = 500  
c_bx = region_side / 3 + region_side / 3 * np.random.random()
c_by = region_side / 3 + region_side / 3 * np.random.random()
r = region_side / 4 + region_side / 30 * np.random.random()
v_l_1 = random_2D_position(c_bx, c_by, r,region_side)
v_l_2 = random_2D_position(c_bx, c_by, r, region_side)
plot_region(c_bx, c_by, r, v_l_1, v_l_2, region_side)
plt.show()

I need to find the portion of the line that is inside the circle if the line intersects the circle. I tried to solve it as follows:

First, find the line equation :

def Find_line_eq(points):
    #https://stackoverflow.com/a/21566184/6214597
    x_coords, y_coords = zip(*points)
    A = np.vstack([x_coords,np.ones(len(x_coords))]).T
    m, c = np.linalg.lstsq(A, y_coords)[0]
    return m,c

Then, I used the following code to get the roots of the quadratic formula:

m,c_l=Find_line_eq(points=[v_l_1,v_l_2])
a=1+m**2
b=2*c_l*m
c=c_l**2-r**2
coeff=[a, b, c]
x1,x2= np.roots(coeff)

Questions:

1- Given the two points of the line, the circle center, and radius (v_l_1, v_l_2, c_bx, c_by,r), how to obtain the portion of the line that is inside the circle? Most of the time x1,x2 are not real, Numpy returns a complex number, even if the line is intersecting the circle?

2- How can we generalize this code to obtain these intersection points for 3D circle and line (the two points of the line are not at the same level )?

\$\endgroup\$
5
  • \$\begingroup\$ Are you sure this code is working as intended? My editor gives me a bunch of indentation errors and missing brackets. \$\endgroup\$ Commented Aug 24, 2021 at 18:13
  • \$\begingroup\$ @N3buchadnezzar it works, I copied it from my Jupyter Notebook . The first function ‘’return" was not in place. I have edited it now \$\endgroup\$
    – Paulo
    Commented Aug 24, 2021 at 18:24
  • \$\begingroup\$ @N3buchadnezzar You can check now \$\endgroup\$
    – Paulo
    Commented Aug 24, 2021 at 18:40
  • 1
    \$\begingroup\$ (@N3buchadnezzar: Did you intend working as indented?;) \$\endgroup\$
    – greybeard
    Commented Aug 25, 2021 at 7:13
  • \$\begingroup\$ With code presented that produces the result/effect required, CodeReview@SE is ready to provide insights about and opinions on the source code. Questions about how to achieve a result are off topic, with those asking how to achieve something additional close to the border between on- and off-topic with Code Review. \$\endgroup\$
    – greybeard
    Commented Aug 25, 2021 at 7:31

1 Answer 1

1
\$\begingroup\$

Are you asking for feedback on your code, or do you just want someone to fix your two bullet points? Asking for extensions to the code outside of (is there a better algorithm/can this be done in a better way), is outside the scope of this site and is better suited for stack overflow.

Lets just quickly address some missing points from your code

Your code is very hard to read

  • Avoid using single letter variable names
  • Structure your code better using more functions and a if name == "__main__" guard.
  • Add comments and docstrings to explain what the heck is going on.
  • Add typing hints to hint at what the variables mean (but first improve their names).
  • Follow PEP 8 and run a linter over your code. Even better get some help from a language server like pyright.

Lack of detail

def random_2D_position(c_bx, c_by, r, region_side):
    while True:
        v_pt = region_side * np.random.random((2, ))
        if np.linalg.norm(v_pt - np.array([c_bx, c_by])) > r:
            return v_pt
    return
def plot_region(c_bx, c_by, r, v_l_1, v_l_2, region_side):
    plt.plot(v_l_1[0], v_l_1[1], 's')

The second return statement in random_2D_position will never be reached. Secondly you should always have two spaces between functions in the outer scope. This line has trailing whitespace

region_side = 500  

For instance the lines

m,c_l=Find_line_eq(points=[v_l_1,v_l_2])
a=1+m**2
b=2*c_l*m
c=c_l**2-r**2
coeff=[a, b, c]
x1,x2= np.roots(coeff)

Is better formated as

m, c_l = Find_line_eq(points=[v_l_1, v_l_2])
a = 1 + m ** 2
b = 2 * c_l * m
c = c_l ** 2 - r ** 2
coeff = [a, b, c]
x1, x2 = np.roots(coeff)

Notice how you have inconsistent spacing around the = symbols for instance. All of this is solved by running a proper formater over your code.

Magic numbers

Your code is littered with magic numbers what does any of this mean?

c_bx = region_side / 3 + region_side / 3 * np.random.random()
c_by = region_side / 3 + region_side / 3 * np.random.random()
r = region_side / 4 + region_side / 30 * np.random.random()

These should ideally be

  1. Extracted to their own function
  2. The magic numbers should be defined as global constants (the standard way of doing this in Python is UPPERCASE)
  3. Some docstrings / typing hints would go a long way explaining what is going on.
  4. As a last resort add some comments.

Nitpicking

There is a new way of getting random numbers in numpy, your method is supported because of legacy reasons, but the newer syntax should be prefered. See https://numpy.org/doc/stable/reference/random/generator.html?highlight=default_rng#numpy.random.default_rng

If you copy code from Stack Overflow you should link to it or at the very least explain it. Now this is just a hunch, but I would never naturally write this

np.linalg.norm(v_pt - np.array([c_bx, c_by]))

for calculating the euclidean distance between two points. Yet it is the first result that shows up when googling "How do I calculate the euclidian distance in numpy", it might be a coincidence but at the very least extract this piece of code into it's own code and source it

def euclidian_distance(A, B):
    """Calculates the euclidian distance between two points

    See https://stackoverflow.com/a/1401828/1048781 for details
    """
    return np.linalg.norm(A - B)

Using add_artist but it's not recommended see here for a more modern way.

Revised code

The following code works for any higher dimension and calculates the orthogonal projection to figure out if the line and the circle intersects. I will leave it to you to figure out how to compute the intersections. Adding typing hints and fleshing out the docstrings is also left as homework =)

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng()

SIDE = 500
DIMENSIONS = 2
EPSILON = 1


def get_random_point(start, stop=None):
    if stop is None:
        start, stop = 0, start
    return (stop - start) * rng.random((DIMENSIONS,)) + start


def get_random_circle():
    pointA = get_random_point(SIDE)
    max_radius = distance_2_closest_side(pointA)
    radius = max_radius * rng.random()
    return (pointA, radius)


def distance_2_closest_side(A):
    return min(min(a, SIDE - a) for a in A)


def euclidian_distance(A, B):
    """Calculates the euclidian distance between two points

    See https://stackoverflow.com/a/1401828/1048781 for details
    """
    return np.linalg.norm(A - B)


def get_line_outside_circle(circle, radius):
    A = B = circle
    while euclidian_distance(A, circle) < radius:
        A = get_random_point(SIDE)
    in_circle = euclidian_distance(B, circle) < radius
    AB_too_close = euclidian_distance(A, B) < EPSILON
    while in_circle and AB_too_close:
        B = get_random_point(SIDE)
        in_circle = euclidian_distance(B, circle) < radius
    line = [A, B]
    return line


def projection(b, a):
    """orthogonal projection of a onto a straight line parallel to b

    See https://en.wikipedia.org/wiki/Vector_projection for further details
    """
    return ((b @ a) / (b @ b)) * b


def orthogonal_projection(A, B, circle):
    """Calculates the orthogonal projection from circle onto the line AB

    See https://stackoverflow.com/a/9368901/1048781 for further details
    """
    AB = B - A
    AC = circle - A
    # projection = AC - AB * ((AB @ AC) / (AB @ AB))
    proj = AC - projection(AB, AC)

    orthogonal_circle_line = circle - proj
    return orthogonal_circle_line


def plot_results(A, B, circle, radius):

    # now make a circle with no fill, which is good for hi-lighting key results
    circle1 = plt.Circle(circle, radius, color="r", fill=False)

    ax = plt.gca()
    ax.cla()  # clear things for fresh plot

    # change default range so that new circles will work
    ax.set_xlim((0, SIDE))
    ax.set_ylim((0, SIDE))
    # some data
    ax.scatter([A[0], B[0]], [A[1], B[1]], color="blue")
    ax.axline(A, B, color="blue")
    ax.scatter([circle[0]], [circle[1]], color="red")
    # key data point that we are encircling
    ax.plot((5), (5), "o", color="y")

    ax.add_patch(circle1)
    plt.show()


def main():
    circle, radius = get_random_circle()
    A, B = get_line_outside_circle(circle, radius)

    proj = orthogonal_projection(A, B, circle)
    if euclidian_distance(circle, proj) < radius:
        print("We have an intersection!")

    plot_results(A, B, circle, radius)


if __name__ == "__main__":

    main()
\$\endgroup\$
2
  • \$\begingroup\$ First, thank you very much. We assume that the line could be inside or outside the circle, in case the line intersects with the circle, what is the length of the line that is inside the circle according to the code above? Does this code return the intersection points (x1,x2),(y1,y2)? \$\endgroup\$
    – Paulo
    Commented Aug 25, 2021 at 1:31
  • \$\begingroup\$ @Paulo "Asking for extensions to the code outside of (is there a better algorithm/can this be done in a better way), is outside the scope of this site and is better suited for stack overflow." Your method works for finding the intersection points fine and can be extended to arbitrary dimensions, see for instance math.stackexchange.com/a/2536095/18908. The distance can be computed from the inner product (or our function euclidian_distance) see math.stackexchange.com/a/2981910/18908. Again, you have to do some work yourself. \$\endgroup\$ Commented Aug 25, 2021 at 6:24

Not the answer you're looking for? Browse other questions tagged or ask your own question.