SurfaceTracerImpl.mesa
Copyright © 1986 by Xerox Corporation. All rights reserved.
James Rauen, August 25, 1986 11:15:51 pm PDT
Last edited by: James Rauen January 13, 1988 4:20:03 pm PST
Running SurfaceTracerImpl adds the command "MakePic" to the Command Tool.
DIRECTORY
AIS USING [CloseFile, CreateFile, FRef, OpenWindow, Raster, RasterPart, WRef, WriteSample],
CADTypes USING [Scad, ScadSequence, VariableRec, VisibleMask, VisibleMaskSequence],
CADIO USING [Error, WhatCellIsItOn],
Commander USING [CommandProc, Register],
ImagerColor USING [RGB],
IO USING [int, PutFR],
LightingModels USING [AmbientSource, Model, PointSource, PointSourceSequence],
MessageWindow USING [Append],
MultiPolynomial USING [Differentiate, EvaluationBindings, MultiPolSequence, Ref, RefFromRope, TotallyEvaluate, TotallySubstituteUnivariates, TSUBindings],
Polynomial USING [Linear, RealRootRec, RealRoots, Ref],
Real USING [Float, InlineFixC],
Rope USING [Concat, ROPE],
ShadingModels USING [Model, ShadingSequence],
SurfaceTracer USING [],
Geometry3dVector USING [Add, Combine, Cross, Dot, Mul, Normalize, Sub, Triple];

