Made Image.voxToWorld and worldToVox methods more flexible - single...

Made Image.voxToWorld and worldToVox methods more flexible - single points/coordinates can now be transformed. Updated bet to work with new orthopanel. A bunch of little cosmetic tweaks to make my new syntax checker stop complaining.

fsl/data/fslimage.py

@@ -16,7 +16,6 @@ import numpy              as np
 import numpy.linalg       as linalg
 import nibabel            as nib
 import matplotlib.cm      as mplcm
-import matplotlib.colors  as mplcolors
 
 import fsl.props            as props
 import fsl.data.imagefile   as imagefile
@@ -69,45 +68,106 @@ class Image(object):
         """
         """
 
-        points = np.zeros((2,3), dtype=np.float32)
-        points[1,axis] = self.shape[axis]-1
+        points = np.zeros((2, 3), dtype=np.float32)
+        points[1, axis] = self.shape[axis] - 1
         
         tx = self.voxToWorld(points)
 
-        lo,hi = sorted(tx[:,axis])
+        lo, hi = sorted(tx[:, axis])
 
-        lo = lo - self.pixdim[axis]*0.5
-        hi = hi + self.pixdim[axis]*0.5
+        lo = lo - self.pixdim[axis] * 0.5
+        hi = hi + self.pixdim[axis] * 0.5
 
         return (lo, hi)
 
 
-    def worldToVox(self, p):
-        return self._transform(p, self.worldToVoxMat)
+    def worldToVox(self, p, axes=None):
+        """
+        Transforms the given set of points in voxel coordinates to
+        points in world coordinates, according to the affine
+        transformation specified in the image file. Parameters:
+        
+          - p:    N*A array, where N is the number of points, and A
+                  is the number of axes to consider (default: 3)
+        
+          - axes: If None, it is assumed that the input p is a N*3
+                  array, with each point being specified by x,y,z
+                  coordinates. If a single value in the range (0-2),
+                  it is assumed that p is a 1D array. Or, if a
+                  sequence of 2 or 3 values, p must be an array of
+                  N*2 or N*3, respectively.
+        """
+
+        voxp = self._transform(p, self.worldToVoxMat, axes)
+        voxp = np.array(voxp, dtype=np.int64)
+        
+        return voxp
 
 
-    def voxToWorld(self, p):
-        return self._transform(p, self.voxToWorldMat)
+    def voxToWorld(self, p, axes=None):
+        """
+        Transforms the given set of points in world coordinates to
+        points in voxel coordinates, according to the affine
+        transformation specified in the image file.  See the
+        worldToVox docstring for more details.
+        
+        """ 
+        return self._transform(p, self.voxToWorldMat, axes)
 
         
-    def _transform(self, p, a):
+    def _transform(self, p, a, axes):
         """
-        Parameters:
-          - p: N*3 numpy array of (x,y,z) coordinates.
-          - a: 4*4 affine transformation matrix.
+        Transforms the given set of points p according to the given
+        affine transformation a. See the worldToVox docstring for
+        more details.
         """
 
-        t = np.zeros(p.shape, dtype=p.dtype)
+        p = self._fillPoints(p, axes)
+        t = np.zeros((len(p), 3), dtype=p.dtype)
+
+        x = p[:, 0]
+        y = p[:, 1]
+        z = p[:, 2]
 
-        x = p[:,0]
-        y = p[:,1]
-        z = p[:,2]
+        t[:, 0] = x * a[0, 0] + y * a[1, 0] + z * a[2, 0] + a[3, 0]
+        t[:, 1] = x * a[0, 1] + y * a[1, 1] + z * a[2, 1] + a[3, 1]
+        t[:, 2] = x * a[0, 2] + y * a[1, 2] + z * a[2, 2] + a[3, 2]
 
-        t[:,0] = x * a[0,0] + y * a[1,0] + z * a[2,0] + a[3,0]
-        t[:,1] = x * a[0,1] + y * a[1,1] + z * a[2,1] + a[3,1]
-        t[:,2] = x * a[0,2] + y * a[1,2] + z * a[2,2] + a[3,2]
+        if axes is None: axes = [0, 1, 2]
 
