From 0319664a2af9576b48b19305f9062dc628edf376 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Wed, 9 Sep 2015 11:06:26 -0700 Subject: [PATCH 01/35] added some initial documentation on the Aligner module. --- docs/align.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/align.rst b/docs/align.rst index 804cd46b9..a54c0d0dd 100644 --- a/docs/align.rst +++ b/docs/align.rst @@ -1,3 +1,5 @@ Aligner ======= +This module contains functions for manual and automatic alignment of brain images and cortical surfaces. For each transform, it saves a transform matrix, which maps pixels to voxels. +The automatic() function does epi-to-anat registration using FLIRT with BBR, then inverts this with Transform.from_fsl() From f244c6044ea07eec80fa94a3d81cb040326234ea Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Wed, 23 Sep 2015 12:16:35 -0700 Subject: [PATCH 02/35] current state of the align_to_mni function, which has all the components in place but is still waiting for WarpMapper and its final form. --- cortex/align.py | 90 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 90 insertions(+) diff --git a/cortex/align.py b/cortex/align.py index 2a49a2435..8e6181072 100644 --- a/cortex/align.py +++ b/cortex/align.py @@ -152,6 +152,96 @@ def automatic(subject, xfmname, reference, noclean=False, bbrtype="signed"): return retval +def anat_to_mni(subject, xfmname, noclean=False): + """Create an automatic alignment of an anatomical image to the MNI standard. + + If `noclean`, intermediate files will not be removed from /tmp. The `reference` image and resulting + transform called `xfmname` will be automatically stored in the database. + + Parameters + ---------- + subject : str + Subject identifier. + xfmname : str + String identifying the transform to be created. + anatimg : str + Path to a nibabel-readable image that will be used as the reference for this transform. + This should be a 3D anatomical volume. + noclean : bool, optional + If True intermediate files will not be removed from /tmp (this is useful for debugging things), + and the returned value will be the name of the temp directory. Default False. + + Returns + ------- + Nothing unless `noclean` is True. + """ + + import shlex + import shutil + import tempfile + import subprocess as sp + + from .database import db + from .xfm import Transform + from .options import config + + fsl_prefix = config.get("basic", "fsl_prefix") + schfile = os.path.join(os.path.split(os.path.abspath(__file__))[0], "bbr.sch") + + print('anat_to_mni, subject: %s, xfmname: %s' % (subject, xfmname)) + + try: + raw_anat = db.get_anat(subject, type='raw').get_filename() + bet_anat = db.get_anat(subject, type='brainmask').get_filename() + anat_dir = os.path.dirname(raw_anat) + odir = '/tmp' + + # stem for the reoriented-into-MNI anatomical images (required by FLIRT/FNIRT) + reorient_anat = 'reorient_anat' + reorient_cmd = '{fslpre}fslreorient2std {raw_anat} {adir}/{ra_raw}'.format(fslpre=fsl_prefix,raw_anat=raw_anat, adir=odir, ra_raw=reorient_anat) + print('Reorienting anatomicals using fslreorient2std') + if sp.call(reorient_cmd, shell=True) != 0: + raise IOError('Error calling fslreorient2std on raw anatomical') + reorient_cmd = '{fslpre}fslreorient2std {bet_anat} {adir}/{ra_raw}_brain'.format(fslpre=fsl_prefix,bet_anat=bet_anat, adir=odir, ra_raw=reorient_anat) + if sp.call(reorient_cmd, shell=True) != 0: + raise IOError('Error calling fslreorient2std on brain-extracted anatomical') + + fsldir = os.environ['FSLDIR'] + standard = '%s/data/standard/MNI152_T1_1mm'%fsldir + bet_standard = '%s_brain'%standard + standardmask = '%s_mask_dil'%bet_standard + cout = 'anat2mni' #stem of the filenames of the transform estimates + + # initial affine anatomical-to-standard registration using FLIRT. required, as the output xfm is used as a start by FNIRT. + flirt_cmd = '{fslpre}flirt -in {adir}/{ra_raw}_brain -ref {bet_standard} -dof 6 -omat /tmp/{cout}_flirt' + flirt_cmd = flirt_cmd.format(fslpre=fsl_prefix, ra_raw=reorient_anat, bet_standard=bet_standard, adir=odir, cout=cout) + print('Running FLIRT to estimate initial affine transform') + #if sp.call(flirt_cmd, shell=True) != 0: + # raise IOError('Error calling FLIRT with command: %s' % flirt_cmd) + + # FNIRT anat-to-mni transform estimation cmd (does not apply any transform, but generates estimate [cout]) + cmd = '{fslpre}fnirt --in={ad}/{ra_raw} --ref={standard} --refmask={standardmask} --aff=/tmp/{cout}_flirt --cout={anat_dir}/{cout}_fnirt --fout=/tmp/{cout}_field --config=T1_2_MNI152_2mm' + cmd = cmd.format(fslpre=fsl_prefix, ra_raw=reorient_anat, standard=standard, standardmask=standardmask, ad=odir, anat_dir=anat_dir, cout=cout) + print('Running FNIRT to estimate transform... this can take a while') + #if sp.call(cmd, shell=True) != 0: + # raise IOError('Error calling fnirt') + + # we now have, in /tmp/cout_fnirt, the warp estimate that should be passed to img2stdcoord. + # let's get all vertex coordinates + cfile = '/tmp/fid_coords' + cfile_warped = '/tmp/mni_coords' + [pts, polys] = db.get_surf(subject,"fiducial",merge=True) + np.savetxt(cfile, pts, fmt='%g') + + xfm_cmd = 'cat {coordfile} | {fslpre}img2stdcoord -mm -img {ad}/{ra_raw} -std {standard} -warp {anat_dir}/{cout}_fnirt > {cfile_warped}' + xfm_cmd = xfm_cmd.format(coordfile=cfile, fslpre=fsl_prefix, ra_raw=reorient_anat, standard=standard, ad=odir, anat_dir=anat_dir, cout=cout, cfile_warped=cfile_warped) + + print('raw anatomical: %s\nbet anatomical: %s\nflirt cmd:%s\nfnirt cmd: %s\nxfm cmd: %s' % (raw_anat,bet_anat,flirt_cmd,cmd,xfm_cmd)) + + finally: + pass + + def autotweak(subject, xfmname): """Tweak an alignment using the FLIRT boundary-based alignment (BBR) from FSL. Ideally this function should actually use a limited search range, but it doesn't. From 57f7bccdb0eb429f40627c4911b1504d15f0a116 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Wed, 23 Sep 2015 12:22:06 -0700 Subject: [PATCH 03/35] stub WarpMapper and modified list of mappers to include it. --- cortex/mapper/__init__.py | 5 +++-- cortex/mapper/warp.py | 32 ++++++++++++++++++++++++++++++++ 2 files changed, 35 insertions(+), 2 deletions(-) create mode 100644 cortex/mapper/warp.py diff --git a/cortex/mapper/__init__.py b/cortex/mapper/__init__.py index da88908aa..0f16cda5d 100644 --- a/cortex/mapper/__init__.py +++ b/cortex/mapper/__init__.py @@ -10,7 +10,7 @@ def get_mapper(subject, xfmname, type='nearest', recache=False, **kwargs): from ..database import db - from . import point, patch, volume, line + from . import point, patch, volume, line, warp mapcls = dict( nearest=point.PointNN, @@ -21,7 +21,8 @@ def get_mapper(subject, xfmname, type='nearest', recache=False, **kwargs): const_patch_trilin=patch.ConstPatchTrilin, const_patch_lanczos=patch.ConstPatchLanczos, line_nearest=line.LineNN, - line_trilinear=line.LineTrilin) + line_trilinear=line.LineTrilin, + warp=warp.WarpMapper) Map = mapcls[type] ptype = Map.__name__.lower() kwds ='_'.join(['%s%s'%(k,str(v)) for k, v in list(kwargs.items())]) diff --git a/cortex/mapper/warp.py b/cortex/mapper/warp.py new file mode 100644 index 000000000..453d04a6a --- /dev/null +++ b/cortex/mapper/warp.py @@ -0,0 +1,32 @@ +import numpy as np +from scipy import sparse + +from . import Mapper +from . import samplers + +class WarpMapper(Mapper): + sampler = staticmethod(samplers.nearest) + @classmethod + def _getmask(cls, coords, polys, shape, **kwargs): + valid = np.unique(polys) + mcoords = np.nan * np.ones_like(coords) + mcoords[valid] = coords[valid] + i, j, data = cls.sampler(mcoords, shape, **kwargs) + csrshape = len(coords), np.prod(shape) + return sparse.csr_matrix((data, np.array([i, j])), shape=csrshape) + + @classmethod + def _cache(cls, filename, subject, xfmname, **kwargs): + print('Caching mapper...') + from ..database import db + masks = [] + # the next line is where the big changes will have to start... + xfm = db.get_xfm(subject, xfmname, xfmtype='coord') + fid = db.get_surf(subject, 'fiducial', merge=False, nudge=False) + flat = db.get_surf(subject, 'flat', merge=False, nudge=False) + + for (pts, _), (_, polys) in zip(fid, flat): + masks.append(cls._getmask(xfm(pts), polys, xfm.shape, **kwargs)) + + _savecache(filename, masks[0], masks[1], xfm.shape) + return cls(masks[0], masks[1], xfm.shape) \ No newline at end of file From 0b961c385c8ff1b2cdfbecb15064befc0c10df74 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Tue, 29 Sep 2015 20:21:22 -0700 Subject: [PATCH 04/35] WarpMapper class to load up warp field and sample it, currently using nearest-neighbor --- cortex/mapper/warp.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/cortex/mapper/warp.py b/cortex/mapper/warp.py index 453d04a6a..bd93a77be 100644 --- a/cortex/mapper/warp.py +++ b/cortex/mapper/warp.py @@ -1,7 +1,7 @@ import numpy as np from scipy import sparse -from . import Mapper +from . import Mapper, _savecache from . import samplers class WarpMapper(Mapper): @@ -17,11 +17,19 @@ def _getmask(cls, coords, polys, shape, **kwargs): @classmethod def _cache(cls, filename, subject, xfmname, **kwargs): - print('Caching mapper...') from ..database import db + from ..xfm import Transform + + print('Caching warpmapper...') + print filename, subject, xfmname + masks = [] # the next line is where the big changes will have to start... - xfm = db.get_xfm(subject, xfmname, xfmtype='coord') + # xfm = db.get_xfm(subject, xfmname, xfmtype='coord') + re_anat = db.get_anat(subject,'reorient_anat') + xfm = Transform(np.linalg.inv(re_anat.get_affine()),re_anat) + warp = db.get_anat(subject,'mni2anat_field') + warpfield = warp.get_data() fid = db.get_surf(subject, 'fiducial', merge=False, nudge=False) flat = db.get_surf(subject, 'flat', merge=False, nudge=False) From 54de55093b399ad3fd4cd9a17c833cd93c5cf7b8 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Wed, 30 Sep 2015 11:18:55 -0700 Subject: [PATCH 05/35] changes to anat_to_mni function, which is still incomplete. the reoriented anatomicals are saved in the anat_dir, not /tmp, and the commands are printed for debugging use. in future the anatomicals may be saved to /tmp again, since the reoriented versions are not necessary for getting the correct coordinates in the identity transform --- cortex/align.py | 26 ++++++++++++++++---------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/cortex/align.py b/cortex/align.py index 8e6181072..fe16d7d5e 100644 --- a/cortex/align.py +++ b/cortex/align.py @@ -193,47 +193,53 @@ def anat_to_mni(subject, xfmname, noclean=False): try: raw_anat = db.get_anat(subject, type='raw').get_filename() bet_anat = db.get_anat(subject, type='brainmask').get_filename() + betmask_anat = db.get_anat(subject, type='brainmask_mask').get_filename() anat_dir = os.path.dirname(raw_anat) - odir = '/tmp' + odir = anat_dir # stem for the reoriented-into-MNI anatomical images (required by FLIRT/FNIRT) reorient_anat = 'reorient_anat' reorient_cmd = '{fslpre}fslreorient2std {raw_anat} {adir}/{ra_raw}'.format(fslpre=fsl_prefix,raw_anat=raw_anat, adir=odir, ra_raw=reorient_anat) - print('Reorienting anatomicals using fslreorient2std') + print('Reorienting anatomicals using fslreorient2std, cmd like: \n%s' % reorient_cmd) if sp.call(reorient_cmd, shell=True) != 0: raise IOError('Error calling fslreorient2std on raw anatomical') reorient_cmd = '{fslpre}fslreorient2std {bet_anat} {adir}/{ra_raw}_brain'.format(fslpre=fsl_prefix,bet_anat=bet_anat, adir=odir, ra_raw=reorient_anat) if sp.call(reorient_cmd, shell=True) != 0: raise IOError('Error calling fslreorient2std on brain-extracted anatomical') + + ra_betmask = reorient_anat + "_brainmask" + reorient_cmd = '{fslpre}fslreorient2std {bet_anat} {adir}/{ra_betmask}'.format(fslpre=fsl_prefix,bet_anat=betmask_anat, adir=odir, ra_betmask=ra_betmask) + if sp.call(reorient_cmd, shell=True) != 0: + raise IOError('Error calling fslreorient2std on brain-extracted mask') fsldir = os.environ['FSLDIR'] standard = '%s/data/standard/MNI152_T1_1mm'%fsldir bet_standard = '%s_brain'%standard standardmask = '%s_mask_dil'%bet_standard - cout = 'anat2mni' #stem of the filenames of the transform estimates + cout = 'mni2anat' #stem of the filenames of the transform estimates # initial affine anatomical-to-standard registration using FLIRT. required, as the output xfm is used as a start by FNIRT. - flirt_cmd = '{fslpre}flirt -in {adir}/{ra_raw}_brain -ref {bet_standard} -dof 6 -omat /tmp/{cout}_flirt' + flirt_cmd = '{fslpre}flirt -in {bet_standard} -ref {adir}/{ra_raw}_brain -dof 6 -omat /tmp/{cout}_flirt' flirt_cmd = flirt_cmd.format(fslpre=fsl_prefix, ra_raw=reorient_anat, bet_standard=bet_standard, adir=odir, cout=cout) print('Running FLIRT to estimate initial affine transform') #if sp.call(flirt_cmd, shell=True) != 0: # raise IOError('Error calling FLIRT with command: %s' % flirt_cmd) - # FNIRT anat-to-mni transform estimation cmd (does not apply any transform, but generates estimate [cout]) - cmd = '{fslpre}fnirt --in={ad}/{ra_raw} --ref={standard} --refmask={standardmask} --aff=/tmp/{cout}_flirt --cout={anat_dir}/{cout}_fnirt --fout=/tmp/{cout}_field --config=T1_2_MNI152_2mm' - cmd = cmd.format(fslpre=fsl_prefix, ra_raw=reorient_anat, standard=standard, standardmask=standardmask, ad=odir, anat_dir=anat_dir, cout=cout) + # FNIRT mni-to-anat transform estimation cmd (does not apply any transform, but generates estimate [cout]) + cmd = '{fslpre}fnirt --in={standard} --ref={ad}/{ra_raw} --refmask={ad}/{refmask} --aff=/tmp/{cout}_flirt --cout={anat_dir}/{cout}_fnirt --fout={anat_dir}/{cout}_field --iout=/tmp/mni2anat_iout --config=T1_2_MNI152_2mm' + cmd = cmd.format(fslpre=fsl_prefix, ra_raw=reorient_anat, standard=standard, refmask=ra_betmask, ad=odir, anat_dir=anat_dir, cout=cout) print('Running FNIRT to estimate transform... this can take a while') #if sp.call(cmd, shell=True) != 0: - # raise IOError('Error calling fnirt') + # raise IOError('Error calling fnirt with cmd: %s'%cmd) # we now have, in /tmp/cout_fnirt, the warp estimate that should be passed to img2stdcoord. # let's get all vertex coordinates cfile = '/tmp/fid_coords' cfile_warped = '/tmp/mni_coords' [pts, polys] = db.get_surf(subject,"fiducial",merge=True) - np.savetxt(cfile, pts, fmt='%g') + # np.savetxt(cfile, pts, fmt='%g') - xfm_cmd = 'cat {coordfile} | {fslpre}img2stdcoord -mm -img {ad}/{ra_raw} -std {standard} -warp {anat_dir}/{cout}_fnirt > {cfile_warped}' + xfm_cmd = 'cat {coordfile} | {fslpre}img2stdcoord -mm -std {ad}/{ra_raw} -img {standard} -warp {anat_dir}/{cout}_fnirt > {cfile_warped}' xfm_cmd = xfm_cmd.format(coordfile=cfile, fslpre=fsl_prefix, ra_raw=reorient_anat, standard=standard, ad=odir, anat_dir=anat_dir, cout=cout, cfile_warped=cfile_warped) print('raw anatomical: %s\nbet anatomical: %s\nflirt cmd:%s\nfnirt cmd: %s\nxfm cmd: %s' % (raw_anat,bet_anat,flirt_cmd,cmd,xfm_cmd)) From 7bde7781833918db717268ad91392cb9b7c149df Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Thu, 1 Oct 2015 14:13:21 -0700 Subject: [PATCH 06/35] deleted warpmapper since it's not actually needed and pointmapper will do, reverted the addition in __init__ also --- cortex/mapper/__init__.py | 3 +-- cortex/mapper/warp.py | 40 --------------------------------------- 2 files changed, 1 insertion(+), 42 deletions(-) delete mode 100644 cortex/mapper/warp.py diff --git a/cortex/mapper/__init__.py b/cortex/mapper/__init__.py index 0f16cda5d..c781aa53e 100644 --- a/cortex/mapper/__init__.py +++ b/cortex/mapper/__init__.py @@ -21,8 +21,7 @@ def get_mapper(subject, xfmname, type='nearest', recache=False, **kwargs): const_patch_trilin=patch.ConstPatchTrilin, const_patch_lanczos=patch.ConstPatchLanczos, line_nearest=line.LineNN, - line_trilinear=line.LineTrilin, - warp=warp.WarpMapper) + line_trilinear=line.LineTrilin) Map = mapcls[type] ptype = Map.__name__.lower() kwds ='_'.join(['%s%s'%(k,str(v)) for k, v in list(kwargs.items())]) diff --git a/cortex/mapper/warp.py b/cortex/mapper/warp.py deleted file mode 100644 index bd93a77be..000000000 --- a/cortex/mapper/warp.py +++ /dev/null @@ -1,40 +0,0 @@ -import numpy as np -from scipy import sparse - -from . import Mapper, _savecache -from . import samplers - -class WarpMapper(Mapper): - sampler = staticmethod(samplers.nearest) - @classmethod - def _getmask(cls, coords, polys, shape, **kwargs): - valid = np.unique(polys) - mcoords = np.nan * np.ones_like(coords) - mcoords[valid] = coords[valid] - i, j, data = cls.sampler(mcoords, shape, **kwargs) - csrshape = len(coords), np.prod(shape) - return sparse.csr_matrix((data, np.array([i, j])), shape=csrshape) - - @classmethod - def _cache(cls, filename, subject, xfmname, **kwargs): - from ..database import db - from ..xfm import Transform - - print('Caching warpmapper...') - print filename, subject, xfmname - - masks = [] - # the next line is where the big changes will have to start... - # xfm = db.get_xfm(subject, xfmname, xfmtype='coord') - re_anat = db.get_anat(subject,'reorient_anat') - xfm = Transform(np.linalg.inv(re_anat.get_affine()),re_anat) - warp = db.get_anat(subject,'mni2anat_field') - warpfield = warp.get_data() - fid = db.get_surf(subject, 'fiducial', merge=False, nudge=False) - flat = db.get_surf(subject, 'flat', merge=False, nudge=False) - - for (pts, _), (_, polys) in zip(fid, flat): - masks.append(cls._getmask(xfm(pts), polys, xfm.shape, **kwargs)) - - _savecache(filename, masks[0], masks[1], xfm.shape) - return cls(masks[0], masks[1], xfm.shape) \ No newline at end of file From 08b339ef28067e88591bf6dfa3c6f64f50faad27 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Wed, 14 Oct 2015 16:27:31 -0700 Subject: [PATCH 07/35] working align function (though some calls are still disabled in this commit). reorients anatomicals, does flirt, feeds this to fnirt, saves the resulting warpfield in the database. then samples the warp field to get a Volume of translations, maps this to vertices, and creates arrays corresponding to the pts array but with the MNI coordinates, not the anatomical ones. saves this as a surfinfo in the database. --- cortex/align.py | 90 +++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 73 insertions(+), 17 deletions(-) diff --git a/cortex/align.py b/cortex/align.py index fe16d7d5e..c96922166 100644 --- a/cortex/align.py +++ b/cortex/align.py @@ -180,10 +180,12 @@ def anat_to_mni(subject, xfmname, noclean=False): import shutil import tempfile import subprocess as sp + import nibabel as nib from .database import db from .xfm import Transform from .options import config + from .dataset import Volume fsl_prefix = config.get("basic", "fsl_prefix") schfile = os.path.join(os.path.split(os.path.abspath(__file__))[0], "bbr.sch") @@ -201,16 +203,16 @@ def anat_to_mni(subject, xfmname, noclean=False): reorient_anat = 'reorient_anat' reorient_cmd = '{fslpre}fslreorient2std {raw_anat} {adir}/{ra_raw}'.format(fslpre=fsl_prefix,raw_anat=raw_anat, adir=odir, ra_raw=reorient_anat) print('Reorienting anatomicals using fslreorient2std, cmd like: \n%s' % reorient_cmd) - if sp.call(reorient_cmd, shell=True) != 0: - raise IOError('Error calling fslreorient2std on raw anatomical') + #if sp.call(reorient_cmd, shell=True) != 0: + # raise IOError('Error calling fslreorient2std on raw anatomical') reorient_cmd = '{fslpre}fslreorient2std {bet_anat} {adir}/{ra_raw}_brain'.format(fslpre=fsl_prefix,bet_anat=bet_anat, adir=odir, ra_raw=reorient_anat) - if sp.call(reorient_cmd, shell=True) != 0: - raise IOError('Error calling fslreorient2std on brain-extracted anatomical') + #if sp.call(reorient_cmd, shell=True) != 0: + # raise IOError('Error calling fslreorient2std on brain-extracted anatomical') ra_betmask = reorient_anat + "_brainmask" reorient_cmd = '{fslpre}fslreorient2std {bet_anat} {adir}/{ra_betmask}'.format(fslpre=fsl_prefix,bet_anat=betmask_anat, adir=odir, ra_betmask=ra_betmask) - if sp.call(reorient_cmd, shell=True) != 0: - raise IOError('Error calling fslreorient2std on brain-extracted mask') + #if sp.call(reorient_cmd, shell=True) != 0: + # raise IOError('Error calling fslreorient2std on brain-extracted mask') fsldir = os.environ['FSLDIR'] standard = '%s/data/standard/MNI152_T1_1mm'%fsldir @@ -221,28 +223,82 @@ def anat_to_mni(subject, xfmname, noclean=False): # initial affine anatomical-to-standard registration using FLIRT. required, as the output xfm is used as a start by FNIRT. flirt_cmd = '{fslpre}flirt -in {bet_standard} -ref {adir}/{ra_raw}_brain -dof 6 -omat /tmp/{cout}_flirt' flirt_cmd = flirt_cmd.format(fslpre=fsl_prefix, ra_raw=reorient_anat, bet_standard=bet_standard, adir=odir, cout=cout) - print('Running FLIRT to estimate initial affine transform') + print('Running FLIRT to estimate initial affine transform with command:\n%s'%flirt_cmd) #if sp.call(flirt_cmd, shell=True) != 0: # raise IOError('Error calling FLIRT with command: %s' % flirt_cmd) # FNIRT mni-to-anat transform estimation cmd (does not apply any transform, but generates estimate [cout]) cmd = '{fslpre}fnirt --in={standard} --ref={ad}/{ra_raw} --refmask={ad}/{refmask} --aff=/tmp/{cout}_flirt --cout={anat_dir}/{cout}_fnirt --fout={anat_dir}/{cout}_field --iout=/tmp/mni2anat_iout --config=T1_2_MNI152_2mm' cmd = cmd.format(fslpre=fsl_prefix, ra_raw=reorient_anat, standard=standard, refmask=ra_betmask, ad=odir, anat_dir=anat_dir, cout=cout) - print('Running FNIRT to estimate transform... this can take a while') + print('Running FNIRT to estimate transform, using the following command... this can take a while:\n%s'%cmd) #if sp.call(cmd, shell=True) != 0: # raise IOError('Error calling fnirt with cmd: %s'%cmd) - # we now have, in /tmp/cout_fnirt, the warp estimate that should be passed to img2stdcoord. - # let's get all vertex coordinates - cfile = '/tmp/fid_coords' - cfile_warped = '/tmp/mni_coords' - [pts, polys] = db.get_surf(subject,"fiducial",merge=True) + [pts, polys] = db.get_surf(subject,"fiducial",merge="True") # np.savetxt(cfile, pts, fmt='%g') - xfm_cmd = 'cat {coordfile} | {fslpre}img2stdcoord -mm -std {ad}/{ra_raw} -img {standard} -warp {anat_dir}/{cout}_fnirt > {cfile_warped}' - xfm_cmd = xfm_cmd.format(coordfile=cfile, fslpre=fsl_prefix, ra_raw=reorient_anat, standard=standard, ad=odir, anat_dir=anat_dir, cout=cout, cfile_warped=cfile_warped) - - print('raw anatomical: %s\nbet anatomical: %s\nflirt cmd:%s\nfnirt cmd: %s\nxfm cmd: %s' % (raw_anat,bet_anat,flirt_cmd,cmd,xfm_cmd)) + #print('raw anatomical: %s\nbet anatomical: %s\nflirt cmd:%s\nfnirt cmd: %s\npts: %s' % (raw_anat,bet_anat,flirt_cmd,cmd,pts)) + + # take the reoriented anatomical, get its affine coord transform, invert this, and save it + reo_xfmnm = 'reorient_inv' + re_anat = db.get_anat(subject,reorient_anat) + reo_xfm = Transform(np.linalg.inv(re_anat.get_affine()),re_anat) + reo_xfm.save(subject,reo_xfmnm,"coord") + + # get the reoriented anatomical's qform and its inverse, they will be needed later + aqf = re_anat.get_qform() + aqfinv = np.linalg.inv(aqf) + + # load the warp field data as a volume + warp = db.get_anat(subject,'%s_field'%cout) + wd = warp.get_data() + # need in (t,z,y,x) order + wd = np.swapaxes(wd,0,3) # x <--> t + wd = np.swapaxes(wd,1,2) # y <--> z + wv = Volume(wd,subject,reo_xfmnm) + + # now do the mapping! this gets the warp field values at the corresponding points + # (uses fiducial surface by default) + warpvd = wv.map(projection="nearest") + + # reshape into something sensible + warpverts_L = [vs for vs in np.swapaxes(warpvd.left,0,1)] + warpverts_R = [vs for vs in np.swapaxes(warpvd.right,0,1)] + warpverts_ordered = np.concatenate((warpverts_L, warpverts_R)) + + # append 1s for matrix multiplication (coordinate transformation) + o = np.ones((len(pts),1)) + pad_pts = np.append(pts, o, axis=1) + + # print pts, len(pts), len(pts[0]), warpverts_ordered, len(warpverts_ordered), pad_pts, len(pad_pts), pad_pts[0] + + # transform vertex coords from mm to vox using the anat's qform + voxcoords = [aqfinv.dot(padpt) for padpt in pad_pts] + # add the offsets specified in the warp at those locations (ignoring the 1s here) + mnivoxcoords = [voxcoords[n][:-1] + warpverts_ordered[n] for n in range(len(voxcoords))] + # re-pad for matrix multiplication + pad_mnivox = np.append(mnivoxcoords, o, axis=1) + + # multiply by the standard's qform to recover mm coords + std = nib.load('%s.nii.gz'%standard) + stdqf = std.get_qform() + mni_coords = np.array([stdqf.dot(padmni)[:-1] for padmni in pad_mnivox]) + + # some debug output + # print pts, mni_coords + print pts[0], mni_coords[0] + print len(pts), len(mni_coords) + print type(pts), type(pts[0][0]), type(mni_coords) + + # now split mni_coords into left and right arrays for saving + nverts_L = len(warpverts_L) + print nverts_L + left = mni_coords[:nverts_L] + right = mni_coords[nverts_L:] + print len(left), len(right) + + mni_surfinfo_fn = db.get_paths(subject)['surfinfo'].format(type='mnicoords',opts='') + np.savez(mni_surfinfo_fn,leftpts=left,rightpts=right) finally: pass From 20cb2ae7c332df598bc54ea2490fdd4d6bd8ab3d Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Wed, 14 Oct 2015 16:33:45 -0700 Subject: [PATCH 08/35] modifications to BrainCTM to store MNI data as an additional attribute of the mesh --- cortex/brainctm.py | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/cortex/brainctm.py b/cortex/brainctm.py index c242b1406..29ac9bbcc 100644 --- a/cortex/brainctm.py +++ b/cortex/brainctm.py @@ -61,6 +61,7 @@ def __init__(self, subject, decimate=False): fidpolys = set(tuple(f) for f in polyutils.sort_polys(hemi.polys)) flatpolys = set(tuple(f) for f in polyutils.sort_polys(ptpoly[1])) hemi.aux[np.array(list(fidpolys - flatpolys)).astype(int), 0] = 1 + hemi.mni[np.array(list(fidpolys - flatpolys)).astype(int), 0] = 1 #Find the flatmap limits if fleft is not None: @@ -89,6 +90,12 @@ def addCurvature(self, **kwargs): self.left.aux[:,1] = npz.left self.right.aux[:,1] = npz.right + def addMNI(self, **kwargs): + print('Adding MNI coords...') + npz = db.get_surfinfo(self.subject, type='mnicoords', **kwargs) + self.left.mni[:,:-1] = npz['leftpts'] + self.right.mni[:,:-1] = npz['rightpts'] + def save(self, path, method='mg2', disp_layers=['rois'], extra_disp=None, **kwargs): """Save CTM file for static html display. @@ -105,6 +112,8 @@ def save(self, path, method='mg2', disp_layers=['rois'], extra_disp=None, **kwar svgname = path+".svg" jsname = path+".json" + print("In method save() of BrainCTM, ctmname: %s"%ctmname) + # Save CTM concatenation (lpts, _, _), lbin = self.left.save(method=method, **kwargs) (rpts, _, _), rbin = self.right.save(method=method, **kwargs) @@ -177,9 +186,11 @@ def __init__(self, pts, polys, norms=None): self.flat = None self.surfs = {} self.aux = np.zeros((len(self.ctm), 4)) + self.mni = np.zeros((len(self.ctm), 4)) def addSurf(self, pts, name=None, renorm=True): '''Scales the in-between surfaces to be same scale as fiducial''' + print('In Hemi.addSurf()...') if name is None: name = 'morphTarget%d'%len(self.surfs) @@ -191,6 +202,7 @@ def addSurf(self, pts, name=None, renorm=True): attrib = np.hstack([rnorm, np.zeros((len(rnorm),1))]) self.surfs[name] = attrib + print('Adding attrib %s'%name) self.ctm.addAttrib(attrib, name) def setFlat(self, pts): @@ -199,6 +211,8 @@ def setFlat(self, pts): def save(self, **kwargs): self.ctm.addAttrib(self.aux, 'auxdat') + print('Adding attrib mnicoords, just added auxdat...') + self.ctm.addAttrib(self.mni, 'mnicoords') self.ctm.save(**kwargs) ctm = CTMfile(self.tf.name) @@ -231,6 +245,7 @@ def __init__(self, pts, polys, fpolys, pia=None): self.aux[idxmap[mwidx], 0] = 1 self.mask = mask self.idxmap = idxmap + self.mni = np.zeros((len(self.ctm), 4)) def setFlat(self, pts): super(DecimatedHemi, self).setFlat(pts[self.mask]) @@ -242,8 +257,10 @@ def make_pack(outfile, subj, types=("inflated",), method='raw', level=0, decimate=False, disp_layers=['rois'],extra_disp=None): """Generates a cached CTM file""" + print 'In make_pack, outfile %s'%outfile ctm = BrainCTM(subj, decimate=decimate) ctm.addCurvature() + ctm.addMNI() for name in types: ctm.addSurf(name) @@ -257,6 +274,7 @@ def make_pack(outfile, subj, types=("inflated",), method='raw', level=0, extra_disp=extra_disp) def read_pack(ctmfile): + print('In read_pack, ctmfile: %s'%ctmfile) fname = os.path.splitext(ctmfile)[0] jsfile = json.load(open(fname+".json")) offset = jsfile['offsets'] @@ -275,5 +293,4 @@ def read_pack(ctmfile): ctm = CTMfile(tf.name, "r") pts, polys, norms = ctm.getMesh() meshes.append((pts, polys)) - return meshes From 2cb70a9249bd8f2b525038ce1cf359a0c5b4d9b4 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Wed, 14 Oct 2015 16:39:22 -0700 Subject: [PATCH 09/35] initial mapper documentation --- docs/mapper.rst | 5 +++++ 1 file changed, 5 insertions(+) create mode 100644 docs/mapper.rst diff --git a/docs/mapper.rst b/docs/mapper.rst new file mode 100644 index 000000000..b3e4f0fb5 --- /dev/null +++ b/docs/mapper.rst @@ -0,0 +1,5 @@ +Mapper +====== +This module contains functions and classes that translate between Volumes of data and the corresponding Vertices of a surface using any of several samplers. That is, if you have some data defined on a per-voxel basis (whether it's BOLD or a warp vector doesn't matter), you can create a Volume and use the map() function to sample from that data at the vertices of a surface. And you can specify the sampling method, be it nearest-neighbor, trilinear, sinc, etc... + +Important Classes: \ No newline at end of file From e614145ae74b22032ec9f44c33be4edea5fbbefa Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Fri, 23 Oct 2015 16:22:40 -0700 Subject: [PATCH 10/35] changes to javascript side to display MNI coordinates associated with a picked vertex. --- cortex/webgl/resources/css/mriview.css | 22 ++++++++++++++ cortex/webgl/resources/js/ctm/CTMLoader.js | 6 ++++ cortex/webgl/resources/js/dataset.js | 4 +-- cortex/webgl/resources/js/facepick.js | 34 +++++++++++++++++++++- cortex/webgl/resources/js/mriview.js | 8 +++++ cortex/webgl/template.html | 3 ++ 6 files changed, 74 insertions(+), 3 deletions(-) diff --git a/cortex/webgl/resources/css/mriview.css b/cortex/webgl/resources/css/mriview.css index 826147f4b..de8ad6774 100644 --- a/cortex/webgl/resources/css/mriview.css +++ b/cortex/webgl/resources/css/mriview.css @@ -136,6 +136,28 @@ a:visited { color:white; } background:rgba(255,255,255,.4); font-weight:bold; } + +/***** MNI Coordinate box, top left, below dataname *********/ +#mnibox { + display:none; + position:absolute; + z-index:12; + float:left; + text-shadow:0px 2px 8px black, 0px 1px 8px black; + color:white; + font-size:24pt; + padding:10px; + margin:20px; + margin-top: 90px ; + border-radius:10px; + background:rgba(255,255,255,.4); + font-weight:bold; +} + +#mnicoords, #mnisubmit { + z-index:13; +} + .datadesc { font-size:10pt; font-weight:normal; diff --git a/cortex/webgl/resources/js/ctm/CTMLoader.js b/cortex/webgl/resources/js/ctm/CTMLoader.js index e8a82e12e..72e6c96b4 100644 --- a/cortex/webgl/resources/js/ctm/CTMLoader.js +++ b/cortex/webgl/resources/js/ctm/CTMLoader.js @@ -93,6 +93,8 @@ THREE.CTMLoader.prototype.load = function( url, callback, useWorker, useBuffers, var length = 0; + console.log(url) + xhr.onreadystatechange = function() { if ( xhr.readyState === 4 ) { @@ -200,6 +202,8 @@ THREE.CTMLoader.prototype.createModelBuffers = function ( file, callback, header var Model = function ( ) { + console.log(file); + console.log('createModelBuffers'); var scope = this; scope.materials = []; @@ -234,6 +238,8 @@ THREE.CTMLoader.prototype.createModelBuffers = function ( file, callback, header attrname = file.body.attrMaps[ i ].name; array = file.body.attrMaps[ i ].attr; + console.log(attrname) + console.log(array) if (attrname.slice(0, 11) == "morphTarget") { diff --git a/cortex/webgl/resources/js/dataset.js b/cortex/webgl/resources/js/dataset.js index 45ffa14aa..b7c64e141 100644 --- a/cortex/webgl/resources/js/dataset.js +++ b/cortex/webgl/resources/js/dataset.js @@ -247,7 +247,7 @@ var dataset = (function(module) { volxfm: { type:'m4', value: new THREE.Matrix4() }, data: { type: 't', value: 0, texture: null }, }, - attributes: { flatpos: true, wm:true, auxdat:true, }, + attributes: { flatpos: true, wm:true, auxdat:true, mnicoords:true}, }), targets = { left: new THREE.WebGLRenderTarget(res, res, { @@ -262,7 +262,7 @@ var dataset = (function(module) { }), } var hemi, geom, target, mesh, name, attr; - var names = ["position", "wm", "auxdat", "index"]; + var names = ["position", "wm", "auxdat", "index", "mnicoords"]; var limits = {top:-1000, bottom:1000}; var xfm = shader.uniforms.volxfm.value; xfm.set.apply(xfm, this.xfm); diff --git a/cortex/webgl/resources/js/facepick.js b/cortex/webgl/resources/js/facepick.js index 6df798239..f29ab632b 100644 --- a/cortex/webgl/resources/js/facepick.js +++ b/cortex/webgl/resources/js/facepick.js @@ -143,6 +143,7 @@ FacePick.prototype = { var pos = world(xbuf, ybuf, zbuf); if (pos && this.lkdt && this.rkdt) { var left = this.lkdt.nearest([pos.x, pos.y, pos.z], 1)[0]; + console.log(left); var right = this.rkdt.nearest([pos.x, pos.y, pos.z], 1)[0]; if (left[1] < right[1]) return {hemi:"left", ptidx:left[0][3], dist:left[1], pos:pos}; @@ -156,8 +157,39 @@ FacePick.prototype = { var p = this._pick(x, y); if (p) { var vec = this.viewer.uniforms.volxfm.value[0].multiplyVector3(p.pos.clone()); - console.log("Picked vertex "+p.ptidx+" in "+p.hemi+" hemisphere, distance="+p.dist+", voxel=["+vec.x+","+vec.y+","+vec.z+"]"); + mniidx = (p.ptidx)*4 ; + if (p.hemi==="left") + hem = this.viewer.meshes.left.geometry ; + if (p.hemi==="right") + hem = this.viewer.meshes.right.geometry ; + + mnix = hem.attributes.mnicoords.array[mniidx] ; + mniy = hem.attributes.mnicoords.array[mniidx+1] ; + mniz = hem.attributes.mnicoords.array[mniidx+2] ; + + console.dir(p); + console.log("Clicked at " + x +"," + y + " and picked vertex "+p.ptidx+" in "+p.hemi+" hemisphere, distance="+p.dist+", pos=["+p.pos.x+","+p.pos.y+","+p.pos.z+"], voxel=["+vec.x+","+vec.y+","+vec.z+"]... \nMNI x index = " + mniidx + "mnix = " + mnix); + console.dir(hem); + this.addMarker(p.hemi, p.ptidx, keep); + $(this.viewer.object).find("#mnibox").show() ; + $(this.viewer.object).find("#mnibox").html("
"); + // store reference to viewer so we can send a message to it within this binding function? + /*var vwr = this ; + $(this.viewer.object).find("#mnisubmit").submit(function() { + console.dir(this) ; + coordstr = this.firstChild.value ; + console.log(this.firstChild.value) ; + console.log(coordstr) ; + coords = coordstr.split(",") ; + sub_x = coords[0] ; + sub_y = coords[1] ; + sub_z = coords[2] ; + vwr.viewer.figure.notify("pick",vwr,vwr.viewer.uniforms.volxfm.value[0].multiplyVector3(THREE.Vector3(sub_x,sub_y,sub_z))) ; + }) ;*/ + //$(this.viewer.object).find("#mnicoords").prop('disabled',false) ; + /*$(this.viewer.object).find("#mnicoords").show() ; + $(this.viewer.object).find("#mnicoords").value = ""+mnix.toFixed(2)+","+mniy.toFixed(2)+","+mniz.toFixed(2) ;*/ this.viewer.figure.notify("pick", this, [vec]); if (this.callback !== undefined) this.callback(vec, p.hemi, p.ptidx); diff --git a/cortex/webgl/resources/js/mriview.js b/cortex/webgl/resources/js/mriview.js index a12b4f5b0..0c0a45a20 100644 --- a/cortex/webgl/resources/js/mriview.js +++ b/cortex/webgl/resources/js/mriview.js @@ -203,6 +203,10 @@ var mriview = (function(module) { module.Viewer.prototype.load = function(ctminfo, callback) { var loader = new THREE.CTMLoader(false); loader.loadParts( ctminfo, function( geometries, materials, header, json ) { + console.dir(geometries) + console.dir(materials) + console.dir(header) + //console.dir(json) geometries[0].computeBoundingBox(); geometries[1].computeBoundingBox(); @@ -246,7 +250,9 @@ var mriview = (function(module) { for (var i = 0, il = geometries[right].morphTargets.length; i < il; i++) { posdata[name].push(geometries[right].morphTargets[i]); } + console.dir(geometries[right].attributes); geometries[right].reorderVertices(); + console.dir(geometries[right].attributes); geometries[right].dynamic = true; //var geom = splitverts(geometries[right], name == "right" ? leftlen : 0); @@ -303,6 +309,8 @@ var mriview = (function(module) { this.picker.undblpick(); }.bind(this)); + console.dir(this); + this.setState("target", this.surfcenter); if (typeof(callback) == "function") diff --git a/cortex/webgl/template.html b/cortex/webgl/template.html index 797f0dd11..30191dbc2 100644 --- a/cortex/webgl/template.html +++ b/cortex/webgl/template.html @@ -286,6 +286,9 @@
Loading brain...
+
+
+
From 12315334bc5a5e7f33fb72de9572f5b38d79a411 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Fri, 23 Oct 2015 16:24:18 -0700 Subject: [PATCH 11/35] modified align.py to return the points and associated mni coords. also removed extraneous argument. --- cortex/align.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/cortex/align.py b/cortex/align.py index c96922166..45e0bf12f 100644 --- a/cortex/align.py +++ b/cortex/align.py @@ -164,9 +164,6 @@ def anat_to_mni(subject, xfmname, noclean=False): Subject identifier. xfmname : str String identifying the transform to be created. - anatimg : str - Path to a nibabel-readable image that will be used as the reference for this transform. - This should be a 3D anatomical volume. noclean : bool, optional If True intermediate files will not be removed from /tmp (this is useful for debugging things), and the returned value will be the name of the temp directory. Default False. @@ -300,6 +297,8 @@ def anat_to_mni(subject, xfmname, noclean=False): mni_surfinfo_fn = db.get_paths(subject)['surfinfo'].format(type='mnicoords',opts='') np.savez(mni_surfinfo_fn,leftpts=left,rightpts=right) + return (pts, mni_coords) + finally: pass From c372264c706ff5f89dae0ac33c96a184e97d04b6 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Tue, 27 Oct 2015 12:05:32 -0700 Subject: [PATCH 12/35] completed MNI coordinate lookup feature. clicking anywhere in the brian brings up the MNI coords box for the selected vertex. Coordinates are editable, and hitting enter or Submit looks up the nearest MNI coordinate using a kdTree of MNI points which was previously created. It then adds a marker there. --- cortex/webgl/resources/css/mriview.css | 8 +++++ cortex/webgl/resources/js/ctm/CTMLoader.js | 6 ---- cortex/webgl/resources/js/facepick.js | 33 ++++++-------------- cortex/webgl/resources/js/facepick_worker.js | 7 +++-- cortex/webgl/resources/js/mriview.js | 25 ++++++++++----- cortex/webgl/template.html | 10 +++++- 6 files changed, 48 insertions(+), 41 deletions(-) diff --git a/cortex/webgl/resources/css/mriview.css b/cortex/webgl/resources/css/mriview.css index de8ad6774..7bd1256ac 100644 --- a/cortex/webgl/resources/css/mriview.css +++ b/cortex/webgl/resources/css/mriview.css @@ -154,6 +154,14 @@ a:visited { color:white; } font-weight:bold; } +#mnibox input { + border: 0; + font-size: 16pt; + font-weight: bold ; + width: 100px ; + margin-right: 5px ; +} + #mnicoords, #mnisubmit { z-index:13; } diff --git a/cortex/webgl/resources/js/ctm/CTMLoader.js b/cortex/webgl/resources/js/ctm/CTMLoader.js index 72e6c96b4..e8a82e12e 100644 --- a/cortex/webgl/resources/js/ctm/CTMLoader.js +++ b/cortex/webgl/resources/js/ctm/CTMLoader.js @@ -93,8 +93,6 @@ THREE.CTMLoader.prototype.load = function( url, callback, useWorker, useBuffers, var length = 0; - console.log(url) - xhr.onreadystatechange = function() { if ( xhr.readyState === 4 ) { @@ -202,8 +200,6 @@ THREE.CTMLoader.prototype.createModelBuffers = function ( file, callback, header var Model = function ( ) { - console.log(file); - console.log('createModelBuffers'); var scope = this; scope.materials = []; @@ -238,8 +234,6 @@ THREE.CTMLoader.prototype.createModelBuffers = function ( file, callback, header attrname = file.body.attrMaps[ i ].name; array = file.body.attrMaps[ i ].attr; - console.log(attrname) - console.log(array) if (attrname.slice(0, 11) == "morphTarget") { diff --git a/cortex/webgl/resources/js/facepick.js b/cortex/webgl/resources/js/facepick.js index f29ab632b..bb9004f2c 100644 --- a/cortex/webgl/resources/js/facepick.js +++ b/cortex/webgl/resources/js/facepick.js @@ -7,11 +7,14 @@ function FacePick(viewer, left, right) { return Math.sqrt((a[0]-b[0])*(a[0]-b[0]) + (a[1]-b[1])*(a[1]-b[1]) + (a[2]-b[2])*(a[2]-b[2])); } kdt = new kdTree([], dist, [0, 1, 2]); + mni_kdt = new kdTree([], dist, [0, 1, 2]); kdt.root = e.data.kdt; + mni_kdt.root = e.data.mnikdt; this[e.data.name] = kdt; + this['mni_' + e.data.name] = mni_kdt; }.bind(this)); - worker.postMessage({pos:left, name:"lkdt"}); - worker.postMessage({pos:right, name:"rkdt"}); + worker.postMessage({pos:left, name:"lkdt", mni:this.viewer.meshes.left.geometry.attributes.mnicoords.array}); + worker.postMessage({pos:right, name:"rkdt", mni:this.viewer.meshes.right.geometry.attributes.mnicoords.array}); this.axes = []; @@ -143,7 +146,6 @@ FacePick.prototype = { var pos = world(xbuf, ybuf, zbuf); if (pos && this.lkdt && this.rkdt) { var left = this.lkdt.nearest([pos.x, pos.y, pos.z], 1)[0]; - console.log(left); var right = this.rkdt.nearest([pos.x, pos.y, pos.z], 1)[0]; if (left[1] < right[1]) return {hemi:"left", ptidx:left[0][3], dist:left[1], pos:pos}; @@ -167,29 +169,11 @@ FacePick.prototype = { mniy = hem.attributes.mnicoords.array[mniidx+1] ; mniz = hem.attributes.mnicoords.array[mniidx+2] ; - console.dir(p); - console.log("Clicked at " + x +"," + y + " and picked vertex "+p.ptidx+" in "+p.hemi+" hemisphere, distance="+p.dist+", pos=["+p.pos.x+","+p.pos.y+","+p.pos.z+"], voxel=["+vec.x+","+vec.y+","+vec.z+"]... \nMNI x index = " + mniidx + "mnix = " + mnix); - console.dir(hem); - this.addMarker(p.hemi, p.ptidx, keep); $(this.viewer.object).find("#mnibox").show() ; - $(this.viewer.object).find("#mnibox").html("
"); - // store reference to viewer so we can send a message to it within this binding function? - /*var vwr = this ; - $(this.viewer.object).find("#mnisubmit").submit(function() { - console.dir(this) ; - coordstr = this.firstChild.value ; - console.log(this.firstChild.value) ; - console.log(coordstr) ; - coords = coordstr.split(",") ; - sub_x = coords[0] ; - sub_y = coords[1] ; - sub_z = coords[2] ; - vwr.viewer.figure.notify("pick",vwr,vwr.viewer.uniforms.volxfm.value[0].multiplyVector3(THREE.Vector3(sub_x,sub_y,sub_z))) ; - }) ;*/ - //$(this.viewer.object).find("#mnicoords").prop('disabled',false) ; - /*$(this.viewer.object).find("#mnicoords").show() ; - $(this.viewer.object).find("#mnicoords").value = ""+mnix.toFixed(2)+","+mniy.toFixed(2)+","+mniz.toFixed(2) ;*/ + $(this.viewer.object).find("#mniX").val(mnix.toFixed(2)) ; + $(this.viewer.object).find("#mniY").val(mniy.toFixed(2)) ; + $(this.viewer.object).find("#mniZ").val(mniz.toFixed(2)) ; this.viewer.figure.notify("pick", this, [vec]); if (this.callback !== undefined) this.callback(vec, p.hemi, p.ptidx); @@ -198,6 +182,7 @@ FacePick.prototype = { this.axes[i].obj.parent.remove(this.axes[i].obj); } this.axes = []; + $(this.viewer.object).find("#mnibox").hide() ; this.viewer.schedule(); } }, diff --git a/cortex/webgl/resources/js/facepick_worker.js b/cortex/webgl/resources/js/facepick_worker.js index 47c3142b9..a13bd98c4 100644 --- a/cortex/webgl/resources/js/facepick_worker.js +++ b/cortex/webgl/resources/js/facepick_worker.js @@ -3,15 +3,18 @@ importScripts( "kdTree-min.js" ); var num = 0; self.onmessage = function( event ) { var pts = []; + var mnipts = []; for (var i = 0, il = event.data.pos.length; i < il; i+= 3) pts.push([event.data.pos[i], event.data.pos[i+1], event.data.pos[i+2], i/3]); - + for (var i = 0, il = event.data.mni.length; i < il; i+= 4) + mnipts.push([event.data.mni[i], event.data.mni[i+1], event.data.mni[i+2], i/4]) var dist = function (a, b) { return Math.sqrt((a[0]-b[0])*(a[0]-b[0]) + (a[1]-b[1])*(a[1]-b[1]) + (a[2]-b[2])*(a[2]-b[2])); } var kdt = new kdTree(pts, dist, [0, 1, 2]); + var mnikdt = new kdTree(mnipts, dist, [0, 1, 2]); - self.postMessage( {kdt:kdt.root, name:event.data.name} ); + self.postMessage( {kdt:kdt.root, name:event.data.name, mnikdt:mnikdt.root} ); if (++num > 1) self.close(); } diff --git a/cortex/webgl/resources/js/mriview.js b/cortex/webgl/resources/js/mriview.js index 0c0a45a20..2d64e44ff 100644 --- a/cortex/webgl/resources/js/mriview.js +++ b/cortex/webgl/resources/js/mriview.js @@ -203,10 +203,6 @@ var mriview = (function(module) { module.Viewer.prototype.load = function(ctminfo, callback) { var loader = new THREE.CTMLoader(false); loader.loadParts( ctminfo, function( geometries, materials, header, json ) { - console.dir(geometries) - console.dir(materials) - console.dir(header) - //console.dir(json) geometries[0].computeBoundingBox(); geometries[1].computeBoundingBox(); @@ -250,9 +246,7 @@ var mriview = (function(module) { for (var i = 0, il = geometries[right].morphTargets.length; i < il; i++) { posdata[name].push(geometries[right].morphTargets[i]); } - console.dir(geometries[right].attributes); geometries[right].reorderVertices(); - console.dir(geometries[right].attributes); geometries[right].dynamic = true; //var geom = splitverts(geometries[right], name == "right" ? leftlen : 0); @@ -309,8 +303,6 @@ var mriview = (function(module) { this.picker.undblpick(); }.bind(this)); - console.dir(this); - this.setState("target", this.surfcenter); if (typeof(callback) == "function") @@ -1111,6 +1103,23 @@ var mriview = (function(module) { }.bind(this)); } + $(this.object).find("#mniform").submit(function() { + x = $(this.object).find("#mniX").val(); + y = $(this.object).find("#mniY").val(); + z = $(this.object).find("#mniZ").val(); + + var left = this.picker.mni_lkdt.nearest([x, y, z], 1)[0]; + var right = this.picker.mni_rkdt.nearest([x, y, z], 1)[0]; + if (left[1] < right[1]) { + this.picker.addMarker("left", left[0][3], false); + } + else { + this.picker.addMarker("right", right[0][3], false); + } + return(0); //do not reload page + }.bind(this)); + + // Setup controls for multiple overlays var updateOverlays = function() { this.roipack.update(this.renderer).done(function(tex){ diff --git a/cortex/webgl/template.html b/cortex/webgl/template.html index 30191dbc2..36eb15c41 100644 --- a/cortex/webgl/template.html +++ b/cortex/webgl/template.html @@ -287,7 +287,15 @@
Loading brain...
-
+
+ MNI coordinates +
+ X: mm
+ Y: mm
+ Z: mm
+ +
+
From e1f323899a097ed6891212235c67e8d75e0b5c1a Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Tue, 27 Oct 2015 12:09:12 -0700 Subject: [PATCH 13/35] change text of Submit button to Go --- cortex/webgl/template.html | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cortex/webgl/template.html b/cortex/webgl/template.html index 36eb15c41..79516eb31 100644 --- a/cortex/webgl/template.html +++ b/cortex/webgl/template.html @@ -293,7 +293,7 @@ X: mm
Y: mm
Z: mm
- + From 944a4ebf6e17adc663c8bc7fc421f2946dac11be Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Tue, 27 Oct 2015 14:34:29 -0700 Subject: [PATCH 14/35] correct docstring, remove debug output, enable actual calling of commands. --- cortex/align.py | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/cortex/align.py b/cortex/align.py index 45e0bf12f..799e5b648 100644 --- a/cortex/align.py +++ b/cortex/align.py @@ -152,7 +152,7 @@ def automatic(subject, xfmname, reference, noclean=False, bbrtype="signed"): return retval -def anat_to_mni(subject, xfmname, noclean=False): +def anat_to_mni(subject, xfmname): """Create an automatic alignment of an anatomical image to the MNI standard. If `noclean`, intermediate files will not be removed from /tmp. The `reference` image and resulting @@ -170,7 +170,8 @@ def anat_to_mni(subject, xfmname, noclean=False): Returns ------- - Nothing unless `noclean` is True. + pts : the vertices of the fiducial surface + mnipts : the mni coordinates of those vertices (same shape as pts) """ import shlex @@ -200,16 +201,16 @@ def anat_to_mni(subject, xfmname, noclean=False): reorient_anat = 'reorient_anat' reorient_cmd = '{fslpre}fslreorient2std {raw_anat} {adir}/{ra_raw}'.format(fslpre=fsl_prefix,raw_anat=raw_anat, adir=odir, ra_raw=reorient_anat) print('Reorienting anatomicals using fslreorient2std, cmd like: \n%s' % reorient_cmd) - #if sp.call(reorient_cmd, shell=True) != 0: - # raise IOError('Error calling fslreorient2std on raw anatomical') + if sp.call(reorient_cmd, shell=True) != 0: + raise IOError('Error calling fslreorient2std on raw anatomical') reorient_cmd = '{fslpre}fslreorient2std {bet_anat} {adir}/{ra_raw}_brain'.format(fslpre=fsl_prefix,bet_anat=bet_anat, adir=odir, ra_raw=reorient_anat) - #if sp.call(reorient_cmd, shell=True) != 0: - # raise IOError('Error calling fslreorient2std on brain-extracted anatomical') + if sp.call(reorient_cmd, shell=True) != 0: + raise IOError('Error calling fslreorient2std on brain-extracted anatomical') ra_betmask = reorient_anat + "_brainmask" reorient_cmd = '{fslpre}fslreorient2std {bet_anat} {adir}/{ra_betmask}'.format(fslpre=fsl_prefix,bet_anat=betmask_anat, adir=odir, ra_betmask=ra_betmask) - #if sp.call(reorient_cmd, shell=True) != 0: - # raise IOError('Error calling fslreorient2std on brain-extracted mask') + if sp.call(reorient_cmd, shell=True) != 0: + raise IOError('Error calling fslreorient2std on brain-extracted mask') fsldir = os.environ['FSLDIR'] standard = '%s/data/standard/MNI152_T1_1mm'%fsldir @@ -221,18 +222,17 @@ def anat_to_mni(subject, xfmname, noclean=False): flirt_cmd = '{fslpre}flirt -in {bet_standard} -ref {adir}/{ra_raw}_brain -dof 6 -omat /tmp/{cout}_flirt' flirt_cmd = flirt_cmd.format(fslpre=fsl_prefix, ra_raw=reorient_anat, bet_standard=bet_standard, adir=odir, cout=cout) print('Running FLIRT to estimate initial affine transform with command:\n%s'%flirt_cmd) - #if sp.call(flirt_cmd, shell=True) != 0: - # raise IOError('Error calling FLIRT with command: %s' % flirt_cmd) + if sp.call(flirt_cmd, shell=True) != 0: + raise IOError('Error calling FLIRT with command: %s' % flirt_cmd) # FNIRT mni-to-anat transform estimation cmd (does not apply any transform, but generates estimate [cout]) cmd = '{fslpre}fnirt --in={standard} --ref={ad}/{ra_raw} --refmask={ad}/{refmask} --aff=/tmp/{cout}_flirt --cout={anat_dir}/{cout}_fnirt --fout={anat_dir}/{cout}_field --iout=/tmp/mni2anat_iout --config=T1_2_MNI152_2mm' cmd = cmd.format(fslpre=fsl_prefix, ra_raw=reorient_anat, standard=standard, refmask=ra_betmask, ad=odir, anat_dir=anat_dir, cout=cout) print('Running FNIRT to estimate transform, using the following command... this can take a while:\n%s'%cmd) - #if sp.call(cmd, shell=True) != 0: - # raise IOError('Error calling fnirt with cmd: %s'%cmd) + if sp.call(cmd, shell=True) != 0: + raise IOError('Error calling fnirt with cmd: %s'%cmd) [pts, polys] = db.get_surf(subject,"fiducial",merge="True") - # np.savetxt(cfile, pts, fmt='%g') #print('raw anatomical: %s\nbet anatomical: %s\nflirt cmd:%s\nfnirt cmd: %s\npts: %s' % (raw_anat,bet_anat,flirt_cmd,cmd,pts)) @@ -283,16 +283,16 @@ def anat_to_mni(subject, xfmname, noclean=False): # some debug output # print pts, mni_coords - print pts[0], mni_coords[0] - print len(pts), len(mni_coords) - print type(pts), type(pts[0][0]), type(mni_coords) + # print pts[0], mni_coords[0] + # print len(pts), len(mni_coords) + # print type(pts), type(pts[0][0]), type(mni_coords) # now split mni_coords into left and right arrays for saving nverts_L = len(warpverts_L) - print nverts_L + #print nverts_L left = mni_coords[:nverts_L] right = mni_coords[nverts_L:] - print len(left), len(right) + #print len(left), len(right) mni_surfinfo_fn = db.get_paths(subject)['surfinfo'].format(type='mnicoords',opts='') np.savez(mni_surfinfo_fn,leftpts=left,rightpts=right) From fb2aa8dd072303e7f22847f3f943d57ec7962375 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Mon, 2 Nov 2015 16:46:19 -0800 Subject: [PATCH 15/35] removed some debug output --- cortex/brainctm.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/cortex/brainctm.py b/cortex/brainctm.py index 29ac9bbcc..39c07a9e8 100644 --- a/cortex/brainctm.py +++ b/cortex/brainctm.py @@ -112,8 +112,6 @@ def save(self, path, method='mg2', disp_layers=['rois'], extra_disp=None, **kwar svgname = path+".svg" jsname = path+".json" - print("In method save() of BrainCTM, ctmname: %s"%ctmname) - # Save CTM concatenation (lpts, _, _), lbin = self.left.save(method=method, **kwargs) (rpts, _, _), rbin = self.right.save(method=method, **kwargs) @@ -190,7 +188,6 @@ def __init__(self, pts, polys, norms=None): def addSurf(self, pts, name=None, renorm=True): '''Scales the in-between surfaces to be same scale as fiducial''' - print('In Hemi.addSurf()...') if name is None: name = 'morphTarget%d'%len(self.surfs) @@ -202,7 +199,6 @@ def addSurf(self, pts, name=None, renorm=True): attrib = np.hstack([rnorm, np.zeros((len(rnorm),1))]) self.surfs[name] = attrib - print('Adding attrib %s'%name) self.ctm.addAttrib(attrib, name) def setFlat(self, pts): @@ -211,7 +207,6 @@ def setFlat(self, pts): def save(self, **kwargs): self.ctm.addAttrib(self.aux, 'auxdat') - print('Adding attrib mnicoords, just added auxdat...') self.ctm.addAttrib(self.mni, 'mnicoords') self.ctm.save(**kwargs) From acb4d62e2210c14fa5020241d5763a3b37afb207 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Mon, 2 Nov 2015 16:47:15 -0800 Subject: [PATCH 16/35] removed some debug output --- cortex/brainctm.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/cortex/brainctm.py b/cortex/brainctm.py index 39c07a9e8..cc266958a 100644 --- a/cortex/brainctm.py +++ b/cortex/brainctm.py @@ -252,7 +252,6 @@ def make_pack(outfile, subj, types=("inflated",), method='raw', level=0, decimate=False, disp_layers=['rois'],extra_disp=None): """Generates a cached CTM file""" - print 'In make_pack, outfile %s'%outfile ctm = BrainCTM(subj, decimate=decimate) ctm.addCurvature() ctm.addMNI() @@ -269,7 +268,6 @@ def make_pack(outfile, subj, types=("inflated",), method='raw', level=0, extra_disp=extra_disp) def read_pack(ctmfile): - print('In read_pack, ctmfile: %s'%ctmfile) fname = os.path.splitext(ctmfile)[0] jsfile = json.load(open(fname+".json")) offset = jsfile['offsets'] From 2e2490ebceec9c79713925c3ca00211935abdc54 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Mon, 2 Nov 2015 16:49:32 -0800 Subject: [PATCH 17/35] removed superfluous reference to warpmapper module, which no longer exists. --- cortex/mapper/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cortex/mapper/__init__.py b/cortex/mapper/__init__.py index c781aa53e..da88908aa 100644 --- a/cortex/mapper/__init__.py +++ b/cortex/mapper/__init__.py @@ -10,7 +10,7 @@ def get_mapper(subject, xfmname, type='nearest', recache=False, **kwargs): from ..database import db - from . import point, patch, volume, line, warp + from . import point, patch, volume, line mapcls = dict( nearest=point.PointNN, From d57d6fd052b224d33fcd55ef44b6acab5d180f05 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Mon, 2 Nov 2015 16:52:06 -0800 Subject: [PATCH 18/35] added stub documentation about anat_to_mni() --- docs/align.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/align.rst b/docs/align.rst index a54c0d0dd..39776297b 100644 --- a/docs/align.rst +++ b/docs/align.rst @@ -3,3 +3,5 @@ Aligner This module contains functions for manual and automatic alignment of brain images and cortical surfaces. For each transform, it saves a transform matrix, which maps pixels to voxels. The automatic() function does epi-to-anat registration using FLIRT with BBR, then inverts this with Transform.from_fsl() + +The anat_to_mni() function uses FNIRT to register the MNI surface to the subject's surface, then estimates the MNI coords corresponding to vertex coords. \ No newline at end of file From fc055ff0d72a588c9297f12fd26ed0b769185c70 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Mon, 9 Nov 2015 11:24:36 -0800 Subject: [PATCH 19/35] styling changes to make the coordinates box less space-filling. --- cortex/webgl/resources/css/mriview.css | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/cortex/webgl/resources/css/mriview.css b/cortex/webgl/resources/css/mriview.css index 7bd1256ac..4d366d567 100644 --- a/cortex/webgl/resources/css/mriview.css +++ b/cortex/webgl/resources/css/mriview.css @@ -145,20 +145,20 @@ a:visited { color:white; } float:left; text-shadow:0px 2px 8px black, 0px 1px 8px black; color:white; - font-size:24pt; + font-size:16pt; padding:10px; margin:20px; - margin-top: 90px ; + margin-top: 80px ; border-radius:10px; background:rgba(255,255,255,.4); font-weight:bold; } #mnibox input { - border: 0; - font-size: 16pt; + /*border: 0;*/ + font-size: 12pt; font-weight: bold ; - width: 100px ; + width: 70px ; margin-right: 5px ; } From f7cba2b8447a7ed30754669f6c9451dc39b563fa Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Mon, 9 Nov 2015 23:57:04 -0800 Subject: [PATCH 20/35] some small but important changes: all files except the final surfinfo are now stored in a per-session tempfile, so they don't permanently clutter the db but are still safe for multiple sessions. some unneeded parameters and imports were removed. additional comments were added, and the do parameter, which allows the function to be run in dummy mode if False. --- cortex/align.py | 57 ++++++++++++++++++++++++++++++------------------- 1 file changed, 35 insertions(+), 22 deletions(-) diff --git a/cortex/align.py b/cortex/align.py index 799e5b648..0504c9683 100644 --- a/cortex/align.py +++ b/cortex/align.py @@ -152,26 +152,27 @@ def automatic(subject, xfmname, reference, noclean=False, bbrtype="signed"): return retval -def anat_to_mni(subject, xfmname): +def anat_to_mni(subject, do=True): """Create an automatic alignment of an anatomical image to the MNI standard. - If `noclean`, intermediate files will not be removed from /tmp. The `reference` image and resulting - transform called `xfmname` will be automatically stored in the database. + This function does the following: + 1) Re-orders orientation labels on anatomical images using fslreorient2std (without modifying the existing files) + 2) Calls FLIRT + 3) Calls FNIRT with the transform estimated by FLIRT as the specified + 4) Gets the resulting warp field, samples it at each vertex location, and calculates MNI coordinates. + 5) Saves these coordinates as a surfinfo file in the db. Parameters ---------- subject : str Subject identifier. - xfmname : str - String identifying the transform to be created. - noclean : bool, optional - If True intermediate files will not be removed from /tmp (this is useful for debugging things), - and the returned value will be the name of the temp directory. Default False. + do : bool + Actually execute the commands (True), or just print them (False, useful for debugging). Returns ------- pts : the vertices of the fiducial surface - mnipts : the mni coordinates of those vertices (same shape as pts) + mnipts : the mni coordinates of those vertices (same shape as pts, corresponding indices) """ import shlex @@ -186,30 +187,30 @@ def anat_to_mni(subject, xfmname): from .dataset import Volume fsl_prefix = config.get("basic", "fsl_prefix") - schfile = os.path.join(os.path.split(os.path.abspath(__file__))[0], "bbr.sch") + cache = tempfile.mkdtemp() - print('anat_to_mni, subject: %s, xfmname: %s' % (subject, xfmname)) + print('anat_to_mni, subject: %s' % subject) try: raw_anat = db.get_anat(subject, type='raw').get_filename() bet_anat = db.get_anat(subject, type='brainmask').get_filename() betmask_anat = db.get_anat(subject, type='brainmask_mask').get_filename() anat_dir = os.path.dirname(raw_anat) - odir = anat_dir + odir = cache # stem for the reoriented-into-MNI anatomical images (required by FLIRT/FNIRT) reorient_anat = 'reorient_anat' reorient_cmd = '{fslpre}fslreorient2std {raw_anat} {adir}/{ra_raw}'.format(fslpre=fsl_prefix,raw_anat=raw_anat, adir=odir, ra_raw=reorient_anat) print('Reorienting anatomicals using fslreorient2std, cmd like: \n%s' % reorient_cmd) - if sp.call(reorient_cmd, shell=True) != 0: + if do and sp.call(reorient_cmd, shell=True) != 0: raise IOError('Error calling fslreorient2std on raw anatomical') reorient_cmd = '{fslpre}fslreorient2std {bet_anat} {adir}/{ra_raw}_brain'.format(fslpre=fsl_prefix,bet_anat=bet_anat, adir=odir, ra_raw=reorient_anat) - if sp.call(reorient_cmd, shell=True) != 0: + if do and sp.call(reorient_cmd, shell=True) != 0: raise IOError('Error calling fslreorient2std on brain-extracted anatomical') ra_betmask = reorient_anat + "_brainmask" reorient_cmd = '{fslpre}fslreorient2std {bet_anat} {adir}/{ra_betmask}'.format(fslpre=fsl_prefix,bet_anat=betmask_anat, adir=odir, ra_betmask=ra_betmask) - if sp.call(reorient_cmd, shell=True) != 0: + if do and sp.call(reorient_cmd, shell=True) != 0: raise IOError('Error calling fslreorient2std on brain-extracted mask') fsldir = os.environ['FSLDIR'] @@ -219,26 +220,33 @@ def anat_to_mni(subject, xfmname): cout = 'mni2anat' #stem of the filenames of the transform estimates # initial affine anatomical-to-standard registration using FLIRT. required, as the output xfm is used as a start by FNIRT. - flirt_cmd = '{fslpre}flirt -in {bet_standard} -ref {adir}/{ra_raw}_brain -dof 6 -omat /tmp/{cout}_flirt' + flirt_cmd = '{fslpre}flirt -in {bet_standard} -ref {adir}/{ra_raw}_brain -dof 6 -omat {adir}/{cout}_flirt' flirt_cmd = flirt_cmd.format(fslpre=fsl_prefix, ra_raw=reorient_anat, bet_standard=bet_standard, adir=odir, cout=cout) print('Running FLIRT to estimate initial affine transform with command:\n%s'%flirt_cmd) - if sp.call(flirt_cmd, shell=True) != 0: + if do and sp.call(flirt_cmd, shell=True) != 0: raise IOError('Error calling FLIRT with command: %s' % flirt_cmd) # FNIRT mni-to-anat transform estimation cmd (does not apply any transform, but generates estimate [cout]) - cmd = '{fslpre}fnirt --in={standard} --ref={ad}/{ra_raw} --refmask={ad}/{refmask} --aff=/tmp/{cout}_flirt --cout={anat_dir}/{cout}_fnirt --fout={anat_dir}/{cout}_field --iout=/tmp/mni2anat_iout --config=T1_2_MNI152_2mm' + # the MNI152 2mm config is used even though we're referencing 1mm, per this FSL list post: + # https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=FSL;d14e5a9d.1105 + cmd = '{fslpre}fnirt --in={standard} --ref={ad}/{ra_raw} --refmask={ad}/{refmask} --aff={ad}/{cout}_flirt --cout={ad}/{cout}_fnirt --fout={ad}/{cout}_field --iout={ad}/{cout}_iout --config=T1_2_MNI152_2mm' cmd = cmd.format(fslpre=fsl_prefix, ra_raw=reorient_anat, standard=standard, refmask=ra_betmask, ad=odir, anat_dir=anat_dir, cout=cout) print('Running FNIRT to estimate transform, using the following command... this can take a while:\n%s'%cmd) - if sp.call(cmd, shell=True) != 0: + if do and sp.call(cmd, shell=True) != 0: raise IOError('Error calling fnirt with cmd: %s'%cmd) [pts, polys] = db.get_surf(subject,"fiducial",merge="True") - #print('raw anatomical: %s\nbet anatomical: %s\nflirt cmd:%s\nfnirt cmd: %s\npts: %s' % (raw_anat,bet_anat,flirt_cmd,cmd,pts)) + print('raw anatomical: %s\nbet anatomical: %s\nflirt cmd:%s\nfnirt cmd: %s\npts: %s' % (raw_anat,bet_anat,flirt_cmd,cmd,pts)) # take the reoriented anatomical, get its affine coord transform, invert this, and save it reo_xfmnm = 'reorient_inv' - re_anat = db.get_anat(subject,reorient_anat) + # need to change this line, as the reoriented anatomical is not in the db but in /tmp now + # re_anat = db.get_anat(subject,reorient_anat) + reo_anat_fn = '{odir}/{reorient_anat}.nii.gz'.format(odir=odir,reorient_anat=reorient_anat) + print(reo_anat_fn) + # since the reoriented anatomicals aren't stored in the db anymore, db.get_anat() will not work (?) + re_anat = nib.load(reo_anat_fn) reo_xfm = Transform(np.linalg.inv(re_anat.get_affine()),re_anat) reo_xfm.save(subject,reo_xfmnm,"coord") @@ -247,7 +255,12 @@ def anat_to_mni(subject, xfmname): aqfinv = np.linalg.inv(aqf) # load the warp field data as a volume - warp = db.get_anat(subject,'%s_field'%cout) + # since it's not in the db anymore but in /tmp instead of: + # warp = db.get_anat(subject,'%s_field'%cout) + # it's this: + warp_fn = '{ad}/{cout}_field.nii.gz'.format(ad=odir,cout=cout) + print warp_fn + warp = nib.load(warp_fn) wd = warp.get_data() # need in (t,z,y,x) order wd = np.swapaxes(wd,0,3) # x <--> t From 5f2927beaf1228e5fdc6bde0f77fc502a64797a4 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Tue, 10 Nov 2015 08:51:48 -0800 Subject: [PATCH 21/35] changed interpolation from nearest-neighbor to lanczos, which is ~sinc --- cortex/align.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cortex/align.py b/cortex/align.py index 0504c9683..279331c79 100644 --- a/cortex/align.py +++ b/cortex/align.py @@ -269,7 +269,7 @@ def anat_to_mni(subject, do=True): # now do the mapping! this gets the warp field values at the corresponding points # (uses fiducial surface by default) - warpvd = wv.map(projection="nearest") + warpvd = wv.map(projection="lanczos") # reshape into something sensible warpverts_L = [vs for vs in np.swapaxes(warpvd.left,0,1)] From cc009a490f3d808f946bf84ecb611bbaf06178dd Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Tue, 10 Nov 2015 09:56:56 -0800 Subject: [PATCH 22/35] added a check to ensure that this still works if there are no MNI coords... --- cortex/brainctm.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/cortex/brainctm.py b/cortex/brainctm.py index cc266958a..7cfffa073 100644 --- a/cortex/brainctm.py +++ b/cortex/brainctm.py @@ -93,8 +93,13 @@ def addCurvature(self, **kwargs): def addMNI(self, **kwargs): print('Adding MNI coords...') npz = db.get_surfinfo(self.subject, type='mnicoords', **kwargs) - self.left.mni[:,:-1] = npz['leftpts'] - self.right.mni[:,:-1] = npz['rightpts'] + try: + self.left.mni[:,:-1] = npz['leftpts'] + self.right.mni[:,:-1] = npz['rightpts'] + except AttributeError: + self.left.mni = [] + self.right.mni = [] + def save(self, path, method='mg2', disp_layers=['rois'], extra_disp=None, **kwargs): """Save CTM file for static html display. From 993f493e5f2654a179eba575c819aa3655b4d456 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Tue, 10 Nov 2015 09:58:09 -0800 Subject: [PATCH 23/35] added coord_disp=True parameter to show(). if False, the MNI coordiante box is left out of the template. I verified that the picker still works in this case. --- cortex/webgl/template.html | 2 ++ cortex/webgl/view.py | 3 ++- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/cortex/webgl/template.html b/cortex/webgl/template.html index 79516eb31..c73430363 100644 --- a/cortex/webgl/template.html +++ b/cortex/webgl/template.html @@ -286,6 +286,7 @@
Loading brain...
+ {% if coord_disp %}
MNI coordinates @@ -297,6 +298,7 @@
+ {% end %}
diff --git a/cortex/webgl/view.py b/cortex/webgl/view.py index 83db638b9..7899a2e90 100644 --- a/cortex/webgl/view.py +++ b/cortex/webgl/view.py @@ -204,7 +204,7 @@ def make_static(outpath, data, types=("inflated",), recache=False, cmap="RdBu_r" def show(data, types=("inflated",), recache=False, cmap='RdBu_r', layout=None, autoclose=True, open_browser=True, port=None, pickerfun=None, - disp_layers=['rois'], extra_disp=None, **kwargs): + disp_layers=['rois'], extra_disp=None, coord_disp=True, **kwargs): """Display a dynamic viewer using the given dataset. See cortex.webgl.make_static for help. """ data = dataset.normalize(data) @@ -314,6 +314,7 @@ def get(self): subjects=subjectjs, disp_layers=disp_layers+dl, disp_defaults=_make_disp_defaults(disp_layers+dl), + coord_disp=coord_disp, **viewopts) self.write(generated) From ab216b6b0b9e68bd06093246e4053570e18870b1 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Tue, 10 Nov 2015 16:59:56 -0800 Subject: [PATCH 24/35] UI updates: reduced size of coord box, added radio button to select MNI or magnet isocenter coords, hooked this button up to a function that toggles the coordinate display, made the picker smarter also so that it displays the selected coordinates and not MNI every time. --- cortex/webgl/resources/css/mriview.css | 6 ++--- cortex/webgl/resources/js/facepick.js | 21 +++++++++++++---- cortex/webgl/resources/js/mriview.js | 31 ++++++++++++++++++++++++++ cortex/webgl/template.html | 7 +++++- 4 files changed, 57 insertions(+), 8 deletions(-) diff --git a/cortex/webgl/resources/css/mriview.css b/cortex/webgl/resources/css/mriview.css index 4d366d567..90ee6e539 100644 --- a/cortex/webgl/resources/css/mriview.css +++ b/cortex/webgl/resources/css/mriview.css @@ -145,7 +145,7 @@ a:visited { color:white; } float:left; text-shadow:0px 2px 8px black, 0px 1px 8px black; color:white; - font-size:16pt; + font-size:12pt; padding:10px; margin:20px; margin-top: 80px ; @@ -154,9 +154,9 @@ a:visited { color:white; } font-weight:bold; } -#mnibox input { +#mnibox input.mniinput { /*border: 0;*/ - font-size: 12pt; + font-size: 10pt; font-weight: bold ; width: 70px ; margin-right: 5px ; diff --git a/cortex/webgl/resources/js/facepick.js b/cortex/webgl/resources/js/facepick.js index bb9004f2c..2ccb68479 100644 --- a/cortex/webgl/resources/js/facepick.js +++ b/cortex/webgl/resources/js/facepick.js @@ -159,21 +159,34 @@ FacePick.prototype = { var p = this._pick(x, y); if (p) { var vec = this.viewer.uniforms.volxfm.value[0].multiplyVector3(p.pos.clone()); - mniidx = (p.ptidx)*4 ; if (p.hemi==="left") hem = this.viewer.meshes.left.geometry ; if (p.hemi==="right") hem = this.viewer.meshes.right.geometry ; - mnix = hem.attributes.mnicoords.array[mniidx] ; - mniy = hem.attributes.mnicoords.array[mniidx+1] ; - mniz = hem.attributes.mnicoords.array[mniidx+2] ; + space = $(this.viewer.object).find(".radio:checked").val(); + if (space==="magnet") { + coordarray = hem.attributes.position ; + $(this.viewer.object).find("#coordsys_mag").prop('checked',true) ; + } + else { //mni or undefined + coordarray = hem.attributes.mnicoords ; + $(this.viewer.object).find("#coordsys_mni").prop('checked',true) ; + } + + mniidx = (p.ptidx)*coordarray.itemSize ; + + mnix = coordarray.array[mniidx] ; + mniy = coordarray.array[mniidx+1] ; + mniz = coordarray.array[mniidx+2] ; this.addMarker(p.hemi, p.ptidx, keep); $(this.viewer.object).find("#mnibox").show() ; $(this.viewer.object).find("#mniX").val(mnix.toFixed(2)) ; $(this.viewer.object).find("#mniY").val(mniy.toFixed(2)) ; $(this.viewer.object).find("#mniZ").val(mniz.toFixed(2)) ; + $(this.viewer.object).find("#ptidx").val(p.ptidx) ; + $(this.viewer.object).find("#pthem").val(p.hemi) ; this.viewer.figure.notify("pick", this, [vec]); if (this.callback !== undefined) this.callback(vec, p.hemi, p.ptidx); diff --git a/cortex/webgl/resources/js/mriview.js b/cortex/webgl/resources/js/mriview.js index 2d64e44ff..a7dd0c86e 100644 --- a/cortex/webgl/resources/js/mriview.js +++ b/cortex/webgl/resources/js/mriview.js @@ -1103,6 +1103,37 @@ var mriview = (function(module) { }.bind(this)); } + $(this.object).find(".radio").change(function() { //the .radio class is specifically for coord space selection, not all radio buttons + space = $(this.object).find(".radio:checked").val(); + ptidx = $(this.object).find("#ptidx").val(); + pthem = $(this.object).find("#pthem").val(); + if (pthem==="left") + h = this.meshes.left.geometry ; + else if (pthem==="right") + h = this.meshes.right.geometry ; + + if (typeof h !== 'undefined') { + if (space==="magnet") { + coordarray = h.attributes.position ; + $(this.object).find("#coordsys_mag").prop('checked',true) ; + } + else { //mni or undefined + coordarray = h.attributes.mnicoords ; + $(this.object).find("#coordsys_mni").prop('checked',true) ; + } + + mniidx = (ptidx)*coordarray.itemSize ; + + mnix = coordarray.array[mniidx] ; + mniy = coordarray.array[mniidx+1] ; + mniz = coordarray.array[mniidx+2] ; + + $(this.object).find("#mniX").val(mnix.toFixed(2)) ; + $(this.object).find("#mniY").val(mniy.toFixed(2)) ; + $(this.object).find("#mniZ").val(mniz.toFixed(2)) ; + } + }.bind(this)); + $(this.object).find("#mniform").submit(function() { x = $(this.object).find("#mniX").val(); y = $(this.object).find("#mniY").val(); diff --git a/cortex/webgl/template.html b/cortex/webgl/template.html index c73430363..a941c4a0b 100644 --- a/cortex/webgl/template.html +++ b/cortex/webgl/template.html @@ -289,8 +289,13 @@ {% if coord_disp %}
- MNI coordinates + Coordinates
+ + MNI + Magnet
 + + X: mm
