ImageFFTImpl.mesa
Copyright (c) 1984, Xerox Corporation.
Michael Plass, June 28, 1984 10:08:53 am PDT
DIRECTORY Basics, Complex, FS, ImageFFT, ImagerPixelMaps, ImagerPixelRow, IO, Real, RealFns, Rope, SafeStorage, Seq;
ImageFFTImpl: CEDAR PROGRAM
IMPORTS Basics, Complex, FS, ImagerPixelMaps, ImagerPixelRow, IO, Real, RealFns, SafeStorage
EXPORTS ImageFFT ~ BEGIN OPEN ImageFFT;
Error: PUBLIC ERROR [code: ATOM] ~ CODE;
Width: PUBLIC PROC [a: Image] RETURNS [width: NAT] ~ {
width ← NAT.LAST;
FOR p: LIST OF ComplexSequence ← a, p.rest UNTIL p=NIL DO
IF width = NAT.LAST THEN width ← p.first.length;
IF width # p.first.length THEN Error[$InconsistentShape];
ENDLOOP;
};
Height: PUBLIC PROC [a: Image] RETURNS [height: NAT] ~ {
height ← 0;
FOR p: LIST OF ComplexSequence ← a, p.rest UNTIL p=NIL DO
height ← height + 1;
ENDLOOP;
};
Lg: PROC [n: LONG CARDINAL] RETURNS [k: NAT ← 0] ~ {
WHILE n # 1 DO
IF n = 0 OR n MOD 2 = 1 THEN Error[$SizeNotPowerOfTwo];
n ← n / 2;
k ← k + 1;
ENDLOOP;
};
ReverseBits: PROC [i: CARDINAL, k: NAT] RETURNS [r: CARDINAL ← 0] ~ {
THROUGH [0..k) DO
r ← 2*r + i MOD 2;
i ← i / 2;
ENDLOOP;
};
Transpose: PUBLIC PROC [a: Image] RETURNS [b: Image] ~ {
n: NAT ← Width[a];
m: NAT ← Height[a];
i: NAT ← 0;
b ← NIL;
THROUGH [0..n) DO
b ← CONS[NEW[Seq.ComplexSequenceRec[m]], b];
ENDLOOP;
FOR p: LIST OF ComplexSequence ← a, p.rest UNTIL p=NIL DO
ai: ComplexSequence ← p.first;
j: NAT ← 0;
FOR p: LIST OF ComplexSequence ← b, p.rest UNTIL p=NIL DO
p.first[i] ← ai[j];
j ← j + 1;
ENDLOOP;
i ← i + 1;
ENDLOOP;
};
InSituTranspose: PUBLIC PROC [a: LIST OF ComplexSequence] ~ {
n: NAT ← Width[a];
i: NAT ← 0;
IF n # Height[a] THEN Error[$ImageNotSquare];
FOR p: LIST OF ComplexSequence ← a, p.rest UNTIL p=NIL DO
ai: ComplexSequence ← p.first;
j: NAT ← i;
FOR q: LIST OF ComplexSequence ← p, q.rest UNTIL q=NIL DO
t: Complex.Vec ← ai[j];
ai[j] ← q.first[i];
q.first[i] ← t;
j ← j + 1;
ENDLOOP;
i ← i + 1;
ENDLOOP;
};
Transform: PUBLIC PROC [a: Image, inverse: BOOLEAN] ~ {
TransformRows[a, inverse];
InSituTranspose[a];
TransformRows[a, inverse];
InSituTranspose[a];
};
TransformRows: PUBLIC PROC [a: Image, inverse: BOOLEAN] ~ {
n: CARDINAL ← Width[a];
k: CARDINAL ← Lg[n];
w: ComplexSequence ← NEW[Seq.ComplexSequenceRec[n]];
rev: Seq.NatSequence ← NEW[Seq.NatSequenceRec[n]];
t: ComplexSequence ← NEW[Seq.ComplexSequenceRec[n]];
FOR i: CARDINAL IN [0..n) DO
theta: REAL ← 360.0*i/n;
IF inverse THEN theta ← -theta;
w[i] ← [RealFns.CosDeg[theta], RealFns.SinDeg[theta]];
rev[i] ← ReverseBits[i, k];
ENDLOOP;
FOR p: LIST OF ComplexSequence ← a, p.rest UNTIL p=NIL DO
FOR i: NAT IN [0..n) DO
t[i] ← p.first[i];
ENDLOOP;
FOR l: NAT IN [0..k) DO
bit: CARDINAL ← Basics.BITSHIFT[1, k-1-l];
FOR i: CARDINAL IN [0..n) DO
IF Basics.BITAND[i, bit] = 0 THEN {
j: CARDINAL ← Basics.BITOR[i, bit];
di: CARDINAL ← Basics.BITAND[n-1, Basics.BITSHIFT[rev[i], k-1-l]];
dj: CARDINAL ← Basics.BITAND[n-1, Basics.BITSHIFT[rev[j], k-1-l]];
ti: Complex.Vec ← t[i];
tj: Complex.Vec ← t[j];
t[i] ← Complex.Add[ti, Complex.Mul[w[di], tj]];
t[j] ← Complex.Add[ti, Complex.Mul[w[dj], tj]];
};
ENDLOOP;
ENDLOOP;
FOR i: NAT IN [0..n) DO
ti: Complex.Vec ← t[i];
IF inverse THEN ti ← [ti.x/n, ti.y/n];
p.first[rev[i]] ← ti;
ENDLOOP;
ENDLOOP;
};
Mul: PUBLIC PROC [dest, source: Image] ~ {
q: LIST OF ComplexSequence ← source;
FOR p: LIST OF ComplexSequence ← dest, p.rest UNTIL p=NIL DO
FOR j: NAT IN [0..p.first.length) DO
p.first[j] ← Complex.Mul[p.first[j], q.first[j]];
ENDLOOP;
q ← q.rest;
ENDLOOP;
};
ScalarMul: PUBLIC PROC [a: Image, scalar: Complex.Vec] ~ {
FOR p: LIST OF ComplexSequence ← a, p.rest UNTIL p=NIL DO
FOR j: NAT IN [0..p.first.length) DO
p.first[j] ← Complex.Mul[p.first[j], scalar];
ENDLOOP;
ENDLOOP;
};
Add: PUBLIC PROC [dest, source: Image] ~ {
q: LIST OF ComplexSequence ← source;
FOR p: LIST OF ComplexSequence ← dest, p.rest UNTIL p=NIL DO
FOR j: NAT IN [0..p.first.length) DO
p.first[j] ← Complex.Add[p.first[j], q.first[j]];
ENDLOOP;
q ← q.rest;
ENDLOOP;
};
Sub: PUBLIC PROC [dest, source: Image] ~ {
q: LIST OF ComplexSequence ← source;
FOR p: LIST OF ComplexSequence ← dest, p.rest UNTIL p=NIL DO
FOR j: NAT IN [0..p.first.length) DO
p.first[j] ← Complex.Sub[p.first[j], q.first[j]];
ENDLOOP;
q ← q.rest;
ENDLOOP;
};
Store: PUBLIC PROC [image: Image, fileName: Rope.ROPE] ~ {
n: INT ← Width[image];
IF Height[image] # n THEN Error[$ImageNotSquare]
ELSE {
stream: IO.STREAMFS.StreamOpen[fileName: fileName, accessOptions: $create, createByteCount: 8*n*n];
bytesPerRow: INT ← Basics.bytesPerWord*SIZE[Complex.Vec]*n;
FOR p: LIST OF ComplexSequence ← image, p.rest UNTIL p=NIL DO
TRUSTED {
dataPointer: LONG POINTER ← @p.first[0];
IO.UnsafePutBlock[stream, [base: dataPointer, count: bytesPerRow]];
};
ENDLOOP;
IO.Close[stream];
};
};
Load: PUBLIC PROC [fileName: Rope.ROPE, scratch: Image ← NIL] RETURNS [Image] ~ {
stream: IO.STREAMFS.StreamOpen[fileName: fileName, accessOptions: $read, streamOptions: [tiogaRead: FALSE, commitAndReopenTransOnFlush: TRUE, truncatePagesOnClose: TRUE, finishTransOnClose: TRUE, closeFSOpenFileOnClose: TRUE]];
b: Image ← NIL;
elementSize: NAT ~ Basics.bytesPerWord*SIZE[Complex.Vec];
BEGIN ENABLE UNWIND => {IO.Close[stream]};
lgnsrq: NAT;
n: NAT;
bytes: LONG CARDINALIO.GetLength[stream];
bytesPerRow: INT;
IF bytes MOD elementSize # 0 THEN Error[$InvalidFileSize];
lgnsrq ← Lg[bytes/elementSize];
IF lgnsrq MOD 2 # 0 THEN Error[$InvalidFileSize];
n ← Basics.BITSHIFT[1, lgnsrq/2];
bytesPerRow ← elementSize*n;
FOR i: NAT IN [0..n) DO
IF scratch#NIL AND scratch.first.length = n THEN {
t: LIST OF ComplexSequence ← scratch;
scratch ← scratch.rest;
t.rest ← b;
b ← t;
}
ELSE b ← CONS[NEW[Seq.ComplexSequenceRec[n]], b];
TRUSTED {
dataPointer: LONG POINTER ← @b.first[0];
IF bytesPerRow # IO.UnsafeGetBlock[stream, [base: dataPointer, count: bytesPerRow]] THEN Error[$CantHappen];
};
ENDLOOP;
END;
IO.Close[stream];
scratch ← Destroy[scratch];
RETURN [RevList[b]]
};
RevList: PROC [a: LIST OF ComplexSequence] RETURNS [b: LIST OF ComplexSequence] ~ {
b ← NIL;
WHILE a # NIL DO
t: LIST OF ComplexSequence ← a;
a ← a.rest;
t.rest ← b;
b ← t;
ENDLOOP;
};
Copy: PUBLIC PROC [a: Image, scratch: Image] RETURNS [b: Image] ~ {
n: NAT ← Width[a];
b ← NIL;
FOR p: LIST OF ComplexSequence ← a, p.rest UNTIL p=NIL DO
IF scratch#NIL AND scratch.first.length = n THEN {
t: LIST OF ComplexSequence ← scratch;
scratch ← scratch.rest;
t.rest ← b;
b ← t;
}
ELSE b ← CONS[NEW[Seq.ComplexSequenceRec[n]], b];
FOR j: NAT IN [0..p.first.length) DO
b.first[j] ← p.first[j];
ENDLOOP;
ENDLOOP;
scratch ← Destroy[scratch];
b ← RevList[b];
};
Destroy: PUBLIC PROC [a: Image] RETURNS [Image ← NIL] ~ {
IF a # NIL THEN {
WHILE a = NIL DO
t: Image ← a;
a ← a.rest;
t.rest ← NIL;
ENDLOOP;
SafeStorage.ReclaimCollectibleObjects[];
};
};
TestSeq: PROC [n: NAT, j0, j1: NAT] RETURNS [LIST OF ComplexSequence] ~ {
a: ComplexSequence ← NEW[Seq.ComplexSequenceRec[n]];
FOR i: NAT IN [0..n) DO
a[i] ← [IF i IN [j0..j1) THEN 1.0 ELSE 0.0, 0];
ENDLOOP;
RETURN [LIST[a]]
};
FromPixelMap: PUBLIC PROC [pixelMap: PixelMap] RETURNS [a: Image ← NIL] ~ {
bb: ImagerPixelMaps.DeviceRectangle ← pixelMap.Window;
rows, columns: NAT ← 1;
maxPixel: REAL ← Basics.BITSHIFT[1, Basics.BITSHIFT[1, pixelMap.refRep.lgBitsPerPixel]]-1;
row: ImagerPixelRow.PixelRow;
WHILE rows < pixelMap.sSize DO rows ← rows * 2 ENDLOOP;
WHILE columns < pixelMap.fSize DO columns ← columns * 2 ENDLOOP;
row ← ImagerPixelRow.CreatePixelRow[bb.sMin+rows, bb.fMin, columns];
FOR s: NAT DECREASING IN [bb.sMin..bb.sMin+rows) DO
a ← CONS[NEW[Seq.ComplexSequenceRec[columns]], a];
row.sOrigin ← s;
ImagerPixelRow.ClearPixelRow[row];
ImagerPixelRow.LoadPixelRow[row, pixelMap];
FOR j: NAT IN [0..columns) DO
a.first[j] ← [row[j]/maxPixel, 0];
ENDLOOP;
ENDLOOP;
};
Range: PUBLIC PROC [a: Image] RETURNS [xmin, ymin, xmax, ymax: REAL] ~ {
xmin ← ymin ← Real.LargestNumber;
xmax ← ymax ← -Real.LargestNumber;
FOR p: LIST OF ComplexSequence ← a, p.rest UNTIL p=NIL DO
ai: ComplexSequence ← p.first;
FOR j: NAT IN [0..p.first.length) DO
aij: Complex.Vec ← ai[j];
xmin ← MIN[xmin, aij.x];
ymin ← MIN[ymin, aij.y];
xmax ← MAX[xmax, aij.x];
ymax ← MAX[ymax, aij.y];
ENDLOOP;
ENDLOOP;
};
MaskedRange: PUBLIC PROC [a: Image, bitmap: PixelMap] RETURNS [xmin, ymin, xmax, ymax: REAL] ~ {
i: NAT ← 0;
IF bitmap.sSize # Height[a]
OR bitmap.fSize # Width[a]
OR bitmap.sOrigin + bitmap.sMin # 0
OR bitmap.fOrigin + bitmap.fMin # 0
OR bitmap.refRep.lgBitsPerPixel # 0
THEN Error[$NonconformingParameters];
xmin ← ymin ← Real.PlusInfinity;
xmax ← ymax ← Real.MinusInfinity;
FOR p: LIST OF ComplexSequence ← a, p.rest UNTIL p=NIL DO
ai: ComplexSequence ← p.first;
FOR j: NAT IN [0..p.first.length) DO
IF bitmap.GetBit[i, j] = 1 THEN {
aij: Complex.Vec ← ai[j];
xmin ← MIN[xmin, aij.x];
ymin ← MIN[ymin, aij.y];
xmax ← MAX[xmax, aij.x];
ymax ← MAX[ymax, aij.y];
};
ENDLOOP;
i ← i + 1;
ENDLOOP;
};
AbsBound: PROC [z: Complex.Vec] RETURNS [REAL] ~ {
x: REALABS[z.x];
y: REALABS[z.y];
RETURN [MIN[MAX[x, y], 1.414214*(x+y)]]
};
TestAbsBound: PROC RETURNS [min, max: REAL] ~ {
min ← 1000;
max ← -1000;
FOR i: NAT IN [0..360) DO
z: Complex.Vec ← [RealFns.CosDeg[i], RealFns.SinDeg[i]];
m: REAL ← AbsBound[z];
min ← MIN[min, m];
max ← MAX[max, m];
ENDLOOP;
};
TransferToPixels: PUBLIC PROC [dest: ImagerPixelMaps.PixelMap, source: Image, scalar: REAL] ~ {
row: ImagerPixelRow.PixelRow ← ImagerPixelRow.CreatePixelRow[0, 0, Width[source]];
maxPixel: REAL ← Basics.BITSHIFT[1, Basics.BITSHIFT[1, dest.refRep.lgBitsPerPixel]]-0.01;
FOR p: LIST OF ComplexSequence ← source, p.rest UNTIL p=NIL DO
FOR j: NAT IN [0..row.fSize) DO
mag: REAL ← AbsBound[p.first[j]];
row[j] ← Real.FixC[MIN[MAX[mag*scalar, 0], maxPixel]];
ENDLOOP;
ImagerPixelRow.StorePixelRow[row, dest];
row.sOrigin ← row.sOrigin + 1;
ENDLOOP;
};
MakeMask: PUBLIC PROC [a: Image, threshold: REAL] RETURNS [bitmap: PixelMap] ~ {
tt: REAL ← threshold*threshold;
i: NAT ← 0;
bitmap ← ImagerPixelMaps.Create[0, [0, 0, Height[a], Width[a]]];
bitmap.Clear;
FOR p: LIST OF ComplexSequence ← a, p.rest UNTIL p=NIL DO
FOR j: NAT IN [0..p.first.length) DO
IF Complex.SqrAbs[p.first[j]] >= tt THEN bitmap.Fill[[i, j, 1, 1], 1];
ENDLOOP;
i ← i+1;
ENDLOOP;
};
ApplyMask: PUBLIC PROC [a: Image, mask: PixelMap] ~ {
i: NAT ← 0;
IF mask.sSize # Height[a]
OR mask.fSize # Width[a]
OR mask.sOrigin + mask.sMin # 0
OR mask.fOrigin + mask.fMin # 0
OR mask.refRep.lgBitsPerPixel # 0
THEN Error[$NonconformingParameters];
FOR p: LIST OF ComplexSequence ← a, p.rest UNTIL p=NIL DO
FOR j: NAT IN [0..p.first.length) DO
IF mask.GetBit[i, j] = 1 THEN p.first[j] ← [0, 0];
ENDLOOP;
i ← i+1;
ENDLOOP;
};
END.