DIRECTORY AISFileFormat, Args, CtBasic, CtDispatch, CtFilter, ImagerSample, IO, RealFns, Rope;

CtFilterCommandsImpl: CEDAR PROGRAM

IMPORTS Args, CtBasic, CtDispatch, CtFilter, IO, RealFns, Rope

~ BEGIN
ROPE:					TYPE ~ Rope.ROPE;
STREAM:				TYPE ~ IO.STREAM;
Box:					TYPE ~ ImagerSample.Box;
SampleMaps:			TYPE ~ CtBasic.SampleMaps;
Filter:					TYPE ~ CtFilter.Filter;
FilterEvalProc:		TYPE ~ CtFilter.FilterEvalProc;
FilterInitProc:		TYPE ~ CtFilter.FilterInitProc;
FilterPrintProc:		TYPE ~ CtFilter.FilterPrintProc;
CtExpand: CtDispatch.CtProc ~ {
sxA, syA, swA, shA, dxA, dyA, dwA, dhA: Args.Arg;
[sxA, syA, swA, shA, dxA, dyA, dwA, dhA] ¬ Args.ArgsGet[cmd, "-s%iiii-d%iiii" !
Args.Error => {error ¬ reason; CONTINUE}];
IF error # NIL
THEN RETURN
ELSE {
sx: NAT ¬ IF sxA.ok	THEN sxA.int ELSE 0;
sy: NAT ¬ IF syA.ok THEN syA.int ELSE 0;
sw: NAT ¬ IF swA.ok THEN swA.int ELSE 640;
sh: NAT ¬ IF shA.ok THEN shA.int ELSE 480;
dx: NAT ¬ IF dxA.ok THEN dxA.int ELSE 0;
dy: NAT ¬ IF dyA.ok THEN dyA.int ELSE 0;
dw: NAT ¬ IF dwA.ok THEN dwA.int ELSE 1024;
dh: NAT ¬ IF dhA.ok THEN dhA.int ELSE 768;
srcBox: Box ¬ [[sy, sx], [sy+sh, sx+sw]];
dstBox: Box ¬ [[dy, dx], [dy+dh, dx+dw]];
IF maps.bpp # 8
THEN RETURN[error: "this works only for 8 bpp, try Ct Resample"]
ELSE CtFilter.Expand[maps[0].map, srcBox, dstBox];
affect ¬ dstBox;
};
};
CtResample: CtDispatch.CtProc ~ {
srcXA, srcYA, srcWA, srcHA, dstXA, dstYA, dstWA, dstHA, sxA, syA, txA, tyA: Args.Arg;
filtNameA, supportA, blurA, par1A, par2A, windowNameA: Args.Arg;
[srcXA, srcYA, srcWA, srcHA, dstXA, dstYA, dstWA, dstHA, sxA, syA, txA, tyA,
 filtNameA, supportA, blurA, par1A, par2A, windowNameA] ¬ Args.ArgsGet[
cmd, "-s%iiii-d%iiii-map%rrrr-filt%s-supp%r-blur%r-param%r[r-window%s"
! Args.Error => {error ¬ reason; GOTO Bad}];
error ¬ Resample[
maps,
srcXA.int, srcYA.int, srcWA.int, srcHA.int, dstXA.int, dstYA.int, dstWA.int, dstHA.int,
supportA.real, blurA.real, par1A.real, par2A.real, sxA.real, syA.real, txA.real, tyA.real,
IF filtNameA.ok THEN filtNameA.rope ELSE "triangle",
windowNameA.rope,
cmd.out];
affect ¬ [[dstYA.int, dstXA.int], [dstYA.int+dstHA.int, dstXA.int+dstWA.int]];
CtBasic.ClipMaps[maps, dstXA.int, dstYA.int, dstWA.int, dstHA.int];
EXITS Bad => RETURN;
};

