<> <> <> 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.