Here is my computation. First of all, I had no idea whatsoever how to integrate implicit trigonometric functions, so I decided to switch to the integration with respect to the sides $x,y,z$ of the triangle $XYZ$. Recall that the conditions for the domain is $xz\le y^2+xy$, $x>y$ plus the triangle inequality $z>x-y$. This triangle inequality is compatible with the first condition if and only if $x< (\sqrt 2+1)y$. So the admissible domain is given by
$$
y< x< (\sqrt 2+1)y,\qquad x-y<z<y+\frac {y^2} x\,.
$$
Now we have to integrate over the points on the circle, i.e., for the case when the circumradius is $1$ and with respect to the angular variables. We relax the circumradius to $r>0$, place $Z$ on the positive semi-axis and the circle center at the origin (so $Z=r$), and take $Y=re^{2i\varphi}, X=re^{2i\psi}$, so
$$
x=2r\sin\varphi,\qquad y=2r\sin\psi,\qquad z=2r\sin(\varphi-\psi)
$$
(the last one up to a sign) where $\varphi,\psi$ run independently over $[0,\pi]$. The condition $x>y$ reduces the $d\varphi\,d\psi$ measure to $\frac{\pi^2}2$ and each admissible triple $x,y,z$ is realized twice (up to the flip about the horizontal axes). We can use any measure on $(0,+\infty)$ we want when integrating over $r$, so we will choose $\frac{dr}{r}$.
Now it is time to compute the Jacobian $J$ of the mapping $(r,\varphi,\psi)\mapsto (x,y,z)$. We have
$$
J=8\left|\begin{matrix}
\sin\varphi & r\cos\varphi & 0
\\
\sin\psi & 0 & r\cos\psi
\\
\sin(\varphi-\psi) &r\cos(\varphi-\psi) & -r\cos(\varphi-\psi)
\end{matrix}\right|
\\
=8r^2[\cos\varphi\cos\psi\sin(\varphi-\psi)-\sin\varphi\cos\psi\cos(\varphi-\psi)+\cos\varphi\sin\psi\cos(\varphi-\psi)]
\\
=8r^2[\cos\varphi\cos\psi\sin(\varphi-\psi)-
\sin(\varphi-\psi)\cos(\varphi-\psi)]
\\
=8r^2\sin(\varphi-\psi)\sin\varphi\sin\psi=\frac{xyz}r\,,
$$
i.e.
$$
\frac{dr}r\,d\varphi\,d\psi=\frac{dx\,dy\,dz}{xyz}\,,
$$
which is pretty neat.
Now let $\Omega$ be the set of all good triangles and $\omega$ be the set of good triangles with $r=1$. Then the back of the envelope computation is the following:
$$
\left[\int_0^\infty\frac{dr}{r}\right]\int_\omega d\varphi\,d\psi=
2\int_{\Omega}\frac{dr}{r}\,d\varphi\,d\psi
=2\int_{\Omega}\frac{dx\,dy\,dz}{xyz}
\\
=
2\int_0^\infty\frac{dy}{y}\left[\int_{y}^{(\sqrt 2+1)y}\frac{dx}x\left(\int_{x-y}^{y+\frac{y^2}x}\frac{dz}z\right)\right]
$$
Substituting $x=ty$, $t\in(1,\sqrt 2+1)$, we evaluate the inner integrals to
$$
\int_1^{\sqrt 2+1}\log\frac{t+1}{t(t-1)}\frac{dt}t\,,
$$
which is independent of $y$, so we finally get
$$2\left[\int_0^\infty\frac{dy}y\right]
\int_1^{\sqrt 2+1}\log\frac{t+1}{t(t-1)}\frac{dt}t\,,
$$
Now it remains to cancel $\int_0^\infty\frac{ds}s$ to get the claimed identity.
The little (but somewhat irritating) problem is that the integral $\int_0^\infty \frac{ds}s$ diverges. However, as Landau (the physicist) used to say, "A chicken is not a bird and a logarithm is not infinity". So we will fix that with a little bit more elaborate argument. ${}{}$
We shall impose one extra condition on the admissible triangles, namely, we will demand that $y$ is comparable to $r$. We always trivially have $r\ge \frac y2$, so we'll fix a big $A>0$ and demand that $r<Ay$. Notice that it is a condition on the shape only, not on the size. Denoting by $\omega_A$ the set of admissible triangles with circumradius $1$ and by $\Omega_A$ the set of all admissible triangles and choosing a huge finite $R>0$, we can write, as before,
$$
\left[\int_1^R\frac{dr}r\right]\int_{\omega_A}d\varphi\,d\psi=2\int_{\Omega_A\cap\{1<r<R\}}\frac{dx\,dy\,dz}{xyz}\,.
$$
Now notice that
$$
\Omega_A\cap\{2<y<\frac RA\}\subset \Omega_A\cap\{1<r<R\}\subset\Omega_A\cap\{\frac 1A<y<2R\}\,,
$$
so instead of one identity, we write two inequalities and conclude that $\log R\int_{\omega_A}d\varphi\,d\psi$ is squeezed between $\log\frac{R}{2A}$ times $2\int_1^{\sqrt 2+1}\left(\int_{E(t,A)}\frac{ds}{s}\right)\frac{dt}t$ and $\log(2AR)$ times the same quantity, where we put $x=ty, z=sy$ and defined $E(t,A)$ to be the set of $s$ for which the triangle with the sides $1,t,s$ is admissible.
Now, dividing by $\log R$, letting $R\to+\infty$ and using the squeeze theorem, we conclude that
$$
\int_{\omega_A}d\varphi\,d\psi=
2\int_1^{\sqrt 2+1}\left(\int_{E(t,A)}\frac{ds}{s}\right)\frac{dt}t
$$
for every fixed $A$. But when $A\to+\infty$, the set $\omega_A$ expands to $\omega$ and for each $t$, the set $E(t,A)$ expands to $[t-1,1+\frac 1t]$, so the monotone convergence theorem finishes the story.
That's it. I tried to find a simple change of variables that would avoid using non-trivial identities for $\mathrm{Li}_2$ to get the final $\frac{\pi^2}8$ answer, but failed so far. Still, we have the result, so, I hope, somebody will eventually provide a simpler proof. The value $\frac 12$ of the probability is too nice to be there without a clear reason. :-)
Best Answer
The conjecture is true, and it can be verified with fedja's method developed for your earlier question. Indeed, this method shows that $$P(\text{$a^p+b^p\ge c^p$ and $c>b>a$})=\left(\frac{\pi^2}{3!}\right)^{-1}2\int_1^{\infty}\left(\int_{E(t)}\frac{ds}{s}\right)\frac{dt}t,$$ where $$E(t):=\{s\in\mathbb{R}:\text{$1+t^p\geq s^p$ and $s>t$}\}.$$ Clearly, $$\int_{E(t)}\frac{ds}{s}=\log((1+t^p)^{1/p})-\log(t)=\frac{\log(1+t^{-p})}{p},$$ hence $$P(\text{$a^p+b^p\ge c^p$ and $c>b>a$})=\frac{12}{\pi^2}\int_1^{\infty}\frac{\log(1+t^{-p})}{pt}\,dt.$$ On the other hand, $$\log(1+t^{-p})=\sum_{k=1}^\infty\frac{(-1)^{k-1}}{k t^{kp}},\qquad t>1,\tag{1}$$ therefore $$P(\text{$a^p+b^p\ge c^p$ and $c>b>a$})=\frac{12}{\pi^2}\sum_{k=1}^\infty\frac{(-1)^{k-1}}{k^2 p^2}=\frac{1}{p^2}.$$
Added. In the last step, we used that $$\sum_{k=1}^\infty\frac{(-1)^{k-1}}{k^2}=\frac{\pi^2}{12},\tag{2}$$ which is equivalent to the solution of the Basel problem. On the other hand, for $p=1$, the probability in question equals $1$ by the triangle inequality, hence in that case the above calculation yields a proof of $(2)$. This was observed by the OP, and in a slightly different context by Nilotpal Sinha. In fact the triangle inequality allows us to eliminate both $(1)$ and $(2)$ from the above calculation. Namely, by a change of variable we can see that $$\int_1^{\infty}\log(1+t^{-p})\frac{dt}{t}$$ is proportional to $p^{-1}$, hence the probability sought is proportional to $p^{-2}$.