Fast ray casting in MATLAB

Problem

I've implemented in MATLAB a ray casting algorithm to detect when a 2D point is inside a generic but non self-intersecting 2D polygon. The problem is that my implementation is not sufficiently fast for my applications, so that I'm wonder if some optimization (or maybe some approximation) to my code is possible to reduce the execution time.

Code

function flag = shapeCasting_MATTEO(y, S)
    p0   = 1.1 * max(S) * ones(2, 1);           % a point outside S
    ray  = [p0 y];
    flag = rayCasting(ray, S);
end

function flag = rayCasting(ray, S)
    V   = [S' S(1:2)']';
    hit = 0;
    for j = 1 : 0.5 * size(S, 1)
        Edgej = [V(2*j-1 : 2*j) V(2*(j+1)-1 : 2*(j+1))];
        hit   = hit + hitOrMiss(ray, Edgej);
    end
    if mod(hit, 2) == 0
        flag = 0;
    else
        flag = 1;
    end
end

function flag = hitOrMiss(s1, s2)
    V1 = s1(:, 1); V2 = s1(:, 2);
    V3 = s2(:, 1); V4 = s2(:, 2);
    
    A = [V2 - V1, V3 - V4];
    if abs(det(A)) < 1e-7
        flag = 0;
        return 
    end

    alpha =  A \ (V3 - V1);

    if 0 <= alpha(1) && alpha(1) <= 1 && ...
       0 <= alpha(2) && alpha(2) <= 1
        flag  = 1;
    else
        flag = 0;
    end
end

shapeCasting_MATTEO: this function takes the 2D test point y (encoded as a 2x1 vector) and the polygon S (encoded as a 2nx1 vector) as input. The polygon is encoded through its 2D vertices V1, ..., Vn, which are stacked inside S. The output is flag, which gets value 1 if y is inside and 0 otherwise. The first line defines a ray, which connects a point p0 outside the polygon and the test point y, which is then pass to my rayCasting function. Here ray is a 2x2 matrix, where the first column is point outside the polygon and the second column is the test point.

rayCasting: this function counts how many times the input ray intersect with the contour of the polygon S. The number of intersections is stored in hit, and so if hit is even, then the point y is outside S, otherwise is outside. To count the number of intersections, the algorithm check each one of the n edges of the polygon. To take into account the last edge which closes the polygon, a virtual vertex Vn+1 =V1 is introduced through the definition of the 2(n+1)x1 vector V. Each edge is identified by two consecutive vertices, which are stored in the 2x2 matrix Edgej. In particular, the first column of Edgej is the vertex Vj, while the second column is the vertices Vj+1. The function hitOrMiss(ray, Edgej) checks if ray intersects Edgej.

hitOrMiss: this function checks if the two input segments s1, s2 intersect with each other. I think that hitOrMiss is the critical point of my shapeCasting_MATTEO algorithm, and I'm not sure if my implementation of hitOrMiss is optimal.

To check if there is an intersection between s1 and s2, hitOrMiss searches for a common point in s1 and s2. If a common point exists, then it can be written as a point over s1 in the form (1 - alpha(1)) * V1 + alpha(1) * V2 or as as point over s2 in the form (1 - alpha(2)) * V3 + alpha(2) * V4, where:

V1, V2 are the vertices of s1;
alpha(1) is a scalar between 0 and 1;
V3, V4 are the vertices of s2;
alpha(2) is a scalar between 0 and 1.

Hence, the objective is to solve the linear equation

(1 - alpha(1)) * V1 + alpha(1) * V2 = (1 - alpha(2)) * V3 + alpha(2) * V4

with respect to alpha(1) and alpha(2) and check if alpha(1) and alpha(2) are in the right ranges. The solution of the previous linear equation is obtained with the line code alpha = A \ (V3 - V1) (where A = [V2 - V1, V3 - V4]), which can proved with simple algebraic considerations. To improve the numerical robustness of the method, a pre-eliminary check over the determinant of A is done. If A is ill-conditioned, which happens when s1 and s2 are almost parallel, then the method says that no intersection occours. Clearly this choice may lead to some errors, but some tests show that overall the method works pretty well.

Minimal workable example

clc
close all
clear

