I have a grid map, constructed using Matlab. I am placing some polygons on it. How can i extract the grid coordinates that are inside these polygons? thanks..
% construct grid
MAX_X=10;
MAX_Y=10;
MAX_VAL=10;
MAP=2*(ones(MAX_X,MAX_Y));
axis([1 MAX_X+1 1 MAX_Y+1])
grid on;
hold on;
%obst 1
x = [1 1 4 4];
y = [1 11 11 1];
plot([x x(1)],[y y(1)],'r-');
% obst 2
x = [7 7 11 11];
y = [11 1 1 11];
plot([x x(1)],[y y(1)],'r-');
first make a closed line:
x = [x x(1)]
y = [y y(1)]
and a matrix of values of X and Y at different points:
Y = repmat((1:MAX_Y)',[1,MAX_X])
X = repmat(1:MAX_X,[MAX_Y,1])
Then you can use inpolygon:
MAP = inpolygon(X,Y,x,y)
That should also work on non-convex polygons.
Dealing with non-convex polygons will be tricky, so perhaps you should start off by dividing your polygon into multiple convex polygons. You can use the interior angle to decide where a corner is non-convex, and apply cuts that involve these corners at the minimum. But perhaps matlab does support some triangulation algorithm you can use.
Once you have your convex polygons, you can think of them as intersections of half spaces. If A and B are two corners which span a line, and P is a grid point, you can use the sign of the determinant
| Ax Bx Px |
| Ay By Py |
| 1 1 1 |
to decide on which side of the line the point lies. Which sign is which depends on the order in which you walk around your polygon, but if you process corners in order, then a point P lies within the convex polygon iff the sign never changes. Horizontal or vertical lines won't be special cases in this formulation, and you won't have divisions either which is good for performance and might help with exactness as well.
If you don't want to iterate over all grid points using that method, you can come up with various optimizations. One would be to precompute the cross product A × B of the two corners spanning each line. The dot product between that cross product and a point P is equal to the above determinant (i.e. det(A,B,P)=(A×B)·P), so instead of the full determinant you now only have to compute three products and two sums for each point-line combination.
If you want to optimize this even further, you'd probably best look at something like Bresenham's line algorithm to compute the points at the border of the polygon, and then simply enumerate all points on a horizontal (or vertical) line between the matching border points.
Unless you perform all your computations with integers or other exact numbers, rounding issues might be a major concern in all of this. You have to decide whether or not to count border points, and also make sure to count points inside the original polygon but on the border of its convex portions exactly once. How much effort this requires depends to a large degree on how your input looks like.
A related post is here: What is an simple way to compute the overlap between an image and a polygon?
inpolygon may work for your application, but will only return true if the center point of a particular grid box is inside the polygon. It won't work if you want to know if any part of a grid box is within the polygon. An example of using inpolygon is shown by the asker in the above post.
