MNI coordinate calculation, storage, and display

cortex/align.py

@@ -152,6 +152,157 @@ def automatic(subject, xfmname, reference, noclean=False, bbrtype="signed"):
 
     return retval
 
+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 
+    transform called `xfmname` will be automatically stored in the database.
+
+    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.
+
+    Returns
+    -------
+    pts : the vertices of the fiducial surface
+    mnipts : the mni coordinates of those vertices (same shape as pts)
+    """
+
+    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")
+    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()
+        betmask_anat = db.get_anat(subject, type='brainmask_mask').get_filename()
+        anat_dir = os.path.dirname(raw_anat)
+        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, 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 = '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 = 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)
+
+        # 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)
+
+        [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'
+        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)
+
+        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.

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. 
 
@@ -177,6 +184,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 +207,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 +240,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 +254,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 +286,4 @@ def read_pack(ctmfile):
             ctm = CTMfile(tf.name, "r")
             pts, polys, norms = ctm.getMesh()
             meshes.append((pts, polys))
-
     return meshes

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:24pt;
+    padding:10px;
+    margin:20px;
+    margin-top: 90px ;
+    border-radius:10px;
+    background:rgba(255,255,255,.4);
+    font-weight:bold;
+}
+
+#mnibox input {
+    border: 0;
+    font-size: 16pt;
+    font-weight: bold ;
+    width: 100px ;
+    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 = [];
 
@@ -156,8 +159,21 @@ 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] ;
+
             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.figure.notify("pick", this, [vec]);
             if (this.callback !== undefined)
                 this.callback(vec, p.hemi, p.ptidx);
@@ -166,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();
         }
     },

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

cortex/webgl/resources/js/mriview.js

@@ -1103,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){

cortex/webgl/template.html

@@ -286,6 +286,17 @@
 <div id='main'>
 	<div id='ctmload' class='loadmsg'><img src='resources/images/loading.gif'>Loading brain...</div>
 	<div id='dataname'></div>
+	<div id='mnibox'>
+		<div id="mnicoords">
+		MNI coordinates
+		<form name='mniform' id='mniform' action="javascript:void(0);">
+		X: <input type='text' name='mniX' id='mniX' value='X' class='mniinput' />mm<br />
+		Y: <input type='text' name='mniY' id='mniY' value='Y' class='mniinput' />mm<br />
+		Z: <input type='text' name='mniZ' id='mniZ' value='Z' class='mniinput' />mm<br />
+		<input type='submit' id='mnisubmit' name='mnisubmit' value='Go' />
+		</form>
+		</div>
+	</div>
 	<div id='braincover'></div>
 	<canvas id='brain'></canvas>
 </div>

docs/align.rst

@@ -1,3 +1,7 @@
 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.

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: