<<>> <> <> <> <<>> 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]; SchemeFFTImpl: CEDAR PROGRAM IMPORTS Basics, ImagerSample, ImagerPixel, Real, RealFns, Scheme, SchemeSys ~ BEGIN VEC: TYPE ~ Vector2.VEC; 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; pi2: REAL ¬ pi+pi; 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]; END.