DIRECTORY Basics, CtBasic, CtFilter, CtMod, ImagerPixel, ImagerSample, ImagerTransformation, IO, Process, Real, RealFns, Rope, SF;

CtFilterImpl: CEDAR PROGRAM
IMPORTS Basics, CtBasic, CtMod, ImagerPixel, ImagerSample, ImagerTransformation, IO, Process, Real, RealFns, Rope, SF
EXPORTS CtFilter

~ BEGIN
Pair:					TYPE ~ CtBasic.RealPair;
SampleMaps:			TYPE ~ CtBasic.SampleMaps;
PixelMap:				TYPE ~ ImagerPixel.PixelMap;
Sample:				TYPE ~ ImagerSample.Sample;
SampleBuffer:		TYPE ~ ImagerSample.SampleBuffer;
SampleMap:			TYPE ~ ImagerSample.SampleMap;
Transformation:		TYPE ~ ImagerTransformation.Transformation;
STREAM:		 		TYPE ~ IO.STREAM;
ROPE:					TYPE ~ Rope.ROPE;
Box:					TYPE ~ SF.Box;
Vec:					TYPE ~ SF.Vec;
Filter:					TYPE ~ CtFilter.Filter;
FilterEvalProc:		TYPE ~ CtFilter.FilterEvalProc;
FilterInitProc:		TYPE ~ CtFilter.FilterInitProc;
FilterPrintProc:		TYPE ~ CtFilter.FilterPrintProc;
ScanlineBuf:			TYPE ~ CtFilter.ScanlineBuf;
ScanlineBufRec:		TYPE ~ CtFilter.ScanlineBufRec;
Scanline32Buf:		TYPE ~ CtFilter.Scanline32Buf;
Scanline32BufRec:	TYPE ~ CtFilter.Scanline32BufRec;
Scanline:				TYPE ~ CtFilter.Scanline;
ScanlineRec:			TYPE ~ CtFilter.ScanlineRec;
Pixel:					TYPE ~ CtFilter.Pixel;
Scanline32:			TYPE ~ CtFilter.Scanline32;
Scanline32Rec:		TYPE ~ CtFilter.Scanline32Rec;
Pixel32:				TYPE ~ CtFilter.Pixel32;
Int16Seq:				TYPE ~ CtFilter.Int16Seq;
Int16SeqRec:			TYPE ~ CtFilter.Int16SeqRec;
FilterBuf:				TYPE ~ CtFilter.FilterBuf;
FilterBufRec:			TYPE ~ CtFilter.FilterBufRec;
WeightBuf:			TYPE ~ CtFilter.WeightBuf;
WeightBufRec:		TYPE ~ CtFilter.WeightBufRec;
PixelMapping:		TYPE ~ CtFilter.PixelMapping;
PixelMappingRec:	TYPE ~ CtFilter.PixelMappingRec;
BoolSeq:				TYPE ~ CtFilter.BoolSeq;
BoolSeqRec:			TYPE ~ CtFilter.BoolSeqRec;
BoolFalse:			TYPE ~ CtFilter.BoolFalse;
Error: PUBLIC ERROR [reason: ROPE] ~ CODE;
Expand: PUBLIC PROC [map: SampleMap, srcBox, dstBox: Box] ~ {
sBox: Box ¬ [SF.Sub[srcBox.min, [1, 1]], SF.Add[srcBox.max, [1, 1]]]; 
dBox: Box ¬ SF.Intersect[ImagerSample.GetBox[map], dstBox];
srcPM: PixelMap ¬ ImagerPixel.MakePixelMap[ImagerSample.Copy[map: map, box: sBox]];
scaleF: REAL ¬ REAL[SF.SizeF[dstBox]]/SF.SizeF[srcBox];
scaleS: REAL ¬ REAL[SF.SizeS[dstBox]]/SF.SizeS[srcBox];
CtMod.TransformMap[map, ImagerTransformation.Cat[
ImagerTransformation.Translate[[-srcBox.min.s, -srcBox.min.f]],
ImagerTransformation.Scale2[[scaleS, scaleF]],
ImagerTransformation.Translate[[dstBox.min.s, dstBox.min.f]]]];
};
CheckForAbort: PROC ~ {Process.CheckForAbort[]};

