This problem can be solved by decomposing into parts and using the properties of a Poisson process.
It helps to recall how to generate a Poisson point process of intensity $\rho$ on a bounded subset $\newcommand{\A}{\mathcal A}$ of $\mathbb R^2$. We first generate a Poisson random variable $N$ with rate $\rho |\mathcal A|$ where $|\cdot|$ denotes Lebesgue measure, and then we sprinkle these $N$ points uniformly at random inside of $\A$.
This immediately tells us that as long as $N \geq 2$, if we choose two points (without replacement) at random, then these two points will be independent and uniformly distributed on $\A$. When $N < 2$, we have to do something and one natural choice is to define the desired probability as zero. Note that this happens with probability
$$\renewcommand{\Pr}{\mathbb P}\Pr(N < 2) = (1+\rho|\A|) e^{-\rho|\A|} \>.
$$
This is the only part of the problem that depends on the Poisson process intensity.
Probability conditional on ${N \geq 2}$
We are interested in the probability
$$
p(A,B,r) := \Pr\Big( d_1^2 \leq \frac{d_2^2}{A(1+B d_2^2)}\Big) \>,
$$
where $A > 0$, $B > 0$ and $\A = \{x : \|x\|_2 \leq r\}$. Here $d_1$ and $d_2$ are the radii of two of our uniformly distributed points that fall in $\A$.
Note that for a point randomly distributed in the disc of radius $r$, the distribution of the distance from the origin is $\Pr(D \leq d) = (d/r)^2$, from which we can see that $D^2$ has the same distribution as $r^2 U$ where $U \sim \mathcal U(0,1)$. From this, we can restate the probability of interest as
$$
p(A,B,r) = \Pr\Big( U_1 \leq \frac{U_2}{A(1+B r^2 U_2)}\Big) = \iint 1_{(0 < x < 1)} 1_{(0 < y < 1)} 1_{(0 < y < x/(A+ABr^2 x))} \,\mathrm dy \, \mathrm dx \>.
$$
This integral splits into two cases. To calculate it, we need the general integral
$$
\int_0^t \frac{x}{a+bx} \,\mathrm d x = \frac{1}{b} (t - \frac{a}{b} \log(1+bt/a)) \>.
$$
Case 1: $A(1+B r^2) \geq 1$.
Here we see that $u \leq A(1+B r^2 u)$ for $u \in [0,1]$, so
$$
p(A,B,r) = \frac{1}{ABr^2}\Big(1 - \frac{\log(1+B r^2)}{B r^2}\Big) \>.
$$
Case 2: $A(1+B r^2) < 1$.
Here the integral for $p(A,B,r)$ splits into two pieces since $u \geq A(1+B r^2 u)$ on $[A/(1-ABr^2),1]$. Hence we integrate up to $t = A/(1-A B r^2)$ using the general integral and then tack on an addition area of $1-A/(1-ABr^2)$ for the second piece. So, we get
$$
\begin{align}
p(A,B,r) &= \frac{1}{B r^2} \Big(\frac{1}{1-A B r^2} + \frac{\log(1-AB r^2)}{A B r^2}\Big) + 1 - \frac{A}{1 - A B r^2} \\
&= 1 + \frac{1}{B r^2} \Big(1 + \frac{\log(1-AB r^2)}{A B r^2}\Big) \>.
\end{align}
$$
Oftentimes a picture helps; here is one that shows an example of the integration region for each case. Note that $U_1$ is on the $y$-axis and $U_2$ on the $x$-axis.
The final probability of interest is then, of course, $(1 - (1+\rho\pi r^2) e^{-\rho\pi r^2} ) p(A,B,r)$.
An easy generalization
We can easily generalize the result to use a different shaped ball. In fact, for any arbitrary norm on $\mathbb R^2$, the conditional probability $p(A,B,r)$ is invariant as long as we use the ball induced by the norm instead of the circle!
This is because no matter what norm we choose, the squared radius is uniformly distributed. To see why, let $\delta(\cdot)$ be a norm on $\mathbb R^2$ and $B_\delta(r) = \{x: \delta(x) \leq r\}$ the ball of radius $r$ under the norm $\delta$. Note that $rx \in B_\delta(r)$ if and only if $x \in B_\delta(1)$. The scaling up or down of the unit ball is a linear transformation and by a standard fact about Lebesgue measure, the measure of a linear transformation $T$ of $B_\delta(1)$ is
$$
|B_\delta(r)| = |T B_\delta(1)| = |\det(T)| |B_{\delta(1)}| = r^2 |B_\delta(1)| \>,
$$
since $T(x) = r x = (r x_1, r x_2)$ in this case.
This shows that if $D = \delta(X)$ for $X$ uniformly distributed in $B_\delta(r)$, then
$$
\Pr(D \leq d) = \frac{|B_\delta(d)|}{|B_\delta(r)|} = (d/r)^2 \>.
$$
The eagle-eyed reader will note that we've only used the homogeneity of the norm here, and so a similar result will hold in general for uniform distributions on classes of sets closed under a homogeneous transformation.
Here is a picture with two points selected. The norms shown are the Euclidean norm, $\ell_1$ norm, $\sup$ norm, and the $\ell^p$ norm for $p = 5$. Each unit ball is outlined in black, and the largest ball within which the two randomly selected points lie is drawn in the corresponding color.
The conditional probability $p(A,B,r)$ is the same for each picture when the distance is measured using the corresponding norm.