Resample: PROC [
maps: SampleMaps,
srcX, srcY, srcW, srcH, dstX, dstY, dstW, dstH: NAT ¬ 0,
support, blur, par1, par2, sx, sy, tx, ty: REAL ¬ 0.0,
filtName: ROPE ¬ "triangle",
windowName: ROPE ¬ NIL,
out: STREAM]
RETURNS [err: ROPE]
~ {
IF Rope.Equal[filtName, "?"]
THEN CtFilter.FilterCatalog[out]
ELSE {
src, dst: Box;
filt: Filter ¬ CtFilter.FilterFind[filtName ! CtFilter.Error => {err ¬ reason; GOTO NoFilt}];
IF srcX = 0 THEN srcX ¬ maps.x;
IF srcY = 0 THEN srcY ¬ maps.y;
IF srcW = 0 THEN srcW ¬ maps.w;
IF srcH = 0 THEN srcH ¬ maps.h;
IF dstX = 0 THEN dstX ¬ maps.x;
IF dstY = 0 THEN dstY ¬ maps.y;
IF dstW = 0 THEN dstW ¬ maps.w;
IF dstH = 0 THEN dstH ¬ maps.h;
src ¬ [[srcY, srcX], [srcY+srcH, srcX+srcW]];
dst ¬ [[dstY, dstX], [dstY+dstH, dstX+dstW]];
IF windowName = NIL THEN windowName ¬
IF filt.windowMe THEN "blackman" ELSE "box";	-- else use blackman if IIR
IF support # 0.0 THEN filt.support ¬ support;
IF blur # 0.0 THEN filt.blur ¬ blur;
IF par1 # 0.0 AND par2 # 0.0 AND filt.initProc # NIL
THEN filt.clientData ¬ filt.initProc[par1, par2];
IF filt.printProc # NIL THEN filt.printProc[out, filt.clientData];
IF NOT Rope.Equal[windowName, "box", FALSE] THEN {	-- window the filter
window: Filter ¬ CtFilter.FilterFind[windowName
! CtFilter.Error => {err ¬ reason; GOTO NoFilt}];
IF par1 # 0.0 AND par2 # 0.0 AND filt.initProc = NIL AND window.initProc # NIL
THEN window.clientData ¬ window.initProc[par1, par2];
IF window.printProc # NIL THEN window.printProc[out, window.clientData];
filt.name ¬ Rope.Cat[filt.name, "*", window.name];
filt.clientData ¬ WindowInit[filt, window];
filt.initProc ¬ NIL;
filt.printProc ¬ NIL;
filt.evalProc ¬ WindowEval;	-- WindowEval[x] = filt.eval[x]*window.eval[x/s]
};
IF sx # 0.0 AND tx # 0.0 AND sy # 0.0 AND ty # 0.0
THEN CtFilter.ResampleContinuous[maps, src, dst, [sx, sy], [tx, ty], filt, out]
ELSE CtFilter.Resample[maps, src, dst, filt, out];
};
EXITS
NoFilt => RETURN[err];
};
pi: REAL ~ 3.1415927;

MitchellData: TYPE ~ RECORD [
p0, p2, p3, q0, q1, q2, q3: REAL
];

MitchellEval: FilterEvalProc ~ {
d: REF MitchellData ¬ NARROW[clientData];
SELECT x FROM
<-2. => RETURN[0.];
<-1. => RETURN[d.q0-x*(d.q1-x*(d.q2-x*d.q3))];
< 0. => RETURN[d.p0+x*x*(d.p2-x*d.p3)];
< 1. => RETURN[d.p0+x*x*(d.p2+x*d.p3)];
< 2. => RETURN[d.q0+x*(d.q1+x*(d.q2+x*d.q3))];
ENDCASE => RETURN[0.];
};

MitchellInit: FilterInitProc ~ {
c ¬ b;		-- use same parameter names as Mitchell and Netravali
b ¬ a;
RETURN[NEW[MitchellData ¬ [
p0: (6.-2.*b)/6.,
p2:  (-18.+12.*b+6.*c)/6.,
p3: (12.-9.*b-6.*c)/6.,
q0: (8.*b+24.*c)/6.,
q1: (-12.*b-48.*c)/6.,
q2: (6.*b+30.*c)/6.,
q3: (-b-6.*c)/6.]]];
};

MitchellPrint: FilterPrintProc ~ {
d: REF MitchellData ¬ NARROW[clientData];
IO.PutF[s, "\tp0=%g p2=%g p3=%g\n", IO.real[d.p0], IO.real[d.p2], IO.real[d.p3]];
IO.PutFL[s, "\tq0=%g q1=%g q2=%g q3=%g\n",
LIST[IO.real[d.q0], IO.real[d.q1], IO.real[d.q2], IO.real[d.q3]]];
};

