> Transform: PUBLIC PROC[ oldImage, newImage: ImagerPixelMaps.PixelMap, m: Matrix3d.Matrix4by4, moveX, moveY: REAL _ 0.0, focalLength: REAL _ 1] ~ { GetPassInfo: PRIVATE PROC[] ~ { newWindow: xyRectangle _ GetSize[ newImage]; infoX, infoY, infoX90, infoY90: PassInfo; m90: Matrix3d.Matrix4by4 _ Rotate90[ m, window]; m _ FixMatrix[ m, focalLength]; -- round to precision of 3 digits past decimal point m90 _ FixMatrix[ m90, focalLength]; m90[4][1] _ m[4][1] _ moveX; -- this isn't where they belong, but it's convenient m90[4][2] _ m[4][2] _ moveY; IF NOT FinalArea[ m, window] THEN RETURN[]; infoX _ InterArea[ m, window, newWindow, xPass]; infoY _ InterArea[ m, window, newWindow, yPass]; infoX90 _ InterArea[ m90, window, newWindow, x90Pass]; infoY90 _ InterArea[ m90, window, newWindow, y90Pass]; SELECT TRUE FROM (infoX.area >= infoY.area) AND (infoX.area >= infoX90.area) AND (infoX.area >= infoY90.area) => { IF infoX.area = 0 THEN RETURN[]; EndPixels _ GetEndPixels[ infoX.window.xSize, infoX.window.ySize]; passInfo _ infoX}; (infoY.area >= infoX90.area) AND (infoY.area >= infoY90.area) => { EndPixels _ GetEndPixels[ infoY.window.ySize, infoY.window.xSize]; passInfo _ infoY}; (infoX90.area >= infoY90.area) => { EndPixels _ GetEndPixels[ infoX90.window.xSize, infoX90.window.ySize]; passInfo _ infoX90}; ENDCASE => { EndPixels _ GetEndPixels[ infoY90.window.ySize, infoY90.window.xSize]; passInfo _ infoY90}; }; <> passInfo: PassInfo; window: xyRectangle _ GetSize[ oldImage]; GetPassInfo[]; IF passInfo = NIL THEN RETURN[]; -- can't do it; intermediate area collapses passInfo.image _ ImagerPixelMaps.Create[ oldImage.refRep.lgBitsPerPixel, [ 0, 0, passInfo.window.ySize, passInfo.window.xSize]]; ImagerPixelMaps.Clear[ passInfo.image]; filterWidth _ Filters.GetFilterWidth[]; <> Image _ oldImage; NewImage _ passInfo.image; Dir _ passInfo.pass1; SecondPass _ FALSE; Pass[ passInfo.matrix]; <> Image _ passInfo.image; NewImage _ newImage; Dir _ passInfo.pass2; SecondPass _ TRUE; Pass[ passInfo.matrix, passInfo.min]; <> <> <> EndPixels _ NIL; }; Pass: PROC[ m: Matrix3d.Matrix4by4, min: INTEGER _ 0] ~ { <> max: NAT; Window _ GetSize[ Image]; NewWindow _ GetSize[ NewImage]; SELECT Dir FROM xPass, y90Pass => max _ Window.ySize; ENDCASE => max _ Window.xSize; FOR i:INTEGER IN [ min..max) DO a, b, c, d, e: REAL; [ a, b, c, d, e] _ CrunchMatrix[ m, i]; LineTransform[ i, a, b, c, d, e]; ENDLOOP; }; LineTransform: PROC[ line: INTEGER, a, b, c, d, e: REAL] ~ { <> newLow, newHigh: INTEGER; <<>> [low: newLow, high: newHigh] _ GetLineBounds[ a, b, c, d, e]; IF NOT SecondPass THEN { -- flag all the pixels that won't get covered FOR j:INTEGER IN [0..newLow) DO EndPixels[ j].unCovered[line] _ TRUE; ENDLOOP; FOR j:INTEGER IN (newHigh..EndPixels.length) DO EndPixels[ j].unCovered[line] _ TRUE; ENDLOOP; }; FOR j:INTEGER IN [ newLow..newHigh] DO -- for each new pixel PixelTransform[ line, j, a, b, c, d, e]; ENDLOOP; }; PixelTransform: PROC[ line, pixel: NAT, a, b, c, d, e: REAL] ~ { <<>> <> b1, b2: REAL; oldLow, oldHigh: INTEGER; total, totalUnWeight, totalWeight, unRatio: REAL _ 0.0; colorTest: INTEGER _ 0; newColor: [0..256) _ 0; <> upperBound: INTEGER; cut1, cut2, partCovered: BOOLEAN; SELECT Dir FROM xPass, y90Pass => upperBound _ Window.xSize - 1; ENDCASE => upperBound _ Window.ySize - 1; [b: b1, cut: cut1] _ TestBound[ ((pixel-filterWidth - e)*d - b), (a - (pixel-filterWidth - e)*c), upperBound]; [b: b2, cut: cut2] _ TestBound[ ((pixel+filterWidth - e)*d - b), (a - (pixel+filterWidth - e)*c), upperBound]; partCovered _ cut1 OR cut2; [low: oldLow, high: oldHigh] _ GetPixelBounds[ b1, b2]; FOR k:INTEGER IN [ oldLow..oldHigh] DO <> weight, unWeight: REAL _ 0.0; pixelValue: INTEGER; lPos: REAL _ (a* (k-1) + b) / (c*(k-1) + d) + e; hPos: REAL _ (a* (k+1) + b) / (c*(k+1) + d) + e; newWidth: REAL _ (hPos - lPos) * filterWidth / 2; newPos: REAL _ lPos + newWidth; IF partCovered THEN IF k = oldLow THEN unRatio _ unRatio + Filters.UncoveredRatio[ newPos, newWidth, pixel, filterWidth, FALSE] ELSE IF k = oldHigh THEN unRatio _ unRatio + Filters.UncoveredRatio[ newPos, newWidth, pixel, filterWidth, TRUE]; weight _ Filters.PixelWeight[ newPos, newWidth, pixel, filterWidth]; <> pixelValue _ GetPixel[ line, k, Image]; IF pixelValue = -1 THEN { -- pixel is outside Image unWeight _ weight; weight _ 0 } ELSE IF SecondPass THEN -- check for uncovered pixels from interImage [weight: weight, unWeight: unWeight] _ TestWeights[ line, k, weight]; total _ total + weight * pixelValue; totalWeight _ totalWeight + weight; -- COVERED weight for this pixel totalUnWeight _ totalUnWeight + unWeight; -- UNCOVERED weight for this pixel ENDLOOP; <> IF totalWeight = 0 THEN { -- totally uncovered IF NOT SecondPass THEN EndPixels[ pixel].unCovered[line] _ TRUE} ELSE { IF totalUnWeight = 0 AND unRatio = 0 THEN -- totally covered colorTest _ Real.RoundI[total / totalWeight] ELSE -- partially covered IF SecondPass THEN { background: INTEGER _ GetPixel[ line, pixel, NewImage]; colorTest _ Real.RoundI[ ((background*totalUnWeight + total) / (totalWeight + totalUnWeight)) * (1-unRatio) + background*unRatio] } ELSE { -- add to list of partially covered pixels SetUnPixel[ line, pixel, (totalUnWeight / (totalWeight + totalUnWeight))*(1-unRatio) + unRatio]; colorTest _ Real.RoundI[total / totalWeight]; }; IF colorTest > 0 THEN -- keep color within valid range IF colorTest < 256 THEN newColor _ colorTest ELSE newColor _ 255 ELSE newColor _ 0; SetPixel[ line, pixel, newColor]; }; }; <> <> TransformPixel: PROC[ p: PixelCoords, m: Matrix3d.Matrix4by4] RETURNS[ point: PixelPos, boundsFault: BOOLEAN _ FALSE] ~ { newZ: REAL _ m[3][1]*p[1] + m[3][2]*p[2] + m[3][4]; <> IF newZ < m[1][1]/10 OR newZ < m[2][2]/10 THEN RETURN[ [0, 0], TRUE]; <> point[1] _ (m[1][1]*p[1] + m[1][2]*p[2] + m[1][4]) / newZ + m[4][1]; point[2] _ (m[2][1]*p[1] + m[2][2]*p[2] + m[2][4]) / newZ + m[4][2]; }; FinalArea: PROC[ m: Matrix3d.Matrix4by4, window: xyRectangle] RETURNS[ BOOLEAN] ~ { corner1, corner2, corner3, corner4: PixelPos; area, longestSideLength: REAL; boundsFault: BOOLEAN _ FALSE; [ point: corner1, boundsFault: boundsFault] _ TransformPixel[ [0, 0], m]; IF boundsFault THEN RETURN[ FALSE]; [point: corner2, boundsFault: boundsFault] _ TransformPixel[ [0, window.ySize-1], m]; IF boundsFault THEN RETURN[ FALSE]; [point: corner3, boundsFault: boundsFault] _ TransformPixel[ [window.xSize-1, 0], m]; IF boundsFault THEN RETURN[ FALSE]; [point: corner4, boundsFault: boundsFault] _ TransformPixel[ [window.xSize-1, window.ySize-1], m]; IF boundsFault THEN RETURN[ FALSE]; [area: area, longestSideLength: longestSideLength] _ GetArea[ corner1, corner2, corner3, corner4]; <> IF area < longestSideLength/512 THEN RETURN[ FALSE] ELSE RETURN[ TRUE]; }; InterArea: PROC[ m: Matrix3d.Matrix4by4, window, newWindow: xyRectangle, dir: Direction] RETURNS[ info: PassInfo] ~ TRUSTED { InterWindow: PROC[] ~ TRUSTED { minA, minB, maxA, maxB: INTEGER; [ low: minA, high: maxA] _ GetPixelBounds[ corner1[1], corner2[1]]; [ low: minB, high: maxB] _ GetPixelBounds[ corner3[1], corner4[1]]; SELECT dir FROM xPass => { size: INTEGER; info.min _ MIN[ newWindow.xSize, MAX[ 0, MIN[ minA, minB]]]; size _ MIN[ newWindow.xSize, MAX[ 0, MAX[ maxA+1, maxB+1]]]; info.window _ [ 0, 0, size, window.ySize]}; x90Pass => { size: INTEGER; info.min _ MIN[ newWindow.xSize, MAX[ 0, MIN[ minA, minB]]]; size _ MIN[ newWindow.xSize, MAX[ 0, MAX[ maxA+1, maxB+1]]]; info.window _ [ 0, 0, size, window.xSize]}; yPass => { size: INTEGER; info.min _ MIN[ newWindow.ySize, MAX[ 0, MIN[ minA, minB]]]; size _ MIN[ newWindow.ySize, MAX[ 0, MAX[ maxA+1, maxB+1]]]; info.window _ [ 0, 0, window.xSize, size]}; ENDCASE => { size: INTEGER; info.min _ MIN[ newWindow.ySize, MAX[ 0, MIN[ minA, minB]]]; size _ MIN[ newWindow.ySize, MAX[ 0, MAX[ maxA+1, maxB+1]]]; info.window _ [ 0, 0, window.ySize, size]}; }; <> <<>> corner1, corner2, corner3, corner4: PixelPos; boundsFault: BOOLEAN _ FALSE; a, b: INTEGER; <> matrix: Matrix3d.Matrix4by4 _ Matrix3d.Identity[]; SELECT dir FROM xPass, x90Pass => { matrix[1] _ m[1]; matrix[3] _ m[3]; matrix[4] _ m[4]; }; yPass, y90Pass => { matrix[1] _ [m[2][2], m[2][1], m[2][3], m[2][4]]; matrix[3] _ [m[3][2], m[3][1], m[3][3], m[3][4]]; matrix[4] _ [m[4][2], m[4][1], m[4][3], m[4][4]]; }; ENDCASE; info _ NEW[ PassInfoRec]; SELECT dir FROM xPass => { a _ window.xSize-1; b _ window.ySize-1}; x90Pass => { a _ window.ySize-1; b _ window.xSize-1}; yPass => { a _ window.ySize-1; b _ window.xSize-1}; ENDCASE => { a _ window.xSize-1; b _ window.ySize-1}; [ point: corner1, boundsFault: boundsFault] _ TransformPixel[ [0, 0], matrix]; IF boundsFault THEN RETURN; [ point: corner2, boundsFault: boundsFault] _ TransformPixel[ [0, b], matrix]; IF boundsFault THEN RETURN; [ point: corner3, boundsFault: boundsFault] _ TransformPixel[ [a, 0], matrix]; IF boundsFault THEN RETURN; [ point: corner4, boundsFault: boundsFault] _ TransformPixel[ [a, b], matrix]; IF boundsFault THEN RETURN; [area: info.area] _ GetArea[ corner1, corner2, corner3, corner4]; InterWindow[]; info.matrix _ m; info.pass1 _ dir; SELECT dir FROM xPass, x90Pass => info.pass2 _ yPass; ENDCASE => info.pass2 _ xPass; }; GetArea: PROC[ corner1, corner2, corner3, corner4: PixelPos] RETURNS[ area, longestSideLength: REAL] ~ { side1, side2, side3, side4: PixelPos; -- these are actually vector coordinates CrossProduct: PROC[ a, b: PixelPos] RETURNS[ aXb: REAL] ~ { RETURN[ a[1]*b[2] - a[2]*b[1]]; }; <> side1[1] _ corner2[1] - corner1[1]; side1[2] _ corner2[2] - corner1[2]; side2[1] _ corner3[1] - corner1[1]; side2[2] _ corner3[2] - corner1[2]; side3[1] _ corner4[1] - corner2[1]; side3[2] _ corner4[2] - corner2[2]; side4[1] _ corner4[1] - corner3[1]; side4[2] _ corner4[2] - corner3[2]; <<| (side1 X side2) + (side3 X side2) / 2 + (side1 X side4) / 2 | is area of quadrilateral>> area _ ABS[ CrossProduct[ side1, side2] + CrossProduct[ side3, side2]/2 + CrossProduct[ side1, side4]/2]; longestSideLength _ MAX[ MAX[ MAX[ (side1[1]*side1[1] + side1[2]*side1[2]), (side2[1]*side2[1] + side2[2]*side2[2])], (side3[1]*side3[1] + side3[2]*side3[2])], (side4[1]*side4[1] + side4[2]*side4[2])]; longestSideLength _ Real.SqRt[ longestSideLength]; }; <> Rotate90: PROC[ m: Matrix3d.Matrix4by4, window: xyRectangle] RETURNS[ m90: Matrix3d.Matrix4by4] ~ TRUSTED { m90 _ Matrix3d.Identity[]; m90 _ Matrix3d.Translate[ m90, 0, -(window.xSize-1), 0]; m90 _ Matrix3d.RotateAboutZAxis[ m90, 90]; -- rotate window 90 degrees m90 _ Matrix3d.MatMult[ m, m90]; -- transform matrix }; FixMatrix: PROC[ m: Matrix3d.Matrix4by4, focalLength: REAL] RETURNS[ fixedM: Matrix3d.Matrix4by4] ~ TRUSTED { FOR i: INTEGER IN [1..4] DO FOR j: INTEGER IN [1..4] DO fixedM[i][j] _ Real.RoundLI[ m[i][j]*1000]/1000.0; -- truncate to e-3 resolution ENDLOOP; ENDLOOP; FOR i: INTEGER IN [1..4] DO fixedM[3][i] _ -fixedM[3][i]; ENDLOOP; fixedM _ Matrix3d.Scale[ fixedM, focalLength, focalLength, 1]; -- perspective transformation fixedM _ Matrix3d.Translate[ fixedM, 0, 0, focalLength]; }; CrunchMatrix: PROC[ m: Matrix3d.Matrix4by4, line: INTEGER] RETURNS[ a, b, c, d, e: REAL] ~ { <> SELECT Dir FROM xPass, x90Pass => { IF SecondPass THEN { linePos: REAL _ line - m[4][2]; den: REAL _ (linePos * m[3][2] - m[2][2]); IF den = 0 THEN a _ b _ c _ d _ 0 -- result is infinitely large ELSE { a _ m[1][1] + m[1][2] * (m[2][1] - linePos * m[3][1]) / den; b _ m[1][4] + m[1][2] * (m[2][4] - linePos * m[3][4]) / den; c _ m[3][1] + m[3][2] * (m[2][1] - linePos * m[3][1]) / den; d _ m[3][4] + m[3][2] * (m[2][4] - linePos * m[3][4]) / den} } ELSE { a _ m[1][1]; b _ m[1][2] * line + m[1][4]; c _ m[3][1]; d _ m[3][2] * line + m[3][4]}; e _ m[4][1]}; yPass, y90Pass => { IF SecondPass THEN { linePos: REAL _ line - m[4][1]; den: REAL _ (linePos * m[3][1] - m[1][1]); IF den = 0 THEN a _ b _ c _ d _ e _ 0 -- result is infinitely large ELSE { a _ m[2][2] + m[2][1] * (m[1][2] - linePos * m[3][2]) / den; b _ m[2][4] + m[2][1] * (m[1][4] - linePos * m[3][4]) / den; c _ m[3][2] + m[3][1] * (m[1][2] - linePos * m[3][2]) / den; d _ m[3][4] + m[3][1] * (m[1][4] - linePos * m[3][4]) / den} } ELSE { a _ m[2][2]; b _ m[2][1] * line + m[2][4]; c _ m[3][2]; d _ m[3][1] * line + m[3][4]}; e _ m[4][2]}; ENDCASE => a _ b _ c _ d _ e _ 0 }; <> GetEndPixels: PROC[ seqLen, lineLen: NAT] RETURNS[ seq: EndPixelSeq] ~ { seq _ NEW[ EndPixelSeqRec[ seqLen]]; FOR i: INTEGER IN [0..seqLen) DO seq[i] _ NEW[ UnPixelInfoRec]; seq[i].unCovered _ NEW[ CoveredSeqRec[ lineLen]]; ENDLOOP; }; SetUnPixel: PROC[ line, pixel: INTEGER, unRatio: REAL] ~ { thisUnPixel: UnPixel _ NEW[ UnPixelRec]; thisUnPixel.pixel _ line; thisUnPixel.unRatio _ unRatio; EndPixels[ pixel].unPixels _ CONS[ thisUnPixel, EndPixels[ pixel].unPixels]; }; TestWeights: PROC[ line, pixel: INTEGER, w: REAL] RETURNS[ weight, unWeight: REAL _ 0.0] ~ { weight _ w; IF EndPixels[ line].unCovered[ pixel] = TRUE THEN { unWeight _ w; weight _ 0} ELSE { -- check for partially covered pixels unPixels: UnPixelList _ EndPixels[ line].unPixels; unPixel: UnPixel; WHILE unPixels # NIL DO unPixel _ unPixels.first; IF unPixel.pixel = pixel THEN { -- got one unWeight _ w * unPixel.unRatio; weight _ w * (1 - unPixel.unRatio)}; unPixels _ unPixels.rest; ENDLOOP }; }; <> GetLineBounds: PROC[ a, b, c, d, e: REAL] RETURNS[ low, high: INTEGER] ~ { <> b1, b2, den1, den2: REAL; lastPixel, lastNewPixel: INTEGER; <> SELECT Dir FROM xPass => { lastPixel _ Window.xSize-1; lastNewPixel _ NewWindow.xSize-1}; x90Pass => { lastPixel _ Window.ySize-1; lastNewPixel _ NewWindow.xSize-1}; yPass => { lastPixel _ Window.ySize-1; lastNewPixel _ NewWindow.ySize-1}; ENDCASE => { lastPixel _ Window.xSize-1; lastNewPixel _ NewWindow.ySize-1}; den1 _ d; den2 _ c*lastPixel+d; IF den1 < 0.1 OR den2 < 0.1 THEN -- z is out of range RETURN[ 0, 0]; <> b1 _ b/den1 + e; b2 _ (a*lastPixel+b)/den2 + e; [low: low, high: high] _ GetPixelBounds[ b1, b2]; low _ MIN[ lastNewPixel, MAX[ 0, low]]; high _ MIN[ lastNewPixel, MAX[ 0, high]]; }; TestBound: PROC[ num, den: REAL _ 0.0, upperBound: INTEGER] RETURNS[ b: REAL, cut: BOOLEAN] ~ { <> IF upperBound = 0 THEN RETURN[ 0, TRUE]; IF ABS[ den] <= ABS[ num/upperBound] THEN { -- out of bounds IF (den >= 0 AND num > 0) OR (den < 0 AND num < 0) THEN b _ upperBound ELSE b _ 0; -- num/den < 0 or num/den = 0/0 = 0 cut _ TRUE} ELSE { b _ num / den; IF b < 0 THEN { b _ 0; cut _ TRUE} ELSE cut _ FALSE}; }; GetPixelBounds: PROC[ b1, b2: REAL] RETURNS[ low, high: INTEGER] ~ { <> rem: REAL; IF b1 < b2 THEN { [i:low] _ SplitToIandR[ b1]; [i:high, rem: rem] _ SplitToIandR[ b2 + 1]} ELSE IF b1 = b2 THEN { -- boundary falls within a single pixel [i:low] _ SplitToIandR[ b1]; [i:high, rem: rem] _ SplitToIandR[ b1 + 1]} ELSE { [i:low] _ SplitToIandR[ b2]; [i:high, rem: rem] _ SplitToIandR[ b1 + 1]}; IF rem = 0.0 THEN high _ high - 1; low _ low - (filterWidth-1); high _ high + (filterWidth-1); }; <> GetPixel: PRIVATE PROC[ line, pixel: INTEGER, image: ImagerPixelMaps.PixelMap] RETURNS[ sample: INTEGER _ -1] ~ { pt, ptL: PixelCoords; window: xyRectangle _ GetSize[ image]; SELECT Dir FROM xPass => { pt[1] _ pixel; pt[2] _ line}; x90Pass => { pt[1] _ (window.xSize - 1) - line; -- turn input -90 degrees pt[2] _ pixel}; -- (x,y) yPass => { pt[1] _ line; pt[2] _ pixel}; ENDCASE => { -- y90Pass pt[1] _ (window.xSize - 1) - pixel; pt[2] _ line}; IF 0 <= pt[1] AND pt[1] < window.xSize -- ensure that pixel is inside window AND 0 <= pt[2] AND pt[2] < window.ySize THEN { ptL _ ConvertRtoLCoords[ window, pt]; sample _ ImagerPixelMaps.GetPixel[ image, ptL[1], ptL[2]] } }; SetPixel: PRIVATE PROC[ line, newPixel: NAT, newColor: [0..256)] ~ { pt, ptL: PixelCoords; SELECT Dir FROM xPass, x90Pass => { pt[1] _ newPixel; pt[2] _ line}; -- (s, f) yPass, y90Pass => { pt[1] _ line; pt[2] _ newPixel}; ENDCASE => RETURN; -- no such thing ptL _ ConvertRtoLCoords[ NewWindow, pt]; ImagerPixelMapsExtras.SetPixel[ NewImage, ptL[1], ptL[2], newColor] }; GetSize: PROC[ image: ImagerPixelMaps.PixelMap] RETURNS[ window: xyRectangle] ~ { imageBounds: ImagerPixelMaps.DeviceRectangle _ image.BoundedWindow[]; window.xMin _ imageBounds.fMin; window.xSize _ imageBounds.fSize; window.yMin _ imageBounds.sMin; window.ySize _ imageBounds.sSize; }; ConvertRtoLCoords: PROC[ window: xyRectangle, ptR: PixelCoords] RETURNS[ ptL: PixelCoords] ~ { <> ptL[2] _ window.xMin + ptR[1]; ptL[1] _ (window.yMin + window.ySize - 1) - ptR[2]; }; SplitToIandR: PUBLIC PROC[ r: REAL] RETURNS[ i: INTEGER, rem: REAL] ~ { IF r >= 0.0 THEN i _ Real.FixI[ r] ELSE { i _ Real.FixI[ r]; -- truncates towards zero IF i # r THEN i _ i - 1}; -- I want to truncate towards negative infinity rem _ r - i; }; END. <<>>