diff --git a/.circleci/config.yml b/.circleci/config.yml
index 33e7527233a..3e06ad2ee4f 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -151,8 +151,8 @@ jobs:
     steps:
       - restore_cache:
           keys:
-            - regression-v2-{{ .Revision }}
-            - regression-v2-
+            - regression-v4-{{ .Revision }}
+            - regression-v4-
       - run:
           name: Get truncated BOLD series
           command: |
@@ -175,7 +175,7 @@ jobs:
               echo "Pre-computed masks were cached"
             fi
       - save_cache:
-         key: regression-v2-{{ .Revision }}-{{ epoch }}
+         key: regression-v4-{{ .Revision }}-{{ epoch }}
          paths:
             - /tmp/data
 
@@ -282,13 +282,13 @@ jobs:
 
       - restore_cache:
           keys:
-            - regression-v2-{{ .Revision }}
+            - regression-v4-{{ .Revision }}
       - restore_cache:
           keys:
-            - masks-workdir-v1-{{ .Branch }}-{{epoch}}
-            - masks-workdir-v1-{{ .Branch }}-
-            - masks-workdir-v1-master-
-            - masks-workdir-v1-
+            - masks-workdir-v2-{{ .Branch }}-{{epoch}}
+            - masks-workdir-v2-{{ .Branch }}-
+            - masks-workdir-v2-master-
+            - masks-workdir-v2-
       - run:
           name: Run regression tests on EPI masks
           no_output_timeout: 2h
@@ -305,14 +305,22 @@ jobs:
               coverage run -p --rcfile=setup.cfg \
                      -m pytest --junit-xml=/tmp/masks/reports/regression.xml \
                      niworkflows/func/tests/
-      - save_cache:
-         key: masks-workdir-v1-{{ .Branch }}-{{ epoch }}
-         paths:
-            - /tmp/masks/workdir
+      - run:
+          name: Clear reports folder & delete plot generator cache
+          command: |
+            pushd reports/
+            tar cvfz fmriprep_bold_mask.tar.gz fmriprep_bold_mask/*/*.nii.gz
+            rm -rf /tmp/masks/reports/fmriprep_bold_mask/
+            popd
+            find workdir/ -name "mask_diff_plot" -exec rm -rf {} +
       - store_artifacts:
           path: /tmp/masks/reports
       - store_test_results:
           path: /tmp/masks/reports
+      - save_cache:
+         key: masks-workdir-v2-{{ .Branch }}-{{ epoch }}
+         paths:
+            - /tmp/masks/workdir
 
       - run:
           name: Coverage preparation
@@ -337,15 +345,6 @@ jobs:
             cp /tmp/masks/reports/coverage.xml .
             sed -i "s+/src/niworkflows+/tmp/src/niworkflows+g" coverage.xml
             python -m codecov --file coverage.xml --flags masks -e CIRCLE_JOB
-      - run:
-          name: Package new masks
-          when: always
-          no_output_timeout: 10m
-          working_directory: /tmp/data
-          command: |
-            tar cfz /tmp/masks/fmriprep_bold_mask.tar.gz fmriprep_bold_mask/*/*.nii.gz
-      - store_artifacts:
-          path: /tmp/masks/fmriprep_bold_mask.tar.gz
 
   test_package:
     machine:
diff --git a/niworkflows/func/tests/test_util.py b/niworkflows/func/tests/test_util.py
index 68351be3191..ef59d3f517e 100755
--- a/niworkflows/func/tests/test_util.py
+++ b/niworkflows/func/tests/test_util.py
@@ -85,7 +85,7 @@ def test_masking(input_fname, expected_fname):
         bold_reference_wf.base_dir = str(base_dir)
 
     out_fname = fname_presuffix(
-        basename, suffix="_masks.svg", use_ext=False, newpath=str(newpath)
+        basename, suffix="_mask.svg", use_ext=False, newpath=str(newpath)
     )
     newpath.mkdir(parents=True, exist_ok=True)
 
diff --git a/niworkflows/func/util.py b/niworkflows/func/util.py
index e8b38249e6e..c2c4496268b 100644
--- a/niworkflows/func/util.py
+++ b/niworkflows/func/util.py
@@ -29,6 +29,7 @@
 def init_bold_reference_wf(
     omp_nthreads,
     bold_file=None,
+    brainmask_thresh=0.85,
     pre_mask=False,
     name="bold_reference_wf",
     gen_report=False,
@@ -50,10 +51,16 @@ def init_bold_reference_wf(
 
     Parameters
     ----------
-    bold_file : str
-        BOLD series NIfTI file
     omp_nthreads : int
         Maximum number of threads an individual process may use
+    bold_file : str
+        BOLD series NIfTI file
+    brainmask_thresh: :obj:`float`
+        Lower threshold for the probabilistic brainmask to obtain
+        the final binary mask (default: 0.85).
+    pre_mask : bool
+        Indicates whether the ``pre_mask`` input will be set (and thus, step 1
+        should be skipped).
     name : str
         Name of workflow (default: ``bold_reference_wf``)
     gen_report : bool
@@ -133,7 +140,9 @@ def init_bold_reference_wf(
         EstimateReferenceImage(), name="gen_ref", mem_gb=1
     )  # OE: 128x128x128x50 * 64 / 8 ~ 900MB.
     enhance_and_skullstrip_bold_wf = init_enhance_and_skullstrip_bold_wf(
-        omp_nthreads=omp_nthreads, pre_mask=pre_mask
+        brainmask_thresh=brainmask_thresh,
+        omp_nthreads=omp_nthreads,
+        pre_mask=pre_mask,
     )
 
     calc_dummy_scans = pe.Node(
@@ -188,7 +197,10 @@ def init_bold_reference_wf(
 
 
 def init_enhance_and_skullstrip_bold_wf(
-    name="enhance_and_skullstrip_bold_wf", pre_mask=False, omp_nthreads=1
+    brainmask_thresh=0.5,
+    name="enhance_and_skullstrip_bold_wf",
+    omp_nthreads=1,
+    pre_mask=False,
 ):
     """
     Enhance and run brain extraction on a BOLD EPI image.
@@ -238,13 +250,16 @@ def init_enhance_and_skullstrip_bold_wf(
 
     Parameters
     ----------
+    brainmask_thresh: :obj:`float`
+        Lower threshold for the probabilistic brainmask to obtain
+        the final binary mask (default: 0.5).
     name : str
         Name of workflow (default: ``enhance_and_skullstrip_bold_wf``)
+    omp_nthreads : int
+        number of threads available to parallel nodes
     pre_mask : bool
         Indicates whether the ``pre_mask`` input will be set (and thus, step 1
         should be skipped).
-    omp_nthreads : int
-        number of threads available to parallel nodes
 
     Inputs
     ------
@@ -345,6 +360,8 @@ def init_enhance_and_skullstrip_bold_wf(
     apply_mask = pe.Node(fsl.ApplyMask(), name="apply_mask")
 
     if not pre_mask:
+        from ..interfaces.nibabel import Binarize
+
         bold_template = get_template(
             "MNI152NLin2009cAsym", resolution=2, desc="fMRIPrep", suffix="boldref"
         )
@@ -383,10 +400,22 @@ def init_enhance_and_skullstrip_bold_wf(
         norm.inputs.fixed_image = str(bold_template)
         map_brainmask = pe.Node(
             ApplyTransforms(
-                interpolation="MultiLabel", input_image=str(brain_mask)
+                interpolation="BSpline",
+                float=True,
+                # Use the higher resolution and probseg for numerical stability in rounding
+                input_image=str(
+                    get_template(
+                        "MNI152NLin2009cAsym",
+                        resolution=1,
+                        label="brain",
+                        suffix="probseg",
+                    )
+                ),
             ),
             name="map_brainmask",
         )
+        binarize_mask = pe.Node(Binarize(thresh_low=brainmask_thresh), name="binarize_mask")
+
         # fmt: off
         workflow.connect([
             (inputnode, init_aff, [("in_file", "moving_image")]),
@@ -397,7 +426,8 @@ def init_enhance_and_skullstrip_bold_wf(
                 ("reverse_invert_flags", "invert_transform_flags"),
                 ("reverse_transforms", "transforms"),
             ]),
-            (map_brainmask, pre_dilate, [("output_image", "in_file")]),
+            (map_brainmask, binarize_mask, [("output_image", "in_file")]),
+            (binarize_mask, pre_dilate, [("out_mask", "in_file")]),
         ])
         # fmt: on
     else:
diff --git a/setup.cfg b/setup.cfg
index 6ca1bf5a1cb..382f97a3075 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -40,7 +40,7 @@ install_requires =
     seaborn
     svgutils
     transforms3d
-    templateflow >= 0.4.2
+    templateflow >= 0.6
 test_requires =
     coverage < 5
     pytest >= 4.4