DIRECTORY
ImagerPixelMaps			USING [BoundedWindow, Clear, Create, DeviceRectangle,
									 GetPixel, PixelMap],
ImagerPixelMapsExtras	USING [SetPixel],
Matrix3d,
Real						USING [FixI, RoundI, RoundLI, SqRt],
Filters					USING [GetFilterWidth, PixelWeight, UncoveredRatio],
FilteredTransforms;

FilteredTransformsImpl:  CEDAR PROGRAM
IMPORTS Filters, ImagerPixelMaps, ImagerPixelMapsExtras, Matrix3d, Real
EXPORTS FilteredTransforms

= BEGIN

filterWidth: INTEGER;
Direction: TYPE = {xPass, x90Pass, yPass, y90Pass};
PixelCoords: TYPE ~ ARRAY [1..2] OF INTEGER;
PixelPos: TYPE ~ ARRAY [1..2] OF REAL;
xyRectangle: TYPE ~ FilteredTransforms.xyRectangle;

PassInfoRec: TYPE ~ RECORD[ 
			area: REAL, 
			image: ImagerPixelMaps.PixelMap,
			window: xyRectangle, 
			min: INTEGER, 
			pass1, pass2: Direction, 
			matrix: Matrix3d.Matrix4by4];
PassInfo: TYPE ~ REF PassInfoRec;

EndPixelSeq: TYPE ~ REF EndPixelSeqRec;
EndPixelSeqRec: TYPE ~ RECORD[ element: SEQUENCE length: NAT OF UnPixelInfo];

UnPixelInfo: TYPE ~ REF UnPixelInfoRec;
UnPixelInfoRec: TYPE ~ RECORD[ unCovered: CoveredSeq, unPixels: UnPixelList];

CoveredSeq: TYPE ~ REF CoveredSeqRec;
CoveredSeqRec: TYPE ~ RECORD[ element: SEQUENCE length: NAT OF BOOLEAN];

UnPixelList: TYPE ~ LIST OF UnPixel;
UnPixel: TYPE ~ REF UnPixelRec;
UnPixelRec: TYPE ~ RECORD[ pixel: INTEGER, unRatio: REAL];

Image, NewImage: ImagerPixelMaps.PixelMap;
Window, NewWindow: xyRectangle;
Dir: Direction;
EndPixels: EndPixelSeq;
SecondPass: BOOLEAN;

Transform: PUBLIC PROC[
		oldImage, newImage: ImagerPixelMaps.PixelMap,
		m: Matrix3d.Matrix4by4,
		moveX, moveY: REAL _ 0.0,
		focalLength: REAL _ 1] ~ {

GetPassInfo: PRIVATE PROC[] ~ {
newWindow: xyRectangle _ GetSize[ newImage];
infoX, infoY, infoX90, infoY90: PassInfo;

m90: Matrix3d.Matrix4by4 _ Rotate90[ m, window];
m _ FixMatrix[ m, focalLength];		-- round to precision of 3 digits past decimal point
m90 _ FixMatrix[ m90, focalLength];
m90[4][1] _ m[4][1] _ moveX;		-- this isn't where they belong, but it's convenient
m90[4][2] _ m[4][2] _ moveY;

IF NOT FinalArea[ m, window] THEN RETURN[];
infoX _ InterArea[ m, window, newWindow, xPass];
infoY _ InterArea[ m, window, newWindow, yPass];
infoX90 _ InterArea[ m90, window, newWindow, x90Pass];
infoY90 _ InterArea[ m90, window, newWindow, y90Pass];

SELECT TRUE FROM
(infoX.area >= infoY.area)
AND (infoX.area >= infoX90.area)
AND (infoX.area >= infoY90.area) => {
IF infoX.area = 0 THEN RETURN[];
EndPixels _ GetEndPixels[ infoX.window.xSize, infoX.window.ySize];
passInfo _ infoX};
(infoY.area >= infoX90.area)
AND (infoY.area >= infoY90.area) => {
EndPixels _ GetEndPixels[ infoY.window.ySize, infoY.window.xSize];
passInfo _ infoY};
(infoX90.area >= infoY90.area) => {	
EndPixels _ GetEndPixels[ infoX90.window.xSize, infoX90.window.ySize];
passInfo _ infoX90};
ENDCASE  => {	
EndPixels _ GetEndPixels[ infoY90.window.ySize, infoY90.window.xSize];
passInfo _ infoY90};
};


passInfo: PassInfo;
window: xyRectangle _ GetSize[ oldImage];
GetPassInfo[];
IF passInfo = NIL THEN RETURN[];		-- can't do it; intermediate area collapses
passInfo.image _ ImagerPixelMaps.Create[ oldImage.refRep.lgBitsPerPixel, 
		[ 0, 0, passInfo.window.ySize, passInfo.window.xSize]];
ImagerPixelMaps.Clear[ passInfo.image];
filterWidth _ Filters.GetFilterWidth[];

Image _ oldImage;
NewImage _ passInfo.image;
Dir _ passInfo.pass1;
SecondPass _ FALSE;
Pass[ passInfo.matrix];

Image _ passInfo.image;
NewImage _ newImage;
Dir _ passInfo.pass2;
SecondPass _ TRUE;
Pass[ passInfo.matrix, passInfo.min];

EndPixels _ NIL;
};

