Tomography in Haskell
My task was to make an computer tomography in haskell. First, some explenations how tomographs works:
Yeah, now we will put some Xrays through it. What we will see distribution of material blocking Xrays in our object (you can think about it as a “density” for simplicity, but it’s not that). For example:
Here we’ve got a projection of our ball in angle of 0 radians. We must iterate that method from 0 to π radians. Because our object is a ball, and it’s in the center, all our projections will look like this. Whole table of projections is called a “sinogram”. With that we can recreate the inside of our patient without actually cut him! Let’s look how our reconstruction will looks like:
From two projections we can figure out that our object looks something like a square and is placed in center of our image. Imagine having 500 projections! Making this projections is done by Radon Transform. Unfortunately this will not works so good. Thats because pixels in center are in “focus”, things outside the center will tend to be blury. Look at the explenation:
That will be the “classic tomography”. But that is the relict of the past. I will focus on “computed tomography”. Using some filters we will try to fix that. In this case I use the ramp filter. How it works?
First, we take our projection and by fft we move it to frequency domain. Then we multiply that by our filter. It makes the lower frequencies less significant in our analysis and also removes the highest. Then, thankfully to inverse FFT we go back to our “time” domain. You can do it only in time domain by convolution, but doing it with fft is much faster (but harder, and in real world tomography as far as I now convolution is used).
Let’s jump into implementation. I use here mainly 2 libraries:
Data.Array.Repa for all array work and
Codec.Picture for loading and saving images.
Load the image
1 2 3 4 5 6 7 8 9 10 11 12 rgbToGrey :: PixelRGB8 -> Double rgbToGrey (PixelRGB8 r g b) = let rgb = P.zipWith (*) [0.2126, 0.7152, 0.0722] $ P.map (fromIntegral) [r, g, b] in fromInteger . round . sum $ rgb :: Double imageToArray :: String -> IO (Int, Int, Array U DIM2 Double) imageToArray fname = do Right imgRGB <- readImage fname let (imgGray@(Image w h _)) = convertRGB8 imgRGB arr = fromFunction (Z :. w :. h) (\(Z :. x :. y) -> rgbToGrey $ pixelAt imgGray x y) img <- computeUnboxedP arr return (w, h, img)
Yeah, it should be pretty self explanatory. Load image, turn it to grayscale and return the width, height and array with image itself.
After that we have to create our sinogram. I wrote these functions:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 getY a p r x = (x, round ((p - (fromIntegral (x - r) * (cos a))) / (sin a)) + r) getX a p r y = (round ((p - (fromIntegral (y - r) * (sin a))) / (cos a)) + r, y) filterCoords r (x, y) = ((fromIntegral x - r)**2 + (fromIntegral y - r)**2) < r**2 getLineAvg :: Double -> Double -> Int -> (Double -> Double -> Int -> Int -> (Int, Int)) -> Array U DIM2 Double -> Double getLineAvg a p w getCoord img = let pixelList = P.map (getCoord a p (div w 2)) [0..w-1] r = (fromIntegral w) / 2 pixelList' = filter (filterCoords r) pixelList len = length pixelList' ret = (sum $ P.map (\(x, y) -> img ! (Z :. x :. y)) pixelList') / fromIntegral len in if isNaN ret then 0 else ret getDetectorValue :: Double -> Double -> Int -> Array U DIM2 Double -> Double getDetectorValue a p w img | a < pi/4 || a > (3 * pi)/4 = getLineAvg a p w getX img | otherwise = getLineAvg a p w getY img
I operate here on normal form of the equation of a line: ` y * sin a + x * cos a = p’ so:
a- is an angle
p- is distance from (0, 0) - it’s the center of the image
w- the width of the image
I map getDetectorValue across all angles I want to scan my object, and then map it with distance of detector line from center of image. If think that’s interesting how easy is to write in Haskell functions with pattern matching.
After that it’s time to normalize our array.
1 2 3 4 5 6 normalize :: Monad m => Array U DIM2 Double -> m (Array D DIM2 Double) normalize arr = do minn <- foldAllP min 1 arr let arrmin = R.map (+ (-minn)) arr maxx <- foldAllP max 0 arrmin return $ R.map (/maxx) arrmin
Fold’s in haskell are like reduce. Give it an array and function by which it will combine all elements.
Now we will apply filter to our sinogram. We do it row by row. Here is our implementation of ramp filter:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 myfilter tresh v | v < tresh = fromIntegral v * 0.6 | otherwise = 0 rowFilter :: Monad m => Array D DIM1 Double -> m (Array D DIM1 Double) rowFilter row = do rowComplex <- computeP $ R.map (\x -> x :+ 0 ) row let fftF = fft $ rowComplex (Z :. w) = extent fftF barier = round (fromIntegral w/2) filtered = fromListUnboxed (Z :. w) $ (P.map (myfilter barier) [1..w]) fftFilt = R.zipWith (*) fftF filtered ifftF <- fmap ifft $ computeP fftFilt return $ R.map realPart ifftF
First we make our array an array of complex numbers. Then perform FFT and multiply it by our filter. After performing IFFT we are done. Unfortunately ‘Repa’ library doesn’t provide function to map through certain dimensions of array (at least I cannot find it), so I have to write my own.
1 2 3 4 5 6 7 8 mapRows :: Monad m => (Array D DIM1 Double -> m (Array D DIM1 Double)) -> Array D DIM2 Double -> m (Array U DIM2 Double) mapRows func array = do let (Z :. w :. h) = extent array rows <- mapM (\num -> do let row = getRow array num func row) [0,1..(w-1)] let hugeRow = toList $ foldr1 append rows return $ fromListUnboxed (Z :. w :. h) hugeRow
After that we’ve got filtered singoram.
Last thing what we’ve got to do is to perform Inverse Radon Transform.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 reconstruct :: Array D DIM2 Double -> Double -> Int -> IO () reconstruct img p orgW = do let (Z :. w :. h) = extent img angleStep = pi / fromIntegral h anglesList = takeWhile (<pi) [a * angleStep | a <- [0..]] orgWNum = fromIntegral orgW listOfIndicies (x, y) = let list = P.map (\a -> round $ (x * sin a + y * cos a)/p + orgWNum/2) anglesList list_zip = zip list [0..(h-1)] in [(a,b) | (a,b) <- list_zip, a >= 0, a < w] render (x, y) = let list = listOfIndicies (x, y) pixelList = P.map (\(p, h) -> img ! (Z :. p :. h)) list pixelSum = sum pixelList avg = pixelSum / (fromIntegral $ length pixelList) in if avg > 0 then avg else 0 let imageIndicies = [orgWNum/(-2)..orgWNum/2-1] img' = fromListUnboxed (Z :. orgW :. orgW) (P.map render [(a,b) | a <- imageIndicies , b <- imageIndicies]) img <- normalize img' writePng "res/reconstruct.png" $ generateImage (\x y -> dToPx (img ! (Z :. x :. y))) orgW orgW
listOfIndicies (x, y) returns an array with indicies of sinogram that we should consider in calulating the value of
(x, y) in our reconstructed image.
render (x, y) returns the value of pixel
(x, y) in this image. Let’s look at the results:
Pretty good, but it’s not the limit of this method. Compared to Python solutions I saw, Haskell implementation is pretty fast! Whole repository you can find on github. That’s all, thanks!