diff --git a/Changes.md b/Changes.md new file mode 100644 --- /dev/null +++ b/Changes.md @@ -0,0 +1,23 @@ +# Change log for the `patch-image` package + +## 0.2: + + * Add new executable that is based on LLVM and `knead`. + + * Add new algorithm for assembling the image from its parts. + The algorithm finds exactly matching part shapes, + such that the border of the shapes is where it hurts least visually. + +## 0.1.0.2: + + * Switch from `accelerate-fft` to `accelerate-cufft`. + +## 0.1: + + * Implement the patching algorithm using `accelerate-cuda`. + +## 0.0: + + * Tests for a weighting algorithm using `GeomAlgLib`. + The goal is to find a reasonable weighting + for mixing arbitrary overlapping polygons. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,196 @@ +This is the workflow: +Scan parts of an image that considerably overlap. +They must all be approximately oriented correctly. +The program uses the overlapping areas for reconstruction +of the layout of the parts. +If all parts are in the directory `part` +then in the best case you can simply run: + + patch-image --output=collage.jpeg part/*.jpeg + +If you get blurred areas, +you might enable an additional rotation correction: + + patch-image --finetune-rotate --output=collage.jpeg part/*.jpeg + +It follows an overview of how the program works. +It implies some things you should care about when using the program. + +The program runs three phases: + +* Orientate each image part individually + +* Find overlapping areas in the parts + and reconstruct the part positions within the big image + +* Blend the parts in order to get the big image + +The first phase orientates each part +such that horizontal structures become perfectly aligned. +Only the brightness channel of the image is analysed. +Horizontal structures can be text or the border of the image. +This also means that you should orientate the parts +horizontally, not vertically. +I also recommend not to mix horizontal and vertical scanned parts +since the horizontal and vertical resolution of your scanner +might differ slightly. +However, it should be fine to rotate the image source by 180° +and rotate it back digitally, +before feeding it to the patch-image program. + +## 1st Phase + +Options for the first phase: + +* `--maximum-absolute-angle`: + Maximum angle to test for. + +* `--number-angles`: + Number of angles minus one + to test between negative and positive maximum angle. + +* `--hint-angle`: + If the search for horizontal structures + does not yield satisfying results for an image part, + you may prepend the `--hint-angle` option with the wanted angle + to the image path. + +## 2nd Phase + +In the second phase the program looks +for overlapping parts between all pairs of images. +For every pair it computes a convolution via a Fourier transform. +Only the brightness channel of the image is analysed. + +* `--pad-size`: + Computing a convolution of two big images may exceed your graphics memory. + To this end, images are shrunk before convolution. + The pad size is the size in pixels after shrinking + that holds 2x2 shrunken image parts. + After determination of the distance between the shrunken parts + the matching is repeated on a non-reduced clip of the original image part, + in order to get precise coordinates. + +* `--minimum-overlap`: + There must be a minimum of overlap + in order to accept that the images actually overlap. + The overlap is measured as a portion of the image part size. + +* `--maximum-difference`: + The maximum allowed mean difference + within an overlapping area of two overlapping images. + If the difference is larger, + then the program assumes that the parts do not overlap. + +* `--smooth`: + It is important to eliminate a brightness offset, + that is, big black and big white areas should be handled equally. + To this end the image is smoothed + and the smoothed image is subtracted from the original one. + This option allows to specify the degree (radius) of the smoothing. + I don't think you ever need to touch this parameter. + +* `--output-overlap`: + Writes images for all pairs of image parts. + These images allow you to diagnose + where the matching algorithm failed. + It may help you to adjust the matching parameters. + In the future we might add an option to ignore problematic pairs. + +Since in the first phase every image part is oriented individually +it may happen that the part orientations don't match. +This would result in blurred areas in the final collage. +In order to correct this, +you can run phase two in an extended mode, +that also re-evaluates the part orientations. +The orientation of the composed image is then determined +by the estimated orientation of the first image. + +Options: + +* `--finetune-rotate`: + Enables the extended overlapping mode. + The option `--output-overlap` will then be ignored. + +* `--number-stamps`: + The extended mode selects many small clips in the overlapping area + and tries to match them. + We call these clips /stamps/. + This option controls the number of stamps per overlapping area minus one. + +* `--stamp-size`: + Size of a square stamp in pixels. + +## 3rd Phase + +The third phase composes a big image from the parts. +The parts are weighted such that the part boundaries cannot be seen anymore +and differences in brightness are faded into another. +The downside is that the superposition may lead to blur. + +Options: + +* `--output`: + Path of the output JPEG image with the weighted collage. + +* `--output-hard`: + Alternative output of a JPEG collage + where the image parts are simply averaged. + You will certainly see bumps in brightness + at the borders of the image parts. + This output may be mostly useful to promote the great weighting algorithm + employed by `--output`. + +* `--output-distance-map`: + The weight for every pixel is chosen according to the distance + to an image part boundary that lies within other parts. + The rationale is that the weight shall become zero + when the pixel is close to a position + that will be affected by a disruption otherwise. + This option allows to emit the distance map for every image part. + +* `--distance-gamma`: + If the distances are used for weighting as they are, + the program fades evenly between the overlapping image parts + over the entire overapping area. + This may mean that the overlapping area is blurred. + Raising the distance to a power greater than one reduces the area of blur. + The downside is that it also reduces the area for adaption + of differing brightness. + +The LLVM implementation provides an additional way to assemble the image parts. +The weighting approach tries to blend across all the overlapping area. +This can equalize differences in brightness. +The downside is that imperfectly matching image parts +lead to blurred content in the overlapping area. +An alternative algorithm tries to make the overlapping as small as possible +and additionally performs blending where it hurts least. +More precisely, parts are blended where they differ least. +However, if the brightness of the image parts differ +then the blending boundaries may become visible. + +Options: + +* `--output-shaped`: + Path of the output JPEG image with smoothly blended image parts + along curves of low image difference. + +* `--output-shaped-hard`: + Like before but the image parts are not smoothly faded. + Instead, every pixel belongs to exactly one original image part. + This is more for debugging purposes than of practical use. + +* `--output-shape`: + Emit the smoothed mask of each image part used for blending. + +* `--output-shape-hard`: + Emit the non-smoothed masks. + +* `--shape-smooth`: + Smooth radius of the image masks. + The higher, the smoother is the blending between parts. + +General options: + +* `--quality`: + JPEG quality percentage for writing the images. diff --git a/patch-image.cabal b/patch-image.cabal --- a/patch-image.cabal +++ b/patch-image.cabal @@ -1,5 +1,5 @@ Name: patch-image -Version: 0.1.0.2 +Version: 0.2 License: BSD3 License-File: LICENSE Author: Henning Thielemann @@ -16,163 +16,6 @@ It solves the problem that your scanner may be too small to capture a certain image as a whole. . - This is the workflow: - Scan parts of an image that considerably overlap. - They must all be approximately oriented correctly. - The program uses the overlapping areas for reconstruction - of the layout of the parts. - If all parts are in the directory @part@ - then in the best case you can simply run: - . - > patch-image --output=collage.jpeg part/*.jpeg - . - If you get blurred areas, - you might enable an additional rotation correction: - . - > patch-image --finetune-rotate --output=collage.jpeg part/*.jpeg - . - It follows an overview of how the program works. - It implies some things you should care about when using the program. - . - The program runs three phases: - . - * Orientate each image part individually - . - * Find overlapping areas in the parts - and reconstruct the part positions within the big image - . - * Blend the parts in order to get the big image - . - The first phase orientates each part - such that horizontal structures become perfectly aligned. - Only the brightness channel of the image is analysed. - Horizontal structures can be text or the border of the image. - This also means that you should orientate the parts - horizontally, not vertically. - I also recommend not to mix horizontal and vertical scanned parts - since the horizontal and vertical resolution of your scanner - might differ slightly. - However, it should be fine to rotate the image source by 180° - and rotate it back digitally, - before feeding it to the patch-image program. - . - Options for the first phase: - . - * @--maximum-absolute-angle@: - Maximum angle to test for. - . - * @--number-angles@: - Number of angles minus one - to test between negative and positive maximum angle. - . - * @--hint-angle@: - If the search for horizontal structures - does not yield satisfying results for an image part, - you may prepend the @--hint-angle@ option with the wanted angle - to the image path. - . - In the second phase the program looks - for overlapping parts between all pairs of images. - For every pair it computes a convolution via a Fourier transform. - Only the brightness channel of the image is analysed. - . - * @--pad-size@: - Computing a convolution of two big images may exceed your graphics memory. - To this end, images are shrunk before convolution. - The pad size is the size in pixels after shrinking - that holds 2x2 shrunken image parts. - After determination of the distance between the shrunken parts - the matching is repeated on a non-reduced clip of the original image part, - in order to get precise coordinates. - . - * @--minimum-overlap@: - There must be a minimum of overlap - in order to accept that the images actually overlap. - The overlap is measured as a portion of the image part size. - . - * @--maximum-difference@: - The maximum allowed mean difference - within an overlapping area of two overlapping images. - If the difference is larger, - then the program assumes that the parts do not overlap. - . - * @--smooth@: - It is important to eliminate a brightness offset, - that is, big black and big white areas should be handled equally. - To this end the image is smoothed - and the smoothed image is subtracted from the original one. - This option allows to specify the degree (radius) of the smoothing. - I don't think you ever need to touch this parameter. - . - * @--output-overlap@: - Writes images for all pairs of image parts. - These images allow you to diagnose - where the matching algorithm failed. - It may help you to adjust the matching parameters. - In the future we might add an option to ignore problematic pairs. - . - Since in the first phase every image part is oriented individually - it may happen that the part orientations don't match. - This would result in blurred areas in the final collage. - In order to correct this, - you can run phase two in an extended mode, - that also re-evaluates the part orientations. - The orientation of the composed image is then determined - by the estimated orientation of the first image. - . - Options: - . - * @--finetune-rotate@: - Enables the extended overlapping mode. - The option @--output-overlap@ will then be ignored. - . - * @--number-stamps@: - The extended mode selects many small clips in the overlapping area - and tries to match them. - We call these clips /stamps/. - This option controls the number of stamps per overlapping area minus one. - . - * @--stamp-size@: - Size of a square stamp in pixels. - . - The third phase composes a big image from the parts. - The parts are weighted such that the part boundaries cannot be seen anymore - and differences in brightness are faded into another. - The downside is that the superposition may lead to blur. - . - Options: - . - * @--output@: - Path of the output JPEG image with the weighted collage. - . - * @--output-hard@: - Alternative output of a JPEG collage - where the image parts are simply averaged. - You will certainly see bumps in brightness - at the borders of the image parts. - This output may be mostly useful to promote the great weighting algorithm - employed by @--output@. - . - * @--output-distance-map@: - The weight for every pixel is chosen according to the distance - to an image part boundary that lies within other parts. - The rationale is that the weight shall become zero - when the pixel is close to a position - that will be affected by a disruption otherwise. - This option allows to emit the distance map for every image part. - . - * @--distance-gamma@: - If the distances are used for weighting as they are, - the program fades evenly between the overlapping image parts - over the entire overapping area. - This may mean that the overlapping area is blurred. - Raising the distance to a power greater than one reduces the area of blur. - The downside is that it also reduces the area for adaption - of differing brightness. - . - * @--quality@: - JPEG quality percentage for writing the images. - . Restrictions: . * Only supports JPEG format. @@ -183,9 +26,12 @@ Tested-With: GHC==7.8.3 Cabal-Version: >=1.6 Build-Type: Simple +Extra-Source-Files: + Changes.md + README.md Source-Repository this - Tag: 0.1.0.2 + Tag: 0.2 Type: darcs Location: http://hub.darcs.net/thielema/patch-image/ @@ -193,13 +39,61 @@ Type: darcs Location: http://hub.darcs.net/thielema/patch-image/ +Flag llvm + Description: Build program version base on knead + Default: True + Manual: True + +Flag cuda + Description: Build program version base on accelerate-cuda + Manual: True + Flag buildDraft - description: Build draft program - default: False + Description: Build draft program + Default: False + Manual: True -Executable patch-image +Executable patch-image-llvm + Main-Is: Knead.hs + Other-Modules: + MatchImageBorders + KneadShape + LinearAlgebra + Arithmetic + Option.Utility + Option + Hs-Source-Dirs: src + + GHC-Options: -Wall -threaded -fwarn-tabs -fwarn-incomplete-record-updates + GHC-Prof-Options: -fprof-auto -rtsopts + + If flag(llvm) + Build-Depends: + knead >=0.2.1 && <0.3, + llvm-extra >=0.7 && <0.8, + llvm-tf >=3.1 && <3.2, + tfp >=1.0 && <1.1, + JuicyPixels >=2.0 && <3.3, + hmatrix >=0.15 && <0.16, + vector >=0.10 && <0.13, + pqueue >=1.2 && <1.4, + enumset >=0.0.4 && <0.1, + containers >=0.4.2 && <0.6, + fft >=0.1.7 && <0.2, + carray >=0.1.5 && <0.2, + array >=0.5 && <0.6, + Cabal >=1.18 && <2, + filepath >=1.3 && <1.4, + utility-ht >=0.0.1 && <0.1, + base >=4 && <5 + Else + Buildable: False + +Executable patch-image-cuda Main-Is: Accelerate.hs Other-Modules: + LinearAlgebra + Arithmetic Option.Utility Option Hs-Source-Dirs: src @@ -207,22 +101,25 @@ GHC-Options: -Wall -threaded -fwarn-tabs -fwarn-incomplete-record-updates GHC-Prof-Options: -fprof-auto -rtsopts - 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.16, - accelerate-io >=0.15 && <0.16, - accelerate >=0.15 && <0.16, - JuicyPixels >=2.0 && <3.3, - hmatrix >=0.15 && <0.16, - gnuplot >=0.5 && <0.6, - vector >=0.10 && <0.13, - Cabal >=1.18 && <2, - filepath >=1.3 && <1.4, - utility-ht >=0.0.1 && <0.1, - base >=4 && <5 + 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.15.1, + accelerate-io >=0.15 && <0.16, + accelerate >=0.15 && <0.16, + JuicyPixels >=2.0 && <3.3, + hmatrix >=0.15 && <0.16, + gnuplot >=0.5 && <0.6, + vector >=0.10 && <0.13, + Cabal >=1.18 && <2, + filepath >=1.3 && <1.4, + utility-ht >=0.0.1 && <0.1, + base >=4 && <5 + Else + Buildable: False Executable patch-image-draft Main-Is: Draft.hs @@ -232,7 +129,7 @@ If flag(buildDraft) Build-Depends: - JuicyPixels >=2.0 && <3.2, + JuicyPixels >=2.0 && <3.3, GeomAlgLib >=0.2 && <0.3, utility-ht >=0.0.1 && <0.1, base >=4 && <5 diff --git a/src/Accelerate.hs b/src/Accelerate.hs --- a/src/Accelerate.hs +++ b/src/Accelerate.hs @@ -3,6 +3,26 @@ import qualified Option +import qualified Arithmetic as Arith +import LinearAlgebra ( + absolutePositionsFromPairDisplacements, layoutFromPairDisplacements, + ) +import Arithmetic ( + Point2, + rotateStretchMovePoint, + rotateStretchMoveBackPoint, + boundingBoxOfRotated, + linearIp, + cubicIp, + smooth3, + projectPerp, + distance, + linearScale, + divUp, + pairFromComplex, + mapComplex, + ) + 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 @@ -24,12 +44,6 @@ (:.)((:.)), Z(Z), Any(Any), All(All), (<*), (<=*), (>=*), (==*), (&&*), (||*), (?), (!), ) -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 Graphics.Gnuplot.Advanced as GP import qualified Graphics.Gnuplot.LineSpecification as LineSpec @@ -50,17 +64,14 @@ import Text.Printf (printf) import qualified Data.List.Key as Key -import qualified Data.List.HT as ListHT -import qualified Data.List as List -import qualified Data.Bits as Bit import Control.Monad.HT (void) -import Control.Monad (liftM2, zipWithM_, when, guard) +import Control.Monad (liftM2, when) import Data.Maybe.HT (toMaybe) import Data.Maybe (catMaybes) -import Data.List.HT (removeEach, mapAdjacent, tails) +import Data.List.HT (mapAdjacent, tails) import Data.Traversable (forM) import Data.Foldable (forM_, foldMap) -import Data.Tuple.HT (mapPair, mapFst, mapSnd, fst3, thd3) +import Data.Tuple.HT (mapPair, mapFst, mapSnd, mapThd3) import Data.Word (Word8) import System.IO.Unsafe (unsafePerformIO) @@ -106,6 +117,10 @@ Pic.imageData = snd $ AIO.toVectors arr } +colorImageExtent :: Array DIM3 Word8 -> (Int, Int) +colorImageExtent pic = + case A.arrayShape pic of Z:.height:.width:._chans -> (width, height) + imageFloatFromByte :: (A.Shape sh, A.Elt a, A.IsFloating a) => Acc (Array sh Word8) -> Acc (Array sh a) @@ -152,46 +167,6 @@ floatArray = id -rotatePoint :: (Num a) => (a,a) -> (a,a) -> (a,a) -rotatePoint (c,s) (x,y) = (c*x-s*y, s*x+c*y) - -rotateStretchMovePoint :: - (Fractional a) => - (a, a) -> (a, a) -> - (a, a) -> (a, a) -rotateStretchMovePoint rot (mx,my) p = - mapPair ((mx+), (my+)) $ rotatePoint rot p - -rotateStretchMoveBackPoint :: - (Fractional a) => - (a, a) -> (a, a) -> - (a, a) -> (a, a) -rotateStretchMoveBackPoint (rx,ry) (mx,my) = - let corr = recip $ rx*rx + ry*ry - rot = (corr*rx, -corr*ry) - in \(x,y) -> rotatePoint rot (x - mx, y - my) - - -boundingBoxOfRotated :: (Num a, Ord a) => (a,a) -> (a,a) -> ((a,a), (a,a)) -boundingBoxOfRotated rot (w,h) = - let (xs,ys) = - unzip $ - rotatePoint rot (0,0) : - rotatePoint rot (w,0) : - rotatePoint rot (0,h) : - rotatePoint rot (w,h) : - [] - in ((minimum xs, maximum xs), (minimum ys, maximum ys)) - -linearIp :: (Num a) => (a,a) -> a -> a -linearIp (x0,x1) t = (1-t) * x0 + t * x1 - -cubicIp :: (Fractional a) => (a,a,a,a) -> a -> a -cubicIp (xm1, x0, x1, x2) t = - let lipm12 = linearIp (xm1,x2) t - lip01 = linearIp (x0, x1) t - in lip01 + (t*(t-1)/2) * (lipm12 + (x0+x1) - 3 * lip01) - splitFraction :: (A.Elt a, A.IsFloating a) => Exp a -> (Exp Int, Exp a) splitFraction x = let i = A.floor x @@ -252,14 +227,6 @@ let (_z :. ydst :. xdst) = unliftDim2 p in A.lift $ trans (A.fromIntegral xdst, A.fromIntegral ydst) -inBoxPlain :: - (Ord a, Num a) => - (a, a) -> - (a, a) -> - Bool -inBoxPlain (width,height) (x,y) = - 0<=x && x (Exp a, Exp a) -> @@ -279,8 +246,8 @@ replicateChannel :: (A.Slice ix, A.Shape ix, A.Elt a) => - Exp ix -> Acc (Channel Z a) -> Acc (Channel ix a) -replicateChannel = LinAlg.extrudeMatrix + Exp (ix :. Int :. Int) -> Acc (Channel Z a) -> Acc (Channel ix a) +replicateChannel = LinAlg.extrudeMatrix . A.indexTail . A.indexTail {- | @rotateStretchMove rot mov@ @@ -294,7 +261,7 @@ ExpDIM2 ix -> Acc (Channel ix a) -> (Acc (Channel Z Bool), Acc (Channel ix a)) rotateStretchMove rot mov sh arr = - let ( chansDst :. heightDst :. widthDst) = sh + let (_chansDst :. heightDst :. widthDst) = sh (_chansSrc :. heightSrc :. widthSrc) = unliftDim2 $ A.shape arr coords = rotateStretchMoveCoords rot mov (widthDst, heightDst) @@ -304,7 +271,7 @@ let (chan :. _ydst :. _xdst) = unliftDim2 ix (xsrc,ysrc) = A.unlift coord in indexFrac arr (chan :. ysrc :. xsrc)) - (replicateChannel chansDst coords)) + (replicateChannel (A.lift sh) coords)) rotateLeftTop :: @@ -349,13 +316,6 @@ in \angle arr -> rot (cos angle, sin angle) arr -{- -duplicate of Graphics.Gnuplot.Utility.linearScale --} -linearScale :: Fractional a => Int -> (a,a) -> [a] -linearScale n (x0,x1) = - map (\m -> x0 + (x1-x0) * fromIntegral m / fromIntegral n) [0..n] - analyseRotations :: [Float] -> Array DIM3 Word8 -> IO () analyseRotations angles pic = do histograms <- @@ -520,8 +480,8 @@ rot radius (cos angle, sin angle) arr -ceilingPow2Exp :: Exp Int -> Exp Int -ceilingPow2Exp n = +ceilingPow2 :: Exp Int -> Exp Int +ceilingPow2 n = A.setBit 0 $ A.ceiling $ logBase 2 (fromIntegral n :: Exp Double) pad :: @@ -569,18 +529,14 @@ 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 + 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 $ A.zipWith mulConj (forward x) (forward y) -ceilingPow2 :: Int -> Int -ceilingPow2 n = - Bit.setBit 0 $ ceiling $ logBase 2 (fromIntegral n :: Double) - removeDCOffset :: (A.Elt a, A.IsFloating a) => Acc (Channel Z a) -> Acc (Channel Z a) removeDCOffset arr = @@ -604,9 +560,6 @@ in x==*0 ||* y==*0 ? (0, arr!p) -smooth3 :: (A.Elt a, A.IsFloating a) => A.Stencil3 a -> Exp a -smooth3 (l,m,r) = (l+2*m+r)/4 - lowpass, highpass :: (A.Elt a, A.IsFloating a) => Exp Int -> Acc (Channel Z a) -> Acc (Channel Z a) @@ -668,20 +621,6 @@ attachDisplacements xsplit ysplit arr = 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 (a, (Int, Int))) -> - Acc (Channel Z (a, (Int, Int))) -weightOverlapScores minOverlap (widtha,heighta) (widthb,heightb) = - A.map - (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 ((clipWidth >=* minOverlap &&* clipHeight >=* minOverlap) - ? - (v / (A.fromIntegral clipWidth * A.fromIntegral clipHeight), 0), - dp)) {- | Set all scores to zero within a certain border. @@ -690,17 +629,18 @@ -} minimumOverlapScores :: (A.Elt a, A.IsFloating 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))) -> Acc (Channel Z (a, (Int, Int))) -minimumOverlapScores minOverlap (widtha,heighta) (widthb,heightb) = +minimumOverlapScores weight minOverlap (widtha,heighta) (widthb,heightb) = A.map - (Exp.modify (expr,(expr,expr)) $ \(v, dp@(dy,dx)) -> + (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) ? - (v, 0), + (weight (clipWidth, clipHeight) v, 0), dp)) @@ -725,15 +665,11 @@ (min widthb heightb)) weight = if False - then - weightOverlapScores minOverlap - (widtha, heighta) - (widthb, heightb) - else - minimumOverlapScores minOverlap - (widtha, heighta) - (widthb, heightb) - in weight $ + then \(clipWidth, clipHeight) v -> + v / (A.fromIntegral clipWidth * A.fromIntegral clipHeight) + else const id + in minimumOverlapScores weight minOverlap + (widtha, heighta) (widthb, heightb) $ attachDisplacements (half $ A.lift width - widthb + widtha) (half $ A.lift height - heightb + heighta) $ @@ -771,11 +707,7 @@ \(z:.yi:.xi:.yj:.xj) -> z:.yi*yk+yj:.xi*xk+xj) arr --- cf. numeric-prelude -divUp :: (Integral a) => a -> a -> a -divUp a b = - div (-a) b - type GenDIM2 a = Z :. a :. a shrinkFactors :: (Integral a) => DIM2 -> GenDIM2 a -> GenDIM2 a -> GenDIM2 a @@ -983,126 +915,6 @@ in \d a b -> Acc.the $ diff d a b --- we cannot use leastSquaresSelected here, because the right-hand side is not zero -absolutePositionsFromPairDisplacements :: - Int -> [((Int, Int), (Float, Float))] -> - ([(Double,Double)], [(Double,Double)]) -absolutePositionsFromPairDisplacements numPics displacements = - let (is, ds) = unzip displacements - (dxs, dys) = unzip ds - {- - We fix the first image to position (0,0) - in order to make the solution unique. - To this end I drop the first column from matrix. - -} - matrix = Matrix.dropColumns 1 $ PackST.runSTMatrix $ do - mat <- PackST.newMatrix 0 (length is) numPics - zipWithM_ - (\k (ia,ib) -> do - PackST.writeMatrix mat k ia (-1) - PackST.writeMatrix mat k ib 1) - [0..] is - return mat - pxs = matrix <\> Vector.fromList (map realToFrac dxs) - pys = matrix <\> Vector.fromList (map realToFrac dys) - in (zip (0 : Vector.toList pxs) (0 : Vector.toList pys), - zip (Vector.toList $ matrix <> pxs) (Vector.toList $ matrix <> pys)) - - -leastSquaresSelected :: - Matrix.Matrix Double -> [Maybe Double] -> - ([Double], [Double]) -leastSquaresSelected m mas = - let (lhsCols,rhsCols) = - ListHT.unzipEithers $ - zipWith - (\col ma -> - case ma of - Nothing -> Left col - Just a -> Right $ Container.scale a col) - (Matrix.toColumns m) mas - lhs = Matrix.fromColumns lhsCols - rhs = foldl1 Container.add rhsCols - sol = lhs <\> Container.scale (-1) rhs - in (snd $ - List.mapAccumL - (curry $ \x -> - case x of - (as, Just a) -> (as, a) - (a:as, Nothing) -> (as, a) - ([], Nothing) -> error "too few elements in solution vector") - (Vector.toList sol) mas, - Vector.toList $ - Container.add (lhs <> sol) rhs) - -{- -Approximate rotation from point correspondences. -Here (dx, dy) is the displacement with respect to the origin (0,0), -that is, the pair plays the role of the absolute position. - -x1 = dx + c*x0 - s*y0 -y1 = dy + s*x0 + c*y0 - - /dx\ -/1 0 x0 -y0\ . |dy| = /x1\ -\0 1 y0 x0/ |c | \y1/ - \s / - -Maybe, dx and dy should be scaled down. -Otherwise they are weighted much more than the rotation. --} -layoutFromPairDisplacements :: - Int -> [((Int, (Float, Float)), (Int, (Float, Float)))] -> - ([((Double,Double), HComplex.Complex Double)], - [(Double,Double)]) -layoutFromPairDisplacements numPics correspondences = - let {- - The weight will only influence the result - for under-constrained equation systems. - This is usually not the case. - -} - weight = - let xs = - concatMap - (\((_ia,(xai,yai)),(_ib,(xbi,ybi))) -> [xai, yai, xbi, ybi]) - correspondences - in realToFrac $ maximum xs - minimum xs - matrix = PackST.runSTMatrix $ do - mat <- PackST.newMatrix 0 (2 * length correspondences) (4*numPics) - zipWithM_ - (\k ((ia,(xai,yai)),(ib,(xbi,ybi))) -> do - let xa = realToFrac xai - let xb = realToFrac xbi - let ya = realToFrac yai - let yb = realToFrac ybi - PackST.writeMatrix mat (k+0) (4*ia+0) (-weight) - PackST.writeMatrix mat (k+1) (4*ia+1) (-weight) - PackST.writeMatrix mat (k+0) (4*ia+2) (-xa) - PackST.writeMatrix mat (k+0) (4*ia+3) ya - PackST.writeMatrix mat (k+1) (4*ia+2) (-ya) - PackST.writeMatrix mat (k+1) (4*ia+3) (-xa) - PackST.writeMatrix mat (k+0) (4*ib+0) weight - PackST.writeMatrix mat (k+1) (4*ib+1) weight - PackST.writeMatrix mat (k+0) (4*ib+2) xb - PackST.writeMatrix mat (k+0) (4*ib+3) (-yb) - PackST.writeMatrix mat (k+1) (4*ib+2) yb - PackST.writeMatrix mat (k+1) (4*ib+3) xb) - [0,2..] correspondences - return mat - {- - We fix the first image to position (0,0) and rotation (1,0) - in order to make the solution unique. - -} - (solution, projection) = - leastSquaresSelected matrix - (take (4*numPics) $ - map Just [0,0,1,0] ++ repeat Nothing) - in (map (\[dx,dy,rx,ry] -> ((weight*dx,weight*dy), rx HComplex.:+ ry)) $ - ListHT.sliceVertical 4 solution, - map (\[x,y] -> (x,y)) $ - ListHT.sliceVertical 2 projection) - - overlap2 :: (A.Slice ix, A.Shape ix) => (Exp Int, Exp Int) -> @@ -1121,8 +933,8 @@ ya = y + top; yb = ya-dy pa = A.lift $ chan :. ya :. xa pb = A.lift $ chan :. yb :. xb - inPicA = 0<=*xa &&* xa<*widtha &&* 0<=*ya &&* ya<*heighta - inPicB = 0<=*xb &&* xb<*widthb &&* 0<=*yb &&* yb<*heightb + inPicA = inBox (widtha,heighta) (xa,ya) + inPicB = inBox (widthb,heightb) (xb,yb) in inPicA ? (inPicB ? ((a!pa + b!pb)/2, a!pa), inPicB ? (b!pb, 0)) @@ -1141,48 +953,47 @@ overlap2 (dx, dy) (rot a, rot b) -emptyCanvas :: +emptyCountCanvas :: (A.Slice ix, A.Shape ix) => ix :. Int :. Int -> (Channel Z Int, Channel ix Float) -emptyCanvas = +emptyCountCanvas = Run.with CUDA.run1 $ \sh -> let (_ix :. height :. width) = unliftDim2 sh in (A.fill (A.lift $ Z:.height:.width) 0, A.fill sh 0) -addToCanvas :: +addToCountCanvas :: (A.Slice ix, A.Shape ix, A.Elt a, A.IsNum a) => (Acc (Channel Z Bool), Acc (Channel ix a)) -> (Acc (Channel Z Int), Acc (Channel ix a)) -> (Acc (Channel Z Int), Acc (Channel ix a)) -addToCanvas (mask, pic) (count, canvas) = +addToCountCanvas (mask, pic) (count, canvas) = (A.zipWith (+) (A.map A.boolToInt mask) count, A.zipWith (+) canvas $ A.zipWith (*) pic $ - replicateChannel - (A.indexTail $ A.indexTail $ A.shape pic) - (A.map (A.fromIntegral . A.boolToInt) mask)) + replicateChannel (A.shape pic) $ + A.map (A.fromIntegral . A.boolToInt) mask) -updateCanvas :: - (Float,Float) -> (Float,Float) -> Array DIM3 Word8 -> +updateCountCanvas :: + ((Float,Float), (Float,Float), Array DIM3 Word8) -> (Channel Z Int, Channel DIM1 Float) -> (Channel Z Int, Channel DIM1 Float) -updateCanvas = +updateCountCanvas = Run.with CUDA.run1 $ - \rot mov pic (count,canvas) -> - addToCanvas + \(rot, mov, pic) (count,canvas) -> + addToCountCanvas (rotateStretchMove rot mov (unliftDim2 $ A.shape canvas) $ separateChannels $ imageFloatFromByte pic) (count,canvas) -finalizeCanvas :: (Channel Z Int, Channel DIM1 Float) -> Array DIM3 Word8 -finalizeCanvas = +finalizeCountCanvas :: (Channel Z Int, Channel DIM1 Float) -> Array DIM3 Word8 +finalizeCountCanvas = Run.with CUDA.run1 $ \(count, canvas) -> imageByteFromFloat $ interleaveChannels $ A.zipWith (/) canvas $ - replicateChannel (A.indexTail $ A.indexTail $ A.shape canvas) $ + replicateChannel (A.shape canvas) $ A.map A.fromIntegral count @@ -1209,38 +1020,6 @@ maskedMaximum = A.fold1 (maybePlus max) - -type Line2 a = (Point2 a, Point2 a) - -intersect :: - (Ord a, Fractional a) => Line2 a -> Line2 a -> Maybe (Point2 a) -intersect ((xa,ya), (xb,yb)) ((xc,yc), (xd,yd)) = do - let denom = (xb-xa)*(yd-yc)-(xd-xc)*(yb-ya) - r = ((xd-xc)*(ya-yc)-(xa-xc)*(yd-yc)) / denom - s = ((xb-xa)*(ya-yc)-(xa-xc)*(yb-ya)) / denom - guard (denom/=0) - guard (0<=r && r<=1) - guard (0<=s && s<=1) - return (xa + r*(xb-xa), ya + r*(yb-ya)) - -intersections :: - (Fractional a, Ord a) => - [Line2 a] -> [Line2 a] -> [Point2 a] -intersections segments0 segments1 = - catMaybes $ liftM2 intersect segments0 segments1 - - -type Point2 a = (a,a) - -projectPerp :: - (Fractional a) => - Point2 a -> (Point2 a, Point2 a) -> (a, Point2 a) -projectPerp (xc,yc) ((xa,ya), (xb,yb)) = - let dx = xb-xa - dy = yb-ya - r = ((xc-xa)*dx + (yc-ya)*dy) / (dx*dx + dy*dy) - in (r, (xa + r*dx, ya + r*dy)) - project :: (A.Elt a, A.IsFloating a) => Point2 (Exp a) -> @@ -1251,11 +1030,6 @@ in (0<=*r &&* r<=*1, y) -distance :: (Floating a) => Point2 a -> Point2 a -> a -distance (xa,ya) (xb,yb) = - sqrt $ (xa-xb)^(2::Int) + (ya-yb)^(2::Int) - - distanceMapEdges :: (A.Elt a, A.IsFloating a) => Exp DIM2 -> Acc (Array DIM1 ((a,a),(a,a))) -> Acc (Channel Z a) @@ -1274,10 +1048,12 @@ Run.with CUDA.run1 $ \sh -> imageByteFromFloat . A.map (0.01*) . distanceMapEdges sh +type Geometry a = Arith.Geometry Int a + distanceMapBox :: (A.Elt a, A.IsFloating a) => Exp DIM2 -> - Exp ((a,a), (a,a), (Int,Int)) -> + Exp (Geometry a) -> 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)) = @@ -1355,7 +1131,7 @@ containedAnywhere :: (A.Elt a, A.IsFloating a) => - Acc (Array DIM1 ((a,a), (a,a), (Int,Int))) -> + Acc (Array DIM1 (Geometry a)) -> Acc (Array DIM3 (a,a)) -> Acc (Array DIM3 Bool) containedAnywhere geoms arr = @@ -1372,8 +1148,8 @@ distanceMapContained :: (A.IsFloating a, A.Elt a) => Exp DIM2 -> - Exp ((a, a), (a, a), (Int, Int)) -> - Acc (Array DIM1 ((a, a), (a, a), (Int, Int))) -> + Exp (Geometry a) -> + Acc (Array DIM1 (Geometry a)) -> Acc (Channel Z a) distanceMapContained sh this others = let distMap = @@ -1391,10 +1167,7 @@ contained distMap distanceMapContainedRun :: - DIM2 -> - ((Float,Float),(Float,Float),(Int,Int)) -> - [((Float,Float),(Float,Float),(Int,Int))] -> - Channel Z Word8 + DIM2 -> Geometry Float -> [Geometry Float] -> Channel Z Word8 distanceMapContainedRun = let distances = Run.with CUDA.run1 $ @@ -1461,8 +1234,8 @@ distanceMap :: (A.Elt a, A.IsFloating a) => Exp DIM2 -> - Exp ((a, a), (a, a), (Int, Int)) -> - Acc (Array DIM1 ((a, a), (a, a), (Int, Int))) -> + Exp (Geometry a) -> + Acc (Array DIM1 (Geometry a)) -> Acc (Array DIM1 (a, a)) -> Acc (Channel Z a) distanceMap sh this others points = @@ -1472,8 +1245,8 @@ distanceMapRun :: DIM2 -> - ((Float,Float),(Float,Float),(Int,Int)) -> - [((Float,Float),(Float,Float),(Int,Int))] -> + Geometry Float -> + [Geometry Float] -> [Point2 Float] -> Channel Z Word8 distanceMapRun = @@ -1495,8 +1268,8 @@ (A.Elt a, A.IsFloating a) => Exp a -> Exp DIM2 -> - Exp ((a, a), (a, a), (Int, Int)) -> - Acc (Array DIM1 ((a, a), (a, a), (Int, Int))) -> + Exp (Geometry a) -> + Acc (Array DIM1 (Geometry a)) -> Acc (Array DIM1 (a, a)) -> Acc (Channel Z a) distanceMapGamma gamma sh this others points = @@ -1522,14 +1295,12 @@ addToWeightedCanvas (weight, pic) (weightSum, canvas) = (A.zipWith (+) weight weightSum, A.zipWith (+) canvas $ A.zipWith (*) pic $ - replicateChannel - (A.indexTail $ A.indexTail $ A.shape pic) - weight) + replicateChannel (A.shape pic) weight) -- launch timeout updateWeightedCanvasMerged :: - ((Float,Float),(Float,Float),(Int,Int)) -> - [((Float,Float),(Float,Float),(Int,Int))] -> + Geometry Float -> + [Geometry Float] -> [Point2 Float] -> Array DIM3 Word8 -> (Channel Z Float, Channel DIM1 Float) -> @@ -1553,8 +1324,8 @@ updateWeightedCanvas :: Float -> - ((Float,Float),(Float,Float),(Int,Int)) -> - [((Float,Float),(Float,Float),(Int,Int))] -> + Geometry Float -> + [Geometry Float] -> [Point2 Float] -> Array DIM3 Word8 -> (Channel Z Float, Channel DIM1 Float) -> @@ -1580,8 +1351,8 @@ -- launch timeout updateWeightedCanvasSplit :: - ((Float,Float),(Float,Float),(Int,Int)) -> - [((Float,Float),(Float,Float),(Int,Int))] -> + Geometry Float -> + [Geometry Float] -> [Point2 Float] -> Array DIM3 Word8 -> (Channel Z Float, Channel DIM1 Float) -> @@ -1610,7 +1381,7 @@ \(weightSum, canvas) -> imageByteFromFloat $ interleaveChannels $ A.zipWith (/) canvas $ - replicateChannel (A.indexTail $ A.indexTail $ A.shape canvas) weightSum + replicateChannel (A.shape canvas) weightSum processOverlap :: @@ -1630,16 +1401,9 @@ (Nothing, optimalOverlapBigFine padExtent (Option.minimumOverlap opt)) Nothing -> - let (rotHeights, rotWidths) = - unzip $ - map (\(Z:.height:.width:._chans) -> (height, width)) $ - map (A.arrayShape . snd) picAngles - maxSum2 sizes = - case List.sortBy (flip compare) sizes of - size0 : size1 : _ -> size0+size1 - _ -> error "less than one picture - there should be no pairs" - padWidth = ceilingPow2 $ maxSum2 rotWidths - padHeight = ceilingPow2 $ maxSum2 rotHeights + let (padWidth, padHeight) = + Arith.correlationSize (Option.minimumOverlap opt) $ + map (colorImageExtent . snd) picAngles padExtent = Z :. padHeight :. padWidth in (Just $ allOverlapsRun padExtent (Option.minimumOverlap opt), optimalOverlap padExtent (Option.minimumOverlap opt)) @@ -1701,13 +1465,6 @@ return (floatPoss, picRots) -pairFromComplex :: (RealFloat a) => Complex a -> (a,a) -pairFromComplex z = (HComplex.realPart z, HComplex.imagPart z) - -mapComplex :: (a -> b) -> Complex a -> Complex b -mapComplex f (r HComplex.:+ i) = f r HComplex.:+ f i - - processOverlapRotate :: Option.Args -> [(Float, Array DIM3 Word8)] -> @@ -1838,79 +1595,28 @@ args (map snd picAngles) pairs notice "\ncompose all parts" - let bbox (rot, pic) = - case A.arrayShape pic of - Z:.height:.width:._chans -> - boundingBoxOfRotated rot - (fromIntegral width, fromIntegral height) - ((canvasLeft,canvasRight), (canvasTop,canvasBottom)) = - mapPair - (mapPair (minimum, maximum) . unzip, - mapPair (minimum, maximum) . unzip) $ - unzip $ - zipWith - (\(mx,my) -> - mapPair (mapPair ((mx+), (mx+)), mapPair ((my+), (my+))) . bbox) - floatPoss picRots - canvasWidth = ceiling (canvasRight-canvasLeft) - canvasHeight = ceiling (canvasBottom-canvasTop) - canvasShape = Z :. canvasHeight :. canvasWidth - movRotPics = - zipWith - (\(mx,my) (rot, pic) -> ((mx-canvasLeft, my-canvasTop), rot, pic)) - floatPoss picRots - info $ - printf "canvas %f - %f, %f - %f\n" - canvasLeft canvasRight canvasTop canvasBottom - info $ printf "canvas size %d, %d\n" canvasWidth canvasHeight + let ((canvasWidth, canvasHeight), rotMovPics, canvasMsgs) = + Arith.canvasShape colorImageExtent floatPoss picRots + mapM_ info canvasMsgs + forM_ (Option.outputHard opt) $ \path -> writeImage (Option.quality opt) path $ - finalizeCanvas $ + finalizeCountCanvas $ foldl - (\canvas (mov, rot, pic) -> updateCanvas rot mov pic canvas) - (emptyCanvas (Z :. 3 :. canvasHeight :. canvasWidth)) - movRotPics + (flip updateCountCanvas) + (emptyCountCanvas (Z :. 3 :. canvasHeight :. canvasWidth)) + rotMovPics notice "\ndistance maps" - let geometries = - map - (\(mov, rot, pic) -> - let Z:.height:.width:._chans = A.arrayShape pic - trans = rotateStretchMovePoint rot mov - widthf = fromIntegral width - heightf = fromIntegral height - corner00 = trans (0,0) - corner10 = trans (widthf,0) - corner01 = trans (0,heightf) - corner11 = trans (widthf,heightf) - corners = [corner00, corner01, corner10, corner11] - edges = - [(corner00, corner10), (corner10, corner11), - (corner11, corner01), (corner01, corner00)] - in ((rot, mov, (width,height)), corners, edges)) - movRotPics - let geometryRelations = - flip map (removeEach geometries) $ - \((thisGeom, thisCorners, thisEdges), others) -> - let intPoints = intersections thisEdges $ concatMap thd3 others - overlappingCorners = - filter - (\c -> - any (\(rot, mov, (width,height)) -> - inBoxPlain (width,height) $ - mapPair (round, round) $ - rotateStretchMoveBackPoint rot mov c) $ - map fst3 others) - thisCorners - allPoints = intPoints ++ overlappingCorners - otherGeoms = map fst3 others - in (thisGeom, otherGeoms, allPoints) + Arith.geometryRelations $ + map (Arith.geometryFeatures . mapThd3 colorImageExtent) rotMovPics forM_ (zip geometryRelations picAngles) $ \((thisGeom, otherGeoms, allPoints), (path, _)) -> do let stem = FilePath.takeBaseName path + let canvasShape = Z :. canvasHeight :. canvasWidth when False $ do writeGrey (Option.quality opt) (printf "/tmp/%s-distance-box.jpeg" stem) $ @@ -1940,4 +1646,4 @@ (zip geometryRelations picRots) main :: IO () -main = process =<< Option.get +main = process =<< Option.get Option.Accelerate diff --git a/src/Arithmetic.hs b/src/Arithmetic.hs new file mode 100644 --- /dev/null +++ b/src/Arithmetic.hs @@ -0,0 +1,312 @@ +module Arithmetic where + +import qualified Data.Complex as Complex +import Data.Complex (Complex, ) + +import Control.Monad (liftM2, guard) + +import qualified Data.List.HT as ListHT +import qualified Data.List as List +import qualified Data.Bits as Bit +import Data.Maybe (catMaybes) +import Data.Tuple.HT (mapPair, fst3, thd3) + +import Text.Printf (PrintfArg, printf) + + +inBox :: + (Ord a, Num a) => + (a, a) -> + (a, a) -> + Bool +inBox (width,height) (x,y) = + 0<=x && x (a,a) -> (a,a) -> (a,a) +rotatePoint (c,s) (x,y) = (c*x-s*y, s*x+c*y) + +rotateStretchMovePoint :: + (Fractional a) => + (a, a) -> (a, a) -> + (a, a) -> (a, a) +rotateStretchMovePoint rot (mx,my) p = + mapPair ((mx+), (my+)) $ rotatePoint rot p + +rotateStretchMoveBackPoint :: + (Fractional a) => + (a, a) -> (a, a) -> + (a, a) -> (a, a) +rotateStretchMoveBackPoint (rx,ry) (mx,my) = + let corr = recip $ rx*rx + ry*ry + rot = (corr*rx, -corr*ry) + in \(x,y) -> rotatePoint rot (x - mx, y - my) + + +boundingBoxOfRotated :: (Num a, Ord a) => (a,a) -> (a,a) -> ((a,a), (a,a)) +boundingBoxOfRotated = boundingBoxOfRotatedGen (min,max) + +boundingBoxOfRotatedGen :: + (Num a) => (a -> a -> a, a -> a -> a) -> (a,a) -> (a,a) -> ((a,a), (a,a)) +boundingBoxOfRotatedGen (mini,maxi) rot (w,h) = + let (xs,ys) = + unzip $ + rotatePoint rot (0,0) : + rotatePoint rot (w,0) : + rotatePoint rot (0,h) : + rotatePoint rot (w,h) : + [] + in ((foldl1 mini xs, foldl1 maxi xs), (foldl1 mini ys, foldl1 maxi ys)) + +canvasShape :: + (RealFrac a, Integral i, + PrintfArg a, PrintfArg i) => + (image -> (i, i)) -> [Point2 a] -> [((a, a), image)] -> + ((i, i), [((a, a), (a, a), image)], [String]) +canvasShape extent floatPoss picRots = + let bbox (rot, pic) = + case extent pic of + (width, height) -> + boundingBoxOfRotated rot + (fromIntegral width, fromIntegral height) + ((canvasLeft,canvasRight), (canvasTop,canvasBottom)) = + mapPair + (mapPair (minimum, maximum) . unzip, + mapPair (minimum, maximum) . unzip) $ + unzip $ + zipWith + (\(mx,my) -> + mapPair (mapPair ((mx+), (mx+)), mapPair ((my+), (my+))) . bbox) + floatPoss picRots + canvasWidth = ceiling (canvasRight-canvasLeft) + canvasHeight = ceiling (canvasBottom-canvasTop) + rotMovPics = + zipWith + (\(mx,my) (rot, pic) -> (rot, (mx-canvasLeft, my-canvasTop), pic)) + floatPoss picRots + in ((canvasWidth, canvasHeight), rotMovPics, + [printf "canvas %f - %f, %f - %f\n" + canvasLeft canvasRight canvasTop canvasBottom, + printf "canvas size %d, %d\n" canvasWidth canvasHeight]) + + + +linearIp :: (Num a) => (a,a) -> a -> a +linearIp (x0,x1) t = (1-t) * x0 + t * x1 + +cubicIp :: (Fractional a) => (a,a,a,a) -> a -> a +cubicIp (xm1, x0, x1, x2) t = + let lipm12 = linearIp (xm1,x2) t + lip01 = linearIp (x0, x1) t + in lip01 + (t*(t-1)/2) * (lipm12 + (x0+x1) - 3 * lip01) + + +data Vec a v = + Vec {vecZero :: v, vecAdd :: v -> v -> v, vecScale :: a -> v -> v} + +vecScalar :: (Num a) => Vec a a +vecScalar = Vec 0 (+) (*) + + +linearIpVec :: (Num a) => Vec a v -> (v,v) -> a -> v +linearIpVec vec (x0,x1) t = + vecAdd vec (vecScale vec (1-t) x0) (vecScale vec t x1) + +cubicIpVec :: (Fractional a) => Vec a v -> (v,v,v,v) -> a -> v +cubicIpVec vec (xm1, x0, x1, x2) t = + let lipm12 = linearIpVec vec (xm1,x2) t + lip01 = linearIpVec vec (x0, x1) t + in vecAdd vec lip01 $ + vecScale vec (t*(t-1)/2) + (foldl1 (vecAdd vec) [lipm12, x0, x1, vecScale vec (-3) lip01]) + + + +smooth3 :: (Fractional a) => (a,a,a) -> a +smooth3 (l,m,r) = (l+2*m+r)/4 + + + +type Point2 a = (a,a) + +type Line2 a = (Point2 a, Point2 a) + +intersect :: + (Ord a, Fractional a) => Line2 a -> Line2 a -> Maybe (Point2 a) +intersect ((xa,ya), (xb,yb)) ((xc,yc), (xd,yd)) = do + let denom = (xb-xa)*(yd-yc)-(xd-xc)*(yb-ya) + r = ((xd-xc)*(ya-yc)-(xa-xc)*(yd-yc)) / denom + s = ((xb-xa)*(ya-yc)-(xa-xc)*(yb-ya)) / denom + guard (denom/=0) + guard (0<=r && r<=1) + guard (0<=s && s<=1) + return (xa + r*(xb-xa), ya + r*(yb-ya)) + +intersections :: + (Fractional a, Ord a) => + [Line2 a] -> [Line2 a] -> [Point2 a] +intersections segments0 segments1 = + catMaybes $ liftM2 intersect segments0 segments1 + +type Geometry i a = ((a,a), (a,a), (i,i)) + +geometryFeatures :: + (Fractional a, Integral i) => + Geometry i a -> (Geometry i a, [Point2 a], [Line2 a]) +geometryFeatures geom@(rot, mov, (width,height)) = + let trans = rotateStretchMovePoint rot mov + widthf = fromIntegral width + heightf = fromIntegral height + corner00 = trans (0,0) + corner10 = trans (widthf,0) + corner01 = trans (0,heightf) + corner11 = trans (widthf,heightf) + corners = [corner00, corner01, corner10, corner11] + edges = + [(corner00, corner10), (corner10, corner11), + (corner11, corner01), (corner01, corner00)] + in (geom, corners, edges) + +geometryRelations :: + (RealFrac a, Integral i) => + [(Geometry i a, [Point2 a], [Line2 a])] -> + [(Geometry i a, [Geometry i a], [Point2 a])] +geometryRelations geometries = + flip map (ListHT.removeEach geometries) $ + \((thisGeom, thisCorners, thisEdges), others) -> + let intPoints = intersections thisEdges $ concatMap thd3 others + overlappingCorners = + filter + (\c -> + any (\(rot, mov, (width,height)) -> + inBox (width,height) $ + mapPair (round, round) $ + rotateStretchMoveBackPoint rot mov c) $ + map fst3 others) + thisCorners + allPoints = intPoints ++ overlappingCorners + otherGeoms = map fst3 others + in (thisGeom, otherGeoms, allPoints) + + +projectPerp :: + (Fractional a) => + Point2 a -> (Point2 a, Point2 a) -> (a, Point2 a) +projectPerp (xc,yc) ((xa,ya), (xb,yb)) = + let dx = xb-xa + dy = yb-ya + r = ((xc-xa)*dx + (yc-ya)*dy) / (dx*dx + dy*dy) + in (r, (xa + r*dx, ya + r*dy)) + + +distanceSqr :: (Num a) => Point2 a -> Point2 a -> a +distanceSqr (xa,ya) (xb,yb) = (xa-xb)^(2::Int) + (ya-yb)^(2::Int) + +distance :: (Floating a) => Point2 a -> Point2 a -> a +distance a b = sqrt $ distanceSqr a b + + +{- +duplicate of Graphics.Gnuplot.Utility.linearScale +-} +linearScale :: Fractional a => Int -> (a,a) -> [a] +linearScale n (x0,x1) = + map (\m -> x0 + (x1-x0) * fromIntegral m / fromIntegral n) [0..n] + + +minimumOverlapAbsFromPortion :: (Integral i) => Float -> (i,i) -> i +minimumOverlapAbsFromPortion minOverlapPortion (width, height) = + floor $ minOverlapPortion * fromIntegral (min width height) + + +ceilingPow2 :: (Bit.Bits i, Integral i) => i -> i +ceilingPow2 n = + Bit.setBit 0 $ ceiling $ logBase 2 (fromIntegral n :: Double) + +ceilingSmooth7, ceilingSmooth7_10, ceilingSmooth7_100 :: + (Bit.Bits i, Integral i) => i -> i +ceilingSmooth7 = ceilingSmooth7_100 + +{- | +Rounds to the smallest number of the form 2^k*j, with k>=0 and 1<=j<=10 +that is at least as large as @n@. +-} +ceilingSmooth7_10 n = + let maxFac = 10 + m = ceilingPow2 $ divUp n maxFac + in m * divUp n m + +-- cf. synthesizer-core:NumberTheory +divideByMaximumPower :: (Integral i) => i -> i -> i +divideByMaximumPower b = + let go n = + case divMod n b of + (q,0) -> go q + _ -> n + in go + +(^!) :: (Num a) => a -> Int -> a +(^!) = (^) + +isSmooth7NumberReduce, isSmooth7NumberDiv :: (Integral i) => i -> Bool +isSmooth7NumberReduce = + (1==) . flip (foldl (flip divideByMaximumPower)) [2,3,5,7] + +isSmooth7NumberDiv = + let multBig = 2^!6*3^!4*5^!2*7^!2 + mult = fromInteger multBig + in if toInteger mult == multBig + then \k -> mod mult k == 0 + else error "isSmooth7NumberDiv: Integer type too small" + +propIsSmooth7Number :: Bool +propIsSmooth7Number = + all + (\k -> isSmooth7NumberReduce k == isSmooth7NumberDiv k) + [1 .. (124::Integer)] + +ceilingSmooth7_100 n = + let maxFac = 100 + m = ceilingPow2 $ divUp n maxFac + in m * (head $ filter isSmooth7NumberDiv $ iterate (1+) $ divUp n m) + + + +{- +Since we require a minimum overlap of overlapping image pairs, +we can reduce the correlation image size by maxMinOverlap. +It means, that correlation wraps around +and correlation coefficients from both borders interfer, +but only in a stripe that we ignore. +maxMinWidth is the maximal width that the smaller image of each pair can have. +-} +correlationSize :: (Bit.Bits i, Integral i) => Float -> [(i, i)] -> (i, i) +correlationSize minOverlapPortion extents = + let ((maxWidthSum, maxMinWidth), (maxHeightSum, maxMinHeight)) = + mapPair (maxSum2, maxSum2) $ unzip extents + maxSum2 sizes = + case List.sortBy (flip compare) sizes of + size0 : size1 : _ -> (size0+size1, size1) + _ -> error "less than two pictures - there should be no pairs" + maxMinOverlap = + minimumOverlapAbsFromPortion + minOverlapPortion (maxMinWidth, maxMinHeight) + padWidth = ceilingSmooth7 $ maxWidthSum - maxMinOverlap + padHeight = ceilingSmooth7 $ maxHeightSum - maxMinOverlap + in (padWidth, padHeight) + + +-- cf. numeric-prelude +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 diff --git a/src/Knead.hs b/src/Knead.hs new file mode 100644 --- /dev/null +++ b/src/Knead.hs @@ -0,0 +1,1706 @@ +{-# LANGUAGE TypeFamilies #-} +module Main where + +import qualified Option + +import qualified MatchImageBorders +import qualified Arithmetic as Arith +import MatchImageBorders (arrayCFromKnead, arrayKneadFromC) +import LinearAlgebra ( + absolutePositionsFromPairDisplacements, layoutFromPairDisplacements, + ) +import KneadShape + (Size, Vec2(Vec2), Dim1, Dim2, Shape2, Index2, Ix2, + verticalVal, horizontalVal) + +import qualified Math.FFT as FFT +import Math.FFT.Base (FFTWReal) + +import qualified Data.Array.Knead.Parameterized.Render as RenderP +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.Expression as Expr +import Data.Array.Knead.Simple.Symbolic ((!)) +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) + +import qualified LLVM.Extra.ScalarOrVector as SoV +import qualified LLVM.Extra.Arithmetic as LLVMArith +import qualified LLVM.Extra.Multi.Value.Memory as MultiMem +import qualified LLVM.Extra.Multi.Value as MultiValue +import LLVM.Extra.Multi.Value (Atom, atom) + +import qualified LLVM.Core as LLVM + +import qualified Data.Complex as Complex +import Data.Complex (Complex((:+)), conjugate, realPart) + +import qualified Codec.Picture as Pic + +import qualified Data.Vector.Storable as SV +import Foreign.ForeignPtr (ForeignPtr, castForeignPtr) + +import qualified System.FilePath as FilePath +import qualified System.IO as IO + +import qualified Distribution.Simple.Utils as CmdLine +import qualified Distribution.Verbosity as Verbosity +import Distribution.Verbosity (Verbosity) +import Text.Printf (printf) + +import qualified Control.Functor.HT as FuncHT +import Control.Monad (liftM2, when, foldM, (<=<)) +import Control.Applicative (pure, (<$>), (<*>)) + +import qualified Data.List as List +import Data.Maybe.HT (toMaybe) +import Data.Maybe (catMaybes, isJust) +import Data.List.HT (tails) +import Data.Traversable (forM) +import Data.Foldable (forM_) +import Data.Ord.HT (comparing) +import Data.Tuple.HT (mapPair, mapFst, mapSnd, mapTriple, mapThd3, fst3, swap) +import Data.Word (Word8, Word32) + + + +type SmallSize = Word32 + +type Plane = Phys.Array Dim2 +type SymbPlane = Symb.Array Dim2 +type ColorImage a = Phys.Array Dim2 (YUV a) +type ColorImage8 = ColorImage Word8 + +type YUV a = (a,a,a) + +shape2 :: (Integral i) => i -> i -> Dim2 +shape2 height width = Vec2 (fromIntegral height) (fromIntegral width) + + +readImage :: Verbosity -> FilePath -> IO ColorImage8 +readImage verbosity path = do + epic <- Pic.readImage path + case epic of + Left msg -> ioError $ userError msg + Right dynpic -> + case dynpic of + Pic.ImageYCbCr8 pic -> do + let dat = Pic.imageData pic + CmdLine.info verbosity $ + printf "yuv %dx%d, size %d\n" + (Pic.imageWidth pic) + (Pic.imageHeight pic) + (SV.length dat) + return $ + Phys.Array + (shape2 (Pic.imageHeight pic) (Pic.imageWidth pic)) + (castForeignPtr $ fst $ SV.unsafeToForeignPtr0 dat) + _ -> ioError $ userError "unsupported image type" + + +vectorStorableFrom :: + (Shape.C sh, SV.Storable a) => + (ForeignPtr c -> ForeignPtr a) -> + Phys.Array sh c -> SV.Vector a +vectorStorableFrom castArray img = + SV.unsafeFromForeignPtr0 + (castArray $ Phys.buffer img) (fromIntegral $ Shape.size $ Phys.shape img) + +imageFromArray :: + (Pic.PixelBaseComponent c ~ a, SV.Storable a) => + (ForeignPtr b -> ForeignPtr a) -> Phys.Array Dim2 b -> Pic.Image c +imageFromArray castArray img = + let Vec2 height width = Phys.shape img + in Pic.Image { + Pic.imageWidth = fromIntegral width, + Pic.imageHeight = fromIntegral height, + Pic.imageData = vectorStorableFrom castArray img + } + +writeImage :: Int -> FilePath -> ColorImage8 -> IO () +writeImage quality path img = + Pic.saveJpgImage quality path $ Pic.ImageYCbCr8 $ + imageFromArray castForeignPtr img + +writeGrey :: Int -> FilePath -> Plane Word8 -> IO () +writeGrey quality path img = + Pic.saveJpgImage quality path $ Pic.ImageY8 $ imageFromArray id img + + +colorImageExtent :: ColorImage8 -> (Size, Size) +colorImageExtent pic = + case Phys.shape pic of Vec2 height width -> (width, height) + + +fromInt :: + (MultiValue.NativeInteger i ir, MultiValue.NativeFloating a ar) => + Exp i -> Exp a +fromInt = Expr.liftM MultiValue.fromIntegral + +floatFromByte :: + (MultiValue.NativeFloating a ar, + MultiValue.PseudoRing a, MultiValue.Real a, + MultiValue.RationalConstant a) => + Exp Word8 -> Exp a +floatFromByte = (* Expr.fromRational' (recip 255)) . fromInt + +byteFromFloat :: + (MultiValue.NativeFloating a ar, + MultiValue.Field a, MultiValue.Real a, + MultiValue.RationalConstant a) => + Exp a -> Exp Word8 +byteFromFloat = fastRound . (255*) . Expr.max 0 . Expr.min 1 + + +imageFloatFromByte :: + (Symb.C array, Shape.C sh, + MultiValue.NativeFloating a ar, + MultiValue.PseudoRing a, MultiValue.Real a, + MultiValue.RationalConstant a) => + array sh Word8 -> array sh a +imageFloatFromByte = Symb.map floatFromByte + +imageByteFromFloat :: + (Symb.C array, Shape.C sh, + MultiValue.NativeFloating a ar, + MultiValue.Field a, MultiValue.Real a, + MultiValue.RationalConstant a) => + array sh a -> array sh Word8 +imageByteFromFloat = Symb.map byteFromFloat + + +yuvByteFromFloat :: + (MultiValue.NativeFloating a ar, + MultiValue.Field a, MultiValue.Real a, + MultiValue.RationalConstant a) => + Exp (YUV a) -> Exp (YUV Word8) +yuvByteFromFloat = + Expr.modify (atom,atom,atom) $ + mapTriple (byteFromFloat, byteFromFloat, byteFromFloat) + +colorImageFloatFromByte :: + (Symb.C array, Shape.C sh, + MultiValue.NativeFloating a ar, + MultiValue.PseudoRing a, MultiValue.Real a, + MultiValue.RationalConstant a) => + array sh (YUV Word8) -> array sh (YUV a) +colorImageFloatFromByte = + Symb.map $ Expr.modify (atom,atom,atom) $ + mapTriple (floatFromByte, floatFromByte, floatFromByte) + +colorImageByteFromFloat :: + (Symb.C array, Shape.C sh, + MultiValue.NativeFloating a ar, + MultiValue.Field a, MultiValue.Real a, + MultiValue.RationalConstant a) => + array sh (YUV a) -> array sh (YUV Word8) +colorImageByteFromFloat = Symb.map yuvByteFromFloat + + +fastRound :: + (MultiValue.NativeInteger i ir, MultiValue.NativeFloating a ar) => + Exp a -> Exp i +fastRound = Expr.liftM MultiValue.roundToIntFast + +splitFraction :: + (MultiValue.NativeFloating a ar) => + Exp a -> (Exp Size, Exp a) +splitFraction = Expr.unzip . Expr.liftM MultiValue.splitFractionToInt + +ceilingToInt :: + (MultiValue.NativeFloating a ar) => + Exp a -> Exp Size +ceilingToInt = Expr.liftM MultiValue.ceilingToInt + + + + +atomDim2 :: Shape2 (Atom i) +atomDim2 = Vec2 atom atom + +atomIx2 :: Index2 (Atom i) +atomIx2 = Vec2 atom atom + +dim2 :: Exp i -> Exp i -> Exp (Shape2 i) +dim2 y x = Expr.compose (Vec2 y x) + +ix2 :: Exp i -> Exp i -> Exp (Index2 i) +ix2 y x = Expr.compose (Vec2 y x) + + +fromSize2 :: + (MultiValue.NativeFloating a ar) => + (Exp Size, Exp Size) -> (Exp a, Exp a) +fromSize2 (x,y) = (fromInt x, fromInt y) + + + +indexLimit :: SymbPlane a -> Index2 (Exp Size) -> Exp a +indexLimit img (Vec2 y x) = + let (Vec2 height width) = Expr.decompose atomDim2 $ Symb.shape img + xc = Expr.max 0 $ Expr.min (width -1) x + yc = Expr.max 0 $ Expr.min (height-1) y + in img ! ix2 yc xc + +limitIndices :: + (Symb.C array, Shape.C sh) => + Exp Dim2 -> array sh Ix2 -> array sh Ix2 +limitIndices sh = + Symb.map + (case Expr.decompose atomDim2 sh of + (Vec2 height width) -> + Expr.modify atomIx2 $ + \(Vec2 y x) -> + let xc = Expr.max 0 $ Expr.min (width -1) x + yc = Expr.max 0 $ Expr.min (height-1) y + in Vec2 yc xc) + +shiftIndicesHoriz, shiftIndicesVert :: + (Symb.C array, Shape.C sh) => + Exp Size -> array sh Ix2 -> array sh Ix2 +shiftIndicesHoriz dx = + Symb.map $ Expr.modify atomIx2 $ \(Vec2 y x) -> Vec2 y (x+dx) +shiftIndicesVert dy = + Symb.map $ Expr.modify atomIx2 $ \(Vec2 y x) -> Vec2 (y+dy) x + + +type VecExp a v = Arith.Vec (Exp a) (Exp v) + +vecYUV :: (MultiValue.PseudoRing a) => VecExp a (YUV a) +vecYUV = + Arith.Vec { + Arith.vecZero = Expr.compose (Expr.zero, Expr.zero, Expr.zero), + Arith.vecAdd = + Expr.modify2 (atom,atom,atom) (atom,atom,atom) $ + \(ay,au,av) (by,bu,bv) -> + (Expr.add ay by, Expr.add au bu, Expr.add av bv), + Arith.vecScale = + Expr.modify2 atom (atom,atom,atom) $ + \a (by,bu,bv) -> (Expr.mul a by, Expr.mul a bu, Expr.mul a bv) + } + +{- +Generated code becomes too big for LLVM here. We need sharing! +-} +indexFrac :: + (MultiValue.NativeFloating a ar, + MultiValue.Real a, MultiValue.Field a, + MultiValue.RationalConstant a) => + VecExp a v -> SymbPlane v -> Index2 (Exp a) -> Exp v +indexFrac vec img (Vec2 y x) = + let (xi,xf) = splitFraction x + (yi,yf) = splitFraction y + interpolRow yc = + Arith.cubicIpVec vec + (indexLimit img (Vec2 yc (xi-1)), + indexLimit img (Vec2 yc (xi )), + indexLimit img (Vec2 yc (xi+1)), + indexLimit img (Vec2 yc (xi+2))) + xf + in Arith.cubicIpVec vec + (interpolRow (yi-1), + interpolRow yi, + interpolRow (yi+1), + interpolRow (yi+2)) + yf + +indexFrac1 :: + (MultiValue.NativeFloating a ar, + MultiValue.Real a, MultiValue.Field a, + MultiValue.RationalConstant a) => + VecExp a v -> SymbPlane v -> Index2 (Exp a) -> Exp v +indexFrac1 vec img (Vec2 y x) = + let (xi,xf) = splitFraction x + (yi,yf) = splitFraction y + interpolRow yc = + Arith.linearIpVec vec + (indexLimit img (Vec2 yc (xi-1)), + indexLimit img (Vec2 yc (xi+1))) + xf + in Arith.linearIpVec vec + (interpolRow yi, + interpolRow (yi+1)) + yf + +gatherFrac, gatherFrac_ :: + (MultiValue.NativeFloating a ar, + MultiValue.Real a, MultiValue.Field a, + MultiValue.RationalConstant a, + MultiValue.C v) => + VecExp a v -> + SymbPlane v -> + SymbPlane (Index2 a) -> + SymbPlane v +gatherFrac_ vec src = + Symb.map (indexFrac vec src . Expr.decompose atomIx2) + +gatherFrac vec src poss = + let possSplit = + Symb.map + (Expr.modify atomIx2 $ \(Vec2 y x) -> + let (xi,xf) = splitFraction x + (yi,yf) = splitFraction y + in (Vec2 yf xf, Vec2 yi xi)) + poss + possFrac = Symb.map Expr.fst possSplit + possInt = Symb.map Expr.snd possSplit + gather = flip Symb.gather src . limitIndices (Symb.shape src) + interpolateHoriz possIntShifted = + Symb.zipWith (Arith.cubicIpVec vec . Expr.unzip4) + (Symb.zip4 + (gather $ shiftIndicesHoriz (-1) possIntShifted) + (gather $ shiftIndicesHoriz 0 possIntShifted) + (gather $ shiftIndicesHoriz 1 possIntShifted) + (gather $ shiftIndicesHoriz 2 possIntShifted)) + (Symb.map horizontalVal possFrac) + + in Symb.zipWith (Arith.cubicIpVec vec . Expr.unzip4) + (Symb.zip4 + (interpolateHoriz $ shiftIndicesVert (-1) possInt) + (interpolateHoriz $ shiftIndicesVert 0 possInt) + (interpolateHoriz $ shiftIndicesVert 1 possInt) + (interpolateHoriz $ shiftIndicesVert 2 possInt)) + (Symb.map verticalVal possFrac) + + +rotateStretchMoveCoords :: + (SV.Storable a, MultiMem.C a, + MultiValue.Real a, MultiValue.Field a, + MultiValue.RationalConstant a, MultiValue.NativeFloating a ar) => + Exp (a, a) -> + Exp (a, a) -> + Exp Dim2 -> + SymbPlane (a, a) +rotateStretchMoveCoords rot mov = + Symb.map + (let trans = + Arith.rotateStretchMoveBackPoint + (Expr.unzip rot) (Expr.unzip mov) + in Expr.modify atomIx2 $ \(Vec2 y x) -> trans $ fromSize2 (x,y)) + . + Symb.id + +inRange :: (MultiValue.Comparison a) => Exp a -> Exp a -> Exp Bool +inRange = + Expr.liftM2 $ \ size x -> do + lower <- MultiValue.cmp LLVM.CmpLE MultiValue.zero x + upper <- MultiValue.cmp LLVM.CmpLT x size + MultiValue.and lower upper + +inBox :: + (MultiValue.Comparison a) => + (Exp a, Exp a) -> + (Exp a, Exp a) -> + Exp Bool +inBox (width,height) (x,y) = + Expr.liftM2 MultiValue.and (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 +validCoords (width,height) = + Symb.map $ Expr.modify (atom,atom) $ \(x,y) -> + maskFromBool $ inBox (width,height) (fastRound x, fastRound y) + +{- | +@rotateStretchMove rot mov@ +first rotate and stretches the image according to 'rot' +and then moves the picture. +-} +rotateStretchMove :: + (SV.Storable a, MultiMem.C a, + MultiValue.Real a, MultiValue.Field a, + MultiValue.RationalConstant a, MultiValue.NativeFloating a ar, + MultiValue.C v) => + VecExp a v -> + Exp (a, a) -> + Exp (a, a) -> + Exp Dim2 -> + SymbPlane v -> + SymbPlane (MaskBool, v) +rotateStretchMove vec rot mov sh img = + let coords = rotateStretchMoveCoords rot mov sh + (Vec2 heightSrc widthSrc) = Expr.decompose atomDim2 $ Symb.shape img + in Symb.zip + (validCoords (widthSrc, heightSrc) coords) + (gatherFrac vec img $ + Symb.map (Expr.modify (atom,atom) $ \(x,y) -> ix2 y x) coords) + +rotate :: + (SV.Storable a, MultiMem.C a, + MultiValue.Real a, MultiValue.Field a, + MultiValue.RationalConstant a, MultiValue.NativeFloating a ar, + MultiValue.C v) => + VecExp a v -> + Exp (a, a) -> + SymbPlane v -> + SymbPlane v +rotate vec rot img = + let (Vec2 height width) = Expr.decompose atomDim2 $ Symb.shape img + ((left, right), (top, bottom)) = + Arith.boundingBoxOfRotatedGen (Expr.min, Expr.max) + (Expr.unzip rot) (fromSize2 (width, height)) + in Symb.map Expr.snd $ + rotateStretchMove vec rot (Expr.zip (-left) (-top)) + (dim2 (ceilingToInt (bottom-top)) (ceilingToInt (right-left))) + img + + +runRotate :: IO (Float -> ColorImage8 -> IO ColorImage8) +runRotate = do + rot <- + RenderP.run $ \rot -> + colorImageByteFromFloat . rotate vecYUV rot . colorImageFloatFromByte + return $ \ angle img -> rot (cos angle, sin angle) img + + +brightnessValue :: Exp (YUV a) -> Exp a +brightnessValue = Expr.modify (atom,atom,atom) fst3 + +brightnessPlane :: + (Symb.C array, Shape.C size) => + array size (YUV a) -> array size a +brightnessPlane = Symb.map brightnessValue + +rowHistogram :: + (Symb.C array, MultiValue.Additive a) => + array Dim2 (YUV a) -> array Dim1 a +rowHistogram = + Symb.fold1 Expr.add . + ShapeDep.backpermute + (Expr.modify atomDim2 $ \(Vec2 h w) -> (h,w)) + (Expr.modify (atom,atom) $ \(y,x) -> Vec2 y x) . + brightnessPlane + + +tailArr :: (Symb.C array) => array Dim1 a -> array Dim1 a +tailArr = ShapeDep.backpermute (Expr.max 0 . flip Expr.sub 1) (Expr.add 1) + +differentiate :: + (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 + + +runScoreRotation :: IO (Float -> ColorImage8 -> IO Float) +runScoreRotation = do + rot <- + RenderP.run $ \rot -> + rowHistogram . rotate vecYUV rot . colorImageFloatFromByte + score <- RenderP.run scoreHistogram + return $ \ angle img -> score =<< rot (cos angle, sin angle) img + +findOptimalRotation :: IO ([Float] -> ColorImage8 -> IO Float) +findOptimalRotation = do + scoreRotation <- runScoreRotation + return $ \angles pic -> + fmap (fst . List.maximumBy (comparing snd)) $ + forM angles $ \angle -> + (,) angle <$> scoreRotation (angle * (pi/180)) pic + + + +transpose :: SymbPlane a -> SymbPlane a +transpose = + ShapeDep.backpermute + (Expr.modify atomDim2 $ \(Vec2 height width) -> (Vec2 width height)) + (Expr.modify atomIx2 $ \(Vec2 x y) -> (Vec2 y x)) + +lowpassVert, lowpass :: + (MultiValue.Field a, MultiValue.Real a, MultiValue.RationalConstant a) => + SymbPlane a -> SymbPlane a +lowpassVert img = + let height = verticalVal $ Symb.shape img + in generate (Symb.shape img) $ Expr.modify atomIx2 $ \(Vec2 y x) -> + Arith.smooth3 + (img ! ix2 (Expr.max 0 (y-1)) x, + img ! ix2 y x, + img ! ix2 (Expr.min (height-1) (y+1)) x) + +lowpass = transpose . lowpassVert . transpose . lowpassVert + +nestM :: Monad m => Int -> (a -> m a) -> a -> m a +nestM n f x0 = foldM (\x () -> f x) x0 (replicate n ()) + +lowpassMulti :: IO (Int -> Plane Float -> IO (Plane Float)) +lowpassMulti = do + lp <- RenderP.run lowpass + return $ \n -> nestM n lp + + +highpassMulti :: IO (Int -> Plane Float -> IO (Plane Float)) +highpassMulti = do + lp <- lowpassMulti + sub <- RenderP.run $ Symb.zipWith Expr.sub . fixArray + return $ \n img -> sub img =<< lp n img + + + +-- counterpart to 'clip' +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 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 :: + (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))) + +-- 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) +untangleSpectra2d spec = + CArray.liftArray2 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] $ + amap (uncurry Arith.mulConj) $ + untangleSpectra2d $ FFT.dftN [0,1] $ + CArray.liftArray2 (:+) (padCArray 0 sh a) (padCArray 0 sh b) + +{- | +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)) $ + amap (uncurry Arith.mulConj) $ + untangleSpectra2d $ FFT.dftN [0,1] $ + CArray.liftArray2 (:+) (padCArray 0 sh a) (padCArray 0 sh 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) + + +type Id a = a -> a + +fixArray :: Id (Symb.Array sh a) +fixArray = id + +prepareOverlapMatching :: + IO (Int -> (Float, ColorImage8) -> IO ((Float, Float), Plane Float)) +prepareOverlapMatching = do + bright <- RenderP.run $ brightnessPlane . colorImageFloatFromByte . fixArray + hp <- highpassMulti + rotat <- RenderP.run $ rotate Arith.vecScalar + return $ \radius (angle, img) -> + let Vec2 height width = Phys.shape img + rot = (cos angle, sin angle) + ((left, _right), (top, _bottom)) = + Arith.boundingBoxOfRotated rot + (fromIntegral width, fromIntegral height) + in fmap ((,) (left, top)) $ + rotat rot =<< hp radius =<< bright img + + +wrap :: Exp Size -> Exp Size -> Exp Size -> Exp Size +wrap size split c = Expr.select (c<*split) c (c-size) + +displacementMap :: + Exp Size -> Exp Size -> Exp Dim2 -> SymbPlane (Size, Size) +displacementMap xsplit ysplit sh = + let Vec2 height width = Expr.decompose atomDim2 sh + in generate sh $ Expr.modify atomIx2 $ \(Vec2 y x) -> + (wrap width xsplit x, wrap height ysplit y) + +attachDisplacements :: + Exp Size -> Exp Size -> + SymbPlane a -> SymbPlane (a, (Size, Size)) +attachDisplacements xsplit ysplit img = + Symb.zip img $ displacementMap xsplit ysplit (Symb.shape img) + + +{- | +Set all scores to zero within a certain border. +Otherwise the matching algorithm will try to match strong bars at the borders +that are actually digitalization artifacts. +-} +minimumOverlapScores :: + (MultiValue.Select a, MultiValue.PseudoRing a, + MultiValue.IntegerConstant a, MultiValue.Real a) => + ((Exp Size, Exp Size) -> Exp a -> Exp a) -> + Exp Size -> (Exp Size, Exp Size) -> (Exp Size, Exp Size) -> + SymbPlane (a, (Size, Size)) -> + SymbPlane (a, (Size, Size)) +minimumOverlapScores weight minOverlap (widtha,heighta) (widthb,heightb) = + Symb.map + (Expr.modify (atom,(atom,atom)) $ \(v, dp@(dx,dy)) -> + let clipWidth = Expr.min widtha (widthb + dx) - Expr.max 0 dx + clipHeight = Expr.min heighta (heightb + dy) - Expr.max 0 dy + in (Expr.select + (clipWidth >=* minOverlap &&* clipHeight >=* minOverlap) + (weight (clipWidth, clipHeight) v) 0, + dp)) + + +allOverlapsFromCorrelation :: + Dim2 -> + Exp Float -> + Exp Dim2 -> Exp Dim2 -> SymbPlane Float -> + SymbPlane (Float, (Size, Size)) +allOverlapsFromCorrelation (Vec2 height width) minOverlapPortion = + \sha shb correlated -> + let (Vec2 heighta widtha) = Expr.decompose atomDim2 sha + (Vec2 heightb widthb) = Expr.decompose atomDim2 shb + half = flip Expr.idiv 2 + minOverlap = + fastRound $ + minOverlapPortion + * + fromInt + (Expr.min + (Expr.min widtha heighta) + (Expr.min widthb heightb)) + weight = + if False + then \(clipWidth, clipHeight) v -> + v / (fromInt clipWidth * fromInt clipHeight) + else const id + in minimumOverlapScores weight minOverlap + (widtha, heighta) (widthb, heightb) $ + attachDisplacements + (half $ Expr.fromInteger' (toInteger width) - widthb + widtha) + (half $ Expr.fromInteger' (toInteger height) - heightb + heighta) $ + correlated + + +allOverlapsRun :: + Dim2 -> IO (Float -> Plane Float -> Plane Float -> IO (Plane Word8)) +allOverlapsRun padExtent@(Vec2 height width) = do + run <- + RenderP.run $ \minOverlapPortion sha shb img -> + imageByteFromFloat $ + Symb.map (0.0001*) $ + Symb.map Expr.fst $ + allOverlapsFromCorrelation padExtent minOverlapPortion sha shb img + + return $ \overlap a b -> + run overlap (Phys.shape a) (Phys.shape b) + =<< liftCArray2 (correlatePaddedCArray $ mapPairInt (height, width)) a b + + +argmax :: + (MultiValue.Comparison a, MultiValue.Select a, MultiValue.Select b) => + Exp (a, b) -> Exp (a, b) -> Exp (a, b) +argmax x y = Expr.select (Expr.fst x <=* Expr.fst y) y x + +argmaximum :: + (Shape.C sh, + MultiValue.Comparison a, MultiValue.Select a, MultiValue.Select b) => + Symb.Array sh (a, b) -> Exp (a, b) +argmaximum = fold1All argmax + +optimalOverlap :: + Dim2 -> IO (Float -> Plane Float -> Plane Float -> IO (Float, (Size, Size))) +optimalOverlap padExtent@(Vec2 height width) = do + run <- + RenderP.run $ \minOverlapPortion (sha, shb) img -> + argmaximum $ + allOverlapsFromCorrelation padExtent minOverlapPortion sha shb img + + return $ \overlap a b -> + run overlap (Phys.shape a, Phys.shape b) + =<< liftCArray2 (correlatePaddedCArray $ mapPairInt (height, width)) a b + + +shrink :: + (MultiValue.Field a, MultiValue.RationalConstant a, MultiValue.Real a, + MultiValue.NativeFloating a ar) => + Shape2 (Exp Size) -> SymbPlane a -> SymbPlane a +shrink (Vec2 yk xk) = + Symb.map (/ (fromInt xk * fromInt yk)) . + Symb.fold1 Expr.add . + ShapeDep.backpermute + (Expr.modify atomDim2 $ \(Vec2 height width) -> + (Vec2 (Expr.idiv height yk) (Expr.idiv width xk), Vec2 yk xk)) + (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) + (Vec2 heighta widtha) (Vec2 heightb widthb) = + Vec2 + (Arith.divUp (heighta+heightb) $ fromIntegral heightPad) + (Arith.divUp (widtha +widthb) $ fromIntegral widthPad) + + +optimalOverlapBig :: + Dim2 -> IO (Float -> Plane Float -> Plane Float -> IO (Float, (Size, Size))) +optimalOverlapBig padExtent = do + shrnk <- RenderP.run $ shrink . Expr.decompose atomDim2 + optOverlap <- optimalOverlap padExtent + return $ \minimumOverlap a b -> do + let factors@(Vec2 yk xk) = + shrinkFactors padExtent (Phys.shape a) (Phys.shape b) + aSmall <- shrnk factors a + bSmall <- shrnk factors b + mapSnd (mapPair ((*xk), (*yk))) <$> + optOverlap minimumOverlap aSmall bSmall + + +clip :: + (MultiValue.C a) => + (Exp Size, Exp Size) -> + (Exp Size, Exp Size) -> + SymbPlane a -> SymbPlane a +clip (left,top) (width,height) = + Symb.backpermute + (Expr.compose $ Vec2 height width) + (Expr.modify (Vec2 atom atom) $ \(Vec2 y x) -> Vec2 (y+top) (x+left)) + + +overlappingArea :: + (Ord a, Num a) => + Shape2 a -> + Shape2 a -> + (a, a) -> ((a, a), (a, a), (a, a)) +overlappingArea (Vec2 heighta widtha) (Vec2 heightb widthb) (dx, dy) = + let left = max 0 dx + top = max 0 dy + right = min widtha (widthb + dx) + bottom = min heighta (heightb + dy) + width = right - left + height = bottom - top + in ((left, top), (right, bottom), (width, height)) + + +{- +Like 'optimalOverlapBig' +but computes precise distance in a second step +using a part in the overlapping area. +-} +optimalOverlapBigFine :: + 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 + overlapBig <- optimalOverlapBig padExtent + clp <- RenderP.run clip + return $ \minimumOverlap a b -> do + let shapeA = Phys.shape a + let shapeB = Phys.shape b + coarsed@(coarsedx,coarsedy) <- snd <$> overlapBig minimumOverlap a b + let ((leftOverlap, topOverlap), _, + (widthOverlap, heightOverlap)) + = overlappingArea shapeA shapeB coarsed + widthFocus = min widthOverlap $ div widthPad 2 + heightFocus = min heightOverlap $ div heightPad 2 + extentFocus = (widthFocus,heightFocus) + leftFocus = leftOverlap + div (widthOverlap-widthFocus) 2 + topFocus = topOverlap + div (heightOverlap-heightFocus) 2 + addCoarsePos (xm,ym) = (xm+coarsedx, ym+coarsedy) + clipA <- clp (leftFocus,topFocus) extentFocus a + clipB <- clp (leftFocus-coarsedx,topFocus-coarsedy) extentFocus b + mapSnd addCoarsePos <$> overlap minimumOverlap clipA clipB + + +{- +Like 'optimalOverlapBigFine' +but computes precise distances between many point pairs in a second step +using many parts in the overlapping area. +These point correspondences +can be used to compute corrections to rotation angles. +-} +optimalOverlapBigMulti :: + Dim2 -> Dim2 -> Int -> + IO (Float -> Float -> Plane Float -> Plane Float -> + IO [(Float, (Size, Size), (Size, Size))]) +optimalOverlapBigMulti padExtent (Vec2 heightStamp widthStamp) numCorrs = do + shrnk <- RenderP.run $ shrink . Expr.decompose atomDim2 + optOverlap <- optimalOverlap padExtent + overDiff <- overlapDifferenceRun + clp <- RenderP.run clip + + optOverlapFine <- optimalOverlap $ Vec2 (2*heightStamp) (2*widthStamp) + let overlapFine minimumOverlap a b + anchorA@(leftA, topA) anchorB@(leftB, topB) extent@(width,height) = do + let addCoarsePos (score, (xm,ym)) = + let xc = div (width+xm) 2 + yc = div (height+ym) 2 + in (score, + (leftA+xc, topA+yc), + (leftB+xc-xm, topB+yc-ym)) + clipA <- clp anchorA extent a + clipB <- clp anchorB extent b + addCoarsePos <$> optOverlapFine minimumOverlap clipA clipB + + return $ \maximumDiff minimumOverlap a b -> do + let factors@(Vec2 yk xk) = + shrinkFactors padExtent (Phys.shape a) (Phys.shape b) + aSmall <- shrnk factors a + bSmall <- shrnk factors b + + shrunkd@(shrunkdx, shrunkdy) + <- snd <$> optOverlap minimumOverlap aSmall bSmall + let coarsedx = shrunkdx * xk + let coarsedy = shrunkdy * yk + let coarsed = (coarsedx,coarsedy) + + diff <- overDiff shrunkd aSmall bSmall + + let ((leftOverlap, topOverlap), + (rightOverlap, bottomOverlap), + (widthOverlap, heightOverlap)) + = overlappingArea (Phys.shape a) (Phys.shape b) coarsed + + let widthStampClip = min widthOverlap widthStamp + heightStampClip = min heightOverlap heightStamp + + (if diff < maximumDiff then id else const $ return []) $ + mapM + (\(x,y) -> + overlapFine minimumOverlap a b + (x, y) (x-coarsedx, y-coarsedy) + (widthStampClip, heightStampClip)) $ + zip + (map round $ tail $ init $ + Arith.linearScale (numCorrs+1) + (fromIntegral leftOverlap :: Double, + fromIntegral $ rightOverlap - widthStampClip)) + (map round $ tail $ init $ + Arith.linearScale (numCorrs+1) + (fromIntegral topOverlap :: Double, + fromIntegral $ bottomOverlap - heightStampClip)) + + +overlapDifference :: + (MultiValue.Algebraic a, MultiValue.RationalConstant a, + MultiValue.Real a, MultiValue.NativeFloating a ar) => + (Exp Size, Exp Size) -> + SymbPlane a -> SymbPlane a -> Exp a +overlapDifference (dx,dy) a b = + let (Vec2 heighta widtha) = Expr.decompose atomDim2 $ Symb.shape a + (Vec2 heightb widthb) = Expr.decompose atomDim2 $ Symb.shape b + leftOverlap = Expr.max 0 dx + topOverlap = Expr.max 0 dy + rightOverlap = Expr.min widtha (widthb + dx) + bottomOverlap = Expr.min heighta (heightb + dy) + widthOverlap = rightOverlap - leftOverlap + heightOverlap = bottomOverlap - topOverlap + extentOverlap = (widthOverlap,heightOverlap) + in Expr.sqrt $ + (/(fromInt widthOverlap * fromInt heightOverlap)) $ + fold1All (+) $ + Symb.map Expr.sqr $ + Symb.zipWith (-) + (clip (leftOverlap,topOverlap) extentOverlap a) + (clip (leftOverlap-dx,topOverlap-dy) extentOverlap b) + +overlapDifferenceRun :: + IO ((Size, Size) -> Plane Float -> Plane Float -> IO Float) +overlapDifferenceRun = RenderP.run overlapDifference + + +overlap2 :: + (MultiValue.Field a, MultiValue.Real a, MultiValue.RationalConstant a, + MultiValue.C v) => + VecExp a v -> + (Exp Size, Exp Size) -> + (SymbPlane v, SymbPlane v) -> SymbPlane v +overlap2 vec (dx,dy) (a,b) = + let (Vec2 heighta widtha) = Expr.decompose atomDim2 $ Symb.shape a + (Vec2 heightb widthb) = Expr.decompose atomDim2 $ Symb.shape b + left = Expr.min 0 dx; right = Expr.max widtha (widthb + dx) + top = Expr.min 0 dy; bottom = Expr.max heighta (heightb + dy) + width = right - left + height = bottom - top + in generate (dim2 height width) $ Expr.modify atomIx2 $ \(Vec2 y x) -> + let xa = x + left; xb = xa-dx + ya = y + top; yb = ya-dy + pa = ix2 ya xa + pb = ix2 yb xb + inPicA = inBox (widtha,heighta) (xa,ya) + inPicB = inBox (widthb,heightb) (xb,yb) + in Expr.ifThenElse inPicA + (Expr.ifThenElse inPicB + (Arith.vecScale vec (1/2) $ Arith.vecAdd vec (a!pa) (b!pb)) + (a!pa)) + (Expr.ifThenElse inPicB (b!pb) (Arith.vecZero vec)) + +composeOverlap :: + IO ((Size, Size) -> + ((Float, ColorImage8), (Float, ColorImage8)) -> + IO ColorImage8) +composeOverlap = do + over <- + RenderP.run $ \displacement (ra, picA) (rb, picB) -> + colorImageByteFromFloat $ + overlap2 vecYUV displacement + (rotate vecYUV ra $ colorImageFloatFromByte picA, + rotate vecYUV rb $ colorImageFloatFromByte picB) + let cis angle = (cos angle, sin angle) + return $ \displacement ((angleA,picA), (angleB,picB)) -> + over displacement (cis angleA, picA) (cis angleB, picB) + + + +emptyCountCanvas :: IO (Dim2 -> IO (Plane (Word32, YUV Float))) +emptyCountCanvas = + 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 (Word32, v) -> + SymbPlane (Word32, v) +addToCountCanvas vec = + Symb.zipWith + (Expr.modify2 (atom,atom) (atom,atom) $ \(mask, pic) (count, canvas) -> + (Expr.add (intFromBool mask) count, + Arith.vecAdd vec canvas $ + Arith.vecScale vec (fromInt $ intFromBool mask) pic)) + +updateCountCanvas :: + IO (RotatedImage -> Plane (Word32, YUV Float) -> + IO (Plane (Word32, YUV Float))) +updateCountCanvas = + RenderP.run $ \(rot, mov, pic) countCanvas -> + addToCountCanvas vecYUV + (rotateStretchMove vecYUV rot mov (Symb.shape countCanvas) $ + colorImageFloatFromByte pic) + countCanvas + +finalizeCountCanvas :: IO ((Plane (Word32, YUV Float)) -> IO ColorImage8) +finalizeCountCanvas = + RenderP.run $ + colorImageByteFromFloat . + Symb.map + (Expr.modify (atom,atom) $ \(count, pixel) -> + Arith.vecScale vecYUV (recip $ fromInt count) pixel) . + fixArray + + +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))) +diffWithCanvas = + RenderP.run $ \(rot, mov, pic) avg -> + Symb.zipWith + (Expr.modify2 (atom,atom) atom $ \(b,x) y -> + (b, diffAbs (brightnessValue x) (brightnessValue y))) + (rotateStretchMove vecYUV rot mov (Symb.shape avg) $ + colorImageFloatFromByte pic) + avg + +finalizeCountCanvasFloat :: + IO ((Plane (Word32, YUV Float)) -> IO (Plane (YUV Float))) +finalizeCountCanvasFloat = + RenderP.run $ + Symb.map + (Expr.modify (atom,atom) $ \(count, pixel) -> + Arith.vecScale vecYUV (recip $ fromInt count) pixel) + . + fixArray + +emptyCanvas :: IO (Dim2 -> IO ColorImage8) +emptyCanvas = RenderP.run $ \sh -> Symb.fill sh (Expr.zip3 0 0 0) + +addMaskedToCanvas :: + IO (RotatedImage -> + Plane MaskBool -> + 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 (yuvByteFromFloat . Expr.snd) $ + rotateStretchMove vecYUV rot mov (Symb.shape canvas) $ + colorImageFloatFromByte pic) + canvas + +updateShapedCanvas :: + IO (RotatedImage -> + Plane Float -> + Plane (Float, YUV Float) -> + IO (Plane (Float, YUV Float))) +updateShapedCanvas = + RenderP.run $ \(rot, mov, pic) shape weightCanvas -> + addToWeightedCanvas vecYUV + (Symb.zipWith + (Expr.modify2 atom (atom,atom) $ \s (b,x) -> (fromInt b * s, x)) + shape $ + rotateStretchMove vecYUV rot mov (Symb.shape weightCanvas) $ + colorImageFloatFromByte pic) + weightCanvas + + +maybePlus :: + (MultiValue.C a) => + (Exp a -> Exp a -> Exp a) -> + Exp (Bool, a) -> Exp (Bool, a) -> Exp (Bool, a) +maybePlus f x y = + let (xb,xv) = Expr.unzip x + (yb,yv) = Expr.unzip y + in Expr.ifThenElse xb + (Expr.compose (Expr.true, Expr.ifThenElse yb (f xv yv) xv)) y + +maskedMinimum :: + (Shape.C sh, Symb.C array, MultiValue.Real a) => + array (sh, SmallSize) (Bool, a) -> array sh (Bool, a) +maskedMinimum = Symb.fold1 (maybePlus Expr.min) + + +generate :: + (Shape.C sh) => + Exp sh -> (Exp (Shape.Index sh) -> Exp b) -> Symb.Array sh b +generate sh f = Symb.map f $ Symb.id sh + +type Geometry a = Arith.Geometry Size a + +distanceMapBox :: + (MultiValue.Field a, MultiValue.NativeFloating a ar, + MultiValue.Real a, MultiValue.RationalConstant a) => + Exp Dim2 -> + Exp (Geometry a) -> + SymbPlane (Bool, (((a,(a,a)), (a,(a,a))), ((a,(a,a)), (a,(a,a))))) +distanceMapBox sh geom = + let (rot, mov, extent@(width,height)) = + Expr.decompose ((atom,atom),(atom,atom),(atom,atom)) geom + widthf = fromInt width + heightf = fromInt height + back = Arith.rotateStretchMoveBackPoint rot mov + forth = Arith.rotateStretchMovePoint rot mov + in generate sh $ Expr.modify atomIx2 $ \(Vec2 y x) -> + let (xsrc,ysrc) = back $ fromSize2 (x,y) + leftDist = Expr.max 0 xsrc + rightDist = Expr.max 0 $ widthf - xsrc + topDist = Expr.max 0 ysrc + bottomDist = Expr.max 0 $ heightf - ysrc + in (inBox extent (fastRound xsrc, fastRound ysrc), + (((leftDist, forth (0,ysrc)), + (rightDist, forth (widthf,ysrc))), + ((topDist, forth (xsrc,0)), + (bottomDist, forth (xsrc,heightf))))) + +distance :: + (MultiValue.Algebraic a, MultiValue.Real a, + MultiValue.IntegerConstant a) => + Arith.Point2 (Exp a) -> Arith.Point2 (Exp a) -> Exp a +distance a b = Expr.sqrt $ Arith.distanceSqr a b + +outerProduct :: + (Shape.C sha, Shape.C shb, Symb.C array) => + (Exp a -> Exp b -> Exp c) -> + array sha a -> array shb b -> array (sha,shb) c +outerProduct = + ShapeDep.backpermute2 Expr.zip Expr.fst Expr.snd + +isZero :: + (MultiValue.Comparison i, MultiValue.Integral i, + MultiValue.IntegerConstant i) => + Exp i -> Exp Bool +isZero = (==* Expr.zero) + +expEven :: + (MultiValue.Comparison i, MultiValue.Integral i, + MultiValue.IntegerConstant i) => + Exp i -> Exp Bool +expEven = isZero . flip Expr.irem 2 + +separateDistanceMap :: + (Symb.C array, Shape.C sh, MultiValue.C a) => + array sh (bool, ((a, a), (a, a))) -> + array (sh, SmallSize) (bool, a) +separateDistanceMap array = + outerProduct + (Expr.modify2 (atom, ((atom, atom), (atom, atom))) atom $ + \(b,(horiz,vert)) sel -> + (b, + Expr.ifThenElse (expEven $ Expr.idiv sel 2) + (uncurry (Expr.ifThenElse (expEven sel)) horiz) + (uncurry (Expr.ifThenElse (expEven sel)) vert))) + array (Symb.lift0 $ Symb.id 4) + + +containedAnywhere :: + (Symb.C array, Shape.C sh, + MultiValue.Field a, MultiValue.NativeFloating a ar, + MultiValue.Real a, MultiValue.RationalConstant a) => + array SmallSize (Geometry a) -> + array sh (a,a) -> + array sh Bool +containedAnywhere geoms array = + Symb.fold1 (Expr.liftM2 MultiValue.or) $ + outerProduct + (Expr.modify2 (atom,atom) ((atom,atom),(atom,atom),(atom,atom)) $ + \(xdst,ydst) (rot, mov, extent) -> + let (xsrc,ysrc) = Arith.rotateStretchMoveBackPoint rot mov (xdst,ydst) + in inBox extent (fastRound xsrc, fastRound ysrc)) + array geoms + + +distanceMapContained :: + (MultiValue.RationalConstant a, MultiValue.NativeFloating a ar, + MultiValue.PseudoRing a, MultiValue.Field a, MultiValue.Real a) => + Exp Dim2 -> + Exp (Geometry a) -> + Symb.Array SmallSize (Geometry a) -> + SymbPlane a +distanceMapContained sh this others = + let distMap = separateDistanceMap $ distanceMapBox sh this + contained = + containedAnywhere others $ + Symb.map (Expr.snd . Expr.snd) distMap + in Symb.map (Expr.modify (atom,atom) $ + \(valid, dist) -> Expr.ifThenElse valid dist 0) $ + maskedMinimum $ + Symb.zipWith + (Expr.modify2 atom (atom,(atom,atom)) $ \c (b,(dist,_)) -> + (Expr.liftM2 MultiValue.and c b, dist)) + contained distMap + + +pixelCoordinates :: + (MultiValue.NativeFloating a ar) => + Exp Dim2 -> SymbPlane (a,a) +pixelCoordinates sh = + generate sh $ Expr.modify atomIx2 $ \(Vec2 y x) -> fromSize2 (x,y) + +distanceMapPoints :: + (Shape.C sh, Symb.C array, + MultiValue.Real a, MultiValue.Algebraic a, MultiValue.IntegerConstant a) => + array sh (a,a) -> + array SmallSize (a,a) -> + array sh a +distanceMapPoints a b = + Symb.fold1 Expr.min $ + outerProduct (Expr.modify2 (atom,atom) (atom,atom) distance) a b + + +{- | +For every pixel +it computes the distance to the closest point on the image part boundary +which lies in any other image. +The rationale is that we want to fade an image out, +wherever is another image that can take over. +Such a closest point can either be a perpendicular point +at one of the image edges, +or it can be an image corner +or an intersection between this image border and another image border. +The first kind of points is computed by 'distanceMapContained' +and the second kind by 'distanceMapPoints'. +We simply compute the distances to all special points +and chose the minimal distance. +-} +distanceMap :: + (MultiValue.Algebraic a, MultiValue.Real a, + MultiValue.RationalConstant a, + MultiValue.NativeFloating a ar) => + Exp Dim2 -> + Exp (Geometry a) -> + Symb.Array SmallSize (Geometry a) -> + Symb.Array SmallSize (a, a) -> + SymbPlane a +distanceMap sh this others points = + Symb.zipWith Expr.min + (distanceMapContained sh this others) + (distanceMapPoints (pixelCoordinates sh) points) + + +pow :: + (MultiValue.Repr LLVM.Value a ~ LLVM.Value ar, + LLVM.IsFloating ar, SoV.TranscendentalConstant ar) => + Exp a -> Exp a -> Exp a +pow = + flip $ Expr.liftM2 $ \(MultiValue.Cons x) (MultiValue.Cons y) -> + fmap MultiValue.Cons $ LLVMArith.pow x y + +distanceMapGamma :: + (MultiValue.Algebraic a, MultiValue.Real a, + MultiValue.RationalConstant a, + MultiValue.NativeFloating a ar, + SoV.TranscendentalConstant ar) => + Exp a -> + Exp Dim2 -> + Exp (Geometry a) -> + Symb.Array SmallSize (Geometry a) -> + Symb.Array SmallSize (a, a) -> + SymbPlane a +distanceMapGamma gamma sh this others points = + Symb.map (pow gamma) $ distanceMap sh this others points + + +emptyWeightedCanvas :: IO (Dim2 -> IO (Plane (Float, YUV Float))) +emptyWeightedCanvas = + RenderP.run $ \sh -> Symb.fill sh (Expr.zip 0 $ Expr.zip3 0 0 0) + +addToWeightedCanvas :: + (MultiValue.PseudoRing a, MultiValue.NativeFloating a ar) => + VecExp a v -> + SymbPlane (a, v) -> + SymbPlane (a, v) -> + SymbPlane (a, v) +addToWeightedCanvas vec = + Symb.zipWith + (Expr.modify2 (atom,atom) (atom,atom) $ + \(weight, pic) (weightSum, canvas) -> + (Expr.add weight weightSum, + Arith.vecAdd vec canvas $ Arith.vecScale vec weight pic)) + +updateWeightedCanvas :: + IO (Float -> + Geometry Float -> + [Geometry Float] -> + [Arith.Point2 Float] -> + ColorImage8 -> + Plane (Float, YUV Float) -> + IO (Plane (Float, YUV Float))) +updateWeightedCanvas = do + distances <- RenderP.run distanceMapGamma + + update <- + RenderP.run $ \this pic dist weightSumCanvas -> + let (rot, mov, _) = Expr.unzip3 this + in addToWeightedCanvas vecYUV + (Symb.zip dist $ + Symb.map Expr.snd $ + rotateStretchMove vecYUV rot mov + (Symb.shape weightSumCanvas) $ + colorImageFloatFromByte pic) + weightSumCanvas + + return $ \gamma this others points pic weightSumCanvas -> do + othersVec <- Phys.vectorFromList others + pointsVec <- Phys.vectorFromList points + dists <- + distances + gamma (Phys.shape weightSumCanvas) this + othersVec pointsVec + update this pic dists weightSumCanvas + + +finalizeWeightedCanvas :: IO ((Plane (Float, YUV Float)) -> IO ColorImage8) +finalizeWeightedCanvas = + RenderP.run $ + colorImageByteFromFloat . + Symb.map + (Expr.modify (atom,atom) $ \(weightSum, pixel) -> + Arith.vecScale vecYUV (recip weightSum) pixel) . + fixArray + + + +processOverlap :: + Option.Args -> + [(Float, ColorImage8)] -> + [((Int, (FilePath, ((Float, Float), Plane Float))), + (Int, (FilePath, ((Float, Float), Plane Float))))] -> + IO ([(Float, Float)], [((Float, Float), ColorImage8)]) +processOverlap args picAngles pairs = do + let opt = Option.option args + let info = CmdLine.info (Option.verbosity opt) + + let padSize = fromIntegral $ Option.padSize opt + (maybeAllOverlapsShared, optimalOverlapShared) <- + case Just $ Vec2 padSize padSize of + Just padExtent -> do + overlap <- optimalOverlapBigFine padExtent + return (Nothing, overlap (Option.minimumOverlap opt)) + Nothing -> do + let padExtent = + uncurry Vec2 $ swap $ + Arith.correlationSize (Option.minimumOverlap opt) $ + map (colorImageExtent . snd) picAngles + overlap <- optimalOverlap padExtent + allOverlapsIO <- allOverlapsRun padExtent + return + (Just $ allOverlapsIO (Option.minimumOverlap opt), + overlap (Option.minimumOverlap opt)) + + composeOver <- composeOverlap + overlapDiff <- overlapDifferenceRun + displacements <- + fmap catMaybes $ + forM pairs $ \((ia,(pathA,(leftTopA,picA))), (ib,(pathB,(leftTopB,picB)))) -> do + forM_ maybeAllOverlapsShared $ \allOverlapsShared -> when True $ + writeGrey (Option.quality opt) + (printf "/tmp/%s-%s-score.jpeg" + (FilePath.takeBaseName pathA) (FilePath.takeBaseName pathB)) + =<< allOverlapsShared picA picB + + doffset@(dox,doy) <- snd <$> optimalOverlapShared picA picB + diff <- overlapDiff doffset picA picB + let overlapping = diff < Option.maximumDifference opt + let d = (fromIntegral dox + fst leftTopA - fst leftTopB, + fromIntegral doy + snd leftTopA - snd leftTopB) + info $ + printf "%s - %s, %s, difference %f%s\n" pathA pathB (show d) diff + (if overlapping then "" else " unrelated -> ignoring") + forM_ (Option.outputOverlap opt) $ \format -> + writeImage (Option.quality opt) + (printf format + (FilePath.takeBaseName pathA) (FilePath.takeBaseName pathB)) + -- ToDo: avoid (!!) + =<< composeOver doffset (picAngles!!ia, picAngles!!ib) + return $ toMaybe overlapping ((ia,ib), d) + + let (poss, dps) = + absolutePositionsFromPairDisplacements + (length picAngles) displacements + info "\nabsolute positions" + info $ unlines $ map show poss + + info "\ncompare position differences with pair displacements" + info $ unlines $ + zipWith + (\(dpx,dpy) (dx,dy) -> + printf "(%f,%f) (%f,%f)" dpx dpy dx dy) + dps (map snd displacements) + let (errdx,errdy) = + mapPair (maximum,maximum) $ unzip $ + zipWith + (\(dpx,dpy) (dx,dy) -> + (abs $ dpx - realToFrac dx, abs $ dpy - realToFrac dy)) + dps (map snd displacements) + + info $ + "\n" + ++ + printf "maximum horizontal error: %f\n" errdx + ++ + printf "maximum vertical error: %f\n" errdy + + let picRots = + map (mapFst (\angle -> (cos angle, sin angle))) picAngles + floatPoss = map (mapPair (realToFrac, realToFrac)) poss + + return (floatPoss, picRots) + + +processOverlapRotate :: + Option.Args -> + [(Float, ColorImage8)] -> + [((Int, (FilePath, ((Float, Float), Plane Float))), + (Int, (FilePath, ((Float, Float), Plane Float))))] -> + IO ([(Float, Float)], [((Float, Float), ColorImage8)]) +processOverlapRotate args picAngles pairs = do + let opt = Option.option args + let info = CmdLine.info (Option.verbosity opt) + + let padSize = Option.padSize opt + let stampSize = Option.stampSize opt + optimalOverlapShared <- + optimalOverlapBigMulti + (shape2 padSize padSize) + (shape2 stampSize stampSize) + (Option.numberStamps opt) + <*> pure (Option.maximumDifference opt) + <*> pure (Option.minimumOverlap opt) + + displacements <- + fmap concat $ + forM pairs $ \((ia,(pathA,(leftTopA,picA))), (ib,(pathB,(leftTopB,picB)))) -> do + let add (x0,y0) (x1,y1) = (fromIntegral x0 + x1, fromIntegral y0 + y1) + correspondences <- + map + (\(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 $ \(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 snd correspondences + + let (posRots, dps) = + layoutFromPairDisplacements (length picAngles) displacements + info "\nabsolute positions and rotations: place, rotation (magnitude, phase)" + info $ unlines $ + map + (\(d,r) -> + printf "%s, %s (%7.5f, %6.2f)" (show d) (show r) + (Complex.magnitude r) (Complex.phase r * 180/pi)) + posRots + + info "\ncompare position differences with pair displacements" + info $ unlines $ + zipWith + (\(dpx,dpy) ((_ia,pa),(_ib,pb)) -> + printf "(%f,%f) %s ~ %s" dpx dpy (show pa) (show pb)) + dps displacements + + let picRots = + zipWith + (\(angle,pic) rot -> + (Arith.pairFromComplex $ + Complex.cis angle * Arith.mapComplex realToFrac rot, + pic)) + picAngles (map snd posRots) + floatPoss = map (mapPair (realToFrac, realToFrac) . fst) posRots + + return (floatPoss, picRots) + + +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) + let info = CmdLine.info (Option.verbosity opt) + + notice "\nfind rotation angles" + findOptRot <- findOptimalRotation + picAngles <- + forM paths $ \(imageOption, path) -> do + pic <- readImage (Option.verbosity opt) path + let maxAngle = Option.maximumAbsoluteAngle opt + let angles = + Arith.linearScale (Option.numberAngleSteps opt) + (-maxAngle, maxAngle) + angle <- + case Option.angle imageOption of + Just angle -> return angle + Nothing -> findOptRot angles pic + info $ printf "%s %f\176\n" path angle + return (path, (angle*pi/180, pic)) + + notice "\nfind relative placements" + prepOverlapMatching <- prepareOverlapMatching + rotated <- + mapM (FuncHT.mapSnd (prepOverlapMatching (Option.smooth opt))) picAngles + let pairs = do + (a:as) <- tails $ zip [0..] rotated + b <- as + return (a,b) + + (floatPoss, picRots) <- + (if Option.finetuneRotate opt + then processOverlapRotate + else processOverlap) + args (map snd picAngles) pairs + + notice "\ncompose all parts" + let ((canvasWidth, canvasHeight), rotMovPics, canvasMsgs) = + Arith.canvasShape colorImageExtent floatPoss picRots + let canvasShape = shape2 canvasHeight canvasWidth + mapM_ info canvasMsgs + + forM_ (Option.outputHard opt) $ \path -> do + emptyCanv <- emptyCountCanvas + updateCanv <- updateCountCanvas + finalizeCanv <- finalizeCountCanvas + + empty <- emptyCanv canvasShape + writeImage (Option.quality opt) path =<< finalizeCanv =<< + foldM (flip updateCanv) empty rotMovPics + + + notice "\ndistance maps" + let geometryRelations = + Arith.geometryRelations $ + map (Arith.geometryFeatures . mapThd3 colorImageExtent) rotMovPics + + forM_ (Option.output opt) $ \path -> do + notice "\nweighted composition" + emptyCanv <- emptyWeightedCanvas + updateCanv <- updateWeightedCanvas + finalizeCanv <- finalizeWeightedCanvas + + empty <- emptyCanv canvasShape + writeImage (Option.quality opt) path =<< finalizeCanv =<< + foldM + (\canvas ((thisGeom, otherGeoms, allPoints), (_rot, pic)) -> + updateCanv (Option.distanceGamma opt) + thisGeom otherGeoms allPoints pic canvas) + empty (zip geometryRelations picRots) + + when (isJust (Option.outputShaped opt) || isJust (Option.outputShapedHard opt)) $ do + notice "\nmatch shapes" + emptyCanv <- emptyCountCanvas + updateCanv <- updateCountCanvas + finalizeCanv <- finalizeCountCanvasFloat + + empty <- emptyCanv canvasShape + sumImg <- foldM (flip updateCanv) empty rotMovPics + + avg <- finalizeCanv sumImg + diff <- diffWithCanvas + picDiffs <- mapM (flip diff avg) rotMovPics + getSnd <- RenderP.run $ Symb.map Expr.snd . fixArray + lp <- lowpassMulti + masks <- map (amap ((0/=) . fst)) <$> mapM arrayCFromKnead picDiffs + let smoothRadius = Option.shapeSmooth opt + smoothPicDiffs <- + mapM (arrayCFromKnead <=< lp smoothRadius <=< getSnd) picDiffs + (locs, pqueue) <- + MatchImageBorders.prepareShaping $ zip masks smoothPicDiffs + counts <- thaw . amap (fromIntegral . fst) =<< arrayCFromKnead sumImg + shapes <- MatchImageBorders.shapeParts counts locs pqueue + + let names = map (FilePath.takeBaseName . fst) picAngles + forM_ (Option.outputShapedHard opt) $ \path -> do + forM_ (Option.outputShapeHard opt) $ \format -> + forM_ (zip names shapes) $ \(name,shape) -> + writeGrey (Option.quality opt) (printf format name) $ + arrayKneadFromC $ amap (\b -> if b then 255 else 0) shape + + emptyPlainCanv <- emptyCanvas + addMasked <- addMaskedToCanvas + emptyPlain <- emptyPlainCanv canvasShape + writeImage (Option.quality opt) path =<< + foldM + (\canvas (shape, rotMovPic) -> + addMasked rotMovPic + (arrayKneadFromC $ amap (fromIntegral . fromEnum) shape) + canvas) + emptyPlain (zip shapes rotMovPics) + + forM_ (Option.outputShaped opt) $ \path -> do + smoothShapes <- + mapM + (lp smoothRadius . arrayKneadFromC . + amap (fromIntegral . fromEnum)) + shapes + forM_ (Option.outputShape opt) $ \format -> do + makeByteImage <- RenderP.run $ imageByteFromFloat . fixArray + forM_ (zip names smoothShapes) $ \(name,shape) -> + writeGrey (Option.quality opt) (printf format name) + =<< makeByteImage shape + + emptyWeightedCanv <- emptyWeightedCanvas + updateWeightedCanv <- updateShapedCanvas + finalizeWeightedCanv <- finalizeWeightedCanvas + emptyWeighted <- emptyWeightedCanv canvasShape + writeImage (Option.quality opt) path =<< + finalizeWeightedCanv =<< + foldM + (\canvas (shape, rotMovPic) -> + updateWeightedCanv rotMovPic shape canvas) + emptyWeighted (zip smoothShapes rotMovPics) + + +rotateTest :: IO () +rotateTest = do + rot <- runRotate + img <- readImage Verbosity.normal "/tmp/bild/artikel0005.jpeg" + forM_ [0..11] $ \k -> do + let path = printf "/tmp/rotated/%04d.jpeg" k + putStrLn path + writeImage 100 path =<< rot (fromInteger k * pi/6) img + +scoreTest :: IO () +scoreTest = do + score <- runScoreRotation + img <- readImage Verbosity.normal "/tmp/bild/artikel0005.jpeg" + forM_ [-10..10] $ \k -> do + print =<< score (fromInteger k * 2*pi/(360*10)) img + +main :: IO () +main = process =<< Option.get Option.Knead diff --git a/src/KneadShape.hs b/src/KneadShape.hs new file mode 100644 --- /dev/null +++ b/src/KneadShape.hs @@ -0,0 +1,176 @@ +{-# 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) diff --git a/src/LinearAlgebra.hs b/src/LinearAlgebra.hs new file mode 100644 --- /dev/null +++ b/src/LinearAlgebra.hs @@ -0,0 +1,133 @@ +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 Data.Complex as HComplex + +import qualified Data.List.HT as ListHT +import qualified Data.List as List +import Control.Monad (zipWithM_) + + +-- we cannot use leastSquaresSelected here, because the right-hand side is not zero +absolutePositionsFromPairDisplacements :: + Int -> [((Int, Int), (Float, Float))] -> + ([(Double,Double)], [(Double,Double)]) +absolutePositionsFromPairDisplacements numPics displacements = + let (is, ds) = unzip displacements + (dxs, dys) = unzip ds + {- + We fix the first image to position (0,0) + in order to make the solution unique. + To this end I drop the first column from matrix. + -} + matrix = Matrix.dropColumns 1 $ PackST.runSTMatrix $ do + mat <- PackST.newMatrix 0 (length is) numPics + zipWithM_ + (\k (ia,ib) -> do + PackST.writeMatrix mat k ia (-1) + PackST.writeMatrix mat k ib 1) + [0..] is + return mat + pxs = matrix <\> Vector.fromList (map realToFrac dxs) + pys = matrix <\> Vector.fromList (map realToFrac dys) + in (zip (0 : Vector.toList pxs) (0 : Vector.toList pys), + zip (Vector.toList $ matrix <> pxs) (Vector.toList $ matrix <> pys)) + + +leastSquaresSelected :: + Matrix.Matrix Double -> [Maybe Double] -> + ([Double], [Double]) +leastSquaresSelected m mas = + let (lhsCols,rhsCols) = + ListHT.unzipEithers $ + zipWith + (\col ma -> + case ma of + Nothing -> Left col + Just a -> Right $ Container.scale a col) + (Matrix.toColumns m) mas + lhs = Matrix.fromColumns lhsCols + rhs = foldl1 Container.add rhsCols + sol = lhs <\> Container.scale (-1) rhs + in (snd $ + List.mapAccumL + (curry $ \x -> + case x of + (as, Just a) -> (as, a) + (a:as, Nothing) -> (as, a) + ([], Nothing) -> error "too few elements in solution vector") + (Vector.toList sol) mas, + Vector.toList $ + Container.add (lhs <> sol) rhs) + +{- +Approximate rotation from point correspondences. +Here (dx, dy) is the displacement with respect to the origin (0,0), +that is, the pair plays the role of the absolute position. + +x1 = dx + c*x0 - s*y0 +y1 = dy + s*x0 + c*y0 + + /dx\ +/1 0 x0 -y0\ . |dy| = /x1\ +\0 1 y0 x0/ |c | \y1/ + \s / + +Maybe, dx and dy should be scaled down. +Otherwise they are weighted much more than the rotation. +-} +layoutFromPairDisplacements :: + Int -> [((Int, (Float, Float)), (Int, (Float, Float)))] -> + ([((Double,Double), HComplex.Complex Double)], + [(Double,Double)]) +layoutFromPairDisplacements numPics correspondences = + let {- + The weight will only influence the result + for under-constrained equation systems. + This is usually not the case. + -} + weight = + let xs = + concatMap + (\((_ia,(xai,yai)),(_ib,(xbi,ybi))) -> [xai, yai, xbi, ybi]) + correspondences + in realToFrac $ maximum xs - minimum xs + matrix = PackST.runSTMatrix $ do + mat <- PackST.newMatrix 0 (2 * length correspondences) (4*numPics) + zipWithM_ + (\k ((ia,(xai,yai)),(ib,(xbi,ybi))) -> do + let xa = realToFrac xai + let xb = realToFrac xbi + let ya = realToFrac yai + let yb = realToFrac ybi + PackST.writeMatrix mat (k+0) (4*ia+0) (-weight) + PackST.writeMatrix mat (k+1) (4*ia+1) (-weight) + PackST.writeMatrix mat (k+0) (4*ia+2) (-xa) + PackST.writeMatrix mat (k+0) (4*ia+3) ya + PackST.writeMatrix mat (k+1) (4*ia+2) (-ya) + PackST.writeMatrix mat (k+1) (4*ia+3) (-xa) + PackST.writeMatrix mat (k+0) (4*ib+0) weight + PackST.writeMatrix mat (k+1) (4*ib+1) weight + PackST.writeMatrix mat (k+0) (4*ib+2) xb + PackST.writeMatrix mat (k+0) (4*ib+3) (-yb) + PackST.writeMatrix mat (k+1) (4*ib+2) yb + PackST.writeMatrix mat (k+1) (4*ib+3) xb) + [0,2..] correspondences + return mat + {- + We fix the first image to position (0,0) and rotation (1,0) + in order to make the solution unique. + -} + (solution, projection) = + leastSquaresSelected matrix + (take (4*numPics) $ + map Just [0,0,1,0] ++ repeat Nothing) + in (map (\[dx,dy,rx,ry] -> ((weight*dx,weight*dy), rx HComplex.:+ ry)) $ + ListHT.sliceVertical 4 solution, + map (\[x,y] -> (x,y)) $ + ListHT.sliceVertical 2 projection) diff --git a/src/MatchImageBorders.hs b/src/MatchImageBorders.hs new file mode 100644 --- /dev/null +++ b/src/MatchImageBorders.hs @@ -0,0 +1,135 @@ +{- | +This is an approach for stitching images at narrow bands +along lines of small image differences. +We start with rotated and placed rectangular image masks +and then start to remove pixels from the borders of the image masks +in the order of decreasing pixel value differences. +For the sake of simplicity we calculate the difference +of a pixel value to the average of all pixel values at a certain position. +We do not recalculate the average if a pixel is removed from the priority queue. +-} +module MatchImageBorders where + +import KneadShape (Vec2(Vec2), Dim2) + +import qualified Data.Array.Knead.Simple.Physical as Phys + +import qualified Data.PQueue.Prio.Max as PQ +import qualified Data.Set as Set +import Data.PQueue.Prio.Max (MaxPQueue) +import Data.Set (Set) + +import qualified Data.Array.CArray.Base as CArrayPriv +import Data.Array.IOCArray (IOCArray) +import Data.Array.MArray (readArray, writeArray, freeze, thaw) +import Data.Array.CArray (CArray) +import Data.Array.IArray + (Ix, amap, bounds, range, rangeSize, inRange, (!), (//)) + +import Foreign.Storable (Storable) + +import Data.Traversable (forM) +import Data.Tuple.HT (mapSnd) +import Data.Maybe (mapMaybe, listToMaybe) +import Data.Word (Word8) + +import Control.Monad (filterM) +import Control.Applicative ((<$>)) + + +arrayCFromKnead :: Phys.Array Dim2 a -> IO (CArray (Int,Int) a) +arrayCFromKnead (Phys.Array (Vec2 height width) fptr) = + CArrayPriv.unsafeForeignPtrToCArray fptr + ((0,0), (fromIntegral height - 1, fromIntegral width - 1)) + +arrayKneadFromC :: + (Storable a) => CArray (Int,Int) a -> Phys.Array Dim2 a +arrayKneadFromC carray = + case bounds carray of + ((ly,lx), (uy,ux)) -> + Phys.Array + (Vec2 + (fromIntegral (rangeSize (ly,uy))) + (fromIntegral (rangeSize (lx,ux)))) + (snd $ CArrayPriv.toForeignPtr carray) + + +findBorder :: (Ix i, Enum i, Ix j, Enum j) => CArray (i,j) Bool -> Set (i,j) +findBorder mask = + let ((yl,xl), (yu,xu)) = bounds mask + revRange (l,u) = [u, pred u .. l] + findLeft y = + listToMaybe $ dropWhile (\x -> not $ mask!(y,x)) $ range (xl,xu) + findRight y = + listToMaybe $ dropWhile (\x -> not $ mask!(y,x)) $ revRange (xl,xu) + findTop x = + listToMaybe $ dropWhile (\y -> not $ mask!(y,x)) $ range (yl,yu) + findBottom x = + listToMaybe $ dropWhile (\y -> not $ mask!(y,x)) $ revRange (yl,yu) + in Set.fromList $ + mapMaybe (\y -> (,) y <$> findLeft y) (range (yl,yu)) ++ + mapMaybe (\y -> (,) y <$> findRight y) (range (yl,yu)) ++ + mapMaybe (\x -> flip (,) x <$> findTop x) (range (xl,xu)) ++ + mapMaybe (\x -> flip (,) x <$> findBottom x) (range (xl,xu)) + +pqueueFromBorder :: (Ix ix) => CArray ix Float -> Set ix -> MaxPQueue Float ix +pqueueFromBorder weights = + PQ.fromList . map (\pos -> (weights!pos, pos)) . Set.toList + + +type Location = Word8 + +locOutside, locBorder, locInside :: Location +locOutside = 0 +locBorder = 1 +locInside = 2 + +prepareLocations :: (Ix ix) => CArray ix Bool -> Set ix -> CArray ix Location +prepareLocations mask border = + amap (\b -> if b then locInside else locOutside) mask + // + map (flip (,) locBorder) (Set.toList border) + +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))) +prepareShaping maskWeightss = + fmap (mapSnd PQ.unions . unzip) $ + forM maskWeightss $ \(mask, weights) -> do + let border = findBorder mask + locations <- thaw $ prepareLocations mask border + return + (locations, + fmap ((,) (locations, weights)) $ pqueueFromBorder weights border) + +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] +shapeParts count masks = + let loop queue = + case PQ.maxView queue of + Nothing -> mapM (fmap (amap (/=locOutside)) . freeze) masks + Just (((locs, diffs), pos@(y,x)), remQueue) -> do + n <- readArray count pos + if n<=1 + then loop remQueue + else do + writeArray count pos (n-1) + writeArray locs pos locOutside + envPoss <- + filterM (fmap (locInside ==) . readArray locs) $ + filter (inRange (bounds diffs)) $ + map + (\(dy,dx) -> (y+dy, x+dx)) + [(0,1), (1,0), (0,-1), (-1,0)] + mapM_ + (\envPos -> writeArray locs envPos locBorder) + envPoss + loop $ PQ.union remQueue $ PQ.fromList $ + map (\envPos -> (diffs!envPos, ((locs, diffs), envPos))) envPoss + in loop diff --git a/src/Option.hs b/src/Option.hs --- a/src/Option.hs +++ b/src/Option.hs @@ -10,6 +10,11 @@ import Control.Monad (when) +import qualified Data.EnumSet as EnumSet +import Data.Tuple.HT (mapSnd) +import Data.Monoid ((<>)) +import Data.Word (Word8) + import qualified Distribution.Verbosity as Verbosity import qualified Distribution.ReadE as ReadE import Distribution.Verbosity (Verbosity) @@ -32,8 +37,12 @@ verbosity :: Verbosity, output :: Maybe FilePath, outputHard :: Maybe FilePath, - outputOverlap :: Maybe String, -- e.g. "/tmp/%s-%s-overlap.jpeg" - outputDistanceMap :: Maybe String, -- e.g. "/tmp/%s-distance.jpeg" + outputShaped :: Maybe FilePath, + outputShapedHard :: Maybe FilePath, + outputOverlap :: Maybe String, + outputDistanceMap :: Maybe String, + outputShape :: Maybe String, + outputShapeHard :: Maybe String, quality :: Int, maximumAbsoluteAngle :: Float, numberAngleSteps :: Int, @@ -45,7 +54,8 @@ finetuneRotate :: Bool, numberStamps :: Int, stampSize :: Int, - distanceGamma :: Float + distanceGamma :: Float, + shapeSmooth :: Int } defltOption :: Option @@ -54,8 +64,12 @@ verbosity = Verbosity.verbose, output = Nothing, outputHard = Nothing, + outputShaped = Nothing, + outputShapedHard = Nothing, outputOverlap = Nothing, outputDistanceMap = Nothing, + outputShape = Nothing, + outputShapeHard = Nothing, quality = 99, maximumAbsoluteAngle = 1, numberAngleSteps = 40, @@ -67,7 +81,8 @@ finetuneRotate = False, numberStamps = 5, stampSize = 64, - distanceGamma = 2 + distanceGamma = 2, + shapeSmooth = 200 } @@ -81,15 +96,33 @@ defltImage = Image {angle = Nothing} +data Engine = Knead | Accelerate + deriving (Eq, Ord, Enum) + +type EngineSet = EnumSet.T Word8 Engine + +knead, accelerate, generic :: EngineSet +knead = EnumSet.singleton Knead +accelerate = EnumSet.singleton Accelerate +generic = knead <> accelerate + + type Description a = [Opt.OptDescr (a -> IO a)] +type EngineDescription a = [(EngineSet, Opt.OptDescr (a -> IO a))] +opt :: + EngineSet -> [Char] -> [String] -> ArgDescr a -> String -> + (EngineSet, Opt.OptDescr a) +opt engines short long argDescr help = + (engines, Opt.Option short long argDescr help) + {- Guide for common Linux/Unix command-line options: http://www.faqs.org/docs/artu/ch10s05.html -} -optionDescription :: Description a -> Description Option +optionDescription :: Description a -> EngineDescription Option optionDescription desc = - Opt.Option ['h'] ["help"] + opt generic ['h'] ["help"] (NoArg $ \ _flags -> do programName <- Env.getProgName putStrLn $ @@ -100,7 +133,7 @@ Exit.exitSuccess) "show options" : - Opt.Option ['v'] ["verbose"] + opt generic ['v'] ["verbose"] (flip ReqArg "N" $ \str flags -> do case ReadE.runReadE Verbosity.flagToVerbosity str of Right n -> return (flags{verbosity = n}) @@ -108,117 +141,144 @@ (printf "verbosity level: 0..3, default: %d" (fromEnum $ verbosity defltOption)) : - Opt.Option [] ["output"] + opt generic [] ["output"] (flip ReqArg "PATH" $ \str flags -> return $ flags{output = Just str}) ("path to generated collage") : - Opt.Option [] ["output-hard"] + opt generic [] ["output-hard"] (flip ReqArg "PATH" $ \str flags -> return $ flags{outputHard = Just str}) ("path to collage without fading") : - Opt.Option [] ["output-overlap"] + opt knead [] ["output-shaped"] + (flip ReqArg "PATH" $ \str flags -> + return $ flags{outputShaped = Just str}) + ("path to generated collage") : + + opt knead [] ["output-shaped-hard"] + (flip ReqArg "PATH" $ \str flags -> + return $ flags{outputShapedHard = Just str}) + ("path to collage without fading") : + + opt generic [] ["output-overlap"] (flip ReqArg "FORMAT" $ \str flags -> return $ flags{outputOverlap = Just str}) - ("path format for overlapped pairs") : + ("path format for overlapped pairs like '%s-%s-overlap.jpeg'") : - Opt.Option [] ["output-distance-map"] + opt generic [] ["output-distance-map"] (flip ReqArg "FORMAT" $ \str flags -> return $ flags{outputDistanceMap = Just str}) - ("path format for distance maps") : + ("path format for distance maps like '%s-distance.jpeg'") : - Opt.Option [] ["quality"] + opt knead [] ["output-shape"] + (flip ReqArg "FORMAT" $ \str flags -> + return $ flags{outputShape = Just str}) + ("path format for smooth part shape like '%s-shape-soft.jpeg'") : + + opt knead [] ["output-shape-hard"] + (flip ReqArg "FORMAT" $ \str flags -> + return $ flags{outputShapeHard = Just str}) + ("path format for hard part shape like '%s-shape-hard.jpeg'") : + + opt generic [] ["quality"] (flip ReqArg "PERCENTAGE" $ \str flags -> fmap (\x -> flags{quality = x}) $ parseNumber "compression quality" (\q -> 0<=q && q<=100) "a percentage" str) (printf "JPEG compression quality for output, default: %d" (quality defltOption)) : - Opt.Option [] ["maximum-absolute-angle"] -- "max-abs-angle" + opt generic [] ["maximum-absolute-angle"] -- "max-abs-angle" (flip ReqArg "DEGREE" $ \str flags -> fmap (\x -> flags{maximumAbsoluteAngle = x}) $ parseNumber "maximum absolute angle" (0<=) "non-negative" str) (printf "Maximum absolute angle for test rotations, default: %f" (maximumAbsoluteAngle defltOption)) : - Opt.Option [] ["number-angles"] -- "num-angles" + opt generic [] ["number-angles"] -- "num-angles" (flip ReqArg "NATURAL" $ \str flags -> fmap (\x -> flags{numberAngleSteps = x}) $ parseNumber "number of angle steps" (0<=) "non-negative" str) (printf "Number of steps for test rotations, default: %d" (numberAngleSteps defltOption)) : - Opt.Option [] ["radon"] + opt accelerate [] ["radon"] (NoArg $ \flags -> return $ flags{radonTransform = True}) (printf "Use Radon transform for estimating orientation, default: disabled") : - Opt.Option [] ["smooth"] + opt generic [] ["smooth"] (flip ReqArg "NATURAL" $ \str flags -> fmap (\x -> flags{smooth = x}) $ parseNumber "smooth radius" (0<=) "non-negative" str) (printf "Smooth radius for DC elimination, default: %d" (smooth defltOption)) : - Opt.Option [] ["pad-size"] + opt generic [] ["pad-size"] (flip ReqArg "NATURAL" $ \str flags -> fmap (\x -> flags{padSize = x}) $ parseNumber "pad size" (0<=) "non-negative" str) (printf "Pad size for matching convolution, default: %d" (padSize defltOption)) : - Opt.Option [] ["minimum-overlap"] + opt generic [] ["minimum-overlap"] (flip ReqArg "FRACTION" $ \str flags -> fmap (\x -> flags{minimumOverlap = x}) $ parseNumber "minimum overlap" (0<=) "non-negative" str) (printf "Minimum overlap portion between pairs of images, default: %f" (minimumOverlap defltOption)) : - Opt.Option [] ["maximum-difference"] + opt generic [] ["maximum-difference"] (flip ReqArg "FRACTION" $ \str flags -> fmap (\x -> flags{maximumDifference = x}) $ parseNumber "maximum difference" (\x -> 0<=x && x<=1) "between 0 and 1" str) (printf "Maximum average difference between overlapping parts, default: %f" (maximumDifference defltOption)) : - Opt.Option [] ["finetune-rotate"] + opt generic [] ["finetune-rotate"] (NoArg $ \flags -> return $ flags{finetuneRotate = True}) (printf "Fine-tune rotation together with overlapping, default: disabled") : - Opt.Option [] ["number-stamps"] + opt generic [] ["number-stamps"] (flip ReqArg "NATURAL" $ \str flags -> fmap (\x -> flags{numberStamps = x}) $ parseNumber "number of stamps" (0<) "positive" str) (printf "Number of stamps in an overlap area, default: %d" (numberStamps defltOption)) : - Opt.Option [] ["stamp-size"] + opt generic [] ["stamp-size"] (flip ReqArg "NATURAL" $ \str flags -> fmap (\x -> flags{stampSize = x}) $ parseNumber "stamp size" (0<) "positive" str) (printf "Size of a stamp, default: %d" (stampSize defltOption)) : - Opt.Option [] ["distance-gamma"] + opt generic [] ["distance-gamma"] (flip ReqArg "FRACTION" $ \str flags -> fmap (\x -> flags{distanceGamma = x}) $ parseNumber "gamma exponent" (0<) "positive" str) (printf "Distance exponent, default: %f" (distanceGamma defltOption)) : + opt knead [] ["shape-smooth"] + (flip ReqArg "NATURAL" $ \str flags -> + fmap (\x -> flags{shapeSmooth = x}) $ + parseNumber "smooth radius" (0<=) "non-negative" str) + (printf "Smooth radius for part shapes, default: %d" + (shapeSmooth defltOption)) : + [] -description :: Description (Image, Args) -> Description (Image, Args) +description :: Description (Image, Args) -> EngineDescription (Image, Args) description desc = map - (fmapOptDescr $ \update (image, old) -> do + (mapSnd $ fmapOptDescr $ \update (image, old) -> do new <- update $ option old return (image, old {option = new})) (optionDescription desc) ++ - Opt.Option [] ["hint-angle"] + opt generic [] ["hint-angle"] (flip ReqArg "DEGREE" $ \str (image, args) -> fmap (\x -> (image{angle = Just x}, args)) $ parseNumber "angle" (\w -> -1000<=w && w<=1000) "degree" str) @@ -234,9 +294,9 @@ -get :: IO Args -get = do - let desc = description desc +get :: Engine -> IO Args +get engine = do + let desc = map snd $ filter (EnumSet.get engine . fst) $ description desc argv <- Env.getArgs let (args, _files, errors) = getOpt (Opt.ReturnInOrder addFile) desc argv when (not $ null errors) $