DIRECTORY AIS, Basics, BasicTime, Commander, Complex, Convert, FS, ImageFFT, Imager, ImagerPixelMap, InterminalBackdoor, IO, Process, Real, RealFns, Rope, Seq, Terminal;

ImageFFTCommandsImpl: CEDAR PROGRAM
IMPORTS AIS, Basics, Commander, Complex, Convert, FS, ImageFFT, ImagerPixelMap, InterminalBackdoor, IO, Process, Real, RealFns, Rope, Terminal ~ BEGIN

ROPE: TYPE ~ Rope.ROPE;

PixelMapFromAIS: PROC [name: ROPE] RETURNS [map: ImagerPixelMap.PixelMap, rast: AIS.Raster] ~ TRUSTED {
ais: AIS.FRef _ AIS.OpenFile[name];
window: AIS.WRef _ AIS.OpenWindow[ais];
lineMap: ImagerPixelMap.PixelMap;
lineBufferDesc: AIS.Buffer;
rast _ AIS.ReadRaster[ais];
lineMap _ ImagerPixelMap.Create[BitsToLog[rast.bitsPerPixel], [window.GetWindowParams.firstScan, window.GetWindowParams.firstPixel, 1, window.GetWindowParams.lastPixel+1-window.GetWindowParams.firstPixel]];
map _ ImagerPixelMap.Create[BitsToLog[rast.bitsPerPixel], [window.GetWindowParams.firstScan, window.GetWindowParams.firstPixel, window.GetWindowParams.lastScan+1-window.GetWindowParams.firstScan, window.GetWindowParams.lastPixel+1-window.GetWindowParams.firstPixel]];
lineBufferDesc _ [length: lineMap.refRep.words, addr: lineMap.refRep.pointer];
FOR i: NAT IN [window.GetWindowParams.firstScan..window.GetWindowParams.lastScan] DO
AIS.UnsafeReadLine[window, lineBufferDesc, i];
map.Transfer[lineMap];
lineMap.sOrigin _ lineMap.sOrigin + 1;
ENDLOOP;
};

BitsToLog: PROC [bits: CARDINAL] RETURNS [log: [0 .. 4]] = BEGIN
log _ SELECT bits FROM
1 => 0,  2 => 1,  4 => 2, 8 => 3, 16 => 4, ENDCASE => ERROR;
END;

CachedFFTName: PROC [aisFileName: ROPE] RETURNS [fftName: ROPE] ~ TRUSTED {
cp: FS.ComponentPositions;
fullFName: ROPE;
created: BasicTime.GMT;
stream: IO.STREAM _ IO.ROS[];
[fullFName: fullFName, created: created] _ FS.FileInfo[aisFileName];
cp _ FS.ExpandName[fullFName].cp;
stream.PutF["%g-%g.FFT", IO.rope[fullFName.Substr[cp.base.start, cp.base.length]], IO.int[LOOPHOLE[created]]];
fftName _ stream.RopeFromROS[close: TRUE];
};

StorePixelMap: PROC [aisFileName: ROPE, source: ImagerPixelMap.PixelMap, bitmap: BOOLEAN _ TRUE, comment: ROPE _ NIL] ~ TRUSTED {
output: AIS.FRef _ AIS.CreateFile[name: aisFileName, raster: NEW[AIS.RasterPart _ [
scanCount: source.sSize,
scanLength: source.fSize,
scanMode: rd,
bitsPerPixel: IF source.refRep.lgBitsPerPixel = 0 AND bitmap THEN 0 ELSE Basics.BITSHIFT[1, source.refRep.lgBitsPerPixel],
linesPerBlock: -1,
paddingPerBlock: 65535
]]];
outputWindow: AIS.WRef _ AIS.OpenWindow[output];
lineMap: ImagerPixelMap.PixelMap _ ImagerPixelMap.Create[source.refRep.lgBitsPerPixel, [source.sOrigin+source.sMin, source.fOrigin+source.fMin, 1, source.fSize]];
lineBufferDesc: AIS.Buffer _ [length: lineMap.refRep.words, addr: lineMap.refRep.pointer];
AIS.WriteComment[output, comment];
FOR i: NAT IN [0..source.sSize) DO
lineMap.Clear;
lineMap.Transfer[source];
lineMap.sOrigin _ lineMap.sOrigin + 1;
AIS.UnsafeWriteLine[outputWindow, lineBufferDesc, i];
ENDLOOP;
AIS.CloseFile[output];
};

GetAISComment: PROC [aisFileName: ROPE] RETURNS [comment: ROPE] ~ {
ais: AIS.FRef _ AIS.OpenFile[name: aisFileName];
comment _ AIS.ReadComment[ais];
AIS.CloseFile[ais];
};

Padded: PROC [pixelMap: ImagerPixelMap.PixelMap] RETURNS [padded: ImagerPixelMap.PixelMap] ~ {
d: INT _ MAX[pixelMap.fSize, pixelMap.sSize];
n: INT _ 1;
WHILE n < d DO n _ n+n ENDLOOP;
padded _ ImagerPixelMap.Create[pixelMap.refRep.lgBitsPerPixel, [0, 0, n, n]];
padded.Clear;
padded.Transfer[pixelMap];
IF pixelMap.refRep.lgBitsPerPixel = 0 THEN padded.Fill[[0, 0, n, n], 1, [xor, null]];
};

