--- /dev/null
+spline3: parallel implementation of the Sorokina/Zeilfelder spline scheme
+Copyright (C) 2024 Michael Orlitzky
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Affero General Public License as
+published by the Free Software Foundation, either version 3 of the
+License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU Affero General Public License for more details.
+
+You should have received a copy of the GNU Affero General Public License
+along with this program. If not, see <https://www.gnu.org/licenses/>.
-Wno-prepositive-qualified-module \
-Wno-missing-safe-haskell-mode \
-Wno-missing-deriving-strategies \
+ -Wno-missing-kind-signatures \
-rtsopts \
-threaded
cabal-version: 3.0
name: spline3
-version: 1.0.0
+version: 1.0.2
author: Michael Orlitzky
-maintainer: Michael Orlitzky <michael@orlitzky.com>
-homepage: http://michael.orlitzky.com/code/spline3.xhtml
+maintainer: Michael Orlitzky <michael@orlitzky.com>
+homepage: https://michael.orlitzky.com/code/spline3.xhtml
category: Math
-license: AGPL-3.0-only
+license: AGPL-3.0-or-later
license-file: doc/LICENSE
bug-reports: mailto:michael@orlitzky.com
synopsis:
- A parallel implementation of the Sorokina/Zeilfelder spline scheme.
+ Parallel implementation of the Sorokina/Zeilfelder spline scheme
description:
- Interpolate volumetric data according to "Local quasi-interpolation
- by cubic C^1 splines on type-6 tetrahedral partitions." The defaults
- are tailored to the MRI data contained in data/mri.bin from the
- Stanford volume data archive at <http://graphics.stanford.edu/data/voldata/>.
+ Interpolate volumetric data according to "Local
+ quasi-interpolation by cubic C^1 splines on type-6 tetrahedral
+ partitions." The defaults are tailored to the MRI data contained in
+ data/mri.bin from the Stanford volume data archive at
+ <https://graphics.stanford.edu/data/voldata/>.
For more information and examples, please see the project homepage at
- <http://michael.orlitzky.com/code/spline3.xhtml>.
+ <https://michael.orlitzky.com/code/spline3.xhtml>.
build-type: Simple
extra-source-files:
data/mri.bin
data/mri.info
data/announcement.txt
+ doc/COPYING
doc/README
-- Don't ship the references for copyright reasons.
--doc/references/*.pdf
source-repository head
type: git
- location: http://gitweb.michael.orlitzky.com/spline3.git
+ location: https://gitweb.michael.orlitzky.com/spline3.git
branch: master
assertEqual "auto-rotate equals manual rotate" ((ccwz . ccwz . cwy) expr1) expr2
where
expr1 =
- (3/8)*I +
- (1/12)*(T + R + L + D) +
- (1/64)*(FT + FR + FL + FD) +
- (7/48)*F +
- (1/48)*B +
- (1/96)*(RT + LD + LT + RD) +
- (1/192)*(BT + BR + BL + BD)
+ (3 / 8)*I +
+ (1 / 12)*(T + R + L + D) +
+ (1 / 64)*(FT + FR + FL + FD) +
+ (7 / 48)*F +
+ (1 / 48)*B +
+ (1 / 96)*(RT + LD + LT + RD) +
+ (1 / 192)*(BT + BR + BL + BD)
expr2 =
- (3/8)*I +
- (1/12)*(F + L + R + B) +
- (1/64)*(FT + LT + RT + BT) +
- (7/48)*T +
- (1/48)*D +
- (1/96)*(FL + BR + FR + BL) +
- (1/192)*(FD + LD + RD + BD)
+ (3 / 8)*I +
+ (1 / 12)*(F + L + R + B) +
+ (1 / 64)*(FT + LT + RT + BT) +
+ (7 / 48)*T +
+ (1 / 48)*D +
+ (1 / 96)*(FL + BR + FR + BL) +
+ (1 / 192)*(FD + LD + RD + BD)
-- | A list of all directions, sans the interior and composite types.
all_directions :: [Cardinal]
-- | The left-side boundary of the cube. See Sorokina and Zeilfelder,
-- p. 76.
xmin :: Cube -> Double
-xmin cube = (i' - 1/2)
+xmin cube = (i' - 1 / 2)
where
i' = fromIntegral (i cube) :: Double
-- | The right-side boundary of the cube. See Sorokina and Zeilfelder,
-- p. 76.
xmax :: Cube -> Double
-xmax cube = (i' + 1/2)
+xmax cube = (i' + 1 / 2)
where
i' = fromIntegral (i cube) :: Double
-- | The front boundary of the cube. See Sorokina and Zeilfelder,
-- p. 76.
ymin :: Cube -> Double
-ymin cube = (j' - 1/2)
+ymin cube = (j' - 1 / 2)
where
j' = fromIntegral (j cube) :: Double
-- | The back boundary of the cube. See Sorokina and Zeilfelder,
-- p. 76.
ymax :: Cube -> Double
-ymax cube = (j' + 1/2)
+ymax cube = (j' + 1 / 2)
where
j' = fromIntegral (j cube) :: Double
-- | The bottom boundary of the cube. See Sorokina and Zeilfelder,
-- p. 76.
zmin :: Cube -> Double
-zmin cube = (k' - 1/2)
+zmin cube = (k' - 1 / 2)
where
k' = fromIntegral (k cube) :: Double
-- | The top boundary of the cube. See Sorokina and Zeilfelder,
-- p. 76.
zmax :: Cube -> Double
-zmax cube = (k' + 1/2)
+zmax cube = (k' + 1 / 2)
where
k' = fromIntegral (k cube) :: Double
top_face :: Cube -> Face.Face
top_face cube = Face.Face v0' v1' v2' v3'
where
- delta = (1/2) :: Double
+ delta = (1 / 2) :: Double
cc = center cube
v0' = cc + ( Point delta (-delta) delta )
v1' = cc + ( Point delta delta delta )
back_face :: Cube -> Face.Face
back_face cube = Face.Face v0' v1' v2' v3'
where
- delta = (1/2) :: Double
+ delta = (1 / 2) :: Double
cc = center cube
v0' = cc + ( Point delta (-delta) (-delta) )
v1' = cc + ( Point delta delta (-delta) )
down_face :: Cube -> Face.Face
down_face cube = Face.Face v0' v1' v2' v3'
where
- delta = (1/2) :: Double
+ delta = (1 / 2) :: Double
cc = center cube
v0' = cc + ( Point (-delta) (-delta) (-delta) )
v1' = cc + ( Point (-delta) delta (-delta) )
front_face :: Cube -> Face.Face
front_face cube = Face.Face v0' v1' v2' v3'
where
- delta = (1/2) :: Double
+ delta = (1 / 2) :: Double
cc = center cube
v0' = cc + ( Point (-delta) (-delta) delta )
v1' = cc + ( Point (-delta) delta delta )
left_face :: Cube -> Face.Face
left_face cube = Face.Face v0' v1' v2' v3'
where
- delta = (1/2) :: Double
+ delta = (1 / 2) :: Double
cc = center cube
v0' = cc + ( Point delta (-delta) delta )
v1' = cc + ( Point (-delta) (-delta) delta )
right_face :: Cube -> Face.Face
right_face cube = Face.Face v0' v1' v2' v3'
where
- delta = (1/2) :: Double
+ delta = (1 / 2) :: Double
cc = center cube
v0' = cc + ( Point (-delta) delta delta)
v1' = cc + ( Point delta delta delta )
-- we'd expect the volume of each one to be 1/24.
prop_all_volumes_exact :: Cube -> Bool
prop_all_volumes_exact cube =
- and [volume t ~~= 1/24 | t <- tetrahedra cube]
+ and [volume t ~~= 1 / 24 | t <- tetrahedra cube]
-- | All tetrahedron should have their v0 located at the center of the
-- cube.
prop_c3000_identity :: Cube -> Bool
prop_c3000_identity cube =
c t0 3 0 0 0 ~= c t0 2 1 0 0 + c t6 2 1 0 0
- - ((c t0 2 0 1 0 + c t0 2 0 0 1)/ 2)
+ - ((c t0 2 0 1 0 + c t0 2 0 0 1) / 2)
where
t0 = tetrahedron cube 0
t6 = tetrahedron cube 6
prop_c2010_identity :: Cube -> Bool
prop_c2010_identity cube =
c t0 2 0 1 0 ~= c t0 1 1 1 0 + c t6 1 1 0 1
- - ((c t0 1 0 2 0 + c t0 1 0 1 1)/ 2)
+ - ((c t0 1 0 2 0 + c t0 1 0 1 1) / 2)
where
t0 = tetrahedron cube 0
t6 = tetrahedron cube 6
prop_c2001_identity :: Cube -> Bool
prop_c2001_identity cube =
c t0 2 0 0 1 ~= c t0 1 1 0 1 + c t6 1 1 1 0
- - ((c t0 1 0 0 2 + c t0 1 0 1 1)/ 2)
+ - ((c t0 1 0 0 2 + c t0 1 0 1 1) / 2)
where
t0 = tetrahedron cube 0
t6 = tetrahedron cube 6
prop_c1020_identity :: Cube -> Bool
prop_c1020_identity cube =
c t0 1 0 2 0 ~= c t0 0 1 2 0 + c t6 0 1 0 2
- - ((c t0 0 0 3 0 + c t0 0 0 2 1)/ 2)
+ - ((c t0 0 0 3 0 + c t0 0 0 2 1) / 2)
where
t0 = tetrahedron cube 0
t6 = tetrahedron cube 6
prop_c1002_identity :: Cube -> Bool
prop_c1002_identity cube =
c t0 1 0 0 2 ~= c t0 0 1 0 2 + c t6 0 1 2 0
- - ((c t0 0 0 0 3 + c t0 0 0 1 2)/ 2)
+ - ((c t0 0 0 0 3 + c t0 0 0 1 2) / 2)
where
t0 = tetrahedron cube 0
t6 = tetrahedron cube 6
prop_c1011_identity :: Cube -> Bool
prop_c1011_identity cube =
c t0 1 0 1 1 ~= c t0 0 1 1 1 + c t6 0 1 1 1 -
- ((c t0 0 0 1 2 + c t0 0 0 2 1)/ 2)
+ ((c t0 0 0 1 2 + c t0 0 0 2 1) / 2)
where
t0 = tetrahedron cube 0
t6 = tetrahedron cube 6
-- What gets computed for c2100 of t6.
expr1 = eval (function_values t6) $
- (3/8)*I +
- (1/12)*(T + R + L + D) +
- (1/64)*(FT + FR + FL + FD) +
- (7/48)*F +
- (1/48)*B +
- (1/96)*(RT + LD + LT + RD) +
- (1/192)*(BT + BR + BL + BD)
+ (3 / 8)*I +
+ (1 / 12)*(T + R + L + D) +
+ (1 / 64)*(FT + FR + FL + FD) +
+ (7 / 48)*F +
+ (1 / 48)*B +
+ (1 / 96)*(RT + LD + LT + RD) +
+ (1 / 192)*(BT + BR + BL + BD)
-- What should be computed for c2100 of t6.
expr2 = eval (function_values t0) $
- (3/8)*I +
- (1/12)*(F + R + L + B) +
- (1/64)*(FT + RT + LT + BT) +
- (7/48)*T +
- (1/48)*D +
- (1/96)*(FR + FL + BR + BL) +
- (1/192)*(FD + RD + LD + BD)
+ (3 / 8)*I +
+ (1 / 12)*(F + R + L + B) +
+ (1 / 64)*(FT + RT + LT + BT) +
+ (7 / 48)*T +
+ (1 / 48)*D +
+ (1 / 96)*(FR + FL + BR + BL) +
+ (1 / 192)*(FD + RD + LD + BD)
-- | We know what (c t6 2 1 0 0) should be from Sorokina and
t6 = tetrahedron cube 6
fvs = function_values t0
expected = eval fvs $
- (3/8)*I +
- (1/12)*(F + R + L + B) +
- (1/64)*(FT + RT + LT + BT) +
- (7/48)*T +
- (1/48)*D +
- (1/96)*(FR + FL + BR + BL) +
- (1/192)*(FD + RD + LD + BD)
+ (3 / 8)*I +
+ (1 / 12)*(F + R + L + B) +
+ (1 / 64)*(FT + RT + LT + BT) +
+ (7 / 48)*T +
+ (1 / 48)*D +
+ (1 / 96)*(FR + FL + BR + BL) +
+ (1 / 192)*(FD + RD + LD + BD)
-- Tests to check that the correct edges are incidental.
-- points (hi, hj, jk) with h = 0.5. We should be able to reproduce
-- this from splines based on the 3x3x3 trilinear.
trilinear_zoom_2_list :: [[[Double]]]
-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_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 = fromListUnboxed (n_cube 6) $
-- tetrahedron.
center :: Face -> Point
center (Face v0' v1' v2' v3') =
- (v0' + v1' + v2' + v3') `scale` (1/4)
+ (v0' + v1' + v2' + v3') `scale` (1 / 4)
-- object centered at (i,j,k)
make_values :: Values3D -> Int -> Int -> Int -> FunctionValues
make_values values !i !j !k =
- empty_values { front = value_at values (i-1) j k,
- back = value_at values (i+1) j k,
- left = value_at values i (j-1) k,
- right = value_at values i (j+1) k,
- down = value_at values i j (k-1),
- top = value_at values i j (k+1),
- front_left = value_at values (i-1) (j-1) k,
- front_right = value_at values (i-1) (j+1) k,
- front_down =value_at values (i-1) j (k-1),
- front_top = value_at values (i-1) j (k+1),
- back_left = value_at values (i+1) (j-1) k,
- back_right = value_at values (i+1) (j+1) k,
- back_down = value_at values (i+1) j (k-1),
- back_top = value_at values (i+1) j (k+1),
- left_down = value_at values i (j-1) (k-1),
- left_top = value_at values i (j-1) (k+1),
- right_down = value_at values i (j+1) (k-1),
- right_top = value_at values i (j+1) (k+1),
- front_left_down = value_at values (i-1) (j-1) (k-1),
- front_left_top = value_at values (i-1) (j-1) (k+1),
- front_right_down = value_at values (i-1) (j+1) (k-1),
- front_right_top = value_at values (i-1) (j+1) (k+1),
- back_left_down = value_at values (i+1) (j-1) (k-1),
- back_left_top = value_at values (i+1) (j-1) (k+1),
- back_right_down = value_at values (i+1) (j+1) (k-1),
- back_right_top = value_at values (i+1) (j+1) (k+1),
+ empty_values { front = value_at values (i - 1) j k,
+ back = value_at values (i + 1) j k,
+ left = value_at values i (j - 1) k,
+ right = value_at values i (j + 1) k,
+ down = value_at values i j (k - 1),
+ top = value_at values i j (k + 1),
+ front_left = value_at values (i - 1) (j - 1) k,
+ front_right = value_at values (i - 1) (j + 1) k,
+ front_down =value_at values (i - 1) j (k - 1),
+ front_top = value_at values (i - 1) j (k + 1),
+ back_left = value_at values (i + 1) (j - 1) k,
+ back_right = value_at values (i + 1) (j + 1) k,
+ back_down = value_at values (i + 1) j (k - 1),
+ back_top = value_at values (i + 1) j (k + 1),
+ left_down = value_at values i (j - 1) (k - 1),
+ left_top = value_at values i (j - 1) (k + 1),
+ right_down = value_at values i (j + 1) (k - 1),
+ right_top = value_at values i (j + 1) (k + 1),
+ front_left_down = value_at values (i - 1) (j - 1) (k - 1),
+ front_left_top = value_at values (i - 1) (j - 1) (k + 1),
+ front_right_down = value_at values (i - 1) (j + 1) (k - 1),
+ front_right_top = value_at values (i - 1) (j + 1) (k + 1),
+ back_left_down = value_at values (i + 1) (j - 1) (k - 1),
+ back_left_top = value_at values (i + 1) (j - 1) (k + 1),
+ back_right_down = value_at values (i + 1) (j + 1) (k - 1),
+ back_right_top = value_at values (i + 1) (j + 1) (k + 1),
interior = value_at values i j k }
-- | Takes a 'FunctionValues' and a function that transforms one
where
fvs = function_values g
fvs' = make_values fvs i j k
- tet_vol = (1/24) :: Double
+ tet_vol = (1 / 24) :: Double
-- The first cube along any axis covers (-1/2, 1/2). The second
| otherwise = (ceiling (coord + offset)) - 1
where
(xsize, ysize, zsize) = dims (function_values g)
- offset = (1/2) :: Double
+ offset = (1 / 2) :: Double
-- | Takes a 'Grid', and returns a 'Cube' containing the given 'Point'.
f p
where
g = Grid v3d
- offset = (1/2) :: Double
+ offset = (1 / 2) :: Double
m' = (fromIntegral m) / (fromIntegral sfx) - offset
n' = (fromIntegral n) / (fromIntegral sfy) - offset
o' = (fromIntegral o) / (fromIntegral sfz) - offset
test_trilinear_c0030 :: Assertion
test_trilinear_c0030 =
- assertAlmostEqual "c0030 is correct" (c t 0 0 3 0) (17/8)
+ assertAlmostEqual "c0030 is correct" (c t 0 0 3 0) (17 / 8)
test_trilinear_c0003 :: Assertion
test_trilinear_c0003 =
- assertAlmostEqual "c0003 is correct" (c t 0 0 0 3) (27/8)
+ assertAlmostEqual "c0003 is correct" (c t 0 0 0 3) (27 / 8)
test_trilinear_c0021 :: Assertion
test_trilinear_c0021 =
- assertAlmostEqual "c0021 is correct" (c t 0 0 2 1) (61/24)
+ assertAlmostEqual "c0021 is correct" (c t 0 0 2 1) (61 / 24)
test_trilinear_c0012 :: Assertion
test_trilinear_c0012 =
- assertAlmostEqual "c0012 is correct" (c t 0 0 1 2) (71/24)
+ assertAlmostEqual "c0012 is correct" (c t 0 0 1 2) (71 / 24)
test_trilinear_c0120 :: Assertion
test_trilinear_c0120 =
- assertAlmostEqual "c0120 is correct" (c t 0 1 2 0) (55/24)
+ assertAlmostEqual "c0120 is correct" (c t 0 1 2 0) (55 / 24)
test_trilinear_c0102 :: Assertion
test_trilinear_c0102 =
- assertAlmostEqual "c0102 is correct" (c t 0 1 0 2) (73/24)
+ assertAlmostEqual "c0102 is correct" (c t 0 1 0 2) (73 / 24)
test_trilinear_c0111 :: Assertion
test_trilinear_c0111 =
- assertAlmostEqual "c0111 is correct" (c t 0 1 1 1) (8/3)
+ assertAlmostEqual "c0111 is correct" (c t 0 1 1 1) (8 / 3)
test_trilinear_c0210 :: Assertion
test_trilinear_c0210 =
- assertAlmostEqual "c0210 is correct" (c t 0 2 1 0) (29/12)
+ assertAlmostEqual "c0210 is correct" (c t 0 2 1 0) (29 / 12)
test_trilinear_c0201 :: Assertion
test_trilinear_c0201 =
- assertAlmostEqual "c0201 is correct" (c t 0 2 0 1) (11/4)
+ assertAlmostEqual "c0201 is correct" (c t 0 2 0 1) (11 / 4)
test_trilinear_c0300 :: Assertion
test_trilinear_c0300 =
- assertAlmostEqual "c0300 is correct" (c t 0 3 0 0) (5/2)
+ assertAlmostEqual "c0300 is correct" (c t 0 3 0 0) (5 / 2)
test_trilinear_c1020 :: Assertion
test_trilinear_c1020 =
- assertAlmostEqual "c1020 is correct" (c t 1 0 2 0) (8/3)
+ assertAlmostEqual "c1020 is correct" (c t 1 0 2 0) (8 / 3)
test_trilinear_c1002 :: Assertion
test_trilinear_c1002 =
- assertAlmostEqual "c1002 is correct" (c t 1 0 0 2) (23/6)
+ assertAlmostEqual "c1002 is correct" (c t 1 0 0 2) (23 / 6)
test_trilinear_c1011 :: Assertion
test_trilinear_c1011 =
- assertAlmostEqual "c1011 is correct" (c t 1 0 1 1) (13/4)
+ assertAlmostEqual "c1011 is correct" (c t 1 0 1 1) (13 / 4)
test_trilinear_c1110 :: Assertion
test_trilinear_c1110 =
- assertAlmostEqual "c1110 is correct" (c t 1 1 1 0) (23/8)
+ assertAlmostEqual "c1110 is correct" (c t 1 1 1 0) (23 / 8)
test_trilinear_c1101 :: Assertion
test_trilinear_c1101 =
- assertAlmostEqual "c1101 is correct" (c t 1 1 0 1) (27/8)
+ assertAlmostEqual "c1101 is correct" (c t 1 1 0 1) (27 / 8)
test_trilinear_c1200 :: Assertion
test_trilinear_c1200 =
test_trilinear_c2010 :: Assertion
test_trilinear_c2010 =
- assertAlmostEqual "c2010 is correct" (c t 2 0 1 0) (10/3)
+ assertAlmostEqual "c2010 is correct" (c t 2 0 1 0) (10 / 3)
test_trilinear_c2001 :: Assertion
test_trilinear_c2001 =
test_trilinear_c2100 :: Assertion
test_trilinear_c2100 =
- assertAlmostEqual "c2100 is correct" (c t 2 1 0 0) (7/2)
+ assertAlmostEqual "c2100 is correct" (c t 2 1 0 0) (7 / 2)
test_trilinear_c3000 :: Assertion
test_trilinear_c3000 =
prop_cube_indices_never_go_out_of_bounds :: Grid -> Gen Bool
prop_cube_indices_never_go_out_of_bounds g =
do
- let coordmin = negate (1/2) :: Double
+ let coordmin = negate (1 / 2) :: Double
let (xsize, ysize, zsize) = dims $ function_values g
- let xmax = (fromIntegral xsize) - (1/2) :: Double
- let ymax = (fromIntegral ysize) - (1/2) :: Double
- let zmax = (fromIntegral zsize) - (1/2) :: Double
+ let xmax = (fromIntegral xsize) - (1 / 2) :: Double
+ let ymax = (fromIntegral ysize) - (1 / 2) :: Double
+ let zmax = (fromIntegral zsize) - (1 / 2) :: Double
x <- choose (coordmin, xmax)
y <- choose (coordmin, ymax)
instance Num Point where
- (Point x1 y1 z1) + (Point x2 y2 z2) = Point (x1+x2) (y1+y2) (z1+z2)
+ (Point x1 y1 z1) + (Point x2 y2 z2) = Point (x1 + x2) (y1 + y2) (z1 + z2)
(Point x1 y1 z1) - (Point x2 y2 z2) = Point (x1-x2) (y1-y2) (z1-z2)
(Point x1 y1 z1) * (Point x2 y2 z2) = Point (x1*x2) (y1*y2) (z1*z2)
abs (Point x y z) = Point (abs x) (abs y) (abs z)
{-# INLINE dot #-}
dot :: Point -> Point -> Double
dot (Point x1 y1 z1) (Point x2 y2 z2) =
- (x2 - x1)^(2::Int) + (y2 - y1)^(2::Int) + (z2 - z1)^(2::Int)
+ (x2 - x1) ^ (2::Int) + (y2 - y1) ^ (2::Int) + (z2 - z1) ^ (2::Int)
fexp :: (RealFunction a) -> Int -> (RealFunction a)
fexp f n
| n == 0 = const 1
- | otherwise = \x -> (f x)^n
+ | otherwise = \x -> (f x) ^ n
-- We just average the four vertices.
barycenter :: Tetrahedron -> Point
barycenter (Tetrahedron _ v0' v1' v2' v3' _) =
- (v0' + v1' + v2' + v3') `scale` (1/4)
+ (v0' + v1' + v2' + v3') `scale` (1 / 4)
coefficient :: Int -> Int -> Int -> Int -> Double
coefficient 0 0 3 0 =
- (1/8) * (i' + f + l' + t' + lt + fl + ft + flt)
+ (1 / 8) * (i' + f + l' + t' + lt + fl + ft + flt)
coefficient 0 0 0 3 =
- (1/8) * (i' + f + r + t' + rt + fr + ft + frt)
+ (1 / 8) * (i' + f + r + t' + rt + fr + ft + frt)
coefficient 0 0 2 1 =
- (5/24)*(i' + f + t' + ft) + (1/24)*(l' + fl + lt + flt)
+ (5 / 24)*(i' + f + t' + ft) + (1 / 24)*(l' + fl + lt + flt)
coefficient 0 0 1 2 =
- (5/24)*(i' + f + t' + ft) + (1/24)*(r + fr + rt + frt)
+ (5 / 24)*(i' + f + t' + ft) + (1 / 24)*(r + fr + rt + frt)
coefficient 0 1 2 0 =
- (5/24)*(i' + f) + (1/8)*(l' + t' + fl + ft)
- + (1/24)*(lt + flt)
+ (5 / 24)*(i' + f) + (1 / 8)*(l' + t' + fl + ft)
+ + (1 / 24)*(lt + flt)
coefficient 0 1 0 2 =
- (5/24)*(i' + f) + (1/8)*(r + t' + fr + ft)
- + (1/24)*(rt + frt)
+ (5 / 24)*(i' + f) + (1 / 8)*(r + t' + fr + ft)
+ + (1 / 24)*(rt + frt)
coefficient 0 1 1 1 =
- (13/48)*(i' + f) + (7/48)*(t' + ft)
- + (1/32)*(l' + r + fl + fr)
- + (1/96)*(lt + rt + flt + frt)
+ (13 / 48)*(i' + f) + (7 / 48)*(t' + ft)
+ + (1 / 32)*(l' + r + fl + fr)
+ + (1 / 96)*(lt + rt + flt + frt)
coefficient 0 2 1 0 =
- (13/48)*(i' + f) + (17/192)*(l' + t' + fl + ft)
- + (1/96)*(lt + flt)
- + (1/64)*(r + d + fr + fd)
- + (1/192)*(rt + ld + frt + fld)
+ (13 / 48)*(i' + f) + (17 / 192)*(l' + t' + fl + ft)
+ + (1 / 96)*(lt + flt)
+ + (1 / 64)*(r + d + fr + fd)
+ + (1 / 192)*(rt + ld + frt + fld)
coefficient 0 2 0 1 =
- (13/48)*(i' + f) + (17/192)*(r + t' + fr + ft)
- + (1/96)*(rt + frt)
- + (1/64)*(l' + d + fl + fd)
- + (1/192)*(rd + lt + flt + frd)
+ (13 / 48)*(i' + f) + (17 / 192)*(r + t' + fr + ft)
+ + (1 / 96)*(rt + frt)
+ + (1 / 64)*(l' + d + fl + fd)
+ + (1 / 192)*(rd + lt + flt + frd)
coefficient 0 3 0 0 =
- (13/48)*(i' + f) + (5/96)*(l' + r + t' + d + fl + fr + ft + fd)
- + (1/192)*(rt + rd + lt + ld + frt + frd + flt + fld)
+ (13 / 48)*(i' + f) + (5 / 96)*(l' + r + t' + d + fl + fr + ft + fd)
+ + (1 / 192)*(rt + rd + lt + ld + frt + frd + flt + fld)
coefficient 1 0 2 0 =
- (1/4)*i' + (1/6)*(f + l' + t')
- + (1/12)*(lt + fl + ft)
+ (1 / 4)*i' + (1 / 6)*(f + l' + t')
+ + (1 / 12)*(lt + fl + ft)
coefficient 1 0 0 2 =
- (1/4)*i' + (1/6)*(f + r + t')
- + (1/12)*(rt + fr + ft)
+ (1 / 4)*i' + (1 / 6)*(f + r + t')
+ + (1 / 12)*(rt + fr + ft)
coefficient 1 0 1 1 =
- (1/3)*i' + (5/24)*(f + t')
- + (1/12)*ft
- + (1/24)*(l' + r)
- + (1/48)*(lt + rt + fl + fr)
+ (1 / 3)*i' + (5 / 24)*(f + t')
+ + (1 / 12)*ft
+ + (1 / 24)*(l' + r)
+ + (1 / 48)*(lt + rt + fl + fr)
coefficient 1 1 1 0 =
- (1/3)*i' + (5/24)*f
- + (1/8)*(l' + t')
- + (5/96)*(fl + ft)
- + (1/48)*(d + r + lt)
- + (1/96)*(fd + ld + rt + fr)
+ (1 / 3)*i' + (5 / 24)*f
+ + (1 / 8)*(l' + t')
+ + (5 / 96)*(fl + ft)
+ + (1 / 48)*(d + r + lt)
+ + (1 / 96)*(fd + ld + rt + fr)
coefficient 1 1 0 1 =
- (1/3)*i' + (5/24)*f
- + (1/8)*(r + t')
- + (5/96)*(fr + ft)
- + (1/48)*(d + l' + rt)
- + (1/96)*(fd + lt + rd + fl)
+ (1 / 3)*i' + (5 / 24)*f
+ + (1 / 8)*(r + t')
+ + (5 / 96)*(fr + ft)
+ + (1 / 48)*(d + l' + rt)
+ + (1 / 96)*(fd + lt + rd + fl)
coefficient 1 2 0 0 =
- (1/3)*i' + (5/24)*f
- + (7/96)*(l' + r + t' + d)
- + (1/32)*(fl + fr + ft + fd)
- + (1/96)*(rt + rd + lt + ld)
+ (1 / 3)*i' + (5 / 24)*f
+ + (7 / 96)*(l' + r + t' + d)
+ + (1 / 32)*(fl + fr + ft + fd)
+ + (1 / 96)*(rt + rd + lt + ld)
coefficient 2 0 1 0 =
- (3/8)*i' + (7/48)*(f + t' + l')
- + (1/48)*(r + d + b + lt + fl + ft)
- + (1/96)*(rt + bt + fr + fd + ld + bl)
+ (3 / 8)*i' + (7 / 48)*(f + t' + l')
+ + (1 / 48)*(r + d + b + lt + fl + ft)
+ + (1 / 96)*(rt + bt + fr + fd + ld + bl)
coefficient 2 0 0 1 =
- (3/8)*i' + (7/48)*(f + t' + r)
- + (1/48)*(l' + d + b + rt + fr + ft)
- + (1/96)*(lt + bt + fl + fd + rd + br)
+ (3 / 8)*i' + (7 / 48)*(f + t' + r)
+ + (1 / 48)*(l' + d + b + rt + fr + ft)
+ + (1 / 96)*(lt + bt + fl + fd + rd + br)
coefficient 2 1 0 0 =
- (3/8)*i' + (1/12)*(t' + r + l' + d)
- + (1/64)*(ft + fr + fl + fd)
- + (7/48)*f
- + (1/48)*b
- + (1/96)*(rt + ld + lt + rd)
- + (1/192)*(bt + br + bl + bd)
+ (3 / 8)*i' + (1 / 12)*(t' + r + l' + d)
+ + (1 / 64)*(ft + fr + fl + fd)
+ + (7 / 48)*f
+ + (1 / 48)*b
+ + (1 / 96)*(rt + ld + lt + rd)
+ + (1 / 192)*(bt + br + bl + bd)
coefficient 3 0 0 0 =
- (3/8)*i' + (1/12)*(t' + f + l' + r + d + b)
- + (1/96)*(lt + fl + ft + rt + bt + fr)
- + (1/96)*(fd + ld + bd + br + rd + bl)
+ (3 / 8)*i' + (1 / 12)*(t' + f + l' + r + d + b)
+ + (1 / 96)*(lt + fl + ft + rt + bt + fr)
+ + (1 / 96)*(fd + ld + bd + br + rd + bl)
{-# INLINE volume #-}
volume :: Tetrahedron -> Double
volume (Tetrahedron _ v0' v1' v2' v3' _) =
- (1/6)*(det v0' v1' v2' v3')
+ (1 / 6)*(det v0' v1' v2' v3')
-- | The barycentric coordinates of a point with respect to v0.
{-# INLINE b0 #-}
volume1 :: Assertion
volume1 =
- assertEqual "volume is correct" True (vol ~= (-1/3))
+ assertEqual "volume is correct" True (vol ~= (-1 / 3))
where
vol = volume t
precomputed_volume = 0 }
volume1 :: Assertion
- volume1 = assertEqual "volume1 is correct" True (vol ~= (1/3))
+ volume1 = assertEqual "volume1 is correct" True (vol ~= (1 / 3))
where
vol = volume t
-- | Returns the domain point of t with indices i,j,k,l.
domain_point :: Tetrahedron -> Int -> Int -> Int -> Int -> Point
domain_point t i j k l =
- weighted_sum `scale` (1/3)
+ weighted_sum `scale` (1 / 3)
where
v0' = (v0 t) `scale` (fromIntegral i)
v1' = (v1 t) `scale` (fromIntegral j)
(volume t) > 0 ==>
c t 2 1 0 0 ~= (term1 - term2 + term3 - term4)
where
- term1 = (1/3)*(p t 0 3 0 0)
- term2 = (5/6)*(p t 3 0 0 0)
+ term1 = (1 / 3)*(p t 0 3 0 0)
+ term2 = (5 / 6)*(p t 3 0 0 0)
term3 = 3*(p t 2 1 0 0)
- term4 = (3/2)*(p t 1 2 0 0)
+ term4 = (3 / 2)*(p t 1 2 0 0)
-- | Given in Sorokina and Zeilfelder, p. 78.
prop_c1110_identity :: Tetrahedron -> Property
(volume t) > 0 ==>
c t 1 1 1 0 ~= (term1 + term2 - term3 - term4)
where
- term1 = (1/3)*((p t 3 0 0 0) + (p t 0 3 0 0) + (p t 0 0 3 0))
- term2 = (9/2)*(p t 1 1 1 0)
- term3 = (3/4)*((p t 2 1 0 0) + (p t 1 2 0 0) + (p t 2 0 1 0))
- term4 = (3/4)*((p t 1 0 2 0) + (p t 0 2 1 0) + (p t 0 1 2 0))
+ term1 = (1 / 3)*((p t 3 0 0 0) + (p t 0 3 0 0) + (p t 0 0 3 0))
+ term2 = (9 / 2)*(p t 1 1 1 0)
+ term3 = (3 / 4)*((p t 2 1 0 0) + (p t 1 2 0 0) + (p t 2 0 1 0))
+ term4 = (3 / 4)*((p t 1 0 2 0) + (p t 0 2 1 0) + (p t 0 1 2 0))
where
numerator = x - lower_threshold
denominator = upper_threshold - lower_threshold
- r = numerator/denominator
+ r = numerator / denominator
flip16 :: Word16 -> Word16
-#!/usr/bin/python2
+#!/usr/bin/python3
"""
Display volumetric data from the Stanford Volume Data Archive. The
def usage():
binary = sys.argv[0]
- print "Usage:", binary, "<data file>"
+ print(f"Usage: {binary} <data file>")
if len(sys.argv) < 2:
usage()
import numpy as np
data = np.fromfile(sys.argv[1], dtype='>u2')
-original_data_length = 7143424
-multiplier = cube_root(len(data) / original_data_length)
+original_data_length = 109*256*256
+multiplier = cube_root(len(data) // original_data_length)
data.shape = (109*multiplier, 256*multiplier, 256*multiplier)
data = data.T
# Display the data
-mlab.figure(bgcolor=(0, 0, 0), size=(1000, 1000))
+mlab.figure( bgcolor=(0,0,0), size=(1000,1000) )
+# Our data is scalar (grayscale), but not equally spaced in all
+# directions.
src = mlab.pipeline.scalar_field(data)
-
-# Our data is not equally spaced in all directions:
src.spacing = [1, 1, 1.5]
src.update_image_data = True
-v = mlab.pipeline.volume(src, vmax=2500, vmin=1400)
+
+mlab.pipeline.volume(src, vmin=1400, vmax=2500)
mlab.show()
projects / spline3.git / commitdiff