diff --git a/atlas_densities/app/cell_densities.py b/atlas_densities/app/cell_densities.py
index 9d1986a..b54b44d 100644
--- a/atlas_densities/app/cell_densities.py
+++ b/atlas_densities/app/cell_densities.py
@@ -140,8 +140,14 @@ def app(verbose):
     help="Path where to write the output cell density nrrd file."
     "A voxel value is a number of cells per mm^3",
 )
+@click.option(
+    "--root-region-name",
+    type=str,
+    default="root",
+    help="Name of the region in the hierarchy",
+)
 @log_args(L)
-def cell_density(annotation_path, hierarchy_path, nissl_path, output_path):
+def cell_density(annotation_path, hierarchy_path, nissl_path, output_path, root_region_name):
     """Compute and save the overall mouse brain cell density.
 
     The input Nissl stain volume of AIBS is turned into an actual density field complying with
@@ -174,6 +180,7 @@ def cell_density(annotation_path, hierarchy_path, nissl_path, output_path):
         annotation.raw,
         _get_voxel_volume_in_mm3(annotation),
         nissl.raw,
+        root_region_name=root_region_name,
     )
     nissl.with_data(overall_cell_density).save_nrrd(output_path)
 
@@ -380,6 +387,12 @@ def glia_cell_densities(
     help="Path to the directory where to write the output cell density nrrd files."
     " It will be created if it doesn't exist already.",
 )
+@click.option(
+    "--root-region-name",
+    type=str,
+    default="root",
+    help="Name of the region in the hierarchy",
+)
 @log_args(L)
 def inhibitory_and_excitatory_neuron_densities(
     annotation_path,
@@ -389,6 +402,7 @@ def inhibitory_and_excitatory_neuron_densities(
     neuron_density_path,
     inhibitory_neuron_counts_path,
     output_dir,
+    root_region_name,
 ):  # pylint: disable=too-many-arguments
     """Compute and save the inhibitory and excitatory neuron densities.
 
@@ -442,6 +456,7 @@ def inhibitory_and_excitatory_neuron_densities(
         VoxelData.load_nrrd(nrn1_path).raw,
         neuron_density.raw,
         inhibitory_data=inhibitory_data(inhibitory_df),
+        root_region_name=root_region_name,
     )
 
     if not Path(output_dir).exists():
diff --git a/atlas_densities/densities/cell_density.py b/atlas_densities/densities/cell_density.py
index d59d78a..dd1ebc1 100644
--- a/atlas_densities/densities/cell_density.py
+++ b/atlas_densities/densities/cell_density.py
@@ -1,25 +1,19 @@
 """Functions to compute the overall mouse brain cell density.
 """
-
-from typing import Dict
+from __future__ import annotations
 
 import numpy as np
 from atlas_commons.typing import AnnotationT, BoolArray, FloatArray
 from voxcell import RegionMap  # type: ignore
 
+from atlas_densities.densities import utils
 from atlas_densities.densities.cell_counts import cell_counts
-from atlas_densities.densities.utils import (
-    compensate_cell_overlap,
-    get_group_ids,
-    get_region_masks,
-    normalize_intensity,
-)
 
 
 def fix_purkinje_layer_intensity(
-    region_map: "RegionMap",
+    group_ids: dict[str, set[int]],
     annotation: AnnotationT,
-    region_masks: Dict[str, BoolArray],
+    region_masks: dict[str, BoolArray],
     cell_intensity: FloatArray,
 ) -> None:
     """
@@ -29,7 +23,8 @@ def fix_purkinje_layer_intensity(
     The array `cell_intensity` is modified in place.
 
     Args:
-        region_map: object to navigate the mouse brain regions hierarchy.
+        group_ids: a dictionary whose keys are group names and whose values are
+            sets of AIBS structure identifiers.
         annotation: integer array of shape (W, H, D) enclosing the AIBS annotation of
             the whole mouse brain.
         region_masks: A dictionary whose keys are region group names and whose values are
@@ -40,8 +35,8 @@ def fix_purkinje_layer_intensity(
             way that the Purkinje layer has a constant intensity value.
     """
 
-    group_ids = get_group_ids(region_map)
     purkinje_layer_mask = np.isin(annotation, list(group_ids["Purkinje layer"]))
+
     # Force Purkinje Layer regions of the Cerebellum group to have a constant intensity
     # equal to the average intensity of the complement.
     # pylint: disable=fixme
@@ -64,6 +59,7 @@ def compute_cell_density(
     annotation: AnnotationT,
     voxel_volume: float,
     nissl: FloatArray,
+    root_region_name: str,
 ) -> FloatArray:
     """
     Compute the overall cell density based on Nissl staining and cell counts from literature.
@@ -101,12 +97,12 @@ def compute_cell_density(
     """
 
     nissl = np.asarray(nissl, dtype=np.float64)
-    nissl = normalize_intensity(nissl, annotation, threshold_scale_factor=1.0, copy=False)
-    nissl = compensate_cell_overlap(nissl, annotation, gaussian_filter_stdv=-1.0, copy=False)
+    nissl = utils.normalize_intensity(nissl, annotation, threshold_scale_factor=1.0, copy=False)
+    nissl = utils.compensate_cell_overlap(nissl, annotation, gaussian_filter_stdv=-1.0, copy=False)
 
-    group_ids = get_group_ids(region_map)
-    region_masks = get_region_masks(group_ids, annotation)
-    fix_purkinje_layer_intensity(region_map, annotation, region_masks, nissl)
+    group_ids = utils.get_group_ids(region_map, root_region_name=root_region_name)
+    region_masks = utils.get_region_masks(group_ids, annotation)
+    fix_purkinje_layer_intensity(group_ids, annotation, region_masks, nissl)
     non_zero_nissl = nissl > 0
     for group, mask in region_masks.items():
         mask = np.logical_and(non_zero_nissl, mask)