GetTransform: PROC [inputName: ROPE, msg: IO.STREAM _ NIL] RETURNS [ImageFFT.Image] ~ {
fftName: ROPE _ CachedFFTName[inputName];
original: ImagerPixelMap.PixelMap _ PixelMapFromAIS[inputName].map;
padded: ImagerPixelMap.PixelMap _ Padded[original];
result: ImagerPixelMap.PixelMap _ ImagerPixelMap.Create[3, original.Window];
image: ImageFFT.Image _ NIL;
image _ ImageFFT.Load[fftName ! FS.Error => CONTINUE];
IF image = NIL THEN {
image _ ImageFFT.FromPixelMap[padded];
IF msg # NIL THEN msg.PutF["Generating FFT of %g . . . ", IO.rope[inputName]];
ImageFFT.Transform[image, FALSE];
ImageFFT.Store[image, fftName];
IF msg # NIL THEN msg.PutF["%g written.\n", IO.rope[fftName]];
};
RETURN [image]
};

FilterAIS: PROC [outputName, inputName: ROPE, filters: LIST OF FilterRep, msg: IO.STREAM _ NIL, comment: ROPE] ~ {
fftName: ROPE _ CachedFFTName[inputName];
original: ImagerPixelMap.PixelMap _ PixelMapFromAIS[inputName].map;
padded: ImagerPixelMap.PixelMap _ Padded[original];
result: ImagerPixelMap.PixelMap _ ImagerPixelMap.Create[3, original.Window];
image: ImageFFT.Image _ NIL;
aisloc: INT _ Rope.Find[outputName, ".ais", 0, FALSE];
imageName: ROPE _ IF aisloc = -1 THEN outputName.Concat[".image"] ELSE outputName.Replace[aisloc, 4, ".image"];
comment _ GetAISComment[inputName].Cat["; ", comment];
image _ ImageFFT.Load[fftName ! FS.Error => CONTINUE];
IF image = NIL THEN {
image _ ImageFFT.FromPixelMap[padded];
IF msg # NIL THEN msg.PutF["Generating FFT of %g . . . ", IO.rope[inputName]];
ImageFFT.Transform[image, FALSE];
ImageFFT.Store[image, fftName];
IF msg # NIL THEN msg.PutF["%g written.\n", IO.rope[fftName]];
};
IF msg # NIL THEN msg.PutRope["Filtering . . . "];
ApplyFilters[image, filters];
IF msg # NIL THEN msg.PutRope["Applying inverse transform . . . "];
ImageFFT.Transform[image, TRUE];
ImageFFT.TransferToPixels[result, image, 255];
IF msg # NIL THEN msg.PutF["Writing %g . . . ", IO.rope[imageName]];
ImageFFT.Store[image, imageName];
IF msg # NIL THEN msg.PutF["Writing %g . . . ", IO.rope[outputName]];
image _ ImageFFT.Destroy[image];
StorePixelMap[outputName, result, FALSE, comment];
};

Break: PROC [char: CHAR] RETURNS [IO.CharClass] ~ {
IF char = '_ OR char = '; THEN RETURN [break];
IF char = '  OR char = '	 OR char = ', OR char = '\n THEN RETURN [sepr];
RETURN [other];
};

GetToken: PROC [stream: IO.STREAM] RETURNS [rope: ROPE _ NIL] = {
rope _ stream.GetTokenRope[Break ! IO.EndOfStream => CONTINUE].token;
};

GetInt: PROC [stream: IO.STREAM] RETURNS [int: INT] = {
int _ INT.FIRST;
int _ stream.GetInt[! IO.EndOfStream, IO.Error => CONTINUE];
};

ColorDisplay: PROC [sSize, fSize: NAT] RETURNS [colorDisplay: ImagerPixelMap.PixelMap] ~ TRUSTED {
fb: Terminal.FrameBuffer _ Terminal.GetColorFrameBufferA[InterminalBackdoor.terminal];
colorDisplay _ PixelMapFromFrameBuffer[fb];
colorDisplay _ colorDisplay.ShiftMap[(sSize-colorDisplay.sSize)/2, (fSize-colorDisplay.fSize)/2].Clip[[0, 0, sSize+1, fSize+1]];
};

AbsBound: PROC [z: Complex.VEC] RETURNS [REAL] ~ INLINE {
x: REAL _ ABS[z.x];
y: REAL _ ABS[z.y];
RETURN [MIN[MAX[x, y], 1.414214*(x+y)]]
};

dispScale: REAL _ 40;
dispOffset: REAL _ 40;
Show: PROC [a: ImageFFT.Image] ~ {
height: NAT _ ImageFFT.Height[a];
width: NAT _ ImageFFT.Width[a];
cd: ImagerPixelMap.PixelMap _ ColorDisplay[height, width];
row: CARDINAL _ height/2;
Add: PROC [n: NAT] RETURNS [new: NAT] = INLINE BEGIN
new_ n + width/2; IF new > width THEN new_ new - width;
END;
FOR p: ImageFFT.Image _ a, p.rest UNTIL p=NIL DO
FOR j: NAT IN [0..width) DO
mag: REAL _ AbsBound[p.first[j]];
ln: REAL _ IF mag <= 0.0 THEN 0.0 ELSE RealFns.Ln[mag];
pixel: INT _ Real.RoundLI[ln*dispScale+dispOffset];
ImagerPixelMap.PutPixel[cd, row, Add[j], MIN[MAX[pixel,0],255]];
ENDLOOP;
row_ row + 1;  IF row > height THEN row _ row - height;
ENDLOOP;
};

RoundLowPassFilter: FilterProc ~ {
RETURN [IF Sqr[ABS[z.x]] + Sqr[ABS[z.y]] > Sqr[param[0]] THEN 0.0 ELSE 1.0]
};

