Since the outer integral $\int_0^{2\pi} F(\theta) d\theta$ has a periodic integrand, a trapezoidal rule should work nicely, picking up extra accuracy when $F(\theta)$ is smooth.
For the inner integral you want to take advantage of the known "weight function" $r$ that appears in the integrand, $F(\theta) = \int_0^{r(\theta)} f(r,\theta) r\; dr$.
One way to implement this is to set some threshold $r_\epsilon \gt 0$ around the origin and use a weighted Gaussian quadrature rule on the inner radial segment $[0,r_\epsilon]$, combined if necessary with a compound trapezoidal rule on the exterior segment $[r_\epsilon,r(\theta)]$.
The common Gaussian quadrature schemes are either for the constant weight function $w(x) = 1$ (Gauss-Legendre) or for much fancier weight functions than $w(x) = x$. I looked around the web but did not find a quick reference for this case, although it appears as a homework-type exercise on Math.SE.
Perhaps we should be able to modify the computation of the Gauss-Legendre nodes (roots of Legendre polynomials) to get the nodes and weights for this case, but I have done the 1-, 2- and 3-point rules from "scratch" as concisely outlined in this note.
For an $N$-point rule we need the (monic) polynomial $q(x)$ of degree $N$ which is orthogonal to lower degree polynomials on $[0,1]$ when integrated with weight function $w(x) = x$:
$$ \int_0^1 q(x) x^k\; x dx \text{ for } k=0,1,\ldots,N-1 $$
When $N=1$ we get $q(x) = x - \frac{2}{3}$, so the node is $x=\frac{2}{3}$ and weight $\omega = \frac{1}{2}$ (because $\int_0^1 x dx = \frac{1}{2}$ should be exactly integrated). This amounts to a quadrature rule:
$$ \int_0^1 f(x) x dx \approx \frac{1}{2} f\left(\frac{2}{3}\right) $$
Adapting to the case of integrating over $[0,r(\theta)]$ amounts to throwing in factors:
$$ \int_0^{r(\theta)} f(r) \; r dr \approx \frac{r(\theta)}{2} f\left(\frac{2r(\theta)}{3}\right) $$
For $N=2$ we get a quadratic $q(x) = x^2 - 1.2x + 0.3$ whose roots are two points in $[0,1]$, and the respective weights are found by requiring $1,x$ to be exactly integrated:
$$ x_\pm = \frac{1}{10} (6 \pm \sqrt{6}) \; , \; \omega_\pm = \frac{1}{36} (9 \pm \sqrt{6}) $$
For $N=3$ the cubic polynomial $q(x) = x^3 - \frac{12}{7} x^2 + \frac{6}{7} x - \frac{4}{35}$ turns out to one of those irreducible cases whose three real roots cannot be expressed in terms of radicals without using complex constants. So in lieu of symbolic expressions we will note the single precision values for the three roots:
$$ x_0 \approx 0.21234, x_1 \approx 0.59053, x_2 \approx 0.91141 $$
which can be refined if necessary by Newton iterations to any desired precision. The corresponding weights can then be found by solving the linear system:
$$ \omega_0 + \omega_1 + \omega_2 = \frac{1}{2} $$
$$ x_0 \omega_0 + x_1 \omega_1 + x_2 \omega_2 = \frac{1}{3} $$
$$ x_0^2 \omega_0 + x_1^2 \omega_1 + x_2^2 \omega_2 = \frac{1}{4} $$
Using the single-precision values of nodes $x_i$ gives approximate weights:
$$ \omega_0 \approx 0.06983, \omega_1 \approx 0.22924, \omega_2 \approx 0.20093 $$
so that:
$$ \int_0^1 f(x) x dx \approx \omega_0 f(x_0) + \omega_1 f(x_1) + \omega_2 f(x_2) $$
and this quadrature can be scaled to interval $[0,r(\theta)]$ in a similar manner to that shown for the 1-point rule above.
