Coordinate transform of discrete data

let’s say I have a set of discrete data of dimension n by n by n, which may represent a scalar field $f(x,y,z), x,y,zin [-1/2,1/2]$ in Cartesian 3D space. To be specific,

n=64;
data=Table[i+j+k,{i,-n/2,n/2},{j,-n/2,n/2},{k,-n/2,n/2}];

so that $f(i/n,j/n,k/n)simeq data[[i,j,k]]$.
My question is, is there an efficient way to calculate the "spherically integrated" data? i.e., in the continuous form I would like to calculate
$$g(r)=int f(rsinthetacosphi,rsinthetasinphi,rcostheta)r^2sintheta dtheta dphi$$
In a less efficient way I would write something like

dataSph=Table[0,{l,0,2n}];
For[i=1,i<=n,i++,
  For[j=1,j<=n,j++,
    For[k=1,k<=n,k++,
      For[l=1,l<=2n,l++,
        If[l-1/2<=Sqrt[i^2+j^2+k^2]<l+1/2,dataSph[[l]]+=data[[i,j,k]]]
   ]]]]

or even just interpolate the data into a continuous version. Is there any better way to do the job?

idx = Tuples[Range@n, 3];
idxvals = Norm /@ N[idx];
dataSph = ConstantArray[0, 2 n + 1];
For[l = 1, l <= 2 n, l++,
  dataSph[[l]] = 
    Total[Extract[data, 
      Pick[idx, Unitize[Clip[idxvals, l + {-1/2, 1/2}, {0, 0}]], 1]]];
  ];