Pass: PROC[ m: Matrix3d.Matrix4by4,
			min: INTEGER _ 0] ~ {


max: NAT;
Window _ GetSize[ Image];
NewWindow _ GetSize[ NewImage];

SELECT Dir FROM
xPass, y90Pass => max _ Window.ySize;
ENDCASE => max _ Window.xSize;

FOR  i:INTEGER IN [ min..max) DO
a, b, c, d, e: REAL;
[ a, b, c, d, e] _ CrunchMatrix[ m, i];
LineTransform[ i, a, b, c, d, e];
ENDLOOP;
};

LineTransform: PROC[ line: INTEGER, a, b, c, d, e: REAL] ~ {


newLow, newHigh: INTEGER;
[low: newLow, high: newHigh] _ GetLineBounds[ a, b, c, d, e];

IF NOT SecondPass THEN {			-- flag all the pixels that won't get covered
FOR j:INTEGER IN [0..newLow) DO
EndPixels[ j].unCovered[line] _ TRUE;
ENDLOOP;

FOR j:INTEGER IN (newHigh..EndPixels.length) DO
EndPixels[ j].unCovered[line] _ TRUE;
ENDLOOP;
};

FOR  j:INTEGER IN [ newLow..newHigh] DO	-- for each new pixel
PixelTransform[ line, j, a, b, c, d, e];
ENDLOOP;
};

PixelTransform: PROC[ line, pixel: NAT, a, b, c, d, e: REAL] ~ {

b1, b2: REAL;
oldLow, oldHigh: INTEGER;
total, totalUnWeight, totalWeight, unRatio: REAL _ 0.0;
colorTest: INTEGER _ 0;
newColor: [0..256) _ 0;

upperBound: INTEGER;
cut1, cut2, partCovered: BOOLEAN;
SELECT Dir FROM
xPass, y90Pass => upperBound _ Window.xSize - 1;
ENDCASE => upperBound _ Window.ySize - 1;
[b: b1, cut: cut1] _ TestBound[ 
			((pixel-filterWidth - e)*d - b), 
			(a - (pixel-filterWidth - e)*c), 
			upperBound];
[b: b2, cut: cut2] _ TestBound[ 
			((pixel+filterWidth - e)*d - b), 
			(a - (pixel+filterWidth - e)*c), 
			upperBound];
partCovered _ cut1 OR cut2;

[low: oldLow, high: oldHigh] _ GetPixelBounds[ b1, b2];

FOR  k:INTEGER IN [ oldLow..oldHigh] DO

weight, unWeight: REAL _ 0.0;
pixelValue: INTEGER;
lPos: REAL _ (a* (k-1) + b) / (c*(k-1) + d) + e;
hPos: REAL _ (a* (k+1) + b) / (c*(k+1) + d) + e;
newWidth: REAL _ (hPos - lPos) * filterWidth / 2;
newPos: REAL _ lPos + newWidth;

IF partCovered THEN
IF k = oldLow THEN
unRatio _ unRatio + Filters.UncoveredRatio[ 
			newPos, newWidth, 
			pixel, filterWidth, 
			FALSE]
ELSE IF k = oldHigh THEN
unRatio _ unRatio + Filters.UncoveredRatio[ 
			newPos, newWidth, 
			pixel, filterWidth, 
			TRUE];

weight _ Filters.PixelWeight[ 
		newPos, newWidth, 
		pixel, filterWidth];
		

pixelValue _ GetPixel[ line, k, Image];

IF pixelValue = -1 THEN {				-- pixel is outside Image
unWeight _ weight;
weight _ 0
}
ELSE IF SecondPass THEN				-- check for uncovered pixels from interImage
[weight: weight, unWeight: unWeight] _ TestWeights[ line, k, weight];

total _ total + weight * pixelValue;
totalWeight _ totalWeight + weight;	-- COVERED weight for this pixel
totalUnWeight _ totalUnWeight + unWeight;	-- UNCOVERED weight for this pixel
ENDLOOP;


IF totalWeight = 0 THEN {								-- totally uncovered
IF NOT SecondPass THEN EndPixels[ pixel].unCovered[line] _ TRUE}
ELSE {
IF totalUnWeight = 0 AND unRatio = 0 THEN	-- totally covered
colorTest _ Real.RoundI[total / totalWeight]
ELSE 													-- partially covered
IF SecondPass THEN {
background: INTEGER _ GetPixel[ line, pixel, NewImage];
colorTest _ Real.RoundI[ 
	((background*totalUnWeight + total) / 
	(totalWeight + totalUnWeight)) * (1-unRatio) + 
	background*unRatio]
}
ELSE {										-- add to list of partially covered pixels
SetUnPixel[ line, pixel,
	(totalUnWeight / 
	(totalWeight + totalUnWeight))*(1-unRatio) + 
	unRatio];
colorTest _ Real.RoundI[total / totalWeight];
};

IF colorTest > 0 THEN 		-- keep color within valid range
IF colorTest < 256 THEN newColor _ colorTest
ELSE newColor _ 255
ELSE newColor _ 0;
SetPixel[ line, pixel, newColor];
};
};


TransformPixel: PROC[ p: PixelCoords, m: Matrix3d.Matrix4by4]
	RETURNS[ point: PixelPos, boundsFault: BOOLEAN _ FALSE] ~ {
newZ: REAL _ m[3][1]*p[1] + m[3][2]*p[2] + m[3][4];

IF newZ < m[1][1]/10 OR newZ < m[2][2]/10 THEN
RETURN[ [0, 0], TRUE];

point[1] _ (m[1][1]*p[1] + m[1][2]*p[2] + m[1][4]) / newZ + m[4][1];
point[2] _ (m[2][1]*p[1] + m[2][2]*p[2] + m[2][4]) / newZ + m[4][2];
};