diff --git a/atlas_densities/densities/fitting.py b/atlas_densities/densities/fitting.py
index 09de580..b327a54 100644
--- a/atlas_densities/densities/fitting.py
+++ b/atlas_densities/densities/fitting.py
@@ -620,7 +620,7 @@ def linear_fitting(  # pylint: disable=too-many-arguments
 
     L.info("Getting group names ...")
     # We want group region names to be stable under taking descendants
-    groups = get_group_names(region_map, cleanup_rest=True)
+    groups = get_group_names(region_map, cleanup_rest=True, root_region_name=region_name)
 
     L.info("Computing fitting coefficients ...")
     fitting_coefficients = compute_fitting_coefficients(groups, average_intensities, densities)
diff --git a/atlas_densities/densities/glia_densities.py b/atlas_densities/densities/glia_densities.py
index f8aeb4b..b9d4f37 100644
--- a/atlas_densities/densities/glia_densities.py
+++ b/atlas_densities/densities/glia_densities.py
@@ -1,24 +1,20 @@
 """Functions to compute glia cell densities."""
+from __future__ import annotations
+
 import logging
-from typing import Dict, Set
 
 import numpy as np
 import voxcell
 from atlas_commons.typing import AnnotationT, FloatArray
 
-from atlas_densities.densities.utils import (
-    compensate_cell_overlap,
-    constrain_cell_counts_per_voxel,
-    get_group_ids,
-    normalize_intensity,
-)
+from atlas_densities.densities import utils
 
 L = logging.getLogger(__name__)
 
 
 def compute_glia_cell_counts_per_voxel(  # pylint: disable=too-many-arguments
     glia_cell_count: int,
-    group_ids: Dict[str, Set[int]],
+    region_map: voxcell.RegionMap,
     annotation: AnnotationT,
     glia_intensity: FloatArray,
     cell_counts_per_voxel: FloatArray,
@@ -34,9 +30,7 @@ def compute_glia_cell_counts_per_voxel(  # pylint: disable=too-many-arguments
             cell counts.
     Args:
         glia_cell_count: overall glia cell count found in the scientific literature.
-        group_ids: dictionary whose keys are group names and whose values are
-            sets of AIBS structure ids. These groups are used to identify the voxels
-            with maximum density (fiber tracts) and those of zero density (Purkinje layer).
+        region_map: object to navigate the mouse brain regions hierarchy
         annotation: integer array of shape (W, H, D) enclosing the AIBS annotation of the whole
             mouse brain.
         glia_intensity: float array of shape (W, H, D) with non-negative entries. The input
@@ -53,13 +47,13 @@ def compute_glia_cell_counts_per_voxel(  # pylint: disable=too-many-arguments
         layer).
     """
 
-    fiber_tracts_mask = np.isin(annotation, list(group_ids["Fiber tracts group"]))
+    fiber_tracts_mask = np.isin(annotation, list(utils.get_fiber_tract_ids(region_map)))
     fiber_tracts_free_mask = np.isin(
         annotation,
-        list(group_ids["Purkinje layer"]),
+        list(utils.get_purkinje_layer_ids(region_map)),
     )
 
-    return constrain_cell_counts_per_voxel(
+    return utils.constrain_cell_counts_per_voxel(
         glia_cell_count,
         glia_intensity,
         cell_counts_per_voxel,
@@ -74,11 +68,11 @@ def compute_glia_densities(  # pylint: disable=too-many-arguments
     annotation: AnnotationT,
     voxel_volume: float,
     glia_cell_count: int,
-    glia_intensities: Dict[str, FloatArray],
+    glia_intensities: dict[str, FloatArray],
     cell_density: FloatArray,
-    glia_proportions: Dict[str, str],
+    glia_proportions: dict[str, str],
     copy: bool = False,
-) -> Dict[str, FloatArray]:
+) -> dict[str, FloatArray]:
     """
     Compute the overall glia cell density as well as astrocyte, olgidendrocyte and microglia
     densities.
@@ -131,7 +125,7 @@ def compute_glia_densities(  # pylint: disable=too-many-arguments
         glia_densities[glia_type] = np.asarray(glia_densities[glia_type], dtype=np.float64)
     cell_density = np.asarray(cell_density, dtype=np.float64)
 
-    glia_densities["glia"] = compensate_cell_overlap(
+    glia_densities["glia"] = utils.compensate_cell_overlap(
         np.asarray(glia_densities["glia"], dtype=np.float64),
         annotation,
         gaussian_filter_stdv=-1.0,
@@ -144,19 +138,19 @@ def compute_glia_densities(  # pylint: disable=too-many-arguments
 
     glia_densities["glia"] = compute_glia_cell_counts_per_voxel(
         glia_cell_count,
-        get_group_ids(region_map),
+        region_map,
         annotation,
         glia_densities["glia"],
         cell_density * voxel_volume,
     )
     placed_cells = np.zeros_like(glia_densities["glia"])
     for glia_type in ["astrocyte", "oligodendrocyte"]:
-        glia_densities[glia_type] = normalize_intensity(
+        glia_densities[glia_type] = utils.normalize_intensity(
             np.asarray(glia_densities[glia_type], dtype=np.float64),
             annotation,
             copy=copy,
         )
-        glia_densities[glia_type] = compensate_cell_overlap(
+        glia_densities[glia_type] = utils.compensate_cell_overlap(
             glia_densities[glia_type],
             annotation,
             gaussian_filter_stdv=2.0,
@@ -168,7 +162,7 @@ def compute_glia_densities(  # pylint: disable=too-many-arguments
             glia_type,
             cell_count,
         )
-        glia_densities[glia_type] = constrain_cell_counts_per_voxel(
+        glia_densities[glia_type] = utils.constrain_cell_counts_per_voxel(
             cell_count,
             glia_densities[glia_type],
             glia_densities["glia"] - placed_cells,
diff --git a/atlas_densities/densities/inhibitory_neuron_density.py b/atlas_densities/densities/inhibitory_neuron_density.py
index 41cd568..9b98931 100644
--- a/atlas_densities/densities/inhibitory_neuron_density.py
+++ b/atlas_densities/densities/inhibitory_neuron_density.py
@@ -89,6 +89,7 @@ def compute_inhibitory_neuron_density(  # pylint: disable=too-many-arguments
     neuron_density: FloatArray,
     inhibitory_proportion: Optional[float] = None,
     inhibitory_data: Optional[InhibitoryData] = None,
+    root_region_name: Optional[str] = None,
 ) -> FloatArray:
     """
     Compute the inhibitory neuron density using a prescribed neuron count and the overall neuron
@@ -115,6 +116,7 @@ def compute_inhibitory_neuron_density(  # pylint: disable=too-many-arguments
                 assigning the proportion of ihnibitory neurons in each group named by a key string.
             'neuron_count': the total number of inhibitory neurons (float).
             Used only if `inhibitory_proportion` is None.
+        root_region_name(str): Name of the root region in the hierarchy
 
     Returns:
         float64 array of shape (W, H, D) with non-negative entries.
@@ -138,7 +140,7 @@ def compute_inhibitory_neuron_density(  # pylint: disable=too-many-arguments
                 "Either inhibitory_proportion or inhibitory_data should be provided"
                 ". Both are None."
             )
-        group_ids = get_group_ids(region_map)
+        group_ids = get_group_ids(region_map, root_region_name=root_region_name)
         inhibitory_data["region_masks"] = get_region_masks(group_ids, annotation)
     else:
         inhibitory_data = {
diff --git a/atlas_densities/densities/utils.py b/atlas_densities/densities/utils.py
index aa1151c..bf91776 100644
--- a/atlas_densities/densities/utils.py
+++ b/atlas_densities/densities/utils.py
@@ -1,6 +1,7 @@
 """Utility functions for cell density computation."""
+from __future__ import annotations
 
-from typing import TYPE_CHECKING, Dict, Set
+from typing import TYPE_CHECKING
 from warnings import warn
 
 import numpy as np
@@ -323,7 +324,35 @@ def constrain_cell_counts_per_voxel(  # pylint: disable=too-many-arguments, too-
     return cell_counts
 
 
-def get_group_ids(region_map: "RegionMap", cleanup_rest: bool = False) -> Dict[str, Set[int]]:
+def get_fiber_tract_ids(region_map: "RegionMap") -> set[int]:
+    """
+    Args:
+        region_map: object to navigate the mouse brain regions hierarchy
+    """
+    fiber_tracts_ids = (
+        region_map.find("fiber tracts", attr="name", with_descendants=True)
+        | region_map.find("grooves", attr="name", with_descendants=True)
+        | region_map.find("ventricular systems", attr="name", with_descendants=True)
+        | region_map.find("Basic cell groups and regions", attr="name")
+        | region_map.find("Cerebellum", attr="name")
+    )
+    assert fiber_tracts_ids, "Missing ids in Fiber tracts"
+    return fiber_tracts_ids
+
+
+def get_purkinje_layer_ids(region_map: "RegionMap") -> set[int]:
+    """
+    Args:
+        region_map: object to navigate the mouse brain regions hierarchy
+    """
+    purkinje_layer_ids = region_map.find("@.*Purkinje layer", attr="name", with_descendants=True)
+    assert purkinje_layer_ids, "Missing ids in Purkinje layer"
+    return purkinje_layer_ids
+
+
+def get_group_ids(
+    region_map: "RegionMap", cleanup_rest: bool = False, root_region_name: str | None = None
+) -> dict[str, set[int]]:
     """
     Get AIBS structure ids for several region groups of interest.
 
@@ -341,6 +370,7 @@ def get_group_ids(region_map: "RegionMap", cleanup_rest: bool = False) -> Dict[s
         A dictionary whose keys are region group names and whose values are
         sets of structure identifiers.
     """
+    # pylint: disable=too-many-locals
     cerebellum_group_ids = region_map.find(
         "Cerebellum", attr="name", with_descendants=True
     ) | region_map.find("arbor vitae", attr="name", with_descendants=True)
@@ -349,19 +379,14 @@ def get_group_ids(region_map: "RegionMap", cleanup_rest: bool = False) -> Dict[s
         | region_map.find("Entorhinal area", attr="name", with_descendants=True)
         | region_map.find("Piriform area", attr="name", with_descendants=True)
     )
-    purkinje_layer_ids = region_map.find("@.*Purkinje layer", attr="name", with_descendants=True)
-    fiber_tracts_ids = (
-        region_map.find("fiber tracts", attr="name", with_descendants=True)
-        | region_map.find("grooves", attr="name", with_descendants=True)
-        | region_map.find("ventricular systems", attr="name", with_descendants=True)
-        | region_map.find("Basic cell groups and regions", attr="name")
-        | region_map.find("Cerebellum", attr="name")
-    )
+    purkinje_layer_ids = get_purkinje_layer_ids(region_map)
+    fiber_tracts_ids = get_fiber_tract_ids(region_map)
     cerebellar_cortex_ids = region_map.find("Cerebellar cortex", attr="name", with_descendants=True)
     molecular_layer_ids = cerebellar_cortex_ids & region_map.find(
         "@.*molecular layer", attr="name", with_descendants=True
     )
-    rest_ids = region_map.find("root", attr="name", with_descendants=True)
+    rest_ids = region_map.find(root_region_name, attr="name", with_descendants=True)
+    assert rest_ids, f"Did not find any ids in {root_region_name}"
     rest_ids -= cerebellum_group_ids | isocortex_group_ids
 
     if cleanup_rest:
@@ -372,7 +397,7 @@ def get_group_ids(region_map: "RegionMap", cleanup_rest: bool = False) -> Dict[s
         ascendant_ids |= set(region_map.get(cerebellum_id, attr="id", with_ascendants=True))
         rest_ids -= ascendant_ids
 
-    return {
+    ret = {
         "Cerebellum group": cerebellum_group_ids,
         "Isocortex group": isocortex_group_ids,
         "Fiber tracts group": fiber_tracts_ids,
@@ -381,9 +406,14 @@ def get_group_ids(region_map: "RegionMap", cleanup_rest: bool = False) -> Dict[s
         "Cerebellar cortex": cerebellar_cortex_ids,
         "Rest": rest_ids,
     }
+    for name, ids in ret.items():
+        assert ids, f"Missing ids in {name}"
+    return ret
 
 
-def get_group_names(region_map: "RegionMap", cleanup_rest: bool = False) -> Dict[str, Set[str]]:
+def get_group_names(
+    region_map: "RegionMap", cleanup_rest: bool = False, root_region_name: str | None = None
+) -> dict[str, set[str]]:
     """
     Get AIBS names for regions in several region groups of interest.
 
@@ -399,7 +429,7 @@ def get_group_names(region_map: "RegionMap", cleanup_rest: bool = False) -> Dict
         sets of brain region names.
     """
 
-    group_ids = get_group_ids(region_map, cleanup_rest)
+    group_ids = get_group_ids(region_map, cleanup_rest, root_region_name)
 
     return {
         group_name: {region_map.get(id_, attr="name") for id_ in group_ids[group_name]}
@@ -408,8 +438,8 @@ def get_group_names(region_map: "RegionMap", cleanup_rest: bool = False) -> Dict
 
 
 def get_region_masks(
-    group_ids: Dict[str, Set[int]], annotation: FloatArray
-) -> Dict[str, BoolArray]:
+    group_ids: dict[str, set[int]], annotation: FloatArray
+) -> dict[str, BoolArray]:
     """
     Get the boolean masks of several region groups of interest.
 
@@ -434,7 +464,7 @@ def get_region_masks(
     }
 
 
-def get_aibs_region_names(region_map: "RegionMap") -> Set[str]:
+def get_aibs_region_names(region_map: "RegionMap") -> set[str]:
     """
     Retrieve the names of every region in `region_map`.
 
diff --git a/tests/densities/test_cell_density.py b/tests/densities/test_cell_density.py
index e25d46e..b46ba0f 100644
--- a/tests/densities/test_cell_density.py
+++ b/tests/densities/test_cell_density.py
@@ -25,9 +25,15 @@ def test_compute_cell_density():
     nissl = rng.random(annotation.shape)
     nissl[0][0][0] = 1e-5  # the outside voxels' intensity should be low
 
-    cell_density = tested.compute_cell_density(region_map, annotation, voxel_volume, nissl)
+    cell_density = tested.compute_cell_density(
+        region_map,
+        annotation,
+        voxel_volume,
+        nissl,
+        root_region_name="root",
+    )
     # Each group has a prescribed cell count
-    group_ids = get_group_ids(region_map)
+    group_ids = get_group_ids(region_map, root_region_name="root")
     region_masks = get_region_masks(group_ids, annotation)
     for group, mask in region_masks.items():
         npt.assert_array_almost_equal(
@@ -55,9 +61,10 @@ def test_cell_density_with_soma_correction():
         annotation,
         voxel_volume,
         nissl,
+        root_region_name="root",
     )
     # Each group has a prescribed cell count
-    group_ids = get_group_ids(region_map)
+    group_ids = get_group_ids(region_map, root_region_name="root")
     region_masks = get_region_masks(group_ids, annotation)
     for group, mask in region_masks.items():
         npt.assert_array_almost_equal(
@@ -79,17 +86,23 @@ def test_cell_density_options():
     rng = np.random.default_rng(seed=42)
     nissl = rng.random(annotation.shape)
     nissl[0][0][0] = 1e-5  # the outside voxels' intensity should be low
-    group_ids = get_group_ids(region_map)
+    group_ids = get_group_ids(region_map, root_region_name="root")
     region_masks = get_region_masks(group_ids, annotation)
 
     with patch(
-        "atlas_densities.densities.cell_density.compensate_cell_overlap",
+        "atlas_densities.densities.utils.compensate_cell_overlap",
         return_value=nissl,
     ):
-        actual = tested.compute_cell_density(region_map, annotation, voxel_volume, nissl)
+        actual = tested.compute_cell_density(
+            region_map,
+            annotation,
+            voxel_volume,
+            nissl,
+            root_region_name="root",
+        )
         expected_intensity = nissl.copy()
         tested.fix_purkinje_layer_intensity(
-            region_map, annotation, cell_counts(), expected_intensity
+            group_ids, annotation, cell_counts(), expected_intensity
         )
         for group, mask in region_masks.items():
             expected_intensity[mask] = nissl[mask] * (cell_counts()[group] / np.sum(nissl[mask]))
diff --git a/tests/densities/test_fitting.py b/tests/densities/test_fitting.py
index 8d1e8b5..8b7e09d 100644
--- a/tests/densities/test_fitting.py
+++ b/tests/densities/test_fitting.py
@@ -423,6 +423,25 @@ def get_hierarchy():
                 "name": "Hippocampal formation",
                 "children": [],
             },
+            {
+                "id": 42,
+                "acronym": "PKL",
+                "name": "Purkinje layer",
+                "children": [],
+            },
+            {
+                "id": 38,
+                "acronym": "CCTX",
+                "name": "Cerebellar cortex",
+                "children": [
+                    {
+                        "id": 37,
+                        "acronym": "MLL",
+                        "name": "molecular layer",
+                        "children": [],
+                    },
+                ],
+            },
         ],
     }
 
@@ -494,7 +513,7 @@ def test_linear_fitting():
         )
         warnings_ = [w for w in warnings_ if isinstance(w.message, AtlasDensitiesWarning)]
         # Three warnings for recording NaN coefficients, three warnings for using them
-        assert len(warnings_) == 6
+        assert len(warnings_) == 9
 
     assert np.allclose(fitting_maps["Rest"]["pv+"]["coefficient"], 2.0 / 3.0)
 
diff --git a/tests/densities/test_glia_densities.py b/tests/densities/test_glia_densities.py
index 29b330c..ab7319b 100644
--- a/tests/densities/test_glia_densities.py
+++ b/tests/densities/test_glia_densities.py
@@ -2,7 +2,7 @@
 
 import numpy as np
 import numpy.testing as npt
-from voxcell import RegionMap, VoxelData
+from voxcell import RegionMap
 
 import atlas_densities.densities.glia_densities as tested
 
@@ -10,15 +10,37 @@
 
 
 def test_compute_glia_cell_counts_per_voxel():
-    group_ids = {
-        "Purkinje layer": set({1, 2}),
-        "Fiber tracts group": set({3, 6}),
-    }
+    region_map = RegionMap.from_dict(
+        {
+            "id": 997,
+            "acronym": "root",
+            "name": "root",
+            "children": [
+                {
+                    "id": 1,
+                    "name": "1 Purkinje layer",
+                    "acronym": "pl1",
+                    "children": [
+                        {"id": 2, "name": "2 Purkinje layer", "acronym": "pl2", "children": []}
+                    ],
+                },
+                {
+                    "id": 3,
+                    "name": "fiber tracts",
+                    "acronym": "ft3",
+                    "children": [
+                        {"id": 6, "name": "6 fiber tracts group", "acronym": "tf6", "children": []}
+                    ],
+                },
+            ],
+        }
+    )
+
     annotation = np.array([[[1, 10, 10, 2, 3]]], dtype=np.uint32)
     cell_density = np.array([[[0.1, 0.5, 0.75, 0.1, 1.0]]], dtype=float)
     glia_density = np.array([[[0.15, 0.1, 0.8, 0.2, 0.8]]], dtype=float)
     corrected_glia_density = tested.compute_glia_cell_counts_per_voxel(
-        2, group_ids, annotation, glia_density, cell_density, copy=False
+        2, region_map, annotation, glia_density, cell_density, copy=False
     )
     expected = np.array([[[0.0, 0.25, 0.75, 0.0, 1.0]]], dtype=float)
     npt.assert_allclose(corrected_glia_density, expected, rtol=1e-2)
@@ -29,7 +51,7 @@ def test_compute_glia_cell_counts_per_voxel():
     glia_density = np.array([[[0.15, 0.1, 0.2, 0.8, 0.8]]], dtype=float)
     glia_density_copy = glia_density.copy()
     corrected_glia_density = tested.compute_glia_cell_counts_per_voxel(
-        2, group_ids, annotation, glia_density, cell_density
+        2, region_map, annotation, glia_density, cell_density
     )
     expected = np.array([[[0.0, 1.0 / 3.0, 2.0 / 3.0, 0.0, 1.0]]], dtype=float)
     npt.assert_allclose(corrected_glia_density, expected, rtol=1e-2)
@@ -48,14 +70,11 @@ def test_glia_cell_counts_per_voxel_input():
     glia_cell_count = 25000
     glia_density = rng.random(annotation.shape).reshape(shape)
     glia_density[0][0][0] = 0.0  # cell density outside the brain should be null.
-    group_ids = {
-        "Purkinje layer": set({1, 2, 7, 11, 20, 25, 33, 200, 1000, 31, 16}),
-        "Fiber tracts group": set({3, 6, 14, 56, 62, 88, 279, 2200, 5667, 7668}),
-    }
+    region_map = RegionMap.load_json(Path(TESTS_PATH, "1.json"))
 
     output_glia_density = tested.compute_glia_cell_counts_per_voxel(
         glia_cell_count,
-        group_ids,
+        region_map,
         annotation,
         glia_density,
         cell_density,
diff --git a/tests/densities/test_inhibitory_neuron_density.py b/tests/densities/test_inhibitory_neuron_density.py
index 86dba22..05bd351 100644
--- a/tests/densities/test_inhibitory_neuron_density.py
+++ b/tests/densities/test_inhibitory_neuron_density.py
@@ -97,6 +97,7 @@ def test_compute_inhibitory_neuron_density():
                     data["nrn1"],
                     neuron_density,
                     inhibitory_data=data["inhibitory_data"],
+                    root_region_name="root",
                 )
                 expected = np.array([[[0.0, 0.0, 4.0, 25.0 / 62.0, 99.0 / 62.0]]]) / voxel_volume
                 npt.assert_almost_equal(inhibitory_neuron_density, expected)
@@ -112,6 +113,7 @@ def test_compute_inhibitory_neuron_density_exception():
             np.array([[[1]]], dtype=float),
             np.array([[[1]]], dtype=float),
             np.array([[[1]]], dtype=float),
+            root_region_name="root",
         )
     assert "inhibitory_proportion" in str(error_)
     assert "inhibitory_data" in str(error_)
@@ -160,6 +162,7 @@ def test_compute_inhibitory_density_large_input():
         data["nrn1"],
         data["neuron_density"],
         inhibitory_data=data["inhibitory_data"],
+        root_region_name="root",
     )
 
     assert np.all(inhibitory_neuron_density <= data["neuron_density"])
diff --git a/tests/densities/test_utils.py b/tests/densities/test_utils.py
index 13ff32b..fdea383 100644
--- a/tests/densities/test_utils.py
+++ b/tests/densities/test_utils.py
@@ -101,7 +101,7 @@ def test_compensate_cell_overlap():
 
 def test_get_group_ids():
     region_map = RegionMap.load_json(Path(TESTS_PATH, "1.json"))
-    group_ids = tested.get_group_ids(region_map)
+    group_ids = tested.get_group_ids(region_map, root_region_name="root")
     for ids in group_ids.values():
         assert len(ids) > 0
     assert len(group_ids["Molecular layer"] & group_ids["Purkinje layer"]) == 0
@@ -117,7 +117,7 @@ def test_get_group_ids():
 
 def test_get_region_masks():
     region_map = RegionMap.load_json(Path(TESTS_PATH, "1.json"))
-    group_ids = tested.get_group_ids(region_map)
+    group_ids = tested.get_group_ids(region_map, root_region_name="root")
     annotation_raw = np.arange(27000).reshape(30, 30, 30)
     region_masks = tested.get_region_masks(group_ids, annotation_raw)
     brain_mask = np.logical_or(