ImplicitConvolveImpl.mesa
Copyright Ó 1990 by Xerox Corporation. All rights reserved.
Bloomenthal, April 20, 1993 4:40 pm PDT
DIRECTORY Commander, CtBasic, CtFilter, CtMod, G2dBasic, G2dScan, G2dVector, G3dBasic, G3dDraw, G3dMatrix, G3dPlane, G3dPolygon, G3dVector, Imager, ImagerColor, ImagerSample, ImplicitConvolve, ImplicitDebug, ImplicitDesign, IO, Real, RealFns, Rope, ImagerPath, ImagerMaskCapture, ImagerTransformation, SF;
ImplicitConvolveImpl: CEDAR PROGRAM
IMPORTS Commander, CtBasic, CtFilter, CtMod, G2dBasic, G2dScan, G2dVector, G3dDraw, G3dPlane, G3dPolygon, G3dVector, Imager, ImagerSample, ImplicitDebug, ImplicitDesign, IO, Real, RealFns
EXPORTS ImplicitConvolve
~ BEGIN
Types and Constants
Filter:     TYPE ~ CtFilter.Filter;
Box:     TYPE ~ G3dBasic.Box;
Pair:     TYPE ~ G3dBasic.Pair;
IntegerPair:   TYPE ~ G3dBasic.IntegerPair;
Quad:     TYPE ~ G3dBasic.Quad;
RealSequence:  TYPE ~ G3dBasic.RealSequence;
Triple:    TYPE ~ G3dBasic.Triple;
TripleSequence:  TYPE ~ G3dBasic.TripleSequence;
DrawType:   TYPE ~ G3dDraw.DrawType;
Matrix:    TYPE ~ G3dMatrix.Matrix;
Viewport:    TYPE ~ G3dMatrix.Viewport;
Context:    TYPE ~ Imager.Context;
RGB:     TYPE ~ ImagerColor.RGB;
SampleMap:   TYPE ~ ImagerSample.SampleMap;
PairList:    TYPE ~ ImplicitConvolve.PairList;
Primitive:    TYPE ~ ImplicitConvolve.Primitive;
PrimitiveList:  TYPE ~ ImplicitConvolve.PrimitiveList;
PrimitiveRep:  TYPE ~ ImplicitConvolve.PrimitiveRep;
ThicknessRep:  TYPE ~ ImplicitConvolve.ThicknessRep;
Twist:     TYPE ~ ImplicitConvolve.Twist;
ROPE:     TYPE ~ Rope.ROPE;
defaultFilter:   Filter ¬ CtFilter.FilterFind["cubic"];
debug:    BOOL ¬ FALSE;
Public Procedures
MakePrimitive: PUBLIC PROC [
points:  TripleSequence,
extent:  REAL ¬ 0.1,
thickness: REAL ¬ 0.0,
intensity:  REAL ¬ 1.0,
intensities: RealSequence,
indices:  LIST OF INTEGER,
res:   IntegerPair ¬ [50, 50],
color:   RGB ¬ [1.0, 1.0, 1.0],
twist:   Twist ¬ []]
RETURNS [p: Primitive]
~ {
p ¬ NEW[PrimitiveRep ¬ [
points: points,
extent: extent,
intensity: intensity,
intensities: intensities,
res: res,
color: color,
twist: twist,
nVertices: points.length,
filter: defaultFilter
]];
IF thickness # 0 THEN p.thickness ¬ NEW[ThicknessRep ¬ [value: thickness]];
IF p.twist.p0 = p.twist.p1 THEN p.twist.tw0 ¬ p.twist.tw1 ¬ 0.0;
IF indices # NIL THEN {
nVerts: NAT ¬ 0;
FOR l: LIST OF INTEGER ¬ indices, l.rest WHILE l # NIL DO nVerts ¬ nVerts+1; ENDLOOP;
p.indices ¬ NEW[G2dBasic.NatSequenceRep[nVerts]];
FOR n: NAT IN [0..p.indices.length ¬ nVerts) DO
p.indices[n] ¬ indices.first;
indices ¬ indices.rest;
ENDLOOP;
p.nVertices ¬ nVerts;
};
SetGeometry[p];
};
GetBounds: PROC [p: Primitive] RETURNS [b: Box] ~ {
extentF: REAL ¬ p.extent*2.0*1.0*1.8; -- usual extent factor*fudge*sqrt(3) (if on diag)
huge: REAL ~ Real.LargestNumber;
b ¬ [[huge, huge, huge], [-huge, -huge, -huge]];
FOR n: NAT IN [0..p.nVertices) DO
t: Triple ¬ GetPoint[p, n];
b.min ¬ [MIN[b.min.x, t.x], MIN[b.min.y, t.y], MIN[b.min.z, t.z]];
b.max ¬ [MAX[b.max.x, t.x], MAX[b.max.y, t.y], MAX[b.max.z, t.z]];
ENDLOOP;
b.min ¬ G3dVector.Sub[b.min, [extentF, extentF, extentF]];
b.max ¬ G3dVector.Add[b.max, [extentF, extentF, extentF]];
};
GetPoint: PUBLIC PROC [p: Primitive, n: INTEGER] RETURNS [Triple] ~ {
RETURN[p.points[IF p.indices # NIL THEN p.indices[n] ELSE n]];
};
SetGeometry: PUBLIC PROC [p: Primitive] ~ {
NearnessAccelerator: PROC [p0, p1: Triple] RETURNS [v: Triple, w: REAL] ~ {
v ¬ G3dVector.Sub[p1, p0];
IF (w ¬ v.x*v.x+v.y*v.y+v.z*v.z) # 0.0 THEN v ¬ [v.x/w, v.y/w, v.z/w];
};
SELECT p.nVertices FROM
> 2 => {
PairsFromProjection: PROC RETURNS [list: PairList] ~ {
FOR n: NAT IN [0..p.nVertices) DO
q: Triple ¬ GetPoint[p, n];
t: Triple ¬ G3dVector.Sub[p.origin, G3dPlane.ProjectPointToPlane[q, pln]];
t: Triple ¬ G3dPlane.ProjectPointToPlane[q, pln];
pp: Pair ¬ [G3dVector.Dot[t, xAxis], G3dVector.Dot[t, yAxis]];
list ¬ CONS[[G3dVector.Dot[t, xAxis], G3dVector.Dot[t, yAxis]], list];
IF debug THEN {
Out[IO.PutFR["q = (%g, %g, %g\n", IO.real[q.x], IO.real[q.y], IO.real[q.z]]];
Out[IO.PutFLR["t = (%g, %g, %g), p = (%g, %g)\n",
LIST[IO.real[t.x], IO.real[t.y], IO.real[t.z], IO.real[pp.x], IO.real[pp.y]]]];
};
ENDLOOP;
};
normal: Triple ¬ p.normal ¬ G3dPolygon.PolygonNormal[p.points, p.indices, TRUE];
pln: G3dPlane.Plane ¬ G3dPlane.FromPointAndNormal[GetPoint[p, 0], normal];
plane: Quad ¬ p.plane ¬ [pln.x, pln.y, pln.z, pln.w];
xAxis: Triple ¬ p.xAxis ¬ IF p.twist.tw0 # 0 OR p.twist.tw1 # 0 -- x axis in the plane
THEN G3dVector.Unit[G3dVector.Sub[p.twist.p1, p.twist.p0]]
ELSE G3dVector.Ortho[normal];
yAxis: Triple ¬ p.yAxis ¬ G3dVector.Cross[xAxis, normal];
origin: Triple ¬ G3dVector.Mul[normal, -p.plane.w];    -- origin of plane
mm: G2dBasic.Box ¬ MinMaxOfPairs[PairsFromProjection[]];  -- 2d image bounds
min: Pair ¬ G2dVector.Sub[mm.min, [2*p.extent, 2*p.extent]];  -- extent margin
p.size ¬ G2dVector.Add[[4*p.extent, 4*p.extent], G2dVector.Sub[mm.max, mm.min]];
p.scale ¬ MIN[REAL[p.res.x]/p.size.x, REAL[p.res.y]/p.size.y];
p.origin ¬ G3dVector.Add[origin, G3dVector.Combine[xAxis, min.x, yAxis, min.y]];
IF debug THEN {
PrintTriple["xAxis", p.xAxis];
PrintTriple["yAxis", p.yAxis];
Out[IO.PutFLR["size: (%g, %g), scale: %g, extent = %g\n",
LIST[IO.real[p.size.x], IO.real[p.size.y], IO.real[p.scale], IO.real[p.extent]]]];
Out[IO.PutFR["min: (%g, %g)\n", IO.real[min.x], IO.real[min.y]]];
PrintTriple["origin", p.origin];
Out[IO.PutFLR["plane: (%g, %g, %g, %g)\n",
LIST[IO.real[p.plane.x], IO.real[p.plane.y], IO.real[p.plane.z], IO.real[p.plane.w]]]];
};
If obviate GetImageProjection scale:
p.xAxis ¬ G3dVector.Mul[xAxis, p.scale];
p.yAxis ¬ G3dVector.Mul[yAxis, p.scale];
p.normal ¬ G3dVector.Mul[normal, p.scale];
If using major plane projection in 2d:
p.major ¬ G3dPlane.GetMajorPlane[p.plane];
p.xAxis ¬ G3dPlane.ProjectPointToMajorPlane[p.xAxis, p.major];
p.yAxis ¬ G3dPlane.ProjectPointToMajorPlane[p.yAxis, p.major];
IF p.twist.tw0 # 0 OR p.twist.tw1 # 0 THEN {
[Artwork node; type 'Artwork on' to command tool]
p.p0x ¬ G3dVector.Add[
p.twist.p0, G3dVector.Project[G3dVector.Sub[p.origin, p.twist.p0], xAxis]];
p.p1x ¬ G3dVector.Add[p.p0x, G3dVector.Mul[xAxis, p.size.x]];
p.h ¬ G3dVector.Dot[yAxis, G3dVector.Sub[p.twist.p0, p.origin]];
[p.accV, p.accW] ¬ NearnessAccelerator[p.twist.p0, p.twist.p1];
};
};
2 => {
p0: Triple ¬ GetPoint[p, 0];
p1: Triple ¬ GetPoint[p, 1];
axis: Triple ¬ G3dVector.Unit[G3dVector.Sub[p1, p0]];
p.origin ¬ p0 ¬ G3dVector.ScaleRay[[p0, axis], -2.0*p.extent];
p1 ¬ G3dVector.ScaleRay[[p1, axis], 2.0*p.extent];
[p.accV, p.accW] ¬ NearnessAccelerator[p0, p1];
p.res.y ¬ 1;
p.scale ¬ REAL[p.res.x]/G3dVector.Distance[p1, p0];
IF debug THEN
Out[IO.PutFR["segment, scale: %g, extent = %g\n", IO.real[p.scale], IO.real[p.extent]]];
};
1 => NULL;
ENDCASE => ERROR;
p.recipExtent ¬ 1.0/p.extent;
p.filter.support ¬ p.filter.blur ¬ p.scale*p.extent;
p.bounds ¬ GetBounds[p];
};
GetImageProjection: PUBLIC PROC [q: Triple, p: Primitive]
RETURNS [xy: Pair, distance: REAL]
~ {
IF p.twist.tw0 # 0 OR p.twist.tw1 # 0
THEN {
a: REAL ¬ G3dVector.Dot[p.accV, G3dVector.Sub[q, p.twist.p0]];
tw: REAL ¬ p.twist.tw0+a*(p.twist.tw1-p.twist.tw0);
n: Triple ¬ G3dVector.Combine[p.normal, RealFns.Cos[tw], p.yAxis, RealFns.Sin[tw]];
y: Triple ¬ G3dVector.Cross[p.xAxis, n];
o: Triple ¬ G3dVector.Add[p.twist.p0, G3dVector.Mul[p.accV, a*p.accW]];
plane: G3dPlane.Plane ¬ G3dPlane.FromPointAndNormal[o, n, TRUE];
qq: Triple ¬ G3dPlane.ProjectPointToPlane[q, plane];
xy.x ¬ G3dVector.Distance[o, p.p0x];
xy.y ¬ p.h+G3dVector.Dot[G3dVector.Sub[qq, o], y];
IF NOT p.distanceSet THEN p.distance ¬ q.x*plane.x+q.y*plane.y+q.z*plane.z+plane.w;
}
ELSE {
If using projection to major plane:
v: Pair ¬ G3dPlane.ProjectPointToMajorPlane[G3dVector.Sub[q, p.origin], p.major];
xy ¬ [G2dVector.Dot[v, p.Axis], G2dVector.Dot[v, p.yAxis]];
v: Triple ¬ G3dVector.Sub[q, p.origin];
xy ¬ [G3dVector.Dot[v, p.xAxis], G3dVector.Dot[v, p.yAxis]];
IF NOT p.distanceSet
THEN p.distance ¬ q.x*p.plane.x+q.y*p.plane.y+q.z*p.plane.z+p.plane.w;
};
p.distanceSet ¬ TRUE;
distance ¬ p.distance;
xy ¬ G2dVector.Mul[xy, p.scale];
};
MinMaxOfPairs: PROC [pairs: PairList] RETURNS [b: G2dBasic.Box] ~ {
b ¬ [[Real.LargestNumber, Real.LargestNumber], [-Real.LargestNumber, -Real.LargestNumber]];
FOR l: PairList ¬ pairs, l.rest WHILE l # NIL DO
b.min ¬ [MIN[b.min.x, l.first.x], MIN[b.min.y, l.first.y]];
b.max ¬ [MAX[b.max.x, l.first.x], MAX[b.max.y, l.first.y]];
ENDLOOP;
};
ValueOfPrimitive: PUBLIC PROC [q: Triple, p: Primitive] RETURNS [v: REAL ¬ 0.0] ~ {
distance, intensity: REAL ¬ 1.0;
IF NOT p.active THEN RETURN;
IF q.x < p.bounds.min.x OR q.y < p.bounds.min.y OR q.z < p.bounds.min.z OR
q.x > p.bounds.max.x OR q.y > p.bounds.max.y OR q.z > p.bounds.max.z THEN RETURN;
SELECT p.nVertices FROM
> 2 => IF p.image # NIL THEN {
xy: Pair;
ix, iy: NAT;
[xy, distance] ¬ GetImageProjection[q, p];
IF p.oneSided AND distance < 0.0 THEN RETURN[0.0];
IF xy.x < 0 OR xy.y < 0 OR xy.x >= p.res.x OR xy.y >= p.res.y THEN RETURN[0.0];
IF NOT p.interp THEN {ix ¬ Real.Round[xy.x]; iy ¬ Real.Round[xy.y]};
intensity ¬ (1.0/255.0)*(IF p.interp
THEN CtMod.BiLerpReal[p.image, xy.x, xy.y]
ELSE REAL[p.image[iy][ix]]);
IF p.thickness # NIL AND p.thickness.array # NIL THEN
distance ¬ MAX[0.0, ABS[distance]-p.thickness.value*(0.5/255.0)*(IF p.interp
THEN CtMod.BiLerpReal[p.thickness.array, xy.x, xy.y]
ELSE REAL[p.thickness.array[iy][ix]])];
};
2 => IF p.buffer # NIL THEN {
Lerp: PROC [buffer: ImagerSample.SampleBuffer, x: REAL] RETURNS [l: REAL] ~ {
ix0, ix1: INT16 ¬ Real.Floor[x];
l ¬ buffer[ix0];
IF (ix1 ¬ ix0+1) < buffer.length THEN l ¬ l+(x-ix0)*(buffer[ix1]-l);
};
a: REAL ¬ G3dVector.Dot[p.accV, G3dVector.Sub[q, p.origin]];
IF a NOT IN (0.0..1.0) THEN RETURN[0.0];
distance ¬ G3dVector.Distance[q, G3dVector.ScaleRay[[p.origin, p.accV], a*p.accW]];
intensity ¬ (1.0/255.0)*(IF p.interp
THEN Lerp[p.buffer, 0.5+(p.res.x-1)*a]
ELSE REAL[p.buffer[Real.Round[0.5+(p.res.x-1)*a]]]);
};
1 => distance ¬ G3dVector.Distance[q, GetPoint[p, 0]];
ENDCASE => ERROR;
v ¬ p.valueScale*intensity*p.filter.evalProc[distance*p.recipExtent];
Out[IO.PutFR["value = %g\n", IO.real[v]]];
};
ValueOfPrimitives: PUBLIC PROC [q: Triple, primitives: PrimitiveList]
RETURNS [v: REAL ¬ 0.0]
~ {
FOR l: PrimitiveList ¬ primitives, l.rest WHILE l # NIL DO
v ¬ v+ValueOfPrimitive[q, l.first];
ENDLOOP;
};
ConvolvePrimitive: PUBLIC PROC [p: Primitive, filter: BOOL] ~ {
IF p.nVertices <= 2 THEN p.res.y ¬ 1;
SELECT p.nVertices FROM
> 2 => {
pairs: PairList ¬ NIL;
FOR n: NAT IN [0..p.nVertices) DO -- ignore twist when making image
v: Triple ¬ G3dVector.Sub[GetPoint[p, n], p.origin];
pairs ¬ CONS[[G3dVector.Dot[v, p.xAxis], G3dVector.Dot[v, p.yAxis]], pairs];
ENDLOOP;
SetImagePath[p, pairs, filter];
};
2 => {
x0: REAL ¬ 2.0*p.scale*p.extent;
x1: REAL ¬ REAL[p.res.x-1]-2.0*p.scale*p.extent;
SetImagePath[p, LIST[[x0, 0], [x1, 0], [x1, 1], [x0, 1]], filter];
};
ENDCASE;
};
GetTwisted: PUBLIC PROC [p: Primitive, q: Triple] RETURNS [Triple] ~ {
IF p.twist.tw0 # 0 OR p.twist.tw1 # 0
THEN {
a: REAL ¬ G3dVector.Dot[p.accV, G3dVector.Sub[q, p.twist.p0]];
tw: REAL ¬ p.twist.tw0+a*(p.twist.tw1-p.twist.tw0);
o: Triple ¬ G3dVector.Add[p.twist.p0, G3dVector.Mul[p.accV, a*p.accW]];
D: Triple ¬ G3dVector.Sub[q, o];
n: Triple ¬ G3dVector.SameLength[D, G3dVector.Cross[D, p.xAxis]];
v: Triple ¬ G3dVector.Combine[D, RealFns.Cos[tw], n, -RealFns.Sin[tw]];
RETURN[G3dVector.Add[o, v]];
}
ELSE RETURN[q];
};
NActivePrimitives: PUBLIC PROC [primitives: PrimitiveList] RETURNS [i: INTEGER ¬ 0] ~ {
FOR l: PrimitiveList ¬ primitives, l.rest WHILE l # NIL DO IF l.first.active THEN i ¬ i+1; ENDLOOP;
};
SetIDs: PUBLIC PROC [list: LIST OF PrimitiveList, startID: INTEGER ¬ 0] ~ {
FOR l: LIST OF PrimitiveList ¬ list, l.rest WHILE l # NIL DO
FOR p: PrimitiveList ¬ l.first, p.rest WHILE p # NIL DO
IF p.first.id = -1 THEN {p.first.id ¬ startID; startID ¬ startID+1};
ENDLOOP;
ENDLOOP;
};
CombinePrimitives: PUBLIC PROC [list: LIST OF PrimitiveList] RETURNS [ret: PrimitiveList] ~ {
FOR l: LIST OF PrimitiveList ¬ list, l.rest WHILE l # NIL DO
FOR p: PrimitiveList ¬ l.first, p.rest WHILE p # NIL DO
IF p.first # NIL THEN ret ¬ CONS[p.first, ret]
ENDLOOP;
ENDLOOP;
};
InvalidateDistance: PUBLIC PROC [primitives: PrimitiveList] ~ {
FOR l: PrimitiveList ¬ primitives, l.rest WHILE l # NIL DO
l.first.distanceSet ¬ FALSE;
ENDLOOP;
};
InsideSolid: PUBLIC PROC [q: Triple, convexSolid: PrimitiveList] RETURNS [BOOL] ~ {
Surface normals are presumed to point towards the outside of the convex solid.
FOR l: PrimitiveList ¬ convexSolid, l.rest WHILE l # NIL DO
l.first.distance ¬ q.x*l.first.plane.x+q.y*l.first.plane.y+q.z*l.first.plane.z+l.first.plane.w;
l.first.distanceSet ¬ TRUE;
IF l.first.distance > 0.0 THEN RETURN[FALSE];
ENDLOOP;
RETURN[TRUE];
};
ScalePairList: PROC [pairs: PairList, s: REAL] RETURNS [PairList] ~ {
FOR l: PairList ¬ pairs, l.rest WHILE l # NIL DO l.first ¬ G2dVector.Mul[l.first, s]; ENDLOOP;
RETURN[pairs];
};
LengthPairList: PROC [pairs: PairList] RETURNS [len: INTEGER ¬ 0] ~ {
FOR l: PairList ¬ pairs, l.rest WHILE l # NIL DO len ¬ len+1; ENDLOOP;
};
CopyPairList: PROC [pairs: PairList] RETURNS [ret: PairList ¬ NIL] ~ {
tmp: PairList ¬ NIL;
FOR l: PairList ¬ pairs, l.rest WHILE l # NIL DO tmp ¬ CONS[l.first, tmp]; ENDLOOP;
FOR l: PairList ¬ tmp, l.rest WHILE l # NIL DO ret ¬ CONS[l.first, ret]; ENDLOOP;
};
Shaded Display
DisplayPrimitive: PUBLIC PROC [p: Primitive, map: SampleMap] ~ {
ShowReflected: PROC [pa: CtBasic.PixelArray, offset: IntegerPair] ~ {
temp: SampleMap ¬ CtBasic.SampleMapFromPixelArray[pa];
CtMod.ReflectH[temp];
ImagerSample.Transfer[map, temp, [offset.y, offset.x]];
};
IF p.image # NIL THEN ShowReflected[p.image, p.display];
IF p.thickness # NIL AND p.thickness.array # NIL
THEN ShowReflected[p.thickness.array, p.thickness.display];
IF p.buffer # NIL THEN {
line: SampleMap ¬ ImagerSample.NewSampleMap[[[0, 0], [1, p.buffer.length]], 8];
ImagerSample.PutSamples[line, [0, 0],, p.buffer, 0, p.buffer.length];
ImagerSample.Transfer[map, line, [p.display.y, p.display.x]];
};
};
MaxIntensity: PROC [intensities: RealSequence] RETURNS [max: REAL ¬ Real.MinusInfinity] ~ {
IF intensities # NIL THEN FOR i: INTEGER IN [0..intensities.length) DO
max ¬ MAX[intensities[i], max];
ENDLOOP;
};
SetImagePath: PUBLIC PROC [p: Primitive, path: PairList, filter: BOOL ¬ TRUE] ~ {
flat: BOOL ¬ TRUE;
local storage to avoid accumulation of "normalization" if repeated calls to SetImagePath:
intensities: RealSequence ¬ G2dBasic.CopyRealSequence[p.intensities];
p.path ¬ ScalePairList[CopyPairList[path], p.scale];
p.valueScale ¬ p.intensity;
IF intensities # NIL THEN {
max: REAL ¬ MaxIntensity[intensities];
FOR i: INT IN [0..intensities.length) DO IF max # intensities[i] THEN flat ¬ FALSE; ENDLOOP;
p.valueScale ¬ max*p.intensity; -- normalize
FOR i: INT IN [0..intensities.length) DO intensities[i] ¬ intensities[i]/max; ENDLOOP;
};
SELECT p.nVertices FROM
> 2 => {
MakeImage: PROC ~ {
i: INTEGER ¬ -1;
pairs: LIST OF G2dBasic.IntegerPair ¬ NIL;
vertices: LIST OF G2dScan.Vertex ¬ NIL;
ImagerSample.Clear[map];
FOR l: PairList ¬ p.path, l.rest WHILE l # NIL DO
p: IntegerPair ¬ [Real.Round[l.first.x], Real.Round[l.first.y]];
IF flat
THEN pairs ¬ CONS[p, pairs]
ELSE vertices ¬ CONS[[p.x, p.y, Real.Round[255.0*intensities[i ¬ i+1]]], vertices];
ENDLOOP;
IF flat
THEN [] ¬ G2dScan.FlatShade[map, pairs, 255]
ELSE [] ¬ G2dScan.GouraudShade[map, vertices];
};
box: ImagerSample.Box ¬ [[0, 0], [p.res.y, p.res.x]];
map: SampleMap ¬ ImagerSample.ObtainScratchMap[box, 8];
MakeImage[];
IF filter THEN {
max: CARDINAL ¬ 0;
buf: ImagerSample.SampleBuffer ¬ ImagerSample.NewSamples[p.res.x];
CtFilter.ResampleMap[map, box, box, p.filter]; -- unnormalized filtering & quantization
FOR y: INT IN [0..p.res.y) DO
ImagerSample.GetSamples[map, [y, 0],, buf, 0, p.res.x];
FOR x: INT IN [0..p.res.x) DO IF buf[x] > max THEN max ¬ buf[x]; ENDLOOP;
ENDLOOP;
IF (max ¬ max+1) IN (0..255) THEN { -- normalization of image:
MakeImage[]; -- inexpensive to recompute
p.filter.scale ¬ 255.0/REAL[max]; -- scale filtered image so its max pixel = 255
p.valueScale ¬ p.valueScale/p.filter.scale; -- compensate during ValueOfPrimitive
CtFilter.ResampleMap[map, box, box, p.filter]; -- normalized filtering (expensive)
p.filter.scale ¬ 1.0; -- reset
};
};
p.image ¬ CtBasic.PixelArrayFromSampleMap[map];
ImagerSample.ReleaseScratchMap[map];
};
2 => {
pairs: PairList ¬ path;
box: ImagerSample.Box ¬ [[0, 0], [1, p.res.x]];
map: SampleMap ¬ ImagerSample.ObtainScratchMap[box, 8];
line: ImagerSample.SampleBuffer ~ ImagerSample.ObtainScratchSamples[p.res.x];
FOR x: INT IN [0..p.res.x) DO
line.samples[x] ¬ IF x IN [pairs.first.x..pairs.rest.first.x] THEN 255 ELSE 0;
ENDLOOP;
IF NOT flat THEN
FOR x: INT IN [0..p.res.x) DO
a: REAL ¬ REAL[x]/REAL[p.res.x-1];
factor: REAL ¬ p.intensities[0]+a*(p.intensities[1]-p.intensities[0]);
line.samples[x] ¬ Real.Round[factor*REAL[line.samples[x]]];
ENDLOOP;
ImagerSample.PutSamples[map, [0, 0],, line, 0, p.res.x];
IF filter THEN CtFilter.ResampleMap[map, box, box, p.filter];
ImagerSample.GetSamples[map,,, p.buffer ¬ ImagerSample.NewSamples[p.res.x]];
ImagerSample.ReleaseScratchSamples[line];
ImagerSample.ReleaseScratchMap[map];
};
1 => NULL;
ENDCASE;
};
Line Drawing
DrawPrimitive: PUBLIC PROC [
p: Primitive,
context: Context,
view: Matrix,
viewport: Viewport ¬ [],
whatChanged: REF ¬ NIL,
type: DrawType,
log: IO.STREAM ¬ NIL]
~ {
T: PROC [i: INT] RETURNS [t: Triple] ~ {t ¬ GetTwisted[p, GetPoint[p, i]]};
close: BOOL ¬ p.nVertices > 2;
G3dDraw.SetColor[context, [p.color.R, p.color.G, p.color.B]];
IF p.twist.tw0 # 0 OR p.twist.tw1 # 0
THEN
G3dDraw.DrawWithConnectedPoints[context, p.nVertices, T, view, viewport, close, type]
ELSE
G3dDraw.ConnectedPoints[context, p.points, view, viewport, close, type];
};
DrawPrimitives: PUBLIC PROC [
list: LIST OF Primitive,
context: Context,
view: Matrix,
viewport: Viewport ¬ [],
whatChanged: REF ¬ NIL,
type: DrawType,
log: IO.STREAM ¬ NIL]
~ {
IF whatChanged = $IPOut THEN Imager.SetStrokeWidth[context, 3.0];
FOR l: LIST OF Primitive ¬ list, l.rest WHILE l # NIL DO
IF l.first.active
THEN DrawPrimitive[l.first, context, view, viewport, whatChanged, type, log];
ENDLOOP;
};
Debug
Out: PROC [message: ROPE] ~ {ImplicitDebug.Write[message]};
PrintTriple: PROC [name: ROPE, t: Triple] ~ {
Out[IO.PutFLR["%g: (%g, %g, %g)\n",
LIST[IO.rope[name], IO.real[t.x], IO.real[t.y], IO.real[t.z]]]];
};
PrintMatrix: PROC [matrix: G3dMatrix.Matrix, name: ROPE ¬ NIL] ~ {
IF out = NIL OR matrix = NIL THEN RETURN;
IF name # NIL THEN IO.PutF[out, "%g\n", IO.rope[name]];
FOR i: NAT IN [0..3] DO
Out[IO.PutFR[out, "%6.3f\t%6.3f\t%6.3f\t%6.3f\n",
IO.real[matrix[i][0]], IO.real[matrix[i][1]], IO.real[matrix[i][2]], IO.real[matrix[i][3]]]];
ENDLOOP;
};
PrintPrimitive: PROC [p: Primitive] ~ {
Out[IO.PutFR[out, "nVertices = %g\n", IO.int[p.nVertices]]];
};
Start Code
Commander.Register["ImplicitConvolve", ImplicitDesign.DesignCmd, "ImplicitConvolve [option | ?]"];
END...
SetImagePath: PUBLIC PROC [p: Primitive, path: PairList, filter: BOOL ¬ TRUE] ~ {
p.path ¬ ScalePairList[CopyPairList[path], p.scale];
SELECT p.nVertices FROM
> 2 => { -- this is an incredibly awkward way to smooth shade a polygon, but . . .
Trajectory: TYPE ~ ImagerPath.Trajectory;
TrajectoryFromList: PROC [pairs: PairList] RETURNS [t: Trajectory] ~ {
t ¬ ImagerPath.MoveTo[pairs.first];
FOR l: PairList ¬ pairs.rest, l.rest WHILE l # NIL DO
t ¬ ImagerPath.LineTo[t, l.first];
ENDLOOP;
};
Bit: PROC [context: Context] ~ {
Imager.SetGray[context, 1];
Imager.MaskRectangle[context, [0, 0, 1, 1]];     -- grotesque: fix minimum
Imager.MaskRectangle[context, [p.res.x-1, p.res.y-1, 1, 1]]; -- grotesque: fix maximum
Imager.MaskFillTrajectory[context, TrajectoryFromList[p.path]];
};
bit: SampleMap ¬ ImagerSample.ZeroOrigin[
ImagerMaskCapture.CaptureBitmap[Bit, ImagerTransformation.Rotate[-90.0]]];
box: ImagerSample.Box ¬ [[0, 0], [p.res.y, p.res.x]];
map: SampleMap ¬ ImagerSample.ObtainScratchMap[box, 8];
size: ImagerSample.Vec ¬ SF.Min[ImagerSample.GetSize[map], ImagerSample.GetSize[bit]];
line: ImagerSample.SampleBuffer ~ ImagerSample.ObtainScratchSamples[size.f];
ImagerSample.Clear[map];
FOR y: NAT IN [1..size.s-1) DO -- avoid "fix" points
ImagerSample.GetSamples[bit, [y, 0],, line, 0, size.f];
FOR x: NAT IN [0..size.f) DO
line.samples[x] ¬ IF line.samples[x] = 1 THEN 255 ELSE 0;
ENDLOOP;
ImagerSample.FlipSamples[line]; -- ever more ugliness
ImagerSample.PutSamples[map, [y, 0],, line, 0, size.f];
ENDLOOP;
IF p.intensities # NIL THEN {
same: BOOL ¬ TRUE;
max: REAL ¬ p.intensities[0];
FOR i: INT IN [1..p.intensities.length) DO
IF max # p.intensities[i] THEN {max ¬ MAX[max, p.intensities[i]]; same ¬ FALSE};
ENDLOOP;
p.intensity ¬ max*p.intensity; -- normalize
FOR i: INT IN [0..p.intensities.length) DO p.intensities[i]¬p.intensities[i]/max; ENDLOOP;
};
IF filter THEN CtFilter.ResampleMap[map, box, box, p.filter];
p.image ¬ CtBasic.PixelArrayFromSampleMap[map];
ImagerSample.ReleaseScratchSamples[line];
ImagerSample.ReleaseScratchMap[map];
};
2 => {
pairs: PairList ¬ path;
box: ImagerSample.Box ¬ [[0, 0], [1, p.res.x]];
map: SampleMap ¬ ImagerSample.ObtainScratchMap[box, 8];
line: ImagerSample.SampleBuffer ~ ImagerSample.ObtainScratchSamples[p.res.x];
FOR x: INT IN [0..p.res.x) DO
line.samples[x] ¬ IF x IN [pairs.first.x..pairs.rest.first.x] THEN 255 ELSE 0;
ENDLOOP;
IF p.intensities # NIL AND p.intensities.length = 2 THEN
FOR x: INT IN [0..p.res.x) DO
a: REAL ¬ REAL[x]/REAL[p.res.x-1];
factor: REAL ¬ p.intensities[1]+a*(p.intensities[1]-p.intensities[0]);
line.samples[x] ¬ Real.Round[factor*REAL[line.samples[x]]];
ENDLOOP;
ImagerSample.PutSamples[map, [0, 0],, line, 0, p.res.x];
IF filter THEN CtFilter.ResampleMap[map, box, box, p.filter];
ImagerSample.GetSamples[map,,, p.buffer ¬ ImagerSample.NewSamples[p.res.x]];
ImagerSample.ReleaseScratchSamples[line];
ImagerSample.ReleaseScratchMap[map];
};
1 => NULL;
ENDCASE;
};
polygon => {
o: Triple ¬ GetPoint[p, 0];
n: Triple ¬ p.normal ¬ G3dPolygon.PolygonNormal[p.points, p.indices, TRUE];
x: Triple ¬ G3dVector.Unit[G3dVector.Ortho[n]];
y: Triple ¬ G3dVector.Unit[G3dVector.Cross[n, x]];
m: Matrix ¬ p.xform ¬ G3dMatrix.Invert[G3dMatrix.MakeFromTriad[x, y, n, o]];
mm: G2dBasic.Box ¬ MinMaxOfPairs[PairsFromXformedPrimitive[p]];
dif: Pair ¬ G2dVector.Add[[4*p.extent, 4*p.extent], G2dVector.Sub[mm.max, mm.min]];
p.scale ¬ MIN[REAL[p.res.x]/dif.x, REAL[p.res.y]/dif.y];
p.plane ¬ G3dPlane.FromPointAndNormal[o, n];
p.xform ¬ m ¬ G3dMatrix.Scale[p.xform, p.scale, p.xform];
mm ¬ MinMaxOfPairs[PairsFromXformedPrimitive[p]];
p.move ¬ [2*p.filter.support-mm.min.x, 2*p.filter.support-mm.min.y];
p.xform ¬ G3dMatrix.Translate[m, [p.move.x, p.move.y, 0], m];
};
PairsFromXformedPrimitive: PROC [p: Primitive] RETURNS [pairs: PairList] ~ {
FOR n: NAT IN [0..p.nVertices) DO
pairs ¬ CONS[XYFrom3d[GetPoint[p, n], p], pairs];
ENDLOOP;
};
XYFrom3d: PROC [point: Triple, p: Primitive] RETURNS [pp: Pair] ~ {
pp.x ¬ point.x*p.xform[0][0]+point.y*p.xform[1][0]+point.z*p.xform[2][0]+p.xform[3][0];
pp.y ¬ point.x*p.xform[0][1]+point.y*p.xform[1][1]+point.z*p.xform[2][1]+p.xform[3][1];
};
scratchReals: ARRAY [0..9] OF RealSequence ¬ ALL[NIL];
ObtainReals: ENTRY PROC [strength: REAL] RETURNS [reals: RealSequence] ~ {
FOR i: NAT IN [0..9] DO
IF (reals ¬ scratchReals[i]) # NIL THEN {scratchReals[i] ¬ NIL; EXIT};
ENDLOOP;
IF reals = NIL THEN reals ¬ NEW[G3dBasic.RealSequenceRep[4]];
reals[0] ¬ reals[1] ¬ reals[2] ¬ strength;
reals[3] ¬ 0.0;
reals.length ¬ 4;
};
ReleaseReals: ENTRY PROC [reals: RealSequence] ~ {
FOR i: NAT IN [0..9] DO IF scratchReals[i] = NIL THEN {scratchReals[i] ¬ reals; EXIT}; ENDLOOP;
};
context: G3dRender.Context;
ConvolvePrimitive: PUBLIC PROC [primitive: Primitive, intensity: REAL,res:IntegerPair¬[50,50]]~{
For polygons with differing vertex intensities
FancyPatch: TYPE ~ G3dRenderWithPixels.FancyPatch;
MakeMap: PROC RETURNS [SampleMap] ~ {
patch: REF FancyPatch ¬ NEW[G3dRenderWithPixels.FancyPatch[primitive.points.length]];
IF context = NIL THEN context ¬ G3dRender.CreateUndisplayedContext[, res.x, res.y, gray];
G3dRender.SetAntiAliasing[context, TRUE];
G3dSortandDisplay.ValidateContext[context];
FOR n: NAT IN [0..primitive.points.length) DO
pair: Pair ¬ G3dMatrix.TransformD[p.points[n], p.xform];
patch[n] ¬ [pair.x, res.y-pair.y, ObtainReals[p.strengths[n]]];
ENDLOOP;
G3dRenderWithPixels.FillInConstantBackGround[context, [0, 0, 0], 0];
G3dRenderWithPixels.RealFancyTiler[context, patch];
FOR n: NAT IN [0..primitive.points.length) DO ReleaseReals[patch[n].val]; ENDLOOP;
RETURN[context.pixels[0]];
};
m: SampleMap;
SetGeometry[primitive, intensity, res];
m ¬ MakeMap[];
CtFilter.ResampleMap[map,ImagerSample.GetBox[m], ImagerSample.GetBox[m], primitive.filter];
primitive.image ¬ CtBasic.PixelArrayFromSampleMap[m];
};
ProjectToMajorPlane: PROC [p: Triple, plane: Quad, majorPlane: MajorPlane] RETURNS [Pair]
~ {
normal: Triple ~ [plane.x, plane.y, plane.z];
distance: REAL ~ G3dVector.Dot[p, normal]+plane.w;
RETURN[SELECT majorPlane FROM
xy =>   G2dVector.Sub[[p.x, p.y], G3dVector.Mul[[plane.x, plane.y], distance]];
xz =>   G2dVector.Sub[[p.x, p.z], G3dVector.Mul[[plane.x, plane.z], distance]];
ENDCASE => G2dVector.Sub[[p.y, p.z], G3dVector.Mul[[plane.y, plane.z], distance]]];
};
IntersectPolygonWithRay: PROC RETURNS [list: LIST OF REAL] ~ {
normal: Triple ¬ G3dVector.Ortho[p.twist.axis.axis];
plane: Quad ¬ G3dPlane.FromPointAndNormal[p.twist.axisbase, normal];
FOR n: NAT IN [0..p.nVertices) DO
p0: Triple ¬ GetPoint[p, n];
p1: Triple ¬ GetPoint[p, (n+1) MOD p.nVertices];
d0: REAL ¬ G3dPlane.DistanceToPoint[p0, plane];
d1: REAL ¬ G3dPlane.DistanceToPoint[p1, plane];
IF d0*d1 < 0.0 THEN {
a0: REAL ¬ ABS[d0];
a1: REAL ¬ ABS[d1];
t: Triple ¬ G3dVector.Combine[p0, a1/(a0+a1), p1, a0/(a+a1)];
};
ENDLOOP;
};