From 2206df8da7fde205580faec58fffdd9513afe9e8 Mon Sep 17 00:00:00 2001 From: Martin Reuter Date: Mon, 18 Dec 2023 18:51:22 +0100 Subject: [PATCH 1/6] initial checkin sample_parc --- recon_surf/sample_parc.py | 336 +++++++++++++++++++++++++++++++++++++ recon_surf/smooth_aparc.py | 124 ++++++++------ 2 files changed, 412 insertions(+), 48 deletions(-) create mode 100644 recon_surf/sample_parc.py diff --git a/recon_surf/sample_parc.py b/recon_surf/sample_parc.py new file mode 100644 index 00000000..5bb77918 --- /dev/null +++ b/recon_surf/sample_parc.py @@ -0,0 +1,336 @@ +#!/usr/bin/env python3 + + +# Copyright 2024 Image Analysis Lab, German Center for Neurodegenerative Diseases (DZNE), Bonn +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + + +# IMPORTS +import optparse +import sys +import numpy as np +from numpy import typing as npt +import nibabel.freesurfer.io as fs +from nibabel import load as nibload +from scipy import sparse +from lapy import TriaMesh +from smooth_aparc import smooth_aparc + + +HELPTEXT = """ +Script to sample labels from image to surface and clean up. + +USAGE: +sample_parc --inseg --insurf --incort + --seglut --surflut --outaparc + --projmm --radius + + +Dependencies: + Python 3.8 + + Numpy + http://www.numpy.org + + Nibabel to read and write FreeSurfer surface meshes + http://nipy.org/nibabel/ + + +Original Author: Martin Reuter +Date: Dec-18-2023 + +""" + +h_inseg = "path to input segmentation image" +h_incort = "path to input cortex label mask" +h_insurf = "path to input surface" +h_outaparc = "path to output aparc" +h_surflut = "FreeSurfer look-up-table for values on surface" +h_seglut = "Look-up-table for values in segmentation image (rows need to correspond to surflut)" +h_projmm = "Sample along normal at projmm distance (in mm), default 0" +h_radius = "Search around sample location at radius (in mm) for label if 'unknown', default None" + + +def options_parse(): + """Command line option parser. + + Returns + ------- + options + object holding options + + """ + parser = optparse.OptionParser( + version="$Id: smooth_aparc,v 1.0 2018/06/24 11:34:08 mreuter Exp $", + usage=HELPTEXT, + ) + parser.add_option("--inseg", dest="inseg", help=h_inseg) + parser.add_option("--insurf", dest="insurf", help=h_insurf) + parser.add_option("--incort", dest="incort", help=h_incort) + parser.add_option("--surflut", dest="surflut", help=h_surflut) + parser.add_option("--seglut", dest="seglut", help=h_seglut) + parser.add_option("--outaparc", dest="outaparc", help=h_outaparc) + parser.add_option("--projmm", dest="projmm", help=h_projmm, default=0.0, type="float") + parser.add_option("--radius", dest="radius", help=h_radius, default=None, type="float") + (options, args) = parser.parse_args() + + if options.insurf is None or options.inseg is None or options.outaparc is None: + sys.exit("ERROR: Please specify input surface, input image and output aparc!") + + if options.surflut is None or options.seglut is None: + sys.exit("ERROR: Please specify surface and segmentatin image LUT!") + + # maybe later add functionality, to not have a cortex label, e.g. + # like FreeSurfer find largest connected component and fill only + # the other unknown regions + if options.incort is None: + sys.exit("ERROR: Please specify surface cortex label!") + + return options + + + +def sample_nearest_nonzero(img, vox_coords, radius=3.0): + """Sample closest non-zero value in a ball of radius around vox_coords. + + Parameters + ---------- + img : nibabel.image + Image to sample. Voxels need to be isotropic. + vox_coords : ndarray float shape(n,3) + Coordinates in voxel space around which to search. + radius : float default 3.0 + Consider all voxels inside this radius to find a non-zero value. + + Returns + ------- + samples : np.ndarray (n,) + Sampled values. Retruns zero for vertices where values are zero in ball. + """ + # check for isotropic voxels + voxsize = img.header.get_zooms() + print("Check isotropic vox sizes: {}".format(voxsize)) + assert (np.max(np.abs(voxsize - voxsize[0])) < 0.001), 'Voxels not isotropic!' + data = np.asarray(img.dataobj) + + # radius in voxels: + rvox = radius * voxsize[0] + + # sample window around nearest voxel + x_nn = np.rint(vox_coords).astype(int) + # Reason: to always have the same number of voxels that we check + # and to be consistent with FreeSurfer, we center the window at + # the nearest neighbor voxel, instead of at the float vox coordinates + + # create box with 2*rvox+1 side length to fully contain ball + # and get coordiante offsets with zero at center + ri = np.floor(rvox).astype(int) + l = np.arange(-ri,ri+1) + xv, yv, zv = np.meshgrid(l, l, l) + dd = np.sqrt(xv*xv + yv*yv + zv*zv).flatten() + + # flatten and keep only sphere with radius + xv = xv.flatten()[dd<=rvox] + yv = yv.flatten()[dd<=rvox] + zv = zv.flatten()[dd<=rvox] + dd = dd[dd<=rvox] + + # stack to get offset vectors + offsets = np.column_stack((xv, yv, zv)) + + # sort offsets according to distance + # Note: we keep the first zero voxel so we can later + # determine if all voxels are zero with the argmax trick + sortidx = np.argsort(dd) + dd = dd[sortidx] + offsets = offsets[sortidx,:] + + # reshape and tile to add to list of coords + n = x_nn.shape[0] + toffsets = np.tile(offsets.transpose().reshape(1,3,offsets.shape[0]),(n,1,1)) + s_coords = x_nn[:, :, np.newaxis] + toffsets + + # get image data at the s_coords locations + s_data = data[s_coords[:,0], s_coords[:,1], s_coords[:,2]] + + # get first non-zero if possible + nzidx = (s_data!=0).argmax(axis=1) + # the above return index zero if all elements are zero which is OK for us + # as we can then sample there and get a value of zero + samples = s_data[np.arange(s_data.shape[0]),nzidx] + return samples + + +def sample_img(surf, img, cortex=None, projmm=0.0, radius=None): + """Sample volume at a distance from the surface. + + Parameters + ---------- + surf : tuple | str + Surface as returned by nibabel fs.read_geometry, where: + surf[0] is the np.array of (n, 3) vertex coordinates and + surf[1] is the np.array of (m, 3) triangle indices. + If type is str, read surface from file. + img : nibabel.image | str + Image to sample. + If type is str, read image from file. + cortex : np.ndarray | str + Filename of cortex label or np.array with cortex indices + projmm : float + Sample projmm mm along the surface vertex normals (default=0). + radius : float [optional] | None + If given and if the sample is equal to zero, then consider + all voxels inside this radius to find a non-zero value. + + Returns + ------- + samples : np.ndarray (n,) + Sampled values. + """ + if isinstance(surf, str): + surf = fs.read_geometry(surf, read_metadata=True) + if isinstance(img, str): + img = nibload(img) + if isinstance(cortex, str): + cortex = fs.read_label(cortex) + nvert = surf[0].shape[0] + # Compute Cortex Mask + if cortex is not None: + mask = np.zeros(nvert, dtype=bool) + mask[cortex] = True + else: + mask = np.ones(nvert, dtype=bool) + + data = np.asarray(img.dataobj) + # Use LaPy TriaMesh for vertex normal computation + T = TriaMesh(surf[0], surf[1]) + # compute sample coordinates projmm mm along the surface normal + # in surface RAS coordiante system: + x = T.v + projmm * T.vertex_normals() + # mask cortex + xx = x[mask] + + # compute Transformation from surface RAS to voxel space: + Torig = img.header.get_vox2ras_tkr() + Tinv = np.linalg.inv(Torig) + x_vox = np.dot(xx, Tinv[:3, :3].T) + Tinv[:3, 3] + # sample at nearest voxel + x_nn = np.rint(x_vox).astype(int) + samples_nn = data[x_nn[:,0], x_nn[:,1], x_nn[:,2]] + # no search for zeros, done: + if not radius: + samplesfull = np.zeros(nvert, dtype="int") + samplesfull[mask] = samples_nn + return samplesfull + # search for zeros, but no zeros exist, done: + zeros = np.asarray(samples_nn==0).nonzero()[0] + if zeros.size == 0: + samplesfull = np.zeros(nvert, dtype="int") + samplesfull[mask] = samples_nn + return samplesfull + # here we need to do the hard work of searching in a windows + # for non-zero samples + print("sample_img: found {} zero samples, searching radius ...",zeros.size) + z_nn = x_nn[zeros] + z_samples = sample_nearest_nonzero(img, z_nn, radius=radius) + samples_nn[zeros] = z_samples + samplesfull = np.zeros(nvert, dtype="int") + samplesfull[mask] = samples_nn + return samplesfull + + +def replace_labels(img_labels, img_lut, surf_lut): + """Replace image labels with corresponding surface labels or unknown. + + Parameters + ---------- + img_labels : np.ndarray(n,) + Array with imgage label ids. + img_lut : str + Filename for image label look up table. + surf_lut : str + Filename for surface label look up table. + + Returns + ------- + surf_labels : np.ndarray (n,) + Array with surface label ids. + surf_ctab : np.ndarray shape(m,4) + Surface color table (RGBA). + surf_names : np.ndarray[str] shape(m,) + Names of label regions. + """ + surflut = np.loadtxt(surf_lut, usecols=(0,2,3,4,5), dtype="int") + surfids = surflut[:,0] + surf_ctab = surflut[:,1:5] + surf_names = np.loadtxt(surf_lut, usecols=(1), dtype="str") + imglut = np.loadtxt(img_lut, usecols=(0,2,3,4,5), dtype="int") + img_names = np.loadtxt(img_lut, usecols=(1), dtype="str") + assert (np.all(img_names == surf_names)), "Label names in the LUTs do not agree!" + # find labels that have no replacement and map to 0=unknonw + missing = np.logical_not(np.isin(img_labels, surfids)) + img_labels = img_labels.copy() + img_labels[missing] = 0 + # translate the remaining labels + newids = imglut[:,0] + d = dict(zip(surfids, newids)) + surf_labels=np.asarray([d.get(e, e) for e in img_labels]) + return surf_labels, surf_ctab, surf_names + + +def sample_parc (surf, seg, imglut, surflut, outaparc, cortex=None, projmm=0.0, radius=None): + """Replace image labels with corresponding surface labels or unknown. + + Parameters + ---------- + surf : tuple | str + Surface as returned by nibabel fs.read_geometry, where: + surf[0] is the np.array of (n, 3) vertex coordinates and + surf[1] is the np.array of (m, 3) triangle indices. + If type is str, read surface from file. + img : nibabel.image | str + Image to sample. + If type is str, read image from file. + imglut : str + Filename for image label look up table. + surflut : str + Filename for surface label look up table. + outaparc : str + Filename for output surface parcellation. + cortex : np.ndarray | str + Filename of cortex label or np.ndarray with cortex indices + projmm : float + Sample projmm mm along the surface vertex normals (default=0). + radius : float [optional] | None + If given and if the sample is equal to zero, then consider + all voxels inside this radius to find a non-zero value. + """ + if isinstance(cortex, str): + cortex = fs.read_label(cortex) + if isinstance(surf, str): + surf = fs.read_geometry(surf, read_metadata=True) + samples = sample_img(surf, seg, cortex, projmm, radius) + surfsamples, surfctab, surfnames = replace_labels(samples, imglut, surflut) + smooths = smooth_aparc(surf, surfsamples, cortex) + fs.write_annot(outaparc, smooths, ctab=surfctab, names=surfnames) + + +if __name__ == "__main__": + # Command Line options are error checking done here + options = options_parse() + + sample_parc(options.insurf, options.inseg, options.seglut, options.surflut, options.outaparc, options.incort, options.projmm, options.radius) + + sys.exit(0) + diff --git a/recon_surf/smooth_aparc.py b/recon_surf/smooth_aparc.py index fbd50e78..89f91163 100644 --- a/recon_surf/smooth_aparc.py +++ b/recon_surf/smooth_aparc.py @@ -79,19 +79,19 @@ def options_parse(): def get_adjM(trias: npt.NDArray, n: int): - """[MISSING]. + """Create symmetric sparse adjacency matrix of triangle mesh. Parameters ---------- - trias : npt.NDArray + trias : np.ndarray shape (m, 3) int n : int - Shape of tje matrix + Shape of output (n,n) adjaceny matrix, where n>=m. Returns ------- - adjM : np.ndarray - Adjoint matrix + adjM : np.ndarray (bool) shape (n,n) + Symmetric adjacency matrix, true corresponds to an edge. """ I = trias @@ -112,13 +112,13 @@ def bincount2D_vectorized(a: npt.NDArray) -> np.ndarray: Parameters ---------- - a : npt.NDArray - Array + a : np.ndarray + Input 2D array of non-negative ints. Returns ------- np.ndarray - Array of counted values + Array of counted values. """ N = a.max() + 1 @@ -129,26 +129,29 @@ def bincount2D_vectorized(a: npt.NDArray) -> np.ndarray: def mode_filter( adjM: sparse.csr_matrix, labels: npt.NDArray[str], - fillonlylabel: str = "", + fillonlylabel = None, novote: npt.ArrayLike = [] -) -> npt.NDArray[str]: - """[MISSING]. +) -> npt.NDArray[int]: + """Apply mode filter (smoothing) to int labels on mesh vertices. Parameters ---------- - adjM : sparse.csr_matrix - Adjoint matrix - labels : npt.NDArray[str] - List of labels - fillonlylabel : str - Label to fill exclusively. Defaults to "" + adjM : sparse.csr_matrix[bool] + Symmetric adjacency matrix defining edges between vertices. + This determines what edges can vote so usually one adds the + identity to the adjacency matrix so that each vertex is included + in its own vote. + labels : npt.NDArray[int] + List of integer labels. + fillonlylabel : int + Label to fill exclusively. Defaults to None. novote : npt.ArrayLike - Entries that should not vote. Defaults to [] + Entries that should not vote. Defaults to []. Returns ------- - labels_new - New filtered labels + labels_new : npt.NDArray[int] + New smoothed labels. """ # make sure labels lengths equals adjM dimension @@ -231,7 +234,7 @@ def mode_filter( rempty += 1 continue # print(str(rvals)) - mvals = mode(rvals)[0] + mvals = mode(rvals, keepdims=True)[0] # print(str(mvals)) if mvals.size != 0: # print(str(row)+' '+str(ids[row])+' '+str(mvals[0])) @@ -246,35 +249,23 @@ def mode_filter( return labels_new -def smooth_aparc( - insurfname: str, - inaparcname: str, - incortexname: str, - outaparcname: str -) -> None: +def smooth_aparc(surf, labels, cortex = None): """Smoothes aparc. Parameters ---------- - insurfname : str - Suface filepath and name of source - inaparcname : str - Annotation filepath and name of source - incortexname : str - Label filepath and name of source - outaparcname : str - Suface filepath and name of destination + surf : nibabel surface + Suface filepath and name of source. + labels : np.array[int] + Labels at each vertex (int). + cortex : np.array[int] + Vertex ids inside cortex mask. + Returns + ------- + smoothed_labels : np.array[int] + Smoothed labels. """ - # read input files - print("Reading in surface: {} ...".format(insurfname)) - surf = fs.read_geometry(insurfname, read_metadata=True) - print("Reading in annotation: {} ...".format(inaparcname)) - aparc = fs.read_annot(inaparcname) - print("Reading in cortex label: {} ...".format(incortexname)) - cortex = fs.read_label(incortexname) - # set labels (n) and triangles (n x 3) - labels = aparc[0] faces = surf[1] nvert = labels.size if labels.size != surf[0].shape[0]: @@ -286,8 +277,11 @@ def smooth_aparc( ) # Compute Cortex Mask - mask = np.zeros(labels.shape, dtype=bool) - mask[cortex] = True + if cortex is not None: + mask = np.zeros(labels.shape, dtype=bool) + mask[cortex] = True + else: + mask = np.ones(labels.shape, dtype=bool) # check if we have places where non-cortex has some labels noncortnum = np.where(~mask & (labels != -1)) print( @@ -320,6 +314,7 @@ def smooth_aparc( # print("minlab: "+str(np.min(labels))+" maxlab: "+str(np.max(labels))) # set all labels inside cortex that are -1 or 0 to fill label + labels = labels.copy() fillonlylabel = np.max(labels) + 1 labels[mask & (labels == -1)] = fillonlylabel labels[mask & (labels == 0)] = fillonlylabel @@ -356,14 +351,47 @@ def smooth_aparc( labels = mode_filter(adjM, labels) # set labels outside cortex to -1 labels[~mask] = -1 + return labels + + +def smooth_aparc_files( + insurfname: str, + inaparcname: str, + incortexname: str, + outaparcname: str +) -> None: + """Smoothes aparc. + + Parameters + ---------- + insurfname : str + Suface filepath and name of source + inaparcname : str + Annotation filepath and name of source + incortexname : str + Label filepath and name of source + outaparcname : str + Suface filepath and name of destination + + """ + # read input files + print("Reading in surface: {} ...".format(insurfname)) + surf = fs.read_geometry(insurfname, read_metadata=True) + print("Reading in annotation: {} ...".format(inaparcname)) + aparc = fs.read_annot(inaparcname) + print("Reading in cortex label: {} ...".format(incortexname)) + cortex = fs.read_label(incortexname) + # set labels (n) and triangles (n x 3) + labels = aparc[0] + slabels = smooth_aparc(surf, labels, cortex) print("Outputting fixed annot: {}".format(outaparcname)) - fs.write_annot(outaparcname, labels, aparc[1], aparc[2]) + fs.write_annot(outaparcname, slabels, aparc[1], aparc[2]) if __name__ == "__main__": # Command Line options are error checking done here options = options_parse() - smooth_aparc(options.insurf, options.inaparc, options.incort, options.outaparc) + smooth_aparc_files(options.insurf, options.inaparc, options.incort, options.outaparc) sys.exit(0) From 2652a38ca4b95d7f28cf862ccc77cb151dfc9bd2 Mon Sep 17 00:00:00 2001 From: Martin Reuter Date: Mon, 18 Dec 2023 19:02:52 +0100 Subject: [PATCH 2/6] switch to new sample_parc --- recon_surf/recon-surf.sh | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/recon_surf/recon-surf.sh b/recon_surf/recon-surf.sh index 36d31cad..143bd18a 100755 --- a/recon_surf/recon-surf.sh +++ b/recon_surf/recon-surf.sh @@ -818,13 +818,13 @@ echo "echo \" \"" |& tee -a $CMDF # this is dangerous, as some cortices could be < 0.6 mm, but then there is no volume label probably anyway. # Also note that currently we cannot mask non-cortex regions here, should be done in mris_anatomical stats later # the smoothing helps - cmd="mris_sample_parc -ct $FREESURFER_HOME/average/colortable_desikan_killiany.txt -file ${binpath}$hemi.DKTatlaslookup.txt -projmm 0.6 -f 5 -surf white.preaparc $subject $hemi aparc.DKTatlas+aseg.orig.mgz aparc.DKTatlas.mapped.prefix.annot" + #cmd="mris_sample_parc -ct $FREESURFER_HOME/average/colortable_desikan_killiany.txt -file ${binpath}$hemi.DKTatlaslookup.txt -projmm 0.6 -f 5 -surf white.preaparc $subject $hemi aparc.DKTatlas+aseg.orig.mgz aparc.DKTatlas.mapped.prefix.annot" + #RunIt "$cmd" $LF $CMDF + #cmd="$python ${binpath}smooth_aparc.py --insurf $sdir/$hemi.white.preaparc --inaparc $ldir/$hemi.aparc.DKTatlas.mapped.prefix.annot --incort $ldir/$hemi.cortex.label --outaparc $ldir/$hemi.aparc.DKTatlas.mapped.annot" + #RunIt "$cmd" $LF $CMDF + cmd="$python ${binpath}sample_parc.py --inseg $mdir/aparc.DKTatlas+aseg.orig.mgz --insurf $sdir/$hemi.white.preaparc --incort $ldir/$hemi.cortex.label --outaparc $ldir/$hemi.aparc.DKTatlas.mapped.annot --surflut ${binpath}$hemi.DKTatlaslookup.txt --seglut ${binpath}DKTatlaslookup.txt --projmm 0.6 --radius 3" RunIt "$cmd" $LF $CMDF - cmd="$python ${binpath}smooth_aparc.py --insurf $sdir/$hemi.white.preaparc --inaparc $ldir/$hemi.aparc.DKTatlas.mapped.prefix.annot --incort $ldir/$hemi.cortex.label --outaparc $ldir/$hemi.aparc.DKTatlas.mapped.annot" - RunIt "$cmd" $LF $CMDF - - # if we segment with FS or if surface registration is requested do it here: if [ "$fsaparc" == "1" ] || [ "$fssurfreg" == "1" ] ; then echo "echo \" \"" |& tee -a $CMDF From 643cfeecd7b3d1a2ce269d635674d39384eedf84 Mon Sep 17 00:00:00 2001 From: Martin Reuter Date: Mon, 18 Dec 2023 19:03:24 +0100 Subject: [PATCH 3/6] remove rewrite_mc_surface --- recon_surf/recon-surf.sh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/recon_surf/recon-surf.sh b/recon_surf/recon-surf.sh index 143bd18a..001a7b96 100755 --- a/recon_surf/recon-surf.sh +++ b/recon_surf/recon-surf.sh @@ -726,8 +726,8 @@ else RunIt "$cmd" $LF $CMDF # Rewrite surface orig.nofix to fix vertex locs bug (scannerRAS instead of surfaceRAS set with mc) - cmd="$python ${binpath}rewrite_mc_surface.py --input $outmesh --output $outmesh --filename_pretess $mdir/filled-pretess$hemivalue.mgz" - RunIt "$cmd" $LF $CMDF + #cmd="$python ${binpath}rewrite_mc_surface.py --input $outmesh --output $outmesh --filename_pretess $mdir/filled-pretess$hemivalue.mgz" + #RunIt "$cmd" $LF $CMDF # Check if the surfaceRAS was correctly set and exit otherwise (sanity check in case nibabel changes their default header behaviour) cmd="mris_info $outmesh | tr -s ' ' | grep -q 'vertex locs : surfaceRAS'" From 599ab809b90ae5e0fb17b006ee8c897c9ab3a619 Mon Sep 17 00:00:00 2001 From: Martin Reuter Date: Fri, 22 Dec 2023 13:03:27 +0100 Subject: [PATCH 4/6] faster mapping, more smoothing --- recon_surf/recon-surf.sh | 2 +- recon_surf/sample_parc.py | 32 +++++++++++++++++--------------- recon_surf/smooth_aparc.py | 15 +++++++++------ 3 files changed, 27 insertions(+), 22 deletions(-) diff --git a/recon_surf/recon-surf.sh b/recon_surf/recon-surf.sh index 001a7b96..48359d7e 100755 --- a/recon_surf/recon-surf.sh +++ b/recon_surf/recon-surf.sh @@ -822,7 +822,7 @@ echo "echo \" \"" |& tee -a $CMDF #RunIt "$cmd" $LF $CMDF #cmd="$python ${binpath}smooth_aparc.py --insurf $sdir/$hemi.white.preaparc --inaparc $ldir/$hemi.aparc.DKTatlas.mapped.prefix.annot --incort $ldir/$hemi.cortex.label --outaparc $ldir/$hemi.aparc.DKTatlas.mapped.annot" #RunIt "$cmd" $LF $CMDF - cmd="$python ${binpath}sample_parc.py --inseg $mdir/aparc.DKTatlas+aseg.orig.mgz --insurf $sdir/$hemi.white.preaparc --incort $ldir/$hemi.cortex.label --outaparc $ldir/$hemi.aparc.DKTatlas.mapped.annot --surflut ${binpath}$hemi.DKTatlaslookup.txt --seglut ${binpath}DKTatlaslookup.txt --projmm 0.6 --radius 3" + cmd="$python ${binpath}sample_parc.py --inseg $mdir/aparc.DKTatlas+aseg.orig.mgz --insurf $sdir/$hemi.white.preaparc --incort $ldir/$hemi.cortex.label --outaparc $ldir/$hemi.aparc.DKTatlas.mapped.annot --seglut ${binpath}$hemi.DKTatlaslookup.txt --surflut ${binpath}DKTatlaslookup.txt --projmm 0.6 --radius 2" RunIt "$cmd" $LF $CMDF # if we segment with FS or if surface registration is requested do it here: diff --git a/recon_surf/sample_parc.py b/recon_surf/sample_parc.py index 5bb77918..6b825857 100644 --- a/recon_surf/sample_parc.py +++ b/recon_surf/sample_parc.py @@ -22,7 +22,7 @@ import numpy as np from numpy import typing as npt import nibabel.freesurfer.io as fs -from nibabel import load as nibload +import nibabel as nib from scipy import sparse from lapy import TriaMesh from smooth_aparc import smooth_aparc @@ -201,7 +201,7 @@ def sample_img(surf, img, cortex=None, projmm=0.0, radius=None): if isinstance(surf, str): surf = fs.read_geometry(surf, read_metadata=True) if isinstance(img, str): - img = nibload(img) + img = nib.load(img) if isinstance(cortex, str): cortex = fs.read_label(cortex) nvert = surf[0].shape[0] @@ -241,7 +241,7 @@ def sample_img(surf, img, cortex=None, projmm=0.0, radius=None): return samplesfull # here we need to do the hard work of searching in a windows # for non-zero samples - print("sample_img: found {} zero samples, searching radius ...",zeros.size) + print("sample_img: found {} zero samples, searching radius ...".format(zeros.size)) z_nn = x_nn[zeros] z_samples = sample_nearest_nonzero(img, z_nn, radius=radius) samples_nn[zeros] = z_samples @@ -272,20 +272,16 @@ def replace_labels(img_labels, img_lut, surf_lut): Names of label regions. """ surflut = np.loadtxt(surf_lut, usecols=(0,2,3,4,5), dtype="int") - surfids = surflut[:,0] + surf_ids = surflut[:,0] surf_ctab = surflut[:,1:5] surf_names = np.loadtxt(surf_lut, usecols=(1), dtype="str") imglut = np.loadtxt(img_lut, usecols=(0,2,3,4,5), dtype="int") + img_ids = imglut[:,0] img_names = np.loadtxt(img_lut, usecols=(1), dtype="str") assert (np.all(img_names == surf_names)), "Label names in the LUTs do not agree!" - # find labels that have no replacement and map to 0=unknonw - missing = np.logical_not(np.isin(img_labels, surfids)) - img_labels = img_labels.copy() - img_labels[missing] = 0 - # translate the remaining labels - newids = imglut[:,0] - d = dict(zip(surfids, newids)) - surf_labels=np.asarray([d.get(e, e) for e in img_labels]) + lut = np.zeros((img_labels.max()+1,), dtype = img_labels.dtype) + lut[img_ids] = surf_ids + surf_labels = lut[img_labels] return surf_labels, surf_ctab, surf_names @@ -299,7 +295,7 @@ def sample_parc (surf, seg, imglut, surflut, outaparc, cortex=None, projmm=0.0, surf[0] is the np.array of (n, 3) vertex coordinates and surf[1] is the np.array of (m, 3) triangle indices. If type is str, read surface from file. - img : nibabel.image | str + seg : nibabel.image | str Image to sample. If type is str, read image from file. imglut : str @@ -320,8 +316,14 @@ def sample_parc (surf, seg, imglut, surflut, outaparc, cortex=None, projmm=0.0, cortex = fs.read_label(cortex) if isinstance(surf, str): surf = fs.read_geometry(surf, read_metadata=True) - samples = sample_img(surf, seg, cortex, projmm, radius) - surfsamples, surfctab, surfnames = replace_labels(samples, imglut, surflut) + if isinstance(seg, str): + seg = nib.load(seg) + # get rid of unknown labels first and translate the rest (avoids too much filling) + segdata, surfctab, surfnames = replace_labels(np.asarray(seg.dataobj), imglut, surflut) + # create img with new data (needed by sample img) + seg2 = nib.MGHImage(segdata, seg.affine, seg.header) + + surfsamples = sample_img(surf, seg2, cortex, projmm, radius) smooths = smooth_aparc(surf, surfsamples, cortex) fs.write_annot(outaparc, smooths, ctab=surfctab, names=surfnames) diff --git a/recon_surf/smooth_aparc.py b/recon_surf/smooth_aparc.py index 89f91163..5b8f3032 100644 --- a/recon_surf/smooth_aparc.py +++ b/recon_surf/smooth_aparc.py @@ -91,7 +91,7 @@ def get_adjM(trias: npt.NDArray, n: int): Returns ------- adjM : np.ndarray (bool) shape (n,n) - Symmetric adjacency matrix, true corresponds to an edge. + Symmetric sparse CSR adjacency matrix, true corresponds to an edge. """ I = trias @@ -132,7 +132,7 @@ def mode_filter( fillonlylabel = None, novote: npt.ArrayLike = [] ) -> npt.NDArray[int]: - """Apply mode filter (smoothing) to int labels on mesh vertices. + """Apply mode filter (smoothing) to integer labels on mesh vertices. Parameters ---------- @@ -142,11 +142,11 @@ def mode_filter( identity to the adjacency matrix so that each vertex is included in its own vote. labels : npt.NDArray[int] - List of integer labels. + List of integer labels at each vertex of the mesh. fillonlylabel : int - Label to fill exclusively. Defaults to None. + Label to fill exclusively. Defaults to None to smooth all labels. novote : npt.ArrayLike - Entries that should not vote. Defaults to []. + Label ids that should not vote. Defaults to []. Returns ------- @@ -347,7 +347,10 @@ def smooth_aparc(surf, labels, cortex = None): idssize = ids.size counter += 1 # SMOOTH other labels (first with wider kernel then again fine-tune): - labels = mode_filter(adjM * adjM, labels) + adjM2 = adjM * adjM + adjM6 = adjM2 * adjM2 * adjM2 + labels = mode_filter(adjM6, labels) + labels = mode_filter(adjM2, labels) labels = mode_filter(adjM, labels) # set labels outside cortex to -1 labels[~mask] = -1 From 6d4d212a0a6751acc36b35db57182551b5e9aa72 Mon Sep 17 00:00:00 2001 From: Martin Reuter Date: Fri, 22 Dec 2023 15:11:47 +0100 Subject: [PATCH 5/6] avoid some randomness for different radii --- recon_surf/sample_parc.py | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/recon_surf/sample_parc.py b/recon_surf/sample_parc.py index 6b825857..b08849b5 100644 --- a/recon_surf/sample_parc.py +++ b/recon_surf/sample_parc.py @@ -138,13 +138,22 @@ def sample_nearest_nonzero(img, vox_coords, radius=3.0): ri = np.floor(rvox).astype(int) l = np.arange(-ri,ri+1) xv, yv, zv = np.meshgrid(l, l, l) + # modify distances slightly, to avoid randomness when + # sorting with different radius values for voxels that otherwise + # have the same distance to center + xvd = xv+0.001 + yvd = yv+0.002 + zvd = zv+0.003 + ddm = np.sqrt(xvd*xvd + yvd*yvd + zvd*zvd).flatten() + # also compute correct distance for ball mask below dd = np.sqrt(xv*xv + yv*yv + zv*zv).flatten() + ddball = dd<=rvox - # flatten and keep only sphere with radius - xv = xv.flatten()[dd<=rvox] - yv = yv.flatten()[dd<=rvox] - zv = zv.flatten()[dd<=rvox] - dd = dd[dd<=rvox] + # flatten and keep only ball with radius + xv = xv.flatten()[ddball] + yv = yv.flatten()[ddball] + zv = zv.flatten()[ddball] + ddm = ddm[ddball] # stack to get offset vectors offsets = np.column_stack((xv, yv, zv)) @@ -152,8 +161,7 @@ def sample_nearest_nonzero(img, vox_coords, radius=3.0): # sort offsets according to distance # Note: we keep the first zero voxel so we can later # determine if all voxels are zero with the argmax trick - sortidx = np.argsort(dd) - dd = dd[sortidx] + sortidx = np.argsort(ddm) offsets = offsets[sortidx,:] # reshape and tile to add to list of coords From 7cfff9fe53f7533673e1dc00f3fbffaff23fa04b Mon Sep 17 00:00:00 2001 From: Martin Reuter Date: Sat, 13 Jan 2024 13:18:06 +0100 Subject: [PATCH 6/6] add island finding and filling --- recon_surf/sample_parc.py | 86 +++++++++++++++++++++++++++++++++++--- recon_surf/smooth_aparc.py | 11 +++-- 2 files changed, 88 insertions(+), 9 deletions(-) diff --git a/recon_surf/sample_parc.py b/recon_surf/sample_parc.py index b08849b5..322fa2e7 100644 --- a/recon_surf/sample_parc.py +++ b/recon_surf/sample_parc.py @@ -24,6 +24,7 @@ import nibabel.freesurfer.io as fs import nibabel as nib from scipy import sparse +from scipy.sparse.csgraph import connected_components from lapy import TriaMesh from smooth_aparc import smooth_aparc @@ -32,9 +33,9 @@ Script to sample labels from image to surface and clean up. USAGE: -sample_parc --inseg --insurf --incort - --seglut --surflut --outaparc - --projmm --radius +sample_parc --inseg --insurf --incort + --seglut --surflut --outaparc + --projmm --radius Dependencies: @@ -99,7 +100,73 @@ def options_parse(): return options +def construct_adj_cluster(tria, annot): + """Adjacency matrix of edges from same annotation label only. + Operates only on triangles and removes edges that cross annotation + label boundaries. + + Returns + ------- + csc_matrix + The non-directed adjacency matrix + will be symmetric. Each inner edge (i,j) will have + the number of triangles that contain this edge. + Inner edges usually 2, boundary edges 1. Higher + numbers can occur when there are non-manifold triangles. + The sparse matrix can be binarized via: + adj.data = np.ones(adj.data.shape). + """ + t0 = tria[:, 0] + t1 = tria[:, 1] + t2 = tria[:, 2] + i = np.column_stack((t0, t1, t1, t2, t2, t0)).reshape(-1) + j = np.column_stack((t1, t0, t2, t1, t0, t2)).reshape(-1) + ia = annot[i] + ja = annot[j] + keep_edges = (ia == ja) + i = i[keep_edges] + j = j[keep_edges] + dat = np.ones(i.shape) + n = annot.shape[0] + return sparse.csc_matrix((dat, (i, j)), shape=(n, n)) + +def find_all_islands(surf, annot): + """Find vertices in disconnected islands for all labels in surface annotation. + + Parameters + ---------- + surf : tuple + Surface as returned by nibabel fs.read_geometry, where: + surf[0] is the np.array of (n, 3) vertex coordinates and + surf[1] is the np.array of (m, 3) triangle indices. + annot : np.ndarray + Annotation as an int array of (n,) with label ids for each vertex. + This is for example the first element of the tupel returned by + nibabel fs.read_annot. + + Returns + ------- + vidx : np.ndarray (i,) + Arrray listing vertex indices of island vertices, empty if no islands + (components disconnetcted from largest label region) are found. + """ + # construct adjaceny matrix without edges across regions: + adjM = construct_adj_cluster(surf[1], annot) + # compute disconnected components + n_comp, labels = connected_components(csgraph=adjM, directed=False, return_labels=True) + # for each label, get islands that are not connected to main component + lids = np.unique(annot) + vidx = np.array([], dtype = np.int32) + for lid in lids: + ll = labels[annot==lid] + lidx = np.arange(labels.size)[annot==lid] + lmax = np.bincount(ll).argmax() + v = lidx[(ll != lmax)] + if v.size > 0: + print("Found disconnected islands ({} vertices total) for label {}!".format(v.size, lid)) + vidx = np.concatenate((vidx,v)) + return vidx def sample_nearest_nonzero(img, vox_coords, radius=3.0): """Sample closest non-zero value in a ball of radius around vox_coords. @@ -294,7 +361,7 @@ def replace_labels(img_labels, img_lut, surf_lut): def sample_parc (surf, seg, imglut, surflut, outaparc, cortex=None, projmm=0.0, radius=None): - """Replace image labels with corresponding surface labels or unknown. + """Sample cortical GM labels from image to surface and smooth. Parameters ---------- @@ -326,13 +393,20 @@ def sample_parc (surf, seg, imglut, surflut, outaparc, cortex=None, projmm=0.0, surf = fs.read_geometry(surf, read_metadata=True) if isinstance(seg, str): seg = nib.load(seg) - # get rid of unknown labels first and translate the rest (avoids too much filling) + # get rid of unknown labels first and translate the rest (avoids too much filling + # later as sampling will search around sample point if label is zero) segdata, surfctab, surfnames = replace_labels(np.asarray(seg.dataobj), imglut, surflut) # create img with new data (needed by sample img) seg2 = nib.MGHImage(segdata, seg.affine, seg.header) - + # sample from image to surface (and search if zero label) surfsamples = sample_img(surf, seg2, cortex, projmm, radius) + # find label islands + vidx = find_all_islands(surf, surfsamples) + # set islands to zero (to ensure they get smoothed away later) + surfsamples[vidx] = 0 + # smooth boundaries and remove islands inside cortex region smooths = smooth_aparc(surf, surfsamples, cortex) + # write annotation fs.write_annot(outaparc, smooths, ctab=surfctab, names=surfnames) diff --git a/recon_surf/smooth_aparc.py b/recon_surf/smooth_aparc.py index 5b8f3032..23fcd89b 100644 --- a/recon_surf/smooth_aparc.py +++ b/recon_surf/smooth_aparc.py @@ -250,7 +250,12 @@ def mode_filter( def smooth_aparc(surf, labels, cortex = None): - """Smoothes aparc. + """Smoothes aparc and fills holes. + + First all labels with 0 and -1 unside cortex are filled via repeated + mode filtering, then all labels are smoothed first with a wider and + then with smaller filters to produce smooth label boundaries. Labels + outside cortex are set to -1 at the end. Parameters ---------- @@ -348,8 +353,8 @@ def smooth_aparc(surf, labels, cortex = None): counter += 1 # SMOOTH other labels (first with wider kernel then again fine-tune): adjM2 = adjM * adjM - adjM6 = adjM2 * adjM2 * adjM2 - labels = mode_filter(adjM6, labels) + adjM4 = adjM2 * adjM2 + labels = mode_filter(adjM4, labels) labels = mode_filter(adjM2, labels) labels = mode_filter(adjM, labels) # set labels outside cortex to -1