Colorize a 3D distribution of points according to density

Question

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.

Answer 1

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]

Answer 2

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.

Answer 3

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.