Lifting the Contours

by
Allan Hayes
Department of Mathematics
The University
Leicester LE1 7RH
United Kingdom
hay@le.ac.uk

Mathematica in Education
Vol.3 No.3
Summer 1994
(c) TELOS/Springer-Verlag

Introduction

One of my  students wanted to know about finding local minima of a function of two variables and I thought that it would be a good idea to have a  3D picture of the function with the contour lines on the surface and color coded for height. Below are two examples given by the package. The rest of the article outlines the main steps in producing them. The details are available in the annotated package.

<<haypacks`ContourLines3D`

The surface in the first figure is  similar to what you might create using the built-in Plot3D.

fig1 = Plot3DContoured[Sin[Cos[x]+ y]/(3+ x^2+y^2),{x,-4,4},{y, -4, 4}]

Show[fig1, ColorOutput->GrayLevel,
	PlotLabel -> FontForm["Figure 1", {"Helvetica-Bold", 10}]]

Figure 1

The surface in the second figure is given by a set of parametric functions.

fig2 = ParametricPlot3DContoured[
	{(2 + Cos[t])Cos[s], Sin[t], (2+Cos[t])Sin[s]},
	{s,-Pi,Pi},{t,-Pi/2,Pi/2}
]

Show[fig2, ColorOutput->GrayLevel,
	PlotLabel -> FontForm["Figure 2", {"Helvetica-Bold", 10}]]

Figure 2

Clear["`*"]

Surface Plots: Figure 1

The system function ContourPlot gives the usual two dimensional contour map. Since I don't need the contour shading I turn it off with the option ContourShading -> False. It may make some of the explanation clearer if I define a function for the formula

f[x_,y_] := Sin[Cos[x]+ y]/(3+ x^2+y^2);

contourplot = 
	ContourPlot[
			f[x,y], {x,-4,4},{y, -4, 4},
			PlotPoints -> 25,
			ContourShading	 -> False
	]

I'm going to make use of the data in the output  ±ContourGraphics± . This ContourGraphics object includes neither the surface nor the contour lines explicitely ±  it is just a matrix of heights at the mesh points plus options  from which the contours and the contour shading can be calculated. We can look at it by using InputForm ( I cut down the number of lines displayed by using Short Ð you may have to adjust the number )

Short[InputForm[contourplot], 12]

I can get the surface by converting to SurfaceGraphics; the option Mesh -> False turns off the mesh.

surface =
	Show[ 
		SurfaceGraphics[ contourplot ], 
		Axes ->True,
		AxesLabel -> {"x","y","z"},
		Mesh -> False
	]

Show[%,ColorOutput->GrayLevel]

To get the contour lines explicitely we must convert to a Graphics object (here and in several plaes below I use Show merely to look at the display for reassurance Ð it is outside the equality)

Show[lines1 = Graphics[contourplot]]

Here is  the output

Short[InputForm[lines1], 5]

The last list is in fact list of options.. If I had not turned off the contour shading there would have been polygons as well as lines in the first list.I only want the lines .

lines2 = Graphics[Cases[lines1, _Line, {3}]]

Check the output:

Short[InputForm[lines2], 4]

To make the contours three dimensional I need to add the heights. This is done in a series of steps:

First note that each contour line is given by an expression similar to that for the first one:

lines2[[1,1]]

I shall change each pair of real numbers  {x,y} to the triple  {x,y,f[x,y]}by using a replacement rule.  I also need to change the head from Graphics to Graphics3D.

Show[
	lines3 = Graphics3D@@(lines2/.{x_?NumberQ,y_} :> {x,y,f[x,y]})
]

The picture is rather flat because Graphics3D uses by default the actual, unscaled dimensions. The scaling can be altered by means of the option BoxRatio.

Show[%53, BoxRatios -> {1, 1, .4}]

The first contour line code has been changed to.

lines3[[1,1]]

Next, I color the contour at height  z in the hue Hue[scale[z]] where scale is a function whose value runs from 0 to 1 as  z runs from the minimum value, zmin to the maximum value zmax of the z-coordinates used in surface (not in lines).

Here is the scale function:

scale[z_] := (z-zmin)/(zmax - zmin)

and  zmin and zmax are found by

{zmin,zmax} = FullOptions[surface, PlotRange][[-1]]

We get