Resample: PUBLIC PROC [maps: SampleMaps, src, dst: Box, filter: Filter ¬ [], s: STREAM ¬ NIL] ~ {
ax0: NAT ¬ src.min.f;
ay0: NAT ¬ src.min.s;
bx0: NAT ¬ dst.min.f;
by0: NAT ¬ dst.min.s;
sx: REAL ¬ REAL[dst.max.f-bx0]/(src.max.f-ax0);		-- scale factors of resampling
sy: REAL ¬ REAL[dst.max.s-by0]/(src.max.s-ay0);
ResampleContinuous[maps, src, dst, [sx, sy], [bx0-.5-sx*(ax0-.5), by0-.5-sy*(ay0-.5)], filter, s];
};
ResampleContinuous: PUBLIC PROC [maps: SampleMaps, src, dst: Box, scale, tran: Pair, filter: Filter ¬ [], s: STREAM ¬ NIL] ~ {

filt: Filter ¬ IF filter = [] THEN FilterFind["triangle"] ELSE filter;
ax0: NAT ¬ src.min.f;	ax1: NAT ¬ src.max.f;
ay0: NAT ¬ src.min.s;	ay1: NAT ¬ src.max.s;
bx0: NAT ¬ dst.min.f;	bx1: NAT ¬ dst.max.f;
by0: NAT ¬ dst.min.s;	by1: NAT ¬ dst.max.s;
sx: REAL ~ scale.x;												-- scale of warp
sy: REAL ~ scale.y;
tx: REAL ~ tran.x;												-- translate of warp
ty: REAL ~ tran.y;
axrad: REAL ~ filt.blur*MAX[1., 1./sx];						-- scale of filter in a space
ayrad: REAL ~ filt.blur*MAX[1., 1./sy];
axsupp: REAL ~ MAX[.5, axrad*filt.support];					-- radius of scaled filter
aysupp: REAL ~ MAX[.5, ayrad*filt.support];
epsilon: REAL ~ 1e-4;								-- fudge factor to make sure buffers are big enough; note that Round[2.5]=2 and Round[3.5]=4, by IEEE floating point rules.  Without epsilon this rounding rule causes bounds faults on yweight in MakeWeightBuf
bx0c, by0c, bx1c, by1c: INT16;						-- valid dst box, we clip to this

[[[ay0, ax0], [ay1, ax1]]] ¬ SF.Intersect[[[ay0, ax0], [ay1, ax1]], maps.box];
[[[by0, bx0], [by1, bx1]]] ¬ SF.Intersect[[[by0, bx0], [by1, bx1]], maps.box];
bx0c ¬ Real.Ceiling[sx*(ax0-axsupp)+tx+epsilon];
by0c ¬ Real.Ceiling[sy*(ay0-aysupp)+ty+epsilon];
bx1c ¬ Real.Floor[sx*(ax1-1+axsupp)+tx-epsilon]+1;
by1c ¬ Real.Floor[sy*(ay1-1+aysupp)+ty-epsilon]+1;

IF bx0 < bx0c OR bx1 > bx1c OR by0 < by0c OR by1 > by1c THEN {
IF s # NIL THEN {
PrintBox[s, "clipping odst=", bx0, by0, bx1, by1];
PrintBox[s, " to ", bx0c, by0c, bx1c, by1c];
IO.PutRope[s, "\n"];
};
[[[by0, bx0], [by1, bx1]]] ¬
SF.Intersect[[[by0, bx0], [by1, bx1]],  [[by0c, bx0c], [by1c, bx1c]]];
};
{
anx: NAT ~ ax1-ax0;									-- size of source image
any: NAT ~ ay1-ay0;
bnx: NAT ~ bx1-bx0;									-- size of destination image
bny: NAT ~ by1-by0;
overlap: BOOL ¬ ax0<bx1 AND ax1>bx0 AND ay0<by1 AND ay1>by0;	-- do src&dst overlap?
ux: REAL ~ bx0-sx*(ax0-.5)-tx;		-- precompute offsets for Map routines, normally .5
uy: REAL ~ by0-sy*(ay0-.5)-ty;
rgb: BOOL ~ maps.bpp = 24;							-- black&white or RGB?

MapI: PROC [x: NAT, scale, offset: REAL] RETURNS [INT16] ~ INLINE
{RETURN[Trunc[(REAL[x]+offset)/scale]];};
MapR: PROC [x, scale, offset: REAL] RETURNS [REAL] ~ INLINE
{RETURN[(x+offset)/scale];};

MakeMapTable: PROC [scale, tran, asupp: REAL, a0, b0, bn: NAT, overlap: BOOL, map: Int16Seq] ~ {

b, bi: NAT ¬ 0;
split: INT16;

IF overlap											-- src and dst boxes overlap
THEN {
IF scale = 1.
THEN											-- no scale change, translation only
split ¬ IF a0 < b0 THEN 0 ELSE bn
ELSE {										-- magnify or minify
IF scale > 1.								-- magnify
THEN split ¬ Real.Floor[(tran+scale*asupp-.5)/(1.-scale)-b0+1.]
ELSE {								-- minify
i0: INT16;
u, cen: REAL;
split ¬ Real.Ceiling[(tran+scale*asupp)/(1.-scale)-b0];
u ¬ REAL[b0]-scale*(REAL[a0]-.5)-tran;		-- recalculate offset
cen ¬ MapR[split-1, scale, u];					-- cen at b=split-1
i0 ¬ Round[cen-asupp];						-- i0 at b=split-1
IF a0+i0 >= b0+split THEN {split ¬ split-1; IO.PutRope[s, "! "];};
};
IF split < 0 THEN split ¬ 0
ELSE IF split > bn THEN split ¬ bn;
IF s # NIL THEN IO.PutF1[s, "split at y=%g\n", IO.int[b0+split]];
};

IF scale >= 1.				-- magnify: scan in toward the split
THEN {
FOR b				   IN [0..split)	DO map[bi] ¬ b;  bi ¬ bi+1; ENDLOOP;
FOR b DECREASING IN [split..bn)	DO map[bi] ¬ b;  bi ¬ bi+1; ENDLOOP;
}
ELSE {				-- minify: scan out away from the split
FOR b DECREASING IN [0..split)	DO map[bi] ¬ b;  bi ¬ bi+1; ENDLOOP;
FOR b				    IN [split..bn)  DO map[bi] ¬ b;  bi ¬ bi+1; ENDLOOP;
};
}
ELSE							-- no overlap, either order will work, we opt for forward
FOR b IN [0..bn)	DO map[bi] ¬ b;  bi ¬ bi+1; ENDLOOP;

}; -- end of MakeMapTable

ResampleUnfiltered: PROC ~ TRUSTED {
ay, ayold, bx, by, byi: NAT;
xmap1, xmap2, xmap3: PixelMapping;
ymap: Int16Seq ~ NEW[Int16SeqRec[bny]];
abuf1, abuf2, abuf3, bbuf1, bbuf2, bbuf3: SampleBuffer;
xmap1 ¬ NEW[PixelMappingRec[bnx]];
abuf1 ¬ ImagerSample.ObtainScratchSamples[anx];	-- src scan line buffer
bbuf1 ¬ ImagerSample.ObtainScratchSamples[bnx];	-- dst scan line buffer
IF rgb THEN {
xmap2 ¬ NEW[PixelMappingRec[bnx]];
xmap3 ¬ NEW[PixelMappingRec[bnx]];
abuf2 ¬ ImagerSample.ObtainScratchSamples[anx];		-- src scan line buffer
abuf3 ¬ ImagerSample.ObtainScratchSamples[anx];		-- src scan line buffer
bbuf2 ¬ ImagerSample.ObtainScratchSamples[bnx];		-- dst scan line buffer
bbuf3 ¬ ImagerSample.ObtainScratchSamples[bnx];		-- dst scan line buffer
};

IF s # NIL THEN {
IO.PutRope[s, "filt=point\n"];
IO.PutF[s, "X  scale=%g tran=%g\n", IO.real[sx], IO.real[tx]];
IO.PutF[s, "Y  scale=%g tran=%g\n", IO.real[sy], IO.real[ty]];
};

MakeMapTable[sy, ty, 0.5, ay0, by0, bny, overlap, ymap];-- y map will be in place
FOR bx IN [0..bnx) DO								-- x mapping is not (but it could be)
ax: NAT ~ MapI[bx, sx, ux];
xmap1[bx] ¬ @abuf1[ax]; -- address of source pixel (in abuf1) that maps to bbuf1[bx]
IF rgb THEN {
xmap2[bx] ¬ @abuf2[ax]; -- source pixel address (in abuf2) that maps to bbuf2[bx]
xmap3[bx] ¬ @abuf3[ax]; -- source pixel address (in abuf2) that maps to bbuf2[bx]
};
ENDLOOP;

ayold ¬ any;	-- impossible value for ay
FOR byi IN [0..bny) DO								-- loop over dest scanlines
CheckForAbort[];
by ¬ ymap[byi];
ay ¬ MapI[by, sy, uy];
IF ay # ayold THEN {
bbuf1p, bbuf2p, bbuf3p: LONG POINTER TO Sample;
xmap1p, xmap2p, xmap3p: LONG POINTER TO LONG POINTER TO Sample;
bbuf1p ¬ @bbuf1[0];
xmap1p ¬ @xmap1[0];
IF rgb THEN {
bbuf2p ¬ @bbuf2[0];
bbuf3p ¬ @bbuf3[0];
xmap2p ¬ @xmap2[0];
xmap3p ¬ @xmap3[0];
};
ayold ¬ ay;
ImagerSample.GetSamples[maps[0].map, [ay0+ay, ax0],, abuf1,, anx];
IF rgb THEN {
ImagerSample.GetSamples[maps[1].map, [ay0+ay, ax0],, abuf2,, anx]; 
ImagerSample.GetSamples[maps[2].map, [ay0+ay, ax0],, abuf3,, anx]; 
};
IF rgb
THEN
FOR bx IN [0..bnx) DO							-- resample abuf into bbuf
bbuf1p­ ¬ xmap1p­­;						-- copy pixel abuf[ax] to bbuf[bx]
bbuf2p­ ¬ xmap2p­­;						-- copy pixel abuf[ax] to bbuf[bx]
bbuf3p­ ¬ xmap3p­­;						-- copy pixel abuf[ax] to bbuf[bx]
bbuf1p ¬ bbuf1p+SIZE[Sample];
bbuf2p ¬ bbuf2p+SIZE[Sample];
bbuf3p ¬ bbuf3p+SIZE[Sample];
xmap1p ¬ xmap1p+SIZE[LONG POINTER TO Sample];
xmap2p ¬ xmap2p+SIZE[LONG POINTER TO Sample];
xmap3p ¬ xmap3p+SIZE[LONG POINTER TO Sample];
ENDLOOP
ELSE
FOR bx IN [0..bnx) DO							-- resample abuf into bbuf
bbuf1p­ ¬ xmap1p­­;						-- copy pixel abuf[ax] to bbuf[bx]
bbuf1p ¬ bbuf1p+SIZE[Sample];
xmap1p ¬ xmap1p+SIZE[LONG POINTER TO Sample];
ENDLOOP;
};
ImagerSample.PutSamples[maps[0].map, [by0+by, bx0],, bbuf1,, bnx]; -- bbuf1 scan
IF rgb THEN {
ImagerSample.PutSamples[maps[1].map, [by0+by, bx0], , bbuf2, , bnx]; -- bbuf2
ImagerSample.PutSamples[maps[2].map, [by0+by, bx0], , bbuf3, , bnx]; -- bbuf3
};
ENDLOOP;

FOR l: LIST OF SampleBuffer ¬ IF rgb
THEN LIST[abuf1, abuf2, abuf3, bbuf1, bbuf2, bbuf3]
ELSE LIST [abuf1, bbuf1],
l.rest WHILE l # NIL DO
ImagerSample.ReleaseScratchSamples[l.first];
ENDLOOP;
}; -- end of ResampleUnfiltered

ResampleFiltered: PROC [] ~ {

MakeWeightBuf: PROC [b, a0: NAT, cen, rad, supp: REAL, len: NAT, wbuf: WeightBuf ¬ NIL] RETURNS [WeightBuf] ~ {
i0: INT16 ¬ Round[cen-supp];
i1: INT16 ¬ Round[cen+supp];
den: INT32 ¬ 0;
sum: INT16 ¬ 0;
sc: REAL;
IF i0 < 0 THEN i0 ¬ 0;
IF i1 > len THEN i1 ¬ len;
IF i0 >= i1 THEN Error["null horizontal filter"];
IF wbuf = NIL THEN wbuf ¬ NEW[WeightBufRec[i1-i0]];
wbuf.i0 ¬ i0;
wbuf.i1 ¬ i1;
FOR i: NAT IN [i0..i1) DO
den ¬ den+Round[weightone*filt.evalProc[(REAL[i]+.5-cen)/rad, filt.clientData]];
ENDLOOP;
sc ¬ IF den = 0. THEN weightone ELSE REAL[weightone]*weightone/den;
sum ¬ 0;
IF ABS[sc] > LAST[INT16]
THEN FOR i: NAT IN [i0..i1) DO	-- extra care at filter fringes to avoid overflow
t: INT16 ~ Round[MAX[FIRST[INT16], MIN[LAST[INT16], sc*filt.evalProc[(REAL[i]+.5-cen)/rad, filt.clientData]]]];
wbuf[i-i0] ¬ t;
sum ¬ sum+t;
ENDLOOP
ELSE FOR i: NAT IN [i0..i1) DO			-- normal case
t: INT16 ~ Round[sc*filt.evalProc[(REAL[i]+.5-cen)/rad, filt.clientData]];
wbuf[i-i0] ¬ t;
sum ¬ sum+t;
ENDLOOP;
IF sum # weightone THEN {
i: NAT ~ MAX[i0, MIN[i1-1, Round[cen]]];
wbuf[i-i0] ¬ wbuf[i-i0]+weightone-sum;
};
RETURN[wbuf];
}; -- end of MakeWeightBuf

abuf: Scanline ~ NEW[ScanlineRec[anx]];	-- src scanline buffer
bbuf: Scanline ~ NEW[ScanlineRec[bnx]];	-- dst scanline buffer
scratch: SampleBuffer ~ ImagerSample.ObtainScratchSamples[MAX[anx, bnx]];	-- scratch scanline buffer
accum: Scanline32 ~ NEW[Scanline32Rec[bnx]];			-- accumulator buffer
axwid: NAT ~ Real.Ceiling[2.*axsupp+epsilon];			-- diameter of this filter
aywid: NAT ~ Real.Ceiling[2.*aysupp+epsilon];
ymap: Int16Seq ~ NEW[Int16SeqRec[bny]];
linebuf: Scanline32Buf ~ NEW[Scanline32BufRec[aywid]];	-- circular buffer of scanlines
linewritten: BoolSeq ~ NEW[BoolSeqRec[MAX[ay1, by1]]];-- line y been written? (debug)
xfilt: FilterBuf ~ NEW[FilterBufRec[bnx]];		 -- one WeightBuf for each dest x pos
yweight: WeightBuf ~ NEW[WeightBufRec[aywid]]; -- WeightBuf for one dest y pos
weightshift: INT16 ~ 14;
weightone: INT16 ~ Basics.BITSHIFT[1, weightshift];
finalshift: NAT ~ 2*weightshift-8;					-- shift after x and y filter passes

abuf.len ¬ anx;
bbuf.len ¬ bnx;
accum.len ¬ bnx;
FOR byf: NAT IN [0..aywid) DO
linebuf[byf] ¬ NEW[Scanline32Rec[bnx]];
linebuf[byf].len ¬ bnx;
linebuf[byf].y ¬ -1;						-- mark scanline as unread
ENDLOOP;

IF s # NIL THEN {
IO.PutF[s, "filt=%g supp=%g\n", IO.rope[filt.name], IO.real[filt.support]];
IO.PutFL[s, "X  scale=%g tran=%g,  rad=%g supp=%g wid=%g\n",
LIST[IO.real[sx], IO.real[tx], IO.real[axrad], IO.real[axsupp], IO.int[axwid]]];
IO.PutFL[s, "Y  scale=%g tran=%g,  rad=%g supp=%g wid=%g\n",
LIST[IO.real[sy], IO.real[ty], IO.real[ayrad], IO.real[aysupp], IO.int[aywid]]];
};

FOR bx: NAT IN [0..bnx) DO
xfilt[bx] ¬ MakeWeightBuf[bx0+bx, ax0, MapR[bx, sx, ux], axrad, axsupp, anx, NIL];
ENDLOOP;

MakeMapTable[sy, ty, aysupp, ay0, by0, bny, overlap, ymap];

FOR byi: NAT IN [0..bny) DO					-- loop over dest scanlines
by: NAT ~ ymap[byi];
CheckForAbort[];

[] ¬ MakeWeightBuf[by0+by, ay0, MapR[by, sy, uy], ayrad, aysupp, any, yweight];
IF rgb
THEN ScanlineZeroRGB[accum]
ELSE ScanlineZeroBW[accum];

FOR ayf: NAT IN [yweight.i0 .. yweight.i1) DO	-- loop over src scanlines that influence this dest scanline
tbuf: Scanline32 ~ linebuf[ayf MOD aywid];
IF tbuf.y # ayf THEN {						-- scanline needs to be read in
IF overlap AND linewritten[ay0+ayf] THEN
Error[IO.PutFR1["FEEDBACK ON LINE %g\n", IO.int[ay0+ayf]]];
tbuf.y ¬ ayf;
ScanlineRead[maps, ax0, ay0+ayf, scratch, abuf];	-- read scan line into abuf
IF rgb
THEN ScanlineResampleRGB[filt, xfilt, abuf, tbuf]
ELSE ScanlineResampleBW[filt, xfilt, abuf, tbuf];
};

IF rgb
THEN ScanlineAccumRGB[filt, yweight[ayf-yweight.i0], tbuf, accum]
ELSE ScanlineAccumBW[filt, yweight[ayf-yweight.i0], tbuf, accum];

ENDLOOP;

IF rgb
THEN ScanlineShiftRGB[filt, finalshift, accum, bbuf]
ELSE ScanlineShiftBW[filt, finalshift, accum, bbuf];

ScanlineWrite[maps, bx0, by0+by, scratch, bbuf];		 -- write scan line from bbuf
linewritten[by0+by] ¬ TRUE;
ENDLOOP;

ImagerSample.ReleaseScratchSamples[scratch];
}; -- end of ResampleFiltered

IF s # NIL THEN {
PrintBox[s, "src=", ax0, ay0, ax1, ay1];
PrintBox[s, " dst=", bx0, by0, bx1, by1];
IO.PutRope[s, "\n"];
};
IF anx <= 0 OR any <= 0 OR bnx <= 0 OR bny <= 0 THEN RETURN;
IF axsupp <= .5 AND aysupp <= .5 THEN filt ¬ FilterFind["point"];

IF anx = bnx AND any = bny AND filt.cardinal AND filt.blur = 1. AND ABS[tx-Round[tx]] < epsilon AND ABS[ty-Round[ty]] < epsilon
THEN												-- optimize unfiltered translation
CtBasic.MoveMaps[maps, src.min, [by0, bx0], [any, anx]]
ELSE												-- scaling and/or filtering required
IF Rope.Equal[filt.name, "point", FALSE]
THEN ResampleUnfiltered[]
ELSE ResampleFiltered[];
};
}; -- end of ResampleContinuous

ResampleMap: PUBLIC PROC [map: SampleMap, src, dst: Box, filter: Filter ¬ [], s: STREAM ¬ NIL] ~ {
box: Box ¬ ImagerSample.GetBox[map];
size: Vec ¬ ImagerSample.GetSize[map];
maps: SampleMaps ¬ CtBasic.CreateMaps[8, box.min.f, box.min.s, size.f, size.s, FALSE];
maps[0].map ¬ map;
Resample[maps, src, dst, filter, s];
};
FilterSeq: TYPE ~ CtFilter.FilterSeq;
FilterSeqRep: TYPE ~ CtFilter.FilterSeqRep;

filterRegistry: FilterSeq ¬ NIL;							-- the registry

FiltersGet: PUBLIC PROC RETURNS [FilterSeq] ~ {RETURN[filterRegistry]};

FilterFind: PUBLIC PROC [name: ROPE] RETURNS [Filter] ~ {
i: INT16 ¬ FilterFindIndex[name];
IF i < 0
THEN Error[IO.PutFR1["no such filter: %g", IO.rope[name]]]
ELSE RETURN[filterRegistry[i]];
};

FilterCatalog: PUBLIC PROC [s: STREAM] ~ {
IO.PutRope[s, "FILTER    SUPPORT\n"];
FOR i: NAT IN [0..filterRegistry.length) DO
filt: Filter ~ filterRegistry[i];
IO.PutF[s, " %-9g    %4.2f%g\n", IO.rope[filt.name], IO.real[filt.support], IO.rope[IF filt.windowMe THEN " (windowed by default)" ELSE ""]];
IF filt.printProc # NIL THEN filt.printProc[s, filt.clientData];
ENDLOOP;
};
FilterRegister: PUBLIC PROC [
name: ROPE,						-- name of filter
evalProc: FilterEvalProc,		-- the filter function for (usually) unit-spaced samples
support: REAL,					-- radius of nonzero portion, if FIR, or nominal width, if IIR
blur: REAL ¬ 1.,					-- blur factor, 1=normal, <1=sharp, >1=blurry
windowMe: BOOL,				-- should filter be windowed? (is it IIR?)
cardinal: BOOL,					-- is this filter cardinal,
							   	   -- ie, does eval(x) = (if x=0 THEN 1 ELSE 0) for integer x?
unitrange: BOOL ¬ FALSE,		-- does filter stay within the range [0..1]?
initProc: FilterInitProc ¬ NIL,	-- call this once to initialize client data, if any
printProc: FilterPrintProc ¬ NIL,-- print client data, if any
scale: REAL ¬ 1.0,				-- rescale filter result before quantization to pixel value
clientData: REF ANY ¬ NIL		-- client info to be passed to evalProc
] ~ {
IF filterRegistry = NIL THEN filterRegistry ¬ NEW[FilterSeqRep[10]];
{
i: INT16 ¬ FilterFindIndex[name];			-- if already exists, modify
IF i < 0 THEN i ¬ filterRegistry.length;	-- if new, append to sequence
IF i >= filterRegistry.lengthmax
THEN filterRegistry ¬ LengthenFilterSeq[filterRegistry, 1.3];
filterRegistry[i] ¬ [name, evalProc, support, blur, windowMe, cardinal, unitrange, initProc, printProc, scale, clientData];
IF i >= filterRegistry.length THEN filterRegistry.length ¬ i+1;
};
};

