diff --git a/.gitignore b/.gitignore
index 57505380..a181cb67 100644
--- a/.gitignore
+++ b/.gitignore
@@ -60,3 +60,6 @@ docs/_build
 docs/auto_examples
 docs/generated
 docs/colormaps.rst
+
+# Python virtual environment
+.venv
diff --git a/cortex/blender/__init__.py b/cortex/blender/__init__.py
index e4f54020..23e6f375 100644
--- a/cortex/blender/__init__.py
+++ b/cortex/blender/__init__.py
@@ -5,6 +5,7 @@
 from mda_xdrlib import xdrlib
 import tempfile
 import subprocess as sp
+import site
 
 import numpy as np
 
@@ -16,30 +17,76 @@
 default_blender = options.config.get('dependency_paths', 'blender')
 
 _base_imports = """import sys
-sys.path.insert(0, '{path}')
+for site_dir in {site_dirs}:
+    print("Adding python site directory to sys.path:", site_dir)
+    sys.path.insert(0, site_dir)
+
 from mda_xdrlib import xdrlib
 import blendlib
 import bpy.ops
 from bpy import context as C
 from bpy import data as D
-""".format(path=os.path.split(os.path.abspath(__file__))[0])
+""".format(site_dirs=[
+    os.path.split(os.path.abspath(__file__))[0],
+    *site.getsitepackages(),
+])
+
+
+def _wrap_code(code, filename):
+    """
+    Wrap code for running in blender
+
+    Parameters
+    ----------
+    code : str
+        code to run in blender
+    filename : str
+        file path for blender file (must end in ".blend")
+    """
+    wrapped_code = _base_imports
+    if not os.path.exists(filename):
+        wrapped_code += "blendlib.clear_all()\n"
+    wrapped_code += code
+    wrapped_code += "\nbpy.ops.wm.save_mainfile(filepath='{fname}')".format(fname=filename)
+    return wrapped_code
 
-def _call_blender(filename, code, blender_path=default_blender):
-    """Call blender, while running the given code. If the filename doesn't exist, save a new file in that location.
+
+def _call_blender(filename, code=None, background=True, blender_path=default_blender):
+    """
+    Call blender, while running the given code. If the filename doesn't exist, save a new file in that location.
     New files will be initially cleared by deleting all objects.
+
+    Parameters
+    ----------
+    filename : str
+        file path for blender file (must end in ".blend")
+    code : str, optional
+        code to run in blender. If None, blender will be opened without running any code.
+    background : bool, optional
+        If True, blender will be opened in background mode.
+    blender_path : str, optional
+        Path to blender executable. If None, defaults to the path specified in pycortexconfig file.
     """
     with tempfile.NamedTemporaryFile() as tf:
         print("In new named temp file: %s"%tf.name)
-        startcode = _base_imports
-        endcode = "\nbpy.ops.wm.save_mainfile(filepath='{fname}')".format(fname=filename)
-        cmd = "{blender_path} -b {fname} -P {tfname}".format(blender_path=blender_path, fname=filename, tfname=tf.name)
-        if not os.path.exists(filename):
-            startcode += "blendlib.clear_all()\n"
-            cmd = "{blender_path} -b -P {tfname}".format(blender_path=blender_path, tfname=tf.name)
-        else:
+
+        # Backup
+        if os.path.exists(filename):
             _legacy_blender_backup(filename, blender_path=blender_path)
-        tf.write((startcode+code+endcode).encode())
-        tf.flush()
+
+        # Construct command
+        cmd = blender_path
+        if background:
+            cmd += " -b"
+        if os.path.exists(filename):
+            cmd += " " + filename
+        if code is not None:
+            wrapped_code = _wrap_code(code, filename)
+            tf.write(wrapped_code.encode())
+            tf.flush()
+            cmd += " -P {tfname}".format(tfname=tf.name)
+
+        print(f"Calling blender:\n    {cmd}")
         sp.check_call([w.encode() for w in shlex.split(cmd)],)
 
 
@@ -97,7 +144,7 @@ def _legacy_blender_backup(fname, blender_path=default_blender):
                 shutil.copy(fname, fname_bkup)
 
 
-def add_cutdata(fname, braindata, name="retinotopy", projection="nearest", mesh="hemi", blender_path=default_blender):
+def add_cutdata(fname, braindata, name="retinotopy", projection="nearest", mesh="hemi", blender_path=None):
     """Add data as vertex colors to blender mesh
     
     Useful to add localizer data for help in placing flatmap cuts
@@ -117,6 +164,8 @@ def add_cutdata(fname, braindata, name="retinotopy", projection="nearest", mesh=
     mesh : string
         ...
     """
+    blender_path = blender_path or default_blender
+
     if isinstance(braindata, dataset.Dataset):
         for view_name, data in braindata.views.items():
             add_cutdata(fname, data, name=view_name, projection=projection, mesh=mesh)
