patch-image 0.1.0.1 → 0.1.0.2
raw patch · 3 files changed
+264/−173 lines, 3 filesdep +accelerate-cufftdep +accelerate-fourierdep −accelerate-fftdep ~Cabaldep ~JuicyPixelsdep ~accelerate-arithmetic
Dependencies added: accelerate-cufft, accelerate-fourier
Dependencies removed: accelerate-fft
Dependency ranges changed: Cabal, JuicyPixels, accelerate-arithmetic, accelerate-utility, vector
Files
- patch-image.cabal +12/−11
- src/Accelerate.hs +246/−162
- src/Option.hs +6/−0
patch-image.cabal view
@@ -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 <haskell@henning-thielemann.de> Maintainer: Henning Thielemann <haskell@henning-thielemann.de>-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
src/Accelerate.hs view
@@ -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
src/Option.hs view
@@ -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 ->