HanningWindow: FilterEvalProc ~ {RETURN[.5+.5*RealFns.Cos[pi*x]];};

HammingWindow: FilterEvalProc ~ {RETURN[.54+.46*RealFns.Cos[pi*x]];};

BlackmanWindow: FilterEvalProc ~ {RETURN[.42+.50*RealFns.Cos[pi*x]+.08*RealFns.Cos[2.*pi*x]];};


KaiserData: TYPE ~ RECORD [
a: REAL,				-- = w*(N-1)/2 in Oppenheim & Schafer notation, typically 4<a<9
i0a: REAL		
];

KaiserWindow: FilterEvalProc ~ {
d: REF KaiserData ¬ NARROW[clientData];
RETURN[I0[d.a*RealFns.SqRt[1.-x*x]]*d.i0a];
};

KaiserInit: FilterInitProc ~ {
RETURN[NEW[KaiserData ¬ [a, 1./I0[a]]]];
};

KaiserPrint: FilterPrintProc ~ {
d: REF KaiserData ¬ NARROW[clientData];
IO.PutF[s, "\ta=%g 1/I0[a]=%g\n", IO.real[d.a], IO.real[d.i0a]];
};

I0: PROC [x: REAL] RETURNS [REAL] ~ {
sum: REAL ¬ 1.;
y: REAL ¬ x*x/4.;
t: REAL ¬ y;
i: NAT ¬ 2;
epsilon: REAL ~ 1e-7;
WHILE t > epsilon DO
sum ¬ sum+t;
t ¬ t*y/(i*i);
i ¬ i+1;
ENDLOOP;
RETURN[sum];
};
WindowData: TYPE ~ RECORD [
support: REAL,						-- desired width of scaled window
filterEvalProc: FilterEvalProc,		-- filter function (presumably IIR)
windowEvalProc: FilterEvalProc,	-- window function (presumably FIR with support 1)
filtData, windowData: REF ANY	-- pass along the client data to each EvalProc
];

WindowEval: FilterEvalProc ~ {
d: REF WindowData ¬ NARROW[clientData];
RETURN[d.filterEvalProc[x, d.filtData]*d.windowEvalProc[x/d.support, d.windowData]];
};

WindowInit: PROC [filt, window: Filter] RETURNS [REF ANY] ~ {
RETURN[NEW[WindowData ¬
[filt.support, filt.evalProc, window.evalProc, filt.clientData, window.clientData]]];
};
ctExpandUsage:		ROPE ~ "Ct Expand [-s<xywh>] [-d<xywh>]: expand s box to d box.";

ctResampleUsage:	ROPE ~ "Ct Resample [options]
Resample src window to dst window using the specified filter.
The src and dst boxes can overlap!
Options:
	-s <x y w h>		source window
	-d <x y w h>	destination window
	-map <xscale yscale xtran ytran>
						optional src to dst mapping:
						dstx = xscale*srcx+xtran, similar for y.
						Mapping is indep. of src and dst windows.
						Default mapping stretches src to dst.
	-filt <filtName>	set filter name (default=triangle)
						'ct resample -filt ?' prints filter catalog
	-supp <width>	filter support (width), defaults in catalog
	-blur <factor>	blur factor: > 1 is blurry, < 1 is aliased
	-window <windowName>
						window an IIR (infinite) filter.
						Common windows are triangle, hanning,
						hamming, blackman, and kaiser.  IIR filters
						are windowed with blackman by default;
						windowing may be disabled with
						'-window box'.
	-param <a b>		set misc filter parameters
";


CtFilter.FilterRegister["mitchell", MitchellEval, 2.0, 1.0, FALSE, FALSE, FALSE, MitchellInit, MitchellPrint, 1.0, MitchellInit[1./3., 1./3.]];	-- preferred Mitchell Filter

