X-Git-Url: https://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FMain.hs;h=c9165d3e3dc71287cdcee00590a5e35142a6e40b;hb=2032a8190ba101305a0e10300dbb7ec37b7daea1;hp=4467d174b94fa4d31d7273c47d9ef7afcdc5aa90;hpb=39994b55486bc843001a238def99a3e4091c5b37;p=spline3.git diff --git a/src/Main.hs b/src/Main.hs index 4467d17..c9165d3 100644 --- a/src/Main.hs +++ b/src/Main.hs @@ -1,90 +1,73 @@ module Main where -import Cube -import Face -import Grid -import Misc (flatten) -import Point -import RealFunction -import Tetrahedron -import ThreeDimensional +import Data.Array.Repa ( + DIM2, + DIM3, + Z(..), + (:.)(..), + slice, + reshape, + Any(..), + All(..) + ) +import qualified Data.Array.Repa as R (map) -trilinear :: [[[Double]]] -trilinear = [ [ [ 1, 2, 3 ], - [ 1, 3, 5 ], - [ 1, 4, 7 ] ], - [ [ 1, 2, 3 ], - [ 1, 4, 7 ], - [ 1, 6, 11 ] ], - [ [ 1, 2, 3 ], - [ 1, 5, 9 ], - [ 1, 8, 15 ]]] +import Data.Array.Repa.IO.BMP (writeComponentsToBMP) +import Data.Word +import System.Environment (getArgs) -zeros :: [[[Double]]] -zeros = [ [ [ 0, 0, 0 ], - [ 0, 0, 0 ], - [ 0, 0, 0 ] ], - -- - [ [ 0, 0, 0 ], - [ 0, 0, 0 ], - [ 0, 0, 0 ] ], - -- - [ [ 0, 0, 0 ], - [ 0, 0, 0 ], - [ 0, 0, 0 ]]] +import Grid (make_grid, zoom) +import MRI +import Values (read_values_3d, write_values_1d) -dummy :: [[[Double]]] -dummy = [ [ [ 0, 1, 2 ], - [ 3, 4, 5 ], - [ 6, 7, 8 ] ], - -- - [ [ 9, 10, 11 ], - [ 12, 13, 14 ], - [ 15, 16, 17 ] ], - -- - [ [ 18, 19, 20 ], - [ 21, 22, 23 ], - [ 24, 25, 26 ]]] +mri_shape2d :: DIM2 +mri_shape2d = (Z :. 256*2 :. 256*2) - -find_point_value :: RealFunction Point -find_point_value p = poly p - where - g0 = make_grid 1 trilinear - the_cubes = flatten (cubes g0) - good_cubes = filter ((flip contains_point) p) the_cubes - target_cube = head good_cubes - good_tets = filter ((flip contains_point) p) (tetrahedrons target_cube) - target_tetrahedron = head good_tets - poly = polynomial target_tetrahedron +mri_shape3d :: DIM3 +mri_shape3d = (Z :. 256 :. 256 :. 1) main :: IO () main = do - print $ find_point_value (0,0,0) - print $ find_point_value (1,0,0) - print $ find_point_value (2,0,0) - print $ find_point_value (0,1,0) - print $ find_point_value (1,1,0) - print $ find_point_value (2,1,0) - print $ find_point_value (0,2,0) - print $ find_point_value (1,2,0) - print $ find_point_value (2,2,0) - print $ find_point_value (0,0,1) - print $ find_point_value (1,0,1) - print $ find_point_value (2,0,1) - print $ find_point_value (0,1,1) - print $ find_point_value (1,1,1) - print $ find_point_value (2,1,1) - print $ find_point_value (0,2,1) - print $ find_point_value (1,2,1) - print $ find_point_value (2,2,1) - print $ find_point_value (0,0,2) - print $ find_point_value (1,0,2) - print $ find_point_value (2,0,2) - print $ find_point_value (0,1,2) - print $ find_point_value (1,1,2) - print $ find_point_value (2,1,2) - print $ find_point_value (0,2,2) - print $ find_point_value (1,2,2) - print $ find_point_value (2,2,2) +-- 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