Given the cuboid $C:=[x_0,x_1]\times[y_0,y_1]\times[z_0,z_1]$ with vertices $v_{ijk}=(x_i,y_j,z_k)$ and assigned values $u_{ijk}$ there is a unique "trilinear function" $T$ with $T(v_{ijk})=u_{ijk}$.
Enclose your sphere $S$ with center $p$ and radius $r>0$ in a cube $C'\subset C$ of sidelength $2r$. $C'$ has vertices $v'_{ijk}$ and assigned values $u'_{ijk}:=T(v_{ijk}')$.
You cannot expect that $T$ has a unique maximum on $S$. If, e.g., $T$ is $=1$ on four tetrahedrally arranged vertices $v_{ijk}'$ of $C'$ and $=-1$ on the other four vertices of $C'$ then $T$ has four local maxima, four local minima and six saddle-points on $S$. Therefore, if you need the exact value of the maximum, you have to do a lot of work in such a case.
The following is an easy way out and should produce a good approximation, in particular if $S$ is rather small compared to $C$.
I propose to do a linear regression on the data concerning $C'$, where I omit the $'$ in the sequel. This means the following: There is a unique linear function
$$f(x,y,z):=a x+ by + c z + d$$
which makes $\sum_{ijk} \bigl(f(v_{ijk})- u_{ijk}\bigr)^2$ minimal. The unknown coefficients $a$, $b$, $c$, $d$ depend linearly on the data vector $u:=(u_{ijk})$. Symmetry considerations suggest that
$$a={\sum_{jk} (u_{1jk}-u_{0jk})\over 4(x_1-x_0)}\ ,$$
and similarly for $b$ and $c$ ; finally $d$ should be chosen such that at the center $p$ of $C$ one has $f(p)={1\over8}\sum_{ijk} u_{ijk}$.
The maximum $M$ of $f$ on the sphere $S$ is then given by
$$M=f(p)+r\sqrt{a^2+b^2+c^2}\ .$$