diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,31 @@ +Copyright (c) 2014, Henning Thielemann + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * The names of contributors may not be used to endorse or promote + products derived from this software without specific prior + written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,3 @@ +#! /usr/bin/env runhaskell +> import Distribution.Simple +> main = defaultMain diff --git a/patch-image.cabal b/patch-image.cabal new file mode 100644 --- /dev/null +++ b/patch-image.cabal @@ -0,0 +1,237 @@ +Name: patch-image +Version: 0.1 +License: BSD3 +License-File: LICENSE +Author: Henning Thielemann +Maintainer: Henning Thielemann +Homepage: http://code.haskell.org/~thielema/patch-image/ +Category: Graphics +Synopsis: Compose a big image from overlapping parts +Description: + Compose a collage from overlapping image parts. + In contrast to Hugin, + this is not intended for creating panoramas from multiple photographies, + but instead is specialised to creating highly accurate reconstructions + of flat but big image sources, like record covers, posters or newspapers. + 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. + . + * Images must be approximately correctly oriented. + . + * May have problems with unstructured areas in the image. +Tested-With: GHC==7.8.3 +Cabal-Version: >=1.6 +Build-Type: Simple + +Source-Repository this + Tag: 0.1 + Type: darcs + Location: http://code.haskell.org/~thielema/patch-image/ + +Source-Repository head + Type: darcs + Location: http://code.haskell.org/~thielema/patch-image/ + +Flag buildDraft + description: Build draft program + default: False + +Executable patch-image + Main-Is: Accelerate.hs + Other-Modules: + Option.Utility + Option + Hs-Source-Dirs: src + + GHC-Options: -Wall -threaded -fwarn-tabs -fwarn-incomplete-record-updates + GHC-Prof-Options: -fprof-auto -rtsopts + + Build-Depends: + accelerate-arithmetic >=0.0 && <0.1, + accelerate-utility >=0.0 && <0.1, + accelerate-cuda >=0.15 && <0.16, + accelerate-fft >=0.15 && <0.16, + accelerate-io >=0.15 && <0.16, + accelerate >=0.15 && <0.16, + JuicyPixels >=2.0 && <3.2, + hmatrix >=0.15 && <0.16, + gnuplot >=0.5 && <0.6, + vector >=0.10 && <0.11, + Cabal >=1.18 && <1.22, + filepath >=1.3 && <1.4, + utility-ht >=0.0.1 && <0.1, + base >=4 && <5 + +Executable patch-image-draft + Main-Is: Draft.hs + Hs-Source-Dirs: src + + GHC-Options: -Wall -fwarn-tabs -fwarn-incomplete-record-updates + + If flag(buildDraft) + Build-Depends: + JuicyPixels >=2.0 && <3.2, + GeomAlgLib >=0.2 && <0.3, + utility-ht >=0.0.1 && <0.1, + base >=4 && <5 + Else + Buildable: False diff --git a/src/Accelerate.hs b/src/Accelerate.hs new file mode 100644 --- /dev/null +++ b/src/Accelerate.hs @@ -0,0 +1,1859 @@ +{-# LANGUAGE TypeOperators #-} +module Main where + +import qualified Option + +import qualified Data.Array.Accelerate.Math.FFT as FFT +import qualified Data.Array.Accelerate.Data.Complex as Complex +import qualified Data.Array.Accelerate.CUDA as CUDA +import qualified Data.Array.Accelerate.IO as AIO +import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as LinAlg +import qualified Data.Array.Accelerate.Utility.Lift.Run as Run +import qualified Data.Array.Accelerate.Utility.Lift.Acc as Acc +import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp +import qualified Data.Array.Accelerate.Utility.Arrange as Arrange +import qualified Data.Array.Accelerate.Utility.Loop as Loop +import qualified Data.Array.Accelerate as A +import Data.Array.Accelerate.Data.Complex (Complex((:+)), ) +import Data.Array.Accelerate.Utility.Lift.Exp (atom) +import Data.Array.Accelerate + (Acc, Array, Exp, DIM1, DIM2, DIM3, + (:.)((:.)), Z(Z), Any(Any), All(All), + (<*), (<=*), (>=*), (==*), (&&*), (||*), (?), ) + +import qualified Data.Packed.Matrix as Matrix +import qualified Data.Packed.Vector as Vector +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 + +import qualified Graphics.Gnuplot.Plot.TwoDimensional as Plot2D +import qualified Graphics.Gnuplot.Graph.TwoDimensional as Graph2D + +import qualified Data.Complex as HComplex + +import qualified Codec.Picture as Pic + +import qualified Data.Vector.Storable as SV + +import qualified System.FilePath as FilePath + +import qualified Distribution.Simple.Utils as CmdLine +import Distribution.Verbosity (Verbosity) +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 Data.Maybe.HT (toMaybe) +import Data.Maybe (catMaybes) +import Data.List.HT (removeEach, mapAdjacent, tails) +import Data.Traversable (forM) +import Data.Foldable (forM_, foldMap) +import Data.Tuple.HT (mapPair, mapFst, mapSnd, fst3, thd3) +import Data.Word (Word8) + + +readImage :: Verbosity -> FilePath -> IO (Array DIM3 Word8) +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 $ + AIO.fromVectors + (Z :. Pic.imageHeight pic :. Pic.imageWidth pic :. 3) + ((), dat) + _ -> ioError $ userError "unsupported image type" + +writeImage :: Int -> FilePath -> Array DIM3 Word8 -> IO () +writeImage quality path arr = do + let (Z :. height :. width :. 3) = A.arrayShape arr + Pic.saveJpgImage quality path $ Pic.ImageYCbCr8 $ + Pic.Image { + Pic.imageWidth = width, + Pic.imageHeight = height, + Pic.imageData = snd $ AIO.toVectors arr + } + +writeGrey :: Int -> FilePath -> Array DIM2 Word8 -> IO () +writeGrey quality path arr = do + let (Z :. height :. width) = A.arrayShape arr + Pic.saveJpgImage quality path $ Pic.ImageY8 $ + Pic.Image { + Pic.imageWidth = width, + Pic.imageHeight = height, + Pic.imageData = snd $ AIO.toVectors arr + } + +imageFloatFromByte :: + (A.Shape sh, A.Elt a, A.IsFloating a) => + Acc (Array sh Word8) -> Acc (Array sh a) +imageFloatFromByte = A.map ((/255) . A.fromIntegral) + +imageByteFromFloat :: + (A.Shape sh, A.Elt a, A.IsFloating a) => + Acc (Array sh a) -> Acc (Array sh Word8) +imageByteFromFloat = A.map (fastRound . (255*) . max 0 . min 1) + + +cycleLeftDim3 :: Exp DIM3 -> Exp DIM3 +cycleLeftDim3 = + Exp.modify (atom :. atom :. atom :. atom) $ + \(z :. chans :. height :. width) -> + z :. height :. width :. chans + +cycleRightDim3 :: Exp DIM3 -> Exp DIM3 +cycleRightDim3 = + Exp.modify (atom :. atom :. atom :. atom) $ + \(z :. height :. width :. chans) -> + z :. chans :. height :. width + +separateChannels :: (A.Elt a) => Acc (Array DIM3 a) -> Acc (Array DIM3 a) +separateChannels arr = + A.backpermute + (cycleRightDim3 $ A.shape arr) + cycleLeftDim3 + arr + +interleaveChannels :: (A.Elt a) => Acc (Array DIM3 a) -> Acc (Array DIM3 a) +interleaveChannels arr = + A.backpermute + (cycleLeftDim3 $ A.shape arr) + cycleRightDim3 + arr + + +fastRound :: + (A.Elt i, A.IsIntegral i, A.Elt a, A.IsFloating a) => Exp a -> Exp i +fastRound x = A.floor (x+0.5) + +floatArray :: Acc (Array sh Float) -> Acc (Array sh Float) +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 + in (i, x - A.fromIntegral i) + + + +type Channel ix a = Array (ix :. Int :. Int) a + +type ExpDIM2 ix = Exp ix :. Exp Int :. Exp Int +type ExpDIM3 ix = Exp ix :. Exp Int :. Exp Int :. Exp Int + +unliftDim2 :: + (A.Slice ix) => + Exp (ix :. Int :. Int) -> ExpDIM2 ix +unliftDim2 = A.unlift + + +indexLimit :: + (A.Slice ix, A.Shape ix, A.Elt a) => + Acc (Channel ix a) -> ExpDIM2 ix -> Exp a +indexLimit arr (ix:.y:.x) = + let (_ :. height :. width) = unliftDim2 $ A.shape arr + xc = max 0 $ min (width -1) x + yc = max 0 $ min (height-1) y + in arr A.! A.lift (ix :. yc :. xc) + +indexFrac :: + (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) => + Acc (Channel ix a) -> Exp ix :. Exp a :. Exp a -> Exp a +indexFrac arr (ix:.y:.x) = + let (xi,xf) = splitFraction x + (yi,yf) = splitFraction y + interpolRow yc = + cubicIp + (indexLimit arr (ix:.yc:.xi-1), + indexLimit arr (ix:.yc:.xi ), + indexLimit arr (ix:.yc:.xi+1), + indexLimit arr (ix:.yc:.xi+2)) + xf + in cubicIp + (interpolRow (yi-1), + interpolRow yi, + interpolRow (yi+1), + interpolRow (yi+2)) + yf + + +rotateStretchMoveCoords :: + (A.Elt a, A.IsFloating a) => + (Exp a, Exp a) -> + (Exp a, Exp a) -> + (Exp Int, Exp Int) -> + Acc (Channel Z (a, a)) +rotateStretchMoveCoords rot mov (width,height) = + let trans = rotateStretchMoveBackPoint rot mov + in A.generate (A.lift $ Z:.height:.width) $ \p -> + let (_z :. ydst :. xdst) = unliftDim2 p + in A.lift $ trans (A.fromIntegral xdst, A.fromIntegral ydst) + +inBoxPlain :: + (Ord a, Num a) => + (a, a) -> + (a, a) -> + Bool +inBoxPlain (width,height) (x,y) = + 0<=x && x + (Exp a, Exp a) -> + (Exp a, Exp a) -> + Exp Bool +inBox (width,height) (x,y) = + 0<=*x &&* x<*width &&* 0<=*y &&* y<*height + +validCoords :: + (A.Elt a, A.IsFloating a) => + (Exp Int, Exp Int) -> + Acc (Channel Z (a, a)) -> + Acc (Channel Z Bool) +validCoords (width,height) = + A.map $ A.lift1 $ \(x,y) -> + inBox (width,height) (fastRound x, fastRound y) + +replicateChannel :: + (A.Slice ix, A.Shape ix, A.Elt a) => + Exp ix -> Acc (Channel Z a) -> Acc (Channel ix a) +replicateChannel = LinAlg.extrudeMatrix + +{- | +@rotateStretchMove rot mov@ +first rotate and stretches the image according to 'rot' +and then moves the picture. +-} +rotateStretchMove :: + (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) => + (Exp a, Exp a) -> + (Exp a, Exp a) -> + ExpDIM2 ix -> Acc (Channel ix a) -> + (Acc (Channel Z Bool), Acc (Channel ix a)) +rotateStretchMove rot mov sh arr = + let ( chansDst :. heightDst :. widthDst) = sh + (_chansSrc :. heightSrc :. widthSrc) = unliftDim2 $ A.shape arr + coords = rotateStretchMoveCoords rot mov (widthDst, heightDst) + + in (validCoords (widthSrc, heightSrc) coords, + Arrange.mapWithIndex + (\ix coord -> + let (chan :. _ydst :. _xdst) = unliftDim2 ix + (xsrc,ysrc) = A.unlift coord + in indexFrac arr (chan :. ysrc :. xsrc)) + (replicateChannel chansDst coords)) + + +rotateLeftTop :: + (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) => + (Exp a, Exp a) -> Acc (Channel ix a) -> + ((Acc (A.Scalar a), Acc (A.Scalar a)), Acc (Channel ix a)) +rotateLeftTop rot arr = + let (chans :. height :. width) = unliftDim2 $ A.shape arr + ((left, right), (top, bottom)) = + boundingBoxOfRotated rot (A.fromIntegral width, A.fromIntegral height) + in ((A.unit left, A.unit top), + snd $ + rotateStretchMove rot (-left,-top) + (chans :. A.ceiling (bottom-top) :. A.ceiling (right-left)) arr) + +rotate :: + (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) => + (Exp a, Exp a) -> + Acc (Channel ix a) -> Acc (Channel ix a) +rotate rot arr = snd $ rotateLeftTop rot arr + + +brightnessPlane :: + (A.Slice ix, A.Shape ix) => + Acc (Channel (ix:.Int) Float) -> Acc (Channel ix Float) +brightnessPlane = flip A.slice (A.lift (Any :. (0::Int) :. All :. All)) + +rowHistogram :: Acc (Channel DIM1 Float) -> Acc (Array DIM1 Float) +rowHistogram = A.fold (+) 0 . brightnessPlane + + +rotateHistogram :: + Float -> Array DIM3 Word8 -> (Array DIM3 Word8, Array DIM1 Float) +rotateHistogram = + let rot = + Run.with CUDA.run1 $ \orient arr -> + let rotated = + rotate orient $ + separateChannels $ imageFloatFromByte arr + in (imageByteFromFloat $ interleaveChannels rotated, + rowHistogram rotated) + 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 <- + forM angles $ \degree -> do + let (rotated, histogram) = rotateHistogram (degree * pi/180) pic + let stem = printf "rotated%+07.2f" degree + writeImage 90 ("/tmp/" ++ stem ++ ".jpeg") rotated + let diffHistogram = map abs $ mapAdjacent (-) $ A.toList histogram + printf "%s: maxdiff %8.3f, sqrdiff %8.0f\n" + stem (maximum diffHistogram) (sum $ map (^(2::Int)) diffHistogram) + return (stem, histogram) + void $ GP.plotDefault $ + foldMap + (\(label, histogram) -> + fmap (Graph2D.lineSpec (LineSpec.title label LineSpec.deflt)) $ + Plot2D.list Graph2D.listLines $ A.toList histogram) + histograms + void $ GP.plotDefault $ + foldMap + (\(label, histogram) -> + fmap (Graph2D.lineSpec (LineSpec.title label LineSpec.deflt)) $ + Plot2D.list Graph2D.listLines $ + map abs $ mapAdjacent (-) $ A.toList histogram) + histograms + + + +differentiate :: + (A.Elt a, A.IsNum a) => + Acc (Array DIM1 a) -> Acc (Array DIM1 a) +differentiate arr = + let size = A.unindex1 $ A.shape arr + in A.generate (A.index1 (size-1)) $ \i -> + arr A.! (A.index1 $ A.unindex1 i + 1) - arr A.! i + +scoreRotation :: Float -> Array DIM3 Word8 -> Float +scoreRotation = + let rot = + Run.with CUDA.run1 $ \orient arr -> + A.sum $ A.map (^(2::Int)) $ differentiate $ rowHistogram $ + rotate orient $ separateChannels $ imageFloatFromByte arr + in \angle arr -> Acc.the $ rot (cos angle, sin angle) arr + +findOptimalRotation :: [Float] -> Array DIM3 Word8 -> Float +findOptimalRotation angles pic = + Key.maximum (flip scoreRotation pic . (* (pi/180))) angles + + + + + +rotateManifest :: Float -> Array DIM3 Word8 -> Array DIM3 Float +rotateManifest = + let rot = + Run.with CUDA.run1 $ \orient arr -> + rotate orient $ separateChannels $ imageFloatFromByte arr + in \angle arr -> rot (cos angle, sin angle) arr + + +prepareOverlapMatching :: + Int -> (Float, Array DIM3 Word8) -> ((Float,Float), Channel Z Float) +prepareOverlapMatching = + let rot = + Run.with CUDA.run1 $ \radius orient arr -> + rotateLeftTop orient $ + (if True + then highpass radius + else removeDCOffset) $ + brightnessPlane $ separateChannels $ imageFloatFromByte arr + in \radius (angle, arr) -> + mapFst (mapPair (Acc.the, Acc.the)) $ + rot radius (cos angle, sin angle) arr + + +ceilingPow2Exp :: Exp Int -> Exp Int +ceilingPow2Exp n = + A.setBit 0 $ A.ceiling $ logBase 2 (fromIntegral n :: Exp Double) + +pad :: + (A.Elt a) => + Exp a -> Exp DIM2 -> Acc (Channel Z a) -> Acc (Channel Z a) +pad a sh arr = + let (height, width) = A.unlift $ A.unindex2 $ A.shape arr + in A.generate sh $ \p -> + let (y, x) = A.unlift $ A.unindex2 p + in (y<*height &&* x<*width) + ? + (arr A.! A.index2 y x, a) + +convolveImpossible :: + (A.Elt a, A.IsFloating a) => + Acc (Channel Z a) -> Acc (Channel Z a) -> Acc (Channel Z a) +convolveImpossible x y = + let (heightx, widthx) = A.unlift $ A.unindex2 $ A.shape x + (heighty, widthy) = A.unlift $ A.unindex2 $ A.shape y + width = ceilingPow2Exp $ widthx + widthy + height = ceilingPow2Exp $ heightx + heighty + sh = A.index2 height width + forward z = + FFT.fft2D FFT.Forward $ CUDA.run $ + A.map (A.lift . (:+ 0)) $ pad 0 sh z + in A.map Complex.real $ + FFT.fft2D FFT.Inverse $ CUDA.run $ + A.zipWith (*) (forward x) (forward y) + + +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 = + let sh = A.shape arr + (_z :. height :. width) = unliftDim2 sh + s = + A.the (A.fold1All (+) arr) + / (A.fromIntegral width * A.fromIntegral height) + in A.map (subtract s) arr + +{- +We cannot remove DC offset in the spectrum, +because we already padded the images with zeros. +-} +clearDCCoefficient :: + (A.Elt a, A.IsFloating a) => + Acc (Array DIM2 (Complex a)) -> Acc (Array DIM2 (Complex a)) +clearDCCoefficient arr = + A.generate (A.shape arr) $ \p -> + let (_z:.y:.x) = unliftDim2 p + in x==*0 ||* y==*0 ? (0, arr A.! 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) +lowpass count = + Loop.nest count $ + A.stencil (\(a,m,b) -> smooth3 (smooth3 a, smooth3 m, smooth3 b)) A.Clamp + +highpass count arr = + A.zipWith (-) arr $ lowpass count arr + + +convolvePaddedSimple :: + (A.Elt a, A.IsFloating a) => + DIM2 -> Acc (Channel Z a) -> Acc (Channel Z a) -> Acc (Channel Z a) +convolvePaddedSimple sh@(Z :. height :. width) = + let forward = + FFT.fft2D' FFT.Forward width height . + A.map (A.lift . (:+ 0)) . pad 0 (A.lift sh) + inverse = FFT.fft2D' FFT.Inverse width height + in \ x y -> + A.map Complex.real $ inverse $ + A.zipWith (\xi yi -> xi * Complex.conjugate yi) (forward x) (forward y) + + +imagUnit :: (A.Elt a, A.IsNum a) => Exp (Complex a) +imagUnit = Exp.modify2 atom atom (:+) 0 1 + +{- | +Let f and g be two real valued images. +The spectrum of f+i*g is spec f + i * spec g. +Let 'flip' be the spectrum with negated indices modulo image size. +It holds: flip (spec f) = conj (spec f). + +(a + conj b) / 2 + = (spec (f+i*g) + conj (flip (spec (f+i*g)))) / 2 + = (spec f + i*spec g + conj (flip (spec f)) + conj (flip (spec (i*g)))) / 2 + = (2*spec f + i*spec g + conj (i*flip (spec g))) / 2 + = (2*spec f + i*spec g - i * conj (flip (spec g))) / 2 + = spec f + +(a - conj b) * (-i/2) + = (-i*a + conj (-i*b)) / 2 + -> this swaps role of f and g in the proof above +-} +untangleRealSpectra :: + (A.Elt a, A.IsFloating a) => + Acc (Array DIM2 (Complex a)) -> Acc (Array DIM2 (Complex a, Complex a)) +untangleRealSpectra spec = + A.zipWith + (\a b -> + A.lift $ + ((a + Complex.conjugate b) / 2, + (a - Complex.conjugate b) * (-imagUnit / 2))) + spec $ + A.backpermute (A.shape spec) + (Exp.modify (atom:.atom:.atom) $ + \(_z:.y:.x) -> + let (_z:.height:.width) = unliftDim2 $ A.shape spec + in Z :. mod (-y) height :. mod (-x) width) + spec + +{- +This is more efficient than 'convolvePaddedSimple' +since it needs only one forward Fourier transform, +where 'convolvePaddedSimple' needs two of them. +For the analysis part, +perform two real-valued Fourier transforms using one complex-valued transform. +Afterwards we untangle the superposed spectra. +-} +convolvePadded :: + (A.Elt a, A.IsFloating a) => + DIM2 -> Acc (Channel Z a) -> Acc (Channel Z a) -> Acc (Channel Z a) +convolvePadded sh@(Z :. height :. width) = + let forward = FFT.fft2D' FFT.Forward width height + inverse = FFT.fft2D' FFT.Inverse width height + in \ a b -> + A.map Complex.real $ inverse $ + A.map (Exp.modify (atom,atom) $ \(ai,bi) -> ai * Complex.conjugate bi) $ + untangleRealSpectra $ forward $ + pad 0 (A.lift sh) $ + A.zipWith (Exp.modify2 atom atom (:+)) a b + + +attachDisplacements :: + (A.Elt a, A.IsScalar a) => + Exp Int -> Exp Int -> + Acc (Channel Z a) -> Acc (Channel Z ((Int, Int), a)) +attachDisplacements xsplit ysplit arr = + let sh = A.shape arr + (_z :. height :. width) = unliftDim2 sh + in A.generate sh $ \p -> + let (_z:.y:.x) = unliftDim2 p + wrap size split c = c<*split ? (c, c-size) + in A.lift ((wrap width xsplit x, wrap height ysplit y), arr A.! p) + +weightOverlapScores :: + (A.Elt a, A.IsFloating a, A.IsScalar a) => + Exp Int -> (Exp Int, Exp Int) -> (Exp Int, Exp Int) -> + Acc (Channel Z ((Int, Int), a)) -> + Acc (Channel Z ((Int, Int), a)) +weightOverlapScores minOverlap (widtha,heighta) (widthb,heightb) = + A.map + (Exp.modify ((atom,atom),atom) $ \(dp@(dy,dx),v) -> + let clipWidth = min widtha (widthb + dx) - max 0 dx + clipHeight = min heighta (heightb + dy) - max 0 dy + in (dp, + (clipWidth >=* minOverlap &&* clipHeight >=* minOverlap) + ? + (v / (A.fromIntegral clipWidth * A.fromIntegral clipHeight), 0))) + +{- | +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 :: + (A.Elt a, A.IsFloating a, A.IsScalar a) => + Exp Int -> (Exp Int, Exp Int) -> (Exp Int, Exp Int) -> + Acc (Channel Z ((Int, Int), a)) -> + Acc (Channel Z ((Int, Int), a)) +minimumOverlapScores minOverlap (widtha,heighta) (widthb,heightb) = + A.map + (Exp.modify ((atom,atom),atom) $ \(dp@(dy,dx),v) -> + let clipWidth = min widtha (widthb + dx) - max 0 dx + clipHeight = min heighta (heightb + dy) - max 0 dy + in (dp, + (clipWidth >=* minOverlap &&* clipHeight >=* minOverlap) + ? + (v, 0))) + +argmax :: + (A.Elt a, A.Elt b, A.IsScalar b) => + Exp (a, b) -> Exp (a, b) -> Exp (a, b) +argmax x y = A.snd x <* A.snd y ? (y,x) + +argmaximum :: + (A.Elt a, A.Elt b, A.IsScalar b) => + Acc (Channel Z (a, b)) -> Acc (A.Scalar (a, b)) +argmaximum = A.fold1All argmax + + +allOverlaps :: + DIM2 -> + Exp Float -> + Acc (Channel Z Float) -> Acc (Channel Z Float) -> + Acc (Channel Z ((Int, Int), Float)) +allOverlaps size@(Z :. height :. width) minOverlapPortion = + let convolve = convolvePadded size + in \a b -> + let (Z :. heighta :. widtha) = A.unlift $ A.shape a + (Z :. heightb :. widthb) = A.unlift $ A.shape b + half = flip div 2 + minOverlap = + fastRound $ + minOverlapPortion + * + A.fromIntegral + (min + (min widtha heighta) + (min widthb heightb)) + weight = + if False + then + weightOverlapScores minOverlap + (widtha, heighta) + (widthb, heightb) + else + minimumOverlapScores minOverlap + (widtha, heighta) + (widthb, heightb) + in weight $ + attachDisplacements + (half $ A.lift width - widthb + widtha) + (half $ A.lift height - heightb + heighta) $ + convolve a b + + +allOverlapsRun :: + DIM2 -> Float -> Channel Z Float -> Channel Z Float -> Channel Z Word8 +allOverlapsRun padExtent = + Run.with CUDA.run1 $ \minOverlap picA picB -> + imageByteFromFloat $ + -- A.map (2*) $ + A.map (0.0001*) $ + A.map A.snd $ allOverlaps padExtent minOverlap picA picB + +optimalOverlap :: + DIM2 -> Float -> Channel Z Float -> Channel Z Float -> ((Int, Int), Float) +optimalOverlap padExtent = + let run = + Run.with CUDA.run1 $ \minimumOverlap a b -> + argmaximum $ allOverlaps padExtent minimumOverlap a b + in \overlap a b -> Acc.the $ run overlap a b + + +shrink :: + (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) => + GenDIM2 (Exp Int) -> Acc (Channel ix a) -> Acc (Channel ix a) +shrink (_:.yk:.xk) arr = + let (shape:.height:.width) = unliftDim2 $ A.shape arr + in A.map (/ (A.fromIntegral xk * A.fromIntegral yk)) $ + A.fold1 (+) $ A.fold1 (+) $ + A.backpermute + (A.lift $ shape :. div height yk :. div width xk :. yk :. xk) + (Exp.modify (atom:.atom:.atom:.atom:.atom) $ + \(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 +shrinkFactors (Z:.heightPad:.widthPad) + (Z :. heighta :. widtha) (Z :. heightb :. widthb) = + let yk = divUp (heighta+heightb) $ fromIntegral heightPad + xk = divUp (widtha +widthb) $ fromIntegral widthPad + in Z :. yk :. xk + +{- +Reduce image sizes below the padExtent before matching images. +-} +optimalOverlapBig :: + DIM2 -> Float -> Channel Z Float -> Channel Z Float -> ((Int, Int), Float) +optimalOverlapBig padExtent = + let run = + Run.with CUDA.run1 $ \minimumOverlap a b -> + let factors@(_z:.yk:.xk) = + shrinkFactors padExtent + (A.unlift $ A.shape a) (A.unlift $ A.shape b) + scalePos = + Exp.modify ((atom,atom), atom) $ + \((xm,ym), score) -> ((xm*xk, ym*yk), score) + in A.map scalePos $ argmaximum $ + allOverlaps padExtent minimumOverlap + (shrink factors a) (shrink factors b) + in \minimumOverlap a b -> Acc.the $ run minimumOverlap a b + + +clip :: + (A.Slice ix, A.Shape ix, A.Elt a) => + (Exp Int, Exp Int) -> + (Exp Int, Exp Int) -> + Acc (Channel ix a) -> Acc (Channel ix a) +clip (left,top) (width,height) arr = + A.backpermute + (A.lift $ A.indexTail (A.indexTail (A.shape arr)) :. height :. width) + (Exp.modify (atom:.atom:.atom) $ + \(z :. y :. x) -> z :. y+top :. x+left) + arr + + +overlappingArea :: + (Ord a, Num a) => + GenDIM2 a -> + GenDIM2 a -> + (a, a) -> ((a, a), (a, a), (a, a)) +overlappingArea (Z :. heighta :. widtha) (Z :. 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 -> Float -> Channel Z Float -> Channel Z Float -> ((Int, Int), Float) +optimalOverlapBigFine padExtent@(Z:.heightPad:.widthPad) = + let run = + Run.with CUDA.run1 $ \minimumOverlap a b -> + let shapeA = A.unlift $ A.shape a + shapeB = A.unlift $ A.shape b + factors@(_z:.yk:.xk) = shrinkFactors padExtent shapeA shapeB + coarsed@(coarsedx,coarsedy) = + mapPair ((xk*), (yk*)) $ + Exp.unliftPair $ A.fst $ A.the $ argmaximum $ + allOverlaps padExtent minimumOverlap + (shrink factors a) (shrink factors b) + + ((leftOverlap, topOverlap), _, + (widthOverlap, heightOverlap)) + = overlappingArea shapeA shapeB coarsed + + widthFocus = min widthOverlap $ A.lift $ div widthPad 2 + heightFocus = min heightOverlap $ A.lift $ div heightPad 2 + extentFocus = (widthFocus,heightFocus) + leftFocus = leftOverlap + div (widthOverlap-widthFocus) 2 + topFocus = topOverlap + div (heightOverlap-heightFocus) 2 + addCoarsePos = + Exp.modify ((atom,atom), atom) $ + \((xm,ym), score) -> ((xm+coarsedx, ym+coarsedy), score) + in A.map addCoarsePos $ argmaximum $ + allOverlaps padExtent minimumOverlap + (clip (leftFocus,topFocus) extentFocus a) + (clip (leftFocus-coarsedx,topFocus-coarsedy) extentFocus b) + in \minimumOverlap a b -> Acc.the $ run minimumOverlap a b + + +{- +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 -> + Float -> Float -> Channel Z Float -> Channel Z Float -> + [((Int, Int), (Int, Int), Float)] +optimalOverlapBigMulti padExtent (Z:.heightStamp:.widthStamp) numCorrs = + let overlapShrunk = + Run.with CUDA.run1 $ + \minimumOverlap factors a b -> + argmaximum $ + allOverlaps padExtent minimumOverlap + (shrink factors a) (shrink factors b) + diffShrunk = + Run.with CUDA.run1 $ + \shrunkd factors a b -> + overlapDifference shrunkd + (shrink factors a) (shrink factors b) + + allOverlapsFine = allOverlaps (Z :. 2*heightStamp :. 2*widthStamp) + overlapFine = + Run.with CUDA.run1 $ + \minimumOverlap a b anchorA@(leftA, topA) anchorB@(leftB, topB) + extent@(width,height) -> + let addCoarsePos = + Exp.modify ((atom,atom), atom) $ + \((xm,ym), score) -> + let xc = div (width+xm) 2 + yc = div (height+ym) 2 + in ((leftA+xc, topA+yc), + (leftB+xc-xm, topB+yc-ym), + score) + in A.map addCoarsePos $ argmaximum $ + allOverlapsFine minimumOverlap + (clip anchorA extent a) + (clip anchorB extent b) + + in \maximumDiff minimumOverlap a b -> + let factors@(Z:.yk:.xk) = + shrinkFactors padExtent (A.arrayShape a) (A.arrayShape b) + + ((shrunkdx, shrunkdy), _score) = + Acc.the $ overlapShrunk minimumOverlap factors a b + + coarsedx = shrunkdx * xk + coarsedy = shrunkdy * yk + coarsed = (coarsedx,coarsedy) + + diff = Acc.the $ diffShrunk (shrunkdx, shrunkdy) factors a b + + ((leftOverlap, topOverlap), + (rightOverlap, bottomOverlap), + (widthOverlap, heightOverlap)) + = overlappingArea (A.arrayShape a) (A.arrayShape b) coarsed + + widthStampClip = min widthOverlap widthStamp + heightStampClip = min heightOverlap heightStamp + + in (if diff < maximumDiff then id else const []) $ + map + (\(x,y) -> + Acc.the $ + overlapFine minimumOverlap a b + (x, y) (x-coarsedx, y-coarsedy) + (widthStampClip, heightStampClip)) $ + zip + (map round $ tail $ init $ + linearScale (numCorrs+1) + (fromIntegral leftOverlap :: Double, + fromIntegral $ rightOverlap - widthStampClip)) + (map round $ tail $ init $ + linearScale (numCorrs+1) + (fromIntegral topOverlap :: Double, + fromIntegral $ bottomOverlap - heightStampClip)) + + +overlapDifference :: + (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) => + (Exp Int, Exp Int) -> + Acc (Channel ix a) -> Acc (Channel ix a) -> Acc (A.Scalar a) +overlapDifference (dx,dy) a b = + let (_ :. heighta :. widtha) = unliftDim2 $ A.shape a + (_ :. heightb :. widthb) = unliftDim2 $ A.shape b + leftOverlap = max 0 dx + topOverlap = max 0 dy + rightOverlap = min widtha (widthb + dx) + bottomOverlap = min heighta (heightb + dy) + widthOverlap = rightOverlap - leftOverlap + heightOverlap = bottomOverlap - topOverlap + extentOverlap = (widthOverlap,heightOverlap) + in A.map sqrt $ + A.map (/(A.fromIntegral widthOverlap * A.fromIntegral heightOverlap)) $ + A.fold1All (+) $ + A.map (^(2::Int)) $ + A.zipWith (-) + (clip (leftOverlap,topOverlap) extentOverlap a) + (clip (leftOverlap-dx,topOverlap-dy) extentOverlap b) + +overlapDifferenceRun :: + (Int, Int) -> + Channel Z Float -> Channel Z Float -> Float +overlapDifferenceRun = + let diff = Run.with CUDA.run1 overlapDifference + 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) -> + (Acc (Channel ix Float), Acc (Channel ix Float)) -> Acc (Channel ix Float) +overlap2 (dx,dy) (a,b) = + let (chansa :. heighta :. widtha) = unliftDim2 $ A.shape a + (chansb :. heightb :. widthb) = unliftDim2 $ A.shape b + left = min 0 dx; right = max widtha (widthb + dx) + top = min 0 dy; bottom = max heighta (heightb + dy) + width = right - left + height = bottom - top + chans = A.intersect chansa chansb + in A.generate (A.lift (chans :. height :. width)) $ \p -> + let (chan :. y :. x) = unliftDim2 p + xa = x + left; xb = xa-dx + 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 + in inPicA ? + (inPicB ? ((a A.! pa + b A.! pb)/2, a A.! pa), + inPicB ? (b A.! pb, 0)) + +composeOverlap :: + (Int, Int) -> + ((Float, Array DIM3 Word8), (Float, Array DIM3 Word8)) -> + Array DIM3 Word8 +composeOverlap = + let rot (angle,pic) = + rotate (cos angle, sin angle) $ + separateChannels $ imageFloatFromByte pic + in Run.with CUDA.run1 $ + \(dx,dy) (a,b) -> + imageByteFromFloat $ interleaveChannels $ + overlap2 (dx, dy) (rot a, rot b) + + +emptyCanvas :: + (A.Slice ix, A.Shape ix) => + ix :. Int :. Int -> + (Channel Z Int, Channel ix Float) +emptyCanvas = + 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 :: + (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) = + (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)) + +updateCanvas :: + (Float,Float) -> (Float,Float) -> Array DIM3 Word8 -> + (Channel Z Int, Channel DIM1 Float) -> + (Channel Z Int, Channel DIM1 Float) +updateCanvas = + Run.with CUDA.run1 $ + \rot mov pic (count,canvas) -> + addToCanvas + (rotateStretchMove rot mov (unliftDim2 $ A.shape canvas) $ + separateChannels $ imageFloatFromByte pic) + (count,canvas) + +finalizeCanvas :: (Channel Z Int, Channel DIM1 Float) -> Array DIM3 Word8 +finalizeCanvas = + Run.with CUDA.run1 $ + \(count, canvas) -> + imageByteFromFloat $ interleaveChannels $ + A.zipWith (/) canvas $ + replicateChannel (A.indexTail $ A.indexTail $ A.shape canvas) $ + A.map A.fromIntegral count + + + +maybePlus :: + (A.Elt a) => + (Exp a -> Exp a -> Exp a) -> + Exp (Bool, a) -> Exp (Bool, a) -> Exp (Bool, a) +maybePlus f x y = + let (xb,xv) = Exp.unliftPair x + (yb,yv) = Exp.unliftPair y + in A.cond xb (A.lift (True, A.cond yb (f xv yv) xv)) y + +maskedMinimum :: + (A.Shape ix, A.Elt a, A.IsScalar a) => + LinAlg.Vector ix (Bool, a) -> + LinAlg.Scalar ix (Bool, a) +maskedMinimum = A.fold1 (maybePlus min) + +maskedMaximum :: + (A.Shape ix, A.Elt a, A.IsScalar a) => + LinAlg.Vector ix (Bool, a) -> + LinAlg.Scalar ix (Bool, a) +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) -> + (Point2 (Exp a), Point2 (Exp a)) -> + (Exp Bool, Point2 (Exp a)) +project x ab = + let (r, y) = projectPerp x ab + in (0<=*r &&* r<=*1, y) + + +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) +distanceMapEdges sh edges = + A.map (Exp.modify (atom,atom) $ \(valid, dist) -> valid ? (dist, 0)) $ + maskedMinimum $ + outerVector + (Exp.modify2 (atom,atom) ((atom, atom), (atom, atom)) $ \p (q0, q1) -> + mapSnd (distance p) $ project p (q0, q1)) + (pixelCoordinates sh) + edges + +distanceMapEdgesRun :: + DIM2 -> Array DIM1 ((Float,Float),(Float,Float)) -> Channel Z Word8 +distanceMapEdgesRun = + Run.with CUDA.run1 $ \sh -> + imageByteFromFloat . A.map (0.01*) . distanceMapEdges sh + +distanceMapBox :: + (A.Elt a, A.IsFloating a) => + Exp DIM2 -> + Exp ((a,a), (a,a), (Int,Int)) -> + Acc (Channel Z (Bool, (((a,(a,a)), (a,(a,a))), ((a,(a,a)), (a,(a,a)))))) +distanceMapBox sh geom = + let (rot, mov, extent@(width,height)) = + Exp.unlift ((atom,atom),(atom,atom),(atom,atom)) geom + widthf = A.fromIntegral width + heightf = A.fromIntegral height + back = rotateStretchMoveBackPoint rot mov + forth = rotateStretchMovePoint rot mov + in A.generate sh $ \p -> + let _z:.y:.x = unliftDim2 p + (xsrc,ysrc) = back (A.fromIntegral x, A.fromIntegral y) + leftDist = max 0 xsrc + rightDist = max 0 $ widthf - xsrc + topDist = max 0 ysrc + bottomDist = max 0 $ heightf - ysrc + in A.lift $ + (inBox extent (fastRound xsrc, fastRound ysrc), + (((leftDist, forth (0,ysrc)), + (rightDist, forth (widthf,ysrc))), + ((topDist, forth (xsrc,0)), + (bottomDist, forth (xsrc,heightf))))) + + +-- cf. Data.Array.Accelerate.Arithmetic.Interpolation +outerVector :: + (A.Slice ix, A.Shape ix, A.Elt a, A.Elt b, A.Elt c) => + (Exp a -> Exp b -> Exp c) -> + LinAlg.Scalar ix a -> LinAlg.Vector Z b -> LinAlg.Vector ix c +outerVector f x y = + A.zipWith f + (A.replicate (A.lift $ Any :. LinAlg.numElems y) x) + (LinAlg.extrudeVector (A.shape x) y) + +separateDistanceMap :: + (A.Elt a) => + Acc (Channel Z (Bool, ((a, a), (a, a)))) -> + Acc (Array DIM3 (Bool, a)) +separateDistanceMap arr = + outerVector + (Exp.modify2 (atom, ((atom, atom), (atom, atom))) (atom,atom) $ + \(b,(horiz,vert)) (orient,side) -> + (b, orient ? (side ? horiz, side ? vert))) + arr + (A.use $ A.fromList (Z:.(4::Int)) $ + liftM2 (,) [False,True] [False,True]) + +distanceMapBoxRun :: + DIM2 -> ((Float,Float),(Float,Float),(Int,Int)) -> Channel Z Word8 +distanceMapBoxRun = + Run.with CUDA.run1 $ \sh geom -> + let scale = + (4/) $ A.fromIntegral $ uncurry min $ + Exp.unliftPair $ Exp.thd3 geom + in imageByteFromFloat $ + A.map (Exp.modify (atom,atom) $ + \(valid, dist) -> valid ? (scale*dist, 0)) $ + maskedMinimum $ + A.map (Exp.mapSnd A.fst) $ + separateDistanceMap $ + distanceMapBox sh geom + + +-- maybe move to Accelerate.Utility +{- | +We use it as a work-around. +Fusion of 'fold1' and 'replicate' would be very welcome +but it seems to fail with current accelerate version. +-} +breakFusion :: (A.Arrays a) => Acc a -> Acc a +breakFusion = id A.>-> id + +array1FromList :: (A.Elt a) => [a] -> Array DIM1 a +array1FromList xs = A.fromList (Z :. length xs) xs + + +containedAnywhere :: + (A.Elt a, A.IsFloating a) => + Acc (Array DIM1 ((a,a), (a,a), (Int,Int))) -> + Acc (Array DIM3 (a,a)) -> + Acc (Array DIM3 Bool) +containedAnywhere geoms arr = + A.fold1 (||*) $ + breakFusion $ + outerVector + (Exp.modify2 (atom,atom) ((atom,atom),(atom,atom),(atom,atom)) $ + \(xdst,ydst) (rot, mov, extent) -> + let (xsrc,ysrc) = rotateStretchMoveBackPoint rot mov (xdst,ydst) + in inBox extent (fastRound xsrc, fastRound ysrc)) + arr geoms + + +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))) -> + Acc (Channel Z a) +distanceMapContained sh this others = + let distMap = + separateDistanceMap $ + distanceMapBox sh this + contained = + containedAnywhere others $ + A.map (A.snd . A.snd) distMap + in A.map (Exp.modify (atom,atom) $ + \(valid, dist) -> valid ? (dist, 0)) $ + maskedMinimum $ + A.zipWith + (Exp.modify2 atom (atom,(atom,atom)) $ \c (b,(dist,_)) -> + (c&&*b, dist)) + contained distMap + +distanceMapContainedRun :: + DIM2 -> + ((Float,Float),(Float,Float),(Int,Int)) -> + [((Float,Float),(Float,Float),(Int,Int))] -> + Channel Z Word8 +distanceMapContainedRun = + let distances = + Run.with CUDA.run1 $ + \sh this others -> + let scale = + (4/) $ A.fromIntegral $ uncurry min $ + Exp.unliftPair $ Exp.thd3 this + in imageByteFromFloat $ A.map (scale*) $ + distanceMapContained sh this others + in \sh this others -> + distances sh this $ array1FromList others + + +pixelCoordinates :: + (A.Elt a, A.IsFloating a) => + Exp DIM2 -> Acc (Channel Z (a,a)) +pixelCoordinates sh = + A.generate sh $ Exp.modify (atom:.atom:.atom) $ \(_z:.y:.x) -> + (A.fromIntegral x, A.fromIntegral y) + +distanceMapPoints :: + (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) => + Acc (Array ix (a,a)) -> + Acc (Array DIM1 (a,a)) -> + Acc (Array ix a) +distanceMapPoints a b = + A.fold1 min $ + outerVector + (Exp.modify2 (atom,atom) (atom,atom) distance) + a b + +distanceMapPointsRun :: + DIM2 -> + [Point2 Float] -> + Channel Z Word8 +distanceMapPointsRun = + let distances = + Run.with CUDA.run1 $ + \sh points -> + let scale = + case Exp.unlift (atom:.atom:.atom) sh of + _z:.y:.x -> (4/) $ A.fromIntegral $ min x y + in imageByteFromFloat $ A.map (scale*) $ + distanceMapPoints (pixelCoordinates sh) points + in \sh points -> + distances sh $ array1FromList points + + +{- | +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 :: + (A.Elt a, A.IsFloating a) => + Exp DIM2 -> + Exp ((a, a), (a, a), (Int, Int)) -> + Acc (Array DIM1 ((a, a), (a, a), (Int, Int))) -> + Acc (Array DIM1 (a, a)) -> + Acc (Channel Z a) +distanceMap sh this others points = + A.zipWith min + (distanceMapContained sh this others) + (distanceMapPoints (pixelCoordinates sh) points) + +distanceMapRun :: + DIM2 -> + ((Float,Float),(Float,Float),(Int,Int)) -> + [((Float,Float),(Float,Float),(Int,Int))] -> + [Point2 Float] -> + Channel Z Word8 +distanceMapRun = + let distances = + Run.with CUDA.run1 $ + \sh this others points -> + let scale = + case Exp.unlift (atom:.atom:.atom) sh of + _z:.y:.x -> (4/) $ A.fromIntegral $ min x y + in imageByteFromFloat $ A.map (scale*) $ + distanceMap sh this others points + in \sh this others points -> + distances sh this + (array1FromList others) + (array1FromList points) + + +distanceMapGamma :: + (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))) -> + Acc (Array DIM1 (a, a)) -> + Acc (Channel Z a) +distanceMapGamma gamma sh this others points = + A.map (**gamma) $ distanceMap sh this others points + + +emptyWeightedCanvas :: + (A.Slice ix, A.Shape ix) => + ix :. Int :. Int -> + (Channel Z Float, Channel ix Float) +emptyWeightedCanvas = + Run.with CUDA.run1 $ \sh -> + let (_ix :. height :. width) = unliftDim2 sh + in (A.fill (A.lift $ Z:.height:.width) 0, + A.fill sh 0) + + +addToWeightedCanvas :: + (A.Slice ix, A.Shape ix, A.Elt a, A.IsNum a) => + (Acc (Channel Z a), Acc (Channel ix a)) -> + (Acc (Channel Z a), Acc (Channel ix a)) -> + (Acc (Channel Z a), Acc (Channel ix a)) +addToWeightedCanvas (weight, pic) (weightSum, canvas) = + (A.zipWith (+) weight weightSum, + A.zipWith (+) canvas $ A.zipWith (*) pic $ + replicateChannel + (A.indexTail $ A.indexTail $ A.shape pic) + weight) + +-- launch timeout +updateWeightedCanvasMerged :: + ((Float,Float),(Float,Float),(Int,Int)) -> + [((Float,Float),(Float,Float),(Int,Int))] -> + [Point2 Float] -> + Array DIM3 Word8 -> + (Channel Z Float, Channel DIM1 Float) -> + (Channel Z Float, Channel DIM1 Float) +updateWeightedCanvasMerged = + let update = + Run.with CUDA.run1 $ + \this others points pic (weightSum,canvas) -> + let (rot, mov, _) = + Exp.unlift ((atom,atom), (atom,atom), atom) this + in addToWeightedCanvas + (distanceMap (A.shape weightSum) this others points, + snd $ rotateStretchMove rot mov (unliftDim2 $ A.shape canvas) $ + separateChannels $ imageFloatFromByte pic) + (weightSum,canvas) + in \this others points pic canvas -> + update this + (array1FromList others) + (array1FromList points) + pic canvas + +updateWeightedCanvas :: + Float -> + ((Float,Float),(Float,Float),(Int,Int)) -> + [((Float,Float),(Float,Float),(Int,Int))] -> + [Point2 Float] -> + Array DIM3 Word8 -> + (Channel Z Float, Channel DIM1 Float) -> + (Channel Z Float, Channel DIM1 Float) +updateWeightedCanvas = + let distances = Run.with CUDA.run1 distanceMapGamma + update = + Run.with CUDA.run1 $ + \this pic dist (weightSum,canvas) -> + let (rot, mov, _) = + Exp.unlift ((atom,atom), (atom,atom), atom) this + in addToWeightedCanvas + (dist, + snd $ rotateStretchMove rot mov (unliftDim2 $ A.shape canvas) $ + separateChannels $ imageFloatFromByte pic) + (weightSum,canvas) + in \gamma this others points pic (weightSum,canvas) -> + update this pic + (distances gamma (A.arrayShape weightSum) this + (array1FromList others) + (array1FromList points)) + (weightSum,canvas) + +-- launch timeout +updateWeightedCanvasSplit :: + ((Float,Float),(Float,Float),(Int,Int)) -> + [((Float,Float),(Float,Float),(Int,Int))] -> + [Point2 Float] -> + Array DIM3 Word8 -> + (Channel Z Float, Channel DIM1 Float) -> + (Channel Z Float, Channel DIM1 Float) +updateWeightedCanvasSplit = + let update = Run.with CUDA.run1 addToWeightedCanvas + distances = Run.with CUDA.run1 distanceMap + rotated = + Run.with CUDA.run1 $ + \sh rot mov pic -> + snd $ rotateStretchMove rot mov (unliftDim2 sh) $ + separateChannels $ imageFloatFromByte pic + in \this@(rot, mov, _) others points pic (weightSum,canvas) -> + update + (distances (A.arrayShape weightSum) this + (array1FromList others) + (array1FromList points), + rotated (A.arrayShape canvas) rot mov pic) + (weightSum,canvas) + + +finalizeWeightedCanvas :: + (Channel Z Float, Channel DIM1 Float) -> Array DIM3 Word8 +finalizeWeightedCanvas = + Run.with CUDA.run1 $ + \(weightSum, canvas) -> + imageByteFromFloat $ interleaveChannels $ + A.zipWith (/) canvas $ + replicateChannel (A.indexTail $ A.indexTail $ A.shape canvas) weightSum + + +processOverlap :: + Option.Args -> + [(Float, Array DIM3 Word8)] -> + [((Int, (FilePath, ((Float, Float), Channel Z Float))), + (Int, (FilePath, ((Float, Float), Channel Z Float))))] -> + IO ([(Float, Float)], [((Float, Float), Array DIM3 Word8)]) +processOverlap args picAngles pairs = do + let opt = Option.option args + let info = CmdLine.info (Option.verbosity opt) + + let padSize = Option.padSize opt + let (maybeAllOverlapsShared, optimalOverlapShared) = + case Just $ Z :. padSize :. padSize of + Just padExtent -> + (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 + padExtent = Z :. padHeight :. padWidth + in (Just $ allOverlapsRun padExtent (Option.minimumOverlap opt), + optimalOverlap padExtent (Option.minimumOverlap opt)) + + displacements <- + fmap catMaybes $ + forM pairs $ \((ia,(pathA,(leftTopA,picA))), (ib,(pathB,(leftTopB,picB)))) -> do + forM_ maybeAllOverlapsShared $ \allOverlapsShared -> when False $ + writeGrey (Option.quality opt) + (printf "/tmp/%s-%s-score.jpeg" + (FilePath.takeBaseName pathA) (FilePath.takeBaseName pathB)) $ + allOverlapsShared picA picB + + let doffset@(dox,doy) = fst $ optimalOverlapShared picA picB + let diff = overlapDifferenceRun 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)) $ + composeOverlap 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) + + +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)] -> + [((Int, (FilePath, ((Float, Float), Channel Z Float))), + (Int, (FilePath, ((Float, Float), Channel Z Float))))] -> + IO ([(Float, Float)], [((Float, Float), Array DIM3 Word8)]) +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 + let optimalOverlapShared = + optimalOverlapBigMulti + (Z :. padSize :. padSize) + (Z :. stampSize :. stampSize) + (Option.numberStamps opt) + (Option.maximumDifference opt) + (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) + let correspondences = + map + (\(pa,pb,score) -> + (((ia, add pa leftTopA), (ib, add pb leftTopB)), score)) $ + optimalOverlapShared picA picB + info $ printf "left-top: %s, %s" (show leftTopA) (show leftTopB) + info $ printf "%s - %s" pathA pathB + forM_ correspondences $ \(((_ia,pa@(xa,ya)),(_ib,pb@(xb,yb))), score) -> + info $ + printf "%s ~ %s, (%f,%f), %f" + (show pa) (show pb) (xb-xa) (yb-ya) score + return $ map fst 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) + (HComplex.magnitude r) (HComplex.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 -> + (pairFromComplex $ + HComplex.cis angle * 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 + 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" + picAngles <- + forM paths $ \(imageOption, path) -> do + pic <- readImage (Option.verbosity opt) path + let maxAngle = Option.maximumAbsoluteAngle opt + let angles = + linearScale (Option.numberAngleSteps opt) + (-maxAngle, maxAngle) + when False $ analyseRotations angles pic + let angle = + maybe (findOptimalRotation angles pic) id $ + Option.angle imageOption + info $ printf "%s %f\176\n" path angle + return (path, (angle*pi/180, pic)) + + notice "\nfind relative placements" + let rotated = + map (mapSnd (prepareOverlapMatching (Option.smooth opt))) picAngles + let prepared = map (snd . snd) rotated + let pairs = do + (a:as) <- tails $ zip [0..] rotated + b <- as + return (a,b) + + when False $ do + notice "write fft" + let pic0 : pic1 : _ = prepared + size = (Z:.512:.1024 :: DIM2) + writeGrey (Option.quality opt) "/tmp/padded.jpeg" $ + CUDA.run1 + (imageByteFromFloat . + pad 0 (A.lift size)) $ + pic0 + writeGrey (Option.quality opt) "/tmp/spectrum.jpeg" $ + CUDA.run $ imageByteFromFloat $ A.map Complex.real $ + FFT.fft2D FFT.Forward $ + CUDA.run1 + (A.map (A.lift . (:+ 0)) . + pad 0 (A.lift size)) $ + pic0 + writeGrey (Option.quality opt) "/tmp/convolution.jpeg" $ + CUDA.run $ imageByteFromFloat $ A.map (0.000001*) $ + convolvePadded size (A.use pic0) (A.use pic1) + + (floatPoss, picRots) <- + (if Option.finetuneRotate opt + then processOverlapRotate + else processOverlap) + 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 + forM_ (Option.outputHard opt) $ \path -> + writeImage (Option.quality opt) path $ + finalizeCanvas $ + foldl + (\canvas (mov, rot, pic) -> updateCanvas rot mov pic canvas) + (emptyCanvas (Z :. 3 :. canvasHeight :. canvasWidth)) + movRotPics + + 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) + + forM_ (zip geometryRelations picAngles) $ + \((thisGeom, otherGeoms, allPoints), (path, _)) -> do + + let stem = FilePath.takeBaseName path + when False $ do + writeGrey (Option.quality opt) + (printf "/tmp/%s-distance-box.jpeg" stem) $ + distanceMapBoxRun canvasShape thisGeom + + writeGrey (Option.quality opt) + (printf "/tmp/%s-distance-contained.jpeg" stem) $ + distanceMapContainedRun canvasShape thisGeom otherGeoms + + writeGrey (Option.quality opt) + (printf "/tmp/%s-distance-points.jpeg" stem) $ + distanceMapPointsRun canvasShape allPoints + + forM_ (Option.outputDistanceMap opt) $ \format -> + writeGrey (Option.quality opt) (printf format stem) $ + distanceMapRun canvasShape thisGeom otherGeoms allPoints + + forM_ (Option.output opt) $ \path -> do + notice "\nweighted composition" + writeImage (Option.quality opt) path $ + finalizeWeightedCanvas $ + foldl + (\canvas ((thisGeom, otherGeoms, allPoints), (_rot, pic)) -> + updateWeightedCanvas (Option.distanceGamma opt) + thisGeom otherGeoms allPoints pic canvas) + (emptyWeightedCanvas (Z :. 3 :. canvasHeight :. canvasWidth)) + (zip geometryRelations picRots) + +main :: IO () +main = process =<< Option.get diff --git a/src/Draft.hs b/src/Draft.hs new file mode 100644 --- /dev/null +++ b/src/Draft.hs @@ -0,0 +1,164 @@ +module Main where + +import qualified Codec.Picture.Png as PNG +import qualified Codec.Picture as Pic + +import qualified Polygon +import qualified Line +import Polygon (Polygon(PolygonCW)) +import Line (Line2, Line(Segment)) +import Point2 (Point2(Point2), distance) + +import qualified Data.List.Match as Match +import Control.Monad (liftM2) +import Data.List.HT (removeEach) +import Data.Maybe.HT (toMaybe) +import Data.Maybe (catMaybes, mapMaybe) + + +project :: + (Fractional a, Ord a) => + Point2 a -> (Point2 a, Point2 a) -> Maybe (Point2 a) +project x (a, b) = + let (r, _, _, y) = Line.distanceAux a b x + in toMaybe (0<=r && r<=1) y + +segments :: (Eq a, Num a) => [Point2 a] -> [Line2 a] +segments = map (uncurry Segment) . Polygon.edges + +intersections :: + (Fractional a, Ord a) => + [Line2 a] -> [Line2 a] -> [Point2 a] +intersections segments0 segments1 = + catMaybes $ liftM2 Line.intersect segments0 segments1 + +averageRGB :: (Eq a, Fractional a) => [(a, (a,a,a))] -> (a,a,a) +averageRGB weightColors = + let (rs0, cs) = unzip weightColors + rs = if all (0==) rs0 then Match.replicate rs0 1 else rs0 + rsum = sum rs + scale r (red,green,blue) = (r*red, r*green, r*blue) + add (red0,green0,blue0) (red1,green1,blue1) = + (red0+red1, green0+green1, blue0+blue1) + in foldl1 add $ zipWith (\r c -> scale (r/rsum) c) rs cs + +polyMany :: + (Floating a, RealFrac a) => + [((a,a,a), [Point2 a])] -> Pic.Image Pic.PixelRGB8 +polyMany colorPointss = + let pointss = map snd colorPointss + polys = map PolygonCW pointss + intPointss = + map + (\(points, otherPointss) -> + intersections (segments points) + (concatMap segments otherPointss)) $ + removeEach pointss + render xi yi = + let x = fromIntegral xi + y = fromIntegral yi + p = Point2 (x,y) + minDist intPoints others this = + minimum $ map (distance p) $ + (intPoints ++) $ + filter (\q -> any (flip Polygon.contains q) others) $ + (this++) $ mapMaybe (project p) $ Polygon.edges this + minDists = + zipWith + (\intPoints (this, others) -> + minDist intPoints (map PolygonCW others) this) + intPointss (removeEach pointss) + contained = map (flip Polygon.contains p) polys + weightColors = + map snd $ + filter fst $ + zip contained $ + zip minDists (map fst colorPointss) + (red,green,blue) = + case weightColors of + [] -> (1,1,1) + [(_, color)] -> color + _ -> averageRGB weightColors + in Pic.PixelRGB8 + (round (red*255)) (round (green*255)) (round (blue*255)) + in Pic.generateImage render 256 256 + +poly2 :: + (Floating a, RealFrac a) => + [Point2 a] -> [Point2 a] -> Pic.Image Pic.PixelRGB8 +poly2 points0 points1 = + let poly0 = PolygonCW points0 + poly1 = PolygonCW points1 + intPoints = intersections (segments points0) (segments points1) + render xi yi = + let x = fromIntegral xi + y = fromIntegral yi + p = Point2 (x,y) + minDist other this = + minimum $ map (distance p) $ + (intPoints ++) $ filter (Polygon.contains other) $ + (this++) $ mapMaybe (project p) $ Polygon.edges this + minDist0 = minDist poly1 points0 + minDist1 = minDist poly0 points1 + scale = 255 / (minDist0 + minDist1) + in case (Polygon.contains poly0 p, Polygon.contains poly1 p) of + (False, False) -> Pic.PixelRGB8 255 255 255 + (False, True) -> Pic.PixelRGB8 0 0 255 + (True, False) -> Pic.PixelRGB8 255 0 0 + (True, True) -> + Pic.PixelRGB8 + (round (minDist0*scale)) 0 + (round (minDist1*scale)) + in Pic.generateImage render 256 256 + + +triangle :: Pic.Image Pic.PixelRGB8 +triangle = + let poly = PolygonCW [Point2 (10, 10), Point2 (200, 100), Point2 (100, 200)] + render xi yi = + let x, y :: Integer + x = fromIntegral xi + y = fromIntegral yi + in if Polygon.contains poly $ Point2 (x, y) + then Pic.PixelRGB8 0 0 255 + else Pic.PixelRGB8 255 255 255 + in Pic.generateImage render 256 256 + + +gradient :: Pic.Image Pic.PixelRGB8 +gradient = + Pic.generateImage + (\x y -> Pic.PixelRGB8 (fromIntegral x) (fromIntegral y) 128) + 256 256 + + +triA, triB :: [Point2 Double] +triA = [Point2 (10, 10), Point2 (245, 100), Point2 (100, 245)] +triB = map (\(Point2 (x,y)) -> Point2 (255-x,y)) triA + +isoBox :: (Eq a, Num a) => (a, a) -> (a, a) -> [Point2 a] +isoBox (l,t) (r,b) = + [Point2 (l,t), Point2 (r,t), Point2 (r,b), Point2 (l,b)] + +boxA, boxB :: [Point2 Double] +boxA = isoBox (10, 10) (200, 200) +boxB = isoBox (55, 55) (245, 245) + +buxA, buxB, buxC, buxD :: [Point2 Double] +buxA = isoBox ( 10, 10) (150, 100) +buxB = isoBox (110, 10) (250, 100) +buxC = isoBox ( 30, 50) (170, 150) +buxD = isoBox (110, 80) (249, 170) + +main :: IO () +main = do + PNG.writePng "/tmp/mixed.png" $ + polyMany [((1,0,0), boxA), ((1,1,0), boxB), + ((0,1,0), triA), ((0,0,1), triB)] + PNG.writePng "/tmp/boxes.png" $ + polyMany [((1,0,0), buxA), ((1,1,0), buxB), + ((0,1,0), buxC), ((0,0,1), buxD)] + PNG.writePng "/tmp/box2.png" $ poly2 boxA boxB + PNG.writePng "/tmp/triangle2.png" $ poly2 triA triB + PNG.writePng "/tmp/triangle.png" triangle + PNG.writePng "/tmp/test.png" gradient diff --git a/src/Option.hs b/src/Option.hs new file mode 100644 --- /dev/null +++ b/src/Option.hs @@ -0,0 +1,246 @@ +module Option where + +import Option.Utility (exitFailureMsg, parseNumber, fmapOptDescr) + +import qualified System.Console.GetOpt as Opt +import qualified System.Environment as Env +import System.Console.GetOpt (ArgDescr(NoArg, ReqArg), getOpt, usageInfo) + +import qualified System.Exit as Exit + +import Control.Monad (when) + +import qualified Distribution.Verbosity as Verbosity +import qualified Distribution.ReadE as ReadE +import Distribution.Verbosity (Verbosity) + +import Text.Printf (printf) + + +data Args = + Args { + option :: Option, + inputs :: [(Image, FilePath)] + } + +defltArgs :: Args +defltArgs = Args {option = defltOption, inputs = []} + + +data Option = + Option { + 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" + quality :: Int, + maximumAbsoluteAngle :: Float, + numberAngleSteps :: Int, + smooth :: Int, + padSize :: Int, + minimumOverlap :: Float, + maximumDifference :: Float, + finetuneRotate :: Bool, + numberStamps :: Int, + stampSize :: Int, + distanceGamma :: Float + } + +defltOption :: Option +defltOption = + Option { + verbosity = Verbosity.verbose, + output = Nothing, + outputHard = Nothing, + outputOverlap = Nothing, + outputDistanceMap = Nothing, + quality = 99, + maximumAbsoluteAngle = 1, + numberAngleSteps = 40, + smooth = 20, + padSize = 1024, + minimumOverlap = 1/4, + maximumDifference = 0.2, + finetuneRotate = False, + numberStamps = 5, + stampSize = 64, + distanceGamma = 2 + } + + +data Image = + Image { + angle :: Maybe Float + } + deriving (Eq) + +defltImage :: Image +defltImage = Image {angle = Nothing} + + +type Description a = [Opt.OptDescr (a -> IO a)] + +{- +Guide for common Linux/Unix command-line options: + http://www.faqs.org/docs/artu/ch10s05.html +-} +optionDescription :: Description a -> Description Option +optionDescription desc = + Opt.Option ['h'] ["help"] + (NoArg $ \ _flags -> do + programName <- Env.getProgName + putStrLn $ + usageInfo + ("Usage: " ++ programName ++ + " [OPTIONS]... [[INPUTOPTIONS]... INPUT]...") $ + desc + Exit.exitSuccess) + "show options" : + + Opt.Option ['v'] ["verbose"] + (flip ReqArg "N" $ \str flags -> do + case ReadE.runReadE Verbosity.flagToVerbosity str of + Right n -> return (flags{verbosity = n}) + Left msg -> exitFailureMsg msg) + (printf "verbosity level: 0..3, default: %d" + (fromEnum $ verbosity defltOption)) : + + Opt.Option [] ["output"] + (flip ReqArg "PATH" $ \str flags -> + return $ flags{output = Just str}) + ("path to generated collage") : + + Opt.Option [] ["output-hard"] + (flip ReqArg "PATH" $ \str flags -> + return $ flags{outputHard = Just str}) + ("path to collage without fading") : + + Opt.Option [] ["output-overlap"] + (flip ReqArg "FORMAT" $ \str flags -> + return $ flags{outputOverlap = Just str}) + ("path format for overlapped pairs") : + + Opt.Option [] ["output-distance-map"] + (flip ReqArg "FORMAT" $ \str flags -> + return $ flags{outputDistanceMap = Just str}) + ("path format for distance maps") : + + Opt.Option [] ["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" + (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" + (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 [] ["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"] + (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"] + (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"] + (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"] + (NoArg $ \flags -> return $ flags{finetuneRotate = True}) + (printf "Fine-tune rotation together with overlapping, default: disabled") : + + Opt.Option [] ["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"] + (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"] + (flip ReqArg "FRACTION" $ \str flags -> + fmap (\x -> flags{distanceGamma = x}) $ + parseNumber "gamma exponent" (0<) "positive" str) + (printf "Distance exponent, default: %f" + (distanceGamma defltOption)) : + + [] + + +description :: Description (Image, Args) -> Description (Image, Args) +description desc = + map + (fmapOptDescr $ \update (image, old) -> do + new <- update $ option old + return (image, old {option = new})) + (optionDescription desc) + ++ + + Opt.Option [] ["hint-angle"] + (flip ReqArg "DEGREE" $ \str (image, args) -> + fmap (\x -> (image{angle = Just x}, args)) $ + parseNumber "angle" (\w -> -1000<=w && w<=1000) "degree" str) + (printf "Angle of the next image in first phase, default: %s" $ + maybe "automatic estimation" show (angle defltImage)) : + + [] + + +addFile :: FilePath -> ((Image, Args) -> IO (Image, Args)) +addFile path (image, args) = + return (defltImage, args {inputs = (image,path) : inputs args}) + + + +get :: IO Args +get = do + let desc = description desc + argv <- Env.getArgs + let (args, _files, errors) = getOpt (Opt.ReturnInOrder addFile) desc argv + when (not $ null errors) $ + exitFailureMsg (init (concat errors)) + + (lastImage, parsedArgs) <- foldl (>>=) (return (defltImage, defltArgs)) args + + when (lastImage /= defltImage) $ + exitFailureMsg "unused trailing image options" + + case inputs parsedArgs of + [] -> exitFailureMsg "no input files" + images -> return $ parsedArgs {inputs = reverse images} diff --git a/src/Option/Utility.hs b/src/Option/Utility.hs new file mode 100644 --- /dev/null +++ b/src/Option/Utility.hs @@ -0,0 +1,34 @@ +module Option.Utility where + +import qualified System.Console.GetOpt as G +import qualified System.Exit as Exit +import qualified System.IO as IO + + +parseNumber :: + (Read a) => + String -> (a -> Bool) -> String -> String -> IO a +parseNumber name constraint constraintName str = + case reads str of + [(n, "")] -> + if constraint n + then return n + else exitFailureMsg $ name ++ " must be a " ++ constraintName ++ " number" + _ -> + exitFailureMsg $ name ++ " must be a number, but is '" ++ str ++ "'" + +exitFailureMsg :: String -> IO a +exitFailureMsg msg = do + IO.hPutStrLn IO.stderr msg + Exit.exitFailure + +fmapArgDescr :: (a -> b) -> (G.ArgDescr a -> G.ArgDescr b) +fmapArgDescr f d = + case d of + G.NoArg a -> G.NoArg $ f a + G.ReqArg g str -> G.ReqArg (f.g) str + G.OptArg g str -> G.OptArg (f.g) str + +fmapOptDescr :: (a -> b) -> (G.OptDescr a -> G.OptDescr b) +fmapOptDescr f (G.Option short long arg help) = + G.Option short long (fmapArgDescr f arg) help