$$\int_0 ^ {1} \int_{-\sqrt{x-x^2}} ^ {\sqrt {x-x^2}} (x^2+y^2) ~dy~dx$$
My findings are:
$$\int_?^? \int_?^? r^3 ~dr~d\theta$$
Region is the circle of radius $~\frac{1}{2}~$ centered at $~(\frac{1}{2}, 0)~$, but I cannot figure out the limits for r and $\theta$.
Sketch of region of limit is
There are a number of ways you could go about this problem. One way is to do the coordinate transformation as you have: $$\int_{x=0}^1 \int_{y=-\sqrt{x-x^2}}^{\sqrt{x-x^2}} x^2 + y^2 \, dy \, dx = \iint_R r^3 \, dr \, d\theta,$$ where $R$ is the region of integration parametrized in polar coordinates. This is where you got stuck. To determine $R$, we note that the boundary of the region in Cartesian coordinates is given by $$y^2 = x - x^2,$$ or $$x^2 + y^2 = x.$$ Since the Cartesian to polar transformation is $$r^2 = x^2 + y^2, \quad x = r \cos \theta, \quad y = r \sin \theta,$$ we can substitute directly: $$r^2 = r \cos \theta$$ from which we obtain $$r = \cos \theta.$$ Now consider the interval of integration for $\theta$: if $r \ge 0$, then a natural choice is $\theta \in [-\pi/2, \pi/2]$, hence $$\iint_R = \int_{\theta = - \pi/2}^{\pi/2} \int_{r = 0}^{\cos \theta},$$ since we are integrating over the disk, rather than the boundary.
Alternatively, you can be more sophisticated and note that the region of integration is a circle of radius $1/2$ centered at $(1/2, 0)$; therefore, an initial translation of the form $$(x,y) \mapsto (x + 1/2, y)$$ changes the integral to $$\int_{x=-1/2}^{1/2} \int_{y=-\sqrt{1/4 - x^2}}^{\sqrt{1/4 - x^2}} (x+1/2)^2 + y^2 \, dy \, dx.$$ Now the region of integration is simply a circle of radius $1/2$ centered at the origin, and the polar coordinate conversion is straightforward: $$\int_{r=0}^{1/2} \int_{\theta = 0}^{2\pi} \left((r \cos \theta + 1/2)^2 + (r \sin \theta)^2 \right) r \, d\theta \, dr.$$ The integrand simplifies to $$\frac{r}{4} + r^3 + r^2 \cos \theta,$$ and since the variables of integration are not dependent on each other, this computation is easily performed.