(* ETH Oberon, Copyright 2001 ETH Zuerich Institut fuer Computersysteme, ETH Zentrum, CH-8092 Zuerich.
Refer to the "General ETH Oberon System Source License" contract available at: http://www.oberon.ethz.ch/ *)

MODULE GfxMatrix; (** portable *)	(* eos  *)
(** AUTHOR "eos"; PURPOSE "Affine transformations in 2D"; *)

	(*
		21.02.2000 - bugfix in Scaled: didn't recognize downscaling
		15.04.2000 - added Atan2
	*)

	IMPORT
		Streams, Math;


	CONST
		Eps = 1.0E-5;


	TYPE
		(**
			Transformation matrix
				3x2_matrices can represent any combination of affine transformations, i.e. of translation, rotation, scaling and
				shearing.

			Translate by tx, ty:
				[  1  0 ]
				[  0  1 ]
				[ tx ty ]

			Scale by sx, sy:
				[ sx   0 ]
				[  0  sy ]
				[  0   0 ]

			Rotate counter-clockwise by angle phi:
				[  cos(phi)   sin(phi) ]
				[ -sin(phi)  cos(phi) ]
				[       0             0       ]

			Shear along x_axis by factor f:
				[ 1   0 ]
				[ f    1 ]
				[ 0   0 ]
		**)

		Matrix* = ARRAY 3, 2 OF REAL;


	VAR
		Identity*: Matrix;	(** identity matrix (read_only) **)


	(**--- Matrix Computation ---**)

	(** initialize matrix with given values **)
	PROCEDURE Init* (VAR m: Matrix; m00, m01, m10, m11, m20, m21: REAL);
	BEGIN
		m[0, 0] := m00; m[0, 1] := m01;
		m[1, 0] := m10; m[1, 1] := m11;
		m[2, 0] := m20; m[2, 1] := m21
	END Init;

	(**
		Procedures Get3PointTransform, Get2PointTransform and Invert may not be able to find a solution. In that case,
		they return a singular matrix with all elements set to zero.
	**)

	(** calculate matrix that maps p0 to p1, q0 to q1, and r0 to r1 **)
	PROCEDURE Get3PointTransform* (px0, py0, px1, py1, qx0, qy0, qx1, qy1, rx0, ry0, rx1, ry1: REAL; VAR res: Matrix);
		VAR m: ARRAY 6, 7 OF REAL; i, j, k: LONGINT; max, t: REAL; v: ARRAY 6 OF REAL;
	BEGIN
		(* initialize set of linear equations for matrix coefficients *)
		m[0, 0] := px0; m[0, 1] := py0; m[0, 2] := 1.0; m[0, 3] := 0.0; m[0, 4] := 0.0; m[0, 5] := 0.0; m[0, 6] := px1;
		m[1, 0] := qx0; m[1, 1] := qy0; m[1, 2] := 1.0; m[1, 3] := 0.0; m[1, 4] := 0.0; m[1, 5] := 0.0; m[1, 6] := qx1;
		m[2, 0] := rx0; m[2, 1] := ry0; m[2, 2] := 1.0; m[2, 3] := 0.0; m[2, 4] := 0.0;  m[2, 5] := 0.0; m[2, 6] := rx1;
		m[3, 0] := 0.0; m[3, 1] := 0.0; m[3, 2] := 0.0; m[3, 3] := px0; m[3, 4] := py0; m[3, 5] := 1.0; m[3, 6] := py1;
		m[4, 0] := 0.0; m[4, 1] := 0.0; m[4, 2] := 0.0; m[4, 3] := qx0; m[4, 4] := qy0; m[4, 5] := 1.0; m[4, 6] := qy1;
		m[5, 0] := 0.0; m[5, 1] := 0.0; m[5, 2] := 0.0; m[5, 3] := rx0; m[5, 4] := ry0; m[5, 5] := 1.0; m[5, 6] := ry1;

		(* Gaussian elimination with pivoting *)
		FOR i := 0 TO 5 DO
			k := i; max := ABS(m[i, i]);
			FOR j := i+1 TO 5 DO
				IF ABS(m[j, i]) > max THEN
					k := j; max := ABS(m[j, i])
				END
			END;
			IF max < Eps THEN	(* matrix is singular *)
				Init(res, 0, 0, 0, 0, 0, 0);
				RETURN
			END;
			IF k # i THEN	(* swap rows to bring largest element up *)
				FOR j := i TO 6 DO
					t := m[i, j]; m[i, j] := m[k, j]; m[k, j] := t
				END
			END;
			FOR k := i+1 TO 5 DO
				t := m[k, i]/m[i, i];
				FOR j := i+1 TO 6 DO
					m[k, j] := m[k, j] - t * m[i, j]
				END
			END
		END;

		(* solve equations *)
		FOR i := 5 TO 0 BY -1 DO
			t := m[i, 6];
			FOR j := i+1 TO 5 DO
				t := t - v[j] * m[i, j]
			END;
			v[i] := t/m[i, i]
		END;

		Init(res, v[0], v[3], v[1], v[4], v[2], v[5])
	END Get3PointTransform;

	(** calculate matrix that maps p0 to p1 and q0 to q1 **)
	PROCEDURE Get2PointTransform* (px0, py0, px1, py1, qx0, qy0, qx1, qy1: REAL; VAR res: Matrix);
		VAR rx0, ry0, rx1, ry1: REAL;
	BEGIN
		rx0 := px0 + py0 - qy0; ry0 := py0 + qx0 - px0;
		rx1 := px1 + py1 - qy1; ry1 := py1 + qx1 - px1;
		Get3PointTransform(px0, py0, px1, py1, qx0, qy0, qx1, qy1, rx0, ry0, rx1, ry1, res)
	END Get2PointTransform;

	(** calculate inverse of matrix **)
	PROCEDURE Invert* (m: Matrix; VAR res: Matrix);
		VAR det, inv: REAL;
	BEGIN
		det := m[0, 0] * m[1, 1] - m[0, 1] * m[1, 0];
		IF ABS(det) >= Eps THEN	(* matrix can be inverted; use Cramer's rule *)
			inv := 1/det;
			res[0, 0] := +inv * m[1, 1];
			res[0, 1] := -inv * m[0, 1];
			res[1, 0] := -inv * m[1, 0];
			res[1, 1] := +inv * m[0, 0];
			res[2, 0] := +inv * (m[1, 0] * m[2, 1] - m[1, 1] * m[2, 0]);
			res[2, 1] := +inv * (m[0, 1] * m[2, 0] - m[0, 0] * m[2, 1])
		ELSE
			Init(res, 0, 0, 0, 0, 0, 0)
		END
	END Invert;


	(**--- Detection of Special Cases ---**)

	(** return determinant of matrix **)
	PROCEDURE Det* (VAR m: Matrix): REAL;
	BEGIN
		RETURN m[0, 0] * m[1, 1] - m[0, 1] * m[1, 0]
	END Det;

	(** return whether matrix is singular **)
	PROCEDURE Singular* (VAR m: Matrix): BOOLEAN;
	BEGIN
		RETURN ABS(m[0, 0] * m[1, 1] - m[0, 1] * m[1, 0]) < Eps
	END Singular;

	(** return whether matrix changes vector lengths **)
	PROCEDURE Scaled* (VAR m: Matrix): BOOLEAN;
	BEGIN
		RETURN ABS(ABS(m[0, 0] * m[1, 1] - m[0, 1] * m[1, 0]) - 1) > Eps
	END Scaled;

	(** return whether matrix includes rotation, shear, or mirror transformation **)
	PROCEDURE Rotated* (VAR m: Matrix): BOOLEAN;
	BEGIN
		RETURN (m[0, 0] < -Eps) OR (m[1, 1] < -Eps) OR (ABS(m[0, 1]) > Eps) OR (ABS(m[1, 0]) > Eps)
	END Rotated;

	(** return whether matrices should be considered equal **)
	PROCEDURE Equal* (VAR m, n: Matrix): BOOLEAN;
	BEGIN
		RETURN
			(ABS(m[0, 0] - n[0, 0]) < Eps) & (ABS(m[0, 1] - n[0, 1]) < Eps) &
			(ABS(m[1, 0] - n[1, 0]) < Eps) & (ABS(m[1, 1] - n[1, 1]) < Eps) &
			(ABS(m[2, 0] - n[2, 0]) < Eps) & (ABS(m[2, 1] - n[2, 1]) < Eps)
	END Equal;


	(**--- Matrix Concatenation ---**)

	(**
		Combinations of single transformations are evaluated from left to right. Executing Translate, Rotate or Scale
		pre-concatenates a corresponding matrix to the left of the given matrix parameter. This has the effect that
		the new transformation is applied before all previously accumulated transformations. Every transformation is
		therefore executed in the context of the coordinate system defined by the concatenation of all transformations
		to its right.
	**)

	(** translation by (dx, dy) **)
	PROCEDURE Translate* (m: Matrix; dx, dy: REAL; VAR res: Matrix);
	BEGIN
		res[0, 0] := m[0, 0]; res[0, 1] := m[0, 1];
		res[1, 0] := m[1, 0]; res[1, 1] := m[1, 1];
		res[2, 0] := m[2, 0] + dx * m[0, 0] + dy * m[1, 0];
		res[2, 1] := m[2, 1] + dx * m[0, 1] + dy * m[1, 1]
	END Translate;

	(** scale by (sx, sy) **)
	PROCEDURE Scale* (m: Matrix; sx, sy: REAL; VAR res: Matrix);
	BEGIN
		res[0, 0] := sx * m[0, 0]; res[0, 1] := sx * m[0, 1];
		res[1, 0] := sy * m[1, 0]; res[1, 1] := sy * m[1, 1];
		res[2, 0] := m[2, 0]; res[2, 1] := m[2, 1]
	END Scale;

	(** scale at (ox, oy) by (sx, sy) **)
	PROCEDURE ScaleAt* (m: Matrix; ox, oy, sx, sy: REAL; VAR res: Matrix);
		VAR tx, ty: REAL;
	BEGIN
		res[0, 0] := sx * m[0, 0]; res[0, 1] := sx * m[0, 1];
		res[1, 0] := sy * m[1, 0]; res[1, 1] := sy * m[1, 1];
		tx := ox * (1-sx); ty := oy * (1-sy);
		res[2, 0] := tx * m[0, 0] + ty * m[1, 0] + m[2, 0];
		res[2, 1] := tx * m[0, 1] + ty * m[1, 1] + m[2, 1]
	END ScaleAt;

	(** rotate counter-clockwise by angle specified by its sine and cosine **)
	PROCEDURE Rotate* (m: Matrix; sin, cos: REAL; VAR res: Matrix);
	BEGIN
		res[0, 0] := cos * m[0, 0] + sin * m[1, 0]; res[0, 1] := cos * m[0, 1] + sin * m[1, 1];
		res[1, 0] := -sin * m[0, 0] + cos * m[1, 0]; res[1, 1] := -sin * m[0, 1] + cos * m[1, 1];
		res[2, 0] := m[2, 0]; res[2, 1] := m[2, 1]
	END Rotate;

	(** rotate counter-clockwise around (ox, oy) by angle specified by its sine and cosine **)
	PROCEDURE RotateAt* (m: Matrix; ox, oy, sin, cos: REAL; VAR res: Matrix);
		VAR tx, ty: REAL;
	BEGIN
		res[0, 0] := cos * m[0, 0] + sin * m[1, 0]; res[0, 1] := cos * m[0, 1] + sin * m[1, 1];
		res[1, 0] := -sin * m[0, 0] + cos * m[1, 0]; res[1, 1] := -sin * m[0, 1] + cos * m[1, 1];
		tx := ox * (1-cos) + oy * sin; ty := oy * (1-cos) - ox * sin;
		res[2, 0] := tx * m[0, 0] + ty * m[1, 0] + m[2, 0];
		res[2, 1] := tx * m[0, 1] + ty * m[1, 1] + m[2, 1]
	END RotateAt;

	(** concatenate matrices **)
	PROCEDURE Concat* (m, n: Matrix; VAR res: Matrix);
	BEGIN
		res[0, 0] := m[0, 0] * n[0, 0] + m[0, 1] * n[1, 0];
		res[0, 1] := m[0, 0] * n[0, 1] + m[0, 1] * n[1, 1];
		res[1, 0] := m[1, 0] * n[0, 0] + m[1, 1] * n[1, 0];
		res[1, 1] := m[1, 0] * n[0, 1] + m[1, 1] * n[1, 1];
		res[2, 0] := m[2, 0] * n[0, 0] + m[2, 1] * n[1, 0] + n[2, 0];
		res[2, 1] := m[2, 0] * n[0, 1] + m[2, 1] * n[1, 1] + n[2, 1]
	END Concat;


	(**--- Arctan of Vector ---**)

	PROCEDURE Atan2* (x, y: REAL): REAL;
		VAR phi: REAL;
	BEGIN
		IF (ABS(x) < 1.0) & (ABS(y) >= ABS(x * MAX(REAL))) THEN	(* y/x would overflow *)
			IF y >= 0.0 THEN phi := Math.pi/2
			ELSE phi := -Math.pi/2
			END
		ELSIF x > 0.0 THEN	(* 1st or 4th quadrant *)
			phi := Math.arctan(y/x)
		ELSIF x < 0.0 THEN	(* 2nd or 3rd quadrant *)
			phi := Math.arctan(y/x) + Math.pi
		END;
		RETURN phi
	END Atan2;


	(**--- Matrix Application ---**)

	(** apply transformation matrix to point **)
	PROCEDURE Apply* (VAR m: Matrix; xin, yin: REAL; VAR xout, yout: REAL);
	BEGIN
		xout := xin * m[0, 0] + yin * m[1, 0] + m[2, 0];
		yout := xin * m[0, 1] + yin * m[1, 1] + m[2, 1]
	END Apply;

	(** apply transformation matrix to vector (ignoring translation) **)
	PROCEDURE ApplyToVector* (VAR m: Matrix; xin, yin: REAL; VAR xout, yout: REAL);
	BEGIN
		xout := xin * m[0, 0] + yin * m[1, 0];
		yout := xin * m[0, 1] + yin * m[1, 1]
	END ApplyToVector;

	(** apply transformation matrix to distance **)
	PROCEDURE ApplyToDist* (VAR m: Matrix; din: REAL; VAR dout: REAL);
		VAR x, y: REAL;
	BEGIN
		x := din * m[0, 0]; y := din * m[0, 1];
		IF ABS(y) < 1.0E-3 THEN dout := x
		ELSE dout := Math.sqrt(x * x + y * y)
		END
	END ApplyToDist;

	(** apply transformation matrix to axis-aligned rectangle; result is enclosing axis-aligned rectangle **)
	PROCEDURE ApplyToRect* (VAR m: Matrix; ilx, ily, irx, iuy: REAL; VAR olx, oly, orx, ouy: REAL);
		VAR l, h: REAL;
	BEGIN
		olx := m[2, 0]; orx := m[2, 0];
		l := ilx * m[0, 0]; h := irx * m[0, 0];
		IF l <= h THEN olx := olx + l; orx := orx + h ELSE olx := olx + h; orx := orx + l END;
		l := ily * m[1, 0]; h := iuy * m[1, 0];
		IF l <= h THEN olx := olx + l; orx := orx + h ELSE olx := olx + h; orx := orx + l END;
		oly := m[2, 1]; ouy := m[2, 1];
		l := ilx * m[0, 1]; h := irx * m[0, 1];
		IF l <= h THEN oly := oly + l; ouy := ouy + h ELSE oly := oly + h; ouy := ouy + l END;
		l := ily * m[1, 1]; h := iuy * m[1, 1];
		IF l <= h THEN oly := oly + l; ouy := ouy + h ELSE oly := oly + h; ouy := ouy + l END
	END ApplyToRect;

	(** apply inverse of matrix to point **)
	PROCEDURE Solve* (VAR m: Matrix; u, v: REAL; VAR x, y: REAL);
		VAR det: REAL;
	BEGIN
		det := m[0, 0] * m[1, 1] - m[0, 1] * m[1, 0];
		IF ABS(det) >= Eps THEN	(* matrix can be inverted *)
			u := u - m[2, 0]; v := v - m[2, 1];
			x := (m[1, 1] * u - m[1, 0] * v)/det;
			y := (m[0, 0] * v - m[0, 1] * u)/det
		END
	END Solve;


	(**--- Matrix I/O ---**)

	PROCEDURE Write* (VAR w: Streams.Writer; VAR m: Matrix);
		VAR i: LONGINT;
	BEGIN
		FOR i := 0 TO 2 DO
			w.RawReal(m[i, 0]); w.RawReal(m[i, 1])
		END
	END Write;

	PROCEDURE Read* (VAR r: Streams.Reader; VAR m: Matrix);
		VAR i: LONGINT;
	BEGIN
		FOR i := 0 TO 2 DO
			r.RawReal(m[i, 0]); r.RawReal(m[i, 1])
		END
	END Read;


BEGIN
	Init(Identity, 1, 0, 0, 1, 0, 0)
END GfxMatrix.