MNI coordinate calculation, storage, and display

cortex/brainctm.py

@@ -89,6 +89,17 @@ 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='mni_nl', **kwargs)
+        try:
+            self.left.mni[:,:-1] = npz.left.T
+            self.right.mni[:,:-1] = npz.right.T
+        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. 
 
@@ -177,6 +188,7 @@ 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'''
@@ -199,6 +211,7 @@ def setFlat(self, pts):
 
     def save(self, **kwargs):
         self.ctm.addAttrib(self.aux, 'auxdat')
+        self.ctm.addAttrib(self.mni, 'mnicoords')
         self.ctm.save(**kwargs)
 
         ctm = CTMfile(self.tf.name)
@@ -231,6 +244,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])
@@ -244,6 +258,7 @@ def make_pack(outfile, subj, types=("inflated",), method='raw', level=0,
 
     ctm = BrainCTM(subj, decimate=decimate)
     ctm.addCurvature()
+    ctm.addMNI()
     for name in types:
         ctm.addSurf(name)
 
@@ -275,5 +290,4 @@ def read_pack(ctmfile):
             ctm = CTMfile(tf.name, "r")
             pts, polys, norms = ctm.getMesh()
             meshes.append((pts, polys))
-
     return meshes

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

cortex/surfinfo.py

@@ -153,3 +153,177 @@ def make_surface_graph(tris):
             lines.append(pts[pbnd,:2])
 
     np.savez(outfile, lines=lines, ismwalls=ismwalls)
+
+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:
+    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)
+    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'
+        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=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=tmpfile)
+        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=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)
+
+        # 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=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)
+
+        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)
+        # since the reoriented anatomicals aren't stored in the db anymore, db.get_anat() will not work (?)
+        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")
+
+        # 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,left=left.T,right=right.T)
+

cortex/webgl/resources/css/mriview.css

@@ -136,6 +136,36 @@ 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:12pt;
+    padding:10px;
+    margin:20px;
+    margin-top: 80px ;
+    border-radius:10px;
+    background:rgba(255,255,255,.4);
+    font-weight:bold;
+}
+
+#mnibox input.mniinput {
+    /*border: 0;*/
+    font-size: 10pt;
+    font-weight: bold ;
+    width: 70px ;
+    margin-right: 5px ;
+}
+
+#mnicoords, #mnisubmit {
+    z-index:13;
+}
+
 .datadesc {
     font-size:10pt;
     font-weight:normal;

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);

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 = [];
 
@@ -157,7 +160,36 @@ FacePick.prototype = {
         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")
+                hem = this.viewer.meshes.right.geometry ;
+
+            space = $(this.viewer.object).find(".radio:checked").val();
+            if (space==="magnet") {
+                mnix = vec.x ;
+                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 ;
+                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) ;
+                $(this.viewer.object).find(".units").text("mm");
+            }
+
             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);
@@ -166,6 +198,7 @@ FacePick.prototype = {
                 this.axes[i].obj.parent.remove(this.axes[i].obj);
             }
             this.axes = [];
+            $(this.viewer.object).find("#mnibox").hide() ;
             this.viewer.schedule();
         }
     },
@@ -281,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) {

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();
 }