How to improve performance of huge cellfun + inpolygon?

Dear all,

I have these variables:

B: matrix 19,381x31,768, contains my query points (this variable is not included in the 'test.mat')

lat: vector 19,381x1, contains latitudes of the query points

lon: vector 31,768x1, contains longitudes of the query points

PgonsLat: 5x5 cell array, each cell contains 5 latitudes that define the edges of a polygon

PgonsLon: 5x5 cell array, each cell contains 5 longitudes that define the edges of a polygon

I want to find which query points lie inside each of the polygons (the resulting 'inp1' needs to be cell array). So I use inpolygons function (https://uk.mathworks.com/matlabcentral/fileexchange/7187-inpolygons?s_tid=srchtitle_support_results_2_inpolygons), which has the same syntax as MATLAB inpolygon function (but is faster, I guess).

I wrote this code:

coord = NaN(length(lat)*length(lon),2);

Unrecognized function or variable 'lat'.

for i=1:length(lat)
    for j=1:length(lon)
        coord(j+((i-1)*(length(lon)-1)),:)=[lat(i),lon(j)];
    end
end

checkInsidePolygon1=@(llat,llon) inpolygons(coord(:,1),coord(:,2),llat,llon);

inp1=cellfun(checkInsidePolygon1, PgonsLat,PgonsLon,'UniformOutput',false);

It works, but the performance is poor. The line below takes 184 seconds.

>> tic
inp1=cellfun(checkInsidePolygon1, PgonsLat,PgonsLon,'UniformOutput',false);
toc
>> Elapsed time is 184.079170 seconds.

Unfortunately this line is inside a loop, and I need to iterate over 35,000 times. So it would take weeks to finish.

Is there any way to speed this up?

Answer 1

[Vedant Shah]

Hi @kinga kulesza,

To determine which query points, lie inside each polygon using the ‘inpolygons’ function, the provided code takes longer time because it involves checking all coordinate pairs against each polygon making it infeasible to run for 35000 iterations.

To improve performance, the number of points tested against each polygon can be substantially reduced by applying a bounding box filter. A bounding box is the smallest axis-aligned rectangle that fully contains a polygon. By limiting the point-in-polygon checks to only those coordinates that fall within this bounding box, unnecessary computations are avoided.

Below is an optimized code snippet, which demonstrates the above-mentioned approach:

% Generate the full coordinate grid once 
[LonGrid, LatGrid] = meshgrid(lon, lat); 
coord = [LatGrid(:), LonGrid(:)]; 

%  Preallocate the output cell array 
inp1 = cell(size(PgonsLat)); 

% Loop through each polygon 
for i = 1:size(PgonsLat, 1) 
    for j = 1:size(PgonsLat, 2) 
        % Extract polygon vertices 
        polyLat = PgonsLat{i, j}; 
        polyLon = PgonsLon{i, j}; 
        
        % Compute bounding box of the polygon 
        latMin = min(polyLat); latMax = max(polyLat); 
        lonMin = min(polyLon); lonMax = max(polyLon); 
        
        % Filter points within the bounding box 
        mask = coord(:,1) >= latMin & coord(:,1) <= latMax & ... 
            coord(:,2) >= lonMin & coord(:,2) <= lonMax; 
        coordSubset = coord(mask, :); 
        
        % Run inpolygons only on the filtered subset 
        inside = inpolygons(coordSubset(:,1), coordSubset(:,2), polyLat, polyLon); 
        
        % Store result in full-size logical array 
        result = false(size(coord,1),1); 
        result(mask) = inside; 
        inp1{i,j} = result; 
        end 
end 

The time elapsed while executing the above code is significantly reduced as can be seen below:

This method significantly improves efficiency by reducing the number of point-in-polygon checks, thereby lowering computational time and memory usage. It remains scalable across thousands of polygons or iterations.

For more information refer to the following documentations:

• https://www.mathworks.com/help/releases/R2022b/matlab/ref/meshgrid.html
• https://www.mathworks.com/help/releases/R2022b/matlab/ref/cell.html

[kinga kulesza]

It's been a while... but thanks anyway!