Matching red to blue dots

I have two red points, $r_1$ and $r_2$, and two blue points, $b_1$ and $b_2$. They are all placed randomly and uniformly in $[0,1]^2$.

Each dot points to the closest dot from another colour; closest is defined wrt the Euclidean distance. We use $x \to y$ to indicate dot $x$ points to dot $y$.

If $r_1 \to b_1$, what is the probability that $r_2 \to b_1$ too?

NOTE it must be larger than 1/2 because $r_1 \to b_1$ tells us in a way that $b_1$ is likely to have a centric location, and thus is likely than it is closer to $r_2$ too than $b_2$.

From what you have written, J assume that all random dots in the question are distributed independently.

Now, suppose $b_1 = (x_1, y_1)$ and $b_2 = (x_2, y_2)$ are fixed. Thus for a random dot $r_1$, uniformly distributed on $[0; 1]^2$, the probability, that it lies closer to $b_1$, than to $b_2$ is $\mu(\{(x, y) \in [0;1]^2| (x_1 + x_2 - 2x)(x_2 -x_1) + (y_1 + y_2 - 2y)(y_2 -y_1)>0\})$, where $\mu$ stands for Lebesgue measure. Now, as $r_1$ and $r_2$ are independent, then the probability, that both $r_1$ and $r_2$ lie closer to $b_1$, than to $b_2$ is $(\mu(\{(x, y) \in [0;1]^2| (x_1 + x_2 - 2x)(x_2 -x_1) + (y_1 + y_2 - 2y)(y_2 -y_1)>0\}))^2$.

Now as $b_1$ and $b_2$ are also independent and uniformly distributed, we can conclude, that in our initial problem $$P(r_2 \to b_1|r_1 \to b_1) = \frac{\int_0^1 \int_0^1 \int_0^1 \int_0^1 (\mu(\{(x, y) \in [0;1]^2| (x_1 + x_2 - 2x)(x_2 -x_1) + (y_1 + y_2 - 2y)(y_2 -y_1)>0\}))^2dx_1dy_1dx_2dy_2}{\int_0^1 \int_0^1 \int_0^1 \int_0^1 \mu(\{(x, y) \in [0;1]^2| (x_1 + x_2 - 2x)(x_2 -x_1) + (y_1 + y_2 - 2y)(y_2 -y_1)>0\}))dx_1dy_1dx_2dy_2}$$

Starting off in the same manner as Yanior Weg, assume $b_1$ and $b_2$ are fixed. Then $P(r_2 \to b_1|r_1 \to b_1) = \frac{P(r_2 \to b_1 \bigcap r_1 \to b_1)}{P(r_1 \to b_1)} = \frac{P(r_2 \to b_1 \bigcap r_1 \to b_1)}{\frac{1}{2}} = 2 \cdot P(r_2 \to b_1 \bigcap r_1 \to b_1) = 2(P(r_1 \to b_1))^2$

