Slashdot Mirror
Writing Code for Surface Plots?
MySchizoBuddy asks: "In what ways can you code plots of surface charts using a z(x,y) function or a cloud of points? I'm looking for a tutorial that explains this that doesn't use DirectX or OpenGL libraries (the language I'm using cannot use either framework anyway). How is the 3D mess generated and how can the 2D contour plots be generated as well? I'm assuming once I know that I can also use it to make torus plots as well. Remember, I'm asking for the explanation of the underlying math and an example code that does that. The GNUPlot gallery has some examples that I find helpful, but are there similar examples out there? (Remember, I am writing the plotting code as well)? Can anyone help or point me in the right direction?"
first plot
...
hum
http://www.detachedsolutions.com/graph3/
This is software that does what you need, written for the TI-83 (should work fine on an 84 too). Since it's open source, you can peek at the source code to get an idea of how its 3d rendering works.
Unfortunately, most TI calculator apps are coded in assembly...
You could plot graphics that represent levels as colours. I think I invented a kind of R^2->R function plot that represents a Z value as the degree of grey (0 being black, 255 being white) corresponding to Z mod 255. The result is that the more frequent "cuts" (sudden changes from white to black) are, the highest is the derivative.
Looky here, here and here.
I generated solutions to partial differential equations in Mathematica that looked similiar to those plots. You can animate them (Bessel Functions, Heat Equation, and harmonic motion). I personally would use Mathematica. If you're a student you can get that application for cheap.
Yes! I listen to NYC Speedcore and do math at 3AM. I suggest you try it too.
I looked at 3D rendering quite a while back (about 15 years ago). Believe me -- you probably don't really want to built this up from scratch.
It sounds like you just need some sample code that doesn't rely on OpenGL or DirectX. To that end, the following (open-source) projects may be useful:
Gan Family Homepage
Meshes ? OMG data structures to manage lines. My brain just froze.
So what is you're API that lets you turn pixels on and off? Its not that hard to go from a cloud of points in openGL to a cloud of points in another api.
Contour lines can be harder. you need to connect points on the same elevation without drawing a line through higher/lower elevation. I'd say don't bother connecting lines. just draw a lot of points (of certain colors for certain heights) with a fine enough granularity that they bleed into each other.
http://en.wikipedia.org/wiki/Jury_nullification
just a thought, but why not use MATLAB. You can sometimes use it to make stand-alone programs and maybe even DLLs. then add it to your code. I help develop code for a Cray X1 and all visualization (for my project) is done locally using MATLAB. It's pretty beefy, but it'll make a contour plot.
887321 = 337*2633
Good 2D contour plotting is not easy (I know!!!).
If you are not afraid of coding, then the best reference that I know of is:
by Albrecht Preusser, ACM TOMS V10#4
The code (in FORTRAN, unfortunately) is also in TOMS as:
Algorithm 626 TRICP: A Contour Plot Program Triangular Meshes
You can get ACM TOMS (Transaction on Mathematical Software) from any reasonable college library.
This might be a bit of overkill, as it is designed for irregular data sets (versus regular datasets - grids). However, the code can be separated into a separate Delaunay Triangulation step. If you have a grid (or in your case, a function f(x,y)), it is easy to generate fixed triangles.
If anybody knows of a better algorithm, I'm all ears!!!
I know this doesn't help you, but I'm still curious -- under what circumstances can you not use OpenGL for this? Under what circumstances would it be easier to implement it yourself than to find/create a decent GL binding? Or maybe you could just read the software implementation from Mesa?
Don't thank God, thank a doctor!
Check out the R-project at www.r-project.org. It is very good for that sort of thing.
I wrote a program to do that in 1968, for a UNIVAC 1108 with a pen plotter. Even did hidden surface elimination. We even ran it on a display a few times during really slow periods. But I don't have the ALGOL code any more.
The algorithm was in JACM around 1967, but you wouldn't do it that way today.
My advice would be to pick up a decent intro book on graphics. The Foley and Van Dam is a classic, but nearly any good book on the topic tends to have a review of the basic math involved in 3d: linear algebra with 3- or 4-value vectors and 3x3 or 4x4 matrices, the perspective divide, etc. Knowing the basics will help you have a lot better idea of what to do. But for anyone curious about the most basic 3d, here's a simple way to do it that should work decently for simple wireframe plots.
" > /> <!-- First line to draw --> /> <!-- Second line to draw --> /> <!-- Third line to draw, etc... -->
1) Take each 3d point of interest and apply some set of rotations around the X, Y, and Z axis to get the view angle right. You'll probably want to just do one rotation around each axis. e.g.:
x2 = x * cos( z_angle ) + y * sin( z_angle );
y2 = x * -sin( z_angle ) + y * cos( z_angle );
x3 = x2 * cos( y_angle ) + z * sin( y_angle );
z2 = x2 * -sin( y_angle ) + z * cos( y_angle );
y3 = y2 * cos( x_angle ) + z2 * sin( x_angle );
z3 = y2 * -sin( x_angle ) + z2 * cos( x_angle );
At this point, you'll have rotated your original coordinate (x,y,z) by some angle around some axis to get (x3,y3,z3)
(To my fellow graphics geeks: yes, yes, I know full well that matrices and quaternions are better -- I'm trying keep this simple.)
2) Maybe add some distance to the z value, to move ("translate") it away from the "camera" at the origin. You'll need to add enough to make sure that all of your z-values are positive and non-zero before step 3:
z4 = z3 + z_translate
3) Do what's called the "perspective divide" to find the two dimensional coordinate of this point. e.g.:
px = x3 / z4 * scale + center_x;
py = y3 / z4 * scale + center_y;
4) Either plot a point at the (px,py) that you just got, or compute steps 1 and 2 twice for the end points of lines and then use your system's graphics primitives to draw lines between (px1,py1) and (px2,py2). If you're not using a language or a library that can do simple 2d graphics, then I'd suggest writing out a file with the lines in SVG format (and then rasterize to bitmap with a program like Batik, Inkscape, Illustrator, etc.) e.g.:
<?xml version="1.0" standalone="no"?>
<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd
<svg width="4cm" height="4cm" viewBox="0 0 400 400" xmlns="http://www.w3.org/2000/svg" version="1.1">
<path d="M px1 py1 L px2 py2" fill="none" stroke="black" stroke-width="1"
<path d="M px1 py1 L px2 py2" fill="none" stroke="black" stroke-width="1"
<path d="M px1 py1 L px2 py2" fill="none" stroke="black" stroke-width="1"
</svg>
5) If you want more than wireframe, things will be a lot more complicated. At the very least, you'll need to be able to fill polygons and sort back to front by the z-value after the rotation so that everything gets drawn in the correct order. Z-buffers would be even better.
Well, I hope all that helps. Again, if anyone's really interested, I'd suggest getting a real book on the topic. But I hope this's enough to get the curious started.
dim ymax(screenwidth),ymin(screenwidth)
For xtemp = xlower to xupper step xincr
for ytemp = ylower to yupper step yincr
compute z=f(xtemp,ytemp)
(x,y)=transform(xtemp,ytemp,z)
if yymax(x) then plot x,y : ymax(x) = y
next ytemp, xtemp
Adjust your viewpoint for your 3D to 2D matrix so it plots from front to back.
Algorithm comes from an old book on doing 3D on an Apple ][ - it wasn't too hard to convert it to work on a TRS-80 Model I or a CoCo.
I just looked at the GNUPlot page, and regardless of what anyone else says, I don't think it should be too difficult to code at all, at least the grid stuff like near the top of that page. The reason - it's all drawn without perspective, which means a little vector math tells you where to put each point. You might be a bit foolish to try to hand-code something substantial with perspective, but orthogonal projections are straightforward.
Here's what you need to do. First, figure out what a unit in the X, Y, and Z axes corresponds to in screen coordinates. In the examples on the GNUPlot page, I'd estimate that a unit in the X direction is about 10 pixels to the right and 2 pixels down in screen coordinates, or (10, -2). Similarly, a unit in the Y direction is about (3, 3), and a unit in the Z direction is (0, 1). Now, to find the screen location of a point given in 3D, say (x,y,z), we start at the screen coordinate of the origin (wherever you want to put (0,0,0)). To this, we add x times the vector (10, -2), y times the vector (3, 3), and z times the vector (0, 1). For example, the point (10, 3, 6), with an origin at (100, 100), would show up at (100, 100) + (100, -20) + (9, 9) + (0, 6) = (209, 95).
In other words, if we call the unit vectors in the X, Y and Z directions (given in screen coordinates) Xhat, Yhat, and Zhat, then the point (x,y,z) should be plotted at (origin + x*Xhat+y*Yhat+z*Zhat).
I presume that from here you should be able to figure out how to draw the grid lines - just sample a bunch of points, and draw straight lines between neighbors. If you want to fill in the tiles (like the pictures near the bottom of the GNUPlot page), you'll need to be clever about the order you fill them in so that the ones in front are filled last. For this, you can sort the tiles by the appropriate combination of x and y coordinates (in "real" space, not screen space) - in our example, Xhat - 3.333Yhat points straight down, so you would sort by the combination x - 3.333y and draw in that order.
This should at least get you on your way. I'm sure there are subtleties to consider, but hey - it's your project, not mine, so I'll leave it at that!
Hate to disappoint you, but you're probably not going to find a tutorial and code that conveniently meets your requirements. For starters, eliminating OpenGL and D3D is going to remove a large number of tutorials, simply because it doesn't make sense for most people to roll their own rendering code.
Your best bet would be to look through "Computer Graphics: Principles and Practice" by Foley, van Dam, Feiner and Hughes. It covers a large number of topics, in enough detail to get a good idea of how to get started and where to look for more information. Also "Real-Time Rendering" by Tomas Möller and Eric Haines. Despite the name, a large portion of the book is applicable to almost all computer graphics, not just real-time. Also, a book on linear algebra would probably be helpful.
Maybe not
The actual code for such a plot, pixel cloud or such, is easy. 5-10 lines, nicely covered elsewhere in the comments already. Most difficulty lies in putting the first pixel on the screen, the right way.
In times of Amiga, it would be like typing "ScreenOpen 320,256,32", then a line defining the palette in equally simple terms, then you can start the drawing stuff. On PCs there are dozens of graphics libraries and each one more difficult to use than another. Microsoft Visual someshit have this pretty easy to make in a crappy way, not too difficult in a bit harder way. SVGALib, currently considered obsolete, is quite friendly, SDL is a bit of a bitch to initialize, though both require you to roll your own gfx function libraries (another 10 lines of code like #define plot(x,y,c) SCREEN[(y)*WIDTH+(x)]=(c) ). Stay far away from QT, GTK and alikes. They are a true hell to draw arbitrary stuff. ImageMagick is very nice for generating gfx files, but sucks at displaying stuff. Matlab makes ploting a breese, but that's not really a product you can ship to a customer or deploy in production environment. You could actually take GNUplot sources and plug your piece into them, or use some sockets or scripts to launch GNUplot. Or take one of the ancient languages like Turbo C, Turbo Pascal, Fortran, and write that program in them, it's much easier.
Anagram("United States of America") == "Dine out, taste a Mac, fries"
As part of some final year course thing I wrote a 3D rendering engine in Java based on implicit surface rendering, I did it in two stages. The first was simply an engine that converted the surface equation (f(x) + f(y) + f(z) = whatever) and then plotted it everything into a spacial array then dumpped everything into java's 3d library for scaling and rotation. Then I knocked out the Java 3D API and replaced it with my own so that everything worked all the way back to java 1 with no extra libraries. It was God aweful slow but it did everything, translation, rotation, scaling, specular highlighting, Phong and Gourond(sp?) shading as I remember.
I could probably dig that out.
Regards, Phil
The idea behind plotting z = f(x,y) is to sample the function on a grid of points (x,y). You then interpolate linearly between adjacent points to get a polygon. For example, if my grid contains (x,y) = {(0,0), (0,1), (1,0), (1,1)} then my polygon in 3-space has corners (0,0,f(0,0)), (0,1,f(0,1)), etc. Project this polygon onto the screen using the standard 3d transformations.
There are a couple of things to watch out for. First, all polygons should be convex. This is usually not a problem, but you can guarantee it by making sure that the 2D grid is composed of convex polygons - typically rectangles or triangles are used. It might be best to use a grid of rectangles first because that is easier to understand and debug, but when that is working move to a grid of triangles because everything will be much faster.
Second, if the function changes rapidly between the grid points then your plot will probably be wrong. One solution which works in most cases is to use an adaptive grid. For each polygon in 3D, average all the corner points to obtain the center. Compute the actual function value using the x and y coordinates of the center, and compare with the center's z-value. If they differ significantly then you split the polygon up, using the center point (with the real function value) as a new vertex, and repeat. A rectangle will split into 4 smaller rectangles, a triangle into 3 smaller triangles, etc. To determine if the values differ significantly, project them onto the view screen and see if they map to the same pixel. There are ways to optimize all of this significantly.
For plotting a cloud of points you are really screwed unless you know the function is of the form z = f(x,y). In that case, each (x,y) point has two closest neighbours which should give you a triangle. Sort the values by x and then by y, and as you construct each triangle remove the initial point (but not the neighbours) from the list. This will give you a triangular mesh. There are a few weird things which can happen, but you can mostly solve these problems by throwing away bad points.
In Soviet America the banks rob you!
If its a very small part of a larger program, then look for some good examples and try to trudge through getting it to work. Otherwise, start writing code and examine the visual output of the screen. Look for algorithms and methods for producing the output you want. When you begin to understand them, you can integrate the concepts into your solution.
Unfortunately, even if you have an API, you still have to understand what it does and how to interact with it. Many people believe that modular programming simply means using code alrealy written. This is a false premise as there is no substitute for some degree of understanding.
In the end, you may wish to use a public library because it is highly accurate, optimized, well researched, or portable, but by that time you will have the knowledge to understand why you are selecting the particular library.
I actually understand your question and the first place to look is a book Statistics and Data Analysis in Geology" by Davis. I used it in graduate school back in the 80's. It has understandable code for generating this output. It starts with line printer plots to get you running, first with trend fitting the data, and then actually machine contouring the data. I used their examples to program my C64 do draw contour maps on the 1526 plotter. I have tried surface plots with Excel off and on and it is frustrating because the mindset of the Excel developers is in a different universe. Google around for COGS Computer Oriented Geological Society, they may have some of the source as FORTRAN in their archives.
... (0,n) and then start with the next row, (1,0), (1,1),(1,2)...(1,N). Keep going down the rows.
If you have $$$$.00 to through at the problem use IDL from Research Systems Inc. They just got bought by someone so you will have to Google for their new digs. Expect to pay three kilo-bucks.
Gridding is a cool problem. Basically, you start at (0,0) and calculate the distance of each point to (0,0) and then sort by that distance. Then you take the nearest every how many points and fit a trend and get the value at (0,0),. next you do (0,1) and calculate distances, sort, estimate the value, move to (0,2)
Next you have to decide bottom contour line value, top contour line value, and contour interval. Code as if you are a pen plotter. Look at the values of each corner of each grid cell, one at a time, and if a contour line goes through that side connect it to where it crosses the other side. This requires a lot of logic that is beyond this scope but reply to me and we can get deeper if you need.
The Commodore 64 would grid all weekend, but did a decent job of making plots and maps. Core was cheep back then.
One thng I can think of is for each Z value chosen, to make lines between adjacent points and test if any point on that line segment crosses that Z. It's pretty iterative but you'll get your point cloud. That's relatively straightforward, provided you have a fairly regular series of points.
Another way was to use vector math, make a triangle between three terrain points and try to see if it intersects with your Z plane. Maybe ten years ago I made an algorithm that returns a line segment if there was such an intersection.
If you wanted to do curve fitting and see where that curve passes your z plane, that could be done, but that's beyond my education.
R is OSS and has the commands he wants (image, countour, and wireframe) it sounds like exactly what he wants: code that makes these plots that he can look at.
You can use the marching cubes algorithm. You can go from a cloud of points to a 3d object that you can transform like you want. I think they use that in MRI imaging. I used the 2d version (marching squares) to optimize a metablob routine a long time ago.
see: http://www.exaflop.org/docs/marchcubes/index.html
Generic Mapping Tool, will do everything you want, and a lot you have not even thought of yet. GMT produces postscript and encapsuleated postscript, so if you have a postscript viewer your good to go, runs on unix/linux, is GPL'd so there source code to analyse. There are lots of examples with generating scripts such as time-series collected along a track or a 3D perspective, artificialy illuminated, greyscale image; just reading and understanding the accompanying ducumentation should be good for college credit in cartography.
Apocalypse Cancelled, Sorry, No Ticket Refunds
- the 3d to 2d mapping from a viewpoint
- dealing with hidden points
The 3d to 2d mapping is actually pretty easy, with a little Linear Algebra. If you don't know any, grab a Linear algebra textbook and read up on vectors, dot products, cross products, etc.You pick a viewpoint V = (x_v, y_v, z_v) and a direction vector D = (d_x, d_y, d_z) ( D is a unit vector, it's length should be 1). and an "up" vector U (which way is up, also a unit vector, perpendicular to D). Next we need a "right" vector, R, which is always D x U (vector cross product) So now our mapped coordinates for any location P are:
z = (V - P)*D [where * is the dot product]
x = ((V - P)*R) / z
y = ((V - P)*U) / z
So next we need to deal with hidden lines. There's an old trick for this, make a list of line endpoints, and compute z, above, for the center of each line. Sort the lines by increasing z. Now start at the first positive z value, and draw each line that would not cross an existing line.
- "History shows again and again how nature points out the folly of men" -- Blue Oyster Cult, 'Godzilla'
http://www.pacifict.com/Products.html
GCHQ Quantum Insert installed. If only our tongues were made of glass, how much more careful we would be when we speak
The standard way of doing this is to use the Marching Cubes algorithm. Google it --> http://www.exaflop.org/docs/marchcubes/
:wq ~ ~ ~ ~ ~
Pardon me if I missed seeing a reference to the PLOT10 library in this discussion. Used to be THE lib for TEKTRONIX terminals (line perinter then color CRT) for graphics anything. Wrote some stuff for plots in C and fortran for data surfaces way back (Z80 CP/M, and later IBM PC), as it was ported as a general purpose lib, to graph on a flat bed pen plotter and CRT (Hercules monochrome graphics adapter). Also, I seem to recall a few geeks simulation Rubik's Cube and Enterprise (Dammit Jim, the Original Series!) on a Tektronix Color CRT when I was a wee lad in college (pre IBM PC). You should be able to google for the source or for a compatible port.
m l
http://www.gaeinc.com/plot10.htm
http://www.cvrti.utah.edu/~macleod/docs/pscont.ht
http://www.ill.fr/Computing/pgplotSS.html
and then there's http://plplot.sourceforge.net/ as well. (includs the classic "Lena" demo pic)
I've got some code here somwhere's....hmm.....now what did I do with those 8 inch floppies and drive.....
--
"The most exciting phrase to hear in science, the one that heralds new discoveries, is not 'Eureka!' (I found it!) but 'That's funny...'" -Isaac Asimov
Foley and Van Dam have lots of info: http://www.amazon.com/gp/product/0201609215/sr=8-1 /qid=1153828209/ref=sr_1_1/103-0593733-9387027?ie= UTF8/
You could look into Octave, which is an OS version of MATLAB released under GPL. You could either use it directly or look into the source for some helpful ideas.