packages feed

patch-image 0.3.0.1 → 0.3.1

raw patch · 10 files changed

+468/−336 lines, 10 filesdep +storable-tupledep ~knead

Dependencies added: storable-tuple

Dependency ranges changed: knead

Files

Changes.md view
@@ -1,5 +1,9 @@ # Change log for the `patch-image` package +## 0.3.1:++ * Speed up computation by moving more stuff to Knead/LLVM.+ ## 0.3:   * Allow to save program state, modify it manually
patch-image.cabal view
@@ -1,5 +1,5 @@ Name:           patch-image-Version:        0.3.0.1+Version:        0.3.1 License:        BSD3 License-File:   LICENSE Author:         Henning Thielemann <haskell@henning-thielemann.de>@@ -41,7 +41,7 @@   README.md  Source-Repository this-  Tag:         0.3.0.1+  Tag:         0.3.1   Type:        darcs   Location:    http://hub.darcs.net/thielema/patch-image/ @@ -67,7 +67,9 @@   Main-Is: Knead.hs   Other-Modules:     MatchImageBorders-    KneadShape+    Knead.CArray+    Knead.Shape+    Complex     State     LinearAlgebra     Arithmetic@@ -81,7 +83,7 @@    If flag(llvm)     Build-Depends:-      knead >=0.2.1 && <0.3,+      knead >=0.2.2 && <0.3,       llvm-extra >=0.7 && <0.8,       llvm-tf >=3.1 && <3.2,       tfp >=1.0 && <1.1,@@ -92,6 +94,7 @@       enumset >=0.0.4 && <0.1,       containers >=0.4.2 && <0.6,       fft >=0.1.7 && <0.2,+      storable-tuple >=0.0.3 && <0.1,       carray >=0.1.5 && <0.2,       array >=0.4 && <0.6,       cassava >=0.4.5 && <0.5,@@ -112,6 +115,7 @@     State     LinearAlgebra     Arithmetic+    Complex     Degree     Option.Utility     Option
src/Accelerate.hs view
@@ -5,6 +5,7 @@ import qualified State  import qualified Arithmetic as Arith+import qualified Complex as Komplex import qualified Degree import LinearAlgebra (    absolutePositionsFromPairDisplacements, fixAtLeastOnePosition,@@ -724,11 +725,24 @@  type GenDIM2 a = Z :. a :. a -shrinkFactors :: (Integral a) => DIM2 -> GenDIM2 a -> GenDIM2 a -> GenDIM2 a-shrinkFactors (Z:.heightPad:.widthPad)+-- cf. Arithmetic.minimumOverlapAbsFromPortion+minimumOverlapAbsFromPortion ::+   (Num a, Ord i) => (a -> i, i -> a) -> a -> (i, i) -> i+minimumOverlapAbsFromPortion+   (afloor, fromInt) minOverlapPortion (width, height) =+      afloor $ minOverlapPortion * fromInt (min width height)++shrinkFactors ::+   (Num a, Integral i) =>+   (a -> i, i -> a) ->+   DIM2 -> a -> GenDIM2 i -> GenDIM2 i -> GenDIM2 i+shrinkFactors methods (Z:.heightPad:.widthPad) minOverlapPortion    (Z :. heighta :. widtha) (Z :. heightb :. widthb) =-      let yk = divUp (heighta+heightb) $ fromIntegral heightPad-          xk = divUp (widtha +widthb)  $ fromIntegral widthPad+      let minOverlap =+            minimumOverlapAbsFromPortion methods minOverlapPortion+               (min widtha widthb, min heighta heightb)+          yk = divUp (heighta+heightb-minOverlap) $ fromIntegral heightPad+          xk = divUp (widtha +widthb -minOverlap) $ fromIntegral widthPad       in  Z :. yk :. xk  {-@@ -740,7 +754,8 @@    let run =           Run.with CUDA.run1 $ \minimumOverlap a b ->              let factors@(_z:.yk:.xk) =-                    shrinkFactors padExtent+                    shrinkFactors (A.floor, A.fromIntegral) padExtent+                       minimumOverlap                        (A.unlift $ A.shape a) (A.unlift $ A.shape b)                  scalePos =                     Exp.modify (expr, (expr,expr)) $@@ -792,7 +807,9 @@           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+                 factors@(_z:.yk:.xk) =+                    shrinkFactors (A.floor, A.fromIntegral) padExtent+                       minimumOverlap shapeA shapeB                  coarsed@(coarsedx,coarsedy) =                     mapPair ((xk*), (yk*)) $                     Exp.unliftPair $ A.snd $ A.the $ argmaximum $@@ -862,7 +879,8 @@     in  \minimumOverlap mMaximumDiff a b ->           let factors@(Z:.yk:.xk) =-                 shrinkFactors padExtent (A.arrayShape a) (A.arrayShape b)+                 shrinkFactors (floor, fromIntegral) padExtent+                    minimumOverlap (A.arrayShape a) (A.arrayShape b)                (_score, shrunkd@(shrunkdx, shrunkdy)) =                  Acc.the $ overlapShrunk minimumOverlap factors a b@@ -1597,7 +1615,7 @@        map picColored pics,        map          (mapPair-            (mapPair (realToFrac, realToFrac), Arith.mapComplex realToFrac))+            (mapPair (realToFrac, realToFrac), Komplex.map realToFrac))          posRots)
src/Arithmetic.hs view
@@ -1,5 +1,6 @@ module Arithmetic where +import qualified Complex as Komplex import qualified Data.Complex as Complex import Data.Complex (Complex, ) @@ -71,7 +72,7 @@    let posRotPics =          zipWith             (\(angle,pic) (pos,rot) ->-               (pos, (pairFromComplex (Complex.cis angle * rot), pic)))+               (pos, (Komplex.toPair (Complex.cis angle * rot), pic)))             picAngles floatPosRots        bbox (rot, pic) =          case extent pic of@@ -307,16 +308,6 @@ divUp :: (Integral a) => a -> a -> a divUp a b = - div (-a) b ---pairFromComplex :: (RealFloat a) => Complex a -> (a,a)-pairFromComplex z = (Complex.realPart z, Complex.imagPart z)--mapComplex :: (a -> b) -> Complex a -> Complex b-mapComplex f (r Complex.:+ i)  =  f r Complex.:+ f i--mulConj :: (RealFloat a) => Complex a -> Complex a -> Complex a-mulConj x y = x * Complex.conjugate y   -- ToDo: move to a new utility module
+ src/Complex.hs view
@@ -0,0 +1,29 @@+module Complex where++import qualified Data.Complex as HComplex+import Data.Complex (Complex((:+)))+++toPair :: (RealFloat a) => Complex a -> (a,a)+toPair z = (HComplex.realPart z, HComplex.imagPart z)++map :: (a -> b) -> Complex a -> Complex b+map f (r :+ i)  =  f r :+ f i++conjugate :: (Num a) => Complex a -> Complex a+conjugate (r :+ i)  =  r :+ negate i++add :: (Num a) => Complex a -> Complex a -> Complex a+add (xr:+xi) (yr:+yi) = (xr+yr) :+ (xi+yi)++sub :: (Num a) => Complex a -> Complex a -> Complex a+sub (xr:+xi) (yr:+yi) = (xr-yr) :+ (xi-yi)++mul :: (Num a) => Complex a -> Complex a -> Complex a+mul (xr:+xi) (yr:+yi) = (xr*yr-xi*yi) :+ (xr*yi+xi*yr)++mulConj :: (Num a) => Complex a -> Complex a -> Complex a+mulConj x y = mul x $ conjugate y++mulConj_ :: (RealFloat a) => Complex a -> Complex a -> Complex a+mulConj_ x y = x * Complex.conjugate y
src/Knead.hs view
@@ -6,6 +6,8 @@  import qualified MatchImageBorders import qualified Arithmetic as Arith+import qualified Knead.CArray as KneadCArray+import qualified Complex as Komplex import qualified Degree import MatchImageBorders (arrayCFromKnead, arrayKneadFromC) import Arithmetic (guardedPairs, maximum0)@@ -13,7 +15,7 @@    absolutePositionsFromPairDisplacements, fixAtLeastOnePosition,    layoutFromPairDisplacements, fixAtLeastOneAnglePosition,    )-import KneadShape+import Knead.Shape          (Size, Vec2(Vec2), Dim1, Dim2, Shape2, Index2, Ix2,           verticalVal, horizontalVal) import Degree (Degree(Degree), getDegree)@@ -31,7 +33,6 @@ import Data.Array.Knead.Expression          (Exp, (==*), (/=*), (<*), (<=*), (>=*), (&&*)) -import qualified Data.Array.CArray as CArray import Data.Array.IArray (amap) import Data.Array.CArray (CArray) import Data.Array.MArray (thaw)@@ -45,7 +46,7 @@ import qualified LLVM.Core as LLVM  import qualified Data.Complex as Complex-import Data.Complex (Complex((:+)), conjugate, realPart)+import Data.Complex (Complex((:+)))  import qualified Codec.Picture as Pic @@ -61,7 +62,7 @@ import Text.Printf (printf)  import qualified Control.Monad.HT as MonadHT-import Control.Monad (liftM2, when, join, foldM, (<=<))+import Control.Monad (when, join, foldM, (<=<)) import Control.Applicative (pure, (<$>), (<*>))  import qualified Data.Foldable as Fold@@ -71,12 +72,12 @@ import Data.Monoid ((<>)) import Data.Maybe.HT (toMaybe) import Data.Maybe (mapMaybe, isJust, isNothing)+import Data.Bits (Bits) import Data.Traversable (forM) import Data.Foldable (forM_) import Data.Ord.HT (comparing) import Data.Tuple.HT-         (mapPair, mapFst, mapSnd, mapTriple, swap,-          mapThd3, fst3, thd3, uncurry3)+         (mapPair, mapFst, mapSnd, mapTriple, swap, mapThd3, fst3, uncurry3) import Data.Word (Word8, Word32)  @@ -563,125 +564,63 @@   -- counterpart to 'clip'-pad ::-   (MultiValue.C a) =>-   Exp a -> Exp Dim2 -> SymbPlane a -> SymbPlane a+pad :: (MultiValue.C a) => Exp a -> Exp Dim2 -> SymbPlane a -> SymbPlane a pad a sh img =    let Vec2 height width = Expr.decompose atomDim2 $ Symb.shape img    in  generate sh $ \p ->          let Vec2 y x = Expr.decompose atomIx2 p          in  Expr.ifThenElse (y<*height &&* x<*width) (img ! p) a -padCArray ::-   (SV.Storable a) =>-   a -> (Int,Int) -> CArray (Int,Int) a -> CArray (Int,Int) a-padCArray a (height, width) img =-   CArray.listArray ((0,0), (height-1, width-1)) (repeat a)-   CArray.//-   CArray.assocs img--clipCArray ::-   (SV.Storable a) => (Int,Int) -> CArray (Int,Int) a -> CArray (Int,Int) a-clipCArray (height, width) =-   CArray.ixmap ((0,0), (height-1, width-1)) id--mapPairInt :: (Integral i, Integral j) => (i,i) -> (j,j)-mapPairInt = mapPair (fromIntegral, fromIntegral)--correlatePaddedSimpleCArray ::-   (FFTWReal a) =>-   (Int,Int) ->-   CArray (Int,Int) a ->-   CArray (Int,Int) a ->-   CArray (Int,Int) a-correlatePaddedSimpleCArray sh =-   let forward = FFT.dftRCN [0,1] . padCArray 0 sh-       inverse = FFT.dftCRN [0,1]-   in  \ a b ->-         inverse $ CArray.liftArray2 Arith.mulConj (forward a) (forward b)---- expects zero-based arrays-cyclicReverse2d :: (SV.Storable a) => CArray (Int,Int) a -> CArray (Int,Int) a+cyclicReverse2d :: (MultiValue.C a) => SymbPlane a -> SymbPlane a cyclicReverse2d spec =-   let (height, width) = mapPair ((1+), (1+)) $ snd $ CArray.bounds spec-   in  CArray.ixmap (CArray.bounds spec)-         (\(y,x) -> (mod (-y) height, mod (-x) width)) spec+   let (Vec2 height width) = Expr.decompose atomDim2 $ Symb.shape spec+   in  Symb.backpermute (Symb.shape spec)+         (Expr.modify atomIx2 $ \(Vec2 y x) ->+            Vec2+               (wrap height height (height-y))+               (wrap width width (width-x)))+         spec -untangleCoefficient ::-   (RealFloat a) => Complex a -> Complex a -> (Complex a, Complex a)-untangleCoefficient a b =-   let bc = conjugate b-   in  ((a + bc) / 2, (a - bc) * (0 :+ (-1/2)))+atomComplex :: Complex (Atom a)+atomComplex = atom:+atom --- ToDo: could be moved to fft package untangleSpectra2d ::-   (RealFloat a, SV.Storable a) =>-   CArray (Int,Int) (Complex a) -> CArray (Int,Int) (Complex a, Complex a)+   (MultiValue.C a, MultiValue.Field a,+    MultiValue.Real a, MultiValue.RationalConstant a) =>+   SymbPlane (Complex a) -> SymbPlane (Complex a, Complex a) untangleSpectra2d spec =-   CArray.liftArray2 untangleCoefficient spec (cyclicReverse2d spec)--{- |-Equivalent to @amap (uncurry Arith.mulConj) . untangleSpectra2d@-but much faster, since it avoids the slow @instance Storable (a,b)@-based on @storable-tuple:storePair@.--}-mulConjUntangledSpectra2d ::-   (RealFloat a, SV.Storable a) =>-   CArray (Int,Int) (Complex a) -> CArray (Int,Int) (Complex a)-mulConjUntangledSpectra2d spec =-   CArray.liftArray2-      ((uncurry Arith.mulConj .) . untangleCoefficient)+   Symb.zipWith+      (Expr.modify2 atomComplex atomComplex KneadCArray.untangleCoefficient)       spec (cyclicReverse2d spec) --{--This is more efficient than 'correlatePaddedSimpleCArray'-since it needs only one complex forward Fourier transform,-where 'correlatePaddedSimpleCArray' 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.--}-correlatePaddedComplexCArray ::-   (FFTWReal a) =>-   (Int,Int) ->-   CArray (Int,Int) a ->-   CArray (Int,Int) a ->-   CArray (Int,Int) a-correlatePaddedComplexCArray sh a b =-   amap realPart $ FFT.idftN [0,1] $-   mulConjUntangledSpectra2d $ FFT.dftN [0,1] $-   CArray.liftArray2 (:+) (padCArray 0 sh a) (padCArray 0 sh b)+correlatePadded ::+   (FFTWReal a, MultiValue.Real a, MultiMem.C a,+    MultiValue.Field a, MultiValue.RationalConstant a) =>+   Dim2 -> IO (Plane a -> Plane a -> IO (Plane a))+correlatePadded padExtent@(Vec2 height width) = do+   let sh = Expr.cons padExtent+   mergePlanes <-+      RenderP.run $ \a b ->+         Symb.zipWith Expr.consComplex (pad 0 sh a) (pad 0 sh b)+   let (halfWidth,parity) = divMod width 2+   let exprFromInt = Expr.cons . fromIntegral+   mulSpecs <-+      RenderP.run $+         clip (0,0) (exprFromInt $ halfWidth+1, exprFromInt height) .+         Symb.map+            (Expr.modify (atomComplex, atomComplex) $ uncurry Komplex.mulConj) .+         untangleSpectra2d -{- |-Should be yet a little bit more efficient than 'correlatePaddedComplexCArray'-since it uses a real back transform.--}-correlatePaddedCArray ::-   (FFTWReal a) =>-   (Int,Int) ->-   CArray (Int,Int) a ->-   CArray (Int,Int) a ->-   CArray (Int,Int) a-correlatePaddedCArray sh@(height,width) a b =-   (case divMod width 2 of-      (halfWidth,0) -> FFT.dftCRN [0,1] . clipCArray (height,halfWidth+1)-      (halfWidth,_) -> FFT.dftCRON [0,1] . clipCArray (height,halfWidth+1)) $-   mulConjUntangledSpectra2d $ FFT.dftN [0,1] $-   CArray.liftArray2 (:+) (padCArray 0 sh a) (padCArray 0 sh b)+   return $ \ a b ->+      liftCArray (if parity==0 then FFT.dftCRN [0,1] else FFT.dftCRON [0,1]) =<<+      mulSpecs =<< liftCArray (FFT.dftN [0,1]) =<< mergePlanes a b  -liftCArray2 ::-   (SV.Storable a) =>-   (CArray (Int,Int) a -> CArray (Int,Int) a -> CArray (Int,Int) a) ->-   Plane a -> Plane a -> IO (Plane a)-liftCArray2 f a b =-   arrayKneadFromC <$>-   liftM2 f-      (arrayCFromKnead a)-      (arrayCFromKnead b)+liftCArray ::+   (SV.Storable a, SV.Storable b) =>+   (CArray (Int,Int) a -> CArray (Int,Int) b) ->+   Plane a -> IO (Plane b)+liftCArray f a = arrayKneadFromC . f <$> arrayCFromKnead a   fixArray :: Id (Symb.Array sh a)@@ -776,17 +715,17 @@  allOverlapsRun ::    Dim2 -> IO (Float -> Plane Float -> Plane Float -> IO (Plane Word8))-allOverlapsRun padExtent@(Vec2 height width) = do+allOverlapsRun padExtent = do    run <-       RenderP.run $ \minOverlapPortion sha shb img ->          imageByteFromFloat $          Symb.map (0.0001*) $          Symb.map Expr.fst $          allOverlapsFromCorrelation padExtent minOverlapPortion sha shb img+   correlate <- correlatePadded padExtent     return $ \overlap a b ->-      run overlap (Phys.shape a) (Phys.shape b)-         =<< liftCArray2 (correlatePaddedCArray $ mapPairInt (height, width)) a b+      run overlap (Phys.shape a) (Phys.shape b) =<< correlate a b   argmax ::@@ -802,15 +741,15 @@  optimalOverlap ::    Dim2 -> IO (Float -> Plane Float -> Plane Float -> IO (Float, (Size, Size)))-optimalOverlap padExtent@(Vec2 height width) = do+optimalOverlap padExtent = do    run <-       RenderP.run $ \minOverlapPortion (sha, shb) img ->          argmaximum $          allOverlapsFromCorrelation padExtent minOverlapPortion sha shb img+   correlate <- correlatePadded padExtent     return $ \overlap a b ->-      run overlap (Phys.shape a, Phys.shape b)-         =<< liftCArray2 (correlatePaddedCArray $ mapPairInt (height, width)) a b+      run overlap (Phys.shape a, Phys.shape b) =<< correlate a b   shrink ::@@ -826,14 +765,42 @@       (Expr.modify (atomIx2, atomIx2) $          \(Vec2 yi xi, Vec2 yj xj) -> Vec2 (yi*yk+yj) (xi*xk+xj)) -shrinkFactors :: (Integral a) => Dim2 -> Shape2 a -> Shape2 a -> Shape2 a-shrinkFactors (Vec2 heightPad widthPad)+{-+The implementation accepts overlapping of at most minOverlapPortion+of the two shrunken images.+However, in practice this optimization is rarely effective.+In most cases the shrink factors are the same+independent from whether minOverlap is zero or not.+-}+shrinkFactors ::+   (Integral a) => Dim2 -> Float -> Shape2 a -> Shape2 a -> Shape2 a+shrinkFactors (Vec2 heightPad widthPad) minOverlapPortion    (Vec2 heighta widtha) (Vec2 heightb widthb) =-      Vec2-         (Arith.divUp (heighta+heightb) $ fromIntegral heightPad)-         (Arith.divUp (widtha +widthb)  $ fromIntegral widthPad)+      let minOverlap =+            Arith.minimumOverlapAbsFromPortion minOverlapPortion+               (min widtha widthb, min heighta heightb)+      in Vec2+            (Arith.divUp (heighta+heightb-minOverlap) $ fromIntegral heightPad)+            (Arith.divUp (widtha +widthb -minOverlap) $ fromIntegral widthPad) +{-+Should compute almost the same as shrinkFactors+but is less optimized and more idiomatic.+@correlationSize@ has a final @ceilingSmooth7@.+This is not necessary here+since we expect that the user chooses an FFT friendly target size.+-}+shrinkFactorsAlt ::+   (Bits a, Integral a) => Float -> Dim2 -> Shape2 a -> Shape2 a -> Shape2 a+shrinkFactorsAlt minOverlapPortion (Vec2 heightPad widthPad) a b =+   let (widthc,heightc) =+         Arith.correlationSize minOverlapPortion $+         map (\(Vec2 height width) -> (width, height)) [a,b]+   in Vec2+         (Arith.divUp heightc $ fromIntegral heightPad)+         (Arith.divUp widthc $ fromIntegral widthPad) + optimalOverlapBig ::    Dim2 -> IO (Float -> Plane Float -> Plane Float -> IO (Float, (Size, Size))) optimalOverlapBig padExtent = do@@ -841,7 +808,7 @@    optOverlap <- optimalOverlap padExtent    return $ \minimumOverlap a b -> do       let factors@(Vec2 yk xk) =-            shrinkFactors padExtent (Phys.shape a) (Phys.shape b)+            shrinkFactors padExtent minimumOverlap (Phys.shape a) (Phys.shape b)       aSmall <- shrnk factors a       bSmall <- shrnk factors b       mapSnd (mapPair ((*xk), (*yk))) <$>@@ -883,7 +850,7 @@    Dim2 -> IO (Float -> Plane Float -> Plane Float -> IO (Float, (Size, Size))) optimalOverlapBigFine padExtent@(Vec2 heightPad widthPad) = do    overlap <- optimalOverlap padExtent-   -- optimalOverlap is rendered again here+   -- optimalOverlap is compiled again here    overlapBig <- optimalOverlapBig padExtent    clp <- RenderP.run clip    return $ \minimumOverlap a b -> do@@ -936,7 +903,7 @@     return $ \minimumOverlap mMaximumDiff a b -> do       let factors@(Vec2 yk xk) =-            shrinkFactors padExtent (Phys.shape a) (Phys.shape b)+            shrinkFactors padExtent minimumOverlap (Phys.shape a) (Phys.shape b)       aSmall <- shrnk factors a       bSmall <- shrnk factors b @@ -1306,9 +1273,7 @@     MultiValue.NativeFloating a ar) =>    Exp (Geometry b) -> SymbPlane a -> SymbPlane Word8 scaleDistanceMapGeom geom img =-   let scale =-         (4/) $ fromInt $-         Expr.modify (atom,atom,(atom,atom)) (uncurry Expr.min . thd3) geom+   let scale = (4/) $ fromInt $ Expr.uncurry Expr.min $ Expr.thd3 geom    in  imageByteFromFloat $ Symb.map (scale*) img  @@ -1689,7 +1654,7 @@        map picColored pics,        map          (mapPair-            (mapPair (realToFrac, realToFrac), Arith.mapComplex realToFrac))+            (mapPair (realToFrac, realToFrac), Komplex.map realToFrac))          posRots)  @@ -1728,20 +1693,22 @@    when False $ do       notice "write fft"       let pic0 : pic1 : _ = map snd rotated-          size = (1024,768)-      cpic0 <- arrayCFromKnead pic0-      cpic1 <- arrayCFromKnead pic1+          size = Vec2 1024 768       makeByteImage <-          RenderP.run $ \k -> imageByteFromFloat . Symb.map (k*) . fixArray+      runPad <- RenderP.run pad       writeGrey (Option.quality opt) "/tmp/padded.jpeg" =<<-         (makeByteImage 1 $ arrayKneadFromC $ padCArray 0 size cpic0)+         (makeByteImage 1 =<< runPad 0 size pic0)+      runMagnitude <-+         RenderP.run $+         Symb.map (Expr.modify atomComplex $ \(r:+i) -> Expr.sqrt$ r*r+i*i)+            . fixArray       writeGrey (Option.quality opt) "/tmp/spectrum.jpeg" =<<-         (makeByteImage 0.1 $ arrayKneadFromC $-          CArray.liftArray Complex.magnitude $-          FFT.dftRCN [0,1] $ padCArray 0 size cpic0)+         (makeByteImage 0.1 =<< runMagnitude =<<+          liftCArray (FFT.dftRCN [0,1]) =<< runPad 0 size pic0)+      correlate <- correlatePadded size       writeGrey (Option.quality opt) "/tmp/convolution.jpeg" =<<-         (makeByteImage 0.1 $ arrayKneadFromC $-          correlatePaddedCArray size cpic0 cpic1)+         (makeByteImage 0.1 =<< correlate pic0 pic1)     return $       zipWith3
+ src/Knead/CArray.hs view
@@ -0,0 +1,119 @@+module Knead.CArray where++import qualified Complex as Komplex++import qualified Math.FFT as FFT+import Math.FFT.Base (FFTWReal)++import Foreign.Storable (Storable)+import Foreign.Storable.Tuple ()++import qualified Data.Array.CArray as CArray+import Data.Array.IArray (amap)+import Data.Array.CArray (CArray)++import Data.Complex (Complex((:+)), realPart)++import Data.Tuple.HT (mapPair)+++pad ::+   (Storable a) => a -> (Int,Int) -> CArray (Int,Int) a -> CArray (Int,Int) a+pad a (height, width) img =+   CArray.listArray ((0,0), (height-1, width-1)) (repeat a)+   CArray.//+   CArray.assocs img++clip :: (Storable a) => (Int,Int) -> CArray (Int,Int) a -> CArray (Int,Int) a+clip (height, width) = CArray.ixmap ((0,0), (height-1, width-1)) id++correlatePaddedSimple ::+   (FFTWReal a) =>+   (Int,Int) ->+   CArray (Int,Int) a ->+   CArray (Int,Int) a ->+   CArray (Int,Int) a+correlatePaddedSimple sh =+   let forward = FFT.dftRCN [0,1] . pad 0 sh+       inverse = FFT.dftCRN [0,1]+   in  \ a b ->+         inverse $ CArray.liftArray2 Komplex.mulConj (forward a) (forward b)++-- expects zero-based arrays+cyclicReverse2d :: (Storable a) => CArray (Int,Int) a -> CArray (Int,Int) a+cyclicReverse2d spec =+   let (height, width) = mapPair ((1+), (1+)) $ snd $ CArray.bounds spec+   in  CArray.ixmap (CArray.bounds spec)+         (\(y,x) -> (mod (-y) height, mod (-x) width)) spec++untangleCoefficient ::+   (Fractional a) => Complex a -> Complex a -> (Complex a, Complex a)+untangleCoefficient a b =+   let bc = Komplex.conjugate b+   in  (Komplex.mul (Komplex.add a bc) ((1/2) :+ 0),+        Komplex.mul (Komplex.sub a bc) (0 :+ (-1/2)))++untangleCoefficient_ ::+   (RealFloat a) => Complex a -> Complex a -> (Complex a, Complex a)+untangleCoefficient_ a b =+   let bc = Komplex.conjugate b+   in  ((a + bc) / 2, (a - bc) * (0 :+ (-1/2)))++-- ToDo: could be moved to fft package+untangleSpectra2d ::+   (Fractional a, Storable a) =>+   CArray (Int,Int) (Complex a) -> CArray (Int,Int) (Complex a, Complex a)+untangleSpectra2d spec =+   CArray.liftArray2 untangleCoefficient spec (cyclicReverse2d spec)++{- |+Equivalent to @amap (uncurry Komplex.mulConj) . untangleSpectra2d@+but much faster, since it avoids the slow @instance Storable (a,b)@+based on @storable-tuple:storePair@.+-}+mulConjUntangledSpectra2d ::+   (Fractional a, Storable a) =>+   CArray (Int,Int) (Complex a) -> CArray (Int,Int) (Complex a)+mulConjUntangledSpectra2d spec =+   CArray.liftArray2+      ((uncurry Komplex.mulConj .) . untangleCoefficient)+      spec (cyclicReverse2d spec)+++{-+This is more efficient than 'correlatePaddedSimpleCArray'+since it needs only one complex forward Fourier transform,+where 'correlatePaddedSimpleCArray' 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.+-}+correlatePaddedComplex ::+   (FFTWReal a) =>+   (Int,Int) ->+   CArray (Int,Int) a ->+   CArray (Int,Int) a ->+   CArray (Int,Int) a+correlatePaddedComplex sh a b =+   amap realPart $ FFT.idftN [0,1] $+   mulConjUntangledSpectra2d $ FFT.dftN [0,1] $+   CArray.liftArray2 (:+) (pad 0 sh a) (pad 0 sh b)++{- |+Should be yet a little bit more efficient than 'correlatePaddedComplexCArray'+since it uses a real back transform.+-}+correlatePadded ::+   (FFTWReal a) =>+   (Int,Int) ->+   CArray (Int,Int) a ->+   CArray (Int,Int) a ->+   CArray (Int,Int) a+correlatePadded sh@(height,width) a b =+   (case divMod width 2 of+      (halfWidth,0) -> FFT.dftCRN [0,1] . clip (height,halfWidth+1)+      (halfWidth,_) -> FFT.dftCRON [0,1] . clip (height,halfWidth+1)) $+   mulConjUntangledSpectra2d $ FFT.dftN [0,1] $+   CArray.liftArray2 (:+) (pad 0 sh a) (pad 0 sh b)
+ src/Knead/Shape.hs view
@@ -0,0 +1,176 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE EmptyDataDecls #-}+module Knead.Shape where++import qualified Data.Array.Knead.Index.Nested.Shape as Shape+import qualified Data.Array.Knead.Expression as Expr++import qualified LLVM.Extra.Multi.Value.Memory as MultiMem+import qualified LLVM.Extra.Multi.Value as MultiValue+import qualified LLVM.Extra.Arithmetic as A+import LLVM.Extra.Multi.Value (atom)++import qualified LLVM.Core as LLVM++import qualified Type.Data.Num.Decimal as TypeNum++import Foreign.Storable+         (Storable, sizeOf, alignment, poke, pokeElemOff, peek, peekElemOff)+import Foreign.Ptr (Ptr, castPtr)++import qualified Control.Monad.HT as Monad+import Control.Monad (join)++import Data.Int (Int64)+++{- |+I choose a bit complicated Dim2 definition+to make it distinct from size pairs with width and height swapped.+Alternatives would be Index.Linear or intentionally complicated Shape types like:++type Dim0 = ()+type Dim1 = ((), Size)+type Dim2 = ((), Size, Size)++Problems with Index.Linear is that it is fixed to Word32 dimensions+which causes trouble with negative coordinates+that we encounter on rotations.++The custom shape type requires lots of new definitions+but it is certainly the cleanest solution.+-}+type Size = Int64+type Dim0 = ()+type Dim1 = Size+type Dim2 = Shape2 Size+type Ix2  = Index2 Size++data Vec2 tag i = Vec2 {vertical, horizontal :: i}++data ShapeTag+data IndexTag++type Shape2 = Vec2 ShapeTag+type Index2 = Vec2 IndexTag++++squareShape :: n -> Vec2 tag n+squareShape n = Vec2 n n++castToElemPtr :: Ptr (Vec2 tag a) -> Ptr a+castToElemPtr = castPtr++instance (Storable n) => Storable (Vec2 tag n) where+   -- cf. sample-frame:Frame.Stereo+   sizeOf ~(Vec2 n m) =+      sizeOf n + mod (- sizeOf n) (alignment m) + sizeOf m+   alignment ~(Vec2 n _) = alignment n+   poke p (Vec2 n m) =+      let q = castToElemPtr p+      in  poke q n >> pokeElemOff q 1 m+   peek p =+      let q = castToElemPtr p+      in  Monad.lift2 Vec2 (peek q) (peekElemOff q 1)++instance (MultiValue.C n) => MultiValue.C (Vec2 tag n) where+   type Repr f (Vec2 tag n) = Vec2 tag (MultiValue.Repr f n)+   cons (Vec2 n m) =+      MultiValue.compose $ Vec2 (MultiValue.cons n) (MultiValue.cons m)+   undef = MultiValue.compose $ squareShape MultiValue.undef+   zero = MultiValue.compose $ squareShape MultiValue.zero+   phis bb a =+      case MultiValue.decompose (squareShape atom) a of+         Vec2 a0 a1 ->+            fmap MultiValue.compose $+            Monad.lift2 Vec2 (MultiValue.phis bb a0) (MultiValue.phis bb a1)+   addPhis bb a b =+      case (MultiValue.decompose (squareShape atom) a,+            MultiValue.decompose (squareShape atom) b) of+         (Vec2 a0 a1, Vec2 b0 b1) ->+            MultiValue.addPhis bb a0 b0 >>+            MultiValue.addPhis bb a1 b1++type instance+   MultiValue.Decomposed f (Vec2 tag pat) =+      Vec2 tag (MultiValue.Decomposed f pat)+type instance+   MultiValue.PatternTuple (Vec2 tag pat) =+      Vec2 tag (MultiValue.PatternTuple pat)++instance (MultiValue.Compose n) => MultiValue.Compose (Vec2 tag n) where+   type Composed (Vec2 tag n) = Vec2 tag (MultiValue.Composed n)+   compose (Vec2 n m) =+      case (MultiValue.compose n, MultiValue.compose m) of+         (MultiValue.Cons rn, MultiValue.Cons rm) ->+            MultiValue.Cons (Vec2 rn rm)++instance (MultiValue.Decompose pn) => MultiValue.Decompose (Vec2 tag pn) where+   decompose (Vec2 pn pm) (MultiValue.Cons (Vec2 n m)) =+      Vec2+         (MultiValue.decompose pn (MultiValue.Cons n))+         (MultiValue.decompose pm (MultiValue.Cons m))++instance (MultiMem.C i) => MultiMem.C (Vec2 tag i) where+   type Struct (Vec2 tag i) =+         LLVM.Struct (MultiMem.Struct i, (MultiMem.Struct i, ()))+   decompose nm =+      Monad.lift2 zipShape+         (MultiMem.decompose =<< LLVM.extractvalue nm TypeNum.d0)+         (MultiMem.decompose =<< LLVM.extractvalue nm TypeNum.d1)+   compose nm =+      case unzipShape nm of+         Vec2 n m -> do+            sn <- MultiMem.compose n+            sm <- MultiMem.compose m+            rn <- LLVM.insertvalue (LLVM.value LLVM.undef) sn TypeNum.d0+            LLVM.insertvalue rn sm TypeNum.d1+++unzipShape :: MultiValue.T (Vec2 tag n) -> Vec2 tag (MultiValue.T n)+unzipShape = MultiValue.decompose (squareShape atom)++zipShape :: MultiValue.T n -> MultiValue.T n -> MultiValue.T (Vec2 tag n)+zipShape y x = MultiValue.compose $ Vec2 y x++instance (tag ~ ShapeTag, Shape.C i) => Shape.C (Vec2 tag i) where+   type Index (Vec2 tag i) = Index2 (Shape.Index i)+   intersectCode a b =+      case (unzipShape a, unzipShape b) of+         (Vec2 an am, Vec2 bn bm) ->+            Monad.lift2 zipShape+               (Shape.intersectCode an bn)+               (Shape.intersectCode am bm)+   sizeCode nm =+      case unzipShape nm of+         Vec2 n m ->+            join $ Monad.lift2 A.mul (Shape.sizeCode n) (Shape.sizeCode m)+   size (Vec2 n m) = Shape.size n * Shape.size m+   flattenIndexRec nm ij =+      case (unzipShape nm, unzipShape ij) of+         (Vec2 n m, Vec2 i j) -> do+            (ns, il) <- Shape.flattenIndexRec n i+            (ms, jl) <- Shape.flattenIndexRec m j+            Monad.lift2 (,)+               (A.mul ns ms)+               (A.add jl =<< A.mul ms il)+   loop code nm =+      case unzipShape nm of+         Vec2 n m ->+            Shape.loop (\i -> Shape.loop (\j -> code (zipShape i j)) m) n+++instance (Expr.Compose n) => Expr.Compose (Vec2 tag n) where+   type Composed (Vec2 tag n) = Vec2 tag (Expr.Composed n)+   compose (Vec2 n m) = Expr.lift2 zipShape (Expr.compose n) (Expr.compose m)++instance (Expr.Decompose p) => Expr.Decompose (Vec2 tag p) where+   decompose (Vec2 pn pm) vec =+      Vec2+         (Expr.decompose pn (verticalVal vec))+         (Expr.decompose pm (horizontalVal vec))++verticalVal, horizontalVal :: (Expr.Value val) => val (Vec2 tag n) -> val n+verticalVal = Expr.lift1 (MultiValue.lift1 vertical)+horizontalVal = Expr.lift1 (MultiValue.lift1 horizontal)
− src/KneadShape.hs 
@@ -1,176 +0,0 @@-{-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE EmptyDataDecls #-}-module KneadShape where--import qualified Data.Array.Knead.Index.Nested.Shape as Shape-import qualified Data.Array.Knead.Expression as Expr--import qualified LLVM.Extra.Multi.Value.Memory as MultiMem-import qualified LLVM.Extra.Multi.Value as MultiValue-import qualified LLVM.Extra.Arithmetic as A-import LLVM.Extra.Multi.Value (atom)--import qualified LLVM.Core as LLVM--import qualified Type.Data.Num.Decimal as TypeNum--import Foreign.Storable-         (Storable, sizeOf, alignment, poke, pokeElemOff, peek, peekElemOff)-import Foreign.Ptr (Ptr, castPtr)--import qualified Control.Monad.HT as Monad-import Control.Monad (join)--import Data.Int (Int64)---{- |-I choose a bit complicated Dim2 definition-to make it distinct from size pairs with width and height swapped.-Alternatives would be Index.Linear or intentionally complicated Shape types like:--type Dim0 = ()-type Dim1 = ((), Size)-type Dim2 = ((), Size, Size)--Problems with Index.Linear is that it is fixed to Word32 dimensions-which causes trouble with negative coordinates-that we encounter on rotations.--The custom shape type requires lots of new definitions-but it is certainly the cleanest solution.--}-type Size = Int64-type Dim0 = ()-type Dim1 = Size-type Dim2 = Shape2 Size-type Ix2  = Index2 Size--data Vec2 tag i = Vec2 {vertical, horizontal :: i}--data ShapeTag-data IndexTag--type Shape2 = Vec2 ShapeTag-type Index2 = Vec2 IndexTag----squareShape :: n -> Vec2 tag n-squareShape n = Vec2 n n--castToElemPtr :: Ptr (Vec2 tag a) -> Ptr a-castToElemPtr = castPtr--instance (Storable n) => Storable (Vec2 tag n) where-   -- cf. sample-frame:Frame.Stereo-   sizeOf ~(Vec2 n m) =-      sizeOf n + mod (- sizeOf n) (alignment m) + sizeOf m-   alignment ~(Vec2 n _) = alignment n-   poke p (Vec2 n m) =-      let q = castToElemPtr p-      in  poke q n >> pokeElemOff q 1 m-   peek p =-      let q = castToElemPtr p-      in  Monad.lift2 Vec2 (peek q) (peekElemOff q 1)--instance (MultiValue.C n) => MultiValue.C (Vec2 tag n) where-   type Repr f (Vec2 tag n) = Vec2 tag (MultiValue.Repr f n)-   cons (Vec2 n m) =-      MultiValue.compose $ Vec2 (MultiValue.cons n) (MultiValue.cons m)-   undef = MultiValue.compose $ squareShape MultiValue.undef-   zero = MultiValue.compose $ squareShape MultiValue.zero-   phis bb a =-      case MultiValue.decompose (squareShape atom) a of-         Vec2 a0 a1 ->-            fmap MultiValue.compose $-            Monad.lift2 Vec2 (MultiValue.phis bb a0) (MultiValue.phis bb a1)-   addPhis bb a b =-      case (MultiValue.decompose (squareShape atom) a,-            MultiValue.decompose (squareShape atom) b) of-         (Vec2 a0 a1, Vec2 b0 b1) ->-            MultiValue.addPhis bb a0 b0 >>-            MultiValue.addPhis bb a1 b1--type instance-   MultiValue.Decomposed f (Vec2 tag pat) =-      Vec2 tag (MultiValue.Decomposed f pat)-type instance-   MultiValue.PatternTuple (Vec2 tag pat) =-      Vec2 tag (MultiValue.PatternTuple pat)--instance (MultiValue.Compose n) => MultiValue.Compose (Vec2 tag n) where-   type Composed (Vec2 tag n) = Vec2 tag (MultiValue.Composed n)-   compose (Vec2 n m) =-      case (MultiValue.compose n, MultiValue.compose m) of-         (MultiValue.Cons rn, MultiValue.Cons rm) ->-            MultiValue.Cons (Vec2 rn rm)--instance (MultiValue.Decompose pn) => MultiValue.Decompose (Vec2 tag pn) where-   decompose (Vec2 pn pm) (MultiValue.Cons (Vec2 n m)) =-      Vec2-         (MultiValue.decompose pn (MultiValue.Cons n))-         (MultiValue.decompose pm (MultiValue.Cons m))--instance (MultiMem.C i) => MultiMem.C (Vec2 tag i) where-   type Struct (Vec2 tag i) =-         LLVM.Struct (MultiMem.Struct i, (MultiMem.Struct i, ()))-   decompose nm =-      Monad.lift2 zipShape-         (MultiMem.decompose =<< LLVM.extractvalue nm TypeNum.d0)-         (MultiMem.decompose =<< LLVM.extractvalue nm TypeNum.d1)-   compose nm =-      case unzipShape nm of-         Vec2 n m -> do-            sn <- MultiMem.compose n-            sm <- MultiMem.compose m-            rn <- LLVM.insertvalue (LLVM.value LLVM.undef) sn TypeNum.d0-            LLVM.insertvalue rn sm TypeNum.d1---unzipShape :: MultiValue.T (Vec2 tag n) -> Vec2 tag (MultiValue.T n)-unzipShape = MultiValue.decompose (squareShape atom)--zipShape :: MultiValue.T n -> MultiValue.T n -> MultiValue.T (Vec2 tag n)-zipShape y x = MultiValue.compose $ Vec2 y x--instance (tag ~ ShapeTag, Shape.C i) => Shape.C (Vec2 tag i) where-   type Index (Vec2 tag i) = Index2 (Shape.Index i)-   intersectCode a b =-      case (unzipShape a, unzipShape b) of-         (Vec2 an am, Vec2 bn bm) ->-            Monad.lift2 zipShape-               (Shape.intersectCode an bn)-               (Shape.intersectCode am bm)-   sizeCode nm =-      case unzipShape nm of-         Vec2 n m ->-            join $ Monad.lift2 A.mul (Shape.sizeCode n) (Shape.sizeCode m)-   size (Vec2 n m) = Shape.size n * Shape.size m-   flattenIndexRec nm ij =-      case (unzipShape nm, unzipShape ij) of-         (Vec2 n m, Vec2 i j) -> do-            (ns, il) <- Shape.flattenIndexRec n i-            (ms, jl) <- Shape.flattenIndexRec m j-            Monad.lift2 (,)-               (A.mul ns ms)-               (A.add jl =<< A.mul ms il)-   loop code nm =-      case unzipShape nm of-         Vec2 n m ->-            Shape.loop (\i -> Shape.loop (\j -> code (zipShape i j)) m) n---instance (Expr.Compose n) => Expr.Compose (Vec2 tag n) where-   type Composed (Vec2 tag n) = Vec2 tag (Expr.Composed n)-   compose (Vec2 n m) = Expr.lift2 zipShape (Expr.compose n) (Expr.compose m)--instance (Expr.Decompose p) => Expr.Decompose (Vec2 tag p) where-   decompose (Vec2 pn pm) vec =-      Vec2-         (Expr.decompose pn (verticalVal vec))-         (Expr.decompose pm (horizontalVal vec))--verticalVal, horizontalVal :: (Expr.Value val) => val (Vec2 tag n) -> val n-verticalVal = Expr.lift1 (MultiValue.lift1 vertical)-horizontalVal = Expr.lift1 (MultiValue.lift1 horizontal)
src/MatchImageBorders.hs view
@@ -10,7 +10,7 @@ -} module MatchImageBorders where -import KneadShape (Vec2(Vec2), Dim2)+import Knead.Shape (Vec2(Vec2), Dim2)  import qualified Data.Array.Knead.Simple.Physical as Phys