src/MRI.hs

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