Vectorisation of a function solving two symbolic equations with two unknows and using structures

Hello,

I wanna find the two intersection points between a straight line and a circle. Im not considering the case where there is a single intersection point or none. For this i created a function wich works well for one single straight line. But i wanna use my function for many straight lines at the same time and not use a for loop.

For one straight line, the function solves a two equations two unknows system. This system of equation is defined by the circle equation, and the straight line equation in a 2D space. The parameters i used to define this system is the radius of the circle, the coordinates of the center of the circle, a direction vector of the straight line, and a point on this straight line.

How can i vectorize this function please ?

I attached the full script with a matrices of many straight line i wanna find the intersection with a circle. Here is the function i created.

function intersection_point = intersection_cercle_test(vector_ray, point_ray, radius_circle,  point_center_circle)
% Function wich gives us the intersection points of a
% ray and a circle
% Im not consdering the case where there is no intersection or only one. 

% vector_ray is a 2 line one column matrice containing the direction vector
% of the ray 
% point_ray is a 1 line 2 column matrice containing a point wich is part of
% that ray
% radius circle if the radius of the circle. Its a strictly positiv scalar.
% point_center_circle is a 1 line 2 column matrice containing the point at
% the center of the circle


norm_vector_ray = [vecnorm(vector_ray); vecnorm(vector_ray)];
vector_ray =  vector_ray ./ norm_vector_ray;

syms u v
eqns = [(u-point_center_circle(1))^2 + (v-point_center_circle(2))^2 == radius_circle^2, vector_ray(2,:)*(u-point_ray(:,1))+ vector_ray(1,:)*(point_ray(:,2)-v) == 0];
vars = [v u];
[solv, solu] = solve(eqns,vars);
intersection_point = double([solu solv]);


end

[John D'Errico]

I think the important point is to see that you don't want to vectorize a call that involves solve. Solve will be slow, no matter what. And the mathematics of the solve is trivial, involving nothing more than the solution of a quadratic polynomial, so finding the two roots of the quadratic.

Don't just throw stuff at solve and hope it will be efficient, and then look for ways to vectorize that approach. Instead, use mathematics to solve the problem. @Matt J has shown exactly what you can do.

[Camille]

Thank you for your answear. I didnt know solve was slow.

Answer 1

[Matt J]

function [Intersection1,Intersection2] = intersection_cercle_test(vector_rays, point_rays, radius_circle,  point_center_circle)
% vector_rays is a 2xN matrix containing direction vectors of rays.
% point_rays is a 2xN matrix containing starting points of rays.
% radius circle is the scalar radius of the circle.
% point_center_circle is 2x1 vector coorindate of circle center.


    vector_rays= normalize(vector_rays,1,'n');
    
    point_rays = point_rays-point_center_circle;
    
    B=2*sum(vector_rays.*point_rays,1);
    
    C=vecnorm(point_rays,2,1).^2-radius_circle^2;
    
    discrim=sqrt(B.^2-4*C); %quadratic formula vectorized
    
    Intersection1=point_rays+(-B+discrim)/2.*vector_rays + point_center_circle;
    Intersection2=point_rays+(-B-discrim)/2.*vector_rays + point_center_circle;

end

[Camille]

Thank you very much for your answear. It works very well. Can you explain to me what your doing with line B and below please ?

[Matt J]

The parametric equation of each ray is P(t)=t*d+P0 where d is a unit direction vector_ray and P0 is a point_ray. To lie on a circle of radius R centered at the origin, a point on the ray P(t) must satisfy the equation,

which reduces to a quadratic equation in t,

where B and C are as given in the code. This equation is then solved with the quadratic formula.