# 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?

• – a06e
Oct 4, 2012 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. Oct 4, 2012 at 23:31
• Oct 5, 2012 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, 2012 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, 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, 2012 at 8:44
• Wow, cool stuff, - and +1 for reading papers ;) Nov 4, 2012 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[])] xyzCoords[];
yy = Sqrt[2/(1 + 10^-10 - xyzCoords[])] xyzCoords[];
{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,
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,
PlotStyle -> {Thickness[.01], White}, Joined -> True]]
}

• Neat idea to use Lambert! :) Oct 5, 2012 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, 2012 at 3:36
• @J.M. You're right, thanks. I was too lazy to take Normalize[] out from your code.
– VLC
Oct 7, 2012 at 16:35
• +1 Thanks for bringing the Lambert azimuthal projection to my attention. Very useful.
– a06e
Nov 1, 2012 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? Oct 5, 2012 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, 2012 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...
– a06e
Nov 1, 2012 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, 2012 at 1:57
• @Vitaliy, I figured out an alternative; see my other answer. Nov 2, 2012 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[]}, {b, 1, dims[]}], 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[]}, {y, 1, dims[]}, {z, 1, dims[]}];
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! Oct 5, 2012 at 1:30
• @whuber, oops. I missed the unit part. Oh well. Enjoy my ugly pile of green and red boxes. Oct 5, 2012 at 1:40

The following approach constructs a histogram based on the triangulation of a sphere, or better: its dual. It is much faster than the method shown by J.M., and thus it is useful for cases when there's a lot of data. However, it is only usable when there is a lot of data as it does not perform smoothing.

We can take the dataset vecs from J.M.'s answer, but generate 10000 points instead of 150.

Let's make a nice polyhedron to use as the basis of histogramming.

reg = BoundaryDiscretizeRegion[
Ball[], PrecisionGoal -> 1, MaxCellMeasure -> 1/5,
MeshCellStyle -> {1 -> None}
] Visually, it is made of equally sized triangles. That is not exactly true, but luckily it is easy to measure the size of the triangles:

cellSizes = PropertyValue[{reg, 2}, MeshCellMeasure];


If the triangles are small enough, these approximate well the corresponding area of the spherical surface.

Let us take the vector pointing to the centre of each triangle, and construct a NearestFunction based on the CosineDistance.

nf = Nearest[
PropertyValue[{reg, 2}, MeshCellCentroid] -> "Index",
DistanceFunction -> CosineDistance
]


Each unit vector from vecs will fall into the triangle whose centre it is the closest to.

indices = First /@ nf /@ vecs; // AbsoluteTiming
(* {5.52186, Null} *)

counts = Counts[indices];


Once we have the corresponding counts, we can plot them:

normCounts = Lookup[counts, Range@MeshCellCount[reg, 2], 0]/cellSizes;
normCounts = normCounts/Max[normCounts];

SetProperty[{reg, {2, All}},
MeshCellStyle -> ColorData["Rainbow"] /@ normCounts] The problem is that triangles are not a very nice shape for a histogram.

We can do better by taking the dual of this polyhedron by taking the centre of each triangle. It would look like this (using IGraph/M):

IGMeshCellAdjacencyGraph[reg, 2, VertexCoordinates -> Automatic,
VertexSize -> Medium] We won't convert this to a mesh region though, as for the histogramming we only need the centre of each face anyway.

A difficulty here is measuring the area of each polygon, which I will leave for the next version of this answer. Let's proceed without that for now.

nf = Nearest[
MeshCoordinates[reg] -> "Index",
DistanceFunction -> CosineDistance
]

indices = First /@ nf /@ vecs;

counts = Counts[indices];


To plot the result, we use GraphicsComplex with its VertexColors and VertexNormals options. Since this surface is a sphere at the original, the vertex normals will be the same as the vertex vectors. Ambient lighting makes the colouring independent of the rotation of the sphere (there won't be any shadows).

Graphics3D[
GraphicsComplex[MeshCoordinates[reg],
{EdgeForm[None], MeshCells[reg, 2]},
VertexColors -> ColorData["MintColors"] /@ Lookup[(counts/Max[counts]), Range@MeshCellCount[reg, 2], 0],
VertexNormals -> MeshCoordinates[reg]
],
Lighting -> {"Ambient", White}
]
` • neat addition to this. +1 Jun 4, 2018 at 17:13
• "It is much faster than the method shown by J.M." - heh. I need to come back to optimize it someday, but I haven't had the time to read all the relevant references. Oct 5, 2018 at 12:36