patch-image 0.3.1 → 0.3.2
raw patch · 8 files changed
+271/−282 lines, 8 filesdep +accelerate-llvm-ptxdep +bool8dep +dspdep −accelerate-cudadep −hmatrixdep ~Cabaldep ~acceleratedep ~accelerate-arithmetic
Dependencies added: accelerate-llvm-ptx, bool8, dsp, prelude-compat, semigroups, storable-complex
Dependencies removed: accelerate-cuda, hmatrix
Dependency ranges changed: Cabal, accelerate, accelerate-arithmetic, accelerate-cufft, accelerate-fourier, accelerate-io, accelerate-utility, cassava, containers, filepath, knead, llvm-extra, pqueue, unordered-containers
Files
- Changes.md +8/−0
- patch-image.cabal +33/−21
- src/Accelerate.hs +109/−119
- src/Degree.hs +5/−1
- src/Knead.hs +39/−65
- src/Knead/Shape.hs +7/−1
- src/LinearAlgebra.hs +61/−69
- src/MatchImageBorders.hs +9/−6
Changes.md view
@@ -1,5 +1,13 @@ # Change log for the `patch-image` package +## 0.3.2:++ * Use package `dsp` instead of `hmatrix`+ for the solution of the linear least squares problem+ for determining absolute coordinates from image pair displacements.+ This removes dependency from LAPACK and GSL+ and makes the code a bit simpler.+ ## 0.3.1: * Speed up computation by moving more stuff to Knead/LLVM.
patch-image.cabal view
@@ -1,5 +1,5 @@ Name: patch-image-Version: 0.3.1+Version: 0.3.2 License: BSD3 License-File: LICENSE Author: Henning Thielemann <haskell@henning-thielemann.de>@@ -41,7 +41,7 @@ README.md Source-Repository this- Tag: 0.3.1+ Tag: 0.3.2 Type: darcs Location: http://hub.darcs.net/thielema/patch-image/ @@ -83,28 +83,32 @@ If flag(llvm) Build-Depends:- knead >=0.2.2 && <0.3,- llvm-extra >=0.7 && <0.8,+ knead >=0.3 && <0.4,+ llvm-extra >=0.8 && <0.9, llvm-tf >=3.1 && <3.2, tfp >=1.0 && <1.1, JuicyPixels >=2.0 && <3.3,- hmatrix >=0.15 && <0.16,+ dsp >=0.2.4 && <0.3, vector >=0.10 && <0.13,- pqueue >=1.2 && <1.4,+ pqueue >=1.2 && <1.5, enumset >=0.0.4 && <0.1,- containers >=0.4.2 && <0.6,+ containers >=0.4.2 && <0.7,+ semigroups >=0.1 && <1.0, fft >=0.1.7 && <0.2,+ storable-complex >=0.2.2 && <0.3, storable-tuple >=0.0.3 && <0.1,+ bool8 >=0.0 && <0.1, carray >=0.1.5 && <0.2, array >=0.4 && <0.6,- cassava >=0.4.5 && <0.5,- unordered-containers >=0.2.5 && <0.2.9,+ cassava >=0.4.5 && <0.5.2,+ unordered-containers >=0.2.5 && <0.2.10, bytestring >=0.9.2 && <0.11, explicit-exception >=0.1.7 && <0.2,- Cabal >=1.18 && <2,- filepath >=1.3 && <1.4,+ Cabal >=1.18 && <3,+ filepath >=1.3 && <1.5, non-empty >=0.2 && <0.4, utility-ht >=0.0.13 && <0.1,+ prelude-compat ==0.0.*, base >=4 && <5 Else Buildable: False@@ -126,19 +130,27 @@ If flag(cuda) Build-Depends:- 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.17,- accelerate-io >=0.15 && <0.16,- accelerate >=0.15 && <0.16,+ accelerate-fourier >=1.0 && <1.1,+ accelerate-arithmetic >=1.0 && <1.1,+ accelerate-utility >=1.0 && <1.1,+ accelerate-cufft >=1.0 && <1.1,+ accelerate-llvm-ptx >=1.1 && <1.2,+ accelerate-io >=1.0 && <1.1,+ accelerate >=1.1 && <1.2, JuicyPixels >=2.0 && <3.3,- hmatrix >=0.15 && <0.16,+ cassava >=0.4.5 && <0.5.2,+ dsp >=0.2.4 && <0.3, gnuplot >=0.5 && <0.6,+ containers >=0.4.2 && <0.7,+ array >=0.4 && <0.6, vector >=0.10 && <0.13,- Cabal >=1.18 && <2,- filepath >=1.3 && <1.4,+ unordered-containers >=0.2.5 && <0.2.10,+ bytestring >=0.9.2 && <0.11,+ enumset >=0.0.4 && <0.1,+ explicit-exception >=0.1.7 && <0.2,+ Cabal >=1.18 && <3,+ filepath >=1.3 && <1.5,+ non-empty >=0.2 && <0.4, utility-ht >=0.0.1 && <0.1, base >=4 && <5 Else
src/Accelerate.hs view
@@ -1,11 +1,11 @@ {-# LANGUAGE TypeOperators #-}+{-# LANGUAGE FlexibleContexts #-} module Main where import qualified Option import qualified State import qualified Arithmetic as Arith-import qualified Complex as Komplex import qualified Degree import LinearAlgebra ( absolutePositionsFromPairDisplacements, fixAtLeastOnePosition,@@ -31,8 +31,8 @@ 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 AComplex-import qualified Data.Array.Accelerate.CUDA.Foreign as CUDAForeign-import qualified Data.Array.Accelerate.CUDA as CUDA+import qualified Data.Array.Accelerate.Data.Bits as ABits+import qualified Data.Array.Accelerate.LLVM.PTX as CUDA import qualified Data.Array.Accelerate.IO as AIO import qualified Data.Array.Accelerate.LinearAlgebra as LinAlg import qualified Data.Array.Accelerate.Utility.Lift.Run as Run@@ -46,8 +46,7 @@ import Data.Array.Accelerate.Utility.Ord (argmaximum) import Data.Array.Accelerate (Acc, Array, Exp, DIM1, DIM2, DIM3,- (:.)((:.)), Z(Z), Any(Any), All(All),- (<*), (<=*), (>=*), (==*), (&&*), (||*), (?), (!), )+ (:.)((:.)), Z(Z), Any(Any), All(All), (?), (!), ) import qualified Graphics.Gnuplot.Advanced as GP import qualified Graphics.Gnuplot.LineSpecification as LineSpec@@ -106,8 +105,7 @@ (SV.length dat) return $ AIO.fromVectors- (Z :. Pic.imageHeight pic :. Pic.imageWidth pic :. 3)- ((), dat)+ (Z :. Pic.imageHeight pic :. Pic.imageWidth pic :. 3) dat _ -> ioError $ userError "unsupported image type" writeImage :: Int -> FilePath -> ColorImage8 -> IO ()@@ -117,7 +115,7 @@ Pic.Image { Pic.imageWidth = width, Pic.imageHeight = height,- Pic.imageData = snd $ AIO.toVectors arr+ Pic.imageData = AIO.toVectors arr } writeGrey :: Int -> FilePath -> Array DIM2 Word8 -> IO ()@@ -127,7 +125,7 @@ Pic.Image { Pic.imageWidth = width, Pic.imageHeight = height,- Pic.imageData = snd $ AIO.toVectors arr+ Pic.imageData = AIO.toVectors arr } colorImageExtent :: ColorImage8 -> (Int, Int)@@ -135,12 +133,12 @@ case A.arrayShape pic of Z:.height:.width:._chans -> (width, height) imageFloatFromByte ::- (A.Shape sh, A.Elt a, A.IsFloating a) =>+ (A.Shape sh, A.Floating a, A.FromIntegral Word8 a) => Acc (Array sh Word8) -> Acc (Array sh a) imageFloatFromByte = A.map ((/255) . A.fromIntegral) imageByteFromFloat ::- (A.Shape sh, A.Elt a, A.IsFloating a) =>+ (A.Shape sh, A.RealFloat a) => Acc (Array sh a) -> Acc (Array sh Word8) imageByteFromFloat = A.map (fastRound . (255*) . max 0 . min 1) @@ -172,21 +170,30 @@ arr -fastRound ::- (A.Elt i, A.IsIntegral i, A.Elt a, A.IsFloating a) => Exp a -> Exp i+fastRound :: (A.Elt i, A.IsIntegral i, A.RealFloat a) => Exp a -> Exp i fastRound x = A.floor (x+0.5) floatArray :: Acc (Array sh Float) -> Acc (Array sh Float) floatArray = id -splitFraction :: (A.Elt a, A.IsFloating a) => Exp a -> (Exp Int, Exp a)+splitFraction ::+ (A.RealFloat a, A.FromIntegral Int a) => Exp a -> (Exp Int, Exp a) splitFraction x = let i = A.floor x in (i, x - A.fromIntegral i) +target :: CUDA.PTX+target = unsafePerformIO CUFFT.getBestTarget +cudaRun :: (A.Arrays a) => Acc a -> a+cudaRun = CUDA.runWith target++cudaRun1 :: (A.Arrays a, A.Arrays b) => (Acc a -> Acc b) -> a -> b+cudaRun1 = CUDA.run1With target++ type Channel ix = Array (ix :. Int :. Int) type Plane = Channel Z @@ -209,7 +216,7 @@ in arr ! A.lift (ix :. yc :. xc) indexFrac ::- (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) =>+ (A.Slice ix, A.Shape ix, A.RealFloat a, A.FromIntegral Int a) => Acc (Channel ix a) -> Exp ix :. Exp a :. Exp a -> Exp a indexFrac arr (ix:.y:.x) = let (xi,xf) = splitFraction x@@ -230,7 +237,7 @@ rotateStretchMoveCoords ::- (A.Elt a, A.IsFloating a) =>+ (A.Floating a, A.FromIntegral Int a) => (Exp a, Exp a) -> (Exp a, Exp a) -> (Exp Int, Exp Int) ->@@ -242,15 +249,15 @@ in A.lift $ trans (A.fromIntegral xdst, A.fromIntegral ydst) inBox ::- (A.Elt a, A.IsNum a, A.IsScalar a) =>+ (A.Num a, A.Ord a) => (Exp a, Exp a) -> (Exp a, Exp a) -> Exp Bool inBox (width,height) (x,y) =- 0<=*x &&* x<*width &&* 0<=*y &&* y<*height+ 0 A.<= x A.&& x A.< width A.&& 0 A.<= y A.&& y A.< height validCoords ::- (A.Elt a, A.IsFloating a) =>+ (A.RealFloat a) => (Exp Int, Exp Int) -> Acc (Channel Z (a, a)) -> Acc (Channel Z Bool)@@ -269,7 +276,7 @@ and then moves the picture. -} rotateStretchMove ::- (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) =>+ (A.Slice ix, A.Shape ix, A.RealFloat a, A.FromIntegral Int a) => (Exp a, Exp a) -> (Exp a, Exp a) -> ExpDIM2 ix -> Acc (Channel ix a) ->@@ -289,7 +296,7 @@ rotateLeftTop ::- (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) =>+ (A.Slice ix, A.Shape ix, A.RealFloat a, A.FromIntegral Int a) => (Exp a, Exp a) -> Acc (Channel ix a) -> ((Acc (A.Scalar a), Acc (A.Scalar a)), Acc (Channel ix a)) rotateLeftTop rot arr =@@ -302,7 +309,7 @@ (chans :. A.ceiling (bottom-top) :. A.ceiling (right-left)) arr) rotate ::- (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) =>+ (A.Slice ix, A.Shape ix, A.RealFloat a, A.FromIntegral Int a) => (Exp a, Exp a) -> Acc (Channel ix a) -> Acc (Channel ix a) rotate rot arr = snd $ rotateLeftTop rot arr@@ -321,7 +328,7 @@ Degree Float -> ColorImage8 -> (ColorImage8, Array DIM1 Float) rotateHistogram = let rot =- Run.with CUDA.run1 $ \orient arr ->+ Run.with cudaRun1 $ \orient arr -> let rotated = rotate orient $ separateChannels $ imageFloatFromByte arr@@ -357,9 +364,7 @@ -differentiate ::- (A.Elt a, A.IsNum a) =>- Acc (Array DIM1 a) -> Acc (Array DIM1 a)+differentiate :: (A.Num a) => Acc (Array DIM1 a) -> Acc (Array DIM1 a) differentiate arr = let size = A.unindex1 $ A.shape arr in A.generate (A.index1 (size-1)) $ \i ->@@ -368,7 +373,7 @@ scoreRotation :: Degree Float -> ColorImage8 -> Float scoreRotation = let rot =- Run.with CUDA.run1 $ \orient arr ->+ Run.with cudaRun1 $ \orient arr -> A.sum $ A.map (^(2::Int)) $ differentiate $ rowHistogram $ rotate orient $ separateChannels $ imageFloatFromByte arr in \angle arr -> Acc.the $ rot (Degree.cis angle) arr@@ -378,18 +383,16 @@ Key.maximum (flip scoreRotation pic) angles -magnitudeSqr :: (A.Elt a, A.IsNum a) => Exp (Complex a) -> Exp a+magnitudeSqr :: (A.Num a) => Exp (Complex a) -> Exp a magnitudeSqr = Exp.modify (expr:+expr) $ \(r:+i) -> r*r+i*i fourierTransformationRun :: ColorImage8 -> IO (Array DIM2 Word8) fourierTransformationRun pic = do let (shape@(Z:.height:.width):._) = A.arrayShape pic- plan <-- CUDAForeign.inDefaultContext $- CUFFT.plan2D CUFFT.forwardReal shape+ plan <- CUFFT.plan2D target CUFFT.forwardReal shape let trans =- Run.with CUDA.run1 $ \arr ->+ Run.with cudaRun1 $ \arr -> imageByteFromFloat $ A.map (1e-9*) $ A.zipWith (*)@@ -423,7 +426,7 @@ thus I leave it as it is. -} scoreSlopes ::- (A.Elt a, A.IsFloating a) =>+ (A.RealFloat a, A.FromIntegral Int a) => (Exp Int, Exp Int) -> Acc (Channel Z (Complex a)) -> Acc (Array DIM1 a) scoreSlopes (minX, maxX) arr =@@ -453,16 +456,14 @@ (Degree Float, Degree Float) -> ColorImage8 -> IO (Degree Float) radonAngle (minAngle,maxAngle) pic = do let (shape@(Z :. height :. _width):._) = A.arrayShape pic- plan <-- CUDAForeign.inDefaultContext $- CUFFT.plan2D CUFFT.forwardReal shape+ plan <- CUFFT.plan2D target CUFFT.forwardReal shape let height2 = fromIntegral (div height 2) let slope w = tan (Degree.toRadian w) * height2 let minX = floor $ slope minAngle let maxX = ceiling $ slope maxAngle let angle s = Degree.fromRadian $ atan (s/height2) let trans =- Run.with CUDA.run1 $ \arr ->+ Run.with cudaRun1 $ \arr -> A.map A.snd $ argmaximum $ Arrange.mapWithIndex (\ix s -> A.lift (s, A.unindex1 ix)) $ scoreSlopes (A.constant minX, A.constant maxX) $@@ -475,7 +476,7 @@ rotateManifest :: Degree Float -> ColorImage8 -> Array DIM3 Float rotateManifest = let rot =- Run.with CUDA.run1 $ \orient arr ->+ Run.with cudaRun1 $ \orient arr -> rotate orient $ separateChannels $ imageFloatFromByte arr in \angle arr -> rot (Degree.cis angle) arr @@ -484,7 +485,7 @@ Int -> (Degree Float, ColorImage8) -> ((Float,Float), Plane Float) prepareOverlapMatching = let rot =- Run.with CUDA.run1 $ \radius orient arr ->+ Run.with cudaRun1 $ \radius orient arr -> rotateLeftTop orient $ (if True then highpass radius@@ -497,7 +498,7 @@ ceilingPow2 :: Exp Int -> Exp Int ceilingPow2 n =- A.setBit 0 $ A.ceiling $ logBase 2 (fromIntegral n :: Exp Double)+ ABits.setBit 0 $ A.ceiling $ logBase 2 (fromIntegral n :: Exp Double) pad :: (A.Elt a) =>@@ -506,12 +507,12 @@ let (height, width) = A.unlift $ A.unindex2 $ A.shape arr in A.generate sh $ \p -> let (y, x) = A.unlift $ A.unindex2 p- in (y<*height &&* x<*width)+ in (y A.< height A.&& x A.< width) ? (arr ! A.index2 y x, a) mulConj ::- (A.Elt a, A.IsFloating a) =>+ (A.RealFloat a, A.FromIntegral Int a) => Exp (Complex a) -> Exp (Complex a) -> Exp (Complex a) mulConj x y = x * AComplex.conjugate y @@ -520,15 +521,14 @@ (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+ CUFFT.transform $ unsafePerformIO $ CUFFT.plan2D target 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+ A.use $ cudaRun1 (fft2DGen mode $ A.arrayShape arr) arr fft2D :: (A.Elt e, CUFFT.Real e, A.Elt a, A.Elt b) =>@@ -547,13 +547,14 @@ width = ceilingPow2 $ widthx + widthy height = ceilingPow2 $ heightx + heighty sh = A.index2 height width- forward z = fft2DPlain CUFFT.forwardReal $ CUDA.run $ pad 0 sh z- in fft2DPlain CUFFT.inverseReal $ CUDA.run $+ forward z = fft2DPlain CUFFT.forwardReal $ cudaRun $ pad 0 sh z+ in fft2DPlain CUFFT.inverseReal $ cudaRun $ A.zipWith mulConj (forward x) (forward y) removeDCOffset ::- (A.Elt a, A.IsFloating a) => Acc (Channel Z a) -> Acc (Channel Z a)+ (A.Floating a, A.FromIntegral Int a) =>+ Acc (Channel Z a) -> Acc (Channel Z a) removeDCOffset arr = let sh = A.shape arr (_z :. height :. width) = unliftDim2 sh@@ -567,20 +568,19 @@ because we already padded the images with zeros. -} clearDCCoefficient ::- (A.Elt a, A.IsFloating a) =>+ (A.RealFloat a, A.FromIntegral Int a) => Acc (Array DIM2 (Complex a)) -> Acc (Array DIM2 (Complex a)) clearDCCoefficient arr = A.generate (A.shape arr) $ \p -> let (_z:.y:.x) = unliftDim2 p- in x==*0 ||* y==*0 ? (0, arr!p)+ in x A.== 0 A.|| y A.== 0 ? (0, arr!p) lowpass, highpass ::- (A.Elt a, A.IsFloating a) =>- Exp Int -> Acc (Channel Z a) -> Acc (Channel Z a)+ (A.Floating a) => Exp Int -> Acc (Channel Z a) -> Acc (Channel Z a) lowpass count = Loop.nest count $- A.stencil (\(a,m,b) -> smooth3 (smooth3 a, smooth3 m, smooth3 b)) A.Clamp+ A.stencil (\(a,m,b) -> smooth3 (smooth3 a, smooth3 m, smooth3 b)) A.clamp highpass count arr = A.zipWith (-) arr $ lowpass count arr@@ -619,7 +619,7 @@ wrap :: Exp Int -> Exp Int -> Exp Int -> Exp Int-wrap size split c = c<*split ? (c, c-size)+wrap size split c = c A.< split ? (c, c-size) displacementMap :: Exp Int -> Exp Int -> Exp DIM2 -> Acc (Channel Z (Int, Int))@@ -643,7 +643,7 @@ that are actually digitalization artifacts. -} minimumOverlapScores ::- (A.Elt a, A.IsFloating a, A.IsScalar a) =>+ (A.Floating a, A.IsScalar a) => ((Exp Int, Exp Int) -> Exp a -> Exp a) -> Exp Int -> (Exp Int, Exp Int) -> (Exp Int, Exp Int) -> Acc (Channel Z (a, (Int, Int))) ->@@ -653,7 +653,7 @@ (Exp.modify (expr,(expr,expr)) $ \(v, dp@(dx,dy)) -> let clipWidth = min widtha (widthb + dx) - max 0 dx clipHeight = min heighta (heightb + dy) - max 0 dy- in ((clipWidth >=* minOverlap &&* clipHeight >=* minOverlap)+ in ((clipWidth A.>= minOverlap A.&& clipHeight A.>= minOverlap) ? (weight (clipWidth, clipHeight) v, 0), dp))@@ -694,7 +694,7 @@ allOverlapsRun :: DIM2 -> Float -> Plane Float -> Plane Float -> Plane Word8 allOverlapsRun padExtent =- Run.with CUDA.run1 $ \minOverlap picA picB ->+ Run.with cudaRun1 $ \minOverlap picA picB -> imageByteFromFloat $ -- A.map (2*) $ A.map (0.0001*) $@@ -704,13 +704,13 @@ DIM2 -> Float -> Plane Float -> Plane Float -> (Float, (Int, Int)) optimalOverlap padExtent = let run =- Run.with CUDA.run1 $ \minimumOverlap a b ->+ Run.with cudaRun1 $ \minimumOverlap a b -> argmaximum $ allOverlaps padExtent minimumOverlap a b in \overlap a b -> Acc.the $ run overlap a b shrink ::- (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) =>+ (A.Slice ix, A.Shape ix, A.Floating a, A.FromIntegral Int a) => GenDIM2 (Exp Int) -> Acc (Channel ix a) -> Acc (Channel ix a) shrink (_:.yk:.xk) arr = let (shape:.height:.width) = unliftDim2 $ A.shape arr@@ -752,7 +752,7 @@ DIM2 -> Float -> Plane Float -> Plane Float -> (Float, (Int, Int)) optimalOverlapBig padExtent = let run =- Run.with CUDA.run1 $ \minimumOverlap a b ->+ Run.with cudaRun1 $ \minimumOverlap a b -> let factors@(_z:.yk:.xk) = shrinkFactors (A.floor, A.fromIntegral) padExtent minimumOverlap@@ -804,7 +804,7 @@ optimalOverlapBigFine padExtent@(Z:.heightPad:.widthPad) = let overlaps = allOverlaps padExtent run =- Run.with CUDA.run1 $ \minimumOverlap a b ->+ Run.with cudaRun1 $ \minimumOverlap a b -> let shapeA = A.unlift $ A.shape a shapeB = A.unlift $ A.shape b factors@(_z:.yk:.xk) =@@ -848,20 +848,20 @@ [(Float, (Int, Int), (Int, Int))] optimalOverlapBigMulti padExtent (Z:.heightStamp:.widthStamp) numCorrs = let overlapShrunk =- Run.with CUDA.run1 $+ Run.with cudaRun1 $ \minimumOverlap factors a b -> argmaximum $ allOverlaps padExtent minimumOverlap (shrink factors a) (shrink factors b) diffShrunk =- Run.with CUDA.run1 $+ Run.with cudaRun1 $ \shrunkd factors a b -> overlapDifference shrunkd (shrink factors a) (shrink factors b) allOverlapsFine = allOverlaps (Z :. 2*heightStamp :. 2*widthStamp) overlapFine =- Run.with CUDA.run1 $+ Run.with cudaRun1 $ \minimumOverlap a b anchorA@(leftA, topA) anchorB@(leftB, topB) extent@(width,height) -> let addCoarsePos =@@ -922,7 +922,7 @@ overlapDifference ::- (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) =>+ (A.Slice ix, A.Shape ix, A.Floating a, A.FromIntegral Int a) => (Exp Int, Exp Int) -> Acc (Channel ix a) -> Acc (Channel ix a) -> Acc (A.Scalar a) overlapDifference (dx,dy) a b =@@ -947,7 +947,7 @@ (Int, Int) -> Plane Float -> Plane Float -> Float overlapDifferenceRun =- let diff = Run.with CUDA.run1 overlapDifference+ let diff = Run.with cudaRun1 overlapDifference in \d a b -> Acc.the $ diff d a b @@ -983,7 +983,7 @@ let rotat (rot,pic) = rotate rot $ separateChannels $ imageFloatFromByte pic in (\f d (a,b) -> f d (mapFst Degree.cis a, mapFst Degree.cis b)) $- Run.with CUDA.run1 $+ Run.with cudaRun1 $ \(dx,dy) (a,b) -> imageByteFromFloat $ interleaveChannels $ overlap2 (dx, dy) (rotat a, rotat b)@@ -994,14 +994,14 @@ ix :. Int :. Int -> (Plane Int, Channel ix Float) emptyCountCanvas =- Run.with CUDA.run1 $ \sh ->+ Run.with cudaRun1 $ \sh -> let (_ix :. height :. width) = unliftDim2 sh in (A.fill (A.lift $ Z:.height:.width) 0, A.fill sh 0) addToCountCanvas ::- (A.Slice ix, A.Shape ix, A.Elt a, A.IsNum a) =>+ (A.Slice ix, A.Shape ix, A.Num a, A.FromIntegral Int a) => (Acc (Plane Bool), Acc (Channel ix a)) -> (Acc (Plane Int), Acc (Channel ix a)) -> (Acc (Plane Int), Acc (Channel ix a))@@ -1016,7 +1016,7 @@ (Plane Int, Channel DIM1 Float) -> (Plane Int, Channel DIM1 Float) updateCountCanvas =- Run.with CUDA.run1 $+ Run.with cudaRun1 $ \(rot, mov, pic) (count,canvas) -> addToCountCanvas (rotateStretchMove rot mov (unliftDim2 $ A.shape canvas) $@@ -1025,7 +1025,7 @@ finalizeCountCanvas :: (Plane Int, Channel DIM1 Float) -> ColorImage8 finalizeCountCanvas =- Run.with CUDA.run1 $+ Run.with cudaRun1 $ \(count, canvas) -> imageByteFromFloat $ interleaveChannels $ A.zipWith (/) canvas $@@ -1044,30 +1044,30 @@ in A.cond xb (A.lift (True, A.cond yb (f xv yv) xv)) y maskedMinimum ::- (A.Shape ix, A.Elt a, A.IsScalar a) =>+ (A.Shape ix, A.Ord a) => LinAlg.Vector ix (Bool, a) -> LinAlg.Scalar ix (Bool, a) maskedMinimum = A.fold1 (maybePlus min) maskedMaximum ::- (A.Shape ix, A.Elt a, A.IsScalar a) =>+ (A.Shape ix, A.Ord a) => LinAlg.Vector ix (Bool, a) -> LinAlg.Scalar ix (Bool, a) maskedMaximum = A.fold1 (maybePlus max) project ::- (A.Elt a, A.IsFloating a) =>+ (A.RealFloat a, A.FromIntegral Int a) => Point2 (Exp a) -> (Point2 (Exp a), Point2 (Exp a)) -> (Exp Bool, Point2 (Exp a)) project x ab = let (r, y) = projectPerp x ab- in (0<=*r &&* r<=*1, y)+ in (0 A.<= r A.&& r A.<= 1, y) distanceMapEdges ::- (A.Elt a, A.IsFloating a) =>+ (A.RealFloat a, A.FromIntegral Int a) => Exp DIM2 -> Acc (Array DIM1 ((a,a),(a,a))) -> Acc (Channel Z a) distanceMapEdges sh edges = A.map (Exp.modify (expr,expr) $ \(valid, dist) -> valid ? (dist, 0)) $@@ -1081,13 +1081,13 @@ distanceMapEdgesRun :: DIM2 -> Array DIM1 ((Float,Float),(Float,Float)) -> Plane Word8 distanceMapEdgesRun =- Run.with CUDA.run1 $ \sh ->+ Run.with cudaRun1 $ \sh -> imageByteFromFloat . A.map (0.01*) . distanceMapEdges sh type Geometry a = Arith.Geometry Int a distanceMapBox ::- (A.Elt a, A.IsFloating a) =>+ (A.RealFloat a, A.FromIntegral Int a) => Exp DIM2 -> Exp (Geometry a) -> Acc (Channel Z (Bool, (((a,(a,a)), (a,(a,a))), ((a,(a,a)), (a,(a,a))))))@@ -1138,7 +1138,7 @@ distanceMapBoxRun :: DIM2 -> Geometry Float -> Plane Word8 distanceMapBoxRun =- Run.with CUDA.run1 $ \sh geom ->+ Run.with cudaRun1 $ \sh geom -> scaleDistanceMapGeom geom $ A.map (Exp.modify (expr,expr) $ \(valid, dist) -> valid ? (dist, 0)) $ maskedMinimum $@@ -1161,12 +1161,12 @@ containedAnywhere ::- (A.Elt a, A.IsFloating a) =>+ (A.RealFloat a) => Acc (Array DIM1 (Geometry a)) -> Acc (Array DIM3 (a,a)) -> Acc (Array DIM3 Bool) containedAnywhere geoms arr =- A.fold1 (||*) $+ A.fold1 (A.||) $ breakFusion $ outerVector (Exp.modify2 (expr,expr) ((expr,expr),(expr,expr),(expr,expr)) $@@ -1177,7 +1177,7 @@ distanceMapContained ::- (A.IsFloating a, A.Elt a) =>+ (A.RealFloat a, A.FromIntegral Int a) => Exp DIM2 -> Exp (Geometry a) -> Acc (Array DIM1 (Geometry a)) ->@@ -1189,19 +1189,19 @@ maskedMinimum $ A.zipWith (Exp.modify2 expr (expr,(expr,expr)) $ \c (b,(dist,_)) ->- (c&&*b, dist))+ (c A.&& b, dist)) contained distMap distanceMapContainedRun :: DIM2 -> Geometry Float -> [Geometry Float] -> Plane Word8 distanceMapContainedRun = let distances =- Run.with CUDA.run1 $+ Run.with cudaRun1 $ \sh this -> scaleDistanceMapGeom this . distanceMapContained sh this in \sh this others -> distances sh this $ array1FromList others scaleDistanceMapGeom ::- (A.IsFloating a, A.Elt a, A.Elt b, A.Shape ix) =>+ (A.RealFloat a, A.FromIntegral Int a, A.Elt b, A.Shape ix) => Exp (Geometry b) -> Acc (Array ix a) -> Acc (Array ix Word8) scaleDistanceMapGeom this = let scale = (4/) $ A.fromIntegral $ A.uncurry min $ Exp.thd3 this@@ -1209,14 +1209,13 @@ pixelCoordinates ::- (A.Elt a, A.IsFloating a) =>- Exp DIM2 -> Acc (Channel Z (a,a))+ (A.RealFloat a, A.FromIntegral Int a) => Exp DIM2 -> Acc (Channel Z (a,a)) pixelCoordinates sh = A.generate sh $ Exp.modify (expr:.expr:.expr) $ \(_z:.y:.x) -> (A.fromIntegral x, A.fromIntegral y) distanceMapPoints ::- (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) =>+ (A.Slice ix, A.Shape ix, A.RealFloat a) => Acc (Array ix (a,a)) -> Acc (Array DIM1 (a,a)) -> Acc (Array ix a)@@ -1229,12 +1228,13 @@ distanceMapPointsRun :: DIM2 -> [Point2 Float] -> Plane Word8 distanceMapPointsRun = let distances =- Run.with CUDA.run1 $+ Run.with cudaRun1 $ \sh -> scaleDistanceMap . distanceMapPoints (pixelCoordinates sh) in \sh points -> distances sh $ array1FromList points scaleDistanceMap ::- (A.Elt a, A.IsFloating a) => Acc (Channel Z a) -> Acc (Channel Z Word8)+ (A.RealFloat a, A.FromIntegral Int a) =>+ Acc (Channel Z a) -> Acc (Channel Z Word8) scaleDistanceMap arr = let scale = case Exp.unlift (expr:.expr:.expr) (A.shape arr) of@@ -1258,7 +1258,7 @@ and chose the minimal distance. -} distanceMap ::- (A.Elt a, A.IsFloating a) =>+ (A.RealFloat a, A.FromIntegral Int a) => Exp DIM2 -> Exp (Geometry a) -> Acc (Array DIM1 (Geometry a)) ->@@ -1277,7 +1277,7 @@ Plane Word8 distanceMapRun = let distances =- Run.with CUDA.run1 $+ Run.with cudaRun1 $ \sh this others -> scaleDistanceMap . distanceMap sh this others in \sh this others points -> distances sh this@@ -1286,7 +1286,7 @@ distanceMapGamma ::- (A.Elt a, A.IsFloating a) =>+ (A.RealFloat a, A.FromIntegral Int a) => Exp a -> Exp DIM2 -> Exp (Geometry a) ->@@ -1302,14 +1302,14 @@ ix :. Int :. Int -> (Plane Float, Channel ix Float) emptyWeightedCanvas =- Run.with CUDA.run1 $ \sh ->+ Run.with cudaRun1 $ \sh -> let (_ix :. height :. width) = unliftDim2 sh in (A.fill (A.lift $ Z:.height:.width) 0, A.fill sh 0) addToWeightedCanvas ::- (A.Slice ix, A.Shape ix, A.Elt a, A.IsNum a) =>+ (A.Slice ix, A.Shape ix, A.Num a) => (Acc (Channel Z a), Acc (Channel ix a)) -> (Acc (Channel Z a), Acc (Channel ix a)) -> (Acc (Channel Z a), Acc (Channel ix a))@@ -1328,7 +1328,7 @@ (Plane Float, Channel DIM1 Float) updateWeightedCanvasMerged = let update =- Run.with CUDA.run1 $+ Run.with cudaRun1 $ \this others points pic (weightSum,canvas) -> let (rot, mov, _) = Exp.unlift ((expr,expr), (expr,expr), expr) this@@ -1352,9 +1352,9 @@ (Plane Float, Channel DIM1 Float) -> (Plane Float, Channel DIM1 Float) updateWeightedCanvas =- let distances = Run.with CUDA.run1 distanceMapGamma+ let distances = Run.with cudaRun1 distanceMapGamma update =- Run.with CUDA.run1 $+ Run.with cudaRun1 $ \this pic dist (weightSum,canvas) -> let (rot, mov, _) = Exp.unlift ((expr,expr), (expr,expr), expr) this@@ -1379,10 +1379,10 @@ (Plane Float, Channel DIM1 Float) -> (Plane Float, Channel DIM1 Float) updateWeightedCanvasSplit =- let update = Run.with CUDA.run1 addToWeightedCanvas- distances = Run.with CUDA.run1 distanceMap+ let update = Run.with cudaRun1 addToWeightedCanvas+ distances = Run.with cudaRun1 distanceMap rotated =- Run.with CUDA.run1 $+ Run.with cudaRun1 $ \sh rot mov pic -> snd $ rotateStretchMove rot mov (unliftDim2 sh) $ separateChannels $ imageFloatFromByte pic@@ -1398,7 +1398,7 @@ finalizeWeightedCanvas :: (Plane Float, Channel DIM1 Float) -> ColorImage8 finalizeWeightedCanvas =- Run.with CUDA.run1 $+ Run.with cudaRun1 $ \(weightSum, canvas) -> imageByteFromFloat $ interleaveChannels $ A.zipWith (/) canvas $@@ -1505,8 +1505,7 @@ let (errdx,errdy) = mapPair (maximum0, maximum0) $ unzip $ zipWith- (\(dpx,dpy) (dx,dy) ->- (abs $ dpx - realToFrac dx, abs $ dpy - realToFrac dy))+ (\(dpx,dpy) (dx,dy) -> (abs $ dpx - dx, abs $ dpy - dy)) dps (map snd overlaps) info $@@ -1516,10 +1515,7 @@ ++ printf "maximum vertical error: %f\n" errdy - return- (map picPath pics,- map picColored pics,- map (flip (,) 1) $ map (mapPair (realToFrac, realToFrac)) poss)+ return (map picPath pics, map picColored pics, map (flip (,) 1) poss) processOverlapRotate ::@@ -1610,13 +1606,7 @@ dpx dpy xa ya xb yb) dps overlaps - return- (map picPath pics,- map picColored pics,- map- (mapPair- (mapPair (realToFrac, realToFrac), Komplex.map realToFrac))- posRots)+ return (map picPath pics, map picColored pics, posRots) processRotation ::@@ -1659,14 +1649,14 @@ let pic0 : pic1 : _ = map snd rotated size = (Z:.512:.1024 :: DIM2) writeGrey (Option.quality opt) "/tmp/padded.jpeg" $- CUDA.run1 (imageByteFromFloat . pad 0 (A.lift size)) pic0+ cudaRun1 (imageByteFromFloat . pad 0 (A.lift size)) pic0 writeGrey (Option.quality opt) "/tmp/spectrum.jpeg" $- CUDA.run $ imageByteFromFloat $ A.map AComplex.magnitude $+ cudaRun $ imageByteFromFloat $ A.map AComplex.magnitude $ fft2DPlain CUFFT.forwardReal $- CUDA.run1 (pad 0 (A.lift size)) $+ cudaRun1 (pad 0 (A.lift size)) $ pic0 writeGrey (Option.quality opt) "/tmp/convolution.jpeg" $- CUDA.run $ imageByteFromFloat $ A.map (0.000001*) $+ cudaRun $ imageByteFromFloat $ A.map (0.000001*) $ correlatePadded size (A.use pic0) (A.use pic1) return $
src/Degree.hs view
@@ -4,14 +4,18 @@ import qualified Data.Csv as Csv import Data.Monoid (Monoid, mempty, mappend)+import Data.Semigroup (Semigroup, (<>)) newtype Degree a = Degree {getDegree :: a} deriving (Eq, Show) +instance (Num a) => Semigroup (Degree a) where+ Degree x <> Degree y = Degree $ x+y+ instance (Num a) => Monoid (Degree a) where mempty = Degree 0- mappend (Degree x) (Degree y) = Degree $ x+y+ mappend = (<>) instance Functor Degree where fmap f (Degree x) = Degree $ f x
src/Knead.hs view
@@ -27,11 +27,11 @@ import qualified Data.Array.Knead.Simple.Physical as Phys import qualified Data.Array.Knead.Simple.ShapeDependent as ShapeDep import qualified Data.Array.Knead.Simple.Symbolic as Symb-import qualified Data.Array.Knead.Index.Nested.Shape as Shape+import qualified Data.Array.Knead.Shape.Nested as Shape import qualified Data.Array.Knead.Expression as Expr import Data.Array.Knead.Simple.Symbolic ((!)) import Data.Array.Knead.Expression- (Exp, (==*), (/=*), (<*), (<=*), (>=*), (&&*))+ (Exp, (==*), (<*), (<=*), (>=*), (||*), (&&*)) import Data.Array.IArray (amap) import Data.Array.CArray (CArray)@@ -52,6 +52,7 @@ import qualified Data.Vector.Storable as SV import Foreign.ForeignPtr (ForeignPtr, castForeignPtr)+import Foreign.Storable.Complex () import qualified System.FilePath as FilePath import qualified System.IO as IO@@ -68,7 +69,7 @@ import qualified Data.Foldable as Fold import qualified Data.List as List import qualified Data.Map as Map-import Data.Function.HT (Id)+import qualified Data.Bool8 as Bool8 import Data.Monoid ((<>)) import Data.Maybe.HT (toMaybe) import Data.Maybe (mapMaybe, isJust, isNothing)@@ -79,9 +80,13 @@ import Data.Tuple.HT (mapPair, mapFst, mapSnd, mapTriple, swap, mapThd3, fst3, uncurry3) import Data.Word (Word8, Word32)+import Data.Bool8 (Bool8) +import Prelude2010+import Prelude () + type SmallSize = Word32 type Plane = Phys.Array Dim2@@ -411,17 +416,17 @@ (Exp a, Exp a) -> Exp Bool inBox (width,height) (x,y) =- Expr.liftM2 MultiValue.and (inRange width x) (inRange height y)+ inRange width x &&* inRange height y validCoords :: (MultiValue.NativeFloating a ar, MultiValue.Field a, MultiValue.Real a, MultiValue.RationalConstant a) => (Exp Size, Exp Size) ->- SymbPlane (a, a) -> SymbPlane MaskBool+ SymbPlane (a, a) -> SymbPlane Bool8 validCoords (width,height) = Symb.map $ Expr.modify (atom,atom) $ \(x,y) ->- maskFromBool $ inBox (width,height) (fastRound x, fastRound y)+ Expr.bool8FromP $ inBox (width,height) (fastRound x, fastRound y) {- | @rotateStretchMove rot mov@@@ -438,7 +443,7 @@ Exp (a, a) -> Exp Dim2 -> SymbPlane v ->- SymbPlane (MaskBool, v)+ SymbPlane (Bool8, v) rotateStretchMove vec rot mov sh img = let coords = rotateStretchMoveCoords rot mov sh (Vec2 heightSrc widthSrc) = Expr.decompose atomDim2 $ Symb.shape img@@ -501,16 +506,8 @@ (Symb.C array, MultiValue.Additive a) => array Dim1 a -> array Dim1 a differentiate xs = Symb.zipWith Expr.sub (tailArr xs) xs -the :: Symb.Array () a -> Exp a-the = Symb.the--fold1All ::- (Shape.C sh, MultiValue.C a) =>- (Exp a -> Exp a -> Exp a) -> Symb.Array sh a -> Exp a-fold1All f = the . Symb.fold1All f- scoreHistogram :: (MultiValue.PseudoRing a) => Symb.Array Dim1 a -> Exp a-scoreHistogram = fold1All Expr.add . Symb.map Expr.sqr . differentiate+scoreHistogram = Symb.fold1All Expr.add . Symb.map Expr.sqr . differentiate runScoreRotation :: IO (Degree Float -> ColorImage8 -> IO Float)@@ -558,7 +555,7 @@ highpassMulti :: IO (Int -> Plane Float -> IO (Plane Float)) highpassMulti = do lp <- lowpassMulti- sub <- RenderP.run $ Symb.zipWith Expr.sub . fixArray+ sub <- RenderP.run $ Symb.zipWith Expr.sub . Symb.fix return $ \n img -> sub img =<< lp n img @@ -623,13 +620,10 @@ liftCArray f a = arrayKneadFromC . f <$> arrayCFromKnead a -fixArray :: Id (Symb.Array sh a)-fixArray = id- prepareOverlapMatching :: IO (Int -> (Degree Float, ColorImage8) -> IO ((Float, Float), Plane Float)) prepareOverlapMatching = do- bright <- RenderP.run $ brightnessPlane . colorImageFloatFromByte . fixArray+ bright <- RenderP.run $ brightnessPlane . colorImageFloatFromByte . Symb.fix hp <- highpassMulti rotat <- RenderP.run $ rotate Arith.vecScalar return $ \radius (angle, img) ->@@ -737,7 +731,7 @@ (Shape.C sh, MultiValue.Comparison a, MultiValue.Select a, MultiValue.Select b) => Symb.Array sh (a, b) -> Exp (a, b)-argmaximum = fold1All argmax+argmaximum = Symb.fold1All argmax optimalOverlap :: Dim2 -> IO (Float -> Plane Float -> Plane Float -> IO (Float, (Size, Size)))@@ -961,7 +955,7 @@ extentOverlap = (widthOverlap,heightOverlap) in Expr.sqrt $ (/(fromInt widthOverlap * fromInt heightOverlap)) $- fold1All (+) $+ Symb.fold1All (+) $ Symb.map Expr.sqr $ Symb.zipWith (-) (clip (leftOverlap,topOverlap) extentOverlap a)@@ -1020,31 +1014,20 @@ RenderP.run $ \sh -> Symb.fill sh (Expr.zip 0 $ Expr.zip3 0 0 0) -type MaskBool = Word8--maskFromBool :: Exp Bool -> Exp MaskBool-maskFromBool = Expr.liftM $ MultiValue.liftM $ LLVM.zext--boolFromMask :: Exp MaskBool -> Exp Bool-boolFromMask = (/=* 0)--intFromBool :: Exp MaskBool -> Exp Word32-intFromBool = Expr.liftM $ MultiValue.liftM $ LLVM.ext- type RotatedImage = ((Float,Float), (Float,Float), ColorImage8) addToCountCanvas :: (MultiValue.PseudoRing a, MultiValue.NativeFloating a ar) => VecExp a v ->- SymbPlane (MaskBool, v) ->+ SymbPlane (Bool8, v) -> SymbPlane (Word32, v) -> SymbPlane (Word32, v) addToCountCanvas vec = Symb.zipWith (Expr.modify2 (atom,atom) (atom,atom) $ \(mask, pic) (count, canvas) ->- (Expr.add (intFromBool mask) count,+ (Expr.add (Expr.intFromBool8 mask) count, Arith.vecAdd vec canvas $- Arith.vecScale vec (fromInt $ intFromBool mask) pic))+ Arith.vecScale vec (Expr.floatFromBool8 mask) pic)) updateCountCanvas :: IO (RotatedImage -> Plane (Word32, YUV Float) ->@@ -1063,14 +1046,14 @@ Symb.map (Expr.modify (atom,atom) $ \(count, pixel) -> Arith.vecScale vecYUV (recip $ fromInt count) pixel) .- fixArray+ Symb.fix diffAbs :: (MultiValue.Real a) => Exp a -> Exp a -> Exp a diffAbs = Expr.liftM2 $ \x y -> MultiValue.abs =<< MultiValue.sub x y diffWithCanvas ::- IO (RotatedImage -> Plane (YUV Float) -> IO (Plane (MaskBool, Float)))+ IO (RotatedImage -> Plane (YUV Float) -> IO (Plane (Bool8, Float))) diffWithCanvas = RenderP.run $ \(rot, mov, pic) avg -> Symb.zipWith@@ -1088,20 +1071,20 @@ (Expr.modify (atom,atom) $ \(count, pixel) -> Arith.vecScale vecYUV (recip $ fromInt count) pixel) .- fixArray+ Symb.fix emptyCanvas :: IO (Dim2 -> IO ColorImage8) emptyCanvas = RenderP.run $ \sh -> Symb.fill sh (Expr.zip3 0 0 0) addMaskedToCanvas :: IO (RotatedImage ->- Plane MaskBool ->+ Plane Bool8 -> Plane (YUV Word8) -> IO (Plane (YUV Word8))) addMaskedToCanvas = RenderP.run $ \(rot, mov, pic) mask canvas -> Symb.zipWith3 Expr.ifThenElse- (Symb.map boolFromMask mask)+ (Symb.map Expr.boolPFrom8 mask) (Symb.map (yuvByteFromFloat . Expr.snd) $ rotateStretchMove vecYUV rot mov (Symb.shape canvas) $ colorImageFloatFromByte pic)@@ -1116,7 +1099,8 @@ RenderP.run $ \(rot, mov, pic) shape weightCanvas -> addToWeightedCanvas vecYUV (Symb.zipWith- (Expr.modify2 atom (atom,atom) $ \s (b,x) -> (fromInt b * s, x))+ (Expr.modify2 atom (atom,atom) $ \s (b,x) ->+ (Expr.floatFromBool8 b * s, x)) shape $ rotateStretchMove vecYUV rot mov (Symb.shape weightCanvas) $ colorImageFloatFromByte pic)@@ -1231,7 +1215,7 @@ array sh (a,a) -> array sh Bool containedAnywhere geoms array =- Symb.fold1 (Expr.liftM2 MultiValue.or) $+ Symb.fold1 (||*) $ outerProduct (Expr.modify2 (atom,atom) ((atom,atom),(atom,atom),(atom,atom)) $ \(xdst,ydst) (rot, mov, extent) ->@@ -1257,7 +1241,7 @@ maskedMinimum $ Symb.zipWith (Expr.modify2 atom (atom,(atom,atom)) $ \c (b,(dist,_)) ->- (Expr.liftM2 MultiValue.and c b, dist))+ (c &&* b, dist)) contained distMap distanceMapContainedRun ::@@ -1436,7 +1420,7 @@ Symb.map (Expr.modify (atom,atom) $ \(weightSum, pixel) -> Arith.vecScale vecYUV (recip weightSum) pixel) .- fixArray+ Symb.fix data@@ -1544,8 +1528,7 @@ let (errdx,errdy) = mapPair (maximum0, maximum0) $ unzip $ zipWith- (\(dpx,dpy) (dx,dy) ->- (abs $ dpx - realToFrac dx, abs $ dpy - realToFrac dy))+ (\(dpx,dpy) (dx,dy) -> (abs $ dpx - dx, abs $ dpy - dy)) dps (map snd overlaps) info $@@ -1555,10 +1538,7 @@ ++ printf "maximum vertical error: %f\n" errdy - return- (map picPath pics,- map picColored pics,- map (flip (,) 1) $ map (mapPair (realToFrac, realToFrac)) poss)+ return (map picPath pics, map picColored pics, map (flip (,) 1) poss) processOverlapRotate ::@@ -1649,13 +1629,7 @@ dpx dpy xa ya xb yb) dps overlaps - return- (map picPath pics,- map picColored pics,- map- (mapPair- (mapPair (realToFrac, realToFrac), Komplex.map realToFrac))- posRots)+ return (map picPath pics, map picColored pics, posRots) processRotation ::@@ -1695,14 +1669,14 @@ let pic0 : pic1 : _ = map snd rotated size = Vec2 1024 768 makeByteImage <-- RenderP.run $ \k -> imageByteFromFloat . Symb.map (k*) . fixArray+ RenderP.run $ \k -> imageByteFromFloat . Symb.map (k*) . Symb.fix runPad <- RenderP.run pad writeGrey (Option.quality opt) "/tmp/padded.jpeg" =<< (makeByteImage 1 =<< runPad 0 size pic0) runMagnitude <- RenderP.run $ Symb.map (Expr.modify atomComplex $ \(r:+i) -> Expr.sqrt$ r*r+i*i)- . fixArray+ . Symb.fix writeGrey (Option.quality opt) "/tmp/spectrum.jpeg" =<< (makeByteImage 0.1 =<< runMagnitude =<< liftCArray (FFT.dftRCN [0,1]) =<< runPad 0 size pic0)@@ -1815,9 +1789,9 @@ avg <- finalizeCanv sumImg diff <- diffWithCanvas picDiffs <- mapM (flip diff avg) rotMovPics- getSnd <- RenderP.run $ Symb.map Expr.snd . fixArray+ getSnd <- RenderP.run $ Symb.map Expr.snd . Symb.fix lp <- lowpassMulti- masks <- map (amap ((0/=) . fst)) <$> mapM arrayCFromKnead picDiffs+ masks <- map (amap (Bool8.toBool . fst)) <$> mapM arrayCFromKnead picDiffs let smoothRadius = Option.shapeSmooth opt smoothPicDiffs <- mapM (arrayCFromKnead <=< lp smoothRadius <=< getSnd) picDiffs@@ -1840,7 +1814,7 @@ foldM (\canvas (shape, rotMovPic) -> addMasked rotMovPic- (arrayKneadFromC $ amap (fromIntegral . fromEnum) shape)+ (arrayKneadFromC $ amap Bool8.fromBool shape) canvas) emptyPlain (zip shapes rotMovPics) @@ -1851,7 +1825,7 @@ amap (fromIntegral . fromEnum)) shapes forM_ (Option.outputShape opt) $ \format -> do- makeByteImage <- RenderP.run $ imageByteFromFloat . fixArray+ makeByteImage <- RenderP.run $ imageByteFromFloat . Symb.fix forM_ (zip names smoothShapes) $ \(name,shape) -> writeGrey (Option.quality opt) (printf format name) =<< makeByteImage shape
src/Knead/Shape.hs view
@@ -2,11 +2,12 @@ {-# LANGUAGE EmptyDataDecls #-} module Knead.Shape where -import qualified Data.Array.Knead.Index.Nested.Shape as Shape+import qualified Data.Array.Knead.Shape.Nested 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.Iterator as Iter import qualified LLVM.Extra.Arithmetic as A import LLVM.Extra.Multi.Value (atom) @@ -159,6 +160,11 @@ case unzipShape nm of Vec2 n m -> Shape.loop (\i -> Shape.loop (\j -> code (zipShape i j)) m) n+ iterator nm =+ case unzipShape nm of+ Vec2 n m ->+ fmap (uncurry zipShape) $+ Iter.cartesian (Shape.iterator n) (Shape.iterator m) instance (Expr.Compose n) => Expr.Compose (Vec2 tag n) where
src/LinearAlgebra.hs view
@@ -1,19 +1,21 @@ module LinearAlgebra where -import qualified Data.Packed.Matrix as Matrix-import qualified Data.Packed.Vector as Vector-import qualified Data.Packed.ST as PackST-import qualified Numeric.Container as Container-import Numeric.Container ((<\>), (<>))+import qualified Matrix.Sparse as Sparse+import qualified Matrix.Vector as Vector+import qualified Matrix.QR.Givens as QR import Data.Complex (Complex((:+))) import qualified Data.List.HT as ListHT import qualified Data.List as List+import qualified Data.Array as Array+import qualified Data.Map as Map+import Data.Array (Array) import Data.Tuple.HT (mapPair, mapSnd)-import Data.Maybe (isJust, fromMaybe)+import Data.Maybe (isJust)+import Data.Map (Map) -import Control.Applicative ((<$>))+import Control.Applicative ((<$>), (<$)) fixAtLeastOne :: a -> [Maybe a] -> [Maybe a]@@ -47,20 +49,33 @@ parallel fs = uncurry zip . mapPair fs . unzip -sparseMatrix :: Int -> Int -> [((Int, Int), Double)] -> Matrix.Matrix Double-sparseMatrix numRows numCols xs =- PackST.runSTMatrix $ do- mat <- PackST.newMatrix 0 numRows numCols- mapM_ (\((r,c), x) -> PackST.writeMatrix mat r c x) xs- return mat+type Matrix = Sparse.Matrix Int Int+type Vector = Array Int +sparseMatrix :: Int -> Int -> [((Int, Int), Float)] -> Matrix Float+sparseMatrix numRows numCols =+ Sparse.fromMap ((0,0), (numRows-1, numCols-1)) . Map.fromList++-- every column must be non-empty+sparseToColumns :: Matrix a -> [Map Int a]+sparseToColumns = Map.elems . Sparse.toColumns++sparseFromColumns :: (Int,Int) -> [Map Int a] -> Matrix a+sparseFromColumns (m0,m1) cols =+ Sparse.fromColumns ((m0,0), (m1, length cols-1)) $+ Map.fromList $ zip [0..] cols++addSparseColumn :: (Num a) => Vector a -> Map Int a -> Vector a+addSparseColumn v col = Array.accum (+) v $ Map.toList col++ elm :: Int -> Int -> a -> ((Int, Int), a) elm row col x = ((row, col), x) absolutePositionsFromPairDisplacements :: [(Maybe Float, Maybe Float)] -> [((Int, Int), (Float, Float))] ->- ([(Double,Double)], [(Double,Double)])+ ([(Float,Float)], [(Float,Float)]) absolutePositionsFromPairDisplacements mxys displacements = let numPics = length mxys (mxs, mys) = unzip mxys@@ -68,18 +83,14 @@ matrix = sparseMatrix (length is) numPics $ concat $ zipWith (\k (ia,ib) -> [elm k ia (-1), elm k ib 1]) [0..] is- solve ms ds =- leastSquaresSelected matrix- (map (fmap realToFrac) ms)- (Vector.fromList (map realToFrac ds))+ solve ms ds = leastSquaresSelected matrix ms (Vector.fromList ds) (pxs, achievedDxs) = solve mxs dxs (pys, achievedDys) = solve mys dys in (zip pxs pys, zip achievedDxs achievedDys) leastSquaresSelected ::- Matrix.Matrix Double -> [Maybe Double] -> Vector.Vector Double ->- ([Double], [Double])+ Matrix Float -> [Maybe Float] -> Vector Float -> ([Float], [Float]) leastSquaresSelected m mas rhs0 = let (lhsCols,rhsCols) = ListHT.unzipEithers $@@ -87,13 +98,12 @@ (\col ma -> case ma of Nothing -> Left col- Just a -> Right $ Container.scale a col)- (Matrix.toColumns m) mas- lhs = Matrix.fromColumns lhsCols- rhs = foldl Container.add (Container.scale 0 rhs0) rhsCols- sol = lhs <\> Container.sub rhs0 rhs- in if Vector.dim rhs0 == 0 then (map (fromMaybe 0) mas, []) else- (snd $+ Just a -> Right $ fmap (a*) col)+ (sparseToColumns m) mas+ lhs = sparseFromColumns (Array.bounds rhs0) lhsCols+ rhs = foldl addSparseColumn (0 <$ rhs0) rhsCols+ sol = QR.leastSquares lhs $ Vector.sub rhs0 rhs+ in (snd $ List.mapAccumL (curry $ \x -> case x of@@ -101,14 +111,11 @@ (a:as, Nothing) -> (as, a) ([], Nothing) -> error "too few elements in solution vector") (Vector.toList sol) mas,- Vector.toList $- Container.add (lhs <> sol) rhs)+ Vector.toList $ Vector.add (Sparse.mulVector lhs sol) rhs) -zeroVector, _zeroVector :: Int -> Vector.Vector Double+zeroVector :: Int -> Vector Float zeroVector n = Vector.fromList $ replicate n 0--- fails for vectors of size 0-_zeroVector n = Container.constant 0 n {- Approximate rotation from point correspondences.@@ -123,58 +130,43 @@ \0 1 y0 x0/ |c | \y1/ \s / -We try to scale dx and dy down using 'weight'.-Otherwise they are weighted much more than the rotation.-However the weight will only influence the result+Formerly, we tried to scale dx and dy down using a weight,+such that dx and dy had magnitude similar to c and s.+However the weight would only influence the result for under-constrained equation systems.-This is usually not the case.+This is usually not the case and+the solver does not support rank-deficient matrices, anyway. -} layoutFromPairDisplacements :: [(Maybe (Float, Float), (Maybe Float, Maybe Float))] -> [((Int, Int), ((Float, Float), (Float, Float)))] ->- ([((Double,Double), Complex Double)],- [(Double,Double)])+ ([((Float,Float), Complex Float)], [(Float,Float)]) layoutFromPairDisplacements mrxys correspondences = let numPics = length mrxys- weight =- let xs =- concatMap- (\(_i, ((xai,yai),(xbi,ybi))) -> [xai, yai, xbi, ybi])- correspondences- in if null xs then 1 else realToFrac $ maximum xs - minimum xs matrix = sparseMatrix (2 * length correspondences) (4*numPics) $ concat $ zipWith- (\k ((ia,ib), ((xai,yai),(xbi,ybi))) ->- let xa = realToFrac xai- xb = realToFrac xbi- ya = realToFrac yai- yb = realToFrac ybi- in elm (k+0) (4*ia+0) (-weight) :- elm (k+1) (4*ia+1) (-weight) :- elm (k+0) (4*ia+2) (-xa) :- elm (k+0) (4*ia+3) ya :- elm (k+1) (4*ia+2) (-ya) :- elm (k+1) (4*ia+3) (-xa) :- elm (k+0) (4*ib+0) weight :- elm (k+1) (4*ib+1) weight :- elm (k+0) (4*ib+2) xb :- elm (k+0) (4*ib+3) (-yb) :- elm (k+1) (4*ib+2) yb :- elm (k+1) (4*ib+3) xb :- [])+ (\k ((ia,ib), ((xa,ya),(xb,yb))) ->+ elm (k+0) (4*ia+0) (-1) :+ elm (k+1) (4*ia+1) (-1) :+ elm (k+0) (4*ia+2) (-xa) :+ elm (k+0) (4*ia+3) ya :+ elm (k+1) (4*ia+2) (-ya) :+ elm (k+1) (4*ia+3) (-xa) :+ elm (k+0) (4*ib+0) 1 :+ elm (k+1) (4*ib+1) 1 :+ elm (k+0) (4*ib+2) xb :+ elm (k+0) (4*ib+3) (-yb) :+ elm (k+1) (4*ib+2) yb :+ elm (k+1) (4*ib+3) xb :+ []) [0,2..] correspondences (solution, projection) = leastSquaresSelected matrix- (concatMap- (\(mr, (mx,my)) ->- [(/weight) . realToFrac <$> mx,- (/weight) . realToFrac <$> my,- realToFrac . fst <$> mr,- realToFrac . snd <$> mr]) $+ (concatMap (\(mr, (mx,my)) -> [mx, my, fst <$> mr, snd <$> mr]) $ mrxys) (zeroVector (2 * length correspondences))- in (map (\[dx,dy,rx,ry] -> ((weight*dx,weight*dy), rx :+ ry)) $+ in (map (\[dx,dy,rx,ry] -> ((dx,dy), rx :+ ry)) $ ListHT.sliceVertical 4 solution, map (\[x,y] -> (x,y)) $ ListHT.sliceVertical 2 projection)
src/MatchImageBorders.hs view
@@ -90,11 +90,15 @@ // map (flip (,) locBorder) (Set.toList border) ++type+ Queue i j =+ MaxPQueue Float ((IOCArray (i,j) Location, CArray (i,j) Float), (i,j))+ prepareShaping :: (Ix i, Enum i, Ix j, Enum j) => [(CArray (i,j) Bool, CArray (i,j) Float)] ->- IO ([IOCArray (i,j) Location],- MaxPQueue Float ((IOCArray (i,j) Location, CArray (i,j) Float), (i,j)))+ IO ([IOCArray (i,j) Location], Queue i j) prepareShaping maskWeightss = fmap (mapSnd PQ.unions . unzip) $ forM maskWeightss $ \(mask, weights) -> do@@ -107,11 +111,10 @@ shapeParts :: IOCArray (Int, Int) Int -> [IOCArray (Int, Int) Location] ->- MaxPQueue Float- ((IOCArray (Int, Int) Location, CArray (Int, Int) Float), (Int, Int)) ->- IO [CArray (Int, Int) Bool]+ Queue Int Int -> IO [CArray (Int, Int) Bool] shapeParts count masks =- let loop queue =+ let loop :: Queue Int Int -> IO [CArray (Int, Int) Bool]+ loop queue = case PQ.maxView queue of Nothing -> mapM (fmap (amap (/=locOutside)) . freeze) masks Just (((locs, diffs), pos@(y,x)), remQueue) -> do