IDL Surface Tutorial - Two Intersecting Surfaces

[ Lund University ] [ Fysicum ] [ Maxlab ] [ Sljus ] [ STM Group ]

Part of Struan's IDL Surface Plotting Tutorial

Plotting Two Intersecting Surfaces

Often we want to display two surfaces on the same plot in order to compare them directly. IDL provides a built in z-buffer which does hidden surface and line removal on such plots so that intersecting surfaces can be correctly displayed. All the plotting commands that can be used on the screen device can also be used in the z-buffer, so it is easy to produce plots like those above which show wireframe, shaded and coloured views of the same two interlocking surfaces.

However, just as with single surfaces the standard plots do not always give a clear picture of the data. The grid view can be extremely confusing because there are no depth clues for the viewer to use to decide which surface is at the front and which is at the back. The shaded and coloured views give a more accurate presentation of the frontmost surfaces but completely hide the back parts and, as with the single surface plots, the curvature of the surfaces is not always clear.

By combining the images in various ways it is possible to get round these problems and produce a better looking and more informative plot, such as the one below.

This page of the tutorial describes the steps needed to create this image, but assumes that you have already read the tutorial on plotting single surfaces and understand the operation of the utility routines TV_24 and TVRD_24 from the example code page.

Shading Coloured Intersecting Surfaces

This is done in exactly the same way as with the single surface, except that the actual plotting of the individual pictures is done in the z-buffer and only the final image is displayed on the screen device. As before, we create the coloured view of the two surfaces and then multiply by the pixel values of the shaded view to give the impression of shading the coloured surface.

Implementation

We begin by creating the two test surfaces. One is a gaussian hat just as before and the other is the sum of two half-height hats displaced either side of the first with a col or saddle point between the two.

    surfone = shift(dist(20,16), 10, 8)
    surfone = bytscl(exp(-(surfone/5)^2))
    surftwo = bytscl(fix(shift(surfone, 6, 0)) + shift(surfone, -6, 0))/2b