FinalArea: PROC[ m: Matrix3d.Matrix4by4, window: xyRectangle]
	RETURNS[ BOOLEAN] ~ {
corner1, corner2, corner3, corner4: PixelPos;
area, longestSideLength: REAL;
boundsFault: BOOLEAN _ FALSE;

[ point: corner1, boundsFault: boundsFault] _ TransformPixel[ [0, 0], m];
IF boundsFault THEN RETURN[ FALSE];
[point: corner2, boundsFault: boundsFault] _ TransformPixel[ [0, window.ySize-1], m];
IF boundsFault THEN RETURN[ FALSE];
[point: corner3, boundsFault: boundsFault] _ TransformPixel[ [window.xSize-1, 0], m];
IF boundsFault THEN RETURN[ FALSE];
[point: corner4, boundsFault: boundsFault] _ TransformPixel[ [window.xSize-1, window.ySize-1], m];
IF boundsFault THEN RETURN[ FALSE];

[area: area, longestSideLength: longestSideLength] _ GetArea[ corner1, corner2, corner3, corner4];

IF area < longestSideLength/512 THEN RETURN[ FALSE] ELSE RETURN[ TRUE];
};

InterArea: PROC[ m: Matrix3d.Matrix4by4, 
			window, newWindow: xyRectangle, 
			dir: Direction]
	RETURNS[ info: PassInfo] ~ TRUSTED {

InterWindow: PROC[] ~ TRUSTED {
minA, minB, maxA, maxB: INTEGER;

[ low: minA, high: maxA] _ GetPixelBounds[ corner1[1], corner2[1]];
[ low: minB, high: maxB] _ GetPixelBounds[ corner3[1], corner4[1]];

SELECT dir FROM
xPass => {
size: INTEGER;
info.min _ MIN[ newWindow.xSize, MAX[ 0, MIN[ minA, minB]]];
size _ MIN[ newWindow.xSize, MAX[ 0, MAX[ maxA+1, maxB+1]]];
info.window _ [ 0, 0, size, window.ySize]};
x90Pass => {
size: INTEGER;
info.min _ MIN[ newWindow.xSize, MAX[ 0, MIN[ minA, minB]]];
size _ MIN[ newWindow.xSize, MAX[ 0, MAX[ maxA+1, maxB+1]]];
info.window _ [ 0, 0, size, window.xSize]};
yPass => {
size: INTEGER;
info.min _ MIN[ newWindow.ySize, MAX[ 0, MIN[ minA, minB]]];
size _ MIN[ newWindow.ySize, MAX[ 0, MAX[ maxA+1, maxB+1]]];
info.window _ [ 0, 0, window.xSize, size]};
ENDCASE => {
size: INTEGER;
info.min _ MIN[ newWindow.ySize, MAX[ 0, MIN[ minA, minB]]];
size _ MIN[ newWindow.ySize, MAX[ 0, MAX[ maxA+1, maxB+1]]];
info.window _ [ 0, 0, window.ySize, size]};
};

corner1, corner2, corner3, corner4: PixelPos;
boundsFault: BOOLEAN _ FALSE;
a, b: INTEGER;

matrix: Matrix3d.Matrix4by4 _ Matrix3d.Identity[];
SELECT dir FROM
xPass, x90Pass => {
matrix[1] _ m[1];
matrix[3] _ m[3];
matrix[4] _ m[4];
};
yPass, y90Pass => {
matrix[1] _ [m[2][2], m[2][1], m[2][3], m[2][4]];
matrix[3] _ [m[3][2], m[3][1], m[3][3], m[3][4]];
matrix[4] _ [m[4][2], m[4][1], m[4][3], m[4][4]];
};
ENDCASE;

info _ NEW[ PassInfoRec];

SELECT dir FROM
xPass => { a _ window.xSize-1;	b _ window.ySize-1};
x90Pass => { a _ window.ySize-1;	b _ window.xSize-1};
yPass => { a _ window.ySize-1;	b _ window.xSize-1};
ENDCASE => { a _ window.xSize-1;	b _ window.ySize-1};

[ point: corner1, boundsFault: boundsFault] _ TransformPixel[ [0, 0], matrix];
IF boundsFault THEN RETURN;
[ point: corner2, boundsFault: boundsFault] _ TransformPixel[ [0, b], matrix];
IF boundsFault THEN RETURN;
[ point: corner3, boundsFault: boundsFault] _ TransformPixel[ [a, 0], matrix];
IF boundsFault THEN RETURN;
[ point: corner4, boundsFault: boundsFault] _ TransformPixel[ [a, b], matrix];
IF boundsFault THEN RETURN;

[area: info.area] _ GetArea[ corner1, corner2, corner3, corner4];
InterWindow[];
info.matrix _ m;
info.pass1 _ dir;
SELECT dir FROM
xPass, x90Pass => info.pass2 _ yPass;
ENDCASE => info.pass2 _ xPass;
};

