Interpolating values between points constrained to a 3D surface

I have run into a bit of a problem using scatteredInterpolate between measurements on a 3D surface and am hoping somebody here can provide some suggestions.

My data comes from computed tomography, where I have fit a surface around features in the reconstructed volume, and then taken sparse measurements of features on the surface. So each measurement is associated with 3D coordinate - F(X,Y,Z) - which corresponds to a point on the surface.

I now want to interpolate the values from the measured points across the entire surface. The measurements are made on a meshed surface, although the surface I wish to interpolate over has been converted to a voxelized surface. I originally thought this would work well with scatteredInterpolate, but I have found that because this method based on a triangulation, most values on the surface end up outside of the convex hull of the measurement points and need to be extrapolated. Nearest neighbour interpolation/extrapolation produces the expected (coarse) results, however the linear method produces poor results, as there are large local discontinuities between values that result from interpolation compared to those from extrapolation.

The surface is open but slightly curved with some overlap/wraparound (i.e. it could not be described as Z(X,Y) across its full X,Y range). So I suspect that the discontinuities result from nearly all the points on the surface being either close to the inner border of, or outside of, the convex hull of the original measurement points (see https://se.mathworks.com/help/matlab/math/interpolating-scattered-data.html#bspy6dx-1)

I now realize that it would be better to do the interpolation/extrapolation such that it was constrained to the voxelized surface, rather than essentially treating the measurements as unconstrained points in a 3D volume with scatteredInterpolate). I've been searching for a method to solve this problem for a little while (so I don't think I have overlooked any obvious Matlab functionality), but I don't think there is any inbuilt functionality in matlab for this. Maybe one way would be to calculate geodesic distance transforms on the surface from each of the original measurement points, and then use those in some kind of inverse distance weighting calculation. Are there any other suggestions or anything I have overlooked?