To plot bivariate functions, we typically either use Plot3D, whose dynamic effects allow flexible investigation, or DensityPlot, which is better for publications. However, sometimes the output of DensityPlot is a little ... obscure(?), and shadows from ReliefPlot would add a lot perspective.

Example from the ReliefPlot documentation:

ReliefPlot[Table[i + Sin[i^2 + j^2], {i, -4, 4, .03}, {j, -4, 4, .03}], ColorFunction -> "SunsetColors"]

Compare a replication using DensityPlot:

DensityPlot[i + Sin[i^2 + j^2], {i, -4, 4}, {j, -4, 4}, ColorFunction -> "SunsetColors", PlotPoints -> 40]

Where the well-known oddity of ReliefPlot's orientation can be seen.

For years on MMA SE, a lot of members have used ReliefPlot via manual sampling for function plotting questions (e.g., this post), yielding better aesthetic results.
Alternatives include:

  • Using Plot3D and set the viewpoint at infinity above, but it can also hairy to do. (See @MichaelE2's answer which does this elegantly)
  • Manually define normal vectors and do a normal mapping, which I am not sure how to implement in Mathematica right now. (See my answer to get a crude idea)

But why isn't ReliefPlot for functions been officially implemented (i.e., with adaptive sampling like other *Plot* built-ins)? What is the real technical difficulty?

P.S.: I believe extensive about this problem can already be found on MMA SE (like here and here), but they are scattered everywhere so I am raising this question as a port.

For this thread to serve as a guide to others with the same problem, what are the ways to mimick ReliefPlot's texture? (please provide simple code and output examples) From your personal experience, at what point should one bother to do this for the nicer touch?


ReliefPlot is using the 3D lighting model (I think) to highlight changes in the gradient.

Deploy@Plot3D[i + Sin[i^2 + j^2], {i, -4, 4}, {j, -4, 4}, 
  ColorFunction -> (ColorData["SunsetColors"][#3] &), 
  PlotPoints -> 35, MaxRecursion -> 4, Mesh -> None, 
  ViewPoint -> {0, 0, Infinity}, 
  Lighting -> {{"Ambient", GrayLevel[0.75]}, {"Directional", 
     GrayLevel[0.15], ImageScaled[{-2, 0, 2}]}, {"Directional", 
     GrayLevel[0.15], ImageScaled[{-2, -2, 2}]}, {"Directional", 
     GrayLevel[0.15], ImageScaled[{0, -2, 2}]}}, Boxed -> False, 
  Axes -> False]
  • $\begingroup$ Really elegant results. But I do wonder if ReliefPlot actually constructs a 3D surface internally; it is a function often used on large geographic datasets, so actual ray tracing would be quite expensive. I do wonder that if the plot is complicated (e.g., when poles occur), employing Plot3D may lead to performance issues. I suspect ReliefPlot calculates the gradient of the array and do a normal mapping accordingly. $\endgroup$ – Gravifer Mar 9 at 4:52
  • 1
    $\begingroup$ @Gravifer Well, yeah, of course. Plot3D does not do raytracing either, but calculates the normals (symbolically or numerically, depending on the function). Graphics3D uses the normals (not Plot3D) and Lighting to determine the colors of the polygons. (Unless that is done by the GPU. I'm not really up on exactly what GPUs are used to do.) ReliefPlot constructs a Raster of the colors that it determines. $\endgroup$ – Michael E2 Mar 9 at 4:55
  • $\begingroup$ I am a little confused. I just verified that ReliefPlot returns a Graphics object; does it calls upon Graphics3D functionalities internally? I didn't realize that it uses polygon based reflection rendering. $\endgroup$ – Gravifer Mar 9 at 5:12
  • 1
    $\begingroup$ I don't know how ReliefPlot actually does what it does internally, but only inferred it from how well it matches the Plot3D output. It probably uses the same basic algorithm for 3D shading to determine how to color each square in the raster, but it might use a routine optimized for producing a 2D raster, instead of rasterizing a 3D plot. Note that ReliefPlot has an option LightingAngle. That and other language in the docs suggest that it produces a 2D rendering of a 3D surface. The Method option can be used to switching between two shading algorithms. $\endgroup$ – Michael E2 Mar 9 at 5:50
  • 1
    $\begingroup$ @Gravifer This might give a slightly closer reproduction, though the shadow transitions are still not a diffuse: Plot3D[i + Sin[i^2 + j^2], {i, -4, 4}, {j, -4, 4}, ColorFunction -> (ColorData["SunsetColors"][#3] &), PlotPoints -> 265, MaxRecursion -> 0, Mesh -> None, ViewPoint -> {0, 0, Infinity}, Lighting -> {{"Ambient", GrayLevel[0.85]}, {"Directional", GrayLevel[0.25], ImageScaled[{-1, 1, Sqrt[2]}]}}, Boxed -> False, Axes -> False] $\endgroup$ – Michael E2 Mar 9 at 14:42

As this is a guideline question, I'll place my own research in a standalone answer.
The ideas comes from the discussion with @MichaelE2 under his answer following the Plot3D approach.

I mentioned there that it is possible do do artificial normal mapping and adding the shadow as a mask over the original DensityPlot. The following code is an example for how to do this:

theta=3Pi/4; phi=Pi/4;
func = j + Sin[i^2 + j^2];
base = DensityPlot[func, {i,-4, 4}, {j,-4, 4}, ColorFunction->"SunsetColors", PlotPoints->40];
shadow = Insert[
       FromSphericalCoordinates[{1,theta,Pi-phi}](*Default source position of ReliefPlot*)
       . Normalize@Append[-Grad[func, {i, j}], 10](*2nd argument of GrayLevel is Opacity*)
     ], {i,-4, 4}, {j,-4, 4}, ColorFunction->(GrayLevel[#,(2#-1)^2]&), PlotPoints->80],
 Opacity[.5], {1, 1}];
GraphicsRow[{base, shadow,
  Show[{base, shadow}],
  ReliefPlot[Table[i + Sin[i^2 + j^2],
    {i,-4, 4, .03}, {j,-4, 4, .03}], ColorFunction->"SunsetColors"]}]

mask approach result
Note that:
(1) I switched i and j so the orientation matches.
(2) The $z$ component I appended to the -Grad vector is $10$ instead of $1$; this effectively shrinked the surface and smoothed the shadows. You may tweak this to achieve better visual effects.

The output is still isn't really quite the same to the built-in, but can show the gist of it. I used a quadratic opacity scaling in the mask; it can be modified to improve the output. What you need to know is that in a true normal mapping, the shadow mask is multiplicative; a additive mask will never create highlight faithfully.

Note that there are a few artifacts: the direction of the normal is somewhat 'singular' where the gradient vanishes.

Pack this to a single function (can be really buggy)

reliefPlot[f_, {x_, xmin_, xmax_}, {y_, ymin_, ymax_}, 
  ops : OptionsPattern[{LightingAngle -> {(3 \[Pi])/4, \[Pi]/4}, 
     DensityPlot, ReliefPlot}]] := 
 Show[{DensityPlot[f, {x, xmin, xmax}, {y, ymin, ymax}, ops], 
         OptionValue[LightingAngle][[1]], \[Pi] - 
          OptionValue[LightingAngle][[-1]]}] . 
       Normalize@Append[-Grad[f, {x, y}], 10]], {x, xmin, xmax}, {y, 
      ymin, ymax}, ColorFunction -> (GrayLevel[#, (2 # - 1)^2] &), 
     Evaluate[Sequence @@ FilterRules[{ops}, PlotPoints]]], 
    Opacity[.5], {1, 1}]}]

(If I have tried other ideas, I will update them to this answer too.)


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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