<> <> <> <> 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] ~ { <> 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 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.