DIRECTORY
Basics USING [BITAND, BITOR, BITSHIFT],
ImagerPixel USING [NewPixelMap, PixelMap, PixelProc],
ImagerSample USING [BitsPerSample, Get, GetBitsPerSample, ObtainScratchSamples, PutSamples, ReleaseScratchSamples, Sample, SampleBuffer, SampleMap],
Real USING [Round],
RealFns USING [ArcTan, CosDeg, SinDeg, SqRt],
Rope USING [ROPE],
Scheme USING [Any, Car, Cdr, Complain, Cons, DefinePrimitive, Environment, false, Flonum, KCheck, ListLength, MakeFixnum, NumberRep, Primitive, ProperList, RegisterInit, TheREAL, true, unspecified, Vector, VectorLength, VectorRep],
SchemeSys USING [CheckForAbort],
SF USING [Box, Vec],
Vector2 USING [VEC];
~
BEGIN
ComplexVectorRep:
TYPE ~
RECORD [
SEQUENCE size:
NAT
OF
VEC];
inexactZero: Scheme.Flonum ~
NEW[Scheme.NumberRep.flonum ¬ [
FALSE, flonum[0.0]]];
Lg:
PROC [n:
CARD]
RETURNS [k:
NAT ¬ 0] ~ {
arg: INT ~ n;
WHILE n # 1
DO
IF n = 0 OR n MOD 2 = 1 THEN ERROR Scheme.Complain[Scheme.MakeFixnum[arg], "is not a power of two"];
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;
};
MakeFlonum:
PROC [x:
REAL]
RETURNS [Scheme.Flonum] ~ {
RETURN [IF x = 0.0 THEN inexactZero ELSE NEW[Scheme.NumberRep.flonum ¬ [FALSE, flonum[x]]]]
};
ComplexVectorRef:
PROC [self: Scheme.Vector, index:
INT]
RETURNS [Scheme.Any] ~ {
data: REF ComplexVectorRep ~ NARROW[self.data];
z: VEC ~ data[index];
IF z.y = 0.0
THEN RETURN [MakeFlonum[z.x]]
ELSE RETURN [NEW[Scheme.NumberRep.complex ¬ [FALSE, complex[MakeFlonum[z.x], MakeFlonum[z.y]]]]]
};
TheComplex:
PROC [value: Scheme.Any]
RETURNS [
VEC] ~ {
z:
VEC ~
WITH value
SELECT
FROM
c: REF Scheme.NumberRep.complex => [x: Scheme.TheREAL[c.x], y: Scheme.TheREAL[c.y]],
ENDCASE => [x: Scheme.TheREAL[value], y: 0.0];
RETURN [z];
};
TheBOOL:
PROC [a: Scheme.Any]
RETURNS [
BOOL] ~ {
RETURN [SELECT a FROM NIL, Scheme.false => FALSE, Scheme.true => TRUE ENDCASE => ERROR Scheme.Complain[a, "should be #t or #f"]]
};
ComplexVectorSet:
PROC [self: Scheme.Vector, index:
INT, value: Scheme.Any] ~ {
data: REF ComplexVectorRep ~ NARROW[self.data];
data[index] ¬ TheComplex[value];
};
MakeComplexVector:
PROC [size:
NAT, initialValue:
VEC]
RETURNS [Scheme.Vector] ~ {
data: REF ComplexVectorRep ~ NEW[ComplexVectorRep[size]];
FOR i:
NAT
IN [0..size)
DO
data[i] ¬ initialValue;
ENDLOOP;
RETURN [NEW[Scheme.VectorRep ¬ [length: size, ref: ComplexVectorRef, set: ComplexVectorSet, data: data]]]
};
TheComplexVector:
PROC [a: Scheme.Any]
RETURNS [
REF ComplexVectorRep] ~ {
WITH a
SELECT
FROM
v: Scheme.Vector => {
WITH v.data
SELECT
FROM
c: REF ComplexVectorRep => RETURN [c];
ENDCASE => NULL;
};
ENDCASE => NULL;
ERROR Scheme.Complain[a, "is not a complex-vector"];
};
ComplexVectorCopy:
PROC [a: Scheme.Any]
RETURNS [Scheme.Any] ~ {
v: REF ComplexVectorRep ~ TheComplexVector[a];
result: Scheme.Any ~ MakeComplexVector[v.size, [0, 0]];
z: REF ComplexVectorRep ~ TheComplexVector[result];
FOR i:
NAT
IN [0..v.size)
DO
z[i] ¬ v[i];
ENDLOOP;
RETURN [result]
};
ComplexVectorTransferReal:
PROC [dst, src: Scheme.Any, start, end:
INT] ~ {
d: REF ComplexVectorRep ~ TheComplexVector[dst];
s: REF ComplexVectorRep ~ TheComplexVector[src];
FOR i:
INT
IN [start..end)
DO
d[i].x ¬ s[i].x
ENDLOOP;
};
ComplexVectorTransferImaginary:
PROC [dst, src: Scheme.Any, start, end:
INT] ~ {
d: REF ComplexVectorRep ~ TheComplexVector[dst];
s: REF ComplexVectorRep ~ TheComplexVector[src];
FOR i:
INT
IN [start..end)
DO
d[i].y ¬ s[i].y
ENDLOOP;
};
ComplexVectorTransferMagnitude:
PROC [dst, src: Scheme.Any, start, end:
INT] ~ {
d: REF ComplexVectorRep ~ TheComplexVector[dst];
s: REF ComplexVectorRep ~ TheComplexVector[src];
FOR i:
INT
IN [start..end)
DO
di: VEC ~ d[i];
si: VEC ~ s[i];
sqrmagd: REAL ~ di.x*di.x + di.y*di.y;
IF sqrmagd # 0.0
THEN {
scale: REAL ~ RealFns.SqRt[(si.x*si.x + si.y*si.y)/sqrmagd];
d[i] ¬ [x: di.x*scale, y: di.y*scale];
};
ENDLOOP;
};
ComplexVectorTransferAngle:
PROC [dst, src: Scheme.Any, start, end:
INT] ~ {
d: REF ComplexVectorRep ~ TheComplexVector[dst];
s: REF ComplexVectorRep ~ TheComplexVector[src];
FOR i:
INT
IN [start..end)
DO
di: VEC ~ d[i];
si: VEC ~ s[i];
sqrmags: REAL ~ si.x*si.x + si.y*si.y;
IF sqrmags # 0.0
THEN {
scale: REAL ~ RealFns.SqRt[(di.x*di.x + di.y*di.y)/sqrmags];
d[i] ¬ [x: si.x*scale, y: si.y*scale];
};
ENDLOOP;
};
ComplexVectorClip:
PROC [a: Scheme.Any, eps:
REAL] ~ {
v: REF ComplexVectorRep ~ TheComplexVector[a];
FOR i:
NAT
IN [0..v.size)
DO
vi: VEC ~ v[i];
IF vi.x < 0 OR ABS[vi.y] > eps THEN v[i] ¬ [0, 0];
ENDLOOP;
};
ThePixelMap:
PROC [a: Scheme.Any]
RETURNS [ImagerPixel.PixelMap] ~ {
WITH a
SELECT
FROM
pm: ImagerPixel.PixelMap => RETURN [pm];
ENDCASE => ERROR Scheme.Complain[a, "not a pixelmap"];
};
IsComplexVector:
PROC [a: Scheme.Any]
RETURNS [
BOOL] ~ {
WITH a
SELECT
FROM
v: Scheme.Vector => {
WITH v.data
SELECT
FROM
c: REF ComplexVectorRep => RETURN [TRUE];
ENDCASE => RETURN [FALSE];
};
ENDCASE => RETURN [FALSE];
};
NatSequenceRep:
TYPE ~
RECORD [
SEQUENCE size:
NAT
OF
NAT];
CAdd:
PROC [a:
VEC, b:
VEC]
RETURNS [
VEC] =
INLINE
{RETURN[[(a.x+b.x), (a.y+b.y)]]}; -- complex sum
CMul:
PROC [a:
VEC, b:
VEC]
RETURNS [
VEC] =
{RETURN[[(a.x*b.x - a.y*b.y), (a.x*b.y + a.y*b.x)]]}; -- complex product
FFTransformRows:
PROC [a: Scheme.Any, inverse:
BOOLEAN] ~ {
n: CARDINAL ¬ Scheme.VectorLength[Scheme.Car[a]];
k: CARDINAL ¬ Lg[n];
w: REF ComplexVectorRep ~ NEW[ComplexVectorRep[n]];
rev: REF NatSequenceRep ~ NEW[NatSequenceRep[n]];
t: REF ComplexVectorRep ~ NEW[ComplexVectorRep[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: Scheme.Any ¬ a, Scheme.Cdr[p]
UNTIL p=
NIL
DO
v: REF ComplexVectorRep ~ TheComplexVector[Scheme.Car[p]];
SchemeSys.CheckForAbort[];
FOR i:
NAT
IN [0 .. n)
DO
t[i] ¬ v[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: VEC ¬ t[i];
tj: VEC ¬ t[j];
t[i] ¬ CAdd[ti, CMul[w[di], tj]];
t[j] ¬ CAdd[ti, CMul[w[dj], tj]];
};
ENDLOOP;
ENDLOOP;
FOR i:
NAT
IN [0 .. n)
DO
ti: VEC ¬ t[i];
IF inverse THEN ti ¬ [ti.x/n, ti.y/n];
v[rev[i]] ¬ ti;
ENDLOOP;
ENDLOOP;
};
InSituTranspose:
PROC [a: Scheme.Any] ~ {
n: NAT ¬ Scheme.VectorLength[Scheme.Car[a]];
i: NAT ¬ 0;
IF n # Scheme.ListLength[a] THEN ERROR Scheme.Complain[$image, "matrix not square"];
FOR p: Scheme.Any ¬ a, Scheme.Cdr[p]
UNTIL p=
NIL
DO
ai: REF ComplexVectorRep ~ TheComplexVector[Scheme.Car[p]];
j: NAT ¬ i;
SchemeSys.CheckForAbort[];
FOR q: Scheme.ProperList ¬
NARROW[p],
NARROW[q.cdr]
UNTIL q=
NIL
DO
v: REF ComplexVectorRep ~ TheComplexVector[q.car];
t: VEC ¬ ai[j];
ai[j] ¬ v[i];
v[i] ¬ t;
j ¬ j + 1;
ENDLOOP;
i ¬ i + 1;
ENDLOOP;
};
Transpose:
PROC [a: Scheme.Any]
RETURNS [result: Scheme.Any] ~ {
n: NAT ¬ Scheme.ListLength[a];
m: NAT ¬ Scheme.VectorLength[Scheme.Car[a]];
new: Scheme.Any ¬ NIL;
FOR i:
NAT
IN [0..m)
DO
new ¬ Scheme.Cons[MakeComplexVector[size: n, initialValue: [0, 0]], new];
ENDLOOP;
result ¬ new;
FOR i:
NAT
IN [0..m)
DO
j: NAT ¬ 0;
v: REF ComplexVectorRep ~ TheComplexVector[Scheme.Car[new]];
FOR q: Scheme.ProperList ¬
NARROW[a],
NARROW[q.cdr]
UNTIL q=
NIL
DO
u: REF ComplexVectorRep ~ TheComplexVector[q.car];
v[j] ¬ u[i];
j ¬ j + 1;
ENDLOOP;
new ¬ Scheme.Cdr[new];
ENDLOOP;
};
FFTransformImage:
PROC [a: Scheme.Any, inverse:
BOOLEAN] ~ {
FFTransformRows[a, inverse];
InSituTranspose[a];
FFTransformRows[a, inverse];
InSituTranspose[a];
};
ComplexImageFromPixelMap:
PROC [pixelMap: ImagerPixel.PixelMap, k:
NAT]
RETURNS [a: Scheme.ProperList ¬
NIL] ~ {
bb: SF.Box ¬ pixelMap.box;
rows, columns: NAT ¬ 1;
map: ImagerSample.SampleMap ~ pixelMap[k];
bits: ImagerSample.BitsPerSample ~ ImagerSample.GetBitsPerSample[map];
maxPixel: REAL ~ Basics.BITSHIFT[1, bits] - 1;
WHILE rows < bb.max.s DO rows ¬ rows * 2 ENDLOOP;
WHILE columns < bb.max.f DO columns ¬ columns * 2 ENDLOOP;
FOR s:
NAT
DECREASING
IN [bb.min.s .. bb.min.s+rows)
DO
v: Scheme.Vector ~ MakeComplexVector[columns, [0.0, 0.0]];
rep: REF ComplexVectorRep ~ NARROW[v.data];
FOR j:
NAT
IN [0 .. columns)
DO
sample: ImagerSample.Sample ~ ImagerSample.Get[map, [MIN[s, bb.max.s-1], MIN[j+bb.min.f, bb.max.f-1]]];
rep[j] ¬ [sample/maxPixel, 0.0];
ENDLOOP;
a ¬ Scheme.Cons[v, a];
ENDLOOP;
};
MagnitudeBounds:
PROC [image: Scheme.Any]
RETURNS [
REAL] ~ {
b: REAL ¬ 0;
FOR q: Scheme.ProperList ¬
NARROW[image],
NARROW[q.cdr]
UNTIL q=
NIL
DO
v: REF ComplexVectorRep ~ TheComplexVector[q.car];
FOR j:
NAT
IN [0 .. v.size)
DO
t: VEC ¬ v[j];
b ¬ MAX[t.x*t.x+t.y*t.y, b];
ENDLOOP;
ENDLOOP;
RETURN [RealFns.SqRt[b]]
};
pi: REAL ~ 3.14159265;
PixelMapFromComplexImage:
PROC [image: Scheme.Any]
RETURNS [ImagerPixel.PixelMap] ~ {
size: SF.Vec ~ [s: Scheme.ListLength[image], f: Scheme.VectorLength[a: Scheme.Car[image]]];
b: REAL ~ MagnitudeBounds[image];
maxSample: ImagerSample.Sample ¬ 255;
MaxSample: ImagerPixel.PixelProc ~ {RETURN [maxSample]};
scale: REAL ¬ IF b = 0.0 THEN 1 ELSE maxSample/b;
angleScale: REAL ¬ IF b = 0.0 THEN 1 ELSE maxSample/pi2;
pixelMap: ImagerPixel.PixelMap ~ ImagerPixel.NewPixelMap[samplesPerPixel: 2, box: [max: size], maxSample: MaxSample];
tail: Scheme.Any ¬ image;
mag: ImagerSample.SampleBuffer ~ ImagerSample.ObtainScratchSamples[size.f];
phase: ImagerSample.SampleBuffer ~ ImagerSample.ObtainScratchSamples[size.f];
FOR s:
NAT
IN [0..
NAT[size.s])
DO
v: REF ComplexVectorRep ~ TheComplexVector[Scheme.Car[tail]];
FOR f:
NAT
IN [0..
NAT[size.f])
DO
t: VEC ~ v[f];
m: REAL ~ RealFns.SqRt[t.x*t.x + t.y*t.y];
p: REAL ¬ RealFns.ArcTan[y: t.y, x: t.x];
WHILE p < 0 DO p ¬ p + pi2 ENDLOOP;
mag[f] ¬ Real.Round[m*scale];
phase[f] ¬ Real.Round[p*angleScale];
ENDLOOP;
ImagerSample.PutSamples[map: pixelMap[0], initIndex: [s, 0], buffer: mag];
ImagerSample.PutSamples[map: pixelMap[1], initIndex: [s, 0], buffer: phase];
tail ¬ Scheme.Cdr[tail];
ENDLOOP;
ImagerSample.ReleaseScratchSamples[mag];
ImagerSample.ReleaseScratchSamples[phase];
RETURN [pixelMap]
};
Op:
TYPE ~ {makecomplexvector, complexvector, complexvectorcopy, fftransformrows, insitutranspose, transpose, fftransformimage, compleximagefrompixelmap, pixelmapfromcompleximage, transferx, transfery, transfermag, transferangle, clip};
FFTPrim:
PROC [self: Scheme.Primitive, a, b, c: Scheme.Any, rest: Scheme.ProperList]
RETURNS [result: Scheme.Any ¬ Scheme.unspecified] ~ {
op: Op ~ NARROW[self.data, REF Op];
SELECT op
FROM
$makecomplexvector => {
result ¬ MakeComplexVector[Scheme.KCheck[a, NAT.LAST], TheComplex[b]]
};
$complexvector => {
result ¬ IF IsComplexVector[a] THEN Scheme.true ELSE Scheme.false
};
$complexvectorcopy => {
result ¬ ComplexVectorCopy[a]
};
$fftransformrows => {
FFTransformRows[a, TheBOOL[b]];
};
$insitutranspose => {
InSituTranspose[a];
};
$transpose => {
result ¬ Transpose[a];
};
$fftransformimage => {
FFTransformImage[a, TheBOOL[b]];
};
$compleximagefrompixelmap => {
pm: ImagerPixel.PixelMap ~ ThePixelMap[a];
result ¬ ComplexImageFromPixelMap[ThePixelMap[a], Scheme.KCheck[b, pm.samplesPerPixel-1]];
};
$pixelmapfromcompleximage => {
result ¬ PixelMapFromComplexImage[a];
};
$transferx => {
ComplexVectorTransferReal[a, b, Scheme.KCheck[c], Scheme.KCheck[rest.car]];
};
$transfery => {
ComplexVectorTransferImaginary[a, b, Scheme.KCheck[c], Scheme.KCheck[rest.car]];
};
$transfermag => {
ComplexVectorTransferMagnitude[a, b, Scheme.KCheck[c], Scheme.KCheck[rest.car]];
};
$transferangle => {
ComplexVectorTransferAngle[a, b, Scheme.KCheck[c], Scheme.KCheck[rest.car]];
};
$clip => {
ComplexVectorClip[a, Scheme.TheREAL[b]];
};
ENDCASE => ERROR;
};
Init:
PROC [env: Scheme.Environment] ~ {
Def:
PROC [name: Rope.
ROPE, nArgs:
NAT, doc: Rope.
ROPE, op: Op] ~ {
Scheme.DefinePrimitive[name: name, nArgs: nArgs, dotted: FALSE, proc: FFTPrim, doc: doc, env: env, data: NEW[Op ¬ op], optional: 0];
};
Def["make-complex-vector", 2, "(k z) make a vector to hold inexact complex numbers", $makecomplexvector];
Def["complex-vector?", 1, "Test for a complex-vector", $complexvector];
Def["complex-vector-copy", 1, "Copy a complex-vector", $complexvectorcopy];
Def["f-f-transform-rows!", 2, "(list inverse-bool) do in-situ FFT on elements of the list", $fftransformrows];
Def["in-situ-transpose!", 1, "Transpose a square image", $insitutranspose];
Def["transpose", 1, "Transpose an image", $transpose];
Def["f-f-transform-image!", 2, "(square-image inverse-bool) do 2-D in-situ FFT on square image", $fftransformimage];
Def["complex-image-from-pixelmap", 2, "(pixelmap layer) make a list-of-complex-vector representation of an image", $compleximagefrompixelmap];
Def["pixelmap-from-complex-image", 1, "(image) make a pixelmap representation of magnitude and phase of an image", $pixelmapfromcompleximage];
Def["complex-vector-transfer-real!", 4, "(dst src start end) transfer real parts", $transferx];
Def["complex-vector-transfer-imaginary!", 4, "(dst src start end) transfer imaginary parts", $transfery];
Def["complex-vector-transfer-magnitude!", 4, "(dst src start end) transfer magnitudes", $transfermag];
Def["complex-vector-transfer-angle!", 4, "(dst src start end) transfer angles", $transferangle];
Def["complex-vector-clip!", 2, "(v eps) set non-(positive-real)s to 0", $clip];
};
Scheme.RegisterInit[Init];