Original Article Date: Feb 1, 2001

MSM featured this article authored by Dave Brumbaugh which discusses generating contours from survey data and even supplied the necessary code to accomplish this task. The article has now found a new home and reads.....
Originally published in MSM magazine February 2001

Note: The content of this article applies to MicroStaiton /J



My August 1999 article about generating surface models from USGS data and my May 2000 article about generating artificial terrain-like surfaces have generated a lot of E-mail. Most of the messages take the form: “My data isn’t nicely gridded, like USGS data, and I’m not very interested in creating an imaginary surface. But can you help with contours or other data my surveyor is providing me?” Well, here goes…

I’ve added the ability to generate surfaces from contour data to my original landscape application from last May. Actually, your data can be 3D scattered points, connected line strings, curves, ellipses or even shapes. When you key in JSURFACE FROM CONTOURS, you’ll be prompted to fence an area and the surface generated will be based on all the points read from all the elements found in the fence. Your working view “Z” direction is what will be used to calculate the resulting surface. As with the previous application, JSURFACE FROM CONTOURS BSPLINE will generate a B-spline surface and JSURFACE FROM CONTOURS LINESTRINGS will generate a linestring profile representation of the surface.

1. The figure above shows a triangulated surface generated from spot elevations (recorded as circles, level 20) found in the MicroStation example file Bentley\Workspace\projects\ examples\civil\dgn\survey.dgn. You can compare my triangulated surface with a more traditional digital terrain model as found in dtm.dgn in the same directory.

1. The above figure shows the contours in plan of the Old Faithful quad sheet available for free download from http://www.topodepot.com/. I downloaded this “24K Contour Sample Old Faithful, WY” in DWG format. I translated it to MicroStation format using just the default settings.

Inner Workings
The application uses a technique called “weighted interpolation.” You can refer to Paul Bourke’s article, Nearest Neighbor Weighted Interpolation at http://astronomy.swin.edu.au/~pbourke/modelling/weightinterp/ for the relevant equations.

This technique superimposes a rectangular grid over a random sampling of points and assigns a Z value to the grid points by computing a weighted average of the height of all the points in the sample. The weighting is computed using the x,y distance between the sample points and the grid points. The closest points are strongly weighted and points that are further away are weighted using an inverse exponential function. The default setting for the application uses an inverse cubic function for the weighting of the points. Thus, the weighting of points even a short distance away falls off very rapidly. The drawback to this method is that calculating the weighted average for each grid point requires iterating through all the sample points. This can take a while, and the densest grid I had the patience to try was 75 x 75. Of course, if your sample is relatively sparse (e.g., a limited number of rough contours), this can calculate rather quickly, but there is no way to approach the very fine mesh I was able to use in the May article unless you’ve got a fast computer and a lot of time.

Because the calculation of the weighted average is so computation-intensive, I’ve split it off into a dynamic link library called interpolate.dll. The Visual C++ source code for the library is included in the download. The application thus has three parts: a JMDL application containing the user interface and the design file processing code, an MDL shared library for placing the B-spline surface (see the May article for details) and a C++ .dll for the intensive computations. You can see my February 2000 article about using jmdl_import to call dynamic link library functions from JMDL.

1. The same contours of Old Faithful in an isometric view looking toward the northeast.

 

Installing The New Landscape
Download the code.

This file includes LandscapeSource.zip, the source code, and LandscapeClasses.zip, the .mclass and .class files for the application, landscape.ma (the MDL shared library for B-spline placement) and interpolate.dll (the weighted interpolation dynamic link library). InterpolateSource.zip contains the source for the library. A Java .jar file, landscape.jar, is also included. The paths these files were installed in on my computer are included in the zip file. You should put them in the proper corresponding directories on your machine.

Put the .mjava, .java and .mke files in a source file directory, and put the .mclass and .class files in \Bentley\Program\MicroStation\jmdl\landscape\. If you compile the source files, the \landscape directory will be created for you. If you prefer, you can just put landscape.jar in a directory somewhere and add its full path to your CLASSPATH variable.

Put interpolate.dll in your MicroStation root directory and put landscape.ma in your \mdlapps directory. To load the application and use the commands, key in java landscape.Landscape. The letter case matters.

1. Triangulated surface generated from Old Faithful 24K contours using a 50 x 75 grid mesh.

 

Using the new LandscapeLoad the application with java landscape.Landscape (case matters) just like before. In addition to the PLACE LANDSCAPE command, the application now adds the command JSURFACE FROM CONTOURS with two modifiers, LINESTRINGS and BSPLINE. JSURFACE FROM CONTOURS with no modifiers will place a triangulated surface.

Keying in JSURFACE FROM CONTOURS opens up a settings dialog where you can specify the width and height of your grid, the weighting exponent (the default 6.0 works well and results in the inverse cubic weighting), a stroking tolerance (for extracting points from curves, arcs and ellipses) and the colors to apply to your triangulated surface. The color list will apply the last color number to the highest band of triangles, next to last to the next highest band, etc. There will be as many bands of triangles (each on a different level with the lowest band on level 62, next highest on level 61, etc.) as there are color numbers in your colors string. The default string works well with the default color table.

