From 3b44c7ec7069f0953c3356304112cadbe55d36f1 Mon Sep 17 00:00:00 2001
From: daubners <s.daubner@imperial.ac.uk>
Date: Wed, 15 Jan 2025 16:57:26 +0100
Subject: [PATCH 01/12] refactor volume fraction computation

---
 taufactor/metrics.py | 18 ++++++++++--------
 1 file changed, 10 insertions(+), 8 deletions(-)

diff --git a/taufactor/metrics.py b/taufactor/metrics.py
index b84667f..2a5735a 100644
--- a/taufactor/metrics.py
+++ b/taufactor/metrics.py
@@ -2,7 +2,7 @@
 import torch
 import torch.nn.functional as F
 
-def volume_fraction(img, phases={}):
+def volume_fraction(img, phases={}, device=torch.device('cuda')):
     """
     Calculates volume fractions of phases in an image
     :param img: segmented input image with n phases
@@ -11,15 +11,17 @@ def volume_fraction(img, phases={}):
     """
 
     if type(img) is not type(torch.tensor(1)):
-        img = torch.tensor(img)
+        img = torch.tensor(img, device=device)
 
     if phases=={}:
-        phases = torch.unique(img)
-        vf_out = []
-        for p in phases:
-            vf_out.append((img==p).to(torch.float).mean().item())
-        if len(vf_out)==1:
-            vf_out=vf_out[0]
+        volume = torch.numel(img)
+        labels, counts = torch.unique(img, return_counts=True)
+        labels = labels.int()
+        counts = counts.float()
+        counts /= volume
+        vf_out = {}
+        for i, label in enumerate(labels):
+            vf_out[str(label.item())] = counts[i].item()
     else:
         vf_out={}
         for p in phases:

From 1469ab7bc79c18596686e76e1ef059d056e1d01f Mon Sep 17 00:00:00 2001
From: daubners <s.daubner@imperial.ac.uk>
Date: Wed, 15 Jan 2025 17:29:21 +0100
Subject: [PATCH 02/12] add through feature computation

---
 taufactor/metrics.py | 121 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 121 insertions(+)

diff --git a/taufactor/metrics.py b/taufactor/metrics.py
index 2a5735a..58e867a 100644
--- a/taufactor/metrics.py
+++ b/taufactor/metrics.py
@@ -2,6 +2,8 @@
 import torch
 import torch.nn.functional as F
 
