Spherical interpolation in Matlab from Cartesian Grid

Question

I have a 3D image. I need to take spherical wave decomposition of this image hence I need to convert my 3d image in the cartesian grid to spherical coordinates.

If below is too long to read, gist of what I want is that I want to a mapping that will interpolate the spherical coordinates such that, especially around my origin I don't have a lot of missing information.

First, I do the spherical conversion.

    [ydim,xdim,zdim] = size(myimg);

    [x,y,z] = meshgrid(1:xdim, 1:ydim, 1:zdim);
    x = x - median(x(:)); y = y - median(y(:)); z = z - median(z(:)); 
    [phis, thetas, rs] = cart2sph(x,y,z);

From here I am stuck. How can I use my phis and thetas and r's to do interpolation along an arc (I just assume it would be an arc since it's a sphere).

I actually searched about this a bit and brought together this code for interpolation but I can't verify it works. Mainly because most of it is copied and modified from a similar problem.

function [theta0,phi0,rho0] = my_interp(X, Y, Z, nTheta0, nPhi0)
% forget about spherical, let's just interpolate in cartesian then convert 
% to spherical. 
[theta, phi, V] = cart2sph(X, Y, Z);

% X,Y,Z are meshgrid output from above code snippet.
P = [2 1 3];
X = permute(X, P);
Y = permute(Y, P);
Z = permute(Z, P);
V = permute(V, P);

% create a cartesian interpolant, and we'll use it in spherical.
F = griddedInterpolant(X,Y,Z,V);

% prepare grid for meshing (we'll mesh xyz data, not theta,phi)
theta0 = linspace(-pi, pi, nTheta0);
phi0   = linspace(-pi/2, pi/2, nPhi0);
[theta0, phi0] = meshgrid(theta0, phi0);

[x_,y_,z_] = sph2cart(theta0, phi0, 1 ); % !! here is why I get confused. my 
% radius on the sphere is not just 1, I have changing radius. Am I really     
% missing the key insight here? On the other hand, I don't know how to      
% account for all the r's on the image anyways so I can't change this.
rho0 = F(x_,y_,z_);

theta0 = repmat(theta0, 1, 1, size(X, 3));
phi0 = repmat(phi0, 1, 1, size(X, 3));
rho0 = repmat(rho0, 1, 1, size(X, 3));
end

Answer 1

Instead of creating a grid of cartedian coordinates and converting those to spherical, you want to create a grid of spherical coordinates, convert those to cartesian, then use interp3 to find the image values at those coordinates.

Spherical coordinate grid could be something like meshgrid(0:maxR,0:phiStep:2*pi,0:thetaStep:pi).