# 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_]=(#[]-x)^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. – Henry B Mar 14 '12 at 5:32
• Please include an example in your question. Do not leave us guessing. – Mr.Wizard Mar 14 '12 at 5:35
• Also, welcome to Mathematica.SE. – Mr.Wizard Mar 14 '12 at 5:36
• Is there a reason you aren't using EuclideanDistance[#,{x,y}]&/@ data and are instead rolling your own distance function? – Verbeia Mar 14 '12 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. – Matariki Mar 14 '12 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, BinCounts[data, .1, .1]] // ListContourPlot – Szabolcs Mar 14 '12 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 :-(. – whuber Mar 14 '12 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. – whuber Mar 14 '12 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. – Szabolcs Mar 14 '12 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$. – whuber Mar 14 '12 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. – whuber Mar 16 '12 at 14:01
• @whuber True, a reflection in the x==y axis. – Sjoerd C. de Vries Mar 16 '12 at 14:57