CtFilter.FilterRegister["hanning",  HanningWindow,	1.0, 1.0,	FALSE, TRUE, TRUE];
CtFilter.FilterRegister["hamming", HammingWindow,	1.0, 1.0,	FALSE, TRUE, TRUE];
CtFilter.FilterRegister["blackman",	BlackmanWindow,	1.0, 1.0,	FALSE, TRUE, TRUE];
CtFilter.FilterRegister["kaiser",		KaiserWindow,		1.0, 1.0,	FALSE, TRUE, TRUE, KaiserInit, KaiserPrint, 1.0, KaiserInit[6.5]];

CtDispatch.RegisterCtOp["Resampling and Filtering:",	NIL,		NIL];
CtDispatch.RegisterCtOp["Expand",		CtExpand,				ctExpandUsage];
CtDispatch.RegisterCtOp["Resample",		CtResample,				ctResampleUsage];

END.

l
CtFilterCommandsImpl.mesa
Copyright Σ 1988, 1992 by Xerox Corporation.  All rights reserved.
written by Paul Heckbert, August 26, 1988 1:33:39 pm PDT
Bloomenthal, July 5, 1992 9:04 pm PDT
Imported Types
Expand Command
Resample Command
CtResampleFile: CtDispatch.CtProc ~ {
FileBox: PROC [filename: ROPE] RETURNS [box: Box] ~ {
V: PROC [c: AISFileFormat.Card16] RETURNS [CARD] ~ {RETURN[Basics.Card16FromH[c]]};
i: AISIO.Info _ AISIO.ReadInfo[filename];
box _ [[0, 0], [V[i.raster.scanCount], V[i.raster.scanLength]]];
};
srcFileA, dstFileA, srcXA, srcYA, srcWA, srcHA, dstWA, dstHA: Args.Arg;
sxA, syA, txA, tyA, filtNameA, supportA, blurA, par1A, par2A, windowNameA: Args.Arg;
[srcFileA, dstFileA, srcXA, srcYA, srcWA, srcHA, dstWA, dstHA, sxA, syA, txA, tyA,
 filtNameA, supportA, blurA, par1A, par2A, windowNameA] _ Args.ArgsGet[
cmd, "%ss-s%iiii-d%ii-map%rrrr-filt%s-supp%r-blur%r-param%r[r-window%s"
! Args.Error => {msg _ reason; CONTINUE}];
IF msg # NIL THEN RETURN[$Failure, msg];
IF NOT dstWA.ok OR NOT dstHA.ok
THEN RETURN[$Failure, "missing destination size (-d)"]
ELSE {
maps: SampleMaps;
srcFilename: ROPE _ CtFile.Parse[srcFileA.rope].names[0];
type: CtFile.FileType _ CtFile.Parse[srcFileA.rope].type;
srcB: Box _ FileBox[srcFilename ! FS.Error => GOTO BadFile];
srcX: NAT _ IF srcXA.ok THEN srcXA.int ELSE srcB.min.f;
srcY: NAT _ IF srcYA.ok THEN srcYA.int ELSE srcB.min.s;
srcW: NAT _ IF srcWA.ok THEN srcWA.int ELSE srcB.max.f-srcB.min.f;
srcH: NAT _ IF srcHA.ok THEN srcHA.int ELSE srcB.max.s-srcB.min.s;
dstW: NAT _ dstWA.int;
dstH: NAT _ dstHA.int;
bpp: NAT _ SELECT type FROM bw => 8, color => 24, ENDCASE => 0;
IF bpp = 0 THEN RETURN[$Failure, "bad file type"];
srcB _ SF.Intersect[-- file: -- srcB, [[srcY, srcX], -- user: --[srcY+srcH, srcX+srcW]]];
Set maps to max possible size:
maps _ CtBasic.CreateMaps[bpp, 0, 0, MAX[srcB.max.f, dstW], MAX[srcB.max.s, dstH]];
[] _ CtFile.ReadFile[srcFileA.rope, srcX, srcY, srcW, srcH, maps];
msg _ Resample[
maps,
srcX, srcY, srcW, srcH, 0, 0, dstW, dstH,
supportA.real, blurA.real, par1A.real, par2A.real, sxA.real, syA.real, txA.real, tyA.real,
IF filtNameA.ok THEN filtNameA.rope ELSE "triangle",
windowNameA.rope,
cmd.out];
IF msg # NIL THEN RETURN[$Failure, msg];
CtBasic.ClipMaps[maps, 0, 0, dstW, dstH];
CtFile.WriteMapsToAIS[maps, dstFileA.rope, cmd.out];
EXITS BadFile => RETURN[$Failure, "bad AIS file"];
};
};