GetArea: PROC[ corner1, corner2, corner3, corner4: PixelPos]
	RETURNS[ area, longestSideLength: REAL] ~ {
side1, side2, side3, side4: PixelPos;		-- these are actually vector coordinates

CrossProduct: PROC[ a, b: PixelPos]
	RETURNS[ aXb: REAL] ~ {
RETURN[ a[1]*b[2] - a[2]*b[1]];
};

side1[1] _ corner2[1] - corner1[1];		side1[2] _ corner2[2] - corner1[2];
side2[1] _ corner3[1] - corner1[1];		side2[2] _ corner3[2] - corner1[2];
side3[1] _ corner4[1] - corner2[1];		side3[2] _ corner4[2] - corner2[2];
side4[1] _ corner4[1] - corner3[1];		side4[2] _ corner4[2] - corner3[2];

area _ ABS[ CrossProduct[ side1, side2] + 
				CrossProduct[ side3, side2]/2 +
				CrossProduct[ side1, side4]/2];

longestSideLength _ MAX[ MAX[ MAX[ 
			(side1[1]*side1[1] + side1[2]*side1[2]),
			(side2[1]*side2[1] + side2[2]*side2[2])],
			(side3[1]*side3[1] + side3[2]*side3[2])],
			(side4[1]*side4[1] + side4[2]*side4[2])];
longestSideLength _ Real.SqRt[ longestSideLength];
};


Rotate90: PROC[ m: Matrix3d.Matrix4by4, window: xyRectangle] 
	RETURNS[ m90: Matrix3d.Matrix4by4] ~ TRUSTED {
m90 _ Matrix3d.Identity[];
m90 _ Matrix3d.Translate[ m90, 0, -(window.xSize-1), 0];
m90 _ Matrix3d.RotateAboutZAxis[ m90, 90];				-- rotate window 90 degrees
m90 _ Matrix3d.MatMult[ m, m90];							-- transform matrix
};

FixMatrix: PROC[ m: Matrix3d.Matrix4by4, focalLength: REAL] 
	RETURNS[ fixedM: Matrix3d.Matrix4by4] ~ TRUSTED {
FOR i: INTEGER IN [1..4] DO
FOR j: INTEGER IN [1..4] DO
fixedM[i][j] _ Real.RoundLI[ m[i][j]*1000]/1000.0;	-- truncate to e-3 resolution
ENDLOOP;
ENDLOOP;
FOR i: INTEGER IN [1..4] DO
fixedM[3][i] _ -fixedM[3][i];
ENDLOOP;

fixedM _ Matrix3d.Scale[ fixedM, focalLength, focalLength, 1];	-- perspective transformation
fixedM _ Matrix3d.Translate[ fixedM, 0, 0, focalLength];
};

CrunchMatrix: PROC[ m: Matrix3d.Matrix4by4, line: INTEGER] 
	RETURNS[ a, b, c, d, e: REAL] ~ {


SELECT Dir FROM
xPass, x90Pass => {
IF SecondPass THEN {
linePos: REAL _ line - m[4][2];
den: REAL _ (linePos * m[3][2] - m[2][2]);
IF den = 0 THEN a _ b _ c _ d _ 0		-- result is infinitely large
ELSE {
a _ m[1][1] + m[1][2] * (m[2][1] - linePos * m[3][1]) / den;
b _ m[1][4] + m[1][2] * (m[2][4] - linePos * m[3][4]) / den;
c _ m[3][1] + m[3][2] * (m[2][1] - linePos * m[3][1]) / den;
d _ m[3][4] + m[3][2] * (m[2][4] - linePos * m[3][4]) / den}
}
ELSE {
a _ m[1][1];
b _ m[1][2] * line + m[1][4];
c _ m[3][1];
d _ m[3][2] * line + m[3][4]};
e _ m[4][1]};

yPass, y90Pass => {
IF SecondPass THEN {
linePos: REAL _ line - m[4][1];
den: REAL _ (linePos * m[3][1] - m[1][1]);
IF den = 0 THEN a _ b _ c _ d _ e _ 0		-- result is infinitely large
ELSE {
a _ m[2][2] + m[2][1] * (m[1][2] - linePos * m[3][2]) / den;
b _ m[2][4] + m[2][1] * (m[1][4] - linePos * m[3][4]) / den;
c _ m[3][2] + m[3][1] * (m[1][2] - linePos * m[3][2]) / den;
d _ m[3][4] + m[3][1] * (m[1][4] - linePos * m[3][4]) / den}
}
ELSE {
a _ m[2][2];
b _ m[2][1] * line + m[2][4];
c _ m[3][2];
d _ m[3][1] * line + m[3][4]};
e _ m[4][2]};

ENDCASE => a _ b _ c _ d _ e _ 0
};


GetEndPixels: PROC[ seqLen, lineLen: NAT] RETURNS[ seq: EndPixelSeq] ~ {
seq _ NEW[ EndPixelSeqRec[ seqLen]];
FOR i: INTEGER IN [0..seqLen) DO
seq[i] _ NEW[ UnPixelInfoRec];
seq[i].unCovered _ NEW[ CoveredSeqRec[ lineLen]];
ENDLOOP;
};

SetUnPixel: PROC[ line, pixel: INTEGER, unRatio: REAL] ~ {
thisUnPixel: UnPixel _ NEW[ UnPixelRec];
thisUnPixel.pixel _ line;
thisUnPixel.unRatio _ unRatio;
EndPixels[ pixel].unPixels _ CONS[ thisUnPixel, EndPixels[ pixel].unPixels];
};

