From 558d3c039d587790f045713b6b848f5670f701c9 Mon Sep 17 00:00:00 2001 From: Yujing Huang Date: Mon, 29 Jun 2026 16:09:05 -0400 Subject: [PATCH 1/3] implement surfa.mesh.FSLabelFile class - transform coordinates between 'world', 'voxel', and 'surface' spaces - if inverse deformation field is provide, apply the defomration field - the deformation field is expected to be in abs-crs data representation --- surfa/mesh/__init__.py | 1 + surfa/mesh/fslabelfile.py | 151 ++++++++++++++++++++++++++++++++++++++ surfa/transform/space.py | 30 ++++++++ surfa/transform/warp.py | 22 ++++++ 4 files changed, 204 insertions(+) create mode 100644 surfa/mesh/fslabelfile.py diff --git a/surfa/mesh/__init__.py b/surfa/mesh/__init__.py index 33d033c..ddacb94 100644 --- a/surfa/mesh/__init__.py +++ b/surfa/mesh/__init__.py @@ -5,3 +5,4 @@ from .overlay import cast_overlay from .distance import surface_distance from .timeseries import TimeSeries +from .fslabelfile import FSLabelFile diff --git a/surfa/mesh/fslabelfile.py b/surfa/mesh/fslabelfile.py new file mode 100644 index 0000000..38a2e6d --- /dev/null +++ b/surfa/mesh/fslabelfile.py @@ -0,0 +1,151 @@ +import sys +import numpy as np + +from surfa.transform import Space + + +# this class implements the following methods: +# 1. read(): read list of coords in Freesurfer .label format +# 2. write(): write list of coords in Freesurfer .label format +# 3. transform(): transform coordinates between 'world', 'voxel', and 'surface' spaces +# if inverse deformation field is provide, apply the defomration field +# the deformation field is expected to be in abs-crs data representation +class FSLabelFile(): + + def __init__(self, geom=None, coords=None, surface_vertices=None, coordspace=Space.FS_COORDS_UNKNOWN): + # [N, 3] array containing the coordinates + self._coords = coords + # [N] 1-D array containing the surface vertex number + self._surface_vertices = surface_vertices + # coordinates space + self._coordspace = coordspace + # number of coords + self._num_points = 0 + # volume geometry + self._geom = geom + + if (self._surface_vertices is None and self._coords is not None): + self._surface_vertices = np.zeros(self._coords.shape[0], dtype=int) + + + # read list of coords in Freesurfer .label format + # For Freesurfer .label format, + # see https://surfer.nmr.mgh.harvard.edu/fswiki/LabelsClutsAnnotationFiles#Label_file + def read(self, infile): + fp = open(infile, "r") + + # read the comment line, locate coordinate space + line = fp.readline().strip() + index = line.find("vox2ras=") # find the index where 'vox2ras=' starts + if (index != -1): + loc = index + len("vox2ras=") + space = line[loc:] # skip pass 'vox2ras=' + if (space == 'voxel'): + self._coordspace = Space.FS_COORDS_VOXEL + elif (space == 'scanner'): + self._coordspace = Space.FS_COORDS_SCANNER_RAS + elif (space == 'TkReg'): + self._coordspace = Space.FS_COORDS_TKREG_RAS + else: + raise ValueError(f"Invalid `vox2ras`: {space}. It can either 'voxel', 'scanner', or 'TkReg'") + + # read in the count of coords + line = fp.readline().strip() + self._num_points = int(line) + + # read all the coords + coords, surface_vertices = [], [] + line = fp.readline().strip() + while (line): + # .split() with no arguments automatically handles any consecutive whitespace (multiple spaces, tabs, or newlines) + # as a single delimiter. It also automatically ignores empty elements caused by leading or trailing spaces. + v, x, y, z, stat = line.split() + coords.append((float(x), float(y), float(z))) + surface_vertices.append(int(v)) + + # next line + line = fp.readline().strip() + + self._coords = np.array(coords) # [N, 3] + self._surface_vertices = np.array(surface_vertices) # [N] + + fp.close() + + + # write list of coords in Freesurfer .label format + def write(self, outfile=None): + fp = sys.stdout + if (outfile is not None): + fp = open(outfile, "w") + + # write the comment line + if (self._coordspace == Space.FS_COORDS_TKREG_RAS): + space = "TkReg" + elif (self._coordspace == Space.FS_COORDS_SCANNER_RAS): + space = "scanner" + elif (self._coordspace == Space.FS_COORDS_VOXEL): + space = "voxel" + fp.write(f"#!ascii label, vox2ras={space}\n") + + # write number of coodinates + num_points = self._coords.shape[0] + fp.write(f"{num_points}\n") + + # write the coordinates + for v in range(num_points): + fp.write(f"{self._surface_vertices[v]} {self._coords[v, 0]:.3f} {self._coords[v, 1]:.3f} {self._coords[v, 2]:.3f} 0.0000000000\n") + + # can't close stdout + if (outfile is not None): + fp.close() + + + # transform the coordinates into a different space specified as 'coordspace' + # apply the warp field first if 'warp' is specified + def transform(self, warp=None, coordspace=Space.FS_COORDS_TKREG_RAS, geom=None): + if (warp is not None): + self.apply_warp(warp) + self.change_space(coordspace=coordspace, geom=geom) + + + # apply warp field to the coordinates + # 'warp' is the inverse (backward) warp field mapping from source (moving) to target (fixed) + def apply_warp(self, warp): + # the coordinates needs to be FS_COORDS_VOXEL space assuming the warp in abs-crs + if (self._coordspace != Space.FS_COORDS_VOXEL): + self.change_space(coordspace=Space.FS_COORDS_VOXEL) + + from surfa.transform import Warp + self._coords = Warp.warp_pointsets(warp, self._coords) + + + # change coordinates space + def change_space(self, coordspace=Space.FS_COORDS_TKREG_RAS, geom=None): + if (self._coordspace == coordspace): + return + + transform_geom = geom if (geom is not None) else self._geom + space_from = Space.CODE_TO_LETTER[self._coordspace] + space_to = Space.CODE_TO_LETTER[coordspace] + affine = transform_geom.affine(space_from, space_to) + self._coords = affine[:3, :3] @ self._coords.T + affine[:3, 3:] + self._coords = self._coords.T # [ N, 3] + self._coordspace = coordspace + + + @property + def coordspace(self): + return self._coordspace + + @property + def geom(self): + return self._geom + + @property + def coords(self): + return self._coords + + @property + def surface_vertices(self): + return self._surface_vertices + diff --git a/surfa/transform/space.py b/surfa/transform/space.py index 5e73bc5..d023727 100644 --- a/surfa/transform/space.py +++ b/surfa/transform/space.py @@ -1,5 +1,35 @@ class Space: + """ + Here are the 3 spaces supported (see https://github.com/freesurfer/freesurfer/blob/dev/include/mri.h) + FS_COORDS_UNKNOWN + FS_COORDS_TKREG_RAS - Freesurfer surface RAS space + FS_COORDS_SCANNER_RAS - scanner RAS space + FS_COORDS_VOXEL - voxel space + """ + # see defines in https://github.com/freesurfer/freesurfer/blob/dev/include/mri.h + FS_COORDS_UNKNOWN = 0 + FS_COORDS_TKREG_RAS = 1 + FS_COORDS_SCANNER_RAS = 2 + FS_COORDS_VOXEL = 3 + + # to be consistent with sf.transform.Space + # 'voxel' = voxel space + # 'world' = scanner RAS space + # 'surface' = Freesurfer surface RAS space (tkreg space) + LETTER_TO_CODE = { + 'voxel' : FS_COORDS_VOXEL, + 'world' : FS_COORDS_SCANNER_RAS, + 'surface' : FS_COORDS_TKREG_RAS, + } + + # consistent with sf.transform.Space + CODE_TO_LETTER = { + FS_COORDS_VOXEL : 'voxel', + FS_COORDS_SCANNER_RAS : 'world', + FS_COORDS_TKREG_RAS : 'surface', + } + def __init__(self, name): """ Coordinate space representation. Supported spaces are: diff --git a/surfa/transform/warp.py b/surfa/transform/warp.py index 84a4017..f9e12a3 100644 --- a/surfa/transform/warp.py +++ b/surfa/transform/warp.py @@ -268,3 +268,25 @@ def target(self, value): Set target (or fixed) image geometry. Invokes parent setter. """ self.geom = value + + + @staticmethod + def warp_pointsets(warp, pointsets): + from scipy.interpolate import RegularGridInterpolator + + data = warp.data + nx, ny, nz, _ = data.shape + grid = (np.arange(nx), np.arange(ny), np.arange(nz)) + + # interp(source_crs) = target crs + interp = [ RegularGridInterpolator(grid, data[..., k], bounds_error=False, fill_value=np.nan) + for k in range(data.shape[-1]) ] + + def _sample_warpfield(in_crs): + return np.array([f(in_crs) for f in interp]) + + # input pointsets is [N, 3] array + pointsets = np.array([f(pointsets) for f in interp]) + pointsets = pointsets.T # [N, 3] + + return pointsets \ No newline at end of file From 734f5a49c18c44d846d7eb455a37cca545b106aa Mon Sep 17 00:00:00 2001 From: Yujing Huang Date: Mon, 29 Jun 2026 16:11:55 -0400 Subject: [PATCH 2/3] fix missing newline at end of file --- surfa/transform/warp.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/surfa/transform/warp.py b/surfa/transform/warp.py index f9e12a3..5d077d0 100644 --- a/surfa/transform/warp.py +++ b/surfa/transform/warp.py @@ -289,4 +289,4 @@ def _sample_warpfield(in_crs): pointsets = np.array([f(pointsets) for f in interp]) pointsets = pointsets.T # [N, 3] - return pointsets \ No newline at end of file + return pointsets From ce93cea3e5ca1057ee904df75e9d0d17db098772 Mon Sep 17 00:00:00 2001 From: Yujing Huang Date: Tue, 30 Jun 2026 13:43:08 -0400 Subject: [PATCH 3/3] convert warp_pointsets() input Warp to abs-crs if it is necessary --- surfa/mesh/fslabelfile.py | 2 +- surfa/transform/warp.py | 4 ++++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/surfa/mesh/fslabelfile.py b/surfa/mesh/fslabelfile.py index 38a2e6d..ea12d15 100644 --- a/surfa/mesh/fslabelfile.py +++ b/surfa/mesh/fslabelfile.py @@ -109,7 +109,7 @@ def transform(self, warp=None, coordspace=Space.FS_COORDS_TKREG_RAS, geom=None): # apply warp field to the coordinates - # 'warp' is the inverse (backward) warp field mapping from source (moving) to target (fixed) + # 'warp' is the inverse (backward) warp field mapping from target (fixed) to source (moving) def apply_warp(self, warp): # the coordinates needs to be FS_COORDS_VOXEL space assuming the warp in abs-crs if (self._coordspace != Space.FS_COORDS_VOXEL): diff --git a/surfa/transform/warp.py b/surfa/transform/warp.py index 5d077d0..664c273 100644 --- a/surfa/transform/warp.py +++ b/surfa/transform/warp.py @@ -272,6 +272,10 @@ def target(self, value): @staticmethod def warp_pointsets(warp, pointsets): + # the input warp needs to be abs_crs data interpretation + if (warp.format != Warp.Format.abs_crs): + warp = warp.convert(format=Warp.Format.abs_crs) + from scipy.interpolate import RegularGridInterpolator data = warp.data