In this Desmos plot, the depiction of this can be seen (I'll be referring to this later, so it may be useful to have it open). $P(r_1 \to b_1)$ is the area of the shaded region inside the square. By limiting $x_2, y_2$ such that $x_1 < x_2 < 1$, $y_1 < y_2 < 1$, later calculations can be simplified. The final answer is then $$\int_0^1\int_0^1\int_{x_1}^1\int_{y_1}^14\cdot2A^2 dy_2dx_2dy_1dx_1 = 8\int_0^1\int_0^1\int_{x_1}^1\int_{y_1}^1A^2 dy_2dx_2dy_1dx_1$$, where $A$ is the area of the shaded region inside the square. $2A^2$ is multiplied by $4$ to account for $x_2<x_1$ and $y_2<y_1$.

The line separating the shaded region has equation $$y = f\left(x\right)=-\frac{x_{1}-x_{2}}{y_{1}-y_{2}}x+\frac{\left(x_{1}-x_{2}\right)\left(x_{1}+x_{2}\right)}{2\left(y_{1}-y_{2}\right)}+\frac{y_{1}+y_{2}}{2}$$

This intersects $y = 1$ at $$x = I_1 = \frac{\left(y_{1}-y_{2}\right)\left(y_{1}+y_{2}-2\right)}{2\left(x_{1}-x_{2}\right)}+\frac{x_{1}+x_{2}}{2}$$

and $y = 0$ at $$x = I_2 = \frac{\left(y_{1}-y_{2}\right)\left(y_{2}+y_{1}\right)}{2\left(x_{1}-x_{2}\right)}+\frac{x_{1}+x_{2}}{2}$$

There are now four cases: $(1)\ I_1 < 0, I_2 < 1, \ (2)\ I_1 < 0, I_2 > 1, \ (3)\ I_1 > 0, I_2 < 1$, and $(4)\ I_1 > 0, I_2 > 1$.

Letting $$F(x) = \int_0^x f(t)dt = -\frac{x_{1}-x_{2}}{2\left(y_{1}-y_{2}\right)}x^{2}+\left(\frac{\left(x_{1}-x_{2}\right)\left(x_{1}+x_{2}\right)}{2\left(y_{1}-y_{2}\right)}+\frac{y_{1}+y_{2}}{2}\right)x$$ here $A_i$ represents the area of case $i$:

$A_1 = F(I_2)$

$A_2 = F(1)$

$A_3 = I_1 - F(I_1) + F(I_2)$

$A_4 = I_1 - F(I_1) + F(1)$

From here, $$J = \underbrace{\int_{(1)}A_1^2dy_2dx_2dy_1dx_1}_{J_1}+\underbrace{\int_{(2)}A_2^2dy_2dx_2dy_1dx_1}_{J_2}+\underbrace{\int_{(3)}A_3^2dy_2dx_2dy_1dx_1}_{J_3}+\underbrace{\int_{(4)}A_4^2dy_2dx_2dy_1dx_1}_{J_4}$$, where $(1)$ is the region in $x_1, y_1, x_2, y_2$ such that case 1 happens (restricted to $0 < x1 < x2 < 1$ and $0 < y1 < y2 < 1$), etc.

To simplify, making the substitution $x_s = x_2 + x_1, x_d = x_2 - x_1$ and $y_s = y_2 + y_1, y_d = y_2 - y_1$ helps a lot. The integral would then need to be multiplied by the Jacobian of $\frac{1}{4}$.

For case $1$, the integral can be written out as $$J_1 = \frac{1}{4}\int_0^1 \int_{x_d}^{2-x_d} \left(\int_{1-\sqrt{1-x_{d}x_{s}}}^{x_{d}}\int_{0}^{-\frac{x_{d}x_{s}}{y_{d}}+2}A_1^2dy_{s}dy_{d}+\int_{x_{d}}^{\sqrt{2x_{d}-x_{d}x_{s}}}\int_{0}^{\frac{2x_{d}-x_{d}x_{s}}{y_{d}}}A_1^{2}dy_{s}dy_{d}-\int_{1-\sqrt{1-x_{d}x_{s}}}^{\sqrt{2x_{d}-x_{d}x_{s}}}\int_{0}^{y_{d}}A_1^{2}dy_{s}dy_{d}\right)dx_s dx_d = \frac{1}{32}-\frac{1}{24}\ln(2)$$

For case $2$: $$J_2 = \frac{1}{4}\int_0^1 \int_{x_d}^{2-x_d} \left(\int_{x_{d}}^{\sqrt{x_{s}x_{d}}}\int_{0}^{2-\frac{x_{s}x_{d}}{y_{d}}}A_2^{2}dy_{s}dy_{d}+\int_{\sqrt{x_{s}x_{d}}}^{1}\int_{0}^{2-y_{d}}A_2^{2}dy_{s}dy_{d}-\int_{x_{d}}^{\sqrt{2x_{d}-x_{d}x_{s}}}\int_{0}^{\frac{2x_{d}-x_{d}x_{s}}{y_{d}}}A_2^{2}dy_{s}dy_{d}-\int_{\sqrt{2x_{d}-x_{d}x_{s}}}^{1}\int_{0}^{y_{d}}A_2^{2}dy_{s}dy_{d}\right) dx_s dx_d = \frac{2501}{14400}-\frac{3}{80}\pi-\frac{2}{45}\ln\left(2\right)$$

For case $3$: $$J_3 = \frac{1}{4}\int_0^1 \int_{y_d}^{2-y_d} \left(\int_{y_{d}}^{\sqrt{y_{d}y_{s}}}\int_{0}^{2-\frac{y_{d}y_{s}}{x_{d}}}A_{3}^{2}dx_{s}dx_{d}+\int_{\sqrt{y_{d}y_{s}}}^{1}\int_{0}^{2-x_{d}}A_{3}^{2}dx_{s}dx_{d}-\int_{y_{d}}^{\sqrt{y_{d}\left(2-y_{s}\right)}}\int_{0}^{\frac{y_{d}\left(2-y_{s}\right)}{x_{d}}}A_{3}^{2}dx_{s}dx_{d}-\int_{\sqrt{y_{d}\left(2-y_{s}\right)}}^{1}\int_{0}^{x_{d}}A_{3}^{2}dx_{s}dx_{d}\right) dy_s dy_d = \frac{2501}{14400}-\frac{3}{80}\pi-\frac{2}{45}\ln\left(2\right)$$

For case $4$: $$J_4 = \frac{1}{4}\int_{0}^{1}\int_{1-\sqrt{1-x_{d}^{2}}}^{x_{d}}\int_{2+\frac{y_{d}^{2}-2y_{d}}{x_{d}}}^{2-x_{d}}\int_{\frac{2x_{d}-x_{d}x_{s}}{y_{d}}}^{2-y_{d}}A_{4}^{2}dy_{s}dx_{s}dy_{d}dx_{d}+\frac{1}{4}\int_{0}^{1}\int_{x_{d}}^{\sqrt{2x_{d}-x_{d}^{2}}}\int_{\frac{y_{d}^{2}}{x_{d}}}^{2-x_{d}}\int_{2-\frac{x_{d}x_{s}}{y_{d}}}^{2-y_{d}}A_{4}^{2}dy_{s}dx_{s}dy_{d}dx_{d} = -\frac{95}{288}+\frac{1}{12}\pi+\frac{1}{8}\ln(2)$$

Adding these up yields that $J = \frac{39}{800}+\frac{\pi}{120}-\frac{\ln(2)}{180}$. Multiplying by $8$ gives the final answer as $$\frac{39}{100}+\frac{\pi}{15}-\frac{2\ln(2)}{45} \approx 0.569$$, which is in the simulated range Arthur mentioned in the comments.