+from scipy.ndimage import label, generate_binary_structure
+
 def volume_fraction(img, phases={}, device=torch.device('cuda')):
     """
     Calculates volume fractions of phases in an image
@@ -144,3 +146,122 @@ def triple_phase_boundary(img):
             tpb += torch.sum(tpb_map)
 
     return tpb/total_edges
+
+def label_periodic(field, grayscale_value, neighbour_structure, periodic, debug=False):
+    # Initialize phi field whith enlarged dimensions in periodic directions. Boundary values of
+    # array are copied into ghost cells which are necessary to impose boundary conditions.
+    padx = int(periodic[0])
+    pady = int(periodic[1])
+    padz = int(periodic[2])
+    mask = np.pad(field, ((padx, padx), (pady, pady), (padz, padz)), mode='wrap')
+    labeled_mask, num_labels = label(mask==grayscale_value, structure=neighbour_structure)
+    count = 1
+    for k in range(100):
+        # Find indices where labels are different at the boundaries and create swaplist
+        swap_list = np.zeros((1,2))
+        if periodic[0]:
+            # right x
+            indices = np.where((labeled_mask[0,:,:]!=labeled_mask[-2,:,:]) & (labeled_mask[0,:,:]!=0) & (labeled_mask[-2,:,:]!=0))
+            additional_swaps = np.column_stack((labeled_mask[0,:,:][indices], labeled_mask[-2,:,:][indices]))
+            swap_list = np.row_stack((swap_list,additional_swaps))
+            # left x
+            indices = np.where((labeled_mask[1,:,:]!=labeled_mask[-1,:,:]) & (labeled_mask[1,:,:]!=0) & (labeled_mask[-1,:,:]!=0))
+            additional_swaps = np.column_stack((labeled_mask[1,:,:][indices], labeled_mask[-1,:,:][indices]))
+            swap_list = np.row_stack((swap_list,additional_swaps))
+        if periodic[1]:
+            # top y
+            indices = np.where((labeled_mask[:,0,:]!=labeled_mask[:,-2,:]) & (labeled_mask[:,0,:]!=0) & (labeled_mask[:,-2,:]!=0))
+            additional_swaps = np.column_stack((labeled_mask[:,0,:][indices], labeled_mask[:,-2,:][indices]))
+            swap_list = np.row_stack((swap_list,additional_swaps))
+            # bottom y
+            indices = np.where((labeled_mask[:,1,:]!=labeled_mask[:,-1,:]) & (labeled_mask[:,1,:]!=0) & (labeled_mask[:,-1,:]!=0))
+            additional_swaps = np.column_stack((labeled_mask[:,1,:][indices], labeled_mask[:,-1,:][indices]))
+            swap_list = np.row_stack((swap_list,additional_swaps))
+        if periodic[2]:
+            # front z
+            indices = np.where((labeled_mask[:,:,0]!=labeled_mask[:,:,-2]) & (labeled_mask[:,:,0]!=0) & (labeled_mask[:,:,-2]!=0))
+            additional_swaps = np.column_stack((labeled_mask[:,:,0][indices], labeled_mask[:,:,-2][indices]))
+            swap_list = np.row_stack((swap_list,additional_swaps))
+            # back z
+            indices = np.where((labeled_mask[:,:,1]!=labeled_mask[:,:,-1]) & (labeled_mask[:,:,1]!=0) & (labeled_mask[:,:,-1]!=0))
+            additional_swaps = np.column_stack((labeled_mask[:,:,1][indices], labeled_mask[:,:,-1][indices]))
+            swap_list = np.row_stack((swap_list,additional_swaps))
+        swap_list = swap_list[1:,:]
+        # Sort swap list columns to ensure consistent ordering
+        swap_list = np.sort(swap_list, axis=1)
+
+        # Remove duplicates from swap_list
+        swap_list = np.unique(swap_list, axis=0)
+        # print(f"swap_list contains {swap_list.shape[0]} elements.")
+        if (swap_list.shape[0]==0):
+            break
+        for i in range(swap_list.shape[0]):
+            index = swap_list.shape[0] - i -1
+            labeled_mask[labeled_mask == swap_list[index][1]] = swap_list[index][0]
+        count += 1
+    if(debug):
+        print(f"Did {count} iterations for periodic labelling.")
+    dim = labeled_mask.shape
+    return labeled_mask[padx:dim[0]-padx,pady:dim[1]-pady,padz:dim[2]-padz], np.unique(labeled_mask).size-1
+
+def find_spanning_labels(labelled_array, axis):
+    """
+    Find labels that appear on both ends along given axis
+
+    Returns:
+        set: Labels that appear on both ends of the first axis.
+    """
+    if axis == "x":
+        front = np.s_[0,:,:]
+        end   = np.s_[-1,:,:]
+    elif axis == "y":
+        front = np.s_[:,0,:]
+        end   = np.s_[:,-1,:]
+    elif axis == "z":
+        front = np.s_[:,:,0]
+        end   = np.s_[:,:,-1]
+    else:
+        raise ValueError("Axis should be x, y or z!")
+
+    first_slice_labels = np.unique(labelled_array[front])
+    last_slice_labels = np.unique(labelled_array[end])
+    spanning_labels = set(first_slice_labels) & set(last_slice_labels)
+    spanning_labels.discard(0)  # Remove the background label if it exists
+    return spanning_labels
+
+def extract_through_feature(array, grayscale_value, axis, periodic=[False,False,False], connectivity=1, debug=False):
+    if array.ndim != 3:
+        print(f"Expected a 3D array, but got an array with {array.ndim} dimension(s).")
+        return None
+
+    # Compute volume fraction of given grayscale value
+    vol_phase = volume_fraction(array, phases={'1': grayscale_value})['1']
+
+    # Define a list of connectivities to loop over
+    connectivities_to_loop_over = [connectivity] if connectivity else range(1, 4)
+    through_feature = []
+    through_feature_fraction = np.zeros(len(connectivities_to_loop_over))
+
+    # Compute the largest interconnected features depending on given connectivity
+    count = 0
+    for conn in connectivities_to_loop_over:
+        # connectivity 1 = cells connected by sides (6 neighbours)
+        # connectivity 2 = cells connected by sides & edges (14 neighbours)
+        # connectivity 3 = cells connected by sides & edges & corners (26 neighbours)
+        neighbour_structure = generate_binary_structure(3,conn)
+        # Label connected components in the mask with given neighbour structure
+        if any(periodic):
+            labeled_mask, num_labels = label_periodic(array, grayscale_value, neighbour_structure, periodic, debug=debug)
+        else:
+            labeled_mask, num_labels = label(array == grayscale_value, structure=neighbour_structure)
+        if(debug):
+            print(f"Found {num_labels} labelled regions. For connectivity {conn} and grayscale {grayscale_value}.")
+
+        through_labels = find_spanning_labels(labeled_mask,axis)
+        spanning_network = np.isin(labeled_mask, list(through_labels))
+
+        through_feature.append(spanning_network)
+        through_feature_fraction[count] = volume_fraction(spanning_network, phases={'1': 1})['1']/vol_phase
+        count += 1
+
+    return through_feature, through_feature_fraction
\ No newline at end of file

From c53e81846f8b9914028cc5f512d500e34291ec05 Mon Sep 17 00:00:00 2001
From: daubners <s.daubner@imperial.ac.uk>
Date: Wed, 15 Jan 2025 17:30:09 +0100
Subject: [PATCH 03/12] check inf for blocked pores

---
 taufactor/taufactor.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/taufactor/taufactor.py b/taufactor/taufactor.py
index a11322f..85fc192 100644
--- a/taufactor/taufactor.py
+++ b/taufactor/taufactor.py
@@ -7,6 +7,7 @@
 except ImportError:
     raise ImportError("Pytorch is required to use this package. Please install pytorch and try again. More information about TauFactor's requirements can be found at https://taufactor.readthedocs.io/en/latest/")
 import warnings
+from .metrics import extract_through_feature
 
 class BaseSolver:
     def __init__(self, img, bc=(-0.5, 0.5), device=torch.device('cuda')):
@@ -89,7 +90,9 @@ def check_vertical_flux(self, conv_crit):
         fl = torch.sum(vert_flux, (0, 2, 3))
         err = (fl.max() - fl.min())/(fl.max())
         if fl.min() == 0:
-            return 'zero_flux', torch.mean(fl), err
+            _ , frac = extract_through_feature(self.cpu_img[0], 1, 'x')
+            if frac == 0:
+                return 'zero_flux', torch.mean(fl), err
         if err < conv_crit or torch.isnan(err).item():
             return True, torch.mean(fl), err
         return False, torch.mean(fl), err

From a675fe2be0a03f4ff874c65a51ce7d8e3e139d3f Mon Sep 17 00:00:00 2001
From: daubners <s.daubner@imperial.ac.uk>
Date: Wed, 15 Jan 2025 17:45:54 +0100
Subject: [PATCH 04/12] update requirements

---
 requirements_dev.txt | 3 ++-
 setup.py             | 2 +-
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/requirements_dev.txt b/requirements_dev.txt
index 8c71d95..b0a47ee 100644
--- a/requirements_dev.txt
+++ b/requirements_dev.txt
@@ -13,4 +13,5 @@ matplotlib>3.6
 pytest-runner==5.2
 numpy>=1.24.2
 tifffile==2023.2.3
-myst-parser==0.18.1
\ No newline at end of file
+myst-parser==0.18.1
+scikit-image>=0.20.0
\ No newline at end of file
diff --git a/setup.py b/setup.py
index 24c6814..3a78823 100644
--- a/setup.py
+++ b/setup.py
@@ -10,7 +10,7 @@
 with open('HISTORY.md') as history_file:
     history = history_file.read()
 
-requirements = ['Click>=7.0', 'numpy>=1.0', 'matplotlib>=3.4', 'tifffile>=2023.2.3']
+requirements = ['Click>=7.0', 'numpy>=1.0', 'matplotlib>=3.4', 'tifffile>=2023.2.3', 'scikit-image>=0.20.0']
 
 setup_requirements = ['pytest-runner', ]
 

From 1a522bd39143cfc3ad024c16e8b1efbd5af304ae Mon Sep 17 00:00:00 2001
From: daubners <s.daubner@imperial.ac.uk>
Date: Wed, 15 Jan 2025 17:53:26 +0100
Subject: [PATCH 05/12] add scipy requirements

---
 requirements_dev.txt | 3 ++-
 setup.py             | 9 ++++++++-
 2 files changed, 10 insertions(+), 2 deletions(-)

diff --git a/requirements_dev.txt b/requirements_dev.txt
index b0a47ee..54ed682 100644
--- a/requirements_dev.txt
+++ b/requirements_dev.txt
@@ -14,4 +14,5 @@ pytest-runner==5.2
 numpy>=1.24.2
 tifffile==2023.2.3
 myst-parser==0.18.1
-scikit-image>=0.20.0
\ No newline at end of file
+scikit-image>=0.20.0
+scipy>=1.3
\ No newline at end of file
diff --git a/setup.py b/setup.py
index 3a78823..cc136eb 100644
--- a/setup.py
+++ b/setup.py
@@ -10,7 +10,14 @@
 with open('HISTORY.md') as history_file:
     history = history_file.read()
 
-requirements = ['Click>=7.0', 'numpy>=1.0', 'matplotlib>=3.4', 'tifffile>=2023.2.3', 'scikit-image>=0.20.0']
+requirements = [
+    'Click>=7.0',
+    'numpy>=1.0',
+    'matplotlib>=3.4',
+    'tifffile>=2023.2.3',
+    'scikit-image>=0.20.0',
+    'scipy>=1.3'
+]
 
 setup_requirements = ['pytest-runner', ]
 

From 2040e71d112cf6f1d4ed6e35c51e4bcf19b21c2f Mon Sep 17 00:00:00 2001
From: daubners <s.daubner@imperial.ac.uk>
Date: Thu, 16 Jan 2025 11:51:45 +0100
Subject: [PATCH 06/12] update tests for volume fraction

---
 tests/test_taufactor.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/tests/test_taufactor.py b/tests/test_taufactor.py
index b36c4d5..dae7ee2 100644
--- a/tests/test_taufactor.py
+++ b/tests/test_taufactor.py
@@ -125,7 +125,7 @@ def test_volume_fraction_on_uniform_block():
     """Run volume fraction on uniform block"""
     l = 20
     img = np.ones([l, l, l]).reshape(1, l, l, l)
-    vf = volume_fraction(img)
+    vf = volume_fraction(img)['1']
 
     assert np.around(vf, decimals=5) == 1.0
 
@@ -134,7 +134,7 @@ def test_volume_fraction_on_empty_block():
     """Run volume fraction on empty block"""
     l = 20
     img = np.zeros([l, l, l]).reshape(1, l, l, l)
-    vf = volume_fraction(img)
+    vf = volume_fraction(img)['0']
 
     assert np.around(vf, decimals=5) == 1.0
 
@@ -143,9 +143,9 @@ def test_volume_fraction_on_checkerboard():
     """Run volume fraction on checkerboard block"""
     l = 20
     img = generate_checkerboard(l)
-    vf = volume_fraction(img)
+    vf = volume_fraction(img, phases={'zeros': 0, 'ones': 1})
 
-    assert vf == [0.5, 0.5]
+    assert (vf['zeros'], vf['ones']) == [0.5, 0.5]
 
 
 def test_volume_fraction_on_strip_of_ones():

From 278c562695ff9b666a8e5870d3eb25696fb5af41 Mon Sep 17 00:00:00 2001
From: daubners <s.daubner@imperial.ac.uk>
Date: Thu, 16 Jan 2025 12:56:56 +0100
Subject: [PATCH 07/12] remove cuda from volume

---
 taufactor/metrics.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/taufactor/metrics.py b/taufactor/metrics.py
index 58e867a..2176cbe 100644
--- a/taufactor/metrics.py
+++ b/taufactor/metrics.py
@@ -4,7 +4,7 @@
 
 from scipy.ndimage import label, generate_binary_structure
 
-def volume_fraction(img, phases={}, device=torch.device('cuda')):
+def volume_fraction(img, phases={}):
     """
     Calculates volume fractions of phases in an image
     :param img: segmented input image with n phases
