Oh, I should think a bit more before asking as I've now answered my own question. It is simply a case of requiring 2 scattered interpolations instead of 1 grid interpolation.
By treating the 3D volume as a series of 2D slices (each 'page' of the matrix is a new 2D slice), we can first interpolate to find the correct xyz co-ordinates. The plane of the slice dictates the values to be interpolated.
E.G I wish to make a slice in the xy plane. Therefore I require all (unique) y-coordinates, and the x-coord of the plane. By gridding these with meshgrid, I can then interpolate Z of each 2D slice by using TriScatteredInterp (on the original x-y-z grids with the NaNs removed). This gives us the meshgrid results, xi and yi, and zi, the interpolated values. Remember at the end to replicate the n*1 xi-yi grids to the number of 2D slices.
Almost there - I now have a slice with the correct XYZ coordinates. The final step is to interpolate the colour data with another call to TriScatteredInterp, this time a 3D interpolation from the original X-Y-Z-D grids.
A quick call to squeeze to reduce the dimensions to 2 and I am left with the xi,yi,zi and di matrices required by surface.