% generate a test polygon in [0, 1]^2 with a big number of vertices
n = 1e4;
S = NaN(2*n, 1);
for i = 1 : n
    S(2*i - 1) = cos(2 * pi * (i - 1) / (n - 1));
    S(2*i)     = sin(2 * pi * (i - 1) / (n - 1));
end

% generate a test point
y = rand(2, 1);

% run the algorithm multiple times
nTest   = 100;
runTime = NaN(1, nTest);
for i = 1 : nTest
    tic
    shapeCasting_MATTEO(y, S);
    runTime(i) = toc;
end
runTime = round(1e3 * mean(runTime));

% visual check
Sx = NaN(1, n);
Sy = Sx;
for i = 1 : n
    Sx(i) = S(2*i - 1);
    Sy(i) = S(2*i);
end

figure();
hold on
plot([Sx Sx(1)], [Sy Sy(1)], 'r');
plot(y(1), y(2), '.k');
axis equal
if shapeCasting_MATTEO(y, S)
    title(['Point inside - Execution time: ' num2str(runTime) ' ms']);
else
    title(['Point outside - Execution time: ' num2str(runTime) ' ms']);
end

Solution

Using profile on before your test and profile off; profile viewer afterwards, we can see the bottlenecks in your code. The main slow-down was this line:

alpha = A \ (V3 - V1)

You're calling this thousands of times, and it's a slow operation at the best of times, so we need to leverage some MATLAB-esque matrix operations to make this call more efficient.

There was a recent MathWorks blog about using page-wise functions which were expanded in release R2022b, and offer performance benefits for repeated matrix operations. In particular the pagemldivide looks very applicable to your problem

X = pagemldivide(A,B) computes the left matrix divide of each page of N-D array A into each page of N-D array B. Each page of the output array X is given by X(:,:,i) = A(:,:,i) \ B(:,:,i).

To take advantage of this, we have to create A and B where each "page" (or layer in the 3rd dimension) is defined so instead of the current loop which is (in abstract) something like

for j = 1:N
  alpha(j) = A \ (V3 - V1);
end

We can simply do alpha = pagemldivide( A, B );

Further vectorisation can give us additional runtime improvements

You are recreating V1 and V2 for every hitOrMiss test but they're always defined by ray which doesn't change so take this out of the loop
You can create all of the Edgej vectors in one hit as reshapes of the input S instead of doing this with indexing in a loop. Reshapes are very computationally cheap.

On my PC, this reduces the time of your 100-fold test from 19.24 sec to 0.92 sec, or a performance improvement of ~20x.

If you want to test these back-to-back you should use rng(0) at the top of your test script to get a repeatable set of random values for testing.

Here is the resulting code, I've combined the rayCasting and hitOrMiss functions, see the comments for details:

function flag = rayCasting(ray, S)
    % V1 and V2 are unchanging, defined from ray
    V1 = ray(:, 1); V2 = ray(:, 2);
    % V3 is 2x1xN so each pair of neighbours from S is its own "page"
    V3 = reshape(S,2,1,[]);
    % V4 is a wrap-around of V3, move all pairs by one page
    V4 = V3(:,:,[2:end,1]);
    % Compute B in one hit, no loop required now we have paged V1 and V3
    B = V3 - V1;

    % Right-hand side of A is defined as diff of V3 and V4 as before
    AR = V3 - V4;
    % Left-hand side of A was always fixed by ray anyway, quick repmat
    AL = repmat(V2-V1, 1, 1, size(AR,3));
    % Define A, 2x2xN square matrices "paged" for each check
    % where N = 0.5*size(S,1)
    A = [AL, AR];

    % check for ill-conditioned matrices from the abs determinant
    % faster to compute the det manually as it's only 2x2 and we can
    % vectorise it this way
    idx = abs( A(1,1,:).*A(2,2,:) - A(1,2,:).*A(2,1,:) ) >= 1e-7;
    % For all well conditioned indices, combute A\b
    alpha = squeeze( pagemldivide( A(:,:,idx), B(:,:,idx) ) );
    % Count the number of elements where both rows of alpha are 0<=a<=1
    hit = nnz( (0 <= alpha(1,:)) & (alpha(1,:) <= 1) & (0 <= alpha(2,:)) & (alpha(2,:) <= 1) );   

    % Output flag
    if mod(hit, 2) == 0
        flag = 0;
    else
        flag = 1;
    end
end