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, 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[])] 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! :) – J. M. will be back soon 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]. – J. M. will be back soon 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. – J. M. will be back soon 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... – J. M. will be back soon Nov 1 '12 at 1:57
• @Vitaliy, I figured out an alternative; see my other answer. – J. M. will be back soon 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[]}, {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! – 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

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 – rcollyer Jun 4 '18 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. – J. M. will be back soon Oct 5 '18 at 12:36