ImagerTransformationImpl.mesa
Copyright © 1984 Xerox Corporation. All rights reserved.
Michael Plass, February 8, 1984 8:27:05 am PST
Doug Wyatt, October 30, 1984 7:49:36 pm PST
DIRECTORY
ImagerTransformation USING [IntRectangle, FactoredTransformation, Rectangle, Transformation, TransformationRep],
Real USING [FScale, RealException, RoundLI, SqRt],
RealFns USING [ArcTanDeg, CosDeg, SinDeg],
Vector2 USING [InlineAdd, InlineSub, VEC];
ImagerTransformationImpl: CEDAR PROGRAM
IMPORTS Real, RealFns, Vector2
EXPORTS ImagerTransformation
~ BEGIN OPEN ImagerTransformation;
Transformation: TYPE ~ ImagerTransformation.Transformation;
TransformationRep: TYPE ~ ImagerTransformation.TransformationRep;
VEC: TYPE ~ Vector2.VEC;
FromRec:
PUBLIC
PROC[r: TransformationRep]
RETURNS[Transformation] ~ {
RETURN[NEW[TransformationRep ← r]];
};
Create:
PUBLIC
PROC[a, b, c, d, e, f:
REAL]
RETURNS[Transformation] ~ {
RETURN[NEW[TransformationRep ← [a, b, c, d, e, f]]];
};
Copy:
PUBLIC
PROC[m: Transformation]
RETURNS[Transformation] ~ {
RETURN[NEW[TransformationRep ← m^]];
};
Scale:
PUBLIC
PROC[s:
REAL]
RETURNS[Transformation] ~ {
RETURN[NEW[TransformationRep ← [a: s, b: 0, c: 0, d: 0, e: s, f: 0]]];
};
Scale2:
PUBLIC
PROC[s:
VEC]
RETURNS[Transformation] ~ {
RETURN[NEW[TransformationRep ← [a: s.x, b: 0, c: 0, d: 0, e: s.y, f: 0]]];
};
Rotate:
PUBLIC
PROC[a:
REAL]
RETURNS[Transformation] ~ {
cos, sin: REAL;
SELECT a
FROM
0.0 => { cos ← 1; sin ← 0 };
90.0 => { cos ← 0; sin ← 1 };
180.0 => { cos ← -1; sin ← 0 };
-90.0, 270.0 => { cos ← 0; sin ← -1 };
ENDCASE => { cos ← RealFns.CosDeg[a]; sin ← RealFns.SinDeg[a] };
RETURN[NEW[TransformationRep ← [a: cos, b: -sin, c: 0, d: sin, e: cos, f: 0]]];
};
Translate:
PUBLIC
PROC[t:
VEC]
RETURNS[Transformation] ~ {
RETURN[NEW[TransformationRep ← [a: 1, b: 0, c: t.x, d: 0, e: 1, f: t.y]]];
};
Concat:
PUBLIC
PROC[m, n: Transformation]
RETURNS[Transformation] ~ {
a: REAL ~ m.a*n.a + m.d*n.b;
d: REAL ~ m.a*n.d + m.d*n.e;
b: REAL ~ m.b*n.a + m.e*n.b;
e: REAL ~ m.b*n.d + m.e*n.e;
c: REAL ~ m.c*n.a + m.f*n.b + n.c;
f: REAL ~ m.c*n.d + m.f*n.e + n.f;
RETURN[NEW[TransformationRep ← [a, b, c, d, e, f]]];
};
Cat:
PUBLIC
PROC[m1, m2, m3, m4, m5, m6: Transformation ←
NIL]
RETURNS[Transformation] ~ {
m: Transformation ~ NEW[TransformationRep ← m1^];
IF m2#NIL THEN PostMultiply[m, m2];
IF m3#NIL THEN PostMultiply[m, m3];
IF m4#NIL THEN PostMultiply[m, m4];
IF m5#NIL THEN PostMultiply[m, m5];
IF m6#NIL THEN PostMultiply[m, m6];
RETURN[m];
};
Invert:
PUBLIC
PROC[m: Transformation]
RETURNS[Transformation] ~ {
det: REAL ~ m.a*m.e - m.d*m.b;
RETURN[
NEW[TransformationRep ← [
a: m.e / det, d: -m.d / det,
b: -m.b / det, e: m.a / det,
c: (m.b * m.f - m.e * m.c) / det, f: (m.d * m.c - m.a * m.f) / det
]]];
};
PreMultiply:
PUBLIC
PROC[m, pre: Transformation] ~ {
rep: TransformationRep ← m^;
rep.a ← pre.a*m.a + pre.d*m.b;
rep.d ← pre.a*m.d + pre.d*m.e;
rep.b ← pre.b*m.a + pre.e*m.b;
rep.e ← pre.b*m.d + pre.e*m.e;
rep.c ← pre.c*m.a + pre.f*m.b + m.c;
rep.f ← pre.c*m.d + pre.f*m.e + m.f;
m^ ← rep;
};
PreScale:
PUBLIC
PROC[m: Transformation, s:
REAL] ~ {
IF s=1 THEN NULL
ELSE {
rep: TransformationRep ← m^;
rep.a ← s*m.a;
rep.d ← s*m.d;
rep.b ← s*m.b;
rep.e ← s*m.e;
m^ ← rep;
};
};
PreScale2:
PUBLIC
PROC[m: Transformation, s:
VEC] ~ {
rep: TransformationRep ← m^;
rep.a ← s.x*m.a;
rep.d ← s.x*m.d;
rep.b ← s.y*m.b;
rep.e ← s.y*m.e;
m^ ← rep;
};
PreRotate:
PUBLIC
PROC[m: Transformation, a:
REAL] ~ {
IF a=0 THEN NULL
ELSE {
rep: TransformationRep ← m^;
cos: REAL ~ RealFns.CosDeg[a];
sin: REAL ~ RealFns.SinDeg[a];
rep.a ← cos*m.a + sin*m.b;
rep.d ← cos*m.d + sin*m.e;
rep.b ← cos*m.b - sin*m.a;
rep.e ← cos*m.e - sin*m.d;
m^ ← rep;
};
};
PreTranslate:
PUBLIC
PROC[m: Transformation, t:
VEC] ~ {
rep: TransformationRep ← m^;
rep.c ← t.x*m.a + t.y*m.b + m.c;
rep.f ← t.x*m.d + t.y*m.e + m.f;
m^ ← rep;
};
PostMultiply:
PUBLIC
PROC[m, post: Transformation] ~ {
rep: TransformationRep ← m^;
rep.a ← m.a*post.a+m.d*post.b;
rep.b ← m.b*post.a+m.e*post.b;
rep.c ← m.c*post.a+m.f*post.b+post.c;
rep.d ← m.a*post.d+m.d*post.e;
rep.e ← m.b*post.d+m.e*post.e;
rep.f ← m.c*post.d+m.f*post.e+post.f;
m^ ← rep;
};
PostScale:
PUBLIC
PROC[m: Transformation, s:
REAL] ~ {
IF s=1 THEN NULL
ELSE {
rep: TransformationRep ← m^;
rep.a ← m.a*s;
rep.b ← m.b*s;
rep.c ← m.c*s;
rep.d ← m.d*s;
rep.e ← m.e*s;
rep.f ← m.f*s;
m^ ← rep;
};
};
PostScale2:
PUBLIC
PROC[m: Transformation, s:
VEC] ~ {
rep: TransformationRep ← m^;
rep.a ← m.a*s.x;
rep.b ← m.b*s.x;
rep.c ← m.c*s.x;
rep.d ← m.d*s.y;
rep.e ← m.e*s.y;
rep.f ← m.f*s.y;
m^ ← rep;
};
PostRotate:
PUBLIC
PROC[m: Transformation, a:
REAL] ~ {
IF a=0 THEN NULL
ELSE {
rep: TransformationRep ← m^;
cos: REAL ~ RealFns.CosDeg[a];
sin: REAL ~ RealFns.SinDeg[a];
rep.a ← m.a*cos + m.b*sin;
rep.b ← m.b*cos - m.a*sin;
rep.d ← m.d*cos + m.e*sin;
rep.e ← m.e*cos - m.d*sin;
m^ ← rep;
};
};
PostTranslate:
PUBLIC
PROC[m: Transformation, t:
VEC] ~ {
rep: TransformationRep ← m^;
rep.c ← m.c+t.x;
rep.f ← m.f+t.y;
m^ ← rep;
};
Transform:
PUBLIC
PROC[m: Transformation, v:
VEC]
RETURNS[
VEC] ~ {
IF m=NIL THEN RETURN[v]
ELSE RETURN[[v.x*m.a + v.y*m.b + m.c, v.x*m.d + v.y*m.e + m.f]];
};
TransformVec:
PUBLIC
PROC[m: Transformation, v:
VEC]
RETURNS[
VEC] ~ {
IF m=NIL THEN RETURN[v]
ELSE RETURN[[v.x*m.a + v.y*m.b, v.x*m.d + v.y*m.e]];
};
InverseTransform:
PUBLIC
PROC[m: Transformation, v:
VEC]
RETURNS[
VEC] ~ {
IF m=NIL THEN RETURN[v]
ELSE RETURN[InverseTransformVec[m, v.InlineSub[[m.c, m.f]]]];
};
InverseTransformVec:
PUBLIC
PROC[m: Transformation, v:
VEC]
RETURNS[
VEC] ~ {
IF m=NIL THEN RETURN[v]
ELSE {
det: REAL ~ m.a*m.e - m.b*m.d;
RETURN[[x: (v.x*m.e - v.y*m.b)/det, y: (v.y*m.a - v.x*m.d)/det]];
};
};
DRound:
PUBLIC
PROC[v:
VEC]
RETURNS[
VEC] ~ {
Round: PROC[r: REAL] RETURNS[REAL] ~ { RETURN[Real.RoundLI[r]] };
v.x ← Round[v.x ! Real.RealException => CONTINUE];
v.y ← Round[v.y ! Real.RealException => CONTINUE];
RETURN[v];
};
RoundXY:
PUBLIC
PROC[m: Transformation, v:
VEC]
RETURNS[
VEC] ~ {
RETURN[InverseTransform[m, DRound[Transform[m, v]]]];
};
RoundXYVec:
PUBLIC
PROC[m: Transformation, v:
VEC]
RETURNS[
VEC] ~ {
RETURN[InverseTransformVec[m, DRound[TransformVec[m, v]]]];
};
NumericalInstability: SIGNAL[real: REAL] ~ CODE;
maxRelError:
REAL ← 1.0e-5;
RestrictAbsTo:
PROC[limit:
REAL, real:
REAL]
RETURNS[
REAL] ~ {
Like a MOD, but almost symmetric around zero.
limit ← ABS[limit];
WHILE real < -limit DO real ← real + 2*limit ENDLOOP;
WHILE real >= limit DO real ← real - 2*limit ENDLOOP;
RETURN[real]
};
FactoredTransformation: TYPE ~ ImagerTransformation.FactoredTransformation;
Factor:
PUBLIC
PROC[m: Transformation]
RETURNS[FactoredTransformation] ~ {
t: TransformationRep ~ m^;
theta:
REAL ~ RestrictAbsTo[45,
0.5*RealFns.ArcTanDeg[2*(t.a*t.b+t.d*t.e), (t.a*t.a-t.b*t.b+t.d*t.d-t.e*t.e)]];
cosTheta: REAL ~ RealFns.CosDeg[theta];
sinTheta: REAL ~ RealFns.SinDeg[theta];
aa: REAL ~ cosTheta*t.a + sinTheta*t.b;
bb: REAL ~ -sinTheta*t.a + cosTheta*t.b;
dd: REAL ~ cosTheta*t.d + sinTheta*t.e;
ee: REAL ~ -sinTheta*t.d + cosTheta*t.e;
phi:
REAL ~ RestrictAbsTo[90,
IF ABS[aa]>ABS[ee] THEN RealFns.ArcTanDeg[-dd, aa] ELSE RealFns.ArcTanDeg[bb, ee]];
cosPhi: REAL ~ RealFns.CosDeg[phi];
sinPhi: REAL ~ RealFns.SinDeg[phi];
a: REAL ~ aa*cosPhi - dd*sinPhi;
b: REAL ~ bb*cosPhi - ee*sinPhi;
d: REAL ~ aa*sinPhi + dd*cosPhi;
e: REAL ~ bb*sinPhi + ee*cosPhi;
mag: REAL ~ ABS[a]+ABS[d];
IF ABS[b]+ABS[d] > maxRelError*mag THEN SIGNAL NumericalInstability[ABS[b]+ABS[d]];
RETURN[[preRotate: -theta, scale: [a, e], postRotate: -phi, translate: [t.c, t.f]]]
};
Combine:
PUBLIC
PROC[f: FactoredTransformation]
RETURNS[Transformation] ~ {
m: Transformation ~ Translate[f.translate];
PreRotate[m, f.postRotate];
PreScale2[m, f.scale];
PreRotate[m, f.preRotate];
RETURN[m];
};
Rectangle: TYPE ~ ImagerTransformation.Rectangle;
IntRectangle: TYPE ~ ImagerTransformation.IntRectangle;
TransformRectangle:
PUBLIC
PROC[m: Transformation, rect: Rectangle]
RETURNS[Rectangle] ~ {
w: VEC ← [rect.w*m.a, rect.w*m.d];
h: VEC ← [rect.h*m.b, rect.h*m.e];
p0: VEC ← Transform[m, [rect.x, rect.y]];
p1: VEC ← p0.InlineAdd[w];
p2: VEC ← p1.InlineAdd[h];
p3: VEC ← p0.InlineAdd[h];
rect.x ← MIN[p0.x, p1.x, p2.x, p3.x];
rect.y ← MIN[p0.y, p1.y, p2.y, p3.y];
rect.w ← MAX[p0.x, p1.x, p2.x, p3.x] - rect.x;
rect.h ← MAX[p0.y, p1.y, p2.y, p3.y] - rect.y;
RETURN[rect];
};
TransformIntRectangle:
PUBLIC
PROC[m: Transformation, rect: Rectangle]
RETURNS[IntRectangle] ~ {
new: Rectangle ← TransformRectangle[m, rect];
xMax: REAL ← new.x+new.w;
yMax: REAL ← new.y+new.h;
x0: INT ← Real.RoundLI[new.x];
y0: INT ← Real.RoundLI[new.y];
x1: INT ← Real.RoundLI[xMax];
y1: INT ← Real.RoundLI[yMax];
WHILE x0>new.x DO x0 ← x0-1 ENDLOOP;
WHILE y0>new.y DO y0 ← y0-1 ENDLOOP;
WHILE x1<xMax DO x1 ← x1+1 ENDLOOP;
WHILE y1<yMax DO y1 ← y1+1 ENDLOOP;
RETURN [[x0, y0, x1-x0, y1-y0]]
};
WithinQuarterPixel:
PROC[p, q:
VEC]
RETURNS[
BOOLEAN] ~ {
RETURN [ABS[p.x-q.x] <= 0.25 AND ABS[p.y-q.y] <= 0.25]
};
CloseEnough:
PUBLIC
PROC[s, t: Transformation, rangeSize:
REAL]
RETURNS[
BOOLEAN] ~ {
p00: VEC ← Transform[s, InverseTransform[t, [0, 0]]];
p10: VEC ← Transform[s, InverseTransform[t, [rangeSize, 0]]];
p01: VEC ← Transform[s, InverseTransform[t, [0, rangeSize]]];
RETURN [
WithinQuarterPixel[p00, [0, 0]] AND
WithinQuarterPixel[p10, [rangeSize, 0]] AND
WithinQuarterPixel[p10, [0, rangeSize]]
]
};
CloseToTranslation:
PUBLIC
PROC[s, t: Transformation, rangeSize:
REAL]
RETURNS[
BOOLEAN] ~ {
p10: VEC ← TransformVec[s, InverseTransformVec[t, [rangeSize, 0]]];
p01: VEC ← TransformVec[s, InverseTransformVec[t, [0, rangeSize]]];
RETURN [
WithinQuarterPixel[p10, [rangeSize, 0]] AND
WithinQuarterPixel[p10, [0, rangeSize]]
]
};
SingularValues:
PUBLIC
PROC[m: Transformation]
RETURNS[
VEC] ~ {
t1: REAL ← m.a*m.a + m.d*m.d;
t2: REAL ← m.b*m.b + m.e*m.e;
t3: REAL ← t1 - t2;
t4: REAL ← Real.FScale[m.a*m.b + m.d*m.e, 1];
s: REAL ← Real.SqRt[t3*t3+t4*t4];
x: REAL ← Real.SqRt[Real.FScale[t1+t2+s, -1]];
y: REAL ← Real.SqRt[Real.FScale[t1+t2-s, -1]];
RETURN[[x, y]]
};
fuzz: REAL ← 2.0e-6;
Different:
PROC [a, b:
REAL]
RETURNS [
BOOLEAN] ~ {
IF ABS[a-b] < fuzz THEN RETURN [FALSE];
RETURN [ABS[a - b]/MAX[ABS[a], ABS[b]] > fuzz];
};
ShouldBeSameT:
PROC [s, t: Transformation] ~ {
IF Different[s.a, t.a] THEN ERROR;
IF Different[s.b, t.b] THEN ERROR;
IF Different[s.c, t.c] THEN ERROR;
IF Different[s.d, t.d] THEN ERROR;
IF Different[s.e, t.e] THEN ERROR;
IF Different[s.f, t.f] THEN ERROR;
};
ShouldBeSamePt:
PROC [p, q:
VEC] ~ {
IF Different[p.x, q.x] THEN ERROR;
IF Different[p.y, q.y] THEN ERROR;
};
PreTrans:
PROC[t:
VEC, r: Transformation]
RETURNS[Transformation] ~ {
m: Transformation ~ NEW[TransformationRep ← r^];
PreTranslate[m, t];
RETURN[m];
};
PreRot:
PROC[a:
REAL, r: Transformation]
RETURNS[Transformation] ~ {
m: Transformation ~ NEW[TransformationRep ← r^];
PreRotate[m, a];
RETURN[m];
};
PreScl:
PROC[s:
VEC, r: Transformation]
RETURNS[Transformation] ~ {
m: Transformation ~ NEW[TransformationRep ← r^];
PreScale2[m, s];
RETURN[m];
};
SelfTest:
PROC[r: Transformation] ~ {
p: VEC ← [3.14159, -265.35];
ShouldBeSameT[PreTrans[p, r], Concat[Translate[p], r]];
ShouldBeSameT[Rotate[90], Concat[Rotate[45], Rotate[45]]];
ShouldBeSameT[Rotate[0], Concat[r, Invert[r]]];
ShouldBeSamePt[InverseTransform[r, p], Transform[Invert[r], p]];
ShouldBeSamePt[InverseTransformVec[r, p], TransformVec[Invert[r], p]];
ShouldBeSamePt[TransformVec[r, TransformVec[r, p]], TransformVec[Concat[r, r], p]];
ShouldBeSamePt[SingularValues[r], SingularValues[PreRot[33.5, r]]];
};
SelfTest[PreTrans[[12, 6], PreScl[[1, 5], Rotate[61]]]];
SelfTest[Create[1, 2, 4, 9, 3, 1]];
SelfTest[Create[3, 1, 4, 1, 5, 9]];
END.