# Spherical density plot of data set

I have a somewhat large number of 3-dimensional unit vectors ( ~ 100 000 unit vectors), and I am trying to visualize their distribution. First thing I tried was ListPlot3D, but it doesn't look too good.

Now, I would like to make a density plot on the surface of the unit sphere using this data. For example, the sphere would be bluish in the regions with little or no points, and reddish in regions with many points (you can use any colors you want).

Can I do this with Mathematica?

-
–  becko Oct 4 '12 at 23:24
What is the structure of your data? Can you give the code of what you tried, even with a smaller toy-set of vectors? It's better if you upload original data somewhere and give the link, or give the code to how you generated them. –  Vitaliy Kaurov Oct 4 '12 at 23:31
–  s0rce Oct 5 '12 at 2:48
Sorry but I could not ask question because I do not have enough reputation. I am new to mathematica and I am just trying to understand above posted code. I would be glad if you can help me on this. Thanks –  user4992 Dec 13 '12 at 16:23

I've decided to write a new answer instead of editing my previous one, since the following discussion will be rather long, and the method I am about to discuss is similar, but not quite the same.

The following approach makes use of formulae from this paper to determine density estimates for spherical data. Due to being specially adapted for spherical data, it does not suffer from boundary effects as with my previous answer. From here, one could either display the "smooth spherical histogram" as a texture (like in my previous answer), or as a surface (in complete analogy with SmoothHistogram3D[]).

Here's a demonstration:

(* generate a random unit vector following the Dimroth-Watson distribution, girdle case *)
dimrothWatsonRandom[μ_?VectorQ, κ_ /; NumericQ[κ] && Negative[κ]] := Module[{c, d, u, v, w},
c = Sqrt[-κ]; d = ArcTan[c];
While[
{u, v} = RandomReal[1, 2]; w = Tan[d u]/c;
v > (1 - κ w^2) Exp[κ w^2]];
RotationTransform[{{0, 0, 1}, Normalize[μ]}][Append[
Sqrt[1 - w^2] Normalize[RandomVariate[NormalDistribution[], 2]], RandomChoice[{-1, 1}] w]]]

BlockRandom[SeedRandom[2031, Method -> "MersenneTwister"]; (* for reproducibility *)
vecs = Table[dimrothWatsonRandom[Normalize[{-1, 0, 1}], -10], {150}]];

(* smoothing parameter, automatically determined by maximizing "pseudo-log likelihood" *)
cc = With[{n = Length[vecs]},
NArgMax[Sum[Log[
Total[c Csch[c] Exp[c Delete[vecs, k].Extract[vecs, k]]/(4 π (n - 1))]],
{k, n}], c]];