LengthenFilterSeq: PROC [filts: FilterSeq, factor: REAL ¬ 1.3] RETURNS [new: FilterSeq] ~ {
newlen: NAT ¬ MAX[Round[factor*filts.lengthmax], INT[filts.lengthmax+1]];
new ¬ NEW[FilterSeqRep[newlen]];
FOR i: NAT IN [0..filts.length) DO new[i] ¬ filts[i]; ENDLOOP;
new.length ¬ filts.length;
};
FilterUnregisterAll: PUBLIC PROC ~ {
IF filterRegistry # NIL THEN filterRegistry.length ¬ 0;
};
FilterUnregister: PUBLIC PROC [name: ROPE] ~ {
i: INT16 ¬ FilterFindIndex[name];
IF i >= 0 THEN {
FOR j: NAT IN [i+1..filterRegistry.length) DO
filterRegistry[j-1] ¬ filterRegistry[j];
ENDLOOP;
filterRegistry.length ¬ filterRegistry.length-1;
};
};

FilterRegisterDefaults: PUBLIC PROC ~ {
FilterRegister["point",		NIL,				0.0, 1.0,	FALSE,	TRUE,		TRUE];
FilterRegister["box",			BoxEval,			0.5, 1.0,	FALSE,	TRUE,		TRUE];
FilterRegister["triangle",	TriangleEval,	1.0, 1.0,	FALSE,	TRUE,		TRUE];
FilterRegister["quadratic",	QuadraticEval,	1.5, 1.0,	FALSE,	FALSE,	TRUE];
FilterRegister["cubic",		CubicEval,		2.0, 1.0,	FALSE,	FALSE,	TRUE];

FilterRegister["catrom",		CatRomEval,		2.0, 1.0,	FALSE,	TRUE,		FALSE];
FilterRegister["gaussian",	GaussianEval,	1.25,1.0,	FALSE,	FALSE,	TRUE];
FilterRegister["sinc",		SincEval,			4.0, 1.0,	TRUE,		TRUE,		FALSE];
FilterRegister["bessel",		BesselEval,		3.2383,1.0,	TRUE,		FALSE,	FALSE];

FilterRegister["normal",		NormalEval,		2.5, 0.5,	FALSE,	FALSE,	FALSE];

};

FilterFindIndex: PROC [name: ROPE] RETURNS [INT16] ~ {
IF filterRegistry = NIL THEN RETURN[-1];
FOR i: NAT IN [0..filterRegistry.length) DO
IF Rope.Equal[name, filterRegistry[i].name, FALSE] THEN RETURN[i];
ENDLOOP;
RETURN[-1];
};
pi: REAL ~ 3.1415927;

BoxEval: PUBLIC FilterEvalProc ~ {
SELECT x FROM
<-.5 => RETURN[0.];
< .5 => RETURN[1.];
ENDCASE => RETURN[0.];
};

TriangleEval: PUBLIC FilterEvalProc ~ {
SELECT x FROM
<-1. => RETURN[0.];
< 0. => RETURN[1.+x];
< 1. => RETURN[1.-x];
ENDCASE => RETURN[0.];
};

QuadraticEval: PUBLIC FilterEvalProc ~ {
SELECT x FROM
<-1.5 => RETURN[0.];
<- .5 => {t: REAL ~ x+1.5; RETURN[.5*t*t];};
<  .5 => RETURN[.75-x*x];
< 1.5 => {t: REAL ~ x-1.5; RETURN[.5*t*t];};
ENDCASE => RETURN[0.];
};

CubicEval: PUBLIC FilterEvalProc ~ {
SELECT x FROM
<-2. => RETURN[0.];
<-1. => {t: REAL ~ 2.+x; RETURN[t*t*t/6.];};
< 0. => RETURN[(4.+x*x*(-6.+x*-3.))/6.];
< 1. => RETURN[(4.+x*x*(-6.+x*3.))/6.];
< 2. => {t: REAL ~ 2.-x; RETURN[t*t*t/6.];};
ENDCASE => RETURN[0.];
};

CatRomEval: PUBLIC FilterEvalProc ~ {
SELECT x FROM
<-2. => RETURN[0.];
<-1. => RETURN[.5*(4.+x*(8.+x*(5.+x)))];
< 0. => RETURN[.5*(2.+x*x*(-5.+x*-3.))];
< 1. => RETURN[.5*(2.+x*x*(-5.+x*3.))];
< 2. => RETURN[.5*(4.+x*(-8.+x*(5.-x)))];
ENDCASE => RETURN[0.];
};

GaussianEval: PUBLIC FilterEvalProc ~ {
RETURN[RealFns.Exp[-2.*x*x]*RealFns.SqRt[2./pi]];
};

SincEval: PUBLIC FilterEvalProc ~ {
RETURN[IF x = 0. THEN 1. ELSE RealFns.Sin[pi*x]/(pi*x)];
};

BesselEval: PUBLIC FilterEvalProc ~ {
RETURN[IF x = 0. THEN pi/4. ELSE RealFns.J1[pi*x]/(2.*x)];
};


NormalEval: PUBLIC FilterEvalProc ~ {
RETURN[RealFns.Exp[-.5*x*x]/RealFns.SqRt[2.*pi]];
};
ScanlineRead: PROC [maps: SampleMaps, x0, y0: NAT, scratch: SampleBuffer, buf: Scanline] ~ {
npixel: NAT ~ buf.len;
FOR i: NAT IN [0..maps.nChannels) DO
ImagerSample.GetSamples[maps[i].map, [y0, x0], , scratch, , npixel]; -- read scan line
SELECT maps[i].type FROM
$Red =>
FOR x: NAT IN [0..npixel) DO
buf[x].r ¬ scratch[x];
ENDLOOP;
$Green =>
FOR x: NAT IN [0..npixel) DO
buf[x].g ¬ scratch[x];
ENDLOOP;
$Blue =>
FOR x: NAT IN [0..npixel) DO
buf[x].b ¬ scratch[x];
ENDLOOP;
$BW =>
FOR x: NAT IN [0..npixel) DO
buf[x].r ¬ scratch[x];
ENDLOOP;
ENDCASE => ERROR;
ENDLOOP;
};

ScanlineWrite: PROC [maps: SampleMaps, x0, y0: NAT, scratch: SampleBuffer, buf: Scanline] ~ {
npixel: NAT ~ buf.len;
FOR i: NAT IN [0..maps.nChannels) DO
SELECT maps[i].type FROM
$Red		=> FOR x: NAT IN [0..npixel) DO scratch[x] ¬ buf[x].r; ENDLOOP;
$Green	=> FOR x: NAT IN [0..npixel) DO scratch[x] ¬ buf[x].g; ENDLOOP;
$Blue		=> FOR x: NAT IN [0..npixel) DO scratch[x] ¬ buf[x].b; ENDLOOP;
$BW		=> FOR x: NAT IN [0..npixel) DO scratch[x] ¬ buf[x].r; ENDLOOP;
ENDCASE => ERROR;
ImagerSample.PutSamples[maps[i].map, [y0, x0],, scratch,, npixel]; -- write scan line
ENDLOOP;
};

ScanlineResampleBW: PROC [filt: Filter, xfilt: FilterBuf, abuf: Scanline, tbuf: Scanline32] ~ {
FOR bx: NAT IN [0..tbuf.len) DO
TRUSTED {
i0: NAT ~ xfilt[bx].i0;
i1: NAT ~ xfilt[bx].i1;
xfiltp: LONG POINTER TO INT16 ¬ @xfilt[bx][0];
abufp: LONG POINTER TO Pixel ¬ @abuf[i0];
sum: INT32 ¬ 0;
IF filt.unitrange
THEN									-- if filter in [0..1] then use fast CARD16 multiply
FOR axf: NAT IN [i0..i1) DO
sum ¬ sum+LongMult[xfiltp­, abufp.r];
xfiltp ¬ xfiltp+SIZE[INT16];
abufp ¬ abufp+SIZE[Pixel];
ENDLOOP
ELSE									-- else use slower, more general INT32 multiply
FOR axf: NAT IN [i0..i1) DO
sum ¬ sum+INT32[xfiltp­]*abufp.r;
xfiltp ¬ xfiltp+SIZE[INT16];
abufp ¬ abufp+SIZE[Pixel];
ENDLOOP;
tbuf[bx].r ¬ IF filt.scale # 1 THEN Real.Round[REAL[sum]*(filt.scale/256.0)] ELSE sum/256;
};
ENDLOOP;
};

