Your understandig is not appropriate. You use cartesian coordinate and wonder why the result is not spherical. You remain on the cartesion grid and only assign each cube a spherical calculated value of a combined gray level with opacity.
The problem is not changed by
res = 10;
coords = Table[{Sqrt[i^2 + j^2 + k^2], ArcTan[k, Sqrt[i^2 + j^2]],
ArcTan[i, j]}, {i, -res, res}, {j, -res, res}, {k, -res, res}];
dists = Map[N@Norm[# - .5] &, coords[[All, All, All, 1]], {3}];
opa = 1/(dists + 0.3);
opa /= Max@opa;
ImageAdjust@Image3D[opa, ColorFunction -> "GrayOpacity"]

Some small turn and You see the cubes much clearer:
.
res is not the right parameter to get the desired shading!
The very same is in Your solution at work.
The documentation page for DelaunayMesh shows an example how to avoid cubic volume grids in Mathematica:
pts = RandomReal[{-1, 1}, {25, 3}];
Rsc = DelaunayMesh[pts];
HighlightMesh[Rsc, {Style[0,
Directive[PointSize[Medium], Black]], Style[2, Opacity[0.1]]}]

Your procedure path may be Delauney mesh spheres stacked in each other and assign the volumes of this mesh elements gray levels and opactiy. Then increase res again.
This guide list Solid Geometry what is in Mathematica. Cone
has some samples how this looks on the basis of Cone
and Circle
in color.
Sphere show some example: Place Sphere instead of Cube. The main problem is that Mathematica visualizes surfaces instead of volumes by primitives. A volume is represented by the boundary and that is in 3D a surface. The professional term is tessalation. Boundaries are reprensented at most by Polygons
. The shadings are achieved if desires by Your res parameter alike parameters. This relies on a scene graph like structure and that is it.
But as shown in Volume
Mathematica has both MeshRegion and MeshBoundaryRegion. The built-in MeshPrimitives should suit Your question real hidden intent, but it does not of-the-shelf. Sphere meshing is problematic in Mathematica, surface-mesh-of-a-hemisphere. Have a look at this: draw-a-triangulated-sphere This answer how-to-discretize-a-sphere suggests
Needs["PolyhedronOperations`"]
G = Geodesate[PolyhedronData["Icosahedron", "Faces"], 5] // N;
Graphics3D@G
as discrete closed accurate approximation of the Sphere. It is suitable for Your intent. Or more random: voronoi grid on a sphere.
Or use my full example and increase $res$ and optimize the left over parameters. So chose a different distribution with a sharper upper drop. For example Fermi or piecewise, step function and alike.
This is a counterexample:
Manipulate[
Graphics3D[{Opacity[a], Black, Specularity[White, 5], Sphere[]},
Lighting -> "Neutral"], {a, 0, 1, 0.1}]

This shows up clearly the concept of Mathematica Graphics3D and color and lighting and opacity. And it works. The Sphere is special and somewhat cheating the critical human eye. That is all. No second level or higher. The documentation of Specularity
clearly state that Specularity
does not work with flat 3D surface. You are making use of flat 3D surfaces.
Another example:
Manipulate[
Graphics3D[{Opacity[a], GrayLevel[.25], Specularity[White, 10],
Sphere[]}, Lighting -> "Neutral"], {a, 0, 1, 0.1}]

Manipulate[
Graphics3D[{Opacity[a], Glow[Black], White, Sphere[]},
Lighting -> "Neutral"], {a, 0, 1, 0.1}]

DensityPlot[Exp[-.15/(x^2 + y^2)], {x, -1, 1}, {y, -1, 1},
ColorFunction -> GrayLevel, Frame -> None]

Table[DensityPlot3D[
Exp[-.5/(x^2 + y^2 + z^2)], {x, y, z} \[Element] Ball[],
ColorFunction -> GrayLevel,
OpacityFunction -> Function[f, (1 - f)^2], Boxed -> False,
PlotPoints -> pp], {pp, {50, 75, 100}}]

The more plotpoints the smoother.
Or make use of the DensityPlot
as a texture on the Sphere
.


Make use of the image dimensions in the parameter set for ImageSize.
The great trick with Raster
will help it further the last step.