Probability of a relation on the uniform distribution of points over 2D space Assume a set of nodes is scattered over a 2D surface $\mathcal{S}$ so that for any given $\mathcal{A} \subset \mathcal{S}$, the number of nodes inside $\mathcal{A}$ follows a Poisson distribution with parameter $|\mathcal{A}| \rho$, where $|\mathcal{A}|$ shows the area of the subset $\mathcal{A}$ and $\rho$ is the intensity of the points (average number of points per unit area).
We are only interested in the points inside a given circle with radius $r$. The number of nodes inside the circle is a Poisson variable with parameter $\rho \pi r^2$. We pick two nodes from inside the circle at random. Let $d_1$ and $d_2$ show the distance of the first and the second node from the center of the circle.

How can I compute the probability of the event:
$$
{d_1}^2 < \frac{{d_2}^2}{A(1+B{d_2}^2)}
$$
where $A$ and $B$ are constants.
Edit:


*

*Assume $A > 0$ and $B > 0$.

*I am interested in the process itself, not the points generated by the process (as whuber described in his answer below).

*How about the case that ${d_1}^2$ and ${d_2}^2$ is replaced with ${d_1}^\alpha$ and ${d_2}^\alpha$ for $\alpha > 2$ (I guess, this modifies the problem since ${d_1}^\alpha$ and ${d_2}^\alpha$ are not uniformly distributed anymore).
 A: There are at least two interpretations: one concerns the actual points generated by this process and the other concerns the process itself.
If a realization of the Poisson process is given and pairs of points are to be chosen from that realization, then there is nothing to be done except systematically compare all distances to all other distances (a double loop over the points).
Otherwise, if the procedure is intended to consist of (i) creating a realization of the process and then (ii) selecting a pair of points at random, then the assumptions imply the two points are selected uniformly and independently from the circle.  The calculation for this situation can be performed once and for all.
Notice that the squared distances $r_1 = d_1^2$ and $r_2 = d_2^2$ are uniformly distributed, whence the desired probability is
$$p(a,b) = \Pr\left(d_1^2 \lt \frac{d_2^2}{a(1 + b d_2^2)}\right) = \int_0^1 d r_2 \int_0^{\max(0, \min(1, r_2 / (a(1 + b r_2))))} d r_1.$$
The $\max$ and $\min$ can be handled by breaking into cases.  Some special values of $a$ and $b$ have to be handled.  Because the integration is a square window over a region generically bounded by lines and lobes of a hyperbola (with vertical axis at $1/(ab)$ and horizontal axis at $-1/b$), the result is straightforward but messy; it should involve rational expressions in $a$ and $b$ and some inverse hyperbolic functions (that is, natural logarithms).  I had Mathematica write it out:
$$\begin{array}{ll}
 \frac{b+1}{b} & \left(-1\leq a<0\land \frac{1}{a}-b\leq 1\land b<-1\right)\\ &\lor \left(a<-1\land \frac{1}{a}-b<1\land b<-1\right) \\
 -\frac{1}{b (a b-1)} & \frac{1}{a}-b=1\land a<-1 \\
 \frac{a^2 b+2 a b+a-2}{2 (a b-1)} & b=0\land a>0\land \frac{1}{a}-b>1 \\
 \frac{b-\log (b+1)}{a b^2} & a>0\land \frac{1}{a}-b\leq 1\land b>-1 \\
 \frac{a b^2+a b-a b \log (b+1)-b+\log (b+1)}{a b^2 (a b-1)} & a>0\land \frac{1}{a}-b\leq 1\land b\leq -1 \\
 \frac{\log (1-a b)}{a b^2} & a>0\land \frac{1}{a}-b>1\land b\leq -1 \\
 \frac{a b^2+a b+\log (1-a b)}{a b^2} & \left(-1<b<0\land a>0\land \frac{1}{a}-b>1\right) \\
 & \lor \left(b>0\land a>0\land \frac{1}{a}-b>1\right) \\
 \frac{b-\log ((-b-1) (a b-1))}{a b^2} & a<0\land \frac{1}{a}-b>1
\end{array}$$
Numeric integration and simulation over the ranges $-2 \le a \le 2$ and $-5 \le b \le 5$ confirm these results.
Edit
The modified question asks to replace $d_i^2$ by $d_i^\alpha$ and assumes $a$ and $b$ both positive. Upon making a substitution $r_i = d_i^\alpha$, the region of integration remains the same and integrand becomes $(2/\alpha)^2(r_1 r_2)^{2/\alpha-1}$ instead of $1$.  Writing $\theta = \alpha/2$, we obtain
$$\frac{1}{2} a^{-1/\theta } \, _2F_1\left(\frac{1}{\theta },\frac{2}{\theta };\frac{\theta +2}{\theta };-b\right)$$
when $(a>0\land a<1\land a b+a\geq 1)$ or $a\geq 1$ and otherwise the result is
$$-a^{\frac{1}{\theta }} \left(\frac{1}{1-a b}\right)^{\frac{1}{\theta }}+\frac{1}{2} a^{\frac{1}{\theta }} (1-a b)^{-2/\theta } \, _2F_1\left(\frac{1}{\theta },\frac{2}{\theta };\frac{\theta +2}{\theta };1+\frac{1}{a b-1}\right)+1.$$
Here, $_2F_1$ is the hypergeometric function.  The original case of $\alpha=2$ corresponds to $\theta=1$ and then these formulae reduce to the fourth and seventh of the eight previous cases.  I have checked this result with a simulation, letting $\theta$ range from $1$ through $3$ and covering substantial ranges of $a$ and $b$.
A: 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.

