InterpolationOrder for ContourPlot

How do I acheive something like InterpolationOrder in ListContourPlot for ContourPlot? My main problem is I want to smooth the contours in ContourPlot.

I have some points which I'm going to call data, which is 10,000 points and looks like this:

{{24.3,40.0},{29.2,56.0} ... }


I create this distance function:

D[x_,y_]=(#[[1]]-x)^2+(#[[2]]-y)^2&/@data


Now, my function is supposed to count the number of points within a certain distance of any x,y position and divide by the total number of points:

F[x_,y_,r_]:=Length[Select[D[x,y],(#<r^2)&]]/Length[data]


Now I can plot my function,

ContourPlot[F[x,y,5],{x,0,100},{y,0,100}]


The problem is my function is not naturally smooth.

I know about SmoothDensityHistogram and DensityHistogram, but I like the more simple look of smoothed contours. Can anyone help with this?

• The problem is ContourPlot takes a function. I would use ListContourPlot, but I want to superimpose another plot on top of the contour and I lose the (x,y) information with ListContourPlot. Commented Mar 14, 2012 at 5:32
• Please include an example in your question. Do not leave us guessing. Commented Mar 14, 2012 at 5:35
• Also, welcome to Mathematica.SE. Commented Mar 14, 2012 at 5:36
• Is there a reason you aren't using EuclideanDistance[#,{x,y}]&/@ data and are instead rolling your own distance function? Commented Mar 14, 2012 at 6:35
• One general suggestion: Don't define functions/symbols that start with a capital letter as they might shadow or clash with MMa defined functions and symbols. An example is your use of D which MMa is using for the derivative. Commented Mar 14, 2012 at 7:04

The most efficient way to carry out your task, which is to plot a contour map of a kernel density of your points, is by converting the points to raster format and using a Fast Fourier Transform to convolve them with a density kernel. But that takes some work. If you're willing to wait a few seconds, the whole procedure is (less efficiently) built into Mathematica's SmoothKernelDistribution function.

Here is an example taken, with minor changes (to make it more interesting), directly from the help page:

(* Create some data--around 10,000 points--for the illustration *)
data = Join @@ Table[RandomVariate[BinormalDistribution[m, {1/2, 1/2}, 0], 1500],
{m, RandomReal[{1, 9}, {7, 2}]}];

(* Create a rough (D1) and smooth (D2) density for contouring *)
D1 = SmoothKernelDistribution[data, 0.02]; (* Takes a few seconds *)
D2 = SmoothKernelDistribution[data, 0.5];  (* Takes a few more seconds *)

(* Plot the points and their densities *)
points = ListPlot[data, PlotRange -> {{0, 10}, {0, 10}}, AspectRatio -> 1];
TableForm[{
Prepend[
Table[
ContourPlot[
Evaluate@PDF[D, {x, y}], {x, 0, 10}, {y, 0, 10},
PlotRange -> All,
ColorFunction -> "TemperatureMap"],
{D, {D1, D2}}
], points]
}]


• Actually your code runs near instantaneously here. Did you use a lot more points when you found that it takes a few seconds? If you don't want to use SmoothKernelDistribution because of performance concerns (actually I don't know how it's implemented...), and you'd rather use something FFT based, that's easy to do as well: ListConvolve[GaussianMatrix[5], BinCounts[data, .1, .1]] // ListContourPlot Commented Mar 14, 2012 at 16:58
• Thanks. I don't think ListConvolve, as directly applied to the data, does what the OP intends. You first need to create a 2D array, somewhat as in your reply here, but you don't evaluate the density F on a regular grid: you grid the points (as sums of indicator functions) and then you could use ListConvolve et al. with a suitably chosen kernel. SmoothKernelDistribution really bogged down for me when using custom kernels; at one point it ate all my RAM and forced a reboot :-(. Commented Mar 14, 2012 at 17:02
• I see you're using BinCounts to rasterize the points, so that's ok: nice, simple solution. In practice, far more bins would likely be desirable. Commented Mar 14, 2012 at 17:04
• "The most efficient way to carry out your task, which is to plot a contour map of a kernel density of your points, is by converting the points to raster format and using a Fast Fourier Transform to convolve them with a density kernel." <-- Doesn't this require a regular binning? Is there a way to use the FFT to speed up the convolution when we don't have a function sampled on a regular grid? Yes, you are right that the BinCounts + ListConvolve method will get slow quickly as the bin size is decreased. Commented Mar 14, 2012 at 17:24
• It might be worth some testing and benchmarking. In principle, using sparse matrices, BinCounts is just a $O(n)$ operation (time and storage) for $n$ input points. ListConvolve--or at least the FFT equivalent--will need $O(1/h^2)$ storage and $O(-\log(h)/h)$ time for a bin width of $h$. Commented Mar 14, 2012 at 17:46

I'd say you could solve this problem by evaluating your function F on a regular grid and feeding this data to ListContourPlot with the InterpolationOrder set at an appropriate level.

data = RandomVariate[
MultinormalDistribution[{50, 50}, 150 {{1, 0.8}, {0.8, 4}}], 10000];
dist[x_, y_] := EuclideanDistance[#, {x, y}] & /@ data;
f[x_, y_, r_] := Length[Select[dist[x, y], (# < r^2) &]]/Length[data]
pts = Table[f[x, y, 5], {x, 0, 100}, {y, 0, 100}];

ListContourPlot[pts, InterpolationOrder -> 0]


ListContourPlot[pts, InterpolationOrder -> 1]


The job ContourPlot does isn't that bad, by the way:

ContourPlot[f[x, y, 5], {x, 0, 100}, {y, 0, 100}]


(Yes, I know the plot is rotated with respect to the above ones. I forgot to order x and y in the correct order in the Table function)

• +1 The last plot is actually a reflection, not a rotation. Commented Mar 16, 2012 at 14:01
• @whuber True, a reflection in the x==y axis. Commented Mar 16, 2012 at 14:57