Your approach will not work. And I think it is easier to prove the general Pick's theorem for triangles first, and then apply it to the special case where the triangle has no lattice points inside or on the boundary!
To prove the theorem for a general triangle $T$, let $R$ be a smallest-area lattice rectangle that contains $T$. Clearly each side of $R$ must have some vertex of $T$, otherwise we can shift that side inwards to obtain a smaller lattice rectangle that contains $T$. Since $T$ has $3$ vertices, at least one vertex $A$ of $T$ must be on two sides of $R$. By symmetry (of the rectangle) we can assume that $A$ is at the bottom-left corner of $R$, and that the other vertices $B,C$ of $T$ are on the top and right side of $R$ respectively. The area of $T$ can be easily expressed as the area of $R$ minus the total area of the right-angled triangles formed by the sides of $R$ and $T$.
Thus the problem reduces to proving the theorem for a right-angled triangle $ABC$ with two axis-parallel sides $AB$ and $AC$. The area of $ABC$ is half the area of the rectangle $R$ with sides $AB$ and $AC$, and it is trivial to count the number of lattice points inside or on the boundary of $R$. I shall leave the rest of the details to you. =)
I found a cute geometric argument that $\mathbb{E}[X^2]=\frac{\pi^2}2$, and it seemed worth presenting here. (Many statements like "almost surely" and "if points are in general position" and "except for a set of measure $0$" have been omitted for ease of reading.)
Claim: If five points are randomly chosen on the sphere, their convex hull will be a triangle with probability $\frac5{16}$.
Proof: We start with a lemma.
Lemma: Given any three random points on the sphere, adding a random fourth point has a $50\%$ chance of producing a triangular convex hull.
Proof of lemma: Draw the three great circles connecting each pair of the three points, which subdivides the sphere into eight regions. Four of these regions have the property that a fourth point inside them will yield a triangular convex hull, and each such region is opposite a congruent region inside which a fourth point would not yield a triangular convex hull. So for any location $P$ of the fourth point, exactly one of $P$ and $-P$ would work. This completes the proof.
Now, back to the main theorem.
Given five points, consider ways of labeling one point $A$ and another $B$. Obviously, every five points have $20$ such labelings. Say that a collection of labeled points is good if, starting with the three unnamed points, we can add $A$ to get a triangular convex hull, and add $B$ onto the previous four to again produce a triangular convex hull.
For instance, in the following diagram, the left arrangement is good, but the right two are not (the second one has a quadrilateral convex hull when we add $A$, the third when we add $B$).
Now, fixing a labeling, what are the odds that starting with the unnamed points, adding $A$ gives a triangular hull and adding $B$ to those four points does as well? By our lemma, $\frac12\cdot\frac12=\frac14$. So every choice of five points produces an expected $\frac{20}4=5$ good labelings.
Conversely, given $5$ points whose convex hull is a triangle, how many good labelings does it have? Well, any choice of three points gives a valid starting triangle, and by the time we add the fifth point, we'll have a triangle once again, so the only way it could fail is if the first four points form a quadrilateral. How many ways can that happen?
It shouldn't be too hard to convince yourself that every $5$-point arrangement with a triangular convex hull looks like this, in terms of incidences and intersections (I've drawn lines as straight for convenience):
Inspecting this a bit, we can see that only one of the five points, when removed, leaves a quadrilateral convex hull (the top one, in the above diagram). So that has to be our choice of $B$ for things to fail, and then any of the other four could be $A$. So there are $20-4=16$ good labelings that produce a given configuration.
So, a random five-point arrangement produces $5$ good labelings on average, and a working arrangement produces $16$ good labelings. (Obviously, if a five-point arrangement doesn't have a triangular hull, it won't produce any good labelings.) So the fraction of five-point arrangements that work must be $\frac{5}{16}$.
Okay, proof complete. Now what?
Observe that the probability five points have a triangular convex hull is $10$ times the probability that they have a triangular convex hull and the three points on said triangle were the first three chosen (since any of the ten possible three-subsets could have just as easily been that triangle). So there is a $\frac1{32}$ chance that, if we pick three random points, the next two random points we choose will be inside that triangle.
But given a triangle, the odds that two random points lie inside it is just the square of the fraction of the area of the sphere that it takes up! So what we have shown is that $\mathbb{E}[\left(\frac{X}{4\pi}\right)^2]=\frac1{32}$, and hence $\mathbb{E}[X^2]=\frac{\pi^2}{2}$.
Note that this sort of argument can't extend to $6$ points, because the resulting probability of a triangular convex hull will be ${6\choose 3}\cdot \mathbb{E}[\left(\frac{X}{4\pi}\right)^3] = \frac{15}{32} - \frac{15\log(2)}{4 \pi^2}$, which is almost certainly irrational.
What breaks with $6$ points? I think the fundamental problem is that there is more than one point/line arrangement with positive probability, so the number of good labelings isn't constant. Consider the below two configurations:
![[diagram]](https://i.sstatic.net/CBOyo.png)
In the left case, removing any of the outer triangle's vertices will yield a quadrilateral convex hull, but in the right case, any vertex can be removed while preserving triangularity.
Best Answer
What makes this interesting is how to obtain the joint distribution of the points on the circle. Writing three uniformly random points in spherical coordinates and trying to obtain the corresponding angles on the circle leads to a convoluted mess, so let’s try to go the other way and calculate the Jacobian for three points parameterized using the circle they form.
Let the circle be parallel to the $x$-$y$ plane at polar angle $\theta$, with three points at azimuthal angles $\varphi_i$, so their coordinates are $\vec r_i=(\sin\theta\cos\varphi_i,\sin\theta\sin\varphi_i,\cos\theta)$. To move this configuration into general position, we can rotate the $z$ axis into an arbitrary direction, e.g. by applying rotations about the $y$ axis and then the $x$ axis, through $\psi$ and $\xi$, say, for a total of six parameters, as required.
Consider the infinitesimal changes along the six orthogonal unit vectors in the directions $\frac{\partial\vec r_i}{\partial\theta}$ and $\frac{\partial\vec r_i}{\partial\varphi_i}$ when we change these parameters. That’s six changes for six parameters, so we get a $6\times6$ Jacobian. Label the columns by $\varphi_1,\varphi_2,\varphi_3,\theta,\psi,\xi$ and the rows by $\frac{\partial\vec r_1}{\partial\varphi_1},\frac{\partial\vec r_2}{\partial\varphi_2},\frac{\partial\vec r_3}{\partial\varphi_3},\frac{\partial\vec r_1}{\partial\theta},\frac{\partial\vec r_2}{\partial\theta},\frac{\partial\vec r_3}{\partial\theta}$.
Since we only need the dependence of the joint distribution on the $\varphi_i$, I’m going to drop all factors that don’t depend on them. Then the upper left $3\times3$ minor is the identity, and the one below it is zero, since changing $\varphi_i$ only causes a change along the corresponding vector $\frac{\partial\vec r_i}{\partial\varphi_i}$ and the magnitude of that change is the same for all points (namely $\sin\theta$) and doesn’t depend on the $\varphi_i$.
Thus, to get the determinant of the Jacobian, we can use Laplace expansion along the first three columns, which leaves us with the determinant of the lower right $3\times3$ minor.
The first column of that minor, which contains the changes with $\theta$ along $\frac{\partial r_i}{\partial\theta}$, has all entries $1$. To obtain the other two columns, consider the effect of the infinitesimal generators $\pmatrix{0&0&\epsilon\\0&0&0\\-\epsilon&0&0}$ and $\pmatrix{0&0&0\\0&0&\epsilon\\0&-\epsilon&0}$ of rotations about the $y$ and $x$ axis, respectively, on the $\vec r_i$:
$$ \pmatrix{0&0&\epsilon\\0&0&0\\-\epsilon&0&0}\pmatrix{\sin\theta\cos\varphi_i\\\sin\theta\sin\varphi_i\\\cos\theta}=\epsilon\pmatrix{\cos\theta\\0\\-\sin\theta\cos\varphi_i}\;, \\[18pt] \pmatrix{0&0&0\\0&0&\epsilon\\0&-\epsilon&0}\pmatrix{\sin\theta\cos\varphi_i\\\sin\theta\sin\varphi_i\\\cos\theta}=\epsilon\pmatrix{0\\\cos\theta\\-\sin\theta\sin\varphi_i}\;. $$
Taking the scalar product with the unit vector $\frac{\partial r_i}{\partial\theta}=(\cos\theta\cos\varphi_i,\cos\theta\sin\varphi_i,-\sin\theta)$ yields $\cos\varphi_i$ and $\sin\varphi_i$, respectively, so the lower right minor is
$$\begin{vmatrix}1&\cos\varphi_1&\sin\varphi_1\\1&\cos\varphi_2&\sin\varphi_2\\1&\cos\varphi_3&\sin\varphi_3\end{vmatrix}\;.$$
Interestingly, this happens to be proportional to the area of the triangle formed by the points, whose ratio to the area of the circle is, as Blitzer noted in a comment, given by
$$\frac2\pi\left|\sin\frac\alpha2\sin\frac\beta2\sin\frac\gamma2\right|=\frac1{2\pi}|\sin \alpha+\sin\beta+\sin\gamma|\;,$$
where $\alpha=\varphi_2-\varphi_1$, $\beta=\varphi_3-\varphi_2$ and $\gamma=\varphi_1-\varphi_3$.
Thus, the expected value of this ratio is
$$ \mathsf E\left[\frac\triangle\bigcirc\right] = \frac1{2\pi}\cdot\frac{\iiint(\sin\alpha+\sin\beta+\sin\gamma)^2\mathrm d\varphi_1\mathrm d\varphi_2\mathrm d\varphi_3}{\iiint|\sin\alpha+\sin\beta+\sin\gamma|\mathrm d\varphi_1\mathrm d\varphi_2\mathrm d\varphi_3}\;. $$
In the numerator, the integral over the mixed terms vanishes, leaving
$$ \iiint\left(\sin^2\alpha+\sin^2\beta+\sin^2\gamma\right)\mathrm d\varphi_1\mathrm d\varphi_2\mathrm d\varphi_3=\frac32(2\pi)^3\;. $$
In the denominator, by symmetry we can restrict to the $\varphi_i$ ordered counterclockwise (and thus $\alpha$, $\beta$, $\gamma$ positive) to get rid of the absolute value and obtain
\begin{eqnarray*} \iiint|\sin\alpha+\sin\beta+\sin\gamma|\,\mathrm d\varphi_1\mathrm d\varphi_2\mathrm d\varphi_3 &=& 2\cdot3\cdot2\pi\int_0^{2\pi}\int_0^{2\pi-\alpha}\sin\alpha\,\mathrm d\beta\,\mathrm d\alpha \\ &=& 2\cdot3\cdot2\pi\int_0^{2\pi}(2\pi-\alpha)\sin\alpha\,\mathrm d\alpha \\ &=& 2\cdot3\cdot(2\pi)^2\;. \end{eqnarray*}
Thus, with
$$ \frac1{2\pi}\cdot\frac{\frac32(2\pi)^3}{2\cdot3\cdot(2\pi)^2}=\frac14\;, $$
your numerical findings are confirmed.