Use a symmetric mesh in ListDensityPlot

Question

I am trying to produce a ListDensityPlot but I am having issues with aliasing. The problem arises from the way that Mathematica is breaking up the mesh into triangles, such that the orientations of the triangles is always in one direction, such that along one edge where there is a steep change in the function it appears smooth but along another edge is appears jagged. This is purely artificial due to the nature of the mesh, and could be solved if I could use a mesh which was somewhat more isotropic or which was axisymmetric about the center, but I cannot figure out how to do this. I'm surprised the mesh does not adapt where the derivative becomes large.

I am attaching two plots, with and without the mesh to demonstrate the issue. You can see aliasing on the top left and bottom right edges, whereas the top right and bottom left edges appear smooth. I am pretty certain this is caused by the mesh, as you can see in the second image the triangles are oriented differently along these different edges. The data set used to generate the plot is a simple square lattice. All help appreciated.

ListDensityPlot[data, ColorFunction -> "SunsetColors"]

ListDensityPlot[data, ColorFunction -> "SunsetColors", Mesh->All]

Edit:

Changing the interpolation order does not fix the issue, as per Henrik's response:

GraphicsGrid[
  Partition[
    Table[ListDensityPlot[data, ColorFunction -> "SunsetColors", 
      InterpolationOrder -> i], {i, 1, 6}],
      3],
  ImageSize -> Full]

Data used to generate this plot can be found here or here.

Answer 1

Just a random example:

f = {x, y} \[Function] Sin[7 Pi (x + y)] + 0.7 Cos[12 (-1 x + 2 y + 1.)];
n = 25;
{x, y} = Transpose[Outer[List, Subdivide[0., 1., n], Subdivide[0., 1., n], 1], {2, 3, 1}];
data = f[x, y];
GraphicsGrid[
 Partition[
  Table[ListDensityPlot[data, InterpolationOrder -> i], {i, 1, 6}],
  3], ImageSize -> Large
 ]

So the mesh is not necessarily the issue.

Usually, the option MaxPlotPoints should help. But for some reason, it doesn't. So, let's upsample the image by hand: After download, we load the data, create an interpolation function and sample that on a finer grid.

data = Import["example_SF.mx"];
f = Interpolation[data, InterpolationOrder -> 3];
resolution = 400;
upsampledvals = Outer[f, 
 Subdivide[0., 4., resolution], 
 Subdivide[0., 4., resolution], 
 1];

Now we have the upsampled data but how to render it? Fortunately, @Wizard provided us in this answer with an efficient method to use Image with a costumized color pallete:

renderImage[array_?MatrixQ, cf_, q_Integer: 2048, 
  opts : OptionsPattern[Image]] := 
 Module[{tbl}, 
  tbl = List @@@ Array[cf, q, {0`, 1`}] // N // 
    Developer`ToPackedArray;
  Image[tbl[[# + 1]] & /@ Round[(q - 1`) array], opts]]

Applying it to the upsampled date leads to

img = renderImage[Rescale[upsampledvals], ColorData["SunsetColors"]]

So far, this missed the axes. Those can be added like this:

Graphics[
 Inset[img, {0, 0}, {0, 0}, {4, 4}],
 Axes -> True,
 Frame -> True,
 PlotRange -> {{0, 4}, {0, 4}},
 PlotRangePadding -> 0.1
 ]