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.
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.