I was telling one of our developers about screen-space reflections (SSR) recently, and partway through I realized I had seen this before. Not in a GPU paper. In World Construction Set 2, the terrain renderer Gary Huber and I built for the Amiga in 1995.
Modern SSR is commonly associated with Crytek‘s published work around CryENGINE 3, particularly its 2011 presentation. The recipe goes like this. Render the scene once into specialized per-pixel buffers, a G-buffer, that hold color, depth and surface information for every pixel. Then, for each reflective pixel, compute a reflected ray and march it across the screen, checking it against the depth you already stored, until it hits something. Use the color of whatever it hit as the reflection.
That is what WCS 2 did to render reflective water, sixteen years earlier, on Amigas all the way down to a stock 7 MHz Motorola 68000 with no floating-point coprocessor.
So let’s find out whether the WCS 2 reflection code really is SSR, and prove when it existed. The source code is public, so I can show the actual algorithm. Dating it turned out to be the harder part (dating is hard, I hear. We can’t even use Tinder.)
What SSR is
The modern technique has four parts:
- A G-buffer. After the main render, each pixel has its color, its depth, its surface normal and some notion of how reflective it is.
- A reflected ray per pixel, computed from the view direction and the normal.
- A march along that ray in screen space, comparing the ray’s depth against the stored depth at each step. When the ray passes behind the stored surface, that’s a hit.
- A color fetch from the already-rendered frame at the hit, with some fallback (usually an environment map) when the ray leaves the screen.
It caught on because it needs no second render of the scene and no access to the scene geometry. Everything it needs is already sitting in the frame buffer. Its reflection-pass cost scales primarily with pixels and ray-march length, rather than scene geometry complexity.
The usual references are Tiago Sousa, Nickolay Kasyan and Nicolas Schulz’s “Secrets of CryENGINE 3 Graphics Technology” in the SIGGRAPH 2011 Advances in Real-Time Rendering course, where it appears as “Real Time Local Reflections”; Bart Wronski’s 2014 write-up, which says Crytek first mentioned it in the DirectX 11 update for Crysis 2; and McGuire and Mara’s “Efficient GPU Screen-Space Ray Tracing” (JCGT, 2014). McGuire’s companion blog post says “you can see this technique in action in Crysis 3 and Just Cause 2.” Just Cause 2 shipped in March 2010, though I couldn’t find anything from Avalanche describing their implementation.
Why we couldn’t just ray-trace
World Construction Set rendered real-world terrain from USGS Digital Elevation Models, with fractal detail added on top. Our target machines were Amigas with 6 to 8 MB of RAM. WCS 1 and 2 even ran on a plain 68000 without an FPU, doing all of their floating-point math in software. The 2.x SMakeFile has a noaccel target for exactly that case, linked against SAS/C‘s software IEEE math library, scmieee.lib. It worked. It was not fast. Our own feature sheet from the time says WCS “can render an unlimited scene size in about 6 megs of RAM since only a small portion of the terrain is operated on at any given time.”
That’s the rub right there. The WHOLE scene never existed in memory at once. The renderer generated one terrain cell at a time (four points, two triangles), rasterized it, and threw it away. By the time you wanted to know what a patch of water should reflect, the mountain it was reflecting was long gone.
A conventional ray tracer needs the geometry to trace against. KPT Bryce, which came out on the Mac in late 1994, was a real ray tracer over real geometry objects, and it could do reflective water that way. VistaPro and Scenery Animator, our direct competitors on the Amiga, had water but no reflections.
My weird idea
It started with a demo trick. European Amiga demos of the era faked reflections by drawing a squashed, upside-down copy of the sky onto the ground plane, sometimes just by having the Copper re-emit selected scanlines from the sky rows. I’d written one myself: a 16-level grayscale, 640 x 480 Amiga demo with a starry sky reflected in the ground plane, where the stars fell one at a time straight down to the horizon. It looked terrible in motion, because interlace timing made the stars jitter vertically. (Last year I rebuilt that star-fall effect, inefficiently, as a ShaderToy shader to explain something about shaders.)
That trick only works when the camera looks straight at the horizon. WCS cameras could be anywhere, pointed anywhere, so it was out. My next idea was to render the whole scene a second time from a camera flipped to sit under the water surface, looking up. That would have doubled the render time and only worked for a flat mirror, which wasn’t great. Then I chewed on rendering that looking-up reflected scene into an image and using a distortion map to push its pixels around according to the angle of the waves.
The weird-ideas realization was that we didn’t need a second pass at all. The only complete record of the scene was the rendered image itself, plus whatever we chose to save alongside each pixel. If we kept extra per-pixel buffers during the main render, a second pass could treat the finished image as a 2.5D height field and trace reflected rays through it. The trade-off was real: a separate reflection render would have captured things the main view never sees, like the undersides of things. But it would have cost a whole extra render. Gary and I talked through what extra data each pixel would need, I worried about the memory for those extra attributes and about how it would work with segmented renders, and Gary built it into the renderer.
We had reflections on the V2 plan by December 1994. A planning note from our archives, from our Questar Productions days (the 1996-era company name, before 3D Nature and before AlphaPixel), lists them between clouds and waves:
Stereo Rendering
Clouds
Waves
Reflections
2D Trees from scanned or painted imagesThe credits in the code describe the division of labor fairly. Version.c credits Gary with “Concept, initial implementation, algorithms, renderer, good ideas, principal programming” and me with “Interface style, initial GUI rework, optimization, weird ideas, additional programming.” Rendering reflections from the frame buffer definitely falls under “weird ideas.”
I remember being super excited when Gary got it working. The first renders were only so-so, because the blur, the ripple wiggle and the faux Fresnel falloff weren’t there yet. They looked like liquid metal. That looked cool, but Gary didn’t like it because it didn’t look like water, and he knew (or went back and worked out) why. He added the wiggle and the Fresnel-like falloff; I didn’t even know the Fresnel coefficient of water back then. The blur may have been my suggestion, since “what if we just mask-blurred it to get rid of the aliasing?” sounds like my kind of lazy fix, but I don’t trust that memory.
When I sent Gary a draft of this post, he said the idea might have been inspired by some SIGGRAPH article, though he can’t remember which, and that he did a lot of his daydreaming on yellow legal tablets at Boulder’s Dushanbe Teahouse. He gave me some of those notebooks years ago, and I saved them. They’re in a box somewhere. If they turn up, I’ll post what’s in them.
Show me the code
The legit Amiga V2 source for WCS is public on GitHub at AlphaPixel/3DNature, under the OSGPL. The links below are pinned to a specific commit so the line numbers stay valid. The complete text of the two main files is in the appendices.
| File | What it does for reflections |
|---|---|
GlobeMapSupport.c |
Reflection_Render(), lines 719 to 1023. The reflection pass. |
MapTopo.c |
The polygon rasterizer. Writes the per-pixel buffers at lines 1052 to 1059. |
GlobeMap.c |
Allocates the buffers (lines 191 to 207) and runs the reflection pass after terrain, foliage and clouds are done (lines 711 to 717). |
ColorBlends.c |
Sets reflectivity for water and foam (lines 245 to 262). |
One warning if you read the files on GitHub: they’re Latin-1, and a couple of the variables in Reflection_Render() are named æ and ß. GitHub shows them as garbage characters. The appendices here have been converted to UTF-8.
A specialized G-buffer, 1995 edition
When reflections were turned on, GlobeMap.c allocated three extra full-frame buffers next to the RGB image and the z-buffer:
if (settings.reflections)
{
SlopeMap = (float *)get_Memory(bmapsize * sizeof (float), MEMF_ANY);
ElevationMap = (float *)get_Memory(bmapsize * sizeof (float), MEMF_ANY);
if ((ReflectionMap = (UBYTE *)get_Memory(bmapsize, MEMF_ANY)) == NULL
|| ! ElevationMap || ! SlopeMap)
{
if (User_Message("Render Module", "Out of memory allocating Reflection buffer!\n\
Continue without Reflections?", "Continue|Cancel", "oc"))
settings.reflections = 0;The rasterizer in MapTopo.c filled them in at the same moment it wrote the z-buffer:
*(zbuf + zip) = QPt;
if (Reflections && ReflectionMap)
ReflectionMap[zip] = Reflections;
if (ElevationMap)
ElevationMap[zip] = ElPt;
if (SlopeMap)
SlopeMap[zip] = EloQY;So each pixel carried:
| Buffer | Type | Contents |
|---|---|---|
bitmap[0..2] |
3 x UBYTE | Red, green, blue |
zbuf |
float | Distance from the camera |
ElevationMap |
float | World elevation of the surface at that pixel |
SlopeMap |
float | Wave slope along the view direction |
ReflectionMap |
UBYTE | Reflectivity, 0 to 254. 255 marks foam, which never reflects. |
In modern terminology, this is effectively a specialized G-buffer: per-pixel depth, elevation, local slope and reflectivity stored alongside the finished color image. It isn’t a deferred-shading G-buffer in the textbook sense, since the color is already shaded and elevation and slope aren’t the usual channels, but it plays the same role for the reflection pass.
That’s 16 bytes per pixel. Our demo renders were 752 x 480, which works out to 5,775,360 bytes (5.78 MB, or 5.51 MiB) of buffers on a machine that might have 6 MB total. In practice, a 752 x 480 render with reflections turned on normally needed more than 8 MB of RAM once the terrain data, the image textures and the rest of the renderer were loaded. “Textures” in that era meant things like the raster image of the moon and a library of very small cutout sprites for trees and other vegetation. Aerial photo textures of the Google Earth kind generally weren’t available to us. That memory pressure is why there’s an out-of-memory prompt that offers to render without reflections, and why InitCloudMap() pages the reflection buffers out to disk while it builds the 3D cloud model, then reads them back. Reflections were also a per-render option. You could turn them off on any render to save memory and time, and plenty of people did.
The other trick for big images was segmented rendering. WCS could split a render into horizontal bands and render them one at a time, so only one band’s worth of framebuffer (RGB, Z and the extras) had to be in memory. Each finished band was written to disk, and once the whole scene was done we could free the renderer’s resources and stitch the bands into the full output image. Segments and reflections didn’t mix well. In the 1996 code the reflection pass runs inside each band, with only that band’s buffers loaded; we had no virtual memory to fall back on, so a reflection whose source was in the band above had nothing to find. We’d talked about exactly this when we designed the reflection pass, wondering whether vertical segments with a little overscan would fix it, but for some reason vertical segmenting didn’t arrive until much later. Later versions could cut the image into vertical strips (“Tiles in X”), which suits reflections much better, because a water reflection mostly looks straight up the screen and stays inside its own strip.
The reflectivity came from the water material. ColorBlends.c sets it to 254 * ReflectionStrength for water, where the strength is a user parameter that can be animated, and to the 255 sentinel for foam on breaking waves.
The reflection pass
Reflection_Render() runs once the whole frame is finished. Here is what it does, in order.
Scan every pixel. It loops over the frame and handles only pixels with nonzero reflectivity that aren’t foam (lines 740 to 744).
Rebuild the geometry from the buffers. There’s no scene to query, so everything comes from the stored values:
A = PARC_RNDR_MOTION(0) - ElevationMap[zip] * PARC_RNDR_MOTION(14) * .001;
D = sqrt(zbuf[zip] * zbuf[zip] - A * A);
if (D <= 0.0)
continue;
AoD = A / D;PARC_RNDR_MOTION(0) is the camera altitude and PARC_RNDR_MOTION(14) is vertical exaggeration, so A is how far the camera sits above this bit of water. zbuf holds the straight-line distance to the pixel, so D is the horizontal distance. A / D is the tangent of the angle at which the view ray hits the water. On a flat mirror, the reflected ray leaves at the same angle, so AoD is also the slope of the reflected ray.
Tilt the ray by the wave slope. If waves are on, the stored slope rotates the reflected ray. This is the law of reflection with a per-pixel surface angle:
if (SlopeMap)
{
dAodX = SlopeMap[zip] * zbuf[zip] / D;
ß = atan(dAodX);
æ = atan(AoD);
AoD = tan(ß * 2.0 + æ);
} /* if waves */Tilting a mirror by an angle turns the reflected ray by twice that angle, hence ß * 2.0. That’s just the law of reflection applied to a tilted facet.
A note on the names
Coming back to this code after thirty years, some of the names need a decoder ring.
| Name | Meaning |
|---|---|
æ, ß |
Stand-ins for the Greek letters α (alpha) and β (beta). The Amiga used the Latin-1 character set, which has no Greek letters, so these were the closest things on the keyboard. |
æ (α) |
The angle between the view ray and the water plane. By convention α is the first angle in a problem, often the angle of incidence. Optics textbooks measure incidence from the surface normal; this code measures it from the water surface, the grazing angle, which is the complement of the textbook one. |
ß (β) |
The tilt of the local wave facet, from SlopeMap. β is the conventional name for a second angle. |
A, D |
Altitude and Distance: the camera’s height above this bit of water, and the horizontal distance to it. |
AoD |
“A over D”, the tangent of α, and so the slope of the reflected ray. |
Apr, Dpr, Aprpr |
“A prime”, “D prime” and “A double-prime”, after the prime notation mathematicians use for a related quantity: the height and distance of the candidate pixel being tested, and the camera’s height above that candidate. |
zip, sZip |
The pixel’s index into every per-pixel buffer, and the index of the pixel currently being searched. |
dsX, dOffset |
The march’s current screen X as a double (so it can move by fractions of a pixel), and how far it drifts sideways per row. |
PARC_RNDR_MOTION(n) |
Render motion parameter n for the current frame: 0 is camera altitude, 8 is camera bank, 14 is vertical exaggeration. |
zip deserves a word. It’s a precomputed subscript. Early compilers weren’t reliable about noticing that y * width + x was the same value in five places, so we computed the index once, kept it in a variable, and used it for every buffer. The same habit shows up everywhere in the code. We precomputed anything we could and turned repeated divisions into multiplications, because on a 68000 a divide cost several times what a multiply did, and floating point without an FPU was slower still. Instead of dividing by 180 and multiplying by Ï€ to get radians, Defines.h has PiOver180 as a constant, used 94 times across the code. Instead of dividing elevations by 1,000, Reflection_Render() multiplies by .001. When a divisor was going to be reused, we took its reciprocal once and multiplied by that, as in DEM.c, which normalizes a filter kernel with sum = 1.0 / sum; and then weight[i][j] *= sum; in the loop. A render did these operations millions of times, and every cycle added up on a CISC chip like the m68k. We were voracious about it.
Fade with angle. ReflectionMap[zip] *= (1.0 - AoD); (line 779). Reflections get weaker as the reflected ray gets steeper. It’s a crude Fresnel-like angular falloff. In 1996 I still pronounced it frez-nel. My French teacher would be appalled.
Crude or not, it’s a big part of why our water looked like water. Look down into a lake near your feet and you mostly see into it; look across it toward the far shore and it turns into a mirror. The 1 - AoD falloff does the same thing: steep, downward-looking views get a weak reflection and let the water color through, and distant glancing views reflect strongly. Users and reviewers regularly told us they were amazed that WCS water looked like real water instead of liquid metal, which is what a uniformly mirrored surface in a ray tracer like Bryce tended to look like.
March the ray up the screen. This is the core of it (lines 781 to 805):
for (sY=y-1, dsX=sX-dOffset; sY>=0; sY--, dsX-=dOffset)
{
sX = dsX;
if (sX < 0)
{
dOffset *= (-1.0);
sX = 0;
}
else if (sX > wide)
{
dOffset *= (-1.0);
sX = wide;
}
sZip = sY * settings.scrnwidth + sX;
if (ElevationMap[sZip] < ElevationMap[zip] - MaxWaveAmp)
continue;
Aprpr = PARC_RNDR_MOTION(0) - ElevationMap[sZip] * PARC_RNDR_MOTION(14) * .001;
Dpr = sqrt(zbuf[sZip] * zbuf[sZip] - Aprpr * Aprpr) - D;
Apr = (ElevationMap[sZip] - ElevationMap[zip]) * PARC_RNDR_MOTION(14) * .001;
if (Dpr > 0.0)
{
if (Apr / Dpr > AoD)
break;
} /* if */
} /* for sY=... */The loop walks up the screen one row at a time. It drifts sideways by dOffset so the path converges toward a vanishing point, which is computed earlier and corrected for camera bank. If the path runs off the side of the frame, it bounces back.
For each pixel it visits, it rebuilds that surface’s height above the water (Apr) and its horizontal distance beyond the reflecting point (Dpr), again purely from the stored depth and elevation. If Apr / Dpr is greater than AoD, the stored surface rises above the reflected ray, so the ray has hit it. That’s the same depth-buffer-as-height-field test at the center of modern SSR. It also skips anything lower than the water surface minus the wave height, since that can’t be in the reflection.
Handle a miss. If the ray leaves the top of the frame without hitting anything, the code clamps to the top row (line 809), which in a terrain render is nearly always sky. Modern SSR falls back to an environment map for the same reason. Later versions intersected the reflected ray with the sky dome itself: GetSkyColor() finds where the ray meets a sphere 200 km above the planet and evaluates the sky’s gradients in that direction.
Add ripple. A sinusoidal sideways offset, scaled by distance and animated with the frame number (line 815), makes the reflection shimmer. Our feature sheet at the time said reflections “wiggle and dance just like the real thing.” This is the code that does the wiggling.
Fetch and blend. The code takes a bilinear-weighted sample of the frame around the hit point, from a 3 x 3 neighborhood (lines 846 to 970). It scales that to 80% and blends it into the water pixel by the pixel’s reflectivity.
Smooth. A final 2 x 2 box blur runs over every reflective pixel (line 994 onward).
SlopeMap bends the reflected ray.
Side by side
| Modern SSR | WCS 2, 1995 |
|---|---|
| G-buffer with depth, normal, roughness | zbuf, ElevationMap, SlopeMap, ReflectionMap |
| Rebuild position from depth | Rebuild height and horizontal distance from zbuf and ElevationMap |
| Reflect the view vector about the normal | tan(2ß + æ) about the local wave slope |
| Linear or hierarchical march against depth | Linear march, one screen row per step |
| Hit when the ray passes behind stored depth | Hit when Apr / Dpr > AoD |
| Fetch color from the lit frame | Bilinear fetch from the finished RGB frame |
| Fade by angle and roughness | 1 - AoD falloff times user strength |
| Environment map on a miss | Top row of the frame, usually sky |
What it doesn’t do
The 1996 code only handles horizontal, water-like reflectors. That was all we needed for lakes and oceans, but it means the march only ever goes up the screen, and the hit test works with elevation angles instead of general 3D positions. Crytek’s version handles surfaces at any angle. So did ours, a few years later (see below).
The “normal” is a single slope value per pixel, along the view direction. It tilts the reflected ray up and down. The sideways shimmer is faked with the sine offset, not derived from a real normal.
Like every screen-space technique, it can only reflect what’s in the frame. It got squirrely with a heavily banked camera, where “up the screen” stops meaning “up in the world” and the bank correction can only do so much. It also failed when the thing that should be reflected wasn’t on screen because of how the shot was framed. Views from under a bridge are the classic case: the underside of the bridge should be mirrored in the water, but it’s above the top of the frame, so the ray runs off the screen and gets sky. And if something in the foreground covers the object that should appear in the reflection, that object was never stored, so the reflection can’t show it.
What came later
Gary rewrote the reflection pass from scratch for the Windows-era renderer. The header of RenderPost.cpp says “Built from scratch 11/12/99 by Gary R. Huber,” and it is no longer limited to water. Every reflective pixel stores a reflection percentage and a full surface normal. The pass rebuilds the pixel’s position from its depth, reflects the view vector about the normal (R = 2(N·V)N − V) in real 3D, projects that ray back onto the screen, and walks along it one pixel at a time, interpolating depth as it goes. A hit is a stored surface whose depth falls in a narrow band around the ray’s depth at that pixel. A miss falls back to the sky dome, as above. It is still screen-space: the ray can only find what the renderer kept.
A later optimization in the WCS 5 or 6 era by our colleague, the late Frank Weed II (he passed away in 2010), scanned the image by COLUMNS, starting at the top left and going to the bottom left, then moving one column right and repeating. Frank discovered this was MUCH more cache coherent because the reflective pixel below the previous one would re-iterate almost the same input G-buffer pixels, which were likely hot in the cache. You can see the column loop in the high-quality reflection pass, RenderReflections_BTFRAGS().
That kept geometry got deeper in WCS 6. Gary’s PixelManager (“Written from scratch on 7/5/99”) replaced the single z-buffer with an A-buffer, the technique Loren Carpenter described in 1984. Every pixel keeps a depth-sorted list of fragments, each with its own color, depth, alpha, a 64-bit coverage mask (an 8 x 8 grid of subpixels, all sharing that fragment’s single color and depth), and, for reflective surfaces, its reflection percentage and normal. The WCS 6 manual called it “Pixel Fragment Rendering” and said the result “should be as good as doing ray tracing with 64 rays per pixel – at least as far as antialiasing goes.” The high-quality reflection pass walks each pixel’s whole fragment list, not just the front one, so it can reflect things that are hidden from the camera. A “Render Occluded” option told the renderer to keep those hidden fragments instead of culling them to save memory. In the manual’s words, it “will force WCS to render those faces of the object not visible to the camera. This will allow those faces to reflect in water and other reflective surfaces, in circumstances where the reflection of the surface would be visible, even when the surface itself is not.” That fixed the foreground-occlusion problem, at the cost of a deeper fragment stack.
WCS 6 and Visual Nature Studio 2 shipped this in 2002 and 2003. The manual’s HTML files are dated August 2002, and they say 3D object materials “can now Reflect their surroundings and other 3D Objects, as well as parts of their own geometry,” and that terrain materials “can now reflect their surroundings, as well as parts of their own geometry.” That’s arbitrary-surface screen-space reflection, in a shipping product, about nine years before Crytek’s talk.
The fade got a proper treatment too. The later code replaced the 1 - AoD approximation with the actual Fresnel equation, in UsefulMath.cpp: a version specialized for water (index of refraction 1.33, which is where the odd constant 0.7689, 1.33² − 1, comes from) and a general version that takes any refractive index. The WCS 6 manual describes water and foam as using “theoretically correct Fresnel reflectivity.”
Could the graphics cards of the day have done any of this? Not really. WCS 6 shipped in 2002, the year the first DirectX 9 card, ATI’s Radeon 9700, arrived. Its Shader Model 2.0 pixel shaders were limited to 64 arithmetic and 32 texture instructions with no loops or branches, which rules out marching a ray across the screen. Shader Model 3.0, with real loops and dynamic branching, arrived with NVIDIA’s GeForce 6800 in 2004, and shader-based ray marching through a height field showed up soon after as parallax occlusion mapping. Even then, building a full G-buffer and marching every reflective pixel through it at frame rate had to wait for the DirectX 10 and 11 generation. Crytek’s SSR shipped in a DirectX 11 update. We had the luxury of minutes per frame on the CPU.
And it wasn’t real time. WCS was an offline renderer. The readme for the modern Emerald-Caesar build (below) times its canyon example scene (CanyonSet000) at 4 hours 59 minutes on a 25 MHz 68040 A4000T with the original 2.04 build.
What it does have is the structure of the technique: G-buffer-like per-pixel data written during rasterization, a reflected ray for every reflective pixel, and a march through the frame’s own depth to find what that ray hits and fetch its color.
See it run
You don’t have to take my word for any of this. WCS is open source now, and it still runs. worldconstructionset.com has the current Amiga release, links to the source on GitHub, and the original software from the 3D Nature days.
The current Amiga build is WCS 2.031 “Emerald-Caesar,” released in October 2025 in time for the Amiga’s 40th anniversary. Alexander Fritsch recompiled the released source with a modern GCC, fixed bugs, added catalogs for German, Dutch, Italian and Spanish, and added a full-color render display for P96 and CyberGraphX graphics cards. It’s also on Aminet as gfx/3d/wcs. It runs on real 68020-or-better Amigas with an FPU, on AROS, and under WinUAE emulation. You’ll need the original WCS 2 install disks, which 3D Nature gives away as ADF images in the Legacy Software downloads. The new build is faster, too: the canyon example that took 4:59:44 with the original 2.04 build on a 25 MHz 68040 renders in 3:57:40 with his 68040 version. Turn on reflections in the Render Settings Editor, give the water some Reflection Strength, and you can watch Reflection_Render() do its thing on the same code I walked through above.
If you’d rather not set up an Amiga, the AMIGA alive channel posted QC: Terraforming with World Construction Set 2.0 (February 2025), which shows WCS in use and renders this same canyon scene, reflections and all. The reflections start at 5:35.
When did this exist?
Having the code is one thing. Showing when it was written is harder.
Git doesn’t help
The repository was uploaded to GitHub in 2015. Before that we used SVN and even CVS prior to that internally. GlobeMapSupport.c has exactly one commit, d15d569, “Initial source upload,” dated October 5, 2015. Git doesn’t record file modification times, so the original dates were lost on the way in. They might be in some ancient backup.
The binary does
The repo also contains a compiled Amiga executable, Amiga/WCS. Version.h builds the version string from the compiler’s __DATE__ macro:
#define APP_VERSION APP_TLA " " APP_VERS "." APP_REV " " __DATE__So the build date is baked into the binary:
$ strings WCS | grep "WCS 2"
WCS 2.031 Jul 19 1996
$VER: WCS 2.031 Jul 19 1996The second line is the standard Amiga $VER: version cookie, which the Version command reads.
The same binary contains the strings “No Reflections”, “Reflection %” and “Out of memory allocating Reflection buffer!” A build with the reflection code existed by July 19, 1996, going by the build machine’s clock.
strings on the Amiga executable in the repository.Our archives
Our archives from our Questar Productions days fill in the earlier dates:
| Date | Document | What it shows |
|---|---|---|
| Nov 20, 1994 | WCS 1 feature list | No reflections. WCS 1 never had them. |
| Dec 22, 1994 | V2 planning note | Reflections on the V2 list. |
| Sep 12, 1995 | V2 feature sheet | Reflections described in full: “Any water body can reflect the terrain and sky… They wiggle and dance just like the real thing.” |
| Dec 6, 1995 | “Version 2.x Prerelease” manual addendum | A Reflections section: “Reflections are turned on in the Render Settings Editor.” |
| Mar 3, 1996 | GlobeMapSupport.c |
Contains Reflection_Render(). |
| Apr 10, 1996 | Version.h |
Version 2.04, code-named “Ruby”. |
| Apr 11, 1996 | Box label for the V2 prerelease | “This revision includes waves, clouds, reflections…” |
The March 1996 copy of Reflection_Render() is slightly different from the one on GitHub. It guards against a negative square root in the march, and the ripple offset is commented out. It’s a different snapshot of the same function, which is what you’d expect from code that was still being tuned.
Reflection_Render().Those are our own records, though, and anyone can make a file with an old date on it. What I really wanted was something dated that someone else kept.
Light-ROM 3
In September 1995 we sent a Grand Canyon fly-through, rendered in WCS 2, to Michael Meshew for Light-ROM 3. That was a three-CD collection of LightWave and Amiga graphics resources that Graphic Detail Inc. published in late 1995. Somebody archived the discs on the Internet Archive in 2020: archive.org/details/light-rom-3.
The disc images still carry their original file dates:
| Disc | File | Date on disc |
|---|---|---|
| Disc 2 | PROGRAMS/AMIGA/DEMOS/WCS_V2/WCS_V2.IFF |
1995-09-14 13:29:48 |
| Disc 1 | WCS_JPG/CANYON/CANYN000 to CANYN359 |
1995-09-26, 12:41 to 14:50 |
| Disc 1 | WCS.TXT (our V2 feature text) |
1995-10-06 |
WCS_V2.IFF is a 752 x 480, 24-bit Amiga ILBM image. Its header records that it was saved by ASDG’s Art Department Professional on “14.09.95.” It’s the image at the top of this post, brightened for today’s displays: the canyon at night, rippled water, a large moon, the moon’s glowing path across the water, and the dark mesas reflected in it.
The 360 JPEG frames on disc 1 are the full fly-through. It starts at sunset, with the mesas mirrored in the water and broken up by ripples, and ends after moonrise with the moon path on the water.
Here’s the whole thing. All 360 frames play at 15 frames per second, the rate we recommended at the time; watch the water once the moon comes up.
A note on the shape of these images: if they look a little squat, it’s because they are. These renders were made for NTSC video, where each pixel of a 752 x 480 overscan frame is taller than it is wide, about 0.85:1 on a 4:3 screen. On a monitor of the day they looked right. On today’s square pixels they come out roughly 15% too short. The magazine scans further down show them the way they were meant to be seen: printed at the correct aspect, with a round moon.
For contrast, the same disc carries our WCS 1 demo images, dated March 24, 1995. The WCS 1 Grand Canyon image is a sunset over dry canyon walls. No water, no reflections.
The dates on a CD come from whoever mastered it, not from an independent authority. But they’re consistent with everything else: the case says “Copyright © 1995 Graphic Detail Inc.”, and Amiga Report ran a news item on Light-ROM 3 on November 16, 1995, describing “a stunning 360 frame Jpeg image sequence created with World Construction Set Version 2.”
The press
We showed V2 publicly at the Video Toaster User Expo at the Universal City Hilton, November 1 to 4, 1995. Amiga Report’s show coverage (issue 3.20) described our tape of “beaches with waves rolling onto them” with whitecaps and cloud shadows, and said the output “looked LOTS better than Kai’s Power Tools: BRYCE.” MicroTimes also covered our booth.
Once V2 was out, the reviewers noticed the water:
| Date | Publication | What it said |
|---|---|---|
| Oct 1995 | Amiga Shopper #54 | Our “‘To Do’ list already includes adding clouds, reflections, waves…” |
| Dec 1995 | Amiga Computing #94 | WCS 1 “doesn’t yet offer reflective water or clouds.” |
| Feb 1996 | Amiga Magazine (Italy) #75 | “…meravigliose riflessioni, senza per altro le lunghe attese dovute alla resa di tipo ray tracing.” “…marvellous reflections, without the long waits of ray-traced rendering.” |
| Mar 1996 | Amiga Format #82 | “water (complete with realistic reflections and waves!)” |
| Mar 31, 1996 | Amiga Report 4.05 | Our release: “Waves, breakers and reflections on the water.” |
The Italian one is my favorite. Antonio De Lorenzo’s column says the beaches come “complete with marvellous reflections, without the long waits of ray-traced rendering.” That’s the entire point of the technique, noticed by a reviewer in February 1996. Page 12 of that issue prints the same canyon at sunset and at night, captioned “the same scene rendered at two different times of day with the new release of the program.” The scans are on the Internet Archive: page 11, page 12.
“Without ray tracing” was a hard sell with our customers, precisely because they knew graphics. Our salespeople spent long sessions on the phone with prospective and current users who insisted the reflections had to be ray traced, because to their eyes they sure looked it. Sales would eventually transfer the call to one of us on the technical side, and we’d say it again: no, it really isn’t ray traced, and it couldn’t be, because the scene doesn’t exist in memory to trace against. Then we’d send them back to sales, and they usually bought it.
The V2 prerelease shipped to Amiga customers in spring 1996, with a coupon for a free upgrade to the final release.
NAB 1995, as I remember it
There’s one more piece, and I have to be upfront that it’s memory. In April 1995 Gary, Mindy and I drove to Las Vegas for NAB, which ran April 10 to 13, all three of us squeezed into Gary’s little Toyota coupe. We didn’t have a booth. We camped in the desert northeast of town to save money and cleaned up at a truck stop before heading in. I remember sitting outside in the sun showing printouts of WCS 2 canyon renders with reflective water to some of the Amiga press, as a sneak preview before V2 had been announced. One of them, as I recall, was the late Harv Laser (1949–2022), who had reviewed WCS 1 for the January 1995 AmigaWorld.
Our records from the time confirm the trip and that we went without a booth. I haven’t found anything in print about the printouts. If that memory is right, the reflection code was working by April 1995. If not, the dated evidence still puts rendered reflections on a published CD by September 1995.
Prior art
Plenty of the ingredients existed before 1995. Saito and Takahashi explicitly called their per-pixel storage of geometric properties, including depth and normals, “G-buffers” at SIGGRAPH 1990. Musgrave published fast ray tracing of height fields in 1988. Image-based rendering work in the late 1990s (McMillan’s 1997 thesis, Lischinski and Rappoport in 1998) marched rays through depth images.
The closest academic work I found is from 2005 and 2006. Szirmay-Kalos and colleagues searched along reflected rays against distance values stored in the texels of an environment map around an object (Computer Graphics Forum, 2005). Popescu and colleagues traced reflected rays against reflected-scene impostors built from depth maps (Computer Graphics Forum, 2006). Same search, different buffer: in both cases the depth came from separate images, not the main frame.
After that came screen-space ambient occlusion (Crytek, 2007), then Just Cause 2 (2010) and Crytek’s reflections (2011).
What about the demoscene? Amiga and PC demos of the era did plenty of water tricks, but the ones I know of copied rows of the screen upside down as a flat 2D mirror. No depth buffer, no per-pixel ray, no search for a hit. That produces a reflection-like image. It isn’t the same algorithm.
I’ll add the obvious caveat. 1990s Usenet, demoscene productions and trade magazines are poorly indexed. If someone did this before 1995, I’d like to hear about it.
So did we invent it?
Here’s what I can back up. In 1995, World Construction Set 2 rendered water reflections by storing G-buffer-like per-pixel data during rasterization, computing a reflected ray for every reflective pixel, and marching that ray through the frame’s own depth and elevation data to find and fetch what it hit. Rendered frames with those reflections were on a published CD by September 1995 and in print by February 1996, and the renderer shipped to customers in the spring of 1996. I haven’t found an earlier published description of this specific combination: tracing a reflected ray through depth belonging to the rendered view and fetching the reflection from that same rendered frame. The earliest shipped use elsewhere that I know of is Just Cause 2, in 2010.
We never named it, wrote a paper or filed a patent, so it had no influence on anyone. Crytek arrived at SSR independently, and they made it run in real time on a GPU, which is the hard part. We had generalized our own version to arbitrary surfaces by 2002, but it stayed inside an offline renderer, and nobody outside our customers ever heard about it. We were solving a narrower problem: making lakes look right on an Amiga that couldn’t hold the scene in memory.
Still, the code is on GitHub and the renders are on the Internet Archive. I’m comfortable saying Gary and I built screen-space reflections, in a shipping product, about sixteen years before the technique got its name.
Addendum: Gary weighs in
I sent Gary a draft of this post. Here’s our exchange, lightly trimmed.
Gary wrote:
Might have been inspired by some Siggraph article, can’t remember. I do remember writing the code and all the optimizations and enhancements you mention. There might be notes around here. I did a lot of my daydreaming at the Dushanbe Teahouse on yellow legal tablets.
I wrote back:
I am gonna look for the notebooks you gave me. They’re in a box somewhere, but I saved them. Somewhere.
I researched, and there isn’t any academic material on the technique prior to us that I can find. No SIGGRAPH papers on the technique (though G-buffer dates from the 90s, nobody used it for screen-space reflections).
I actually remember coming up with the idea and discussing what extra data it would need with you and worrying about the extra attribute memory and how to make it work with segments. We even talked about whether vertical segments with overscan would fix it but vertical segmenting didn’t happen until much later for some reason.
I also remember being super excited when you wrote it and it actually worked. The first renders were so-so because we didn’t have the blur or the wiggle or faux Fresnel. I remember the first reflection renders looking like liquid metal and it looked cool but you didn’t like it because it didn’t look like water and I didn’t know why but you knew (or you went back and figured it out). I think you added wiggle and Fresnel, because I didn’t even know what the Fresnel coefficient of water was, but I feel like the blur suggestion was me since it feels like my kind of lazy “well, what if we just mask-blurred it to get rid of the aliasing?” technique. I almost remember making that suggestion, but I don’t think I can trust that memory.
My main recollection surrounding the technique is that I was familiar with the Eurodemo technique of drawing a squashed inverted copy of the sky on the ground plane (sometimes just using the copper to re-emit the selected scanlines from the sky rows). I’d made a 16-gray-scale Amiga demo in 640×480 with a sky of stars reflected in the ground plane, and then one at a time, the stars fell straight down to the horizon. It looked terrible in motion because the interlace timing made the stars vertically jitter weirdly. I might even still have that source on a floppy. I had also been trying to do weird things with grayscale drawing like antialiased lines while I was working at PP&S but they weren’t impressed and didn’t hire me for programming. Last year I re-created that star-fall effect (inefficiently) in a ShaderToy shader to explain something to someone about shaders.
At the time we’d talked about how that copy-mirror technique wouldn’t work because of the unconstrained camera pose, and then I had the idea of rendering the whole scene a second time from a camera pose that was flipped over to be under the water surface, looking up. Which would have taken 2x as long, and only worked with a mirror flat surface, which wasn’t great. Then I chewed on the idea of rendering the looking-up reflected scene into an image and using some kind of distortion map to perturb the pixels according to the reflection angle of the waves on the surface. And then came the weird-ideas realization that we could actually do that during the MAIN render, and not need a whole second reflections-scene render pass (though the reflections scene would have captured things the main render wouldn’t, like the undersides of things).
Appendix A: GlobeMapSupport.c
The complete file, converted from Latin-1 to UTF-8. Reflection_Render() is lines 719 to 1023. On GitHub: current version, pinned version. Licensed under the OSGPL.
Show the full source (1,306 lines)
/* GlobeMapSupport.c (ne gisglobemapsupport.c 14 Jan 1994 CXH)
** Some auxiliary calculation and IO routines for mapping in GIS renderer.
** Built (ripped) from gis.c on 24 Jul 1993 by Chris "Xenon" Hanson.
** Original code by Gary R. Huber.
*/
#include "WCS.h"
/*********************************************************************/
double ComputeBanking(short frame)
{
short i, j, ffr, lfr, compfr;
double angledif, nextx = 0.0, nexty = 0.0, lastx = 0.0, lasty = 0.0,
lonscale, Banking;
if (! KT || KT_MaxFrames < 3) /* no auto banking allowed */
return (0.0);
if (frame <= 1)
compfr = 2;
else if (frame >= KT_MaxFrames)
compfr = KT_MaxFrames - 1;
else
compfr = frame;
ffr = compfr - 6;
lfr = compfr + 6;
if (ffr < 2)
ffr = 2;
if (lfr > KT_MaxFrames - 1)
lfr = KT_MaxFrames - 1;
Banking = 0.0;
lonscale = LATSCALE * cos(KT[1].Val[0][frame] * PiOver180);
for (i=ffr; i<=lfr; i++)
{
for (j=1; j<9; j++)
{
if (i > j)
{
if (KT[2].Key)
lastx = (KT[2].Val[0][i] - KT[2].Val[0][i - j]) * lonscale;
if (KT[1].Key)
lasty = (KT[1].Val[0][i] - KT[1].Val[0][i - j]) * LATSCALE;
} /* if */
else
{
if (KT[2].Key)
lastx = (KT[2].Val[0][i] - KT[2].Val[0][1]) * lonscale;
if (KT[1].Key)
lasty = (KT[1].Val[0][i] - KT[1].Val[0][1]) * LATSCALE;
} /* else */
if (i <= KT_MaxFrames - j)
{
if (KT[2].Key)
nextx = (KT[2].Val[0][i + j] - KT[2].Val[0][i]) * lonscale;
if (KT[1].Key)
nexty = (KT[1].Val[0][i + j] - KT[1].Val[0][i]) * LATSCALE;
} /* if */
else
{
if (KT[2].Key)
nextx = (KT[2].Val[0][KT_MaxFrames] - KT[2].Val[0][i]) * lonscale;
if (KT[1].Key)
nexty = (KT[1].Val[0][KT_MaxFrames] - KT[1].Val[0][i]) * LATSCALE;
} /* else */
angledif = findangle2(nextx, nexty) - findangle2(lastx, lasty);
if (angledif > Pi) angledif -= TwoPi;
else if (angledif < -Pi) angledif += TwoPi;
Banking += settings.bankfactor * angledif * PiUnder180;
} /* for j=1... */
Banking /= 8.0;
} /* for i=ParHdr.FirstFrame+1... */
return (Banking);
} /* ComputeBanking() */
/*********************************************************************/
void setview(void)
{
VP.lat = PARC_RNDR_MOTION(1);
VP.lon = PARC_RNDR_MOTION(2);
VP.alt = PARC_RNDR_MOTION(0);
if (settings.bankturn && KT)
{
short i;
double BankFactor = 0.0;
for (i=-5; i<=5; i++)
{
BankFactor += ComputeBanking(frame + i);
}
PARC_RNDR_MOTION(8) += (BankFactor / 11.0);
} /* if automatic turn banking - multiple frames to smooth out the wobblies */
/* set these values for fixfocus = 0 before converting view point */
FP.lat = VP.lat;
FP.lon = VP.lon;
FP.alt = -EARTHRAD;
if ((settings.lookahead && KT) || fixfocus)
{
if (settings.lookahead && KT)
{
short AheadFrame;
if (frame + settings.lookaheadframes > KT_MaxFrames)
AheadFrame = KT_MaxFrames;
else
AheadFrame = frame + settings.lookaheadframes;
if (AheadFrame == frame && frame > 1)
{
if (KT[1].Key)
FP.lat = KT[1].Val[0][frame] + (KT[1].Val[0][frame] - KT[1].Val[0][frame - 1]);
else
FP.lat = PARC_RNDR_MOTION(1);
if (KT[2].Key)
FP.lon = KT[2].Val[0][frame] + (KT[2].Val[0][frame] - KT[2].Val[0][frame - 1]);
else
FP.lon = PARC_RNDR_MOTION(2);
if (KT[0].Key)
FP.alt = KT[0].Val[0][frame] + (KT[0].Val[0][frame] - KT[0].Val[0][frame - 1]);
else
FP.alt = PARC_RNDR_MOTION(0);
} /* if focus and camera same pt */
else
{
if (KT[1].Key)
FP.lat = KT[1].Val[0][AheadFrame];
else
FP.lat = PARC_RNDR_MOTION(1);
if (KT[2].Key)
FP.lon = KT[2].Val[0][AheadFrame];
else
FP.lon = PARC_RNDR_MOTION(2);
if (KT[0].Key)
FP.alt = KT[0].Val[0][AheadFrame];
else
FP.alt = PARC_RNDR_MOTION(0);
} /* else normal look ahead */
} /* if lookahead */
else
{
FP.lat = PARC_RNDR_MOTION(4);
FP.lon = PARC_RNDR_MOTION(5);
FP.alt = PARC_RNDR_MOTION(3);
} /* if fixfocus && !settings.focpath && !lookahead */
} /* if settings.lookahead || settings.focpath */
if ((FP.lat == VP.lat && FP.lon == VP.lon) || ! fixfocus)
{
convertpt(&VP);
convertpt(&FP);
findposvector(&FP, &VP);
yrot = findangle2(FP.x, FP.z);
rotate(&FP.x, &FP.z, yrot, yrot);
xrot = findangle2(FP.y, FP.z);
rotate(&FP.y, &FP.z, xrot, xrot);
zrot = 0.0;
} /* if same camera and focus lat/lon */
else
{
DP.lat = FP.lat;
DP.lon = FP.lon;
DP.alt = FP.alt + 1000;
convertpt(&VP);
convertpt(&FP);
findposvector(&FP, &VP);
yrot = findangle2(FP.x, FP.z);
rotate(&FP.x, &FP.z, yrot, yrot);
xrot = findangle2(FP.y, FP.z);
rotate(&FP.y, &FP.z, xrot, xrot);
convertpt(&DP);
findposvector(&DP, &VP);
rotate(&DP.x, &DP.z, yrot, findangle2(DP.x, DP.z));
rotate(&DP.y, &DP.z, xrot, findangle2(DP.y, DP.z));
zrot=findangle2(DP.x, DP.y);
rotate(&DP.x, &DP.y, zrot, zrot);
} /* else fix focus */
BuildRotationMatrix(-xrot, -yrot, zrot, NoBankMatx);
/* FP.q = sqrt(FP.x * FP.x + FP.y * FP.y + FP.z * FP.z);*/
zrot -= (PARC_RNDR_MOTION(8) * PiOver180);
cosxrot = cos(-xrot);
cosyrot = cos(-yrot);
coszrot = cos(-zrot);
sinxrot = sin(-xrot);
sinyrot = sin(-yrot);
sinzrot = sin(-zrot);
BuildRotationMatrix(-xrot, -yrot, zrot, ScrRotMatx);
/* Sun position */
DP.alt = 149.0E+6;
DP.lat = PARC_RNDR_MOTION(15);
DP.lon = PARC_RNDR_MOTION(16);
/* This call sets the cartesian coordinates of the sun in SP */
getscrnpt(&SP, NULL);
VectorMagnitude(&SP);
UnitVector(&SP);
memcpy(&SNB, &SP, sizeof (struct coords));
RotatePoint(&SNB, NoBankMatx);
NegateVector(&SNB);
sprintf(str, "Frame %hd VP.lat=%f, VP.lon=%f, VP.alt=%f\n",
frame, VP.lat * PiUnder180, VP.lon * PiUnder180, VP.alt);
Log(DTA_NULL, str);
sprintf(str, "FP.lat=%f, FP.lon=%f, FP.alt=%f\n",
FP.lat * PiUnder180, FP.lon * PiUnder180, FP.alt);
Log(MSG_UTIL_TAB, str);
sprintf(str, "SP.x=%f, SP.y=%f, SP.z=%f\n",
SP.x, SP.y, SP.z);
Log(MSG_UTIL_TAB, str);
sprintf(str, "Y rot=%f, X rot=%f, Z rot=%f\n",
yrot * PiUnder180, xrot * PiUnder180, zrot * PiUnder180);
Log(MSG_UTIL_TAB, str);
sprintf(str, "Q max=%f, Q focus=%f, Banking=%f\n",
qmax, FP.q, PARC_RNDR_MOTION(8));
Log(MSG_UTIL_TAB, str);
} /* setview() */
/*********************************************************************/
short setquickview(void)
{
VP.lat = PARC_RNDR_MOTION(1);
VP.lon = PARC_RNDR_MOTION(2);
VP.alt = PARC_RNDR_MOTION(0);
if (settings.bankturn && KT)
{
short i;
double BankFactor = 0.0;
for (i=-5; i<=5; i++)
{
BankFactor += ComputeBanking(frame + i);
}
PARC_RNDR_MOTION(8) += (BankFactor / 11.0);
} /* if automatic turn banking - multiple frames to smooth out the wobblies */
FP.lat = VP.lat;
FP.lon = VP.lon;
FP.alt = -EARTHRAD;
if ((settings.lookahead && KT) || fixfocus)
{
if (settings.lookahead && KT)
{
short AheadFrame;
if (frame + settings.lookaheadframes > KT_MaxFrames)
AheadFrame = KT_MaxFrames;
else
AheadFrame = frame + settings.lookaheadframes;
if (AheadFrame == frame && frame > 1)
{
if (KT[1].Key)
FP.lat = KT[1].Val[0][frame] + (KT[1].Val[0][frame] - KT[1].Val[0][frame - 1]);
else
FP.lat = PARC_RNDR_MOTION(1);
if (KT[2].Key)
FP.lon = KT[2].Val[0][frame] + (KT[2].Val[0][frame] - KT[2].Val[0][frame - 1]);
else
FP.lon = PARC_RNDR_MOTION(2);
if (KT[0].Key)
FP.alt = KT[0].Val[0][frame] + (KT[0].Val[0][frame] - KT[0].Val[0][frame - 1]);
else
FP.alt = PARC_RNDR_MOTION(0);
} /* if focus and camera same pt */
else
{
if (KT[1].Key)
FP.lat = KT[1].Val[0][AheadFrame];
else
FP.lat = PARC_RNDR_MOTION(1);
if (KT[2].Key)
FP.lon = KT[2].Val[0][AheadFrame];
else
FP.lon = PARC_RNDR_MOTION(2);
if (KT[0].Key)
FP.alt = KT[0].Val[0][AheadFrame];
else
FP.alt = PARC_RNDR_MOTION(0);
} /* else normal look ahead */
} /* if lookahead */
else
{
FP.lat = PARC_RNDR_MOTION(4);
FP.lon = PARC_RNDR_MOTION(5);
FP.alt = PARC_RNDR_MOTION(3);
} /* if fixfocus && !lookahead */
} /* if lookahead or fixfocus */
if ((FP.lat == VP.lat && FP.lon == VP.lon) || ! fixfocus)
{
/* Jamie didn't like this
if (fixfocus)
{
if (! User_Message("Parameters Module",
"Warning!\nCamera and focus at same latitude\
and longitude coordinates (possibly a result of \"Look Ahead\" Render Setting\
enabled).\n\
View may be different than expected.", "OK|Cancel", "oc"))
return (0);
}
*/
convertpt(&VP);
convertpt(&FP);
findposvector(&FP, &VP);
yrot = findangle2(FP.x, FP.z);
rotate(&FP.x, &FP.z, yrot, yrot);
xrot = findangle2(FP.y, FP.z);
rotate(&FP.y, &FP.z, xrot, xrot);
zrot = 0.0;
} /* if !fixfocus */
else
{
DP.lat = FP.lat;
DP.lon = FP.lon;
DP.alt = FP.alt + 1000;
convertpt(&VP);
convertpt(&FP);
findposvector(&FP, &VP);
yrot = findangle2(FP.x, FP.z);
rotate(&FP.x, &FP.z, yrot, yrot);
xrot = findangle2(FP.y, FP.z);
rotate(&FP.y, &FP.z, xrot, xrot);
convertpt(&DP);
findposvector(&DP, &VP);
rotate(&DP.x, &DP.z, yrot, findangle2(DP.x, DP.z));
rotate(&DP.y, &DP.z, xrot, findangle2(DP.y, DP.z));
zrot = findangle2(DP.x, DP.y);
rotate(&DP.x, &DP.y, zrot, zrot);
} /* if fixfocus */
BuildRotationMatrix(-xrot, -yrot, zrot, NoBankMatx);
/* FP.q = sqrt(FP.x * FP.x + FP.y * FP.y + FP.z * FP.z);*/
zrot -= PARC_RNDR_MOTION(8) * PiOver180;
cosxrot = cos(-xrot);
cosyrot = cos(-yrot);
coszrot = cos(-zrot);
sinxrot = sin(-xrot);
sinyrot = sin(-yrot);
sinzrot = sin(-zrot);
BuildRotationMatrix(-xrot, -yrot, zrot, ScrRotMatx);
DP.alt = 149.0E+6;
DP.lat = PARC_RNDR_MOTION(15);
DP.lon = PARC_RNDR_MOTION(16);
getscrnpt(&SP, NULL);
VectorMagnitude(&SP);
UnitVector(&SP);
memcpy(&SNB, &SP, sizeof (struct coords));
RotatePoint(&SNB, NoBankMatx);
NegateVector(&SNB);
return (1);
} /* setquickview() */
/*********************************************************************/
/* Obsolete - see Tree.c
void loadtreemodel(void)
{
/* reads standard GENGYM ascii output */
char filename[255];
short error, spp, dbh, i, fractlevel, dataread = 0;
long mapyear, year, species;
double overlap,latitude,polyarea,polyperacre, sum, runsum;
FILE *fdat;
for (spp=0; spp<MAXSPECIES; spp++) {
for (dbh=0; dbh<MAXDBH; dbh++) standdata[spp][dbh].tpa=0.0;
}
gets(str);
printf("Name of file to be read: ");
gets(filename);
printf("Year to read: ");
gets(str);
mapyear=atoi(str);
if ((fdat=fopen(filename,"r"))==NULL) {
printf("Can't open file %s to read\n",filename);
settings.treemodel=0;
return;
}
for (i=1; i<10; i++) fgets(str,119,fdat);
while (1) {
error=fscanf(fdat,"%hd",&year);
if (error==EOF || error==0) break;
if (year>mapyear) break;
if (year==mapyear) {
fscanf(fdat,"%d",&species);
switch (species) {
case 93: /* englemann spruce */
spp=0;
break;
case 108: /* lodgepole pine */
spp=1;
break;
case 19: /* subalpine fir */
spp=2;
break;
default:
spp=0;
break;
}
fscanf(fdat,"%hd",&dbh);
if (dbh>=MAXDBH) {
fgets(str,119,fdat);
continue;
}
fscanf(fdat,"%le %hd %hd",&standdata[spp][dbh].tpa,
&standdata[spp][dbh].ht,&standdata[spp][dbh].pctcrn);
dataread=1;
}
else fgets(str,119,fdat);
}
fclose(fdat);
if (!dataread) {
printf("No tree model data read!\n");
settings.treemodel=0;
return;
}
printf("Tree model read\n");
sum=0.0;
for (dbh=0; dbh<MAXDBH; dbh++) treearea[dbh]=
Pi*((4.344+1.029*(dbh+.5))/2.0)*((4.344+1.029*(dbh+.5))/2.0);
for (spp=0; spp<MAXSPECIES; spp++) {
for (dbh=0; dbh<MAXDBH; dbh++) {
sum +=(standdata[spp][dbh].tpa*treearea[dbh]);
}
}
if (sum<=0.0) {
printf("Error reading model data\n");
settings.treemodel=0;
return;
}
printf("Canopy closure= %f sq ft/acre = %f%%\n",sum,sum/435.60);
EnterStock:
printf("Clustering level in percent (overlap): ");
scanf("%le",&overlap);
printf("Fractal level (level 0 = 1201 pts per degree): ");
scanf("%d",&fractlevel);
printf("Latitude of mapped area (degrees): ");
scanf("%le",&latitude);
polyarea=.5*(LATSCALE*5280.0/1200.0)*(LATSCALE*5280.0/1200.0)*
cos(latitude*PiOver180);
for (i=1; i<=fractlevel; i++) polyarea /=4.0;
polyperacre=43560.0/polyarea;
if (overlap<0.0 || overlap>=100.0) goto EnterStock;
overlap *=.01;
overlap=1.0-overlap;
runsum=0.0;
for (spp=0; spp<MAXSPECIES; spp++) {
for (dbh=0; dbh<MAXDBH; dbh++) {
standdata[spp][dbh].tpa *=overlap;
standdata[spp][dbh].tpa /=(polyarea/treearea[dbh]);
runsum +=standdata[spp][dbh].tpa/polyperacre;
standdata[spp][dbh].percent=runsum;
}
}
printf("Tree model percents computed. Occupied sites = %f\n",runsum);
if (runsum>1.0)
printf("Warning occupatipon > 100% will result in loss of high dbh items\n");
printf("Species 93 (Englemann Spruce) ecosystem: ");
scanf("%hd", &treeeco[0]);
printf("Species 108 (Lodgepole Pine) ecosystem: ");
scanf("%hd", &treeeco[1]);
printf("Species 19 (Subalpine Fir) ecosystem: ");
scanf("%hd", &treeeco[2]);
printf("Default ecosystem (no trees): ");
scanf("%hd", &settings.defaulteco);
} /* loadtreemodel() */
*/
/*********************************************************************/
short makesky(short renderseg, struct Window *win)
{
short rotatesky = 0, zenithline = 0, MergePts, red, green, blue, error = 0;
long x, y, zip, colval;
double angle, xxx, yyy, halfsky, skyfact, maxskyfact, lonscale;
struct BusyWindow *BWDE;
BWDE = BusyWin_New("Sky", settings.scrnheight, 0, 'BWDE');
if (! settings.horfix)
{
lonscale = cos(FP.lat) * LATSCALE;
xxx = (FP.lon - VP.lon) * PiUnder180 * lonscale;
yyy = (FP.lat - VP.lat) * PiUnder180 * LATSCALE;
angle = findangle2(xxx, yyy);
xxx = qmax * sin(angle);
yyy = qmax * cos(angle);
DP.lat = VP.lat * PiUnder180 + yyy / LATSCALE;
DP.lon = VP.lon * PiUnder180 + xxx / lonscale;
DP.alt = 0.0;
getscrnpt(NULL, NULL);
horline = PARC_RNDR_MOTION(7) - (DP.y * PROJPLANE / DP.z) / vertscale;
DP.lat = VP.lat * PiUnder180 + yyy / LATSCALE;
DP.lon = VP.lon * PiUnder180 + xxx / lonscale;
DP.alt = settings.zenith; /* qmax; */
convertpt(&DP);
findposvector(&DP, &VP);
rotate(&DP.x, &DP.z, yrot, findangle2(DP.x, DP.z));
rotate(&DP.y, &DP.z, xrot, findangle2(DP.y, DP.z));
rotate(&DP.x, &DP.y, zrot - PARC_RNDR_MOTION(8) * PiOver180,
findangle2(DP.x, DP.y));
zenithline = PARC_RNDR_MOTION(7) - (DP.y * PROJPLANE / DP.z) / vertscale;
} /* if ! horfix, horizon and zenith computed at qmax from viewer */
sprintf(str, "horline=%d, zenithline=%d\n",horline,zenithline);
Log(MSG_UTIL_TAB, str);
zenithline = horline - zenithline;
if (horpt >= settings.scrnwidth / 2)
{
maxskyfact = sqrt(zenithline * zenithline +
(horpt * horstretch) * (horpt * horstretch));
} /* if horpt is right of center */
else
{
maxskyfact = sqrt(zenithline * zenithline +
((wide - horpt) * horstretch) * ((wide - horpt) * horstretch));
} /* else */
redsky = (PARC_RNDR_COLOR(3, 0) - (float)PARC_RNDR_COLOR(4, 0)) / maxskyfact;
greensky = (PARC_RNDR_COLOR(3, 1) - (float)PARC_RNDR_COLOR(4, 1)) / maxskyfact;
bluesky = (PARC_RNDR_COLOR(3, 2) - (float)PARC_RNDR_COLOR(4, 2)) / maxskyfact;
srand48(1010);
halfsky = settings.skyalias * .5;
if (fabs(PARC_RNDR_MOTION(8)) > .0001)
rotatesky = 1;
renderseg *= settings.scrnheight;
zip = 0;
for (y=0; y<settings.scrnheight; y++)
{
for (x=0; x<settings.scrnwidth; x++, zip++)
{
if (((unsigned int)bytemap[zip]) < 100)
{
DP.y = (double)y + renderseg - horline;
DP.x = (double)x - horpt;
if (rotatesky) getskypt();
skyfact = sqrt(DP.y * DP.y + (DP.x * horstretch) * (DP.x * horstretch));
skyfact = skyfact < maxskyfact ? skyfact: maxskyfact;
flred = PARC_RNDR_COLOR(3, 0) - redsky * skyfact;
flgreen = PARC_RNDR_COLOR(3, 1) - greensky * skyfact;
flblue = PARC_RNDR_COLOR(3, 2) - bluesky * skyfact;
random = settings.skyalias * drand48() - halfsky;
aliasred = flred + random * flred / 255.0;
aliasgreen = flgreen + random * flgreen / 255.0;
aliasblue = flblue + random * flblue / 255.0;
if (aliasred > 255) aliasred = 255;
if (aliasgreen > 255) aliasgreen = 255;
if (aliasblue > 255) aliasblue = 255;
red = aliasred < 0 ? 0: aliasred;
green = aliasgreen < 0 ? 0: aliasgreen;
blue = aliasblue < 0 ? 0: aliasblue;
if (((unsigned int)bytemap[zip]) > 0)
{
MergePts = 100 - (unsigned int)bytemap[zip];
colval = (*(bitmap[0] + zip) * (unsigned int)bytemap[zip] + red * MergePts) / 100;
*(bitmap[0] + zip) = (UBYTE)colval;
colval = (*(bitmap[1] + zip) * (unsigned int)bytemap[zip] + green * MergePts) / 100;
*(bitmap[1] + zip) = (UBYTE)colval;
colval = (*(bitmap[2] + zip) * (unsigned int)bytemap[zip] + blue * MergePts) / 100;
*(bitmap[2] + zip) = (UBYTE)colval;
} /* values already present */
else
{
*(bitmap[0] + zip) = (UBYTE)red;
*(bitmap[1] + zip) = (UBYTE)green;
*(bitmap[2] + zip) = (UBYTE)blue;
} /* no values present */
if ((render & 0x11) == 0x11)
{
ScreenPixelPlot(win, bitmap, x + drawoffsetX, y + drawoffsetY, zip);
}
if (y > 2)
if (zbuf[zip - 3 * settings.scrnwidth] == FLT_MAX)
{
zbuf[zip] = FLT_MAX; /* this is so partial sky pixels can receive
clouds and aren't picked up in
reflections where they don't belong */
if (ReflectionMap)
ReflectionMap[zip] = 0;
} /* if */
} /* if bytemap[zip] not full */
} /* for x=0... */
if (CheckInput_ID() == ID_BW_CLOSE)
{
error = 1;
break;
}
BusyWin_Update(BWDE, y + 1);
} /* for y=0... */
if (BWDE) BusyWin_Del(BWDE);
return (error);
} /* makesky() */
/*********************************************************************/
void getskypt(void)
{
double angle,newangle,length;
if (DP.x == 0.0) {
angle = DP.y > 0.0 ? HalfPi: OneAndHalfPi;
}
else {
angle = atan(DP.y / DP.x);
if (DP.x < 0.0) angle += Pi;
}
length = sqrt(DP.x * DP.x + DP.y * DP.y);
newangle = angle - PARC_RNDR_MOTION(8) * PiOver180;
DP.x = length * cos(newangle);
DP.y = length * sin(newangle);
}
/*********************************************************************/
void antialias(void)
{
long x, y, j, pts, zip;
if (settings.zbufalias) {
for (j=0; j<settings.aliasfactor; j++) {
for (y=0; y<high; y++) {
for (x=pts=0; x<wide; x++,pts=0) {
zip=y*settings.scrnwidth+x;
red=(UBYTE)(*(bitmap[0]+zip));
green=(UBYTE)(*(bitmap[1]+zip));
blue=(UBYTE)(*(bitmap[2]+zip));
if (fabs( *(zbuf + zip) - *(zbuf + zip + 1) ) > settings.zalias) {
altred=(UBYTE)(*(bitmap[0]+zip+1));
altgreen=(UBYTE)(*(bitmap[1]+zip+1));
altblue=(UBYTE)(*(bitmap[2]+zip+1));
pts=1;
}
if (fabs( *(zbuf + zip) - *(zbuf + zip + settings.scrnwidth)) >
settings.zalias) {
if (pts) {
aliasred=((long)red*2+altred+(UBYTE)(*(bitmap[0]+zip+settings.scrnwidth)))/4;
aliasgreen=((long)green*2+altgreen+(UBYTE)(*(bitmap[1]+zip+settings.scrnwidth)))/4;
aliasblue=((long)blue*2+altblue+(UBYTE)(*(bitmap[2]+zip+settings.scrnwidth)))/4;
}
else {
aliasred=((long)red*2+(UBYTE)(*(bitmap[0]+zip+settings.scrnwidth)))/3;
aliasgreen=((long)green*2+(UBYTE)(*(bitmap[1]+zip+settings.scrnwidth)))/3;
aliasblue=((long)blue*2+(UBYTE)(*(bitmap[2]+zip+settings.scrnwidth)))/3;
pts=1;
}
}
else {
if (pts) {
aliasred=((long)red*2+altred)/3;
aliasgreen=((long)green*2+altgreen)/3;
aliasblue=((long)blue*2+altblue)/3;
}
}
if (pts) {
*(bitmap[0]+zip)=(UBYTE)aliasred;
*(bitmap[1]+zip)=(UBYTE)aliasgreen;
*(bitmap[2]+zip)=(UBYTE)aliasblue;
}
} /* for x=0... */
} /* for y=0... */
} /* for j=0... */
} /* if settings.zbufalias */
else {
for (j=0; j<settings.aliasfactor; j++) {
for (y=0; y<high; y++) {
for (x=0; x<wide; x++) {
zip=y*settings.scrnwidth+x;
red=(UBYTE)(*(bitmap[0]+zip));
altred=((UBYTE)(*(bitmap[0]+zip+1))+(UBYTE)(*(bitmap[0]+zip+settings.scrnwidth)))/2;
green=(UBYTE)(*(bitmap[1]+zip));
altgreen=((UBYTE)(*(bitmap[1]+zip+1))+(UBYTE)(*(bitmap[1]+zip+settings.scrnwidth)))/2;
blue=(UBYTE)(*(bitmap[2]+zip));
altblue=((UBYTE)(*(bitmap[2]+zip+1))+(UBYTE)(*(bitmap[2]+zip+settings.scrnwidth)))/2;
aliasred=((int)red+(int)altred)/2;
aliasgreen=((int)green+(int)altgreen)/2;
aliasblue=((int)blue+(int)altblue)/2;
*(bitmap[0]+zip)=(UBYTE)aliasred;
*(bitmap[1]+zip)=(UBYTE)aliasgreen;
*(bitmap[2]+zip)=(UBYTE)aliasblue;
} /* for x=0... */
} /* for y=0... */
} /* for j=0 */
} /* if !settings.zbufalias */
} /* antialias() */
/**********************************************************************/
#define HOR_JITTER_AMP 1.0
#define VERT_JITTER_AMP 30.0
short Reflection_Render(struct Window *win)
{
short error = 0, Red, Grn, Blu, Banking = 0;
long x, y, sX, sY, zip, sZip, csZip;
double A, D, æ, ß, dAodX, Aprpr, Apr, Dpr, AoD, dsX, dsY, oX, oY,
Wx[3], Wy[3], dOffset, HalfWidth, BankAngle;
struct BusyWindow *BWDE;
if (! ReflectionMap || ! ElevationMap || ! bitmap[0] || ! bitmap[1] || ! bitmap[2])
return (0);
BWDE = BusyWin_New("Reflections", settings.scrnheight, 0, 'BWDE');
HalfWidth = wide / 2.0;
if (fabs(PARC_RNDR_MOTION(8)) > .0001)
{
BankAngle = PARC_RNDR_MOTION(8) * PiOver180;
Banking = 1;
} /* if gotta bank those reflections */
for (y=1, zip=settings.scrnwidth; y<settings.scrnheight; y++)
{
for (x=0; x<settings.scrnwidth; x++, zip++)
{
if (ReflectionMap[zip] && ReflectionMap[zip] != 255)
{
sX = x;
dOffset = (double)(x - HalfWidth) / (double)(y - high - 4 * settings.scrnwidth);
if (Banking)
{
dOffset = atan(dOffset) - BankAngle;
if (fabs(dOffset - HalfPi) < .0001)
continue;
dOffset = tan(dOffset);
} /* if Banking */
A = PARC_RNDR_MOTION(0) - ElevationMap[zip] * PARC_RNDR_MOTION(14) * .001;
/* Uncomment this if you don't want reflections when below the water surface
if (A < 0.0)
continue;
*/
D = sqrt(zbuf[zip] * zbuf[zip] - A * A);
if (D <= 0.0)
continue;
AoD = A / D;
if (SlopeMap)
{
dAodX = SlopeMap[zip] * zbuf[zip] / D;
ß = atan(dAodX);
æ = atan(AoD);
AoD = tan(ß * 2.0 + æ);
} /* if waves */
if (AoD > 1.0)
AoD = 1.0;
if (AoD > 0.0)
ReflectionMap[zip] *= (1.0 - AoD);
for (sY=y-1, dsX=sX-dOffset; sY>=0; sY--, dsX-=dOffset)
{
sX = dsX;
if (sX < 0)
{
dOffset *= (-1.0);
sX = 0;
}
else if (sX > wide)
{
dOffset *= (-1.0);
sX = wide;
}
sZip = sY * settings.scrnwidth + sX;
if (ElevationMap[sZip] < ElevationMap[zip] - MaxWaveAmp)
continue;
Aprpr = PARC_RNDR_MOTION(0) - ElevationMap[sZip] * PARC_RNDR_MOTION(14) * .001;
Dpr = sqrt(zbuf[sZip] * zbuf[sZip] - Aprpr * Aprpr) - D;
Apr = (ElevationMap[sZip] - ElevationMap[zip]) * PARC_RNDR_MOTION(14) * .001;
if (Dpr > 0.0)
{
if (Apr / Dpr > AoD)
break;
} /* if */
} /* for sY=... */
/* horizontal perturbation */
if (sY < 0)
{
sY = 0;
} /* if */
dsY = sY + .5;
dsX += ((HOR_JITTER_AMP / (2.0 * zbuf[zip]))
* sin((((x - frame / 5.0) * zbuf[zip]) * TwoPi) / 10.0));
sX = dsX;
if (sX <= 0)
{
sX = 0;
dsX = 0.5;
}
else if (sX >= wide)
{
sX = wide;
dsX = wide + .5;
}
/* vertical perturbation */
sY = dsY;
if (sY <= 0)
{
sY = 0;
dsY = .5;
}
else if (sY >= high)
{
sY = high;
dsY = high + .5;
}
sZip = sY * settings.scrnwidth + sX;
/* antialias the source pixel */
oY = dsY - sY;
oX = dsX - sX;
if (oY >= .5)
{
Wy[0] = 0.0;
if (sY < high)
{
Wy[1] = .5 + 1.0 - oY;
Wy[2] = 1.0 - Wy[1];
}
else
{
Wy[1] = 1.0;
Wy[2] = 0.0;
}
}
else
{
Wy[2] = 0.0;
if (sY > 0)
{
Wy[1] = .5 + oY;
Wy[0] = 1.0 - Wy[1];
}
else
{
Wy[1] = 1.0;
Wy[0] = 0.0;
}
}
if (oX >= .5)
{
Wx[0] = 0.0;
if (sX < wide)
{
Wx[1] = .5 + 1.0 - oX;
Wx[2] = 1.0 - Wx[1];
}
else
{
Wx[1] = 1.0;
Wx[2] = 0.0;
}
}
else
{
Wx[2] = 0.0;
if (sX > 0)
{
Wx[1] = .5 + oX;
Wx[0] = 1.0 - Wx[1];
}
else
{
Wx[1] = 1.0;
Wx[0] = 0.0;
}
}
Red = Grn = Blu = 0;
csZip = sZip - 1 - settings.scrnwidth;
if (csZip >= 0)
{
Red += *(bitmap[0] + csZip) * Wx[0] * Wy[0];
Grn += *(bitmap[1] + csZip) * Wx[0] * Wy[0];
Blu += *(bitmap[2] + csZip) * Wx[0] * Wy[0];
}
csZip ++;
if (csZip >= 0)
{
Red += *(bitmap[0] + csZip) * Wx[1] * Wy[0];
Grn += *(bitmap[1] + csZip) * Wx[1] * Wy[0];
Blu += *(bitmap[2] + csZip) * Wx[1] * Wy[0];
}
csZip ++;
if (csZip >= 0)
{
Red += *(bitmap[0] + csZip) * Wx[2] * Wy[0];
Grn += *(bitmap[1] + csZip) * Wx[2] * Wy[0];
Blu += *(bitmap[2] + csZip) * Wx[2] * Wy[0];
}
csZip += (settings.scrnwidth - 2);
if (csZip >= 0)
{
Red += *(bitmap[0] + csZip) * Wx[0] * Wy[1];
Grn += *(bitmap[1] + csZip) * Wx[0] * Wy[1];
Blu += *(bitmap[2] + csZip) * Wx[0] * Wy[1];
}
csZip ++;
Red += *(bitmap[0] + csZip) * Wx[1] * Wy[1];
Grn += *(bitmap[1] + csZip) * Wx[1] * Wy[1];
Blu += *(bitmap[2] + csZip) * Wx[1] * Wy[1];
csZip ++;
if (csZip < bmapsize)
{
Red += *(bitmap[0] + csZip) * Wx[2] * Wy[1];
Grn += *(bitmap[1] + csZip) * Wx[2] * Wy[1];
Blu += *(bitmap[2] + csZip) * Wx[2] * Wy[1];
}
csZip += (settings.scrnwidth - 2);
if (csZip < bmapsize)
{
Red += *(bitmap[0] + csZip) * Wx[0] * Wy[2];
Grn += *(bitmap[1] + csZip) * Wx[0] * Wy[2];
Blu += *(bitmap[2] + csZip) * Wx[0] * Wy[2];
}
csZip ++;
if (csZip < bmapsize)
{
Red += *(bitmap[0] + csZip) * Wx[1] * Wy[2];
Grn += *(bitmap[1] + csZip) * Wx[1] * Wy[2];
Blu += *(bitmap[2] + csZip) * Wx[1] * Wy[2];
}
csZip ++;
if (csZip < bmapsize)
{
Red += *(bitmap[0] + csZip) * Wx[2] * Wy[2];
Grn += *(bitmap[1] + csZip) * Wx[2] * Wy[2];
Blu += *(bitmap[2] + csZip) * Wx[2] * Wy[2];
}
/* bitmap averaging */
*(bitmap[0] + zip) = (*(bitmap[0] + zip) * (254 - ReflectionMap[zip]) +
((Red * 8) / 10) * ReflectionMap[zip]) / 254;
*(bitmap[1] + zip) = (*(bitmap[1] + zip) * (254 - ReflectionMap[zip]) +
((Grn * 8) / 10) * ReflectionMap[zip]) / 254;
*(bitmap[2] + zip) = (*(bitmap[2] + zip) * (254 - ReflectionMap[zip]) +
((Blu * 8) / 10) * ReflectionMap[zip]) / 254;
if ((render & 0x11) == 0x11)
{
ScreenPixelPlot(win, bitmap, x + drawoffsetX, y + drawoffsetY, zip);
}
} /* if reflection */
} /* for x=0... */
if (CheckInput_ID() == ID_BW_CLOSE)
{
error = 1;
break;
}
BusyWin_Update(BWDE, y + 1);
} /* for y=0... */
/* do a little blurring to smooth out reflections */
for (y=0; y<high; y++)
{
zip = y * settings.scrnwidth;
for (x=0; x<wide; x++, zip++)
{
if (ReflectionMap[zip])
{
*(bitmap[0] + zip) = ((unsigned int)(*(bitmap[0] + zip)) +
(unsigned int)(*(bitmap[0] + zip + 1)) +
(unsigned int)(*(bitmap[0] + zip + settings.scrnwidth)) +
(unsigned int)(*(bitmap[0] + zip + settings.scrnwidth + 1))) / 4;
*(bitmap[1] + zip) = ((unsigned int)(*(bitmap[1] + zip)) +
(unsigned int)(*(bitmap[1] + zip + 1)) +
(unsigned int)(*(bitmap[1] + zip + settings.scrnwidth)) +
(unsigned int)(*(bitmap[1] + zip + settings.scrnwidth + 1))) / 4;
*(bitmap[2] + zip) = ((unsigned int)(*(bitmap[2] + zip)) +
(unsigned int)(*(bitmap[2] + zip + 1)) +
(unsigned int)(*(bitmap[2] + zip + settings.scrnwidth)) +
(unsigned int)(*(bitmap[2] + zip + settings.scrnwidth + 1))) / 4;
} /* if reflection */
} /* for x=0... */
} /* for y=0... */
if (BWDE) BusyWin_Del(BWDE);
return (error);
} /* Reflection_Render() */
/***********************************************************************/
short Celestial_Bodies(UBYTE **Bitmap, long Width, long Height, struct Window *win)
{
char filename[256];
short error = 0, Colors, Sw, Sh;
float fDx, fDy, fDq;
double Dr, Luminosity, Visible;
UBYTE *Image[3];
struct ColorComponents CC;
/* Sun */
if (settings.sun)
{
Luminosity = Visible = 0.0;
Image[0] = Image[1] = Image[2] = NULL;
Sw = Sh = 0;
Dr = PARC_RNDR_MOTION(28);
CC.Red = PARC_RNDR_COLOR(19, 0);
CC.Grn = PARC_RNDR_COLOR(19, 1);
CC.Blu = PARC_RNDR_COLOR(19, 2);
DP.alt = 149.0E+6;
DP.lat = PARC_RNDR_MOTION(26);
DP.lon = PARC_RNDR_MOTION(27);
getscrncoords(&fDx, &fDy, &fDq, NULL, NULL);
strmfp(filename, imagepath, sunfile);
if (LoadImage(filename, 0, Image, Sw, Sh, 0, &Sw, &Sh, &Colors))
{
if (! (error = Image_Composite(Bitmap, Image, Width, Height,
Sw, Sh, Dr, (double)fDx, (double)fDy, (double)fDq,
&Visible, &Luminosity, &CC, "Sun", win)))
{
if (settings.sunhalo)
error = HaloEffect(Bitmap, Width, Height, Dr * 10.0, (double)fDx, (double)fDy,
Luminosity * Visible, (double)fDq, &CC, "Sun Halo", win);
} /* if */
} /* if image loaded */
else
{
error = 1;
User_Message("Render Module", "Error loading Sun Image!\n\
Operation terminated.", "OK", "o");
} /* else */
if (Image[0])
free_Memory(Image[0], Sw * Sh);
if (Image[1])
free_Memory(Image[1], Sw * Sh);
if (Image[2])
free_Memory(Image[2], Sw * Sh);
} /* if sun */
if (error)
goto EndCelest;
/* Moon */
if (settings.moon)
{
Luminosity = Visible = 0.0;
Image[0] = Image[1] = Image[2] = NULL;
Sw = Sh = 0;
Dr = PARC_RNDR_MOTION(31);
CC.Red = PARC_RNDR_COLOR(20, 0);
CC.Grn = PARC_RNDR_COLOR(20, 1);
CC.Blu = PARC_RNDR_COLOR(20, 2);
DP.alt = 38.0E+4;
DP.lat = PARC_RNDR_MOTION(29);
DP.lon = PARC_RNDR_MOTION(30);
getscrncoords(&fDx, &fDy, &fDq, NULL, NULL);
strmfp(filename, imagepath, moonfile);
if (LoadImage(filename, 0, Image, Sw, Sh, 0, &Sw, &Sh, &Colors))
{
if (! (error = Image_Composite(Bitmap, Image, Width, Height,
Sw, Sh, Dr, (double)fDx, (double)fDy, (double)fDq,
&Visible, &Luminosity, &CC, "Moon", win)))
{
if (settings.moonhalo)
error = HaloEffect(Bitmap, Width, Height, Dr * 3.0, (double)fDx, (double)fDy,
Luminosity * Visible, (double)fDq, &CC, "Moon Halo", win);
} /* if */
} /* if image loaded */
else
{
error = 1;
User_Message("Render Module", "Error loading Moon Image!\n\
Operation terminated.", "OK", "o");
} /* else */
if (Image[0])
free_Memory(Image[0], Sw * Sh);
if (Image[1])
free_Memory(Image[1], Sw * Sh);
if (Image[2])
free_Memory(Image[2], Sw * Sh);
} /* if moon */
EndCelest:
return (error);
} /* Celestial_Bodies() */
/***********************************************************************/
short HaloEffect(UBYTE **Bitmap, long Width, long Height,
double Dr, double Dx, double Dy, double Intensity, double NoHaloDist,
struct ColorComponents *CC, char *NameStr, struct Window *win)
{
short error = 0;
long x, y, zip, CheckByteMap;
double HaloStr, BGStr, Dist, Distx, Disty, DistQMax;
struct BusyWindow *BWDE;
DistQMax = qmax / 4.0;
BWDE = BusyWin_New(NameStr, Height, 0, 'BWDE');
zip = CheckByteMap = 0;
for (y=0; y<Height; y++)
{
for (x=0; x<Width; x++, zip++, CheckByteMap=0)
{
if (fabs(zbuf[zip] - NoHaloDist) < 1.0)
{
if (bytemap[zip] >= 99)
continue;
CheckByteMap = 1;
} /* if this is the object we're haloing */
Distx = Dx - x;
Disty = (Dy - y) * settings.picaspect;
Dist = sqrt(Distx * Distx + Disty * Disty);
if (Dr > Dist)
{
HaloStr = 1.0 - Dist / Dr;
HaloStr *= HaloStr;
HaloStr *= Intensity;
if (CheckByteMap)
HaloStr *= ((100.0 - (unsigned int)bytemap[zip]) / 100.0);
if (zbuf[zip] < DistQMax)
HaloStr -= (HaloStr * (DistQMax - zbuf[zip]) / DistQMax);
BGStr = 1.0 - HaloStr;
Bitmap[0][zip] = (unsigned int)Bitmap[0][zip] * BGStr + CC->Red * HaloStr;
Bitmap[1][zip] = (unsigned int)Bitmap[1][zip] * BGStr + CC->Grn * HaloStr;
Bitmap[2][zip] = (unsigned int)Bitmap[2][zip] * BGStr + CC->Blu * HaloStr;
if ((render & 0x11) == 0x11)
{
ScreenPixelPlot(win, bitmap, x + drawoffsetX, y + drawoffsetY, zip);
}
} /* if in area of effect */
} /* for x=0... */
if (CheckInput_ID() == ID_BW_CLOSE)
{
error = 1;
break;
} /* if user abort */
BusyWin_Update(BWDE, y + 1);
} /* for y=0... */
if (BWDE) BusyWin_Del(BWDE);
return (error);
} /* HaloEffect() */
/***********************************************************************/
short Image_Composite(UBYTE **Bitmap, UBYTE **Source, long Iw, long Ih,
short Sw, short Sh, double Dr, double Dx, double Dy, double Distance,
double *Visible, double *Luminosity, struct ColorComponents *CC,
char *NameStr, struct Window *win)
{
short error = 0;
double Dox, Doy, Dex, Dey, dX, dY, Sox, Soy, Cox, Coy, Cex, Cey,
wtys, wtye, wty, wtxs, wtxe, wt, PixWt, MaxWt;
long Px, Py, Pxp1, Pyp1, x, y, DRows, DCols, DxStart, DyStart, PixVal,
zip, SourceZip, i, j, PtsSeen = 0, PtsExist = 0, Wt;
struct BusyWindow *BWDE;
if (Dr <= 0.0 || Sw == 0 || Sh == 0)
{
*Visible = *Luminosity = 0.0;
return (0);
} /* if something not kosher */
Dox = Dx - Dr;
Dex = Dx + Dr;
Doy = Dy - Dr / settings.picaspect;
Dey = Dy + Dr / settings.picaspect;
dX = Sw / (Dex - Dox);
dY = Sh / (Dey - Doy);
MaxWt = dX * dY;
Sox = ((int)Dox - Dox) * dX;
Soy = ((int)Doy - Doy) * dY;
DxStart = Dox;
DyStart = Doy;
DCols = Dex - DxStart + 1.0;
DRows = Dey - DyStart + 1.0;
BWDE = BusyWin_New(NameStr, DRows, 0, 'BWDE');
for (y=DyStart, Coy=Soy, Cey=Soy+dY, j=0, PixWt = 0, PixVal=0.0; j<DRows;
j++, y++, Coy+=dY, Cey+=dY)
{
if (y < 0)
continue;
if (y >= Ih)
break;
zip = y * Iw + DxStart;
for (x=DxStart, Cox=Sox, Cex=Sox+dX, i=0; i<DCols;
i++, x++, Cox+=dX, Cex+=dX, PixVal=0, PixWt=0.0, zip++)
{
if (x < 0)
continue;
if (x >= Iw)
break;
for (Py=Coy, Pyp1=Coy+1; Py<Cey && Py<Sh; Py++, Pyp1++)
{
if (Py < 0)
continue;
wtys = Py > Coy ? 1.0: Pyp1 - Coy;
wtye = Pyp1 < Cey ? 1.0: Cey - Py;
wty = wtys * wtye;
for (Px=Cox, Pxp1=Cox+1; Px<Cex && Px<Sw; Px++, Pxp1++)
{
if (Px < 0)
continue;
wtxs = Px > Cox ? 1.0: Pxp1 - Cox;
wtxe = Pxp1 < Cex ? 1.0: Cex - Px;
wt = wty * wtxs * wtxe;
SourceZip = Py * Sw + Px;
if (Source[0][SourceZip])
{
PixWt += wt;
PixVal += wt * (unsigned int)Source[0][SourceZip];
}
} /* for Px=... */
} /* for Py=... */
if (PixVal && PixWt > 0.0)
{
PtsExist ++;
if (Distance < zbuf[zip])
{
PtsSeen ++;
PixVal /= PixWt;
*Luminosity += PixVal;
PixWt /= MaxWt;
zbuf[zip] = Distance;
bytemap[zip] = 100 * PixWt;
Wt = PixWt * PixVal;
Bitmap[0][zip] += ((Wt * (CC->Red - Bitmap[0][zip])) / 255);
Bitmap[1][zip] += ((Wt * (CC->Grn - Bitmap[1][zip])) / 255);
Bitmap[2][zip] += ((Wt * (CC->Blu - Bitmap[2][zip])) / 255);
if ((render & 0x11) == 0x11)
{
ScreenPixelPlot(win, bitmap, x + drawoffsetX, y + drawoffsetY, zip);
} /* if */
} /* if */
} /* if */
} /* for Tox=... */
if (CheckInput_ID() == ID_BW_CLOSE)
{
error = 1;
break;
} /* if user abort */
BusyWin_Update(BWDE, j + 1);
} /* for Toy=... */
*Visible = PtsExist ? (double)PtsSeen / (double)PtsExist: 0.0;
*Luminosity = PtsSeen ? (*Luminosity / PtsSeen) / 255.0: 0.0;
if (BWDE) BusyWin_Del(BWDE);
return (error);
} /* Image_Composite() */Appendix B: MapTopo.c
The complete file, converted from Latin-1 to UTF-8. The wave slope is computed at lines 891 to 894 and the per-pixel buffers are written at lines 1052 to 1059. On GitHub: current version, pinned version. Licensed under the OSGPL.
Show the full source (1,646 lines)
/* MapTopo.c (ne gismaptopo.c 14 Jan 1994 CXH)
** Maps a Topo polygon
** Grabbed/Built from gisam.c on 24 Jul 1993 by Chris "Xenon" Hanson.
** Original code by Gary R. Huber.
*/
#include "WCS.h"
#define ZBLEND
#define MODE_REPLACE 0
#define MODE_AVERAGE 1
void MapTopo(struct elmapheaderV101 *map, struct Window *win, short MapAsSFC,
short MakeWater, short Visible, double *Elev)
{
short eco = 0, understory = 0, notsnow = 0, NoAlias = 0, VolumeTexture = 0,
Texture, QInterp = 1, ElInterp = 1, i, j, k, l, sum, Mode;
long NewSubPixArraySize, SubPixWidth, SubPixHeight, MinX, MaxX, StartY, EndY,
colval, x, y, zip, WtAvg;
/* changed from float to double 10/19/95 */
double AAX[3], AAY[3], Xoffset, dist;
double ecoline, snowline, m, h, d;
struct QCvalues QC;
UBYTE *PixPtr;
struct ColorComponents CC[4];
double fred, fgrn, fblu;
short PolyMinXPt, PolyMaxXPt, Tex, MinPt, MidPt, MaxPt, MinQPt, MaxQPt;
double ElX, ElY, LatX, LatY, LonX, LonY, QX, QY,
ElStart, LatStart, LonStart, QStart, ElPt, LatPt, LonPt, QPt,
Temp, X5, Y5, Y4, X3, El3, EloQY;
/* initialize */
Reflections = 0; /* enabled in WaterEco_Set() */
treedraw = undertreedraw = 0; /* enabled in ecoset() */
if (OldLightingModel)
{
h = sunlong - facelong * PiOver180 - diplong;
if (h > Pi) h -= TwoPi;
else if (h < -Pi) h += TwoPi;
d = sunlat - facelat * PiOver180 - diplat;
if (d > Pi) d -= TwoPi;
else if (d < -Pi) d += TwoPi;
sunangle = sqrt(h * h + d * d) * .7071;
if (sunangle > HalfPi) sunangle = HalfPi;
sunfactor = 1.0 - cos(sunangle);
sunshade = sunfactor * PARC_RNDR_MOTION(22);
sunshade += ((1.0 - sunshade) * cloudcover);
} /* if */
/* compute fade and fog factors */
fade = (qqq - PARC_RNDR_MOTION(20)) / PARC_RNDR_MOTION(21);
if (fade > 1.0) fade = 1.0;
else if (fade < 0.0) fade = 0.0;
if (fogrange == 0.0) fog = 0.0;
else
{
fog = (el - PARC_RNDR_MOTION(23)) / fogrange;
if (fog > 1.0) fog = 1.0;
else if (fog < 0.0) fog = 0.0;
} /* else */
/* if surface */
if (MapAsSFC)
{
short range;
double gradfact;
if (settings.reliefshade || (render & 0x01))
{
FloatCol = 8.0 + 7.99 * sunfactor + ROUNDING_KLUDGE;
ColMax = COL_SFC_MAX;
} /* if shading by sun angle */
if ((! settings.reliefshade) || (render & 0x01))
{
if (el < settings.surfel[1])
{
gradfact = ((double)settings.surfel[1] - el)
/ (settings.surfel[1] - settings.surfel[0]);
range = 0;
} /* low gradient */
else if (el >= settings.surfel[1] && el < settings.surfel[2])
{
gradfact = ((double)settings.surfel[2] - el)
/ (settings.surfel[2] - settings.surfel[1]);
range = 1;
} /* else if middle gradient */
else
{
gradfact = ((double)settings.surfel[3] - el)
/ (settings.surfel[3] - settings.surfel[2]);
range = 2;
} /* else upper gradient */
if (gradfact > 1.0) gradfact = 1.0;
else if (gradfact < 0.0) gradfact = 0.0;
FloatCol = 8.0 + 7.99 * gradfact + ROUNDING_KLUDGE;
ColMax = COL_SFC_MAX;
} /* if */
if (render & 0x01)
{
switch (range)
{
case 0:
{
CC[0].Red = PARC_RNDR_COLOR(7, 0);
CC[0].Red += (PARC_RNDR_COLOR(6, 0) - CC[0].Red) * gradfact;
CC[0].Grn = PARC_RNDR_COLOR(7, 1);
CC[0].Grn += (PARC_RNDR_COLOR(6, 1) - CC[0].Grn) * gradfact;
CC[0].Blu = PARC_RNDR_COLOR(7, 2);
CC[0].Blu += (PARC_RNDR_COLOR(6, 2) - CC[0].Blu) * gradfact;
break;
} /* low range */
case 1:
{
CC[0].Red = PARC_RNDR_COLOR(8, 0);
CC[0].Red += (PARC_RNDR_COLOR(7, 0) - CC[0].Red) * gradfact;
CC[0].Grn = PARC_RNDR_COLOR(8, 1);
CC[0].Grn += (PARC_RNDR_COLOR(7, 1) - CC[0].Grn) * gradfact;
CC[0].Blu = PARC_RNDR_COLOR(8, 2);
CC[0].Blu += (PARC_RNDR_COLOR(7, 2) - CC[0].Blu) * gradfact;
break;
} /* middle range */
case 2:
{
CC[0].Red = PARC_RNDR_COLOR(9, 0);
CC[0].Red += (PARC_RNDR_COLOR(8, 0) - CC[0].Red) * gradfact;
CC[0].Grn = PARC_RNDR_COLOR(9, 1);
CC[0].Grn += (PARC_RNDR_COLOR(8, 1) - CC[0].Grn) * gradfact;
CC[0].Blu = PARC_RNDR_COLOR(9, 2);
CC[0].Blu += (PARC_RNDR_COLOR(8, 2) - CC[0].Blu) * gradfact;
break;
} /* high range */
} /* switch color range */
CC[0].Red -= CC[0].Red * sunshade;
CC[0].Grn -= CC[0].Grn * sunshade;
CC[0].Blu -= CC[0].Blu * sunshade;
} /* if render to RGB */
eco = understory = 50;
goto EndDrawFace;
} /* if surface */
/* compute the shading factors for normal (not water) ecosystems */
relfactor = relel + 100.0 * random;
/* see if there is snow on the ground */
snowline = PARC_RNDRLN_ECO(1) +
diplong * PARC_SKLN_ECO(1) + diplat * PARC_SKLT_ECO(1)
+ PARC_RNDRRE_ECO(1) * relfactor;
if (settings.worldmap)
snowline -= ((abs(facelat) - abs(settings.globreflat)) * settings.globecograd);
if (settings.flatteneco)
snowline += ( (PARC_RNDR_MOTION(13) - snowline) * PARC_RNDR_MOTION(12) );
if (el >= snowline)
{
if (slope >= PARC_MNSL_ECO(1) && slope <= PARC_MXSL_ECO(1))
{
if (relfactor >= PARC_RNDRNR_ECO(1)
&& relfactor <= PARC_RNDRXR_ECO(1))
{
FloatCol = 1.0 + 4.99 * sunfactor;
ColMax = COL_SNOW_MAX;
CC[0].Red = PARC_MCOL_ECO(1, 0);
CC[0].Red -= pow(sunshade, 3.3) * CC[0].Red * PARC_RNDR_MOTION(22);
CC[0].Grn = PARC_MCOL_ECO(1, 1);
CC[0].Grn -= pow(sunshade, 3.4) * CC[0].Grn * PARC_RNDR_MOTION(22);
CC[0].Blu = PARC_MCOL_ECO(1, 2);
CC[0].Blu -= pow(sunshade, 3.8) * CC[0].Blu * PARC_RNDR_MOTION(22);
} /* if relel matchup */
else notsnow = 1;
} /* if slope matchup */
else notsnow = 1;
} /* if el > snowline */
else notsnow = 1;
/* We'll add some randomness to our colors for other than snow.
** Since "random" varies from -.1 to +.1, this will give us a variation of
** ±20 */
switch (dir)
{
case 0:
redrand = random * 200.0;
greenrand = random * 200.0;
bluerand = 0;
break;
case 1:
redrand = 0;
greenrand = random * 200.0;
bluerand = random * 200.0;
break;
case 2:
redrand = random * 200.0;
greenrand = 0;
bluerand = random * 200.0;
break;
case 3:
redrand = 0;
greenrand = random * 200.0;
bluerand = 0;
break;
} /* switch */
/* if below water level, make it wet */
if (MakeWater || el <= SeaLevel)
{
eco = WaterEco_Set(MakeWater, CC);
understory = PAR_UNDER_ECO(eco);
ecocount[eco] ++;
goto EndDrawFace;
} /* if make water */
/* Do the colormap thing. If it finds a result >= 0, we're outa here! */
if (cmap)
{
if ((eco = colormap(map, notsnow, CC, &understory)) >= 0)
{
ecocount[eco] ++;
goto EndDrawFace;
}
} /* if cmap */
/* We now know whether there is snow on the ground and that this is not
** a color-mapped polygon. Now search for an ecosystem match between
** the physical terrain conditions and the ecosystem parameter list */
for (eco=12; eco<ECOPARAMS; eco++)
{
ecoline = PARC_RNDRLN_ECO(eco) +
diplong * PARC_SKLN_ECO(eco) + diplat * PARC_SKLT_ECO(eco)
+ PARC_RNDRRE_ECO(eco) * relfactor;
if (settings.worldmap)
ecoline -= ((abs(facelat) - abs(settings.globreflat)) * settings.globecograd);
if (settings.flatteneco)
ecoline += ( (PARC_RNDR_MOTION(13) - ecoline) * PARC_RNDR_MOTION(12) );
if (el <= ecoline)
{
if (slope >= PARC_MNSL_ECO(eco) && slope <= PARC_MXSL_ECO(eco))
{
if (relfactor >= PARC_RNDRNR_ECO(eco) &&
relfactor <= PARC_RNDRXR_ECO(eco))
{
understory = ecoset(eco, notsnow, CC);
break;
} /* if relel matchup */
} /* if slope matchup */
} /* if elevation matchup */
} /* for eco=12... */
if (eco >= ECOPARAMS)
{
Log(WNG_ILL_VAL, "Ecosystem out of range.");
eco = settings.defaulteco;
understory = ecoset(eco, notsnow, CC);
} /* if */
ecocount[eco] ++;
/* Now hopefully, we have colors in all the CC arrays we need:
** CC[0] = the color of the bare ground (the understory's understory color
** unless there is snow).
** CC[1] = the understory texture color (only used if "undertreedraw"
** is 1 or more. It could be trees or rock texture, grass or whatever.
** CC[2] = the overstory color if there are trees or textures in the overstory.
** If the understory and overstory are the same ecosystem, this one will not
** be drawn since it would be rather redundant.
** We will now draw the plain polygon, unadorned by any bitmapped textures.
** If the ecosystem understory class is strata rock then the polygons will
** be given the strata treatment.
*/
EndDrawFace:
if (Visible)
{
if (render & 0x100)
{
long ct;
QC.compval1 = (long)(min(el, 32767.0)) * 65536 + eco * 256 + understory;
ct = aspect < 0.0 ? 0: aspect * PiUnder180 - 90.0;
if (ct < 0) ct += 360;
QC.compval2 = ((long)relfactor + 1000) * 65536 + ct;
QC.compval3 = (short)(slope * PiUnder180) * 256 + (short)(sunangle * PiUnder180);
} /* if render to QC buffer */
if (render & 0x01)
{
CC[3].Red = sunshade * PARC_RNDR_COLOR(1, 0); /* ambient pre-calc */
CC[3].Grn = sunshade * PARC_RNDR_COLOR(1, 1);
CC[3].Blu = sunshade * PARC_RNDR_COLOR(1, 2);
/* This is now done on a pixel basis to facilitate texture mapping
CC[0].Red += sunshade * PARC_RNDR_COLOR(1, 0); /* ambient */
CC[0].Grn += sunshade * PARC_RNDR_COLOR(1, 1);
CC[0].Blu += sunshade * PARC_RNDR_COLOR(1, 2);
CC[0].Red += (PARC_RNDR_COLOR(2, 0) - CC[0].Red) * fade;
CC[0].Red += (PARC_RNDR_COLOR(2, 0) - CC[0].Red) * fog;
CC[0].Red *= redsun;
CC[0].Grn += (PARC_RNDR_COLOR(2, 1) - CC[0].Grn) * fade;
CC[0].Grn += (PARC_RNDR_COLOR(2, 1) - CC[0].Grn) * fog;
CC[0].Grn *= greensun;
CC[0].Blu += (PARC_RNDR_COLOR(2, 2) - CC[0].Blu) * fade;
CC[0].Blu += (PARC_RNDR_COLOR(2, 2) - CC[0].Blu) * fog;
CC[0].Blu *= bluesun;
if (CC[0].Red > 255) red = 255;
else red = (CC[0].Red < 0) ? 0: CC[0].Red;
if (CC[0].Grn > 255) green = 255;
else green = CC[0].Grn < 0 ? 0: CC[0].Grn;
if (CC[0].Blu > 255) blue = 255;
else blue = CC[0].Blu < 0 ? 0: CC[0].Blu;
*/
} /* if render to RGB */
/* determine which poly corner is largest and smallest along x axis */
/* temporary switch - should also be enabled in settings editor */
Texture = PAR_TYPE_ECO(PAR_UNDER_ECO(understory));
if (! MapAsSFC)
{
if (((Texture & 0x00ff) >= 100)
&& ((Texture & 0x00ff) < 150)
&& (PARC_RNDRDN_ECO(PAR_UNDER_ECO(understory)) > treerand2 * 100.0)) /* if procedural texture */
VolumeTexture = 1;
}
if (polyx[b][2] > polyx[b][0])
{
if (polyx[b][2] > polyx[b][1])
{
PolyMaxXPt = 2;
PolyMinXPt = polyx[b][1] < polyx[b][0] ? 1: 0;
} /* if */
else
{
PolyMaxXPt = 1;
PolyMinXPt = polyx[b][2] < polyx[b][0] ? 2: 0;
} /* else */
} /* if */
else
{
if (polyx[b][0] > polyx[b][1])
{
PolyMaxXPt = 0;
PolyMinXPt = polyx[b][2] < polyx[b][1] ? 2: 1;
} /* if */
else
{
PolyMaxXPt = 1;
PolyMinXPt = polyx[b][2] < polyx[b][0] ? 2: 0;
} /* else */
} /* else */
if (polyy[b][0] < 0.0) yy[0] --;
if (polyy[b][1] < 0.0) yy[1] --;
if (polyy[b][2] < 0.0) yy[2] --;
if (polyx[b][0] < 0.0) xx[0] --;
if (polyx[b][1] < 0.0) xx[1] --;
if (polyx[b][2] < 0.0) xx[2] --;
MinX = xx[PolyMinXPt];
MaxX = xx[PolyMaxXPt];
if (MaxX < 0 || MinX > wide || yy[0] > high || yy[2] < 0)
goto DrawTrees;
SubPixWidth = 10 * (1 + MaxX - MinX);
SubPixHeight = 10 * (1 + yy[2] - yy[0]);
NewSubPixArraySize = SubPixWidth * SubPixHeight;
if (NewSubPixArraySize > SubPixArraySize)
{
free_Memory(SubPix, SubPixArraySize);
SubPixArraySize = 0;
if ((SubPix = get_Memory(NewSubPixArraySize, MEMF_CLEAR)) == NULL)
{
NoAlias = 1;
}
else
{
SubPixArraySize = NewSubPixArraySize;
}
} /* if need larger array */
else
{
memset(SubPix, 0, NewSubPixArraySize);
} /* else clear as much of array as needed */
if (EdgeSize < SubPixHeight * sizeof (short))
{
free_Memory(Edge1, EdgeSize);
free_Memory(Edge2, EdgeSize);
EdgeSize = 0;
Edge1 = (short *)get_Memory(SubPixHeight * sizeof (short), MEMF_CLEAR);
Edge2 = (short *)get_Memory(SubPixHeight * sizeof (short), MEMF_CLEAR);
if (! Edge1 || ! Edge2)
{
NoAlias = 1;
}
else
{
EdgeSize = SubPixHeight * sizeof (short);
}
} /* if need new edge arrays */
if (! NoAlias)
{
AAX[0] = 10.0 * (polyx[b][0] - MinX);
AAX[1] = 10.0 * (polyx[b][1] - MinX);
AAY[0] = 10.0 * (polyy[b][0] - yy[0]);
AAY[1] = 10.0 * (polyy[b][1] - yy[0]);
AAY[2] = 10.0 * (polyy[b][2] - yy[0]);
if (polyy[b][0] == polyy[b][2]) m = 0.0;
else m = (polyx[b][2] - polyx[b][0]) / (polyy[b][2] - polyy[b][0]);
StartY = AAY[0] - (int)AAY[0] < .5 ? (int)AAY[0]: 1 + (int)AAY[0];
EndY = AAY[2] - (int)AAY[2] > .5 ? (int)AAY[2]: (int)AAY[2] - 1;
Xoffset = AAX[0] + (StartY + .5 - AAY[0]) * m;
for (y=StartY; y<=EndY; y++)
{
Edge1[y] = (int)Xoffset;
Xoffset += m;
}
if (polyy[b][0] == polyy[b][1]) m = 0.0;
else m = (polyx[b][1] - polyx[b][0]) / (polyy[b][1] - polyy[b][0]);
EndY = AAY[1] - (int)AAY[1] > .5 ? (int)AAY[1]: (int)AAY[1] - 1;
Xoffset = AAX[0] + (StartY + .5 - AAY[0]) * m;
for (y=StartY; y<=EndY; y++)
{
Edge2[y] = (int)Xoffset;
Xoffset += m;
}
PixPtr = &SubPix[StartY * SubPixWidth];
for (y=StartY; y<=EndY; y++)
{
if (Edge2[y] >= Edge1[y])
{
for (x=Edge1[y]; x<=Edge2[y]; x++)
{
PixPtr[x] = 0x01;
}
}
else
{
for (x=Edge2[y]; x<=Edge1[y]; x++)
{
PixPtr[x] = 0x01;
}
}
PixPtr += SubPixWidth;
}
if (polyy[b][2] == polyy[b][1]) m = 0.0;
else m = (polyx[b][2] - polyx[b][1]) / (polyy[b][2] - polyy[b][1]);
StartY = AAY[1] - (int)AAY[1] < .5 ? (int)AAY[1]: 1 + (int)AAY[1];
EndY = AAY[2] - (int)AAY[2] > .5 ? (int)AAY[2]: (int)AAY[2] - 1;
Xoffset = AAX[1] + (StartY + .5 - AAY[1]) * m;
for (y=StartY; y<=EndY; y++)
{
Edge2[y] = (int)Xoffset;
Xoffset += m;
}
PixPtr = &SubPix[StartY * SubPixWidth];
for (y=StartY; y<=EndY; y++)
{
if (Edge2[y] >= Edge1[y])
{
for (x=Edge1[y]; x<=Edge2[y]; x++)
{
PixPtr[x] = 0x01;
}
}
else
{
for (x=Edge2[y]; x<=Edge1[y]; x++)
{
PixPtr[x] = 0x01;
}
}
PixPtr += SubPixWidth;
}
/* figure x and y gradients across polygon for lat, lon, Z and elev */
/* latitude */
if (polylat[b][2] > polylat[b][0])
{
if (polylat[b][2] > polylat[b][1])
{
MaxPt = 2;
MinPt = polylat[b][1] < polylat[b][0] ? 1: 0;
MidPt = MinPt == 1 ? 0: 1;
} /* if */
else
{
MaxPt = 1;
MinPt = polylat[b][2] < polylat[b][0] ? 2: 0;
MidPt = MinPt == 2 ? 0: 2;
} /* else */
} /* if */
else
{
if (polylat[b][0] > polylat[b][1])
{
MaxPt = 0;
MinPt = polylat[b][2] < polylat[b][1] ? 2: 1;
MidPt = MinPt == 2 ? 1: 2;
} /* if */
else
{
MaxPt = 1;
MinPt = polylat[b][2] < polylat[b][0] ? 2: 0;
MidPt = MinPt == 2 ? 0: 2;
} /* else */
} /* else */
if (polylat[b][MaxPt] == polylat[b][MinPt])
{
LatX = 0.0;
LatY = 0.0;
} /* that was easy */
else
{
Temp = (polylat[b][MidPt] - polylat[b][MinPt]) /
(polylat[b][MaxPt] - polylat[b][MinPt]);
X5 = Temp * (polyx[b][MaxPt] - polyx[b][MinPt]) + polyx[b][MinPt];
Y5 = Temp * (polyy[b][MaxPt] - polyy[b][MinPt]) + polyy[b][MinPt];
if (polyy[b][MaxPt] == polyy[b][MinPt])
{
LatY = LatX = 0.0;
LatStart = facelat;
goto StartLon;
/* X3 = Y4 = infinity;*/
} /* if */
X3 = polyx[b][MinPt] + (polyy[b][MidPt] - polyy[b][MinPt]) *
(polyx[b][MaxPt] - polyx[b][MinPt]) /
(polyy[b][MaxPt] - polyy[b][MinPt]);
if (X5 == polyx[b][MidPt])
{
LatY = LatX = 0.0;
LatStart = facelat;
goto StartLon;
/* Y4 = infinity;*/
} /* if */
Y4 = polyy[b][MidPt] + (X3 - polyx[b][MidPt]) *
(Y5 - polyy[b][MidPt]) / (X5 - polyx[b][MidPt]);
if (polyx[b][MaxPt] != polyx[b][MinPt])
El3 = polylat[b][MinPt] + (X3 - polyx[b][MinPt]) *
(polylat[b][MaxPt] - polylat[b][MinPt]) /
(polyx[b][MaxPt] - polyx[b][MinPt]);
else if (polyy[b][MaxPt] != polyy[b][MinPt])
El3 = polylat[b][MinPt] + (polyy[b][MidPt] - polyy[b][MinPt]) *
(polylat[b][MaxPt] - polylat[b][MinPt]) /
(polyy[b][MaxPt] - polyy[b][MinPt]);
else
{
LatY = LatX = 0.0;
LatStart = facelat;
goto StartLon;
} /* else */
if (X3 != polyx[b][MidPt] && Y4 != polyy[b][MidPt])
{
LatX = (El3 - polylat[b][MidPt]) / (X3 - polyx[b][MidPt]);
LatY = (El3 - polylat[b][MidPt]) / (polyy[b][MidPt] - Y4);
} /* if */
else
{
LatY = LatX = 0.0;
LatStart = facelat;
goto StartLon;
} /* else */
} /* else not so easy */
LatStart = polylat[b][MidPt] + (polyx[b][PolyMinXPt] - polyx[b][MidPt])
* LatX;
LatStart += (polyy[b][0] - polyy[b][MidPt]) * LatY;
/* longitude */
StartLon:
if (polylon[b][2] > polylon[b][0])
{
if (polylon[b][2] > polylon[b][1])
{
MaxPt = 2;
MinPt = polylon[b][1] < polylon[b][0] ? 1: 0;
MidPt = MinPt == 1 ? 0: 1;
} /* if */
else
{
MaxPt = 1;
MinPt = polylon[b][2] < polylon[b][0] ? 2: 0;
MidPt = MinPt == 2 ? 0: 2;
} /* else */
} /* if */
else
{
if (polylon[b][0] > polylon[b][1])
{
MaxPt = 0;
MinPt = polylon[b][2] < polylon[b][1] ? 2: 1;
MidPt = MinPt == 2 ? 1: 2;
} /* if */
else
{
MaxPt = 1;
MinPt = polylon[b][2] < polylon[b][0] ? 2: 0;
MidPt = MinPt == 2 ? 0: 2;
} /* else */
} /* else */
if (polylon[b][MaxPt] == polylon[b][MinPt])
{
LonX = 0.0;
LonY = 0.0;
} /* that was easy */
else
{
Temp = (polylon[b][MidPt] - polylon[b][MinPt]) /
(polylon[b][MaxPt] - polylon[b][MinPt]);
X5 = Temp * (polyx[b][MaxPt] - polyx[b][MinPt]) + polyx[b][MinPt];
Y5 = Temp * (polyy[b][MaxPt] - polyy[b][MinPt]) + polyy[b][MinPt];
if (polyy[b][MaxPt] == polyy[b][MinPt])
{
LonY = LonX = 0.0;
LonStart = facelong;
goto StartZ;
/* X3 = Y4 = infinity;*/
} /* if */
X3 = polyx[b][MinPt] + (polyy[b][MidPt] - polyy[b][MinPt]) *
(polyx[b][MaxPt] - polyx[b][MinPt]) /
(polyy[b][MaxPt] - polyy[b][MinPt]);
if (X5 == polyx[b][MidPt])
{
LonY = LonX = 0.0;
LonStart = facelong;
goto StartZ;
/* Y4 = infinity;*/
} /* if */
Y4 = polyy[b][MidPt] + (X3 - polyx[b][MidPt]) *
(Y5 - polyy[b][MidPt]) / (X5 - polyx[b][MidPt]);
if (polyx[b][MaxPt] != polyx[b][MinPt])
El3 = polylon[b][MinPt] + (X3 - polyx[b][MinPt]) *
(polylon[b][MaxPt] - polylon[b][MinPt]) /
(polyx[b][MaxPt] - polyx[b][MinPt]);
else if (polyy[b][MaxPt] != polyy[b][MinPt])
El3 = polylon[b][MinPt] + (polyy[b][MidPt] - polyy[b][MinPt]) *
(polylon[b][MaxPt] - polylon[b][MinPt]) /
(polyy[b][MaxPt] - polyy[b][MinPt]);
else
{
LonY = LonX = 0.0;
LonStart = facelong;
goto StartZ;
} /* else */
if (X3 != polyx[b][MidPt] && Y4 != polyy[b][MidPt])
{
LonX = (El3 - polylon[b][MidPt]) / (X3 - polyx[b][MidPt]);
LonY = (El3 - polylon[b][MidPt]) / (polyy[b][MidPt] - Y4);
} /* if */
else
{
LonY = LonX = 0.0;
LonStart = facelong;
goto StartZ;
} /* else */
} /* else not so easy */
LonStart = polylon[b][MidPt] + (polyx[b][PolyMinXPt] - polyx[b][MidPt])
* LonX;
LonStart += (polyy[b][0] - polyy[b][MidPt]) * LonY;
StartZ:
/* z buffer */
if (polyq[b][2] > polyq[b][0])
{
if (polyq[b][2] > polyq[b][1])
{
MaxQPt = 2;
MinQPt = polyq[b][1] < polyq[b][0] ? 1: 0;
MidPt = MinQPt == 1 ? 0: 1;
} /* if */
else
{
MaxQPt = 1;
MinQPt = polyq[b][2] < polyq[b][0] ? 2: 0;
MidPt = MinQPt == 2 ? 0: 2;
} /* else */
} /* if */
else
{
if (polyq[b][0] > polyq[b][1])
{
MaxQPt = 0;
MinQPt = polyq[b][2] < polyq[b][1] ? 2: 1;
MidPt = MinQPt == 2 ? 1: 2;
} /* if */
else
{
MaxQPt = 1;
MinQPt = polyq[b][2] < polyq[b][0] ? 2: 0;
MidPt = MinQPt == 2 ? 0: 2;
} /* else */
} /* else */
if (polyq[b][MaxQPt] == polyq[b][MinQPt])
{
QX = 0.0;
QY = 0.0;
} /* that was easy */
else
{
Temp = (polyq[b][MidPt] - polyq[b][MinQPt]) /
(polyq[b][MaxQPt] - polyq[b][MinQPt]);
X5 = Temp * (polyx[b][MaxQPt] - polyx[b][MinQPt]) + polyx[b][MinQPt];
Y5 = Temp * (polyy[b][MaxQPt] - polyy[b][MinQPt]) + polyy[b][MinQPt];
if (polyy[b][MaxQPt] == polyy[b][MinQPt])
{
QInterp = 0;
QY = QX = 0.0;
QStart = qqq;
goto StartEl;
/* X3 = Y4 = infinity;*/
} /* if */
X3 = polyx[b][MinQPt] + (polyy[b][MidPt] - polyy[b][MinQPt]) *
(polyx[b][MaxQPt] - polyx[b][MinQPt]) /
(polyy[b][MaxQPt] - polyy[b][MinQPt]);
if (X5 == polyx[b][MidPt])
{
QInterp = 0;
QY = QX = 0.0;
QStart = qqq;
goto StartEl;
/* Y4 = infinity;*/
} /* if */
Y4 = polyy[b][MidPt] + (X3 - polyx[b][MidPt]) *
(Y5 - polyy[b][MidPt]) / (X5 - polyx[b][MidPt]);
if (polyx[b][MaxQPt] != polyx[b][MinQPt])
El3 = polyq[b][MinQPt] + (X3 - polyx[b][MinQPt]) *
(polyq[b][MaxQPt] - polyq[b][MinQPt]) /
(polyx[b][MaxQPt] - polyx[b][MinQPt]);
else if (polyy[b][MaxQPt] != polyy[b][MinQPt])
El3 = polyq[b][MinQPt] + (polyy[b][MidPt] - polyy[b][MinQPt]) *
(polyq[b][MaxQPt] - polyq[b][MinQPt]) /
(polyy[b][MaxQPt] - polyy[b][MinQPt]);
else
{
QInterp = 0;
QY = QX = 0.0;
QStart = qqq;
goto StartEl;
} /* else */
if (X3 != polyx[b][MidPt] && Y4 != polyy[b][MidPt])
{
QX = (El3 - polyq[b][MidPt]) / (X3 - polyx[b][MidPt]);
QY = (El3 - polyq[b][MidPt]) / (polyy[b][MidPt] - Y4);
} /* if */
else
{
QInterp = 0;
QY = QX = 0.0;
QStart = qqq;
goto StartEl;
} /* else */
} /* else not so easy */
QStart = polyq[b][MidPt] + (polyx[b][PolyMinXPt] - polyx[b][MidPt])
* QX;
QStart += (polyy[b][0] - polyy[b][MidPt]) * QY;
/* elevation */
StartEl:
if (Elev[2] > Elev[0])
{
if (Elev[2] > Elev[1])
{
MaxPt = 2;
MinPt = Elev[1] < Elev[0] ? 1: 0;
MidPt = MinPt == 1 ? 0: 1;
} /* if */
else
{
MaxPt = 1;
MinPt = Elev[2] < Elev[0] ? 2: 0;
MidPt = MinPt == 2 ? 0: 2;
} /* else */
} /* if */
else
{
if (Elev[0] > Elev[1])
{
MaxPt = 0;
MinPt = Elev[2] < Elev[1] ? 2: 1;
MidPt = MinPt == 2 ? 1: 2;
} /* if */
else
{
MaxPt = 1;
MinPt = Elev[2] < Elev[0] ? 2: 0;
MidPt = MinPt == 2 ? 0: 2;
} /* else */
} /* else */
if (Elev[MaxPt] == Elev[MinPt])
{
ElX = 0.0;
ElY = 0.0;
} /* that was easy */
else
{
Temp = (Elev[MidPt] - Elev[MinPt]) /
(Elev[MaxPt] - Elev[MinPt]);
X5 = Temp * (polyx[b][MaxPt] - polyx[b][MinPt]) + polyx[b][MinPt];
Y5 = Temp * (polyy[b][MaxPt] - polyy[b][MinPt]) + polyy[b][MinPt];
if (polyy[b][MaxPt] == polyy[b][MinPt])
{
ElInterp = 0;
ElY = ElX = 0.0;
ElStart = el;
goto StartDraw;
/* X3 = Y4 = infinity;*/
} /* if */
X3 = polyx[b][MinPt] + (polyy[b][MidPt] - polyy[b][MinPt]) *
(polyx[b][MaxPt] - polyx[b][MinPt]) /
(polyy[b][MaxPt] - polyy[b][MinPt]);
if (X5 == polyx[b][MidPt])
{
ElInterp = 0;
ElY = ElX = 0.0;
ElStart = el;
goto StartDraw;
/* Y4 = infinity;*/
} /* if */
Y4 = polyy[b][MidPt] + (X3 - polyx[b][MidPt]) *
(Y5 - polyy[b][MidPt]) / (X5 - polyx[b][MidPt]);
if (polyx[b][MaxPt] != polyx[b][MinPt])
El3 = Elev[MinPt] + (X3 - polyx[b][MinPt]) *
(Elev[MaxPt] - Elev[MinPt]) /
(polyx[b][MaxPt] - polyx[b][MinPt]);
else if (polyy[b][MaxPt] != polyy[b][MinPt])
El3 = Elev[MinPt] + (polyy[b][MidPt] - polyy[b][MinPt]) *
(Elev[MaxPt] - Elev[MinPt]) /
(polyy[b][MaxPt] - polyy[b][MinPt]);
else
{
ElInterp = 0;
ElY = ElX = 0.0;
ElStart = el;
goto StartDraw;
} /* else */
if (X3 != polyx[b][MidPt] && Y4 != polyy[b][MidPt])
{
ElX = (El3 - Elev[MidPt]) / (X3 - polyx[b][MidPt]);
ElY = (El3 - Elev[MidPt]) / (polyy[b][MidPt] - Y4);
} /* if */
else
{
ElInterp = 0;
ElY = ElX = 0.0;
ElStart = el;
goto StartDraw;
} /* else */
} /* else not so easy */
ElStart = Elev[MidPt] + (polyx[b][PolyMinXPt] - polyx[b][MidPt]) * ElX;
ElStart += (polyy[b][0] - polyy[b][MidPt]) * ElY;
StartDraw:
if (QY != 0.0)
EloQY = .001 * ElY / QY;
else
EloQY = 0.0;
for (y=yy[0], k=0; y<=yy[2]; y++,k++,QStart+=QY,ElStart+=ElY,LatStart+=LatY,LonStart+=LonY)
{
if (y < 0) continue;
if (y > oshigh) break;
zip = scrnrowzip[y] + MinX;
for (x=MinX, l=0, QPt=QStart, ElPt=ElStart, LatPt=LatStart, LonPt=LonStart;
x<=MaxX; x++, l++, zip++, QPt+=QX, ElPt+=ElX, LatPt+=LatX, LonPt+=LonX)
{
if (x < 0)
continue;
if (x > wide)
break;
dist = *(zbuf + zip);
if (! QInterp)
QPt = qqq;
else
{
if (QPt < polyq[b][MinQPt])
QPt = polyq[b][MinQPt];
else if (QPt > polyq[b][MaxQPt])
QPt = polyq[b][MaxQPt];
} /* else */
if (! ElInterp)
ElPt = el;
else
{
if (ElPt < Elev[MinPt])
ElPt = Elev[MinPt];
else if (ElPt > Elev[MaxPt])
ElPt = Elev[MaxPt];
} /* else */
if (QPt <= dist + PtrnOffset || (unsigned int)bytemap[zip] < 100)
{
sum = 0;
PixPtr = &SubPix[k * 10 * SubPixWidth + l * 10];
for (i=0; i<10; i++)
{
for (j=0; j<10; j++)
{
sum += (unsigned int)PixPtr[j];
} /* for j=0... */
PixPtr += SubPixWidth;
} /* for i=0... */
if (sum)
{
if (render & 0x01)
{
/* changed this from ptrnoffset to HalfPtrnOffset 1/2/96 to elliminate
spider webs */
if (bytemap[zip]
&& (QPt >= dist || ((sum < 100) && (dist - QPt < HalfPtrnOffset))))
Mode = MODE_AVERAGE;
else
Mode = MODE_REPLACE;
if (VolumeTexture)
{
if (Texture & 0x0200) /* if strata colors */
{
Tex = ComputeTextureColor(ElPt, LatPt, LonPt, ElY, &CC[0]);
CC[0].Red -= sunshade * CC[0].Red; /* shading */
CC[0].Grn -= sunshade * CC[0].Grn;
CC[0].Blu -= sunshade * CC[0].Blu;
} /* if */
else if (Texture & 0x0100) /* if plain strata */
{
Tex = ComputeTexture(ElPt, LatPt, LonPt, ElY);
} /* else if */
else Tex = 255;
/* This is a test - replace with optional fracture texture
Tex *= ((355 - MakeNoise(NoiseMap, 200, LatPt * 480.0, LonPt * 480.0)) / 255.0);
if (Tex > 255)
Tex = 255;
*/
} /* if strata */
/*
else if (NoiseTexture)
{
Tex = ((255 - MakeNoise(NoiseMap, 75, LatPt * 480.0, LonPt * 480.0)) / 255.0);
/* Tex = ComputeBumpMapTexture(LatPt, LonPt);*/
} /* else no strata */
*/
else
Tex = 255;
fred = (Tex * CC[0].Red) / 255.;
fgrn = (Tex * CC[0].Grn) / 255.;
fblu = (Tex * CC[0].Blu) / 255.;
fred += CC[3].Red; /* sunshade * PARC_RNDR_COLOR(1, 0); ambient */
fgrn += CC[3].Grn; /* sunshade * PARC_RNDR_COLOR(1, 1); */
fblu += CC[3].Blu; /* sunshade * PARC_RNDR_COLOR(1, 2); */
fred += (PARC_RNDR_COLOR(2, 0) - fred) * fade;
fgrn += (PARC_RNDR_COLOR(2, 1) - fgrn) * fade;
fblu += (PARC_RNDR_COLOR(2, 2) - fblu) * fade;
fred += (PARC_RNDR_COLOR(2, 0) - fred) * fog;
fgrn += (PARC_RNDR_COLOR(2, 1) - fgrn) * fog;
fblu += (PARC_RNDR_COLOR(2, 2) - fblu) * fog;
fred *= redsun;
fgrn *= greensun;
fblu *= bluesun;
if (fred > 255.) fred = 255.;
else if (fred < 0.) fred = 0.;
if (fgrn > 255.) fgrn = 255.;
else if (fgrn < 0.) fgrn = 0.;
if (fblu > 255.) fblu = 255.;
else if (fblu < 0.) fblu = 0.;
if (Mode)
{
if ((unsigned int)bytemap[zip] >= 110)
continue; /* added 9/25/95 */
WtAvg = (unsigned int)bytemap[zip] + sum >= 110 ? 110:
(unsigned int)bytemap[zip] + sum; /* changed 9/25/95 */
sum = WtAvg - bytemap[zip]; /* added 9/25/95 */
colval = (*(bitmap[0] + zip) * (unsigned int)bytemap[zip] + fred * sum)
/ WtAvg;
*(bitmap[0] + zip) = (UBYTE)colval;
colval = (*(bitmap[1] + zip) * (unsigned int)bytemap[zip] + fgrn * sum)
/ WtAvg;
*(bitmap[1] + zip) = (UBYTE)colval;
colval = (*(bitmap[2] + zip) * (unsigned int)bytemap[zip] + fblu * sum)
/ WtAvg;
*(bitmap[2] + zip) = (UBYTE)colval;
} /* if already a value */
else
{
*(bitmap[0] + zip) = (UBYTE)fred;
*(bitmap[1] + zip) = (UBYTE)fgrn;
*(bitmap[2] + zip) = (UBYTE)fblu;
} /* else first value */
} /* if render to bitmaps */
if (QPt < dist)
{
if (render & 0x10)
{
if (render & 0x01)
{
ScreenPixelPlot(win, bitmap, x + drawoffsetX, y + drawoffsetY, zip);
} /* if render to bitmaps */
else
{
NoRGBScreenPixelPlot(win, FloatCol, ColMax, x + drawoffsetX, y + drawoffsetY);
} /* else no bitmaps */
} /* if render to screen */
if (render & 0x100)
{
*(QCmap[0] + zip) = QC.compval1;
*(QCmap[1] + zip) = QC.compval2;
*(QCmap[2] + zip) = QC.compval3;
*(QCcoords[0] + zip) = facelat;
*(QCcoords[1] + zip) = facelong;
} /* if render & 0x100 */
*(zbuf + zip) = QPt;
if (Reflections && ReflectionMap)
ReflectionMap[zip] = Reflections;
if (ElevationMap)
ElevationMap[zip] = ElPt;
if (SlopeMap)
SlopeMap[zip] = EloQY;
} /* if lower QPt value */
bytemap[zip] += sum;
} /* if at least one point */
} /* if QPt */
} /* for x=... */
} /* for y=... */
} /* if array is large enough to hold current polygon */
} /* if polygon is visible */
/* Now we draw the understory layer if undertreedraw > 0 */
DrawTrees:
if (undertreedraw)
{
DetailTree(win, understory, &QC, 1, &CC[1], &CC[3], Elev);
} /* if understory tree draw */
/* Now we draw the overstory layer if treedraw > 0.
** treedraw could be 2 if forest model
*/
if (treedraw)
{
DetailTree(win, eco, &QC, (short)(treedraw == 1), &CC[2], &CC[3], Elev);
} /* if treedraw */
} /* MapTopo() */
/***********************************************************************/
short colormap(struct elmapheaderV101 *map, short notsnow,
struct ColorComponents *CC, short *understory)
{
UBYTE matchred, matchgreen, matchblue;
short colpts = 0, cpts = 0, eco;
long Row, Col, facept[3];
#ifdef ENABLE_STATISTICS
long ct, elct;
#endif /* ENABLE_STATISTICS */
double ecoline;
/* determine number of points falling on color map */
if (settings.mastercmap)
{
if (facelat >= CMap->South && facelat <= CMap->North &&
facelong <= CMap->West && facelong >= CMap->East)
{
if (settings.cmaporientation)
{
Row = (CMap->North - facelat) / CMap->StepLat;
Col = (CMap->West - facelong) / CMap->StepLon;
facept[0] = facept[1] = facept[2] = Row * CMap->Cols + Col;
} /* if image orientation */
else
{
Col = (facelat - CMap->South) / CMap->StepLat;
Row = (CMap->West - facelong) / CMap->StepLon;
facept[0] = facept[1] = facept[2] = Row * CMap->Cols + Col;
} /* else WCS DEM orientation */
if (Row < CMap->Rows && Col < CMap->Cols)
{
if ( colmap[0][facept[0]] || colmap[1][facept[0]] || colmap[2][facept[0]])
{
colpts = 7;
cpts = 3;
} /* if */
else
return (-1);
} /* if */
else
{
return (-1);
}
} /* if in bounds */
else
return (-1);
} /* if mastercmap */
else if (settings.borderandom)
{
randomint = treerand * 2.99 + ROUNDING_KLUDGE;
facept[0] = map->facept[randomint];
facept[1] = map->facept[randomint];
facept[2] = map->facept[randomint];
if ( colmap[0][facept[0]] || colmap[1][facept[0]] || colmap[2][facept[0]] )
{
colpts = 7;
cpts = 3;
} /* if */
} /* if settings.borderandom */
else
{
facept[0] = map->facept[0];
facept[1] = map->facept[1];
facept[2] = map->facept[2];
if ( colmap[0][facept[0]] || colmap[1][facept[0]] || colmap[2][facept[0]] )
{
cpts = colpts = 1;
} /* if */
if ( colmap[0][facept[1]] || colmap[1][facept[1]]
|| colmap[2][facept[1]] )
{
colpts += 2;
cpts += 1;
} /* if */
if ( colmap[0][facept[2]] || colmap[1][facept[2]]
|| colmap[2][facept[2]] )
{
colpts += 4;
cpts += 1;
} /* if */
} /* else not random borders */
/* if at least one color map point */
if (cpts)
{
#ifdef ENABLE_STATISTICS
if (settings.statistics && cpts > 1)
{
ct = aspect * 5.72957795 + 1.0;
elct = (el - settings.surfel[3]) * 20
/ (settings.surfel[3] - settings.surfel[0]);
if (elct < 0) elct = 0;
else if (elct > 19) elct = 19;
statcount[ct][elct] ++;
meanel += el;
meanaspect += aspect;
meanslope +=slope;
samples ++;
} /* if settings.statistics */
#endif /* ENABLE_STATISTICS */
/* if match colors to ecosystems */
if (settings.ecomatch)
{
switch (colpts)
{
case 1:
{
matchred = colmap[0][facept[0]];
matchgreen = colmap[1][facept[0]];
matchblue = colmap[2][facept[0]];
break;
} /* case 1 */
case 2:
{
matchred = (UBYTE)colmap[0][facept[1]];
matchgreen = (UBYTE)colmap[1][facept[1]];
matchblue = (UBYTE)colmap[2][facept[1]];
break;
} /* case 2 */
case 3:
{
if (treerand > .5)
{
matchred = (UBYTE)colmap[0][facept[0]];
matchgreen = (UBYTE)colmap[1][facept[0]];
matchblue = (UBYTE)colmap[2][facept[0]];
} /* if */
else
{
matchred = (UBYTE)colmap[0][facept[1]];
matchgreen = (UBYTE)colmap[1][facept[1]];
matchblue = (UBYTE)colmap[2][facept[1]];
} /* else */
break;
} /* case 3 */
case 4:
{
matchred = (UBYTE)colmap[0][facept[2]];
matchgreen = (UBYTE)colmap[1][facept[2]];
matchblue = (UBYTE)colmap[2][facept[2]];
break;
} /* case 4 */
case 5:
{
if (treerand > .5)
{
matchred = (UBYTE)colmap[0][facept[0]];
matchgreen = (UBYTE)colmap[1][facept[0]];
matchblue = (UBYTE)colmap[2][facept[0]];
} /* if */
else
{
matchred = (UBYTE)colmap[0][facept[2]];
matchgreen = (UBYTE)colmap[1][facept[2]];
matchblue = (UBYTE)colmap[2][facept[2]];
} /* else */
break;
} /* case 5 */
case 6:
{
if (treerand > .5)
{
matchred = (UBYTE)colmap[0][facept[1]];
matchgreen = (UBYTE)colmap[1][facept[1]];
matchblue = (UBYTE)colmap[2][facept[1]];
} /* if */
else
{
matchred = (UBYTE)colmap[0][facept[2]];
matchgreen = (UBYTE)colmap[1][facept[2]];
matchblue = (UBYTE)colmap[2][facept[2]];
} /* else */
break;
} /* case 6 */
case 7:
{
if (treerand > .67)
{
matchred = colmap[0][facept[0]];
matchgreen = colmap[1][facept[0]];
matchblue = colmap[2][facept[0]];
} /* if */
else if (treerand > .33)
{
matchred = colmap[0][facept[1]];
matchgreen = colmap[1][facept[1]];
matchblue = colmap[2][facept[1]];
} /* else if */
else
{
matchred = colmap[0][facept[2]];
matchgreen = colmap[1][facept[2]];
matchblue = colmap[2][facept[2]];
} /* else */
break;
} /* case 7 */
} /* switch colpts */
/* find ecosystem match if there is one and slope is appropriate */
for (eco=0; eco<ECOPARAMS; eco++)
{
if (matchred == PAR_MTCH_ECO(eco, 0))
{
if (matchgreen == PAR_MTCH_ECO(eco, 1))
{
if (matchblue == PAR_MTCH_ECO(eco, 2))
{
if (cpts > 1 && slope <= PARC_MXSL_ECO(eco))
{
if (eco != 1) *understory = ecoset(eco, notsnow, CC);
return (eco);
} /* if cpts > 1 */
else
{
return (-1);
} /* else cpts=1 */
} /* match blue */
} /* match green */
} /* match red */
} /* for eco=0... */
} /* if settings.ecomatch */
/* normal color mapping */
/* if all three points color mapped set color to average */
if (cpts == 3)
{
CC[0].Red = ((UBYTE)colmap[0][facept[0]] + (UBYTE)colmap[0][facept[1]] +
(UBYTE)colmap[0][facept[2]]) / 3.0;
CC[0].Grn = ((UBYTE)colmap[1][facept[0]] + (UBYTE)colmap[1][facept[1]] +
(UBYTE)colmap[1][facept[2]]) / 3.0;
CC[0].Blu = ((UBYTE)colmap[2][facept[0]] + (UBYTE)colmap[2][facept[1]] +
(UBYTE)colmap[2][facept[2]]) / 3.0;
if (! settings.cmapluminous)
{
CC[0].Red -= sunshade * CC[0].Red; /* shading */
CC[0].Red += sunshade * PARC_RNDR_COLOR(1, 0); /* ambient */
CC[0].Grn -= sunshade * CC[0].Grn;
CC[0].Grn += sunshade * PARC_RNDR_COLOR(1, 1);
CC[0].Blu -= sunshade * CC[0].Blu;
CC[0].Blu += sunshade * PARC_RNDR_COLOR(1, 2);
} /* if not luminous color mapping */
eco = settings.defaulteco;
} /* if cpts = 3 */
/* otherwise find normal ecosystem */
if (cpts < 3 || settings.cmaptrees)
{
if (notsnow)
{
for (eco=12; eco<ECOPARAMS; eco++)
{
ecoline = PARC_RNDRLN_ECO(eco) +
diplong * PARC_SKLN_ECO(eco) + diplat * PARC_SKLT_ECO(eco)
+ PARC_RNDRRE_ECO(eco) * relfactor;
if (settings.worldmap)
ecoline -= ((abs(facelat) - abs(settings.globreflat)) * settings.globecograd);
if (settings.flatteneco)
ecoline += ( (PARC_RNDR_MOTION(13) - ecoline) * PARC_RNDR_MOTION(12) );
if (el <= ecoline)
{
if (slope >= PARC_MNSL_ECO(eco) && slope <= PARC_MXSL_ECO(eco))
{
if (relfactor >= PARC_RNDRNR_ECO(eco) &&
relfactor <= PARC_RNDRXR_ECO(eco))
{
break;
} /* if relel matchup */
} /* if slope matchup */
} /* if elevation matchup */
} /* for eco=12... */
if (eco >= ECOPARAMS)
{
Log(WNG_ILL_VAL, "Ecosystem out of range.");
eco = settings.defaulteco;
} /* if */
CC[0].Red = PARC_SCOL_ECO(eco, 0);
CC[0].Grn = PARC_SCOL_ECO(eco, 1);
CC[0].Blu = PARC_SCOL_ECO(eco, 2);
if (! settings.cmapluminous)
{
CC[0].Red -= sunshade * CC[0].Red; /* shading */
CC[0].Red += sunshade * PARC_RNDR_COLOR(1, 0); /* ambient */
CC[0].Grn -= sunshade * CC[0].Grn;
CC[0].Grn += sunshade * PARC_RNDR_COLOR(1, 1);
CC[0].Blu -= sunshade * CC[0].Blu;
CC[0].Blu += sunshade * PARC_RNDR_COLOR(1, 2);
} /* if not luminous */
} /* if not snow */
colmapavg(map, colpts, &CC[0]);
} /* if cpts < 3 || cmaptrees */
else
eco = settings.defaulteco;
/* set colors for rendering */
if (settings.cmaptrees && PARC_RNDRDN_ECO(eco) > treerand2 * 100)
treedraw = 1;
CC[2].Red = CC[1].Red = CC[0].Red;
CC[2].Grn = CC[1].Grn = CC[0].Grn;
CC[2].Blu = CC[1].Blu = CC[0].Blu;
return (eco);
} /* if cpts */
else
{
return (-1);
} /* else ! cpts */
} /* colormap() */
/********************************************************************/
#define SCAN_NORMAL 1
#define SCAN_REVERSE 0
void MapCloud(struct Window *win, short *CloudVal, short *IllumVal, long Elev)
{
short x, y, ScanOrder = SCAN_NORMAL, Edge2aHt, Edge2bHt, Edge1Ht;
long scrnrow, zip;
double m, EdgeVal, ValCX, ValCY, dCX, dCY, ValIX, ValIY, dIX, dIY;
if (! CloudVal)
return;
Edge2aHt = yy[1] - yy[0];
Edge2bHt = yy[2] - yy[1];
Edge1Ht = yy[2] - yy[0];
if (EdgeSize < Edge1Ht * sizeof (short))
{
free_Memory(Edge1, EdgeSize);
free_Memory(Edge2, EdgeSize);
EdgeSize = Edge1Ht * sizeof (short);
Edge1 = (short *)get_Memory(EdgeSize, MEMF_CLEAR);
Edge2 = (short *)get_Memory(EdgeSize, MEMF_CLEAR);
if (! Edge1 || ! Edge2)
{
return;
}
} /* if need new edge arrays */
if (Edge1Ht)
{
m = ((float)xx[2] - xx[0]) / (Edge1Ht);
dCY = ((float)CloudVal[2] - CloudVal[0]) / (float)Edge1Ht;
dIY = ((float)IllumVal[2] - IllumVal[0]) / (float)Edge1Ht;
} /* if */
else
{
m = 0.0;
dCY = 0.0;
dIY = 0.0;
}
EdgeVal = xx[0];
for (y=0; y<=Edge1Ht; y++, EdgeVal += m)
Edge1[y] = EdgeVal;
if (Edge2aHt)
m = ((float)xx[1] - xx[0]) / (Edge2aHt);
else
m = 0.0;
EdgeVal = xx[0];
for (y=0; y<=Edge2aHt; y++, EdgeVal += m)
Edge2[y] = EdgeVal;
if (Edge2bHt)
m = ((float)xx[2] - xx[1]) / (Edge2bHt);
else
m = 0.0;
EdgeVal = xx[1];
for (y=Edge2aHt; y<=Edge1Ht; y++, EdgeVal += m)
Edge2[y] = EdgeVal;
if (Edge1[Edge2aHt] > Edge2[Edge2aHt])
ScanOrder = SCAN_REVERSE;
if (Edge2[Edge2aHt] != Edge1[Edge2aHt])
{
ValCY = CloudVal[0] + dCY * Edge2aHt;
ValIY = IllumVal[0] + dIY * Edge2aHt;
dCX = ((float)CloudVal[1] - ValCY) / ((float)Edge2[Edge2aHt] - Edge1[Edge2aHt]);
dIX = ((float)IllumVal[1] - ValIY) / ((float)Edge2[Edge2aHt] - Edge1[Edge2aHt]);
} /* if */
else
{
dCX = 0.0;
dIX = 0.0;
}
ValCY = CloudVal[0];
ValIY = IllumVal[0];
if (render & 0x010)
SetAPen(win->RPort, AltPen[1]);
if (ScanOrder)
{
for (y=0, scrnrow=yy[0]; y<=Edge1Ht; y++, scrnrow++, ValCY += dCY, ValIY += dIY)
{
if (scrnrow < 0)
continue;
if (scrnrow > high)
break;
zip = scrnrowzip[scrnrow] + Edge1[y];
ValCX = ValCY;
ValIX = ValIY;
for (x=Edge1[y]; x<=Edge2[y]; x++, zip++, ValCX += dCX, ValIX += dIX)
{
if (x < 0)
continue;
if (x > wide)
break;
if (qqq < *(zbuf + zip))
{
if (ValCX < 0.0)
ValCX = 0.0;
else if (ValCX > 255.0)
ValCX = 255.0;
if (ValIX < 0.0)
ValIX = 0.0;
else if (ValIX > 255.0)
ValIX = 255.0;
*(bytemap + zip) = (USHORT)ValCX;
*(bytemap + zip) += (((USHORT)ValIX) << 8);
*(zbuf + zip) = qqq;
if (ElevationMap)
ElevationMap[zip] = Elev;
} /* if */
} /* for x=... */
} /* for y=0... */
} /* if normal scan - left to right */
else
{
for (y=0, scrnrow=yy[0]; y<=Edge1Ht; y++, scrnrow++, ValCY += dCY, ValIY += dIY)
{
if (scrnrow < 0)
continue;
if (scrnrow > high)
break;
zip = scrnrowzip[scrnrow] + Edge1[y];
ValCX = ValCY;
ValIX = ValIY;
for (x=Edge1[y]; x>=Edge2[y]; x--, zip--, ValCX -= dCX, ValIX -= dIX)
{
if (x < 0)
continue;
if (x > wide)
break;
if (qqq < *(zbuf + zip))
{
if (ValCX < 0.0)
ValCX = 0.0;
else if (ValCX > 255.0)
ValCX = 255.0;
if (ValIX < 0.0)
ValIX = 0.0;
else if (ValIX > 255.0)
ValIX = 255.0;
*(bytemap + zip) = (USHORT)ValCX;
*(bytemap + zip) += (((USHORT)ValIX) << 8);
*(zbuf + zip) = qqq;
if (ElevationMap)
ElevationMap[zip] = Elev;
} /* if */
} /* for x=... */
} /* for y=0... */
} /* else scan reverse */
} /* MapCloud() */
/***********************************************************************/
double CloudCover_Set(struct CloudData *CD, double Lat, double Lon)
{
return(DEM_InterpPt(&CD->Map, Lat, Lon) / 510.0);
/*
long Row, Col;
double CloudCover = 0.0, LonOff, LatOff, LonInvOff, LatInvOff, wt[4], val[4], Area;
if (Lat >= CD->Map.lolat &&
Lat <= CD->Map.lolat + (CD->Map.columns - 1) * CD->Map.steplat &&
Lon <= CD->Map.lolong &&
Lon >= CD->Map.lolong - CD->Map.rows * CD->Map.steplong)
{
LatOff = Lat - CD->Map.lolat;
LonOff = CD->Map.lolong - Lon;
Col = LatOff / CD->Map.steplat;
Row = LonOff / CD->Map.steplong;
LatOff -= (Col * CD->Map.steplat);
LonOff -= (Row * CD->Map.steplong);
LonInvOff = CD->Map.steplong - LonOff;
LatInvOff = CD->Map.steplat - LatOff;
if (Col >= CD->Map.columns - 1)
{
if (Row >= CD->Map.rows)
{
CloudCover = CD->Map.map[Row * CD->Map.columns + Col] / 510.0;
}
else
{
wt[0] = LonInvOff / CD->Map.steplong;
val[0] = CD->Map.map[Row * CD->Map.columns + Col];
wt[3] = 1.0 - wt[0];
val[3] = CD->Map.map[(Row + 1) * CD->Map.columns + Col];
CloudCover = (wt[0] * val[0] + wt[3] * val[3]) / 510.0;
}
}
else if (Row >= CD->Map.rows)
{
wt[0] = LatInvOff / CD->Map.steplat;
val[0] = CD->Map.map[Row * CD->Map.columns + Col];
wt[1] = 1.0 - wt[0];
val[1] = CD->Map.map[Row * CD->Map.columns + Col + 1];
CloudCover = (wt[0] * val[0] + wt[1] * val[1]) / 510.0;
}
else
{
Area = CD->Map.steplat * CD->Map.steplong;
wt[0] = LatInvOff * LonInvOff / Area;
val[0] = CD->Map.map[Row * CD->Map.columns + Col];
wt[1] = LatOff * LonInvOff / Area;
val[1] = CD->Map.map[Row * CD->Map.columns + Col + 1];
wt[2] = LatOff * LonOff / Area;
val[2] = CD->Map.map[(Row + 1) * CD->Map.columns + Col + 1];
wt[3] = LonOff * LatInvOff / Area;
val[3] = CD->Map.map[(Row + 1) * CD->Map.columns + Col];
CloudCover = (wt[0] * val[0] + wt[1] * val[1] + wt[2] * val[2] + wt[3] * val[3]) / 510.0;
}
/* max value of cloud map is 255 so dividing by 510 allows cloudcover
to reach 50% */
} /* if in bounds */
return (CloudCover);
*/
} /* CloudCover_Set() */