@@ -13,7 +13,7 @@ def volume_fraction(img, phases={}, device=torch.device('cuda')):
     """
 
     if type(img) is not type(torch.tensor(1)):
-        img = torch.tensor(img, device=device)
+        img = torch.tensor(img)
 
     if phases=={}:
         volume = torch.numel(img)

From b9da23c934984e64cec0e58ec4d05aee61d329c8 Mon Sep 17 00:00:00 2001
From: daubners <s.daubner@imperial.ac.uk>
Date: Thu, 16 Jan 2025 13:11:24 +0100
Subject: [PATCH 08/12] fix tests

---
 taufactor/metrics.py    | 2 ++
 tests/test_taufactor.py | 2 +-
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/taufactor/metrics.py b/taufactor/metrics.py
index 2176cbe..7e37706 100644
--- a/taufactor/metrics.py
+++ b/taufactor/metrics.py
@@ -236,6 +236,8 @@ def extract_through_feature(array, grayscale_value, axis, periodic=[False,False,
 
     # Compute volume fraction of given grayscale value
     vol_phase = volume_fraction(array, phases={'1': grayscale_value})['1']
+    if vol_phase == 0:
+        return 0, 0
 
     # Define a list of connectivities to loop over
     connectivities_to_loop_over = [connectivity] if connectivity else range(1, 4)
diff --git a/tests/test_taufactor.py b/tests/test_taufactor.py
index dae7ee2..46dade0 100644
--- a/tests/test_taufactor.py
+++ b/tests/test_taufactor.py
@@ -145,7 +145,7 @@ def test_volume_fraction_on_checkerboard():
     img = generate_checkerboard(l)
     vf = volume_fraction(img, phases={'zeros': 0, 'ones': 1})
 
-    assert (vf['zeros'], vf['ones']) == [0.5, 0.5]
+    assert (vf['zeros'], vf['ones']) == (0.5, 0.5)
 
 
 def test_volume_fraction_on_strip_of_ones():

From 96b6e94d5747ec99dfe5664ff8a6ee4db01f0852 Mon Sep 17 00:00:00 2001
From: daubners <s.daubner@imperial.ac.uk>
Date: Tue, 21 Jan 2025 20:46:51 +0000
Subject: [PATCH 09/12] add new surface area computation plus testing

---
 taufactor/metrics.py    | 208 +++++++++++++++++++++++++++++-----------
 tests/test_metrics.py   | 179 ++++++++++++++++++++++++++++++++++
 tests/test_taufactor.py | 153 ++---------------------------
 3 files changed, 339 insertions(+), 201 deletions(-)
 create mode 100644 tests/test_metrics.py

diff --git a/taufactor/metrics.py b/taufactor/metrics.py
index 7e37706..b4eca0c 100644
--- a/taufactor/metrics.py
+++ b/taufactor/metrics.py
@@ -1,8 +1,10 @@
 import numpy as np
 import torch
 import torch.nn.functional as F
+import warnings
 
-from scipy.ndimage import label, generate_binary_structure
+from scipy.ndimage import label, generate_binary_structure, convolve
+from skimage import measure
 
 def volume_fraction(img, phases={}):
     """
