diff --git a/Changes.md b/Changes.md --- a/Changes.md +++ b/Changes.md @@ -1,5 +1,10 @@ # Change log for the `patch-image` package +## 0.3: + + * Allow to save program state, modify it manually + and re-run with adapted parameters using CSV files. + ## 0.2: * Add new executable that is based on LLVM and `knead`. diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,3 +1,5 @@ +# Automatic mode + This is the workflow: Scan parts of an image that considerably overlap. They must all be approximately oriented correctly. @@ -97,7 +99,7 @@ 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 +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, @@ -124,7 +126,7 @@ ## 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 +These 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. @@ -158,8 +160,9 @@ 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. +Our LLVM implementation provides an additional way to assemble the image parts. +The already known 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. @@ -194,3 +197,187 @@ * `--quality`: JPEG quality percentage for writing the images. + + +# Semi-Automatic mode: inspection and correction + +If the program does not work correctly on a particular set of images +or if you want to archive what it has done, +then the following options may help you. + +* `--output-state=state-%s.csv`: + Use of this option let the program generate three files + that correspond to the stages of the image processing. + For each file the program replaces the placeholder `%s` + by an identifier for the particular stage. + For the example above the generated files are: + + * `state-angle.csv`: + After the rotation detection + the program emits this file containing the image paths + and their angles in degrees. + + * `state-relation.csv`: + This file contains for all possible pairs of image parts + whether the program found that they overlap and, if yes, + what their relative offset is. + + * `state-position.csv`: + From the relative offsets the program + computes the absolute part positions. + This CSV file is essentially an extension of `state-angle.csv` + by the absolute positions. + However, in `--finetune-rotate` mode + the angles may differ from the ones in `state-angle.csv`. + +* `--state=state-position.csv`: + You can feed the program with various parameters + such that the program does not need to compute them itself. + This can be useful for repeating a previous computation, + speedup a computation or improve a previous computation. + The idea is to take either `state-angle.csv` or `state-position.csv` + as generated by `--output-state`. + You may omit any column or re-order them, or you may clear individual cells. + At least the Image column must be present. + The most simple CSV file would thus be a text file + where the first line is `Image` and + each of the following lines contains one image path. + The images listed in the CSV file are loaded + additionally to the files given at the command-line. + A complete CSV file would contain the columns `Image`, `Angle`, `X`, `Y`, + and in `--finetune-rotate` mode additionally `DAngle`. + + A missing column or empty cell means + that the program computes the parameter itself. + You may even fix say an X coordinate + but let the program compute the corresponding Y coordinate. + I do not know whether there is a useful application of this feature, + but the feature is a natural result of how the computation works. + + You may use the CSV file to set the rotation angle of the assembled image. + To this end, choose an image part where you know the angle for sure + and add the `Angle` to the image's row + while omitting the `Angle` for all other images. + + If an image is grossly turned away from horizontal alignment, + that is, rotated by, say, more than 2 degree, + or the image does not contain a horizontal structure to align to, + then the automatic angle detection will fail. + In this case you might use the Angle column to specify the angle directly. + + In `--finetune-rotate` mode the `Angle` column means the following: + If an `Angle` cell is set + then the rotation angle is fixed both in the first rotation detection phase, + and in the subsequent rotation correction. + If an `Angle` cell is empty + then the rotation angle is both determined by horizontal structures + in the first phase + and then adjusted in the subsequent correction phase. + Further, if you want to say that all angles are only hints + that may be adjusted in the second phase, + then you should add an empty `DAngle` column. + If you want to exclude some angles from correction, + then enter a `0` at the corresponding cells in the `DAngle` column. + + Summarized: + If both `Angle` and `DAngle` columns are present, + then the `Angle` column specifies the angle as determined in the first phase + and the `DAngle` column specifies the angle difference + as determined in the second phase. + You may enter values different from `0` in the `DAngle` column. + E.g., setting `1` means that any angle determined in the first phase + is unconditionally increased by 1 degree in the second phase. + I do not know, whether this can ever be useful, + but providing the columns `Angle` and `DAngle` this way + seemed to be the most natural way to support all sensible use cases. + If only the `Angle` column is available + this is equivalent to adding a `DAngle` column + and set `DAngle=0` if `Angle` is set + and leave `DAngle` empty if `Angle` is empty. + + If you list images both on the command-line and in a CSV file, + then both lists are concatenated. + I recommend that you use either command-line arguments or a CSV file + to specify the image list. + You can easily get this wrong + when running the program with a command-line file list first + and then rerunning the computation starting with a saved `state.csv` file. + +* `--relations=state-relation.csv`: + The table `relation.csv` contains pairs of images + and their relative offsets. + The most simple use case is to take the `state-relation.csv` file + as generated by `--output-state` + and alter or delete rows that you find inappropriate. + The common columns are `ImageA`, `ImageB`, and `Rel`. + `Rel` can be: + + * empty: + The user wants the program to determine whether the image parts overlap. + + * `X`: + The user knows that the image parts overlap. + + * `-`: + The user knows that the image parts do not overlap. + + Without `--finetune-rotate` the columns `DX` and `DY` follow. + They tell the offset of the image parts to each other. + In `--finetune-rotate` mode + there are the additional columns `XA`, `YA`, `XB`, `YB`. + They contain coordinates of corresponding points in the image pairs. + There are usually many such correspondences for one image pair + and they are listed below header lines for the image pairs. + + It is uncommon to have both a `state.csv` file with absolute positions + and a `relation.csv` file, + though the program supports it. + If there is a pair of image parts + where both images have fixed positions in `state.csv` + and there is also a fixed offset registered in `relation.csv`, + then the fixed positions win. + The differing offset between the image parts + are seen as an unfortunate impossibility + to optimally match the optimization goals. + + There is another difficulty in mixing fixed positions + together with fixed offsets: + For a pair of image parts with fixed positions + the program does not match images + in order to compute the relative offsets of the parts + because it would not affect the final positions, anyway. + This optimization makes the program skip the phase + of determining image relations completely, + if all part positions are fixed. + This is the wanted effect. + The unwanted effect is that since no relative offsets are computed, + no relative offsets can be written to `state-relation.csv` + as requested by `--output-state`. + That is, re-running the program with `--state` using a CSV file + as previously emitted by `--output-state` + while using `--output-state` anew + will generate a new `state-relation.csv` file + with less information than before or no information at all. + + Files listed in `relation.csv` are not loaded. + Instead of this, the program checks consistency + with `state.csv` and files listed as command-line arguments. + + +# Trouble shooting + +## Unrelated images are recognized as overlapping or vice versa + +If the program misses related image pairs +or recognizes unrelated images as overlapping, +you may first watch the overlap differences the program prints. +You may then adapt the threshold via the `--maximum-difference` option. + +The total image might contain repetitive elements +that make unrelated image parts look like they overlap. +So, if there is no clear threshold +or you find the above procedure inappropriate, +you may alternatively let the program emit a `relation.csv` file +using the `--output-state` option. +You can then edit `relation.csv` file +and re-run the program with the option `--relations=relation.csv`. diff --git a/patch-image.cabal b/patch-image.cabal --- a/patch-image.cabal +++ b/patch-image.cabal @@ -1,5 +1,5 @@ Name: patch-image -Version: 0.2 +Version: 0.3 License: BSD3 License-File: LICENSE Author: Henning Thielemann @@ -16,6 +16,16 @@ It solves the problem that your scanner may be too small to capture a certain image as a whole. . + There are two versions of the program: + . + * A version using LLVM via the @knead@ package: + Enable it with the Cabal flag @-fllvm@. + . + * A version using CUDA via @accelerate-cuda@: + Enable it with the Cabal flag @-fcuda@. + . + Both of them have their individual installation difficulties. + . Restrictions: . * Only supports JPEG format. @@ -31,7 +41,7 @@ README.md Source-Repository this - Tag: 0.2 + Tag: 0.3 Type: darcs Location: http://hub.darcs.net/thielema/patch-image/ @@ -58,8 +68,10 @@ Other-Modules: MatchImageBorders KneadShape + State LinearAlgebra Arithmetic + Degree Option.Utility Option Hs-Source-Dirs: src @@ -81,10 +93,15 @@ containers >=0.4.2 && <0.6, fft >=0.1.7 && <0.2, carray >=0.1.5 && <0.2, - array >=0.5 && <0.6, + array >=0.4 && <0.6, + cassava >=0.4.5 && <0.5, + unordered-containers >=0.2.5 && <0.2.9, + bytestring >=0.9.2 && <0.11, + explicit-exception >=0.1.7 && <0.2, Cabal >=1.18 && <2, filepath >=1.3 && <1.4, - utility-ht >=0.0.1 && <0.1, + non-empty >=0.2 && <0.4, + utility-ht >=0.0.13 && <0.1, base >=4 && <5 Else Buildable: False @@ -92,8 +109,10 @@ Executable patch-image-cuda Main-Is: Accelerate.hs Other-Modules: + State LinearAlgebra Arithmetic + Degree Option.Utility Option Hs-Source-Dirs: src diff --git a/src/Accelerate.hs b/src/Accelerate.hs --- a/src/Accelerate.hs +++ b/src/Accelerate.hs @@ -2,11 +2,15 @@ module Main where import qualified Option +import qualified State import qualified Arithmetic as Arith +import qualified Degree import LinearAlgebra ( - absolutePositionsFromPairDisplacements, layoutFromPairDisplacements, + absolutePositionsFromPairDisplacements, fixAtLeastOnePosition, + layoutFromPairDisplacements, fixAtLeastOneAnglePosition, ) +import Degree (Degree(Degree), getDegree) import Arithmetic ( Point2, rotateStretchMovePoint, @@ -19,13 +23,13 @@ distance, linearScale, divUp, - pairFromComplex, - mapComplex, + guardedPairs, + maximum0, ) import qualified Data.Array.Accelerate.Fourier.Real as FourierReal import qualified Data.Array.Accelerate.CUFFT.Single as CUFFT -import qualified Data.Array.Accelerate.Data.Complex as Complex +import qualified Data.Array.Accelerate.Data.Complex as AComplex import qualified Data.Array.Accelerate.CUDA.Foreign as CUDAForeign import qualified Data.Array.Accelerate.CUDA as CUDA import qualified Data.Array.Accelerate.IO as AIO @@ -50,7 +54,7 @@ import qualified Graphics.Gnuplot.Plot.TwoDimensional as Plot2D import qualified Graphics.Gnuplot.Graph.TwoDimensional as Graph2D -import qualified Data.Complex as HComplex +import qualified Data.Complex as Complex import qualified Codec.Picture as Pic @@ -60,15 +64,21 @@ 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 Data.List.Key as Key import Control.Monad.HT (void) -import Control.Monad (liftM2, when) +import Control.Monad (liftM2, when, join) +import Control.Applicative ((<$>)) + +import qualified Data.Foldable as Fold +import qualified Data.List.Key as Key +import qualified Data.Map as Map +import Data.Monoid ((<>)) import Data.Maybe.HT (toMaybe) -import Data.Maybe (catMaybes) -import Data.List.HT (mapAdjacent, tails) +import Data.Maybe (mapMaybe, isNothing) +import Data.List.HT (mapAdjacent) import Data.Traversable (forM) import Data.Foldable (forM_, foldMap) import Data.Tuple.HT (mapPair, mapFst, mapSnd, mapThd3) @@ -77,7 +87,9 @@ import System.IO.Unsafe (unsafePerformIO) -readImage :: Verbosity -> FilePath -> IO (Array DIM3 Word8) +type ColorImage8 = Array DIM3 Word8 + +readImage :: Verbosity -> FilePath -> IO ColorImage8 readImage verbosity path = do epic <- Pic.readImage path case epic of @@ -97,7 +109,7 @@ ((), dat) _ -> ioError $ userError "unsupported image type" -writeImage :: Int -> FilePath -> Array DIM3 Word8 -> IO () +writeImage :: Int -> FilePath -> ColorImage8 -> IO () writeImage quality path arr = do let (Z :. height :. width :. 3) = A.arrayShape arr Pic.saveJpgImage quality path $ Pic.ImageYCbCr8 $ @@ -117,7 +129,7 @@ Pic.imageData = snd $ AIO.toVectors arr } -colorImageExtent :: Array DIM3 Word8 -> (Int, Int) +colorImageExtent :: ColorImage8 -> (Int, Int) colorImageExtent pic = case A.arrayShape pic of Z:.height:.width:._chans -> (width, height) @@ -174,7 +186,8 @@ -type Channel ix a = Array (ix :. Int :. Int) a +type Channel ix = Array (ix :. Int :. Int) +type Plane = Channel Z type ExpDIM2 ix = Exp ix :. Exp Int :. Exp Int type ExpDIM3 ix = Exp ix :. Exp Int :. Exp Int :. Exp Int @@ -304,7 +317,7 @@ rotateHistogram :: - Float -> Array DIM3 Word8 -> (Array DIM3 Word8, Array DIM1 Float) + Degree Float -> ColorImage8 -> (ColorImage8, Array DIM1 Float) rotateHistogram = let rot = Run.with CUDA.run1 $ \orient arr -> @@ -313,15 +326,15 @@ separateChannels $ imageFloatFromByte arr in (imageByteFromFloat $ interleaveChannels rotated, rowHistogram rotated) - in \angle arr -> rot (cos angle, sin angle) arr + in \angle arr -> rot (Degree.cis angle) arr -analyseRotations :: [Float] -> Array DIM3 Word8 -> IO () +analyseRotations :: [Degree Float] -> ColorImage8 -> 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 + let (rotated, histogram) = rotateHistogram degree pic + let stem = printf "rotated%+07.2f" $ getDegree degree writeImage 90 ("/tmp/" ++ stem ++ ".jpeg") rotated let diffHistogram = map abs $ mapAdjacent (-) $ A.toList histogram printf "%s: maxdiff %8.3f, sqrdiff %8.0f\n" @@ -351,24 +364,24 @@ in A.generate (A.index1 (size-1)) $ \i -> arr ! (A.index1 $ A.unindex1 i + 1) - arr ! i -scoreRotation :: Float -> Array DIM3 Word8 -> Float +scoreRotation :: Degree Float -> ColorImage8 -> 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 + in \angle arr -> Acc.the $ rot (Degree.cis angle) arr -findOptimalRotation :: [Float] -> Array DIM3 Word8 -> Float +findOptimalRotation :: [Degree Float] -> ColorImage8 -> Degree Float findOptimalRotation angles pic = - Key.maximum (flip scoreRotation pic . (* (pi/180))) angles + Key.maximum (flip scoreRotation pic) angles magnitudeSqr :: (A.Elt a, A.IsNum a) => Exp (Complex a) -> Exp a magnitudeSqr = Exp.modify (expr:+expr) $ \(r:+i) -> r*r+i*i -fourierTransformationRun :: Array DIM3 Word8 -> IO (Array DIM2 Word8) +fourierTransformationRun :: ColorImage8 -> IO (Array DIM2 Word8) fourierTransformationRun pic = do let (shape@(Z:.height:.width):._) = A.arrayShape pic plan <- @@ -393,7 +406,7 @@ imageFloatFromByte arr return $ trans pic -fourierTransformation :: Option.Option -> FilePath -> Array DIM3 Word8 -> IO () +fourierTransformation :: Option.Option -> FilePath -> ColorImage8 -> IO () fourierTransformation opt path pic = do let stem = FilePath.takeBaseName path spec <- fourierTransformationRun pic @@ -435,17 +448,18 @@ z1 = weighted ! A.lift (Z :. y :. mod (xi+1) width) in linearIp (z0,z1) frac -radonAngle :: (Float, Float) -> Array DIM3 Word8 -> IO Float +radonAngle :: + (Degree Float, Degree Float) -> ColorImage8 -> IO (Degree Float) radonAngle (minAngle,maxAngle) pic = do let (shape@(Z :. height :. _width):._) = A.arrayShape pic plan <- CUDAForeign.inDefaultContext $ CUFFT.plan2D CUFFT.forwardReal shape let height2 = fromIntegral (div height 2) - let slope w = tan (w*pi/180) * height2 + let slope w = tan (Degree.toRadian w) * height2 let minX = floor $ slope minAngle let maxX = ceiling $ slope maxAngle - let angle s = atan (s/height2) * 180/pi + let angle s = Degree.fromRadian $ atan (s/height2) let trans = Run.with CUDA.run1 $ \arr -> A.map A.snd $ argmaximum $ @@ -457,16 +471,16 @@ return $ angle $ fromIntegral $ Acc.the (trans pic) + minX -rotateManifest :: Float -> Array DIM3 Word8 -> Array DIM3 Float +rotateManifest :: Degree Float -> ColorImage8 -> 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 + in \angle arr -> rot (Degree.cis angle) arr prepareOverlapMatching :: - Int -> (Float, Array DIM3 Word8) -> ((Float,Float), Channel Z Float) + Int -> (Degree Float, ColorImage8) -> ((Float,Float), Plane Float) prepareOverlapMatching = let rot = Run.with CUDA.run1 $ \radius orient arr -> @@ -477,7 +491,7 @@ brightnessPlane $ separateChannels $ imageFloatFromByte arr in \radius (angle, arr) -> mapFst (mapPair (Acc.the, Acc.the)) $ - rot radius (cos angle, sin angle) arr + rot radius (Degree.cis angle) arr ceilingPow2 :: Exp Int -> Exp Int @@ -498,7 +512,7 @@ mulConj :: (A.Elt a, A.IsFloating a) => Exp (Complex a) -> Exp (Complex a) -> Exp (Complex a) -mulConj x y = x * Complex.conjugate y +mulConj x y = x * AComplex.conjugate y fft2DGen :: @@ -568,7 +582,7 @@ A.stencil (\(a,m,b) -> smooth3 (smooth3 a, smooth3 m, smooth3 b)) A.Clamp highpass count arr = - A.zipWith (-) arr $ lowpass count arr + A.zipWith (-) arr $ lowpass count arr correlatePaddedSimple :: @@ -677,7 +691,7 @@ allOverlapsRun :: - DIM2 -> Float -> Channel Z Float -> Channel Z Float -> Channel Z Word8 + DIM2 -> Float -> Plane Float -> Plane Float -> Plane Word8 allOverlapsRun padExtent = Run.with CUDA.run1 $ \minOverlap picA picB -> imageByteFromFloat $ @@ -686,7 +700,7 @@ A.map A.fst $ allOverlaps padExtent minOverlap picA picB optimalOverlap :: - DIM2 -> Float -> Channel Z Float -> Channel Z Float -> (Float, (Int, Int)) + DIM2 -> Float -> Plane Float -> Plane Float -> (Float, (Int, Int)) optimalOverlap padExtent = let run = Run.with CUDA.run1 $ \minimumOverlap a b -> @@ -721,7 +735,7 @@ Reduce image sizes below the padExtent before matching images. -} optimalOverlapBig :: - DIM2 -> Float -> Channel Z Float -> Channel Z Float -> (Float, (Int, Int)) + DIM2 -> Float -> Plane Float -> Plane Float -> (Float, (Int, Int)) optimalOverlapBig padExtent = let run = Run.with CUDA.run1 $ \minimumOverlap a b -> @@ -771,7 +785,7 @@ using a part in the overlapping area. -} optimalOverlapBigFine :: - DIM2 -> Float -> Channel Z Float -> Channel Z Float -> (Float, (Int, Int)) + DIM2 -> Float -> Plane Float -> Plane Float -> (Float, (Int, Int)) optimalOverlapBigFine padExtent@(Z:.heightPad:.widthPad) = let overlaps = allOverlaps padExtent run = @@ -813,7 +827,7 @@ -} optimalOverlapBigMulti :: DIM2 -> DIM2 -> Int -> - Float -> Float -> Channel Z Float -> Channel Z Float -> + Float -> Maybe Float -> Plane Float -> Plane Float -> [(Float, (Int, Int), (Int, Int))] optimalOverlapBigMulti padExtent (Z:.heightStamp:.widthStamp) numCorrs = let overlapShrunk = @@ -846,18 +860,22 @@ (clip anchorA extent a) (clip anchorB extent b) - in \maximumDiff minimumOverlap a b -> + in \minimumOverlap mMaximumDiff a b -> let factors@(Z:.yk:.xk) = shrinkFactors padExtent (A.arrayShape a) (A.arrayShape b) - (_score, (shrunkdx, shrunkdy)) = + (_score, shrunkd@(shrunkdx, shrunkdy)) = 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 + doesOverlap = + case mMaximumDiff of + Just maximumDiff -> + maximumDiff > Acc.the (diffShrunk shrunkd factors a b) + Nothing -> True ((leftOverlap, topOverlap), (rightOverlap, bottomOverlap), @@ -867,7 +885,7 @@ widthStampClip = min widthOverlap widthStamp heightStampClip = min heightOverlap heightStamp - in (if diff < maximumDiff then id else const []) $ + in (if doesOverlap then id else const []) $ map (\(x,y) -> Acc.the $ @@ -909,7 +927,7 @@ overlapDifferenceRun :: (Int, Int) -> - Channel Z Float -> Channel Z Float -> Float + Plane Float -> Plane Float -> Float overlapDifferenceRun = let diff = Run.with CUDA.run1 overlapDifference in \d a b -> Acc.the $ diff d a b @@ -941,22 +959,22 @@ composeOverlap :: (Int, Int) -> - ((Float, Array DIM3 Word8), (Float, Array DIM3 Word8)) -> - Array DIM3 Word8 + ((Degree Float, ColorImage8), (Degree Float, ColorImage8)) -> + ColorImage8 composeOverlap = - let rot (angle,pic) = - rotate (cos angle, sin angle) $ - separateChannels $ imageFloatFromByte pic - in Run.with CUDA.run1 $ + let rotat (rot,pic) = + rotate rot $ separateChannels $ imageFloatFromByte pic + in (\f d (a,b) -> f d (mapFst Degree.cis a, mapFst Degree.cis b)) $ + Run.with CUDA.run1 $ \(dx,dy) (a,b) -> imageByteFromFloat $ interleaveChannels $ - overlap2 (dx, dy) (rot a, rot b) + overlap2 (dx, dy) (rotat a, rotat b) emptyCountCanvas :: (A.Slice ix, A.Shape ix) => ix :. Int :. Int -> - (Channel Z Int, Channel ix Float) + (Plane Int, Channel ix Float) emptyCountCanvas = Run.with CUDA.run1 $ \sh -> let (_ix :. height :. width) = unliftDim2 sh @@ -966,9 +984,9 @@ 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)) + (Acc (Plane Bool), Acc (Channel ix a)) -> + (Acc (Plane Int), Acc (Channel ix a)) -> + (Acc (Plane Int), Acc (Channel ix a)) addToCountCanvas (mask, pic) (count, canvas) = (A.zipWith (+) (A.map A.boolToInt mask) count, A.zipWith (+) canvas $ A.zipWith (*) pic $ @@ -976,9 +994,9 @@ A.map (A.fromIntegral . A.boolToInt) mask) updateCountCanvas :: - ((Float,Float), (Float,Float), Array DIM3 Word8) -> - (Channel Z Int, Channel DIM1 Float) -> - (Channel Z Int, Channel DIM1 Float) + ((Float,Float), (Float,Float), ColorImage8) -> + (Plane Int, Channel DIM1 Float) -> + (Plane Int, Channel DIM1 Float) updateCountCanvas = Run.with CUDA.run1 $ \(rot, mov, pic) (count,canvas) -> @@ -987,7 +1005,7 @@ separateChannels $ imageFloatFromByte pic) (count,canvas) -finalizeCountCanvas :: (Channel Z Int, Channel DIM1 Float) -> Array DIM3 Word8 +finalizeCountCanvas :: (Plane Int, Channel DIM1 Float) -> ColorImage8 finalizeCountCanvas = Run.with CUDA.run1 $ \(count, canvas) -> @@ -1043,7 +1061,7 @@ edges distanceMapEdgesRun :: - DIM2 -> Array DIM1 ((Float,Float),(Float,Float)) -> Channel Z Word8 + DIM2 -> Array DIM1 ((Float,Float),(Float,Float)) -> Plane Word8 distanceMapEdgesRun = Run.with CUDA.run1 $ \sh -> imageByteFromFloat . A.map (0.01*) . distanceMapEdges sh @@ -1100,20 +1118,15 @@ (A.use $ A.fromList (Z:.(4::Int)) $ liftM2 (,) [False,True] [False,True]) -distanceMapBoxRun :: - DIM2 -> ((Float,Float),(Float,Float),(Int,Int)) -> Channel Z Word8 +distanceMapBoxRun :: DIM2 -> Geometry Float -> Plane 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 (expr,expr) $ - \(valid, dist) -> valid ? (scale*dist, 0)) $ - maskedMinimum $ - A.map (Exp.mapSnd A.fst) $ - separateDistanceMap $ - distanceMapBox sh geom + scaleDistanceMapGeom geom $ + A.map (Exp.modify (expr,expr) $ \(valid, dist) -> valid ? (dist, 0)) $ + maskedMinimum $ + A.map (Exp.mapSnd A.fst) $ + separateDistanceMap $ + distanceMapBox sh geom -- maybe move to Accelerate.Utility @@ -1152,14 +1165,9 @@ Acc (Array DIM1 (Geometry a)) -> 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 (expr,expr) $ - \(valid, dist) -> valid ? (dist, 0)) $ + let distMap = separateDistanceMap $ distanceMapBox sh this + contained = containedAnywhere others $ A.map (A.snd . A.snd) distMap + in A.map (Exp.modify (expr,expr) $ \(valid, dist) -> valid ? (dist, 0)) $ maskedMinimum $ A.zipWith (Exp.modify2 expr (expr,(expr,expr)) $ \c (b,(dist,_)) -> @@ -1167,20 +1175,21 @@ contained distMap distanceMapContainedRun :: - DIM2 -> Geometry Float -> [Geometry Float] -> Channel Z Word8 + DIM2 -> Geometry Float -> [Geometry Float] -> Plane 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 + \sh this -> scaleDistanceMapGeom this . distanceMapContained sh this + in \sh this others -> distances sh this $ array1FromList others +scaleDistanceMapGeom :: + (A.IsFloating a, A.Elt a, A.Elt b, A.Shape ix) => + Exp (Geometry b) -> Acc (Array ix a) -> Acc (Array ix Word8) +scaleDistanceMapGeom this = + let scale = (4/) $ A.fromIntegral $ A.uncurry min $ Exp.thd3 this + in imageByteFromFloat . A.map (scale*) + pixelCoordinates :: (A.Elt a, A.IsFloating a) => Exp DIM2 -> Acc (Channel Z (a,a)) @@ -1199,23 +1208,22 @@ (Exp.modify2 (expr,expr) (expr,expr) distance) a b -distanceMapPointsRun :: - DIM2 -> - [Point2 Float] -> - Channel Z Word8 +distanceMapPointsRun :: DIM2 -> [Point2 Float] -> Plane Word8 distanceMapPointsRun = let distances = Run.with CUDA.run1 $ - \sh points -> - let scale = - case Exp.unlift (expr:.expr:.expr) sh of - _z:.y:.x -> (4/) $ A.fromIntegral $ min x y - in imageByteFromFloat $ A.map (scale*) $ - distanceMapPoints (pixelCoordinates sh) points - in \sh points -> - distances sh $ array1FromList points + \sh -> scaleDistanceMap . distanceMapPoints (pixelCoordinates sh) + in \sh points -> distances sh $ array1FromList points +scaleDistanceMap :: + (A.Elt a, A.IsFloating a) => Acc (Channel Z a) -> Acc (Channel Z Word8) +scaleDistanceMap arr = + let scale = + case Exp.unlift (expr:.expr:.expr) (A.shape arr) of + _z:.y:.x -> 4 / A.fromIntegral (min x y) + in imageByteFromFloat $ A.map (scale*) arr + {- | For every pixel it computes the distance to the closest point on the image part boundary @@ -1248,16 +1256,11 @@ Geometry Float -> [Geometry Float] -> [Point2 Float] -> - Channel Z Word8 + Plane Word8 distanceMapRun = let distances = Run.with CUDA.run1 $ - \sh this others points -> - let scale = - case Exp.unlift (expr:.expr:.expr) sh of - _z:.y:.x -> (4/) $ A.fromIntegral $ min x y - in imageByteFromFloat $ A.map (scale*) $ - distanceMap sh this others points + \sh this others -> scaleDistanceMap . distanceMap sh this others in \sh this others points -> distances sh this (array1FromList others) @@ -1279,7 +1282,7 @@ emptyWeightedCanvas :: (A.Slice ix, A.Shape ix) => ix :. Int :. Int -> - (Channel Z Float, Channel ix Float) + (Plane Float, Channel ix Float) emptyWeightedCanvas = Run.with CUDA.run1 $ \sh -> let (_ix :. height :. width) = unliftDim2 sh @@ -1302,9 +1305,9 @@ Geometry Float -> [Geometry Float] -> [Point2 Float] -> - Array DIM3 Word8 -> - (Channel Z Float, Channel DIM1 Float) -> - (Channel Z Float, Channel DIM1 Float) + ColorImage8 -> + (Plane Float, Channel DIM1 Float) -> + (Plane Float, Channel DIM1 Float) updateWeightedCanvasMerged = let update = Run.with CUDA.run1 $ @@ -1327,9 +1330,9 @@ Geometry Float -> [Geometry Float] -> [Point2 Float] -> - Array DIM3 Word8 -> - (Channel Z Float, Channel DIM1 Float) -> - (Channel Z Float, Channel DIM1 Float) + ColorImage8 -> + (Plane Float, Channel DIM1 Float) -> + (Plane Float, Channel DIM1 Float) updateWeightedCanvas = let distances = Run.with CUDA.run1 distanceMapGamma update = @@ -1354,9 +1357,9 @@ Geometry Float -> [Geometry Float] -> [Point2 Float] -> - Array DIM3 Word8 -> - (Channel Z Float, Channel DIM1 Float) -> - (Channel Z Float, Channel DIM1 Float) + ColorImage8 -> + (Plane Float, Channel DIM1 Float) -> + (Plane Float, Channel DIM1 Float) updateWeightedCanvasSplit = let update = Run.with CUDA.run1 addToWeightedCanvas distances = Run.with CUDA.run1 distanceMap @@ -1375,7 +1378,7 @@ finalizeWeightedCanvas :: - (Channel Z Float, Channel DIM1 Float) -> Array DIM3 Word8 + (Plane Float, Channel DIM1 Float) -> ColorImage8 finalizeWeightedCanvas = Run.with CUDA.run1 $ \(weightSum, canvas) -> @@ -1384,16 +1387,30 @@ replicateChannel (A.shape canvas) weightSum +data + Picture param = + Picture { + picPath :: FilePath, + picParam :: param, + picColored :: (Degree Float, ColorImage8), + picPlane :: ((Float, Float), Plane Float) + } + +mapPicParam :: (a -> b) -> Picture a -> Picture b +mapPicParam f pic = pic{picParam = f $ picParam pic} + + 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 + IO ([FilePath], [(Degree Float, ColorImage8)], + [((Float, Float), Complex Float)]) +processOverlap args = do let opt = Option.option args let info = CmdLine.info (Option.verbosity opt) + pics <- + map (mapPicParam (\(State.NoAngleCorrection, pos) -> pos)) <$> + processRotation args let padSize = Option.padSize opt let (maybeAllOverlapsShared, optimalOverlapShared) = case Just $ Z :. padSize :. padSize of @@ -1403,38 +1420,61 @@ Nothing -> let (padWidth, padHeight) = Arith.correlationSize (Option.minimumOverlap opt) $ - map (colorImageExtent . snd) picAngles + map (colorImageExtent . snd . picColored) pics padExtent = Z :. padHeight :. padWidth in (Just $ allOverlapsRun padExtent (Option.minimumOverlap opt), optimalOverlap padExtent (Option.minimumOverlap opt)) + relations <- + maybe (return Map.empty) + (State.readDisplacement (map picPath pics)) + (Option.relations opt) + + let open = map ((\(mx,my) -> isNothing mx || isNothing my) . picParam) pics displacements <- - fmap catMaybes $ - forM pairs $ \((ia,(pathA,(leftTopA,picA))), (ib,(pathB,(leftTopB,picB)))) -> do + forM (guardedPairs open $ zip [0..] pics) $ + \((ia, Picture pathA _ origA (leftTopA,picA)), + (ib, Picture pathB _ origB (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) = snd $ 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 relation = Map.lookup (pathA,pathB) relations + md <- + case (join $ fmap fst relation, join $ fmap snd relation) of + (Just State.NonOverlapping, _) -> return Nothing + (Just State.Overlapping, Just d) -> return $ Just d + (related, _) -> do + let doffset@(dox,doy) = snd $ optimalOverlapShared picA picB + let diff = overlapDifferenceRun doffset picA picB + let overlapping = + related == Just State.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 (origA, origB) + return $ toMaybe overlapping d + return ((ia,ib), (pathA,pathB), md) + forM_ (Option.outputState opt) $ \format -> + State.writeDisplacement (printf format "relation") displacements + + let overlaps = mapMaybe (\(i,_paths,md) -> (,) i <$> md) displacements let (poss, dps) = absolutePositionsFromPairDisplacements - (length picAngles) displacements + (fixAtLeastOnePosition (0,0) $ map picParam pics) overlaps info "\nabsolute positions" info $ unlines $ map show poss @@ -1443,13 +1483,13 @@ zipWith (\(dpx,dpy) (dx,dy) -> printf "(%f,%f) (%f,%f)" dpx dpy dx dy) - dps (map snd displacements) + dps (map snd overlaps) let (errdx,errdy) = - mapPair (maximum,maximum) $ unzip $ + mapPair (maximum0, maximum0) $ unzip $ zipWith (\(dpx,dpy) (dx,dy) -> (abs $ dpx - realToFrac dx, abs $ dpy - realToFrac dy)) - dps (map snd displacements) + dps (map snd overlaps) info $ "\n" @@ -1458,23 +1498,24 @@ ++ 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) + return + (map picPath pics, + map picColored pics, + map (flip (,) 1) $ map (mapPair (realToFrac, realToFrac)) poss) 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 + IO ([FilePath], [(Degree Float, ColorImage8)], + [((Float, Float), Complex Float)]) +processOverlapRotate args = do let opt = Option.option args let info = CmdLine.info (Option.verbosity opt) + let infoPlain = when (Option.verbosity opt >= Verbosity.verbose) . putStr + pics <- + map (mapPicParam (mapFst State.getAngleCorrection)) <$> + processRotation args let padSize = Option.padSize opt let stampSize = Option.stampSize opt let optimalOverlapShared = @@ -1482,105 +1523,127 @@ (Z :. padSize :. padSize) (Z :. stampSize :. stampSize) (Option.numberStamps opt) - (Option.maximumDifference opt) (Option.minimumOverlap opt) + relations <- + maybe (return Map.empty) + (State.readRotated (map picPath pics)) + (Option.relations opt) + + let open = + map + ((\(ma, (mx,my)) -> isNothing ma || isNothing mx || isNothing my) + . picParam) + pics 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 - (\(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 + forM (guardedPairs open $ zip [0..] pics) $ + \((ia, Picture pathA _ _ (leftTopA,picA)), + (ib, Picture pathB _ _ (leftTopB,picB))) -> do + let relation = Map.lookup (pathA,pathB) relations + correspondences <- + case (join $ fmap fst relation, Fold.fold $ fmap snd relation) of + (Just State.NonOverlapping, _) -> return [] + (Just State.Overlapping, corrs@(_:_)) -> return corrs + (related, _) -> do + let add (x0,y0) (x1,y1) = + (fromIntegral x0 + x1, fromIntegral y0 + y1) + let mMaxDiff = + toMaybe (related /= Just State.Overlapping) $ + Option.maximumDifference opt + let corrs = + map + (\(score,pa,pb) -> + (score, (add pa leftTopA, add pb leftTopB))) $ + optimalOverlapShared mMaxDiff picA picB + info $ printf "left-top: %s, %s" (show leftTopA) (show leftTopB) + info $ printf "%s - %s" pathA pathB + forM_ corrs $ \(score, (pa@(xa,ya),pb@(xb,yb))) -> + info $ + printf "%s ~ %s, (%f,%f), %f" + (show pa) (show pb) (xb-xa) (yb-ya) score + return $ map snd corrs + return ((ia,ib), (pathA,pathB), correspondences) + forM_ (Option.outputState opt) $ \format -> + State.writeRotated (printf format "relation") displacements + + let overlaps = concatMap (\(i,_paths,ps) -> (,) i <$> ps) displacements let (posRots, dps) = - layoutFromPairDisplacements (length picAngles) displacements + layoutFromPairDisplacements + (map (mapFst (fmap Degree.cis)) $ + fixAtLeastOneAnglePosition (Degree 0, (0,0)) $ + map picParam pics) + overlaps info "\nabsolute positions and rotations: place, rotation (magnitude, phase)" - info $ unlines $ + infoPlain $ unlines $ map - (\(d,r) -> - printf "%s, %s (%7.5f, %6.2f)" (show d) (show r) - (HComplex.magnitude r) (HComplex.phase r * 180/pi)) + (\((dx,dy),r) -> + printf "(%8.2f,%8.2f), %8.6f :+ %9.6f (%8.6f, %7.3f)" dx dy + (Complex.realPart r) (Complex.imagPart r) + (Complex.magnitude r) + (getDegree $ Degree.fromRadian $ Complex.phase r)) posRots info "\ncompare position differences with pair displacements" - info $ unlines $ + infoPlain $ 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 + (\(dpx,dpy) (_i, ((xa,ya),(xb,yb))) -> + printf "(%8.5f,%8.5f) (%7.2f,%7.2f) ~ (%7.2f,%7.2f)" + dpx dpy xa ya xb yb) + dps overlaps - return (floatPoss, picRots) + return + (map picPath pics, + map picColored pics, + map + (mapPair + (mapPair (realToFrac, realToFrac), Arith.mapComplex realToFrac)) + posRots) -process :: Option.Args -> IO () -process args = do - IO.hSetBuffering IO.stdout IO.LineBuffering - IO.hSetBuffering IO.stderr IO.LineBuffering - - let paths = Option.inputs args +processRotation :: + (State.AngleCorrected angleCorr) => + Option.Args -> IO [Picture (angleCorr, (Maybe Float, Maybe Float))] +processRotation args = do let opt = Option.option args let notice = CmdLine.notice (Option.verbosity opt) let info = CmdLine.info (Option.verbosity opt) + inputs <- Option.images args + notice "\nfind rotation angles" picAngles <- - forM paths $ \(imageOption, path) -> do + forM inputs $ \(State.Proposed path (maybeAngle, _) _) -> do pic <- readImage (Option.verbosity opt) path let maxAngle = Option.maximumAbsoluteAngle opt - let angles = - linearScale (Option.numberAngleSteps opt) - (-maxAngle, maxAngle) + let angles = Degree.linearScale (Option.numberAngleSteps opt) maxAngle when False $ analyseRotations angles pic when False $ fourierTransformation opt path pic angle <- - case Option.angle imageOption of + case maybeAngle of Just angle -> return angle Nothing -> if Option.radonTransform opt - then radonAngle (-maxAngle, maxAngle) pic + then radonAngle (fmap negate maxAngle, maxAngle) pic else return $ findOptimalRotation angles pic - info $ printf "%s %f\176\n" path angle - return (path, (angle*pi/180, pic)) + info $ printf "%s %f\176\n" path (getDegree angle) + return (angle, pic) + forM_ (Option.outputState opt) $ \format -> + State.write (printf format "angle") $ + zipWith State.Angle (map State.propPath inputs) (map fst picAngles) + 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) + let rotated = map (prepareOverlapMatching (Option.smooth opt)) picAngles when False $ do notice "write fft" - let pic0 : pic1 : _ = prepared + let pic0 : pic1 : _ = map snd rotated size = (Z:.512:.1024 :: DIM2) writeGrey (Option.quality opt) "/tmp/padded.jpeg" $ - CUDA.run1 - (imageByteFromFloat . - pad 0 (A.lift size)) $ - pic0 + CUDA.run1 (imageByteFromFloat . pad 0 (A.lift size)) pic0 writeGrey (Option.quality opt) "/tmp/spectrum.jpeg" $ - CUDA.run $ imageByteFromFloat $ A.map Complex.magnitude $ + CUDA.run $ imageByteFromFloat $ A.map AComplex.magnitude $ fft2DPlain CUFFT.forwardReal $ CUDA.run1 (pad 0 (A.lift size)) $ pic0 @@ -1588,15 +1651,38 @@ CUDA.run $ imageByteFromFloat $ A.map (0.000001*) $ correlatePadded size (A.use pic0) (A.use pic1) - (floatPoss, picRots) <- - (if Option.finetuneRotate opt - then processOverlapRotate - else processOverlap) - args (map snd picAngles) pairs + return $ + zipWith3 + (\(State.Proposed path (_,angleCorr) maybePos) colored plane -> + Picture path (angleCorr, maybePos) colored plane) + inputs picAngles rotated +process :: Option.Args -> IO () +process args = do + IO.hSetBuffering IO.stdout IO.LineBuffering + IO.hSetBuffering IO.stderr IO.LineBuffering + + let opt = Option.option args + let notice = CmdLine.notice (Option.verbosity opt) + let info = CmdLine.info (Option.verbosity opt) + + (paths, picAngles, posRots) <- + if Option.finetuneRotate opt + then processOverlapRotate args + else processOverlap args + + forM_ (Option.outputState opt) $ \format -> + State.write (printf format "position") $ + zipWith3 + (\path (angle, _) (pos, rot) -> + State.Position path + (angle <> Degree.fromRadian (Complex.phase rot)) pos) + paths picAngles posRots + notice "\ncompose all parts" let ((canvasWidth, canvasHeight), rotMovPics, canvasMsgs) = - Arith.canvasShape colorImageExtent floatPoss picRots + Arith.canvasShape colorImageExtent + (map (mapFst Degree.toRadian) picAngles) posRots mapM_ info canvasMsgs forM_ (Option.outputHard opt) $ \path -> @@ -1612,8 +1698,8 @@ Arith.geometryRelations $ map (Arith.geometryFeatures . mapThd3 colorImageExtent) rotMovPics - forM_ (zip geometryRelations picAngles) $ - \((thisGeom, otherGeoms, allPoints), (path, _)) -> do + forM_ (zip geometryRelations paths) $ + \((thisGeom, otherGeoms, allPoints), path) -> do let stem = FilePath.takeBaseName path let canvasShape = Z :. canvasHeight :. canvasWidth @@ -1643,7 +1729,7 @@ updateWeightedCanvas (Option.distanceGamma opt) thisGeom otherGeoms allPoints pic canvas) (emptyWeightedCanvas (Z :. 3 :. canvasHeight :. canvasWidth)) - (zip geometryRelations picRots) + (zip geometryRelations picAngles) main :: IO () main = process =<< Option.get Option.Accelerate diff --git a/src/Arithmetic.hs b/src/Arithmetic.hs --- a/src/Arithmetic.hs +++ b/src/Arithmetic.hs @@ -3,11 +3,14 @@ import qualified Data.Complex as Complex import Data.Complex (Complex, ) +import qualified Control.Functor.HT as FuncHT import Control.Monad (liftM2, guard) +import qualified Data.NonEmpty as NonEmpty import qualified Data.List.HT as ListHT import qualified Data.List as List import qualified Data.Bits as Bit +import Data.NonEmpty ((!:)) import Data.Maybe (catMaybes) import Data.Tuple.HT (mapPair, fst3, thd3) @@ -50,40 +53,44 @@ (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) : + FuncHT.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)) + in ((NonEmpty.foldBalanced mini xs, NonEmpty.foldBalanced maxi xs), + (NonEmpty.foldBalanced mini ys, NonEmpty.foldBalanced maxi ys)) canvasShape :: - (RealFrac a, Integral i, + (RealFloat a, Integral i, PrintfArg a, PrintfArg i) => - (image -> (i, i)) -> [Point2 a] -> [((a, a), image)] -> + (image -> (i, i)) -> [(a,image)] -> [(Point2 a, Complex a)] -> ((i, i), [((a, a), (a, a), image)], [String]) -canvasShape extent floatPoss picRots = - let bbox (rot, pic) = +canvasShape extent picAngles floatPosRots = + let posRotPics = + zipWith + (\(angle,pic) (pos,rot) -> + (pos, (pairFromComplex (Complex.cis angle * rot), pic))) + picAngles floatPosRots + bbox (rot, pic) = case extent pic of (width, height) -> boundingBoxOfRotated rot (fromIntegral width, fromIntegral height) + appPairs fs = mapPair fs . unzip + minimax = appPairs (minimum, maximum) + shift d = mapPair ((d+), (d+)) ((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 + appPairs (minimax, minimax) $ + map (uncurry $ \(mx,my) -> mapPair (shift mx, shift my) . bbox) $ + posRotPics canvasWidth = ceiling (canvasRight-canvasLeft) canvasHeight = ceiling (canvasBottom-canvasTop) rotMovPics = - zipWith - (\(mx,my) (rot, pic) -> (rot, (mx-canvasLeft, my-canvasTop), pic)) - floatPoss picRots + map + (\((mx,my), (rot,pic)) -> (rot, (mx-canvasLeft, my-canvasTop), pic)) + posRotPics in ((canvasWidth, canvasHeight), rotMovPics, [printf "canvas %f - %f, %f - %f\n" canvasLeft canvasRight canvasTop canvasBottom, @@ -310,3 +317,22 @@ mulConj :: (RealFloat a) => Complex a -> Complex a -> Complex a mulConj x y = x * Complex.conjugate y + + +-- ToDo: move to a new utility module +{- +This helps to compare only those pairs of images, +where at least one parameter is open to be determined. +Omitting pairs with completely fixed parameters +only accelerates the computation, +but the optimum will be the same. +-} +guardedPairs :: [Bool] -> [a] -> [(a,a)] +guardedPairs cs xs = do + ((takeA,a):as) <- ListHT.tails $ zip cs xs + (takeB,b) <- as + guard $ takeA || takeB + return (a,b) + +maximum0 :: (Num a, Ord a) => [a] -> a +maximum0 = NonEmpty.maximum . NonEmpty.cons 0 diff --git a/src/Degree.hs b/src/Degree.hs new file mode 100644 --- /dev/null +++ b/src/Degree.hs @@ -0,0 +1,35 @@ +module Degree where + +import qualified Arithmetic as Arith + +import qualified Data.Csv as Csv +import Data.Monoid (Monoid, mempty, mappend) + + +newtype Degree a = Degree {getDegree :: a} + deriving (Eq, Show) + +instance (Num a) => Monoid (Degree a) where + mempty = Degree 0 + mappend (Degree x) (Degree y) = Degree $ x+y + +instance Functor Degree where + fmap f (Degree x) = Degree $ f x + +instance Csv.FromField a => Csv.FromField (Degree a) where + parseField s = fmap Degree $ Csv.parseField s + +toRadian :: (Floating a) => Degree a -> a +toRadian (Degree angle) = angle*pi/180 + +fromRadian :: (Floating a) => a -> Degree a +fromRadian angle = Degree $ angle*180/pi + +linearScale :: (Fractional a) => Int -> Degree a -> [Degree a] +linearScale num (Degree maxAngle) = + map Degree $ Arith.linearScale num (-maxAngle, maxAngle) + +cis :: (Floating a) => Degree a -> (a,a) +cis deg = + case toRadian deg of + rad -> (cos rad, sin rad) diff --git a/src/Knead.hs b/src/Knead.hs --- a/src/Knead.hs +++ b/src/Knead.hs @@ -2,16 +2,21 @@ module Main where import qualified Option +import qualified State import qualified MatchImageBorders import qualified Arithmetic as Arith +import qualified Degree import MatchImageBorders (arrayCFromKnead, arrayKneadFromC) +import Arithmetic (guardedPairs, maximum0) import LinearAlgebra ( - absolutePositionsFromPairDisplacements, layoutFromPairDisplacements, + absolutePositionsFromPairDisplacements, fixAtLeastOnePosition, + layoutFromPairDisplacements, fixAtLeastOneAnglePosition, ) import KneadShape (Size, Vec2(Vec2), Dim1, Dim2, Shape2, Index2, Ix2, verticalVal, horizontalVal) +import Degree (Degree(Degree), getDegree) import qualified Math.FFT as FFT import Math.FFT.Base (FFTWReal) @@ -55,18 +60,23 @@ import Distribution.Verbosity (Verbosity) import Text.Printf (printf) -import qualified Control.Functor.HT as FuncHT -import Control.Monad (liftM2, when, foldM, (<=<)) +import qualified Control.Monad.HT as MonadHT +import Control.Monad (liftM2, when, join, foldM, (<=<)) import Control.Applicative (pure, (<$>), (<*>)) +import qualified Data.Foldable as Fold import qualified Data.List as List +import qualified Data.Map as Map +import Data.Function.HT (Id) +import Data.Monoid ((<>)) import Data.Maybe.HT (toMaybe) -import Data.Maybe (catMaybes, isJust) -import Data.List.HT (tails) +import Data.Maybe (mapMaybe, isJust, isNothing) 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.Tuple.HT + (mapPair, mapFst, mapSnd, mapTriple, swap, + mapThd3, fst3, thd3, uncurry3) import Data.Word (Word8, Word32) @@ -456,12 +466,12 @@ img -runRotate :: IO (Float -> ColorImage8 -> IO ColorImage8) +runRotate :: IO (Degree Float -> ColorImage8 -> IO ColorImage8) runRotate = do rot <- RenderP.run $ \rot -> colorImageByteFromFloat . rotate vecYUV rot . colorImageFloatFromByte - return $ \ angle img -> rot (cos angle, sin angle) img + return $ \ angle img -> rot (Degree.cis angle) img brightnessValue :: Exp (YUV a) -> Exp a @@ -502,21 +512,20 @@ scoreHistogram = fold1All Expr.add . Symb.map Expr.sqr . differentiate -runScoreRotation :: IO (Float -> ColorImage8 -> IO Float) +runScoreRotation :: IO (Degree 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 + return $ \ angle img -> score =<< rot (Degree.cis angle) img -findOptimalRotation :: IO ([Float] -> ColorImage8 -> IO Float) +findOptimalRotation :: IO ([Degree Float] -> ColorImage8 -> IO (Degree Float)) findOptimalRotation = do scoreRotation <- runScoreRotation return $ \angles pic -> fmap (fst . List.maximumBy (comparing snd)) $ - forM angles $ \angle -> - (,) angle <$> scoreRotation (angle * (pi/180)) pic + forM angles $ \angle -> (,) angle <$> scoreRotation angle pic @@ -539,13 +548,10 @@ 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 + return $ \n -> MonadHT.nest n lp highpassMulti :: IO (Int -> Plane Float -> IO (Plane Float)) @@ -614,6 +620,20 @@ untangleSpectra2d spec = CArray.liftArray2 untangleCoefficient spec (cyclicReverse2d spec) +{- | +Equivalent to @amap (uncurry Arith.mulConj) . untangleSpectra2d@ +but much faster, since it avoids the slow @instance Storable (a,b)@ +based on @storable-tuple:storePair@. +-} +mulConjUntangledSpectra2d :: + (RealFloat a, SV.Storable a) => + CArray (Int,Int) (Complex a) -> CArray (Int,Int) (Complex a) +mulConjUntangledSpectra2d spec = + CArray.liftArray2 + ((uncurry Arith.mulConj .) . untangleCoefficient) + spec (cyclicReverse2d spec) + + {- This is more efficient than 'correlatePaddedSimpleCArray' since it needs only one complex forward Fourier transform, @@ -632,8 +652,7 @@ CArray (Int,Int) a correlatePaddedComplexCArray sh a b = amap realPart $ FFT.idftN [0,1] $ - amap (uncurry Arith.mulConj) $ - untangleSpectra2d $ FFT.dftN [0,1] $ + mulConjUntangledSpectra2d $ FFT.dftN [0,1] $ CArray.liftArray2 (:+) (padCArray 0 sh a) (padCArray 0 sh b) {- | @@ -650,8 +669,7 @@ (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] $ + mulConjUntangledSpectra2d $ FFT.dftN [0,1] $ CArray.liftArray2 (:+) (padCArray 0 sh a) (padCArray 0 sh b) @@ -666,20 +684,18 @@ (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)) + IO (Int -> (Degree 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) + rot = Degree.cis angle ((left, _right), (top, _bottom)) = Arith.boundingBoxOfRotated rot (fromIntegral width, fromIntegral height) @@ -897,7 +913,7 @@ -} optimalOverlapBigMulti :: Dim2 -> Dim2 -> Int -> - IO (Float -> Float -> Plane Float -> Plane Float -> + IO (Float -> Maybe 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 @@ -918,7 +934,7 @@ clipB <- clp anchorB extent b addCoarsePos <$> optOverlapFine minimumOverlap clipA clipB - return $ \maximumDiff minimumOverlap a b -> do + return $ \minimumOverlap mMaximumDiff a b -> do let factors@(Vec2 yk xk) = shrinkFactors padExtent (Phys.shape a) (Phys.shape b) aSmall <- shrnk factors a @@ -930,7 +946,11 @@ let coarsedy = shrunkdy * yk let coarsed = (coarsedx,coarsedy) - diff <- overDiff shrunkd aSmall bSmall + doesOverlap <- + case mMaximumDiff of + Just maximumDiff -> + (maximumDiff>) <$> overDiff shrunkd aSmall bSmall + Nothing -> return True let ((leftOverlap, topOverlap), (rightOverlap, bottomOverlap), @@ -940,7 +960,7 @@ let widthStampClip = min widthOverlap widthStamp heightStampClip = min heightOverlap heightStamp - (if diff < maximumDiff then id else const $ return []) $ + (if doesOverlap then id else const $ return []) $ mapM (\(x,y) -> overlapFine minimumOverlap a b @@ -1013,7 +1033,7 @@ composeOverlap :: IO ((Size, Size) -> - ((Float, ColorImage8), (Float, ColorImage8)) -> + ((Degree Float, ColorImage8), (Degree Float, ColorImage8)) -> IO ColorImage8) composeOverlap = do over <- @@ -1022,9 +1042,9 @@ 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) + over displacement + (Degree.cis angleA, picA) (Degree.cis angleB, picB) @@ -1223,7 +1243,19 @@ (uncurry (Expr.ifThenElse (expEven sel)) vert))) array (Symb.lift0 $ Symb.id 4) +distanceMapBoxRun :: + IO (Dim2 -> Geometry Float -> IO (Plane Word8)) +distanceMapBoxRun = + RenderP.run $ \sh geom -> + scaleDistanceMapGeom geom $ + Symb.map + (Expr.modify (atom,atom) $ \(valid, dist) -> Expr.select valid dist 0) $ + maskedMinimum $ + Symb.map (Expr.mapSnd Expr.fst) $ + separateDistanceMap $ + distanceMapBox sh geom + containedAnywhere :: (Symb.C array, Shape.C sh, MultiValue.Field a, MultiValue.NativeFloating a ar, @@ -1261,10 +1293,27 @@ (Expr.liftM2 MultiValue.and c b, dist)) contained distMap +distanceMapContainedRun :: + IO (Dim2 -> Geometry Float -> [Geometry Float] -> IO (Plane Word8)) +distanceMapContainedRun = do + distances <- + RenderP.run $ + \sh this -> scaleDistanceMapGeom this . distanceMapContained sh this + return $ \sh this others -> distances sh this =<< Phys.vectorFromList others +scaleDistanceMapGeom :: + (MultiValue.Field a, MultiValue.Real a, MultiValue.RationalConstant a, + MultiValue.NativeFloating a ar) => + Exp (Geometry b) -> SymbPlane a -> SymbPlane Word8 +scaleDistanceMapGeom geom img = + let scale = + (4/) $ fromInt $ + Expr.modify (atom,atom,(atom,atom)) (uncurry Expr.min . thd3) geom + in imageByteFromFloat $ Symb.map (scale*) img + + pixelCoordinates :: - (MultiValue.NativeFloating a ar) => - Exp Dim2 -> SymbPlane (a,a) + (MultiValue.NativeFloating a ar) => Exp Dim2 -> SymbPlane (a,a) pixelCoordinates sh = generate sh $ Expr.modify atomIx2 $ \(Vec2 y x) -> fromSize2 (x,y) @@ -1278,7 +1327,26 @@ Symb.fold1 Expr.min $ outerProduct (Expr.modify2 (atom,atom) (atom,atom) distance) a b +distanceMapPointsRun :: + IO (Dim2 -> [Arith.Point2 Float] -> IO (Plane Word8)) +distanceMapPointsRun = do + distances <- + RenderP.run $ + \sh -> scaleDistanceMap . distanceMapPoints (pixelCoordinates sh) + return $ \sh points -> distances sh =<< Phys.vectorFromList points + +scaleDistanceMap :: + (MultiValue.Field a, MultiValue.Real a, MultiValue.RationalConstant a, + MultiValue.NativeFloating a ar) => + SymbPlane a -> SymbPlane Word8 +scaleDistanceMap img = + let scale = + case Expr.decompose atomDim2 $ Symb.shape img of + Vec2 y x -> 4 / fromInt (Expr.min x y) + in imageByteFromFloat $ Symb.map (scale*) img + + {- | For every pixel it computes the distance to the closest point on the image part boundary @@ -1308,7 +1376,22 @@ (distanceMapContained sh this others) (distanceMapPoints (pixelCoordinates sh) points) +distanceMapRun :: + IO (Dim2 -> + Geometry Float -> + [Geometry Float] -> + [Arith.Point2 Float] -> + IO (Plane Word8)) +distanceMapRun = do + distances <- + RenderP.run $ + \sh this others -> scaleDistanceMap . distanceMap sh this others + return $ \sh this others points -> do + othersVec <- Phys.vectorFromList others + pointsVec <- Phys.vectorFromList points + distances sh this othersVec pointsVec + pow :: (MultiValue.Repr LLVM.Value a ~ LLVM.Value ar, LLVM.IsFloating ar, SoV.TranscendentalConstant ar) => @@ -1381,7 +1464,7 @@ update this pic dists weightSumCanvas -finalizeWeightedCanvas :: IO ((Plane (Float, YUV Float)) -> IO ColorImage8) +finalizeWeightedCanvas :: IO (Plane (Float, YUV Float) -> IO ColorImage8) finalizeWeightedCanvas = RenderP.run $ colorImageByteFromFloat . @@ -1391,17 +1474,30 @@ fixArray +data + Picture param = + Picture { + picPath :: FilePath, + picParam :: param, + picColored :: (Degree Float, ColorImage8), + picPlane :: ((Float, Float), Plane Float) + } +mapPicParam :: (a -> b) -> Picture a -> Picture b +mapPicParam f pic = pic{picParam = f $ picParam pic} + + 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 + IO ([FilePath], [(Degree Float, ColorImage8)], + [((Float, Float), Complex Float)]) +processOverlap args = do let opt = Option.option args let info = CmdLine.info (Option.verbosity opt) + pics <- + map (mapPicParam (\(State.NoAngleCorrection, pos) -> pos)) <$> + processRotation args let padSize = fromIntegral $ Option.padSize opt (maybeAllOverlapsShared, optimalOverlapShared) <- case Just $ Vec2 padSize padSize of @@ -1412,43 +1508,65 @@ let padExtent = uncurry Vec2 $ swap $ Arith.correlationSize (Option.minimumOverlap opt) $ - map (colorImageExtent . snd) picAngles + map (colorImageExtent . snd . picColored) pics overlap <- optimalOverlap padExtent allOverlapsIO <- allOverlapsRun padExtent return (Just $ allOverlapsIO (Option.minimumOverlap opt), overlap (Option.minimumOverlap opt)) + relations <- + maybe (return Map.empty) + (State.readDisplacement (map picPath pics)) + (Option.relations opt) + composeOver <- composeOverlap overlapDiff <- overlapDifferenceRun + let open = map ((\(mx,my) -> isNothing mx || isNothing my) . picParam) pics displacements <- - fmap catMaybes $ - forM pairs $ \((ia,(pathA,(leftTopA,picA))), (ib,(pathB,(leftTopB,picB)))) -> do - forM_ maybeAllOverlapsShared $ \allOverlapsShared -> when True $ + forM (guardedPairs open $ zip [0..] pics) $ + \((ia, Picture pathA _ origA (leftTopA,picA)), + (ib, Picture pathB _ origB (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 - 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 relation = Map.lookup (pathA,pathB) relations + md <- + case (join $ fmap fst relation, join $ fmap snd relation) of + (Just State.NonOverlapping, _) -> return Nothing + (Just State.Overlapping, Just d) -> return $ Just d + (related, _) -> do + doffset@(dox,doy) <- snd <$> optimalOverlapShared picA picB + diff <- overlapDiff doffset picA picB + let overlapping = + related == Just State.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)) + =<< composeOver doffset (origA, origB) + return $ toMaybe overlapping d + return ((ia,ib), (pathA,pathB), md) + forM_ (Option.outputState opt) $ \format -> + State.writeDisplacement (printf format "relation") displacements + + let overlaps = mapMaybe (\(i,_paths,md) -> (,) i <$> md) displacements let (poss, dps) = absolutePositionsFromPairDisplacements - (length picAngles) displacements + (fixAtLeastOnePosition (0,0) $ map picParam pics) overlaps info "\nabsolute positions" info $ unlines $ map show poss @@ -1457,13 +1575,13 @@ zipWith (\(dpx,dpy) (dx,dy) -> printf "(%f,%f) (%f,%f)" dpx dpy dx dy) - dps (map snd displacements) + dps (map snd overlaps) let (errdx,errdy) = - mapPair (maximum,maximum) $ unzip $ + mapPair (maximum0, maximum0) $ unzip $ zipWith (\(dpx,dpy) (dx,dy) -> (abs $ dpx - realToFrac dx, abs $ dpy - realToFrac dy)) - dps (map snd displacements) + dps (map snd overlaps) info $ "\n" @@ -1472,23 +1590,24 @@ ++ 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) + return + (map picPath pics, + map picColored pics, + map (flip (,) 1) $ map (mapPair (realToFrac, realToFrac)) poss) 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 + IO ([FilePath], [(Degree Float, ColorImage8)], + [((Float, Float), Complex Float)]) +processOverlapRotate args = do let opt = Option.option args let info = CmdLine.info (Option.verbosity opt) + let infoPlain = when (Option.verbosity opt >= Verbosity.verbose) . putStr + pics <- + map (mapPicParam (mapFst State.getAngleCorrection)) <$> + processRotation args let padSize = Option.padSize opt let stampSize = Option.stampSize opt optimalOverlapShared <- @@ -1496,99 +1615,166 @@ (shape2 padSize padSize) (shape2 stampSize stampSize) (Option.numberStamps opt) - <*> pure (Option.maximumDifference opt) <*> pure (Option.minimumOverlap opt) + relations <- + maybe (return Map.empty) + (State.readRotated (map picPath pics)) + (Option.relations opt) + + let open = + map + ((\(ma, (mx,my)) -> isNothing ma || isNothing mx || isNothing my) + . picParam) + pics 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) + forM (guardedPairs open $ zip [0..] pics) $ + \((ia, Picture pathA _ _ (leftTopA,picA)), + (ib, Picture pathB _ _ (leftTopB,picB))) -> do + let relation = Map.lookup (pathA,pathB) relations 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 + case (join $ fmap fst relation, Fold.fold $ fmap snd relation) of + (Just State.NonOverlapping, _) -> return [] + (Just State.Overlapping, corrs@(_:_)) -> return corrs + (related, _) -> do + let add (x0,y0) (x1,y1) = + (fromIntegral x0 + x1, fromIntegral y0 + y1) + let mMaxDiff = + toMaybe (related /= Just State.Overlapping) $ + Option.maximumDifference opt + corrs <- + map + (\(score,pa,pb) -> + (score, (add pa leftTopA, add pb leftTopB))) <$> + optimalOverlapShared mMaxDiff picA picB + info $ printf "left-top: %s, %s" (show leftTopA) (show leftTopB) + info $ printf "%s - %s" pathA pathB + forM_ corrs $ \(score, (pa@(xa,ya),pb@(xb,yb))) -> + info $ + printf "%s ~ %s, (%f,%f), %f" + (show pa) (show pb) (xb-xa) (yb-ya) score + return $ map snd corrs + return ((ia,ib), (pathA,pathB), correspondences) + forM_ (Option.outputState opt) $ \format -> + State.writeRotated (printf format "relation") displacements + + let overlaps = concatMap (\(i,_paths,ps) -> (,) i <$> ps) displacements let (posRots, dps) = - layoutFromPairDisplacements (length picAngles) displacements + layoutFromPairDisplacements + (map (mapFst (fmap Degree.cis)) $ + fixAtLeastOneAnglePosition (Degree 0, (0,0)) $ + map picParam pics) + overlaps info "\nabsolute positions and rotations: place, rotation (magnitude, phase)" - info $ unlines $ + infoPlain $ unlines $ map - (\(d,r) -> - printf "%s, %s (%7.5f, %6.2f)" (show d) (show r) - (Complex.magnitude r) (Complex.phase r * 180/pi)) + (\((dx,dy),r) -> + printf "(%8.2f,%8.2f), %8.6f :+ %9.6f (%8.6f, %7.3f)" dx dy + (Complex.realPart r) (Complex.imagPart r) + (Complex.magnitude r) + (getDegree $ Degree.fromRadian $ Complex.phase r)) posRots info "\ncompare position differences with pair displacements" - info $ unlines $ + infoPlain $ 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 + (\(dpx,dpy) (_i, ((xa,ya),(xb,yb))) -> + printf "(%8.5f,%8.5f) (%7.2f,%7.2f) ~ (%7.2f,%7.2f)" + dpx dpy xa ya xb yb) + dps overlaps - return (floatPoss, picRots) + return + (map picPath pics, + map picColored pics, + map + (mapPair + (mapPair (realToFrac, realToFrac), Arith.mapComplex realToFrac)) + posRots) -process :: Option.Args -> IO () -process args = do - IO.hSetBuffering IO.stdout IO.LineBuffering - IO.hSetBuffering IO.stderr IO.LineBuffering - - let paths = Option.inputs args +processRotation :: + (State.AngleCorrected angleCorr) => + Option.Args -> IO [Picture (angleCorr, (Maybe Float, Maybe Float))] +processRotation args = do let opt = Option.option args let notice = CmdLine.notice (Option.verbosity opt) let info = CmdLine.info (Option.verbosity opt) + inputs <- Option.images args + notice "\nfind rotation angles" findOptRot <- findOptimalRotation picAngles <- - forM paths $ \(imageOption, path) -> do + forM inputs $ \(State.Proposed path (maybeAngle, _) _) -> do pic <- readImage (Option.verbosity opt) path let maxAngle = Option.maximumAbsoluteAngle opt - let angles = - Arith.linearScale (Option.numberAngleSteps opt) - (-maxAngle, maxAngle) + let angles = Degree.linearScale (Option.numberAngleSteps opt) maxAngle angle <- - case Option.angle imageOption of + case maybeAngle of Just angle -> return angle Nothing -> findOptRot angles pic - info $ printf "%s %f\176\n" path angle - return (path, (angle*pi/180, pic)) + info $ printf "%s %f\176\n" path (getDegree angle) + return (angle, pic) + forM_ (Option.outputState opt) $ \format -> + State.write (printf format "angle") $ + zipWith State.Angle (map State.propPath inputs) (map fst picAngles) + 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) + rotated <- mapM (prepOverlapMatching (Option.smooth opt)) picAngles - (floatPoss, picRots) <- - (if Option.finetuneRotate opt - then processOverlapRotate - else processOverlap) - args (map snd picAngles) pairs + when False $ do + notice "write fft" + let pic0 : pic1 : _ = map snd rotated + size = (1024,768) + cpic0 <- arrayCFromKnead pic0 + cpic1 <- arrayCFromKnead pic1 + makeByteImage <- + RenderP.run $ \k -> imageByteFromFloat . Symb.map (k*) . fixArray + writeGrey (Option.quality opt) "/tmp/padded.jpeg" =<< + (makeByteImage 1 $ arrayKneadFromC $ padCArray 0 size cpic0) + writeGrey (Option.quality opt) "/tmp/spectrum.jpeg" =<< + (makeByteImage 0.1 $ arrayKneadFromC $ + CArray.liftArray Complex.magnitude $ + FFT.dftRCN [0,1] $ padCArray 0 size cpic0) + writeGrey (Option.quality opt) "/tmp/convolution.jpeg" =<< + (makeByteImage 0.1 $ arrayKneadFromC $ + correlatePaddedCArray size cpic0 cpic1) + return $ + zipWith3 + (\(State.Proposed path (_,angleCorr) maybePos) colored plane -> + Picture path (angleCorr, maybePos) colored plane) + inputs picAngles rotated + +process :: Option.Args -> IO () +process args = do + IO.hSetBuffering IO.stdout IO.LineBuffering + IO.hSetBuffering IO.stderr IO.LineBuffering + + let opt = Option.option args + let notice = CmdLine.notice (Option.verbosity opt) + let info = CmdLine.info (Option.verbosity opt) + + (paths, picAngles, posRots) <- + if Option.finetuneRotate opt + then processOverlapRotate args + else processOverlap args + + forM_ (Option.outputState opt) $ \format -> + State.write (printf format "position") $ + zipWith3 + (\path (angle, _) (pos, rot) -> + State.Position path + (angle <> Degree.fromRadian (Complex.phase rot)) pos) + paths picAngles posRots + notice "\ncompose all parts" let ((canvasWidth, canvasHeight), rotMovPics, canvasMsgs) = - Arith.canvasShape colorImageExtent floatPoss picRots + Arith.canvasShape colorImageExtent + (map (mapFst Degree.toRadian) picAngles) posRots let canvasShape = shape2 canvasHeight canvasWidth mapM_ info canvasMsgs @@ -1607,6 +1793,35 @@ Arith.geometryRelations $ map (Arith.geometryFeatures . mapThd3 colorImageExtent) rotMovPics + forM_ (Option.outputDistanceMap opt) $ \format -> do + debug <- + if True + then return $ \ _stem _geoms -> return () + else do + distMapBox <- distanceMapBoxRun + distMapContained <- distanceMapContainedRun + distMapPoints <- distanceMapPointsRun + + return $ \stem (thisGeom, otherGeoms, allPoints) -> do + writeGrey (Option.quality opt) + (printf "/tmp/%s-distance-box.jpeg" stem) + =<< distMapBox canvasShape thisGeom + + writeGrey (Option.quality opt) + (printf "/tmp/%s-distance-contained.jpeg" stem) + =<< distMapContained canvasShape thisGeom otherGeoms + + writeGrey (Option.quality opt) + (printf "/tmp/%s-distance-points.jpeg" stem) + =<< distMapPoints canvasShape allPoints + + distMap <- distanceMapRun + forM_ (zip geometryRelations paths) $ \(geoms, path) -> do + let stem = FilePath.takeBaseName path + debug stem geoms + writeGrey (Option.quality opt) (printf format stem) =<< + uncurry3 (distMap canvasShape) geoms + forM_ (Option.output opt) $ \path -> do notice "\nweighted composition" emptyCanv <- emptyWeightedCanvas @@ -1619,7 +1834,7 @@ (\canvas ((thisGeom, otherGeoms, allPoints), (_rot, pic)) -> updateCanv (Option.distanceGamma opt) thisGeom otherGeoms allPoints pic canvas) - empty (zip geometryRelations picRots) + empty (zip geometryRelations picAngles) when (isJust (Option.outputShaped opt) || isJust (Option.outputShapedHard opt)) $ do notice "\nmatch shapes" @@ -1644,7 +1859,7 @@ counts <- thaw . amap (fromIntegral . fst) =<< arrayCFromKnead sumImg shapes <- MatchImageBorders.shapeParts counts locs pqueue - let names = map (FilePath.takeBaseName . fst) picAngles + let names = map FilePath.takeBaseName paths forM_ (Option.outputShapedHard opt) $ \path -> do forM_ (Option.outputShapeHard opt) $ \format -> forM_ (zip names shapes) $ \(name,shape) -> @@ -1693,14 +1908,14 @@ forM_ [0..11] $ \k -> do let path = printf "/tmp/rotated/%04d.jpeg" k putStrLn path - writeImage 100 path =<< rot (fromInteger k * pi/6) img + writeImage 100 path =<< rot (Degree $ fromInteger k * 30) 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 + print =<< score (Degree $ fromInteger k / 10) img main :: IO () main = process =<< Option.get Option.Knead diff --git a/src/LinearAlgebra.hs b/src/LinearAlgebra.hs --- a/src/LinearAlgebra.hs +++ b/src/LinearAlgebra.hs @@ -6,43 +6,81 @@ import qualified Numeric.Container as Container import Numeric.Container ((<\>), (<>)) -import qualified Data.Complex as HComplex +import Data.Complex (Complex((:+))) import qualified Data.List.HT as ListHT import qualified Data.List as List -import Control.Monad (zipWithM_) +import Data.Tuple.HT (mapPair, mapSnd) +import Data.Maybe (isJust, fromMaybe) +import Control.Applicative ((<$>)) --- we cannot use leastSquaresSelected here, because the right-hand side is not zero + +fixAtLeastOne :: a -> [Maybe a] -> [Maybe a] +fixAtLeastOne zero ms = + case (any isJust ms, ms) of + (True, _) -> ms + (False, _:nothings) -> Just zero : nothings + (False, []) -> error "fixAtLeastOne: empty image list" + +{- | +If no coordinate is fixed, then the first one will be fixed to the given value. +This is not strictly necessary. +Without a fixed coordinate, +the solver will center all solutions around zero. +However, there will not necessarily be an image with a zero coordinate, +which is somehow ugly. +-} +fixAtLeastOnePosition :: + (a,b) -> [(Maybe a, Maybe b)] -> [(Maybe a, Maybe b)] +fixAtLeastOnePosition (a,b) = + parallel (fixAtLeastOne a, fixAtLeastOne b) + +fixAtLeastOneAnglePosition :: + (angle, (a,b)) -> + [(Maybe angle, (Maybe a, Maybe b))] -> + [(Maybe angle, (Maybe a, Maybe b))] +fixAtLeastOneAnglePosition (angle, ab) = + parallel (fixAtLeastOne angle, fixAtLeastOnePosition ab) + +parallel :: ([a0] -> [a1], [b0] -> [b1]) -> ([(a0,b0)] -> [(a1,b1)]) +parallel fs = uncurry zip . mapPair fs . unzip + + +sparseMatrix :: Int -> Int -> [((Int, Int), Double)] -> Matrix.Matrix Double +sparseMatrix numRows numCols xs = + PackST.runSTMatrix $ do + mat <- PackST.newMatrix 0 numRows numCols + mapM_ (\((r,c), x) -> PackST.writeMatrix mat r c x) xs + return mat + +elm :: Int -> Int -> a -> ((Int, Int), a) +elm row col x = ((row, col), x) + + absolutePositionsFromPairDisplacements :: - Int -> [((Int, Int), (Float, Float))] -> + [(Maybe Float, Maybe Float)] -> [((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)) +absolutePositionsFromPairDisplacements mxys displacements = + let numPics = length mxys + (mxs, mys) = unzip mxys + (is, (dxs, dys)) = mapSnd unzip $ unzip displacements + matrix = + sparseMatrix (length is) numPics $ concat $ + zipWith (\k (ia,ib) -> [elm k ia (-1), elm k ib 1]) [0..] is + solve ms ds = + leastSquaresSelected matrix + (map (fmap realToFrac) ms) + (Vector.fromList (map realToFrac ds)) + (pxs, achievedDxs) = solve mxs dxs + (pys, achievedDys) = solve mys dys + in (zip pxs pys, zip achievedDxs achievedDys) leastSquaresSelected :: - Matrix.Matrix Double -> [Maybe Double] -> + Matrix.Matrix Double -> [Maybe Double] -> Vector.Vector Double -> ([Double], [Double]) -leastSquaresSelected m mas = +leastSquaresSelected m mas rhs0 = let (lhsCols,rhsCols) = ListHT.unzipEithers $ zipWith @@ -52,9 +90,10 @@ 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 $ + rhs = foldl Container.add (Container.scale 0 rhs0) rhsCols + sol = lhs <\> Container.sub rhs0 rhs + in if Vector.dim rhs0 == 0 then (map (fromMaybe 0) mas, []) else + (snd $ List.mapAccumL (curry $ \x -> case x of @@ -65,6 +104,12 @@ Vector.toList $ Container.add (lhs <> sol) rhs) + +zeroVector, _zeroVector :: Int -> Vector.Vector Double +zeroVector n = Vector.fromList $ replicate n 0 +-- fails for vectors of size 0 +_zeroVector n = Container.constant 0 n + {- Approximate rotation from point correspondences. Here (dx, dy) is the displacement with respect to the origin (0,0), @@ -78,56 +123,58 @@ \0 1 y0 x0/ |c | \y1/ \s / -Maybe, dx and dy should be scaled down. +We try to scale dx and dy down using 'weight'. Otherwise they are weighted much more than the rotation. +However the weight will only influence the result +for under-constrained equation systems. +This is usually not the case. -} layoutFromPairDisplacements :: - Int -> [((Int, (Float, Float)), (Int, (Float, Float)))] -> - ([((Double,Double), HComplex.Complex Double)], + [(Maybe (Float, Float), (Maybe Float, Maybe Float))] -> + [((Int, Int), ((Float, Float), (Float, Float)))] -> + ([((Double,Double), 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. - -} +layoutFromPairDisplacements mrxys correspondences = + let numPics = length mrxys weight = let xs = concatMap - (\((_ia,(xai,yai)),(_ib,(xbi,ybi))) -> [xai, yai, xbi, ybi]) + (\(_i, ((xai,yai),(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 + in if null xs then 1 else realToFrac $ maximum xs - minimum xs + matrix = + sparseMatrix (2 * length correspondences) (4*numPics) $ concat $ + zipWith + (\k ((ia,ib), ((xai,yai),(xbi,ybi))) -> let xa = realToFrac xai - 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) + xb = realToFrac xbi + ya = realToFrac yai + yb = realToFrac ybi + in elm (k+0) (4*ia+0) (-weight) : + elm (k+1) (4*ia+1) (-weight) : + elm (k+0) (4*ia+2) (-xa) : + elm (k+0) (4*ia+3) ya : + elm (k+1) (4*ia+2) (-ya) : + elm (k+1) (4*ia+3) (-xa) : + elm (k+0) (4*ib+0) weight : + elm (k+1) (4*ib+1) weight : + elm (k+0) (4*ib+2) xb : + elm (k+0) (4*ib+3) (-yb) : + elm (k+1) (4*ib+2) yb : + elm (k+1) (4*ib+3) xb : + []) [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)) $ + (concatMap + (\(mr, (mx,my)) -> + [(/weight) . realToFrac <$> mx, + (/weight) . realToFrac <$> my, + realToFrac . fst <$> mr, + realToFrac . snd <$> mr]) $ + mrxys) + (zeroVector (2 * length correspondences)) + in (map (\[dx,dy,rx,ry] -> ((weight*dx,weight*dy), rx :+ ry)) $ ListHT.sliceVertical 4 solution, map (\[x,y] -> (x,y)) $ ListHT.sliceVertical 2 projection) diff --git a/src/Option.hs b/src/Option.hs --- a/src/Option.hs +++ b/src/Option.hs @@ -1,16 +1,23 @@ module Option where +import qualified State import Option.Utility (exitFailureMsg, parseNumber, fmapOptDescr) +import Degree (Degree(Degree, getDegree)) 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 qualified System.IO as IO import Control.Monad (when) import qualified Data.EnumSet as EnumSet +import qualified Data.Vector as Vector +import qualified Data.Map as Map +import qualified Data.List.HT as ListHT +import qualified Data.List as List import Data.Tuple.HT (mapSnd) import Data.Monoid ((<>)) import Data.Word (Word8) @@ -35,6 +42,8 @@ data Option = Option { verbosity :: Verbosity, + state :: Maybe FilePath, + relations :: Maybe FilePath, output :: Maybe FilePath, outputHard :: Maybe FilePath, outputShaped :: Maybe FilePath, @@ -43,8 +52,9 @@ outputDistanceMap :: Maybe String, outputShape :: Maybe String, outputShapeHard :: Maybe String, + outputState :: Maybe String, quality :: Int, - maximumAbsoluteAngle :: Float, + maximumAbsoluteAngle :: Degree Float, numberAngleSteps :: Int, radonTransform :: Bool, smooth :: Int, @@ -62,6 +72,8 @@ defltOption = Option { verbosity = Verbosity.verbose, + state = Nothing, + relations = Nothing, output = Nothing, outputHard = Nothing, outputShaped = Nothing, @@ -70,8 +82,9 @@ outputDistanceMap = Nothing, outputShape = Nothing, outputShapeHard = Nothing, + outputState = Nothing, quality = 99, - maximumAbsoluteAngle = 1, + maximumAbsoluteAngle = Degree 1, numberAngleSteps = 40, radonTransform = False, smooth = 20, @@ -88,14 +101,20 @@ data Image = Image { - angle :: Maybe Float + angle :: Maybe (Degree Float) } deriving (Eq) defltImage :: Image defltImage = Image {angle = Nothing} +proposedFromImage :: + (State.AngleCorrected angleCorr) => + (Image, FilePath) -> State.Proposed angleCorr +proposedFromImage (Image ang, path) = + State.Proposed path (ang, State.autoAngleCorrection) (Nothing, Nothing) + data Engine = Knead | Accelerate deriving (Eq, Ord, Enum) @@ -141,6 +160,16 @@ (printf "verbosity level: 0..3, default: %d" (fromEnum $ verbosity defltOption)) : + opt generic [] ["state"] + (flip ReqArg "PATH" $ \path flags -> + return $ flags{state = Just path}) + ("CSV file with predefined parameters") : + + opt generic [] ["relations"] + (flip ReqArg "PATH" $ \str flags -> + return $ flags{relations = Just str}) + ("CSV file with image pairs") : + opt generic [] ["output"] (flip ReqArg "PATH" $ \str flags -> return $ flags{output = Just str}) @@ -181,6 +210,11 @@ return $ flags{outputShapeHard = Just str}) ("path format for hard part shape like '%s-shape-hard.jpeg'") : + opt generic [] ["output-state"] + (flip ReqArg "FORMAT" $ \str flags -> + return $ flags{outputState = Just str}) + ("path format for program states like 'patch-state-%s.csv'") : + opt generic [] ["quality"] (flip ReqArg "PERCENTAGE" $ \str flags -> fmap (\x -> flags{quality = x}) $ @@ -190,10 +224,10 @@ opt generic [] ["maximum-absolute-angle"] -- "max-abs-angle" (flip ReqArg "DEGREE" $ \str flags -> - fmap (\x -> flags{maximumAbsoluteAngle = x}) $ + fmap (\x -> flags{maximumAbsoluteAngle = Degree x}) $ parseNumber "maximum absolute angle" (0<=) "non-negative" str) (printf "Maximum absolute angle for test rotations, default: %f" - (maximumAbsoluteAngle defltOption)) : + (getDegree $ maximumAbsoluteAngle defltOption)) : opt generic [] ["number-angles"] -- "num-angles" (flip ReqArg "NATURAL" $ \str flags -> @@ -280,10 +314,10 @@ opt generic [] ["hint-angle"] (flip ReqArg "DEGREE" $ \str (image, args) -> - fmap (\x -> (image{angle = Just x}, args)) $ + fmap (\x -> (image{angle = Just (Degree 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)) : + maybe "automatic estimation" (show . getDegree) (angle defltImage)) : [] @@ -293,7 +327,6 @@ return (defltImage, args {inputs = (image,path) : inputs args}) - get :: Engine -> IO Args get engine = do let desc = map snd $ filter (EnumSet.get engine . fst) $ description desc @@ -307,6 +340,20 @@ when (lastImage /= defltImage) $ exitFailureMsg "unused trailing image options" - case inputs parsedArgs of + return parsedArgs + +images :: + (State.AngleCorrected angleCorr) => Args -> IO [State.Proposed angleCorr] +images args = do + xs <- maybe (return Vector.empty) State.read (state $ option args) + case Vector.toList xs ++ (reverse $ map proposedFromImage $ inputs args) of [] -> exitFailureMsg "no input files" - images -> return $ parsedArgs {inputs = reverse images} + imgs -> do + let duplicates = + Map.keys $ + Map.filter (ListHT.lengthAtLeast 2) $ Map.fromListWith (++) $ + map (\img -> (State.propPath img, [img])) imgs + when (not $ null duplicates) $ + IO.hPutStrLn IO.stderr $ + "duplicate file paths: " ++ List.intercalate ", " duplicates + return imgs diff --git a/src/State.hs b/src/State.hs new file mode 100644 --- /dev/null +++ b/src/State.hs @@ -0,0 +1,396 @@ +module State where + +import Degree (Degree(Degree)) + +import qualified Data.Csv as Csv +import qualified Data.HashMap.Strict as HashMap +import qualified Data.ByteString.Lazy.Char8 as BL +import qualified Data.ByteString.Char8 as B +import qualified Data.Traversable as Trav +import qualified Data.Foldable as Fold +import qualified Data.Vector as Vector +import qualified Data.List.HT as ListHT +import qualified Data.List as List +import qualified Data.Map as Map +import qualified Data.Set as Set +import Data.Traversable (traverse) +import Data.Foldable (foldMap) +import Data.Csv ((.=), (.:)) +import Data.Vector (Vector) +import Data.Monoid ((<>)) +import Data.Maybe (fromMaybe, catMaybes, isNothing) +import Data.Tuple.HT (mapFst, mapSnd, mapPair, swap) +import Data.Map (Map) + +import qualified Control.Monad.Exception.Synchronous as ME +import qualified Control.Functor.HT as FuncHT +import Control.Monad (when, join) +import Control.Applicative (pure, liftA2, (<$>), (<$), (<*>), empty) + +import Text.Printf (printf) + +import qualified System.IO as IO +import Prelude hiding (read) + + +newtype File = File FilePath + +data Angle = Angle FilePath (Degree Float) + +data Position = Position FilePath (Degree Float) (Float, Float) + +data Displacement = + Displacement FilePath FilePath (Maybe Relation) (Maybe (Float, Float)) + +data Rotated = + Rotated + (Maybe ((FilePath, FilePath), Maybe Relation)) + (Maybe ((Float, Float), (Float, Float))) + + +imageId, angleId, dAngleId, xId, yId :: B.ByteString +imageId = B.pack "Image" +angleId = B.pack "Angle" +dAngleId = B.pack "DAngle" +xId = B.pack "X" +yId = B.pack "Y" + +imageAId, imageBId, relationId, dxId, dyId :: B.ByteString +imageAId = B.pack "ImageA" +imageBId = B.pack "ImageB" +relationId = B.pack "Rel" +dxId = B.pack "DX" +dyId = B.pack "DY" + +xaId, yaId, xbId, ybId :: B.ByteString +xaId = B.pack "XA" +yaId = B.pack "YA" +xbId = B.pack "XB" +ybId = B.pack "YB" + + +instance Csv.ToNamedRecord File where + toNamedRecord (File path) = Csv.namedRecord [imageId .= path] +instance Csv.DefaultOrdered File where + headerOrder _ = Csv.header [imageId] + +instance Csv.ToNamedRecord Angle where + toNamedRecord (Angle path (Degree angle)) = + Csv.namedRecord [imageId .= path, angleId .= angle] +instance Csv.DefaultOrdered Angle where + headerOrder _ = Csv.header [imageId, angleId] + +data Relation = NonOverlapping | Overlapping + deriving (Eq, Ord, Enum) + +instance Csv.ToField Relation where + toField NonOverlapping = B.singleton '-' + toField Overlapping = B.singleton 'X' + +instance Csv.FromField Relation where + parseField bstr = + case B.unpack bstr of + "-" -> pure NonOverlapping + "X" -> pure Overlapping + _ -> empty + +unrelated :: Bool -> Relation +unrelated b = + if b + then State.NonOverlapping + else State.Overlapping + + +instance Csv.ToNamedRecord Position where + toNamedRecord (Position path (Degree angle) (x,y)) = + Csv.namedRecord [imageId .= path, angleId .= angle, xId .= x, yId .= y] +instance Csv.DefaultOrdered Position where + headerOrder _ = Csv.header [imageId, angleId, xId, yId] + +instance Csv.ToNamedRecord Displacement where + toNamedRecord (Displacement pathA pathB rel disp) = + Csv.namedRecord + [imageAId .= pathA, imageBId .= pathB, + relationId .= rel, dxId .= fmap fst disp, dyId .= fmap snd disp] +instance Csv.DefaultOrdered Displacement where + headerOrder _ = Csv.header [imageAId, imageBId, relationId, dxId, dyId] + +instance Csv.FromNamedRecord Displacement where + parseNamedRecord m = + Displacement + <$> m .: imageAId + <*> m .: imageBId + <*> m .: relationId + <*> runCombinedParser (parseMaybePair (m .:# dxId) (m .:# dyId)) + + +instance Csv.ToNamedRecord Rotated where + toNamedRecord (Rotated pathsRel rot) = + let ((pathA, pathB), rel) = mapFst FuncHT.unzip $ FuncHT.unzip pathsRel + ((xa,ya),(xb,yb)) = + mapPair (FuncHT.unzip, FuncHT.unzip) $ FuncHT.unzip rot + in Csv.namedRecord $ + [imageAId .= pathA, imageBId .= pathB, relationId .= join rel, + xaId .= xa, yaId .= ya, xbId .= xb, ybId .= yb] +instance Csv.DefaultOrdered Rotated where + headerOrder _ = + Csv.header [imageAId, imageBId, relationId, xaId, yaId, xbId, ybId] + + +{- | +We use it like a @Compose Writer Csv.Parser@, +we do not have much use for the 'Applicative' combinators. +-} +data NamedParser a = NamedParser [B.ByteString] (Csv.Parser a) + +(.:#) :: Csv.FromField a => Csv.NamedRecord -> B.ByteString -> NamedParser a +(.:#) m name = NamedParser [name] (m .: name) + +enumerateNames :: [B.ByteString] -> String +enumerateNames = B.unpack . B.intercalate (B.pack ", ") + +parseMaybePair :: + NamedParser (Maybe a) -> NamedParser (Maybe b) -> + NamedParser (Maybe (a, b)) +parseMaybePair (NamedParser na pa) (NamedParser nb pb) = + let n = na++nb + in NamedParser n $ do + m <- liftA2 (,) pa pb + case m of + (Nothing, Nothing) -> return Nothing + (Just a, Just b) -> return $ Just (a,b) + _ -> + fail $ + printf "The columns %s must be all set or all empty." $ + enumerateNames n + +parseImpliedMaybe :: + NamedParser (Maybe a) -> NamedParser (Maybe b) -> + NamedParser (Maybe (a, Maybe b)) +parseImpliedMaybe (NamedParser na pa) (NamedParser nb pb) = + NamedParser (na++nb) $ do + m <- liftA2 (,) pa pb + case m of + (Nothing, Just _) -> + fail $ + printf "If the columns %s are set, then %s must be set as well." + (enumerateNames nb) (enumerateNames na) + (ma,mb) -> return $ flip (,) mb <$> ma + +runCombinedParser :: NamedParser a -> Csv.Parser a +runCombinedParser (NamedParser _ p) = p + + +instance Csv.FromNamedRecord Rotated where + parseNamedRecord m = + Rotated + <$> + (runCombinedParser $ + parseImpliedMaybe + (parseMaybePair (m .:# imageAId) (m .:# imageBId)) + (m .:# relationId)) + <*> + (runCombinedParser $ + parseMaybePair + (parseMaybePair (m .:# xaId) (m .:# yaId)) + (parseMaybePair (m .:# xbId) (m .:# ybId))) + +segmentRotated :: + [Rotated] -> + ME.Exceptional String + [((FilePath, FilePath), + (Maybe Relation, [((Float, Float), (Float, Float))]))] +segmentRotated = + (\(prefix, blocks) -> + ME.assert "leading relations without image header" (null prefix) + >> + (ME.Success $ + map + (\((mrot0,(paths,mrel)), rots) -> + (paths, + (mrel, + catMaybes $ mrot0 : map (\(Rotated _ mrot) -> mrot) rots))) + blocks)) + . + ListHT.segmentBeforeMaybe + (\(Rotated mPathRel mrot) -> (,) mrot <$> mPathRel) + +imagePairMap :: + [((FilePath, FilePath), a)] -> + ME.Exceptional String (Map (FilePath, FilePath) a) +imagePairMap = + Trav.sequence . + Map.fromListWithKey + (\(pathA,pathB) _ _ -> + ME.Exception $ printf "duplicate image pair: %s, %s" pathA pathB) . + map (mapSnd ME.Success) + +unmatchedImages :: + [FilePath] -> Map (FilePath, FilePath) a -> [FilePath] +unmatchedImages paths pairs = + Set.toList $ + Set.difference + (Set.fromList $ concatMap (\(pathA,pathB) -> [pathA,pathB]) $ + Map.keys pairs) + (Set.fromList paths) + +warnUnmatchedImages :: [FilePath] -> Map (FilePath, FilePath) a -> IO () +warnUnmatchedImages paths relations = do + let unmatched = unmatchedImages paths relations + when (not $ null unmatched) $ IO.hPutStrLn IO.stderr $ + "image in relations but not in the image list: " ++ + List.intercalate ", " unmatched + + + +write :: (Csv.ToNamedRecord a, Csv.DefaultOrdered a) => FilePath -> [a] -> IO () +write path = BL.writeFile path . Csv.encodeDefaultOrderedByName + + +writeDisplacement :: + FilePath -> [(ix, (FilePath, FilePath), Maybe (Float, Float))] -> IO () +writeDisplacement path = + write path . + map + (\(_, (pathA,pathB), md) -> + Displacement pathA pathB (Just $ unrelated $ isNothing md) md) + +writeRotated :: + FilePath -> + [(ix, (FilePath, FilePath), [((Float, Float), (Float, Float))])] -> IO () +writeRotated path = + write path . + concatMap + (\(_, paths, rots) -> + Rotated (Just (paths, Just $ unrelated $ null rots)) Nothing + : + map (\rot -> Rotated Nothing (Just rot)) rots) + + + +{- +(.:) accepts missing fields and (.:?) instead accepts missing columns. + +Warning: This implementation would also accept ill-formated cells: +(.:?) m field = (m.:field) <|> pure Nothing +-} +(.:?) :: + Csv.FromField a => Csv.NamedRecord -> B.ByteString -> Csv.Parser (Maybe a) +(.:?) m field = traverse Csv.parseField $ HashMap.lookup field m + +maybeCorrection :: + (Csv.FromField a) => a -> Maybe a -> Maybe (Maybe a) -> Maybe a +maybeCorrection zero a mda = fromMaybe (zero <$ a) mda + +parseAngle :: Csv.FromField a => Csv.NamedRecord -> Csv.Parser (Maybe a) +parseAngle m = join <$> m .:? angleId + +class Csv.DefaultOrdered angleCorr => AngleCorrected angleCorr where + autoAngleCorrection :: angleCorr + parseAngles :: + Csv.NamedRecord -> Csv.Parser (Maybe (Degree Float), angleCorr) + +newtype + AngleCorrection = + AngleCorrection {getAngleCorrection :: Maybe (Degree Float)} + +instance AngleCorrected AngleCorrection where + autoAngleCorrection = AngleCorrection Nothing + parseAngles m = + liftA2 + (\a mda -> (a, AngleCorrection $ maybeCorrection (Degree 0) a mda)) + (parseAngle m) (m .:? dAngleId) + +data NoAngleCorrection = NoAngleCorrection + +instance AngleCorrected NoAngleCorrection where + autoAngleCorrection = NoAngleCorrection + parseAngles m = liftA2 (,) (parseAngle m) (pure NoAngleCorrection) + + +data + Proposed angleCorr = + Proposed FilePath + (Maybe (Degree Float), angleCorr) + (Maybe Float, Maybe Float) + +propPath :: Proposed angleCorr -> FilePath +propPath (Proposed path _ _) = path + +propAngleCorr :: Proposed angleCorr -> angleCorr +propAngleCorr (Proposed _ (_,angleCorr) _) = angleCorr + +instance Csv.DefaultOrdered AngleCorrection where + headerOrder _ = Csv.header [dAngleId] + +instance Csv.DefaultOrdered NoAngleCorrection where + headerOrder _ = Csv.header [] + +instance + Csv.DefaultOrdered angleCorr => + Csv.DefaultOrdered (Proposed angleCorr) where + headerOrder p = + Csv.header [imageId, angleId] <> + Csv.headerOrder (propAngleCorr p) <> + Csv.header [xId, yId] + +instance + (AngleCorrected angleCorr) => + Csv.FromNamedRecord (Proposed angleCorr) where + parseNamedRecord m = + Proposed <$> m .: imageId + <*> parseAngles m + <*> liftA2 (,) (join <$> m .:? xId) (join <$> m .:? yId) + +read :: + (Csv.FromNamedRecord a, Csv.DefaultOrdered a) => + FilePath -> IO (Vector a) +read path = do + (headers, body) <- + either (ioError . userError) return . Csv.decodeByName + =<< BL.readFile path + let ignored = + Set.toList $ + Set.difference + (Set.fromList $ Fold.toList headers) + (Set.fromList $ Fold.toList $ + Csv.headerOrder $ Vector.head body) + when (not $ null ignored) $ IO.hPutStrLn IO.stderr $ + "ignore unknown columns: " ++ List.intercalate ", " (map B.unpack ignored) + return body + + +readDisplacement :: + [FilePath] -> FilePath -> + IO (Map (FilePath, FilePath) (Maybe Relation, Maybe (Float, Float))) +readDisplacement imagePaths path = do + relations <- + ME.switch (ioError . userError) return . + imagePairMap . + foldMap + (\(Displacement pathA pathB rel d) -> + ((pathA,pathB), (rel,d)) : + ((pathB,pathA), (rel, mapPair (negate,negate) <$> d)) : + []) + =<< read path + warnUnmatchedImages imagePaths relations + return relations + +readRotated :: + [FilePath] -> FilePath -> + IO (Map (FilePath, FilePath) + (Maybe Relation, [((Float, Float), (Float, Float))])) +readRotated imagePaths path = do + relationsPlain <- read path + relations <- + ME.switch (ioError . userError) return $ + imagePairMap . + concatMap + (\((pathA,pathB), (rel,ds)) -> + ((pathA,pathB), (rel,ds)) : + ((pathB,pathA), (rel, map swap ds)) : + []) + =<< + segmentRotated (Vector.toList relationsPlain) + warnUnmatchedImages imagePaths relations + return relations