Show[
	lines4 = 
		(lines3/.
			ln:Line[{ {__,z_},__}] :> 
				{Hue[scale[z]], ln}
		) 
]

Show[%,ColorOutput->GrayLevel]

Put the the surface and the lines together. Notice that we automatically change to Graphics3D and that the default box ratios {1,1,0.4} from surface are used, because surface comes  before lines5 in Show (try them the other way round).

Show[surface,lines4]

Show[%,ColorOutput->GrayLevel]

Well, that's not very good, parts of the contour lines are hidden behind the surface.
So, we lift them up a bit.

lines5 = lines4/.{x_,y_,z_?NumberQ} :> {x,y,z + .02(zmax-zmin)}

Show[surface, lines5]

Show[%,ColorOutput->GrayLevel]

Show[Graphics[%%], 
	PlotRange->{{0.405, 0.65}, {0.525, 0.68}}]

Show[%,ColorOutput->GrayLevel]

Here is a clip from the previous figure showing the lifted contour (the steps are: copy the preceding picture and paste it below; click on the new figure to select it  ± the bounding box will show ± ; with the Command key held down drag the bounding box handles so that the box encloses only the region that you want to look at; release the Command key and drag  the bounding box handles to make the box the size that you want; if necessary use the menu command to make lines thin).

Show[Graphics[Graphics3D[fig1 ]], 
	PlotRange->{{0.405, 0.65}, {0.525, 0.68}}]

Show[%,ColorOutput->GrayLevel]

This clearly leaves something to be desired: the package function gives the corresponding clip below. You might like to spot the differences and work out an improvement of your own Ð I will come back to this problem after dealing with parametric plots (simply reducing the height that the contours are lifted is not the answer).

Parametric Plots: Figure 2

ParametricPlot3D can, of course, give much more complicated graphics than Plot3D but, happily, the basic idea is very much the same before.  I will go through the common steps rather more quickly.
Some descriptive notation might make things clearer.

{x, y, z} = {(R+r Cos[t])Cos[s],r Sin[t], (R+r Cos[t])Sin[s]};
R = 2; r = 1;

The first problem is that I want the surface without the edges  of the constituent polygons. Unfortunately Mesh -> False  is not an option for ParametricPlot3D, so I have to insert the directive  EdgeForm[] after first generating the graphics  (with option DisplayFunction -> Identity to avoid showing it prematurely).

surface1 = 
	Show[
		Insert[
			ParametricPlot3D[
				{x,y,z},{s,-Pi,Pi},{t,-Pi/2,Pi/2},
				DisplayFunction -> Identity
			],
			EdgeForm[], 
			{1,1}
		],
		DisplayFunction -> $DisplayFunction
	]

Show[%, ColorOutput->GrayLevel]

The contours can again be found by using ContourPlot,  but for parametric plots this has to be done in two stages:

First we find the contours for the formula z in the variables  s and t.

Show[
	lines1 =
		Graphics[Cases[ 
			(*the list found by Cases needs to be wrapped in Graphics*) 
			Graphics[	
				ContourPlot[
					z, {s,-Pi,Pi},{t,-Pi/2,Pi/2},
					DisplayFunction -> Identity,
					ContourShading -> False
				],
			],
			_Line, {3}
		]]	 
]

Next we need the real contours ± in (x,y,z) -sapce. Each sub- expression  Line[{ ... {a,b} ...}] of lines1 defines one of the contours shown above; all that is needed is to replace each {a,b} with the result of replacing each s and t  in {x,y,z) with  a and b respectively, and then change the head of the result to Graphics3D.

Show[lines2 = Graphics3D@@(lines1/.{s_?NumberQ,t_} ->{x,y,z})]

Now we color the contours. This time I will just use a formula for the scaling.

{zmin,zmax} = FullOptions[surface1, PlotRange][[-1]]

Show[
	lines3 = 
		lines2/.{
			ln:Line[{{__,z_},__}] :> 
				{Hue[(z-zmin)/(zmax - zmin)], ln}
		} 
]

Show[%, ColorOutput -> GrayLevel]

and put the the surface and the lines together.

Show[surface1, lines3]

Show[%, ColorOutput -> GrayLevel]

As before, the contours need to be lifted.

lines4 = lines3/.{x_,y_,z_?NumberQ} :> {x,y,z + .01(zmax-zmin)}

Show[surface1, lines4]

Show[%, ColorOutput -> GrayLevel]

