How to get the estimator? They introduce a new correlation. For $\pi\in \Pi(\mu,\nu)$ the set of coupling of two probability measures $\mu$ and $\nu$ on a Polish space $(X,d)$. The author introduces a plugin estimator.
Question: How to understand $\hat{\pi}$?

For the estimator of the denominator, I use the same notation as the answer.
\begin{equation*}
    D(\hat{\pi})=\sum_{y,z}d(y,z)\hat{\pi}_2(\{y\})\hat{\pi}_2(\{z\})
    =\sum_{y,z}d(y,z)\frac{1}{N^2}\sum 1[\phi(Y_n)\in\{y\}]\sum 1[\phi(Y_m)\in\{z\}]
\end{equation*}
Expand the product of these two summation,
$$
=\frac{1}{N^2}\sum_{y,z}d(y,z)(\sum_{i=1}^N 1[\phi(Y_i)\in\{y\}]1[\phi(Y_i)\in\{z\}]+\sum_{i\neq j}1[\phi(Y_i)\in\{y\}]1[\phi(Y_j)\in\{z\}] )
$$
 A: $\newcommand{\de}{\delta}$We have
\begin{equation*}
    W(\pi):=\frac{N(\pi)}{D(\pi)}, \tag{1}\label{1}
\end{equation*}
where
\begin{equation*}
    N(\pi):=\int W_1(\pi_{x_1},\nu)\mu(dx_1)=\int W_1(\pi_{x_1},\pi_2)\pi_1(dx_1),  
\end{equation*}
\begin{equation*}
        D(\pi):=\int d(y,z)\nu(dy)\nu(dz)=\int d(y,z)\pi_2(dy)\pi_2(dz), 
\end{equation*}
$\pi_{x_1}$ is the conditional distribution of $\eta$ given $\xi=x_1$ if the joint distribution of a pair $(\xi,\eta)$ of random variables is the measure $\pi$ with marginals $\mu=\pi_1$ and $\nu=\pi_2$.
The standard notation $\de_a$ is for the Dirac delta measure supported on the singleton set $\{a\}$. So,
\begin{equation*}
    \hat{\pi}:=\frac1N\,\sum_{i=1}^N \de_{(\phi(X_i),\,\phi(Y_i))}
\end{equation*}
is the (random) probability measure such that
\begin{equation*}
    \hat{\pi}(A)=\frac1N\,\sum_{i=1}^N 1((\phi(X_i),\,\phi(Y_i))\in A) \\ 
    =\frac1N\,\sum_{i=1}^N |\{i\in\{1,\dots,N\}\colon (\phi(X_i),\,\phi(Y_i))\in A\}|
\end{equation*}
for any Borel set $A\in[0,1]^2$. Here $1(\mathcal A)$ is the indicator of an assertion $\mathcal A$, so that $1(\mathcal A)=1$ if $\mathcal A$ is true and $1(\mathcal A)=0$ if $\mathcal A$ is false.
Then for the first marginal $\hat\pi_1$ of $\hat\pi$ and each $x_1\in X=[0,1]$ we have
\begin{equation*}
    \hat\pi_1(\{x_1\})=\frac1N\,\sum_{i=1}^N 1(\phi(X_i)\in \{x_1\}) \\ 
    =\frac1N\,\sum_{i=1}^N 1(X_i\in\phi^{-1}(\{x_1\}))
    =\frac1N\,\sum_{i=1}^N 1(X_i\in G)
\end{equation*}
for $G=\phi^{-1}(\{x_1\})$,
and for the conditional distribution $\hat\pi_{x_1}$ of $\eta$ given $\xi=x_1$ if the joint distribution of a pair $(\xi,\eta)$ is $\hat\pi$ we have
\begin{equation*}
\hat\pi_{x_1}(B)=\frac{\hat\pi(\{x_1\}\times B)}{\hat\pi(\{x_1\}\times X)}, \tag{2}\label{2}
\end{equation*}
where $B$ is any Borel subset of $X=[0,1]$; if the denominator of the ratio in \eqref{2} is $0$, let $\hat\pi_{x_1}$ be an arbitrary probability measure. Next,
\begin{equation*}
    \hat\pi(\{x_1\}\times B)=\frac1N\,\sum_{i=1}^N 1(\phi(X_i)=x_1,\,\phi(Y_i)\in B) \\ 
    =\frac1N\,\sum_{i=1}^N 1(X_i\in\phi^{-1}(\{x_1\}),\,\phi(Y_i)\in B). 
\end{equation*}
So, for $G=\phi^{-1}(\{x_1\})$,
\begin{equation*}
\hat\pi_{x_1}(B)=\hat\pi_G(B):=\frac{\frac1N\,\sum_{i=1}^N 1(X_i\in G,\,\phi(Y_i)\in B)}
    {\frac1N\,\sum_{i=1}^N 1(X_i\in G)}.  \tag{3}\label{3}
\end{equation*}
So,
\begin{equation*} 
\begin{aligned}
    N(\hat\pi)&=\int W_1(\hat\pi_{x_1},\hat\pi_2)\hat\pi_1(dx_1) \\ 
    &=\sum_{x_1} W_1(\hat\pi_{x_1},\hat\pi_2)\hat\pi_1(\{x_1\}) \\ 
&   =\sum_{G\in \Phi} W_1(\hat\pi_G,\hat\pi_2)\frac1N\,\sum_{i=1}^N 1(X_i\in G) \\ 
&   =\sum_{G\in \Phi} \frac{1}{N}|i\in\{1,\dots,N\}\ \text{s.t.}\ X_i\in G|\, 
W_1(\hat{\pi}_G,\hat{\pi}_2).  
\end{aligned}
\end{equation*}
Also, if $d$ is the standard metric over $[0,1]$, then, as you have it in your "new edition",
\begin{equation}
\begin{aligned}
    D(\hat\pi)&=\int|y-z|\hat\pi_2(dz)\hat\pi_2(dz) \\ 
    &=\sum_{y,z}|y-z|\hat\pi_2(\{y\})\hat\pi_2(\{z\}) \\ 
    &=\sum_{y,z}|y-z|\frac{1}{N^2}\sum_{n,m=1}^N 1(\phi(Y_n)=y,\phi(Y_m)=z) \\ 
    &=\frac{1}{N^2}\sum_{n,m=1}^N \sum_{y,z}|y-z|1(\phi(Y_n)=y,\phi(Y_m)=z) \\ 
    &=\frac{1}{N^2}\sum_{n,m=1}^N |\phi(Y_n)-\phi(Y_m)|. 
\end{aligned}
\end{equation}
So,
\begin{equation*}
     W(\hat\pi):=\frac{N(\hat\pi)}{D(\hat\pi)}
    =\frac{\sum_{G\in \Phi} \frac{1}{N}|i\in\{1,\dots,N\}\ \text{s.t.}\ X_i\in G|\, W(\hat{\pi}_G, \hat{\pi}_2)}{\frac{1}{N^2}\sum_{n,m=1}^N|\phi(Y_n)-\phi(Y_m)|} 
\end{equation*}
is a plug-in estimator of $W(\pi)$.