At this point it is worth creating some utility variables: a string to hold the name of the screen device (this line will only work if the screen is your current plotting device, ie if you haven't already typed 'set_plot, 'z''), and two variables to hold the size of the plotting window (we make these long integers because we'll often be multiplying them together and we don't want to have to worry about overflow errors). Finally, we need to set the default z-range of our plots so that the built in surface plotting routines don't scale the second surface to be the same height as the first.

    screen = !D.name
    xwin = 400L
    ywin = 400L
    !z.range = [0,255]

Now we are ready to use the z-buffer to plot a shaded version of the image:

    set_plot, 'z'
    device, set_resolution=[xwin, ywin]
    shade_surf, surfone
    shade_surf, surftwo, /noerase
    shade_array = tvrd()

Then we create the colour arrays for the colour coded image, plot it and read it back with the TVRD_24 procedure which automatically converts it to a pixel interleaved, 24-bit, RGB image.

    surfone_col = 128b + surfone/2b
    surftwo_col = 127b - surftwo
    shade_surf, surfone, shades=surfone_col
    shade_surf, surftwo, shades=surftwo_col, /noerase
    loadct, 25
    colour_surf = tvrd_24()

We have used a small trick to plot the two surfaces with different colours. The z-buffer is only an 8-bit colour device (something I fervently hope will be changed in IDL 5.0 hint, hint) so it is necessary to share the 256 colours between the two surfaces. A simple way of doing this is to use the bottom half of the colour table for one surface and the top half for the other, and if we also invert the values for one surface we can ensure that they both have the same colour at their minimum heights. This is done by creating an array for the SHADES keyword which runs from 128 up to 255 for surfone and from 127 down to 0 for surftwo. Naturally, a multitude of other colour schemes are also possible.

All that remains is to multiply the individual RGB planes of the colour image by the pixel intensities of the shaded image, rescale to the range [0,255] and display the image on the screen device.

    shade_col = long(colour_surf)
    for i=0,2 do shade_col(i,*,*) = shade_col(i,*,*) * shade_array
    shade_col = byte(shade_col/255)

    set_plot, screen
    window, /free, xsize=xwin, ysize=ywin
    tv_24, shade_col

Grids Revisited

If we want to overlay a grid onto an image created in the z-buffer we have to use a masking or overlaying technique because of the submarining discussed on the single surface plotting page.

Because of the way the z-buffer works, producing a suitable grid to manipulate a shaded or coloured image is not a trivial exercise. The image at the top of the page shows the standard output from two calls to SURFACE in the z-buffer, which includes lines behind what would be the frontmost surface of the coloured or shaded plot and will give the wrong results if used to mask them. We might be tempted to plot a surface the same colour as the background to hide the lines behind it, but then we are back to the submarining problem.

The trick is to create the visible grids for each surface separately and then add them together to make a grid for the whole plot. For each surface we plot a normal grid and then plot the other surface as a colour plot where all the entries in the colour array passed to the SHADES keyword are set to the background colour. The idea is illustrated below:

In the left hand two images, which show the trick for each surface in turn, we have left the second surface grey so that it can be seen. Plotting it using the background colour would leave only the visible portions of the grids in the image, and these then can be combined to produce the image on the right which shows only the parts of the grids which lie on surfaces visible in a shaded or coloured version of the plot.

We can then take the completed grid and use it to mask or overlay an image as before. Better still, we can take the individual grids and use them to overlay a different colour grid onto each surface, as has been done here:

Implementation

The construction of the grids is very simple. To make the grid for the first surface we use the normal SURFACE procedure and then overplot with a coloured version of the second surface using a dummy array full of zeros to force the colour of the surface to be whatever the background is.

    loadct, 0
    set_plot, 'z'
    device, set_resolution=[xwin,ywin]
    surface, surfone
    shade_surf, surftwo, shades=bytarr(xwin, ywin), /noerase
    grid_one = tvrd()

We then repeat this for the second surface:

    surface, surftwo
    shade_surf, surfone, shades=bytarr(xwin, ywin), /noerase
    grid_two = tvrd()

We can then combine the two grids into a single one for both interlocking surfaces:

    surf_grid = grid_one OR grid_two

Alternatively, we can use the grids separately to give each surface a different coloured grid. For the 24-bit shaded, coloured image created above the code would look something like this:

    shade_col_grid = reform(shade_col, 3, xwin*ywin)
    shade_col_grid(*,where(grid_one ne 0)) = 20b
    shade_col_grid(*,where(grid_two ne 0)) = 230b
    shade_col_grid = reform(shade_col_grid, 3, xwin, ywin, /overwrite)

    set_plot, screen
    tv_24, shade_col_grid

Transparency

The gridded, shaded, coloured plot produced above is a much better representation of the two intersecting surfaces than the output of the standard IDL surface plotting commands. However, it still suffers from the problem that the viewer has no impression of the shape of the hidden parts of the surfaces. It turns out that we can use a trick to give the impression of transparency to the frontmost surfaces without incurring the computational or programming overhead of a full transparency algorithm.

The idea is very simple. If a surface is semi-transparent we can partly see what is behind it, so to give the impression of transparency all we have to do is make a separate plot of the hidden surfaces and partially add them to the image of the visible surfaces. In RGB colour space the addition is trivial, and if we divide the image of the hidden surfaces by a suitable factor before performing the addition we can control the degree of transparency too.

We don't even need to calculate which parts of the individual surfaces are visible and which are not. The parts of any single surface that are visible on the combined plot will look the same when only that surface is plotted, so if we add an image of the complete single surface to that of the two intersecting surfaces the visible part remains the same and a partial image of the hidden part appears.

If we add images of both individual surfaces to the combined plot we get our final image, where both surfaces are gridded, coloured, shaded and semi-transparent.

Note that this is not true transparency. For each individual surface only the non-visible parts of the other surface are made visible but not its own hidden parts: for example, the lines on the back of the large central peak are not plotted at all. However, the plot undoubtedly conveys more information than the non-transparent version, and so serves a useful purpose even if not perfect.

Implementation

We already have an image of the combined surface, shade_col_grid, created above, so our first step is to create individual plots of the single surfaces using the techniques described in the single surface tutorial:

; Create a gridded, shaded, coloured plot of the first surface

    set_plot, 'z'
    device, set_resolution=[xwin, ywin]
    
    surface, surfone                                  ; create grid image
    onegrid = tvrd()
    
    shade_surf, surfone                               ; create shaded image
    oneshades = tvrd()
    
    loadct, 25                                        ; create coloured image
    shade_surf, surfone, shades=surfone_col
    onepic = long(tvrd_24())
    
    for i=0,2 do onepic(i,*,*) = onepic(i,*,*) * oneshades       ; combine
    onepic = reform(byte(onepic/255), 3, xwin*ywin)              ; them to
    onepic(*,where(onegrid ne 0)) = 20b                          ; form the
    onepic = reform(onepic, 3, xwin, ywin, /overwrite)           ; composite
    
    
; Repeat for the second surface

    surface, surftwo                                  ; create grid image
    twogrid = tvrd()
    
    shade_surf, surftwo                               ; create shaded image
    twoshades = tvrd()
    
    shade_surf, surftwo, shades=surftwo_col           ; create coloured image
    twopic = long(tvrd_24())
    
    for i=0,2 do twopic(i,*,*) = twopic(i,*,*) * twoshades       ; combine 
    twopic = reform(byte(twopic/255), 3, xwin*ywin)              ; them to
    twopic(*,where(twogrid ne 0)) = 230b                         ; form the
    twopic = reform(twopic, 3, xwin, ywin, /overwrite)           ; composite

Then all that remains is to add the three images together in the right proportions and plot on the screen device.

    transp = 1.0/3.0    
    finalpic = (1.0-transp)*shade_col_grid + transp*onepic + transp*twopic
    finalpic = byte(finalpic/(1.0+transp))
   
    set_plot, screen
    tv_24, finalpic

Notice that altering the value of transp allows us to vary the transparency of the surfaces. It would of course be possible to have different degrees of transparency for the two surfaces by using different transp values for each picture in the final summation, but care needs to be taken that the total contribution to all the visible portions is the same.

[ Lund University ] [ Fysicum ] [ Maxlab ] [ Sljus ] [ STM Group ]

Maintained by: Sljus Webmaster
970220