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

Inverse problems lecture

parent 3216ae96
...@@ -401,8 +401,10 @@ ...@@ -401,8 +401,10 @@
"name": "haskell" "name": "haskell"
}, },
"language_info": { "language_info": {
"codemirror_mode": "ihaskell",
"file_extension": ".hs",
"name": "haskell", "name": "haskell",
"version": "7.10.2" "version": "7.10.3"
} }
}, },
"nbformat": 4, "nbformat": 4,
......
{-# LANGUAGE BangPatterns #-}
import Data.Array.Repa
import Data.Array.Repa.Helpers (runCGreg, saveArrayAsBMP, loadArrayFromBMP)
import Data.Array.Repa.Stencil
import Data.Array.Repa.Stencil.Dim2
import Data.Random.Normal (normalsIO')
blur :: Stencil DIM2 Double -> Array U DIM2 Double -> Array PC5 DIM2 Double
blur stencil = mapStencil2 (BoundConst 0) stencil
blurAdjoint :: Stencil DIM2 Double -> Array U DIM2 Double -> Array PC5 DIM2 Double
blurAdjoint = blur
mkKernel :: Int -> Stencil DIM2 Double
mkKernel n = makeStencil2 n n getElem
where getElem (Z:.(!j):.(!k))
| max (abs j) (abs k) <= n = Just (1 / fromIntegral (n^(2::Int)))
| otherwise = Nothing
-- blur' :: Int -> Array U DIM2 Double -> Array D DIM2 Double
-- blur' n x =
main :: IO ()
main = do
img <- loadArrayFromBMP "grumpy.bmp"
saveArrayAsBMP "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
let reconstruct = runCGreg 1e-3 (blur kernel) (blurAdjoint kernel)
blurry <- computeP $ blur kernel img
saveArrayAsBMP "blurry.bmp" blurry
putStrLn "blurry reconstructed"
reconstruct 0 blurry >>= saveArrayAsBMP "blurry_reconstructed.bmp"
putStrLn "blurry regularized"
reconstruct 0.05 blurry >>= saveArrayAsBMP "blurry_regularized.bmp"
mean <- (/ fromIntegral (size sh)) <$> sumAllP blurry
rnd <- normalsIO' (0, 0.1*mean)
let noiseArr = fromListUnboxed sh $ take (size sh) rnd
noisy <- computeUnboxedP $ blurry +^ noiseArr
saveArrayAsBMP "noisy.bmp" noisy
putStrLn "noisy reconstructed"
reconstruct 0 noisy >>= saveArrayAsBMP "noisy_reconstructed.bmp"
putStrLn "noisy regularized"
reconstruct 0.5 noisy >>= saveArrayAsBMP "noisy_regularized.bmp"
Supports Markdown
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