OPTIMIZATIONS += -fexcess-precision
OPTIMIZATIONS += -fno-spec-constr-count
+REPA_INCLUDES := -ivendor/repa-head/repa/dist/build/
+REPA_INCLUDES += -ivendor/repa-head/repa-io/dist/build
+
GHC_OPTS += $(OPTIMIZATIONS)
GHC_OPTS += $(GHC_WARNINGS)
+GHC_OPTS += $(REPA_INCLUDES)
GHC_OPTS += -odir $(TMPDIR)
GHC_OPTS += -hidir $(TMPDIR)
GHC_OPTS += --make
GHC_OPTS += -threaded
GHC_OPTS += -o bin/${BIN}
-
.PHONY : test publish_doc doc src_html hlint
$(BIN): src/*.hs
n_cube n = (Repa.Z Repa.:. n Repa.:. n Repa.:. n)
trilinear :: Values3D
-trilinear = Repa.fromList (n_cube 3) $
+trilinear = Repa.fromListUnboxed (n_cube 3) $
flatten $
transpose_xz
trilinear_list
trilinear_zoom_2_list = [[[1, 3/2, 2, 5/2, 3], [1, 7/4, 5/2, 13/4, 4], [1, 2, 3, 4, 5], [1, 9/4, 7/2, 19/4, 6], [1, 5/2, 4, 11/2, 7]], [[1, 3/2, 2, 5/2, 3], [1, 15/8, 11/4, 29/8, 9/2], [1, 9/4, 7/2, 19/4, 6], [1, 21/8, 17/4, 47/8, 15/2], [1, 3, 5, 7, 9]], [[1, 3/2, 2, 5/2, 3], [1, 2, 3, 4, 5], [1, 5/2, 4, 11/2, 7], [1, 3, 5, 7, 9], [1, 7/2, 6, 17/2, 11]], [[1, 3/2, 2, 5/2, 3], [1, 17/8, 13/4, 35/8, 11/2], [1, 11/4, 9/2, 25/4, 8], [1, 27/8, 23/4, 65/8, 21/2], [1, 4, 7, 10, 13]], [[1, 3/2, 2, 5/2, 3], [1, 9/4, 7/2, 19/4, 6], [1, 3, 5, 7, 9], [1, 15/4, 13/2, 37/4, 12], [1, 9/2, 8, 23/2, 15]]]
trilinear_zoom_2 :: Values3D
-trilinear_zoom_2 = Repa.fromList (n_cube 6) $
+trilinear_zoom_2 = Repa.fromListUnboxed (n_cube 6) $
flatten $
transpose_xz
trilinear_zoom_2_list
trilinear9x9x9 :: Values3D
-trilinear9x9x9 = Repa.fromList (n_cube 9) $
+trilinear9x9x9 = Repa.fromListUnboxed (n_cube 9) $
flatten $
transpose_xz
trilinear9x9x9_list
-- No need to transpose_xz this one.
zeros :: Values3D
-zeros = Repa.fromList (n_cube 3) $ flatten zeros_list
+zeros = Repa.fromListUnboxed (n_cube 3) $ flatten zeros_list
-- | A 3x3x3 array of numbers, starting at (0,0,0) == 0 and counting
[ 24, 25, 26 ]]]
naturals :: Values3D
-naturals = Repa.fromList (n_cube 3) $
+naturals = Repa.fromListUnboxed (n_cube 3) $
flatten $
transpose_xz
naturals_list
-- | Used in at least one test where we need a 1x20x1 array.
naturals_1d :: Values3D
-naturals_1d = Repa.fromList twenty_vector (flatten naturals_1d_list)
+naturals_1d = Repa.fromListUnboxed twenty_vector (flatten naturals_1d_list)
zoom v3d scale_factor
| xsize == 0 || ysize == 0 || zsize == 0 = empty3d
| otherwise =
- R.force $ R.unsafeTraverse v3d transExtent f
+ R.compute $ R.unsafeTraverse v3d transExtent f
where
(xsize, ysize, zsize) = dims v3d
transExtent = zoom_shape scale_factor
+{-# LANGUAGE FlexibleContexts #-}
-- | The MRI module contains functionsd and definitions relevant (and
-- specific) to the MRI data files found at,
--
import Data.Word
import Data.Bits
import Data.Array.Repa as R
+import Data.Array.Repa.Repr.Unboxed as R
import Data.Array.Repa.IO.Binary as R
-import Data.Array.Repa.IO.ColorRamp as R
-import Data.Array.Repa.IO.BMP as R (writeComponentsToBMP)
+import Data.Array.Repa.Algorithms.ColorRamp as R
+import Data.Array.Repa.IO.BMP as R (writeImageToBMP)
import Values
-- | RawData is an array of words (16 bits), as contained in the MRI
-- data files.
-type RawData sh = Array sh Word16
+type RawData sh = Array U sh Word16
-- | A specialization of the 'RawData' type, to three dimensions.
type RawData3D = RawData DIM3
type RGB = (Word8, Word8, Word8)
-type ColorData sh = Array sh RGB
+type ColorData sh = Array U sh RGB
+{-# INLINE read_word16s #-}
read_word16s :: FilePath -> IO RawData3D
read_word16s path = do
arr <- R.readArrayFromStorableFile path mri_shape
- arr `deepSeqArray` return ()
- return arr
+ arr' <- now $ R.copy arr
+ return arr'
bracket :: Double -> Word16
shift xx 8 .|. (shift xx (-8) .&. 0x00ff)
-swap_bytes :: (Shape sh) => (RawData sh) -> (RawData sh)
+{-# INLINE swap_bytes #-}
+swap_bytes :: (Shape sh, Repr r Word16) => Array r sh Word16
+ -> Array D sh Word16
swap_bytes arr =
- R.force $ R.map flip16 arr
+ R.map flip16 arr
-bracket_array :: (Shape sh) => (Values sh) -> (RawData sh)
+bracket_array :: Shape sh => Values sh -> Array D sh Word16
bracket_array arr =
- R.force $ R.map bracket arr
+ R.map bracket arr
-round_array :: (Shape sh) => (Values sh) -> (RawData sh)
+{-# INLINE round_array #-}
+round_array :: Shape sh => Values sh -> Array D sh Word16
round_array arr =
- R.force $ R.map round arr
+ R.map round arr
-flip_y :: RawData3D -> RawData3D
+flip_y :: Repr r Word16 => Array r DIM3 Word16 -> Array D DIM3 Word16
flip_y arr =
- R.force $ R.traverse arr id
+ R.unsafeTraverse arr id
(\get (Z :. z :. y :. x) ->
get (Z :. z :. (mri_height - 1) - y :. x))
-flip_x :: RawData3D -> RawData3D
+flip_x :: Repr r Word16 => Array r DIM3 Word16 -> Array D DIM3 Word16
flip_x arr =
- R.force $ R.traverse arr id
+ R.unsafeTraverse arr id
(\get (Z :. z :. y :. x) ->
get (Z :. z :. y :. (mri_width - 1) - x))
+
+{-# INLINE write_word16s #-}
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
+ R.compute $ R.map (truncate_rgb . ramp_it) arr
where
ramp_it :: Double -> (Double, Double, Double)
ramp_it x =
b' = truncate (b * 255)
-
-z_slice :: Elt a => Int -> Array DIM3 a -> Array DIM2 a
+z_slice :: (R.Unbox a, Repr r a) => Int -> Array r DIM3 a -> Array D DIM2 a
z_slice n arr =
slice arr (Any :. n :. All :. All)
-
write_values_slice_to_bitmap :: Values2D -> FilePath -> IO ()
write_values_slice_to_bitmap v3d path =
- R.writeComponentsToBMP path routput goutput boutput
+ R.writeImageToBMP path colors
where
arr_bracketed = bracket_array v3d
- colors = values_to_colors $ R.map fromIntegral arr_bracketed
- routput = R.map (\(red, _, _) -> red) colors
- goutput = R.map (\(_, green, _) -> green) colors
- boutput = R.map (\(_, _, blue) -> blue) colors
+ colors = values_to_colors $ R.compute $ R.map fromIntegral arr_bracketed
let zoom_factor = (scale, scale, scale)
let out_file = "output.bin"
arr <- read_word16s in_file
- let arr' = swap_bytes arr
- let arrMRI = R.reshape mri_shape arr'
- let dbl_data = R.force $ R.map fromIntegral arrMRI
- let output = zoom dbl_data zoom_factor
- let word16_output = round_array output
+ let arr' = swap_bytes arr
+ let arrMRI = R.reshape mri_shape arr'
+ let dbl_data = R.compute $ R.map fromIntegral arrMRI
+ let output = zoom dbl_data zoom_factor
+ let word16_output = R.compute $ round_array output
write_word16s out_file word16_output
let scale = read s :: Int
let zoom_factor = (1, scale, scale)
let out_file = "output.bmp"
- arr <- read_word16s in_file
- let arr' = swap_bytes arr
- let arrInv = flip_x $ flip_y arr'
- let arrSlice = z_slice 50 arrInv
+ arr <- read_word16s in_file
+ let arrSlice = R.computeUnboxed $ z_slice 50 $ flip_x $ flip_y $ swap_bytes arr
let arrSlice' = R.reshape mri_slice3d arrSlice
- let dbl_data = R.map fromIntegral arrSlice'
- let output = zoom dbl_data zoom_factor
- write_values_slice_to_bitmap (z_slice 0 output) out_file
+
+ -- If zoom isn't being inlined we need to extract the slice before hand,
+ -- and convert it to the require formed.
+ let dbl_data = R.compute $ R.map fromIntegral arrSlice'
+ let output = zoom dbl_data zoom_factor
+ let arrSlice0 = R.computeUnboxed $ z_slice 0 output
+
+ write_values_slice_to_bitmap arrSlice0 out_file
where
import Data.Array.Repa (
- Array,
+ Array, U,
Z(..),
(:.)(..),
DIM2,
DIM3,
extent,
- fromList,
+ fromListUnboxed,
unsafeIndex,
)
import ScaleFactor (ScaleFactor)
-type Values sh = Array sh Double
-type Values2D = Values DIM2
-type Values3D = Values DIM3
+type Values sh = Array U sh Double
+type Values2D = Values DIM2
+type Values3D = Values DIM3
instance Arbitrary Values3D where
z_dim <- choose (1, 27)
elements <- vectorOf (x_dim * y_dim * z_dim) (arbitrary :: Gen Double)
let new_shape = (Z :. x_dim :. y_dim :. z_dim)
- let three_d = Data.Array.Repa.fromList new_shape elements
+ let three_d = Data.Array.Repa.fromListUnboxed new_shape elements
return three_d
empty3d :: Values3D
-empty3d = Data.Array.Repa.fromList (Z :. 0 :. 0 :. 0) []
+empty3d = Data.Array.Repa.fromListUnboxed (Z :. 0 :. 0 :. 0) []
dims :: Values3D -> (Int, Int, Int)
projects / spline3.git / commitdiff