From 050d9c654c5a4fd15bbade6a2baafc789bc2d94b Mon Sep 17 00:00:00 2001 From: Declan Gaylo Date: Tue, 29 Jul 2025 11:52:10 -0400 Subject: [PATCH 1/6] Moved different reconstruction methods to own files --- blobid/reconstruction/central_difference.py | 19 +++++ blobid/reconstruction/no_inversion.py | 88 ++++++++++++++++++++ blobid/reconstruction/reconstruction.py | 91 +-------------------- tests/reconstruction/test_reconstruction.py | 3 +- 4 files changed, 113 insertions(+), 88 deletions(-) create mode 100644 blobid/reconstruction/central_difference.py create mode 100644 blobid/reconstruction/no_inversion.py diff --git a/blobid/reconstruction/central_difference.py b/blobid/reconstruction/central_difference.py new file mode 100644 index 0000000..832a803 --- /dev/null +++ b/blobid/reconstruction/central_difference.py @@ -0,0 +1,19 @@ +""" +Implementation of simple central differencing for interface reconstruction +""" +import numpy as np +import numpy.typing as npt + + +def get_normals(f: np.ndarray, normals_dtype: npt.DTypeLike) -> np.ndarray: + """ + Calculate the sign of the interface normal using central difference + """ + + n = np.empty((3, f.shape[0]-2, f.shape[1]-2, f.shape[2]-2), dtype=normals_dtype) + + n[0] = np.sign(f[:-2, 1:-1, 1:-1] - f[2:, 1:-1, 1:-1]) + n[1] = np.sign(f[1:-1, :-2, 1:-1] - f[1:-1, 2:, 1:-1]) + n[2] = np.sign(f[1:-1, 1:-1, :-2] - f[1:-1, 1:-1, 2:]) + + return n diff --git a/blobid/reconstruction/no_inversion.py b/blobid/reconstruction/no_inversion.py new file mode 100644 index 0000000..fe9f648 --- /dev/null +++ b/blobid/reconstruction/no_inversion.py @@ -0,0 +1,88 @@ +""" +Implementation of No Inversion VOF Interface Reconstruction Algorithm (NIVIRA) + +Reference +--------- + +G. D. Weymouth & D. K.-P. Yue, "Conservative Volume-of-Fluid method for free-surface simulations on +Cartesian-grids," Journal of Computational Physics, vol. 229, pp. 2853--2865, 2010. +[10.1016/j.jcp.2009.12.018](https://doi.org/10.1016/j.jcp.2009.12.018) +""" +import warnings +import numpy as np +import numpy.typing as npt + +from .._numba_support import njit, numba_availible + + +def get_normals(f: np.ndarray, normals_dtype: npt.DTypeLike) -> np.ndarray: + """ + Calculate the sign of the interface normal using NIVIRA + """ + + # Numba only supports float32 and float64 + if numba_availible and (f.dtype not in (np.float32, np.float64)): + warnings.warn(f"void_fraction type {f.dtype.name} not supported, using float32") + f = f.astype(np.float32) + + # estimate dominant direction + d_dom = _get_dominant_direction(f) + + # calculate normals using 3x3x3 stencil described by Weymouth and Yue + n = np.empty((3, f.shape[0]-2, f.shape[1]-2, f.shape[2]-2), dtype=normals_dtype) + + for d in range(3): + n[d] = _normals_WY(f, d_dom, d, normals_dtype) + + return n + + +def _get_dominant_direction(f): + # estimate the normal using central difference + n_approx = np.empty((3, f.shape[0]-2, f.shape[1]-2, f.shape[2]-2), dtype=f.dtype) + n_approx[0] = f[:-2, 1:-1, 1:-1] - f[2:, 1:-1, 1:-1] + n_approx[1] = f[1:-1, :-2, 1:-1] - f[1:-1, 2:, 1:-1] + n_approx[2] = f[1:-1, 1:-1, :-2] - f[1:-1, 1:-1, 2:] + + # cells do not have a normal if empty, full, or undefined CD + invalid = (f[1:-1, 1:-1, 1:-1] == 0) | (f[1:-1, 1:-1, 1:-1] == 1) | np.all(n_approx == 0, axis=0) + + # calculate the dominant direction, use d_dom=-1 to indicate no valid normal + d_dom = np.where(invalid, -1, np.argmax(np.abs(n_approx), axis=0)) + + return d_dom + + +@njit +def _normals_WY(f, d_dom, d, dtype) -> np.ndarray: + n = np.zeros(d_dom.shape, dtype=dtype) + + for i, j, k in np.argwhere(d_dom != -1): + if d_dom[i, j, k] == d: + # extract column in direction d + block = np.swapaxes(f[i:i+3, j:j+3, k:k+3], 0, d) + + # use central difference + n[i, j, k] = np.sign(block[0, 1, 1] - block[2, 1, 1]) + + else: + # skip the dimension we don't care about + d_drop = 3 - d - d_dom[i, j, k] + block = np.take(f[i:i+3, j:j+3, k:k+3], indices=1, axis=d_drop) + + # sum over d_dom + height = np.sum(block, axis=1 if (d_dom[i, j, k] > d) else 0) + + # central difference of summed columns + n_tmp = (height[0] - height[2])/2 + + if abs(n_tmp) < 0.5: + n[i, j, k] = np.sign(n_tmp) + else: + # for steep interface, use one-sided differences + if n_tmp * f[i+1, j+1, k+1] >= 0: + n[i, j, k] = np.sign(height[1] - height[2]) + else: + n[i, j, k] = np.sign(height[0] - height[1]) + + return n diff --git a/blobid/reconstruction/reconstruction.py b/blobid/reconstruction/reconstruction.py index 46269a8..c07f32a 100644 --- a/blobid/reconstruction/reconstruction.py +++ b/blobid/reconstruction/reconstruction.py @@ -1,10 +1,10 @@ """ Calculate the sign of the interface normal based on a 3x3(x3) stencil """ -import warnings import numpy as np -from .._numba_support import njit, numba_availible +from .central_difference import get_normals as get_normals_CD +from .no_inversion import get_normals as get_normals_WY NORMALS_TYPE = np.int8 @@ -15,91 +15,8 @@ def normals(f: np.ndarray, normals_method: str) -> np.ndarray: """ match normals_method: case 'CD': - return normals_CD(f) + return get_normals_CD(f, NORMALS_TYPE) case 'WY': - return normals_WY(f) + return get_normals_WY(f, NORMALS_TYPE) case _: raise ValueError(f"normals_method '{normals_method}' is not supported") - - -def normals_CD(f: np.ndarray) -> np.ndarray: - """ - Calculate the sign of the interface normal using central difference - """ - - n = np.empty((3, f.shape[0]-2, f.shape[1]-2, f.shape[2]-2), dtype=NORMALS_TYPE) - - n[0] = np.sign(f[:-2, 1:-1, 1:-1] - f[2:, 1:-1, 1:-1]) - n[1] = np.sign(f[1:-1, :-2, 1:-1] - f[1:-1, 2:, 1:-1]) - n[2] = np.sign(f[1:-1, 1:-1, :-2] - f[1:-1, 1:-1, 2:]) - - return n - - -def normals_WY(f: np.ndarray) -> np.ndarray: - # Numba only supports float32 and float64 - if numba_availible and (f.dtype not in (np.float32, np.float64)): - warnings.warn(f"void_fraction type {f.dtype.name} not supported, using float32") - f = f.astype(np.float32) - - # estimate dominant direction - d_dom = _get_dominant_direction(f) - - # calculate normals using 3x3x3 stencil described by Weymouth and Yue - n = np.empty((3, f.shape[0]-2, f.shape[1]-2, f.shape[2]-2), dtype=NORMALS_TYPE) - - for d in range(3): - n[d] = _normals_WY(f, d_dom, d) - - return n - - -def _get_dominant_direction(f): - # estimate the normal using central difference - n_approx = np.empty((3, f.shape[0]-2, f.shape[1]-2, f.shape[2]-2), dtype=f.dtype) - n_approx[0] = f[:-2, 1:-1, 1:-1] - f[2:, 1:-1, 1:-1] - n_approx[1] = f[1:-1, :-2, 1:-1] - f[1:-1, 2:, 1:-1] - n_approx[2] = f[1:-1, 1:-1, :-2] - f[1:-1, 1:-1, 2:] - - # cells do not have a normal if empty, full, or undefined CD - invalid = (f[1:-1, 1:-1, 1:-1] == 0) | (f[1:-1, 1:-1, 1:-1] == 1) | np.all(n_approx == 0, axis=0) - - # calculate the dominant direction, use d_dom=-1 to indicate no valid normal - d_dom = np.where(invalid, -1, np.argmax(np.abs(n_approx), axis=0)) - - return d_dom - - -@njit -def _normals_WY(f, d_dom, d) -> np.ndarray: - n = np.zeros(d_dom.shape, dtype=NORMALS_TYPE) - - for i, j, k in np.argwhere(d_dom != -1): - if d_dom[i, j, k] == d: - # extract column in direction d - block = np.swapaxes(f[i:i+3, j:j+3, k:k+3], 0, d) - - # use central difference - n[i, j, k] = np.sign(block[0, 1, 1] - block[2, 1, 1]) - - else: - # skip the dimension we don't care about - d_drop = 3 - d - d_dom[i, j, k] - block = np.take(f[i:i+3, j:j+3, k:k+3], indices=1, axis=d_drop) - - # sum over d_dom - height = np.sum(block, axis=1 if (d_dom[i, j, k] > d) else 0) - - # central difference of summed columns - n_tmp = (height[0] - height[2])/2 - - if abs(n_tmp) < 0.5: - n[i, j, k] = np.sign(n_tmp) - else: - # for steep interface, use one-sided differences - if n_tmp * f[i+1, j+1, k+1] >= 0: - n[i, j, k] = np.sign(height[1] - height[2]) - else: - n[i, j, k] = np.sign(height[0] - height[1]) - - return n diff --git a/tests/reconstruction/test_reconstruction.py b/tests/reconstruction/test_reconstruction.py index d1f03fc..20b4bc3 100644 --- a/tests/reconstruction/test_reconstruction.py +++ b/tests/reconstruction/test_reconstruction.py @@ -4,7 +4,8 @@ import numpy as np from blobid._numba_support import numba_availible -from blobid.reconstruction.reconstruction import normals, NORMALS_TYPE, _get_dominant_direction +from blobid.reconstruction.reconstruction import normals, NORMALS_TYPE +from blobid.reconstruction.no_inversion import _get_dominant_direction def test_reconstruction_CD(): From 283e7fb815d57cac9587842899b5ff00d8262dbb Mon Sep 17 00:00:00 2001 From: Declan Gaylo Date: Tue, 29 Jul 2025 11:58:55 -0400 Subject: [PATCH 2/6] Rename normals -> get_normals --- blobid/__init__.py | 3 ++- blobid/reconstruction/__init__.py | 4 ++-- blobid/reconstruction/reconstruction.py | 2 +- blobid/solver.py | 2 +- tests/reconstruction/test_reconstruction.py | 20 ++++++++++---------- tests/test_boundaries.py | 6 +++--- 6 files changed, 19 insertions(+), 18 deletions(-) diff --git a/blobid/__init__.py b/blobid/__init__.py index f13abce..cf94de3 100644 --- a/blobid/__init__.py +++ b/blobid/__init__.py @@ -1,4 +1,5 @@ from .solver import get_labels from . import labeling +from . import reconstruction -__all__ = ['get_labels', 'labeling'] +__all__ = ['get_labels', 'labeling', 'reconstruction'] diff --git a/blobid/reconstruction/__init__.py b/blobid/reconstruction/__init__.py index bd984a5..9710e13 100644 --- a/blobid/reconstruction/__init__.py +++ b/blobid/reconstruction/__init__.py @@ -1,3 +1,3 @@ -from .reconstruction import normals +from .reconstruction import get_normals -__all__ = ['normals'] +__all__ = ['get_normals'] diff --git a/blobid/reconstruction/reconstruction.py b/blobid/reconstruction/reconstruction.py index c07f32a..08dce17 100644 --- a/blobid/reconstruction/reconstruction.py +++ b/blobid/reconstruction/reconstruction.py @@ -9,7 +9,7 @@ NORMALS_TYPE = np.int8 -def normals(f: np.ndarray, normals_method: str) -> np.ndarray: +def get_normals(f: np.ndarray, normals_method: str) -> np.ndarray: """ Calculate the sign of the interface normals using method specified by normals_method """ diff --git a/blobid/solver.py b/blobid/solver.py index 662e436..213d585 100644 --- a/blobid/solver.py +++ b/blobid/solver.py @@ -160,7 +160,7 @@ def get_labels( # calculate connectivity if use_normals: - normals = reconstruction.normals(domain.vof(padding=1), normals_method) + normals = reconstruction.get_normals(domain.vof(padding=1), normals_method) else: normals = None diff --git a/tests/reconstruction/test_reconstruction.py b/tests/reconstruction/test_reconstruction.py index 20b4bc3..026a1bd 100644 --- a/tests/reconstruction/test_reconstruction.py +++ b/tests/reconstruction/test_reconstruction.py @@ -4,14 +4,14 @@ import numpy as np from blobid._numba_support import numba_availible -from blobid.reconstruction.reconstruction import normals, NORMALS_TYPE +from blobid.reconstruction.reconstruction import get_normals, NORMALS_TYPE from blobid.reconstruction.no_inversion import _get_dominant_direction def test_reconstruction_CD(): # all zeros f = np.zeros((3, 3, 3)) - n = normals(f, 'CD').squeeze() + n = get_normals(f, 'CD').squeeze() assert n[0] == 0 assert n[1] == 0 assert n[2] == 0 @@ -21,14 +21,14 @@ def test_reconstruction_CD(): # all ones f = np.zeros((3, 3, 3)) - n = normals(f, 'CD').squeeze() + n = get_normals(f, 'CD').squeeze() assert n[0] == 0 assert n[1] == 0 assert n[2] == 0 for _ in range(1000): f = np.random.rand(3, 3, 3) - n = normals(f, 'CD').squeeze() + n = get_normals(f, 'CD').squeeze() assert n[0] == -np.sign(f[2, 1, 1]-f[0, 1, 1]) assert n[1] == -np.sign(f[1, 2, 1]-f[1, 0, 1]) @@ -38,7 +38,7 @@ def test_reconstruction_CD(): def test_reconstruction_WY(): # all zeros f = np.zeros((3, 3, 3)) - n = normals(f, 'WY').squeeze() + n = get_normals(f, 'WY').squeeze() assert n[0] == 0 assert n[1] == 0 assert n[2] == 0 @@ -48,7 +48,7 @@ def test_reconstruction_WY(): # all ones f = np.zeros((3, 3, 3)) - n = normals(f, 'WY').squeeze() + n = get_normals(f, 'WY').squeeze() assert n[0] == 0 assert n[1] == 0 assert n[2] == 0 @@ -62,7 +62,7 @@ def test_reconstruction_WY(): assert _get_dominant_direction(f)[0, 0, 0] == 0 - n = normals(f, 'WY').squeeze() + n = get_normals(f, 'WY').squeeze() assert n[0] > 0 assert n[1] > 0 assert n[2] == 0 @@ -74,10 +74,10 @@ def test_vof_type(): # float16 should give a warning if numba_availible: with pytest.warns(): - normals(f.astype(np.float16), 'WY') + get_normals(f.astype(np.float16), 'WY') # float32 and float64 should work without warning with warnings.catch_warnings(): warnings.simplefilter("error") - normals(f.astype(np.float32), 'WY') - normals(f.astype(np.float64), 'WY') + get_normals(f.astype(np.float32), 'WY') + get_normals(f.astype(np.float64), 'WY') diff --git a/tests/test_boundaries.py b/tests/test_boundaries.py index 3f4c47a..7254ada 100644 --- a/tests/test_boundaries.py +++ b/tests/test_boundaries.py @@ -4,7 +4,7 @@ from blobid._domain import _pad_array from blobid.labeling.database import LabelDatabase from blobid.labeling.labeling import _stitch_boundaries -from blobid.reconstruction import normals +from blobid.reconstruction import get_normals @pytest.fixture @@ -138,8 +138,8 @@ def test_normal_calculation_at_boundaries(vof_3d): f = _pad_array(vof_3d.copy(), per, 1, 1) # normals should be the same on periodic sides - n = normals(f, normals_method='CD') + n = get_normals(f, normals_method='CD') check_periodicity(np.moveaxis(n, 0, 3), per, 1) - n = normals(f, normals_method='WY') + n = get_normals(f, normals_method='WY') check_periodicity(np.moveaxis(n, 0, 3), per, 1) From 97829ee258f7ea69bab64999030542595a189388 Mon Sep 17 00:00:00 2001 From: Declan Gaylo Date: Tue, 29 Jul 2025 12:20:54 -0400 Subject: [PATCH 3/6] Documentation update --- blobid/labeling/labeling.py | 8 +++++ blobid/reconstruction/reconstruction.py | 44 ++++++++++++++++++++++--- blobid/solver.py | 23 +++---------- 3 files changed, 52 insertions(+), 23 deletions(-) diff --git a/blobid/labeling/labeling.py b/blobid/labeling/labeling.py index 18f462f..df4289d 100644 --- a/blobid/labeling/labeling.py +++ b/blobid/labeling/labeling.py @@ -33,6 +33,14 @@ def apply_ccl( labels : ndarray[ni, nj, nk] An array of type `label_type` with the same shape as `is_object`. + Notes + ----- + The labeling algorithm used is based on work by He, Chao & Suzuki.[1]_ + + + .. [1] L. He, Y. Chao & K. Suzuki, "A Linear-Time Two-Scan Labeling Algorithm," IEEE International Conference + on Image Processing, 2007. [10.1109/ICIP.2007.4379810](https://doi.org/10.1109/ICIP.2007.4379810) + """ # checks assert is_object.ndim == 3 diff --git a/blobid/reconstruction/reconstruction.py b/blobid/reconstruction/reconstruction.py index 08dce17..2287a71 100644 --- a/blobid/reconstruction/reconstruction.py +++ b/blobid/reconstruction/reconstruction.py @@ -9,14 +9,48 @@ NORMALS_TYPE = np.int8 -def get_normals(f: np.ndarray, normals_method: str) -> np.ndarray: - """ - Calculate the sign of the interface normals using method specified by normals_method +def get_normals(void_fraction: np.ndarray, normals_method: str = 'CD') -> np.ndarray: + r""" + Calculate the sign of the interface normals using method specified by `normals_method` + + Parameters + ---------- + void_fraction : ndarray[ni+2, nj+2, nk+2] + The void fractions for each grid cell. + normals_method: {'CD', 'WY'}, optional + Method used to calculate interface normals. Defaults to 'CD'. + + Returns + ------- + normals : ndarray[3, ni, nj, nk] + An integer array which returns the sign of the interface normal. + + Notes + ----- + + The signs of the interface normal :math:`\langle n_x, n_y, n_z\rangle` are calculated using the method set by + `normals_method`: + + - If `normals_method` is `CD`, central differencing is used. For example + .. math:: + \mathrm{sign}(n_x) = - \mathrm{sign}(f_{i+1} - f_{i-1}) + + - If `normals_method` is `WY`, normals are calculated using the No Inversion VOF Interface Reconstruction Algorithm + (NIVIRA) described by Weymouth and Yue.[1]_ `WY` is more accurate than `CD` but can be slower + (especially if not using Numba). + + For simplicity, all methods assume constant grid spacing. + + + .. [1] G. D. Weymouth & D. K.-P. Yue, "Conservative Volume-of-Fluid method for free-surface simulations on + Cartesian-grids," Journal of Computational Physics, vol. 229, pp. 2853--2865, 2010. + [10.1016/j.jcp.2009.12.018](https://doi.org/10.1016/j.jcp.2009.12.018) + """ match normals_method: case 'CD': - return get_normals_CD(f, NORMALS_TYPE) + return get_normals_CD(void_fraction, NORMALS_TYPE) case 'WY': - return get_normals_WY(f, NORMALS_TYPE) + return get_normals_WY(void_fraction, NORMALS_TYPE) case _: raise ValueError(f"normals_method '{normals_method}' is not supported") diff --git a/blobid/solver.py b/blobid/solver.py index 213d585..a58d278 100644 --- a/blobid/solver.py +++ b/blobid/solver.py @@ -37,7 +37,8 @@ def get_labels( use_normals : bool, optional Whether to use interface normals to determine connectedness. Defaults to True. normals_method: {'CD', 'WY'}, optional - Method used to calculate interface normals. Defaults to 'CD'. + Passed to [reconstruction.get_normals()](blobid/reconstruction.html#get_normals) for calculate interface + normals. label_type : dtype, optional Determines the integer data-type of `labels`. Defaults to `numpy.uint32` @@ -79,20 +80,13 @@ def get_labels( In 3-D, each cell has 6 neighbor cells it can be connected to. The first requirement for two neighbor cells to be connected is that they are both object cells. In addition (unless `use_normals` is False), the interface normals of each cell must not be opposed.[2]_ - - The sign of the interface normal :math:`\langle n_x, n_y, n_z\rangle` using the method set by `normals_method`. - - - If `normals_method` is `CD`, central differencing is used. For example - .. math:: - \mathrm{sign}(n_x) = - \mathrm{sign}(f_{i+1} - f_{i-1}) - - - If `normals_method` is `WY`, normals are calculated based on a 3x3x3 stencil using a second-order interface - reconstruction method.[3]_ + Interface normals are calculated using the method set by `normals_method` + (see [reconstruction.get_normals()](blobid/reconstruction.html#get_normals)). Labeling -------- Based on this object and connectivity criteria, labels are assigned to each blob (defined as a connected set - of object cells) with a two-scan labeling algorithm.[4]_ + of object cells) with a two-scan labeling algorithm.[3]_ Examples ======== @@ -136,13 +130,6 @@ def get_labels( of connected components with volume-of-fluid method," Computers & Fluids, vol. 197, pp. 104373, 2020. [10.1016/j.compfluid.2019.104373](https://doi.org/10.1016/j.compfluid.2019.104373) - .. [3] G. D. Weymouth & D. K.-P. Yue, "Conservative Volume-of-Fluid method for free-surface simulations on - Cartesian-grids," Journal of Computational Physics, vol. 229, pp. 2853--2865, 2010. - [10.1016/j.jcp.2009.12.018](https://doi.org/10.1016/j.jcp.2009.12.018) - - .. [4] L. He, Y. Chao & K. Suzuki, "A Linear-Time Two-Scan Labeling Algorithm," IEEE International Conference - on Image Processing, 2007. [10.1109/ICIP.2007.4379810](https://doi.org/10.1109/ICIP.2007.4379810) - """ # Defaults From eb27bece11756f556dc973db8657295d82db6777 Mon Sep 17 00:00:00 2001 From: Declan Gaylo Date: Tue, 29 Jul 2025 13:00:35 -0400 Subject: [PATCH 4/6] Rewrite reconstruction testing --- blobid/reconstruction/__init__.py | 4 +- blobid/reconstruction/reconstruction.py | 8 ++ tests/reconstruction/test_get_normals.py | 43 +++++++++++ tests/reconstruction/test_no_inversion.py | 40 ++++++++++ tests/reconstruction/test_reconstruction.py | 83 --------------------- 5 files changed, 93 insertions(+), 85 deletions(-) create mode 100644 tests/reconstruction/test_get_normals.py create mode 100644 tests/reconstruction/test_no_inversion.py delete mode 100644 tests/reconstruction/test_reconstruction.py diff --git a/blobid/reconstruction/__init__.py b/blobid/reconstruction/__init__.py index 9710e13..458fe28 100644 --- a/blobid/reconstruction/__init__.py +++ b/blobid/reconstruction/__init__.py @@ -1,3 +1,3 @@ -from .reconstruction import get_normals +from .reconstruction import get_normals, SUPPORTED_METHODS -__all__ = ['get_normals'] +__all__ = ['get_normals', 'SUPPORTED_METHODS'] diff --git a/blobid/reconstruction/reconstruction.py b/blobid/reconstruction/reconstruction.py index 2287a71..37d8ec3 100644 --- a/blobid/reconstruction/reconstruction.py +++ b/blobid/reconstruction/reconstruction.py @@ -8,6 +8,9 @@ NORMALS_TYPE = np.int8 +SUPPORTED_METHODS = ['CD', 'WY'] +r"""Supported methods for calculating interface normals""" + def get_normals(void_fraction: np.ndarray, normals_method: str = 'CD') -> np.ndarray: r""" @@ -24,6 +27,7 @@ def get_normals(void_fraction: np.ndarray, normals_method: str = 'CD') -> np.nda ------- normals : ndarray[3, ni, nj, nk] An integer array which returns the sign of the interface normal. + If the interface in a grid cell is undefined, :math:`\langle 0, 0, 0\rangle` will be returned. Notes ----- @@ -47,6 +51,10 @@ def get_normals(void_fraction: np.ndarray, normals_method: str = 'CD') -> np.nda [10.1016/j.jcp.2009.12.018](https://doi.org/10.1016/j.jcp.2009.12.018) """ + # checks + assert void_fraction.ndim == 3 + assert all(dim > 2 for dim in void_fraction.shape) + match normals_method: case 'CD': return get_normals_CD(void_fraction, NORMALS_TYPE) diff --git a/tests/reconstruction/test_get_normals.py b/tests/reconstruction/test_get_normals.py new file mode 100644 index 0000000..f2d76c1 --- /dev/null +++ b/tests/reconstruction/test_get_normals.py @@ -0,0 +1,43 @@ +import hashlib + +import numpy as np +import blobid.reconstruction as rec + + +def test_undefined_normals(): + """Tests where normal should be all zeros (undefined)""" + + def assert_undefined(norm): + n = norm.squeeze() + assert n[0] == 0 + assert n[1] == 0 + assert n[2] == 0 + + # all zeros + f = np.zeros((3, 3, 3)) + for method in rec.SUPPORTED_METHODS: + assert_undefined(rec.get_normals(f, method)) + + # all ones + f = np.ones((3, 3, 3)) + for method in rec.SUPPORTED_METHODS: + assert_undefined(rec.get_normals(f, method)) + + # any value in a range between 0 and 1 + for _ in range(1000): + f = np.ones((3, 3, 3)) * np.random.rand() + for method in rec.SUPPORTED_METHODS: + assert_undefined(rec.get_normals(f, method)) + + +def test_end_to_end(fs_vof): + expected_result = [ + ['CD', 'c90626929711b744621692f3561700a1'], + ['WY', '9ac131e45e11e9f9457327c01ad3fb8e'] + ] + + for method, signature in expected_result: + normals = rec.get_normals(fs_vof, method) + + hash_object = hashlib.md5(normals.tobytes(order='C')) + assert hash_object.hexdigest() == signature diff --git a/tests/reconstruction/test_no_inversion.py b/tests/reconstruction/test_no_inversion.py new file mode 100644 index 0000000..0bd53c7 --- /dev/null +++ b/tests/reconstruction/test_no_inversion.py @@ -0,0 +1,40 @@ +import warnings + +import pytest +import numpy as np +import blobid.reconstruction.no_inversion as nivira + +from blobid._numba_support import numba_availible +from blobid.reconstruction.reconstruction import NORMALS_TYPE + + +def test_WY_Fig2(): + """Figure 2 from Weymouth and Yue""" + + f = np.empty((3, 3, 3)) + for k in range(3): + f[:, :, k] = np.array([[1.0, 1.0, 0.95], + [1.0, 0.9, 0.3], + [0.7, 0.15, 0.0]]) + + assert nivira._get_dominant_direction(f)[0, 0, 0] == 0 + + n = nivira.get_normals(f, NORMALS_TYPE).squeeze() + assert n[0] > 0 + assert n[1] > 0 + assert n[2] == 0 + + +def test_vof_type(): + f = np.random.rand(3, 3, 3) + + # float16 should give a warning + if numba_availible: + with pytest.warns(): + nivira.get_normals(f.astype(np.float16), NORMALS_TYPE) + + # float32 and float64 should work without warning + with warnings.catch_warnings(): + warnings.simplefilter("error") + nivira.get_normals(f.astype(np.float32), NORMALS_TYPE) + nivira.get_normals(f.astype(np.float64), NORMALS_TYPE) diff --git a/tests/reconstruction/test_reconstruction.py b/tests/reconstruction/test_reconstruction.py deleted file mode 100644 index 026a1bd..0000000 --- a/tests/reconstruction/test_reconstruction.py +++ /dev/null @@ -1,83 +0,0 @@ -import warnings - -import pytest -import numpy as np - -from blobid._numba_support import numba_availible -from blobid.reconstruction.reconstruction import get_normals, NORMALS_TYPE -from blobid.reconstruction.no_inversion import _get_dominant_direction - - -def test_reconstruction_CD(): - # all zeros - f = np.zeros((3, 3, 3)) - n = get_normals(f, 'CD').squeeze() - assert n[0] == 0 - assert n[1] == 0 - assert n[2] == 0 - - assert isinstance(n, np.ndarray) - assert isinstance(n[0], NORMALS_TYPE) - - # all ones - f = np.zeros((3, 3, 3)) - n = get_normals(f, 'CD').squeeze() - assert n[0] == 0 - assert n[1] == 0 - assert n[2] == 0 - - for _ in range(1000): - f = np.random.rand(3, 3, 3) - n = get_normals(f, 'CD').squeeze() - - assert n[0] == -np.sign(f[2, 1, 1]-f[0, 1, 1]) - assert n[1] == -np.sign(f[1, 2, 1]-f[1, 0, 1]) - assert n[2] == -np.sign(f[1, 1, 2]-f[1, 1, 0]) - - -def test_reconstruction_WY(): - # all zeros - f = np.zeros((3, 3, 3)) - n = get_normals(f, 'WY').squeeze() - assert n[0] == 0 - assert n[1] == 0 - assert n[2] == 0 - - assert isinstance(n, np.ndarray) - assert isinstance(n[0], NORMALS_TYPE) - - # all ones - f = np.zeros((3, 3, 3)) - n = get_normals(f, 'WY').squeeze() - assert n[0] == 0 - assert n[1] == 0 - assert n[2] == 0 - - # Figure 2 from Weymouth and Yue - f = np.empty((3, 3, 3)) - for k in range(3): - f[:, :, k] = np.array([[1.0, 1.0, 0.95], - [1.0, 0.9, 0.3], - [0.7, 0.15, 0.0]]) - - assert _get_dominant_direction(f)[0, 0, 0] == 0 - - n = get_normals(f, 'WY').squeeze() - assert n[0] > 0 - assert n[1] > 0 - assert n[2] == 0 - - -def test_vof_type(): - f = np.random.rand(3, 3, 3) - - # float16 should give a warning - if numba_availible: - with pytest.warns(): - get_normals(f.astype(np.float16), 'WY') - - # float32 and float64 should work without warning - with warnings.catch_warnings(): - warnings.simplefilter("error") - get_normals(f.astype(np.float32), 'WY') - get_normals(f.astype(np.float64), 'WY') From 8d009c9cb6928a285ce570a8fc4d5797e3fa15e5 Mon Sep 17 00:00:00 2001 From: Declan Gaylo Date: Tue, 29 Jul 2025 13:26:09 -0400 Subject: [PATCH 5/6] Rename files --- blobid/__init__.py | 4 ++-- blobid/plic/__init__.py | 8 ++++++++ .../central_difference.py | 0 .../no_inversion.py => plic/nivira.py} | 0 .../reconstruction.py => plic/plic.py} | 4 ++-- blobid/reconstruction/__init__.py | 3 --- blobid/solver.py | 6 +++--- .../{reconstruction => plic}/test_get_normals.py | 16 ++++++++-------- .../test_no_inversion.py => plic/test_nivira.py} | 4 ++-- tests/test_boundaries.py | 7 ++++--- 10 files changed, 29 insertions(+), 23 deletions(-) create mode 100644 blobid/plic/__init__.py rename blobid/{reconstruction => plic}/central_difference.py (100%) rename blobid/{reconstruction/no_inversion.py => plic/nivira.py} (100%) rename blobid/{reconstruction/reconstruction.py => plic/plic.py} (95%) delete mode 100644 blobid/reconstruction/__init__.py rename tests/{reconstruction => plic}/test_get_normals.py (67%) rename tests/{reconstruction/test_no_inversion.py => plic/test_nivira.py} (89%) diff --git a/blobid/__init__.py b/blobid/__init__.py index cf94de3..b858f86 100644 --- a/blobid/__init__.py +++ b/blobid/__init__.py @@ -1,5 +1,5 @@ from .solver import get_labels from . import labeling -from . import reconstruction +from . import plic -__all__ = ['get_labels', 'labeling', 'reconstruction'] +__all__ = ['get_labels', 'labeling', 'plic'] diff --git a/blobid/plic/__init__.py b/blobid/plic/__init__.py new file mode 100644 index 0000000..75c7461 --- /dev/null +++ b/blobid/plic/__init__.py @@ -0,0 +1,8 @@ +r""" +Implements piecewise linear interface calculation (PLIC), where the interface in each cell can be represented by +:math:`\vec{n}\cdot\vec{x}=\alpha` +""" + +from .plic import get_normals, SUPPORTED_METHODS + +__all__ = ['get_normals', 'SUPPORTED_METHODS'] diff --git a/blobid/reconstruction/central_difference.py b/blobid/plic/central_difference.py similarity index 100% rename from blobid/reconstruction/central_difference.py rename to blobid/plic/central_difference.py diff --git a/blobid/reconstruction/no_inversion.py b/blobid/plic/nivira.py similarity index 100% rename from blobid/reconstruction/no_inversion.py rename to blobid/plic/nivira.py diff --git a/blobid/reconstruction/reconstruction.py b/blobid/plic/plic.py similarity index 95% rename from blobid/reconstruction/reconstruction.py rename to blobid/plic/plic.py index 37d8ec3..34be0f1 100644 --- a/blobid/reconstruction/reconstruction.py +++ b/blobid/plic/plic.py @@ -4,12 +4,12 @@ import numpy as np from .central_difference import get_normals as get_normals_CD -from .no_inversion import get_normals as get_normals_WY +from .nivira import get_normals as get_normals_WY NORMALS_TYPE = np.int8 SUPPORTED_METHODS = ['CD', 'WY'] -r"""Supported methods for calculating interface normals""" +"""Supported methods for calculating interface normals""" def get_normals(void_fraction: np.ndarray, normals_method: str = 'CD') -> np.ndarray: diff --git a/blobid/reconstruction/__init__.py b/blobid/reconstruction/__init__.py deleted file mode 100644 index 458fe28..0000000 --- a/blobid/reconstruction/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ -from .reconstruction import get_normals, SUPPORTED_METHODS - -__all__ = ['get_normals', 'SUPPORTED_METHODS'] diff --git a/blobid/solver.py b/blobid/solver.py index a58d278..9ec98d8 100644 --- a/blobid/solver.py +++ b/blobid/solver.py @@ -8,7 +8,7 @@ from ._domain import VOFDomain from . import labeling -from . import reconstruction +from . import plic def get_labels( @@ -37,7 +37,7 @@ def get_labels( use_normals : bool, optional Whether to use interface normals to determine connectedness. Defaults to True. normals_method: {'CD', 'WY'}, optional - Passed to [reconstruction.get_normals()](blobid/reconstruction.html#get_normals) for calculate interface + Passed to [plic.get_normals()](blobid/plic.html#get_normals) for calculating interface normals. label_type : dtype, optional Determines the integer data-type of `labels`. Defaults to `numpy.uint32` @@ -147,7 +147,7 @@ def get_labels( # calculate connectivity if use_normals: - normals = reconstruction.get_normals(domain.vof(padding=1), normals_method) + normals = plic.get_normals(domain.vof(padding=1), normals_method) else: normals = None diff --git a/tests/reconstruction/test_get_normals.py b/tests/plic/test_get_normals.py similarity index 67% rename from tests/reconstruction/test_get_normals.py rename to tests/plic/test_get_normals.py index f2d76c1..56fa964 100644 --- a/tests/reconstruction/test_get_normals.py +++ b/tests/plic/test_get_normals.py @@ -1,7 +1,7 @@ import hashlib import numpy as np -import blobid.reconstruction as rec +import blobid.plic as plic def test_undefined_normals(): @@ -15,19 +15,19 @@ def assert_undefined(norm): # all zeros f = np.zeros((3, 3, 3)) - for method in rec.SUPPORTED_METHODS: - assert_undefined(rec.get_normals(f, method)) + for method in plic.SUPPORTED_METHODS: + assert_undefined(plic.get_normals(f, method)) # all ones f = np.ones((3, 3, 3)) - for method in rec.SUPPORTED_METHODS: - assert_undefined(rec.get_normals(f, method)) + for method in plic.SUPPORTED_METHODS: + assert_undefined(plic.get_normals(f, method)) # any value in a range between 0 and 1 for _ in range(1000): f = np.ones((3, 3, 3)) * np.random.rand() - for method in rec.SUPPORTED_METHODS: - assert_undefined(rec.get_normals(f, method)) + for method in plic.SUPPORTED_METHODS: + assert_undefined(plic.get_normals(f, method)) def test_end_to_end(fs_vof): @@ -37,7 +37,7 @@ def test_end_to_end(fs_vof): ] for method, signature in expected_result: - normals = rec.get_normals(fs_vof, method) + normals = plic.get_normals(fs_vof, method) hash_object = hashlib.md5(normals.tobytes(order='C')) assert hash_object.hexdigest() == signature diff --git a/tests/reconstruction/test_no_inversion.py b/tests/plic/test_nivira.py similarity index 89% rename from tests/reconstruction/test_no_inversion.py rename to tests/plic/test_nivira.py index 0bd53c7..8616ae0 100644 --- a/tests/reconstruction/test_no_inversion.py +++ b/tests/plic/test_nivira.py @@ -2,10 +2,10 @@ import pytest import numpy as np -import blobid.reconstruction.no_inversion as nivira +import blobid.plic.nivira as nivira from blobid._numba_support import numba_availible -from blobid.reconstruction.reconstruction import NORMALS_TYPE +from blobid.plic.plic import NORMALS_TYPE def test_WY_Fig2(): diff --git a/tests/test_boundaries.py b/tests/test_boundaries.py index 7254ada..a7ab173 100644 --- a/tests/test_boundaries.py +++ b/tests/test_boundaries.py @@ -1,10 +1,11 @@ import pytest import numpy as np +from blobid import plic + from blobid._domain import _pad_array from blobid.labeling.database import LabelDatabase from blobid.labeling.labeling import _stitch_boundaries -from blobid.reconstruction import get_normals @pytest.fixture @@ -138,8 +139,8 @@ def test_normal_calculation_at_boundaries(vof_3d): f = _pad_array(vof_3d.copy(), per, 1, 1) # normals should be the same on periodic sides - n = get_normals(f, normals_method='CD') + n = plic.get_normals(f, normals_method='CD') check_periodicity(np.moveaxis(n, 0, 3), per, 1) - n = get_normals(f, normals_method='WY') + n = plic.get_normals(f, normals_method='WY') check_periodicity(np.moveaxis(n, 0, 3), per, 1) From e56aec827b6b90aa1985f04971afef3b6bbf3db7 Mon Sep 17 00:00:00 2001 From: Declan Gaylo Date: Tue, 29 Jul 2025 13:46:22 -0400 Subject: [PATCH 6/6] Fixed numba-dependent test results --- tests/plic/test_get_normals.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/plic/test_get_normals.py b/tests/plic/test_get_normals.py index 56fa964..9b80407 100644 --- a/tests/plic/test_get_normals.py +++ b/tests/plic/test_get_normals.py @@ -37,7 +37,7 @@ def test_end_to_end(fs_vof): ] for method, signature in expected_result: - normals = plic.get_normals(fs_vof, method) + normals = plic.get_normals(fs_vof.astype(np.float32), method) hash_object = hashlib.md5(normals.tobytes(order='C')) assert hash_object.hexdigest() == signature