TestWeights: PROC[ line, pixel: INTEGER, w: REAL] 
	RETURNS[ weight, unWeight: REAL _ 0.0] ~ {
weight _ w;
IF EndPixels[ line].unCovered[ pixel] = TRUE THEN {
unWeight _ w;
weight _ 0}
ELSE {					-- check for partially covered pixels
unPixels: UnPixelList _ EndPixels[ line].unPixels;
unPixel: UnPixel;
WHILE unPixels # NIL DO
unPixel _ unPixels.first;
IF unPixel.pixel = pixel THEN {	-- got one
unWeight _ w * unPixel.unRatio;
weight _ w * (1 - unPixel.unRatio)};
unPixels _ unPixels.rest;
ENDLOOP
};
};


GetLineBounds: PROC[ a, b, c, d, e: REAL]
	RETURNS[ low, high: INTEGER] ~ {


b1, b2, den1, den2: REAL;
lastPixel, lastNewPixel: INTEGER;

SELECT Dir FROM
xPass => { lastPixel _ Window.xSize-1; lastNewPixel _ NewWindow.xSize-1}; 
x90Pass => { lastPixel _ Window.ySize-1; lastNewPixel _ NewWindow.xSize-1};
yPass => { lastPixel _ Window.ySize-1; lastNewPixel _ NewWindow.ySize-1};
ENDCASE => { lastPixel _ Window.xSize-1; lastNewPixel _ NewWindow.ySize-1}; 

den1 _ d;
den2 _ c*lastPixel+d;
IF den1 < 0.1 OR den2 < 0.1 THEN 	-- z is out of range
RETURN[ 0, 0];

b1 _ b/den1 + e;
b2 _ (a*lastPixel+b)/den2 + e;

[low: low, high: high] _ GetPixelBounds[ b1, b2];

low _ MIN[ lastNewPixel, MAX[ 0, low]];
high _ MIN[ lastNewPixel, MAX[ 0, high]];
};

TestBound: PROC[ num, den: REAL _ 0.0, upperBound: INTEGER] 
	RETURNS[ b: REAL, cut: BOOLEAN] ~ {


IF upperBound = 0 THEN RETURN[ 0, TRUE];
IF ABS[ den] <= ABS[ num/upperBound] THEN {	-- out of bounds
IF (den >= 0 AND num > 0) OR (den < 0 AND num < 0) THEN
b _ upperBound
ELSE b _ 0;				-- num/den < 0 or num/den = 0/0 = 0
cut _ TRUE}
ELSE {
b _ num / den;
IF b < 0 THEN {
b _ 0;
cut _ TRUE}
ELSE cut _ FALSE};
};

GetPixelBounds: PROC[ b1, b2: REAL]
	RETURNS[ low, high: INTEGER] ~ {

rem: REAL;
IF b1 < b2 THEN {
[i:low] _ SplitToIandR[ b1];
[i:high, rem: rem] _ SplitToIandR[ b2 + 1]}
ELSE IF b1 = b2 THEN {					-- boundary falls within a single pixel
[i:low] _ SplitToIandR[ b1];
[i:high, rem: rem] _ SplitToIandR[ b1 + 1]}
ELSE {
[i:low] _ SplitToIandR[ b2];
[i:high, rem: rem] _ SplitToIandR[ b1 + 1]};
IF rem = 0.0 THEN high _ high - 1;

low _ low - (filterWidth-1);
high _ high + (filterWidth-1);
};


GetPixel: PRIVATE PROC[ line, pixel: INTEGER, image: ImagerPixelMaps.PixelMap]
RETURNS[ sample: INTEGER _ -1] ~ {

pt, ptL: PixelCoords;
window: xyRectangle _ GetSize[ image];
SELECT Dir FROM
xPass => {
pt[1] _ pixel;
pt[2] _ line};
x90Pass => {
pt[1] _ (window.xSize - 1) - line;		-- turn input -90 degrees
pt[2] _ pixel};								-- (x,y)
yPass => {
pt[1] _ line;
pt[2] _ pixel};
ENDCASE => {									-- y90Pass
pt[1] _ (window.xSize - 1) - pixel;
pt[2] _ line};

IF 0 <= pt[1] AND pt[1] < window.xSize 			-- ensure that pixel is inside window
AND 0 <= pt[2] AND pt[2] < window.ySize THEN {

ptL _ ConvertRtoLCoords[ window, pt];
sample _ ImagerPixelMaps.GetPixel[ image, ptL[1], ptL[2]]
}
};

SetPixel: PRIVATE PROC[
			line, newPixel: NAT,
			newColor: [0..256)] ~ {

pt, ptL: PixelCoords;
SELECT Dir FROM
xPass, x90Pass => {
pt[1] _ newPixel;
pt[2] _ line};			-- (s, f)
yPass, y90Pass => {
pt[1] _ line;
pt[2] _ newPixel};
ENDCASE => RETURN;		-- no such thing

ptL _ ConvertRtoLCoords[ NewWindow, pt];
ImagerPixelMapsExtras.SetPixel[ NewImage, ptL[1], ptL[2], newColor]
};

GetSize: PROC[ image: ImagerPixelMaps.PixelMap] RETURNS[ window: xyRectangle] ~ {
imageBounds: ImagerPixelMaps.DeviceRectangle _ image.BoundedWindow[];
window.xMin _ imageBounds.fMin;
window.xSize _ imageBounds.fSize;
window.yMin _ imageBounds.sMin;
window.ySize _ imageBounds.sSize;
};

ConvertRtoLCoords: PROC[ window: xyRectangle, ptR: PixelCoords] RETURNS[ ptL: PixelCoords] ~ {


ptL[2] _ window.xMin + ptR[1];
ptL[1] _ (window.yMin + window.ySize - 1) - ptR[2];
};

