Commit f30a3228 authored by Christoph Ruegge's avatar Christoph Ruegge
Browse files

move inverse to lecture/ and solutions/, created exercise

parent beb640fe
%% Cell type:markdown id: tags:
# Image deconvolution
## Task
Given a blurred image $y$, find the original image $x$, where
$$T x = y,$$
and $T$ is a blurring operator, e.g. an operator assigns to each pixel an
average of all pixels in some neighbourhood.
Depending on the blurring operator, it is not clear in general whether this
problem has a solution at all. A frequent approach is to replace the problem by
$$\text{Minimize } {\lVert T x - y \rVert}^2 \text{ with respect to } x,$$
or equivalently
$$T^* T x = T^* y.$$
This equation can e.g. be solved using CG, but it turns out to be numerically
unstable: if the blurred image $y$ is distorted with some noise, the
reconstruction $x$ has significantly amplified noise.
## Regularization
A simple way to stablilize the problem is Tikhonov regularization, which
replaces the problem by
$$\text{Minimize } {\lVert T x - y \rVert}^2 + \alpha {\lVert x \rVert}^2 \text{
with respect to } x,$$
for some *regularization parameter* $\alpha > 0$ or equivalently
$$(T^* T + \alpha \mathbb{1}) x = T^* y.$$
$\alpha$ has to be chosen large enough to avoid the instabilities of the
problem, but small enough in order not to change the problem too much.
There are several ways to choose the parameter (semi-)automatically depending on
the blurred and distorted image $y$ and the noise level. Mathematically, one is
interested in the question under what circumstances and how fast the regularized
reconstruction converge to the true solution.
For this example, we will simply choose a regularization parameter manually.
## Imports and basics
%% Cell type:code id: tags:
``` haskell
{-# OPTIONS_GHC -O3 #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE ScopedTypeVariables #-}
import Control.Monad.Identity
import Data.Array.Repa as R hiding ((++))
import Data.Array.Repa.Eval
import Data.Array.Repa.Stencil
import Data.Array.Repa.Stencil.Dim2
import Data.Array.Repa.Algorithms.Pixel (doubleLuminanceOfRGB8, rgb8OfGreyDouble)
import Data.Array.Repa.IO.BMP (readImageFromBMP, writeImageToBMP)
import Data.Random.Normal (normalsIO')
```
%% Cell type:markdown id: tags:
Alias, for convenience.
%% Cell type:code id: tags:
``` haskell
type Image r = Array r DIM2 Double
```
%% Cell type:markdown id: tags:
Loading and saving of images -- rather uninteresting.
%% Cell type:code id: tags:
``` haskell
loadImage :: FilePath -> IO (Image 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 file img = do
!m <- foldAllP max 0 img
!img' <- computeP $ R.map (rgb8OfGreyDouble . (/ m)) img
writeImageToBMP file img'
```
%% Cell type:markdown id: tags:
## Operator type class
We need functions for both the blurring operator and its adjoint. This can be
done nicely using a typeclass.
%% Cell type:code id: tags:
``` haskell
class Operator a where
evalOp :: a -> Image U -> Image D
adjointOp :: a -> Image U -> Image D
```
%% Cell type:markdown id: tags:
## CG solver
The Repa CG solver for the regularized normal equation
$$(T^* T + \alpha \mathbb{1}) x = T^* y.$$
The code tries to use as few operator evaluations as possible.
%% Cell type:code id: tags:
``` haskell
data CGState = CGState { cgx :: Image U
, cgp :: Image U
, cgr :: Image U
, cgr2 :: Double
}
cgreg :: Operator a => a -> Double -> Image U -> Image U -> [CGState]
cgreg op reg rhs initial =
runIdentity $ do
(res :: Image U) <- computeP $ rhs -^ evalOp op initial
rInit <- computeP $ adjointOp op res
r2Init <- normSquaredP rInit
return $ iterate cgStep (CGState initial rInit rInit r2Init)
where normSquaredP = sumAllP . R.map (^(2::Int))
scale a = R.map (* a)
cgStep :: CGState -> CGState
cgStep (CGState x p r r2) =
runIdentity $ do
!(q :: Image 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
!r' <- computeP $ r -^ scale alpha (s +^ scale reg p)
!r2' <- normSquaredP r'
let beta = r2' / r2
!p' <- computeP $ r' +^ scale beta p
return $ CGState x' p' r' r2'
```
%% Cell type:markdown id: tags:
`cgreg` returns a lazy list of all iterates. `takeUntil` is uses to implement
the stopping rule; it is similar to `takeWhile`, but also returns the final
iterate.
%% Cell type:code id: tags:
``` haskell
takeUntil :: (a -> Bool) -> [a] -> [a]
takeUntil _ [] = []
takeUntil predicate (x:xs)
| predicate x = [x]
| otherwise = x : takeUntil predicate xs
```
%% Cell type:markdown id: tags:
We want to be able to output some information about the iterates, i.e. run an IO
action over the list of iterates, before returning the final one. One option
would be something like
```haskell
fmap last . forM iterates $ \it -> do [...]
```
However, this would perform the `fmap last` only *after* the IO actions are
done, therfore retaining the entire list of iterates in memory.
We want all but the final iterate to be garbage collected directly after
printing information. This can be done using a monadic fold, where we do not
actually accumulate a result but only return the last value.
%% Cell type:code id: tags:
``` haskell
process :: Monad m => [a] -> (a -> m ()) -> m a
process xs f = foldM (\_ x -> f x >> return x) undefined xs
```
%% Cell type:markdown id: tags:
`runCG` wrapper that calls `cg`, print information and implements the
stopping rule based on relative residuals.
%% Cell type:code id: tags:
``` haskell
runCG :: Operator a => Double -> a -> Double -> Image U -> IO (Image U)
runCG tol op reg rhs = do
let initial = computeS $ fromFunction (extent rhs) (const 0)
let steps' = cgreg op reg rhs initial
let r20 = cgr2 $ head steps'
let steps = takeUntil (\x -> sqrt (cgr2 x / r20) < tol) steps'
result <- process (zip [(1::Int)..] steps) $ \(n, cgs) ->
putStrLn $ show n ++ " " ++ show (sqrt $ cgr2 cgs / r20)
return $ cgx . snd $ result
```
%% Cell type:markdown id: tags:
## Blurring operator
The implementation of the averaging operator uses repa stencils.
%% Cell type:code id: tags:
``` haskell
data StencilOp = StencilOp { getStencil :: Stencil DIM2 Double }
instance Operator StencilOp where
evalOp = (delay .) . mapStencil2 (BoundConst 0) . getStencil
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)
```
%% Cell type:markdown id: tags:
## Main
The main function loads an image, blurs and distorts it, and performs
regularized and unregularized reconstructions.
%% Cell type:code id: tags:
``` haskell
main :: IO ()
main = do
img <- loadImage "grumpy.bmp"
saveImage "exact.bmp" img
let sh = extent img
-- Repa only supports stencils of size up to 7. Efficient implementations of
-- more general blurring (or other convolution) operators can be done using e.g.
-- FFTs.
let kernel = mkKernel 7
blurry <- computeP $ evalOp kernel img
saveImage "blurry.bmp" blurry
putStrLn "blurry reconstructed"
runCG 1e-3 kernel 0 blurry >>= saveImage "blurry_reconstructed.bmp"
putStrLn "blurry regularized"
runCG 1e-3 kernel 0.1 blurry >>= saveImage "blurry_regularized.bmp"
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
saveImage "noisy.bmp" noisy
putStrLn "noisy reconstructed"
runCG 1e-3 kernel 0 noisy >>= saveImage "noisy_reconstructed.bmp"
putStrLn "noisy regularized"
runCG 1e-3 kernel 0.1 noisy >>= saveImage "noisy_regularized.bmp"
```
%% Cell type:markdown id: tags:
# Exercise: Inverse Problems
The *Volterra operator* $V$ is defined by
$$V f(x) = \int_{x_0}^x f(x) \,dx.$$
A simple discrete approximation is
$$(V f)_i = h \sum_{j=0}^i f_j,$$
where $f_j = f(x_j)$ are the values of $f$ at equidistant sampling points $(x_j)_{j=1}^{N-1}$ with distance $h = x_1 - x_0$. As a matrix, this can simply be written as
$$V = h \begin{pmatrix}
1 & 0 & 0 & 0 & \ldots & 0 \\
1 & 1 & 0 & 0 & \ldots & 0 \\
1 & 1 & 1 & 0 & \ldots & 0 \\
\vdots & & & \ddots & & \\
1 & 1 & 1 & 1 & \ldots & 1
\end{pmatrix}$$
- Implement the Volterra operator with a parameter `h :: Double` as an instance of the `Operator` type class from the lecture, but with the `Arr` type alias replaced by
```haskell
type Arr r = Array r DIM1 Double
```
Though it is possible to do this using a matrix formulation and e.g. `mmultS`, this unneccessarily inefficient. A better possibility is to note the connection between the Volterra operator and Haskell's `scanl` function. Moreover, there is an implementation of `scanl` for unboxed `Vector`s from `Data.Vector.Unboxed` (*not* HMatrix's `Numeric.Container`) as well as conversion functions `toUnboxed` and `fromUnboxed` for Repa `U` arrays to and from these `Vector`s. Use these to implement the operator.
*Note*: The adjoint (i.e. transpose) of the operator can be found most easily by noting the similarity between the Volterra matrix above and its transpose: the adjoint simply proceeds in the opposite direction.
- Use the regularized CG implementation from the lecture to compute regularized *derivatives* of the functions
$$\begin{align}
f_1(x) &= e^{-x^2} \\
f_2(x) &= \text{sgn}(x) = \begin{cases} -1 & x < 0 \\ 1 & x \geq 0 \end{cases}
\end{align}$$
sampled at 200 points between -2 and 2 and distorted by random noise between -0.1 and 0.1. The CG tolerance should be set suffieciently low (e.g. 1e-6).
- Plot the results together with the noisy functions and, in the case of $f_1$, the analytical derivative.
{"nbformat_minor":0,"nbformat":4,"cells":[{"metadata":{"hidden":false},"source":[" # Image deconvolution\n","\n"," ## Task\n","\n","Given a blurred image $y$, find the original image $x$, where\n","\n","$$T x = y,$$\n","\n","and $T$ is a blurring operator, e.g. an operator assigns to each pixel an\n","average of all pixels in some neighbourhood.\n","\n","Depending on the blurring operator, it is not clear in general whether this\n","problem has a solution at all. A frequent approach is to replace the problem by\n","\n","$$\\text{Minimize } {\\lVert T x - y \\rVert}^2 \\text{ with respect to } x,$$\n","\n","or equivalently\n","\n","$$T^* T x = T^* y.$$\n","\n","This equation can e.g. be solved using CG, but it turns out to be numerically\n","unstable: if the blurred image $y$ is distorted with some noise, the\n","reconstruction $x$ has significantly amplified noise.\n","\n"," ## Regularization\n","\n","A simple way to stablilize the problem is Tikhonov regularization, which\n","replaces the problem by\n","\n","$$\\text{Minimize } {\\lVert T x - y \\rVert}^2 + \\alpha {\\lVert x \\rVert}^2 \\text{\n","with respect to } x,$$\n","\n","for some *regularization parameter* $\\alpha > 0$ or equivalently\n","\n","$$(T^* T + \\alpha \\mathbb{1}) x = T^* y.$$\n","\n","$\\alpha$ has to be chosen large enough to avoid the instabilities of the\n","problem, but small enough in order not to change the problem too much.\n","\n","There are several ways to choose the parameter (semi-)automatically depending on\n","the blurred and distorted image $y$ and the noise level. Mathematically, one is\n","interested in the question under what circumstances and how fast the regularized\n","reconstruction converge to the true solution.\n","\n","For this example, we will simply choose a regularization parameter manually.\n","\n"," ## Imports and basics\n","\n"],"cell_type":"markdown"},{"execution_count":null,"outputs":[],"metadata":{"collapsed":false},"source":["{-# OPTIONS_GHC -O3 #-}\n","\n","{-# LANGUAGE BangPatterns #-}\n","{-# LANGUAGE ScopedTypeVariables #-}\n","\n","import Control.Monad.Identity\n","import Data.Array.Repa as R hiding ((++))\n","import Data.Array.Repa.Eval\n","import Data.Array.Repa.Stencil\n","import Data.Array.Repa.Stencil.Dim2\n","import Data.Array.Repa.Algorithms.Pixel (doubleLuminanceOfRGB8, rgb8OfGreyDouble)\n","import Data.Array.Repa.IO.BMP (readImageFromBMP, writeImageToBMP)\n","import Data.Random.Normal (normalsIO')\n"],"cell_type":"code"},{"metadata":{"hidden":false},"source":["Alias, for convenience.\n","\n"],"cell_type":"markdown"},{"execution_count":null,"outputs":[],"metadata":{"collapsed":false},"source":["type Arr r = Array r DIM2 Double\n"],"cell_type":"code"},{"metadata":{"hidden":false},"source":["Loading and saving of images -- rather uninteresting.\n","\n"],"cell_type":"markdown"},{"execution_count":null,"outputs":[],"metadata":{"collapsed":false},"source":["loadImage :: FilePath -> IO (Arr U)\n","loadImage file = do\n"," loaded <- readImageFromBMP file\n"," case loaded of\n"," Left err -> error $ show err\n"," Right img -> computeP $ R.map doubleLuminanceOfRGB8 img\n","\n","saveImage :: FilePath -> Arr U -> IO ()\n","saveImage file img = do\n"," !m <- foldAllP max 0 img\n"," !img' <- computeP $ R.map (rgb8OfGreyDouble . (/ m)) img\n"," writeImageToBMP file img'\n"],"cell_type":"code"},{"metadata":{"hidden":false},"source":[" ## Operator type class\n","\n","We need functions for both the blurring operator and its adjoint. This can be\n","done nicely using a typeclass.\n","\n"],"cell_type":"markdown"},{"execution_count":null,"outputs":[],"metadata":{"collapsed":false},"source":["class Operator a where\n"," evalOp :: a -> Arr U -> Arr D\n"," adjointOp :: a -> Arr U -> Arr D\n"],"cell_type":"code"},{"metadata":{"hidden":false},"source":[" ## CG solver\n","\n","The Repa CG solver for the regularized normal equation\n","\n","$$(T^* T + \\alpha \\mathbb{1}) x = T^* y.$$\n","\n","The code tries to use as few operator evaluations as possible.\n","\n"],"cell_type":"markdown"},{"execution_count":null,"outputs":[],"metadata":{"collapsed":false},"source":["data CGState = CGState { cgx :: Arr U\n"," , cgp :: Arr U\n"," , cgr :: Arr U\n"," , cgr2 :: Double\n"," }\n","\n","cgreg :: Operator a => a -> Double -> Arr U -> Arr U -> [CGState]\n","cgreg op reg rhs initial =\n"," runIdentity $ do\n"," (res :: Arr U) <- computeP $ rhs -^ evalOp op initial\n"," rInit <- computeP $ adjointOp op res\n"," r2Init <- normSquaredP rInit\n"," return $ iterate cgStep (CGState initial rInit rInit r2Init)\n"," where normSquaredP = sumAllP . R.map (^(2::Int))\n"," scale a = R.map (* a)\n"," cgStep :: CGState -> CGState\n"," cgStep (CGState x p r r2) =\n"," runIdentity $ do\n"," !(q :: Arr U) <- computeP $ evalOp op p\n"," !p2 <- normSquaredP p\n"," !q2 <- normSquaredP q\n"," let alpha = r2 / (q2 + reg*p2)\n"," !x' <- computeP $ x +^ scale alpha p\n"," !(s :: Arr U) <- computeP $ adjointOp op q\n"," !r' <- computeP $ r -^ scale alpha (s +^ scale reg p)\n"," !r2' <- normSquaredP r'\n"," let beta = r2' / r2\n"," !p' <- computeP $ r' +^ scale beta p\n"," return $ CGState x' p' r' r2'\n"],"cell_type":"code"},{"metadata":{"hidden":false},"source":["`cgreg` returns a lazy list of all iterates. `takeUntil` is uses to implement\n","the stopping rule; it is similar to `takeWhile`, but also returns the final\n","iterate.\n","\n"],"cell_type":"markdown"},{"execution_count":null,"outputs":[],"metadata":{"collapsed":false},"source":["takeUntil :: (a -> Bool) -> [a] -> [a]\n","takeUntil _ [] = []\n","takeUntil predicate (x:xs)\n"," | predicate x = [x]\n"," | otherwise = x : takeUntil predicate xs\n"],"cell_type":"code"},{"metadata":{"hidden":false},"source":["We want to be able to output some information about the iterates, i.e. run an IO\n","action over the list of iterates, before returning the final one. One option\n","would be something like\n","\n","```haskell\n","fmap last . forM iterates $ \\it -> do [...]\n","```\n","\n","However, this would perform the `fmap last` only *after* the IO actions are\n","done, therfore retaining the entire list of iterates in memory.\n","\n","We want all but the final iterate to be garbage collected directly after\n","printing information. This can be done using a monadic fold, where we do not\n","actually accumulate a result but only return the last value.\n","\n"],"cell_type":"markdown"},{"execution_count":null,"outputs":[],"metadata":{"collapsed":false},"source":["process :: Monad m => [a] -> (a -> m ()) -> m a\n","process xs f = foldM (\\_ x -> f x >> return x) undefined xs\n"],"cell_type":"code"},{"metadata":{"hidden":false},"source":["`runCG` wrapper that calls `cg`, print information and implements the stopping\n","rule based on relative residuals.\n","\n"],"cell_type":"markdown"},{"execution_count":null,"outputs":[],"metadata":{"collapsed":false},"source":["runCG :: Operator a => Double -> a -> Double -> Arr U -> IO (Arr U)\n","runCG tol op reg rhs = do\n"," let initial = computeS $ fromFunction (extent rhs) (const 0)\n"," let steps' = cgreg op reg rhs initial\n"," let r20 = cgr2 $ head steps'\n"," let steps = takeUntil (\\x -> sqrt (cgr2 x / r20) < tol) steps'\n"," result <- process (zip [(1::Int)..] steps) $ \\(n, cgs) ->\n"," putStrLn $ show n ++ \" \" ++ show (sqrt $ cgr2 cgs / r20)\n"," return $ cgx . snd $ result\n"],"cell_type":"code"},{"metadata":{"hidden":false},"source":[" ## Blurring operator\n","\n","The implementation of the averaging operator uses repa stencils.\n","\n"],"cell_type":"markdown"},{"execution_count":null,"outputs":[],"metadata":{"collapsed":false},"source":["data StencilOp = StencilOp { getStencil :: Stencil DIM2 Double }\n","\n","instance Operator StencilOp where\n"," evalOp = (delay .) . mapStencil2 (BoundConst 0) . getStencil\n"," adjointOp = evalOp\n","\n","mkKernel :: Int -> StencilOp\n","mkKernel n = StencilOp $ makeStencil2 n n getElem\n"," where getElem (Z:.i:.j)\n"," | max (abs i) (abs j) <= n = Just (1 / fromIntegral (n*n))\n"," | otherwise = Nothing\n"],"cell_type":"code"},{"metadata":{"hidden":false},"source":[" ## Main\n","\n","The main function loads an image, blurs and distorts it, and performs\n","regularized and unregularized reconstructions.\n","\n"],"cell_type":"markdown"},{"execution_count":null,"outputs":[],"metadata":{"collapsed":false},"source":["main :: IO ()\n","main = do\n"," img <- loadImage \"grumpy.bmp\"\n"," saveImage \"exact.bmp\" img\n","\n"," let sh = extent img\n","\n"," -- Repa only supports stencils of size up to 7. Efficient implementations of\n"," -- more general blurring (or other convolution) operators can be done using e.g.\n"," -- FFTs.\n"," let kernel = mkKernel 7\n","\n"," blurry <- computeP $ evalOp kernel img\n"," saveImage \"blurry.bmp\" blurry\n","\n"," putStrLn \"blurry reconstructed\"\n"," runCG 1e-3 kernel 0 blurry >>= saveImage \"blurry_reconstructed.bmp\"\n","\n"," putStrLn \"blurry regularized\"\n"," runCG 1e-3 kernel 0.1 blurry >>= saveImage \"blurry_regularized.bmp\"\n","\n"," mean <- (/ fromIntegral (size sh)) <$> sumAllP blurry\n"," rnd <- normalsIO' (0, 0.15*mean)\n"," let noiseArr = fromList sh $ take (size sh) rnd :: Arr U\n"," (noisy :: Arr U) <- computeP $ blurry +^ noiseArr\n"," saveImage \"noisy.bmp\" noisy\n","\n"," putStrLn \"noisy reconstructed\"\n"," runCG 1e-3 kernel 0 noisy >>= saveImage \"noisy_reconstructed.bmp\"\n","\n"," putStrLn \"noisy regularized\"\n"," runCG 1e-3 kernel 0.1 noisy >>= saveImage \"noisy_regularized.bmp\"\n"],"cell_type":"code"}],"metadata":{"language_info":{"name":"haskell","version":"7.10.2"},"kernelspec":{"display_name":"Haskell","name":"haskell","language":"haskell"}}}
\ No newline at end of file
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