# Colorize a 3D distribution of points according to density

Suppose we have some random 3D distribution of points. I'll use the spherical ball defined below as an example and a starting point for the discussion here :

Ball[num_]:=Table[
{
#1 Sqrt[1-#2^2]Cos[#3],
#1 Sqrt[1-#2^2]Sin[#3],
#1 #2
}

&[
Random[NormalDistribution[1, 0.5]],
Random[Real,{-1,1}],
Random[Real,{0,2Pi}]
],{num}]

Graphics3D[{AbsolutePointSize[2],Point[Ball[10000]]},Boxed->True,BoxRatios->{1,1,1},ImageSize->800,SphericalRegion->True]


This code produces the ball shown here :

I would like to colorize the whole distribution with a color shade defined on the density of points, from color1 (highest density of points) to color2 (lowest density of points). Adding more colors to the shade is a desirable option (color1 to color2 to color3 ... to colorN). How can I do that ? What Mathematica code could do this ?

Please, the suggestions should be compatible with Mathematica 7.0.

-
Unfortunately SmoothKernelDistribution is new in 8.0. – Rahul Mar 5 '14 at 19:08
Does this work in v7? arr = Array[Exp[-0.01 {##}.{##}] &, {20, 20, 20}]; ListInterpolation[arr]. (I don't have v7 so I don't know.) – Szabolcs Mar 5 '14 at 19:28
Apparently, this small code works in v7.0, but it doesn't ouput anything. Not sure if it's really working. – Cham Mar 5 '14 at 19:49

There are two ways:

1. Calculate the density analytically. For the distribution you use this is difficult but since this for producing something pretty, and not for accuracy, you can consider using a different distribution.

2. Approximate the distribution numerically.

I'm going to do no. 2. below. I don't have version 7, so it is just a guess that these functions will work. In later versions SmoothKernelDistribution is going to be a better an easier way to estimate the density, but in v7 we can't use that so I went for simple binning.

Summary of the idea:

1. Estimate the density by binning. The result will be sensitive to the bin size. A large bin size might smooth out the data more than desired while a small bin size will increase fluctuations.

2. Interpolate the histogram to obtain a smooth function. Linear interpolation might be the best choice if the histogram has a lot of fluctuations. Higher order interpolation tends to exacerbate the fluctuations.

3. Colour the points based on this interpolated estimate of the density.

Here's code that implements this:

This data will only be used for estimating the density. More data points give a smoother histogram.

In[74]:= data = Ball[500000];


Figure out reasonable bounds for the binning:

In[75]:= Max /@ Transpose[data]
Out[75]= {2.8324, 2.78793, 2.65889}

In[76]:= Min /@ Transpose[data]
Out[76]= {-2.85683, -3.18866, -2.73881}


-3..3 will do, but you might need to increase it if points tend to fall outside of this range. Also make sure to divide by the maximum value so that we have numbers between 0 and 1 (that can be used in Mathematica's colour functions).

In[77]:= bins = N@BinCounts[data, {-3, 3, .1}, {-3, 3, .1}, {-3, 3, .1}];
bins /= Max[bins];

In[79]:= if =
ListInterpolation[bins, {{-3, 3}, {-3, 3}, {-3, 3}},
InterpolationOrder -> 1]


Now we have a density function to colour by:

Graphics3D[{AbsolutePointSize[2],
{ColorData["Rainbow"][if[##]], Point[{##}]} & @@@ Ball[10000]},
Boxed -> True, BoxRatios -> {1, 1, 1}, ImageSize -> 800, SphericalRegion -> True]

-
This code works in v7.0, but the colour ouptput is pretty weird. I'm also getting some warnings about values lying outside the range of data in the interpolating function. Extrapolation will be used. – Cham Mar 5 '14 at 19:58
@Cham This is meant as a basic example that you can adapt to your needs, not a complete solution. All the necessary parts are there, but I expect you'll want to tweak it. You can change the colour function. If points tend to fall outside the range {-3,3}, just increase it. You'll need to look up the functions I used and understand what they do (they're not complicated). – Szabolcs Mar 5 '14 at 20:02
How do you change the ColorData["Rainbow"] to a blend from Color1 to Color2 ? – Cham Mar 5 '14 at 20:09
@Cham Search the docs for "Blend". – Szabolcs Mar 5 '14 at 20:10
@Cham It seems to me that the problem is that you are just taking the code and trying to use it without understanding what it does. Please read the code, and look up the functions you are not familiar with. After that it should be clear to you why it produces a certain result and how you can modify it. If after looking up the functions you still have questions about some details, feel free to ask, but please do make an effort do understand it on your own. The choice of binning does clearly have an effect on the result (because it effectively does some smoothing). – Szabolcs Mar 5 '14 at 20:23

Here is a brute force approach, approximating the density by simply counting near neighbors:

 all = Ball[10000];
blz = {Count[ all, p_ /; (Norm[p - #] < .5)], #} & /@ all ;


..Go get lunch..

Then directly color each point..

 Graphics3D[{AbsolutePointSize[3],
{Hue[N[(2/3) Log[#[[1]]]/7]], Point[#[[2]]]} & /@ blz},
Boxed -> True, BoxRatios -> {1, 1, 1}, ImageSize -> 800,
SphericalRegion -> True]


You can speed this up quite a lot using FindClusters first:

 blz = Flatten[
Function[cc, ({Count[cc, p_ /; (Norm[p - #] < .5)], #} & /@ cc )] /@
FindClusters[Ball[10000] , 10 ], 1];


Be aware this is introducing an additional approximation, as you are only finding neighbors within a cluster. Here is a look at what FindClusters does..looping over each of 10 clusters. As you can see points on the boundary of their cluster are going to get an artificially low density count.

-
+1 for the go get lunch line :D – Sektor Mar 5 '14 at 22:33

This is not particularly clever or efficient, but you could use FindClusters to break the data into pieces and then do a rough density calculation on each piece.

Ball[num_] :=
Table[{#1 Sqrt[1 - #2^2] Cos[#3], #1 Sqrt[
1 - #2^2] Sin[#3], #1 #2} &[Random[NormalDistribution[1, 0.5]],
Random[Real, {-1, 1}], Random[Real, {0, 2 Pi}]], {num}]

colorize[list_List] := Module[{d, densities, min, max},
d = If[Length[#] > 1, Log[Length[#]/Times @@ StandardDeviation[#]], 0] &;
densities = Map[d, list];
{min, max} = {Min[#], Max[#]} &@densities;
{ColorData[{"AvocadoColors", {max, min}}][d[#]], Point[#]} & /@ list]

Graphics3D[colorize[FindClusters[Ball[10000],50]], PlotRange-> {Automatic,{0,Automatic},Automatic}]


you will need to adjust the number of clusters for your preferences.

-
Hmm, this code gives me an error message : AvocadoColors is not a known entity, class or tag for ColorData. I guess it's the FindClusters which isn't working correctly. – Cham Mar 5 '14 at 21:53