src/MRI.hs

   1 {-# LANGUAGE FlexibleContexts #-}
   2 -- | The MRI module contains functionsd and definitions relevant (and
   3 --   specific) to the MRI data files found at,
   4 --
   5 --   <http://graphics.stanford.edu/data/voldata/>
   6 --
   7 module MRI (
   8   flip_x,
   9   flip_y,
  10   mri_shape,
  11   mri_slice3d,
  12   read_word16s,
  13   round_array,
  14   swap_bytes,
  15   write_values_slice_to_bitmap,
  16   write_word16s,
  17   z_slice
  18   )
  19 where
  20
  21 import Data.Word
  22 import Data.Bits
  23 import Data.Array.Repa                  as R
  24 import Data.Array.Repa.Repr.Unboxed     as R
  25 import Data.Array.Repa.IO.Binary        as R
  26 import Data.Array.Repa.Algorithms.ColorRamp     as R
  27 import Data.Array.Repa.IO.BMP           as R (writeImageToBMP)
  28
  29 import Values
  30
  31 mri_depth :: Int
  32 mri_depth = 109
  33
  34 mri_width :: Int
  35 mri_width = 256
  36
  37 mri_height :: Int
  38 mri_height = 256
  39
  40 mri_shape :: DIM3
  41 mri_shape = (Z :. mri_depth :. mri_height :. mri_width)
  42
  43 mri_lower_threshold :: Double
  44 mri_lower_threshold = 1400
  45
  46 mri_upper_threshold :: Double
  47 mri_upper_threshold = 2500
  48
  49 mri_slice3d :: DIM3
  50 mri_slice3d = (Z :. 1 :. mri_height  :. mri_width)
  51
  52 -- | RawData is an array of words (16 bits), as contained in the MRI
  53 --   data files.
  54 type RawData sh = Array U sh Word16
  55
  56 -- | A specialization of the 'RawData' type, to three dimensions.
  57 type RawData3D = RawData DIM3
  58
  59 type RGB = (Word8, Word8, Word8)
  60 type ColorData sh = Array U sh RGB
  61
  62
  63 {-# INLINE read_word16s #-}
  64 read_word16s :: FilePath -> IO RawData3D
  65 read_word16s path = do
  66   arr <- R.readArrayFromStorableFile path mri_shape
  67   arr'   <- now $ R.copy arr
  68   return arr'
  69
  70
  71 bracket :: Double -> Word16
  72 bracket x
  73         | x < mri_lower_threshold      = 0
  74         | x > mri_upper_threshold      = 255
  75         | otherwise                    = truncate (r * 255)
  76             where
  77               numerator = x - mri_lower_threshold
  78               denominator = mri_upper_threshold - mri_lower_threshold
  79               r = numerator/denominator
  80
  81
  82 flip16 :: Word16 -> Word16
  83 flip16 xx =
  84   shift xx 8 .|. (shift xx (-8) .&. 0x00ff)
  85
  86
  87 {-# INLINE swap_bytes #-}
  88 swap_bytes :: (Shape sh, Repr r Word16) => Array r sh Word16
  89              -> Array D sh Word16
  90 swap_bytes arr =
  91   R.map flip16 arr
  92
  93
  94 bracket_array :: Shape sh => Values sh -> Array D sh Word16
  95 bracket_array arr =
  96   R.map bracket arr
  97
  98
  99 {-# INLINE round_array #-}
 100 round_array :: Shape sh => Values sh -> Array D sh Word16
 101 round_array arr =
 102   R.map round arr
 103
 104
 105 flip_y :: Repr r Word16 => Array r DIM3 Word16 -> Array D DIM3 Word16
 106 flip_y arr =
 107   R.unsafeTraverse arr id
 108               (\get (Z :. z :. y :. x) ->
 109                    get (Z :. z :. (mri_height - 1) - y :. x))
 110
 111 flip_x :: Repr r Word16 => Array r DIM3 Word16 -> Array D DIM3 Word16
 112 flip_x arr =
 113   R.unsafeTraverse arr id
 114               (\get (Z :. z :. y :. x) ->
 115                    get (Z :. z :. y :. (mri_width - 1) - x))
 116
 117
 118 {-# INLINE write_word16s #-}
 119 write_word16s :: (Shape sh) => FilePath -> (RawData sh) -> IO ()
 120 write_word16s = R.writeArrayToStorableFile
 121
 122
 123
 124 values_to_colors :: (Shape sh) => (Values sh) -> (ColorData sh)
 125 values_to_colors arr =
 126   R.compute $ R.map (truncate_rgb . ramp_it) arr
 127   where
 128     ramp_it :: Double -> (Double, Double, Double)
 129     ramp_it x =
 130       if x == 0 then
 131         (0, 0, 0)
 132       else
 133         rampColorHotToCold 0 255 x
 134
 135     truncate_rgb :: (Double, Double, Double) -> (Word8, Word8, Word8)
 136     truncate_rgb (r, g, b) =
 137       (r', g', b')
 138       where
 139         r' = truncate (r * 255)
 140         g' = truncate (g * 255)
 141         b' = truncate (b * 255)
 142
 143
 144 z_slice :: (R.Unbox a, Repr r a) => Int -> Array r DIM3 a -> Array D DIM2 a
 145 z_slice n arr =
 146   slice arr (Any :. n :. All :. All)
 147
 148
 149 write_values_slice_to_bitmap :: Values2D -> FilePath -> IO ()
 150 write_values_slice_to_bitmap v3d path =
 151   R.writeImageToBMP path colors
 152   where
 153     arr_bracketed = bracket_array v3d
 154     colors = values_to_colors $ R.compute $ R.map fromIntegral arr_bracketed