SpotFilter: FilterProc ~ {
radius: REAL ~ param[2];
m: REAL ~ param[3];
diam: REAL ~ 2*radius;
z _ [ABS[z.x-param[0]], ABS[z.y-param[1]]];
IF z.x > diam OR z.y > diam THEN RETURN [1.0]
ELSE {
d: REAL _ Complex.Abs[z] / radius;
IF d > 2 THEN RETURN [1.0];
d _ (d-2)*(d+2);
d _ d*d;
d _ d/16*(m-1) + 1;
RETURN [d]
};
};

Sqr: PROC [r: REAL] RETURNS [REAL] ~ INLINE {RETURN [r*r]};

BandBoostFilter: FilterProc ~ {
fuzz: REAL ~ param[2]/2;
m: REAL ~ param[3];
d: REAL _ Complex.SqrAbs[z];
IF d <= Sqr[param[0]-fuzz] OR d >= Sqr[param[1]+fuzz] THEN RETURN [1.0];
IF d >= Sqr[param[0]+fuzz] AND d <= Sqr[param[1]-fuzz] THEN RETURN [m];
d _ Real.SqRt[d];
IF d < param[0] + fuzz THEN d _ (d - param[0])/fuzz
ELSE d _ (param[1] - d)/fuzz;
d _ (d+1)/2;
RETURN [1 + d*(m-1)];
};

FilterProc: TYPE ~ PROC [param: ARRAY [0..maxParam) OF REAL, z: Complex.VEC] RETURNS [REAL];

FilterRep: TYPE ~ RECORD [
filterProc: FilterProc,
param: ARRAY [0..maxParam) OF REAL _ ALL[0]
];

FilterEntry: TYPE ~ RECORD [
filterName: ROPE,
nParam: NAT _ 0,
filterProc: FilterProc,
doc: ROPE
];

FindFilterEntry: PROC [name: ROPE] RETURNS [filterEntry: FilterEntry] ~ {
FOR p: LIST OF FilterEntry _ filterList, p.rest UNTIL p = NIL DO
IF name.Equal[p.first.filterName, FALSE] THEN RETURN [p.first]
ENDLOOP;
RETURN [[NIL, 0, NIL, NIL]]
};

ParseFilterSpec: PROC [stream: IO.STREAM, msg: IO.STREAM] RETURNS [filters: LIST OF FilterRep] ~ {
p: LIST OF FilterRep _ NIL;
token: ROPE _ GetToken[stream];
param: ARRAY [0..maxParam) OF REAL _ ALL[0];
nParam: NAT _ 0;
WHILE token.Length > 0 DO
r: REAL _ -99999999;
r _ Convert.RealFromRope[token ! Convert.Error => CONTINUE];
IF r = -99999999 THEN r _ Convert.IntFromRope[token ! Convert.Error => CONTINUE];
IF r = -99999999 THEN {
filterEntry: FilterEntry _ FindFilterEntry[token];
IF filterEntry.filterName = NIL THEN {
msg.PutF["%g: Unknown filter", IO.rope[token]];
ERROR ABORTED;
};
IF nParam # filterEntry.nParam THEN {
msg.PutF["%g: Wrong number of parameters", IO.rope[token]];
ERROR ABORTED;
};
p _ CONS[[filterEntry.filterProc, param], p];
nParam _ 0;
param _ ALL[0];
}
ELSE {
param[nParam] _ r;
nParam _ nParam + 1;
IF nParam > maxParam THEN {
msg.PutF["%g: Too many parameters", IO.rope[token]];
ERROR ABORTED;
};
};
token _ GetToken[stream];
ENDLOOP;
WHILE p # NIL DO
t: LIST OF FilterRep _ p;
p _ p.rest;
t.rest _ filters;
filters _ t;
ENDLOOP;
};

ApplyFilters: PROC [image: ImageFFT.Image, filters: LIST OF FilterRep] ~ {
n: NAT _ ImageFFT.Height[image];
ninv2: REAL _ 2.0/n;
s: REAL _ 0;
FOR p: ImageFFT.Image _ image, p.rest UNTIL p = NIL DO
row: Seq.ComplexSequence _ p.first;
FOR f: NAT IN [0..n) DO
fr: REAL _ IF 2*f > n THEN f-n ELSE f;
z: Complex.VEC _ [ABS[fr]*ninv2, ABS[s]*ninv2];
aij: Complex.VEC _ row[f];
m: REAL _ 1.0;
FOR f: LIST OF FilterRep _ filters, f.rest UNTIL f = NIL DO
m _ m*f.first.filterProc[f.first.param, z];
ENDLOOP;
row[f] _ [aij.x*m, aij.y*m];
ENDLOOP;
s _ s + 1.0;
IF s > n/2 THEN s _ s - n;
ENDLOOP;
};

FFTShowCommand: Commander.CommandProc ~ {
stream: IO.STREAM _ IO.RIS[cmd.commandLine];
inputName: ROPE _ GetToken[stream];
filters: LIST OF FilterRep _ ParseFilterSpec[stream, cmd.out];
image: ImageFFT.Image _ GetTransform[inputName, cmd.out];
ApplyFilters[image, filters];
Show[image];
image _ image.Destroy;
};

FFTFilterCommand: Commander.CommandProc ~ {
stream: IO.STREAM _ IO.RIS[cmd.commandLine];
outputName: ROPE _ GetToken[stream];
gets: ROPE _ GetToken[stream];
inputName: ROPE _ GetToken[stream];
oldPriority: Process.Priority _ Process.GetPriority[];
filters: LIST OF FilterRep _ ParseFilterSpec[stream, cmd.out];
IF NOT gets.Equal["_"] THEN {
cmd.out.PutRope["Specify output _ input <filterspec>, please.\n"];
RETURN;
};
Process.SetPriority[Process.priorityBackground];
FilterAIS[outputName, inputName, filters, cmd.out, cmd.commandLine ! UNWIND => {Process.SetPriority[oldPriority]}];
Process.SetPriority[oldPriority];
cmd.out.PutRope["Done.\n"];
};

