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.STREAM ← FS.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.STREAM ← FS.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 CARDINAL ← IO.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: REAL ← ABS[z.x];
y: REAL ← ABS[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.