DIRECTORY List, Basics, ImagerTransformation, Real, RealFns, SF, Vector2, ImagerSample, HBrick;
~
BEGIN
Transformation:
TYPE ~ ImagerTransformation.Transformation;
BrickShape:
TYPE ~ HBrick.BrickShape;
ExtendedGCD:
TYPE ~
RECORD [p, q:
INTEGER, g:
NAT];
ExtendedEuclid:
PROC [u, v:
NAT]
RETURNS [ExtendedGCD] ~ {
Finds p, q, g such that p*u + q*v = g = GCD(u, v)
u1: INTEGER ← 1;
u2: INTEGER ← 0;
u3: NAT ← u;
v1: INTEGER ← 0;
v2: INTEGER ← 1;
v3: NAT ← v;
WHILE v3 # 0
DO
quotient: NAT ~ u3/v3;
t1: INTEGER ~ u1-(v1*quotient);
t2: INTEGER ~ u2-(v2*quotient);
t3: INTEGER ~ u3-(v3*quotient);
u1 ← v1; u2 ← v2; u3 ← v3;
v1 ← t1; v2 ← t2; v3 ← t3;
ENDLOOP;
RETURN [[p: u1, q: u2, g: u3]]
};
Sig:
PROC [a, b:
INT]
RETURNS [
INT] ~ {
RETURN [
IF a < 0
THEN -b
ELSE b]};
BrickShapeFromDeltas:
PUBLIC
PROC [s0, f0, s1, f1:
INT]
RETURNS [BrickShape] ~ {
area: CARD ~ ABS[INT[s0]*INT[f1] - INT[s1]*INT[f0]]; -- cross product to get area
gcd: ExtendedGCD ~ ExtendedEuclid[ABS[s0], ABS[s1]];
sSize: CARD ~ gcd.g;
fSize: INT ~ area/sSize;
phase: INT ← Sig[s0, f0]*INT[gcd.p] + Sig[s1, f1]*INT[gcd.q];
WHILE phase < 0 DO phase ← phase + fSize ENDLOOP;
WHILE phase >= fSize DO phase ← phase - fSize ENDLOOP;
RETURN [[sSize, fSize, phase]];
};
SqrAbs: PROC [v: VEC] RETURNS [REAL] ~ INLINE {RETURN [v.x*v.x+v.y*v.y]};
SqrLen:
PROC [u, v:
VEC]
RETURNS [
REAL] ~
INLINE {
RETURN [SqrAbs[Vector2.Sub[u, v]]]};
BrickSpec:
TYPE ~ HBrick.BrickSpec;
BrickSpecFromTransformedRectangle:
PUBLIC
PROC [w, h:
REAL, clientToDevice: Transformation, allowedRelativeError:
REAL, minLevels:
CARD]
RETURNS [BrickSpec] ~ {
v:
ARRAY [0..1]
OF
VEC ~ [
ImagerTransformation.TransformVec[clientToDevice, [w, 0]],
ImagerTransformation.TransformVec[clientToDevice, [0, h]]
];
vMul: ARRAY [0..1] OF VEC ← v;
allowSqr: REAL ~ allowedRelativeError*allowedRelativeError;
n0: INT ← 1;
n1: INT ← 1;
inc: BOOL ← FALSE;
DO
v0x: INT ← Real.Round[vMul[0].x];
v0y: INT ← Real.Round[vMul[0].y];
v1x: INT ← Real.Round[vMul[1].x];
v1y: INT ← Real.Round[vMul[1].y];
SELECT
TRUE
FROM
SqrLen[vMul[0], [v0x, v0y]] > SqrAbs[vMul[0]]*allowSqr
OR (inc
AND n0 < n1) => {
vMul[0] ← Vector2.Add[v[0], vMul[0]];
n0 ← n0 + 1;
inc ← FALSE;
};
SqrLen[vMul[1], [v1x, v1y]] > SqrAbs[vMul[1]]*allowSqr
OR inc => {
vMul[1] ← Vector2.Add[v[1], vMul[1]];
n1 ← n1 + 1;
inc ← FALSE;
};
ENDCASE => {
n0w: REAL ~ REAL[n0]*w;
n1h: REAL ~ REAL[n1*h];
brickShape: BrickShape ~ BrickShapeFromDeltas[v0x, v0y, v1x, v1y];
IF brickShape.sSize*brickShape.fSize < minLevels
THEN { inc ← TRUE }
ELSE {
m: Transformation ~ ImagerTransformation.Create[a: v0x/n0w, b: v1x/n1h, c: clientToDevice.c, d: v0y/n0w, e: v1y/n1h, f: clientToDevice.f];
m should:
map [0, 0] to [b, f],
map vector [n0*w, 0] to vector [v0x, v0y],
map vector [0, n1*h] to vector [v1x, v1y]
relativeError ← Real.SqRt[MAX[SqrLen[vMul[0], [v0x, v0y]]/SqrAbs[vMul[0]], SqrLen[vMul[1], [v1x, v1y]]/SqrAbs[vMul[1]]]];
RETURN [[brickShape, m]];
};
};
ENDLOOP;
};
Brick:
TYPE ~ HBrick.Brick;
FilterProc:
TYPE ~
PROC [x, y:
REAL]
RETURNS [
REAL];
A filter function defined on x in [-1..+1], y in [-1..+1], returning values in [0..1]
eps:
VEC ← [0.314156, 0.271828];
BrickFromDotScreen:
PUBLIC
PROC [pixelsPerDot:
REAL, degrees:
REAL, shape:
REAL, allowedRelativeError:
REAL, minLevels:
CARD, maxSample:
CARDINAL, pixelToDevice: Transformation, trc:
PROC [
REAL]
RETURNS [
REAL]]
RETURNS [Brick] ~ {
Filter: FilterProc ~ {
tx: REAL ← (x-1)*(x-1)*(x+1)*(x+1);
ty: REAL ← (y-1)*(y-1)*(y+1)*(y+1);
tx: REAL ← RealFns.CosDeg[x*180+eps.x];
ty: REAL ← RealFns.CosDeg[y*180+eps.y];
RETURN [shape*tx+(1.0-shape)*ty]
};
m: Transformation ~ ImagerTransformation.Cat[ImagerTransformation.Scale[pixelsPerDot*0.5], ImagerTransformation.Rotate[degrees], pixelToDevice];
brickSpec: BrickSpec ~ BrickSpecFromTransformedRectangle[2, 2, m, allowedRelativeError, minLevels];
brick: Brick ~ BrickFromFilter[brickSpec, Filter, maxSample, trc];
RETURN [brick];
};
LORA:
TYPE ~
LIST
OF
REF
ANY;
NodeRep:
TYPE ~
RECORD [pos:
SF.Vec, value:
REAL];
BrickFromFilter:
PUBLIC
PROC [brickSpec: BrickSpec, filter: FilterProc, maxSample:
CARDINAL, trc:
PROC [
REAL]
RETURNS [
REAL]]
RETURNS [Brick] ~ {
brickShape: BrickShape ~ brickSpec.brickShape;
Encode:
PROC [pos:
SF.Vec]
RETURNS [
CARD] ~ {
RETURN [Basics.LongMult[pos.s, brickShape.fSize]+NAT[pos.f]]
};
n: CARD ~ brickShape.fSize*brickShape.sSize;
m: CARD ~ IF n = 101 THEN 53 ELSE 101;
Equidistribute:
PROC [i:
CARD]
RETURNS [
CARD] ~ {
Assert i IN [0..n)
This should return p[i], where p is a permutation of [0..n) that scatters things around pretty well.
RETURN [(i*m+13) MOD n]
};
list: LORA ← NIL;
Compare: List.CompareProc ~ {
a: REF NodeRep ~ NARROW[ref1];
b: REF NodeRep ~ NARROW[ref2];
IF a.value < b.value THEN RETURN [less];
IF a.value > b.value THEN RETURN [greater];
IF Equidistribute[Encode[a.pos]] < Equidistribute[Encode[b.pos]] THEN RETURN [less];
RETURN [greater];
};
bitsPerSample: NAT ~ IF maxSample > 255 THEN BITS[CARDINAL] ELSE 8;
sampleMap: ImagerSample.RasterSampleMap ~ ImagerSample.NewSampleMap[box: [max: [brickShape.sSize, brickShape.fSize]], bitsPerSample: bitsPerSample];
scale: REAL ~ 1.0/REAL[n-1];
realMaxSample: REAL ~ REAL[maxSample];
ImagerSample.Clear[sampleMap];
FOR s:
NAT
IN [0..
NAT[brickShape.sSize])
DO
FOR f:
NAT
IN [0..
NAT[brickShape.fSize])
DO
p: VEC ← ImagerTransformation.InverseTransformVec[m: brickSpec.m, v: [s, f]];
WHILE p.x > 1.0 DO p.x ← p.x - 2.0 ENDLOOP;
WHILE p.x < -1.0 DO p.x ← p.x + 2.0 ENDLOOP;
WHILE p.y > 1.0 DO p.y ← p.y - 2.0 ENDLOOP;
WHILE p.y < -1.0 DO p.y ← p.y + 2.0 ENDLOOP;
list ← CONS[NEW[NodeRep ← [[s, f], filter[p.x, p.y]]], list];
ENDLOOP;
ENDLOOP;
list ← List.Sort[list: list, compareProc: Compare];
FOR i:
CARD
IN [0..n)
DO
node: REF NodeRep ~ NARROW[list.first];
rest: LIST OF REF ~ list.rest;
valReal: REAL ← i*scale;
IF trc # NIL THEN { valReal ← MIN[MAX[trc[valReal], 0.0], 1.0] };
ImagerSample.Put[map: sampleMap, index: node.pos, value: Real.Round[valReal*realMaxSample]];
list.rest ← NIL;
list ← rest;
ENDLOOP;
RETURN [[maxSample, sampleMap, brickShape.phase]]
};