ClipAISCommand: Commander.CommandProc ~ {
stream: IO.STREAM _ IO.RIS[cmd.commandLine];
outputName: ROPE _ GetToken[stream];
gets: ROPE _ GetToken[stream];
inputName: ROPE _ GetToken[stream];
skipLines: INT _ GetInt[stream];
skipDots: INT _ GetInt[stream];
lineCount: INT _ GetInt[stream];
dotCount: INT _ GetInt[stream];
pixelMap: ImagerPixelMap.PixelMap;
raster: AIS.Raster;
IF NOT (gets.Equal["_"] AND skipLines >= 0 AND skipDots >= 0 AND lineCount > 0 AND dotCount > 0)  THEN {
cmd.out.PutRope["Specify output _ input <skipLines> <skipDots> <lineCount> <dotCount>, please.\n"];
RETURN;
};
[pixelMap, raster] _ PixelMapFromAIS[inputName];
StorePixelMap[outputName, pixelMap.Clip[[skipLines, skipDots, lineCount, dotCount]], raster.bitsPerPixel = 0, GetAISComment[inputName].Cat["; ", cmd.commandLine]];
cmd.out.PutRope["Done.\n"];
};

FFTPeekCommand: Commander.CommandProc ~ {
stream: IO.STREAM _ IO.RIS[cmd.commandLine];
inputName: ROPE _ GetToken[stream];
kind: ROPE _ GetToken[stream];
skipLines: NAT _ GetInt[stream];
skipDots: NAT _ GetInt[stream];
lineCount: NAT _ GetInt[stream];
dotCount: NAT _ GetInt[stream];
realPart: BOOLEAN _ kind.Length = 0 OR kind.Substr[0, 1].Equal["R", FALSE];
elementSize: NAT ~ Basics.bytesPerWord*SIZE[Complex.VEC];
stream _ FS.StreamOpen[fileName: inputName, accessOptions: $read, streamOptions: [tiogaRead: FALSE, commitAndReopenTransOnFlush: TRUE, truncatePagesOnClose: TRUE, finishTransOnClose: TRUE, closeFSOpenFileOnClose: TRUE]];
BEGIN ENABLE UNWIND => {IO.Close[stream]};
Lg: PROC [n: LONG CARDINAL] RETURNS [k: NAT _ 0] ~ {
WHILE n # 1 DO
IF n = 0 OR n MOD 2 = 1 THEN ERROR;
n _ n / 2;
k _ k + 1;
ENDLOOP;
};
lgnsrq: NAT;
n: NAT;
bytes: LONG CARDINAL _ IO.GetLength[stream];
IF bytes MOD elementSize # 0 THEN ERROR;
lgnsrq _ Lg[bytes/elementSize];
IF lgnsrq MOD 2 # 0 THEN ERROR;
n _ Basics.BITSHIFT[1, lgnsrq/2];
FOR i: NAT IN [0..skipDots+dotCount) DO
inrange: BOOLEAN _ i IN [skipLines..skipLines+lineCount);
IF inrange THEN cmd.out.PutRope["["];
FOR j: NAT IN [0..n) DO
complex: Complex.VEC;
TRUSTED {
dataPointer: LONG POINTER _ @complex;
IF elementSize # IO.UnsafeGetBlock[stream, [base: dataPointer, count: elementSize]] THEN ERROR;
};
IF inrange AND j IN [skipDots..skipDots+dotCount) THEN cmd.out.PutF["\t%g", IO.real[IF realPart THEN complex.x ELSE complex.y]];
IF inrange AND j IN [skipDots..skipDots+dotCount-1) THEN cmd.out.PutRope[","];
ENDLOOP;
IF inrange THEN cmd.out.PutRope["]"];
IF i IN [skipLines..skipLines+lineCount-1) THEN cmd.out.PutRope[",\n"];
ENDLOOP;
END;
IO.Close[stream];
};

FFTRankPixelsCommand: Commander.CommandProc ~ {
stream: IO.STREAM _ IO.RIS[cmd.commandLine];
inputName: ROPE _ GetToken[stream];
skipLines: NAT _ GetInt[stream];
skipDots: NAT _ GetInt[stream];
lineCount: NAT _ GetInt[stream];
dotCount: NAT _ GetInt[stream];
elementSize: NAT ~ Basics.bytesPerWord*SIZE[Complex.VEC];
element: REF ARRAY [0..32000) OF REAL _ NEW[ARRAY [0..32000) OF REAL];
position: REF ARRAY [0..32000) OF REAL _ NEW[ARRAY [0..32000) OF REAL];
k: NAT _ 0;
stream _ FS.StreamOpen[fileName: inputName, accessOptions: $read, streamOptions: [tiogaRead: FALSE, commitAndReopenTransOnFlush: TRUE, truncatePagesOnClose: TRUE, finishTransOnClose: TRUE, closeFSOpenFileOnClose: TRUE]];
BEGIN ENABLE UNWIND => {IO.Close[stream]};
Lg: PROC [n: LONG CARDINAL] RETURNS [k: NAT _ 0] ~ {
WHILE n # 1 DO
IF n = 0 OR n MOD 2 = 1 THEN ERROR;
n _ n / 2;
k _ k + 1;
ENDLOOP;
};
lgnsrq: NAT;
n: NAT;
bytes: LONG CARDINAL _ IO.GetLength[stream];
IF bytes MOD elementSize # 0 THEN ERROR;
lgnsrq _ Lg[bytes/elementSize];
IF lgnsrq MOD 2 # 0 THEN ERROR;
n _ Basics.BITSHIFT[1, lgnsrq/2];
FOR i: NAT IN [0..skipDots+dotCount) DO
inrange: BOOLEAN _ i IN [skipLines..skipLines+lineCount);
FOR j: NAT IN [0..n) DO
complex: Complex.VEC;
TRUSTED {
dataPointer: LONG POINTER _ @complex;
IF elementSize # IO.UnsafeGetBlock[stream, [base: dataPointer, count: elementSize]] THEN ERROR;
};
IF inrange AND j IN [skipDots..skipDots+dotCount) THEN {
position[k] _ k;
element[k] _ complex.x;
k _ k + 1;
};
ENDLOOP;
ENDLOOP;
END;
IO.Close[stream];
ShellSort[element, position, k];
FOR i: NAT IN [0..k) DO
element[i] _ REAL[i]/k;
ENDLOOP;
ShellSort[position, element, k];
FOR i: NAT IN [0..k) DO
cmd.out.PutF["\t%g", IO.real[element[i]]];
IF (i+1) MOD dotCount = 0 THEN cmd.out.PutRope["\n"];
ENDLOOP;
};

