Add the MRI code as it is after its first successful output.

diff --git a/src/MRI.hs b/src/MRI.hs
new file mode 100644 (file)
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 8a1efdb76de21e6b94421e63e6ebd07630ff5990..c9165d3e3dc71287cdcee00590a5e35142a6e40b 100644 (file)
--- 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