Copyright (C) 2005 by John
Coulthard. Permission is granted to copy, distribute
modify this document under the terms of the GNU Free Documentation
License, Version 1.1 or any later version published by the Free
Foundation; A copy of the license is included in the file fdl.txt
can be obtained from the Free Software Internet site at
Revised: June 22, 2005
I initially wrote QuikGrid as an exercise to learn C++ programming under Windows. It was implemented and has been maintained using Borland C++ compilers, most recently version 5.0B (1996). My primary source of information for the Windows programming environment (API) has been from the books published by Petzold. My use of classes is primarily to simplify memory management. The code does not make use of any Windows Object Oriented classes like the MS MFC or Borland's OWL. QuikGrid consists of approximately 37 separate source (.cpp) files with associated .h files, etc. Originally the program was called “Surface” and you will occasionally find references to surface in the code (e.g. surface.ico is the icon for QuikGrid).
The programming style varies as I have been learning C++ as I go along. My many years of experience with the FORTRAN language shows through. (The contouring and grid expansion routines are originally translations of FORTRAN code.) Also I tended to experiment with different techniques as I developed the code.
The contouring and grid expansion routines have been distributed as free software under the GNU Lesser GPL license. They are available from my web site www.perspectiveedge.com . That is the more appropriate source code to download if you are only interested in one of those routines. The code has been compiled and tested under Linux and I have a report that it compiles under MS VC.
quikgrid.cpp implements Winmain. The dialog boxes are implemented in dialogbx.cpp and files named dlg*.cpp. The routine paintcon.cpp is the core code for all drawing (think of “paint contours”, which is all it did originally). The core code for data input is in loaddata.cpp. Data output tends to be handled by a host of different files, for example, vrml.cpp for vrml output, datawrit.cpp to output scattered data points, dxfwrite.cpp to output DXF files, etc. Command files are all handled in cmdfile.cpp. The saving and restoring of preferences is all in prefer.cpp. There are lots of little support routines implemented in files like utilitys.cpp and grid.cpp. If I can't remember where a function is implemented I usually find it by doing a “grep” like search on the source files.
At startup WM_CREATE does initializations then posts an IDM_STARTUP message. IDM_STARTUP checks to see if there are any command line arguments. This will happen if QuikGrid is invoked from a DOS command of the form “QuikGrid input.* output.*” or is linked to from another program using such a command. This is also what it looks like to QuikGrid if a data file is “dragged and dropped” onto the QuikGrid Icon or if a file association is made between some file types and QuikGrid and a user double clicks on the data file. .
If there are no command line arguments QuikGrid goes into a normal Windows wait state.
If a command line argument is present QuikGrid assumes the argument is a file name and processes the file as if a person had started QuikGrid and then selected an “File... Input scattered data points...” or File... Process command file” menu operation. Which command is determined from the file Suffix as described in the documentation. Similarly if there is a second argument it is assumed to be a file name and QuikGrid generates a grid and does an appropriate output command. If a second argument is present QuikGrid then terminates. If only one argument is present it processes the input file, generates a grid if it was not a command file, then goes into a normal Windows wait state.
The storage of the input data is handled by the class scatdata.cpp . Speed of access to the x, y and z data points is one of the critical areas with respect to CPU usage. The points are stored as three arrays x, y and z, each of which is one block of contiguous memory. Other storage schemes seem to produce higher overhead. Initially the code allocates enough space for 16,000 data points. If it runs out it allocates new space which is 2.5 times bigger, copies the data over to the new space and releases the old. This will be repeated as often as necessary until all the data points are stored. The upper limit of 128 million is just a sanity check. It also maintains an array of pointers to optional comments attached to a data point and an array of type char which is used to store flags (Of which only one is used - namely the flag that says this data point should be ignored during grid generation).
The Scatdata class keeps track of the max and mins for x, y and z. It also maintains a value called zadjust which is used to asure that all the z coordinates are positive. Negative values are used to flag unevaluated grid intersections throughout the program. The x and y values are normalized by subtracting the x,y value of the first data point from all data points. The normalization values are stored as dxnormalize and dynormalize. They may be retrieved using the function LoadNormalization. This allows all the computations to proceed using single precision arithmetic. Any time an internal value is to be displayed it must be modified by zadjust, dxnormalize or dynormalize. Similarly numbers entered in a dialog box must also be modified appropriately. This was not implemented very neatly and I have never taken the time to clean it up. Be careful - it is a constant source of error - although one easily dealt with.
All the data input are read as strings and the conversion to floating point is done by QuikGrid, not by the IO routines. There is a story behind that! Years ago Microsoft released a version of Microsoft Plus for Windows 95 that broke the floating point input routines for programs compiled using the Borland Compiler. (The bug didn't affect MS C++ compiled programs). The fix was a service pack that was 10 times the size of QuikGrid and I rebelled at the thought of expecting my customers to have to download and install a service pack if they were using Microsoft Plus. Rather than spend money on the Microsoft Compiler (and being more than a little annoyed at MS) I elected to read the data as strings and convert them to floating point internally.
The focus of the display is always the generated grid so LoadData will create an appropriate grid and initialize it (see function InitializeGridingrid.cpp). If a grid has been loaded QuikGrid will create a dummy set of 2 scattered data points which spans the space. (Any set of data points that contains less than 3 points is considered a “dummy” set internally, but it will be set so that the max and min's are appropriate to the grid input so an attempt to display just the data set in cases like this will not fail).
QuikGrid handles a WM_PAINT call differently depending on whether the program is in 2d or 3d mode. In either case there are three steps
1. Determine the scaling
required for the image (Scale2dRectangle or Scale3dCube).
2. Centre the image on the screen (CentreDisplay)
3. Paint the screen (Paint2dSurface, Paint3dSurface).
All those entry points are in Paintcon.
QuikGrid runs in mapmode MM_ISOTROPIC. The data is scaled to fit into the Windows world of 0 to 32k and all scaling, zooming, panning etc. is handled by changing the Windows window and viewport (SetWindow* and SetViewport*) and redrawing the screen. I think of it as “vector” mode as opposed to a “bitmap” mode. The code only rarely has to dip down and deal with the world of bits (but it must occasionally to determine things like line thickness in the Isotropic world – which is not specified in bits, and so the thickness changes depending on the amount of zooming that is being done).
The job of Scale2dRectangle and Scale3dCube is to determine the maximum and minimum values for the image and set up the scaling parameters for the inline functions ScaleX and ScaleY (and UnScaleX and UnScaleY – their inverse). Those inline functions map the data from the QuikGrid floating point world into the Windows integer Isotropic world. This is a trivial exercise for Scale2dRectangle. The scaling is always determined by the grid, not the scattered data points, which is why a grid must always be present, even if only the scattered data points are to be displayed.
Things are a little more complex in the 3d world. Scale3dCube rotates the extremes of the 3d cube to determine the maximum extents of the rotated 2d image.
The value Margin controls the white space at the edge of the image. MaxScale is a constant set so that the image can nicely be both zoomed and unzoomed (set MaxScale too high and the image can only be zoomed in on). Once this is set up all you have to do is remember to scale everything via ScaleX or ScaleY before handing them to any Windows drawing function.
The idea here is that the image is probably being redrawn because the user has clicked somewhere in the image, indicating it should be redrawn, zoomed/unzoomed/panned with the mouse location to be the centre of the new image. CentreDisplay handles the calculation of the new Window and Viewport. By the time CentreDisplay is called the zooming and panning information has already been handled via Paintcon functions like ZoomViewDouble,ZoomIn,ZoomOut and PanView.
For example, when the user left clicks the mouse, case WM_LBUTTONDOWN is invoked in quikgrid.cpp. It sets the mapmode to MM_ISOTROPIC, determines the position of the mouse, converts it to the Windows Isotropic world, then calls ZoomViewDouble, which updates the QuikGrid internal parameters for the desired Window. It then invalidates the current display (done by function PictureChanged) which triggers the WM_PAINT call to redraw the screen. WM_PAINT calls Scale2dRectangle (or Scale3dCube) then calls CentreDisplay to set up the viewport and window, then finally calls Paint2dSurface or Paint3dSurface to do the actual drawing. For the most part you can ignore the fact that scaling and zooming is being done while you work in the Paint2dRectangle or Paint3dCube code.
There is a lot of code in Paintcon that is soley there to reduce
the time taken to paint the image. The code would be a fraction of
it's current size and much simpler to follow if the code just
the image in a straightforward fashion.
I reduce the time taken to
display the image in two main ways.
First the code attempts to remove lines and features which are
visible, i.e. not within the viewing window. The functions Visible,
VisibleY and InTheGrid in paintcon.cpp are used
for this. For example the routine that displays the 3d hidden
grid, Hatch3d (in hatch.cpp),
routines to determine if a grid rectangle is in the viewing window
doesn't draw it if it is not eligible to be viewed.
Secondly the code assembles contiguous line segments into Windows
polylines (See the code in drawstuff.cpp in particular function
These two mechanisms greatly
reduce the number of calls across the Operating System Interface
reduce the time taken to draw an image, especially when a lot of
zooming is going
In other situations I pre-compute scaled and rotated values and save them for later use.
Life is relatively easy in the 2d world. You can take any QuikGrid x,y position, scale it with ScaleX and ScaleY and pass it to a Windows drawing API function like LineTo or MoveToExt and it will work fine. However most of the routines do not make heavy use of those functions. Instead they call internal QuikGrid functions like PlotTo and PolyBuffer which are defined in the module drawstuf.cpp.
PlotTo removes line segments that are not within the viewed window (i.e. Visible). It does not attempt to do any clipping though. If the line pierces the window or potentially passes through it the line is handed to Windows to clip. This decreases the time it takes to display the image, especially if any zooming is going on. The calling sequence is PlotTo(x,y,pen) where x,y are the scaled x,y coordinates and pen is zero for a move and 1 or 2 for a draw. A pen of 1 selects a normal black line and a pen of 2 selects a bold line. A negative pen value tells PlotTo to flush the polyline buffer.
Polybuffer assembles individual line segments together into Windows “polylines”. Polybuffer(x,y,draw) is primarily called by PlotTo. The arguments are the scaled x,y coordinates and draw is 0 for a move or 1 for a draw (PlotTo selects the bold face pen before calling PolyBuffer when a bold line is required). A negative draw tells PolyBuffer to flush the polyline buffer.
Drawing the grid: In an ideal world a grid line could be draw straight across from min to max in one draw. But we have undefined grid intersections that can force a blank segment into the line at any time. Function Draw2dGridline takes care of this complexity for non-coloured grid squares. Function Hatch2d from hatch.cpp is used to paint coloured grids.
Contour lines: Black and bold contour lines are only computed the first time they are defined for a given image. The generated and scaled line coordinates are saved in Save2dContours which is defined by the class SaveLineType. If the user zooms in on an image the precomputed and scaled coordinates are retrieved and passed on quite rapidly. 3d contour lines are handled similarly using Save3dContours.
In the 3d world a lot more has to be kept in mind.
The points all have to be
rotated and optionally a perspective projection applied.
The rotation math is all in rotate.cpp . When the viewer changes the angle of view RotateInitialize is called to set the rotation matrix constants. Thereafter any set of x,y,z coordinates can be rotated by calling the function rotate. Perspective projection is handled separately by function Project. Separating out the perspective projection saves cpu time if an orthographic projection only is wanted. (This was important in the early days).
Typically the z-axis is
scaled to be some fraction of the size of the x-axis (see zratio
in the documentation)
The grid is stored in Zgrid which is defined by the class SurfaceGrid. The implementation of this class is in surfgrid.cpp. The application of the zratio is done within the class. The scaling of the z grid coordinate is initially done by Scale3dCube through the function ApplyzRatio. Depending on what is going on the zratio scaling may be turned off for some computations, for example to obtain the true value of a grid coordinate to determine the colour of the grid square. (There is some real messy code in this area).
Hidden line removal is
achieved by drawing stuff the farthest away from the viewer
I believe this is called "the poor mans hidden surface algorithm". The code draws the xy axes first, and the z axis first or last depending on the angle of the view. The grid surface is drawn from the back to the front. This logic is all handled by function Hatch3d in hatch.cpp. Scattered data points and the outline are displayed last and are never hidden. If hidden contour lines are desired (the default) they are generated on a grid by grid basis by Hatch3d (which is slow!). Non hidden contour lines are rendered by CNTOUR in the same fashion as the 2d case. Non hidden contour lines draw very fast as polylines are generated. Labels can be applied if the 3d contour lines are not hidden.
The entire grid is rotated, scaled and stored in the variable xyGrid early on in the procedure using the Paintcon function RotateGrid. xyGrid is defined by the class xyGridClasswhich is located in xygrid.cpp .
QuikGrid can display every nth grid line (which I call grid resolution) or a section of a grid. This can greatly speed up the viewing of a large grid. This is a relatively recent innovation and so far I have only made use of the grid resolution feature. The nitty gritty is in the class SurfaceGrid and also in the module gridres.cpp.
Source files: XPAND.CPP,
void Xpand( SurfaceGrid &Zgrid, ScatData &RandomData)
Undefined grid intersections are flagged, by default, as -99999. QuikGrid adjusts the z coordinates so they are always positive and then corrects the values for display purposes. Contour makes use of this convention as well and does not contour negative areas of the grid.
XPAND is basically a Nearest Neighbour algorithm with eight sectors. The algorithm lends itself to an efficient implementation, which is why XPAND is fast.
Xpand is organized so that each grid intersection can be evaluated in a separate function call. This feature was included to allow a grid to be evaluated in the background under Windows 3.1. Under Windows 3.1 timesharing only worked if each program voluntarily relinquished control back to the OS after a small time period. With Xpand this was accomplished by only computing one Grid Intersection at a time. The functions involved include XpandInit - to initialize everything and XpandPoint, which evaluated one grid intersection only. It is very unlikely that these functions have any external use any more.
Routine Xpand calls XpandInit to set everything up and then calls XpandPoint for each grid intersection until the grid is evaluated.
Initialization is handled by the functions XpandInit and LocateInit. XpandInit initializes a lot of local variables and then calls LocateInit. LocateInit makes use of the Class GridXtype which sets up an array called Locator, which contains all the scattered data points and associates with each data point the grid intersection it is closest to. The array is then sorted by grid intersection. It then creates a lookup table that makes it very fast to find all the data points that are close to a given grid intersection.
The procedure is then identical for the calculation of each grid intersection.
First the closest data points in each octant are determined. This is done by looking at all the data points close to the grid intersection to be evaluated and then scanning the data points for nearby grid intersections by shelling out from the intersection being evaluated. The process will eventually stop due to:
- An edge of the grid is
- The distance from the grid intersection being evaluated becomes too large.
- The distance becomes far enough that any new data points can not contribute appreciably to the weighted value of the intersection.
Function PutInOctant keeps track of the closest data point in each octant.
Once the scanning process is complete the weighted average of the selected points are used to determine the value for the grid intersection. The grid intersection may be left unevaluated because there are no points nearby or because the points are all on one side of the grid intersection.
Speed of access to the scattered data points is the most critical part of the program. Storing the data points as single dimensioned contiguous memory arrays for X, Y and Z seemed to produce the best results for me. My tests showed that storing them as a structure of x,y,z was slower and storing them as any sort of chained blocks of memory (which I did for data sets larger than 16000 points under Windows 3.1) was much slower. All this was done back in the Windows 3.1 days. Perhaps the compilers are smarter optimizers now?
The rest of the time critical code pretty well cascades back from there to the table lookup for data points close to a given grid intersection (ScanOneGrid) and the shelling out process (SelectPoints).
The sorting of the array of data points and associated grid intersections is another time eater. Xpand uses the standard Unix math function Qsort for this purpose. I don't know what has been going on in the world of sorting these days but back in the 1970's the Quick Sort algorithm was about as good as you could get and the one provided in the library should be written in assembler and highly optimized (but maybe not?).
Much of this code dates back to 1993 when I was using a 25Mhz cpu with no floating point processor. I suppose you might say it was fine tuned for that type of platform.
Xpand suffers from the same artifacts that all grid generation schemes share. For example the contour lines may not "honour" the scattered data points (i.e. the contour line may go on the wrong side of the data point - this is because it is the generated grid being contoured, not the original data points).
The algorithm works best with scattered data points that are more or less evenly distributed. The most common artifact is typically a "ripple" in the generated grid that goes at 45 or 90 degree angles. This may become pronounced if your scattered data points tend to be oriented in rows or columns. The effect is generated as a column or row of data points moves from one octant to another as grid intersections are evaluated.
The algorithm is interpolative, not extrapolative. If you ask it to extrapolate to the edges when there are no points "out there" the results may be quite strange (play with the Edge Sensitivity and Distance cutoff to experiment with this - different data sets will behave in different ways).
If the data points are clustered, leaving large areas of the grid with no data points nearby, the grid generation times will suffer. By default these grid intersections may be left unevaluated but you may over-ride this by used the Distance cutoff parameter (and increasing the grid generation time). If the scattered data points are very sparse in comparison to the grid (a dense grid), grid generation times will suffer. Basically anything that causes Xpand to shell out a long way to find data points will cause the grid generation time to suffer.
Sourcefiles : CONTOUR.CPP,
void Contour( SurfaceGrid &Zgrid, float ContourValue )
Will contour a grid for a single contour value. Call it repeatedly to contour more than one contour value. The function void DoLineTo(float x, float y, int draw), which is defined in Paintcon, will be called repeatedly as the contour line is traced. Contour lines which intercept the edge of the grid are contoured first, then interior closed contour lines are traced.
CONTOUR assigns the average of the four corners of a grid to the centre of the grid and then contours the 4 triangles. There is only one way a contour line can traverse a triangle. Because each contour line is traced continuously from edge to edge or back on itself (in the case of closed contours) the lines can be “painted” using Windows polylines.