shellD: ARRAY [1..11] OF CARDINAL = [1, 4, 13, 40, 121, 364, 1093, 3280, 9841, 29524, 65535];
ShellSort: PROCEDURE [key, data: REF ARRAY [0..32000) OF REAL, length: NAT] = {
passes: NAT _ 1;
UNTIL shellD[passes+2] >= length DO passes _ passes+1 ENDLOOP;
FOR pass: NAT DECREASING IN [1..passes] DO
h: NAT _ shellD[pass];
FOR i: NAT IN [0..length) DO
aKey: REAL _ key[i];
aData: REAL _ data[i];
j: NAT _ i;
WHILE j>=h AND aKey < key[j-h] DO
key[j] _ key[j-h];
data[j] _ data[j-h];
j _ j-h;
ENDLOOP;
key[j] _ aKey;
data[j] _ aData;
ENDLOOP;
ENDLOOP;
};

FFTLocalScaleCommand: Commander.CommandProc ~ {
stream: IO.STREAM _ IO.RIS[cmd.commandLine];
outputName: ROPE _ GetToken[stream];
gets: ROPE _ GetToken[stream];
inputName: ROPE _ GetToken[stream];
neighborRadius: INT _ GetInt[stream];
radiusSqr: INT _ neighborRadius*neighborRadius;
image: ImageFFT.Image _ ImageFFT.Load[inputName];
size: INT _ image.Height;
pixelMap: ImagerPixelMap.PixelMap _ ImagerPixelMap.Create[3, [0, 0, size, size]];
a: ARRAY [0..1024) OF Seq.ComplexSequence;
n: NAT _ 0;
Sqr: PROC [a: NAT] RETURNS [INT] ~ INLINE {RETURN [Basics.LongMult[a, a]]};
FOR p: ImageFFT.Image _ image, p.rest UNTIL p=NIL DO
a[n] _ p.first;
n _ n + 1;
ENDLOOP;
FOR i: INT IN [0..n) DO
FOR j: INT IN [0..n) DO
aij: REAL _ a[i][j].x;
max: REAL _ aij;
min: REAL _ aij;
FOR di: INT IN [-neighborRadius..neighborRadius] DO
FOR dj: INT IN [-neighborRadius..neighborRadius] DO
IF Sqr[ABS[di]] + Sqr[ABS[dj]] <= radiusSqr THEN {
aiijj: REAL;
ii: INT _ i + di;
jj: INT _ j + dj;
WHILE ii >= n DO ii _ ii - n ENDLOOP;
WHILE jj >= n DO jj _ jj - n ENDLOOP;
WHILE ii < 0 DO ii _ ii + n ENDLOOP;
WHILE jj < 0 DO jj _ jj + n ENDLOOP;
aiijj _ a[ii][jj].x;
min _ MIN[min, aiijj];
max _ MAX[max, aiijj];
};
ENDLOOP;
ENDLOOP;
IF min < max THEN aij _ (aij-min)/(max-min);
pixelMap.Fill[[i, j, 1, 1], MIN[MAX[Real.RoundLI[255*aij], 0], 255]];
ENDLOOP;
ENDLOOP;
StorePixelMap[outputName, pixelMap, FALSE, Rope.Concat["FFTLocalScale ", cmd.commandLine]];
FOR i: NAT IN [0..n) DO
a[i] _ NIL;
ENDLOOP;
image _ image.Destroy;
};

PixelMapFromFrameBuffer: PROC [frameBuffer: Terminal.FrameBuffer] 
RETURNS [ImagerPixelMap.PixelMap] ~ {
refRep: REF ImagerPixelMap.PixelMapRep ~ NEW[ImagerPixelMap.PixelMapRep _ [
ref: frameBuffer,
pointer: frameBuffer.base,
words: INT[frameBuffer.wordsPerLine]*INT[frameBuffer.height],
lgBitsPerPixel: BitsToLog[frameBuffer.bitsPerPixel],
rast: frameBuffer.wordsPerLine,
lines: frameBuffer.height
]];
frame: ImagerPixelMap.PixelMap ~ [
sOrigin: 0, fOrigin: 0, sMin: 0, fMin: 0,
sSize: frameBuffer.height, fSize: frameBuffer.width, refRep: refRep
];
RETURN[frame];
};


maxParam: NAT ~ 4;
filterList: LIST OF FilterEntry ~ LIST [
["RoundLowPass", 1, RoundLowPassFilter, "Round-shaped low-pass filter"],
["Spot", 4, SpotFilter, "[x, y, radius, m]  scales fuzzy spot around [x, y] by m"],
["BandBoost", 4, BandBoostFilter, "[innerRadius, outerRadius, edgeFuzz, m]  scales annulus centered on the origin by m"]
];