@@ -162,10 +211,12 @@ def add_cutdata(fname, braindata, name="retinotopy", projection="nearest", mesh=
     return 
 
 
-def gii_cut(fname, subject, hemi, blender_path=default_blender):
+def gii_cut(fname, subject, hemi, blender_path=None):
     '''
     Add gifti surface to blender
     '''
+    blender_path = blender_path or default_blender
+
     from ..database import db
     hemis = dict(lh='left',
                  rh='right')
@@ -194,16 +245,24 @@ def gii_cut(fname, subject, hemi, blender_path=default_blender):
         _call_blender(fname, code, blender_path=blender_path)
 
 
-def fs_cut(fname, subject, hemi, freesurfer_subject_dir=None, blender_path=default_blender):
-    """Cut freesurfer surface using blender interface
+def fs_cut_init(fname, subject, hemi, freesurfer_subject_dir=None, blender_path=None):
+    """Initialize a blender object from a freesurfer volume.
 
     Parameters
     ----------
     fname : str
         file path for new .blend file (must end in ".blend")
-
-    if `freesurfer_subject_dir` is None, it defaults to SUBJECTS_DIR environment variable
+    subject : str
+        subject name
+    hemi : str
+        hemisphere name (lh or rh)
+    freesurfer_subject_dir : str
+        path to freesurfer subject directory. If None, it defaults to SUBJECTS_DIR environment variable
+    blender_path : str
+        path to blender executable. If None, it defaults to the path specified in pycortexconfig file.
     """
+    blender_path = blender_path or default_blender
+    
     wpts, polys, curv = freesurfer.get_surf(subject, hemi, 'smoothwm', freesurfer_subject_dir=freesurfer_subject_dir)
     ipts, _, _ = freesurfer.get_surf(subject, hemi, 'inflated', freesurfer_subject_dir=freesurfer_subject_dir)
     rcurv = np.clip(((-curv + .6) / 1.2), 0, 1)
@@ -225,8 +284,29 @@ def fs_cut(fname, subject, hemi, freesurfer_subject_dir=None, blender_path=defau
         """.format(tfname=tf.name)
         _call_blender(fname, code, blender_path=blender_path)
 
+
+def fs_cut_open(fname, blender_path=None):
+    """Open a blender file in blender for the manual cut
+
+    Parameters
+    ----------
+    fname : str
+        file path for blender file (must end in ".blend")
+    blender_path : str
+        path to blender executable. If None, it defaults to the path specified in pycortexconfig file.
+    """
+    blender_path = blender_path or default_blender
+    
+    _call_blender(fname, background=False, blender_path=blender_path)
+
+
 def write_patch(bname, pname, mesh="hemi", blender_path=default_blender):
-    """Write out the mesh 'mesh' in the blender file 'bname' into patch file 'pname'
+    """Deprecated: please use write_volume_patch instead"""
+    return write_volume_patch(bname, pname, "hemi", mesh, blender_path)
+
+
+def write_volume_patch(bname, pname, hemi, mesh="hemi", blender_path=None):
+    """Write volume patch in freesurfer format.
     This is a necessary step for flattening the surface in freesurfer
 
     Parameters
@@ -235,11 +315,18 @@ def write_patch(bname, pname, mesh="hemi", blender_path=default_blender):
         blender file name that contains the mesh
     pname : str
         name of patch file to be saved
+    hemi : str
+        hemisphere name (lh or rh)
     mesh : str
         name of mesh in blender file
+    blender_path : str, optional
+        path to blender executable. If None, it defaults to the path specified in pycortexconfig file.
     """
+    blender_path = blender_path or default_blender
+
     p = xdrlib.Packer()
     p.pack_string(pname.encode())
+    p.pack_string(hemi.encode())
     p.pack_string(mesh.encode())
     with tempfile.NamedTemporaryFile() as tf:
         tf.write(p.get_buffer())
@@ -247,8 +334,49 @@ def write_patch(bname, pname, mesh="hemi", blender_path=default_blender):
         code = """with open('{tfname}', 'rb') as fp:
             u = xdrlib.Unpacker(fp.read())
             pname = u.unpack_string().decode('utf-8')
+            hemi = u.unpack_string().decode('utf-8')
             mesh = u.unpack_string().decode('utf-8')
-            blendlib.save_patch(pname, mesh)
+            blendlib.write_volume_patch(pname, hemi, mesh)
         """.format(tfname=tf.name)
         _call_blender(bname, code, blender_path=blender_path)
+    return True
+
+def write_flat_patch(bname, pname, hemi, mesh="hemi", method="MINIMUM_STRETCH", blender_path=None):
+    """Write flat patch in freesurfer format.
+    This is a necessary step for flattening the surface in freesurfer
 
+    Parameters
+    ----------
+    bname : str
+        blender file name that contains the mesh
+    pname : str
+        name of patch file to be saved
+    hemi : str
+        hemisphere name (lh or rh)
+    mesh : str
+        name of mesh in blender file
+    method : str
+        method to use for UV unwrap. One of 'CONFORMAL', 'ANGLE_BASED', 'MINIMUM_STRETCH'.
+    blender_path : str, optional
+        path to blender executable. If None, it defaults to the path specified in pycortexconfig file.
+    """
+    blender_path = blender_path or default_blender
+
+    p = xdrlib.Packer()
+    p.pack_string(pname.encode())
+    p.pack_string(hemi.encode())
+    p.pack_string(mesh.encode())
+    p.pack_string(method.encode())
+    with tempfile.NamedTemporaryFile() as tf:
+        tf.write(p.get_buffer())
+        tf.flush()
+        code = """with open('{tfname}', 'rb') as fp:
+            u = xdrlib.Unpacker(fp.read())
+            pname = u.unpack_string().decode('utf-8')
+            hemi = u.unpack_string().decode('utf-8')
+            mesh = u.unpack_string().decode('utf-8')
+            method = u.unpack_string().decode('utf-8')
+            blendlib.write_flat_patch(pname, hemi, mesh, method)
+        """.format(tfname=tf.name)
+        _call_blender(bname, code, blender_path=blender_path)
+    return True
\ No newline at end of file
diff --git a/cortex/blender/blendlib.py b/cortex/blender/blendlib.py
index c81b1b86..22ad05b1 100644
--- a/cortex/blender/blendlib.py
+++ b/cortex/blender/blendlib.py
@@ -1,9 +1,14 @@
-"""This module is intended to be imported directly by blender.
-It provides utility functions for adding meshes and saving them to communicate with the rest of pycortex
+"""
+This module is intended to be imported directly by blender.
+It provides utility functions for adding meshes and saving them to communicate with the rest of pycortex.
+
+Read more about Blender Python API here: https://docs.blender.org/api/current/index.html.
 """
 import struct
 from mda_xdrlib import xdrlib
 import tempfile
+import time
+import math
 
 import bpy.ops
 from bpy import context as C
@@ -131,7 +136,7 @@ def add_shapekey(shape, name=None):
         key.data[i].co = shape[i]
     return key
 
-def write_patch(filename, pts, edges=None):
+def _write_patch(filename, pts, edges=None):
     """Writes a patch file that is readable by freesurfer.
     
     Parameters
@@ -154,8 +159,56 @@ def write_patch(filename, pts, edges=None):
                 fp.write(struct.pack('>i3f', -i-1, *pt))
             else:
                 fp.write(struct.pack('>i3f', i+1, *pt))
+    print("Wrote freesurfer patch to %s"%filename)
+
+def _circularize_uv_coords(pts, u_min, u_max, v_min, v_max):
+    """Transform UV coordinates into a circular shape while preserving relative positions.
+    
+    Parameters
+    ----------
+    pts : dict
+        Dictionary mapping vertex indices to (u, v, z) coordinates
+    u_min, u_max, v_min, v_max : float
+        Original bounds of the UV coordinates
+        
+    Returns
+    -------
+    dict
+        Dictionary mapping vertex indices to new (u, v, z) coordinates
+    """
+    # Convert to normalized coordinates in [-1, 1] range
+    u_center = (u_max + u_min) / 2
+    v_center = (v_max + v_min) / 2
+    u_scale = (u_max - u_min) / 2
+    v_scale = (v_max - v_min) / 2
+    
+    new_pts = {}
+    for idx, (u, v, z) in pts.items():
+        # Normalize coordinates
+        u_norm = (u - u_center) / u_scale
+        v_norm = (v - v_center) / v_scale
+        
+        # Convert to polar coordinates
+        r = math.sqrt(u_norm**2 + v_norm**2)
+        theta = math.atan2(v_norm, u_norm)
+        
+        # Normalize radius to create perfect circle
+        # Use square root to preserve area/density
+        r = math.sqrt(r)
+        
+        # Convert back to Cartesian coordinates
+        u_new = r * math.cos(theta)
+        v_new = r * math.sin(theta)
+        
+        # Scale back to original range
+        u_new = u_new * u_scale + u_center
+        v_new = v_new * v_scale + v_center
+        
+        new_pts[idx] = (u_new, v_new, z)
+    
+    return new_pts
 
-def _get_pts_edges(mesh):
+def _get_geometry(mesh, hemi, flatten, method=None):
     """Function called within blender to get non-cut vertices & edges
 
     Operates on a mesh object within an open instance of blender. 
@@ -164,10 +217,26 @@ def _get_pts_edges(mesh):
     ----------
     mesh : str
         name of mesh to cut
+    hemi : str
+        hemisphere name (lh or rh)
+    flatten : bool
+        if True, returns flattened coordinates using UV unwrap
+    method : str
+        method to use for UV unwrap. One of 'CONFORMAL', 'ANGLE_BASED', 'MINIMUM_STRETCH'.
+
+    Returns
+    -------
+    verts : set
+        set of vertex indices
+    pts : list
+        list of (vertex_index, flattened_coordinates) tuples
+    edges : set
+        set of edge vertex indices
     """
     if isinstance(mesh, str):
         mesh = D.meshes[mesh]
 
+    # Collect edge vertex indices
     bpy.ops.object.mode_set(mode='OBJECT')
     bpy.ops.object.mode_set(mode='EDIT')
     bpy.ops.mesh.select_all(action='DESELECT')
@@ -175,30 +244,33 @@ def _get_pts_edges(mesh):
     bpy.ops.mesh.select_non_manifold()
     bpy.ops.object.mode_set(mode='OBJECT')
 
-    mwall_edge = set()
+    edge_vertex_idxs = set()  # Medial wall in standard case
     for edge in mesh.edges:
         if edge.select:
-            mwall_edge.add(edge.vertices[0])
-            mwall_edge.add(edge.vertices[1])
+            edge_vertex_idxs.add(edge.vertices[0])
+            edge_vertex_idxs.add(edge.vertices[1])
 
+    # Collect seam vertex indices & select seams
     bpy.ops.object.mode_set(mode='EDIT') 
     C.tool_settings.mesh_select_mode = True, False, False
     bpy.ops.mesh.select_all(action='DESELECT')
     bpy.ops.object.mode_set(mode='OBJECT')
-    seam = set()
+    seam_vertex_idxs = set()
     for edge in mesh.edges:
         if edge.use_seam:
-            seam.add(edge.vertices[0])
-            seam.add(edge.vertices[1])
+            seam_vertex_idxs.add(edge.vertices[0])
+            seam_vertex_idxs.add(edge.vertices[1])
             edge.select = True
 
+    # Expand seam selection & collect expanded vertex indices
     bpy.ops.object.mode_set(mode='EDIT') 
     bpy.ops.mesh.select_more()
     bpy.ops.object.mode_set(mode='OBJECT')
-    smore = set()
+    expanded_seam_vertex_idxs = set()
     for i, vert in enumerate(mesh.vertices):
         if vert.select:
-            smore.add(i)
+            expanded_seam_vertex_idxs.add(i)
+
     # Leave cuts (+ area around them) selected.
     # Uncomment the next lines to revert to previous behavior
     # (deselecting everything)
@@ -206,26 +278,93 @@ def _get_pts_edges(mesh):
     # bpy.ops.mesh.select_all(action='DESELECT')
     # bpy.ops.object.mode_set(mode='OBJECT')
 
-    fverts = set()
-    if hasattr(mesh, "polygons"):
-        faces = mesh.polygons
-    else:
-        faces = mesh.faces
-    for face in faces:
-        fverts.add(face.vertices[0])
-        fverts.add(face.vertices[1])
-        fverts.add(face.vertices[2])
-
-    print("exported %d faces"%len(fverts))
-    edges = mwall_edge | (smore - seam)
-    verts = fverts - seam
+    face_vertices = set()
+    for face in getattr(mesh, "polygons", getattr(mesh, "faces", None)):
+        face_vertices.add(face.vertices[0])
+        face_vertices.add(face.vertices[1])
+        face_vertices.add(face.vertices[2])
+
+    verts = face_vertices - seam_vertex_idxs
     pts = [(v, D.shape_keys['Key'].key_blocks['inflated'].data[v].co) for v in verts]
+    edges = edge_vertex_idxs | (expanded_seam_vertex_idxs - seam_vertex_idxs)
+
+    if flatten:
+        # Scales
+        u_coords, v_coords = [u for _, (u, _, _) in pts], [v for _, (_, v, _) in pts]
+        u_min, u_max, v_min, v_max = min(u_coords), max(u_coords), min(v_coords), max(v_coords)
+        print("u_min: %f, u_max: %f, v_min: %f, v_max: %f"%(u_min, u_max, v_min, v_max))
+
+        if not mesh.uv_layers:
+            mesh.uv_layers.new(name="FlattenUV")
+        
+        print("UV unwrapping mesh with method %s (may take a few minutes)..."%method)
+        start = time.time()
+        bpy.ops.object.mode_set(mode='EDIT')
+        bpy.ops.mesh.select_all(action='SELECT')
+        bpy.ops.uv.unwrap(method=method, margin=0.001)
+        bpy.ops.object.mode_set(mode='OBJECT')
+        end = time.time()
+        print("UV unwrapping mesh took %.1f seconds" % (end - start))
+        
+        print("Collecting coordinates for %d verts..."%len(verts))
+        start = time.time()
+        pts = {}
+        for loop in mesh.loops:
+            if loop.vertex_index in verts:
+                coords2d = mesh.uv_layers.active.data[loop.index].uv
+                u_scaled = coords2d[0] * (u_max - u_min) + u_min
+                v_scaled = coords2d[1] * (v_max - v_min) + v_min
+                
+                if hemi == "rh":
+                    # Rotate 180 degrees clockwise
+                    u_center = (u_max + u_min) / 2
+                    v_center = (v_max + v_min) / 2
+                    u_rotated = 2 * u_center - u_scaled
+                    v_rotated = 2 * v_center - v_scaled
+                    pts[loop.vertex_index] = (u_rotated, v_rotated, 0.0)
+                else:
+                    pts[loop.vertex_index] = (u_scaled, v_scaled, 0.0)
+        
+        # Circularize the UV coordinates
+        print("Circularizing UV coordinates...")
+        pts = _circularize_uv_coords(pts, u_min, u_max, v_min, v_max)
+        
+        pts = sorted(pts.items(), key=lambda x: x[0])
+        end = time.time()
+        print("Collecting coordinates took %.1f seconds" % (end - start))
+
+    print("Collected geometry. verts: %d, pts: %d, edges: %d"%(len(verts), len(pts), len(edges)))
     return verts, pts, edges
 
-def save_patch(fname, mesh='hemi'):
-    """Saves patch to file that can be read by freesurfer"""
-    verts, pts, edges = _get_pts_edges(mesh)
-    write_patch(fname, pts, edges)
+
+def save_patch(fname, mesh="hemi"):
+    """Deprecated: please use write_volume_patch instead"""
+    return write_volume_patch(fname, "lh", mesh)
+
+
+def write_volume_patch(fname, hemi, mesh="hemi"):
+    """Write mesh patch in freesurfer format"""
+    _, pts, edges = _get_geometry(mesh, hemi, flatten=False)
+    _write_patch(fname, pts, edges)
+
+
+def write_flat_patch(fname, hemi, mesh="hemi", method="MINIMUM_STRETCH"):
+    """Write flat patch in freesurfer format
+    
+    Parameters
+    ----------
+    fname : str
+        Output filename
+    hemi : str
+        hemisphere name (lh or rh)
+    mesh : str
+        Name of the mesh to flatten
+    method : str
+        UV unwrapping method to use
+    """
+    _, pts, edges = _get_geometry(mesh, hemi, flatten=True, method=method)
+    _write_patch(fname, pts, edges)
+
 
 def read_xdr(filename):
     with open(filename, "rb") as fp:
diff --git a/cortex/freesurfer.py b/cortex/freesurfer.py
index 4d21d901..f384b0ff 100644
--- a/cortex/freesurfer.py
+++ b/cortex/freesurfer.py
@@ -283,6 +283,9 @@ def import_flat(fs_subject, patch, hemis=['lh', 'rh'], cx_subject=None,
         List of hemispheres to import. Defaults to both hemispheres.
     cx_subject : str
         Pycortex subject name
+    flat_type : str
+        Type of flatmap to import. Defaults to 'freesurfer'.
+        Can be 'freesurfer', 'slim', or 'blender'.
     freesurfer_subject_dir : str
         directory for freesurfer subjects. None defaults to environment variable
         $SUBJECTS_DIR
@@ -308,15 +311,19 @@ def import_flat(fs_subject, patch, hemis=['lh', 'rh'], cx_subject=None,
 
     from . import formats
     for hemi in hemis:
-        if flat_type == 'freesurfer':
-            pts, polys, _ = get_surf(fs_subject, hemi, "patch", patch+".flat", freesurfer_subject_dir=freesurfer_subject_dir)
-            # Reorder axes: X, Y, Z instead of Y, X, Z
-            flat = pts[:, [1, 0, 2]]
-            # Flip Y axis upside down
-            flat[:, 1] = -flat[:, 1]
+        if flat_type in ['freesurfer', 'blender']:
+            surf_path = (patch + ".flat") if (flat_type == 'freesurfer') else (patch + ".flat.blender")
+            pts, polys, _ = get_surf(fs_subject, hemi, "patch", surf_path, freesurfer_subject_dir=freesurfer_subject_dir)
+
+            if flat_type == 'freesurfer':
+                # Reorder axes: X, Y, Z instead of Y, X, Z
+                flat = pts[:, [1, 0, 2]]
+                # Flip Y axis upside down
+                flat[:, 1] = -flat[:, 1]
+            else:
+                flat = pts
         elif flat_type == 'slim':
-            flat_file = get_paths(fs_subject, hemi, type='slim',
-                                  freesurfer_subject_dir=freesurfer_subject_dir)
+            flat_file = get_paths(fs_subject, hemi, type='slim', freesurfer_subject_dir=freesurfer_subject_dir)
             flat_file = flat_file.format(name=patch + ".flat")
             flat, polys = formats.read_obj(flat_file)
 
diff --git a/cortex/segment.py b/cortex/segment.py
index da3ea925..0981d313 100644
--- a/cortex/segment.py
+++ b/cortex/segment.py
@@ -6,8 +6,8 @@
 import warnings
 import numpy as np
 import subprocess as sp
-from builtins import input
 import multiprocessing as mp
+from builtins import input
 
 from . import formats
 from . import blender
@@ -18,7 +18,8 @@
 from .freesurfer import autorecon as run_freesurfer_recon
 from .freesurfer import import_subj as import_freesurfer_subject
 
-slim_path = options.config.get('dependency_paths', 'slim')
+default_blender_path = options.config.get('dependency_paths', 'blender')
+default_slim_path = options.config.get('dependency_paths', 'slim')
 
 
 def init_subject(subject, filenames, do_import_subject=False, **kwargs):
@@ -125,7 +126,8 @@ def edit_segmentation(subject,
 
 def cut_surface(cx_subject, hemi, name='flatten', fs_subject=None, data=None,
                 freesurfer_subject_dir=None, flatten_with='freesurfer', 
-                do_import_subject=True, blender_cmd=None, **kwargs):
+                method=None, do_import_subject=True, blender_path=default_blender_path,
+                recache=True, auto_overwrite=False, **kwargs):
     """Initializes an interface to cut the segmented surface for flatmapping.
     This function creates or opens a blend file in your filestore which allows
     surfaces to be cut along hand-defined seams. Blender will automatically
@@ -154,7 +156,7 @@ def cut_surface(cx_subject, hemi, name='flatten', fs_subject=None, data=None,
         Name of Freesurfer subject directory. None defaults to SUBJECTS_DIR
         environment variable
     flatten_with : str
-        'freesurfer' or 'SLIM' - 'freesurfer' (default) uses freesurfer's 
+        One of 'freesurfer','SLIM', or 'blender' - 'freesurfer' (default) uses freesurfer's 
         `mris_flatten` function to flatten the cut surface. 'SLIM' uses
         the SLIM algorithm, which takes much less time but tends to leave
         more distortions in the flatmap. SLIM is an optional dependency, and 
@@ -162,84 +164,125 @@ def cut_surface(cx_subject, hemi, name='flatten', fs_subject=None, data=None,
         (https://github.com/MichaelRabinovich/Scalable-Locally-Injective-Mappings) 
         to your computer and set the slim dependency path in your pycortex config 
         file to point to </path/to/your/slim/install>/ReweightedARAP
+    method : str
+        method to use for UV unwrap. When using Blender, it must be present and
+        can be one of 'CONFORMAL', 'ANGLE_BASED', 'MINIMUM_STRETCH'.
     do_import_subject : bool
         set option to automatically import flatmaps when both are completed 
         (if set to false, you must import later with `cortex.freesurfer.import_flat()`)
+    blender_path : str
+        Path to blender executable. If None, defaults to path specified in pycortexconfig file.
+    recache : boolean
+        Whether or not to recache intermediate files. Takes longer to plot this way, potentially
+        resolves some errors. Useful if you've made changes to the alignment
+    auto_overwrite : bool
+        Whether to overwrite existing flatmaps. If True, the flatmap will be
+        overwritten without asking for confirmation.
     """
+
     if fs_subject is None:
         fs_subject = cx_subject
-    opts = "[hemi=%s,name=%s]"%(hemi, name)
-    fname = db.get_paths(cx_subject)['anats'].format(type='cutsurf', opts=opts, ext='blend')
+
     # Double-check that fiducial and inflated vertex counts match
-    # (these may not match if a subject is initially imported from freesurfer to pycortex, 
+    # (these may not match if a subject is initially imported from freesurfer to pycortex,
     # and then edited further for a better segmentation and not re-imported)
-    ipt, ipoly, inrm = freesurfer.get_surf(fs_subject, hemi, 'inflated')
-    fpt, fpoly, fnrm = freesurfer.get_surf(fs_subject, hemi, 'fiducial')
+    ipt, ipoly, inrm = freesurfer.get_surf(fs_subject, hemi, "inflated")
+    fpt, fpoly, fnrm = freesurfer.get_surf(fs_subject, hemi, "fiducial")
     if ipt.shape[0] != fpt.shape[0]:
-        raise ValueError("Please re-import subject - fiducial and inflated vertex counts don't match!")
+        raise ValueError(
+            "Please re-import subject - fiducial and inflated vertex counts don't match!"
+        )
     else:
-        print('Vert check ok!')
-    if not os.path.exists(fname):
-        blender.fs_cut(fname, fs_subject, hemi, freesurfer_subject_dir)
+        print("Vert check ok!")
+
+    # Create blender file with cuts
+    opts = "[hemi=%s,name=%s]" % (hemi, name)
+    fname = db.get_paths(cx_subject)["anats"].format(
+        type="cutsurf", opts=opts, ext="blend"
+    )
+    if not os.path.exists(fname) or recache:
+        if os.path.exists(fname):
+            os.remove(fname)
+        print("Initializing blender file %s..."%fname)
+        blender.fs_cut_init(fname, fs_subject, hemi, freesurfer_subject_dir, blender_path=blender_path)
+
     # Add localizer data to facilitate cutting
     if data is not None:
-        if isinstance(data, list):
-            for d in data:
-                blender.add_cutdata(fname, d, name=d.description)
-        else:
-            blender.add_cutdata(fname, data, name=data.description)
-    if blender_cmd is None:
-        blender_cmd = options.config.get('dependency_paths', 'blender')
-    # May be redundant after blender.fs_cut above...
-    if os.path.exists(fname):
-        blender._legacy_blender_backup(fname, blender_path=blender_cmd)
-    sp.call([blender_cmd, fname])
-    patchpath = freesurfer.get_paths(fs_subject, hemi,
-                                     freesurfer_subject_dir=freesurfer_subject_dir)
-    patchpath = patchpath.format(name=name)
-    blender.write_patch(fname, patchpath, blender_path=blender_cmd)
-    if flatten_with == 'freesurfer':
-        done = freesurfer.flatten(fs_subject, hemi, patch=name,
-                           freesurfer_subject_dir=freesurfer_subject_dir,
-                           **kwargs)
-        if not done:
-            # If flattening is aborted, skip the rest of this function
-            # (Do not attempt to import completed flatmaps)
-            return
-        if do_import_subject:
-            # Check to see if both hemispheres have been flattened
-            other = freesurfer.get_paths(fs_subject, "lh" if hemi == "rh" else "rh",
-                                         freesurfer_subject_dir=freesurfer_subject_dir)
-            other = other.format(name=name+".flat")
-            # If so, go ahead and import subject
-            if os.path.exists(other):
-                freesurfer.import_flat(fs_subject, name, cx_subject=cx_subject,
-                            flat_type='freesurfer',
-                            freesurfer_subject_dir=freesurfer_subject_dir)
+        data = data if isinstance(data, list) else [data]
+        for d in data:
+            blender.add_cutdata(fname, d, name=d.description, blender_path=blender_path)
+            
+    # Open blender for user to manually do the cuts
+    print("Opening blender file %s..."%fname)
+    blender.fs_cut_open(fname, blender_path=blender_path)
+    
+    # Generate 3D base freesurfer patch to flatten
+    base_patch_path = freesurfer.get_paths(
+        fs_subject, hemi, freesurfer_subject_dir=freesurfer_subject_dir
+    ).format(name=name)
+    if not os.path.exists(base_patch_path) or recache:
+        if os.path.exists(base_patch_path):
+            os.remove(base_patch_path)
+        print("Writing base patch to %s..."%base_patch_path)
+        blender.write_volume_patch(fname, base_patch_path, hemi, blender_path=blender_path)
+    """
+    pts_v, polys_v, _  = freesurfer.get_surf(fs_subject, hemi, "patch", name, freesurfer_subject_dir=freesurfer_subject_dir)
+    pts_f, polys_f, _  = freesurfer.get_surf(fs_subject, hemi, "patch", name+".flat", freesurfer_subject_dir=freesurfer_subject_dir)
+    """
+
+    # Flatten
+    if flatten_with == 'blender':
+        assert method is not None, "method must be provided when using blender"
+
+        flat_patch_path = freesurfer.get_paths(
+            fs_subject, hemi, freesurfer_subject_dir=freesurfer_subject_dir
+        ).format(name=name + ".flat.blender")
+        if os.path.exists(flat_patch_path):
+            os.remove(flat_patch_path)
+
+        print("Generating flat patch via blender method %s to %s..."%(method, flat_patch_path))
+        done = blender.write_flat_patch(fname, flat_patch_path, hemi, method=method, blender_path=blender_path)
+        path_type, flat_type = "patch", "blender"
+    elif flatten_with == 'freesurfer':
+        print("Generating flat patch via freesurfer...")
+        done = freesurfer.flatten(
+            fs_subject, hemi, patch=name, freesurfer_subject_dir=freesurfer_subject_dir, **kwargs
+        )
+        path_type, flat_type = "patch", "freesurfer"
     elif flatten_with == 'SLIM':
-        done = flatten_slim(fs_subject, hemi, patch=name,
-                            freesurfer_subject_dir=freesurfer_subject_dir,
-                            **kwargs)
-        if not done:
-            # If flattening is aborted, skip the rest of this function
-            # (Do not attempt to import completed flatmaps)
-            return
-        if do_import_subject:
-            other = freesurfer.get_paths(fs_subject, "lh" if hemi == "rh" else "rh",
-                                         type='slim',
-                                         freesurfer_subject_dir=freesurfer_subject_dir)
-            other = other.format(name=name)
-            # If so, go ahead and import subject
-            if os.path.exists(other):
-                freesurfer.import_flat(fs_subject, name, cx_subject=cx_subject,
-                            flat_type='slim',
-                            freesurfer_subject_dir=freesurfer_subject_dir)
+        print("Generating flat patch via SLIM...")
+        done = flatten_slim(
+            fs_subject, hemi, patch=name, freesurfer_subject_dir=freesurfer_subject_dir, **kwargs
+        )
+        path_type, flat_type = "slip", "slim"
+    else:
+        raise ValueError(f"Invalid flatten_with: {flatten_with}")
 
-    return
+    # If flattening is aborted, skip the rest of this function
+    # (Do not attempt to import completed flatmaps)
+    if not done:
+        return
+    
+    # Import from freesurfer to pycortex DB
+    if do_import_subject:
+        hemi = "lh" if hemi == "rh" else "rh"
+        other = freesurfer.get_paths(
+            fs_subject, hemi, path_type, freesurfer_subject_dir=freesurfer_subject_dir
+        ).format(name=name)
+
+        if os.path.exists(other):  # Only when both hemi flats are present
+            freesurfer.import_flat(
+                fs_subject,
+                name,
+                cx_subject=cx_subject,
+                flat_type=flat_type,
+                auto_overwrite=auto_overwrite,
+                freesurfer_subject_dir=freesurfer_subject_dir,
+            )
 
 
 def flatten_slim(subject, hemi, patch, n_iterations=20, freesurfer_subject_dir=None,
-                 slim_path=slim_path, do_flatten=None):
+                 slim_path=default_slim_path, do_flatten=None):
     """Flatten brain w/ slim object flattening
 
     Parameters
diff --git a/docs/segmentation_guide.rst b/docs/segmentation_guide.rst
index 4f32d5af..eaf746b3 100644
--- a/docs/segmentation_guide.rst
+++ b/docs/segmentation_guide.rst
@@ -17,7 +17,7 @@ The brain model is exported to a 3D modeling program called Blender, where you w
 
 **3. Labeling ROIs**
 
-Here you will project functional data (semantic betas or localizer data, as well as retinotopic) onto the flatmaps, allowing you to label Regions of Interest on the brain - areas responsive to faces, scenes, or whatever else we’re analyzing.
+Here you will project functional data (semantic betas or localizer data, as well as retinotopic) onto the flatmaps, allowing you to label Regions of Interest on the brain - areas responsive to faces, scenes, or whatever else we're analyzing.
 
 In this guide, we will go over the first two steps.
 
@@ -46,7 +46,7 @@ To open freesurfer:
 For example:
     ``/auto/myfolder/freesurfer/SetUpFreeSurfer.sh``
 
-Create a “subjects” directory. If a "subjects" directory doesn't exist, make one in the FreeSurfer directory. Freesurfer is finicky about directories, so this step is crucial.
+Create a "subjects" directory. If a "subjects" directory doesn't exist, make one in the FreeSurfer directory. Freesurfer is finicky about directories, so this step is crucial.
 
 
 
@@ -69,7 +69,7 @@ In this case, you just give Freesurfer the name of the very first dicom file in
 For example:
     ``/auto/myfolder/anatomy/Subject/Subject_t1_nii -s Subject``
 
-The ‘-s Subject’ portion creates a folder, in this case a folder titled "Subject". The folder should be named for the subject.
+The '-s Subject' portion creates a folder, in this case a folder titled "Subject". The folder should be named for the subject.
 
 
 
@@ -102,14 +102,14 @@ At this stage, you just want to make sure that autorecon1 ran successfully and t
 of non-brain anatomy were not left behind. If big chunks of eye or skull were left behind, it is
 good to manually delete them yourself. If autorecon1 ran successfully, you can probably skip
 manual editing even if some anatomy was left behind since the next step, autorecon2, is quite
-accurate at determining brain surfaces even if non-brain anatomy was left behind. However, it’s good to double check that everything worked out.
+accurate at determining brain surfaces even if non-brain anatomy was left behind. However, it's good to double check that everything worked out.
 
 To pull up the newly stripped brains and make manual edits, type in your terminal:
     ``ipython``
 
     ``import cortex``
 
-    ``cortex.segment.fix_wm(‘Subject’)``
+    ``cortex.segment.fix_wm('Subject')``
 
 This should cause three windows to pop up: a mayavi viewer with the 3D brain, one of the brain in 2D, and one of a tool bar. At this point, you want to edit individual voxels. This mostly consists of getting rid of remaining skull and eyes. To do this, click the edit voxels tool on the toolbox bar or press A on your keyboard as a shortcut. After this, to delete voxels, simply right click the areas you wish to delete. If you erase something by accident and want to undo it, press CTRL + Z (this only works for the last thing you erased so be careful).
 
@@ -151,12 +151,12 @@ Tools > Configure volume brush
     Set Clone Source to Aux Volume
 
 This lets you paint from the aux volume to the mask. 
-Set Mode back to New Value if you’re done.
+Set Mode back to New Value if you're done.
 
 To change brush size:
 Tools > Configure brush info > Change Radius
 
-To change the size of the "paintbrush”, in the tool bar, go to: tools > configure brush info and
+To change the size of the "paintbrush", in the tool bar, go to: tools > configure brush info and
 change the radius. A shortcut to do the same thing is to press the numbers on the keypad of
 your keyboard (where 1 is 1x1, 4 is 4x4, etc).
 Generally you should just work with a 1-pixel radius, though.
@@ -206,9 +206,9 @@ not labeled as white matter when they should be. The command to make these edits
     
     ``import cortex``
     
-    ``cortex.segment.fix_wm(“subject”)``
+    ``cortex.segment.fix_wm("subject")``
 
-We’ll look through the results of autorecon2, examining the white matter curve and masks, and then the pial (gray matter) curve. This can be a lengthy process; because it’s an entirely nonverbal task, I recommend listening to podcasts as you go.    
+We'll look through the results of autorecon2, examining the white matter curve and masks, and then the pial (gray matter) curve. This can be a lengthy process; because it's an entirely nonverbal task, I recommend listening to podcasts as you go.    
 
 |
 
@@ -221,10 +221,94 @@ it shouldn't (such as gray matter and/or leftover pieces of eye or skull) as wel
 green/yellow surfaces. Make sure to hit "A" to switch to edit mode.
 
 
-Autorecon on the white matter surface should take about 2 hours. These manual edits are an iterative process; when it’s done, go back and look over the 3D surface, and make any changes that seem necessary. New spikes can appear in unexpected places, so three or four iterations may be needed, probably more if you are just starting to learn how to do it.
+Autorecon on the white matter surface should take about 2 hours. These manual edits are an iterative process; when it's done, go back and look over the 3D surface, and make any changes that seem necessary. New spikes can appear in unexpected places, so three or four iterations may be needed, probably more if you are just starting to learn how to do it.
 
 
+Making cuts
+##################################
+
+After completing the segmentation phase, the next step is to make cuts in the brain surface to prepare it for flattening. This process involves creating cuts along the brain's sulci to transform the 3D surface into a 2D flatmap with minimal distortion.
+
+PyCortex provides three different methods for cutting and flattening brain surfaces:
+
+**1. Freesurfer (Recommended)**
+The traditional and most reliable method that uses Freesurfer's `mris_flatten` command. This method produces high-quality flatmaps with minimal distortion but takes approximately 2 hours per hemisphere.
+
+**2. SLIM**
+An experimental method using the SLIM algorithm that is very fast but tends to leave more distortions in the flatmap. Requires additional installation of the SLIM dependency.
+
+**3. Blender**
+A newer method that uses Blender's UV unwrapping capabilities for faster flattening (typically 5-15 minutes per hemisphere). While faster, it may introduce more distortion compared to Freesurfer.
+
+The complete process begins with manual cutting in Blender, where you'll make cuts to prepare the surface for flattening. Once the cuts are complete, the cut surface is automatically flattened using your chosen method. Finally, the resulting flatmap is imported into PyCortex for visualization and analysis.
+
+You may follow the steps below or a `Python notebook <https://colab.research.google.com/github/dmitry-mli/pycortex/blob/blender-flattening-support/examples/quickstart/fmri_flattening.ipynb>`_.
+
+Step 1: Manual Cutting in Blender
+***************************************************
+
+Start the cutting process by calling `cortex.segment.cut_surface()`. This function will create a Blender file with your brain surface, open Blender automatically, and allow you to make manual cuts for the left hemisphere.
+
+.. code-block:: python
+
+    import cortex
+    
+    cortex.segment.cut_surface(
+        "sub-01",                   # Your subject ID
+        "lh",                       # Left hemisphere
+        name="flatten",             # Name for this flattening attempt
+        flatten_with="freesurfer",  # Or "SLIM" or "blender"
+        recache=True,               # Force recache of the subject
+        do_import_subject=False,    # Don't import until both hemispheres are done
+    )
+
+To make the cuts please watch the `cutting tutorial video <https://www.youtube.com/watch?v=D4tylQ_mMuM>`_.
+
+Step 2: Repeat for Right Hemisphere
+***************************************************
+
+After completing the left hemisphere, repeat the process for the right hemisphere.
+
+.. code-block:: python
 
+    cortex.segment.cut_surface(
+        "sub-01",                   # Your subject ID
+        "rh",                       # Right hemisphere
+        name="flatten",             # Name for this flattening attempt
+        flatten_with="freesurfer",  # Or "SLIM" or "blender"
+        recache=True,               # Force recache of the subject
+        auto_overwrite=True,        # Overwrite PyCortex record
+        do_import_subject=True,     # Import both hemispheres when done
+    )
 
+After completing both hemispheres, your flatmap will be automatically imported into PyCortex and ready for visualization and analysis.
+
+
+Step 3: Verify the cuts
+***************************************************
 
+After completing both hemispheres, your flatmap will be automatically imported into PyCortex and ready for visualization and analysis.
+
+To verify that your cuts and flattening worked correctly, you can visualize the results using PyCortex's visualization tools. Here's a verification script:
+
+.. code-block:: python
+
+    import cortex
+    import numpy as np
+    from matplotlib import pyplot as plt
+    
+    test_data = np.random.rand(1000)  # Random data
+    
+    vol = cortex.Volume(
+        test_data,
+        subject="sub-01",
+        xfmname="full",
+        vmin=0,
+        vmax=1
+    )
+    
+    # Display the visualization on the flatmap
+    cortex.quickshow(vol, with_colorbar=True, recache=True)
+    plt.show()
 
+Alternatively, you may follow one of the examples from the gallery.
diff --git a/examples/quickstart/fmri_flattening.ipynb b/examples/quickstart/fmri_flattening.ipynb
new file mode 100644
index 00000000..84be61a9
--- /dev/null
+++ b/examples/quickstart/fmri_flattening.ipynb
@@ -0,0 +1,220 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Functional MRI pial flattening in 20 minutes\n",
+    "\n",
+    "This notebook demonstrates how to use PyCortex to:\n",
+    "1. Cut and flatten brain surfaces in under 20 minutes;\n",
+    "2. Preview flattened surface in under 10 seconds;\n",
+    "3. Visualize BOLD fMRI data on the flattened surfaces"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Prerequisites\n",
+    "\n",
+    "You will need to run auto"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### 1. Reconcile subject's anatomy with Freesurfer\n",
+    "\n",
+    "Make sure you have fully reconciled the subject's. Detailed guidance is out of scope for this notebook, but roughly you should've run these lines:\n",
+    "\n",
+    "```\n",
+    "recon-all -autorecon1 -s sub-01\n",
+    "recon-all -autorecon2 -s sub-01\n",
+    "recon-all -autorecon3 -s sub-01\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### 2. Install Python dependencies"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import cortex\n",
+    "import nibabel as nib\n",
+    "import numpy as np\n",
+    "from matplotlib import pyplot as plt"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Step 1: Align Function MRI to Anatomical\n",
+    "\n",
+    "`cortex.align.automatic` perform automatic alignment using Freesurfer's boundary-based registration. It is fully automated. This should take ~1 minute."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "bold_nii_gz_path = \"path to bold.nii.gz\"\n",
+    "cortex.align.automatic(\n",
+    "    subject=\"sub-01\",\n",
+    "    xfmname=\"full\",\n",
+    "    reference=bold_nii_gz_path,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Step 2: Surface Cutting and Flattening\n",
+    "\n",
+    "`cortex.segment.cut_surface` performs surface cutting. It is semi-automatic, meaning it requires a bit of manual work. It will start Blender for us to manually specify the cuts. Make sure Blender is installed on your system.\n",
+    "\n",
+    "#### A. Left hemisphere\n",
+    "\n",
+    "Let's start with the left hemisphere.\n",
+    "\n",
+    "Please follow [this guide](https://www.youtube.com/watch?v=D4tylQ_mMuM) on cuts. This will take ~15 minutes."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cortex.segment.cut_surface(\n",
+    "    \"sub-01\",\n",
+    "    \"lh\",\n",
+    "    name=\"exp\",\n",
+    "    flatten_with=\"blender\",\n",
+    "    method=\"CONFORMAL\",\n",
+    "    recache=True,\n",
+    "    do_import_subject=False,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### B. Right hemisphere\n",
+    "\n",
+    "Then repeat the same for the right hemisphere.\n",
+    "\n",
+    "This will take another ~5 minutes.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Right hemisphere\n",
+    "cortex.segment.cut_surface(\n",
+    "    \"sub-01\",\n",
+    "    \"rh\",\n",
+    "    name=\"exp\",\n",
+    "    flatten_with=\"blender\",\n",
+    "    method=\"CONFORMAL\",\n",
+    "    recache=True,\n",
+    "    auto_overwrite=True,\n",
+    "    do_import_subject=True,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Step 3: Visualize fMRI\n",
+    "\n",
+    "Now we'll load the BOLD fMRI data and prepare it for visualization. We'll compute the temporal mean of the BOLD signal to create a static visualization.\n",
+    "\n",
+    "Then we can use `cortex.quickshow` to visualize the BOLD data on the flattened surfaces using pycortex.\n",
+    "\n",
+    "This will take 30 seconds."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Load functional data\n",
+    "bold_path = \"path to bold.nii.gz\"\n",
+    "img = nib.load(bold_path)\n",
+    "data = img.get_fdata()\n",
+    "\n",
+    "# Compute temporal mean\n",
+    "mean_vol = np.mean(data, axis=3)\n",
+    "\n",
+    "# Transpose to match expected orientation (96,96,76) → (76,96,96)\n",
+    "mean_vol = np.transpose(mean_vol, (2, 1, 0))\n",
+    "\n",
+    "# Select cuboid\n",
+    "cuboid = np.array([[0, 20, 0], [50, 70, 70]])  # (2, 3) as a pair of (min, max) voxels\n",
+    "print(mean_vol.shape)\n",
+    "cuboid_slice = tuple(slice(*span) for span in cuboid.transpose())\n",
+    "mean_vol_filtered = np.zeros_like(mean_vol)\n",
+    "mean_vol_filtered[cuboid_slice] = mean_vol[cuboid_slice]\n",
+    "mean_vol = mean_vol_filtered\n",
+    "\n",
+    "\n",
+    "# Create Volume object for visualization\n",
+    "vol = cortex.Volume(\n",
+    "    mean_vol,\n",
+    "    subject=\"sub-01\",\n",
+    "    depth=1,\n",
+    "    height=512,\n",
+    "    xfmname=\"full\",\n",
+    "    vmin=np.min(mean_vol),\n",
+    "    vmax=np.max(mean_vol),\n",
+    ")\n",
+    "\n",
+    "# Display the visualization\n",
+    "cortex.quickshow(vol, with_colorbar=True, recache=True)\n",
+    "plt.show()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.17"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}