From 2032a8190ba101305a0e10300dbb7ec37b7daea1 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Wed, 7 Sep 2011 15:58:40 -0400 Subject: [PATCH] Add the MRI code as it is after its first successful output. --- src/MRI.hs | 127 ++++++++++++++++++++++++++++++++++++++++++++++++++++ src/Main.hs | 68 ++++++++++++++++++++++------ 2 files changed, 182 insertions(+), 13 deletions(-) create mode 100644 src/MRI.hs diff --git a/src/MRI.hs b/src/MRI.hs new file mode 100644 index 0000000..9f97b08 --- /dev/null +++ b/src/MRI.hs @@ -0,0 +1,127 @@ +{-# LANGUAGE ScopedTypeVariables #-} + +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 + + +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_width :. mri_height) + +mri_lower_threshold :: Int +mri_lower_threshold = 1400 + +mri_upper_threshold :: Int +mri_upper_threshold = 2500 + + +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 :: IO RawData3D +read_word16s = do + arr <- R.readArrayFromStorableFile "../data/mri.bin" mri_shape + arr `deepSeqArray` return () + return arr + + +{-# INLINE bracket #-} +bracket :: Int -> Int -> Int -> Word16 +bracket low high x + | x < low = 0 + | x > high = 255 + | otherwise = truncate (r * 255) + where + numerator = fromIntegral (x - low) :: Double + denominator = fromIntegral (high - low) :: Double + r = numerator/denominator + + +{-# INLINE flip16 #-} +flip16 :: Word16 -> Word16 +flip16 xx = + shift xx 8 .|. (shift xx (-8) .&. 0x00ff) + + +bracket_array :: (Shape sh) => (RawData sh) -> (RawData sh) +bracket_array arr = + R.map (bracket mri_lower_threshold mri_upper_threshold . fromIntegral . flip16) arr + +flip_y :: RawData3D -> RawData3D +flip_y arr = + 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.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 = + -- dump the slice back as word16 + R.writeArrayToStorableFile + + +raw_data_to_color :: (Shape sh) => (RawData sh) -> (ColorData sh) +raw_data_to_color arr = + R.force $ R.map (truncate_rgb . ramp_it . fromIntegral) 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 + + diff --git a/src/Main.hs b/src/Main.hs index 8a1efdb..c9165d3 100644 --- a/src/Main.hs +++ b/src/Main.hs @@ -2,30 +2,72 @@ module Main where import Data.Array.Repa ( + DIM2, DIM3, Z(..), (:.)(..), + slice, + reshape, + Any(..), + All(..) ) +import qualified Data.Array.Repa as R (map) +import Data.Array.Repa.IO.BMP (writeComponentsToBMP) +import Data.Word import System.Environment (getArgs) import Grid (make_grid, zoom) +import MRI import Values (read_values_3d, write_values_1d) -mri_shape :: DIM3 -mri_shape = (Z :. 256 :. 256 :. 1) - - +mri_shape2d :: DIM2 +mri_shape2d = (Z :. 256*2 :. 256*2) +mri_shape3d :: DIM3 +mri_shape3d = (Z :. 256 :. 256 :. 1) main :: IO () main = do - args <- getArgs - let color = head args - let in_file = "./data/MRbrain.40." ++ color - let out_file = "MRbrain.40." ++ color ++ ".out" - mridata <- read_values_3d mri_shape in_file - - let g = make_grid 1 mridata - let output = zoom g (4,4,1) - write_values_1d output out_file +-- args <- getArgs +-- let color = head args +-- let in_file = "./data/MRbrain.40." ++ color + let out_file = "MRbrain.50.red.out" + arr <- read_word16s + let arrBrack = bracket_array arr + let arrInv = flip_x $ flip_y arrBrack + let arrSlice = slice arrInv (Any :. (50 :: Int) :. All :. All) + let arrColor = raw_data_to_color arrSlice + + let arrColor' = reshape mri_shape3d arrColor + let rdata = red_dbl_data arrColor' + let gdata = green_dbl_data arrColor' + let bdata = blue_dbl_data arrColor' + +-- mridata <- read_values_3d mri_shape in_file + + let gr = make_grid 1 rdata + let routput = zoom gr (2,2,1) + + let gg = make_grid 1 gdata + let goutput = zoom gg (2,2,1) + + let gb = make_grid 1 bdata + let boutput = zoom gb (2,2,1) + + let routput' = R.map double_to_word8 (reshape mri_shape2d routput) + let goutput' = R.map double_to_word8 (reshape mri_shape2d goutput) + let boutput' = R.map double_to_word8 (reshape mri_shape2d boutput) + +-- write_values_1d output out_file + writeComponentsToBMP out_file routput' goutput' boutput' + where + double_to_word8 :: Double -> Word8 + double_to_word8 x = + if x > 255 then + 255 :: Word8 + else + if x < 0 then + 0 :: Word8 + else + fromIntegral $ truncate x -- 2.44.2