SurfaceTracerImpl: CEDAR PROGRAM
IMPORTS AIS, CADIO, Commander, IO, MessageWindow, MultiPolynomial, Polynomial, Real, Rope, Geometry3dVector
EXPORTS SurfaceTracer
~ BEGIN
Type declarations
ALPO: TYPE ~ MultiPolynomial.Ref;
RootSequence: TYPE ~ RECORD [
roots: SEQUENCE length: NAT OF REF Polynomial.RealRootRec];
Triple: TYPE ~ Geometry3dVector.Triple;
r: REAL ← 1.0;
RGB: TYPE ~ ImagerColor.RGB;
SurfaceSequence: TYPE ~ SurfaceTracer.SurfaceSequence;
Global constants
infinity: REAL ← 1E30;
epsilon: REAL ← 1E-4;
Public procedures
TraceSurfaces: PUBLIC PROC [
surfaces: REF MultiPolynomial.MultiPolSequence,
variables: CADTypes.VariableRec,
colors: REF ShadingModels.ShadingSequence,
screenCenter, screenU, screenV: Triple,
pixelsU, pixelsV: NAT,
lightingModel: LightingModels.Model,
filename: Rope.ROPE] ~ BEGIN
Construct the ST surface sequence representation.
surfaceSeq: REF STSurfaceSequence ← NEW[STSurfaceSequence[surfaces.length]];
FOR i: NAT IN [0..surfaces.length) DO
substitute for $x, $y, $z
surfaceSeq[i] ← [
scad: [
surface: surfaces[i],
named: FALSE, --garbage
name: "",  --garbage
color: [1, 1, 1], --garbage
cells: NIL  --garbage
],
gradientX: MultiPolynomial.Differentiate[surfaces[i], $x],
gradientY: MultiPolynomial.Differentiate[surfaces[i], $y],
gradientZ: MultiPolynomial.Differentiate[surfaces[i], $z],
mask: NIL,
shading: colors[i]
];
ENDLOOP;
Draw the surface.
DrawSurface [
surfaceSeq: surfaceSeq,
variables: variables,
screenCenter: screenCenter,
screenU: screenU,
screenV: screenV,
pixelsU: pixelsU,
pixelsV: pixelsV,
lightingModel: lightingModel,
filename: filename,
n: 1,
doCells: FALSE];
END;
TraceCells: PUBLIC PROC [
surfaces: REF CADTypes.ScadSequence,
variables: CADTypes.VariableRec,
colors: REF ShadingModels.ShadingSequence,
masks: REF CADTypes.VisibleMaskSequence,
screenCenter, screenU, screenV: Triple,
pixelsU, pixelsV: NAT,
lightingModel: LightingModels.Model,
filename: Rope.ROPE] ~ BEGIN
Construct the ST surface sequence representation.
surfaceSeq: REF STSurfaceSequence ← NEW[STSurfaceSequence[surfaces.length]];
FOR i: NAT IN [0..surfaces.length) DO
substitute for $x, $y, $z
surfaceSeq[i] ← [
scad: surfaces[i],
gradientX: MultiPolynomial.Differentiate[surfaces[i].surface, variables.x],
gradientY: MultiPolynomial.Differentiate[surfaces[i].surface, variables.y],
gradientZ: MultiPolynomial.Differentiate[surfaces[i].surface, variables.z],
mask: masks[i],
shading: colors[i]
];
ENDLOOP;
Draw the surface.
DrawSurface [
surfaceSeq: surfaceSeq,
variables: variables,
screenCenter: screenCenter,
screenU: screenU,
screenV: screenV,
pixelsU: pixelsU,
pixelsV: pixelsV,
lightingModel: lightingModel,
filename: filename,
n: 1,
doCells: TRUE];
END;
Private surface representation
Polynomials are in (x, y, z) coordinates. Other information kept besides the Scad are the gradients, the shading of the surface, and the cell mask.
STSurface: TYPE ~ RECORD [
scad: CADTypes.Scad,
gradientX: MultiPolynomial.Ref,
gradientY: MultiPolynomial.Ref,
gradientZ: MultiPolynomial.Ref,
mask: REF CADTypes.VisibleMask,
shading: ShadingModels.Model];
STSurfaceSequence: TYPE ~ RECORD [
surfaces: SEQUENCE length: NAT OF STSurface];
Image generation
DrawSurface: PROC [surfaceSeq: REF STSurfaceSequence, variables: CADTypes.VariableRec, screenCenter, screenU, screenV: Triple, pixelsU, pixelsV: CARDINAL, lightingModel: LightingModels.Model, n: NAT ← 0, filename: Rope.ROPE ← "Surface", doCells: BOOLEAN] ~ BEGIN
Create the AIS mechanism
raster: AIS.Raster ← NEW[AIS.RasterPart ← [pixelsV, pixelsU, rd, 8, -1, 65535]];
redFile: AIS.FRef ← AIS.CreateFile[Rope.Concat[filename, "-red.ais"], raster];
greenFile: AIS.FRef ← AIS.CreateFile[Rope.Concat[filename, "-green.ais"], raster];
blueFile: AIS.FRef ← AIS.CreateFile[Rope.Concat[filename, "-blue.ais"], raster];
redWindow: AIS.WRef ← AIS.OpenWindow[redFile, 0, pixelsV-1, 0, pixelsU-1];
greenWindow: AIS.WRef ← AIS.OpenWindow[greenFile, 0, pixelsV-1, 0, pixelsU-1];
blueWindow: AIS.WRef ← AIS.OpenWindow[blueFile, 0, pixelsV-1, 0, pixelsU-1];
break: SIGNAL = CODE;
Create Andrew's Adaptive Anti-Aliasing Algorithm mechanism.
GotItSeq: TYPE ~ RECORD [yesOrNo: SEQUENCE length: NAT OF BOOLEAN];
gotIt: REF GotItSeq ← NEW[GotItSeq[9]];
ColorsTracedSeq: TYPE ~ RECORD [colorsTraced: SEQUENCE length: NAT OF RGB];
colorsTraced: REF ColorsTracedSeq ← NEW[ColorsTracedSeq[9]];
TracePixel: PROC[ul, ur, ll, lr: NAT] RETURNS[RGB] ~ BEGIN
ulColor, urColor, llColor, lrColor: RGB;
IF gotIt[ul] THEN ulColor ← colorsTraced[ul]
ELSE BEGIN
ulColor ← colorsTraced[ul] ← Trace[];
gotIt[ul] ← TRUE;
END;
IF gotIt[ur] THEN urColor ← colorsTraced[ur]
ELSE BEGIN
urColor ← colorsTraced[ur] ← Trace[];
gotIt[ur] ← TRUE;
END;
IF gotIt[ll] THEN llColor ← colorsTraced[ll]
ELSE BEGIN
llColor ← colorsTraced[ll] ← Trace[];
gotIt[ll] ← TRUE;
END;
IF gotIt[lr] THEN lrColor ← colorsTraced[lr]
ELSE BEGIN
lrColor ← colorsTraced[lr] ← Trace[];
gotIt[lr] ← TRUE;
END;
IF ur-ul = 1 OR CloseEnough[ulColor, urColor, llColor, lrColor] THEN
RETURN [AverageRGB[ulColor, urColor, llColor, lrColor]]
ELSE BEGIN
um, lm, ml, mr, mm: NAT;
v1, v2, v3, v4: RGB;
um ← (ul + ur)/2;
lm ← (ll + lr)/2;
ml ← (ul + ll)/2;
mr ← (ur + lr)/2;
mm ← (ul + lr)/2;
v1 ← TracePixel [ul, um, ml, mm];
v2 ← TracePixel [um, ur, mm, mr];
v3 ← TracePixel [ml, mm, ll, lm];
v4 ← TracePixel [mm, mr, lm, lr];
RETURN [AverageRGB[v1, v2, v3, v4]];
END;
END;
CloseEnough: PROC[c1, c2, c3, c4: RGB] RETURNS[BOOLEAN] ~ BEGIN
RETURN[TRUE];
END;
AverageRGB: PROC[c1, c2, c3, c4: RGB] RETURNS[RGB] ~ BEGIN
RETURN[[
(c1.R+c2.R+c3.R+c4.R)/4.0,
(c1.G+c2.G+c3.G+c4.G)/4.0,
(c1.B+c2.B+c3.B+c4.B)/4.0]];
END;
UVOffsets: PROC[subPixelIndex: NAT] RETURNS[uOffset, vOffset: REAL] ~ BEGIN
subRow: NAT ←
Stuff
RoughSampleGrid: TYPE ~ RECORD[
rows: SEQUENCE nRows: NAT OF REF RoughSampleRow];
RoughSampleRow: TYPE ~ RECORD[
columns: SEQUENCE nCols: NAT OF RGB];
roughSampleGrid: REF RoughSampleGrid;
Calculate the viewer direction vector.
direction: Triple ← Geometry3dVector.Normalize[Geometry3dVector.Cross[screenV, screenU]];
Generate the first grid of samples.
roughSampleGrid ← NEW[RoughSampleGrid[pixelsV + 1]];
FOR row: NAT IN [0..pixelsV] DO
roughSampleGrid[row] ← NEW[RoughSampleRow[pixelsU + 1]];
ENDLOOP;
FOR row: NAT IN [0..pixelsV] DO
FOR column: NAT IN [0..pixelsU] DO
Calculate the (u, v) of the upper left corner of the pixel at (row, column). The range of u and v is -1 <= u <= 1, -1 <= v <= 1.
u: REAL ← ((2 * Real.Float[column]) / Real.Float[pixelsU]) - 1;
v: REAL ← ((2 * Real.Float[row]) / Real.Float[pixelsV]) - 1;
Calculate the (x, y, z) of the upper left corner of the pixel.
screenOffset: Triple ← Geometry3dVector.Combine[screenU, u, screenV, v];
startingPoint: Triple ← Geometry3dVector.Add[screenCenter, screenOffset];
Trace the upper left corner of the pixel.
color: RGB ← Trace[surfaceSeq, variables, startingPoint, direction, lightingModel, doCells];
roughSampleGrid[column][row] ← color;
Do the next sample. Continue for one extra row and one extra column of pixels, so that tracing all the upper left corners will generate a complete set of samples.
ENDLOOP;
MessageWindow.Append [IO.PutFR ["Sampled %g of %g rows", IO.int[row + 1], IO.int[pixelsV + 1]], TRUE];
ENDLOOP;
Iterate over the pixels and calculate each.
FOR row: CARDINAL IN [0..pixelsV) DO
FOR column: CARDINAL IN [0..pixelsU) DO
red: REAL ← (
roughSampleGrid[row][column].R +
roughSampleGrid[row + 1][column].R +
roughSampleGrid[row][column + 1].R +
roughSampleGrid[row + 1][column + 1].R) / 4;
green: REAL ← (
roughSampleGrid[row][column].G +
roughSampleGrid[row + 1][column].G +
roughSampleGrid[row][column + 1].G +
roughSampleGrid[row + 1][column + 1].G) / 4;
blue: REAL ← (
roughSampleGrid[row][column].B +
roughSampleGrid[row + 1][column].B +
roughSampleGrid[row][column + 1].B +
roughSampleGrid[row + 1][column + 1].B) / 4;
redIntensityCard: CARDINAL ← Real.InlineFixC[red*255.0];
greenIntensityCard: CARDINAL ← Real.InlineFixC[green*255.0];
blueIntensityCard: CARDINAL ← Real.InlineFixC[blue*255.0];
AIS.WriteSample [redWindow, redIntensityCard, row, column];
AIS.WriteSample [greenWindow, greenIntensityCard, row, column];
AIS.WriteSample [blueWindow, blueIntensityCard, row, column];
ENDLOOP;
ENDLOOP;
Terminate gracefully.
AIS.CloseFile[redFile];
AIS.CloseFile[greenFile];
AIS.CloseFile[blueFile];
END;
Trace: PROC [surfaces: REF STSurfaceSequence, variables: CADTypes.VariableRec, fromPoint, direction: Triple, lighting: LightingModels.Model, doCells: BOOLEAN] RETURNS [color: RGB] ~ BEGIN
Declarations
t, xGrad, yGrad, zGrad: REAL;
intersection, gradient, normal, surfaceToEye: Triple;
somethingToDisplay: BOOLEAN;
break: SIGNAL = CODE;
gradientBindings: MultiPolynomial.EvaluationBindings;
roots: REF RootSequence ← NEW[RootSequence[surfaces.length]];
surfaceToDisplay: NAT;
reflectedLight: RGB;
theSurface: STSurface;
theShading: ShadingModels.Model;
Eventually, the following should be passed as parameters
alpha: NAT ← 30;
Define the ray equations.
xRayEquation: Polynomial.Ref ← Polynomial.Linear[[fromPoint.x, direction.x]];
yRayEquation: Polynomial.Ref ← Polynomial.Linear[[fromPoint.y, direction.y]];
zRayEquation: Polynomial.Ref ← Polynomial.Linear[[fromPoint.z, direction.z]];
substitutionBindings: MultiPolynomial.TSUBindings ← LIST[ [$x, xRayEquation], [$y, yRayEquation], [$z, zRayEquation] ];
Iterate over the surfaces. Substitute the ray equations into each surface and find the roots.
FOR i: NAT IN [0..surfaces.length) DO
poly: Polynomial.Ref ← MultiPolynomial.TotallySubstituteUnivariates [surfaces[i].scad.surface, substitutionBindings];
roots[i] ← Polynomial.RealRoots[poly];
ENDLOOP;
Examine the roots in order of increasing value of t, beginning at zero, until one is found that should be displayed. If no such roots are found, return the background color.
t ← 0;
somethingToDisplay ← FALSE;
UNTIL somethingToDisplay DO
Look for all occurrences of the next root.
indicesForNextRoot: LIST OF RECORD[i: NAT, j: NAT] ← NIL;
nextRoot: REAL ← infinity;
FOR i: NAT IN [0..roots.length) DO
FOR j: NAT IN [0..roots[i].nRoots) DO
IF t < roots[i][j] AND roots[i][j] < nextRoot THEN BEGIN
nextRoot ← roots[i][j];
indicesForNextRoot ← CONS[[i,j], NIL];
END
ELSE IF t < roots[i][j] AND roots[i][j] = nextRoot THEN BEGIN
indicesForNextRoot ← CONS[[i,j], indicesForNextRoot];
END;
ENDLOOP;
ENDLOOP;
If no next root was found, return the background.
IF nextRoot = infinity THEN BEGIN
color ← lighting.background;
RETURN;
END;
Otherwise, calculate the point of intersection and see if any occurrence of the root should be displayed.
t ← nextRoot;
intersection ← Geometry3dVector.Combine[fromPoint, 1, direction, t];
IF intersection.x > 0 AND intersection.y > 0 AND intersection.z > 0 THEN ERROR;
UNTIL indicesForNextRoot = NIL OR somethingToDisplay DO
thisSurface: NAT ← indicesForNextRoot.first.i;
cellIsVisible: BOOLEANTRUE;
IF doCells THEN BEGIN
thisCell: NAT;
thisCell ← CADIO.WhatCellIsItOn[intersection, variables, surfaces[thisSurface].scad
! CADIO.Error => {cellIsVisible ← FALSE; CONTINUE}
];
IF cellIsVisible THEN cellIsVisible ← surfaces[thisSurface].mask[thisCell];
END;
IF cellIsVisible THEN BEGIN
somethingToDisplay ← TRUE;
surfaceToDisplay ← thisSurface;
theSurface ← surfaces[surfaceToDisplay];
theShading ← theSurface.shading;
END;
indicesForNextRoot ← indicesForNextRoot.rest;
ENDLOOP;
ENDLOOP;
Calculate the unit vector normal to the surface at the point of intersection.
gradientBindings ← LIST[ [$x, intersection.x], [$y, intersection.y], [$z, intersection.z] ];
xGrad ← MultiPolynomial.TotallyEvaluate[theSurface.gradientX, gradientBindings];
yGrad ← MultiPolynomial.TotallyEvaluate[theSurface.gradientY, gradientBindings];
zGrad ← MultiPolynomial.TotallyEvaluate[theSurface.gradientZ, gradientBindings];
gradient ← [x: xGrad, y: yGrad, z: zGrad];
normal ← Geometry3dVector.Normalize[gradient];
Calculate the total reflected light, beginning with the ambient source.
reflectedLight ← ScaleRGB[lighting.ambientSource, theShading.ambientReflectionCoefficient];
Add in the reflection from each point source.
surfaceToEye ← Geometry3dVector.Sub[fromPoint, intersection];
FOR source: NAT IN [0..lighting.pointSources.length) DO
See if light can reach the surface.
surfaceToSource: Geometry3dVector.Triple ← Geometry3dVector.Sub[lighting.pointSources[source].position, intersection];
xRayEquation: Polynomial.Ref ← Polynomial.Linear[[intersection.x, surfaceToSource.x]];
yRayEquation: Polynomial.Ref ← Polynomial.Linear[[intersection.y, surfaceToSource.y]];
zRayEquation: Polynomial.Ref ← Polynomial.Linear[[intersection.z, surfaceToSource.z]];
substitutionBindings: MultiPolynomial.TSUBindings ← LIST[ [$x, xRayEquation], [$y, yRayEquation], [$z, zRayEquation] ];
poly: Polynomial.Ref ← MultiPolynomial.TotallySubstituteUnivariates [theSurface.scad.surface, substitutionBindings];
t1: REAL;
positiveRootFound, lightCanReachSurface: BOOLEAN;
[t1, positiveRootFound] ← LeastPositiveRoot[poly];
lightCanReachSurface ← (Geometry3dVector.Dot[normal, surfaceToSource] * Geometry3dVector.Dot[normal, surfaceToEye] >= 0) AND ~(positiveRootFound AND (t1 < 1));
If it can, add in the diffuse and specular reflection from that source.
IF lightCanReachSurface THEN BEGIN
cosine: REALABS[Geometry3dVector.Dot[normal, Geometry3dVector.Normalize[surfaceToSource]]];
diffuseReflectance: RGB ← ScaleRGB [
MultiplyCMY [lighting.pointSources[source].color, theShading.surfaceColor],
cosine * theShading.diffuseReflectionCoefficient];
specularReflectance: RGB ← ScaleRGB [
lighting.pointSources[source].color,
IntegerPower[cosine, alpha] * theShading.specularReflectionCoefficient];
totalReflectance: RGB ← AddRGB [
diffuseReflectance,
specularReflectance];
reflectedLight ← AddRGB [
reflectedLight,
totalReflectance];
IF cosine < 0 THEN SIGNAL break;
END;
Do the next source.
ENDLOOP;
Scale down the reflected light if there's too much.
IF reflectedLight.R > 1.0 THEN reflectedLight.R ← 1.0;
IF reflectedLight.G > 1.0 THEN reflectedLight.G ← 1.0;
IF reflectedLight.B > 1.0 THEN reflectedLight.B ← 1.0;
Return the color of the reflected light.
color ← reflectedLight;
END;
Internal procedures
LeastPositiveRoot: PROC[poly: Polynomial.Ref] RETURNS [leastPositiveRoot: REAL, positiveRootFound: BOOLEAN] ~ BEGIN
roots: REF Polynomial.RealRootRec ← Polynomial.RealRoots[poly];
positiveRootFound ← FALSE;
FOR i: INTEGER IN [0..roots.nRoots) DO
IF roots.realRoot[i] > 0.001 THEN
IF positiveRootFound THEN
leastPositiveRoot ← MIN[leastPositiveRoot, roots[i]]
ELSE BEGIN
positiveRootFound ← TRUE;
leastPositiveRoot ← roots[i];
END;
ENDLOOP;
END;
IntegerPower: PROC [base: REAL, exponent: INTEGER] RETURNS [result: REAL] ~ BEGIN
SELECT exponent FROM
0 => result ← 1;
1 => result ← base;
2 => result ← base*base;
3 => result ← base*base*base;
4 => {foo: REAL ← base*base; result ← foo*foo};
5 => {foo: REAL ← base*base; result ← foo*foo*base};
6 => {foo: REAL ← base*base; result ← foo*foo*foo};
>6 => {foo: REAL ← base*base; result ← foo*foo*foo*IntegerPower[base, exponent-6]};
ENDCASE;
END;
AddRGB: PROC [a, b: RGB] RETURNS [sum: RGB] ~ BEGIN
sum ← [
R: a.R + b.R,
G: a.G + b.G,
B: a.B + b.B];
END;
MultiplyRGB: PROC [a, b: RGB] RETURNS [product: RGB] ~ BEGIN
product ← [
R: a.R * b.R,
G: a.G * b.G,
B: a.B * b.B];
END;
ScaleRGB: PROC [a: RGB, s: REAL] RETURNS [product: RGB] ~ BEGIN
product ← [
R: a.R * s,
G: a.G * s,
B: a.B * s];
END;
MultiplyCMY: PROC [a, b: RGB] RETURNS [product: RGB] ~ BEGIN
product ← [
R: 1 - ((1 - a.R) * (1 - b.R)),
G: 1 - ((1 - a.G) * (1 - b.G)),
B: 1 - ((1 - a.B) * (1 - b.B))];
END;
MultiplyCMY: PROC [a, b: RGB] RETURNS [product: RGB] ~ BEGIN
product ← [
R: a.R * b.R,
G: a.G * b.G,
B: a.B * b.B];
END;
MakePicProc: Commander.CommandProc ~ BEGIN
surfaces: REF MultiPolynomial.MultiPolSequence ← NEW[MultiPolynomial.MultiPolSequence[1]];
colors: REF ShadingModels.ShadingSequence ← NEW[ShadingModels.ShadingSequence[1]];
lightingModel: LightingModels.Model ← [
background: [1, 1, 1],
ambientSource: [0.1, 0.1, 0.1],
pointSources: NEW[LightingModels.PointSourceSequence[1]]];
surfaces[0] ← MultiPolynomial.RefFromRope["x^2y^2 + y^2z^2 + x^2z^2 + xyz"];
surfaces[0] ← MultiPolynomial.RefFromRope["x^2 + y^2 + z^2 - 0.11111111"];
surfaces[0] ← MultiPolynomial.RefFromRope["256 z^3 - 128 x^2 z^2 + 144 x y^2 z + 16 x^4 z - 27 y^4 - 4 x^3 y^2"];
colors[0] ← [
surfaceColor: [1, 1, 0],
ambientReflectionCoefficient: 1.0,
diffuseReflectionCoefficient: 0.4,
specularReflectionCoefficient: 0.8];
colors[1] ← [
surfaceColor: [0, 1, 0],
ambientReflectionCoefficient: 1.0,
diffuseReflectionCoefficient: 0.4,
specularReflectionCoefficient: 0.8];
lightingModel.pointSources[0] ← [
position: [2, 0, 0],
color: [1, 1, 1]];
TraceSurfaces[
surfaces: surfaces,
variables: [$x, $y, $z],
colors: colors,
screenCenter: [2, 0, 0],
screenU: Geometry3dVector.Mul[Geometry3dVector.Normalize[[0, 1, 0]], 1],
screenV: Geometry3dVector.Mul[Geometry3dVector.Normalize[[0, 0, 1]], 1],
pixelsU: 100,
pixelsV: 100,
lightingModel: lightingModel,
filename: "FooSurface"];
END;
Commander.Register[key: "MakePic", proc: MakePicProc, doc: "makes an AIS file."];
END.