-        return t
+        return t[:, axes]
+
+
+    def _fillPoints(self, p, axes):
+        """
+        Turns the given array p into a N*3 array of x,y,z coordinates.
+        The array p may be a 1D array, or an N*2 or N*3 array. 
+        """
+
+        if not isinstance(p, collections.Iterable): p = [p]
+        p = np.array(p)
+
+        if axes is None: return p
+
+        if not isinstance(axes, collections.Iterable): axes = [axes]
+
+        if p.ndim == 1:
+            p = p.reshape((len(p), 1))
+
+        if p.ndim != 2:
+            raise ValueError('Points array must be either one or two '
+                             'dimensions')
+
+        if len(axes) != p.shape[1]:
+            raise ValueError('Points array shape does not match specified '
+                             'number of axes')
+
+        newp = np.zeros((len(p), 3), dtype=p.dtype)
+
+        for i, ax in enumerate(axes):
+            newp[:, ax] = p[:, i]
+
+        return newp
 
         
     def getAttribute(self, name):
@@ -163,12 +223,12 @@ class ImageDisplay(props.HasProperties):
     cmap       = props.ColourMap(default=mplcm.Greys_r,
                                  preNotifyFunc=updateColourMap)
     
-    _view   = props.VGroup(('enabled',
-                            'displayMin',
-                            'displayMax',
-                            'alpha',
-                            'rangeClip',
-                            'cmap'))
+    _view = props.VGroup(('enabled',
+                          'displayMin',
+                          'displayMax',
+                          'alpha',
+                          'rangeClip',
+                          'cmap'))
     _labels = {
         'enabled'    : 'Enabled',
         'displayMin' : 'Min.',
@@ -176,7 +236,7 @@ class ImageDisplay(props.HasProperties):
         'alpha'      : 'Opacity',
         'rangeClip'  : 'Clipping',
         'cmap'       : 'Colour map'
-        }
+    }
 
 
     def __init__(self, image):
@@ -252,13 +312,13 @@ class ImageList(object):
             raise TypeError('item must be a fsl.data.fslimage.Image') 
 
             
-    def __len__     (self):        return self._items.__len__()
-    def __getitem__ (self, key):   return self._items.__getitem__(key)
-    def __iter__    (self):        return self._items.__iter__()
+    def __len__(     self):        return self._items.__len__()
+    def __getitem__( self, key):   return self._items.__getitem__(key)
+    def __iter__(    self):        return self._items.__iter__()
     def __contains__(self, item):  return self._items.__contains__(item)
-    def __eq__      (self, other): return self._items.__eq__(other)
-    def __str__     (self):        return self._items.__str__()
-    def __repr__    (self):        return self._items.__repr__()
+    def __eq__(      self, other): return self._items.__eq__(other)
+    def __str__(     self):        return self._items.__str__()
+    def __repr__(    self):        return self._items.__repr__()
 
  
     def append(self, item):
@@ -288,7 +348,8 @@ class ImageList(object):
     def extend(self, items):
         map(self._validate, items)
         self._items.extend(items)
-        log.debug('List extended: {}'.format(', '.join([str(i) for i in item])))
+        log.debug('List extended: {}'.format(
+            ', '.join([str(i) for i in items])))
         self._updateImageAttributes()
         self._notify()
 
@@ -304,12 +365,11 @@ class ImageList(object):
         self._notify()
 
         
-    def addListener   (self, listener): self._listeners.append(listener)
+    def addListener(   self, listener): self._listeners.append(listener)
     def removeListener(self, listener): self._listeners.remove(listener)
-    def _notify       (self):
+    def _notify(       self):
         for listener in self._listeners:
             try:
                 listener(self)
-            except e:
+            except Exception as e:
                 log.debug('Listener raised exception: {}'.format(e.message))
- 

fsl/fslview/orthopanel.py

@@ -12,10 +12,10 @@ import sys
 if True:
     import logging
     logging.basicConfig(
-        format='%(levelname)8s '\
-               '%(filename)20s '\
-               '%(lineno)4d: '\
-               '%(funcName)s - '\
+        format='%(levelname)8s '
+               '%(filename)20s '
+               '%(lineno)4d: '
+               '%(funcName)s - '
                '%(message)s',
         level=logging.DEBUG)
 