Commander.Register["ClipAIS", ClipAISCommand, "Clip a sampled image (output _ input <skipLines> <skipDots> <lineCount> <dotCount>)"];

Commander.Register["FFTFilter", FFTFilterCommand, "Apply spectral filters to (output _ input <params1> filter1 <params2> filter2 . . . )"];

Commander.Register["FFTShow", FFTShowCommand, "Show the filtered FFT of an image on the color display (input <params1> filter1 <params2> filter2 . . . ) (log intensity is displayed)"];

Commander.Register["FFTPeek", FFTPeekCommand, "Type the pixel values of a portion of a .image file (input <REAL | IMAGINARY> <firstScan> <firstPixel> <scanCount> <pixelCount>)"];

Commander.Register["FFTRankPixels", FFTRankPixelsCommand, "Type the ranks of the pixel values of a portion of a .image file (input <firstScan> <firstPixel> <scanCount> <pixelCount>)"];

Commander.Register["FFTLocalScale", FFTLocalScaleCommand, "Rescale each pixel according to the min and max of its local neighborhood (output.ais _ input.image <neighborRadius>)"];

END.
TImageFFTCommandsImpl.mesa
Copyright (c) 1984, 1985, Xerox Corporation. All right reserved.
Michael Plass, July 19, 1984 9:33:16 pm PDT
Tim Diebert: August 12, 1985 3:33:25 pm PDT
row: ImagerPixelRow.PixelRow _ ImagerPixelRow.CreatePixelRow[height/2, width/2, width];
row[j] _ MIN[MAX[pixel,0],255];
row.sOrigin _ row.sOrigin - height;
row.fOrigin _ row.fOrigin - width;
ImagerPixelRow.StorePixelRow[row, cd];
row.fOrigin _ row.fOrigin + width;
ImagerPixelRow.StorePixelRow[row, cd];
row.sOrigin _ row.sOrigin + height;
row.fOrigin _ row.fOrigin - width;
ImagerPixelRow.StorePixelRow[row, cd];
row.fOrigin _ row.fOrigin + width;
ImagerPixelRow.StorePixelRow[row, cd];
row.sOrigin _ row.sOrigin + 1;
ShowAISCommand: Commander.CommandProc ~ {
stream: IO.STREAM _ IO.RIS[cmd.commandLine];
inputName: ROPE _ GetToken[stream];
param: ROPE _ GetToken[stream];
context: Imager.Context _ Imager.Create[$Gray8bpp];
color: ImagerBasic.SampledColor _ ImagerAISUtil.AISToColor[inputName];
horizontalSlotNumber: INT _ 0;
horizontalSlotNumber _ Convert.IntFromRope[param ! Convert.Error => CONTINUE];
context.state.T _ Imager.Scale[1.0];
context.SetColor[color];
context.MaskRectangle[horizontalSlotNumber*color.pa.yPixels, 0, color.pa.yPixels, color.pa.xPixels];
cmd.out.PutRope["AIS comment: "];
cmd.out.PutRope[GetAISComment[inputName]];
cmd.out.PutRope["\n"];
cmd.out.PutRope["Done.\n"];
};

Commander.Register["ShowAIS", ShowAISCommand, "Show a sampled image on the color display (input <horizontalSlotNumber: INTEGER>)"];

Commander.Register["FFTHelp", FFTHelpCommand, "List the names of the available spectral filters"];

ΚΡ˜™J™@J™+Icode™+J˜—šΟk	œ
œ2œ8œ.˜©
J˜—šœ
˜#š	œ
œ'œ0œ+˜–
J˜—Jšœœœ
˜J˜šΟnœœœœ&œœ˜gJšœœœ˜#Jšœœœ˜'Jš
œ!˜!Jšœœ˜Jšœœ˜Jš
œΞ
˜Ξ
Jš
œ‹˜‹Jš
œN˜Nšœœ
œE˜TJšœ+˜.Jš
œ˜Jš
œ&˜&Jšœ
˜—Jš
œ˜J˜—š	ž	œœœœ˜@šœœ˜Jšœ+œœ
˜<—Jšœ
˜J˜—šž
œœœœœœ˜KJšœœ˜Jšœœ
˜Jšœœ
˜Jš	œœ
œœ
œ˜Jšœ+œ˜DJšœœ˜!Jšœœ8œœ˜nJšœ$œ˜*Jš
œ˜J˜—šž
œœœ+œœ
œœ˜
š	œœœ'œ
œ˜SJš
œ˜Jš
œ˜Jš
œ
˜
Jšœœ"œœœœ"˜zJš
œ˜Jš
œ˜Jš
œ˜—Jš
œ0˜0Jš
œ’
˜’
JšœœG˜ZJšœ˜"šœœ
œ˜"Jš
œ˜Jš
œ˜Jš
œ&˜&Jšœ2˜5Jš
˜—Jšœ˜Jš
œ˜J˜—š
ž
œœœœœ˜CJšœœœ˜0Jšœ
œ˜Jšœ˜Jš
œ˜J˜—–

"cedar" stylešžœœ%œ&˜^Jšœœœ!˜-Jšœœ˜Jšœœ	œ
˜Jš
œM˜MJš
œ
˜
Jš
œ˜Jšœ$œ+˜UJš
œ˜J˜—–

"cedar" stylešžœœ
œœ
œœœ˜WJšœ	œ˜)Jš
œC˜CJš
œ3˜3Jš
œL˜LJšœœ
˜Jšœ œ
œ˜6šœ	œ
œ˜Jš
œ&˜&Jšœœ
œ)œ˜NJšœœ˜!Jš
œ˜Jšœœ
œœ˜>Jš
œ˜—Jšœ˜Jš
œ˜J˜—–