hist = With[{n = Length[vecs]}, Image[DensityPlot[
Total[cc Csch[cc] Exp[cc vecs.{Cos[θ] Sin[φ], Sin[θ] Sin[φ], Cos[φ]}]/(4 π n)],
{θ, -π, π}, {φ, 0, π}, AspectRatio -> Automatic,
ColorFunction -> "ThermometerColors", Frame -> False, Mesh -> True,
MeshFunctions -> {#3 &}, ImagePadding -> None, PerformanceGoal -> "Quality",
PlotPoints -> 75, PlotRange -> All, PlotRangePadding -> None],
ImageResolution -> 256]];

(* spherical smooth density histogram as texture *)
ParametricPlot3D[{Cos[θ] Sin[φ], Sin[θ] Sin[φ], Cos[φ]}, {θ, -π, π}, {φ, 0, π},
Lighting -> "Neutral", Mesh -> None,
PlotStyle -> Texture[hist], TextureCoordinateFunction -> ({#4, #5} &)]


(* spherical smooth histogram *)
With[{h = 2/3 (* scaling parameter *), n = Length[vecs]},
SphericalPlot3D[
1 + h Total[cc Csch[cc] Exp[cc vecs.{Cos[θ] Sin[φ], Sin[θ] Sin[φ], Cos[φ]}]/(4 π n)],
{φ, 0, π}, {θ, -π, π}, BoundaryStyle -> None, MeshFunctions -> {#6 &},
MeshStyle -> AbsoluteThickness[1], PlotPoints -> 85]]


I should mention that anyone who has to deal with spherical data should take a look at the book Statistical Analysis of Spherical Data; the book mentions this approach, as well as a number of other methods for visualizing spherical data.

-
Very nice, especially the spherical smooth histogram. –  VLC Nov 3 '12 at 8:44
Wow, cool stuff, - and +1 for reading papers ;) –  Vitaliy Kaurov Nov 4 '12 at 1:07

3D figures always look nice, but sometimes, they are limited in the way they can present data. Let's suppose that we have some data like those generated in the answer that J. M. gave, but they are dense around two points that are on opposite poles. In this case, a 3D figure would always hide half of the data. To show all the data we can transform the data using the Lambert azimuthal equal-area projection and then plot the data in 2D.

Let's generate our data using J. M.'s function. Here, I set one set of points to cluster around $(1,0,0)$ and the other one on the opposite side around $(-1,0,0)$:

BlockRandom[SeedRandom[42, Method -> "MersenneTwister"];
vecs1 = Table[vonMisesFisherRandom[{1, 0, 0}, 3], {5000}]];

BlockRandom[SeedRandom[42, Method -> "MersenneTwister"];
vecs2 = Table[vonMisesFisherRandom[{-1, 0, 0}, 3], {5000}]];

data = Join[vecs1, vecs2];


We create a function for data transformation, transform the data and plot the projected data in 2D:

xyzToLambert[xyzCoords_List] := Module[{xx, yy},
xx = Sqrt[2/(1 + 10^-10 - xyzCoords[[3]])] xyzCoords[[1]];
yy = Sqrt[2/(1 + 10^-10 - xyzCoords[[3]])] xyzCoords[[2]];
{xx, yy}
]

projdata = Table[xyzToLambert[data[[i]]], {i, Length[data]}];

ListPlot[projdata, AspectRatio -> 1, PlotRange -> {{-3, 3}, {-3, 3}}]


This already looks nice, but with many more data points, it is unlikely that we could see anything at all. We can thus proceed and make some other plots:

opts = {ColorFunction -> Function[{height}, ColorData["Rainbow"][height]],
PlotRange -> {{-3, 3}, {-3, 3}}, ImageSize -> Medium};

Row@{
DensityHistogram[projdata, 50, opts],
SmoothDensityHistogram[projdata, opts]
}


## Update

It might be useful to add some spatial information by superimposing a grid on our plots. The lines in light gray correspond to the latitude/longitude grid in steps of 30 degrees, the thick white line corresponds to the equator line.

<< VectorAnalysis
phiLines[theta_] := Module[{long, projlong},
long = Table[
CoordinatesToCartesian[{1, theta Degree, phi}, Spherical], {phi,
0, 2 Pi, Pi/32}];
projlong = Table[xyzToLambert[long[[i]]], {i, Length[long]}]]
thetaLines[phi_] := Module[{long, projlong},
long = Table[
CoordinatesToCartesian[{1, theta, phi Degree},
Spherical], {theta, -Pi, Pi, Pi/32}];
projlong = Table[xyzToLambert[long[[i]]], {i, Length[long]}]]

Row@{
Show[
DensityHistogram[projdata, opts],
Graphics[{LightGray, Circle[{0, 0}, 2]}, AspectRatio -> 1,
PlotRange -> {{-2.5, 2.5}, {-2.5, 2.5}}],
ListPlot[Table[phiLines[phi], {phi, Range[0, 180, 30]}],
PlotStyle -> LightGray, Joined -> True],
ListPlot[Table[thetaLines[theta], {theta, Range[0, 150, 30]}],
PlotStyle -> LightGray, Joined -> True],
ListPlot[phiLines[90],
PlotStyle -> {Thickness[.01], White}, Joined -> True]],
Show[
SmoothDensityHistogram[projdata, opts],
Graphics[{LightGray, Circle[{0, 0}, 2]}, AspectRatio -> 1,
PlotRange -> {{-2.5, 2.5}, {-2.5, 2.5}}],
ListPlot[Table[phiLines[phi], {phi, Range[0, 180, 30]}],
PlotStyle -> LightGray, Joined -> True],
ListPlot[Table[thetaLines[theta], {theta, Range[0, 150, 30]}],
PlotStyle -> LightGray, Joined -> True],
ListPlot[phiLines[90],
PlotStyle -> {Thickness[.01], White}, Joined -> True]]
}

-
Neat idea to use Lambert! :) –  Ｊ. Ｍ. Oct 5 '12 at 14:18
BTW: no need to normalize unit vectors; you can just write it as vonMisesFisherRandom[{1, 0, 0}, 3]. I put Normalize[] in my answer since I was too lazy to write the corresponding unit vector. Also: Join[vecs1, vecs2] is easier to do than Partition[Flatten[{vecs1, vecs2}], 3]. –  Ｊ. Ｍ. Oct 6 '12 at 3:36
@J.M. You're right, thanks. I was too lazy to take Normalize[] out from your code. –  VLC Oct 7 '12 at 16:35
+1 Thanks for bringing the Lambert azimuthal projection to my attention. Very useful. –  becko Nov 1 '12 at 0:55

One could map a SmoothDensityHistogram[] of the unit vectors onto a sphere as a Texture[]; the procedure is completely analogous to what I did in this answer. To illustrate:

(* generate a random unit vector following the von Mises-Fisher distribution *)
vonMisesFisherRandom[μ_?VectorQ, κ_?NumericQ] := Module[{ξ = RandomReal[], w},
w = 1 + (Log[ξ] + Log[1 + (1 - ξ) Exp[-2 κ]/ξ])/κ;
RotationTransform[{{0, 0, 1}, μ}][
Append[Sqrt[1 - w^2] Normalize[RandomVariate[NormalDistribution[], 2]], w]]]

BlockRandom[SeedRandom[42, Method -> "MersenneTwister"]; (* for reproducibility *)
vecs = Table[vonMisesFisherRandom[Normalize[{1, -2, 3}], 10], {100}]];

hist = Image[
SmoothDensityHistogram[{ArcTan[#1, #2], ArcCos[#3]} & @@@ vecs,
AspectRatio -> Automatic, ColorFunction -> "ThermometerColors",
Frame -> False, ImagePadding -> None, PerformanceGoal -> "Quality",
PlotRange -> {{-π, π}, {0, π}}, PlotRangePadding -> None],
ImageResolution -> 256];

(* you can use SphericalPlot3D[] instead *)
ParametricPlot3D[{Cos[θ] Sin[φ], Sin[θ] Sin[φ], Cos[φ]}, {θ, -π, π}, {φ, 0, π},
Lighting -> "Neutral", Mesh -> None,
PlotStyle -> Texture[hist], TextureCoordinateFunction -> ({#4, #5} &)]


Here's a version where the unit vectors are marked by tiny green spheres:

-
This is good "points are away from boundary" case. But this has seam which breaks periodicity along the azimuthal angle. Points on the different side of the seam should be "aware" of each other. What do you think? –  Vitaliy Kaurov Oct 5 '12 at 15:40
@Vitaliy, yes, since SmoothDensityHistogram[] does not know about the periodicity, there is an unsightly seam on the longitude; I haven't fully thought about how to solve this... an overlay might be one method, but this seems a bit wasteful to me. –  Ｊ. Ｍ. Oct 5 '12 at 15:52
+1 @VitaliyKaurov. Luckily, the points of my data are sufficiently far from the boundary, but I can see how it can be an issue in general... –  becko Nov 1 '12 at 0:57
@becko, I'm working on a different answer that solves the boundary problem. For now, may I recommend looking at Statistical Analysis of Spherical Data; it might be useful for your work... –  Ｊ. Ｍ. Nov 1 '12 at 1:57
@Vitaliy, I figured out an alternative; see my other answer. –  Ｊ. Ｍ. Nov 2 '12 at 21:22

I've been working at this a bit, but I've obviously done something stupid. I'm not sure why my $\phi$ only goes from $[-\pi/2, \pi/2]$, giving me a half-spherical shell. I just generated some boring uniform data which doesn't make a very interesting plot.

vectors = Normalize /@ RandomReal[{-1, 1}, {10000, 3}];
sphericalvectors = {ArcCos[#3], ArcTan[#2/#1]} & @@@ vectors;
{Min@#, Max@#} & /@ Transpose[sphericalvectors]

(*
{{0.0504573, 3.10853}, {-1.57067, 1.57018}}
*)

binned = BinCounts[sphericalvectors, {0, Pi, 0.01}, {-Pi/2, Pi/2, 0.01}];
normbinned = N[binned/Max[binned]];
pts = Flatten[Table[{a, b, normbinned[[a, b]]}, {a, 1, dims[[1]]}, {b, 1, dims[[2]]}], 1];
fint = Interpolation[pts];
SphericalPlot3D[1, {theta, 0, Pi}, {phi, -Pi/2, Pi/2},
ColorFunction ->
Function[{x, y, z, theta, pri, r},
ColorData["DarkRainbow"][fint[theta,pri]]],
ColorFunctionScaling -> False, Mesh -> False, Boxed -> False, Axes -> False]


**

old 3D voxel rendering below:

This definitely isn't a great answer, but maybe a starting point.

I just made up some data, I'm not sure if it matches yours well and then I binned the data in a 3D histogram and rendered as translucent cuboids:

vectors = RandomVariate[NormalDistribution[], {100000, 3}];
binning = {-2, 2, 0.3};
binned = BinCounts[vectors, binning, binning, binning];
dims = Dimensions@binned;
normbinned = N[binned/Max[binned]];
coordswithdata = Table[{normbinned[[x, y, z]], {x, y, z}},
{x, 1, dims[[1]]}, {y, 1, dims[[2]]}, {z, 1, dims[[3]]}];
cubes = {Hue@#1, Opacity@#1, EdgeForm[], Cuboid@#2} &
output = ParallelMap[cubes @@ # &, coordswithdata, {3}];
Graphics3D[output]
`

-
Because the data are unit vectors, the density is supported on the two-dimensional sphere $S^2$, not the 3D ball: no voxel rendering will be needed! –  whuber Oct 5 '12 at 1:30
@whuber, oops. I missed the unit part. Oh well. Enjoy my ugly pile of green and red boxes. –  s0rce Oct 5 '12 at 1:40