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?
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.
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]
}
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]]
}
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].
$\endgroup$
Commented
Oct 6, 2012 at 3:36
Normalize[] out from your code.
$\endgroup$
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:
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.
$\endgroup$
Commented
Oct 5, 2012 at 15:52
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]
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}
]
