Acceleration
Digital Image Synthesis g g y Yung-Yu Chuang
10/8/2009 10/8/2009
with slides by Mario Costa Sousa, Gordon Stoll and Pat Hanrahan
Classes
• Primitive (in core/primitive.*)
G t i P i iti – GeometricPrimitive – InstancePrimitive
A t
– Aggregate
• Two types of accelerators are provided (in accelerators/*.cpp)
– GridAccel – KdTreeAccel
Hierarchy
P i iti Primitive
I t A t
Geometric Primitive
Instance Primitive
Aggregate
T
Material Shape
Primitive
class Primitive : public ReferenceCounted {
<Primitive interface>
<Primitive interface>
}
class InstancePrimitive : public Primitive {
…
Reference<Primitive> instance;
} }
Interface
BBox WorldBound();
bool CanIntersect();
geometry bool CanIntersect();
bool Intersect(const Ray &r,
Intersection *in);
// update maxt Intersection *in);
bool IntersectP(const Ray &r);
void Refine(vector<Reference<Primitive>>
void Refine(vector<Reference<Primitive>>
&refined);
id F ll R fi ( t <R f <P i iti >>
void FullyRefine(vector<Reference<Primitive>>
&refined);
A Li ht *G tA Li ht() t i l
AreaLight *GetAreaLight();
BSDF *GetBSDF(const DifferentialGeometry &dg, ) material const Transform &WorldToObject);
Intersection
• primitive stores the actual intersecting primitive hence Primitive >GetAreaLight and primitive, hence Primitive->GetAreaLight and GetBSDF can only be called for
GeometricPrimitive GeometricPrimitive
More information than DifferentialGeometry;
struct Intersection {
y;
also contains material information
<Intersection interface>
DifferentialGeometry dg;
const Primitive *primitive;
Transform WorldToObject;
};
GeometricPrimitive
• represents a single shape
h ld f t h d it i l
• holds a reference to a Shape and its Material, and a pointer to an AreaLight
Reference<Shape> shape;
Reference<Material> material; // BRDF AreaLight *areaLight; // emittance
• Most operations are forwarded to shape
Object instancing
61 unique plant models, 1.1M triangles, 300MBq p , g , 4000 individual plants, 19.5M triangles
InstancePrimitive
R f <P i iti > i t
Reference<Primitive> instance;
Transform InstanceToWorld, WorldToInstance;
Ray ray = WorldToInstance(r);
if (!instance->Intersect(ray, isect)) return false;
r.maxt = ray.maxt;
isect->WorldToObject = isect->WorldToObject
*WorldToInstance;
Aggregates
• Acceleration is a heart component of a ray tracer because ray/scene intersection accounts tracer because ray/scene intersection accounts for the majority of execution time
G l d th b f / i iti
• Goal: reduce the number of ray/primitive
intersections by quick simultaneous rejection of f i iti d th f t th t b groups of primitives and the fact that nearby intersections are likely to be found first
• Two main approaches: spatial subdivision, object subdivision
• No clear winner
Acceleration techniques Bounding volume hierarchy
Bounding volume hierarchy
1) Find bounding box of objects
Bounding volume hierarchy
1) Find bounding box of objects 2) S li bj i
2) Split objects into two groups
Bounding volume hierarchy
1) Find bounding box of objects 2) S li bj i
2) Split objects into two groups 3) Recurse
Bounding volume hierarchy
1) Find bounding box of objects 2) S li bj i
2) Split objects into two groups 3) Recurse
Bounding volume hierarchy
1) Find bounding box of objects 2) S li bj i
2) Split objects into two groups 3) Recurse
Bounding volume hierarchy
1) Find bounding box of objects 2) S li bj i
2) Split objects into two groups 3) Recurse
Where to split?
• At midpoint
• Sort and put half of the objects on each side
• Sort, and put half of the objects on each side
• Use modeling hierarchy
BVH traversal
• If hit parent, then check all children
BVH traversal
• Don't return intersection immediately because the other subvolumes may have a closer
the other subvolumes may have a closer intersection
Bounding volume hierarchy
Bounding volume hierarchy Space subdivision approaches
Quadtree (2D) Unifrom grid Quadtree (2D)
Octree (3D)
Unifrom grid
Space subdivision approaches
KD tree BSP tree
KD tree BSP tree
Uniform grid
Uniform grid
P
Preprocess scene 1. Find bounding box
Uniform grid
P
Preprocess scene 1. Find bounding box 2 Determine grid resolution 2. Determine grid resolution
Uniform grid
P
Preprocess scene 1. Find bounding box 2 Determine grid resolution 2. Determine grid resolution 3. Place object in cell if its
bounding box overlaps the cell
Uniform grid
P
Preprocess scene 1. Find bounding box 2 Determine grid resolution 2. Determine grid resolution 3. Place object in cell if its
bounding box overlaps the cell 4. Check that object overlaps cell
(expensive!)
Uniform grid traversal
P
Preprocess scene Traverse grid
3D li 3D DDA 3D line = 3D-DDA (Digital Differential Analyzer)
Analyzer)
1 2
x x
y m y b
mx y
1
2
x
x
1
1
i
i
x
x
b mx
y
i1
i1
naivem
y
y
i1
i
DDAoctree
Octree
K-d tree
A
A
Leaf nodes correspond to unique regions in space
K-d tree
A
B
A
Leaf nodes correspond to unique regions in space
K-d tree
A
B B
B
A
Leaf nodes correspond to unique regions in space
K-d tree
C A
B B
B
A
K-d tree
C A
B B
C B
C
A
K-d tree
C A
D BB
C B
C
A
K-d tree
C A
D BB
C B
C
A
D
K-d tree
A
D B
B CC CC
D A
Leaf nodes correspond to unique regions in space
K-d tree traversal
A
D B
B CC CC
D
A
Leaf nodes correspond to unique regions in space
BSP tree
6 55
9 7
9 10 8
1 2
11
4 2
3
BSP tree
6 5
5 11
inside outside 9
7
inside ones
outside ones 9
10 8 1
2
11
4 2
3
BSP tree
6 55 11
2 5
9 6 7
3 4
6 7 8 9
10 8 1
2
9 10 11 11
4 2
3
BSP tree
6 5 1
9 5
9b 5
1
9 10 8
7
9a
6 7
11b 8
9b 1
2
11
9a 10 11a 11a
11b
4
2 11a
3
BSP tree
6 5
9 9b
2 1
5
1
9 10 8
7
2 3
5
6 8
9a
11b 1
2
11a 11 4 7 9b
4 2
9a 11b
3 10
11a
BSP tree traversal
6 5
9 5
9b 2
1 5
1
9 10 8
7
2 3
5
6 8
9a
11b 9a
11b 1
2
11a 11 4 7 9b
4 2
i t
9a 11b
3 point
10 11a
BSP tree traversal
6 5
9 5
9b 2
1 5
1
9 10 8
7
2 3
5
6 8
9a
11b 9a
11b 1
2
11a 11 4 7 9b
4 2
i t
9a 11b
3 point
10 11a
BSP tree traversal
6 5
9 9b
2 1
5
1
9 10 8
7
2 3
5
6 8
9a
11b 9a
11b 1
2
11a 11 4 7 9b
2
point
9a 11b 4
point 3 10
11a
Ray-Box intersections
• Both GridAccel and KdTreeAccel require it Q i k j i d i i
• Quick rejection, use enter and exit point to traverse the hierarchy
• AABB is the intersection of three slabs
Ray-Box intersections
1 Dx
O
xt x
1t
1D
xt
1 O
xx
1
x=x0 x=x1
Ray-Box intersections
bool BBox::IntersectP(const Ray &ray,
float *hitt0, float *hitt1) {
float t0 = ray.mint, t1 = ray.maxt;
for (int i = 0; i < 3; ++i) {
float invRayDir = 1.f / ray.d[i];
float tNear = (pMin[i] - ray.o[i]) * invRayDir;
float tFar = (pMax[i] - ray.o[i]) * invRayDir;(p [ ] y [ ]) y ; if (tNear > tFar) swap(tNear, tFar);
t0 = tNear > t0 ? tNear : t0;
i i
t0 = tNear > t0 ? tNear : t0;
t1 = tFar < t1 ? tFar : t1;
if (t0 > t1) return false;
}
segment intersection intersection is empty
}
if (hitt0) *hitt0 = t0;
if (hitt1) *hitt1 = t1;
return true;
}
Grid accelerator
• Uniform grid
Teapot in a stadium problem
• Not adaptive to distribution of primitives.
Have to determine the number of voxels
• Have to determine the number of voxels.
(problem with too many or too few)
GridAccel
Class GridAccel:public Aggregate {
<GridAccel methods>
<GridAccel methods>
u_int nMailboxes;
MailboxPrim *mailboxes;
MailboxPrim *mailboxes;
int NVoxels[3];
BBox bounds;
BBox bounds;
Vector Width, InvWidth;
V l ** l
Voxel **voxels;
ObjectArena<Voxel> voxelArena;
t ti i t M ilb Id static int curMailboxId;
}
mailbox
struct MailboxPrim {
Reference<Primitive> primitive;
Reference<Primitive> primitive;
Int lastMailboxId;
}
GridAccel
GridAccel(vector<Reference<Primitive> > &p, bool forRefined bool refineImmediately) bool forRefined, bool refineImmediately) : gridForRefined(forRefined) {
// Initialize with primitives for grid // Initialize with primitives for grid vector<Reference<Primitive> > prims;
if (refineImmediately) if (refineImmediately)
for (u_int i = 0; i < p.size(); ++i) [i] >F ll R fi ( i )
p[i]->FullyRefine(prims);
else i
prims = p;
GridAccel
// Initialize mailboxes for grid nMailboxes = prims size();
nMailboxes = prims.size();
mailboxes
= (MailboxPrim*)AllocAligned(
= (MailboxPrim*)AllocAligned(
nMailboxes * sizeof(MailboxPrim));
for (u int i 0; i < nMailboxes; ++i) for (u_int i = 0; i < nMailboxes; ++i)
new (&mailboxes[i])MailboxPrim(prims[i]);
Determine number of voxels
• Too many voxels → slow traverse, large
memory consumption (bad cache performance) memory consumption (bad cache performance)
• Too few voxels → too many primitives in a l
voxel
• Let the axis with the largest extent have
3 N3partitions (N:number of primitives)
Vector delta = bounds.pMax - bounds.pMin;
int maxAxis=bounds.MaximumExtent();
float invMaxWidth=1.f/delta[maxAxis];
float cubeRoot=3.f*powf(float(prims.size()),1.f/3.f);
float voxelsPerUnitDist=cubeRoot * invMaxWidth;
float voxelsPerUnitDist=cubeRoot * invMaxWidth;
Calculate voxel size and allocate voxels
for (int axis=0; axis<3; ++axis) {
NVoxels[axis]=Round2Int(delta[axis]*voxelsPerUnitDist);
NVoxels[axis]=Clamp(NVoxels[axis], 1, 64);
}
for (int axis=0; axis<3; ++axis) {
Width[axis]=delta[axis]/NVoxels[axis];
Width[axis]=delta[axis]/NVoxels[axis];
InvWidth[axis]=
(Width[axis]==0.f)?0.f:1.f/Width[axis];
}
int nVoxels = NVoxels[0] * NVoxels[1] * NVoxels[2];
voxels=(Voxel **)AllocAligned(nVoxels*sizeof(Voxel *));
memset(voxels 0 nVoxels * sizeof(Voxel *));
memset(voxels, 0, nVoxels * sizeof(Voxel *));
Conversion between voxel and position
int PosToVoxel(const Point &P, int axis) { int v=Float2Int(
(P[axis]-bounds.pMin[axis])*InvWidth[axis]);
return Clamp(v, 0, NVoxels[axis]-1);
}
float VoxelToPos(int p, int axis) const { return bounds pMin[axis]+p*Width[axis];
return bounds.pMin[axis]+p Width[axis];
}
Point VoxelToPos(int x, int y, int z) const { return bounds pMin+
return bounds.pMin+
Vector(x*Width[0], y*Width[1], z*Width[2]);
}
inline int Offset(int x, int y, int z) {
return z*NVoxels[0]*NVoxels[1] + y*NVoxels[0] + x;
}
Add primitives into voxels
for (u_int i=0; i<prims.size(); ++i) {
<Find voxel extent of primitive>
<Find voxel extent of primitive>
<Add primitive to overlapping voxels>
} }
<Find voxel extent of primitive>
BBox pb = prims[i]->WorldBound();
int vmin[3] vmax[3];
int vmin[3], vmax[3];
for (int axis = 0; axis < 3; ++axis) { vmin[axis] = PosToVoxel(pb pMin axis);
vmin[axis] = PosToVoxel(pb.pMin, axis);
vmax[axis] = PosToVoxel(pb.pMax, axis);
} }
<Add primitive to overlapping voxels>
for (int z = vmin[2]; z <= vmax[2]; ++z) for (int y = vmin[1]; y <= vmax[1]; ++y)
for (int x = vmin[0]; x <= vmax[0]; ++x) { int offset = Offset(x, y, z);
if (!voxels[offset]) { voxels[offset]
= new (voxelArena)Voxel(&mailboxes[i]);
} } else {
// Add primitive to already-allocated voxel voxels[offset]->AddPrimitive(&mailboxes[i]);
voxels[offset] >AddPrimitive(&mailboxes[i]);
} }
Voxel structure
struct Voxel {
<Voxel methods>
<Voxel methods>
union {
MailboxPrim *onePrimitive;
MailboxPrim *onePrimitive;
MailboxPrim **primitives;
};
};
u_int allCanIntersect:1;
i t P i iti 31 u_int nPrimitives:31;
}
Packed into 64 bits
AddPrimitive
void AddPrimitive(MailboxPrim *prim) { if (nPrimitives == 1) {
MailboxPrim **p = new MailboxPrim *[2];
p[0] = onePrimitive;
primitives = p;
}
else if (IsPowerOf2(nPrimitives)) { int nAlloc = 2 * nPrimitives;;
MailboxPrim **p = new MailboxPrim *[nAlloc];
for (u_int i = 0; i < nPrimitives; ++i) p[i] = primitives[i];
p[i] = primitives[i];
delete[] primitives;
primitives = p;
} }
primitives[nPrimitives] = prim;
++nPrimitives;
} }
GridAccel traversal
bool GridAccel::Intersect(
Ray &ray Intersection *isect) { Ray &ray, Intersection *isect) {
<Check ray against overall grid bounds>
<Get ray mailbox id>
<Get ray mailbox id>
<Set up 3D DDA for ray>
<Walk ray through voxel grid>
<Walk ray through voxel grid>
}
Check against overall bound
float rayT;
if (bounds Inside(ray(ray mint))) if (bounds.Inside(ray(ray.mint)))
rayT = ray.mint;
else if (!bounds IntersectP(ray &rayT)) else if (!bounds.IntersectP(ray, &rayT))
return false;
Point gridIntersect ray(rayT);
Point gridIntersect = ray(rayT);
Set up 3D DDA (Digital Differential Analyzer)
• Similar to Bresenhm’s line drawing algorithm
Set up 3D DDA (Digital Differential Analyzer)
blue values changes along the traversal
NextCrossingT[1] Out
voxel
g g
voxel index
DeltaT[0]
rayT
Step[0]=1 DeltaT[0]
NextCrossingT[0]
Pos
DeltaT: the distance change when voxel changes 1 in that direction
Set up 3D DDA
for (int axis=0; axis<3; ++axis) {
Pos[axis]=PosToVoxel(gridIntersect, axis);
if (ray.d[axis]>=0) {
NextCrossingT[axis] = rayT+
(VoxelToPos(Pos[axis]+1 axis)-gridIntersect[axis]) (VoxelToPos(Pos[axis]+1,axis)-gridIntersect[axis]) /ray.d[axis];
DeltaT[axis] = Width[axis] / ray.d[axis];
Step[axis] = 1;
Out[axis] = NVoxels[axis]; 1
Out[axis] NVoxels[axis];
} else { ...
Dx
Step[axis] = -1;
Out[axis] = -1;
} }
} Width[0]
Walk through grid
for (;;) {
*voxel=voxels[Offset(Pos[0] Pos[1] Pos[2])];
*voxel=voxels[Offset(Pos[0],Pos[1],Pos[2])];
if (voxel != NULL) hitSomething |=
hitSomething |=
voxel->Intersect(ray,isect,rayId);
<Advance to next voxel>
<Advance to next voxel>
}
t hitS thi
return hitSomething;
Do not return; cut tmax instead Do not return; cut tmax instead Return when entering a voxel
that is beyond the closest found intersection.
Advance to next voxel
int bits=((NextCrossingT[0]<NextCrossingT[1])<<2) + ((NextCrossingT[0]<NextCrossingT[2])<<1) + ((NextCrossingT[1]<NextCrossingT[2]));
const int cmpToAxis[8] = { 2, 1, 2, 1, 2, 2, 0, 0 };
int stepAxis=cmpToAxis[bits];
if (ray.maxt < NextCrossingT[stepAxis]) break;
Pos[stepAxis]+=Step[stepAxis];
Pos[stepAxis]+ Step[stepAxis];
if (Pos[stepAxis] == Out[stepAxis]) break;
NextCrossingT[stepAxis] += DeltaT[stepAxis];
conditions
x<y x<z y<z
0 0 0 x≥y≥z 2
0 0 1 x≥z>y 1
0 1 0 -
0 1 1 z>x≥y 1 1 0 0 y>x≥z 2 1 0 0 y>x≥z 2
1 0 1 -
1 1 0 y≥z>x 0 1 1 0 y≥z>x 0 1 1 1 z>y>x 0
KD-Tree accelerator
• Non-uniform space subdivision (for example, kd tree and octree) is better than uniform grid kd-tree and octree) is better than uniform grid if the scene is irregularly distributed.
Spatial hierarchies
A A
A
Letters correspond to planes (A) Point Location by recursive search
Spatial hierarchies
A A
B B
A
Letters correspond to planes (A, B) Point Location by recursive search
Spatial hierarchies
A
D
A
B
B C
C
D
A
Letters correspond to planes (A, B, C, D) Point Location by recursive search
Variations
octree
kd tree octree bsp tree
kd-tree bsp-tree
“Hack” kd-tree building
• Split axis
R d bi l t t t
– Round-robin; largest extent
• Split location
– Middle of extent; median of geometry (balanced tree)
• Termination
– Target # of primitives, limited tree depth
• All of these techniques stink.
Building good kd-trees
• What split do we really want?
Cl Id th th t k t i h – Clever Idea: the one that makes ray tracing cheap – Write down an expression of cost and minimize it
G d t ti i ti
– Greedy cost optimization
• What is the cost of tracing a ray through a cell?
Cost(cell) = C_trav + Prob(hit L) * Cost(L) + Prob(hit R) * Cost(R)
Splitting with cost in mind
Split in the middle
To get through this part of empty space, you need to test all triangles on the right.
• Makes the L & R probabilities equal P tt ti t th L & R t
• Pays no attention to the L & R costs
Split at the median
• Makes the L & R costs equal
P tt ti t th L & R b biliti
• Pays no attention to the L & R probabilities
Cost-optimized split
Since Cost(R) is much higher, make it as small as possible
• Automatically and rapidly isolates complexity
P d l h k f t
• Produces large chunks of empty space
Building good kd-trees
• Need the probabilities
T t t b ti l t f
– Turns out to be proportional to surface area
• Need the child cell costs
– Simple triangle count works great (very rough approx.)
ll “b ” – Empty cell “boost”
Cost(cell) = C_trav + Prob(hit L) * Cost(L) + Prob(hit R) * Cost(R)
= C_trav + SA(L) * TriCount(L) + SA(R) * TriCount(R)
i th ti f th t t t t th t t C_trav is the ratio of the cost to traverse to the cost to intersect
C_trav= 1:80 in pbrt (found by experiments)
Surface area heuristic
2n splits;
must coincides with object boundary. Why?
S
aS
ba b
a
p
a S
bb
p S
S
Termination criteria
• When should we stop splitting?
B d d th li it b f t i l – Bad: depth limit, number of triangles – Good: when split does not help any more.
h h ld f
• Threshold of cost improvement
– Stretch over multiple levels
– For example, if cost does not go down after three splits in a row, terminate
• Threshold of cell size
– Absolute probability SA(node)/SA(scene) small
Basic building algorithm
1. Pick an axis, or optimize across all three 2 B ild f did li l i ( 2. Build a set of candidate split locations (cost
extrema must be at bbox vertices) 3. Sort or bin the triangles
4. Sweep to incrementally track L/R counts, cost p y , 5. Output position of minimum cost split
Running time: T ( N ) N log N 2 T ( N / 2 ) Running time:
N N
N T
N T N N N T
log
2) (
) 2 / ( 2 log )
(
• Characteristics of highly optimized tree
– very deep, very small leaves, big empty cells
Ray traversal algorithm
• Recursive inorder traversal
tmax t*
* t
t
t
tmin
*
t t tmin t* tmax
* t
* min
t t
tmaxt min max min
Intersect(L,tmin,tmax) Intersect(L,tmin,t*)Intersect(R,tmin,tmax) Intersect(R,t*,tmax)( , , )
a video for kdtree
Tree representation
8-byte (reduced from 16-byte, 20% gain)
struct KdAccelNode { interior
struct KdAccelNode { ...
union {
u_int flags; // Both float split; // Interior u int nPrims; // Leaf
leaf u_int nPrims; // Leaf
};
union {
n
union {
u_int aboveChild; // Interior MailboxPrim *onePrimitive; // Leaf MailboxPrim **primitives; // Leaf };
} }
Tree representation
1 8 23
S E M
flags
2 n
Flag: 0,1,2 (interior x, y, z) 3 (leaf)
KdTreeAccel construction
• Recursive top-down algorithm d h
8 13l ( )• max depth =
81.3log(N)If (nPrims <= maxPrims || depth==0) {
<create leaf>
}
Interior node
• Choose split axis position
M d i t – Medpoint – Medium cut
A h i ti – Area heuristic
• Create leaf if no good splits were found
• Classify primitives with respect to split
Choose split axis position
cost of no split: N
k i k t
1
) (
cost of split:
k 1
B NA
k k i A N
k k i B
t P t b P t a
t
1 1
) ( )
(
assumptions:
1. t
iiis the same for all primitives p
2. t
i: t
t= 80 : 1 (determined by experiments, main factor for the performance)
cost of split:
cost of no split:
tiN) )(
1
( e B B A A
i
t t b p N p N
t
p )
cost of split:
t i( e)(pB B pA A)s
BA B
p ( | )
B C
A
s
Ap ( | )
B C
A
Choose split axis position
Start from the axis with maximum extent, sort all edge events and process them in order all edge events and process them in order
A C
B
C
a0 b0 a1 b1 c0 c1
Choose split axis position
If there is no split along this axis, try other axes.
When all fail, create a leaf.
When all fail, create a leaf.
KdTreeAccel traversal
KdTreeAccel traversal
tmax
ToDo stack tplane
tmin
far near
KdTreeAccel traversal
ToDo stack tmax
tmin
far near
KdTreeAccel traversal
tmax
far
ToDo stack tplane
tmin
near
KdTreeAccel traversal
ToDo stack tmax
tmin
KdTreeAccel traversal
t
tmin
ToDo stack tmax
KdTreeAccel traversal
bool KdTreeAccel::Intersect
(const Ray &ray, Intersection *isect) (const Ray &ray, Intersection isect) {
if (!bounds.IntersectP(ray, &tmin, &tmax))( ( y, , )) return false;
KdAccelNode *node=&nodes[0];
while (node!=NULL) {
if (ray.maxt<tmin) break;
if (!node->IsLeaf()) <Interior>
else <Leaf>
}
} ToDo stack
} ToDo stack
(max depth)
Leaf node
1. Check whether ray intersects primitive(s) inside the node; update ray’s maxt
inside the node; update ray s maxt 2. Grab next node from ToDo queue
Interior node
1. Determine near and far (by testing which side O is)
O is)
below above below above
node+1 &(nodes[node->aboveChild])
2. Determine whether we can skip a node
node+1 &(nodes[node >aboveChild])
t tplane tplane ttmin
tmax tmin
tmax
near far near far
ttplane
Acceleration techniques Best efficiency scheme
References
• J. Goldsmith and J. Salmon, Automatic Creation of Object Hierarchies for Ray Tracing IEEE CG&A 1987 Object Hierarchies for Ray Tracing, IEEE CG&A, 1987.
• Brian Smits, Efficiency Issues for Ray Tracing, Journal of Graphics Tools, 1998. p ,
• K. Klimaszewski and T. Sederberg, Faster Ray Tracing K. Klimaszewski and T. Sederberg, Faster Ray Tracing Using Adaptive Grids, IEEE CG&A Jan/Feb 1999.
• Whang et. al., Octree-R: An Adaptive Octree for efficient ray tracing, IEEE TVCG 1(4), 1995.
• A. Glassner, Space subdivision for fast ray tracing. IEEE CG&A, 4(10), 1984