@@ -31,69 +33,165 @@ def volume_fraction(img, phases={}):
 
     return vf_out
 
-def surface_area(img, phases, periodic=False):
+def crop_area_of_interest_torch(tensor, labels):
+    indices = torch.nonzero(torch.isin(tensor, labels), as_tuple=True)
+    min_idx = [torch.min(idx).item() for idx in indices]
+    max_idx = [torch.max(idx).item() for idx in indices]
+
+    # Slice the tensor to the bounding box
+    # Make sure to stay inside the bounds of total array
+    sub_tensor = tensor[max(min_idx[0]-3,0):min(max_idx[0]+4,tensor.shape[0]),
+                        max(min_idx[1]-3,0):min(max_idx[1]+4,tensor.shape[1]),
+                        max(min_idx[2]-3,0):min(max_idx[2]+4,tensor.shape[2])]
+    return sub_tensor
+
+def crop_area_of_interest_numpy(array, labels):
+    indices = np.nonzero(np.isin(array, labels))
+    min_idx = [np.min(idx) for idx in indices]
+    max_idx = [np.max(idx) for idx in indices]
+
+    # Slice the array to the bounding box
+    # Make sure to stay inside the bounds of the total array
+    sub_array = array[max(min_idx[0]-3, 0):min(max_idx[0]+4, array.shape[0]),
+                      max(min_idx[1]-3, 0):min(max_idx[1]+4, array.shape[1]),
+                      max(min_idx[2]-3, 0):min(max_idx[2]+4, array.shape[2])]
+    return sub_array
+
+def gaussian_kernel_3d_torch(size=3, sigma=1.0, device=torch.device('cuda')):
+    """Creates a 3D Gaussian kernel using PyTorch"""
+    ax = torch.linspace(-(size // 2), size // 2, size)
+    xx, yy, zz = torch.meshgrid(ax, ax, ax, indexing="ij")
+
+    # Calculate Gaussian function for each point in the grid
+    kernel = torch.exp(-(xx**2 + yy**2 + zz**2) / (2 * sigma**2))
+    kernel /= kernel.sum()
+    kernel = kernel.to(device)
+    return kernel.unsqueeze(0).unsqueeze(0)
+
+def gaussian_kernel_3d_numpy(size=3, sigma=1.0):
+    """Creates a 3D Gaussian kernel using NumPy"""
+    ax = np.linspace(-(size // 2), size // 2, size)
+    xx, yy, zz = np.meshgrid(ax, ax, ax)
+
+    # Calculate Gaussian function for each point in the grid
+    kernel = np.exp(-(xx**2 + yy**2 + zz**2) / (2 * sigma**2))
+    kernel /= np.sum(kernel)
+    return kernel
+
+def specific_surface_area(img, spacing=(1,1,1), phases={}, method='gradient', device=torch.device('cuda'), smoothing=True):
     """
-    Calculate interfacial surface area between two phases or the total surface area of one phase
-    :param img:
-    :param phases: list of phases to calculate SA, if lenght 1 calculate total SA, if length 2 calculate inerfacial SA
-    :param periodic: list of bools indicating if the image is periodic in each dimension
-    :return: the surface area in faces per unit volume
+    Calculate the specific surface area of all (specified) phases
+    :param img: labelled microstructure where each integer value represents a phase
+    :param spacing: voxel size in each dimension [dx,dy,dz]
+    :param phases: dictionary of phases {'name': label, ...}. If empty do all by default.
+    :param method: string to indicate preferred method (face_counting, marching_cubes or gradient)
+    :return: the surface area per unit volume
     """
-    shape = img.shape
-    int_not_in_img = int(np.unique(img).max()+1)
+    [dx,dy,dz] = spacing
+    surface_areas = {}
 
-    dim = len(shape)
-    img = torch.tensor(img)
-    # finding an int that is not in the img for padding:
-    
-    if periodic:
-        periodic.reverse()
-        pad = ()
-        for x in periodic:
-            pad += tuple((int(not x),)*2)
-        img = F.pad(img, pad, 'constant', value=int_not_in_img)
-        periodic.reverse()
-    else:
-        img = F.pad(img, (1,)*dim*2, 'constant', value=int_not_in_img)
-        periodic=[0]*dim
-
-    SA_map = torch.zeros_like(img)
-    if not isinstance(phases, list):
-        raise TypeError('phases should be a list')
-    for i in range(dim):
-        for j in [1, -1]:
-            i_rolled = torch.roll(img, j, i)
-            if len(phases)==2:
-                SA_map[(i_rolled == phases[0]) & (img == phases[1])] += 1
+    if (method == 'gradient') | (method == 'face_counting'):
+        if type(img) is not type(torch.tensor(1)):
+            tensor = torch.tensor(img)
+        else:
+            tensor = img
+        tensor = tensor.to(device)
+
+    if method == 'gradient':
+        if phases=={}:
+            labels = torch.unique(tensor)
+            labels = labels.int()
+            phases = {str(label.item()): label.item() for label in labels}
+        gaussian = gaussian_kernel_3d_torch(device=device)
+
+        volume = torch.numel(tensor)
+        for name, label in phases.items():
+            sub_tensor = crop_area_of_interest_torch(tensor, label)
+            # Create binary mask for the label within the slice
+            mask = (sub_tensor == label).float()
+            if smoothing:
+                mask = mask.unsqueeze(0).unsqueeze(0)
+                mask = F.pad(mask, (1,1,1,1,1,1), mode='reflect')
+                mask = F.conv3d(mask, gaussian, padding='valid')
+                mask = mask.squeeze()
+
+            grad = torch.gradient(mask, spacing=(dx,dy,dz))
+            norm2 = grad[0].pow(2) + grad[1].pow(2) + grad[2].pow(2)
+            surface_area = torch.sum(torch.sqrt(norm2)).item()
+
+            surface_areas[name] = surface_area / volume
+
+    elif method == 'face_counting':
+        # TODO: treat dimensions such that dx!=dz is accounted for
+        tensor = tensor.to(torch.int32)
+        if len(torch.unique(tensor)) == 1:
+            surface_areas = {str(tensor[0,0,0].int().item()): 0.0}
+        else:
+            volume = torch.numel(tensor)*dx
+            phasepairs = torch.tensor([[0,0]], device=device)
+
+            neighbour_idx = torch.nonzero(tensor[:-1,:,:] != tensor[1:,:,:], as_tuple=True)
+            neighbour_list = torch.stack([tensor[:-1,:,:][neighbour_idx], tensor[1:,:,:][neighbour_idx]])
+            phasepairs = torch.cat((phasepairs,torch.transpose(neighbour_list,0,1)), 0)
+            neighbour_idx = torch.nonzero(tensor[:,:-1,:] != tensor[:,1:,:], as_tuple=True)
+            neighbour_list = torch.stack([tensor[:,:-1,:][neighbour_idx], tensor[:,1:,:][neighbour_idx]])
+            phasepairs = torch.cat((phasepairs,torch.transpose(neighbour_list,0,1)), 0)
+            neighbour_idx = torch.nonzero(tensor[:,:,:-1] != tensor[:,:,1:], as_tuple=True)
+            neighbour_list = torch.stack([tensor[:,:,:-1][neighbour_idx], tensor[:,:,1:][neighbour_idx]])
+            phasepairs = torch.cat((phasepairs,torch.transpose(neighbour_list,0,1)), 0)
+
+            # Crop initial dummy values
+            phasepairs = phasepairs[1:]
+
+            if phases=={}:
+                if phasepairs == {}:
+                    surface_areas
+                labels, counts = torch.unique(phasepairs, return_counts=True)
+                labels = labels.int()
+                counts = counts.float()
+                counts /= volume
+                surface_areas = {str(label.item()): counts[i].item() for i, label in enumerate(labels)}
             else:
-                SA_map[(i_rolled == phases[0]) & (img != phases[0])] += 1
-    # remove padding
-    if not periodic[0]:
-        SA_map = SA_map[1:-1, :]
-    if not periodic[1]:
-        SA_map = SA_map[:, 1:-1]
-    x, y = shape[0], shape[1]
-    # scale factor calculated by taking into account edges
-    periodic_mask=[not x for x in periodic]
-    if dim == 3:
-        z = shape[2]
-        if not periodic[2]:
-            SA_map = SA_map[:, :, 1:-1]
-        sf = torch.sum(torch.tensor([x,y,z])[periodic_mask]*torch.roll(torch.tensor([x,y,z])[periodic_mask],1))
-        total_faces = 3*(x*y*z)-sf
-    elif dim == 2:
-        sf = torch.sum(torch.tensor([x,y])[periodic_mask])
-        total_faces = 2*(x+1)*(y+1)-(x+1)-(y+1)-2*sf
+                for name, label in phases.items():
+                    count = torch.sum((phasepairs == label).int()).item()
+                    surface_areas[name] = count / volume
+
+    elif method == 'marching_cubes':
+        if device != 'cpu':
+            warnings.warn("The marching cubes algorithm is performed on the CPU based on scikit-image package.")
+        if dx != dy | dx!= dz | dy!=dz:
+            raise ValueError("Surface area computation based on marching cubes assumes dx=dy=dz.")
+
+        if type(img) is type(torch.tensor(1)):
+            array = np.array(img.cpu())
+        else:
+            array = img
+
+        if phases=={}:
+            labels = np.unique(array).astype(int)
+            phases = {str(label): label for label in labels}
+
+        volume = array.size*dx
+        gaussian = gaussian_kernel_3d_numpy(size=3, sigma=1.0)
+        for name, label in phases.items():
+            sub_array = crop_area_of_interest_numpy(array, label)
+            sub_array = (sub_array == label).astype(float)
+            if smoothing:
+                sub_array = convolve(sub_array, gaussian, mode='nearest')
+            vertices, faces, _, _ = measure.marching_cubes(sub_array, 0.5, method='lewiner')
+            surface_area = measure.mesh_surface_area(vertices, faces)
+            surface_areas[name] = surface_area/volume
+
     else:
-        total_faces=SA_map.size
-    sa = torch.sum(SA_map)/total_faces
-    return sa
+        raise ValueError("Choose method\n 'gradient' for fast phase-field approach\n 'face_counting' for face counting or\n 'marching_cubes' for marching cubes method.")
+
+    return surface_areas
 
 def triple_phase_boundary(img):
     """Calculate triple phase boundary density i.e. fraction of voxel verticies that touch at least 3 phases
 
     Args:
-        img (numpy array): image to calculate metric on     
+        img (numpy array): image to calculate metric on
     Returns:
         float: triple phase boundary density 
     """
@@ -266,4 +364,4 @@ def extract_through_feature(array, grayscale_value, axis, periodic=[False,False,
         through_feature_fraction[count] = volume_fraction(spanning_network, phases={'1': 1})['1']/vol_phase
         count += 1
 
-    return through_feature, through_feature_fraction
\ No newline at end of file
+    return through_feature, through_feature_fraction
diff --git a/tests/test_metrics.py b/tests/test_metrics.py
new file mode 100644
index 0000000..e2072e2
--- /dev/null
+++ b/tests/test_metrics.py
@@ -0,0 +1,179 @@
+#!/usr/bin/env python
+
+"""Tests for `taufactor` package."""
+
+from taufactor.metrics import volume_fraction, specific_surface_area, triple_phase_boundary
+from tests.utils import *
+import numpy as np
+
+
+# Volume fraction
+
+def test_volume_fraction_on_uniform_block():
+    """Run volume fraction on uniform block"""
+    l = 20
+    img = np.ones([l, l, l]).reshape(1, l, l, l)
+    vf = volume_fraction(img)['1']
+
+    assert np.around(vf, decimals=5) == 1.0
+
+
+def test_volume_fraction_on_empty_block():
+    """Run volume fraction on empty block"""
+    l = 20
+    img = np.zeros([l, l, l]).reshape(1, l, l, l)
+    vf = volume_fraction(img)['0']
+
+    assert np.around(vf, decimals=5) == 1.0
+
+
+def test_volume_fraction_on_checkerboard():
+    """Run volume fraction on checkerboard block"""
+    l = 20
+    img = generate_checkerboard(l)
+    vf = volume_fraction(img, phases={'zeros': 0, 'ones': 1})
+
+    assert (vf['zeros'], vf['ones']) == (0.5, 0.5)
+
+
+def test_volume_fraction_on_strip_of_ones():
+    """Run volume fraction on strip of ones"""
+    l = 20
+    img = np.zeros([l, l, l])
+    t = 10
+    img[:, 0:t, 0:t] = 1
+    vf = volume_fraction(img, phases={'zeros': 0, 'ones': 1})
+
+    assert (vf['zeros'], vf['ones']) == (0.75, 0.25)
+
+
+def test_volume_fraction_on_multi_cubes():
+    """Run surface area on multiple cubes"""
+    l = 20
+    img = np.zeros([l, l, l])
+    img[0:10, 0:10, 0:5] = 1
+    img[5:-5, 5:-5, 5:-5] = 2
+    img[0:10, 0:10, 15:] = 1
+    img[10:, 10:, 15:] = 3
+    vf = volume_fraction(img)
+    sum = vf['0'] + vf['1'] + vf['2'] + vf['3']
+
+    assert (vf['1'], vf['2'], vf['3'], sum) == (0.125, 0.125, 0.0625, 1.0)
+
+# Surface area
+
+def test_surface_area_on_uniform_block():
+    """Run surface area on uniform block"""
+    l = 20
+    img = np.ones([l, l, l])
+    sa_f = specific_surface_area(img, method='face_counting')['1']
+    sa_g = specific_surface_area(img, method='gradient')['1']
+    # Marchin cubes will not work unless there are non-uniform values
+    # sa_m = specific_surface_area(img, method='marching_cubes')
+
+    assert (sa_f, sa_g) == (0, 0)
+
+
+def test_surface_area_on_floating_cube():
+    """Run surface area on floating cube"""
+    l = 20
+    img = np.zeros([l, l, l])
+    x1, x2 = 5, 15
+    img[x1:x2, x1:x2, x1:x2] = 1
+    sa_f = specific_surface_area(img, method='face_counting')['1']
+    sa_m = specific_surface_area(img, method='marching_cubes', smoothing=False, device='cpu')['1']
+    sa_g = specific_surface_area(img, method='gradient', smoothing=False)['1']
+
+    # All six sides should be taken into account
+    # True value is 0.075
+    assert (np.around(sa_f, 5), np.around(sa_m, 5), np.around(sa_g, 5)) == (0.075, 0.07051, 0.07085)
+
+
+def test_surface_area_on_corner_cube():
+    """Run surface area on corner cube"""
+    l = 20
+    img = np.zeros([l, l, l])
+    img[0:10, 0:10, 0:10] = 1
+    sa_f = specific_surface_area(img, method='face_counting')['1']
+
+    # Only three inner sides should be taken into account
+    # True value is 0.0375
+    assert np.around(sa_f, 5) == 0.0375
+
+
+def test_surface_area_on_sphere():
+    """Run surface area on sphere"""
+    l = 20
+    img = np.zeros([l, l, l])
+    radius = l*0.5-3
+    x, y, z = np.ogrid[:l, :l, :l]
+    distance_squared = (x - l/2 + 0.5)**2 + (y - l/2 + 0.5)**2 + (z - l/2 + 0.5)**2
+    mask = distance_squared <= radius**2
+    img[mask] = 1
+    a_theo = 4*np.pi*radius**2/img.size
+    sa_f = np.abs(specific_surface_area(img, method='face_counting')['1']-a_theo)/a_theo*100
+    sa_m = np.abs(specific_surface_area(img, method='marching_cubes', device='cpu')['1']-a_theo)/a_theo*100
+    sa_g = np.abs(specific_surface_area(img, method='gradient')['1']-a_theo)/a_theo*100
+
+    # Relative errors should be
+    # - face_counting: 52.01 %,
+    # - marching_cubes: 0.82 %,
+    # - gradient:       0.95 %
+    assert (np.around(sa_f, 2), np.around(sa_m, 2), np.around(sa_g, 2)) == (52.01, 0.82, 0.95)
+
+
+def test_surface_area_on_multi_cubes():
+    """Run surface area on multiple cubes"""
+    l = 20
+    img = np.zeros([l, l, l])
+    img[0:10, 0:10, 0:5] = 1
+    img[5:-5, 5:-5, 5:-5] = 2
+    img[0:10, 0:10, 15:] = 1
+    img[10:, 10:, 15:] = 3
+    sa_f = specific_surface_area(img, method='face_counting')
+    sa_m = specific_surface_area(img, method='marching_cubes', smoothing=False, device='cpu')
+    sa_g = specific_surface_area(img, method='gradient', smoothing=False)
+    results = []
+    for sa in [sa_f, sa_m, sa_g]:
+        for phase in ['1', '2', '3']:
+            results.append(np.around(sa[phase], 4))
+    # All six sides should be taken into account
+    # True value is 0.075
+    reference = [0.0500, 0.0750, 0.0250, \
+                 0.0422, 0.0705, 0.0211, \
+                 0.0482, 0.0709, 0.0241]
+    assert results == reference
+
+
+def test_tpb_2d():
+    """Run tpb on 3x3"""
+    l = 3
+    img = np.zeros([l, l])
+    img[0] = 1
+    img[:, 0] = 2
+    tpb = triple_phase_boundary(img)
+    assert tpb == 0.25
+
+
+def test_tpb_3d_corners():
+    """Run tpb on 2x2"""
+
+    l = 2
+    img = np.zeros([l, l, l])
+    img[0, 0, 0] = 1
+    img[1, 1, 1] = 1
+    img[0, 1, 1] = 2
+    img[1, 0, 0] = 2
+    tpb = triple_phase_boundary(img)
+    assert tpb == 1
+
+
+def test_tpb_3d_corners():
+    """Run tpb on 2x2 corners"""
+
+    l = 2
+    img = np.zeros([l, l, l])
+    img[0] = 1
+    img[:, 0] = 2
+    tpb = triple_phase_boundary(img)
+    assert tpb == 1/3
\ No newline at end of file
diff --git a/tests/test_taufactor.py b/tests/test_taufactor.py
index 46dade0..3eff302 100644
--- a/tests/test_taufactor.py
+++ b/tests/test_taufactor.py
@@ -2,17 +2,13 @@
 
 """Tests for `taufactor` package."""
 
-import pytest
 import taufactor as tau
-from taufactor.metrics import volume_fraction, surface_area, triple_phase_boundary
 import torch as pt
-from tests.utils import *
 import numpy as np
-import matplotlib.pyplot as plt
-import torch
-#  Testing the main solver
 
 
+#  Testing the main solver
+
 def test_solver_on_uniform_block():
     """Run solver on a block of ones."""
     l = 20
@@ -74,9 +70,10 @@ def test_solver_on_slanted_strip_of_ones():
     S = tau.Solver(img, device=pt.device('cpu'))
     S.solve()
     assert np.around(S.tau, decimals=5) == 4
-#  Testing the periodic solver
 
 
+#  Testing the periodic solver
+
 def test_deadend():
     """Test deadend pore"""
     solid  = np.zeros((10,50,50))
@@ -117,146 +114,8 @@ def test_periodic_solver_on_strip_of_ones():
     S.solve()
     assert np.around(S.tau, decimals=5) == 1
 
-# Testing the metrics
-# Volume fraction
-
-
-def test_volume_fraction_on_uniform_block():
-    """Run volume fraction on uniform block"""
-    l = 20
-    img = np.ones([l, l, l]).reshape(1, l, l, l)
-    vf = volume_fraction(img)['1']
-
-    assert np.around(vf, decimals=5) == 1.0
-
-
-def test_volume_fraction_on_empty_block():
-    """Run volume fraction on empty block"""
-    l = 20
-    img = np.zeros([l, l, l]).reshape(1, l, l, l)
-    vf = volume_fraction(img)['0']
-
-    assert np.around(vf, decimals=5) == 1.0
-
-
-def test_volume_fraction_on_checkerboard():
-    """Run volume fraction on checkerboard block"""
-    l = 20
-    img = generate_checkerboard(l)
-    vf = volume_fraction(img, phases={'zeros': 0, 'ones': 1})
-
-    assert (vf['zeros'], vf['ones']) == (0.5, 0.5)
-
-
-def test_volume_fraction_on_strip_of_ones():
-    """Run volume fraction on strip of ones"""
-    l = 20
-    img = np.zeros([l, l, l])
-    t = 10
-    img[:, 0:t, 0:t] = 1
-    vf = volume_fraction(img, phases={'zeros': 0, 'ones': 1})
-
-    assert (vf['zeros'], vf['ones']) == (0.75, 0.25)
-
-# Surface area
-
-
-def test_surface_area_on_uniform_block():
-    """Run surface area on uniform block"""
-    l = 20
-    img = np.ones([l, l, l])
-    sa = surface_area(img, phases=[1])
-
-    assert sa == 0
-
-
-def test_surface_area_on_empty_block():
-    """Run surface area on empty block"""
-    l = 20
-    img = np.zeros([l, l, l])
-    sa = surface_area(img, phases=[0])
-
-    assert sa == 0
-
-def test_surface_area_on_checkerboard():
-    """Run surface area on checkerboard block"""
-    l = 50
-    img = generate_checkerboard(l)
-    sa = surface_area(img, phases=[0, 1])
-
-    assert sa == 1
-
-
-def test_surface_area_on_strip_of_ones():
-    """Run surface area on single one in small 2x2z2 cube"""
-    l = 2
-    img = np.zeros([l, l, l])
-    t = 1
-    img[0, 0, 0] = 1
-    sa = surface_area(img, phases=[0, 1])
-
-    assert sa == 0.25
-
-
-def test_surface_area_on_non_periodic_2d():
-    """Run surface area on a pair of one in small 3x3 square"""
-    img = np.array([[0, 0, 0], [1, 1, 0], [0, 0, 0]])
-    sa = surface_area(img, phases=[0, 1])
-
-    assert sa == 5/12
-
-
-def test_surface_area_on_periodic_2d():
-    """Run surface area on a pair of one in small 3x3 square"""
-    img = np.array([[0, 0, 0], [1, 1, 0], [0, 0, 0]])
-    sa = surface_area(img, phases=[0, 1], periodic=[0, 1])
-
-    assert sa == 6/18
-
-
-def test_surface_area_interfactial_3ph():
-    """Run surface area 3 phases """
-    l = 3
-    img = np.zeros([l, l, l])
-    img[1] = 1
-    img[2] = 2
-    sa = surface_area(img, phases=[1, 2])
-    assert sa == 1/6
-
-
-def test_tpb_2d():
-    """Run tpb on 3x3"""
-    l = 3
-    img = np.zeros([l, l])
-    img[0] = 1
-    img[:, 0] = 2
-    tpb = triple_phase_boundary(img)
-    assert tpb == 0.25
-
-
-def test_tpb_3d_corners():
-    """Run tpb on 2x2"""
-
-    l = 2
-    img = np.zeros([l, l, l])
-    img[0, 0, 0] = 1
-    img[1, 1, 1] = 1
-    img[0, 1, 1] = 2
-    img[1, 0, 0] = 2
-    tpb = triple_phase_boundary(img)
-    assert tpb == 1
-
-
-def test_tpb_3d_corners():
-    """Run tpb on 2x2 corners"""
-
-    l = 2
-    img = np.zeros([l, l, l])
-    img[0] = 1
-    img[:, 0] = 2
-    tpb = triple_phase_boundary(img)
-    assert tpb == 1/3
 
+#  Testing the multiphase solver
 
 def test_multiphase_and_solver_agree():
     """test mph and solver agree when Ds are the same"""
@@ -342,6 +201,8 @@ def test_mphsolver_on_strip_of_ones_and_twos_and_threes():
     assert np.around(S.tau, 4) == 1
 
 
+#  Testing the tau_e solver
+
 def test_taue_deadend():
     """Run solver on a deadend strip of ones"""
 

From b9fa69ed6b477debd9704a02245fc9f830a0e599 Mon Sep 17 00:00:00 2001
From: daubners <s.daubner@imperial.ac.uk>
Date: Tue, 21 Jan 2025 21:05:17 +0000
Subject: [PATCH 10/12] fix cpu default

---
 taufactor/metrics.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/taufactor/metrics.py b/taufactor/metrics.py
index b4eca0c..ad86c50 100644
--- a/taufactor/metrics.py
+++ b/taufactor/metrics.py
@@ -57,7 +57,7 @@ def crop_area_of_interest_numpy(array, labels):
                       max(min_idx[2]-3, 0):min(max_idx[2]+4, array.shape[2])]
     return sub_array
 
-def gaussian_kernel_3d_torch(size=3, sigma=1.0, device=torch.device('cuda')):
+def gaussian_kernel_3d_torch(size=3, sigma=1.0, device):
     """Creates a 3D Gaussian kernel using PyTorch"""
     ax = torch.linspace(-(size // 2), size // 2, size)
     xx, yy, zz = torch.meshgrid(ax, ax, ax, indexing="ij")
@@ -90,6 +90,10 @@ def specific_surface_area(img, spacing=(1,1,1), phases={}, method='gradient', de
     [dx,dy,dz] = spacing
     surface_areas = {}
 
+    if torch.device(device).type.startswith('cuda') and not torch.cuda.is_available():
+        device = torch.device('cpu')
+        warnings.warn("CUDA not available, defaulting device to cpu.")
+
     if (method == 'gradient') | (method == 'face_counting'):
         if type(img) is not type(torch.tensor(1)):
             tensor = torch.tensor(img)
@@ -102,7 +106,7 @@ def specific_surface_area(img, spacing=(1,1,1), phases={}, method='gradient', de
             labels = torch.unique(tensor)
             labels = labels.int()
             phases = {str(label.item()): label.item() for label in labels}
-        gaussian = gaussian_kernel_3d_torch(device=device)
+        gaussian = gaussian_kernel_3d_torch(device)
 
         volume = torch.numel(tensor)
         for name, label in phases.items():

From f4ca7abf4789e112f0ea98bc0d70a1de20bd4499 Mon Sep 17 00:00:00 2001
From: daubners <s.daubner@imperial.ac.uk>
Date: Tue, 21 Jan 2025 21:08:00 +0000
Subject: [PATCH 11/12] argument order

---
 taufactor/metrics.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/taufactor/metrics.py b/taufactor/metrics.py
index ad86c50..8dad6fa 100644
--- a/taufactor/metrics.py
+++ b/taufactor/metrics.py
@@ -57,7 +57,7 @@ def crop_area_of_interest_numpy(array, labels):
                       max(min_idx[2]-3, 0):min(max_idx[2]+4, array.shape[2])]
     return sub_array
 
-def gaussian_kernel_3d_torch(size=3, sigma=1.0, device):
+def gaussian_kernel_3d_torch(device, size=3, sigma=1.0):
     """Creates a 3D Gaussian kernel using PyTorch"""
     ax = torch.linspace(-(size // 2), size // 2, size)
     xx, yy, zz = torch.meshgrid(ax, ax, ax, indexing="ij")

From 56a281e5a39df8cf20690721930ec65d59004722 Mon Sep 17 00:00:00 2001
From: daubners <s.daubner@imperial.ac.uk>
Date: Tue, 21 Jan 2025 21:19:09 +0000
Subject: [PATCH 12/12] add skimage

---
 docs/requirements.txt | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/requirements.txt b/docs/requirements.txt
index 41b7a9b..abcd00d 100644
--- a/docs/requirements.txt
+++ b/docs/requirements.txt
@@ -1,6 +1,7 @@
 numpy
 torch
 scipy
+scikit-image
 matplotlib
 tifffile
 myst_parser