Commit 323bf30e authored by Christoph Ruegge's avatar Christoph Ruegge
Browse files

Rename Image to Arr for reusability in exercise

parent 361b24a6
......@@ -65,20 +65,20 @@ import Data.Random.Normal (normalsIO')
Alias, for convenience.
\begin{code}
type Image r = Array r DIM2 Double
type Arr r = Array r DIM2 Double
\end{code}
Loading and saving of images -- rather uninteresting.
\begin{code}
loadImage :: FilePath -> IO (Image U)
loadImage :: FilePath -> IO (Arr U)
loadImage file = do
loaded <- readImageFromBMP file
case loaded of
Left err -> error $ show err
Right img -> computeP $ R.map doubleLuminanceOfRGB8 img
saveImage :: FilePath -> Image U -> IO ()
saveImage :: FilePath -> Arr U -> IO ()
saveImage file img = do
!m <- foldAllP max 0 img
!img' <- computeP $ R.map (rgb8OfGreyDouble . (/ m)) img
......@@ -92,8 +92,8 @@ done nicely using a typeclass.
\begin{code}
class Operator a where
evalOp :: a -> Image U -> Image D
adjointOp :: a -> Image U -> Image D
evalOp :: a -> Arr U -> Arr D
adjointOp :: a -> Arr U -> Arr D
\end{code}
## CG solver
......@@ -105,16 +105,16 @@ $$(T^* T + \alpha \mathbb{1}) x = T^* y.$$
The code tries to use as few operator evaluations as possible.
\begin{code}
data CGState = CGState { cgx :: Image U
, cgp :: Image U
, cgr :: Image U
data CGState = CGState { cgx :: Arr U
, cgp :: Arr U
, cgr :: Arr U
, cgr2 :: Double
}
cgreg :: Operator a => a -> Double -> Image U -> Image U -> [CGState]
cgreg :: Operator a => a -> Double -> Arr U -> Arr U -> [CGState]
cgreg op reg rhs initial =
runIdentity $ do
(res :: Image U) <- computeP $ rhs -^ evalOp op initial
(res :: Arr U) <- computeP $ rhs -^ evalOp op initial
rInit <- computeP $ adjointOp op res
r2Init <- normSquaredP rInit
return $ iterate cgStep (CGState initial rInit rInit r2Init)
......@@ -123,12 +123,12 @@ cgreg op reg rhs initial =
cgStep :: CGState -> CGState
cgStep (CGState x p r r2) =
runIdentity $ do
!(q :: Image U) <- computeP $ evalOp op p
!(q :: Arr U) <- computeP $ evalOp op p
!p2 <- normSquaredP p
!q2 <- normSquaredP q
let alpha = r2 / (q2 + reg*p2)
!x' <- computeP $ x +^ scale alpha p
!(s :: Image U) <- computeP $ adjointOp op q
!(s :: Arr U) <- computeP $ adjointOp op q
!r' <- computeP $ r -^ scale alpha (s +^ scale reg p)
!r2' <- normSquaredP r'
let beta = r2' / r2
......@@ -172,7 +172,7 @@ process xs f = foldM (\_ x -> f x >> return x) undefined xs
rule based on relative residuals.
\begin{code}
runCG :: Operator a => Double -> a -> Double -> Image U -> IO (Image U)
runCG :: Operator a => Double -> a -> Double -> Arr U -> IO (Arr U)
runCG tol op reg rhs = do
let initial = computeS $ fromFunction (extent rhs) (const 0)
let steps' = cgreg op reg rhs initial
......@@ -195,9 +195,10 @@ instance Operator StencilOp where
adjointOp = evalOp
mkKernel :: Int -> StencilOp
mkKernel n = StencilOp $ makeStencil2 n n $ \(Z:.i:.j) ->
if max (abs i) (abs j) <= n then Just x else Nothing
where x = 1 / fromIntegral (n*n)
mkKernel n = StencilOp $ makeStencil2 n n getElem
where getElem (Z:.i:.j)
| max (abs i) (abs j) <= n = Just (1 / fromIntegral (n*n))
| otherwise = Nothing
\end{code}
## Main
......@@ -229,8 +230,8 @@ main = do
mean <- (/ fromIntegral (size sh)) <$> sumAllP blurry
rnd <- normalsIO' (0, 0.15*mean)
let noiseArr = fromList sh $ take (size sh) rnd :: Image U
(noisy :: Image U) <- computeP $ blurry +^ noiseArr
let noiseArr = fromList sh $ take (size sh) rnd :: Arr U
(noisy :: Arr U) <- computeP $ blurry +^ noiseArr
saveImage "noisy.bmp" noisy
putStrLn "noisy reconstructed"
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment