diff --git a/osl_ephys/report/src_report.py b/osl_ephys/report/src_report.py
index a0f2075..07e5445 100644
--- a/osl_ephys/report/src_report.py
+++ b/osl_ephys/report/src_report.py
@@ -106,6 +106,10 @@ def gen_html_data(config, outdir, subject, reportdir, logger=None, extra_funcs=N
         data["plt_filters_svd"] = f"{subject}/filters_svd.png"
         copy("{}/{}".format(outdir, data["filters_svd_plot"]), "{}/{}/filters_svd.png".format(reportdir, subject))
 
+    if "dipole_locations_plot" in data:
+        data["plt_dipole_locations"] = f"{subject}/dipole_locations.html" # this might not have been rendered yet (if using dask), so potentially need to copy it over later
+        copy("{}/{}".format(outdir, data["dipole_locations_plot"]), "{}/{}/dipole_locations.html".format(reportdir, subject))
+
     if "parc_psd_plot" in data:
         data["plt_parc_psd"] = f"{subject}/parc_psd.png"
         copy("{}/{}".format(outdir, data["parc_psd_plot"]), "{}/{}/parc_psd.png".format(reportdir, subject))
@@ -157,6 +161,7 @@ def gen_html_page(reportdir):
     # Load HTML data
     data = []
     for subdir in subdirs:
+        
         subdir = Path(subdir)
         # Just generate the html page with the successful runs
         try:
@@ -464,6 +469,7 @@ def add_to_data(data_file, info):
         Info to add.
     """
     data_file = Path(data_file)
+
     if data_file.exists():
         data = pickle.load(open(data_file, "rb"))
     else:
diff --git a/osl_ephys/report/templates/src_subject_panel.html b/osl_ephys/report/templates/src_subject_panel.html
index 37752c1..d1466be 100644
--- a/osl_ephys/report/templates/src_subject_panel.html
+++ b/osl_ephys/report/templates/src_subject_panel.html
@@ -105,18 +105,32 @@ <h3>Coregistration
       {% endif %}
   
       {% if data.beamform or data.beamform_and_parcellate %}
-          <div class="tabpage" style='width: 100%; display: none' id={{ data.fif_id }}_beamforming>
-              <h3>Beamforming</h3>
-              <span style="margin-right: 10px;"></span>
-              <button style="cursor: help; border: 1px solid; background-color: lightblue; color: black; padding: 2px 5px; font-size: 12px;" 
-                      title="Does the covariance look sensible (i.e., not dominated by a single channel pair)? Do the eigenvalues drop off gradually before suddenly dropping to zero? Is the rank what you'd expect?">
-                  ?
-              </button>
-              <div>LCMV Filter</div>
-              <div style="width: 70%; margin: auto">
-                  <img src="{{ data.plt_filters_cov }}" alt="" style='max-width: 100%'/>
-                  <img src="{{ data.plt_filters_svd }}" alt="" style='max-width: 100%'/>
-              </div>
+          <div class="tabpage" style='width: 100%' id={{ data.fif_id }}_beamforming>
+                <h3>Beamforming</h3>
+                <span style="margin-right: 10px;"></span>
+                <button style="cursor: help; border: 1px solid; background-color: lightblue; color: black; padding: 2px 5px; font-size: 12px;" 
+                        title="Does the covariance look sensible (i.e., not dominated by a single channel pair)? Do the eigenvalues drop off gradually before suddenly dropping to zero? Is the rank what you'd expect?">
+                    ?
+                </button>
+                <div>LCMV Filter</div>
+                <div style="width: 70%; margin: auto">
+                    <img src="{{ data.plt_filters_cov }}" alt="" style='max-width: 100%'/>
+                    <img src="{{ data.plt_filters_svd }}" alt="" style='max-width: 100%'/>
+                </div>
+                <div>Dipole types and locations</div>
+                <iframe data-src="{{ data.plt_dipole_locations }}" width="525" height="525" id="{{ data.fif_id}}_beamforming_iframe"></iframe>
+                <div>
+                    Click and drag to rotate and scroll to zoom.</br>
+                    Blue dots are dipoles that will be treated as single dipoles.</br>
+                    </br>
+                    If a (non-patched) bilateral beamformer is being used:</br>
+                    - Red lines connect dipole "pairs" that will be beamformed together.</br>
+                    - Green dots are "midline" dipoles that will be treated as single dipoles.</br>
+                    </br>
+                    If a patch (multi-dipole) beamformer is being used:</br>
+                    - Multi-coloured dots are dipoles that belong to multi-dipole patches,</br>
+                    the color indicates the patch each dipole belongs to.
+                </div>
           </div>
       {% endif %}
   
@@ -167,7 +181,6 @@ <h4>Error Log</h4>
 </style>
 
 
-
 <script>
 (function () {
   // initialize only once even if template rendered many times
@@ -517,3 +530,4 @@ <h4>Error Log</h4>
 })();
 </script>
 <!-- END: osl-fif stacking fix -->
+
diff --git a/osl_ephys/source_recon/beamforming.py b/osl_ephys/source_recon/beamforming.py
index 2e58959..b78c908 100644
--- a/osl_ephys/source_recon/beamforming.py
+++ b/osl_ephys/source_recon/beamforming.py
@@ -10,6 +10,9 @@
 
 import numpy as np
 import matplotlib.pyplot as plt
+
+from scipy.spatial import KDTree
+
 import mne
 from mne import read_forward_solution, Covariance, compute_covariance, compute_raw_covariance
 from mne.rank import compute_rank
@@ -19,7 +22,7 @@
 from mne.forward.forward import is_fixed_orient
 from mne.beamformer import read_beamformer
 from mne.beamformer._lcmv import _apply_lcmv
-from mne.beamformer._compute_beamformer import _reduce_leadfield_rank, _sym_inv_sm, Beamformer
+from mne.beamformer._compute_beamformer import _reduce_leadfield_rank, Beamformer
 from mne.minimum_norm.inverse import _check_reference
 from mne.utils import (
     _check_channels_spatial_filter,
@@ -35,6 +38,7 @@
     warn,
 )
 from mne.utils import logger as mne_logger
+from osl_ephys.source_recon import parcellation
 
 try:
     from mne._fiff.meas_info import _simplify_info
@@ -46,10 +50,53 @@
     from mne.io.pick import pick_channels_cov, pick_info
     from mne.io.proj import make_projector
 
-from osl_ephys.source_recon import rhino
+import nibabel as nib
+from osl_ephys.source_recon.rhino import coreg
 from osl_ephys.source_recon.rhino import utils as rhino_utils
+from osl_ephys.source_recon.utils import find_file
 from osl_ephys.utils.logger import log_or_print
+from osl_ephys.source_recon.parcellation import assign_voxels_to_binary_parcellation  
 
+def _sym_inv_sm(x, reduce_rank, inversion, sk):
+    """Symmetric inversion with single- or matrix-style inversion.
+    osl-ephys version of the same fn in MNE
+    Modified to cope with bilateral beamformer, which creates the situation 
+    where x.shape[1:] == (2, 2).
+    """
+    if x.shape[1:] == (1, 1):
+        with np.errstate(divide="ignore", invalid="ignore"):
+            x_inv = 1.0 / x
+        x_inv[~np.isfinite(x_inv)] = 1.0
+
+    elif x.shape[1:] == (3, 3):     
+        #assert x.shape[1:] == (3, 3)
+        if inversion == "matrix":
+            x_inv = _sym_mat_pow(x, -1, reduce_rank=reduce_rank)
+            # Reapply source covariance after inversion
+            x_inv *= sk[:, :, np.newaxis]
+            x_inv *= sk[:, np.newaxis, :]
+            
+        else:
+            # Invert for each dipole separately using plain division
+            diags = np.diagonal(x, axis1=1, axis2=2)
+            assert not reduce_rank  # guaranteed earlier
+            with np.errstate(divide="ignore"):
+                diags = 1.0 / diags
+            # set the diagonal of each 3x3
+            x_inv = np.zeros_like(x)
+            for k in range(x.shape[0]):
+                this = diags[k]
+                # Reapply source covariance after inversion
+                this *= sk[k] * sk[k]
+                x_inv[k].flat[::4] = this
+    # MWW
+    else:
+        x_inv = _sym_mat_pow(x, -1, reduce_rank=reduce_rank)
+        if inversion != "matrix":
+            raise ValueError("inversion must be 'matrix' when x.shape[1:] != (1, 1) or (3, 3)")
+    # MWW end
+    
+    return x_inv
 
 def get_beamforming_filenames(subjects_dir, subject):
     """Get beamforming filenames.
@@ -98,6 +145,11 @@ def make_lcmv(
     reduce_rank=True,
     depth=None,
     inversion="matrix",
+    use_bilateral_pairs=None,
+    bilateral_tol=None,
+    bilateral_tol_midline=None,
+    dipole_patches_file=None,
+    group_patches=False,
     verbose=None,
     save_filters=False,
 ):
@@ -153,6 +205,41 @@ def make_lcmv(
         How to weight (or normalize) the forward using a depth prior (see MNE docs).
     inversion : 'matrix' | 'single'
         The inversion scheme to compute the weights.
+    use_bilateral_pairs : bool
+        If True, find bilateral pairs of dipoles and use a bilateral beamformer.
+    bilateral_tol : float
+        The distance threshold for finding pairs, in mm.
+        Only used if use_bilateral_pairs is True.
+        Recommended to be ~ gridstep / 2
+    bilateral_tol_midline : float
+        The minimum distance in mm from the midline for a dipole to be
+        considered for bilateral pairing. Only used if use_bilateral_pairs
+        is True. Recommended value is ~gridstep / 2  mm.
+        If None then set to bilateral_tol.
+    dipole_patches_file : str | None
+        Path to a 4D nifti binary parcellation file where 
+        each volume should be a mask with 0 for the background 
+        and 1 for the parcel.
+        Dipoles within a parcel will be considered as a group
+        of dipoles that will be beamformed together as a multi-dipole.
+        If None, then no patches are used in the beamformer.
+        File must be in MNI space.
+        This only works with volumetric (RHINO) source spaces.
+        If group_patches=True is also specified,
+        then the dipole_patches_file needs to be
+        accompanied by a text file with the same name but a .xml
+        extension that contains a "group" field
+        with a unique integer value shared by patches that are 
+        bilateral pairs (if group=0 then the patch 
+        is assumed to not be in a bilateral pair).
+        (see also parcellation.load_parcellation_description).
+        Using a dipole_patches_file only works with volumetric (RHINO)
+        source spaces.
+    group_patches : bool
+        If True, group patches using the "group" field in the dipole patches description (
+        xml) file.
+        If False, do not group patches.
+        Only required if dipole_patches_file is not None.
     save_filters : bool
         Should we save the LCMV beamforming filter?
 
@@ -163,12 +250,16 @@ def make_lcmv(
     """
     log_or_print("*** RUNNING OSL MAKE LCMV ***")
 
-    rhino_files = rhino_utils.get_rhino_filenames(subjects_dir, subject)
     beamforming_files = get_beamforming_filenames(subjects_dir, subject)
+    rhino_files = rhino_utils.get_rhino_filenames(subjects_dir, subject)
 
-    # load forward solution
-    fwd_fname = rhino_files["fwd_model"]
-    fwd = read_forward_solution(fwd_fname)
+    fwd = read_forward_solution(rhino_files["fwd_model"])
+
+    src_pnts_inuse_mni = None
+    if use_bilateral_pairs or (dipole_patches_file is not None):
+        src_pnts_inuse_mni = _get_src_pnts_inuse_mni(subjects_dir, 
+                                                 subject, 
+                                                 src_index=0)
 
     if data_cov is None:
         # Note that if chantypes are meg, eeg; and meg includes mag, grad then compute_covariance will project data separately for meg and eeg to reduced
@@ -225,6 +316,12 @@ def make_lcmv(
         rank=rank,
         noise_rank=noise_rank,
         reduce_rank=reduce_rank,
+        use_bilateral_pairs=use_bilateral_pairs,
+        bilateral_tol=bilateral_tol,
+        bilateral_tol_midline=bilateral_tol_midline,
+        dipole_patches_file=dipole_patches_file,
+        group_patches=group_patches,
+        src_pnts_inuse_mni=src_pnts_inuse_mni,
         verbose=verbose,
     )
 
@@ -238,7 +335,7 @@ def make_lcmv(
 
 
 def make_plots(subjects_dir, subject, filters, data):
-    """Plot LCMV beamforming filters.
+    """Plot LCMV beamforming info.
 
     Parameters
     ----------
@@ -262,6 +359,7 @@ def make_plots(subjects_dir, subject, filters, data):
     fig_cov, fig_svd = filters["data_cov"].plot(data.info, show=False, verbose=False)
     fig_cov.savefig(beamforming_files["filters_plot_cov"], dpi=150)
     fig_svd.savefig(beamforming_files["filters_plot_svd"], dpi=150)
+
     plt.close("all")
     return beamforming_files["filters_plot_cov"], beamforming_files["filters_plot_svd"]
 
@@ -346,7 +444,7 @@ def apply_lcmv_raw(raw, filters, reject_by_annotation="omit"):
     return next(stc)
 
 
-def get_recon_timeseries(subjects_dir, subject, coord_mni, recon_timeseries_head):
+def get_recon_timeseries(subjects_dir, subject, coord_mni, recon_timeseries_head, src_index=0):
     """Gets the reconstructed time series nearest to the passed coordinates in MNI space.
 
     Parameters
@@ -359,7 +457,8 @@ def get_recon_timeseries(subjects_dir, subject, coord_mni, recon_timeseries_head
         3D coordinate in MNI space to get timeseries for
     recon_timeseries_head : (ndipoles, ntpts) np.array
         Reconstructed time courses in head (polhemus) space. Assumes that the dipoles are the same (and in the same order) as those in the forward model, rhino_files['fwd_model'].
-
+    src_index : int
+        Index of source space in forward model to use.
     Returns
     -------
     recon_timeseries : numpy.ndarray
@@ -376,8 +475,8 @@ def get_recon_timeseries(subjects_dir, subject, coord_mni, recon_timeseries_head
 
     # Get hold of coords of points reconstructed to. Note, MNE forward model is done in head space in metres. RHINO does everything in mm
     fwd = read_forward_solution(rhino_files["fwd_model"])
-    vs = fwd["src"][0]
-    recon_coords_head = vs["rr"][vs["vertno"]] * 1000  # in mm
+    vs = fwd["src"][src_index]
+    recon_coords_head = vs["rr"][vs["vertno"]] * 1000  # RHINO is in mm
 
     # Convert coords_head from head to mri space to get index of reconstructed coordinate nearest to coord_mni
     head_scaledmri_t = rhino_utils.read_trans(coreg_filenames["head_scaledmri_t_file"])
@@ -396,6 +495,7 @@ def transform_recon_timeseries(
     recon_timeseries,
     spatial_resolution=None,
     reference_brain="mni",
+    src_index=0,
 ):
     """Spatially resamples a (ndipoles x ntpts) array of reconstructed time courses (in head/polhemus space) to dipoles on the brain grid of the specified reference brain.
 
@@ -416,6 +516,8 @@ def transform_recon_timeseries(
         'mri' indicates that the reference_brain is the subject's sMRI in the scaled native/mri space.
         'unscaled_mri' indicates that the reference_brain is the subject's sMRI in unscaled native/mri space.
         Note that scaled/unscaled relates to the allow_smri_scaling option in coreg. If allow_scaling was False, then the unscaled MRI will be the same as the scaled MRI.
+    src_index : int
+        Index of source space in forward model to use.
 
     Returns
     -------
@@ -434,15 +536,16 @@ def transform_recon_timeseries(
     # -------------------------------------------------------------------------------------
     # Get hold of coords of points reconstructed to
     #
-    # Note, MNE forward model is done in head space in metres. RHINO does everything in mm.
+    # Note, MNE forward model is done in head space in metres
+    # RHINO does everything in mm.
     fwd = read_forward_solution(rhino_files["fwd_model"])
-    vs = fwd["src"][0]
+    vs = fwd["src"][src_index]
     recon_coords_head = vs["rr"][vs["vertno"]] * 1000  # in mm
 
     # ----------------------------
     if spatial_resolution is None:
         # Estimate gridstep from forward model
-        rr = fwd["src"][0]["rr"]
+        rr = fwd["src"][src_index]["rr"]
 
         store = []
         for ii in range(rr.shape[0]):
@@ -668,6 +771,356 @@ def get_leadfields(
     return leadfields.T, coords
 
 
+def _find_patches(src_pnts_inuse_mni, 
+                  dipole_patches_file,
+                  group_patches=False):   
+
+    '''
+    
+    Parameters
+    ----------
+    src_pnts_inuse_mni: np.array
+        The 3D coords on the source/dipole points being used.
+        In MNI space in mm.
+        [ndipoles, 3]
+    dipole_patches_file : str
+        Path to a 4D nifti binary parcellation file where 
+        each volume should be a mask with 0 for the background 
+        and 1 for the parcel.
+        Dipoles within a parcel will be considered as a group
+        of dipoles that will be beamformed together as a multi-dipole.
+        File must be in MNI space.
+        This only works with volumetric (RHINO) source spaces.
+    group_patches : bool
+        If True, group patches using the "group" field in the dipole patches description (
+        xml) file.
+        If False, do not group patches.
+    '''
+
+    parcellation_data = nib.load(find_file(dipole_patches_file)).get_fdata()
+
+    # Check if parcellation is binary 
+    binary_mask = np.array_equal(np.unique(parcellation_data), [0, 1])
+    if not binary_mask:
+        raise ValueError("dipole_patches_file must be a binary mask (0 for background, 1 for parcel)")  
+
+    # Check that each voxel in the parcellation belongs to only one parcel
+    nparcels = parcellation_data.shape[3]
+    overlapping_volumes = []
+    for vv in range(nparcels):
+        for ww in range(vv+1, nparcels):
+            overlapping_voxels = np.where((parcellation_data[:, :, :, vv] > 0) & (parcellation_data[:, :, :, ww] > 0))
+            if len(overlapping_voxels[0]) > 0:
+                overlapping_volumes.append((vv, ww))
+                log_or_print(f"Found overlapping volumes: {vv}, {ww}")
+    if len(overlapping_volumes) > 0:
+        raise ValueError("dipole_patches_file must be a parcellation where each voxel belongs to only one parcel.")
+    
+    parcellation_assignment = assign_voxels_to_binary_parcellation(src_pnts_inuse_mni, dipole_patches_file)
+
+    multi_dipoles = []
+    n_dipoles_in_patches = 0
+
+    for ppp in range(nparcels):
+        dipoles_in_this_patch = np.where(np.array(parcellation_assignment) == ppp)[0]
+        n_dipoles_in_patches += len(dipoles_in_this_patch)
+        if len(dipoles_in_this_patch) > 0:
+            multi_dipoles.append(dipoles_in_this_patch.tolist())
+
+    single_dipoles = np.setdiff1d(np.arange(src_pnts_inuse_mni.shape[0]), np.unique(np.concatenate(multi_dipoles))) 
+
+    log_or_print(f"Number of patches: {len(multi_dipoles)}")
+    log_or_print(f"Total number of dipoles in patches: {n_dipoles_in_patches} out of {src_pnts_inuse_mni.shape[0]}")    
+    log_or_print(f"Number of single dipoles (not in patches): {len(single_dipoles)}")
+    log_or_print(f"Total number of dipoles (in patches + single_dipoles): {n_dipoles_in_patches + len(single_dipoles)} out of {src_pnts_inuse_mni.shape[0]}")
+
+    if group_patches:
+
+        log_or_print(f"Loading {dipole_patches_file} to find group field")
+        parcellation_description = parcellation.load_parcellation_description(dipole_patches_file)
+        group_indices = []
+        for index in parcellation_description.keys():
+            if parcellation_description[index]["group"] > 0:
+                group_indices.append(parcellation_description[index]["group"])
+            else:
+                group_indices.append(0)
+
+        group_indices = np.array(group_indices)
+
+        unique_group_indices = np.unique(group_indices)
+        unique_group_indices = unique_group_indices[unique_group_indices != 0]
+
+        log_or_print(f"Found {len(unique_group_indices)} unique group indices")
+
+        new_multi_dipoles = []
+        for gpi in unique_group_indices:
+            patches_with_this_gpi = np.where(group_indices == gpi)[0]
+
+            # Find the dipoles in the patches in this group
+            dipoles_in_this_group = []
+            for pp in patches_with_this_gpi:
+                dipoles_in_this_group += multi_dipoles[pp]
+
+            new_multi_dipoles.append(dipoles_in_this_group)
+
+        multi_dipoles = new_multi_dipoles
+
+        log_or_print(f"After grouping patches, number of multi-dipoles is now {len(multi_dipoles)}")
+
+    return multi_dipoles, single_dipoles
+
+
+def _find_bilateral_pairs(src_pnts_inuse_mni, 
+                          bilateral_tol=None, 
+                          bilateral_tol_midline=None):
+
+    """
+
+    Find pairs of bilateral dipoles in forward["src"][src_index]
+    
+    Parameters
+    ----------
+    src_pnts_inuse_mni: np.array
+        The 3D coords on the source/dipole points being used.
+        In MNI space in mm.
+        [ndipoles, 3]
+    bilateral_tol : float
+        The distance threshold for finding pairs, in mm
+        If None then attempt to estimate gridstep from src_pnts_inuse_mni
+        and set bilateral_tol = gridstep / 2.
+    bilateral_tol_midline : float
+        The distance threshold for finding midline points, in mm
+        If None then bilateral_tol_midline = bilateral_tol.
+
+    Returns
+    -------
+    multi_dipoles : list of lists
+        Each entry in the list is a list of indices (wrt to src_pnts_inuse_mni) for dipoles 
+        that should be treated as a multi-dipole.
+        I.e. each entry is a list of length 2, containing the indices for the two dipoles 
+        that make up a bilateral pair.
+    single_dipoles : (nsingle_dipoles,) numpy.ndarray
+        List of indices (wrt to src_pnts_inuse_mni) for single points that could not be paired.
+    midline_points : (nmidline,) numpy.ndarray
+        List of indices (wrt to src_pnts_inuse_mni) for points close to the midline.
+
+    """
+
+    # Assumes src_pnts_inuse_mni are in MNI space in mm
+    # x-axis indexes left vs right hemisphere
+    #
+    # Aim is to identify pairs of matching points that are in opposite hemispheres 
+    # (I.e. x values are the same magnitude, but one is positive , and the other is negative - within a tolerance)
+
+    if bilateral_tol is None:
+        gridstep = int(rhino_utils.get_gridstep(src_pnts_inuse_mni)/1000)
+        bilateral_tol = gridstep / 2
+        log_or_print(f"gridstep = {gridstep} mm")
+        log_or_print(f"Setting bilateral_tol = {bilateral_tol} mm")
+
+    if bilateral_tol <= 0:
+        raise ValueError("bilateral_tol must be a positive scalar.")
+
+    if bilateral_tol_midline is None:
+        bilateral_tol_midline = bilateral_tol
+        log_or_print(f"Setting bilateral_tol_midline = {bilateral_tol_midline} mm")
+
+    pairs = []
+    single_dipoles = []
+
+    # find all src points that are close to the midline
+    midline_points = np.where(np.abs(src_pnts_inuse_mni[:, 0]) < bilateral_tol_midline)[0]
+    non_midline_points = np.where(np.abs(src_pnts_inuse_mni[:, 0]) >= bilateral_tol_midline)[0]
+    src_inuse_nomidline = src_pnts_inuse_mni[non_midline_points, :]
+
+    # compute distances between all remaining points cross-hemisphere
+    # i.e. between points with positive x and points with negative x
+    left_points = np.where(src_inuse_nomidline[:, 0] < 0)[0]
+    right_points = np.where(src_inuse_nomidline[:, 0] >= 0)[0]
+    left_coords = src_inuse_nomidline[left_points, :]
+    right_coords = src_inuse_nomidline[right_points, :]
+    right_coords[:, 0] = -right_coords[:, 0]  # mirror right hemisphere points to left hemisphere
+
+    # compute distance matrix
+    dist_matrix = np.zeros((left_coords.shape[0], right_coords.shape[0]))
+    for ii in range(left_coords.shape[0]):
+        for jj in range(right_coords.shape[0]):
+            dist_matrix[ii, jj] = np.abs(left_coords[ii, :] - right_coords[jj, :]).mean()
+
+    # find best matching pairs using the distance matrix
+    while True:
+        min_dist = np.min(dist_matrix)
+        if min_dist > bilateral_tol:
+            break
+        min_idx = np.unravel_index(np.argmin(dist_matrix), dist_matrix.shape)
+        pairs.append((min_idx[0], min_idx[1]))
+        dist_matrix[min_idx[0], :] = np.inf
+        dist_matrix[:, min_idx[1]] = np.inf
+        #print(f"Found pair: {left_points[min_idx[0]]} <-> {right_points[min_idx[1]]} with distance {min_dist}") 
+        
+    # find single_dipoles
+    used_left = np.zeros(left_coords.shape[0])
+    used_right = np.zeros(right_coords.shape[0])
+    for pp in pairs:
+        used_left[pp[0]] = 1
+        used_right[pp[1]] = 1
+
+    # midline points already in src_inuse indexing
+
+    # put single_dipoles into src_inuse indexing
+    for ii in range(left_coords.shape[0]):
+        if used_left[ii] == 0:
+            single_dipoles.append(non_midline_points[left_points[ii]])
+    for ii in range(right_coords.shape[0]):
+        if used_right[ii] == 0:
+            single_dipoles.append(non_midline_points[right_points[ii]])
+
+    # put pairs into src_inuse indexing
+    new_pairs = []
+    for pp in pairs:
+        new_pairs.append((non_midline_points[left_points[pp[0]]], 
+                          non_midline_points[right_points[pp[1]]]))                  
+    pairs = new_pairs
+
+    # convert pairs to numpy array
+    pairs = np.array(pairs) 
+    single_dipoles = np.array(single_dipoles)
+    midline_points = np.array(midline_points)
+    
+    log_or_print(f"Found {len(pairs)} pairs and {len(single_dipoles)} single-dipoles and {len(midline_points)} midline points")
+    log_or_print(f"Total points used: {len(pairs)*2 + len(single_dipoles) + len(midline_points)} out of {src_pnts_inuse_mni.shape[0]}")
+
+    # put pairs into list of lists
+    multi_dipoles = []
+    for pp in pairs:
+        multi_dipoles.append([pp[0], pp[1]])
+        
+    return multi_dipoles, single_dipoles, midline_points
+
+def _get_src_pnts_inuse_mni(outdir, subject, src_index):
+
+    '''
+    Get source points that are in use in MNI coordinates.
+    
+    Parameters
+    ----------
+    subjects_dir : string
+        Directory containing subject data.
+    subject : string
+        Subject identifier.
+    head_pnts : array
+        Points in the head coordinate system in mm
+
+    Returns
+    -------
+    mni_pnts : array
+        Points in the MNI coordinate system in mm
+    '''
+
+    # load forward solution
+    rhino_files = rhino_utils.get_rhino_filenames(outdir, subject)
+    fwd_fname = rhino_files["fwd_model"]
+    fwd = read_forward_solution(fwd_fname)
+    if fwd["coord_frame"] != 4:
+        raise ValueError("Source space must be in head coordinates (coord_frame=4).")        
+    src = fwd["src"][src_index]
+    if src['type'] != 'vol':
+        ValueError('Forward model, src[''type''], needs to be vol')
+    src_pnts_inuse = src['rr'][src['inuse']==1, :] * 1000  # convert to mm (MNE works in metres, RHINO works in mm)
+
+    # move to MNI coordinates to find pairs
+    src_pnts_inuse_mni = coreg.apply_head2mni_xform(outdir, subject, src_pnts_inuse) 
+
+    return src_pnts_inuse_mni
+
+def plot_dipole_locations(outdir, 
+                          subject, 
+                          use_bilateral_pairs, 
+                          bilateral_tol, 
+                          bilateral_tol_midline, 
+                          dipole_patches_file,
+                          group_patches,
+                          filename,
+                          src_index=0):
+
+    '''    
+    Parameters
+    ----------
+    outdir : string
+        Directory containing the subject directories.
+    subject : string
+        Subject name.
+    use_bilateral_pairs : bool
+        If True, display bilateral pairs of dipoles.
+        Overridden by dipole_patches_file.
+    bilateral_tol : float
+        The distance threshold for finding pairs, in mm
+        Recommended to be ~ gridstep / 2
+    bilateral_tol_midline : float
+        The distance threshold for finding midline points, in mm
+        Recommended to be ~ gridstep / 2
+        If None then bilateral_tol_midline = bilateral_tol.
+    dipole_patches_file : string
+        File name of dipole patches description (xml) file.
+        If None, do not use dipole patches.
+        Overides use_bilateral_pairs.
+    group_patches : bool
+        If True, group patches using the "group" field in the dipole patches description (xml) file.
+        If False, do not group patches.
+    filename : string
+        file to save plot to
+    src_index : int
+        The index of the source space to use in forward['src'].        
+    '''
+
+    if dipole_patches_file is not None:
+        log_or_print("Plotting dipole patches from {}".format(dipole_patches_file))
+        
+        src_pnts_inuse_mni = _get_src_pnts_inuse_mni(outdir, subject, src_index)
+        multi_dipoles, single_dipoles = _find_patches(src_pnts_inuse_mni, 
+                                          dipole_patches_file,
+                                          group_patches=group_patches,
+                                          )       
+        midline_points = None
+
+        multi_dipoles_display_mode = "multicoloured_points" 
+
+    elif use_bilateral_pairs:   
+
+        log_or_print("Plotting bilateral pairs")
+
+        src_pnts_inuse_mni = _get_src_pnts_inuse_mni(outdir, subject, src_index)
+        multi_dipoles, single_dipoles, midline_points = _find_bilateral_pairs(src_pnts_inuse_mni, 
+                                                                  bilateral_tol=bilateral_tol, 
+                                                                  bilateral_tol_midline=bilateral_tol_midline,
+                                                                 )      
+        multi_dipoles_display_mode = "linepairs"
+
+    else: 
+        multi_dipoles, single_dipoles, midline_points = None, None, None
+
+        multi_dipoles_display_mode = None
+
+    coreg.coreg_display(
+        outdir,
+        subject,
+        plot_type="surf",
+        display_outskin=False,
+        display_sensors=False,
+        display_inskull=True,
+        display_dipoles=True,
+        display_sensor_oris=False,
+        display_fiducials=False,
+        display_headshape_pnts=False,        
+        multi_dipoles_display_mode=multi_dipoles_display_mode,
+        multi_dipoles=multi_dipoles,
+        single_dipoles=single_dipoles,
+        midline_points=midline_points,
+        filename=filename,
+        src_index=src_index,
+    )
+
+
 @verbose
 def _make_lcmv(
     info,
@@ -683,6 +1136,12 @@ def _make_lcmv(
     reduce_rank=False,
     depth=None,
     inversion="matrix",
+    use_bilateral_pairs=False,
+    bilateral_tol=None,
+    bilateral_tol_midline=None,
+    dipole_patches_file=None,
+    group_patches=False,
+    src_pnts_inuse_mni=None,
     verbose=None,
 ):
     """Compute LCMV spatial filter.
@@ -733,6 +1192,51 @@ def _make_lcmv(
         How to weight (or normalize) the forward using a depth prior (see Notes).
     inversion : 'matrix' | 'single'
         The inversion scheme to compute the weights.
+    use_bilateral_pairs : bool
+        If True, use bilateral pairs of dipoles in the beamformer computation.
+        If True, then src_pnts_inuse_mni, bilateral_tol and bilateral_tol_mid
+        must also be specified.
+        If dipole_patches_file is also specified, then use_bilateral_pairs is ignored.
+    bilateral_tol : float
+        The distance threshold for finding pairs, in mm                     
+        Recommended to be ~ gridstep / 2
+        Only required if use_bilateral_pairs is True.
+    bilateral_tol_midline : float
+        The distance threshold for finding midline points, in mm
+        Recommended to be ~ gridstep / 2
+        If None then bilateral_tol_midline = bilateral_tol.
+        Only required if use_bilateral_pairs is True.
+    dipole_patches_file : str | None
+        Path to a 4D nifti binary, mutually exclusive parcellation file where
+        each volume should be a mask with 0 for the background
+        and 1 for the parcel, and each voxel should belong to only one parcel.
+        Dipoles within a parcel will be considered as a group
+        of dipoles that will be beamformed together as a multi-dipole.
+        If None, then no patches are used in the beamformer.
+        File must be in MNI space.
+        This only works with volumetric (RHINO) source spaces.
+        Overrides use_bilateral_pairs.
+
+        If group_patches=True is also specified,
+        then the dipole_patches_file needs to be
+        accompanied by a text file with the same name but a .xml
+        extension that contains a "group" field
+        with a unique integer value shared by patches that are 
+        bilateral pairs (if group=0 then the patch 
+        is assumed to not be in a bilateral pair).
+        (see also parcellation.load_parcellation_description).
+        Using a dipole_patches_file only works with volumetric (RHINO)
+        source spaces.
+    group_patches : bool
+        If True, group patches using the "group" field in the dipole patches description (
+        xml) file.
+        If False, do not group patches.
+        Only required if dipole_patches_file is not None.
+    src_pnts_inuse_mni : (ndipoles, 3) numpy.ndarray
+        The 3D coordinates (in mm) of the dipoles in the source space
+        that are being used in the forward model (i.e. where forward['src'][src_index]['inuse']==1).
+        These coordinates must be in MNI space in mm.
+        Required if use_bilateral_pairs is True or dipole_patches_file is not None. 
     verbose : bool, str, int, or None
         If not None, override default verbose level (see mne.verbose).
         
@@ -790,6 +1294,39 @@ def _make_lcmv(
 
     # Compute spatial filter
     n_orient = 3 if is_free_ori else 1
+ 
+    #### MWW 
+    if dipole_patches_file is not None:
+
+        log_or_print(f"Using patches from {dipole_patches_file}")
+        if group_patches:
+            log_or_print("Using group_patches")
+
+        multi_dipoles, single_dipoles = _find_patches(src_pnts_inuse_mni, 
+                                                    dipole_patches_file,
+                                                    group_patches=group_patches,
+                                                    )
+    
+    elif use_bilateral_pairs:
+
+        log_or_print("Using bilateral pairs")
+        log_or_print(f"bilateral_tol={bilateral_tol}mm")
+        log_or_print(f"bilateral_tol_midline={bilateral_tol_midline}mm")
+
+        multi_dipoles, single_dipoles, midline_points = _find_bilateral_pairs(src_pnts_inuse_mni, 
+                                                                  bilateral_tol=bilateral_tol, 
+                                                                  bilateral_tol_midline=bilateral_tol_midline,
+                                                                  )             
+
+        # combine single_dipoles and midline points
+        single_dipoles = np.concatenate([single_dipoles, midline_points])
+
+    else:
+        single_dipoles = None
+        multi_dipoles = None
+
+    #### end MWW
+
     W, max_power_ori = _compute_beamformer(
         G,
         Cm,
@@ -803,6 +1340,8 @@ def _make_lcmv(
         nn=nn,
         orient_std=orient_std,
         whitener=whitener,
+        multi_dipoles=multi_dipoles, # MWW
+        single_dipoles=single_dipoles # MWW
     )
 
     # Get src type to store with filters for _make_stc
@@ -840,7 +1379,12 @@ def _make_lcmv(
     return filters
 
 
-def _compute_beamformer(G, Cm, reg, n_orient, weight_norm, pick_ori, reduce_rank, rank, inversion, nn, orient_std, whitener):
+def _compute_beamformer(G, Cm, reg, n_orient, 
+                        weight_norm, pick_ori, 
+                        reduce_rank, rank, inversion, 
+                        nn, orient_std, whitener,
+                        multi_dipoles=None, single_dipoles=None):
+    
     """Compute a spatial beamformer filter (LCMV or DICS).
 
     For more detailed information on the parameters, see the docstrings of `make_lcmv` and `make_dics`.
@@ -873,6 +1417,12 @@ def _compute_beamformer(G, Cm, reg, n_orient, weight_norm, pick_ori, reduce_rank
         The std of the orientation prior used in weighting the lead fields.
     whitener : (n_channels, n_channels) numpy.ndarray
         The whitener.
+    multi_dipoles : list | None
+        List of length num_multi_dipole, where each item is a list of dipole indices (wrt to G)
+        in each multi-dipole. If None, all dipoles are treated as single-dipoles.
+    single_dipoles : list | None
+        List of dipole indices (wrt to G) that should be treated as
+        single dipoles. If None, all dipoles are treated as multi-dipoles.
 
     Returns
     -------
@@ -907,6 +1457,7 @@ def _compute_beamformer(G, Cm, reg, n_orient, weight_norm, pick_ori, reduce_rank
     assert nn.shape == (n_sources, 3)
 
     mne_logger.info("Computing beamformer filters for %d source%s" % (n_sources, _pl(n_sources)))
+
     n_channels = G.shape[0]
 
     assert n_orient in (3, 1)
@@ -951,7 +1502,10 @@ def _compute_beamformer(G, Cm, reg, n_orient, weight_norm, pick_ori, reduce_rank
 
     def _compute_bf_terms(Gk, Cm_inv):
         bf_numer = np.matmul(Gk.swapaxes(-2, -1).conj(), Cm_inv)
-        bf_denom = np.matmul(bf_numer, Gk)
+
+        # bf_numer is (nsources x norients x nsensors)
+        # Gk is (nsources x nsensors x norients)
+        bf_denom = np.matmul(bf_numer, Gk) # (nsources x nsensors x norient)
         return bf_numer, bf_denom
 
     # ----------------------------------------------------------
@@ -998,7 +1552,7 @@ def _compute_bf_terms(Gk, Cm_inv):
         max_power_ori = eig_vecs[np.arange(len(eig_vecs)), :, order[:, -1]]
 
         assert max_power_ori.shape == (n_sources, n_orient)
-
+        
         # Set the (otherwise arbitrary) sign to match the normal
         signs = np.sign(np.sum(max_power_ori * nn, axis=1, keepdims=True))
         signs[signs == 0] = 1.0
@@ -1014,29 +1568,74 @@ def _compute_bf_terms(Gk, Cm_inv):
             Gk = Gk[..., 2:3]
             n_orient = 1
 
-    # ----------------------------------------------------------------------
-    # 3. Compute numerator and denominator of beamformer formula (unit-gain)
+    # Note Gk is (n_sources, n_channels, n_orient)
+
+    if multi_dipoles is not None:
+
+        W_multi = []
+        for mm, multi_dipole in enumerate(multi_dipoles):
+            # multi_dipole is (ndipoles,)
+            Gk_multi = np.zeros([n_channels, n_orient*len(multi_dipole)])
+            for dd in range(len(multi_dipole)):
+                Gk_multi[:, dd*n_orient:(dd+1)*n_orient] = Gk[multi_dipole[dd], :, :]
+
+            # 3. Compute numerator and denominator of beamformer formula (unit-gain)
+            bf_numer_multi, bf_denom_multi = _compute_bf_terms(Gk_multi, Cm_inv)
+            del Gk_multi
+
+            # 4. Invert the denominator
+            bf_denom_multi_inv = _sym_inv_sm(bf_denom_multi, reduce_rank, inversion='matrix', sk=sk)
+            W_multi.append(np.matmul(bf_denom_multi_inv, bf_numer_multi))
+    
+        if single_dipoles is not None:
+
+            Gk_single = Gk[single_dipoles, :, :]
 
-    bf_numer, bf_denom = _compute_bf_terms(Gk, Cm_inv)
+            # 3. Compute numerator and denominator of beamformer formula (unit-gain)
+            bf_numer_single, bf_denom_single = _compute_bf_terms(Gk_single, Cm_inv)
+            del Gk_single
+
+            # 4. Invert the denominator
+            # W is W_ug, i.e.: G.T @ Cm_inv / (G.T @ Cm_inv @ G)
+            bf_denom_single_inv = _sym_inv_sm(bf_denom_single, reduce_rank, inversion='matrix', sk=sk)
+            W_single = np.matmul(bf_denom_single_inv, bf_numer_single)
+
+        # Build single W from W_multi and W_single
+        W = np.zeros((n_sources, n_orient, n_channels))
+        
+        for mm, multi_dipole in enumerate(multi_dipoles):
+            for dd in range(len(multi_dipole)):
+                W[multi_dipole[dd], :, :] = W_multi[mm][dd*n_orient:(dd+1)*n_orient, :]
+
+        if single_dipoles is not None:
+            W[single_dipoles, :, :] = W_single
+                    
+        assert W.shape == (n_sources, n_orient, n_channels)
+
+    else:
+            
+        # ----------------------------------------------------------------------
+        # 3. Compute numerator and denominator of beamformer formula (unit-gain)
 
-    assert bf_denom.shape == (n_sources,) + (n_orient,) * 2
-    assert bf_numer.shape == (n_sources, n_orient, n_channels)
+        bf_numer, bf_denom = _compute_bf_terms(Gk, Cm_inv)
 
-    del Gk  # lead field has been adjusted and should not be used anymore
+        assert bf_denom.shape == (n_sources,) + (n_orient,) * 2
+        assert bf_numer.shape == (n_sources, n_orient, n_channels)
 
-    # ----------------------------------------------------------------------
-    # 4. Invert the denominator
+        del Gk  # lead field has been adjusted and should not be used anymore
 
-    # Here W is W_ug, i.e.: G.T @ Cm_inv / (G.T @ Cm_inv @ G)
-    bf_denom_inv = _sym_inv_sm(bf_denom, reduce_rank, inversion, sk)
+        # ----------------------------------------------------------------------
+        # 4. Invert the denominator
 
-    assert bf_denom_inv.shape == (n_sources, n_orient, n_orient)
+        # Here W is W_ug, i.e.: G.T @ Cm_inv / (G.T @ Cm_inv @ G)
+        bf_denom_inv = _sym_inv_sm(bf_denom, reduce_rank, inversion='matrix', sk=sk)
 
-    W = np.matmul(bf_denom_inv, bf_numer)
+        assert bf_denom_inv.shape == (n_sources, n_orient, n_orient)
+        W = np.matmul(bf_denom_inv, bf_numer)
 
-    assert W.shape == (n_sources, n_orient, n_channels)
+        assert W.shape == (n_sources, n_orient, n_channels)
 
-    del bf_denom_inv, sk
+        del bf_denom_inv, sk
 
     # ----------------------------------------------------------------
     # 5. Re-scale filter weights according to the selected weight_norm
diff --git a/osl_ephys/source_recon/files/Glasser52_binary_space-MNI152NLin6_res-8x8x8.txt b/osl_ephys/source_recon/files/Glasser52_binary_space-MNI152NLin6_res-8x8x8.txt
deleted file mode 100644
index 5363af3..0000000
--- a/osl_ephys/source_recon/files/Glasser52_binary_space-MNI152NLin6_res-8x8x8.txt
+++ /dev/null
@@ -1 +0,0 @@
-<data>
diff --git a/osl_ephys/source_recon/files/Glasser52_binary_space-MNI152NLin6_res-8x8x8.xml b/osl_ephys/source_recon/files/Glasser52_binary_space-MNI152NLin6_res-8x8x8.xml
index 2c1f242..497a808 100644
--- a/osl_ephys/source_recon/files/Glasser52_binary_space-MNI152NLin6_res-8x8x8.xml
+++ b/osl_ephys/source_recon/files/Glasser52_binary_space-MNI152NLin6_res-8x8x8.xml
@@ -1,54 +1,54 @@
 <data>
-<label index="1" x="14.5" y="-80.5" z="-0.4">Primary and Early Visual Cortex_rh</label>
-<label index="2" x="19.4" y="-81.4" z="30.6">Dorsal Stream Visual Cortex_rh</label>
-<label index="3" x="29.9" y="-57.6" z="-17.7">Ventral Stream Visual Cortex_rh</label>
-<label index="4" x="42.2" y="-71.5" z="0.0">MT+ Complex and Neighboring Visual Areas_rh</label>
-<label index="5" x="23.6" y="-29.5" z="61.0">Superior Somatosnesory and Motor Cortex_rh</label>
-<label index="6" x="49.1" y="-13.4" z="38.3">Inferior Somatosnesory and Motor Cortex_rh</label>
-<label index="7" x="12.7" y="-2.2" z="62.9">Supplementory Motor Area_rh</label>
-<label index="8" x="10.0" y="-31.2" z="53.6">Cingulate Motor Areas & Area 5_rh</label>
-<label index="9" x="40.9" y="0.9" z="41.1">Premotor Cortex_rh</label>
-<label index="10" x="37.7" y="1.6" z="4.0">Insular & Frontoparietal Operculum_rh</label>
-<label index="11" x="40.2" y="-27.5" z="13.5">Early Auditory Cortex_rh</label>
-<label index="12" x="55.0" y="-14.0" z="-6.5">Auditory Association Cortex_rh</label>
-<label index="13" x="25.0" y="-22.1" z="-22.1">Medial Temporal Cortex_rh</label>
-<label index="14" x="48.5" y="-16.3" z="-25.3">Lateral Temporal Cortex_rh</label>
-<label index="15" x="52.1" y="-48.1" z="12.3">Temporal-Parieto-Occipital Junction_rh</label>
-<label index="16" x="29.5" y="-51.0" z="43.0">Medial Bank of the Intra-parietal Sulcus_rh</label>
-<label index="17" x="19.1" y="-58.4" z="59.6">Superior Medial Parietal Cortex_rh</label>
-<label index="18" x="56.7" y="-29.7" z="35.7">Inferior Parietal Cortex Task-Positive Network_rh</label>
-<label index="19" x="47.4" y="-56.1" z="35.8">Inferior Parietal Cortex Task-Negative Network_rh</label>
-<label index="20" x="36.5" y="-71.1" z="32.6">Intraparietal Sulcus & PGP_rh</label>
-<label index="21" x="10.3" y="-55.1" z="26.8">Posterior Cingulate Cortex_rh</label>
-<label index="22" x="6.1" y="33.8" z="14.2">Anterior Cingulate and Medial Prefrontal Cortex_rh</label>
-<label index="23" x="17.1" y="47.0" z="-11.6">Orbital and Polar Frontal Cortex_rh</label>
-<label index="24" x="45.4" y="32.5" z="3.8">Inferior Frontal Cortex_rh</label>
-<label index="25" x="34.8" y="37.1" z="25.7">Inferior Dorso-Lateral Prefrontal Cortex_rh</label>
-<label index="26" x="21.0" y="31.9" z="45.9">Superior Dorso-Lateral Prefronta lCortex_rh</label>
-<label index="27" x="-16.9" y="-82.0" z="-0.9">Primary and Early Visual Cortex_lh</label>
-<label index="28" x="-21.7" y="-84.3" z="27.0">Dorsal Stream Visual Cortex_lh</label>
-<label index="29" x="-34.5" y="-56.0" z="-18.6">Ventral Stream Visual Cortex_lh</label>
-<label index="30" x="-45.2" y="-70.3" z="-1.4">MT+ Complex and Neighboring Visual Areas_lh</label>
-<label index="31" x="-22.6" y="-30.4" z="62.8">Superior Somatosnesory and Motor Cortex_lh</label>
-<label index="32" x="-48.3" y="-18.7" z="41.7">Inferior Somatosnesory and Motor Cortex_lh</label>
-<label index="33" x="-14.0" y="-2.7" z="62.7">Supplementory Motor Area_lh</label>
-<label index="34" x="-14.4" y="-30.6" z="49.7">Cingulate Motor Areas & Area 5_lh</label>
-<label index="35" x="-40.3" y="-1.7" z="43.8">Premotor Cortex_lh</label>
-<label index="36" x="-41.3" y="0.1" z="3.7">Insular & Frontoparietal Operculum_lh</label>
-<label index="37" x="-45.4" y="-29.4" z="12.1">Early Auditory Cortex_lh</label>
-<label index="38" x="-55.3" y="-17.3" z="-7.8">Auditory Association Cortex_lh</label>
-<label index="39" x="-27.9" y="-22.2" z="-22.6">Medial Temporal Cortex_lh</label>
-<label index="40" x="-49.6" y="-17.9" z="-25.1">Lateral Temporal Cortex_lh</label>
-<label index="41" x="-53.1" y="-54.0" z="14.5">Temporal-Parieto-Occipital Junction_lh</label>
-<label index="42" x="-31.1" y="-50.0" z="41.1">Medial Bank of the Intra-parietal Sulcus_lh</label>
-<label index="43" x="-21.3" y="-60.5" z="58.2">Superior Medial Parietal Cortex_lh</label>
-<label index="44" x="-56.6" y="-36.3" z="36.6">Inferior Parietal Cortex Task-Positive Network_lh</label>
-<label index="45" x="-45.1" y="-63.4" z="35.2">Inferior Parietal Cortex Task-Negative Network_lh</label>
-<label index="46" x="-36.7" y="-70.0" z="30.4">Intraparietal Sulcus & PGP_lh</label>
-<label index="47" x="-11.6" y="-49.5" z="28.0">Posterior Cingulate Cortex_lh</label>
-<label index="48" x="-8.7" y="32.0" z="11.4">Anterior Cingulate and Medial Prefrontal Cortex_lh</label>
-<label index="49" x="-21.9" y="46.3" z="-8.1">Orbital and Polar Frontal Cortex_lh</label>
-<label index="50" x="-45.1" y="28.2" z="7.7">Inferior Frontal Cortex_lh</label>
-<label index="51" x="-37.2" y="36.2" z="25.2">Inferior Dorso-Lateral Prefrontal Cortex_lh</label>
-<label index="52" x="-23.5" y="28.6" z="45.7">Superior Dorso-Lateral Prefrontal Cortex_lh</label>
+<label index="1" group="1" x="14.5" y="-80.5" z="-0.4">Primary and Early Visual Cortex_rh</label>
+<label index="2" group="2" x="19.4" y="-81.4" z="30.6">Dorsal Stream Visual Cortex_rh</label>
+<label index="3" group="3" x="29.9" y="-57.6" z="-17.7">Ventral Stream Visual Cortex_rh</label>
+<label index="4" group="4" x="42.2" y="-71.5" z="0.0">MT+ Complex and Neighboring Visual Areas_rh</label>
+<label index="5" group="5" x="23.6" y="-29.5" z="61.0">Superior Somatosnesory and Motor Cortex_rh</label>
+<label index="6" group="6" x="49.1" y="-13.4" z="38.3">Inferior Somatosnesory and Motor Cortex_rh</label>
+<label index="7" group="7" x="12.7" y="-2.2" z="62.9">Supplementory Motor Area_rh</label>
+<label index="8" group="8" x="10.0" y="-31.2" z="53.6">Cingulate Motor Areas and Area 5_rh</label>
+<label index="9" group="9" x="40.9" y="0.9" z="41.1">Premotor Cortex_rh</label>
+<label index="10" group="10" x="37.7" y="1.6" z="4.0">Insular and Frontoparietal Operculum_rh</label>
+<label index="11" group="11" x="40.2" y="-27.5" z="13.5">Early Auditory Cortex_rh</label>
+<label index="12" group="12" x="55.0" y="-14.0" z="-6.5">Auditory Association Cortex_rh</label>
+<label index="13" group="13" x="25.0" y="-22.1" z="-22.1">Medial Temporal Cortex_rh</label>
+<label index="14" group="14" x="48.5" y="-16.3" z="-25.3">Lateral Temporal Cortex_rh</label>
+<label index="15" group="15" x="52.1" y="-48.1" z="12.3">Temporal-Parieto-Occipital Junction_rh</label>
+<label index="16" group="16" x="29.5" y="-51.0" z="43.0">Medial Bank of the Intra-parietal Sulcus_rh</label>
+<label index="17" group="17" x="19.1" y="-58.4" z="59.6">Superior Medial Parietal Cortex_rh</label>
+<label index="18" group="18" x="56.7" y="-29.7" z="35.7">Inferior Parietal Cortex Task-Positive Network_rh</label>
+<label index="19" group="19" x="47.4" y="-56.1" z="35.8">Inferior Parietal Cortex Task-Negative Network_rh</label>
+<label index="20" group="20" x="36.5" y="-71.1" z="32.6">Intraparietal Sulcus and PGP_rh</label>
+<label index="21" group="21" x="10.3" y="-55.1" z="26.8">Posterior Cingulate Cortex_rh</label>
+<label index="22" group="22" x="6.1" y="33.8" z="14.2">Anterior Cingulate and Medial Prefrontal Cortex_rh</label>
+<label index="23" group="23" x="17.1" y="47.0" z="-11.6">Orbital and Polar Frontal Cortex_rh</label>
+<label index="24" group="24" x="45.4" y="32.5" z="3.8">Inferior Frontal Cortex_rh</label>
+<label index="25" group="25" x="34.8" y="37.1" z="25.7">Inferior Dorso-Lateral Prefrontal Cortex_rh</label>
+<label index="26" group="26" x="21.0" y="31.9" z="45.9">Superior Dorso-Lateral Prefronta lCortex_rh</label>
+<label index="27" group="1" x="-16.9" y="-82.0" z="-0.9">Primary and Early Visual Cortex_lh</label>
+<label index="28" group="2" x="-21.7" y="-84.3" z="27.0">Dorsal Stream Visual Cortex_lh</label>
+<label index="29" group="3" x="-34.5" y="-56.0" z="-18.6">Ventral Stream Visual Cortex_lh</label>
+<label index="30" group="4" x="-45.2" y="-70.3" z="-1.4">MT+ Complex and Neighboring Visual Areas_lh</label>
+<label index="31" group="5" x="-22.6" y="-30.4" z="62.8">Superior Somatosnesory and Motor Cortex_lh</label>
+<label index="32" group="6" x="-48.3" y="-18.7" z="41.7">Inferior Somatosnesory and Motor Cortex_lh</label>
+<label index="33" group="7" x="-14.0" y="-2.7" z="62.7">Supplementory Motor Area_lh</label>
+<label index="34" group="8" x="-14.4" y="-30.6" z="49.7">Cingulate Motor Areas and Area 5_lh</label>
+<label index="35" group="9" x="-40.3" y="-1.7" z="43.8">Premotor Cortex_lh</label>
+<label index="36" group="10" x="-41.3" y="0.1" z="3.7">Insular and Frontoparietal Operculum_lh</label>
+<label index="37" group="11" x="-45.4" y="-29.4" z="12.1">Early Auditory Cortex_lh</label>
+<label index="38" group="12" x="-55.3" y="-17.3" z="-7.8">Auditory Association Cortex_lh</label>
+<label index="39" group="13" x="-27.9" y="-22.2" z="-22.6">Medial Temporal Cortex_lh</label>
+<label index="40" group="14" x="-49.6" y="-17.9" z="-25.1">Lateral Temporal Cortex_lh</label>
+<label index="41" group="15" x="-53.1" y="-54.0" z="14.5">Temporal-Parieto-Occipital Junction_lh</label>
+<label index="42" group="16" x="-31.1" y="-50.0" z="41.1">Medial Bank of the Intra-parietal Sulcus_lh</label>
+<label index="43" group="17" x="-21.3" y="-60.5" z="58.2">Superior Medial Parietal Cortex_lh</label>
+<label index="44" group="18" x="-56.6" y="-36.3" z="36.6">Inferior Parietal Cortex Task-Positive Network_lh</label>
+<label index="45" group="19" x="-45.1" y="-63.4" z="35.2">Inferior Parietal Cortex Task-Negative Network_lh</label>
+<label index="46" group="20" x="-36.7" y="-70.0" z="30.4">Intraparietal Sulcus and PGP_lh</label>
+<label index="47" group="21" x="-11.6" y="-49.5" z="28.0">Posterior Cingulate Cortex_lh</label>
+<label index="48" group="22" x="-8.7" y="32.0" z="11.4">Anterior Cingulate and Medial Prefrontal Cortex_lh</label>
+<label index="49" group="23" x="-21.9" y="46.3" z="-8.1">Orbital and Polar Frontal Cortex_lh</label>
+<label index="50" group="24" x="-45.1" y="28.2" z="7.7">Inferior Frontal Cortex_lh</label>
+<label index="51" group="25" x="-37.2" y="36.2" z="25.2">Inferior Dorso-Lateral Prefrontal Cortex_lh</label>
+<label index="52" group="26" x="-23.5" y="28.6" z="45.7">Superior Dorso-Lateral Prefrontal Cortex_lh</label>
 </data>
diff --git a/osl_ephys/source_recon/files/Glasser52_binary_space-MNI152NLin6_res-8x8x8_10mm.nii.gz b/osl_ephys/source_recon/files/Glasser52_binary_space-MNI152NLin6_res-8x8x8_10mm.nii.gz
new file mode 100644
index 0000000..a9cfb83
Binary files /dev/null and b/osl_ephys/source_recon/files/Glasser52_binary_space-MNI152NLin6_res-8x8x8_10mm.nii.gz differ
diff --git a/osl_ephys/source_recon/files/HarvOxf-sub-Schaefer100-combined-2mm_4d.xml b/osl_ephys/source_recon/files/HarvOxf-sub-Schaefer100-combined-2mm_4d.xml
index 4ba08ce..2b32c4b 100644
--- a/osl_ephys/source_recon/files/HarvOxf-sub-Schaefer100-combined-2mm_4d.xml
+++ b/osl_ephys/source_recon/files/HarvOxf-sub-Schaefer100-combined-2mm_4d.xml
@@ -1,116 +1,116 @@
 <data>
-<label index="1">7Networks_LH_Vis_1</label>
-<label index="2">7Networks_LH_Vis_2</label>
-<label index="3">7Networks_LH_Vis_3</label>
-<label index="4">7Networks_LH_Vis_4</label>
-<label index="5">7Networks_LH_Vis_5</label>
-<label index="6">7Networks_LH_Vis_6</label>
-<label index="7">7Networks_LH_Vis_7</label>
-<label index="8">7Networks_LH_Vis_8</label>
-<label index="9">7Networks_LH_Vis_9</label>
-<label index="10">7Networks_LH_SomMot_1</label>
-<label index="11">7Networks_LH_SomMot_2</label>
-<label index="12">7Networks_LH_SomMot_3</label>
-<label index="13">7Networks_LH_SomMot_4</label>
-<label index="14">7Networks_LH_SomMot_5</label>
-<label index="15">7Networks_LH_SomMot_6</label>
-<label index="16">7Networks_LH_DorsAttn_Post_1</label>
-<label index="17">7Networks_LH_DorsAttn_Post_2</label>
-<label index="18">7Networks_LH_DorsAttn_Post_3</label>
-<label index="19">7Networks_LH_DorsAttn_Post_4</label>
-<label index="20">7Networks_LH_DorsAttn_Post_5</label>
-<label index="21">7Networks_LH_DorsAttn_Post_6</label>
-<label index="22">7Networks_LH_DorsAttn_PrCv_1</label>
-<label index="23">7Networks_LH_DorsAttn_FEF_1</label>
-<label index="24">7Networks_LH_SalVentAttn_ParOper_1</label>
-<label index="25">7Networks_LH_SalVentAttn_FrOperIns_1</label>
-<label index="26">7Networks_LH_SalVentAttn_FrOperIns_2</label>
-<label index="27">7Networks_LH_SalVentAttn_PFCl_1</label>
-<label index="28">7Networks_LH_SalVentAttn_Med_1</label>
-<label index="29">7Networks_LH_SalVentAttn_Med_2</label>
-<label index="30">7Networks_LH_SalVentAttn_Med_3</label>
-<label index="31">7Networks_LH_Limbic_OFC_1</label>
-<label index="32">7Networks_LH_Limbic_TempPole_1</label>
-<label index="33">7Networks_LH_Limbic_TempPole_2</label>
-<label index="34">7Networks_LH_Cont_Par_1</label>
-<label index="35">7Networks_LH_Cont_PFCl_1</label>
-<label index="36">7Networks_LH_Cont_pCun_1</label>
-<label index="37">7Networks_LH_Cont_Cing_1</label>
-<label index="38">7Networks_LH_Default_Temp_1</label>
-<label index="39">7Networks_LH_Default_Temp_2</label>
-<label index="40">7Networks_LH_Default_Par_1</label>
-<label index="41">7Networks_LH_Default_Par_2</label>
-<label index="42">7Networks_LH_Default_PFC_1</label>
-<label index="43">7Networks_LH_Default_PFC_2</label>
-<label index="44">7Networks_LH_Default_PFC_3</label>
-<label index="45">7Networks_LH_Default_PFC_4</label>
-<label index="46">7Networks_LH_Default_PFC_5</label>
-<label index="47">7Networks_LH_Default_PFC_6</label>
-<label index="48">7Networks_LH_Default_PFC_7</label>
-<label index="49">7Networks_LH_Default_pCunPCC_1</label>
-<label index="50">7Networks_LH_Default_pCunPCC_2</label>
-<label index="51">7Networks_RH_Vis_1</label>
-<label index="52">7Networks_RH_Vis_2</label>
-<label index="53">7Networks_RH_Vis_3</label>
-<label index="54">7Networks_RH_Vis_4</label>
-<label index="55">7Networks_RH_Vis_5</label>
-<label index="56">7Networks_RH_Vis_6</label>
-<label index="57">7Networks_RH_Vis_7</label>
-<label index="58">7Networks_RH_Vis_8</label>
-<label index="59">7Networks_RH_SomMot_1</label>
-<label index="60">7Networks_RH_SomMot_2</label>
-<label index="61">7Networks_RH_SomMot_3</label>
-<label index="62">7Networks_RH_SomMot_4</label>
-<label index="63">7Networks_RH_SomMot_5</label>
-<label index="64">7Networks_RH_SomMot_6</label>
-<label index="65">7Networks_RH_SomMot_7</label>
-<label index="66">7Networks_RH_SomMot_8</label>
-<label index="67">7Networks_RH_DorsAttn_Post_1</label>
-<label index="68">7Networks_RH_DorsAttn_Post_2</label>
-<label index="69">7Networks_RH_DorsAttn_Post_3</label>
-<label index="70">7Networks_RH_DorsAttn_Post_4</label>
-<label index="71">7Networks_RH_DorsAttn_Post_5</label>
-<label index="72">7Networks_RH_DorsAttn_PrCv_1</label>
-<label index="73">7Networks_RH_DorsAttn_FEF_1</label>
-<label index="74">7Networks_RH_SalVentAttn_TempOccPar_1</label>
-<label index="75">7Networks_RH_SalVentAttn_TempOccPar_2</label>
-<label index="76">7Networks_RH_SalVentAttn_FrOperIns_1</label>
-<label index="77">7Networks_RH_SalVentAttn_Med_1</label>
-<label index="78">7Networks_RH_SalVentAttn_Med_2</label>
-<label index="79">7Networks_RH_Limbic_OFC_1</label>
-<label index="80">7Networks_RH_Limbic_TempPole_1</label>
-<label index="81">7Networks_RH_Cont_Par_1</label>
-<label index="82">7Networks_RH_Cont_Par_2</label>
-<label index="83">7Networks_RH_Cont_PFCl_1</label>
-<label index="84">7Networks_RH_Cont_PFCl_2</label>
-<label index="85">7Networks_RH_Cont_PFCl_3</label>
-<label index="86">7Networks_RH_Cont_PFCl_4</label>
-<label index="87">7Networks_RH_Cont_Cing_1</label>
-<label index="88">7Networks_RH_Cont_PFCmp_1</label>
-<label index="89">7Networks_RH_Cont_pCun_1</label>
-<label index="90">7Networks_RH_Default_Par_1</label>
-<label index="91">7Networks_RH_Default_Temp_1</label>
-<label index="92">7Networks_RH_Default_Temp_2</label>
-<label index="93">7Networks_RH_Default_Temp_3</label>
-<label index="94">7Networks_RH_Default_PFCv_1</label>
-<label index="95">7Networks_RH_Default_PFCv_2</label>
-<label index="96">7Networks_RH_Default_PFCdPFCm_1</label>
-<label index="97">7Networks_RH_Default_PFCdPFCm_2</label>
-<label index="98">7Networks_RH_Default_PFCdPFCm_3</label>
-<label index="99">7Networks_RH_Default_pCunPCC_1</label>
-<label index="100">7Networks_RH_Default_pCunPCC_2</label>
-<label index="101">LeftThalamus</label>
-<label index="102">LeftCaudate</label>
-<label index="103">LeftPutamen</label>
-<label index="104">LeftPallidum</label>
-<label index="105">LeftHippocampus</label>
-<label index="106">LeftAmygdala</label>
-<label index="107">LeftAccumbens</label>
-<label index="108">RightThalamus</label>
-<label index="109">RightCaudate</label>
-<label index="110">RightPutamen</label>
-<label index="111">RightPallidum</label>
-<label index="112">RightHippocampus</label>
-<label index="113">RightAmygdala</label>
-<label index="114">RightAccumbens</label>
+<label group="0" index="1">7Networks_LH_Vis_1</label>
+<label group="0" index="2">7Networks_LH_Vis_2</label>
+<label group="0" index="3">7Networks_LH_Vis_3</label>
+<label group="0" index="4">7Networks_LH_Vis_4</label>
+<label group="0" index="5">7Networks_LH_Vis_5</label>
+<label group="0" index="6">7Networks_LH_Vis_6</label>
+<label group="0" index="7">7Networks_LH_Vis_7</label>
+<label group="0" index="8">7Networks_LH_Vis_8</label>
+<label group="0" index="9">7Networks_LH_Vis_9</label>
+<label group="0" index="10">7Networks_LH_SomMot_1</label>
+<label group="0" index="11">7Networks_LH_SomMot_2</label>
+<label group="0" index="12">7Networks_LH_SomMot_3</label>
+<label group="0" index="13">7Networks_LH_SomMot_4</label>
+<label group="0" index="14">7Networks_LH_SomMot_5</label>
+<label group="0" index="15">7Networks_LH_SomMot_6</label>
+<label group="0" index="16">7Networks_LH_DorsAttn_Post_1</label>
+<label group="0" index="17">7Networks_LH_DorsAttn_Post_2</label>
+<label group="0" index="18">7Networks_LH_DorsAttn_Post_3</label>
+<label group="0" index="19">7Networks_LH_DorsAttn_Post_4</label>
+<label group="0" index="20">7Networks_LH_DorsAttn_Post_5</label>
+<label group="0" index="21">7Networks_LH_DorsAttn_Post_6</label>
+<label group="0" index="22">7Networks_LH_DorsAttn_PrCv_1</label>
+<label group="0" index="23">7Networks_LH_DorsAttn_FEF_1</label>
+<label group="0" index="24">7Networks_LH_SalVentAttn_ParOper_1</label>
+<label group="0" index="25">7Networks_LH_SalVentAttn_FrOperIns_1</label>
+<label group="0" index="26">7Networks_LH_SalVentAttn_FrOperIns_2</label>
+<label group="0" index="27">7Networks_LH_SalVentAttn_PFCl_1</label>
+<label group="0" index="28">7Networks_LH_SalVentAttn_Med_1</label>
+<label group="0" index="29">7Networks_LH_SalVentAttn_Med_2</label>
+<label group="0" index="30">7Networks_LH_SalVentAttn_Med_3</label>
+<label group="0" index="31">7Networks_LH_Limbic_OFC_1</label>
+<label group="0" index="32">7Networks_LH_Limbic_TempPole_1</label>
+<label group="0" index="33">7Networks_LH_Limbic_TempPole_2</label>
+<label group="0" index="34">7Networks_LH_Cont_Par_1</label>
+<label group="0" index="35">7Networks_LH_Cont_PFCl_1</label>
+<label group="0" index="36">7Networks_LH_Cont_pCun_1</label>
+<label group="0" index="37">7Networks_LH_Cont_Cing_1</label>
+<label group="0" index="38">7Networks_LH_Default_Temp_1</label>
+<label group="0" index="39">7Networks_LH_Default_Temp_2</label>
+<label group="0" index="40">7Networks_LH_Default_Par_1</label>
+<label group="0" index="41">7Networks_LH_Default_Par_2</label>
+<label group="0" index="42">7Networks_LH_Default_PFC_1</label>
+<label group="0" index="43">7Networks_LH_Default_PFC_2</label>
+<label group="0" index="44">7Networks_LH_Default_PFC_3</label>
+<label group="0" index="45">7Networks_LH_Default_PFC_4</label>
+<label group="0" index="46">7Networks_LH_Default_PFC_5</label>
+<label group="0" index="47">7Networks_LH_Default_PFC_6</label>
+<label group="0" index="48">7Networks_LH_Default_PFC_7</label>
+<label group="0" index="49">7Networks_LH_Default_pCunPCC_1</label>
+<label group="0" index="50">7Networks_LH_Default_pCunPCC_2</label>
+<label group="0" index="51">7Networks_RH_Vis_1</label>
+<label group="0" index="52">7Networks_RH_Vis_2</label>
+<label group="0" index="53">7Networks_RH_Vis_3</label>
+<label group="0" index="54">7Networks_RH_Vis_4</label>
+<label group="0" index="55">7Networks_RH_Vis_5</label>
+<label group="0" index="56">7Networks_RH_Vis_6</label>
+<label group="0" index="57">7Networks_RH_Vis_7</label>
+<label group="0" index="58">7Networks_RH_Vis_8</label>
+<label group="0" index="59">7Networks_RH_SomMot_1</label>
+<label group="0" index="60">7Networks_RH_SomMot_2</label>
+<label group="0" index="61">7Networks_RH_SomMot_3</label>
+<label group="0" index="62">7Networks_RH_SomMot_4</label>
+<label group="0" index="63">7Networks_RH_SomMot_5</label>
+<label group="0" index="64">7Networks_RH_SomMot_6</label>
+<label group="0" index="65">7Networks_RH_SomMot_7</label>
+<label group="0" index="66">7Networks_RH_SomMot_8</label>
+<label group="0" index="67">7Networks_RH_DorsAttn_Post_1</label>
+<label group="0" index="68">7Networks_RH_DorsAttn_Post_2</label>
+<label group="0" index="69">7Networks_RH_DorsAttn_Post_3</label>
+<label group="0" index="70">7Networks_RH_DorsAttn_Post_4</label>
+<label group="0" index="71">7Networks_RH_DorsAttn_Post_5</label>
+<label group="0" index="72">7Networks_RH_DorsAttn_PrCv_1</label>
+<label group="0" index="73">7Networks_RH_DorsAttn_FEF_1</label>
+<label group="0" index="74">7Networks_RH_SalVentAttn_TempOccPar_1</label>
+<label group="0" index="75">7Networks_RH_SalVentAttn_TempOccPar_2</label>
+<label group="0" index="76">7Networks_RH_SalVentAttn_FrOperIns_1</label>
+<label group="0" index="77">7Networks_RH_SalVentAttn_Med_1</label>
+<label group="0" index="78">7Networks_RH_SalVentAttn_Med_2</label>
+<label group="0" index="79">7Networks_RH_Limbic_OFC_1</label>
+<label group="0" index="80">7Networks_RH_Limbic_TempPole_1</label>
+<label group="0" index="81">7Networks_RH_Cont_Par_1</label>
+<label group="0" index="82">7Networks_RH_Cont_Par_2</label>
+<label group="0" index="83">7Networks_RH_Cont_PFCl_1</label>
+<label group="0" index="84">7Networks_RH_Cont_PFCl_2</label>
+<label group="0" index="85">7Networks_RH_Cont_PFCl_3</label>
+<label group="0" index="86">7Networks_RH_Cont_PFCl_4</label>
+<label group="0" index="87">7Networks_RH_Cont_Cing_1</label>
+<label group="0" index="88">7Networks_RH_Cont_PFCmp_1</label>
+<label group="0" index="89">7Networks_RH_Cont_pCun_1</label>
+<label group="0" index="90">7Networks_RH_Default_Par_1</label>
+<label group="0" index="91">7Networks_RH_Default_Temp_1</label>
+<label group="0" index="92">7Networks_RH_Default_Temp_2</label>
+<label group="0" index="93">7Networks_RH_Default_Temp_3</label>
+<label group="0" index="94">7Networks_RH_Default_PFCv_1</label>
+<label group="0" index="95">7Networks_RH_Default_PFCv_2</label>
+<label group="0" index="96">7Networks_RH_Default_PFCdPFCm_1</label>
+<label group="0" index="97">7Networks_RH_Default_PFCdPFCm_2</label>
+<label group="0" index="98">7Networks_RH_Default_PFCdPFCm_3</label>
+<label group="0" index="99">7Networks_RH_Default_pCunPCC_1</label>
+<label group="0" index="100">7Networks_RH_Default_pCunPCC_2</label>
+<label group="0" index="101">LeftThalamus</label>
+<label group="0" index="102">LeftCaudate</label>
+<label group="0" index="103">LeftPutamen</label>
+<label group="0" index="104">LeftPallidum</label>
+<label group="0" index="105">LeftHippocampus</label>
+<label group="0" index="106">LeftAmygdala</label>
+<label group="0" index="107">LeftAccumbens</label>
+<label group="0" index="108">RightThalamus</label>
+<label group="0" index="109">RightCaudate</label>
+<label group="0" index="110">RightPutamen</label>
+<label group="0" index="111">RightPallidum</label>
+<label group="0" index="112">RightHippocampus</label>
+<label group="0" index="113">RightAmygdala</label>
+<label group="0" index="114">RightAccumbens</label>
 </data>
diff --git a/osl_ephys/source_recon/files/HarvOxf-sub-Schaefer100-combined-2mm_4d_ds8.xml b/osl_ephys/source_recon/files/HarvOxf-sub-Schaefer100-combined-2mm_4d_ds8.xml
index 4ba08ce..2b32c4b 100644
--- a/osl_ephys/source_recon/files/HarvOxf-sub-Schaefer100-combined-2mm_4d_ds8.xml
+++ b/osl_ephys/source_recon/files/HarvOxf-sub-Schaefer100-combined-2mm_4d_ds8.xml
@@ -1,116 +1,116 @@
 <data>
-<label index="1">7Networks_LH_Vis_1</label>
-<label index="2">7Networks_LH_Vis_2</label>
-<label index="3">7Networks_LH_Vis_3</label>
-<label index="4">7Networks_LH_Vis_4</label>
-<label index="5">7Networks_LH_Vis_5</label>
-<label index="6">7Networks_LH_Vis_6</label>
-<label index="7">7Networks_LH_Vis_7</label>
-<label index="8">7Networks_LH_Vis_8</label>
-<label index="9">7Networks_LH_Vis_9</label>
-<label index="10">7Networks_LH_SomMot_1</label>
-<label index="11">7Networks_LH_SomMot_2</label>
-<label index="12">7Networks_LH_SomMot_3</label>
-<label index="13">7Networks_LH_SomMot_4</label>
-<label index="14">7Networks_LH_SomMot_5</label>
-<label index="15">7Networks_LH_SomMot_6</label>
-<label index="16">7Networks_LH_DorsAttn_Post_1</label>
-<label index="17">7Networks_LH_DorsAttn_Post_2</label>
-<label index="18">7Networks_LH_DorsAttn_Post_3</label>
-<label index="19">7Networks_LH_DorsAttn_Post_4</label>
-<label index="20">7Networks_LH_DorsAttn_Post_5</label>
-<label index="21">7Networks_LH_DorsAttn_Post_6</label>
-<label index="22">7Networks_LH_DorsAttn_PrCv_1</label>
-<label index="23">7Networks_LH_DorsAttn_FEF_1</label>
-<label index="24">7Networks_LH_SalVentAttn_ParOper_1</label>
-<label index="25">7Networks_LH_SalVentAttn_FrOperIns_1</label>
-<label index="26">7Networks_LH_SalVentAttn_FrOperIns_2</label>
-<label index="27">7Networks_LH_SalVentAttn_PFCl_1</label>
-<label index="28">7Networks_LH_SalVentAttn_Med_1</label>
-<label index="29">7Networks_LH_SalVentAttn_Med_2</label>
-<label index="30">7Networks_LH_SalVentAttn_Med_3</label>
-<label index="31">7Networks_LH_Limbic_OFC_1</label>
-<label index="32">7Networks_LH_Limbic_TempPole_1</label>
-<label index="33">7Networks_LH_Limbic_TempPole_2</label>
-<label index="34">7Networks_LH_Cont_Par_1</label>
-<label index="35">7Networks_LH_Cont_PFCl_1</label>
-<label index="36">7Networks_LH_Cont_pCun_1</label>
-<label index="37">7Networks_LH_Cont_Cing_1</label>
-<label index="38">7Networks_LH_Default_Temp_1</label>
-<label index="39">7Networks_LH_Default_Temp_2</label>
-<label index="40">7Networks_LH_Default_Par_1</label>
-<label index="41">7Networks_LH_Default_Par_2</label>
-<label index="42">7Networks_LH_Default_PFC_1</label>
-<label index="43">7Networks_LH_Default_PFC_2</label>
-<label index="44">7Networks_LH_Default_PFC_3</label>
-<label index="45">7Networks_LH_Default_PFC_4</label>
-<label index="46">7Networks_LH_Default_PFC_5</label>
-<label index="47">7Networks_LH_Default_PFC_6</label>
-<label index="48">7Networks_LH_Default_PFC_7</label>
-<label index="49">7Networks_LH_Default_pCunPCC_1</label>
-<label index="50">7Networks_LH_Default_pCunPCC_2</label>
-<label index="51">7Networks_RH_Vis_1</label>
-<label index="52">7Networks_RH_Vis_2</label>
-<label index="53">7Networks_RH_Vis_3</label>
-<label index="54">7Networks_RH_Vis_4</label>
-<label index="55">7Networks_RH_Vis_5</label>
-<label index="56">7Networks_RH_Vis_6</label>
-<label index="57">7Networks_RH_Vis_7</label>
-<label index="58">7Networks_RH_Vis_8</label>
-<label index="59">7Networks_RH_SomMot_1</label>
-<label index="60">7Networks_RH_SomMot_2</label>
-<label index="61">7Networks_RH_SomMot_3</label>
-<label index="62">7Networks_RH_SomMot_4</label>
-<label index="63">7Networks_RH_SomMot_5</label>
-<label index="64">7Networks_RH_SomMot_6</label>
-<label index="65">7Networks_RH_SomMot_7</label>
-<label index="66">7Networks_RH_SomMot_8</label>
-<label index="67">7Networks_RH_DorsAttn_Post_1</label>
-<label index="68">7Networks_RH_DorsAttn_Post_2</label>
-<label index="69">7Networks_RH_DorsAttn_Post_3</label>
-<label index="70">7Networks_RH_DorsAttn_Post_4</label>
-<label index="71">7Networks_RH_DorsAttn_Post_5</label>
-<label index="72">7Networks_RH_DorsAttn_PrCv_1</label>
-<label index="73">7Networks_RH_DorsAttn_FEF_1</label>
-<label index="74">7Networks_RH_SalVentAttn_TempOccPar_1</label>
-<label index="75">7Networks_RH_SalVentAttn_TempOccPar_2</label>
-<label index="76">7Networks_RH_SalVentAttn_FrOperIns_1</label>
-<label index="77">7Networks_RH_SalVentAttn_Med_1</label>
-<label index="78">7Networks_RH_SalVentAttn_Med_2</label>
-<label index="79">7Networks_RH_Limbic_OFC_1</label>
-<label index="80">7Networks_RH_Limbic_TempPole_1</label>
-<label index="81">7Networks_RH_Cont_Par_1</label>
-<label index="82">7Networks_RH_Cont_Par_2</label>
-<label index="83">7Networks_RH_Cont_PFCl_1</label>
-<label index="84">7Networks_RH_Cont_PFCl_2</label>
-<label index="85">7Networks_RH_Cont_PFCl_3</label>
-<label index="86">7Networks_RH_Cont_PFCl_4</label>
-<label index="87">7Networks_RH_Cont_Cing_1</label>
-<label index="88">7Networks_RH_Cont_PFCmp_1</label>
-<label index="89">7Networks_RH_Cont_pCun_1</label>
-<label index="90">7Networks_RH_Default_Par_1</label>
-<label index="91">7Networks_RH_Default_Temp_1</label>
-<label index="92">7Networks_RH_Default_Temp_2</label>
-<label index="93">7Networks_RH_Default_Temp_3</label>
-<label index="94">7Networks_RH_Default_PFCv_1</label>
-<label index="95">7Networks_RH_Default_PFCv_2</label>
-<label index="96">7Networks_RH_Default_PFCdPFCm_1</label>
-<label index="97">7Networks_RH_Default_PFCdPFCm_2</label>
-<label index="98">7Networks_RH_Default_PFCdPFCm_3</label>
-<label index="99">7Networks_RH_Default_pCunPCC_1</label>
-<label index="100">7Networks_RH_Default_pCunPCC_2</label>
-<label index="101">LeftThalamus</label>
-<label index="102">LeftCaudate</label>
-<label index="103">LeftPutamen</label>
-<label index="104">LeftPallidum</label>
-<label index="105">LeftHippocampus</label>
-<label index="106">LeftAmygdala</label>
-<label index="107">LeftAccumbens</label>
-<label index="108">RightThalamus</label>
-<label index="109">RightCaudate</label>
-<label index="110">RightPutamen</label>
-<label index="111">RightPallidum</label>
-<label index="112">RightHippocampus</label>
-<label index="113">RightAmygdala</label>
-<label index="114">RightAccumbens</label>
+<label group="0" index="1">7Networks_LH_Vis_1</label>
+<label group="0" index="2">7Networks_LH_Vis_2</label>
+<label group="0" index="3">7Networks_LH_Vis_3</label>
+<label group="0" index="4">7Networks_LH_Vis_4</label>
+<label group="0" index="5">7Networks_LH_Vis_5</label>
+<label group="0" index="6">7Networks_LH_Vis_6</label>
+<label group="0" index="7">7Networks_LH_Vis_7</label>
+<label group="0" index="8">7Networks_LH_Vis_8</label>
+<label group="0" index="9">7Networks_LH_Vis_9</label>
+<label group="0" index="10">7Networks_LH_SomMot_1</label>
+<label group="0" index="11">7Networks_LH_SomMot_2</label>
+<label group="0" index="12">7Networks_LH_SomMot_3</label>
+<label group="0" index="13">7Networks_LH_SomMot_4</label>
+<label group="0" index="14">7Networks_LH_SomMot_5</label>
+<label group="0" index="15">7Networks_LH_SomMot_6</label>
+<label group="0" index="16">7Networks_LH_DorsAttn_Post_1</label>
+<label group="0" index="17">7Networks_LH_DorsAttn_Post_2</label>
+<label group="0" index="18">7Networks_LH_DorsAttn_Post_3</label>
+<label group="0" index="19">7Networks_LH_DorsAttn_Post_4</label>
+<label group="0" index="20">7Networks_LH_DorsAttn_Post_5</label>
+<label group="0" index="21">7Networks_LH_DorsAttn_Post_6</label>
+<label group="0" index="22">7Networks_LH_DorsAttn_PrCv_1</label>
+<label group="0" index="23">7Networks_LH_DorsAttn_FEF_1</label>
+<label group="0" index="24">7Networks_LH_SalVentAttn_ParOper_1</label>
+<label group="0" index="25">7Networks_LH_SalVentAttn_FrOperIns_1</label>
+<label group="0" index="26">7Networks_LH_SalVentAttn_FrOperIns_2</label>
+<label group="0" index="27">7Networks_LH_SalVentAttn_PFCl_1</label>
+<label group="0" index="28">7Networks_LH_SalVentAttn_Med_1</label>
+<label group="0" index="29">7Networks_LH_SalVentAttn_Med_2</label>
+<label group="0" index="30">7Networks_LH_SalVentAttn_Med_3</label>
+<label group="0" index="31">7Networks_LH_Limbic_OFC_1</label>
+<label group="0" index="32">7Networks_LH_Limbic_TempPole_1</label>
+<label group="0" index="33">7Networks_LH_Limbic_TempPole_2</label>
+<label group="0" index="34">7Networks_LH_Cont_Par_1</label>
+<label group="0" index="35">7Networks_LH_Cont_PFCl_1</label>
+<label group="0" index="36">7Networks_LH_Cont_pCun_1</label>
+<label group="0" index="37">7Networks_LH_Cont_Cing_1</label>
+<label group="0" index="38">7Networks_LH_Default_Temp_1</label>
+<label group="0" index="39">7Networks_LH_Default_Temp_2</label>
+<label group="0" index="40">7Networks_LH_Default_Par_1</label>
+<label group="0" index="41">7Networks_LH_Default_Par_2</label>
+<label group="0" index="42">7Networks_LH_Default_PFC_1</label>
+<label group="0" index="43">7Networks_LH_Default_PFC_2</label>
+<label group="0" index="44">7Networks_LH_Default_PFC_3</label>
+<label group="0" index="45">7Networks_LH_Default_PFC_4</label>
+<label group="0" index="46">7Networks_LH_Default_PFC_5</label>
+<label group="0" index="47">7Networks_LH_Default_PFC_6</label>
+<label group="0" index="48">7Networks_LH_Default_PFC_7</label>
+<label group="0" index="49">7Networks_LH_Default_pCunPCC_1</label>
+<label group="0" index="50">7Networks_LH_Default_pCunPCC_2</label>
+<label group="0" index="51">7Networks_RH_Vis_1</label>
+<label group="0" index="52">7Networks_RH_Vis_2</label>
+<label group="0" index="53">7Networks_RH_Vis_3</label>
+<label group="0" index="54">7Networks_RH_Vis_4</label>
+<label group="0" index="55">7Networks_RH_Vis_5</label>
+<label group="0" index="56">7Networks_RH_Vis_6</label>
+<label group="0" index="57">7Networks_RH_Vis_7</label>
+<label group="0" index="58">7Networks_RH_Vis_8</label>
+<label group="0" index="59">7Networks_RH_SomMot_1</label>
+<label group="0" index="60">7Networks_RH_SomMot_2</label>
+<label group="0" index="61">7Networks_RH_SomMot_3</label>
+<label group="0" index="62">7Networks_RH_SomMot_4</label>
+<label group="0" index="63">7Networks_RH_SomMot_5</label>
+<label group="0" index="64">7Networks_RH_SomMot_6</label>
+<label group="0" index="65">7Networks_RH_SomMot_7</label>
+<label group="0" index="66">7Networks_RH_SomMot_8</label>
+<label group="0" index="67">7Networks_RH_DorsAttn_Post_1</label>
+<label group="0" index="68">7Networks_RH_DorsAttn_Post_2</label>
+<label group="0" index="69">7Networks_RH_DorsAttn_Post_3</label>
+<label group="0" index="70">7Networks_RH_DorsAttn_Post_4</label>
+<label group="0" index="71">7Networks_RH_DorsAttn_Post_5</label>
+<label group="0" index="72">7Networks_RH_DorsAttn_PrCv_1</label>
+<label group="0" index="73">7Networks_RH_DorsAttn_FEF_1</label>
+<label group="0" index="74">7Networks_RH_SalVentAttn_TempOccPar_1</label>
+<label group="0" index="75">7Networks_RH_SalVentAttn_TempOccPar_2</label>
+<label group="0" index="76">7Networks_RH_SalVentAttn_FrOperIns_1</label>
+<label group="0" index="77">7Networks_RH_SalVentAttn_Med_1</label>
+<label group="0" index="78">7Networks_RH_SalVentAttn_Med_2</label>
+<label group="0" index="79">7Networks_RH_Limbic_OFC_1</label>
+<label group="0" index="80">7Networks_RH_Limbic_TempPole_1</label>
+<label group="0" index="81">7Networks_RH_Cont_Par_1</label>
+<label group="0" index="82">7Networks_RH_Cont_Par_2</label>
+<label group="0" index="83">7Networks_RH_Cont_PFCl_1</label>
+<label group="0" index="84">7Networks_RH_Cont_PFCl_2</label>
+<label group="0" index="85">7Networks_RH_Cont_PFCl_3</label>
+<label group="0" index="86">7Networks_RH_Cont_PFCl_4</label>
+<label group="0" index="87">7Networks_RH_Cont_Cing_1</label>
+<label group="0" index="88">7Networks_RH_Cont_PFCmp_1</label>
+<label group="0" index="89">7Networks_RH_Cont_pCun_1</label>
+<label group="0" index="90">7Networks_RH_Default_Par_1</label>
+<label group="0" index="91">7Networks_RH_Default_Temp_1</label>
+<label group="0" index="92">7Networks_RH_Default_Temp_2</label>
+<label group="0" index="93">7Networks_RH_Default_Temp_3</label>
+<label group="0" index="94">7Networks_RH_Default_PFCv_1</label>
+<label group="0" index="95">7Networks_RH_Default_PFCv_2</label>
+<label group="0" index="96">7Networks_RH_Default_PFCdPFCm_1</label>
+<label group="0" index="97">7Networks_RH_Default_PFCdPFCm_2</label>
+<label group="0" index="98">7Networks_RH_Default_PFCdPFCm_3</label>
+<label group="0" index="99">7Networks_RH_Default_pCunPCC_1</label>
+<label group="0" index="100">7Networks_RH_Default_pCunPCC_2</label>
+<label group="0" index="101">LeftThalamus</label>
+<label group="0" index="102">LeftCaudate</label>
+<label group="0" index="103">LeftPutamen</label>
+<label group="0" index="104">LeftPallidum</label>
+<label group="0" index="105">LeftHippocampus</label>
+<label group="0" index="106">LeftAmygdala</label>
+<label group="0" index="107">LeftAccumbens</label>
+<label group="0" index="108">RightThalamus</label>
+<label group="0" index="109">RightCaudate</label>
+<label group="0" index="110">RightPutamen</label>
+<label group="0" index="111">RightPallidum</label>
+<label group="0" index="112">RightHippocampus</label>
+<label group="0" index="113">RightAmygdala</label>
+<label group="0" index="114">RightAccumbens</label>
 </data>
diff --git a/osl_ephys/source_recon/files/Schaefer2018_100Parcels_7Networks_order_FSLMNI152_2mm_4d.xml b/osl_ephys/source_recon/files/Schaefer2018_100Parcels_7Networks_order_FSLMNI152_2mm_4d.xml
index b2e7110..45889ee 100644
--- a/osl_ephys/source_recon/files/Schaefer2018_100Parcels_7Networks_order_FSLMNI152_2mm_4d.xml
+++ b/osl_ephys/source_recon/files/Schaefer2018_100Parcels_7Networks_order_FSLMNI152_2mm_4d.xml
@@ -1,102 +1,102 @@
 <data>
-<label index="0">7Networks_LH_Vis_1</label>
-<label index="1">7Networks_LH_Vis_2</label>
-<label index="2">7Networks_LH_Vis_3</label>
-<label index="3">7Networks_LH_Vis_4</label>
-<label index="4">7Networks_LH_Vis_5</label>
-<label index="5">7Networks_LH_Vis_6</label>
-<label index="6">7Networks_LH_Vis_7</label>
-<label index="7">7Networks_LH_Vis_8</label>
-<label index="8">7Networks_LH_Vis_9</label>
-<label index="9">7Networks_LH_SomMot_1</label>
-<label index="10">7Networks_LH_SomMot_2</label>
-<label index="11">7Networks_LH_SomMot_3</label>
-<label index="12">7Networks_LH_SomMot_4</label>
-<label index="13">7Networks_LH_SomMot_5</label>
-<label index="14">7Networks_LH_SomMot_6</label>
-<label index="15">7Networks_LH_DorsAttn_Post_1</label>
-<label index="16">7Networks_LH_DorsAttn_Post_2</label>
-<label index="17">7Networks_LH_DorsAttn_Post_3</label>
-<label index="18">7Networks_LH_DorsAttn_Post_4</label>
-<label index="19">7Networks_LH_DorsAttn_Post_5</label>
-<label index="20">7Networks_LH_DorsAttn_Post_6</label>
-<label index="21">7Networks_LH_DorsAttn_PrCv_1</label>
-<label index="22">7Networks_LH_DorsAttn_FEF_1</label>
-<label index="23">7Networks_LH_SalVentAttn_ParOper_1</label>
-<label index="24">7Networks_LH_SalVentAttn_FrOperIns_1</label>
-<label index="25">7Networks_LH_SalVentAttn_FrOperIns_2</label>
-<label index="26">7Networks_LH_SalVentAttn_PFCl_1</label>
-<label index="27">7Networks_LH_SalVentAttn_Med_1</label>
-<label index="28">7Networks_LH_SalVentAttn_Med_2</label>
-<label index="29">7Networks_LH_SalVentAttn_Med_3</label>
-<label index="30">7Networks_LH_Limbic_OFC_1</label>
-<label index="31">7Networks_LH_Limbic_TempPole_1</label>
-<label index="32">7Networks_LH_Limbic_TempPole_2</label>
-<label index="33">7Networks_LH_Cont_Par_1</label>
-<label index="34">7Networks_LH_Cont_PFCl_1</label>
-<label index="35">7Networks_LH_Cont_pCun_1</label>
-<label index="36">7Networks_LH_Cont_Cing_1</label>
-<label index="37">7Networks_LH_Default_Temp_1</label>
-<label index="38">7Networks_LH_Default_Temp_2</label>
-<label index="39">7Networks_LH_Default_Par_1</label>
-<label index="40">7Networks_LH_Default_Par_2</label>
-<label index="41">7Networks_LH_Default_PFC_1</label>
-<label index="42">7Networks_LH_Default_PFC_2</label>
-<label index="43">7Networks_LH_Default_PFC_3</label>
-<label index="44">7Networks_LH_Default_PFC_4</label>
-<label index="45">7Networks_LH_Default_PFC_5</label>
-<label index="46">7Networks_LH_Default_PFC_6</label>
-<label index="47">7Networks_LH_Default_PFC_7</label>
-<label index="48">7Networks_LH_Default_pCunPCC_1</label>
-<label index="49">7Networks_LH_Default_pCunPCC_2</label>
-<label index="50">7Networks_RH_Vis_1</label>
-<label index="51">7Networks_RH_Vis_2</label>
-<label index="52">7Networks_RH_Vis_3</label>
-<label index="53">7Networks_RH_Vis_4</label>
-<label index="54">7Networks_RH_Vis_5</label>
-<label index="55">7Networks_RH_Vis_6</label>
-<label index="56">7Networks_RH_Vis_7</label>
-<label index="57">7Networks_RH_Vis_8</label>
-<label index="58">7Networks_RH_SomMot_1</label>
-<label index="59">7Networks_RH_SomMot_2</label>
-<label index="60">7Networks_RH_SomMot_3</label>
-<label index="61">7Networks_RH_SomMot_4</label>
-<label index="62">7Networks_RH_SomMot_5</label>
-<label index="63">7Networks_RH_SomMot_6</label>
-<label index="64">7Networks_RH_SomMot_7</label>
-<label index="65">7Networks_RH_SomMot_8</label>
-<label index="66">7Networks_RH_DorsAttn_Post_1</label>
-<label index="67">7Networks_RH_DorsAttn_Post_2</label>
-<label index="68">7Networks_RH_DorsAttn_Post_3</label>
-<label index="69">7Networks_RH_DorsAttn_Post_4</label>
-<label index="70">7Networks_RH_DorsAttn_Post_5</label>
-<label index="71">7Networks_RH_DorsAttn_PrCv_1</label>
-<label index="72">7Networks_RH_DorsAttn_FEF_1</label>
-<label index="73">7Networks_RH_SalVentAttn_TempOccPar_1</label>
-<label index="74">7Networks_RH_SalVentAttn_TempOccPar_2</label>
-<label index="75">7Networks_RH_SalVentAttn_FrOperIns_1</label>
-<label index="76">7Networks_RH_SalVentAttn_Med_1</label>
-<label index="77">7Networks_RH_SalVentAttn_Med_2</label>
-<label index="78">7Networks_RH_Limbic_OFC_1</label>
-<label index="79">7Networks_RH_Limbic_TempPole_1</label>
-<label index="80">7Networks_RH_Cont_Par_1</label>
-<label index="81">7Networks_RH_Cont_Par_2</label>
-<label index="82">7Networks_RH_Cont_PFCl_1</label>
-<label index="83">7Networks_RH_Cont_PFCl_2</label>
-<label index="84">7Networks_RH_Cont_PFCl_3</label>
-<label index="85">7Networks_RH_Cont_PFCl_4</label>
-<label index="86">7Networks_RH_Cont_Cing_1</label>
-<label index="87">7Networks_RH_Cont_PFCmp_1</label>
-<label index="88">7Networks_RH_Cont_pCun_1</label>
-<label index="89">7Networks_RH_Default_Par_1</label>
-<label index="90">7Networks_RH_Default_Temp_1</label>
-<label index="91">7Networks_RH_Default_Temp_2</label>
-<label index="92">7Networks_RH_Default_Temp_3</label>
-<label index="93">7Networks_RH_Default_PFCv_1</label>
-<label index="94">7Networks_RH_Default_PFCv_2</label>
-<label index="95">7Networks_RH_Default_PFCdPFCm_1</label>
-<label index="96">7Networks_RH_Default_PFCdPFCm_2</label>
-<label index="97">7Networks_RH_Default_PFCdPFCm_3</label>
-<label index="98">7Networks_RH_Default_pCunPCC_1</label>
-<label index="99">7Networks_RH_Default_pCunPCC_2</label>
+<label group="0" index="0">7Networks_LH_Vis_1</label>
+<label group="0" index="1">7Networks_LH_Vis_2</label>
+<label group="0" index="2">7Networks_LH_Vis_3</label>
+<label group="0" index="3">7Networks_LH_Vis_4</label>
+<label group="0" index="4">7Networks_LH_Vis_5</label>
+<label group="0" index="5">7Networks_LH_Vis_6</label>
+<label group="0" index="6">7Networks_LH_Vis_7</label>
+<label group="0" index="7">7Networks_LH_Vis_8</label>
+<label group="0" index="8">7Networks_LH_Vis_9</label>
+<label group="0" index="9">7Networks_LH_SomMot_1</label>
+<label group="0" index="10">7Networks_LH_SomMot_2</label>
+<label group="0" index="11">7Networks_LH_SomMot_3</label>
+<label group="0" index="12">7Networks_LH_SomMot_4</label>
+<label group="0" index="13">7Networks_LH_SomMot_5</label>
+<label group="0" index="14">7Networks_LH_SomMot_6</label>
+<label group="0" index="15">7Networks_LH_DorsAttn_Post_1</label>
+<label group="0" index="16">7Networks_LH_DorsAttn_Post_2</label>
+<label group="0" index="17">7Networks_LH_DorsAttn_Post_3</label>
+<label group="0" index="18">7Networks_LH_DorsAttn_Post_4</label>
+<label group="0" index="19">7Networks_LH_DorsAttn_Post_5</label>
+<label group="0" index="20">7Networks_LH_DorsAttn_Post_6</label>
+<label group="0" index="21">7Networks_LH_DorsAttn_PrCv_1</label>
+<label group="0" index="22">7Networks_LH_DorsAttn_FEF_1</label>
+<label group="0" index="23">7Networks_LH_SalVentAttn_ParOper_1</label>
+<label group="0" index="24">7Networks_LH_SalVentAttn_FrOperIns_1</label>
+<label group="0" index="25">7Networks_LH_SalVentAttn_FrOperIns_2</label>
+<label group="0" index="26">7Networks_LH_SalVentAttn_PFCl_1</label>
+<label group="0" index="27">7Networks_LH_SalVentAttn_Med_1</label>
+<label group="0" index="28">7Networks_LH_SalVentAttn_Med_2</label>
+<label group="0" index="29">7Networks_LH_SalVentAttn_Med_3</label>
+<label group="0" index="30">7Networks_LH_Limbic_OFC_1</label>
+<label group="0" index="31">7Networks_LH_Limbic_TempPole_1</label>
+<label group="0" index="32">7Networks_LH_Limbic_TempPole_2</label>
+<label group="0" index="33">7Networks_LH_Cont_Par_1</label>
+<label group="0" index="34">7Networks_LH_Cont_PFCl_1</label>
+<label group="0" index="35">7Networks_LH_Cont_pCun_1</label>
+<label group="0" index="36">7Networks_LH_Cont_Cing_1</label>
+<label group="0" index="37">7Networks_LH_Default_Temp_1</label>
+<label group="0" index="38">7Networks_LH_Default_Temp_2</label>
+<label group="0" index="39">7Networks_LH_Default_Par_1</label>
+<label group="0" index="40">7Networks_LH_Default_Par_2</label>
+<label group="0" index="41">7Networks_LH_Default_PFC_1</label>
+<label group="0" index="42">7Networks_LH_Default_PFC_2</label>
+<label group="0" index="43">7Networks_LH_Default_PFC_3</label>
+<label group="0" index="44">7Networks_LH_Default_PFC_4</label>
+<label group="0" index="45">7Networks_LH_Default_PFC_5</label>
+<label group="0" index="46">7Networks_LH_Default_PFC_6</label>
+<label group="0" index="47">7Networks_LH_Default_PFC_7</label>
+<label group="0" index="48">7Networks_LH_Default_pCunPCC_1</label>
+<label group="0" index="49">7Networks_LH_Default_pCunPCC_2</label>
+<label group="0" index="50">7Networks_RH_Vis_1</label>
+<label group="0" index="51">7Networks_RH_Vis_2</label>
+<label group="0" index="52">7Networks_RH_Vis_3</label>
+<label group="0" index="53">7Networks_RH_Vis_4</label>
+<label group="0" index="54">7Networks_RH_Vis_5</label>
+<label group="0" index="55">7Networks_RH_Vis_6</label>
+<label group="0" index="56">7Networks_RH_Vis_7</label>
+<label group="0" index="57">7Networks_RH_Vis_8</label>
+<label group="0" index="58">7Networks_RH_SomMot_1</label>
+<label group="0" index="59">7Networks_RH_SomMot_2</label>
+<label group="0" index="60">7Networks_RH_SomMot_3</label>
+<label group="0" index="61">7Networks_RH_SomMot_4</label>
+<label group="0" index="62">7Networks_RH_SomMot_5</label>
+<label group="0" index="63">7Networks_RH_SomMot_6</label>
+<label group="0" index="64">7Networks_RH_SomMot_7</label>
+<label group="0" index="65">7Networks_RH_SomMot_8</label>
+<label group="0" index="66">7Networks_RH_DorsAttn_Post_1</label>
+<label group="0" index="67">7Networks_RH_DorsAttn_Post_2</label>
+<label group="0" index="68">7Networks_RH_DorsAttn_Post_3</label>
+<label group="0" index="69">7Networks_RH_DorsAttn_Post_4</label>
+<label group="0" index="70">7Networks_RH_DorsAttn_Post_5</label>
+<label group="0" index="71">7Networks_RH_DorsAttn_PrCv_1</label>
+<label group="0" index="72">7Networks_RH_DorsAttn_FEF_1</label>
+<label group="0" index="73">7Networks_RH_SalVentAttn_TempOccPar_1</label>
+<label group="0" index="74">7Networks_RH_SalVentAttn_TempOccPar_2</label>
+<label group="0" index="75">7Networks_RH_SalVentAttn_FrOperIns_1</label>
+<label group="0" index="76">7Networks_RH_SalVentAttn_Med_1</label>
+<label group="0" index="77">7Networks_RH_SalVentAttn_Med_2</label>
+<label group="0" index="78">7Networks_RH_Limbic_OFC_1</label>
+<label group="0" index="79">7Networks_RH_Limbic_TempPole_1</label>
+<label group="0" index="80">7Networks_RH_Cont_Par_1</label>
+<label group="0" index="81">7Networks_RH_Cont_Par_2</label>
+<label group="0" index="82">7Networks_RH_Cont_PFCl_1</label>
+<label group="0" index="83">7Networks_RH_Cont_PFCl_2</label>
+<label group="0" index="84">7Networks_RH_Cont_PFCl_3</label>
+<label group="0" index="85">7Networks_RH_Cont_PFCl_4</label>
+<label group="0" index="86">7Networks_RH_Cont_Cing_1</label>
+<label group="0" index="87">7Networks_RH_Cont_PFCmp_1</label>
+<label group="0" index="88">7Networks_RH_Cont_pCun_1</label>
+<label group="0" index="89">7Networks_RH_Default_Par_1</label>
+<label group="0" index="90">7Networks_RH_Default_Temp_1</label>
+<label group="0" index="91">7Networks_RH_Default_Temp_2</label>
+<label group="0" index="92">7Networks_RH_Default_Temp_3</label>
+<label group="0" index="93">7Networks_RH_Default_PFCv_1</label>
+<label group="0" index="94">7Networks_RH_Default_PFCv_2</label>
+<label group="0" index="95">7Networks_RH_Default_PFCdPFCm_1</label>
+<label group="0" index="96">7Networks_RH_Default_PFCdPFCm_2</label>
+<label group="0" index="97">7Networks_RH_Default_PFCdPFCm_3</label>
+<label group="0" index="98">7Networks_RH_Default_pCunPCC_1</label>
+<label group="0" index="99">7Networks_RH_Default_pCunPCC_2</label>
 </data>
diff --git a/osl_ephys/source_recon/files/Schaefer2018_100Parcels_7Networks_order_FSLMNI152_2mm_4d_ds8.xml b/osl_ephys/source_recon/files/Schaefer2018_100Parcels_7Networks_order_FSLMNI152_2mm_4d_ds8.xml
index b2e7110..2b08ea0 100644
--- a/osl_ephys/source_recon/files/Schaefer2018_100Parcels_7Networks_order_FSLMNI152_2mm_4d_ds8.xml
+++ b/osl_ephys/source_recon/files/Schaefer2018_100Parcels_7Networks_order_FSLMNI152_2mm_4d_ds8.xml
@@ -1,102 +1,102 @@
 <data>
-<label index="0">7Networks_LH_Vis_1</label>
-<label index="1">7Networks_LH_Vis_2</label>
-<label index="2">7Networks_LH_Vis_3</label>
-<label index="3">7Networks_LH_Vis_4</label>
-<label index="4">7Networks_LH_Vis_5</label>
-<label index="5">7Networks_LH_Vis_6</label>
-<label index="6">7Networks_LH_Vis_7</label>
-<label index="7">7Networks_LH_Vis_8</label>
-<label index="8">7Networks_LH_Vis_9</label>
-<label index="9">7Networks_LH_SomMot_1</label>
-<label index="10">7Networks_LH_SomMot_2</label>
-<label index="11">7Networks_LH_SomMot_3</label>
-<label index="12">7Networks_LH_SomMot_4</label>
-<label index="13">7Networks_LH_SomMot_5</label>
-<label index="14">7Networks_LH_SomMot_6</label>
-<label index="15">7Networks_LH_DorsAttn_Post_1</label>
-<label index="16">7Networks_LH_DorsAttn_Post_2</label>
-<label index="17">7Networks_LH_DorsAttn_Post_3</label>
-<label index="18">7Networks_LH_DorsAttn_Post_4</label>
-<label index="19">7Networks_LH_DorsAttn_Post_5</label>
-<label index="20">7Networks_LH_DorsAttn_Post_6</label>
-<label index="21">7Networks_LH_DorsAttn_PrCv_1</label>
-<label index="22">7Networks_LH_DorsAttn_FEF_1</label>
-<label index="23">7Networks_LH_SalVentAttn_ParOper_1</label>
-<label index="24">7Networks_LH_SalVentAttn_FrOperIns_1</label>
-<label index="25">7Networks_LH_SalVentAttn_FrOperIns_2</label>
-<label index="26">7Networks_LH_SalVentAttn_PFCl_1</label>
-<label index="27">7Networks_LH_SalVentAttn_Med_1</label>
-<label index="28">7Networks_LH_SalVentAttn_Med_2</label>
-<label index="29">7Networks_LH_SalVentAttn_Med_3</label>
-<label index="30">7Networks_LH_Limbic_OFC_1</label>
-<label index="31">7Networks_LH_Limbic_TempPole_1</label>
-<label index="32">7Networks_LH_Limbic_TempPole_2</label>
-<label index="33">7Networks_LH_Cont_Par_1</label>
-<label index="34">7Networks_LH_Cont_PFCl_1</label>
-<label index="35">7Networks_LH_Cont_pCun_1</label>
-<label index="36">7Networks_LH_Cont_Cing_1</label>
-<label index="37">7Networks_LH_Default_Temp_1</label>
-<label index="38">7Networks_LH_Default_Temp_2</label>
-<label index="39">7Networks_LH_Default_Par_1</label>
-<label index="40">7Networks_LH_Default_Par_2</label>
-<label index="41">7Networks_LH_Default_PFC_1</label>
-<label index="42">7Networks_LH_Default_PFC_2</label>
-<label index="43">7Networks_LH_Default_PFC_3</label>
-<label index="44">7Networks_LH_Default_PFC_4</label>
-<label index="45">7Networks_LH_Default_PFC_5</label>
-<label index="46">7Networks_LH_Default_PFC_6</label>
-<label index="47">7Networks_LH_Default_PFC_7</label>
-<label index="48">7Networks_LH_Default_pCunPCC_1</label>
-<label index="49">7Networks_LH_Default_pCunPCC_2</label>
-<label index="50">7Networks_RH_Vis_1</label>
-<label index="51">7Networks_RH_Vis_2</label>
-<label index="52">7Networks_RH_Vis_3</label>
-<label index="53">7Networks_RH_Vis_4</label>
-<label index="54">7Networks_RH_Vis_5</label>
-<label index="55">7Networks_RH_Vis_6</label>
-<label index="56">7Networks_RH_Vis_7</label>
-<label index="57">7Networks_RH_Vis_8</label>
-<label index="58">7Networks_RH_SomMot_1</label>
-<label index="59">7Networks_RH_SomMot_2</label>
-<label index="60">7Networks_RH_SomMot_3</label>
-<label index="61">7Networks_RH_SomMot_4</label>
-<label index="62">7Networks_RH_SomMot_5</label>
-<label index="63">7Networks_RH_SomMot_6</label>
-<label index="64">7Networks_RH_SomMot_7</label>
-<label index="65">7Networks_RH_SomMot_8</label>
-<label index="66">7Networks_RH_DorsAttn_Post_1</label>
-<label index="67">7Networks_RH_DorsAttn_Post_2</label>
-<label index="68">7Networks_RH_DorsAttn_Post_3</label>
-<label index="69">7Networks_RH_DorsAttn_Post_4</label>
-<label index="70">7Networks_RH_DorsAttn_Post_5</label>
-<label index="71">7Networks_RH_DorsAttn_PrCv_1</label>
-<label index="72">7Networks_RH_DorsAttn_FEF_1</label>
-<label index="73">7Networks_RH_SalVentAttn_TempOccPar_1</label>
-<label index="74">7Networks_RH_SalVentAttn_TempOccPar_2</label>
-<label index="75">7Networks_RH_SalVentAttn_FrOperIns_1</label>
-<label index="76">7Networks_RH_SalVentAttn_Med_1</label>
-<label index="77">7Networks_RH_SalVentAttn_Med_2</label>
-<label index="78">7Networks_RH_Limbic_OFC_1</label>
-<label index="79">7Networks_RH_Limbic_TempPole_1</label>
-<label index="80">7Networks_RH_Cont_Par_1</label>
-<label index="81">7Networks_RH_Cont_Par_2</label>
-<label index="82">7Networks_RH_Cont_PFCl_1</label>
-<label index="83">7Networks_RH_Cont_PFCl_2</label>
-<label index="84">7Networks_RH_Cont_PFCl_3</label>
-<label index="85">7Networks_RH_Cont_PFCl_4</label>
-<label index="86">7Networks_RH_Cont_Cing_1</label>
-<label index="87">7Networks_RH_Cont_PFCmp_1</label>
-<label index="88">7Networks_RH_Cont_pCun_1</label>
-<label index="89">7Networks_RH_Default_Par_1</label>
-<label index="90">7Networks_RH_Default_Temp_1</label>
-<label index="91">7Networks_RH_Default_Temp_2</label>
-<label index="92">7Networks_RH_Default_Temp_3</label>
-<label index="93">7Networks_RH_Default_PFCv_1</label>
-<label index="94">7Networks_RH_Default_PFCv_2</label>
-<label index="95">7Networks_RH_Default_PFCdPFCm_1</label>
-<label index="96">7Networks_RH_Default_PFCdPFCm_2</label>
-<label index="97">7Networks_RH_Default_PFCdPFCm_3</label>
-<label index="98">7Networks_RH_Default_pCunPCC_1</label>
-<label index="99">7Networks_RH_Default_pCunPCC_2</label>
+<label group="0" group="0" index="0">7Networks_LH_Vis_1</label>
+<label group="0" group="0" index="1">7Networks_LH_Vis_2</label>
+<label group="0" group="0" index="2">7Networks_LH_Vis_3</label>
+<label group="0" group="0" index="3">7Networks_LH_Vis_4</label>
+<label group="0" index="4">7Networks_LH_Vis_5</label>
+<label group="0" index="5">7Networks_LH_Vis_6</label>
+<label group="0" index="6">7Networks_LH_Vis_7</label>
+<label group="0" index="7">7Networks_LH_Vis_8</label>
+<label group="0" index="8">7Networks_LH_Vis_9</label>
+<label group="0" index="9">7Networks_LH_SomMot_1</label>
+<label group="0" index="10">7Networks_LH_SomMot_2</label>
+<label group="0" index="11">7Networks_LH_SomMot_3</label>
+<label group="0" index="12">7Networks_LH_SomMot_4</label>
+<label group="0" index="13">7Networks_LH_SomMot_5</label>
+<label group="0" index="14">7Networks_LH_SomMot_6</label>
+<label group="0" index="15">7Networks_LH_DorsAttn_Post_1</label>
+<label group="0" index="16">7Networks_LH_DorsAttn_Post_2</label>
+<label group="0" index="17">7Networks_LH_DorsAttn_Post_3</label>
+<label group="0" index="18">7Networks_LH_DorsAttn_Post_4</label>
+<label group="0" index="19">7Networks_LH_DorsAttn_Post_5</label>
+<label group="0" index="20">7Networks_LH_DorsAttn_Post_6</label>
+<label group="0" index="21">7Networks_LH_DorsAttn_PrCv_1</label>
+<label group="0" index="22">7Networks_LH_DorsAttn_FEF_1</label>
+<label group="0" index="23">7Networks_LH_SalVentAttn_ParOper_1</label>
+<label group="0" index="24">7Networks_LH_SalVentAttn_FrOperIns_1</label>
+<label group="0" index="25">7Networks_LH_SalVentAttn_FrOperIns_2</label>
+<label group="0" index="26">7Networks_LH_SalVentAttn_PFCl_1</label>
+<label group="0" index="27">7Networks_LH_SalVentAttn_Med_1</label>
+<label group="0" index="28">7Networks_LH_SalVentAttn_Med_2</label>
+<label group="0" index="29">7Networks_LH_SalVentAttn_Med_3</label>
+<label group="0" index="30">7Networks_LH_Limbic_OFC_1</label>
+<label group="0" index="31">7Networks_LH_Limbic_TempPole_1</label>
+<label group="0" index="32">7Networks_LH_Limbic_TempPole_2</label>
+<label group="0" index="33">7Networks_LH_Cont_Par_1</label>
+<label group="0" index="34">7Networks_LH_Cont_PFCl_1</label>
+<label group="0" index="35">7Networks_LH_Cont_pCun_1</label>
+<label group="0" index="36">7Networks_LH_Cont_Cing_1</label>
+<label group="0" index="37">7Networks_LH_Default_Temp_1</label>
+<label group="0" index="38">7Networks_LH_Default_Temp_2</label>
+<label group="0" index="39">7Networks_LH_Default_Par_1</label>
+<label group="0" index="40">7Networks_LH_Default_Par_2</label>
+<label group="0" index="41">7Networks_LH_Default_PFC_1</label>
+<label group="0" index="42">7Networks_LH_Default_PFC_2</label>
+<label group="0" index="43">7Networks_LH_Default_PFC_3</label>
+<label group="0" index="44">7Networks_LH_Default_PFC_4</label>
+<label group="0" index="45">7Networks_LH_Default_PFC_5</label>
+<label group="0" index="46">7Networks_LH_Default_PFC_6</label>
+<label group="0" index="47">7Networks_LH_Default_PFC_7</label>
+<label group="0" index="48">7Networks_LH_Default_pCunPCC_1</label>
+<label group="0" index="49">7Networks_LH_Default_pCunPCC_2</label>
+<label group="0" index="50">7Networks_RH_Vis_1</label>
+<label group="0" index="51">7Networks_RH_Vis_2</label>
+<label group="0" index="52">7Networks_RH_Vis_3</label>
+<label group="0" index="53">7Networks_RH_Vis_4</label>
+<label group="0" index="54">7Networks_RH_Vis_5</label>
+<label group="0" index="55">7Networks_RH_Vis_6</label>
+<label group="0" index="56">7Networks_RH_Vis_7</label>
+<label group="0" index="57">7Networks_RH_Vis_8</label>
+<label group="0" index="58">7Networks_RH_SomMot_1</label>
+<label group="0" index="59">7Networks_RH_SomMot_2</label>
+<label group="0" index="60">7Networks_RH_SomMot_3</label>
+<label group="0" index="61">7Networks_RH_SomMot_4</label>
+<label group="0" index="62">7Networks_RH_SomMot_5</label>
+<label group="0" index="63">7Networks_RH_SomMot_6</label>
+<label group="0" index="64">7Networks_RH_SomMot_7</label>
+<label group="0" index="65">7Networks_RH_SomMot_8</label>
+<label group="0" index="66">7Networks_RH_DorsAttn_Post_1</label>
+<label group="0" index="67">7Networks_RH_DorsAttn_Post_2</label>
+<label group="0" index="68">7Networks_RH_DorsAttn_Post_3</label>
+<label group="0" index="69">7Networks_RH_DorsAttn_Post_4</label>
+<label group="0" index="70">7Networks_RH_DorsAttn_Post_5</label>
+<label group="0" index="71">7Networks_RH_DorsAttn_PrCv_1</label>
+<label group="0" index="72">7Networks_RH_DorsAttn_FEF_1</label>
+<label group="0" index="73">7Networks_RH_SalVentAttn_TempOccPar_1</label>
+<label group="0" index="74">7Networks_RH_SalVentAttn_TempOccPar_2</label>
+<label group="0" index="75">7Networks_RH_SalVentAttn_FrOperIns_1</label>
+<label group="0" index="76">7Networks_RH_SalVentAttn_Med_1</label>
+<label group="0" index="77">7Networks_RH_SalVentAttn_Med_2</label>
+<label group="0" index="78">7Networks_RH_Limbic_OFC_1</label>
+<label group="0" index="79">7Networks_RH_Limbic_TempPole_1</label>
+<label group="0" index="80">7Networks_RH_Cont_Par_1</label>
+<label group="0" index="81">7Networks_RH_Cont_Par_2</label>
+<label group="0" index="82">7Networks_RH_Cont_PFCl_1</label>
+<label group="0" index="83">7Networks_RH_Cont_PFCl_2</label>
+<label group="0" index="84">7Networks_RH_Cont_PFCl_3</label>
+<label group="0" index="85">7Networks_RH_Cont_PFCl_4</label>
+<label group="0" index="86">7Networks_RH_Cont_Cing_1</label>
+<label group="0" index="87">7Networks_RH_Cont_PFCmp_1</label>
+<label group="0" index="88">7Networks_RH_Cont_pCun_1</label>
+<label group="0" index="89">7Networks_RH_Default_Par_1</label>
+<label group="0" index="90">7Networks_RH_Default_Temp_1</label>
+<label group="0" index="91">7Networks_RH_Default_Temp_2</label>
+<label group="0" index="92">7Networks_RH_Default_Temp_3</label>
+<label group="0" index="93">7Networks_RH_Default_PFCv_1</label>
+<label group="0" index="94">7Networks_RH_Default_PFCv_2</label>
+<label group="0" index="95">7Networks_RH_Default_PFCdPFCm_1</label>
+<label group="0" index="96">7Networks_RH_Default_PFCdPFCm_2</label>
+<label group="0" index="97">7Networks_RH_Default_PFCdPFCm_3</label>
+<label group="0" index="98">7Networks_RH_Default_pCunPCC_1</label>
+<label group="0" index="99">7Networks_RH_Default_pCunPCC_2</label>
 </data>
diff --git a/osl_ephys/source_recon/parcellation.py b/osl_ephys/source_recon/parcellation.py
index 3ea265e..72273be 100644
--- a/osl_ephys/source_recon/parcellation.py
+++ b/osl_ephys/source_recon/parcellation.py
@@ -10,7 +10,8 @@
 import os.path as op
 from pathlib import Path
 
-import fsl
+from fsl import wrappers as fsl_wrappers
+
 import mne
 import numpy as np
 import nibabel as nib
@@ -23,8 +24,11 @@
 from mpl_toolkits.axes_grid1 import make_axes_locatable
 
 import osl_ephys.source_recon.rhino.utils as rhino_utils
+
 from osl_ephys.utils.logger import log_or_print
 from osl_ephys.source_recon import freesurfer_utils
+from osl_ephys.source_recon.utils import find_file
+
 import pandas as pd
 
 def load_parcellation_description(parcellation_file, freesurfer=False, subject=None):
@@ -44,7 +48,7 @@ def load_parcellation_description(parcellation_file, freesurfer=False, subject=N
     Returns
     -------
     df : dict
-        Dictionary with keys 'index', 'label', 'x', 'y', 'z' for each parcel.
+        Dictionary with keys 'index', 'group', 'label', 'x', 'y', 'z' for each parcel.
     
     """
 
@@ -52,18 +56,17 @@ def load_parcellation_description(parcellation_file, freesurfer=False, subject=N
     parcellation_file = find_file(parcellation_file, freesurfer=False)
 
     if parcellation_file is not None:
-        # look for xml file with same name
+        # if file extension is .nii or .nii.gz replace with .xml
         xml_file = parcellation_file.replace(".nii.gz", ".xml").replace(".nii", ".xml")
-        print(xml_file)
+        log_or_print(f"Looking for {xml_file}")
         if op.exists(xml_file):
             # load xml_file without parcellation library
             df = pd.read_xml(xml_file)
-
             # if no index column, use ordinal index
             df = df.set_index('index') if 'index' in df.columns else df.reset_index(drop=True)
 
             # keep only relevant columns
-            df = df.loc[:, df.columns.str.lower().isin(["index", "label", "x", "y", "z"])]
+            df = df.loc[:, df.columns.str.lower().isin(["index", "group", "label", "x", "y", "z"])]
 
             # turn df into a dictionary using the index as keys
             parcellation_dict = df.to_dict(orient='index')
@@ -118,50 +121,6 @@ def load_parcellation(parcellation_file, freesurfer=False, subject=None):
     return None
 
 
-def _find_package_file(filename):
-    files_dir = str(Path(__file__).parent) + "/files/"
-    if op.exists(files_dir + filename):
-        return files_dir + filename
-
-
-def _find_freesurfer_file(filename):
-    avail = mne.label._read_annot_cands(
-        os.path.join(os.environ["SUBJECTS_DIR"], 'fsaverage', 'label')
-    )
-    if filename in avail:
-        filename, hemis = mne.label._get_annot_fname(
-            None, 'fsaverage', 'both', filename, os.environ['SUBJECTS_DIR']
-        )
-        return filename
-
-
-def find_file(filename, freesurfer=False):
-    """Look for a parcellation file within the package.
-
-    Parameters
-    ----------
-    filename : str
-        Path to parcellation file to look for.
-    freesurfer : bool, optional
-        Should we look in the freesurfer directory?
-
-    Returns
-    -------
-    filename : str
-        Path to parcellation file found.
-    """
-    if not op.exists(filename):
-        if freesurfer:
-            filename = _find_freesurfer_file(filename)
-        else:
-            filename = _find_package_file(filename)
-        
-    if filename is None:
-        raise FileNotFoundError(filename)
-
-    return filename
-
-
 def guess_parcellation(data, return_path=False):
     """Guess parcellation file from data.
     
@@ -509,7 +468,18 @@ def resample_parcellation(parcellation_file, voxel_coords, working_dir=None, fre
     #   but this functionality is deprecated and will be removed in a future release.
     #
     # However, it doesn't look like the input be being truncated, the resampled parcellation appears to be a 4D volume.
-    fsl.wrappers.flirt(parcellation_file, parcellation_file, out=parcellation_resampled, applyisoxfm=gridstep)
+
+    # Check if parcellation is binary 
+    binary_mask = np.array_equal(np.unique(nib.load(parcellation_file).get_fdata()), [0, 1])
+    if binary_mask:
+        interp = "nearestneighbour"
+    else:
+        interp = "trilinear" 
+
+    fsl_wrappers.flirt(parcellation_file, parcellation_file, out=parcellation_resampled, applyisoxfm=gridstep, interp=interp)
+
+    if nib.load(parcellation_resampled).ndim != 4:
+        raise ValueError(f"Resampled parcellation {parcellation_resampled} should be a 4D NIfTI file.")
 
     nparcels = nib.load(parcellation_resampled).get_fdata().shape[3]
 
@@ -531,7 +501,6 @@ def resample_parcellation(parcellation_file, voxel_coords, working_dir=None, fre
 
     return parcellation_asmatrix
 
-
 def local_orthogonalise(timeseries, parcellation_file=None, dist=None, adjacency=None, ):
     """Returns a local orthogonalisation of the timeseries, where the time series of local neighbours are regressed out with multiple linear regression.
     
@@ -1372,3 +1341,60 @@ def plot_source_topo(
     # despite the options of vmin, vmax, the colorbar is always set to -vmax to vmax. correct this
     # plt.gca().get_images()[0].set_clim(vmin, vmax)
     return plt.gca().get_images()[0]
+
+
+def assign_voxels_to_binary_parcellation(voxel_coords, parcel_file):
+    """
+    Assigns voxel to parcels in a 4D parcellation file.
+    Voxel coords and parcel file must be in same coordinate space.
+    Parcel file must be a binary parcellation (0 for background, 1 for parcel)
+    Output is a binary parcellation with each voxel assigned to at most one parcel.
+
+    Assigns voxel to the closest parcel if within gridstep distance.
+    If no parcel is within gridstep distance, then do not assign to any parcel.
+    If multiple parcels are within gridstep distance, then assign to the closest parcel.
+
+    Parameters
+    ----------
+    voxel_coords : (nvoxels, 3) np.array
+        The 3D coords of the voxels to assign to parcels.
+    parcel_file : str
+        Path to a 4D nifti binary parcellation file where 
+        each volume should be a mask with 0 for the background 
+        and 1 for the parcel.
+        File must be in same coordinate space as voxel_coords.
+    Returns
+    -------
+    parcellation_assignment : list[int]
+        List of nvoxels parcel indices assigned to each voxel (-1 for no parcel assigned).
+    """
+
+    parcel_file = find_file(parcel_file)
+    nparcels = nib.load(parcel_file).get_fdata().shape[3]
+
+    kdtrees = []
+    for parcel_ind in range(nparcels):
+        coords_mni, _ = rhino_utils.niimask2mmpointcloud(parcel_file, volindex=parcel_ind)
+        kdtrees.append(KDTree(coords_mni.T))
+
+    gridstep = int(rhino_utils.get_gridstep(coords_mni.T) / 1000)
+
+    parcellation_assignment = []
+    for ind in range(voxel_coords.shape[0]):
+
+        # Assigns voxel to the closest parcel if within gridstep distance.
+        # If no parcel is within gridstep distance, then do not assign to any parcel.
+        # If multiple parcels are within gridstep distance, then assign to the closest parcel.
+        best_distance = np.inf
+        best_index = -1
+        for parcel_ind, kdtree in enumerate(kdtrees):
+            distance, _ = kdtree.query(voxel_coords[ind, :])
+            if distance < gridstep and distance < best_distance:
+                best_distance = distance
+                best_index = parcel_ind    
+        if best_index >= 0:
+            parcellation_assignment.append(best_index)
+        else:
+            parcellation_assignment.append(-1)  # -1 for no parcel assigned
+
+    return parcellation_assignment
diff --git a/osl_ephys/source_recon/rhino/coreg.py b/osl_ephys/source_recon/rhino/coreg.py
index 0329037..35ac50e 100644
--- a/osl_ephys/source_recon/rhino/coreg.py
+++ b/osl_ephys/source_recon/rhino/coreg.py
@@ -135,6 +135,7 @@ def coreg(
     # MRI (native) -- mri_mrivoxel_t (native2nativeindex) --> MRI (native) voxel indices
     #
     # RHINO does everthing in mm
+    # (MNE does everthing in metres)
 
     log_or_print("*** RUNNING OSL-EPHYS RHINO COREGISTRATION ***")
 
@@ -461,22 +462,29 @@ def coreg_display(
     subject,
     plot_type="surf",
     display_outskin=True,
-    display_outskin_with_nose=True,
+    display_outskin_with_nose=False,
     display_sensors=True,
     display_sensor_oris=True,
     display_fiducials=True,
     display_headshape_pnts=True,
+    display_inskull=False,    
+    display_dipoles=False,
+    multi_dipoles_display_mode="multicoloured_points",
+    multi_dipoles=None,
+    single_dipoles=None,
+    midline_points=None,
     filename=None,
+    src_index=0,    
 ):
     """Display coregistration.
 
-    Displays the coregistered RHINO scalp surface and polhemus/sensor locations.
+    3D display of the coregistered RHINO surfaces, polhemus/sensor locations etc.
 
     Display is done in MEG (device) space (in mm).
 
     Purple dots are the polhemus derived fiducials (these only get used to initialse the coreg, if headshape points are being used).
 
-    Yellow diamonds are the MNI standard space derived fiducials (these are the ones that matter).
+    Yellow diamonds are the MNI standard space derived fiducials.
 
     Parameters
     ----------
@@ -499,7 +507,38 @@ def coreg_display(
     display_fiducials - bool
         Whether to include fiducials in the display.
     display_headshape_pnts - bool
-        Whether to include headshape points in the display.
+        Whether to include headshape points in the display.    
+    display_inskull : bool
+        Whether to show in-skull surface in the display.
+    display_dipoles : bool  
+        Whether to show dipoles in the display. Dipoles are only shown if a forward model has been computed.
+    multi_dipoles_display_mode : string
+        Either:
+            'multicoloured_points' to show dipoles as colored points with color indicating which dipoles are part of the same multidipole set.
+            'points' to show all dipoles that are part of any patch as small red points
+            'linepairs' to show pairs of dipoles as red line pairs. Only works when number
+            of dipoles in each multidipole set is 2.
+    multi_dipoles : list of np.ndarray
+        Each element of the list is a set of indices corresponding to the dipoles to plot that
+        will be beamformed together as a multi-dipole set.
+        E.g. for bilateral dipoles, each element of the list will contain 2 indices.
+        E.g. for a patch, each element of the list will contain all the indices of the dipoles in the patch.
+        Indices correspond to the order of the single dipoles in the forward model source space.
+        Set to None to not plot any multidipole sets.
+    single_dipoles : list[int]
+        Each element of the list is a set of indices corresponding to the dipoles to plot as single dipoles.
+        Indices correspond to the order of the single dipoles in the forward model source space.
+        Set to None to not plot any single dipoles -  
+        unless multi_dipoles is also None, in which case all dipoles
+        in the forward model source space will be plotted as single dipoles.
+    midline_points : list[int]
+        Dipoles that lie on the midline of the head.
+        E.g. will be set if doing a bilateral beamformer (where they are treated as single dipoles).
+        Only ever used for plotting purposes.
+        Set to None to not plot any midline dipoles.
+    src_index : int
+        Index of the source space to use in the forward model, if there is more than one source space.
+        Default is 0, which is the first source space.
     filename : str
         Filename to save display to (as an interactive html).
         Must have extension .html.
@@ -519,6 +558,8 @@ def coreg_display(
     bet_outskin_mesh_vtk_file = coreg_filenames["bet_outskin_mesh_vtk_file"]
     bet_outskin_surf_file = coreg_filenames["bet_outskin_surf_file"]
     bet_outskin_plus_nose_surf_file = coreg_filenames["bet_outskin_plus_nose_surf_file"]
+    bet_inskull_mesh_file = coreg_filenames["bet_inskull_mesh_file"]
+    bet_inskull_surf_file = coreg_filenames["bet_inskull_surf_file"]
 
     head_scaledmri_t_file = coreg_filenames["head_scaledmri_t_file"]
     mrivoxel_scaledmri_t_file = coreg_filenames["mrivoxel_scaledmri_t_file"]
@@ -540,6 +581,18 @@ def coreg_display(
         outskin_mesh_4surf_file = bet_outskin_mesh_vtk_file
         outskin_surf_file = bet_outskin_surf_file
 
+    if display_dipoles:
+        rhino_files = rhino_utils.get_rhino_filenames(subjects_dir, subject)
+
+        fwd_fname = rhino_files["fwd_model"]
+        if Path(fwd_fname).exists():
+            forward = read_forward_solution(fwd_fname)
+            src = forward["src"]
+        else:
+            log_or_print("No forward model computed. Not displaying dipoles.")
+            src = None
+            display_dipoles = False
+
     # ------------
     # Setup xforms
 
@@ -669,6 +722,31 @@ def coreg_display(
         meg_sensor_locs = meg_sensor_locs * 1000
         meg_sensor_oris = meg_sensor_oris * 1000
 
+    # ----------------------------
+    # Setup vol source grid points
+
+    if display_dipoles:
+        # stored points are in metres, convert to mm
+        src_pnts = src[src_index]['rr'][src[src_index]['inuse']==1, :] * 1000
+
+        # Move from head space to MEG (device) space
+        src_pnts = rhino_utils.xform_points(head_trans["trans"], src_pnts.T).T
+
+        log_or_print("BEM surface: number of dipoles = {}".format(src_pnts.shape[0]))
+
+        if single_dipoles is not None:
+            single_dipoles = src_pnts[single_dipoles, :]
+        else:
+            if multi_dipoles is None:
+                single_dipoles = src_pnts
+                log_or_print("No dipole types specified, displaying all dipoles as single dipoles")
+            else:
+                single_dipoles = None
+                log_or_print("Will not display any single dipoles")
+
+        if midline_points is not None:
+            midline_points = src_pnts[midline_points, :]
+
     # --------
     # Do plots
     with warnings.catch_warnings():
@@ -771,6 +849,93 @@ def coreg_display(
 
                 renderer.surface(surface=surf_smri, color=(0, 0.7, 0.7), opacity=0.4, backface_culling=False)
 
+            if display_inskull:
+                # Inner skull surface
+                # Load in surface, this is in mm
+                coords_native, faces = nib.freesurfer.read_geometry(bet_inskull_surf_file)
+
+                # Move to MEG (device) space
+                coords_meg = rhino_utils.xform_points(mri_trans["trans"], coords_native.T).T
+
+                surf_smri = dict(rr=coords_meg, tris=faces)
+
+                renderer.surface(surface=surf_smri, color=(0.25, 0.25, 0.25), opacity=0.1, backface_culling=False)
+
+            # dipole (vol source) grid points
+            if display_dipoles:
+
+                if single_dipoles is not None and len(single_dipoles.T) > 0:
+
+                    color, alpha = "blue", 0.5
+                    scale = 0.0015
+
+                    if multi_dipoles is not None and len(multi_dipoles) > 0: 
+                        # make single dipoles a bit bigger if also plotting multi_dipoles
+                        scale = 0.002
+
+                    renderer.sphere(center=single_dipoles, 
+                                        color=color, 
+                                        scale=scale * 1000, 
+                                        opacity=alpha, 
+                                        backface_culling=True)
+
+                if midline_points is not None and len(midline_points.T) > 0:
+
+                    color, alpha = "green", 1
+                    scale = 0.0025
+                    renderer.sphere(center=midline_points, 
+                                    color=color, 
+                                    scale=scale * 1000, 
+                                    opacity=alpha, 
+                                    backface_culling=True)
+                
+                if multi_dipoles is not None and len(multi_dipoles) > 0:
+
+                    if multi_dipoles_display_mode == "multicoloured_points":
+                        colors_in = plt.cm.get_cmap("tab20", len(multi_dipoles)).colors
+                    else:
+                        colors_in = ["red"] * len(multi_dipoles)
+
+                    for index, pp in enumerate(multi_dipoles):
+                        if multi_dipoles_display_mode == "linepairs":
+                            # only show lines if number of dipoles in multi_dipoles is 2
+                            if len(pp) != 2:
+                                multi_dipoles_display_mode = "points"
+                                log_or_print("multi_dipoles_display_mode set to 'points' as number of dipoles is not 2")
+
+                        if multi_dipoles_display_mode == "points" or multi_dipoles_display_mode == "multicoloured_points":
+                            alpha = 0.5
+                            scale = 0.002
+
+                            for i in range(len(pp)):
+                                # spheres
+                                renderer.sphere(center=np.array([[src_pnts[pp[i], 0], 
+                                                                    src_pnts[pp[i], 1],
+                                                                    src_pnts[pp[i], 2]]]), 
+                                                    color=colors_in[index], 
+                                                    scale=scale * 1000, 
+                                                    opacity=alpha, 
+                                                    backface_culling=True)
+                            
+                            
+                        elif multi_dipoles_display_mode == "linepairs":
+                            # show a line connecting each pair
+                            color = "red"
+                            alpha = 0.4
+                            radius = 0.2
+                            # mne/viz/backends/_abstract.py
+                            renderer.tube(origin = np.array([[src_pnts[pp[0], 0], 
+                                                    src_pnts[pp[0], 1],
+                                                    src_pnts[pp[0], 2]]]), 
+                                            destination = np.array([[src_pnts[pp[1], 0], 
+                                                        src_pnts[pp[1], 1],
+                                                        src_pnts[pp[1], 2]]]),
+                                            radius=radius,
+                                            opacity=alpha,
+                                            color=color)
+                        else:
+                            ValueError(f'{multi_dipoles_display_mode} is an invalid multi_dipoles_display_mode')
+
             renderer.set_camera(azimuth=90, elevation=90, distance=600, focalpoint=(0.0, 0.0, 0.0))
 
             # Save or show
@@ -782,15 +947,31 @@ def coreg_display(
             # -------------------
             # Setup scalp surface
 
-            # Load in scalp surface
-            # And turn the nvoxx x nvoxy x nvoxz volume into a 3 x npoints point cloud
-            smri_headshape_nativeindex = rhino_utils.niimask2indexpointcloud(outskin_mesh_file)
+            if display_outskin:
+                # Load in scalp surface
+                # And turn the nvoxx x nvoxy x nvoxz volume into a 3 x npoints point cloud
+                smri_headshape_nativeindex = rhino_utils.niimask2indexpointcloud(outskin_mesh_file)
 
-            # Move from native voxel indices to native space coordinates (in mm)
-            smri_headshape_native = rhino_utils.xform_points(mrivoxel_scaledmri_t["trans"], smri_headshape_nativeindex)
+                # Move from native voxel indices to native space coordinates (in mm)
+                smri_headshape_native = rhino_utils.xform_points(mrivoxel_scaledmri_t["trans"], smri_headshape_nativeindex)
 
-            # Move to MEG (device) space
-            smri_headshape_meg = rhino_utils.xform_points(mri_trans["trans"], smri_headshape_native)
+                # Move to MEG (device) space
+                smri_headshape_meg = rhino_utils.xform_points(mri_trans["trans"], smri_headshape_native)
+
+            # -------------------------
+            # Setup inner skull surface
+
+            if display_inskull:
+                # Load in inner skull surface and turn the nvoxx x nvoxy x nvoxz volume into a 3 x npoints point cloud
+                inner_skull_nativeindex = rhino_utils.niimask2indexpointcloud(bet_inskull_mesh_file)
+
+                # Move from native voxel indices to native space coordinates (in mm)
+                inner_skull_native = rhino_utils.xform_points(mrivoxel_scaledmri_t["trans"], inner_skull_nativeindex)
+
+                # Move to MEG (device) space
+                inner_skull_meg = rhino_utils.xform_points(mri_trans["trans"], inner_skull_native)
+
+            # -------------------
 
             plt.figure()
             ax = plt.axes(projection="3d")
@@ -820,6 +1001,12 @@ def coreg_display(
                     smri_headshape_megt = smri_headshape_meg
                     ax.scatter(smri_headshape_megt[0, 0:-1:10], smri_headshape_megt[1, 0:-1:10], smri_headshape_megt[2, 0:-1:10], color=color, marker=marker, s=scale, alpha=alpha)
 
+            if display_inskull:
+                # Inner skull
+                inner_skull_megt = inner_skull_meg
+                color, scale, alpha, marker = (0.5, 0.5, 0.5), 6, 0.1, "."
+                ax.scatter(inner_skull_megt[0, 0:-1:20], inner_skull_megt[1, 0:-1:20], inner_skull_megt[2, 0:-1:20], color=color, marker=marker, s=scale, alpha=alpha)
+
             if display_headshape_pnts:
                 color, scale, alpha, marker = "red", 8, 0.7, "o"
                 if polhemus_headshape_meg is not None and len(polhemus_headshape_meg) > 0:
@@ -846,6 +1033,41 @@ def coreg_display(
                 else:
                     log_or_print("There are no polhemus derived fiducials to plot")
 
+            if display_dipoles:
+
+                if single_dipoles is not None and len(single_dipoles.T) > 0:
+                    color, scale, alpha, marker = "blue", 1, 0.5, "."
+                    single_pntst = single_dipoles.T
+                    ax.scatter(
+                        single_pntst[0, :],
+                        single_pntst[1, :],
+                        single_pntst[2, :],
+                        color=color,
+                        marker=marker,
+                        s=scale,
+                        alpha=alpha,
+                    )
+
+                if midline_points is not None and len(midline_points.T) > 0:
+                    color, scale, alpha, marker = "green", 1, 0.5, "."
+                    midline_pntst = midline_points.T
+                    ax.scatter(
+                        midline_pntst[0, :],
+                        midline_pntst[1, :],
+                        midline_pntst[2, :],
+                        color=color,
+                        marker=marker,
+                        s=scale,
+                        alpha=alpha,
+                    )
+
+                if multi_dipoles is not None and len(multi_dipoles) > 0:
+                    color = "red"
+                    for pp in multi_dipoles:
+                        ax.plot([src_pnts[pp[0], 0], src_pnts[pp[1], 0]], 
+                                    [src_pnts[pp[0], 1], src_pnts[pp[1], 1]], 
+                                    [src_pnts[pp[0], 2], src_pnts[pp[1], 2]], c=color)
+
             if filename is None:
                 plt.show()
             else:
@@ -855,256 +1077,82 @@ def coreg_display(
         else:
             raise ValueError("invalid plot_type.")
 
-
 def bem_display(
     subjects_dir,
     subject,
     plot_type="surf",
-    display_outskin_with_nose=True,
+    display_outskin=True,
+    display_outskin_with_nose=False,
     display_sensors=False,
+    display_sensor_oris=False,
+    display_fiducials=False,
+    display_headshape_pnts=False,    
+    display_inskull=True,    
+    display_dipoles=True,
+    multi_dipoles=None,
+    single_dipoles=None,
+    midline_points=None,
     filename=None,
+    src_index=0,
 ):
-    """Displays the coregistered RHINO scalp surface and inner skull surface.
-
-    Display is done in MEG (device) space (in mm).
-
-    Parameters
-    ----------
-    subjects_dir : string
-        Directory to find RHINO subject dirs in.
-    subject : string
-        Subject name dir to find RHINO files in.
-    plot_type : string
-        Either:
-            'surf' to do a 3D surface plot using surface meshes.
-            'scatter' to do a scatter plot using just point clouds.
-    display_outskin_with_nose : bool
-        Whether to include nose with scalp surface in the display.
-    display_sensors : bool
-        Whether to include sensor locations in the display.
-    filename : str
-        Filename to save display to (as an interactive html). Must have extension .html.
-    """
+    '''
+    3D display of the Boundary Element Model, including
+    source/dipole locations.
+
+    bem_display calls coreg_display, but has different defaults for the inputs.
+
+    See coreg_display help for more.
+    '''
+
+    coreg_display(    
+        subjects_dir,
+        subject,
+        plot_type=plot_type,
+        display_outskin=display_outskin,
+        display_outskin_with_nose=display_outskin_with_nose,
+        display_sensors=display_sensors,
+        display_sensor_oris=display_sensor_oris,
+        display_fiducials=display_fiducials,
+        display_headshape_pnts=display_headshape_pnts,          
+        display_inskull=display_inskull,    
+        display_dipoles=display_dipoles,
+        multi_dipoles=multi_dipoles,
+        single_dipoles=single_dipoles,
+        midline_points=midline_points,
+        filename=filename,
+        src_index=src_index,
+    )
+
+
+def apply_head2mni_xform(subjects_dir, 
+                         subject, 
+                         head_pnts):
+    '''
+    Apply the head to MNI transform to the given points.
+    This function assumes that the points are in the head coordinate system.
+
+    Inputs:
+        subjects_dir: string
+            Directory containing subject data.
+        subject: string
+            Subject identifier.
+        head_pnts: array
+            Points in the head coordinate system in mm.
+
+    Outputs:
+        mni_pnts: array
+            Points in the MNI coordinate system in mm.
+
+    '''
 
-    # Note the jargon used varies for xforms and coord spaces:
-    # MEG (device) -- dev_head_t --> HEAD (polhemus)
-    # HEAD (polhemus)-- head_mri_t (polhemus2native) --> MRI (native)
-    # MRI (native) -- mri_mrivoxel_t (native2nativeindex) --> MRI (native) voxel indices
-    #
-    # RHINO does everthing in mm
-
-    rhino_files = rhino_utils.get_rhino_filenames(subjects_dir, subject)
-    filenames = rhino_files["coreg"]
-
-    bet_outskin_plus_nose_mesh_file = filenames["bet_outskin_plus_nose_mesh_file"]
-    bet_outskin_plus_nose_surf_file = filenames["bet_outskin_plus_nose_surf_file"]
-    bet_outskin_mesh_file = filenames["bet_outskin_mesh_file"]
-    bet_outskin_mesh_vtk_file = filenames["bet_outskin_mesh_vtk_file"]
-    bet_outskin_surf_file = filenames["bet_outskin_surf_file"]
-    bet_inskull_mesh_file = filenames["bet_inskull_mesh_file"]
-    bet_inskull_surf_file = filenames["bet_inskull_surf_file"]
-
-    head_scaledmri_t_file = filenames["head_scaledmri_t_file"]
-    mrivoxel_scaledmri_t_file = filenames["mrivoxel_scaledmri_t_file"]
-
-    info_fif_file = filenames["info_fif_file"]
-
-    if display_outskin_with_nose:
-        outskin_mesh_file = bet_outskin_plus_nose_mesh_file
-        outskin_mesh_4surf_file = bet_outskin_plus_nose_mesh_file
-        outskin_surf_file = bet_outskin_plus_nose_surf_file
-    else:
-        outskin_mesh_file = bet_outskin_mesh_file
-        outskin_mesh_4surf_file = bet_outskin_mesh_vtk_file
-        outskin_surf_file = bet_outskin_surf_file
-
-    fwd_fname = rhino_files["fwd_model"]
-    if Path(fwd_fname).exists():
-        forward = read_forward_solution(fwd_fname)
-        src = forward["src"]
-    else:
-        src = None
-
-    # ------------
-    # Setup xforms
-
-    info = read_info(info_fif_file)
-
-    mrivoxel_scaledmri_t = read_trans(mrivoxel_scaledmri_t_file)
-
-    # get meg to head xform in metres from info
-    head_scaledmri_t = read_trans(head_scaledmri_t_file)
-    dev_head_t, _ = _get_trans(info["dev_head_t"], "meg", "head")
-
-    # Change xform from metres to mm.
-    # Note that MNE xform in fif.info assume metres, whereas we want it
-    # in mm. To change units on an xform, just need to change the translation
-    # part and leave the rotation alone
-    dev_head_t["trans"][0:3, -1] = dev_head_t["trans"][0:3, -1] * 1000
-
-    # We are going to display everything in MEG (device) coord frame in mm
-    meg_trans = Transform("meg", "meg")
-    mri_trans = invert_transform(combine_transforms(dev_head_t, head_scaledmri_t, "meg", "mri"))
-    head_trans = invert_transform(dev_head_t)
-
-    # -----------------
-    # Setup MEG sensors
-
-    if display_sensors:
-        meg_picks = pick_types(info, meg=True, ref_meg=False, exclude=())
-
-        coil_transs = [_loc_to_coil_trans(info["chs"][pick]["loc"]) for pick in meg_picks]
-        coils = _create_meg_coils([info["chs"][pick] for pick in meg_picks], acc="normal")
-
-        meg_rrs, meg_tris = list(), list()
-        offset = 0
-        for coil, coil_trans in zip(coils, coil_transs):
-            rrs, tris = _sensor_shape(coil)
-            rrs = apply_trans(coil_trans, rrs)
-            meg_rrs.append(rrs)
-            meg_tris.append(tris + offset)
-            offset += len(meg_rrs[-1])
-        if len(meg_rrs) == 0:
-            log_or_print("MEG sensors not found. Cannot plot MEG locations.")
-        else:
-            meg_rrs = apply_trans(meg_trans, np.concatenate(meg_rrs, axis=0))
-            meg_tris = np.concatenate(meg_tris, axis=0)
-
-        # convert to mm
-        meg_rrs = meg_rrs * 1000
-
-    # ----------------------------
-    # Setup vol source grid points
-
-    if src is not None:
-        # stored points are in metres, convert to mm
-        src_pnts = src[0]["rr"][src[0]["vertno"], :] * 1000
-
-        # Move from head space to MEG (device) space
-        src_pnts = rhino_utils.xform_points(head_trans["trans"], src_pnts.T).T
-
-        log_or_print("BEM surface: number of dipoles = {}".format(src_pnts.shape[0]))
-
-    # --------
-    # Do plots
-
-    if plot_type == "surf":
-        # Initialize figure
-        renderer = _get_renderer(None, bgcolor=(0.5, 0.5, 0.5), size=(500, 500))
-
-        # Sensors
-        if display_sensors:
-            if len(meg_rrs) > 0:
-                color, alpha = (0.0, 0.25, 0.5), 0.2
-                surf = dict(rr=meg_rrs, tris=meg_tris)
-                renderer.surface(surface=surf, color=color, opacity=alpha, backface_culling=True)
-
-        # sMRI-derived scalp surface
-        rhino_utils.create_freesurfer_mesh_from_bet_surface(
-            infile=outskin_mesh_4surf_file,
-            surf_outfile=outskin_surf_file,
-            nii_mesh_file=outskin_mesh_file,
-            xform_mri_voxel2mri=mrivoxel_scaledmri_t["trans"],
-        )
-
-        coords_native, faces = nib.freesurfer.read_geometry(outskin_surf_file)
-
-        # Move to MEG (device) space
-        coords_meg = rhino_utils.xform_points(mri_trans["trans"], coords_native.T).T
-
-        surf_smri = dict(rr=coords_meg, tris=faces)
-
-        # plot surface
-        renderer.surface(surface=surf_smri, color=(0.85, 0.85, 0.85), opacity=0.3, backface_culling=False)
-
-        # Inner skull surface
-        # Load in surface, this is in mm
-        coords_native, faces = nib.freesurfer.read_geometry(bet_inskull_surf_file)
-
-        # Move to MEG (device) space
-        coords_meg = rhino_utils.xform_points(mri_trans["trans"], coords_native.T).T
-
-        surf_smri = dict(rr=coords_meg, tris=faces)
-
-        # Plot surface
-        renderer.surface(surface=surf_smri, color=(0.25, 0.25, 0.25), opacity=0.25, backface_culling=False)
-
-        # vol source grid points
-        if src is not None and len(src_pnts.T) > 0:
-            color, scale, alpha = (1, 0, 0), 0.001, 1
-            renderer.sphere(center=src_pnts, color=color, scale=scale * 1000, opacity=alpha, backface_culling=True)
-
-        renderer.set_camera(azimuth=90, elevation=90, distance=600, focalpoint=(0.0, 0.0, 0.0))
-
-        # Save or show
-        rhino_utils.save_or_show_renderer(renderer, filename)
-
-    # --------------------------
-    elif plot_type == "scatter":
+    coreg_filenames = get_coreg_filenames(subjects_dir, subject)
+    surfaces_filenames = get_surfaces_filenames(subjects_dir, subject)
 
-        # -------------------
-        # Setup scalp surface
+    head_mri_t = read_trans(coreg_filenames["head_mri_t_file"])
+    mni_mri_t = read_trans(surfaces_filenames["mni_mri_t_file"])
+    head_mni_t = combine_transforms(head_mri_t, invert_transform(mni_mri_t), 'head', 'mni_tal')
 
-        # Load in scalp surface and turn the nvoxx x nvoxy x nvoxz volume into a 3 x npoints point cloud
-        smri_headshape_nativeindex = rhino_utils.niimask2indexpointcloud(outskin_mesh_file)
+    # Apply this xform to the mni fids to get what we call the "sMRI-derived fids" in native space
+    mni_pnts = rhino_utils.xform_points(head_mni_t["trans"], head_pnts.T).T
 
-        # Move from native voxel indices to native space coordinates (in mm)
-        smri_headshape_native = rhino_utils.xform_points(mrivoxel_scaledmri_t["trans"], smri_headshape_nativeindex)
-
-        # Move to MEG (device) space
-        smri_headshape_meg = rhino_utils.xform_points(mri_trans["trans"], smri_headshape_native)
-
-        # -------------------------
-        # Setup inner skull surface
-
-        # Load in inner skull surface and turn the nvoxx x nvoxy x nvoxz volume into a 3 x npoints point cloud
-        inner_skull_nativeindex = rhino_utils.niimask2indexpointcloud(bet_inskull_mesh_file)
-
-        # Move from native voxel indices to native space coordinates (in mm)
-        inner_skull_native = rhino_utils.xform_points(mrivoxel_scaledmri_t["trans"], inner_skull_nativeindex)
-
-        # Move to MEG (device) space
-        inner_skull_meg = rhino_utils.xform_points(mri_trans["trans"], inner_skull_native)
-
-        ax = plt.axes(projection="3d")
-
-        # Sensors
-        if display_sensors:
-            color, scale, alpha, marker = (0.0, 0.25, 0.5), 2, 0.2, "."
-            if len(meg_rrs) > 0:
-                meg_rrst = meg_rrs.T  # do plot in mm
-                ax.scatter(meg_rrst[0, :], meg_rrst[1, :], meg_rrst[2, :], color=color, marker=marker, s=scale, alpha=alpha)
-
-        # Scalp
-        color, scale, alpha, marker = (0.75, 0.75, 0.75), 6, 0.2, "."
-        if len(smri_headshape_meg) > 0:
-            smri_headshape_megt = smri_headshape_meg
-            ax.scatter(smri_headshape_megt[0, 0:-1:20], smri_headshape_megt[1, 0:-1:20], smri_headshape_megt[2, 0:-1:20], color=color, marker=marker, s=scale, alpha=alpha)
-
-        # Inner skull
-        inner_skull_megt = inner_skull_meg
-        color, scale, alpha, marker = (0.5, 0.5, 0.5), 6, 0.2, "."
-        ax.scatter(inner_skull_megt[0, 0:-1:20], inner_skull_megt[1, 0:-1:20], inner_skull_megt[2, 0:-1:20], color=color, marker=marker, s=scale, alpha=alpha)
-
-        # vol source grid points
-        if src is not None and len(src_pnts.T) > 0:
-            color, scale, alpha, marker = (1, 0, 0), 1, 0.5, "."
-            src_pntst = src_pnts.T
-            ax.scatter(
-                src_pntst[0, :],
-                src_pntst[1, :],
-                src_pntst[2, :],
-                color=color,
-                marker=marker,
-                s=scale,
-                alpha=alpha,
-            )
-
-        if filename is None:
-            plt.show()
-        else:
-            log_or_print(f"saving {filename}")
-            plt.savefig(filename)
-            plt.close()
-    else:
-        raise ValueError("invalid plot_type")
+    return mni_pnts
diff --git a/osl_ephys/source_recon/rhino/utils.py b/osl_ephys/source_recon/rhino/utils.py
index db4f5dd..4e205e6 100644
--- a/osl_ephys/source_recon/rhino/utils.py
+++ b/osl_ephys/source_recon/rhino/utils.py
@@ -1516,7 +1516,6 @@ def rotate_pointcloud(points, angle_degrees, axis='x'):
     
     return np.dot(points, rotation_matrix.T)
 
-import numpy as np
 
 def grid_average_downsample(point_cloud, voxel_size):
     """
@@ -1624,4 +1623,5 @@ def replace_headshape(raw, ds_headshape):
     raw_copy = raw.copy()
     raw_copy.set_montage(montage)
 
-    return raw_copy
\ No newline at end of file
+    return raw_copy
+
diff --git a/osl_ephys/source_recon/utils.py b/osl_ephys/source_recon/utils.py
new file mode 100644
index 0000000..d423534
--- /dev/null
+++ b/osl_ephys/source_recon/utils.py
@@ -0,0 +1,58 @@
+"""source recon utilities.
+
+"""
+
+# Authors: Mark Woolrich <mark.woolrich@ohba.ox.ac.uk>
+#          Chetan Gohil <chetan.gohil@psych.ox.ac.uk>
+#          Rob Seymour <rob.seymour@psych.ox.ac.uk>
+#          Mats van Es <mats.vanes@psych.ox.ac.uk>
+
+from mne import label as mne_label
+
+import os
+import os.path as op
+from pathlib import Path
+
+
+def _find_package_file(filename):
+    files_dir = str(Path(__file__).parent) + "/files/"
+    if op.exists(files_dir + filename):
+        return files_dir + filename
+    
+
+def _find_freesurfer_file(filename):
+    avail = mne_label._read_annot_cands(
+        os.path.join(os.environ["SUBJECTS_DIR"], 'fsaverage', 'label')
+    )
+    if filename in avail:
+        filename, hemis = mne_label._get_annot_fname(
+            None, 'fsaverage', 'both', filename, os.environ['SUBJECTS_DIR']
+        )
+        return filename
+
+
+def find_file(filename, freesurfer=False):
+    """Look for a file within the package.
+
+    Parameters
+    ----------
+    filename : str
+        Path to file to look for.
+    freesurfer : bool, optional
+        Should we look in the freesurfer directory?
+
+    Returns
+    -------
+    filename : str
+        Path to file found.
+    """
+    if not op.exists(filename):
+        if freesurfer:
+            filename = _find_freesurfer_file(filename)
+        else:
+            filename = _find_package_file(filename)
+        
+    if filename is None:
+        raise FileNotFoundError(filename)
+
+    return filename
diff --git a/osl_ephys/source_recon/wrappers.py b/osl_ephys/source_recon/wrappers.py
index 1ee6a5e..b1c66a7 100644
--- a/osl_ephys/source_recon/wrappers.py
+++ b/osl_ephys/source_recon/wrappers.py
@@ -661,6 +661,11 @@ def beamform(
     weight_norm="nai",
     pick_ori="max-power-pre-weight-norm",
     reg=0,
+    use_bilateral_pairs=False,
+    bilateral_tol=5,
+    bilateral_tol_midline=None,
+    dipole_patches_file=None,
+    group_patches=False,
     reportdir=None,
 ):
     """Wrapper function for beamforming.
@@ -689,6 +694,45 @@ def beamform(
         Orientation of the dipoles.
     reg : float, optional
         The regularization for the whitened data covariance.
+    use_bilateral_pairs : bool, optional
+        Whether to use bilateral dipole pairs and a bilateral beamformer.
+        Otherwise, a standard single dipole beamformer is used.
+        Only works with volumetric (RHINO) source spaces.
+        Overridden by dipole_patches_file, i.e. if dipole_patches_file
+        is specified, then use_bilateral_pairs is ignored.
+    bilateral_tol : float, optional
+        The maximum distance in mm between two dipoles to be considered
+        a bilateral pair. Only used if use_bilateral_pairs is True.
+        Recommended value is ~gridstep/2 mm.
+    bilateral_tol_midline : float, optional
+        The minimum distance in mm from the midline for a dipole to be
+        considered for bilateral pairing. Only used if use_bilateral_pairs
+        is True. Recommended value is ~gridstep/2 mm.
+        If None, is set to bilateral_tol. 
+    dipole_patches_file : str, optional
+        Path to a 4D nifti binary, mutually exclusive parcellation file, where
+        each volume should be a mask with 0 for the background  
+        and 1 for the parcel and each parcel should be non-overlapping.
+        Dipoles within a parcel will be considered as a group
+        of dipoles that will be beamformed together as a multi-dipole.
+        If None, then no patches are used in the beamformer.
+        File must be in MNI space.
+        This only works with volumetric (RHINO) source spaces.
+        Overrides use_bilateral_pairs.
+        If group_patches=True is also specified,
+        then the dipole_patches_file needs to be
+        accompanied by a text file with the same name but a .xml
+        extension that contains a "group" field
+        with a unique integer value shared by patches that are 
+        to be grouped into a multi-dipole (if group=0 then the patch 
+        is assumed to not be in a group).
+        (see also parcellation.load_parcellation_description).
+        Using a dipole_patches_file only works with volumetric (RHINO)
+        source spaces.
+    group_patches : bool, optional
+        Whether to group the dipoles within each patch together
+        to form a multi-dipole beamformer.
+        Only used when dipole_patches_file is specified.
     reportdir : str, optional
         Path to report directory
     """
@@ -715,7 +759,7 @@ def beamform(
         )
 
     # Pick channels
-    data.pick(chantypes)
+    chantype_data = data.copy().pick(chantypes)
 
     # Create beamforming filters
     log_or_print("beamforming.make_lcmv")
@@ -724,22 +768,37 @@ def beamform(
     filters = beamforming.make_lcmv(
         subjects_dir=outdir,
         subject=subject,
-        data=data,
+        data=chantype_data,
         chantypes=chantypes,
         reg=reg,
         weight_norm=weight_norm,
         pick_ori=pick_ori,
         rank=rank,
         save_filters=True,
+        use_bilateral_pairs=use_bilateral_pairs,
+        bilateral_tol=bilateral_tol,
+        bilateral_tol_midline=bilateral_tol_midline,
+        dipole_patches_file=dipole_patches_file,  
+        group_patches=group_patches,      
     )
 
     # Make plots
     filters_cov_plot, filters_svd_plot = beamforming.make_plots(
-        outdir, subject, filters, data
+        outdir, subject, filters, chantype_data
     )
     filters_cov_plot = filters_cov_plot.replace(f"{outdir}/", "")
     filters_svd_plot = filters_svd_plot.replace(f"{outdir}/", "")
 
+    dipole_locations_plot = f"{subject}/dipole_locations_plot.html"
+    beamforming.plot_dipole_locations(outdir, 
+                                    subject, 
+                                    use_bilateral_pairs, 
+                                    bilateral_tol, 
+                                    bilateral_tol_midline, 
+                                    dipole_patches_file,
+                                    group_patches,
+                                    f"{outdir}/{dipole_locations_plot}")
+    
     if reportdir is not None:
         # Save info for the report
         src_report.add_to_data(
@@ -752,6 +811,7 @@ def beamform(
                 "freq_range": freq_range,
                 "filters_cov_plot": filters_cov_plot,
                 "filters_svd_plot": filters_svd_plot,
+                "dipole_locations_plot": dipole_locations_plot,
             },
         )
 
@@ -896,7 +956,10 @@ def parcellate(
     epoch_file : str
         Path to epoched preprocessed fif file.
     parcellation_file : str
-        Path to the parcellation file to use.
+        If surface_extraction_method='fsl':
+            Path to a 4D nifti file providing the parcellation to use, 
+            with each volume correspond to a parcel.
+            Must be in MNI space.
     method : str
         Method to use in the parcellation.
     orthogonalisation : str, None
@@ -1115,13 +1178,13 @@ def parcellate(
     # Save plots
     parc_psd_plot = f"{subject}/parc/psd.png"
     parcellation.plot_psd(
-            parcel_data,
-            fs=data.info["sfreq"],
-            freq_range=freq_range,
-            parcellation_file=parcellation_file,
-            filename=f"{outdir}/{parc_psd_plot}",
-            freesurfer=surface_extraction_method=='freesurfer',
-        )
+        parcel_data,
+        fs=data.info["sfreq"],
+        freq_range=freq_range,
+        parcellation_file=parcellation_file,
+        filename=f"{outdir}/{parc_psd_plot}",
+        freesurfer=surface_extraction_method=='freesurfer',
+    )
     parc_corr_plot = f"{subject}/parc/corr.png"
     parcellation.plot_correlation(parcel_data, filename=f"{outdir}/{parc_corr_plot}")
 
@@ -1175,6 +1238,11 @@ def beamform_and_parcellate(
     reference_brain="mni",
     extra_chans="stim",
     neighbour_distance=None,
+    use_bilateral_pairs=False,
+    bilateral_tol=5,
+    bilateral_tol_midline=None, 
+    dipole_patches_file=None,  
+    group_patches=False, 
     reportdir=None,
 ):
     """Wrapper function for beamforming and parcellation.
@@ -1231,6 +1299,45 @@ def beamform_and_parcellate(
     neighbour_distance : float, optional
         Distance in mm between parcel centers to consider neighbours 
         for orthogonalisation='local'.
+    use_bilateral_pairs : bool, optional
+        Whether to use bilateral dipole pairs and a bilateral beamformer.
+        Otherwise, a standard single dipole beamformer is used.
+        This only works with volumetric (RHINO) source spaces.
+        Overridden by dipole_patches_file, i.e. if dipole_patches_file
+        is specified, then use_bilateral_pairs is ignored.
+    bilateral_tol : float, optional
+        The maximum distance in mm between two dipoles to be considered
+        a bilateral pair. Only used if use_bilateral_pairs is True.
+        Recommended value is ~gridstep/2 mm.
+    bilateral_tol_midline : float, optional
+        The minimum distance in mm from the midline for a dipole to be
+        considered for bilateral pairing. Only used if use_bilateral_pairs
+        is True. Recommended value is ~gridstep/2 mm.
+        If None, is set to bilateral_tol.
+    dipole_patches_file : str, optional
+        Path to a 4D nifti binary, mutually exclusive parcellation file, where
+        each volume should be a mask with 0 for the background  
+        and 1 for the parcel and each parcel should be non-overlapping.
+        Dipoles within a parcel will be considered as a group
+        of dipoles that will be beamformed together as a multi-dipole.
+        If None, then no patches are used in the beamformer.
+        File must be in MNI space.
+        This only works with volumetric (RHINO) source spaces.
+        Overrides use_bilateral_pairs.
+        If group_patches=True is also specified,
+        then the dipole_patches_file needs to be
+        accompanied by a text file with the same name but a .xml
+        extension that contains a "group" field
+        with a unique integer value shared by patches that are 
+        to be grouped into a multi-dipole (if group=0 then the patch 
+        is assumed to not be in a group).
+        (see also parcellation.load_parcellation_description).
+        Using a dipole_patches_file only works with volumetric (RHINO)
+        source spaces.
+    group_patches : bool, optional
+        Whether to group the dipoles within each patch together
+        to form a multi-dipole beamformer.
+        Only used when dipole_patches_file is specified.
     reportdir : str, optional
         Path to report directory.
     """
@@ -1266,13 +1373,18 @@ def beamform_and_parcellate(
     filters = beamforming.make_lcmv(
         subjects_dir=outdir,
         subject=subject,
-        data=data,
+        data=chantype_data,
         chantypes=chantypes,
         reg=reg,
         weight_norm=weight_norm,
         pick_ori=pick_ori,
         rank=rank,
         save_filters=True,
+        use_bilateral_pairs=use_bilateral_pairs,
+        bilateral_tol=bilateral_tol,
+        bilateral_tol_midline=bilateral_tol_midline,    
+        dipole_patches_file=dipole_patches_file,
+        group_patches=group_patches,
     )
 
     # Make plots
@@ -1282,6 +1394,16 @@ def beamform_and_parcellate(
     filters_cov_plot = filters_cov_plot.replace(f"{outdir}/", "")
     filters_svd_plot = filters_svd_plot.replace(f"{outdir}/", "")
 
+    dipole_locations_plot = f"{subject}/dipole_locations_plot.html"
+    beamforming.plot_dipole_locations(outdir, 
+                                    subject, 
+                                    use_bilateral_pairs, 
+                                    bilateral_tol, 
+                                    bilateral_tol_midline, 
+                                    dipole_patches_file,
+                                    group_patches,
+                                    f"{outdir}/{dipole_locations_plot}")
+
     # Apply beamforming
     bf_data = beamforming.apply_lcmv(chantype_data, filters)
 
@@ -1289,6 +1411,7 @@ def beamform_and_parcellate(
         bf_data = np.transpose([bf.data for bf in bf_data], axes=[1, 2, 0])
     else:
         bf_data = bf_data.data
+
     bf_data_mni, _, coords_mni, _ = beamforming.transform_recon_timeseries(
         subjects_dir=outdir,
         subject=subject,
@@ -1374,6 +1497,7 @@ def beamform_and_parcellate(
                 "freq_range": freq_range,
                 "filters_cov_plot": filters_cov_plot,
                 "filters_svd_plot": filters_svd_plot,
+                "dipole_locations_plot": dipole_locations_plot,
                 "parcellation_file": parcellation_file,
                 "method": method,
                 "reference_brain": reference_brain,