ScanlineResampleRGB: PROC [filt: Filter, xfilt: FilterBuf, abuf: Scanline, tbuf: Scanline32] ~ {
FOR bx: NAT IN [0..tbuf.len) DO
TRUSTED {
i0: NAT ~ xfilt[bx].i0;
i1: NAT ~ xfilt[bx].i1;
xfiltp: LONG POINTER TO INT16 ¬ @xfilt[bx][0];
abufp: LONG POINTER TO Pixel ¬ @abuf[i0];
sumr: INT32 ¬ 0;
sumg: INT32 ¬ 0;
sumb: INT32 ¬ 0;
IF filt.unitrange
THEN									-- if filter in [0..1] then use fast CARD16 multiply
FOR axf: NAT IN [i0..i1) DO
sumr ¬ sumr+LongMult[xfiltp­, abufp.r];
sumg ¬ sumg+LongMult[xfiltp­, abufp.g];
sumb ¬ sumb+LongMult[xfiltp­, abufp.b];
xfiltp ¬ xfiltp+SIZE[INT16];
abufp ¬ abufp+SIZE[Pixel];
ENDLOOP
ELSE									-- else use slower, more general INT32 multiply
FOR axf: NAT IN [i0..i1) DO
weight: INT32 ¬ xfiltp­;
sumr ¬ sumr+weight*abufp.r;
sumg ¬ sumg+weight*abufp.g;
sumb ¬ sumb+weight*abufp.b;
xfiltp ¬ xfiltp+SIZE[INT16];
abufp ¬ abufp+SIZE[Pixel];
ENDLOOP;
tbuf[bx].r ¬ sumr/256;
tbuf[bx].g ¬ sumg/256;
tbuf[bx].b ¬ sumb/256;
};
ENDLOOP;
};

ScanlineZeroBW: PROC [accum: Scanline32] ~ {
FOR bx: NAT IN [0..accum.len) DO
accum[bx].r ¬ 0;
ENDLOOP
};

ScanlineZeroRGB: PROC [accum: Scanline32] ~ {
FOR bx: NAT IN [0..accum.len) DO
accum[bx].r ¬ 0;
accum[bx].g ¬ 0;
accum[bx].b ¬ 0;
ENDLOOP;
};

ScanlineAccumBW: PROC [filt: Filter, weight: INT16, tbuf: Scanline32, accum: Scanline32]
~ TRUSTED {
accump: LONG POINTER TO Pixel32 ¬ @accum[0];
tbufp: LONG POINTER TO Pixel32 ¬ @tbuf[0];
bnx: NAT ~ tbuf.len;
IF filt.unitrange
THEN									-- if filter in [0..1] then use fast CARD16 multiply
FOR bx: NAT IN [0..bnx) DO
accump.r ¬ accump.r+LongMult32[weight, tbufp.r];
accump ¬ accump+SIZE[Pixel32];
tbufp ¬ tbufp+SIZE[Pixel];
ENDLOOP
ELSE									-- else use slower, more general INT32 multiply
FOR bx: NAT IN [0..bnx) DO
accump.r ¬ accump.r+INT32[weight]*tbufp.r;
accump ¬ accump+SIZE[Pixel32];
tbufp ¬ tbufp+SIZE[Pixel];
ENDLOOP;
};

ScanlineAccumRGB: PROC [filt: Filter, weight: INT16, tbuf: Scanline32, accum: Scanline32]
~ TRUSTED {
accump: LONG POINTER TO Pixel32 ¬ @accum[0];
tbufp: LONG POINTER TO Pixel32 ¬ @tbuf[0];
bnx: NAT ~ tbuf.len;
IF filt.unitrange
THEN									-- if filter in [0..1] then use fast CARD16 multiply
FOR bx: NAT IN [0..bnx) DO
accump.r ¬ accump.r+LongMult[weight, tbufp.r];
accump.g ¬ accump.g+LongMult[weight, tbufp.g];
accump.b ¬ accump.b+LongMult[weight, tbufp.b];
accump ¬ accump+SIZE[Pixel32];
tbufp ¬ tbufp+SIZE[Pixel];
ENDLOOP
ELSE {								-- else use slower, more general INT32 multiply
w: INT32 ¬ weight;
FOR bx: NAT IN [0..bnx) DO
accump.r ¬ accump.r+w*tbufp.r;
accump.g ¬ accump.g+w*tbufp.g;
accump.b ¬ accump.b+w*tbufp.b;
accump ¬ accump+SIZE[Pixel32];
tbufp ¬ tbufp+SIZE[Pixel];
ENDLOOP;
};
};

ScanlineShiftBW: PROC [filt: Filter, shift: NAT, accum: Scanline32, bbuf: Scanline] ~ {
bnx: NAT ~ bbuf.len;
white: INT32 ~ LOOPHOLE[Basics.BITLSHIFT[LOOPHOLE[INT32[1]], 8+shift]];
IF filt.unitrange
THEN
FOR bx: NAT IN [0..bnx) DO
bbuf[bx].r ¬ LOOPHOLE[Basics.BITRSHIFT[LOOPHOLE[accum[bx].r], shift], INT32];
ENDLOOP
ELSE
FOR bx: NAT IN [0..bnx) DO
t: INT32 ¬ accum[bx].r;
bbuf[bx].r ¬ IF t < 0 THEN 0 ELSE IF t >= white THEN 255
ELSE INT16[LOOPHOLE[Basics.BITRSHIFT[LOOPHOLE[t], shift], INT32]];
ENDLOOP;
};

ScanlineShiftRGB: PROC [filt: Filter, shift: NAT, accum: Scanline32, bbuf: Scanline] ~ {
bnx: NAT ~ bbuf.len;
white: INT32 ~ LOOPHOLE[Basics.BITLSHIFT[LOOPHOLE[INT32[1]], 8+shift]];
IF filt.unitrange
THEN
FOR bx: NAT IN [0..bnx) DO
bbuf[bx].r ¬ LOOPHOLE[Basics.BITRSHIFT[LOOPHOLE[accum[bx].r], shift], INT32];
bbuf[bx].g ¬ LOOPHOLE[Basics.BITRSHIFT[LOOPHOLE[accum[bx].g], shift], INT32];
bbuf[bx].b ¬ LOOPHOLE[Basics.BITRSHIFT[LOOPHOLE[accum[bx].b], shift], INT32];
ENDLOOP
ELSE
FOR bx: NAT IN [0..bnx) DO
t: INT32 ¬ accum[bx].r;
bbuf[bx].r ¬ IF t < 0 THEN 0 ELSE IF t >= white THEN 255
ELSE INT16[LOOPHOLE[Basics.BITRSHIFT[LOOPHOLE[t], shift], INT32]];
t ¬ accum[bx].g;
bbuf[bx].g ¬ IF t < 0 THEN 0 ELSE IF t >= white THEN 255
ELSE INT16[LOOPHOLE[Basics.BITRSHIFT[LOOPHOLE[t], shift], INT32]];
t ¬ accum[bx].b;
bbuf[bx].b ¬ IF t < 0 THEN 0 ELSE IF t >= white THEN 255
ELSE INT16[LOOPHOLE[Basics.BITRSHIFT[LOOPHOLE[t], shift], INT32]];
ENDLOOP;
};
LongMult: PROC [a, b: CARD16] RETURNS [INT32] = INLINE {RETURN[a*b]};
LongMult32: PROC [a, b: CARD32] RETURNS [INT32] = INLINE {RETURN[a*b]};
Trunc: PROC [x: REAL] RETURNS [INT16] ~ INLINE
{RETURN[Real.Fix[x]];};

Round: PROC [x: REAL] RETURNS [INT16] ~ INLINE {RETURN[Real.Round[x]]};

PrintBox: PROC [s: STREAM, str: ROPE, x0, y0, x1, y1: INT16] ~ {
IO.PutFL[s, "%g[%g,%g)[%g,%g)",
LIST[IO.rope[str], IO.int[x0], IO.int[x1], IO.int[y0], IO.int[y1]]];
IO.PutF[s, "%gx%g", IO.int[x1-x0], IO.int[y1-y0]];
};
FilterRegisterDefaults[];			-- register the standard set of filters