"cedar" stylešž	œœœœ
œœ
œœœ˜rJšœ	œ˜)Jš
œC˜CJš
œ3˜3Jš
œL˜LJšœœ
˜Jšœœ$œ˜6Jš	œœœ
œœ)˜oJš
œ6˜6Jšœ œ
œ˜6šœ	œ
œ˜Jš
œ&˜&Jšœœ
œ)œ˜NJšœœ˜!Jš
œ˜Jšœœ
œœ˜>Jš
œ˜—Jšœœ
œ!˜2Jš
œ˜Jšœœ
œ2˜CJšœœ˜ Jš
œ.˜.Jšœœ
œœ˜DJš
œ!˜!Jšœœ
œœ˜EJš
œ ˜ Jšœ"œ˜2Jš
œ˜J˜—–

"cedar" styleš
žœœœœœ˜3Jšœœœ
œ	˜.Jšœœœœœ
œ˜HJšœ	˜Jš
œ˜J˜—–

"cedar" stylešžœœ
œ
œœœœ˜AJšœ#œœ˜EJš
œ˜J˜—–

"cedar" stylešžœœ
œ
œœœ˜7Jšœœ
œ
˜Jšœœœ
œ˜<Jš
œ˜J˜—š
žœœœœ+œ˜bJš
œV˜VJš
œ+˜+Jš
œ€
˜€
Jš
œ˜J˜—šžœœ
œœœœ˜9Jšœœœ˜Jšœœœ˜Jšœœ
œ˜'Jš
œ˜J˜—Jšœœ˜Jšœœ˜šžœœ˜"Jšœœ˜!Jšœœ˜Jš
œ:˜:Jš
œW™WJšœœ˜š
žœœœœœœ
˜4Jšœœ
œ˜7Jšœ
˜—šœœœ
˜0šœœ
œ˜Jšœœ˜!Jš	œœœœœ˜7Jšœœ)˜3Jšœ	œ
œ™Jšœ)œ
œ˜@Jšœ
˜—Jšœœœ˜7Jš
œ#™#Jš
œ"™"Jš
œ&™&Jš
œ"™"Jš
œ&™&Jš
œ#™#Jš
œ"™"Jš
œ&™&Jš
œ"™"Jš
œ&™&Jš
œ™Jšœ
˜—Jš
œ˜J˜—šžœ˜"Jšœœœ
œœœ˜KJ˜J˜—šž
œ˜Jšœœ˜Jšœœ˜Jšœœ˜Jšœœœ˜+Jšœœœ
œ˜-šœ˜Jšœœ˜"Jšœœ
œ˜J˜J˜J˜Jšœ˜
Jš
œ˜—J˜J˜—šžœœœœœœœ˜;J˜—šžœ˜Jšœœ˜Jšœœ˜Jšœœ˜Jšœœœ
œ˜HJšœœœ
œ˜GJš
œ˜Jšœœ˜3Jšœ˜J˜Jšœ˜J˜J˜—šž
œœœ	œœ
œ
œœœ˜\J˜—šœœœ˜Jš
œ˜Jš	œœœ
œœ˜+Jš
œ˜J˜—šœ
œœ˜Jšœœ
˜Jšœœ˜Jš
œ˜Jšœ˜	Jš
œ˜J˜—šžœœœœ˜Išœœ
œ"œœ
˜@Jšœ œœ
œ
˜>Jšœ
˜—Jšœœœœ˜Jš
œ˜J˜—šžœœ
œ
œœ
œœœ
œ˜bJšœœ
œ
œ
˜Jšœœ˜Jš	œœœ
œœ˜,Jšœœ˜šœ˜Jšœœ
˜Jšœ2œ˜<Jšœœ2œ˜Qšœœ˜J˜2šœœ
œ˜&Jšœœ˜/Jšœ
œ
˜Jš
œ˜—šœœ˜%Jšœ+œ˜;Jšœ
œ
˜Jš
œ˜—Jšœœ%˜-Jš
œ˜Jšœœ˜Jš
œ
˜
—šœ˜Jš
œ˜Jš
œ˜šœœ˜Jšœ$œ˜4Jšœ
œ
˜Jš
œ˜—Jš
œ˜—Jš
œ˜Jšœ
˜—šœœ
˜Jšœœ
œ˜J˜Jš
œ˜Jš
œ˜Jšœ
˜—Jš
œ˜J˜—šžœœ"œ
œ˜JJšœœ˜ Jšœœ	˜Jšœœ˜šœ#œœ
˜6J˜#šœœ
œ˜Jš	œœœ	œœ˜&Jšœœœœ˜/Jšœ
œ
˜Jšœœ˜šœœ
œœœ
˜;J˜+Jš
˜—Jš
œ˜Jšœ
˜—J˜Jšœ	œ˜Jš
˜—Jš
œ˜J˜—šžœ˜)Jš	œœ
œœ
œ˜,Jšœœ˜#Jšœ	œ
œ.˜>Jš
œ9˜9Jš
œ˜Jš
œ˜Jš
œ˜Jš
œ˜J˜—šžœ˜+Jš	œœ
œœ
œ˜,Jšœœ˜$Jšœœ˜Jšœœ˜#Jš
œ6˜6Jšœ	œ
œ.˜>šœ
œœ˜Jš
œB˜BJšœ
˜Jš
œ˜—Jš
œ0˜0JšœEœ(˜sJš
œ!˜!Jš
œ˜Jš
œ˜J˜—šžœ™)Jš	œœ
œœ
œ™,Jšœœ™#Jšœœ™J™3Jš
œF™FJšœœ™JšœDœ™NJšœ
œ™$Jš
œ™Jš
œd™dJš
œ!™!Jš
œ*™*Jš
œ™Jš
œ™Jš
œ™J™—šžœ˜)Jš	œœ
œœ
œ˜,Jšœœ˜$Jšœœ˜Jšœœ˜#Jšœœ˜ Jšœ
œ˜Jšœœ˜ Jšœ
œ˜Jš
œ"˜"Jšœœ˜šœ
œœœœœœ˜hJš
œc˜cJšœ
˜Jš
œ˜—Jš
œ0˜0Jš
ϣ
˜£
Jš
œ˜Jš
œ˜J˜—šžœ˜)Jš	œœ
œœ
œ˜,Jšœœ˜#Jšœœ˜Jšœœ˜ Jšœ
œ˜Jšœœ˜ Jšœ
œ˜Jšœ
œœœ˜KJšœ
œœ	œ˜9Jšœ]œœœœœ˜ά
šœ
œ
œœ˜*š
žœœ
œœœ	˜4šœ˜Jš
œœœœ
œ
˜#J˜
J˜
Jš
˜—Jš
œ˜—Jšœœ
˜Jšœœ
˜Jšœœ
œœ˜,Jšœœœ
œ
˜(Jš
œ˜Jšœœœ
œ
˜Jšœœ˜!šœœ
œ˜'Jšœ	œœ"˜9Jšœ	œ˜%šœœ
œ˜Jšœœ
˜šœ˜	Jšœ
œ
œ˜%JšœœAœ
œ
˜_Jš
œ˜—Jšœ	œœœœœ
œœ
˜€
Jšœ	œœ!œ˜NJšœ
˜—Jšœ	œ˜%Jšœœ$œ˜GJš
˜—Jšœ
˜—Jšœ˜Jš
œ˜J˜—šžœ˜/Jš	œœ
œœ
œ˜,Jšœœ˜#Jšœœ˜ Jšœ
œ˜Jšœœ˜ Jšœ
œ˜Jšœ
œœ	œ˜9Jšœ	œ
œœ
œœ
œœ
œ˜FJšœ
œ
œœ
œœ
œœ
œ˜GJšœœ˜Jšœ]œœœœœ˜ά
šœ
œ
œœ˜*š
žœœ
œœœ	˜4šœ˜Jš
œœœœ
œ
˜#J˜
J˜
Jš
˜—Jš
œ˜—Jšœœ
˜Jšœœ
˜Jšœœ
œœ˜,Jšœœœ
œ
˜(Jš
œ˜Jšœœœ
œ
˜Jšœœ˜!šœœ
œ˜'Jšœ	œœ"˜9šœœ
œ˜Jšœœ
˜šœ˜	Jšœ
œ
œ˜%JšœœAœ
œ
˜_Jš
œ˜—šœ	œœœ˜8Jš
œ˜Jš
œ˜J˜
Jš
œ˜—Jšœ
˜—Jšœ
˜—Jšœ
˜—Jšœ˜Jš
œ ˜ šœœ
œ˜Jšœ
œ˜Jšœ
˜—Jš
œ ˜ šœœ
œ˜Jšœœ˜*Jšœœœ˜5Jšœ
˜—Jš
œ˜J˜—Iunitšœœ	œ
œ<˜]šž	œ	œ
œ
œœ
œ
œ˜OJšœœ˜Jšœœœ
˜>š	œœ