@@ -31,7 +31,12 @@ import fsl.fslview.slicecanvas as slicecanvas
 # the current cursort location in the image space.
 LocationEvent, EVT_LOCATION_EVENT = wxevent.NewEvent()
 
-class OrthoPanel(wx.Panel):
+
+class OrthoPanel(wx.Panel, props.HasProperties):
+
+    displayXCanvas = props.Boolean(default=True)
+    displayYCanvas = props.Boolean(default=True)
+    displayZCanvas = props.Boolean(default=True)
 
     def __init__(self, parent, imageList):
         """
@@ -46,7 +51,7 @@ class OrthoPanel(wx.Panel):
         self.imageList = imageList
 
         wx.Panel.__init__(self, parent)
-        self.SetMinSize((300,100))
+        self.SetMinSize((300, 100))
 
         self.xcanvas = slicecanvas.SliceCanvas(self, imageList, zax=0)
         self.ycanvas = slicecanvas.SliceCanvas(self, imageList, zax=1,
@@ -115,9 +120,9 @@ class OrthoPanel(wx.Panel):
 
         if not ev.LeftIsDown(): return
 
-        mx,my  = ev.GetPositionTuple()
-        source = ev.GetEventObject()
-        w,h = source.GetClientSize()
+        mx, my  = ev.GetPositionTuple()
+        source  = ev.GetEventObject()
+        w, h    = source.GetClientSize()
 
         my = h - my
 
@@ -128,13 +133,13 @@ class OrthoPanel(wx.Panel):
         mx = mx * (source.xmax - source.xmin) / float(w) + source.xmin
         my = my * (source.ymax - source.ymin) / float(h) + source.ymin
 
-        if   source == self.xcanvas: y,z = mx,my
-        elif source == self.ycanvas: x,z = mx,my
-        elif source == self.zcanvas: x,y = mx,my
+        if   source == self.xcanvas: y, z = mx, my
+        elif source == self.ycanvas: x, z = mx, my
+        elif source == self.zcanvas: x, y = mx, my
 
-        self.setLocation(x,y,z)
+        self.setLocation(x, y, z)
 
-        evt = LocationEvent(x=x,y=y,z=z)
+        evt = LocationEvent(x=x, y=y, z=z)
         wx.PostEvent(self, evt)
 
 

fsl/fslview/slicecanvas.py

@@ -6,11 +6,7 @@
 # Author: Paul McCarthy <pauldmccarthy@gmail.com>
 #
 
-import sys
-
 import numpy             as np
-import matplotlib.colors as mplcolors
-import matplotlib.cm     as mplcm
 
 import                      wx
 import wx.glcanvas       as wxgl
@@ -106,7 +102,7 @@ class GLImageData(object):
         xax   = self.canvas.xax
         yax   = self.canvas.yax
         halfx = self.xlen / 2.0
-        halfy = self.ylen / 2.0 
+        halfy = self.ylen / 2.0
 
         # Data stored in the geometry buffer. Defines
         # the geometry of a single voxel, rendered as
@@ -122,22 +118,22 @@ class GLImageData(object):
         # every voxel in one slice.
         xidxs        = np.arange(self.xdim, dtype=np.float32)
         yidxs        = np.arange(self.ydim, dtype=np.float32)
-        yidxs,xidxs  = np.meshgrid(yidxs, xidxs, indexing='ij')
+        yidxs, xidxs = np.meshgrid(yidxs, xidxs, indexing='ij')
 
         # And put them into a single array (the image.voxToWorld
         # method needs xyz coordinates, hence the N*3 shape here)
-        positionData = np.zeros((self.xdim*self.ydim, 3), dtype=np.float32)
-        positionData[:,xax] = xidxs.ravel(order='C')
-        positionData[:,yax] = yidxs.ravel(order='C')
+        positionData = np.zeros((self.xdim * self.ydim, 3), dtype=np.float32)
+        positionData[:, xax] = xidxs.ravel(order='C')
+        positionData[:, yax] = yidxs.ravel(order='C')
 
         # Then we transform them from voxel
         # coordinates to world coordinates
-        positionData = image.voxToWorld(positionData)[:,(xax, yax)]
+        positionData = image.voxToWorld(positionData)[:, (xax, yax)]
 
         # GL buffers for the geometry and position data
         geomData       = geomData    .ravel(order='C')
         positionData   = positionData.ravel(order='C')
-        geomBuffer     = vbo.VBO(geomData    , gl.GL_STATIC_DRAW)
+        geomBuffer     = vbo.VBO(geomData,     gl.GL_STATIC_DRAW)
         positionBuffer = vbo.VBO(positionData, gl.GL_STATIC_DRAW)
 
         # The image buffer, containing the image data itself
@@ -179,8 +175,8 @@ class GLImageData(object):
         # The image data is normalised to lie
         # between 0 and 255, and cast to uint8
         imageData = np.array(image.data, dtype=np.float32)
-        imageData = 255.0*(imageData       - imageData.min()) / \
-                          (imageData.max() - imageData.min())
+        imageData = 255.0 * (imageData       - imageData.min()) / \
+                            (imageData.max() - imageData.min())
         imageData = np.array(imageData, dtype=np.uint8)
 
         # Then flattened, with fortran dimension ordering,
@@ -295,13 +291,20 @@ class GLImageData(object):
 vertex_shader = """
 #version 120
 