SplitToIandR: PUBLIC PROC[ r: REAL] RETURNS[ i: INTEGER, rem: REAL] ~ {
IF r >= 0.0 THEN i _ Real.FixI[ r]
ELSE {
i _ Real.FixI[ r];					-- truncates towards zero
IF i # r THEN i _ i - 1};		-- I want to truncate towards negative infinity
rem _  r - i;
};

END.

	\FilteredTransformsImpl.mesa
Last Edited by: Hiller, September 7, 1984 0:02:00 am PDT
Copyright (C) 1984 by Xerox Corporation.  All rights reserved.


Main Program Loop
Beginning of Transform
first pass
second pass
clean up global storage before exiting
Image _ NIL;
NewImage _ NIL;
loop through all lines in the image and transform each
transform a line


transform a pixel
back-transform to find the old pixels that will fall under this new pixel
find distance between the last pixel and the next.  model this pixel as a symmetrical one that is centered between the last and the next.
find values for this pixel and add to totals
average the values and blend with background to get final pixel value
Area Calculation
runs an individual pixel through the entire matrix transform (including final translation)
if image scaling is more than 10x in either direction, don't do it (you can't get a good picture anyhow at this point)
m[4] contains the translations that should occur after the matrix transformation has been performed.  This is not where they really belong, but the space isn't being used anyhow, and it seemed easier than passing around the new variables.
if this is true, then the area is so small it won't change the value of any of the pixels it overlaps
Beginning of InterArea

transform first but not second variable
convert corner points into vectors
| (side1 X side2) + (side3 X side2) / 2 + (side1 X side4) / 2 | is area of quadrilateral
Matrix Operations
this procedure returns the relevant scaling and skew coefficients for the given line.  It also does a reverse-transform to determine the values using the already transformed value from the first pass if necessary.  c and d returned are the z scale and skew coefficients, for doing perspective.  e is the translation to give the entire line after it has been transformed, to move it into its final position in the final image.
Edge Blending
Bounds Checking
finds lowest and highest value of newPixel for this line, using the size of Window as the old line length.  The values are constrained by the size of NewWindow.
in the '90' passes, the input is rotated so x bounds become y bounds and vice versa
get real values of edges of image
assumes lower bound is 0
convert real boundaries to integers to find pixels which overlap the boundary
Type Conversion
returns location of the pixel in screen coordinates (s, f)

