src/MRI.hs

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