diff --git a/patch-image.cabal b/patch-image.cabal --- a/patch-image.cabal +++ b/patch-image.cabal @@ -1,10 +1,10 @@ Name: patch-image -Version: 0.1.0.1 +Version: 0.1.0.2 License: BSD3 License-File: LICENSE Author: Henning Thielemann Maintainer: Henning Thielemann -Homepage: http://code.haskell.org/~thielema/patch-image/ +Homepage: http://hub.darcs.net/thielema/patch-image/ Category: Graphics Synopsis: Compose a big image from overlapping parts Description: @@ -185,13 +185,13 @@ Build-Type: Simple Source-Repository this - Tag: 0.1.0.1 + Tag: 0.1.0.2 Type: darcs - Location: http://code.haskell.org/~thielema/patch-image/ + Location: http://hub.darcs.net/thielema/patch-image/ Source-Repository head Type: darcs - Location: http://code.haskell.org/~thielema/patch-image/ + Location: http://hub.darcs.net/thielema/patch-image/ Flag buildDraft description: Build draft program @@ -208,17 +208,18 @@ GHC-Prof-Options: -fprof-auto -rtsopts Build-Depends: - accelerate-arithmetic >=0.0 && <0.1, - accelerate-utility >=0.0 && <0.1, + accelerate-fourier >=0.0 && <0.1, + accelerate-arithmetic >=0.1 && <0.2, + accelerate-utility >=0.1 && <0.2, + accelerate-cufft >=0.0 && <0.1, accelerate-cuda >=0.15 && <0.16, - accelerate-fft >=0.15 && <0.16, accelerate-io >=0.15 && <0.16, accelerate >=0.15 && <0.16, - JuicyPixels >=2.0 && <3.2, + JuicyPixels >=2.0 && <3.3, hmatrix >=0.15 && <0.16, gnuplot >=0.5 && <0.6, - vector >=0.10 && <0.11, - Cabal >=1.18 && <1.22, + vector >=0.10 && <0.13, + Cabal >=1.18 && <2, filepath >=1.3 && <1.4, utility-ht >=0.0.1 && <0.1, base >=4 && <5 diff --git a/src/Accelerate.hs b/src/Accelerate.hs --- a/src/Accelerate.hs +++ b/src/Accelerate.hs @@ -3,11 +3,13 @@ import qualified Option -import qualified Data.Array.Accelerate.Math.FFT as FFT +import qualified Data.Array.Accelerate.Fourier.Real as FourierReal +import qualified Data.Array.Accelerate.CUFFT.Single as CUFFT import qualified Data.Array.Accelerate.Data.Complex as Complex +import qualified Data.Array.Accelerate.CUDA.Foreign as CUDAForeign import qualified Data.Array.Accelerate.CUDA as CUDA import qualified Data.Array.Accelerate.IO as AIO -import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as LinAlg +import qualified Data.Array.Accelerate.LinearAlgebra as LinAlg import qualified Data.Array.Accelerate.Utility.Lift.Run as Run import qualified Data.Array.Accelerate.Utility.Lift.Acc as Acc import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp @@ -15,11 +17,12 @@ import qualified Data.Array.Accelerate.Utility.Loop as Loop import qualified Data.Array.Accelerate as A import Data.Array.Accelerate.Data.Complex (Complex((:+)), ) -import Data.Array.Accelerate.Utility.Lift.Exp (atom) +import Data.Array.Accelerate.Utility.Lift.Exp (expr) +import Data.Array.Accelerate.Utility.Ord (argmaximum) import Data.Array.Accelerate (Acc, Array, Exp, DIM1, DIM2, DIM3, (:.)((:.)), Z(Z), Any(Any), All(All), - (<*), (<=*), (>=*), (==*), (&&*), (||*), (?), ) + (<*), (<=*), (>=*), (==*), (&&*), (||*), (?), (!), ) import qualified Data.Packed.Matrix as Matrix import qualified Data.Packed.Vector as Vector @@ -40,6 +43,7 @@ import qualified Data.Vector.Storable as SV import qualified System.FilePath as FilePath +import qualified System.IO as IO import qualified Distribution.Simple.Utils as CmdLine import Distribution.Verbosity (Verbosity) @@ -59,7 +63,9 @@ import Data.Tuple.HT (mapPair, mapFst, mapSnd, fst3, thd3) import Data.Word (Word8) +import System.IO.Unsafe (unsafePerformIO) + readImage :: Verbosity -> FilePath -> IO (Array DIM3 Word8) readImage verbosity path = do epic <- Pic.readImage path @@ -113,13 +119,13 @@ cycleLeftDim3 :: Exp DIM3 -> Exp DIM3 cycleLeftDim3 = - Exp.modify (atom :. atom :. atom :. atom) $ + Exp.modify (expr :. expr :. expr :. expr) $ \(z :. chans :. height :. width) -> z :. height :. width :. chans cycleRightDim3 :: Exp DIM3 -> Exp DIM3 cycleRightDim3 = - Exp.modify (atom :. atom :. atom :. atom) $ + Exp.modify (expr :. expr :. expr :. expr) $ \(z :. height :. width :. chans) -> z :. chans :. height :. width @@ -211,7 +217,7 @@ let (_ :. height :. width) = unliftDim2 $ A.shape arr xc = max 0 $ min (width -1) x yc = max 0 $ min (height-1) y - in arr A.! A.lift (ix :. yc :. xc) + in arr ! A.lift (ix :. yc :. xc) indexFrac :: (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) => @@ -383,7 +389,7 @@ differentiate arr = let size = A.unindex1 $ A.shape arr in A.generate (A.index1 (size-1)) $ \i -> - arr A.! (A.index1 $ A.unindex1 i + 1) - arr A.! i + arr ! (A.index1 $ A.unindex1 i + 1) - arr ! i scoreRotation :: Float -> Array DIM3 Word8 -> Float scoreRotation = @@ -398,9 +404,99 @@ Key.maximum (flip scoreRotation pic . (* (pi/180))) angles +magnitudeSqr :: (A.Elt a, A.IsNum a) => Exp (Complex a) -> Exp a +magnitudeSqr = + Exp.modify (expr:+expr) $ \(r:+i) -> r*r+i*i +fourierTransformationRun :: Array DIM3 Word8 -> IO (Array DIM2 Word8) +fourierTransformationRun pic = do + let (shape@(Z:.height:.width):._) = A.arrayShape pic + plan <- + CUDAForeign.inDefaultContext $ + CUFFT.plan2D CUFFT.forwardReal shape + let trans = + Run.with CUDA.run1 $ \arr -> + imageByteFromFloat $ + A.map (1e-9*) $ + A.zipWith (*) + (A.map + (Exp.modify (expr,expr) $ \(k,j) -> + magnitudeSqr $ A.lift $ + A.fromIntegral k :+ A.fromIntegral j) $ + displacementMap + (A.constant (div width 2)) + (A.constant (div height 2)) + (A.constant shape)) $ + A.map magnitudeSqr $ + CUFFT.transform plan $ + brightnessPlane $ separateChannels $ + imageFloatFromByte arr + return $ trans pic +fourierTransformation :: Option.Option -> FilePath -> Array DIM3 Word8 -> IO () +fourierTransformation opt path pic = do + let stem = FilePath.takeBaseName path + spec <- fourierTransformationRun pic + writeGrey (Option.quality opt) + (printf "/tmp/%s-spectrum.jpeg" stem) spec +{- +Spectrum coefficients with high indices +are disturbed by periodization artifacts. +That is, we might fade out the picture at the borders +or soften the disturbed spectrum coefficients. +However, I do not know for a universally good solution, +thus I leave it as it is. +-} +scoreSlopes :: + (A.Elt a, A.IsFloating a) => + (Exp Int, Exp Int) -> + Acc (Channel Z (Complex a)) -> Acc (Array DIM1 a) +scoreSlopes (minX, maxX) arr = + let shape = A.shape arr + (_z:.height:.width) = Exp.unlift (expr:.expr:.expr) shape + width2 = div width 2 + height2 = div height 2 + weighted = + A.zipWith (*) (A.map magnitudeSqr arr) $ + A.map + (Exp.modify (expr,expr) $ \(k,j) -> + magnitudeSqr $ A.lift $ + A.fromIntegral k :+ A.fromIntegral j) $ + displacementMap width2 height2 shape + in A.fold (+) 0 $ + A.generate (A.lift (Z:.maxX-minX+1:.height2)) $ + Exp.modify (expr:.expr:.expr) $ \(_z:.x:.y) -> + let (xi,frac) = + splitFraction $ + A.fromIntegral (x + minX) * + A.fromIntegral y / A.fromIntegral height2 + z0 = weighted ! A.lift (Z :. y :. mod xi width) + z1 = weighted ! A.lift (Z :. y :. mod (xi+1) width) + in linearIp (z0,z1) frac + +radonAngle :: (Float, Float) -> Array DIM3 Word8 -> IO Float +radonAngle (minAngle,maxAngle) pic = do + let (shape@(Z :. height :. _width):._) = A.arrayShape pic + plan <- + CUDAForeign.inDefaultContext $ + CUFFT.plan2D CUFFT.forwardReal shape + let height2 = fromIntegral (div height 2) + let slope w = tan (w*pi/180) * height2 + let minX = floor $ slope minAngle + let maxX = ceiling $ slope maxAngle + let angle s = atan (s/height2) * 180/pi + let trans = + Run.with CUDA.run1 $ \arr -> + A.map A.snd $ argmaximum $ + Arrange.mapWithIndex (\ix s -> A.lift (s, A.unindex1 ix)) $ + scoreSlopes (A.constant minX, A.constant maxX) $ + CUFFT.transform plan $ + brightnessPlane $ separateChannels $ + imageFloatFromByte arr + return $ angle $ fromIntegral $ Acc.the (trans pic) + minX + + rotateManifest :: Float -> Array DIM3 Word8 -> Array DIM3 Float rotateManifest = let rot = @@ -437,23 +533,48 @@ let (y, x) = A.unlift $ A.unindex2 p in (y<*height &&* x<*width) ? - (arr A.! A.index2 y x, a) + (arr ! A.index2 y x, a) -convolveImpossible :: +mulConj :: (A.Elt a, A.IsFloating a) => + Exp (Complex a) -> Exp (Complex a) -> Exp (Complex a) +mulConj x y = x * Complex.conjugate y + + +fft2DGen :: + (A.Elt e, CUFFT.Real e) => + CUFFT.Mode DIM2 e a b -> DIM2 -> CUFFT.Transform DIM2 a b +fft2DGen mode sh = + CUFFT.transform $ unsafePerformIO $ + CUDAForeign.inDefaultContext $ CUFFT.plan2D mode sh + +fft2DPlain :: + (A.Elt e, CUFFT.Real e, A.Elt a, A.Elt b) => + CUFFT.Mode DIM2 e a b -> + Channel Z a -> Acc (Channel Z b) +fft2DPlain mode arr = + A.use $ CUDA.run1 (fft2DGen mode $ A.arrayShape arr) arr + +fft2D :: + (A.Elt e, CUFFT.Real e, A.Elt a, A.Elt b) => + CUFFT.Mode DIM2 e a b -> + Int -> Int -> + Acc (Channel Z a) -> Acc (Channel Z b) +fft2D mode width height = fft2DGen mode (Z:.height:.width) + + +correlateImpossible :: + (A.Elt a, CUFFT.Real a) => Acc (Channel Z a) -> Acc (Channel Z a) -> Acc (Channel Z a) -convolveImpossible x y = +correlateImpossible x y = let (heightx, widthx) = A.unlift $ A.unindex2 $ A.shape x (heighty, widthy) = A.unlift $ A.unindex2 $ A.shape y width = ceilingPow2Exp $ widthx + widthy height = ceilingPow2Exp $ heightx + heighty sh = A.index2 height width - forward z = - FFT.fft2D FFT.Forward $ CUDA.run $ - A.map (A.lift . (:+ 0)) $ pad 0 sh z - in A.map Complex.real $ - FFT.fft2D FFT.Inverse $ CUDA.run $ - A.zipWith (*) (forward x) (forward y) + forward z = fft2DPlain CUFFT.forwardReal $ CUDA.run $ pad 0 sh z + in fft2DPlain CUFFT.inverseReal $ CUDA.run $ + A.zipWith mulConj (forward x) (forward y) ceilingPow2 :: Int -> Int @@ -480,7 +601,7 @@ clearDCCoefficient arr = A.generate (A.shape arr) $ \p -> let (_z:.y:.x) = unliftDim2 p - in x==*0 ||* y==*0 ? (0, arr A.! p) + in x==*0 ||* y==*0 ? (0, arr!p) smooth3 :: (A.Elt a, A.IsFloating a) => A.Stencil3 a -> Exp a @@ -497,104 +618,70 @@ A.zipWith (-) arr $ lowpass count arr -convolvePaddedSimple :: - (A.Elt a, A.IsFloating a) => +correlatePaddedSimple :: + (A.Elt a, CUFFT.Real a) => DIM2 -> Acc (Channel Z a) -> Acc (Channel Z a) -> Acc (Channel Z a) -convolvePaddedSimple sh@(Z :. height :. width) = - let forward = - FFT.fft2D' FFT.Forward width height . - A.map (A.lift . (:+ 0)) . pad 0 (A.lift sh) - inverse = FFT.fft2D' FFT.Inverse width height - in \ x y -> - A.map Complex.real $ inverse $ - A.zipWith (\xi yi -> xi * Complex.conjugate yi) (forward x) (forward y) +correlatePaddedSimple sh@(Z :. height :. width) = + let forward = fft2D CUFFT.forwardReal width height . pad 0 (A.lift sh) + inverse = fft2D CUFFT.inverseReal width height + in \ x y -> inverse $ A.zipWith mulConj (forward x) (forward y) -imagUnit :: (A.Elt a, A.IsNum a) => Exp (Complex a) -imagUnit = Exp.modify2 atom atom (:+) 0 1 - -{- | -Let f and g be two real valued images. -The spectrum of f+i*g is spec f + i * spec g. -Let 'flip' be the spectrum with negated indices modulo image size. -It holds: flip (spec f) = conj (spec f). - -(a + conj b) / 2 - = (spec (f+i*g) + conj (flip (spec (f+i*g)))) / 2 - = (spec f + i*spec g + conj (flip (spec f)) + conj (flip (spec (i*g)))) / 2 - = (2*spec f + i*spec g + conj (i*flip (spec g))) / 2 - = (2*spec f + i*spec g - i * conj (flip (spec g))) / 2 - = spec f - -(a - conj b) * (-i/2) - = (-i*a + conj (-i*b)) / 2 - -> this swaps role of f and g in the proof above --} -untangleRealSpectra :: - (A.Elt a, A.IsFloating a) => - Acc (Array DIM2 (Complex a)) -> Acc (Array DIM2 (Complex a, Complex a)) -untangleRealSpectra spec = - A.zipWith - (\a b -> - A.lift $ - ((a + Complex.conjugate b) / 2, - (a - Complex.conjugate b) * (-imagUnit / 2))) - spec $ - A.backpermute (A.shape spec) - (Exp.modify (atom:.atom:.atom) $ - \(_z:.y:.x) -> - let (_z:.height:.width) = unliftDim2 $ A.shape spec - in Z :. mod (-y) height :. mod (-x) width) - spec - {- -This is more efficient than 'convolvePaddedSimple' -since it needs only one forward Fourier transform, -where 'convolvePaddedSimple' needs two of them. +This is more efficient than 'correlatePaddedSimple' +since it needs only one complex forward Fourier transform, +where 'correlatePaddedSimple' needs two real transforms. +Especially for odd sizes +two real transforms are slower than a complex transform. For the analysis part, perform two real-valued Fourier transforms using one complex-valued transform. Afterwards we untangle the superposed spectra. -} -convolvePadded :: - (A.Elt a, A.IsFloating a) => +correlatePadded :: + (A.Elt a, CUFFT.Real a) => DIM2 -> Acc (Channel Z a) -> Acc (Channel Z a) -> Acc (Channel Z a) -convolvePadded sh@(Z :. height :. width) = - let forward = FFT.fft2D' FFT.Forward width height - inverse = FFT.fft2D' FFT.Inverse width height +correlatePadded sh@(Z :. height :. width) = + let forward = fft2D CUFFT.forwardComplex width height + inverse = fft2D CUFFT.inverseReal width height in \ a b -> - A.map Complex.real $ inverse $ - A.map (Exp.modify (atom,atom) $ \(ai,bi) -> ai * Complex.conjugate bi) $ - untangleRealSpectra $ forward $ + inverse $ A.map (A.uncurry mulConj) $ + FourierReal.untangleSpectra2d $ forward $ pad 0 (A.lift sh) $ - A.zipWith (Exp.modify2 atom atom (:+)) a b + A.zipWith (Exp.modify2 expr expr (:+)) a b +wrap :: Exp Int -> Exp Int -> Exp Int -> Exp Int +wrap size split c = c<*split ? (c, c-size) + +displacementMap :: + Exp Int -> Exp Int -> Exp DIM2 -> Acc (Channel Z (Int, Int)) +displacementMap xsplit ysplit sh = + let (_z :. height :. width) = unliftDim2 sh + in A.generate sh $ \p -> + let (_z:.y:.x) = unliftDim2 p + in A.lift (wrap width xsplit x, wrap height ysplit y) + attachDisplacements :: (A.Elt a, A.IsScalar a) => Exp Int -> Exp Int -> - Acc (Channel Z a) -> Acc (Channel Z ((Int, Int), a)) + Acc (Channel Z a) -> Acc (Channel Z (a, (Int, Int))) attachDisplacements xsplit ysplit arr = - let sh = A.shape arr - (_z :. height :. width) = unliftDim2 sh - in A.generate sh $ \p -> - let (_z:.y:.x) = unliftDim2 p - wrap size split c = c<*split ? (c, c-size) - in A.lift ((wrap width xsplit x, wrap height ysplit y), arr A.! p) + A.zip arr $ displacementMap xsplit ysplit (A.shape arr) weightOverlapScores :: (A.Elt a, A.IsFloating a, A.IsScalar a) => Exp Int -> (Exp Int, Exp Int) -> (Exp Int, Exp Int) -> - Acc (Channel Z ((Int, Int), a)) -> - Acc (Channel Z ((Int, Int), a)) + Acc (Channel Z (a, (Int, Int))) -> + Acc (Channel Z (a, (Int, Int))) weightOverlapScores minOverlap (widtha,heighta) (widthb,heightb) = A.map - (Exp.modify ((atom,atom),atom) $ \(dp@(dy,dx),v) -> + (Exp.modify (expr,(expr,expr)) $ \(v, dp@(dy,dx)) -> let clipWidth = min widtha (widthb + dx) - max 0 dx clipHeight = min heighta (heightb + dy) - max 0 dy - in (dp, - (clipWidth >=* minOverlap &&* clipHeight >=* minOverlap) - ? - (v / (A.fromIntegral clipWidth * A.fromIntegral clipHeight), 0))) + in ((clipWidth >=* minOverlap &&* clipHeight >=* minOverlap) + ? + (v / (A.fromIntegral clipWidth * A.fromIntegral clipHeight), 0), + dp)) {- | Set all scores to zero within a certain border. @@ -604,36 +691,26 @@ minimumOverlapScores :: (A.Elt a, A.IsFloating a, A.IsScalar a) => Exp Int -> (Exp Int, Exp Int) -> (Exp Int, Exp Int) -> - Acc (Channel Z ((Int, Int), a)) -> - Acc (Channel Z ((Int, Int), a)) + Acc (Channel Z (a, (Int, Int))) -> + Acc (Channel Z (a, (Int, Int))) minimumOverlapScores minOverlap (widtha,heighta) (widthb,heightb) = A.map - (Exp.modify ((atom,atom),atom) $ \(dp@(dy,dx),v) -> + (Exp.modify (expr,(expr,expr)) $ \(v, dp@(dy,dx)) -> let clipWidth = min widtha (widthb + dx) - max 0 dx clipHeight = min heighta (heightb + dy) - max 0 dy - in (dp, - (clipWidth >=* minOverlap &&* clipHeight >=* minOverlap) - ? - (v, 0))) - -argmax :: - (A.Elt a, A.Elt b, A.IsScalar b) => - Exp (a, b) -> Exp (a, b) -> Exp (a, b) -argmax x y = A.snd x <* A.snd y ? (y,x) - -argmaximum :: - (A.Elt a, A.Elt b, A.IsScalar b) => - Acc (Channel Z (a, b)) -> Acc (A.Scalar (a, b)) -argmaximum = A.fold1All argmax + in ((clipWidth >=* minOverlap &&* clipHeight >=* minOverlap) + ? + (v, 0), + dp)) allOverlaps :: DIM2 -> Exp Float -> Acc (Channel Z Float) -> Acc (Channel Z Float) -> - Acc (Channel Z ((Int, Int), Float)) + Acc (Channel Z (Float, (Int, Int))) allOverlaps size@(Z :. height :. width) minOverlapPortion = - let convolve = convolvePadded size + let correlate = correlatePadded size in \a b -> let (Z :. heighta :. widtha) = A.unlift $ A.shape a (Z :. heightb :. widthb) = A.unlift $ A.shape b @@ -660,7 +737,7 @@ attachDisplacements (half $ A.lift width - widthb + widtha) (half $ A.lift height - heightb + heighta) $ - convolve a b + correlate a b allOverlapsRun :: @@ -670,10 +747,10 @@ imageByteFromFloat $ -- A.map (2*) $ A.map (0.0001*) $ - A.map A.snd $ allOverlaps padExtent minOverlap picA picB + A.map A.fst $ allOverlaps padExtent minOverlap picA picB optimalOverlap :: - DIM2 -> Float -> Channel Z Float -> Channel Z Float -> ((Int, Int), Float) + DIM2 -> Float -> Channel Z Float -> Channel Z Float -> (Float, (Int, Int)) optimalOverlap padExtent = let run = Run.with CUDA.run1 $ \minimumOverlap a b -> @@ -690,7 +767,7 @@ A.fold1 (+) $ A.fold1 (+) $ A.backpermute (A.lift $ shape :. div height yk :. div width xk :. yk :. xk) - (Exp.modify (atom:.atom:.atom:.atom:.atom) $ + (Exp.modify (expr:.expr:.expr:.expr:.expr) $ \(z:.yi:.xi:.yj:.xj) -> z:.yi*yk+yj:.xi*xk+xj) arr @@ -712,7 +789,7 @@ Reduce image sizes below the padExtent before matching images. -} optimalOverlapBig :: - DIM2 -> Float -> Channel Z Float -> Channel Z Float -> ((Int, Int), Float) + DIM2 -> Float -> Channel Z Float -> Channel Z Float -> (Float, (Int, Int)) optimalOverlapBig padExtent = let run = Run.with CUDA.run1 $ \minimumOverlap a b -> @@ -720,8 +797,8 @@ shrinkFactors padExtent (A.unlift $ A.shape a) (A.unlift $ A.shape b) scalePos = - Exp.modify ((atom,atom), atom) $ - \((xm,ym), score) -> ((xm*xk, ym*yk), score) + Exp.modify (expr, (expr,expr)) $ + \(score, (xm,ym)) -> (score, (xm*xk, ym*yk)) in A.map scalePos $ argmaximum $ allOverlaps padExtent minimumOverlap (shrink factors a) (shrink factors b) @@ -736,7 +813,7 @@ clip (left,top) (width,height) arr = A.backpermute (A.lift $ A.indexTail (A.indexTail (A.shape arr)) :. height :. width) - (Exp.modify (atom:.atom:.atom) $ + (Exp.modify (expr:.expr:.expr) $ \(z :. y :. x) -> z :. y+top :. x+left) arr @@ -762,17 +839,18 @@ using a part in the overlapping area. -} optimalOverlapBigFine :: - DIM2 -> Float -> Channel Z Float -> Channel Z Float -> ((Int, Int), Float) + DIM2 -> Float -> Channel Z Float -> Channel Z Float -> (Float, (Int, Int)) optimalOverlapBigFine padExtent@(Z:.heightPad:.widthPad) = - let run = + let overlaps = allOverlaps padExtent + run = Run.with CUDA.run1 $ \minimumOverlap a b -> let shapeA = A.unlift $ A.shape a shapeB = A.unlift $ A.shape b factors@(_z:.yk:.xk) = shrinkFactors padExtent shapeA shapeB coarsed@(coarsedx,coarsedy) = mapPair ((xk*), (yk*)) $ - Exp.unliftPair $ A.fst $ A.the $ argmaximum $ - allOverlaps padExtent minimumOverlap + Exp.unliftPair $ A.snd $ A.the $ argmaximum $ + overlaps minimumOverlap (shrink factors a) (shrink factors b) ((leftOverlap, topOverlap), _, @@ -785,10 +863,10 @@ leftFocus = leftOverlap + div (widthOverlap-widthFocus) 2 topFocus = topOverlap + div (heightOverlap-heightFocus) 2 addCoarsePos = - Exp.modify ((atom,atom), atom) $ - \((xm,ym), score) -> ((xm+coarsedx, ym+coarsedy), score) + Exp.modify (expr, (expr,expr)) $ + \(score, (xm,ym)) -> (score, (xm+coarsedx, ym+coarsedy)) in A.map addCoarsePos $ argmaximum $ - allOverlaps padExtent minimumOverlap + overlaps minimumOverlap (clip (leftFocus,topFocus) extentFocus a) (clip (leftFocus-coarsedx,topFocus-coarsedy) extentFocus b) in \minimumOverlap a b -> Acc.the $ run minimumOverlap a b @@ -804,7 +882,7 @@ optimalOverlapBigMulti :: DIM2 -> DIM2 -> Int -> Float -> Float -> Channel Z Float -> Channel Z Float -> - [((Int, Int), (Int, Int), Float)] + [(Float, (Int, Int), (Int, Int))] optimalOverlapBigMulti padExtent (Z:.heightStamp:.widthStamp) numCorrs = let overlapShrunk = Run.with CUDA.run1 $ @@ -824,13 +902,13 @@ \minimumOverlap a b anchorA@(leftA, topA) anchorB@(leftB, topB) extent@(width,height) -> let addCoarsePos = - Exp.modify ((atom,atom), atom) $ - \((xm,ym), score) -> + Exp.modify (expr, (expr,expr)) $ + \(score, (xm,ym)) -> let xc = div (width+xm) 2 yc = div (height+ym) 2 - in ((leftA+xc, topA+yc), - (leftB+xc-xm, topB+yc-ym), - score) + in (score, + (leftA+xc, topA+yc), + (leftB+xc-xm, topB+yc-ym)) in A.map addCoarsePos $ argmaximum $ allOverlapsFine minimumOverlap (clip anchorA extent a) @@ -840,7 +918,7 @@ let factors@(Z:.yk:.xk) = shrinkFactors padExtent (A.arrayShape a) (A.arrayShape b) - ((shrunkdx, shrunkdy), _score) = + (_score, (shrunkdx, shrunkdy)) = Acc.the $ overlapShrunk minimumOverlap factors a b coarsedx = shrunkdx * xk @@ -1046,8 +1124,8 @@ inPicA = 0<=*xa &&* xa<*widtha &&* 0<=*ya &&* ya<*heighta inPicB = 0<=*xb &&* xb<*widthb &&* 0<=*yb &&* yb<*heightb in inPicA ? - (inPicB ? ((a A.! pa + b A.! pb)/2, a A.! pa), - inPicB ? (b A.! pb, 0)) + (inPicB ? ((a!pa + b!pb)/2, a!pa), + inPicB ? (b!pb, 0)) composeOverlap :: (Int, Int) -> @@ -1182,10 +1260,10 @@ (A.Elt a, A.IsFloating a) => Exp DIM2 -> Acc (Array DIM1 ((a,a),(a,a))) -> Acc (Channel Z a) distanceMapEdges sh edges = - A.map (Exp.modify (atom,atom) $ \(valid, dist) -> valid ? (dist, 0)) $ + A.map (Exp.modify (expr,expr) $ \(valid, dist) -> valid ? (dist, 0)) $ maskedMinimum $ outerVector - (Exp.modify2 (atom,atom) ((atom, atom), (atom, atom)) $ \p (q0, q1) -> + (Exp.modify2 (expr,expr) ((expr, expr), (expr, expr)) $ \p (q0, q1) -> mapSnd (distance p) $ project p (q0, q1)) (pixelCoordinates sh) edges @@ -1203,7 +1281,7 @@ Acc (Channel Z (Bool, (((a,(a,a)), (a,(a,a))), ((a,(a,a)), (a,(a,a)))))) distanceMapBox sh geom = let (rot, mov, extent@(width,height)) = - Exp.unlift ((atom,atom),(atom,atom),(atom,atom)) geom + Exp.unlift ((expr,expr),(expr,expr),(expr,expr)) geom widthf = A.fromIntegral width heightf = A.fromIntegral height back = rotateStretchMoveBackPoint rot mov @@ -1239,7 +1317,7 @@ Acc (Array DIM3 (Bool, a)) separateDistanceMap arr = outerVector - (Exp.modify2 (atom, ((atom, atom), (atom, atom))) (atom,atom) $ + (Exp.modify2 (expr, ((expr, expr), (expr, expr))) (expr,expr) $ \(b,(horiz,vert)) (orient,side) -> (b, orient ? (side ? horiz, side ? vert))) arr @@ -1254,7 +1332,7 @@ (4/) $ A.fromIntegral $ uncurry min $ Exp.unliftPair $ Exp.thd3 geom in imageByteFromFloat $ - A.map (Exp.modify (atom,atom) $ + A.map (Exp.modify (expr,expr) $ \(valid, dist) -> valid ? (scale*dist, 0)) $ maskedMinimum $ A.map (Exp.mapSnd A.fst) $ @@ -1284,7 +1362,7 @@ A.fold1 (||*) $ breakFusion $ outerVector - (Exp.modify2 (atom,atom) ((atom,atom),(atom,atom),(atom,atom)) $ + (Exp.modify2 (expr,expr) ((expr,expr),(expr,expr),(expr,expr)) $ \(xdst,ydst) (rot, mov, extent) -> let (xsrc,ysrc) = rotateStretchMoveBackPoint rot mov (xdst,ydst) in inBox extent (fastRound xsrc, fastRound ysrc)) @@ -1304,11 +1382,11 @@ contained = containedAnywhere others $ A.map (A.snd . A.snd) distMap - in A.map (Exp.modify (atom,atom) $ + in A.map (Exp.modify (expr,expr) $ \(valid, dist) -> valid ? (dist, 0)) $ maskedMinimum $ A.zipWith - (Exp.modify2 atom (atom,(atom,atom)) $ \c (b,(dist,_)) -> + (Exp.modify2 expr (expr,(expr,expr)) $ \c (b,(dist,_)) -> (c&&*b, dist)) contained distMap @@ -1334,7 +1412,7 @@ (A.Elt a, A.IsFloating a) => Exp DIM2 -> Acc (Channel Z (a,a)) pixelCoordinates sh = - A.generate sh $ Exp.modify (atom:.atom:.atom) $ \(_z:.y:.x) -> + A.generate sh $ Exp.modify (expr:.expr:.expr) $ \(_z:.y:.x) -> (A.fromIntegral x, A.fromIntegral y) distanceMapPoints :: @@ -1345,7 +1423,7 @@ distanceMapPoints a b = A.fold1 min $ outerVector - (Exp.modify2 (atom,atom) (atom,atom) distance) + (Exp.modify2 (expr,expr) (expr,expr) distance) a b distanceMapPointsRun :: @@ -1357,7 +1435,7 @@ Run.with CUDA.run1 $ \sh points -> let scale = - case Exp.unlift (atom:.atom:.atom) sh of + case Exp.unlift (expr:.expr:.expr) sh of _z:.y:.x -> (4/) $ A.fromIntegral $ min x y in imageByteFromFloat $ A.map (scale*) $ distanceMapPoints (pixelCoordinates sh) points @@ -1403,7 +1481,7 @@ Run.with CUDA.run1 $ \sh this others points -> let scale = - case Exp.unlift (atom:.atom:.atom) sh of + case Exp.unlift (expr:.expr:.expr) sh of _z:.y:.x -> (4/) $ A.fromIntegral $ min x y in imageByteFromFloat $ A.map (scale*) $ distanceMap sh this others points @@ -1461,7 +1539,7 @@ Run.with CUDA.run1 $ \this others points pic (weightSum,canvas) -> let (rot, mov, _) = - Exp.unlift ((atom,atom), (atom,atom), atom) this + Exp.unlift ((expr,expr), (expr,expr), expr) this in addToWeightedCanvas (distanceMap (A.shape weightSum) this others points, snd $ rotateStretchMove rot mov (unliftDim2 $ A.shape canvas) $ @@ -1487,7 +1565,7 @@ Run.with CUDA.run1 $ \this pic dist (weightSum,canvas) -> let (rot, mov, _) = - Exp.unlift ((atom,atom), (atom,atom), atom) this + Exp.unlift ((expr,expr), (expr,expr), expr) this in addToWeightedCanvas (dist, snd $ rotateStretchMove rot mov (unliftDim2 $ A.shape canvas) $ @@ -1575,7 +1653,7 @@ (FilePath.takeBaseName pathA) (FilePath.takeBaseName pathB)) $ allOverlapsShared picA picB - let doffset@(dox,doy) = fst $ optimalOverlapShared picA picB + let doffset@(dox,doy) = snd $ optimalOverlapShared picA picB let diff = overlapDifferenceRun doffset picA picB let overlapping = diff < Option.maximumDifference opt let d = (fromIntegral dox + fst leftTopA - fst leftTopB, @@ -1656,16 +1734,16 @@ let add (x0,y0) (x1,y1) = (fromIntegral x0 + x1, fromIntegral y0 + y1) let correspondences = map - (\(pa,pb,score) -> - (((ia, add pa leftTopA), (ib, add pb leftTopB)), score)) $ + (\(score,pa,pb) -> + (score, ((ia, add pa leftTopA), (ib, add pb leftTopB)))) $ optimalOverlapShared picA picB info $ printf "left-top: %s, %s" (show leftTopA) (show leftTopB) info $ printf "%s - %s" pathA pathB - forM_ correspondences $ \(((_ia,pa@(xa,ya)),(_ib,pb@(xb,yb))), score) -> + forM_ correspondences $ \(score, ((_ia,pa@(xa,ya)),(_ib,pb@(xb,yb)))) -> info $ printf "%s ~ %s, (%f,%f), %f" (show pa) (show pb) (xb-xa) (yb-ya) score - return $ map fst correspondences + return $ map snd correspondences let (posRots, dps) = layoutFromPairDisplacements (length picAngles) displacements @@ -1698,6 +1776,9 @@ process :: Option.Args -> IO () process args = do + IO.hSetBuffering IO.stdout IO.LineBuffering + IO.hSetBuffering IO.stderr IO.LineBuffering + let paths = Option.inputs args let opt = Option.option args let notice = CmdLine.notice (Option.verbosity opt) @@ -1712,9 +1793,14 @@ linearScale (Option.numberAngleSteps opt) (-maxAngle, maxAngle) when False $ analyseRotations angles pic - let angle = - maybe (findOptimalRotation angles pic) id $ - Option.angle imageOption + when False $ fourierTransformation opt path pic + angle <- + case Option.angle imageOption of + Just angle -> return angle + Nothing -> + if Option.radonTransform opt + then radonAngle (-maxAngle, maxAngle) pic + else return $ findOptimalRotation angles pic info $ printf "%s %f\176\n" path angle return (path, (angle*pi/180, pic)) @@ -1737,15 +1823,13 @@ pad 0 (A.lift size)) $ pic0 writeGrey (Option.quality opt) "/tmp/spectrum.jpeg" $ - CUDA.run $ imageByteFromFloat $ A.map Complex.real $ - FFT.fft2D FFT.Forward $ - CUDA.run1 - (A.map (A.lift . (:+ 0)) . - pad 0 (A.lift size)) $ + CUDA.run $ imageByteFromFloat $ A.map Complex.magnitude $ + fft2DPlain CUFFT.forwardReal $ + CUDA.run1 (pad 0 (A.lift size)) $ pic0 writeGrey (Option.quality opt) "/tmp/convolution.jpeg" $ CUDA.run $ imageByteFromFloat $ A.map (0.000001*) $ - convolvePadded size (A.use pic0) (A.use pic1) + correlatePadded size (A.use pic0) (A.use pic1) (floatPoss, picRots) <- (if Option.finetuneRotate opt diff --git a/src/Option.hs b/src/Option.hs --- a/src/Option.hs +++ b/src/Option.hs @@ -37,6 +37,7 @@ quality :: Int, maximumAbsoluteAngle :: Float, numberAngleSteps :: Int, + radonTransform :: Bool, smooth :: Int, padSize :: Int, minimumOverlap :: Float, @@ -58,6 +59,7 @@ quality = 99, maximumAbsoluteAngle = 1, numberAngleSteps = 40, + radonTransform = False, smooth = 20, padSize = 1024, minimumOverlap = 1/4, @@ -146,6 +148,10 @@ parseNumber "number of angle steps" (0<=) "non-negative" str) (printf "Number of steps for test rotations, default: %d" (numberAngleSteps defltOption)) : + + Opt.Option [] ["radon"] + (NoArg $ \flags -> return $ flags{radonTransform = True}) + (printf "Use Radon transform for estimating orientation, default: disabled") : Opt.Option [] ["smooth"] (flip ReqArg "NATURAL" $ \str flags ->