<> <> <<>> DIRECTORY <> <> 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]]]; }; <> <> 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]) }; }; <> <> << RETURNS[ filter: FilterType _ NEW[ FilterTypeRec[ filterRecLength+1]]] ~ {>> <> <> <> <> <> <> <> <> <> <> <> <> <<>> <> <> <> <<};>> <<>> <> <> << RETURNS[ filter: FilterType _ NEW[ FilterTypeRec[ filterRecLength+1]]] ~ {>> <> <> <> <> <> <> <> <> <> <> <<};>> <> 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; }; <> <> << RETURNS[ filter: FilterType _ NEW[ FilterTypeRec[ filterRecLength+1]]] ~ {>> <> <> <> <> <> <> <> <> <> <> <> <> <> <<};>> <<>> <> <> << RETURNS[ filter: FilterType _ NEW[ FilterTypeRec[ filterRecLength+1]]] ~ {>> <> <> <> <> <> <<>> <> <> <> <<};>> <> 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; }; <> <> << RETURNS[ filter: FilterType _ NEW[ FilterTypeRec[ filterRecLength+1]]] ~ {>> <> <> <> <> <> <<>> <> <> <> <<};>> <<>> <> <> << RETURNS[ filter: FilterType _ NEW[ FilterTypeRec[ filterRecLength+1]]] ~ {>> <> <> <> <> <> <> <> <> <> <> <> <> <> <<};>> <<>> 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 => {>> <> <> 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.