Additional Filters
two-parameter cubic filter from
Mitchell&Netravali "Reconstruction Filters in Computer Graphics", SIGGRAPH 88
Window Functions
Window functions are functions that are multiplied by IIR filters to make them FIR.  These are adapted from A.V. Oppenheim, R.W. Schafer, Digital Signal Processing, Prentice-Hall, 1975 scaled and translated to the domain [-1..1] by substituting n=(x+1)*(N-1)/2 into their formulas.
Note that these window functions do not integrate to 1, as do most of the filters in CtFilter.mesa.  That's not a problem for Resample, since it normalizes its filters.
param a trades off main lobe width (sharpness) for side lobe amplitude (ringing)
from Oppenheim & Schafer
Modified zeroth order Bessel function of the first kind.
Are there better ways to compute this than the power series?
Template For a Windowed Filter
a windowed filter function; this routine is used for all windowing
Start Code
ctResampleFileUsage:	ROPE ~ "Ct ResampleFile <srcFile> <dstFile> [options]
Resample from AIS file to AIS file, rather than in color display.
Options:
	-d <w h>			destination window (note: no origin)
	-s <x y w h>		source window (defaults to entire file)
	-map <xscale yscale xtran ytran>
						optional src to dst mapping:
						dstx = xscale*srcx+xtran, similar for y.
						Mapping is indep. of src and dst windows.
						Default mapping stretches src to dst.
	-filt <filtName>	set filter name (default=triangle)
						'ct resample -filt ?' prints filter catalog
	-supp <width>	filter support (width), defaults in catalog
	-blur <factor>	blur factor: > 1 is blurry, < 1 is aliased
	-window <windowName>
						window an IIR (infinite) filter.
						Common windows are triangle, hanning,
						hamming, blackman, and kaiser.  IIR filters
						are windowed with blackman by default;
						windowing may be disabled with
						'-window box'.
	-param <a b>		set misc filter parameters
";
						   name	 	eval	  			support blur window card. unitrange
CtDispatch.RegisterCtOp["ResampleFile",	CtResampleFile,			ctResampleFileUsage, FALSE];
Κ–
"cedarcode" style•NewlineDelimiter

