patch-image 0.1.0.2 → 0.2
raw patch · 10 files changed
+2944/−600 lines, 10 filesdep +arraydep +carraydep +containersdep ~JuicyPixelsdep ~accelerate-cudanew-component:exe:patch-image-cudanew-component:exe:patch-image-llvm
Dependencies added: array, carray, containers, enumset, fft, knead, llvm-extra, llvm-tf, pqueue, tfp
Dependency ranges changed: JuicyPixels, accelerate-cuda
Files
- Changes.md +23/−0
- README.md +196/−0
- patch-image.cabal +76/−179
- src/Accelerate.hs +96/−390
- src/Arithmetic.hs +312/−0
- src/Knead.hs +1706/−0
- src/KneadShape.hs +176/−0
- src/LinearAlgebra.hs +133/−0
- src/MatchImageBorders.hs +135/−0
- src/Option.hs +91/−31
+ Changes.md view
@@ -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.
+ README.md view
@@ -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.
patch-image.cabal view
@@ -1,5 +1,5 @@ Name: patch-image-Version: 0.1.0.2+Version: 0.2 License: BSD3 License-File: LICENSE Author: Henning Thielemann <haskell@henning-thielemann.de>@@ -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
src/Accelerate.hs view
@@ -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<width && 0<=y && y<height- inBox :: (A.Elt a, A.IsNum a, A.IsScalar a) => (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
+ src/Arithmetic.hs view
@@ -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<width && 0<=y && y<height+++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 = 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
+ src/Knead.hs view
@@ -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
+ src/KneadShape.hs view
@@ -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)
+ src/LinearAlgebra.hs view
@@ -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)
+ src/MatchImageBorders.hs view
@@ -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
src/Option.hs view
@@ -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) $