�Êœ˜š
œ™J™8Jš
œ>™>—J™š
Ïk	˜	Jšœ�œO˜fJšœ�œ˜'Jš
œ	˜	Jšœ
�œ˜.Jšœ�œ/˜@Jš
œ˜J˜—šœ�œ
�˜&Jš�œ@˜GJš�œ˜—J˜Jšœ�˜J˜Jšœ
�œ
˜Jšœ�œ$˜3Jš	œ
�œ�œ�œ
�œ
˜,Jš	œ
�œ�œ�œ
�œ
˜&Jšœ
�œ"˜3J˜Jš	œ
�œ�œ�œH�œ@˜¹
Jšœ
�œ�œ
˜!J˜Jšœ
�œ�œ˜'Jšœ�œ�œ�œ	�œ
�œ˜MJ˜Jšœ
�œ�œ˜'Jšœ�œ�œ0˜MJ˜Jšœ�œ�œ˜%Jš
œ�œ�œ�œ	�œ
�œ
�œ˜HJ˜Jšœ
�œ�œ
�œ	˜$Jšœ	�œ�œ˜Jš	œ�œ�œ	�œ�œ˜:J˜Jš
œ*˜*Jš
œ˜Jš
œ˜Jš
œ˜Jšœ�œ
˜J™š
œ™J˜šÏn	œ�œ
�œ
˜JšœZ�œ�œ	˜‚
J˜šžœ�œ
�œ˜Jš
œ,˜,Jš
œ)˜)J˜Jš
œ0˜0Jšœ!Ïc4˜UJ˜#JšœŸ4˜RJš
œ˜J˜Jš�œ
�œ�œ
�œ˜+I
aš
œ0˜0Jš
œ0˜0Kš
œ6˜6Kš
œ6˜6J˜š�œ
�œ
�˜Jš
œ˜Jš�œ˜ š�œ"˜%Jš�œ�œ
�œ˜ Jš
œB˜BJš
œ˜—Jš
œ˜š�œ"˜%Jš
œB˜BJš
œ˜—š
œ$˜$Jš
œF˜FJš
œ˜—š�œ˜Jš
œF˜FJš
œ˜——J˜J˜—Kš
œ™K˜Kš
œ˜Kš
œ)˜)J˜Jš	�œ�œ
�œ
�œŸ+˜MJš
œƒ
˜ƒ
Jš
œ'˜'J˜'J˜Jš
œ
™
Jš
œ˜Jš
œ˜Jš
œ˜Jšœ
�œ
˜Jš
œ˜J˜Jš
œ™Jš
œ˜Jš
œ˜Jš
œ˜Jšœ
�œ
˜Jš
œ%˜%J˜Jš
œ&™&Jš
œ™Jš
œ™Jšœ�œ
˜J˜—J˜šžœ�œ"�œ	˜<J˜Kš
œ6™6K˜Kšœ�œ
˜	Kš
œ˜Jš
œ˜J˜š�œ�˜Jš
œ%˜%Kš�œ˜—K˜š�œ�œ
�œ
�˜ Jšœ�œ
˜Jš
œ'˜'Jš
œ!˜!Kš�œ
˜—K˜—J˜šž
œ�œ�œ�œ˜<K˜Kš
œ™K˜Kšœ�œ
˜J™š
œ=˜=K˜—š�œ
�œ�œŸ-˜Hš�œ�œ
�œ
�˜Kšœ �œ
˜%Kš�œ
˜K˜—š�œ�œ
�œ�˜/Kšœ �œ
˜%Kš�œ
˜—K˜K˜—š	�œ�œ
�œ�œ
Ÿ˜=Kš
œ(˜(Jš�œ
˜—J˜—J˜šžœ�œ�œ�œ˜@K™Kš
œ™K˜Kšœ�œ
˜
Kšœ�œ
˜Kšœ,�œ˜7Kšœ�œ˜K˜K˜Kš
œI™IKšœ�œ
˜Kšœ�œ
˜!š�œ�˜Kš
œ0˜0Kš�œ"˜)—Kš
œz˜zKš
œz˜zšœ�œ˜K˜—Kš
œ7˜7K˜Kš�œ�œ
�œ�˜'˜Kš
œ‰
™‰
Kšœ�œ˜Kšœ�œ
˜Kšœ�œ&˜0Kšœ�œ&˜0Kšœ
�œ#˜1Kšœ�œ˜K˜š�œ
�˜š�œ�˜Jšœ^�œ
˜d—š�œ
�œ
�˜Jšœ^�œ˜dK˜——Kš
œM˜MKš
œ,™,K˜Kš
œ'˜'K˜š�œ�œŸ˜6Kš
œ˜Kš
œ
˜
K˜
—š�œ
�œ�œŸ-˜HKš
œE˜EK˜—Kš
œ$˜$Kšœ$Ÿ ˜DKšœ*Ÿ"˜LKš�œ
˜—K˜Kš
œE™EK˜š�œ�œ
Ÿ˜5Kš�œ
�œ�œ%�œ
˜@—š�œ˜š�œ�œ
�œ
Ÿ˜<Kš
œ,˜,—š�œŸ˜&š�œ�œ˜Kšœ�œ$˜7Kš
œ‡
˜‡
Kš
œ
˜
—š�œŸ*˜:Kš
œe˜eKš
œ-˜-K˜K˜——š�œ�œŸ ˜8Jš�œ�œ˜,Jš�œ˜—Jš�œ˜Jš
œ!˜!J˜—J˜——J˜š
œ™J˜Kš
œZ™Zšžœ�œ)˜=Kšœ
�œ �œ�œ˜<Kšœ�œ)˜3K˜Kš
œv™vš�œ�œ�˜.Kš�œ
�œ˜—K˜Kš
ϔ
™î
Kš
œD˜DKš
œD˜DK˜—K˜šž	œ�œ.˜=Kšœ
�œ�œ˜K˜-Kšœ�œ
˜Kšœ
�œ�œ
˜K˜Kš
œI˜IKš�œ
�œ
�œ�œ˜#Kš
œU˜UKš�œ
�œ
�œ�œ˜#Kš
œU˜UKš�œ
�œ
�œ�œ˜#Kš
œb˜bKš�œ
�œ
�œ�œ˜#K˜Kš
œb˜bK˜Kš
œe™eKš�œ�œ
�œ�œ�œ
�œ�œ˜GK˜—K˜šž	œ�œQ˜`Kšœ
�œ�œ˜%K˜šžœ�œ�œ˜Kšœ�œ
˜ J˜Jš
œC˜CJš
œC˜CJ˜š�œ�˜š
œ
˜
Jšœ�œ
˜Jšœ�œ�œ�œ˜<Jšœ�œ�œ�œ˜<Jš
œ+˜+—š
œ˜Jšœ�œ
˜Jšœ�œ�œ�œ˜<Jšœ�œ�œ�œ˜<Jš
œ+˜+—š
œ
˜
Jšœ�œ
˜Jšœ�œ�œ�œ˜<Jšœ�œ�œ�œ˜<Jš
œ+˜+—š�œ˜Jšœ�œ
˜Jšœ�œ�œ�œ˜<Jšœ�œ�œ�œ˜<Jš
œ+˜+——K˜K˜—Kš
œ™K™Kš
œ-˜-Kšœ
�œ�œ
˜Jšœ�œ
˜K˜Kš
œ'™'Kš
œ2˜2š�œ�˜š
œ˜Kš
œ˜Kš
œ˜Kš
œ˜K˜—š
œ˜Kš
œ1˜1Kš
œ1˜1Kš
œ1˜1K˜—Kš�œ
˜K˜—Kšœ�œ˜K˜š�œ�˜Jš
œ3˜3Jš
œ5˜5Jš
œ3˜3Jš�œ.˜5—K˜Kš
œN˜NKš�œ
�œ
�œ
˜Kš
œN˜NKš�œ
�œ
�œ
˜Kš
œN˜NKš�œ
�œ
�œ
˜Kš
œN˜NKš�œ
�œ
�œ
˜K˜Kš
œA˜AKš
œ˜Kš
œ˜Kš
œ˜š�œ�˜Jš
œ%˜%Jš�œ˜—K˜—K˜šžœ�œ/˜<Jšœ
�œ�œ˜,Kšœ'Ÿ(˜OK˜šžœ�œ˜#Kšœ
�œ�œ˜Kš�œ˜K˜—K˜Kš
œ"™"Kš
œH˜HKš
œH˜HKš
œH˜HKš
œH˜HK˜Kš
œX™XKšœ�œh˜rK˜Kšœ�œ�œ�œµ
˜Ö
Kš
œ2˜2K˜——J˜š
œ™J˜šžœ�œ/˜=Jšœ
�œ�œ˜/Jš
œ˜Jš
œ8˜8Jšœ.Ÿ˜IJšœ'Ÿ˜:K˜—K˜šž	œ�œ'�œ˜<Kšœ
�œ!�œ˜2š�œ�œ
�œ�˜š�œ�œ
�œ�˜Kšœ3Ÿ˜PKš�œ
˜—Kš�œ
˜—š�œ�œ
�œ�˜Kš
œ˜Kš�œ
˜—K˜Kšœ?Ÿ˜\Kš
œ8˜8K˜—K˜šžœ�œ �œ˜;Jšœ
�œ�œ˜"J˜Jš
œ©™©J˜š�œ�˜˜š�œ�œ˜Jšœ	�œ˜Jšœ�œ!˜*Jš�œ	�œŸ˜@š�œ˜Jš
œ<˜<Jš
œ<˜<Jš
œ<˜<Jš
œ<˜<—J˜
—š�œ˜J˜Jš
œ˜J˜Jš
œ˜—J˜
J˜—˜š�œ�œ˜Jšœ	�œ˜Jšœ�œ!˜*Jš�œ	�œŸ˜Dš�œ˜Jš
œ<˜<Jš
œ<˜<Jš
œ<˜<Jš
œ<˜<—J˜
—š�œ˜J˜Jš
œ˜J˜Jš
œ˜—J˜
J˜—Jš�œ˜ —J˜——J˜š
œ
™
J˜šžœ�œ�œ�œ˜HKšœ�œ˜$š�œ�œ
�œ
�˜ Kšœ	�œ˜Kšœ�œ˜1Kš�œ
˜—K˜—J˜šž
œ�œ�œ�œ˜:Kšœ�œ˜(K˜Kš
œ˜Kšœ�œ+˜LK˜—K˜šžœ�œ�œ�œ˜2Kšœ
�œ�œ˜+Kš
œ˜š�œ&�œ
�œ˜3Kš
œ
˜
Kš
œ˜—š�œŸ%˜0Kš
œ2˜2Kš
œ˜š�œ�œ
�˜Kš
œ˜š�œ�œŸ
˜*Kš
œ˜Kš
œ$˜$—Kš
œ˜Kš
�˜—Kš
œ˜—K˜——J˜š
œ™J˜šž
œ�œ�œ
˜)Jšœ
�œ
�œ˜!J˜Jš
œ 
™ 
J˜Jšœ�œ
˜Jšœ�œ
˜!J˜Jš
œS™Sš�œ�˜Jš
œJ˜JJš
œK˜KJš
œI˜IJš�œE˜LJ˜—Jš
œ	˜	Jš
œ˜š�œ�œ�œŸ˜6Kš�œ	˜—Kš
œ!™!Kš
œ˜Kš
œ˜K˜Kš
œ1˜1K˜Kšœ�œ�œ˜'Kšœ�œ�œ˜)J˜—J˜šž	œ�œ�œ�œ˜<Jšœ
�œ�œ�œ˜$J˜Jš
œ™J˜Jš�œ�œ
�œ�œ˜(š	�œ
�œ
�œ�œŸ˜<š	�œ�œ
�œ
�œ
�˜7J˜—Kš�œŸœ˜2Kšœ�œ
˜—š�œ˜Kš
œ˜š�œ�œ˜K˜Kšœ�œ
˜—Kš�œ�œ˜—K˜—K˜šžœ�œ
�œ
˜#Jšœ
�œ
�œ˜!J˜Jš
œM™MJšœ�œ
˜
š�œ	�œ˜Kš
œ˜Kš
œ+˜+—š�œ
�œ	�œŸ'˜BKš
œ˜Kš
œ+˜+—š�œ˜Kš
œ˜Kš
œ,˜,—Kš�œ�œ˜"K˜Kš
œ˜Kš
œ˜K˜——J˜š
œ™J˜šžœ�œ
�œ�œ"˜NJš�œ
�œ
˜"J˜Jš
œ˜Jš
œ&˜&š�œ�˜š
œ
˜
Jš
œ˜Jš
œ˜—š
œ˜Jšœ$Ÿ˜=JšœŸ˜—š
œ
˜
Jš
œ
˜
Jš
œ˜—š�œŸ
˜Jš
œ#˜#Jš
œ˜—J˜—Jš�œ�œŸ&�œ�œ�œ˜~˜Jš
œ%˜%Jš
œ9˜9Jš
œ
˜
—J˜—J˜šžœ�œ
�œ�œ˜JJ˜Jš
œ˜š�œ�˜š
œ˜Jš
œ˜JšœŸ	˜—š
œ˜Jš
œ
˜
Jš
œ˜—Jš�œ�œŸ˜%—Jš
œ(˜(Jš
œC˜CJ˜J˜—šžœ�œ#�œ˜QJš
œE˜EJš
œ˜Jš
œ!˜!Jš
œ˜Jš
œ!˜!J˜J˜—šžœ�œ)�œ˜^J˜Jš
œ:™:J˜Jš
œ˜Jš
œ3˜3J˜J˜—šžœ�œ
�œ�œ�œ�œ�œ˜GJš�œ
�œ˜"š�œ˜JšœŸ˜0Jš�œ�œŸ/˜J—Jš
œ
˜
J˜——J˜Kš�œ
˜J˜J™—…—GbfZ