™™Jšœ	Οeœ6™BJ™8J™%J˜—JšΟk	œCžœ˜^J˜šΠblœžœ
ž˜#J˜Jšžœ&žœ˜>—J˜Jšœž˜
headš
Οl™Jšžœžœžœ
˜Jšžœžœžœ
žœ
˜Jšœžœ˜!Jšœžœ˜(Jšœžœ˜#Jšœžœ˜0Jšœžœ˜0Jšœžœ˜2—š
 ™šΠbnœ˜J˜1˜OJšœžœ˜*—šžœ	ž˜Jšžœ
ž˜šžœ˜Jš	œžœžœžœ	žœ˜(Jš	œžœžœžœ	žœ˜(Jš	œžœžœžœ	žœ˜*Jš	œžœžœžœ	žœ˜*Jš	œžœžœžœ	žœ˜(Jš	œžœžœžœ	žœ˜(Jš	œžœžœžœ	žœ˜+Jš	œžœžœžœ	žœ˜*J˜)J˜)šžœ
˜Jšžœ
žœ5˜@Jšžœ.˜2—J˜J˜——J˜——š
 ™š‘
œ˜!J˜UJ˜@J˜L˜G˜FJšœ!žœ˜,——˜J˜J˜WJ˜ZJšžœžœžœ˜4Jš
œ˜J˜	—J˜NJ˜CJšžœžœ
˜˜J˜——š‘œ™%šΟnœžœžœžœ™5Jš’
œ
Οs
ž£
ϣ
ϣ
ž£
œ
žœžœ™SJšœžœžœ™)J™@J™—J™GJ™TJ™R™GJ™GJšœžœ™*—Jšžœžœ
žœ™(šžœ
žœ
žœ
žœ	™Jšžœ
žœ+™6šžœ™J™Jšœ
žœ%™9J™9Jšœ"žœ
žœ
™<Jš	œžœžœ	žœžœ™7Jš	œžœžœ
žœžœ™7Jš	œžœžœ
žœžœ™BJš	œžœžœ
žœžœ™BJšœžœ
™Jšœžœ
™Jš	œžœžœžœžœ™?Jšžœ	žœ
žœ™2JšœžœΟcœ€œ™YJ™Jšœ%žœžœ™SJ™B™J™J™)J™ZJšžœžœžœ™4J™J™	—Jšžœžœ
žœ
žœ™(J™)J™4Jšžœžœ™2J™——J™J™—š’œžœ˜Jš
œ˜Jšœ0žœ˜8Jšœ+žœ˜6Jšœ
žœ˜Jšœžœžœ
˜Jšœžœ
˜Jšžœžœ
˜Jš
œ˜šžœ˜Jšžœ˜ šžœ˜Jš
œ˜JšœOžœ
˜]Jšžœ
žœ˜Jšžœ
žœ˜Jšžœ
žœ˜Jšžœ
žœ˜Jšžœ
žœ˜Jšžœ
žœ˜Jšžœ
žœ˜Jšžœ
žœ˜J˜-J˜-šžœžœ
žœ
˜%Jšžœžœžœ€˜H—Jšžœžœ˜-Jšžœžœ˜$šžœžœžœž˜4Jšžœ-˜1—Jšžœžœ
žœ&˜Bš	žœ
žœžœžœ€˜G˜/Jšœ#žœ
˜1—šžœžœžœžœ
žœž˜NJšžœ1˜5—Jšžœžœ
žœ*˜HJ˜2J˜+Jšœžœ
˜Jšœžœ
˜Jšœ€0˜LJ˜—šžœ
žœ
žœ
žœ	˜2JšžœK˜OJšžœ.˜2—J˜——š
ž˜Jšœ
žœ˜—J˜——š
 ™šœžœ
˜J˜—šœžœžœ˜Jšœž˜ J˜J˜—š’œ˜ ™J™M—Jšœžœžœ
˜)šžœž˜
Jšœžœ˜Jšœžœ ˜.Jšœžœ˜'Jšœžœ˜'Jšœžœ ˜.Jšžœžœ˜—J˜J˜—š’œ˜ Jšœ€5˜=J˜šžœ
žœ˜J˜J˜J˜J˜J˜J˜J˜—J˜J˜—š’
œ˜"Jšœžœžœ
˜)Jšžœ"žœ
žœ
žœ
˜Qšžœ(˜*Jš
žœ
žœ
žœ
žœ
žœ˜B—J˜——š
 ™JšœŠ
Οzœv™™™¨
J˜—š’
œžœ˜CJ˜—š’
œžœ˜EJ˜—š’œžœ7˜_J˜J˜—šœžœžœ˜Jšœžœ€?˜KJ™PJšœž˜J˜J˜—š’œ˜ J™Jšœžœžœ
˜'Jšžœ
žœ"˜+J˜J˜—š’
œ˜Jšžœ
žœžœ˜(J˜J˜—š’œ˜ Jšœžœžœ
˜'Jšžœ žœžœ˜@J˜J˜—š
’œžœžœžœžœ˜%J™8J™<Jšœžœ˜Jšœžœ
˜Jšœžœ˜Jšœžœ˜Jšœ	žœ˜šžœ
ž˜J˜J˜J˜Jšžœ
˜—Jšžœ˜J˜——š
 ™šœžœžœ˜Jšœ	žœ€!˜5Jšœ!€#˜DJšœ €2˜RJšœžœ
žœ
€.˜LJ˜J˜—š’
œ˜J™BJšœžœžœ
˜'JšžœN˜TJ˜J˜—š
’
œžœžœžœ
žœ˜=šžœ
žœ
˜Jš
œU˜U—J˜——š
 
™
šœžœ=˜QJ˜—šœžœ˜²J˜—Jšœžœΐ™ΩJ˜Jšœ<žœžœžœA€˜¬
J˜Jšœ	₯4™=Jšœ=žœžœžœ˜PJšœ<žœ
žœžœ˜OJšœ>žœ
žœžœ˜QJšœ<žœžœžœ1˜~J˜Jšœ5žœžœ˜@J˜?J˜EJšœOžœ™VJ˜—Jš
ž˜J˜—…—!=š