-uniform   float alpha;          /* Opacity - constant for a whole image          */
+/* Opacity - constant for a whole image */
+uniform float alpha;
+
+/* Current vertex */
+attribute vec2 inVertex;
 
-attribute vec2  inVertex;       /* Current vertex                                */
-attribute vec2  inPosition;     /* Position of the current voxel                 */
-attribute float voxelValue;     /* Value of the current voxel (in range [0,1])   */
+/* Position of the current voxel */
+attribute vec2 inPos;
 
-varying   float fragVoxelValue; /* Voxel value passed through to fragment shader */ 
+/* Value of the current voxel (in range [0,1]) */
+attribute float voxValue;
+
+/* Voxel value passed through to fragment shader */ 
+varying float fragVoxValue;
 
 void main(void) {
 
@@ -311,10 +314,10 @@ void main(void) {
      * coordinates to screen coordinates).
      */
     gl_Position = gl_ModelViewProjectionMatrix * \
-        vec4(inVertex+inPosition, 0.0, 1.0);
+        vec4(inVertex+inPos, 0.0, 1.0);
 
     /* Pass the voxel value through to the shader. */
-    fragVoxelValue = voxelValue;
+    fragVoxValue = voxValue;
 }
 """
 
@@ -326,11 +329,11 @@ fragment_shader = """
 
 uniform float     alpha; 
 uniform sampler1D colourMap;      /* RGB colour map, stored as a 1D texture */
