DIRECTORY Real USING [RoundI], RealFns USING [SinDeg], Filters; FiltersImpl: CEDAR PROGRAM IMPORTS Real, RealFns 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]]]; }; 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; 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 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]]; 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; 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; }; UncoveredRatio: PUBLIC PROC[ centerA, widthA, centerB, widthB: REAL, afterPixelA: BOOLEAN] RETURNS[ weight: REAL _ 0.0] ~ { 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 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]) }; }; 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; }; 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; }; 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]]]}; ENDCASE; }; GetFilterWidth: PUBLIC PROC[] RETURNS[ INTEGER] ~ { RETURN[ filterWidth]}; GetFilterRecLength: PUBLIC PROC[] RETURNS[ INTEGER] ~ { RETURN[ filterRecLength]}; 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* SetFilter[ 0]; -- initial default END. âFiltersImpl.mesa Last Edited By: Marty Hiller, September 6, 1984 11:57:09 pm PDT Copyright (C) 1984 by Xerox Corporation. All rights reserved. ImagerPixelMaps USING [PixelMap], ImagerPixelMapsExtras USING [SetPixel], -- , ImagerPixelMapsExtras Beginning of PixelWeight find how much of pixel A overlaps (or doesn't overlap) the pixel B find relative beginning and end of pixels if the end of one pixel filter is inside the other pixel filter, index the beginning (or end) of filter subdivide filter to get more accurate areas multiply to get overlap between the pixels 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 find how much of pixel A overlaps (or doesn't overlap) the pixel B if the point is outside the pixel, index the beginning or end of the filter 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 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) 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; }; 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]]]]]]]]]}; 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]; }; Initialize all filters; initial default filter is sin x / x to 1 pixel over halfFilterRecLength _ 5*filterPrecision; filterRecLength _ 2*halfFilterRecLength; Filter[ 7] _ Filter7[]; -- sin x / x to an absurd number of pixels over Ê §˜šœR™RJšœ>™>J™—šÏk ˜ Jšœœ™#Jšœœ™'Jšœ œ ˜Jšœœ ˜Jšœ˜J˜—šœœ˜Jšœ˜Jšœ™Jšœ˜—J˜Jšœ˜J˜Jšœ&œ˜.Jšœœ˜Jšœ œ˜Jšœœœ˜%Jš œœœœ œœœ˜EJšœ˜Jšœœœ˜#J˜Jšœœ˜Jš œœœœœ˜$Jšœœœœ˜,Jšœœœœ˜,J˜šÏnœœœ$œ˜CJšœœ œ˜!J˜Jšœœ˜Jšœœœœ˜"J˜šžœœœ+œ˜[Jšœœœœ ˜6Jšœœœœ˜/Jšœœœ)˜:J˜—J˜Jšœ™J˜JšœB™BJšœ-œ˜2Jšœ$œ˜,J˜Jšœœœœ˜,Jšœ œ ˜Jšœ œ ˜Jšœ˜Jšœ˜J˜Jšœ)™)Jšœ-˜-Jšœ-˜-J˜š œœœœ˜PJšœ Ïc4˜C—J˜Jšœg™gJšœ"˜"Jšœ"˜"šœœ˜Jšœ˜Jšœ!˜!—šœ˜Jšœ˜Jšœ!˜!J˜—šœœ˜Jšœ˜Jšœ,˜,—šœ˜Jšœ˜Jšœ-˜-J˜—Jšœœœœ˜FJ˜Jšœ+™+Jšœ˜šœœ˜šœœ˜/šœd˜dJ˜——šœœ˜/Jšœd˜d—Jšœ˜Jšœ˜J˜—Jšœ*™*šœœ˜Jšœ$˜$Jšœ)˜)Jšœ”˜”Jšœ˜Jšœ˜—J˜—˜Jšœ¶™¶—š žœœœ(œœ˜dJšœœ œ˜!J˜JšœB™BJšœœ˜Jšœœ˜#J˜Jšœœœ˜Jšœœœ ˜1Jšœœœ ˜1J˜Jšœ œŸ˜AJšœ œ ˜Jšœ˜J˜Jšœ"˜"J˜šœ ˜Jšœ,˜,—š˜Jšœ-˜-—J˜JšœK™Kšœ˜Jšœ˜—šœœ˜Jšœ˜—š˜Jšœ:˜:J˜—šœ˜Jšœ˜—šœœ˜Jšœ˜—š˜Jšœ:˜:—J˜šœ œ˜Jšœ˜Jšœœ˜œJ˜—šœ˜Jšœ˜Jšœœ˜œJ˜—J˜J˜—J˜Jšœ1™1šžœœ™Jšœœœ)™KJšœœ™#Jšœœ™%Jšœœ™&šœœœœ™Jšœ™Jšœ™—šœœœœ™Jšœ™Jšœ™—šœœœœ™*Jšœ™Jšœ™J™—šœœœ™)Jšœ2™2Jšœ™—Jšœ™—J™Jšœ ™ šžœœ™Jšœœœ)™KJšœ™šœœœœ™.Jšœ™Jšœ™—šœœœ(œ™Jšœ5™5Jšœ™—J™šœœœœ™*Jšœ2™2Jšœ™—Jšœ™—J™Jšœ™šžœœœ™Jšœœœ)™KJšœ™š œœœœŸ™QJšœœ9™@Jšœ;œ™CJšœ™—JšœB™Bšœœœ(œ™