I'm not sure if this answer is correct. Please feel free to modify it...
In a first step, we perform a variable substitution $r=Rx$ in order to eliminate $R$ from the integration bounds. We have $f(R) = R \int_0^1\!dx\,e^{2\pi i (a R^2 x^2 + b R x)} = R \,\tilde f(R).$
As the integrand is holomorphic, we are allowed to deform the integration contour in the complex plain $z = x + i y$. The stationary point is at $z= -b/2aR$.
As the integral starts at $z=0$, we choose a contour on which the phase does not change. This path is (approximately) given by $x=\pm y$ where $\pm$ corresponds to the sign of $ab$.
Discussing the case with $ab> 0$:
Case 1: (a>0, b>0) From the contribution close to $z=0$, we thus have with $z= e^{i \pi/4} t$ (for $a>0$) $\tilde f(R) \sim e^{i\pi/4} \int_0^c\;dt\,e^{-2\pi a R^2 t^2} \sim\frac{e^{i\pi/4}}{2\sqrt{2a}}. $
Similar, for $z=1$, we choose $z=1+e^{i\pi/4} t$ and we have a contribution from the end of the contour with the value $\tilde f(R) \sim -e^{i\pi/4} e^{2\pi i a R^2} \int_c^0\;dt\,e^{-2\pi a R^2 t^2} \sim-\frac{e^{i\pi/4}e^{2\pi i a R^2}}{2\sqrt{2a}}.$
In total, we have $f(R) \sim R \frac{e^{i\pi/4}}{2\sqrt{2a}} (1-e^{2\pi i a R^2}). $
Case 2: (a<0, b<0)
Everything goes along the line of Case 1, we just need to integrate in the direction with $z=e^{-i\pi/4}t$ with the end-result $f(R) \sim R \frac{e^{-i\pi/4}}{2\sqrt{2|a|}} (1-e^{2\pi i a R^2}). $
Case 3: $a b \leq 0$
In this case, we have to take the saddle point at $z=-b/2aR$ into account (as the integration contour passes through it. The integral is dominated from the stationary point and the asymptotic result reads using $z=-b/2aR + e^{\mathop{\rm sgn}(a)i\pi/4} t$ $ f(R) \sim R e^{i b^2\pi/2a} e^{\mathop{\rm sgn}(a)i\pi/4} \int_{-c}^c e^{-2\pi|a| R^2 t^2} \sim e^{i b^2\pi/2a} \frac{e^{\mathop{\rm sgn}(a)i\pi/4} }{2 \sqrt{|a|}}.$