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}\ .$