END.

¶
CtFilterImpl.mesa
Copyright Ó 1988, 1992 by Xerox Corporation.  All rights reserved.
written by Paul Heckbert, August 18, 1988 5:26:54 pm PDT
Bloomenthal, November 22, 1992 2:23 pm PST (converted to 24 bpp: 3 maps, not RG, B maps)
Andrew Glassner November 14, 1990 4:12 pm PST
Imported Types
Local `Public' Types
Local `Private' Types
Errors
Simple Resampling
Fancy Filtered Resampling
Resample maps in place, mapping src to dst

Resample maps in place according to scale and tran, reading from src and writing to dst.
Notation: prefix 'a' denotes source coordinates, prefix 'b' denotes destination coordinates.

Boxes are inclusive at low end, exclusive at high end,
i.e. src box is [ax0..ax1)x[ay0..ay1) and dst is [bx0..bx1)x[by0..by1)
If axsupp and aysupp are both <= .5 then we've got point sampling.
Point sampling is essentially a special filter whose width is fixed at one src pixel:
Clip the src and dst windows to maps:
Some of the requested dst window lies outside the warped src window (plus support margin), so clip the dst window.
MapI and MapR convert discrete dst pixel coordinates to discrete src pixel coordinates
construct a table for the mapping a=(b-t)/s for b-b0 in [0..bn)
ordered so that the buffer resampling operation
FOR bi: IN [0..bn) DO
b _ map[bi];
a _ MapI(b, scale, offset);
buf[b] _ buf[a];		-- for point sampling
(for filtering, use weighted average of buf[a-asupp]..buf[a+asupp])
ENDLOOP
can work in place without feedback

find fixed point of mapping: let split=b-b0 at the point where a=b
	-- if moving left then scan left to right, else scan right to left
THE MAGIC SPLIT FORMULA:
Choose integer nearest midpoint of valid interval:
	x < split <= x+scale/(scale-1)
where x = (tran+scale*asupp)/(1-scale)-b0
Only one valid split point in this case:
	x <= split < x+1
so we take extra care (x as above)
The above formula is perfect for exact arithmetic, but hardware roundoff can cause split to be one too big.  Check for roundoff by simulating precisely the calculation of i0 in MakeWeightBuf.
if feedback would occur at b=split-1 then fix split point
Do unfiltered resampling (point sampling).
If rgb there are three image buffers (buf1=R, buf2=G, buf3=B), else one (buf1=BW)

scan line ay0+ay goes to by0+by
Read scan line(s):
Find scale factor sc to normalize the filter:
So that sum of sc*eval() is approx weightone:
IF sum = 0 AND s # NIL THEN IO.PutF[s, "x=%g sum=0\n", IO.int[b]];
This happens occasionally near fringes of sinc support
Try weightshift 10 if INT16 and filter returns <= 16; try 14 (or 20 or 24?) if use INT32
construct circular buffer of aywid intermediate scanlines
prepare a WeightBuf (a sampled filter for source pixels) for each dest x position
y mapping will be done in place
prepare a WeightBuf for this dest y position
If unscaled, integer translation, and filter is cardinal then we can really cruise:
Filter Registration



				name			evalFunc	   support  blur	window	cardinal	unitrange

note: we window Sinc and Bessel by default, but not Gaussian or Normal
Filters
unit area filters for unit spaced samples

box filter, a.k.a. pulse, Fourier window, or 1st order (constant) b-spline basis
triangle filter, a.k.a. Bartlett window, or 2nd order (linear) b-spline basis
3rd order (quadratic) b-spline basis
4th order (cubic) b-spline basis
Catmull-Rom spline basis, a.k.a. Overhauser and Parabolic Blend
Gaussian (infinite)
Sinc filter, a perfect lowpass filter (infinite)
Bessel filter (for circularly symmetric 2-d filtering, infinite)
See William Pratt, Digital Image Processing, p. 97
zeros are at approx x = 1.2197, 2.2331, 3.2383, 4.2411, 5.2428, 6.2439, 7.2448, 8.2454
filters for non-unit spaced samples
normal distribution: infinite, and has unit area, but it's not for unit spaced samples
Normal(x) = Gaussian(x/2)/2
Scanline Operations
Read scanline y0 for x in [x0..x0+buf.len) from maps, splitting the red, green, and blue components into buf, using supplied scratch space.  Writes to only buf[i].r if maps is B&W.
Write scanline y0 for x in [x0..x0+npixel) to maps, packing the red, green, and blue components from buf, using given scratch space.  Reads from only buf[i].r if maps is B&W.
Resample abuf into tbuf, which has length tbuf.len, using the specified filter and filter coefficient table xfilt.  tbuf and xfilt both have length tbuf.len.  B&W version.
sum _ sum+LongMult[xfilt[bx][axf-i0], abuf[axf]];
sum _ sum+INT32[xfilt[bx][axf-i0]]*abuf[axf];
Resample abuf, into tbuf, which has length tbuf.len, using the specified filter and filter coefficient table xfilt.  tbuf and xfilt both have length tbuf.len.  RGB version.
Zero out the B&W accum, which has length accum.len.
Zero out the RGB accum, which has length accum.len.
Perform the array operation: accum _ accum+weight*tbuf, where the B&W accum and tbuf have length tbuf.len.
accum[bx] _ accum[bx]+LongMult[weight, tbuf[bx]];
accum[bx] _ accum[bx]+INT32[weight]*tbuf[bx];
Perform the array operation: accum _ accum+weight*tbuf, where the B&W accum and tbuf have length tbuf.len.
Perform the array operation: bbuf _ rightshift[accum, shift], where the B&W accum and bbuf have length bbuf.len.
Perform the array operation: bbuf _ rightshift[accum, shift], where the RGB accum and bbuf have length bbuf.len.
Miscellaneous Support Routines


Start Code
�Ê$´–
"cedarcode" style•NewlineDelimiter

™™Jšœ	Ïeœ6™BJ™8J™XJ™-—J˜JšÏk	œTžœ žœ
˜‚
J˜šÑblnœžœ
ž˜JšžœJžœ ž˜uJšžœ	˜—J˜Jšœž˜
headš
Ïl™Jšœ
žœ˜"Jšœžœ˜(Jšœ
žœ˜)Jšœžœ˜&Jšœžœ˜0Jšœ
žœ˜+Jšœžœ'˜<Jšžœžœžœ
žœ
˜Jšžœžœžœ
˜Jšœ	žœžœ˜Jšœ	žœžœ˜—š
 ™Jšœžœ˜#Jšœžœ˜0Jšœžœ˜0Jšœžœ˜2—š
 ™Jšœžœ˜+Jšœžœ˜0Jšœžœ˜.Jšœžœ˜3Jšœ
