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;
Main Program Loop
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};
};
Beginning of Transform
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[];
first pass
Image ← oldImage;
NewImage ← passInfo.image;
Dir ← passInfo.pass1;
SecondPass ← FALSE;
Pass[ passInfo.matrix];
second pass
Image ← passInfo.image;
NewImage ← newImage;
Dir ← passInfo.pass2;
SecondPass ← TRUE;
Pass[ passInfo.matrix, passInfo.min];
clean up global storage before exiting
Image ← NIL;
NewImage ← NIL;
EndPixels ← NIL;
};
Pass:
PROC[ m: Matrix3d.Matrix4by4,
min:
INTEGER ← 0] ~ {
loop through all lines in the image and transform each
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] ~ {
transform a line
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] ~ {
transform a pixel
b1, b2: REAL;
oldLow, oldHigh: INTEGER;
total, totalUnWeight, totalWeight, unRatio: REAL ← 0.0;
colorTest: INTEGER ← 0;
newColor: [0..256) ← 0;
back-transform to find the old pixels that will fall under this new pixel
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
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.
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];
find values for this pixel and add to totals
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;
average the values and blend with background to get final pixel value
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];
};
};
Area Calculation
runs an individual pixel through the entire matrix transform (including final translation)
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 image scaling is more than 10x in either direction, don't do it (you can't get a good picture anyhow at this point)
IF newZ < m[1][1]/10
OR newZ < m[2][2]/10
THEN
RETURN[ [0, 0], TRUE];
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.
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 this is true, then the area is so small it won't change the value of any of the pixels it overlaps
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]};
};
Beginning of InterArea
corner1, corner2, corner3, corner4: PixelPos;
boundsFault: BOOLEAN ← FALSE;
a, b: INTEGER;
transform first but not second variable
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]];
};
convert corner points into vectors
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];
| (side1 X side2) + (side3 X side2) / 2 + (side1 X side4) / 2 | is area of quadrilateral
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];
};
Matrix Operations
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] ~ {
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.
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
};
Bounds Checking
GetLineBounds:
PROC[ a, b, c, d, e:
REAL]
RETURNS[ low, high: INTEGER] ~ {
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.
b1, b2, den1, den2: REAL;
lastPixel, lastNewPixel: INTEGER;
in the '90' passes, the input is rotated so x bounds become y bounds and vice versa
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];
get real values of edges of image
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] ~ {
assumes lower bound is 0
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] ~ {
convert real boundaries to integers to find pixels which overlap the boundary
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);
};
Type Conversion
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] ~ {
returns location of the pixel in screen coordinates (s, f)
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;
};