10
$\begingroup$

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?

Improve this question
$\endgroup$
8
  • $\begingroup$ 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. $\endgroup$
    – Henry B
    Mar 14, 2012 at 5:32
  • 1
    $\begingroup$ Please include an example in your question. Do not leave us guessing. $\endgroup$
    – Mr.Wizard
    Mar 14, 2012 at 5:35
  • $\begingroup$ Also, welcome to Mathematica.SE. $\endgroup$
    – Mr.Wizard
    Mar 14, 2012 at 5:36
  • $\begingroup$ Is there a reason you aren't using EuclideanDistance[#,{x,y}]&/@ data and are instead rolling your own distance function? $\endgroup$
    – Verbeia
    Mar 14, 2012 at 6:35
  • 4
    $\begingroup$ 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. $\endgroup$
    – Matariki
    Mar 14, 2012 at 7:04

2 Answers 2

Reset to default
8
$\begingroup$

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]
}]

Figure

Improve this answer
$\endgroup$
8
  • $\begingroup$ 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 $\endgroup$
    – Szabolcs
    Mar 14, 2012 at 16:58
  • $\begingroup$ 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 :-(. $\endgroup$
    – whuber
    Mar 14, 2012 at 17:02
  • $\begingroup$ 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. $\endgroup$
    – whuber
    Mar 14, 2012 at 17:04
  • $\begingroup$ "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. $\endgroup$
    – Szabolcs
    Mar 14, 2012 at 17:24
  • $\begingroup$ 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$. $\endgroup$
    – whuber
    Mar 14, 2012 at 17:46
6
$\begingroup$

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]

Mathematica graphics

ListContourPlot[pts, InterpolationOrder -> 1]

Mathematica graphics

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

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

Mathematica graphics

(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)

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ +1 The last plot is actually a reflection, not a rotation. $\endgroup$
    – whuber
    Mar 16, 2012 at 14:01
  • $\begingroup$ @whuber True, a reflection in the x==y axis. $\endgroup$
    – Sjoerd C. de Vries
    Mar 16, 2012 at 14:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.