diff --git a/blobid/__init__.py b/blobid/__init__.py
index f13abce..b858f86 100644
--- a/blobid/__init__.py
+++ b/blobid/__init__.py
@@ -1,4 +1,5 @@
 from .solver import get_labels
 from . import labeling
+from . import plic
 
-__all__ = ['get_labels', 'labeling']
+__all__ = ['get_labels', 'labeling', 'plic']
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/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/plic/central_difference.py b/blobid/plic/central_difference.py
new file mode 100644
index 0000000..832a803
--- /dev/null
+++ b/blobid/plic/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/reconstruction.py b/blobid/plic/nivira.py
similarity index 68%
rename from blobid/reconstruction/reconstruction.py
rename to blobid/plic/nivira.py
index 46269a8..fe9f648 100644
--- a/blobid/reconstruction/reconstruction.py
+++ b/blobid/plic/nivira.py
@@ -1,42 +1,25 @@
 """
-Calculate the sign of the interface normal based on a 3x3(x3) stencil
+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
 
-NORMALS_TYPE = np.int8
-
-
-def normals(f: np.ndarray, normals_method: str) -> np.ndarray:
-    """
-    Calculate the sign of the interface normals using method specified by normals_method
-    """
-    match normals_method:
-        case 'CD':
-            return normals_CD(f)
-        case 'WY':
-            return normals_WY(f)
-        case _:
-            raise ValueError(f"normals_method '{normals_method}' is not supported")
-
 
-def normals_CD(f: np.ndarray) -> np.ndarray:
+def get_normals(f: np.ndarray, normals_dtype: npt.DTypeLike) -> np.ndarray:
     """
-    Calculate the sign of the interface normal using central difference
+    Calculate the sign of the interface normal using NIVIRA
     """
 
-    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")
@@ -46,10 +29,10 @@ def normals_WY(f: np.ndarray) -> np.ndarray:
     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)
+    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)
+        n[d] = _normals_WY(f, d_dom, d, normals_dtype)
 
     return n
 
@@ -71,8 +54,8 @@ def _get_dominant_direction(f):
 
 
 @njit
-def _normals_WY(f, d_dom, d) -> np.ndarray:
-    n = np.zeros(d_dom.shape, dtype=NORMALS_TYPE)
+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:
diff --git a/blobid/plic/plic.py b/blobid/plic/plic.py
new file mode 100644
index 0000000..34be0f1
--- /dev/null
+++ b/blobid/plic/plic.py
@@ -0,0 +1,64 @@
+"""
+Calculate the sign of the interface normal based on a 3x3(x3) stencil
+"""
+import numpy as np
+
+from .central_difference import get_normals as get_normals_CD
+from .nivira import get_normals as get_normals_WY
+
+NORMALS_TYPE = np.int8
+
+SUPPORTED_METHODS = ['CD', 'WY']
+"""Supported methods for calculating interface normals"""
+
+
+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.
+            If the interface in a grid cell is undefined, :math:`\langle 0, 0, 0\rangle` will be returned.
+
+    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)
+
+    """
+    # 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)
+        case 'WY':
+            return get_normals_WY(void_fraction, NORMALS_TYPE)
+        case _:
+            raise ValueError(f"normals_method '{normals_method}' is not supported")
diff --git a/blobid/reconstruction/__init__.py b/blobid/reconstruction/__init__.py
deleted file mode 100644
index bd984a5..0000000
--- a/blobid/reconstruction/__init__.py
+++ /dev/null
@@ -1,3 +0,0 @@
-from .reconstruction import normals
-
-__all__ = ['normals']
diff --git a/blobid/solver.py b/blobid/solver.py
index 662e436..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,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 [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`
 
@@ -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
@@ -160,7 +147,7 @@ def get_labels(
 
     # calculate connectivity
     if use_normals:
-        normals = reconstruction.normals(domain.vof(padding=1), normals_method)
+        normals = plic.get_normals(domain.vof(padding=1), normals_method)
     else:
         normals = None
 
diff --git a/tests/plic/test_get_normals.py b/tests/plic/test_get_normals.py
new file mode 100644
index 0000000..9b80407
--- /dev/null
+++ b/tests/plic/test_get_normals.py
@@ -0,0 +1,43 @@
+import hashlib
+
+import numpy as np
+import blobid.plic as plic
+
+
+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 plic.SUPPORTED_METHODS:
+        assert_undefined(plic.get_normals(f, method))
+
+    # all ones
+    f = np.ones((3, 3, 3))
+    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 plic.SUPPORTED_METHODS:
+            assert_undefined(plic.get_normals(f, method))
+
+
+def test_end_to_end(fs_vof):
+    expected_result = [
+        ['CD', 'c90626929711b744621692f3561700a1'],
+        ['WY', '9ac131e45e11e9f9457327c01ad3fb8e']
+    ]
+
+    for method, signature in expected_result:
+        normals = plic.get_normals(fs_vof.astype(np.float32), method)
+
+        hash_object = hashlib.md5(normals.tobytes(order='C'))
+        assert hash_object.hexdigest() == signature
diff --git a/tests/plic/test_nivira.py b/tests/plic/test_nivira.py
new file mode 100644
index 0000000..8616ae0
--- /dev/null
+++ b/tests/plic/test_nivira.py
@@ -0,0 +1,40 @@
+import warnings
+
+import pytest
+import numpy as np
+import blobid.plic.nivira as nivira
+
+from blobid._numba_support import numba_availible
+from blobid.plic.plic 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 d1f03fc..0000000
--- a/tests/reconstruction/test_reconstruction.py
+++ /dev/null
@@ -1,82 +0,0 @@
-import warnings
-
-import pytest
-import numpy as np
-
-from blobid._numba_support import numba_availible
-from blobid.reconstruction.reconstruction import normals, NORMALS_TYPE, _get_dominant_direction
-
-
-def test_reconstruction_CD():
-    # all zeros
-    f = np.zeros((3, 3, 3))
-    n = 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 = 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()
-
-        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 = 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 = 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 = 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():
-            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')
diff --git a/tests/test_boundaries.py b/tests/test_boundaries.py
index 3f4c47a..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 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 = normals(f, normals_method='CD')
+        n = plic.get_normals(f, normals_method='CD')
         check_periodicity(np.moveaxis(n, 0, 3), per, 1)
 
-        n = normals(f, normals_method='WY')
+        n = plic.get_normals(f, normals_method='WY')
         check_periodicity(np.moveaxis(n, 0, 3), per, 1)