Y: mm
Z: mm
From da69a597dab7932af6fb24b6f64cedf4c08fe86f Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Tue, 10 Nov 2015 17:54:55 -0800 Subject: [PATCH 25/35] bugfix: form submission did not respect coord system radio button, always searched for nearest mni coordinate. that has been fixed now, meaning one can search by MNI and magnet iso-center coords. --- cortex/webgl/resources/js/mriview.js | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/cortex/webgl/resources/js/mriview.js b/cortex/webgl/resources/js/mriview.js index a7dd0c86e..f5626d748 100644 --- a/cortex/webgl/resources/js/mriview.js +++ b/cortex/webgl/resources/js/mriview.js @@ -1138,9 +1138,16 @@ var mriview = (function(module) { x = $(this.object).find("#mniX").val(); y = $(this.object).find("#mniY").val(); z = $(this.object).find("#mniZ").val(); + space = $(this.object).find(".radio:checked").val(); + if (space==="magnet") { + var left = this.picker.lkdt.nearest([x, y, z], 1)[0]; + var right = this.picker.rkdt.nearest([x, y, z], 1)[0]; + } + else { //mni or undefined + var left = this.picker.mni_lkdt.nearest([x, y, z], 1)[0]; + var right = this.picker.mni_rkdt.nearest([x, y, z], 1)[0]; + } - var left = this.picker.mni_lkdt.nearest([x, y, z], 1)[0]; - var right = this.picker.mni_rkdt.nearest([x, y, z], 1)[0]; if (left[1] < right[1]) { this.picker.addMarker("left", left[0][3], false); } From 783e214e556ab7780e8181ef00811f92f97ec3f1 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Wed, 11 Nov 2015 14:16:11 -0800 Subject: [PATCH 26/35] removed some extraneous output, and added one more useful line. --- cortex/align.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/cortex/align.py b/cortex/align.py index 279331c79..9627a2f5e 100644 --- a/cortex/align.py +++ b/cortex/align.py @@ -237,14 +237,14 @@ def anat_to_mni(subject, do=True): [pts, polys] = db.get_surf(subject,"fiducial",merge="True") - print('raw anatomical: %s\nbet anatomical: %s\nflirt cmd:%s\nfnirt cmd: %s\npts: %s' % (raw_anat,bet_anat,flirt_cmd,cmd,pts)) + #print('raw anatomical: %s\nbet anatomical: %s\nflirt cmd:%s\nfnirt cmd: %s\npts: %s' % (raw_anat,bet_anat,flirt_cmd,cmd,pts)) # take the reoriented anatomical, get its affine coord transform, invert this, and save it reo_xfmnm = 'reorient_inv' # need to change this line, as the reoriented anatomical is not in the db but in /tmp now # re_anat = db.get_anat(subject,reorient_anat) reo_anat_fn = '{odir}/{reorient_anat}.nii.gz'.format(odir=odir,reorient_anat=reorient_anat) - print(reo_anat_fn) + # print(reo_anat_fn) # since the reoriented anatomicals aren't stored in the db anymore, db.get_anat() will not work (?) re_anat = nib.load(reo_anat_fn) reo_xfm = Transform(np.linalg.inv(re_anat.get_affine()),re_anat) @@ -259,7 +259,7 @@ def anat_to_mni(subject, do=True): # warp = db.get_anat(subject,'%s_field'%cout) # it's this: warp_fn = '{ad}/{cout}_field.nii.gz'.format(ad=odir,cout=cout) - print warp_fn + # print warp_fn warp = nib.load(warp_fn) wd = warp.get_data() # need in (t,z,y,x) order @@ -308,6 +308,7 @@ def anat_to_mni(subject, do=True): #print len(left), len(right) mni_surfinfo_fn = db.get_paths(subject)['surfinfo'].format(type='mnicoords',opts='') + print('Saving mni coordinates as a surfinfo...') np.savez(mni_surfinfo_fn,leftpts=left,rightpts=right) return (pts, mni_coords) From 3479c53a94746ee6749608bc556cef9779086567 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Fri, 13 Nov 2015 15:24:44 -0800 Subject: [PATCH 27/35] moved align.py/anat_to_mni() to surfinfo.py/mni_nl() and removed its return values --- cortex/align.py | 165 --------------------------------------------- cortex/surfinfo.py | 150 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 150 insertions(+), 165 deletions(-) diff --git a/cortex/align.py b/cortex/align.py index 9627a2f5e..2a49a2435 100644 --- a/cortex/align.py +++ b/cortex/align.py @@ -152,171 +152,6 @@ def automatic(subject, xfmname, reference, noclean=False, bbrtype="signed"): return retval -def anat_to_mni(subject, do=True): - """Create an automatic alignment of an anatomical image to the MNI standard. - - This function does the following: - 1) Re-orders orientation labels on anatomical images using fslreorient2std (without modifying the existing files) - 2) Calls FLIRT - 3) Calls FNIRT with the transform estimated by FLIRT as the specified - 4) Gets the resulting warp field, samples it at each vertex location, and calculates MNI coordinates. - 5) Saves these coordinates as a surfinfo file in the db. - - Parameters - ---------- - subject : str - Subject identifier. - do : bool - Actually execute the commands (True), or just print them (False, useful for debugging). - - Returns - ------- - pts : the vertices of the fiducial surface - mnipts : the mni coordinates of those vertices (same shape as pts, corresponding indices) - """ - - import shlex - import shutil - import tempfile - import subprocess as sp - import nibabel as nib - - from .database import db - from .xfm import Transform - from .options import config - from .dataset import Volume - - fsl_prefix = config.get("basic", "fsl_prefix") - cache = tempfile.mkdtemp() - - print('anat_to_mni, subject: %s' % subject) - - try: - raw_anat = db.get_anat(subject, type='raw').get_filename() - bet_anat = db.get_anat(subject, type='brainmask').get_filename() - betmask_anat = db.get_anat(subject, type='brainmask_mask').get_filename() - anat_dir = os.path.dirname(raw_anat) - odir = cache - - # stem for the reoriented-into-MNI anatomical images (required by FLIRT/FNIRT) - reorient_anat = 'reorient_anat' - reorient_cmd = '{fslpre}fslreorient2std {raw_anat} {adir}/{ra_raw}'.format(fslpre=fsl_prefix,raw_anat=raw_anat, adir=odir, ra_raw=reorient_anat) - print('Reorienting anatomicals using fslreorient2std, cmd like: \n%s' % reorient_cmd) - if do and sp.call(reorient_cmd, shell=True) != 0: - raise IOError('Error calling fslreorient2std on raw anatomical') - reorient_cmd = '{fslpre}fslreorient2std {bet_anat} {adir}/{ra_raw}_brain'.format(fslpre=fsl_prefix,bet_anat=bet_anat, adir=odir, ra_raw=reorient_anat) - if do and sp.call(reorient_cmd, shell=True) != 0: - raise IOError('Error calling fslreorient2std on brain-extracted anatomical') - - ra_betmask = reorient_anat + "_brainmask" - reorient_cmd = '{fslpre}fslreorient2std {bet_anat} {adir}/{ra_betmask}'.format(fslpre=fsl_prefix,bet_anat=betmask_anat, adir=odir, ra_betmask=ra_betmask) - if do and sp.call(reorient_cmd, shell=True) != 0: - raise IOError('Error calling fslreorient2std on brain-extracted mask') - - fsldir = os.environ['FSLDIR'] - standard = '%s/data/standard/MNI152_T1_1mm'%fsldir - bet_standard = '%s_brain'%standard - standardmask = '%s_mask_dil'%bet_standard - cout = 'mni2anat' #stem of the filenames of the transform estimates - - # initial affine anatomical-to-standard registration using FLIRT. required, as the output xfm is used as a start by FNIRT. - flirt_cmd = '{fslpre}flirt -in {bet_standard} -ref {adir}/{ra_raw}_brain -dof 6 -omat {adir}/{cout}_flirt' - flirt_cmd = flirt_cmd.format(fslpre=fsl_prefix, ra_raw=reorient_anat, bet_standard=bet_standard, adir=odir, cout=cout) - print('Running FLIRT to estimate initial affine transform with command:\n%s'%flirt_cmd) - if do and sp.call(flirt_cmd, shell=True) != 0: - raise IOError('Error calling FLIRT with command: %s' % flirt_cmd) - - # FNIRT mni-to-anat transform estimation cmd (does not apply any transform, but generates estimate [cout]) - # the MNI152 2mm config is used even though we're referencing 1mm, per this FSL list post: - # https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=FSL;d14e5a9d.1105 - cmd = '{fslpre}fnirt --in={standard} --ref={ad}/{ra_raw} --refmask={ad}/{refmask} --aff={ad}/{cout}_flirt --cout={ad}/{cout}_fnirt --fout={ad}/{cout}_field --iout={ad}/{cout}_iout --config=T1_2_MNI152_2mm' - cmd = cmd.format(fslpre=fsl_prefix, ra_raw=reorient_anat, standard=standard, refmask=ra_betmask, ad=odir, anat_dir=anat_dir, cout=cout) - print('Running FNIRT to estimate transform, using the following command... this can take a while:\n%s'%cmd) - if do and sp.call(cmd, shell=True) != 0: - raise IOError('Error calling fnirt with cmd: %s'%cmd) - - [pts, polys] = db.get_surf(subject,"fiducial",merge="True") - - #print('raw anatomical: %s\nbet anatomical: %s\nflirt cmd:%s\nfnirt cmd: %s\npts: %s' % (raw_anat,bet_anat,flirt_cmd,cmd,pts)) - - # take the reoriented anatomical, get its affine coord transform, invert this, and save it - reo_xfmnm = 'reorient_inv' - # need to change this line, as the reoriented anatomical is not in the db but in /tmp now - # re_anat = db.get_anat(subject,reorient_anat) - reo_anat_fn = '{odir}/{reorient_anat}.nii.gz'.format(odir=odir,reorient_anat=reorient_anat) - # print(reo_anat_fn) - # since the reoriented anatomicals aren't stored in the db anymore, db.get_anat() will not work (?) - re_anat = nib.load(reo_anat_fn) - reo_xfm = Transform(np.linalg.inv(re_anat.get_affine()),re_anat) - reo_xfm.save(subject,reo_xfmnm,"coord") - - # get the reoriented anatomical's qform and its inverse, they will be needed later - aqf = re_anat.get_qform() - aqfinv = np.linalg.inv(aqf) - - # load the warp field data as a volume - # since it's not in the db anymore but in /tmp instead of: - # warp = db.get_anat(subject,'%s_field'%cout) - # it's this: - warp_fn = '{ad}/{cout}_field.nii.gz'.format(ad=odir,cout=cout) - # print warp_fn - warp = nib.load(warp_fn) - wd = warp.get_data() - # need in (t,z,y,x) order - wd = np.swapaxes(wd,0,3) # x <--> t - wd = np.swapaxes(wd,1,2) # y <--> z - wv = Volume(wd,subject,reo_xfmnm) - - # now do the mapping! this gets the warp field values at the corresponding points - # (uses fiducial surface by default) - warpvd = wv.map(projection="lanczos") - - # reshape into something sensible - warpverts_L = [vs for vs in np.swapaxes(warpvd.left,0,1)] - warpverts_R = [vs for vs in np.swapaxes(warpvd.right,0,1)] - warpverts_ordered = np.concatenate((warpverts_L, warpverts_R)) - - # append 1s for matrix multiplication (coordinate transformation) - o = np.ones((len(pts),1)) - pad_pts = np.append(pts, o, axis=1) - - # print pts, len(pts), len(pts[0]), warpverts_ordered, len(warpverts_ordered), pad_pts, len(pad_pts), pad_pts[0] - - # transform vertex coords from mm to vox using the anat's qform - voxcoords = [aqfinv.dot(padpt) for padpt in pad_pts] - # add the offsets specified in the warp at those locations (ignoring the 1s here) - mnivoxcoords = [voxcoords[n][:-1] + warpverts_ordered[n] for n in range(len(voxcoords))] - # re-pad for matrix multiplication - pad_mnivox = np.append(mnivoxcoords, o, axis=1) - - # multiply by the standard's qform to recover mm coords - std = nib.load('%s.nii.gz'%standard) - stdqf = std.get_qform() - mni_coords = np.array([stdqf.dot(padmni)[:-1] for padmni in pad_mnivox]) - - # some debug output - # print pts, mni_coords - # print pts[0], mni_coords[0] - # print len(pts), len(mni_coords) - # print type(pts), type(pts[0][0]), type(mni_coords) - - # now split mni_coords into left and right arrays for saving - nverts_L = len(warpverts_L) - #print nverts_L - left = mni_coords[:nverts_L] - right = mni_coords[nverts_L:] - #print len(left), len(right) - - mni_surfinfo_fn = db.get_paths(subject)['surfinfo'].format(type='mnicoords',opts='') - print('Saving mni coordinates as a surfinfo...') - np.savez(mni_surfinfo_fn,leftpts=left,rightpts=right) - - return (pts, mni_coords) - - finally: - pass - - def autotweak(subject, xfmname): """Tweak an alignment using the FLIRT boundary-based alignment (BBR) from FSL. Ideally this function should actually use a limited search range, but it doesn't. diff --git a/cortex/surfinfo.py b/cortex/surfinfo.py index cb9ffca95..3ec3825dd 100644 --- a/cortex/surfinfo.py +++ b/cortex/surfinfo.py @@ -153,3 +153,153 @@ def make_surface_graph(tris): lines.append(pts[pbnd,:2]) np.savez(outfile, lines=lines, ismwalls=ismwalls) + +def mni_nl(subject, do=True): + """Create an automatic alignment of an anatomical image to the MNI standard. + + This function does the following: + 1) Re-orders orientation labels on anatomical images using fslreorient2std (without modifying the existing files) + 2) Calls FLIRT + 3) Calls FNIRT with the transform estimated by FLIRT as the specified + 4) Gets the resulting warp field, samples it at each vertex location, and calculates MNI coordinates. + 5) Saves these coordinates as a surfinfo file in the db. + + Parameters + ---------- + subject : str + Subject identifier. + do : bool + Actually execute the commands (True), or just print them (False, useful for debugging). + """ + + import nibabel as nib + + from .options import config + from .dataset import Volume + + fsl_prefix = config.get("basic", "fsl_prefix") + cache = tempfile.mkdtemp() + + print('anat_to_mni, subject: %s' % subject) + + raw_anat = db.get_anat(subject, type='raw').get_filename() + bet_anat = db.get_anat(subject, type='brainmask').get_filename() + betmask_anat = db.get_anat(subject, type='brainmask_mask').get_filename() + anat_dir = os.path.dirname(raw_anat) + odir = cache + + # stem for the reoriented-into-MNI anatomical images (required by FLIRT/FNIRT) + reorient_anat = 'reorient_anat' + reorient_cmd = '{fslpre}fslreorient2std {raw_anat} {adir}/{ra_raw}'.format(fslpre=fsl_prefix,raw_anat=raw_anat, adir=odir, ra_raw=reorient_anat) + print('Reorienting anatomicals using fslreorient2std, cmd like: \n%s' % reorient_cmd) + if do and sp.call(reorient_cmd, shell=True) != 0: + raise IOError('Error calling fslreorient2std on raw anatomical') + + reorient_cmd = '{fslpre}fslreorient2std {bet_anat} {adir}/{ra_raw}_brain'.format(fslpre=fsl_prefix,bet_anat=bet_anat, adir=odir, ra_raw=reorient_anat) + if do and sp.call(reorient_cmd, shell=True) != 0: + raise IOError('Error calling fslreorient2std on brain-extracted anatomical') + + ra_betmask = reorient_anat + "_brainmask" + reorient_cmd = '{fslpre}fslreorient2std {bet_anat} {adir}/{ra_betmask}'.format(fslpre=fsl_prefix,bet_anat=betmask_anat, adir=odir, ra_betmask=ra_betmask) + + if do and sp.call(reorient_cmd, shell=True) != 0: + raise IOError('Error calling fslreorient2std on brain-extracted mask') + + fsldir = os.environ['FSLDIR'] + standard = '%s/data/standard/MNI152_T1_1mm'%fsldir + bet_standard = '%s_brain'%standard + standardmask = '%s_mask_dil'%bet_standard + cout = 'mni2anat' #stem of the filenames of the transform estimates + + # initial affine anatomical-to-standard registration using FLIRT. required, as the output xfm is used as a start by FNIRT. + flirt_cmd = '{fslpre}flirt -in {bet_standard} -ref {adir}/{ra_raw}_brain -dof 6 -omat {adir}/{cout}_flirt' + flirt_cmd = flirt_cmd.format(fslpre=fsl_prefix, ra_raw=reorient_anat, bet_standard=bet_standard, adir=odir, cout=cout) + print('Running FLIRT to estimate initial affine transform with command:\n%s'%flirt_cmd) + if do and sp.call(flirt_cmd, shell=True) != 0: + raise IOError('Error calling FLIRT with command: %s' % flirt_cmd) + + # FNIRT mni-to-anat transform estimation cmd (does not apply any transform, but generates estimate [cout]) + # the MNI152 2mm config is used even though we're referencing 1mm, per this FSL list post: + # https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=FSL;d14e5a9d.1105 + cmd = '{fslpre}fnirt --in={standard} --ref={ad}/{ra_raw} --refmask={ad}/{refmask} --aff={ad}/{cout}_flirt --cout={ad}/{cout}_fnirt --fout={ad}/{cout}_field --iout={ad}/{cout}_iout --config=T1_2_MNI152_2mm' + cmd = cmd.format(fslpre=fsl_prefix, ra_raw=reorient_anat, standard=standard, refmask=ra_betmask, ad=odir, anat_dir=anat_dir, cout=cout) + print('Running FNIRT to estimate transform, using the following command... this can take a while:\n%s'%cmd) + if do and sp.call(cmd, shell=True) != 0: + raise IOError('Error calling fnirt with cmd: %s'%cmd) + + pts, polys = db.get_surf(subject,"fiducial",merge="True") + + #print('raw anatomical: %s\nbet anatomical: %s\nflirt cmd:%s\nfnirt cmd: %s\npts: %s' % (raw_anat,bet_anat,flirt_cmd,cmd,pts)) + + # take the reoriented anatomical, get its affine coord transform, invert this, and save it + reo_xfmnm = 'reorient_inv' + # need to change this line, as the reoriented anatomical is not in the db but in /tmp now + # re_anat = db.get_anat(subject,reorient_anat) + reo_anat_fn = '{odir}/{reorient_anat}.nii.gz'.format(odir=odir,reorient_anat=reorient_anat) + # print(reo_anat_fn) + # since the reoriented anatomicals aren't stored in the db anymore, db.get_anat() will not work (?) + re_anat = nib.load(reo_anat_fn) + reo_xfm = Transform(np.linalg.inv(re_anat.get_affine()),re_anat) + reo_xfm.save(subject,reo_xfmnm,"coord") + + # get the reoriented anatomical's qform and its inverse, they will be needed later + aqf = re_anat.get_qform() + aqfinv = np.linalg.inv(aqf) + + # load the warp field data as a volume + # since it's not in the db anymore but in /tmp instead of: + # warp = db.get_anat(subject,'%s_field'%cout) + # it's this: + warp_fn = '{ad}/{cout}_field.nii.gz'.format(ad=odir,cout=cout) + # print warp_fn + warp = nib.load(warp_fn) + wd = warp.get_data() + # need in (t,z,y,x) order + wd = np.swapaxes(wd,0,3) # x <--> t + wd = np.swapaxes(wd,1,2) # y <--> z + wv = Volume(wd,subject,reo_xfmnm) + + # now do the mapping! this gets the warp field values at the corresponding points + # (uses fiducial surface by default) + warpvd = wv.map(projection="lanczos") + + # reshape into something sensible + warpverts_L = [vs for vs in np.swapaxes(warpvd.left,0,1)] + warpverts_R = [vs for vs in np.swapaxes(warpvd.right,0,1)] + warpverts_ordered = np.concatenate((warpverts_L, warpverts_R)) + + # append 1s for matrix multiplication (coordinate transformation) + o = np.ones((len(pts),1)) + pad_pts = np.append(pts, o, axis=1) + + # print pts, len(pts), len(pts[0]), warpverts_ordered, len(warpverts_ordered), pad_pts, len(pad_pts), pad_pts[0] + + # transform vertex coords from mm to vox using the anat's qform + voxcoords = [aqfinv.dot(padpt) for padpt in pad_pts] + # add the offsets specified in the warp at those locations (ignoring the 1s here) + mnivoxcoords = [voxcoords[n][:-1] + warpverts_ordered[n] for n in range(len(voxcoords))] + # re-pad for matrix multiplication + pad_mnivox = np.append(mnivoxcoords, o, axis=1) + + # multiply by the standard's qform to recover mm coords + std = nib.load('%s.nii.gz'%standard) + stdqf = std.get_qform() + mni_coords = np.array([stdqf.dot(padmni)[:-1] for padmni in pad_mnivox]) + + # some debug output + # print pts, mni_coords + # print pts[0], mni_coords[0] + # print len(pts), len(mni_coords) + # print type(pts), type(pts[0][0]), type(mni_coords) + + # now split mni_coords into left and right arrays for saving + nverts_L = len(warpverts_L) + #print nverts_L + left = mni_coords[:nverts_L] + right = mni_coords[nverts_L:] + #print len(left), len(right) + + mni_surfinfo_fn = db.get_paths(subject)['surfinfo'].format(type='mnicoords',opts='') + print('Saving mni coordinates as a surfinfo...') + np.savez(mni_surfinfo_fn,leftpts=left,rightpts=right) + From cb588edb8ebc74ca2d2e21ea08c5ce16d94d4740 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Fri, 13 Nov 2015 18:12:15 -0800 Subject: [PATCH 28/35] changed the native space display to be voxels, not mm. changes to both the picker and the function that reacts to coord space radio button. --- cortex/webgl/resources/js/facepick.js | 15 ++++++++------- cortex/webgl/resources/js/mriview.js | 19 +++++++++++++------ 2 files changed, 21 insertions(+), 13 deletions(-) diff --git a/cortex/webgl/resources/js/facepick.js b/cortex/webgl/resources/js/facepick.js index 2ccb68479..f995073b5 100644 --- a/cortex/webgl/resources/js/facepick.js +++ b/cortex/webgl/resources/js/facepick.js @@ -159,6 +159,7 @@ FacePick.prototype = { var p = this._pick(x, y); if (p) { var vec = this.viewer.uniforms.volxfm.value[0].multiplyVector3(p.pos.clone()); + console.log("Picked vertex "+p.ptidx+" in "+p.hemi+" hemisphere, distance="+p.dist+", voxel=["+vec.x+","+vec.y+","+vec.z+"]"); if (p.hemi==="left") hem = this.viewer.meshes.left.geometry ; if (p.hemi==="right") @@ -166,20 +167,20 @@ FacePick.prototype = { space = $(this.viewer.object).find(".radio:checked").val(); if (space==="magnet") { - coordarray = hem.attributes.position ; + mnix = vec.x ; + mniy = vec.y ; + mniz = vec.z ; $(this.viewer.object).find("#coordsys_mag").prop('checked',true) ; } else { //mni or undefined coordarray = hem.attributes.mnicoords ; + mniidx = (p.ptidx)*coordarray.itemSize ; + mnix = coordarray.array[mniidx] ; + mniy = coordarray.array[mniidx+1] ; + mniz = coordarray.array[mniidx+2] ; $(this.viewer.object).find("#coordsys_mni").prop('checked',true) ; } - mniidx = (p.ptidx)*coordarray.itemSize ; - - mnix = coordarray.array[mniidx] ; - mniy = coordarray.array[mniidx+1] ; - mniz = coordarray.array[mniidx+2] ; - this.addMarker(p.hemi, p.ptidx, keep); $(this.viewer.object).find("#mnibox").show() ; $(this.viewer.object).find("#mniX").val(mnix.toFixed(2)) ; diff --git a/cortex/webgl/resources/js/mriview.js b/cortex/webgl/resources/js/mriview.js index f5626d748..cb9ae3f89 100644 --- a/cortex/webgl/resources/js/mriview.js +++ b/cortex/webgl/resources/js/mriview.js @@ -1116,18 +1116,25 @@ var mriview = (function(module) { if (space==="magnet") { coordarray = h.attributes.position ; $(this.object).find("#coordsys_mag").prop('checked',true) ; + mniidx = (ptidx)*coordarray.itemSize ; + px = coordarray.array[mniidx] ; + py = coordarray.array[mniidx+1] ; + pz = coordarray.array[mniidx+2] ; + var coord = new THREE.Vector3(px, py, pz); + var vec = this.uniforms.volxfm.value[0].multiplyVector3(coord); + mnix = vec.x ; + mniy = vec.y ; + mniz = vec.z ; } else { //mni or undefined coordarray = h.attributes.mnicoords ; $(this.object).find("#coordsys_mni").prop('checked',true) ; + mniidx = (ptidx)*coordarray.itemSize ; + mnix = coordarray.array[mniidx] ; + mniy = coordarray.array[mniidx+1] ; + mniz = coordarray.array[mniidx+2] ; } - mniidx = (ptidx)*coordarray.itemSize ; - - mnix = coordarray.array[mniidx] ; - mniy = coordarray.array[mniidx+1] ; - mniz = coordarray.array[mniidx+2] ; - $(this.object).find("#mniX").val(mnix.toFixed(2)) ; $(this.object).find("#mniY").val(mniy.toFixed(2)) ; $(this.object).find("#mniZ").val(mniz.toFixed(2)) ; From 7fa1317510d435b85e24feaf72d891f311b00c33 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Tue, 17 Nov 2015 11:07:02 -0800 Subject: [PATCH 29/35] removed extraneous manipulation of mni coords --- cortex/brainctm.py | 1 - 1 file changed, 1 deletion(-) diff --git a/cortex/brainctm.py b/cortex/brainctm.py index 7cfffa073..57808b3b5 100644 --- a/cortex/brainctm.py +++ b/cortex/brainctm.py @@ -61,7 +61,6 @@ def __init__(self, subject, decimate=False): fidpolys = set(tuple(f) for f in polyutils.sort_polys(hemi.polys)) flatpolys = set(tuple(f) for f in polyutils.sort_polys(ptpoly[1])) hemi.aux[np.array(list(fidpolys - flatpolys)).astype(int), 0] = 1 - hemi.mni[np.array(list(fidpolys - flatpolys)).astype(int), 0] = 1 #Find the flatmap limits if fleft is not None: From 352b42ad96008a0d0720f8bc867421f25ad63d9c Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Tue, 17 Nov 2015 11:42:10 -0800 Subject: [PATCH 30/35] added fsl_dir option to defaults.cfg --- cortex/defaults.cfg | 1 + 1 file changed, 1 insertion(+) diff --git a/cortex/defaults.cfg b/cortex/defaults.cfg index f0d4aa0db..75c94e9f1 100644 --- a/cortex/defaults.cfg +++ b/cortex/defaults.cfg @@ -2,6 +2,7 @@ default_cmap = RdBu_r default_cmap2D = RdBu_covar fsl_prefix = fsl5.0- +fsl_dir = 'usr/local/fsl' [mayavi_aligner] line_width = 1 From 1e3c52f2ab26719d98360c135a510df598d04ddb Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Tue, 17 Nov 2015 11:47:27 -0800 Subject: [PATCH 31/35] added the standardfile and projection (interpolation method) arguments to mni_nl, wrote some logic to use them, also correctly used the outfile so that get_surfinfo can specify the filename --- cortex/surfinfo.py | 262 +++++++++++++++++++++++++-------------------- 1 file changed, 143 insertions(+), 119 deletions(-) diff --git a/cortex/surfinfo.py b/cortex/surfinfo.py index 3ec3825dd..fcfda3ef5 100644 --- a/cortex/surfinfo.py +++ b/cortex/surfinfo.py @@ -154,7 +154,7 @@ def make_surface_graph(tris): np.savez(outfile, lines=lines, ismwalls=ismwalls) -def mni_nl(subject, do=True): +def mni_nl(outfile, subject, do=True, standardfile='', projection="lanczos"): """Create an automatic alignment of an anatomical image to the MNI standard. This function does the following: @@ -181,125 +181,149 @@ def mni_nl(subject, do=True): cache = tempfile.mkdtemp() print('anat_to_mni, subject: %s' % subject) - - raw_anat = db.get_anat(subject, type='raw').get_filename() - bet_anat = db.get_anat(subject, type='brainmask').get_filename() - betmask_anat = db.get_anat(subject, type='brainmask_mask').get_filename() - anat_dir = os.path.dirname(raw_anat) - odir = cache - - # stem for the reoriented-into-MNI anatomical images (required by FLIRT/FNIRT) - reorient_anat = 'reorient_anat' - reorient_cmd = '{fslpre}fslreorient2std {raw_anat} {adir}/{ra_raw}'.format(fslpre=fsl_prefix,raw_anat=raw_anat, adir=odir, ra_raw=reorient_anat) - print('Reorienting anatomicals using fslreorient2std, cmd like: \n%s' % reorient_cmd) - if do and sp.call(reorient_cmd, shell=True) != 0: - raise IOError('Error calling fslreorient2std on raw anatomical') - - reorient_cmd = '{fslpre}fslreorient2std {bet_anat} {adir}/{ra_raw}_brain'.format(fslpre=fsl_prefix,bet_anat=bet_anat, adir=odir, ra_raw=reorient_anat) - if do and sp.call(reorient_cmd, shell=True) != 0: - raise IOError('Error calling fslreorient2std on brain-extracted anatomical') - - ra_betmask = reorient_anat + "_brainmask" - reorient_cmd = '{fslpre}fslreorient2std {bet_anat} {adir}/{ra_betmask}'.format(fslpre=fsl_prefix,bet_anat=betmask_anat, adir=odir, ra_betmask=ra_betmask) - - if do and sp.call(reorient_cmd, shell=True) != 0: - raise IOError('Error calling fslreorient2std on brain-extracted mask') - - fsldir = os.environ['FSLDIR'] - standard = '%s/data/standard/MNI152_T1_1mm'%fsldir - bet_standard = '%s_brain'%standard - standardmask = '%s_mask_dil'%bet_standard - cout = 'mni2anat' #stem of the filenames of the transform estimates - - # initial affine anatomical-to-standard registration using FLIRT. required, as the output xfm is used as a start by FNIRT. - flirt_cmd = '{fslpre}flirt -in {bet_standard} -ref {adir}/{ra_raw}_brain -dof 6 -omat {adir}/{cout}_flirt' - flirt_cmd = flirt_cmd.format(fslpre=fsl_prefix, ra_raw=reorient_anat, bet_standard=bet_standard, adir=odir, cout=cout) - print('Running FLIRT to estimate initial affine transform with command:\n%s'%flirt_cmd) - if do and sp.call(flirt_cmd, shell=True) != 0: - raise IOError('Error calling FLIRT with command: %s' % flirt_cmd) - - # FNIRT mni-to-anat transform estimation cmd (does not apply any transform, but generates estimate [cout]) - # the MNI152 2mm config is used even though we're referencing 1mm, per this FSL list post: - # https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=FSL;d14e5a9d.1105 - cmd = '{fslpre}fnirt --in={standard} --ref={ad}/{ra_raw} --refmask={ad}/{refmask} --aff={ad}/{cout}_flirt --cout={ad}/{cout}_fnirt --fout={ad}/{cout}_field --iout={ad}/{cout}_iout --config=T1_2_MNI152_2mm' - cmd = cmd.format(fslpre=fsl_prefix, ra_raw=reorient_anat, standard=standard, refmask=ra_betmask, ad=odir, anat_dir=anat_dir, cout=cout) - print('Running FNIRT to estimate transform, using the following command... this can take a while:\n%s'%cmd) - if do and sp.call(cmd, shell=True) != 0: - raise IOError('Error calling fnirt with cmd: %s'%cmd) + resp = raw_input("This function can take up to two hours to run. Is that okay? [y/n]") + if resp == "y": + raw_anat = db.get_anat(subject, type='raw').get_filename() + bet_anat = db.get_anat(subject, type='brainmask').get_filename() + betmask_anat = db.get_anat(subject, type='brainmask_mask').get_filename() + anat_dir = os.path.dirname(raw_anat) + odir = cache + ext = '.nii.gz' + + if standardfile == '': #user did not specify the standard file, so we will try to determine... + fsl_dir = config.get("basic", "fsl_dir") + if not os.path.isdir(fsl_dir): + if 'FSLDIR' in os.environ.keys() and os.path.isdir(os.environ['FSLDIR']): + fsl_dir = os.environ['FSLDIR'] + else: # no (valid) FSLDIR env var and the config var is missing + print('Could not find FSL directory. Either specify the standard you wish to use, set the fsl_dir config variable, or set $FSLDIR...') + return + # by the time we get here, we know fsl_dir is a directory. + standard = os.path.join(fsl_dir,'data','standard', 'MNI152_T1_1mm') + else: + standard = standardfile + + if not os.path.isfile(standard) and not os.path.isfile(standard+ext): # we don't know if they passed in foo.nii.gz or just foo + print('The standard {sfile} does not exist! Aborting...'.format(sfile=standard)) + return + + bet_standard = '%s_brain'%standard + standardmask = '%s_mask_dil'%bet_standard + cout = 'mni2anat' #stem of the filenames of the transform estimates + full_cout = os.path.join(odir,cout) + + # stem for the reoriented-into-MNI anatomical images (required by FLIRT/FNIRT) + reorient_anat = 'reorient_anat' + outfile = os.path.join(odir, reorient_anat) + outfile_ext = outfile + ext + ra_betmask = outfile + "_brainmask" + + # START DOING THINGS + reorient_cmd = '{fslpre}fslreorient2std {raw_anat} {outfile}'.format(fslpre=fsl_prefix,raw_anat=raw_anat, outfile=outfile) + print('Reorienting anatomicals using fslreorient2std, cmd like: \n%s' % reorient_cmd) + if do and sp.call(reorient_cmd, shell=True) != 0: + raise IOError('Error calling fslreorient2std on raw anatomical') + + reorient_cmd = '{fslpre}fslreorient2std {bet_anat} {outfile}_brain'.format(fslpre=fsl_prefix,bet_anat=bet_anat, outfile=outfile) + if do and sp.call(reorient_cmd, shell=True) != 0: + raise IOError('Error calling fslreorient2std on brain-extracted anatomical') + + reorient_cmd = '{fslpre}fslreorient2std {bet_anat} {ra_betmask}'.format(fslpre=fsl_prefix,bet_anat=betmask_anat, ra_betmask=ra_betmask) + if do and sp.call(reorient_cmd, shell=True) != 0: + raise IOError('Error calling fslreorient2std on brain-extracted mask') + + # initial affine anatomical-to-standard registration using FLIRT. required, as the output xfm is used as a start by FNIRT. + flirt_cmd = '{fslpre}flirt -in {bet_standard} -ref {outfile}_brain -dof 6 -omat {full_cout}_flirt' + flirt_cmd = flirt_cmd.format(fslpre=fsl_prefix, bet_standard=bet_standard, outfile=outfile, full_cout=full_cout) + print('Running FLIRT to estimate initial affine transform with command:\n%s'%flirt_cmd) + if do and sp.call(flirt_cmd, shell=True) != 0: + raise IOError('Error calling FLIRT with command: %s' % flirt_cmd) + + # FNIRT mni-to-anat transform estimation cmd (does not apply any transform, but generates estimate [cout]) + # the MNI152 2mm config is used even though we're referencing 1mm, per this FSL list post: + # https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=FSL;d14e5a9d.1105 + cmd = '{fslpre}fnirt --in={standard} --ref={outfile} --refmask={ra_betmask} --aff={full_cout}_flirt --cout={full_cout}_fnirt --fout={full_cout}_field --iout={full_cout}_iout --config=T1_2_MNI152_2mm' + cmd = cmd.format(fslpre=fsl_prefix, outfile=outfile, standard=standard, ra_betmask=ra_betmask, full_cout=full_cout) + print('Running FNIRT to estimate transform, using the following command... this can take a while:\n%s'%cmd) + if do and sp.call(cmd, shell=True) != 0: + raise IOError('Error calling fnirt with cmd: %s'%cmd) pts, polys = db.get_surf(subject,"fiducial",merge="True") - #print('raw anatomical: %s\nbet anatomical: %s\nflirt cmd:%s\nfnirt cmd: %s\npts: %s' % (raw_anat,bet_anat,flirt_cmd,cmd,pts)) - - # take the reoriented anatomical, get its affine coord transform, invert this, and save it - reo_xfmnm = 'reorient_inv' - # need to change this line, as the reoriented anatomical is not in the db but in /tmp now - # re_anat = db.get_anat(subject,reorient_anat) - reo_anat_fn = '{odir}/{reorient_anat}.nii.gz'.format(odir=odir,reorient_anat=reorient_anat) - # print(reo_anat_fn) - # since the reoriented anatomicals aren't stored in the db anymore, db.get_anat() will not work (?) - re_anat = nib.load(reo_anat_fn) - reo_xfm = Transform(np.linalg.inv(re_anat.get_affine()),re_anat) - reo_xfm.save(subject,reo_xfmnm,"coord") - - # get the reoriented anatomical's qform and its inverse, they will be needed later - aqf = re_anat.get_qform() - aqfinv = np.linalg.inv(aqf) - - # load the warp field data as a volume - # since it's not in the db anymore but in /tmp instead of: - # warp = db.get_anat(subject,'%s_field'%cout) - # it's this: - warp_fn = '{ad}/{cout}_field.nii.gz'.format(ad=odir,cout=cout) - # print warp_fn - warp = nib.load(warp_fn) - wd = warp.get_data() - # need in (t,z,y,x) order - wd = np.swapaxes(wd,0,3) # x <--> t - wd = np.swapaxes(wd,1,2) # y <--> z - wv = Volume(wd,subject,reo_xfmnm) - - # now do the mapping! this gets the warp field values at the corresponding points - # (uses fiducial surface by default) - warpvd = wv.map(projection="lanczos") - - # reshape into something sensible - warpverts_L = [vs for vs in np.swapaxes(warpvd.left,0,1)] - warpverts_R = [vs for vs in np.swapaxes(warpvd.right,0,1)] - warpverts_ordered = np.concatenate((warpverts_L, warpverts_R)) - - # append 1s for matrix multiplication (coordinate transformation) - o = np.ones((len(pts),1)) - pad_pts = np.append(pts, o, axis=1) - - # print pts, len(pts), len(pts[0]), warpverts_ordered, len(warpverts_ordered), pad_pts, len(pad_pts), pad_pts[0] - - # transform vertex coords from mm to vox using the anat's qform - voxcoords = [aqfinv.dot(padpt) for padpt in pad_pts] - # add the offsets specified in the warp at those locations (ignoring the 1s here) - mnivoxcoords = [voxcoords[n][:-1] + warpverts_ordered[n] for n in range(len(voxcoords))] - # re-pad for matrix multiplication - pad_mnivox = np.append(mnivoxcoords, o, axis=1) - - # multiply by the standard's qform to recover mm coords - std = nib.load('%s.nii.gz'%standard) - stdqf = std.get_qform() - mni_coords = np.array([stdqf.dot(padmni)[:-1] for padmni in pad_mnivox]) - - # some debug output - # print pts, mni_coords - # print pts[0], mni_coords[0] - # print len(pts), len(mni_coords) - # print type(pts), type(pts[0][0]), type(mni_coords) - - # now split mni_coords into left and right arrays for saving - nverts_L = len(warpverts_L) - #print nverts_L - left = mni_coords[:nverts_L] - right = mni_coords[nverts_L:] - #print len(left), len(right) - - mni_surfinfo_fn = db.get_paths(subject)['surfinfo'].format(type='mnicoords',opts='') - print('Saving mni coordinates as a surfinfo...') - np.savez(mni_surfinfo_fn,leftpts=left,rightpts=right) + #print('raw anatomical: %s\nbet anatomical: %s\nflirt cmd:%s\nfnirt cmd: %s\npts: %s' % (raw_anat,bet_anat,flirt_cmd,cmd,pts)) + + # take the reoriented anatomical, get its affine coord transform, invert this, and save it + reo_xfmnm = 'reorient_inv' + # need to change this line, as the reoriented anatomical is not in the db but in /tmp now + # re_anat = db.get_anat(subject,reorient_anat) + # since the reoriented anatomicals aren't stored in the db anymore, db.get_anat() will not work (?) + re_anat = nib.load(outfile_ext) # remember, we set outfile to be the full path to the reoriented anat (and outfile_ext to that+.nii.gz) + reo_xfm = Transform(np.linalg.inv(re_anat.get_affine()),re_anat) + reo_xfm.save(subject,reo_xfmnm,"coord") + + # get the reoriented anatomical's qform and its inverse, they will be needed later + aqf = re_anat.get_qform() + aqfinv = np.linalg.inv(aqf) + + # load the warp field data as a volume + # since it's not in the db anymore but in /tmp instead of: + # warp = db.get_anat(subject,'%s_field'%cout) + # it's this: + warp_fn = '{full_cout}_field.nii.gz'.format(full_cout=full_cout) + # print warp_fn + warp = nib.load(warp_fn) + wd = warp.get_data() + # need in (t,z,y,x) order + wd = wd.T + wv = Volume(wd,subject,reo_xfmnm) + # need to delete the reo_xfm + reoxfmpath = os.path.split(reo_xfm.reference.get_filename())[0] + print reoxfmpath + + # now do the mapping! this gets the warp field values at the corresponding points + # (uses fiducial surface by default) + warpvd = wv.map(projection) # passed in via kwargs, defaults to lanczos + + # reshape into something sensible + warpverts_L = [vs for vs in np.swapaxes(warpvd.left,0,1)] + warpverts_R = [vs for vs in np.swapaxes(warpvd.right,0,1)] + warpverts_ordered = np.concatenate((warpverts_L, warpverts_R)) + + # append 1s for matrix multiplication (coordinate transformation) + o = np.ones((len(pts),1)) + pad_pts = np.append(pts, o, axis=1) + + # print pts, len(pts), len(pts[0]), warpverts_ordered, len(warpverts_ordered), pad_pts, len(pad_pts), pad_pts[0] + + # transform vertex coords from mm to vox using the anat's qform + voxcoords = [aqfinv.dot(padpt) for padpt in pad_pts] + # add the offsets specified in the warp at those locations (ignoring the 1s here) + mnivoxcoords = [voxcoords[n][:-1] + warpverts_ordered[n] for n in range(len(voxcoords))] + # re-pad for matrix multiplication + pad_mnivox = np.append(mnivoxcoords, o, axis=1) + + # multiply by the standard's qform to recover mm coords + + if not os.path.splitext(standard)[1]: #will be True iff no extension specified in standard file name + std = nib.load('%s.nii.gz'%standard) + else: + std = nib.load(standard) + stdqf = std.get_qform() + mni_coords = np.array([stdqf.dot(padmni)[:-1] for padmni in pad_mnivox]) + + # some debug output + # print pts, mni_coords + # print pts[0], mni_coords[0] + # print len(pts), len(mni_coords) + # print type(pts), type(pts[0][0]), type(mni_coords) + + # now split mni_coords into left and right arrays for saving + nverts_L = len(warpverts_L) + #print nverts_L + left = mni_coords[:nverts_L] + right = mni_coords[nverts_L:] + #print len(left), len(right) + + print('Saving MNI coordinates as a surfinfo ({outfile})...'.format(outfile=outfile)) + np.savez(outfile,leftpts=left,rightpts=right) From 329f92613c260357d5be4caa6c60fadebf5ce392 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Wed, 18 Nov 2015 11:06:59 -0800 Subject: [PATCH 32/35] fixed bug where outfile was being overwritten and thus the surfinfo was not being saved in teh appropriate place --- cortex/surfinfo.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cortex/surfinfo.py b/cortex/surfinfo.py index fcfda3ef5..eca8461e3 100644 --- a/cortex/surfinfo.py +++ b/cortex/surfinfo.py @@ -154,7 +154,7 @@ def make_surface_graph(tris): np.savez(outfile, lines=lines, ismwalls=ismwalls) -def mni_nl(outfile, subject, do=True, standardfile='', projection="lanczos"): +def mni_nl(si_outfile, subject, do=True, standardfile='', projection="lanczos"): """Create an automatic alignment of an anatomical image to the MNI standard. This function does the following: @@ -324,6 +324,6 @@ def mni_nl(outfile, subject, do=True, standardfile='', projection="lanczos"): right = mni_coords[nverts_L:] #print len(left), len(right) - print('Saving MNI coordinates as a surfinfo ({outfile})...'.format(outfile=outfile)) - np.savez(outfile,leftpts=left,rightpts=right) + print('Saving MNI coordinates as a surfinfo ({si_outfile})...'.format(si_outfile=si_outfile)) + np.savez(si_outfile,leftpts=left,rightpts=right) From dfefa587e5fb6e0324c18d2b40bd708f29469d59 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Wed, 18 Nov 2015 14:59:26 -0800 Subject: [PATCH 33/35] changed first param of mni_nl() back to outfile, and used tmpfile for the intermediate file instead. Also fixed the output to match what other surfinfo functions do: now the hemispheres' data are called left and right (instead of leftpts and rightpts), and I transposed them (plus re-transpose in ctm load) --- cortex/brainctm.py | 4 ++-- cortex/surfinfo.py | 22 +++++++++++----------- 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/cortex/brainctm.py b/cortex/brainctm.py index 57808b3b5..5644149cc 100644 --- a/cortex/brainctm.py +++ b/cortex/brainctm.py @@ -93,8 +93,8 @@ def addMNI(self, **kwargs): print('Adding MNI coords...') npz = db.get_surfinfo(self.subject, type='mnicoords', **kwargs) try: - self.left.mni[:,:-1] = npz['leftpts'] - self.right.mni[:,:-1] = npz['rightpts'] + self.left.mni[:,:-1] = npz['left'].T + self.right.mni[:,:-1] = npz['right'].T except AttributeError: self.left.mni = [] self.right.mni = [] diff --git a/cortex/surfinfo.py b/cortex/surfinfo.py index eca8461e3..04c66a7ec 100644 --- a/cortex/surfinfo.py +++ b/cortex/surfinfo.py @@ -154,7 +154,7 @@ def make_surface_graph(tris): np.savez(outfile, lines=lines, ismwalls=ismwalls) -def mni_nl(si_outfile, subject, do=True, standardfile='', projection="lanczos"): +def mni_nl(outfile, subject, do=True, standardfile='', projection="lanczos"): """Create an automatic alignment of an anatomical image to the MNI standard. This function does the following: @@ -214,17 +214,17 @@ def mni_nl(si_outfile, subject, do=True, standardfile='', projection="lanczos"): # stem for the reoriented-into-MNI anatomical images (required by FLIRT/FNIRT) reorient_anat = 'reorient_anat' - outfile = os.path.join(odir, reorient_anat) - outfile_ext = outfile + ext - ra_betmask = outfile + "_brainmask" + tmpfile = os.path.join(odir, reorient_anat) + tmpfile_ext = tmpfile + ext + ra_betmask = tmpfile + "_brainmask" # START DOING THINGS - reorient_cmd = '{fslpre}fslreorient2std {raw_anat} {outfile}'.format(fslpre=fsl_prefix,raw_anat=raw_anat, outfile=outfile) + reorient_cmd = '{fslpre}fslreorient2std {raw_anat} {outfile}'.format(fslpre=fsl_prefix,raw_anat=raw_anat, outfile=tmpfile) print('Reorienting anatomicals using fslreorient2std, cmd like: \n%s' % reorient_cmd) if do and sp.call(reorient_cmd, shell=True) != 0: raise IOError('Error calling fslreorient2std on raw anatomical') - reorient_cmd = '{fslpre}fslreorient2std {bet_anat} {outfile}_brain'.format(fslpre=fsl_prefix,bet_anat=bet_anat, outfile=outfile) + reorient_cmd = '{fslpre}fslreorient2std {bet_anat} {outfile}_brain'.format(fslpre=fsl_prefix,bet_anat=bet_anat, outfile=tmpfile) if do and sp.call(reorient_cmd, shell=True) != 0: raise IOError('Error calling fslreorient2std on brain-extracted anatomical') @@ -234,7 +234,7 @@ def mni_nl(si_outfile, subject, do=True, standardfile='', projection="lanczos"): # initial affine anatomical-to-standard registration using FLIRT. required, as the output xfm is used as a start by FNIRT. flirt_cmd = '{fslpre}flirt -in {bet_standard} -ref {outfile}_brain -dof 6 -omat {full_cout}_flirt' - flirt_cmd = flirt_cmd.format(fslpre=fsl_prefix, bet_standard=bet_standard, outfile=outfile, full_cout=full_cout) + flirt_cmd = flirt_cmd.format(fslpre=fsl_prefix, bet_standard=bet_standard, outfile=tmpfile, full_cout=full_cout) print('Running FLIRT to estimate initial affine transform with command:\n%s'%flirt_cmd) if do and sp.call(flirt_cmd, shell=True) != 0: raise IOError('Error calling FLIRT with command: %s' % flirt_cmd) @@ -243,7 +243,7 @@ def mni_nl(si_outfile, subject, do=True, standardfile='', projection="lanczos"): # the MNI152 2mm config is used even though we're referencing 1mm, per this FSL list post: # https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=FSL;d14e5a9d.1105 cmd = '{fslpre}fnirt --in={standard} --ref={outfile} --refmask={ra_betmask} --aff={full_cout}_flirt --cout={full_cout}_fnirt --fout={full_cout}_field --iout={full_cout}_iout --config=T1_2_MNI152_2mm' - cmd = cmd.format(fslpre=fsl_prefix, outfile=outfile, standard=standard, ra_betmask=ra_betmask, full_cout=full_cout) + cmd = cmd.format(fslpre=fsl_prefix, outfile=tmpfile, standard=standard, ra_betmask=ra_betmask, full_cout=full_cout) print('Running FNIRT to estimate transform, using the following command... this can take a while:\n%s'%cmd) if do and sp.call(cmd, shell=True) != 0: raise IOError('Error calling fnirt with cmd: %s'%cmd) @@ -257,7 +257,7 @@ def mni_nl(si_outfile, subject, do=True, standardfile='', projection="lanczos"): # need to change this line, as the reoriented anatomical is not in the db but in /tmp now # re_anat = db.get_anat(subject,reorient_anat) # since the reoriented anatomicals aren't stored in the db anymore, db.get_anat() will not work (?) - re_anat = nib.load(outfile_ext) # remember, we set outfile to be the full path to the reoriented anat (and outfile_ext to that+.nii.gz) + re_anat = nib.load(tmpfile_ext) # remember, we set tmpfile to be the full path to the reoriented anat (and tmpfile_ext to that+.nii.gz) reo_xfm = Transform(np.linalg.inv(re_anat.get_affine()),re_anat) reo_xfm.save(subject,reo_xfmnm,"coord") @@ -324,6 +324,6 @@ def mni_nl(si_outfile, subject, do=True, standardfile='', projection="lanczos"): right = mni_coords[nverts_L:] #print len(left), len(right) - print('Saving MNI coordinates as a surfinfo ({si_outfile})...'.format(si_outfile=si_outfile)) - np.savez(si_outfile,leftpts=left,rightpts=right) + print('Saving MNI coordinates as a surfinfo ({outfile})...'.format(outfile=outfile)) + np.savez(outfile,left=left.T,right=right.T) From c31c4e3f82d09b3d9ab900871280d9758c5494b1 Mon Sep 17 00:00:00 2001 From: Arjun Mukerji Date: Wed, 18 Nov 2015 16:51:50 -0800 Subject: [PATCH 34/35] fixed coordinate display and search so that MNI coords are in mm, the other set (which i call 'magnet' in the code) is in voxel indices. both of these are correctly displayed and can be searched against and toggled between. also, the selected vertex's coords are displayed in the box after search. --- cortex/webgl/resources/js/facepick.js | 3 ++ cortex/webgl/resources/js/mriview.js | 54 ++++++++++++++++++++++++--- cortex/webgl/template.html | 6 +-- 3 files changed, 55 insertions(+), 8 deletions(-) diff --git a/cortex/webgl/resources/js/facepick.js b/cortex/webgl/resources/js/facepick.js index f995073b5..ead1a9731 100644 --- a/cortex/webgl/resources/js/facepick.js +++ b/cortex/webgl/resources/js/facepick.js @@ -171,6 +171,7 @@ FacePick.prototype = { mniy = vec.y ; mniz = vec.z ; $(this.viewer.object).find("#coordsys_mag").prop('checked',true) ; + $(this.viewer.object).find(".units").text("vox"); } else { //mni or undefined coordarray = hem.attributes.mnicoords ; @@ -179,6 +180,7 @@ FacePick.prototype = { mniy = coordarray.array[mniidx+1] ; mniz = coordarray.array[mniidx+2] ; $(this.viewer.object).find("#coordsys_mni").prop('checked',true) ; + $(this.viewer.object).find(".units").text("mm"); } this.addMarker(p.hemi, p.ptidx, keep); @@ -312,6 +314,7 @@ FacePick.prototype = { }, addMarker: function(hemi, ptidx, keep) { + console.log('add marker: '+hemi+' idx: '+ptidx) ; var vert = this._getPos(hemi, ptidx); for (var i = 0; i < this.axes.length; i++) { if (keep === true) { diff --git a/cortex/webgl/resources/js/mriview.js b/cortex/webgl/resources/js/mriview.js index cb9ae3f89..b61460588 100644 --- a/cortex/webgl/resources/js/mriview.js +++ b/cortex/webgl/resources/js/mriview.js @@ -1103,6 +1103,7 @@ var mriview = (function(module) { }.bind(this)); } + //this fires when the radio button for selecting coord systems is changed $(this.object).find(".radio").change(function() { //the .radio class is specifically for coord space selection, not all radio buttons space = $(this.object).find(".radio:checked").val(); ptidx = $(this.object).find("#ptidx").val(); @@ -1120,11 +1121,13 @@ var mriview = (function(module) { px = coordarray.array[mniidx] ; py = coordarray.array[mniidx+1] ; pz = coordarray.array[mniidx+2] ; + console.log('mm coords: '+px+', '+py+', '+pz) ; var coord = new THREE.Vector3(px, py, pz); var vec = this.uniforms.volxfm.value[0].multiplyVector3(coord); mnix = vec.x ; mniy = vec.y ; mniz = vec.z ; + $(this.object).find(".units").text("vox"); } else { //mni or undefined coordarray = h.attributes.mnicoords ; @@ -1133,22 +1136,29 @@ var mriview = (function(module) { mnix = coordarray.array[mniidx] ; mniy = coordarray.array[mniidx+1] ; mniz = coordarray.array[mniidx+2] ; + $(this.object).find(".units").text("mm"); } - $(this.object).find("#mniX").val(mnix.toFixed(2)) ; $(this.object).find("#mniY").val(mniy.toFixed(2)) ; $(this.object).find("#mniZ").val(mniz.toFixed(2)) ; } }.bind(this)); + //this fires when the coordinate entry form is submitted -- it's time to look up the nearest vertex $(this.object).find("#mniform").submit(function() { x = $(this.object).find("#mniX").val(); y = $(this.object).find("#mniY").val(); z = $(this.object).find("#mniZ").val(); space = $(this.object).find(".radio:checked").val(); if (space==="magnet") { - var left = this.picker.lkdt.nearest([x, y, z], 1)[0]; - var right = this.picker.rkdt.nearest([x, y, z], 1)[0]; + //do some conversions to eventually get mm coords, which are what's stored in the kdtree + var coord = new THREE.Vector3(x, y, z); + var xfm = this.uniforms.volxfm.value[0] ; + var xfm_inv = new THREE.Matrix4() ; + xfm_inv = xfm_inv.getInverse(xfm); + xfm_inv.multiplyVector3(coord) ; + var left = this.picker.lkdt.nearest([coord.x, coord.y, coord.z], 1)[0]; + var right = this.picker.rkdt.nearest([coord.x, coord.y, coord.z], 1)[0]; } else { //mni or undefined var left = this.picker.mni_lkdt.nearest([x, y, z], 1)[0]; @@ -1156,11 +1166,45 @@ var mriview = (function(module) { } if (left[1] < right[1]) { - this.picker.addMarker("left", left[0][3], false); + ptidx = left[0][3] ; + h = this.meshes.left.geometry ; + hemi = "left" ; } else { - this.picker.addMarker("right", right[0][3], false); + ptidx = right[0][3] ; + h = this.meshes.right.geometry ; + hemi = "right" ; + } + + this.picker.addMarker(hemi, ptidx, false); + + if (space==="magnet") + coordarray = h.attributes.position ; + else //mni or undefined + coordarray = h.attributes.mnicoords ; + + mniidx = (ptidx)*coordarray.itemSize ; + mnix = coordarray.array[mniidx] ; + mniy = coordarray.array[mniidx+1] ; + mniz = coordarray.array[mniidx+2] ; + + //if we're in MNI space, we can just update the box with the coords. if not though, we need to get vox (not mm) coords... + if (space==="magnet") { + var coord2 = new THREE.Vector3(mnix, mniy, mniz); + var vec2 = this.uniforms.volxfm.value[0].multiplyVector3(coord2); + mnix = vec2.x ; + mniy = vec2.y ; + mniz = vec2.z ; } + + //update the hidden fields ptidx and pthem with the picked vertex's info + $(this.object).find("#ptidx").val(ptidx) ; + $(this.object).find("#pthem").val(hemi) ; + + //update coords + $(this.object).find("#mniX").val(mnix.toFixed(2)) ; + $(this.object).find("#mniY").val(mniy.toFixed(2)) ; + $(this.object).find("#mniZ").val(mniz.toFixed(2)) ; return(0); //do not reload page }.bind(this)); diff --git a/cortex/webgl/template.html b/cortex/webgl/template.html index a941c4a0b..d1fba00a0 100644 --- a/cortex/webgl/template.html +++ b/cortex/webgl/template.html @@ -296,9 +296,9 @@ Magnet
- X: mm
- Y: mm
- Z: mm
+ X:
+ Y:
+ Z:
From 388f86bc581e2a16fb8ccdc07ee98296a2bfb3ad Mon Sep 17 00:00:00 2001 From: MarkLescroart Date: Tue, 24 Nov 2015 14:39:59 -0800 Subject: [PATCH 35/35] Small bug fixes in brainctm.py --- cortex/brainctm.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cortex/brainctm.py b/cortex/brainctm.py index 5644149cc..3fdea7097 100644 --- a/cortex/brainctm.py +++ b/cortex/brainctm.py @@ -91,10 +91,10 @@ def addCurvature(self, **kwargs): def addMNI(self, **kwargs): print('Adding MNI coords...') - npz = db.get_surfinfo(self.subject, type='mnicoords', **kwargs) + npz = db.get_surfinfo(self.subject, type='mni_nl', **kwargs) try: - self.left.mni[:,:-1] = npz['left'].T - self.right.mni[:,:-1] = npz['right'].T + self.left.mni[:,:-1] = npz.left.T + self.right.mni[:,:-1] = npz.right.T except AttributeError: self.left.mni = [] self.right.mni = []