œ
œ
˜*Jšœœ˜šœœ
œ
˜Jšœœ
˜Jšœœ˜Jšœœ˜šœœ˜!Jš
œ˜Jš
œ˜Jš
œ˜Jšœ
˜—Jš
œ˜Jš
œ˜Jšœ
˜—Jšœ
˜—Jš
œ˜J˜—šžœ˜/Jš	œœ
œœ
œ˜,Jšœœ˜$Jšœœ˜Jšœœ˜#Jšœœ˜%Jšœœ!˜/Jš
œ1˜1Jšœœ˜Jš
œQ˜QJšœœœ˜*J˜Jšžœœœœœœœ˜Kšœ#œœ
˜4Jš
œ˜Jš
œ
˜
Jš
˜—šœœ
œ˜šœœ
œ˜J˜J˜J˜šœœ
œ#˜3šœœ
œ#˜3šœœœœ˜2Jšœœ
˜Jšœœ
˜Jšœœ
˜Jšœ	œ
œ
˜%Jšœ	œ
œ
˜%Jšœœ
œ
˜$Jšœœ
œ
˜$Jš
œ˜Jš
œ˜Jš
œ˜Jš
œ˜—Jš
˜—Jš
˜—Jš
œ,˜,Jšœœ
œ"˜EJš
˜—Jš
˜—Jšœ$œ2˜[šœœ
œ˜Jšœœ
˜Jš
˜—Jš
œ˜Jš
œ˜—K˜šžœœ%œ˜hšœœœ˜KKš
œ˜Kš
œ˜Kšœœœ˜=Kš
œ4˜4Kš
œ˜Kš
œ˜Kš
œ˜—š
œ"˜"Kš
œ)˜)Kš
œC˜CKš
œ˜—Kšœ˜K˜K˜—J˜Jšœ
œ˜šœœ
œœ˜(Jš
œH˜HJš
œS˜SJš
œx˜xJš
œ˜J˜—š
œƒ
™ƒ
J™—š
œ…
˜…
J˜—š
œb™bJ™—š
œ‹
˜‹
J˜—š
œΈ
˜Έ
J˜—š
œ²
˜²
J˜—š
œΈ
˜Έ
J˜—š
œ³
˜³
J˜—Jšœ
˜——…—Mn9