1. Rendered triangulated surface generated from Old Faithful 24K contours using a 50 x 75 grid mesh.

 

You’ll be prompted to place two corners of a fence block. The fence should contain the elements you want processed for your surface. Any shapes, linestrings, lines, point elements, curves, arcs, ellipses, complex strings or complex shapes will be processed to build a large array of points. The block you placed will define the limits of your surface and the extent of your grid points.

Key in JSURFACE FROM CONTOURS BSPLINE if you want B-spline surfaces generated instead of a triangulated surface. Key in JSURFACE FROM CONTOURS LINESTRINGS to place a line string representation of the surface. As opposed to PLACE LANDSCAPE, all options run in essentially equivalent time.

Building The Application
The key to compiling and building JMDL (and MDL) applications is having your environment set up properly. My environment is set up about like this with the following environment variables defined:

Where you see %MS%, this is expanded to include the string from your MS environment variable. Obviously, you need to set MS to point to your installation of MicroStation/J. You should just append %MS%\mdl\bin;%MS%\jmdl\bin; %MS% to your existing PATH variable rather then replacing it.

When you’ve got your environment set up properly, you should be able to change directories to the source directory, and type BMAKE –A LANDSCAPE.MKE at a system prompt. This should build the .mclass, .class, and .ma files and put them in %MS%\jmdl\landscape\ and %MS%\mdlapps\.

If you get errors when running BMAKE, don’t panic. The error messages are informative, and you should try to parse through them and figure out what is missing from your environment.

You may also want to look for the batch file, mstndevvars.bat, in ...\MicroStation\jmdl\bin. Running this file may set your environment appropriately.

To compile interpolate.dll, you’ll need Visual C++. Put the source in a working directory and open interpolate.dsw.

1. Rendered B-spline surface generated from Old Faithful 24K contours using a 50 x 75 grid mesh.

 

Development DetailsThe logic for the surface creation is contained in the class ContoursCmd (defined in ContoursCmd.mjava).

After the second data point is processed (see the dataPoint() method), the block we’ve been placing in dynamics is turned into a fence with a call to DgnKernel.session.setFenceFromShape(). With the fence in place, an iterator for all the elements in the fence is obtained from the method DgnKernel.session.getFenceIterator(). Note that the fence is set to process in inside clip mode. Note also that getFenceIterator() starts a new command that causes the current command’s stop() method to be called. I dealt with this using some global flags so my dialog box wouldn’t close.

Each element found in the fence is passed to the method scavengePoints(). Depending on the element type (note use of instanceof operator), the element is stroked and the element’s points are added to the class’s vector member variable m_vectorOfPoints.

When all the elements are processed, the x, y and z components of each point in the vector are added into the double array, xyzArray[]. This array is passed to one of the methods buildSurface(), buildBSplineSurface(), or buildLineStrings(), depending on the options you used to start the command. These three methods all have essentially the same structure. Using the points that defined the fence, the methods figure the width and height of the cells that make up the grid. Each x, y grid point is then passed with the double array of x, y, z values to getElevation(), which, in turn, calls interpolate_getWeightedElevation() from interpolate.dll. See the jmdl_import statement for the syntax required for calling this function.

interpolate_getWeightedElevation() performs the distance weighted averaging of all the points in the array to determine the height to use for the grid point. The functions buildSurface() and buildBSplineSurface() save the grid points (x, y and interpolated z value) into the DPoint array, pointsGrid[][]. This array is then used to place either the triangles or the B-spline surface. The B-spline surface is placed with a call to the function landscpe_writeBsplineSurface() in the MDL shared library landscpe.ma.

Hints and enhancements
The most obvious addition is to add processing logic for text nodes and text strings, since spot elevations are often recorded this way in design files. Ideally, you would want your contours to be in a reference file and your surfaces to be in another file. You’d have to do a little more work in the fence processing portion of the application to make that work.

I’ve added some logic that allows this application to be called from the command line after a fence has already been placed (see the start() and doProcessing() methods). This should allow the application to be called from macros and user commands. After placing the fence block, the syntax for calling the application is JSURFACE FROM CONTOURS[TRIANGLES|LINESTRINGS|BSPLINE][,width,height, power,stroke tolerance]. To place a 50 x 50 B-spline surface with inverse cubic distance interpolation and one unit stroke tolerance, place a fence and send the key-in: JSURFACE FROM CONTOURS BSPLINE, 50, 50, 6.0, 1.0 . If you don’t include a width, the default is 50. Height defaults to 50, interpolation constant defaults to 6.0 and stroke tolerance defaults to 1.0.

Enjoy, and make sure you try to compile this. I’m always happy to hear your feedback, particularly if you are having problems. I would also welcome suggestions for future applications articles.


AskInga Article #274