module MRI where import Data.Word import Data.Bits import Data.Array.Repa as R import Data.Array.Repa.IO.Binary as R import Data.Array.Repa.IO.ColorRamp as R import Values mri_depth :: Int mri_depth = 109 mri_width :: Int mri_width = 256 mri_height :: Int mri_height = 256 mri_shape :: DIM3 mri_shape = (Z :. mri_depth :. mri_height :. mri_width) mri_lower_threshold :: Double mri_lower_threshold = 1400 mri_upper_threshold :: Double mri_upper_threshold = 2500 mri_slice3d :: DIM3 mri_slice3d = (Z :. 1 :. mri_height :. mri_width) type RawData sh = Array sh Word16 type RawData3D = RawData DIM3 type RGB = (Word8, Word8, Word8) type ColorData sh = Array sh RGB rgb_to_dbl :: RGB -> (Double, Double, Double) rgb_to_dbl (x,y,z) = (fromIntegral x, fromIntegral y, fromIntegral z) read_word16s :: FilePath -> IO RawData3D read_word16s path = do arr <- R.readArrayFromStorableFile path mri_shape arr `deepSeqArray` return () return arr {-# INLINE bracket #-} bracket :: Double -> Word16 bracket x | x < mri_lower_threshold = 0 | x > mri_upper_threshold = 255 | otherwise = truncate (r * 255) where numerator = x - mri_lower_threshold denominator = mri_upper_threshold - mri_lower_threshold r = numerator/denominator {-# INLINE flip16 #-} flip16 :: Word16 -> Word16 flip16 xx = shift xx 8 .|. (shift xx (-8) .&. 0x00ff) swap_bytes :: (Shape sh) => (RawData sh) -> (RawData sh) swap_bytes arr = R.force $ R.map flip16 arr bracket_array :: (Shape sh) => (Values sh) -> (RawData sh) bracket_array arr = R.force $ R.map f arr where f = bracket flip_y :: RawData3D -> RawData3D flip_y arr = R.force $ R.traverse arr id (\get (Z :. z :. y :. x) -> get (Z :. z :. (mri_height - 1) - y :. x)) flip_x :: RawData3D -> RawData3D flip_x arr = R.force $ R.traverse arr id (\get (Z :. z :. y :. x) -> get (Z :. z :. y :. (mri_width - 1) - x)) write_word16s :: (Shape sh) => FilePath -> (RawData sh) -> IO () write_word16s = R.writeArrayToStorableFile values_to_colors :: (Shape sh) => (Values sh) -> (ColorData sh) values_to_colors arr = R.force $ R.map (truncate_rgb . ramp_it) arr where ramp_it :: Double -> (Double, Double, Double) ramp_it x = if x == 0 then (0, 0, 0) else rampColorHotToCold 0 255 x truncate_rgb :: (Double, Double, Double) -> (Word8, Word8, Word8) truncate_rgb (r, g, b) = (r', g', b') where r' = truncate (r * 255) g' = truncate (g * 255) b' = truncate (b * 255) red_dbl_data :: (Shape sh) => (ColorData sh) -> Array sh Double red_dbl_data = R.map (get_r . rgb_to_dbl) where get_r :: (Double, Double, Double) -> Double get_r (r, _, _) = r green_dbl_data :: (Shape sh) => (ColorData sh) -> Array sh Double green_dbl_data = R.map (get_g . rgb_to_dbl) where get_g :: (Double, Double, Double) -> Double get_g (_, g, _) = g blue_dbl_data :: (Shape sh) => (ColorData sh) -> Array sh Double blue_dbl_data = R.map (get_b . rgb_to_dbl) where get_b :: (Double, Double, Double) -> Double get_b (_, _, b) = b z_slice :: Elt a => Int -> Array DIM3 a -> Array DIM2 a z_slice n arr = slice arr (Any :. n :. All :. All)