This is a rasterization problem. For linear inequalities, the regions are polygonal. The classical approach to polygon rasterization is scan conversion; you may learn about it in any textbook on computer graphics.
Preparatory to scan conversion, you will need to convert your linear inequalities (defining half-planes) into an oriented list of edges. One way is to compute the convex hull of the dualized inequalities.
For general polynomial inequalities, the problem is considerably harder. Scan conversion with plane sweeping and active edge lists will do the trick. Detecting sweep events requires you to solve systems of polynomial equations in two variables, which may be done with a combination of numerical and symbolic methods. One way uses resultants and root finding to determine the intersection events up front. Another approach is more incremental, less efficient but significantly simpler:
Assume a bounding box is given. For each y going from the maximum to the minimum, incrementally evaluate the coefficients of the polynomials as functions of y using forward differencing, yielding polynomials depending only on x for each scan-line. Their roots will serve as the rasterization bounds for that scan-line, and you can learn how to find them in a book on numerical methods (Numerical Recipes is popular).
Edit: When I wrote my answer, the last inequality had subscripts rather than superscripts, so it appeared linear. Now I see it is quadratic. That still makes it a rasterization problem although a non-polygonal one. For this specific case you can use circle rasterization to cut out the annulus. In the most general case, you will need polynomial techniques along the lines of my last paragraph. I added some more detail on this.