
<<ImageFFTImpl.mesa>>
<<Copyright (c) 1984, Xerox Corporation.>>
<<Michael Plass, July 30, 1985 2:28:15 pm PDT>>
    <<Tim Diebert: June 7, 1985 2:13:33 pm PDT>>


DIRECTORY Basics, Complex, FS, ImageFFT, ImagerPixelMap, IO, Real, RealFns, Rope, SafeStorage, Seq, Process;

ImageFFTImpl: CEDAR PROGRAM
IMPORTS Basics, Complex, FS, ImagerPixelMap, IO, Real, RealFns, SafeStorage, Process
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;
Process.CheckForAbort[];
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;
Process.CheckForAbort[];
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
Process.CheckForAbort[];
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: ImagerPixelMap.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
sample: CARDINAL = SELECT pixelMap.refRep.lgBitsPerPixel FROM
0 => ImagerPixelMap.GetBit[pixelMap, s, j],
1 => ImagerPixelMap.Get2Bits[pixelMap, s, j],
2 => ImagerPixelMap.Get4Bits[pixelMap, s, j],
3 => ImagerPixelMap.Get8Bits[pixelMap, s, j],
4 => ImagerPixelMap.Get16Bits[pixelMap, s, j],
ENDCASE => ERROR;
a.first[j] _ [sample/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: ImagerPixelMap.PixelMap, source: Image, scalar: REAL] ~ {
<<row: ImagerPixelRow.PixelRow _ ImagerPixelRow.CreatePixelRow[0, 0, Width[source]];>>
row: CARDINAL_ 0;
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 .. Width[source]) DO
mag: REAL _ AbsBound[p.first[j]];
<<row[j] _ Real.FixC[MIN[MAX[mag*scalar, 0], maxPixel]];>>
ImagerPixelMap.PutPixel[dest, row, j, Real.FixC[MIN[MAX[mag*scalar, 0], maxPixel]]];
ENDLOOP;
row_ row + 1;
<<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 _ ImagerPixelMap.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.