src/MRI.hs

   1 {-# LANGUAGE ScopedTypeVariables #-}
   2
   3 module MRI
   4 where
   5
   6 import Data.Word
   7 import Data.Bits
   8 import Data.Array.Repa                  as R
   9 import Data.Array.Repa.IO.Binary        as R
  10 import Data.Array.Repa.IO.ColorRamp     as R
  11
  12
  13 mri_depth :: Int
  14 mri_depth = 109
  15
  16 mri_width :: Int
  17 mri_width = 256
  18
  19 mri_height :: Int
  20 mri_height = 256
  21
  22 mri_shape :: DIM3
  23 mri_shape = (Z :. mri_depth :. mri_width  :. mri_height)
  24
  25 mri_lower_threshold :: Int
  26 mri_lower_threshold = 1400
  27
  28 mri_upper_threshold :: Int
  29 mri_upper_threshold = 2500
  30
  31
  32 type RawData sh = Array sh Word16
  33 type RawData3D = RawData DIM3
  34
  35 type RGB = (Word8, Word8, Word8)
  36 type ColorData sh = Array sh RGB
  37
  38 rgb_to_dbl :: RGB -> (Double, Double, Double)
  39 rgb_to_dbl (x,y,z) = (fromIntegral x, fromIntegral y, fromIntegral z)
  40
  41
  42 read_word16s :: IO RawData3D
  43 read_word16s = do
  44   arr <- R.readArrayFromStorableFile "../data/mri.bin" mri_shape
  45   arr `deepSeqArray` return ()
  46   return arr
  47
  48
  49 {-# INLINE bracket #-}
  50 bracket :: Int -> Int -> Int -> Word16
  51 bracket low high x
  52         | x < low      = 0
  53         | x > high     = 255
  54         | otherwise    = truncate (r * 255)
  55             where
  56               numerator = fromIntegral (x - low) :: Double
  57               denominator = fromIntegral (high - low) :: Double
  58               r = numerator/denominator
  59
  60
  61 {-# INLINE flip16 #-}
  62 flip16 :: Word16 -> Word16
  63 flip16 xx =
  64   shift xx 8 .|. (shift xx (-8) .&. 0x00ff)
  65
  66
  67 bracket_array :: (Shape sh) => (RawData sh) -> (RawData sh)
  68 bracket_array arr =
  69   R.map (bracket mri_lower_threshold mri_upper_threshold . fromIntegral . flip16) arr
  70
  71 flip_y :: RawData3D -> RawData3D
  72 flip_y arr =
  73   R.traverse arr id (\get (Z :. z :. y :. x) -> get (Z :. z :. (mri_height - 1) - y :. x))
  74
  75 flip_x :: RawData3D -> RawData3D
  76 flip_x arr =
  77   R.traverse arr id (\get (Z :. z :. y :. x) -> get (Z :. z :. y :. (mri_width - 1) - x))
  78
  79 write_word16s :: (Shape sh) => FilePath -> (RawData sh) -> IO ()
  80 write_word16s =
  81   -- dump the slice back as word16
  82   R.writeArrayToStorableFile
  83
  84
  85 raw_data_to_color :: (Shape sh) => (RawData sh) -> (ColorData sh)
  86 raw_data_to_color arr =
  87   R.force $ R.map (truncate_rgb . ramp_it . fromIntegral) arr
  88   where
  89     ramp_it :: Double -> (Double, Double, Double)
  90     ramp_it x =
  91       if x == 0 then
  92         (0, 0, 0)
  93       else
  94         rampColorHotToCold 0 255 x
  95
  96     truncate_rgb :: (Double, Double, Double) -> (Word8, Word8, Word8)
  97     truncate_rgb (r, g, b) =
  98       (r', g', b')
  99       where
 100         r' = truncate (r * 255)
 101         g' = truncate (g * 255)
 102         b' = truncate (b * 255)
 103
 104
 105 red_dbl_data :: (Shape sh) => (ColorData sh) -> Array sh Double
 106 red_dbl_data =
 107   R.map (get_r . rgb_to_dbl)
 108   where
 109     get_r :: (Double, Double, Double) -> Double
 110     get_r (r, _, _) = r
 111
 112 green_dbl_data :: (Shape sh) => (ColorData sh) -> Array sh Double
 113 green_dbl_data =
 114   R.map (get_g . rgb_to_dbl)
 115   where
 116     get_g :: (Double, Double, Double) -> Double
 117     get_g (_, g, _) = g
 118
 119
 120 blue_dbl_data :: (Shape sh) => (ColorData sh) -> Array sh Double
 121 blue_dbl_data =
 122   R.map (get_b . rgb_to_dbl)
 123   where
 124     get_b :: (Double, Double, Double) -> Double
 125     get_b (_, _, b) = b
 126
 127