žœ˜&Jšœžœ˜+Jšœžœ˜!Jšœžœ˜)Jšœžœ˜.Jšœžœ˜$Jšœ
žœ˜&Jšœžœ˜+Jšœžœ˜(Jšœžœ˜-Jšœ
žœ˜'Jšœžœ˜,Jšœžœ˜,Jšœžœ˜1Jšœžœ˜$Jšœžœ˜)Jšœ
žœ˜'—š
 ™Jš
Ïnœžœ
žœ
žœžœ
˜*—š
 ™š¡œžœ
žœ*˜=Jšœ
žœžœ˜FJšœžœ-˜;J˜SJš	œžœžœ
žœžœ˜7Jš	œžœžœ
žœžœ˜7˜1J˜?J˜.J˜?—J˜——š
 ™š¡
œžœ˜0J˜—š
¡œžœ
žœ;žœžœ˜aJ™*Jšœžœ
˜Jšœžœ
˜Jšœžœ
˜Jšœžœ
˜Jšœžœžœ"Ïc˜OJšœžœžœ ˜/Jš	œ,Ïs
ϣ
ϣ
ϣ
œ˜bJ˜J™—š
¡œžœ
žœNžœžœ˜~J˜J™X™\J™—šœžœ
žœžœ˜FJ™6J™F—Jšœžœžœ
˜+Jšœžœžœ
˜+Jšœžœžœ
˜+Jšœžœžœ
˜+Jšœžœ¢˜/Jšœžœ˜Jšœžœ¢˜2Jšœžœ
˜Jšœžœ
žœ¢˜Jšœžœ
žœ˜'J™BJ™U—Jšœžœžœ¢˜JJšœžœžœ˜+Jšœ	žœ¢Ð
˜í
Jšœžœ¢!˜EJ˜J™%Jšœžœ/˜NJšœžœ/˜NJ˜0J˜0J˜2J˜2J˜š
žœžœžœžœžœ˜>J™ršžœžœ
žœ˜J˜2J˜,Jšžœ˜J˜—˜JšžœD˜F—J˜—˜
Jšœžœ¢˜3Jšœžœ˜Jšœžœ¢˜8Jšœžœ˜Jš
œ	žœžœ	žœ	žœ
¢˜SJšœžœ¢3˜SJšœžœ˜Jšœžœ¢˜7J˜J™Vš
¡œžœžœžœžœžœž˜AJšœ
žœžœ˜)—š¡œžœžœžœžœž˜;Jšœ
žœ˜J˜—š
¡œžœžœžœžœ˜`J˜J™?™/™J™J™J™'J™CJ™——Jšœ	Ïbœ™"J™Jšœžœ˜Jšœžœ
˜
J˜J™Bšžœ¢˜1šžœ˜šžœ˜
šžœ¢$˜3Jšœžœ	žœžœ˜!J™C—šžœ¢˜$Jš
¤™šžœ¢
˜šžœ;˜?Jš
œ|™|—šžœ
¢	˜Jšœžœ
˜
Jšœžœ
˜
Jš
œ]™]˜7Jš
œ¿
™¿
—Jšœžœžœ¢˜=Jšœ#¢˜6Jšœ¢˜.šžœžœžœ˜BJ™9—J˜——Jšžœžœ
˜Jšžœ
žœžœ˜#Jš
žœžœ
žœ
žœžœ˜AJ˜J˜——šžœ¢$˜6šžœ˜Jšžœ	žœžœžœ
˜>Jš
žœž
œ
žœ
žœžœ
˜DJ˜
—šžœ¢'˜1Jš
žœž
œ
žœžœžœ
˜CJšžœ
žœžœžœ
˜AJ˜——J˜
—šžœ¢9˜Dšžœžœ	žœžœ
˜5J˜———Jšœ¢Ðbc˜J˜—š¡œžœžœ˜$J™*J™QJ™Jšœžœ
˜Jš
œ"˜"Jšœžœ˜'J˜7Jšœžœ˜"Jšœ0¢˜GJšœ0¢˜Gšžœžœ˜
Jšœžœ˜"Jšœžœ˜"Jšœ1¢˜HJšœ1¢˜HJšœ1¢˜HJšœ1¢˜HJ˜—J˜šžœžœ
žœ˜Jšžœ˜Jšžœ"žœžœ˜>Jšžœ"žœžœ˜>J˜—J˜Jšœ8¢˜Qšžœžœ
žœ¢%˜BJšœžœ˜Jšœ¢<˜Tšžœžœ˜
Jšœ¢9˜QJšœ¢9˜QJ˜—Jšžœ
˜—J˜Jšœ
¢˜'šžœžœ
žœ¢˜9J˜J˜J˜J™šžœžœ˜Jšœžœ
žœ
žœ˜/Jš
œžœ
žœ
žœ
žœ
žœ
žœ˜?J˜J˜šžœžœ˜
J˜J˜J˜J˜J˜—J˜š
œ™Jš
œB˜Bšžœžœ˜
Jš
œC˜CJš
œC˜CJ˜——šžœ˜š
ž˜šžœžœ
žœ¢˜6Jšœ¢"˜;Jšœ¢"˜;Jšœ¢"˜;Jšœžœ	˜Jšœžœ	˜Jšœžœ	˜Jš	œžœ
žœ
žœ
žœ	˜-Jš	œžœ
žœ
žœ
žœ	˜-Jš	œžœ
žœ
žœ
žœ	˜-Jš
ž˜——š
ž˜šžœžœ
žœ¢˜6Jšœ¢"˜;Jšœžœ	˜Jš	œžœ
žœ
žœ
žœ	˜-Jšžœ
˜———J˜—JšœC¢
˜Pšžœžœ˜
JšœE¢˜MJšœE¢˜MJ˜—Jšžœ
˜J˜—šžœžœ
žœžœ˜$Jšžœ
žœ*˜3šžœ
žœ˜šœžœžœ
ž˜J˜,Jšžœ
˜———Jšœ¢¥˜J˜—š¡œžœ˜J˜š¡
œžœ	žœžœžœžœžœ˜oJšœžœ˜Jšœžœ˜Jšœžœ˜Jšœžœ˜Jšœžœ
˜	Jšžœžœ˜Jšžœ
žœ
˜Jšžœ
žœ!˜1Jšžœžœ
žœžœ˜3J˜
J˜
J™-šžœžœ
žœ
ž˜Jšœ)žœ#˜PJšžœ
˜—Jš
œ-™-Jš	œžœ
žœžœ
žœ˜CJ˜šžœ
žœžœ
žœ
˜šžœ
žœžœ
žœ
žœ
¢1˜PJšœžœ	žœ
žœ
žœžœ
žœ
žœžœ%˜oJ˜J˜Jš
ž˜—šžœ
žœžœ
žœ
žœ¢˜/Jšœžœžœ#˜JJ˜J˜Jšžœ
˜——šžœ	žœžœ
žœ
žœžœ	™BJ™6—šžœžœ˜Jšœžœžœžœ˜(J˜&J˜—Jšžœ˜
Jšœ¢¥˜J˜—Jšœžœ¢˜>Jšœžœ¢˜>Jšœ:žœ
¢˜dJšœžœ¢˜DJšœžœ&¢˜JJšœžœ#˜-Jšœžœ˜'Jšœžœ¢˜VJšœžœžœ
¢˜UJšœžœ¢$˜QJšœžœ¢˜NJšœžœ8ž™XJšœ
žœ˜Jšœžœ
žœ˜3Jšœžœ¢$˜KJ˜J˜J˜J˜šžœžœ
žœž˜Jš
œ9™9Jšœžœ˜'J˜Jšœ¢˜4Jšžœ
˜—J˜šžœžœ
žœ˜Jšžœžœžœ˜Kšžœ:˜<Jšžœ
žœžœžœžœžœ˜P—šžœ:˜<Jšžœ
žœžœžœžœžœ˜P—J˜—J˜J™Qšžœžœ
žœ
ž˜JšœMžœ˜RJšžœ
˜J˜—˜;Jš
œ™—J˜š	žœžœ
žœ
žœ¢˜;Jšœžœ
˜J˜J˜J™,Jšœ£
œ
£
œJ˜Ošžœ˜Jšžœ˜Jšžœ˜J˜—š	žœžœ
žœžœ
¢<˜jJšœžœ˜*šžœžœ¢˜;šžœ	žœž˜(Jšœžœ!žœ˜;—J˜
Jšœ1¢˜Lšžœ˜Jšžœ-˜1Jšžœ-˜1—Jš
œ˜J˜—šžœ˜Jšžœ=˜AJšžœ=˜AJ˜—Jšžœ
˜J˜—šžœ˜Jšžœ0˜4Jšžœ0˜4J˜—Jšœ3¢˜OJšœžœ
˜Jšžœ
˜J˜—J˜,Jšœ¢¥˜—J˜šžœžœ
žœ˜J˜(J˜)Jšžœ˜J˜—Jšžœ
žœ
žœ
žœ
žœ
žœ
˜<Jšžœžœžœ˜AJ˜J™Sšžœžœžœžœžœ
žœžœ
žœ˜šžœ¢"˜2J•StartOfExpansionÜ
[map: 
ImagerSample.SampleMap, dstMin: 
SF.Vec _ [s: 
0, f: 
0], srcMin: 
SF.Vec _ [s: 
0, f: 
0], size: 
SF.Vec _ [s: 
32767, f: 
32767], function: 
ImagerSample.Function _ [dstFunc: 
null, srcFunc: 
null]]˜7—šžœ¢$˜4šžœ žœ
˜(Jšžœ˜Jšžœ˜———J˜—Jšœ¢¥˜J˜—š
¡œžœ
žœ9žœžœ˜bJ˜$J˜&JšœOžœ˜VJ˜J˜$J˜——š
 ™Jšœžœ˜%šœžœ˜+J˜—Jšœžœ¢˜6J˜š
¡
œžœ
žœ
žœžœ˜GJ˜—š
¡
œžœ
žœžœžœ
˜9Jšœžœ˜!šžœ˜Jšžœžœžœ
˜:Jšžœ
žœ˜—J˜J˜—š¡
œžœ
žœžœ˜*Jšžœ#˜%šžœžœ
žœž˜+J˜!Jšžœžœžœžœžœžœžœ˜�
Jšžœžœ
žœ$˜@Jšžœ
˜—J˜J™—š¡œžœ
žœ˜Jšœžœ¢˜"Jšœ¢8˜SJšœ	žœ¢>˜QJšœžœ¢-˜BJšœ
žœ¢*˜=Jšœ
žœ¢œ¢;˜yJšœžœžœ¢,˜FJšœžœ¢3˜SJšœžœ
¢˜=Jšœžœ¢;˜QJšœžœ
žœžœ¢'˜BJ˜Jšžœžœ
žœžœ˜D˜
Jšœžœ¢˜@Jšžœžœ¢˜Fšžœ˜ Jšžœ9˜=—J˜{Jšžœžœ˜?J˜—J˜J˜—š¡œžœžœžœ˜[Jšœžœžœ žœ˜IJšœžœ˜ Jš
žœžœ
žœžœžœ
˜>J˜J˜—J™š¡œžœ
žœ˜$Jšžœžœ
žœ˜7J˜J™—š¡œžœ
žœžœ˜.Jšœžœ˜!šžœžœ˜šžœžœ
žœž˜-J˜(Jšžœ
˜—J˜0J˜—J˜J˜—š¡œžœ
žœ˜'JšœÏz:™>J™Jš	œžœžœžœžœ˜?Jšœ-žœžœžœ˜AJšœ3žœžœžœ˜GJšœ5žœžœžœ˜IJšœ/žœžœžœ˜CJ˜Jšœ1žœžœžœ˜FJšœ3žœžœžœ˜GJšœ.žœžœžœ˜CJšœ3žœžœžœ˜HJ˜Jšœ1žœžœžœ˜FJ˜J™FJ˜J˜—š
¡œžœžœžœžœ˜6Jšžœžœ
žœ
žœ˜(šžœžœ
žœž˜+Jšžœ*žœžœ
žœ˜BJšžœ
˜—Jšžœ˜J˜——š
 ™šœžœ
˜J˜—š
 )™)J™š¡œžœ˜"J™Pšžœž˜
Jšœžœ˜Jšœžœ˜Jšžœžœ˜—˜J˜——š¡œžœ˜'J™Mšžœž˜
Jšœžœ˜Jšœžœ˜Jšœžœ˜Jšžœžœ˜—˜J˜——š¡
œžœ˜(J™$šžœž˜
Jšœ	žœ˜Jšœ
žœ
žœ˜,Jšœ	žœ
˜Jšœ
žœ
žœ˜,Jšžœžœ˜—˜J˜——š¡	œžœ˜$J™ šžœž˜
Jšœžœ˜Jšœžœ	žœ
˜,Jšœžœ˜(Jšœžœ˜'Jšœžœ	žœ
˜,Jšžœžœ˜—˜J˜——š¡
œžœ˜%J™?šžœž˜
Jšœžœ˜Jšœžœ˜(Jšœžœ˜(Jšœžœ˜'Jšœžœ˜)Jšžœžœ˜—˜J˜——š¡œžœ˜'J™Jšžœ+˜1˜J˜——š¡œžœ˜#J™0Jšžœ
žœžœžœ˜8˜J˜——š¡
œžœ˜%J™@Jšœ¦œ™2J™VJš
žœ
žœžœžœ	žœ˜:J˜—J˜—š
 #™#J˜š¡
œžœ˜%J™VJ™Jšžœ+˜1J˜———š
 ™š¡œ
£
ž£
ϣ
ϣ
ϣ
ϣ
žœ
£
ϣ
œ
£
ϣ
œ	£
œ
£
œ
˜\J™´
Jšœžœ˜šžœžœ
žœž˜$JšœE¢˜Všžœž˜˜šžœžœ
žœ
ž˜J˜Jšžœ
˜——˜	šžœžœ
žœ
ž˜J˜Jšžœ
˜——˜šžœžœ
žœ
ž˜J˜Jšžœ
˜——˜šžœžœ
žœ
ž˜J˜Jšžœ
˜——Jšžœžœ
˜—Jšžœ
˜—J˜J˜—š¡
œ
£
ž£
ϣ
ϣ
ϣ
ϣ
žœ
£
ϣ
œ
£
ϣ
œ	£
œ
£
œ
˜]J™®
Jšœžœ˜šžœžœ
žœž˜$šžœž˜Jšœžœžœ
žœ
žœžœ
˜EJšœ	žœžœ
žœ
žœžœ
˜FJšœ	žœžœ
žœ
žœžœ
˜FJšœžœžœ
žœ
žœžœ
˜DJšžœžœ
˜—JšœC¢˜UJšžœ
˜—J˜J˜—š¡œžœG˜_J™«
šžœžœ
žœž˜šžœ˜	Jšœžœ˜Jšœžœ˜Jš	œžœ
žœ
žœ
žœ˜.Jšœžœ
žœ
žœ˜)Jšœžœ˜šžœ˜šžœ	¢4˜Ašžœžœ
žœ
ž˜J™1J˜%Jšœžœ
žœ˜Jšœžœ˜Jš
ž˜——šžœ	¢/˜<šžœžœ
žœ
ž˜Jšœ
žœ™-Jšœ
žœ˜!Jšœžœ
žœ˜Jšœžœ˜Jšžœ
˜———Jšœ
£
œ
£
žœ£
œ
£
œ
£
ž£
œžœ£
ž£
œ˜ZJ˜—Jšžœ
˜—J˜J˜—š¡œ
£
žœG˜`J™¬
šžœžœ
žœž˜šžœ˜	Jšœžœ˜Jšœžœ˜Jš	œžœ
žœ
žœ
žœ˜.Jšœžœ
žœ
žœ˜)Jšœžœ˜Jšœžœ˜Jšœžœ˜šžœ˜Jšžœ	¢4˜Ašžœžœ
žœ
ž˜J˜'J˜'J˜'Jšœžœ
žœ˜Jšœžœ˜Jš
ž˜—šžœ	¢/˜<šžœžœ
žœ
ž˜Jšœžœ˜J˜J˜J˜Jšœžœ
žœ˜Jšœžœ˜Jšžœ
˜———J˜J˜J˜J˜—Jšžœ
˜—J˜J˜—š¡œžœ˜,J™3šžœžœ
žœž˜ J˜Jš
ž˜—J˜J˜—š¡œžœ˜-J™3šžœžœ
žœž˜ J˜J˜J˜Jšžœ
˜—J˜J˜—š¡œžœžœ&˜XJšœžœ˜J™jJšœžœ
žœ
žœ˜,Jšœžœ
žœ
žœ˜*Jšœžœ˜šžœ˜šžœ	¢4˜Ašžœžœ
žœ
ž˜J™1J˜0Jšœžœ
˜Jšœžœ˜Jš
ž˜——šžœ	¢/˜<šžœžœ
žœ
ž˜Jšœžœ™-Jšœžœ˜*Jšœžœ
˜Jšœžœ˜Jšžœ
˜———J˜J˜—š¡œžœžœ&˜YJšœžœ˜J™jJšœžœ
žœ
žœ˜,Jšœžœ
žœ
žœ˜*Jšœžœ˜šžœ˜šžœ	¢4˜Ašžœžœ
žœ
ž˜J˜.J˜.J˜.Jšœžœ
˜Jšœžœ˜Jš
ž˜——šžœ
¢/˜=Jšœžœ
˜šžœžœ
žœ
ž˜J˜J˜J˜Jšœžœ
˜Jšœžœ˜Jšžœ
˜—J˜——J˜J˜—š¡œžœžœ(˜WJ™pJšœžœ˜Jšœžœžœž	œ
žœ
žœ˜Gšžœ˜š
ž˜šžœžœ
žœ
ž˜Jš	œ
žœž	œ
žœžœ˜MJš
ž˜——š
ž˜šžœžœ
žœ
ž˜Jšœžœ˜Jšœ
žœžœžœ
žœžœžœ
žœ
žœž	œ
žœ
žœ˜{Jšžœ
˜———˜J˜——š¡œžœžœ(˜XJ™pJšœžœ˜Jšœžœžœž	œ
žœ
žœ˜Gšžœ˜š
ž˜šžœžœ
žœ
ž˜Jš	œ
žœž	œ
žœžœ˜MJš	œ
žœž	œ
žœžœ˜MJš	œ
žœž	œ
žœžœ˜MJš
ž˜——š
ž˜šžœžœ
žœ
ž˜Jšœžœ˜Jšœ
žœžœžœ
žœžœžœ
žœ
žœž	œ
žœ
žœ˜{J˜Jšœ
žœžœžœ
žœžœžœ
žœ
žœž	œ
žœ
žœ˜{J˜Jšœ
žœžœžœ
žœžœžœ
žœ
žœž	œ
žœ
žœ˜{Jšžœ
˜———J˜——š
 ™š¡œžœžœžœžœžœžœ˜EJ™—š¡
œžœžœžœžœžœžœ˜GJ™—š¡œžœžœžœžœž˜.šœ
žœ˜J˜——š¡œžœžœžœžœžœžœ˜GJ˜—š
¡œžœžœžœžœ˜@šžœ˜Jšžœ
žœžœ
žœ
žœ
žœ˜D—Jšžœžœ
žœ
˜2J˜——š
 
™
Jšœ¢'˜C—J˜Jšžœ
˜J˜—…—hØ¢B