# Find grid points interior to a closed curve

Consider a simply structured grid $G$ represent a domain $\Omega$, e.g. Cartesian (2D/3D) or Polar\Spherical etc. For simplicity the following discussion is about 2D case.

Let $\Gamma\subset \Omega$ be a simply closed curve.

I'm looking for an efficient way to find all points of the grid that lie inside the $\Gamma$.

If $\Gamma$ is a circle or an ellipse then the task is quite easy: represent the grid points in polar\elliptical coordinates and compare the radius\"elliptical radii" and decide.

For a more general case, consider a parametric representation of a starshaped $\Gamma$ and find a matching between the parameter and the radius at each point of the grid and then compare.

However, this algorithm is very slow for a very fine grid.

A possible optimization is to start from a coarser grid and to distribute the result to a finer one, so that only points that are near $\Gamma$ will be tested on a finer grid.

Is there more sophisticated algorithm or another optimization for the one above? Any ideas?

Starshaped curves are not a limitation, any other not too limited (like a polygon) class of curves is acceptable. Another representation of the curve, e.g. bezier curves, level set curves or something else is in general fine. Please explain how to construct such a curve and keep in mind that it should be possible to construct local coordinates at each point. I'm also differentiating the function along this curve, so the representation of the curve should be differentiable.

Optionally, consider, if you can\know, also a non smooth $\Gamma$. The ambiguity of representation and non-differentiability at finitely many points of where the curve isn't smooth enough should not be a problem for my needs.

• Your 'more general' approach seems to assume the curve is convex or bounds a starshaped region? Commented Dec 6, 2015 at 16:01
• yes, starshaped. Commented Dec 6, 2015 at 16:30

First approximate your curve by a piecewise linear one (i.e. by a polygon). There are several ways to do this. One approach is to observe that a polyline is a spline with degree = 1, so all the results of spline approximation theory apply. Carl deBoor's book "A Practical Guide to Splines" has a chapter on polyline approximation. Another approach is to first approximate with biarcs, and then (optionally) approximate the biarcs with a polyline.

After doing this, you have reduced your problem to "point in polygon" testing, which is a topic that has been widely studied. You can get a start by reading the Wikipedia page on the subject.

There are two basic approaches: one based on winding number calculations, and the other based on ray casting.

There is code here. I haven't used it myself, but I've used other stuff from the same site and found it to be very high quality.

Your case is slightly different because you want to test a large number of points lying on a grid, not just one. So, for you, ray casting might be a good approach.If a given ray contains a large number of your points, then, once you've intersected this ray with the polygon, you can classify all the points that lie on the ray, all at once.

If you want to do the calculations really fast (but not very accurately), use GPU code. Assign a color to a test point $P$, scan convert the polygon, and see if $P$ changes color. If it does, then it was inside the polygon.

• I've heard about ray-casting or something very similar. I never understood how to use it in my case. Say you have $r(t)=(x(t),y(t))$ with $t\in[0,1]$ and a Cartesian $n\times n$ grid on $[x_1,x_n]\times[y_1,y_n]$. Take $(x_j,y_1)\to(x_j,y_n)$ ray, how you find all the points where the ray crosses $r(t)$? The accuracy is important. Also, how do you efficiently approximate $r(t)$ as a polygon when $r(t)$ may doesn't conform the grid so that all inside points are still inside. Let ask it even simpler, how do you make a staircase approximation of general shaped curve? Commented Dec 9, 2015 at 12:09
• Regarding intersections ... if you have a PWL approximation, then it's just intersecting lines, so obviously easy. Or, in some cases, you can use the outline curve $r$ in its original form. If you have Bezier curves, for example, you only have to solve cubic equations to find intersections. Commented Dec 11, 2015 at 13:03
• Regarding construction of PWL approximations, see additions to the answer. The PWL approximation consists of arbitrarily sloped line segments, not just horizontal or vertical ones, so it doesn't look like a staircase. It's piecewise linear, not piecewise constant. Commented Dec 11, 2015 at 13:12
• The question is not how you do construct piecewise linear approximation of a curve, but how do you do it so that all interior points remains interior. Consider not necessary convex curves or curves that has points with high eccentricity. Note also that some points of the grid may be extremely close to the curve. Commented Dec 11, 2015 at 13:40
• Rootfinding: chebfun.org/examples/roots/index.html Commented Dec 14, 2015 at 13:28