vox_bool = mgz_data == vol_id
// Get the 3 dimensional indices in voxel space
vox_ijk = np.array(np.where(vox_bool)).T
// Transform to MRI coordinates (where our surfaces live)
_, vox_mri_t, _, _, _ = _read_mri_info(mri)
rr_voi = apply_trans(vox_mri_t, vox_ijk) // mri voxels -> MRI surface RAS
// Filter out points too far from volume region voxels
dists = _compute_nearest(rr_voi, rr, return_dists=True)[1]
// Maximum distance from center of mass of a voxel to any of its corners
maxdist = linalg.norm(vox_mri_t["trans"][:3, :3].sum(0) / 2.)
return dists <= maxdist