FiltersImpl.mesa
Last Edited By: Marty Hiller, September 6, 1984 11:57:09 pm PDT
Copyright (C) 1984 by Xerox Corporation. All rights reserved.
DIRECTORY
ImagerPixelMaps USING [PixelMap],
ImagerPixelMapsExtras USING [SetPixel],
Real USING [RoundI],
RealFns USING [SinDeg],
Filters;
FiltersImpl:
CEDAR
PROGRAM
IMPORTS Real, RealFns
-- , ImagerPixelMapsExtras
EXPORTS Filters
= BEGIN
halfFilterRecLength, filterRecLength: INTEGER;
filterPrecision: INTEGER ← 250;
filterWidth: INTEGER;
FilterType: TYPE ~ REF FilterTypeRec;
FilterTypeRec: TYPE ~ RECORD[ element: SEQUENCE length: NAT OF REAL];
defaultFilter: [0..2);
Filter: ARRAY [0..2) OF FilterType;
PI: REAL = 3.14159265;
splitIndices: LIST OF INTEGER ← NIL;
splitListItem: TYPE ~ RECORD[ x,y: INTEGER];
splitListType: TYPE ~ LIST OF splitListItem;
PixelWeight:
PUBLIC
PROC[ centerA, widthA, centerB, widthB:
REAL]
RETURNS[ weight: REAL ← 0.0] ~ {
splitList: splitListType ← NIL;
s: LIST OF INTEGER ← splitIndices;
SplitAt:
PRIVATE
PROC[ list: splitListType, add: splitListItem]
RETURNS[ splitListType] ~ {
IF list.first.x > add.x THEN RETURN[ CONS[ add, list]]
ELSE IF list.first.x = add.x THEN RETURN[ list]
ELSE RETURN[ CONS[ list.first, SplitAt[ list.rest, add]]];
};
Beginning of PixelWeight
find how much of pixel A overlaps (or doesn't overlap) the pixel B
bA, eA, kA, kB, sizeRatio, width2A, width2B: REAL;
bIndexA, eIndexA, bIndexB, eIndexB: INTEGER;
IF widthA = 0 OR widthB = 0 THEN RETURN[ 0];
widthA ← ABS[ widthA];
widthB ← ABS[ widthB];
width2A ← 2*widthA;
width2B ← 2*widthB;
find relative beginning and end of pixels
bA ← (centerB - widthB) - (centerA - widthA);
eA ← (centerA + widthA) - (centerB + widthB);
IF eA > width2A
OR
(-eA > width2B
OR
(bA > width2A
OR
-bA > width2B))
THEN
RETURN[ 0]; -- distance between them is greater than their width
if the end of one pixel filter is inside the other pixel filter, index the beginning (or end) of filter
kA ← halfFilterRecLength / widthA;
kB ← halfFilterRecLength / widthB;
IF bA < 0
THEN {
bIndexA ← 0;
bIndexB ← Real.RoundI[ -bA * kB]}
ELSE {
bIndexB ← 0;
bIndexA ← Real.RoundI[ bA * kA]};
IF eA < 0
THEN {
eIndexA ← filterRecLength;
eIndexB ← Real.RoundI[ (width2B + eA) * kB]}
ELSE {
eIndexB ← filterRecLength;
eIndexA ← Real.RoundI[ (width2A - eA) * kA]};
splitList ← CONS[ [bIndexA, bIndexB], CONS[ [eIndexA, eIndexB], NIL]];
subdivide filter to get more accurate areas
sizeRatio ← widthA / widthB;
WHILE s #
NIL
DO
IF bIndexA < s.first
AND eIndexA > s.first
THEN
splitList ← SplitAt[ splitList, [s.first, Real.RoundI[ (s.first - bIndexA) * sizeRatio] + bIndexB]];
IF bIndexB < s.first
AND eIndexB > s.first
THEN
splitList ← SplitAt[ splitList, [Real.RoundI[ (s.first - bIndexB) / sizeRatio] + bIndexA, s.first]];
s ← s.rest;
ENDLOOP;
multiply to get overlap between the pixels
WHILE splitList.rest #
NIL
DO
p1: splitListItem ← splitList.first;
p2: splitListItem ← splitList.rest.first;
weight ← weight + (Filter[defaultFilter][ p2.x] - Filter[defaultFilter][ p1.x]) *
(Filter[defaultFilter][ p2.y] - Filter[defaultFilter][ p1.y]);
splitList ← splitList.rest;
ENDLOOP;
};
Pixel A is the covering pixel, pixel B is the covered one. This tells you how much of the pixel is left to cover, either before or after the covering pixel, depending on afterPixelA
UncoveredRatio:
PUBLIC
PROC[
centerA, widthA, centerB, widthB:
REAL,
afterPixelA:
BOOLEAN]
RETURNS[ weight: REAL ← 0.0] ~ {
find how much of pixel A overlaps (or doesn't overlap) the pixel B
cA, eA, width2B: REAL;
cIndexA, eIndexA, eIndexB: INTEGER;
IF widthB = 0 THEN RETURN[ 0];
IF widthA < 0 THEN afterPixelA ← NOT afterPixelA;
IF widthB < 0 THEN afterPixelA ← NOT afterPixelA;
widthA ← ABS[ widthA/filterWidth]; -- distance to the next pixel
widthB ← ABS[ widthB];
width2B ← 2*widthB;
cA ← centerA - (centerB - widthB);
IF afterPixelA
THEN
eA ← (centerA + widthA) - (centerB - widthB)
ELSE
eA ← (centerA - widthA) - (centerB - widthB);
if the point is outside the pixel, index the beginning or end of the filter
IF cA < 0
THEN
cIndexA ← 0
ELSE
IF cA > width2B
THEN
cIndexA ← filterRecLength
ELSE
cIndexA ← Real.RoundI[ cA * halfFilterRecLength / widthB];
IF eA < 0
THEN
eIndexA ← 0
ELSE
IF eA > width2B
THEN
eIndexA ← filterRecLength
ELSE
eIndexA ← Real.RoundI[ eA * halfFilterRecLength / widthB];
IF afterPixelA
THEN {
eIndexB ← filterRecLength;
weight ← (Filter[defaultFilter][ eIndexA] - Filter[defaultFilter][ cIndexA]) / 2 +
(Filter[defaultFilter][ eIndexB] - Filter[defaultFilter][ eIndexA])
}
ELSE {
eIndexB ← 0;
weight ← (Filter[defaultFilter][ cIndexA] - Filter[defaultFilter][ eIndexA]) / 2 +
(Filter[defaultFilter][ eIndexA] - Filter[defaultFilter][ eIndexB])
};
};
square; extends with zero value to 1st pixel over
Filter0: PRIVATE PROC
RETURNS[ filter: FilterType ← NEW[ FilterTypeRec[ filterRecLength+1]]] ~ {
j: INTEGER ← halfFilterRecLength/2;
k: INTEGER ← 3*halfFilterRecLength/2;
inc: REAL ← 1.0 / halfFilterRecLength;
FOR i: INTEGER IN [0..j] DO
filter[i] ← 0;
ENDLOOP;
FOR i: INTEGER IN (j..k] DO
filter[i] ← filter[ i-1] + inc;
ENDLOOP;
FOR i: INTEGER IN (k..filterRecLength] DO
filter[i] ← filter[ k];
ENDLOOP;
FOR i: INTEGER IN [0..filterRecLength] DO
filter[i] ← filter[ i] / filter[ filterRecLength];
ENDLOOP;
};
triangular
Filter1: PRIVATE PROC
RETURNS[ filter: FilterType ← NEW[ FilterTypeRec[ filterRecLength+1]]] ~ {
filter[0] ← 0;
FOR i: INTEGER IN [1..halfFilterRecLength] DO
filter[i] ← filter[ i-1] + i;
ENDLOOP;
FOR i: INTEGER IN (halfFilterRecLength..filterRecLength] DO
filter[i] ← filter[ i-1] + (filterRecLength - i);
ENDLOOP;
FOR i: INTEGER IN [0..filterRecLength] DO
filter[i] ← filter[ i] / filter[ filterRecLength];
ENDLOOP;
};
sin x / x to first pixel
Filter2:
PRIVATE
PROC
RETURNS[ filter: FilterType ← NEW[ FilterTypeRec[ filterRecLength+1]]] ~ {
filter[0] ← 0;
FOR i:
INTEGER
IN [1..halfFilterRecLength)
DO
-- calculate integral of pixel
x: REAL ← (i-halfFilterRecLength) / (halfFilterRecLength * 1.0);
filter[i] ← filter[ i-1] + RealFns.SinDeg[ 180 * x] / (PI * x);
ENDLOOP;
filter[halfFilterRecLength] ← filter[halfFilterRecLength - 1] + 1;
FOR i:
INTEGER
IN (halfFilterRecLength..filterRecLength]
DO
x: REAL ← (i-halfFilterRecLength) / (halfFilterRecLength * 1.0);
filter[i] ← filter[ i-1] + RealFns.SinDeg[ 180 * x] / (PI * x);
ENDLOOP;
FOR i:
INTEGER
IN [0..filterRecLength]
DO
filter[i] ← filter[ i] / filter[ filterRecLength];
ENDLOOP;
};
sin x / x to second pixel
Filter3: PRIVATE PROC
RETURNS[ filter: FilterType ← NEW[ FilterTypeRec[ filterRecLength+1]]] ~ {
filter[0] ← 0;
FOR i: INTEGER IN [1..halfFilterRecLength) DO -- calculate integral of pixel
x: REAL ← (i-halfFilterRecLength) / (halfFilterRecLength * 1.0);
filter[i] ← filter[ i-1] + RealFns.SinDeg[ 360 * x] / (2 * PI * x);
ENDLOOP;
filter[halfFilterRecLength] ← filter[halfFilterRecLength - 1] + 1;
FOR i: INTEGER IN (halfFilterRecLength..filterRecLength] DO
x: REAL ← (i-halfFilterRecLength) / (halfFilterRecLength * 1.0);
filter[i] ← filter[ i-1] + RealFns.SinDeg[ 360 * x] / (2 * PI * x);
ENDLOOP;
FOR i: INTEGER IN [0..filterRecLength] DO
filter[i] ← filter[ i] / filter[ filterRecLength];
ENDLOOP;
};
gaussian to second pixel
Filter4: PRIVATE PROC
RETURNS[ filter: FilterType ← NEW[ FilterTypeRec[ filterRecLength+1]]] ~ {
filter[0] ← 1 / RealFns.Exp[ 4];
FOR i: INTEGER IN [1..filterRecLength] DO -- calculate integral of pixel (to e**4 at 2nd pixel)
x: REAL ← (i-halfFilterRecLength) / (halfFilterRecLength * .5);
filter[i] ← filter[ i-1] + 1 / RealFns.Exp[ (x * x)];
ENDLOOP;
FOR i: INTEGER IN [0..filterRecLength] DO
filter[i] ← filter[ i] / filter[ filterRecLength];
ENDLOOP;
};
sin x / x * sin (x/2) / (x/2)
Filter5:
PRIVATE
PROC
RETURNS[ filter: FilterType ← NEW[ FilterTypeRec[ filterRecLength+1]]] ~ {
filter[0] ← 0;
FOR i:
INTEGER
IN [1..halfFilterRecLength)
DO
-- calculate integral of pixel
x: REAL ← (i-halfFilterRecLength) / (halfFilterRecLength * 1.0);
filter[i] ← filter[ i-1] + (RealFns.SinDeg[ 180 * x] / (PI * x)) *
(RealFns.SinDeg[ 360 * x] / (2 * PI * x));
ENDLOOP;
filter[halfFilterRecLength] ← filter[halfFilterRecLength - 1] + 1;
FOR i:
INTEGER
IN (halfFilterRecLength..filterRecLength]
DO
x: REAL ← (i-halfFilterRecLength) / (halfFilterRecLength * 1.0);
filter[i] ← filter[ i-1] + (RealFns.SinDeg[ 180 * x] / (PI * x)) *
(RealFns.SinDeg[ 360 * x] / (2 * PI * x));
ENDLOOP;
FOR i:
INTEGER
IN [0..filterRecLength]
DO
filter[i] ← filter[ i] / filter[ filterRecLength];
ENDLOOP;
};
pointy gaussian to second pixel
Filter6: PRIVATE PROC
RETURNS[ filter: FilterType ← NEW[ FilterTypeRec[ filterRecLength+1]]] ~ {
filter[0] ← 1 / RealFns.Exp[ 9];
FOR i: INTEGER IN [1..filterRecLength] DO -- calculate integral of pixel (to e**9 at 2nd pixel)
x: REAL ← (i-halfFilterRecLength) * 3.0 / halfFilterRecLength;
filter[i] ← filter[ i-1] + 1 / RealFns.Exp[ (x * x)];
ENDLOOP;
FOR i: INTEGER IN [0..filterRecLength] DO
filter[i] ← filter[ i] / filter[ filterRecLength];
ENDLOOP;
};
sin x / x to fifth pixel
Filter7: PRIVATE PROC
RETURNS[ filter: FilterType ← NEW[ FilterTypeRec[ filterRecLength+1]]] ~ {
filter[0] ← 0;
FOR i: INTEGER IN [1..halfFilterRecLength) DO -- calculate integral of pixel
x: REAL ← (i-halfFilterRecLength) / (halfFilterRecLength * 1.0);
filter[i] ← filter[ i-1] + RealFns.SinDeg[ 5*180 * x] / (5*PI * x);
ENDLOOP;
filter[halfFilterRecLength] ← filter[halfFilterRecLength - 1] + 1;
FOR i: INTEGER IN (halfFilterRecLength..filterRecLength] DO
x: REAL ← (i-halfFilterRecLength) / (halfFilterRecLength * 1.0);
filter[i] ← filter[ i-1] + RealFns.SinDeg[ 5*180 * x] / (5*PI * x);
ENDLOOP;
FOR i: INTEGER IN [0..filterRecLength] DO
filter[i] ← filter[ i] / filter[ filterRecLength];
ENDLOOP;
};
SetFilter:
PUBLIC
PROC[ i: [0..2)] ~ {
SELECT i
FROM
0 => filterWidth ← 1;
1 => filterWidth ← 2;
ENDCASE;
defaultFilter ← i;
halfFilterRecLength ← filterWidth*filterPrecision;
filterRecLength ← 2*halfFilterRecLength;
SELECT i
FROM
0 => { -- sinc (to 1 pixel over)
splitIndices ← CONS[ halfFilterRecLength, NIL]};
1 => {
splitIndices ← CONS[ halfFilterRecLength/2,
CONS[ halfFilterRecLength,
CONS[ 3*halfFilterRecLength/2, NIL]]]};
7 => {
spacing: INTEGER ← halfFilterRecLength/5;
splitIndices ← CONS[ spacing,
CONS[ 2*spacing,
CONS[ 3*spacing,
CONS[ 4*spacing,
CONS[ 5*spacing,
CONS[ 6*spacing,
CONS[ 7*spacing,
CONS[ 8*spacing,
CONS[ 9*spacing,
NIL]]]]]]]]]};
ENDCASE;
};
ShowFilter: PUBLIC PROC[ i: [0..7), pix: ImagerPixelMaps.PixelMap] ~ {
p: INTEGER ← 0;
save: INTEGER ← defaultFilter;
offset: INTEGER;
SetFilter[ i];
offset ← ((pix.fSize-pix.fMin)-halfFilterRecLength)/2+pix.fMin;
FOR x: INTEGER IN [0..halfFilterRecLength] DO
ImagerPixelMapsExtras.SetPixel[ pix, -- axis
(pix.sSize - 1) - 50,
offset + x,
75];
ImagerPixelMapsExtras.SetPixel[ pix, -- integral of filter shape
(pix.sSize - 1) - (100 + Filter[i][ x*2]/50),
x,
125];
ImagerPixelMapsExtras.SetPixel[ pix, -- filter shape
(pix.sSize - 1) - (50 + (Filter[i][ x*2] - p)),
offset + x,
175];
p ← Filter[i][ x*2];
ENDLOOP;
SetFilter[ save];
};
GetFilterWidth: PUBLIC PROC[] RETURNS[ INTEGER] ~ { RETURN[ filterWidth]};
GetFilterRecLength: PUBLIC PROC[] RETURNS[ INTEGER] ~ { RETURN[ filterRecLength]};
Initialize all filters; initial default filter is sin x / x to 1 pixel over
halfFilterRecLength ← filterPrecision;
filterRecLength ← 2*halfFilterRecLength;
Filter[ 0] ← Filter2[]; -- sin x / x to 1 pixel over *good*
halfFilterRecLength ← 2*filterPrecision;
filterRecLength ← 2*halfFilterRecLength;
Filter[ 1] ← Filter5[]; -- sin x / x * sin (x/2) / (x/2) *very good* *slow*
halfFilterRecLength ← 5*filterPrecision;
filterRecLength ← 2*halfFilterRecLength;
Filter[ 7] ← Filter7[]; -- sin x / x to an absurd number of pixels over
SetFilter[ 0]; -- initial default
END.