X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FMain.hs;h=951c9c36bf7a819b51f758de42a158484e00d90b;hb=fc0d3c47103269ed75788a87bb5f28ee70408c89;hp=8376ce954ff9f840d2b88f11b65a2393a283cd5f;hpb=89b8b6e94fcc944a1f4611811265f3c6217af850;p=spline3.git diff --git a/src/Main.hs b/src/Main.hs index 8376ce9..951c9c3 100644 --- a/src/Main.hs +++ b/src/Main.hs @@ -1,104 +1,115 @@ +{-# LANGUAGE RecordWildCards, DoAndIfThenElse #-} + module Main where -import Cube -import Face -import Grid -import Misc (flatten) -import Point -import RealFunction -import Tetrahedron -import ThreeDimensional - -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 ]]] - -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 ]]] - -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 ]]] - - -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 = good_cubes !! 0 - good_tets = filter ((flip contains_point) p) (tetrahedrons target_cube) - target_tetrahedron = good_tets !! 0 - poly = polynomial target_tetrahedron +import Data.Maybe (fromJust) +import Control.Monad (when) +import qualified Data.Array.Repa as R +import Data.Maybe (isJust) +import GHC.Conc (getNumProcessors, setNumCapabilities) +import System.IO (hPutStrLn, stderr) +import System.Exit (exitSuccess, exitWith, ExitCode(..)) + +import CommandLine (Args(..), apply_args) +import ExitCodes +import Grid (zoom) +import MRI ( + flip_x, + flip_y, + read_word16s, + round_array, + swap_bytes, + write_values_slice_to_bitmap, + write_word16s, + z_slice + ) + + +validate_args :: Args -> IO () +validate_args Args{..} = do + when (scale <= 0) $ do + hPutStrLn stderr "ERROR: scale must be greater than zero." + exitWith (ExitFailure exit_arg_not_positive) + + when (width <= 0) $ do + hPutStrLn stderr "ERROR: width must be greater than zero." + exitWith (ExitFailure exit_arg_not_positive) + + when (height <= 0) $ do + hPutStrLn stderr "ERROR: height must be greater than zero." + exitWith (ExitFailure exit_arg_not_positive) + + when (depth <= 0) $ do + hPutStrLn stderr "ERROR: depth must be greater than zero." + exitWith (ExitFailure exit_arg_not_positive) + + case slice of + Just s -> + when (s < 0 || s > depth) $ do + hPutStrLn stderr "ERROR: slice must be between zero and depth." + exitWith (ExitFailure exit_arg_out_of_bounds) + Nothing -> return () + main :: IO () main = do - putStrLn $ show $ find_point_value (0,0,0) - putStrLn $ show $ find_point_value (1,0,0) - putStrLn $ show $ find_point_value (2,0,0) - putStrLn $ show $ find_point_value (0,1,0) - putStrLn $ show $ find_point_value (1,1,0) - putStrLn $ show $ find_point_value (2,1,0) - putStrLn $ show $ find_point_value (0,2,0) - putStrLn $ show $ find_point_value (1,2,0) - putStrLn $ show $ find_point_value (2,2,0) - putStrLn $ show $ find_point_value (0,0,1) - putStrLn $ show $ find_point_value (1,0,1) - putStrLn $ show $ find_point_value (2,0,1) - putStrLn $ show $ find_point_value (0,1,1) - putStrLn $ show $ find_point_value (1,1,1) - putStrLn $ show $ find_point_value (2,1,1) - putStrLn $ show $ find_point_value (0,2,1) - putStrLn $ show $ find_point_value (1,2,1) - putStrLn $ show $ find_point_value (2,2,1) - putStrLn $ show $ find_point_value (0,0,2) - putStrLn $ show $ find_point_value (1,0,2) - putStrLn $ show $ find_point_value (2,0,2) - putStrLn $ show $ find_point_value (0,1,2) - putStrLn $ show $ find_point_value (1,1,2) - putStrLn $ show $ find_point_value (2,1,2) - putStrLn $ show $ find_point_value (0,2,2) - putStrLn $ show $ find_point_value (1,2,2) - putStrLn $ show $ find_point_value (2,2,2) - -- let g0 = make_grid 1 trilinear - -- let the_cubes = flatten (cubes g0) - -- putStrLn $ show $ the_cubes - -- let p = (2, 0, 0) - -- let target_cubes = filter ((flip contains_point) p) the_cubes - -- putStrLn $ show $ target_cubes - -- let target_cube = (take 1 target_cubes) !! 0 - -- putStrLn $ show $ target_cube - -- let target_tetrahedra = filter ((flip contains_point) p) (tetrahedrons target_cube) - -- let target_tetrahedron = (take 1 target_tetrahedra) !! 0 - -- putStrLn $ show $ target_tetrahedron - -- let poly = polynomial target_tetrahedron - -- putStrLn $ show $ poly - -- putStrLn $ show $ poly p + args@Args{..} <- apply_args + -- validate_args will simply exit if there's a problem. + validate_args args + + -- The first thing we do is set the number of processors. We get the + -- number of processors (cores) in the machine with + -- getNumProcessors, and set it with setNumCapabilities. This is so + -- we don't have to pass +RTS -Nfoo on the command line every time. + num_procs <- getNumProcessors + setNumCapabilities num_procs + + -- Determine whether we're doing 2d or 3d. If we're given a slice, + -- assume 2d. + let mri_shape = (R.Z R.:. depth R.:. height R.:. width) :: R.DIM3 + + if (isJust slice) then + main2d args mri_shape + else + main3d args mri_shape + + exitSuccess + + where + + + +main3d :: Args -> R.DIM3 -> IO () +main3d Args{..} mri_shape = do + let zoom_factor = (scale, scale, scale) + arr <- read_word16s input mri_shape + let arr' = swap_bytes arr + let arrMRI = R.reshape mri_shape arr' + dbl_data <- R.computeUnboxedP $ R.map fromIntegral arrMRI + raw_output <- zoom dbl_data zoom_factor + word16_output <- R.computeUnboxedP $ round_array raw_output + write_word16s output word16_output + + +main2d :: Args -> R.DIM3 -> IO () +main2d Args{..} mri_shape = do + let zoom_factor = (1, scale, scale) + arr <- read_word16s input mri_shape + arrSlice <- R.computeUnboxedP + $ z_slice (fromJust slice) + $ flip_x width + $ flip_y height + $ swap_bytes arr + let arrSlice' = R.reshape mri_slice3d arrSlice + + -- If zoom isn't being inlined we need to extract the slice before hand, + -- and convert it to the require formed. + dbl_data <- R.computeUnboxedP $ R.map fromIntegral arrSlice' + raw_output <- zoom dbl_data zoom_factor + arrSlice0 <- R.computeUnboxedP $ z_slice 0 raw_output + + write_values_slice_to_bitmap arrSlice0 output + where + mri_slice3d :: R.DIM3 + mri_slice3d = (R.Z R.:. 1 R.:. height R.:. width)