-varying float     fragVoxelValue;
+varying float     fragVoxValue;
 
 void main(void) {
 
-    vec4  voxTexture = texture1D(colourMap, fragVoxelValue);
+    vec4  voxTexture = texture1D(colourMap, fragVoxValue);
     vec3  voxColour  = voxTexture.rgb;
     float voxAlpha   = voxTexture.a;
 
@@ -342,6 +345,7 @@ void main(void) {
 }
 """
 
+
 class SliceCanvas(wxgl.GLCanvas):
     """
     A wx.glcanvas.GLCanvas which may be used to display a single 2D slice from
@@ -549,8 +553,8 @@ class SliceCanvas(wxgl.GLCanvas):
 
         # Indices of all vertex/fragment shader parameters 
         self.inVertexPos   = gl.glGetAttribLocation( self.shaders, 'inVertex')
-        self.voxelValuePos = gl.glGetAttribLocation( self.shaders, 'voxelValue')
-        self.inPositionPos = gl.glGetAttribLocation( self.shaders, 'inPosition')
+        self.voxelValuePos = gl.glGetAttribLocation( self.shaders, 'voxValue')
+        self.inPositionPos = gl.glGetAttribLocation( self.shaders, 'inPos')
         self.alphaPos      = gl.glGetUniformLocation(self.shaders, 'alpha')
         self.colourMapPos  = gl.glGetUniformLocation(self.shaders, 'colourMap')
 
@@ -626,15 +630,13 @@ class SliceCanvas(wxgl.GLCanvas):
             zstride = glImageData.zstride
 
             # Figure out which slice we are drawing
-            point = np.zeros((1,3))
-            point[0,self.zax] = self.zpos
-            zi = int(image.worldToVox(point)[0,self.zax])
+            zi = int(image.worldToVox(self.zpos, self.zax))
 
-            if not imageDisplay.enabled:
-                continue 
+            # and if it's out of range, don't draw it
+            if zi < 0 or zi >= zdim: continue 
 
-            if zi < 0 or zi >= zdim:
-                continue
+            # nor if the display is disabled
+            if not imageDisplay.enabled: continue 
 
             # Set up the colour buffer
             gl.glEnable(gl.GL_TEXTURE_1D)

fsl/tools/bet.py

@@ -5,20 +5,17 @@
 # Author: Paul McCarthy <pauldmccarthy@gmail.com>
 #
 
-
-import os
-import sys
-
 from collections import OrderedDict
 
 import wx
 
 import nibabel as nb
 
-import fsl.props           as props
-import fsl.data.imagefile  as imagefile
-import fsl.utils.runwindow as runwindow
-import fsl.utils.imageview as imageview
+import fsl.props              as props
+import fsl.data.imagefile     as imagefile
+import fsl.data.fslimage      as fslimage
+import fsl.utils.runwindow    as runwindow
+import fsl.fslview.orthopanel as orthopanel
 
 runChoices = OrderedDict((
 
@@ -31,15 +28,22 @@ runChoices = OrderedDict((
     ('-B',  'Bias field & neck cleanup (can be useful in SIENA)'),
     ('-Z',  'Improve BET if FOV is very small in Z'),
     ('-F',  'Apply to 4D FMRI data'),
-    ('-A',  'Run bet2 and then betsurf to get additional skull and scalp surfaces'),
+    ('-A',  'Run bet2 and then betsurf to get additional skull and scalp '
+            'surfaces'),
     ('-A2', 'As above, when also feeding in non-brain extracted T2')))
 
 
 class Options(props.HasProperties):
 
-    inputImage           = props.FilePath(exists=True, suffixes=imagefile._allowedExts, required=True)
-    outputImage          = props.FilePath(                                              required=True)
-    t2Image              = props.FilePath(exists=True, suffixes=imagefile._allowedExts, required=lambda i: i.runChoice == '-A2')
+    inputImage           = props.FilePath(
+        exists=True,
+        suffixes=imagefile._allowedExts,
+        required=True)
+    outputImage          = props.FilePath(required=True)
+    t2Image              = props.FilePath(
+        exists=True,
+        suffixes=imagefile._allowedExts,
+        required=lambda i: i.runChoice == '-A2')
 
     runChoice            = props.Choice(runChoices)
 
@@ -157,7 +161,8 @@ optLabels = {
     'thresholdImages'      : 'Apply thresholding to brain and mask image',
     'outputSkull'          : 'Output exterior skull surface image',
     'outputMesh'           : 'Generate brain surface as mesh in .vtk format',
-    'outputSurfaceOverlay' : 'Output brain surface overlaid onto original image',
+    'outputSurfaceOverlay' : 'Output brain surface overlaid onto original '
+                             'image',
     'fractionalIntensity'  : 'Fractional intensity threshold',
     'thresholdGradient'    : 'Threshold gradient',
     'headRadius'           : 'Head radius (mm)',
@@ -169,10 +174,13 @@ optLabels = {
 
 
 optTooltips = {
-    'fractionalIntensity' : 'Smaller values give larger brain outline estimates.',
-    'thresholdGradient'   : 'Positive values give larger brain outline at bottom, smaller at top.',
+    'fractionalIntensity' : 'Smaller values give larger brain outline '
+                            'estimates.',
+    'thresholdGradient'   : 'Positive values give larger brain outline at '
+                            'bottom, smaller at top.',
     'headRadius'          : 'Initial surface sphere is set to half of this.',
-    'centreCoords'        : 'Coordinates (voxels) for centre of initial brain surface sphere.'
+    'centreCoords'        : 'Coordinates (voxels) for centre of initial '
+                            'brain surface sphere.'
 }
 
 
@@ -182,21 +190,23 @@ def selectHeadCentre(opts, button):
     select the head centre location.
     """
 
-    image  = nb.load(opts.inputImage)
-    parent = button.GetTopLevelParent()
-    frame  = imageview.ImageFrame(parent, image.get_data(), opts.inputImage)
-    panel  = frame.panel
+    image     = fslimage.Image(opts.inputImage)
+    imageList = fslimage.ImageList([image])
+    parent    = button.GetTopLevelParent()
+    frame     = orthopanel.OrthoFrame(parent, imageList, opts.inputImage)
+    panel     = frame.panel
 
-    panel.setLocation(
-        opts.xCoordinate,
-        opts.yCoordinate,
-        opts.zCoordinate)
+    voxCoords   = [opts.xCoordinate, opts.yCoordinate, opts.zCoordinate]
+    worldCoords = image.voxToWorld([voxCoords])[0]
+
+    panel.setLocation(*worldCoords)
 
     # Whenever the x/y/z coordinates change on
-    # the Options object,update the dialog view. 
-    def updateViewX(val, *a): panel.setXLocation(val)
-    def updateViewY(val, *a): panel.setYLocation(val)
-    def updateViewZ(val, *a): panel.setZLocation(val)
+    # the Options object,update the orthopanel
+    # location
+    def updateViewX(val, *a): panel.setXLocation(image.voxToWorld(val, axes=0))
+    def updateViewY(val, *a): panel.setYLocation(image.voxToWorld(val, axes=1))
+    def updateViewZ(val, *a): panel.setZLocation(image.voxToWorld(val, axes=2))
 
     optListeners = (
         ('xCoordinate', 'updateViewX_{}'.format(id(panel)), updateViewX),
@@ -218,11 +228,24 @@ def selectHeadCentre(opts, button):
     # And whenever the x/y/z coordinates change
     # on the dialog, update the option values.
     def updateOpts(ev):
-        opts.xCoordinate = ev.x
-        opts.yCoordinate = ev.y
-        opts.zCoordinate = ev.z
+        x = image.worldToVox(ev.x, axes=0)
+        y = image.worldToVox(ev.y, axes=1)
+        z = image.worldToVox(ev.z, axes=2)
+
+        if   x >= image.shape[0]: x = image.shape[0] - 1
+        elif x <  0:              x = 0
+        
+        if   y >= image.shape[1]: y = image.shape[1] - 1
+        elif y <  0:              y = 0
+        
+        if   z >= image.shape[2]: z = image.shape[2] - 1
+        elif z <  0:              z = 0
 
-    panel.Bind(imageview.EVT_LOCATION_EVENT, updateOpts)
+        opts.xCoordinate = x
+        opts.yCoordinate = y
+        opts.zCoordinate = z
+
+    panel.Bind(orthopanel.EVT_LOCATION_EVENT, updateOpts)
 
     # Position the dialog by the button that was clicked
     pos = button.GetScreenPosition()
@@ -254,9 +277,10 @@ betView = props.NotebookGroup((
             props.HGroup(
                 key='centreCoords',
                 children=(
-                    props.Button(text='Select',
-                                 callback=selectHeadCentre,
-                                 enabledWhen=lambda i: i.isValid('inputImage')),
+                    props.Button(
+                        text='Select',
+                        callback=selectHeadCentre,
+                        enabledWhen=lambda i: i.isValid('inputImage')),
                     'xCoordinate',
                     'yCoordinate',
                     'zCoordinate'))))))
@@ -274,9 +298,9 @@ def runBet(parent, opts):
         if exitCode != 0: return
 
         image = nb.load(imagefile.addExt(opts.outputImage))
-        frame = imageview.ImageFrame(window,
-                                     image.get_data(),
-                                     title=opts.outputImage)
+        frame = orthopanel.OrthoFrame(window,
+                                      image.get_data(),
+                                      title=opts.outputImage)
         frame.Show()
         
     runwindow.checkAndRun('BET', opts, parent, Options.genBetCmd,
@@ -290,4 +314,3 @@ FSL_HELPPAGE  = 'bet'
 FSL_OPTIONS   = Options
 FSL_INTERFACE = interface
 FSL_RUNTOOL   = runBet
-