We have serious problems: the appearance is still patchy and the contours are evidently separated from the surface . Things are even more serious if we want to be able to change the viewpoint.  This is worse than for the Plot3D case; what can be done?

Why not just trace the round the polygons of the surface and construct the contour lines ? I did try this and immediately ran into two problems :(1) the resulting function was a lot slower (at least the way that I did it); (2) there was still a problem with polygons covering lines ± because, when a polygon is  rendered, it is first divided into two triangles using one of its diagonals, and a line from edge to edge of the polygon may lie behind these triangles (see the following figure where the diagonal line is added in the second graphics); still more calculation would be needed to allow for this.

Show[
	GraphicsArray[{
		Graphics3D[ gr =
		{
			Polygon[{{0,0,.4},{1,0,0},{1,1,.4},{0,1,0}}],
			Thickness[.01],
			Hue[0], Line[{{-.2,.5,.2}, {1.2,.5,.2}}]
		}
		],
		Graphics3D[{gr,Line[{{0,0,.4},{1,1,.4}}]}]	
	}]
]

Show[%, ColorOutput -> GrayLevel]

However, a reasonably efficient , high-level, answer is almost literally staring us in the face Ð  it is to pull the contour lines towards the view point. A great virtue of this is that, provided we don't overdo it, then it cannot be seen (although we could, for example, pull a contour line right through a vertical cylinder so that it would be seen as a complete circle  in front of it). The main difficulty that I found was sorting out the different coordinate sytems, but his taught me quite a lot.

The location of the view point in user coordinates depends on the BlockRatios, PlotRange and of course the value of ViewPoint in its own special coordinate system .  Let's find these for surface1

The option ViewPoint -> {vx,vy,vz} implies that these numbers are the coordinates of the viewpoint  relative to a system with origin at the center of the box, and with dimensions assigned as follows:
Suppose that the box ratios are {bx, by, bz},  where  by  is the largest, and that the dimensions of the box are {dimx, dimy, dimz}. In viewpoint coordinates the box is given the dimensions  {bx, by, bz}/by,  =  {bx/by, 1, bz/by}. So, the units for the view point coordinates in the x  y and z directions have lengths  in the user coordinates of { by/bx dimx , dimy, by/bz dimz }. Therefore, relative to the center of the box, the position of the  viewpoint in user coordinates is
		  	{vx,vy,vz}  {dimx,dimy,dimz}  by/{bx,by,bz} .

This may seem simpler in Mathematica notation.

The values of BoxRatios,PlotRange,ViewPoint used for surface are easily found.

{br,pr,vp} = FullOptions[surface1,{ BoxRatios,PlotRange,ViewPoint}]

The dimensions of the box in user coordinates are

dims = pr.{-1,1}

The position of the maximum dimension is

pmd = Position[dims, Max[dims],{1},1][[1,1]]

The center in user coordinates is

center = pr.{1,1}/2

So, the viewpoint in user coordiantes is

uvp =  center + vp dims br[[pmd]]/br

Now, define a function for moving a point a distance d towards the viewpoint.

lift[pt_,d_] :=

pt + d (( #/Sqrt[#.#])&[uvp - pt]);

Use this function to lift the contour lines

lines5 = lines3/. Line[pts_] :> Line[lift[#,.05]&/@pts]

Add the lifted contour lines to the surface.

Show[surface1, lines5]

Show[%,ColorOutput->GrayLevel]

There are, in fact, some further refinements needed, since changes in PlotRange or BoxRatios, either by user choice or by the addition of the contour lines, can change the viewpoint. This can lead to the lifted contours being being partly, or entirely, outside the box and so not completely displayed. If you want to go into these details, I have annotated the code in the package quite fully.

In the package, I have aimed to store up data so as to avoid needing to recalculate, to allow alterations by using options and Show and to provide conversion to a Graphics3D object.

It can be quite tedious to package code properly but I strongly advise it if you want to be able to benefit from your work in the future and be able to share it with others.

We are, of course, trying to approximate to curved lines and surfaces by straight line segments and triangles, so a basic limitation is the number of plot points used. Also,there is an unavoidable disadvantage in drawing the contour line this way -- the illusion requires a particular viewpoint . The prefered solution must be one in which the contours are traced on the polygons (or triangles ) used for the surface.

I should be grateful for any suggestions that are offered for improvement of these functions.