A: 

I always check places like http://mathworld.wolfram.com before going to deep on my own :)

Maybe their curves section would help? Or maybe the entry on maps.

warren
A: 

compare what you have rendered with a real-world topo map - they look identical to me! i wouldn't change a thing...

Steven A. Lowe
Well, the *positions* of the lines are pretty good right now. But the lines themselves are all broken-up and muddy and sloppy looking. Drawing the topo lines as (...cont'd...)
benjismith
...bezier curves would (I think) provide a very nice cleanup. But it'd also give some bonus benefits, like making elevation regions selectable with a mouse click (which would be very very cool!)
benjismith
+6  A: 

The gradient is a mathematical operator that may help you.

If you can turn your interpolation into a differentiable function, the gradient of the height will always point in the direction of steepest ascent. All curves of equal height are perpendicular to the gradient of height evaluated at that point.

Your idea about starting from the highest point is sensible, but might miss features if there is more than one local maximum.

I'd suggest 1. pick height values at which you will draw lines 2. create a bunch of points on a fine, regularly spaced grid, then walk each point in small steps in the gradient direction towards the nearest height at which you want to draw a line 3. create curves by stepping each point perpendicular to the gradient; eliminate excess points by killing a point when another curve comes too close to it-- but to avoid destroying the center of hourglass like figures, you might need to check the angle between the oriented vector perpendicular to the gradient for both of the points. (When I say oriented, I mean make sure that the angle between the gradient and the perpendicular value you calculate is always 90 degrees in the same direction.)

ellisbben
Very nice. See my most recent edit for additional clarifying questions.
benjismith
+3  A: 

Alternately, there is the marching squares algorithm which seems appropriate to your problem, although you may want to smooth the results if you use a coarse grid.

The topo curves you want to draw are isosurfaces of a scalar field over 2 dimensions. For isosurfaces in 3 dimensions, there is the marching cubes algorithm.

ellisbben
+2  A: 

I've wanted something like this myself, but haven't found a vector-based solution.

A raster-based solution isn't that bad, though, especially if your data is raster-based. If your data is vector-based too (in other words, you have a 3D model of your surface), you should be able to do some real math to find the intersection curves with horizontal planes at varying elevations.

For a raster-based approach, I look at each pair of neighboring pixels. If one is above a contour level, and one is below, obviously a contour line runs between them. The trick I used to anti-alias the contour line is to mix the contour line color into both pixels, proportional to their closeness to the idealized contour line.

Maybe some examples will help. Suppose that the current pixel is at an "elevation" of 12 ft, a neighbor is at an elevation of 8 ft, and contour lines are every 10 ft. Then, there is a contour line half way between; paint the current pixel with the contour line color at 50% opacity. Another pixel is at 11 feet and has a neighbor at 6 feet. Color the current pixel at 80% opacity.

alpha = (contour - neighbor) / (current - neighbor)

Unfortunately, I don't have the code handy, and there might have been a bit more to it (I vaguely recall looking at diagonal neighbors too, and adjusting by sqrt(2) / 2). I hope this enough to give you the gist.

erickson
I don't have a 3D model, but any surface in my system is representable as an equation, so it certainly qualifies as an arbitrary-precision vector representation. But the equation is biiiiiiig and nasty. For example, the image on this page uses an equation with more than 300 terms. It's waaaaay ugly.
benjismith
Yes, I remember your other question about prepping your data set for interpolation. You could use that to do something more accurate than linear interpolation, but since your target output is still one pixel, the difference may not be visible.
erickson
I like your approach of determining whether the contour line runs *between* two pixels. The approach that I took in my rendering (above) was to see whether the pixels value was *close* to the real contour line, given some arbitrary epsilon. (...continued...)
benjismith
I think your technique would yield better results though at the expense of 4x increase in processing time. On the other hand, you could combine the two approaches, only performing the multi-pixel routine if the pixel's value was within epsilon of the contour's elevation. Hmmmmmmmm.....
benjismith
Of course, having said all that, I still want vector splines!!!
benjismith
Yes, I still would like vector splines too. I'm actually trying to make orienteering maps for my Boy Scouts, and would love to develop an SVG output.
erickson
+3  A: 

In response to your comment to @erickson and to answer the point about calculating the gradient of your function. Instead of calculating the derivatives of your 300 term function you could do a numeric differentiation as follows.

Given a point [x,y] in your image you could calculate the gradient (direction of steepest decent)

g={  ( f(x+dx,y)-f(x-dx,y) )/(2*dx), 
  {  ( f(x,y+dy)-f(x,y-dy) )/(2*dy)

where dx and dy could be the spacing in your grid. The contour line will run perpendicular to the gradient. So, to get the contour direction, c, we can multiply g=[v,w] by matrix, A=[0 -1, 1 0] giving

c = [-w,v]
Azim
Yeah, this was my first thought, to just get an empirical measurement of the x-slope and y-slope. But I didn't know how to transform that into the gradient. So I started working on calculating the partial derivatives for the whole huge function. Thanks for the matrix multiplication trick!
benjismith
your welcome. glad it helped.
Azim
+1  A: 

It occurred to me that what you're trying to do would be pretty easy to do in MATLAB, using the contour function. Doing things like making low-density approximations to your contours can probably be done with some fairly simple post-processing of the contours.

Fortunately, GNU Octave, a MATLAB clone, has implementations of the various contour plotting functions. You could look at that code for an algorithm and implementation that's almost certainly mathematically sound. Or, you might just be able to offload the processing to Octave. Check out the page on interfacing with other languages to see if that would be easier.

Disclosure: I haven't used Octave very much, and I haven't actually tested it's contour plotting. However, from my experience with MATLAB, I can say that it will give you almost everything you're asking for in just a few lines of code, provided you get your data into MATLAB.

Also, congratulations on making a very VanGough-esque slopefield plot.

A: 

Have you tried these ?

http://local.wasp.uwa.edu.au/~pbourke/papers/conrec/

and www.galiander.ca/quikgrid/index.html

anon