FEAT model fit plotting sorta working, but I need to think about the

design a bit. Partial model fits are not being calculated correctly.

fsl/data/featimage.py

@@ -10,15 +10,126 @@ analysis.
 """
 
 import os.path as op
+import            glob
 
-import image as fslimage
+import numpy   as np
+
+import nibabel as nib
+
+import image   as fslimage
+
+
+def loadDesignMat(designmat):
+    """Loads a FEAT ``design.mat`` file. Returns a ``numpy`` array
+    containing the design matrix data, where the first dimension
+    corresponds to the data points, and the second to the EVs.
+    """
+
+    matrix = None
+    with open(designmat, 'rt') as f:
+
+        while True:
+            line = f.readline()
+            if line.strip() == '/Matrix':
+                break
+
+        matrix = np.loadtxt(f)
+
+    if matrix is None or matrix.size == 0:
+        raise RuntimeError('{} does not appear to be a '
+                           'valid design.mat file'.format(designmat))
+
+    return matrix
+
+
+def loadDesignCon(designcon):
+    """Loads a FEAT ``design.con`` file. Returns a tuple containing:
+    
+      - A dictionary of ``{contrastnum : name}`` mappings
+    
+      - A list of contrast vectors (each of which is a list itself).
+    """
+
+    matrix       = None
+    numContrasts = 0
+    names        = {}
+    with open(designcon, 'rt') as f:
+
+        while True:
+            line = f.readline().strip()
+
+            if line.startswith('/ContrastName'):
+                tkns       = line.split(None, 1)
+                num        = [c for c in tkns[0] if c.isdigit()]
+                num        = int(''.join(num))
+                name       = tkns[1].strip()
+                names[num] = name
+
+            elif line.startswith('/NumContrasts'):
+                numContrasts = int(line.split()[1])
+
+            elif line == '/Matrix':
+                break
+
+        matrix = np.loadtxt(f)
+
+    if matrix       is None             or \
+       numContrasts != matrix.shape[0]:
+        raise RuntimeError('{} does not appear to be a '
+                           'valid design.con file'.format(designcon))
+
+    # Fill in any missing contrast names
+    if len(names) != numContrasts:
+        for i in range(numContrasts):
+            if i + 1 not in names:
+                names[i + 1] = str(i + 1)
+
+    names     = [names[c + 1] for c in range(numContrasts)]
+    contrasts = []
+
+    for row in matrix:
+        contrasts.append(list(row))
+
+    return names, contrasts
+
+
+def loadDesignFsf(designfsf):
+    """
+    """
+
+    settings = {}
+
+    with open(designfsf, 'rt') as f:
+
+        for line in f.readlines():
+            line = line.strip()
+
+            if not line.startswith('set '):
+                continue
+
+            tkns = line.split(None, 2)
+
+            settings[tkns[1].strip()] = tkns[2]
+    
+    return settings
 
 
 def isFEATData(path):
+    
     keys = ['.feat{}filtered_func_data' .format(op.sep),
             '.gfeat{}filtered_func_data'.format(op.sep)]
 
-    isfeat = any([k in path for k in keys])
+    isfeatdir = any([k in path for k in keys])
+
+    dirname   = op.dirname(path)
+    hasdesfsf = op.exists(op.join(dirname, 'design.fsf'))
+    hasdesmat = op.exists(op.join(dirname, 'design.mat'))
+    hasdescon = op.exists(op.join(dirname, 'design.con'))
+
+    isfeat    = (isfeatdir and
+                 hasdesmat and
+                 hasdescon and
+                 hasdesfsf)
     
     return isfeat
 
@@ -32,36 +143,87 @@ class FEATImage(fslimage.Image):
             raise ValueError('{} does not appear to be data from a '
                              'FEAT analysis'.format(self.dataSource))
 
+        featDir     = op.dirname(self.dataSource)
+        settings    = loadDesignFsf(op.join(featDir, 'design.fsf'))
+        design      = loadDesignMat(op.join(featDir, 'design.mat'))
+        names, cons = loadDesignCon(op.join(featDir, 'design.con'))
 
-    # A FEATImage is an Image which has 
-    # some extra utility methods, something 
-    # like all of the below things:
+        self.__featDir       = featDir
+        self.__design        = design
+        self.__contrastNames = names
+        self.__contrasts     = cons
+        self.__settings      = settings
 
-    # def numParameterEstimates(self):
-    #     return 0
+        self.__pes           = [None] * self.numEVs()
+        self.__copes         = [None] * self.numContrasts()
+
+        
+
+    def getDesign(self):
+        return np.array(self.__design)
+        
     
-    # def numCOPEs(self):
-    #     return 0
+    def numPoints(self):
+        return self.__design.shape[0] 
 
     
-    # def getParameterEstimate(self, num):
-    #     pass
+    def numEVs(self):
+        return self.__design.shape[1]
 
     
-    # def getFullModelFit(self):
-    #     pass
+    def numContrasts(self):
+        return len(self.__contrasts)
 
     
-    # def getPartialModelFIt(self, contrast):
-    #     pass
+    def contrastNames(self):
+        return list(self.__contrastNames)
+
+
+    def contrasts(self):
+        return [list(c) for c in self.__contrasts]
+
+
+    def __getPEFile(self, prefix, ev):
+        prefix = op.join(self.__featDir, 'stats', '{}{}'.format(
+            prefix, ev + 1))
+        return glob.glob('{}.*'.format(prefix))[0]
+
+
+    def getPE(self, ev):
+
+        if self.__pes[ev] is None:
+            pefile = self.__getPEFile('pe', ev)
+            self.__pes[ev] = nib.load(pefile).get_data()
+
+        return self.__pes[ev]
 
     
-    # def getCOPEs(self):
-    #     pass
+    def getCOPE(self, num):
+        if self.__copes[num] is None:
+            copefile = self.__getPEFile('cope', num)
+            self.__copes[num] = nib.load(copefile).get_data()
 
+        return self.__copes[num] 
+        
 
-    # def getZStats(self):
-    #     pass
+    def fit(self, contrast, xyz):
+
+        x, y, z = xyz
+        numEVs  = self.numEVs()
+
+        if len(contrast) != numEVs:
+            raise ValueError('Contrast is wrong length')
+
+        X        = self.__design
+        data     = self.data[x, y, z, :]
+        modelfit = np.zeros(len(data))
+
+        for i in range(numEVs):
+
+            pe        = self.getPE(i)[x, y, z]
+            modelfit += np.dot(X[:, i], pe) * contrast[i]
+
+        return modelfit + data.mean()
 
     
     # def getThresholdedZStats(self):

fsl/data/strings.py

@@ -213,7 +213,8 @@ labels = TypeDict({
     'PlotPanel.xlabel'          : 'X',
     'PlotPanel.ylabel'          : 'Y',
 
-    'TimeSeriesControlPanel.currentFEATSettings' : 'FEAT settings for {}',
+    'TimeSeriesControlPanel.currentFEATSettings' : 'FEAT settings for '
+                                                   'selected overlay ({})',
     
 })
 
@@ -276,6 +277,8 @@ properties = TypeDict({
     'HistogramSeries.showOverlay'     : 'Show 3D histogram overlay',
 
     'FEATTimeSeries.plotFullModelFit' : 'Plot full model fit',
+    'FEATTimeSeries.plotPEFits'       : 'Plot {} fit',
+    'FEATTimeSeries.plotCOPEFits'     : 'Plot COPE {} ({}) fit',
 
     'OrthoEditProfile.selectionSize'          : 'Selection size',
     'OrthoEditProfile.selectionIs3D'          : '3D selection',

fsl/fslview/controls/timeseriescontrolpanel.py

@@ -153,8 +153,27 @@ class TimeSeriesControlPanel(fslpanel.FSLViewPanel):
             displayName=strings.labels[self, 'currentFEATSettings'].format(
                 display.name))
 
-        widg = props.makeWidget(self.__widgets, ts, 'plotFullModelFit')
+        full    = props.makeWidget(     self.__widgets, ts, 'plotFullModelFit')
+        pes     = props.makeListWidgets(self.__widgets, ts, 'plotPEFits')
+        copes   = props.makeListWidgets(self.__widgets, ts, 'plotCOPEFits')
         self.__widgets.AddWidget(
-            widg,
+            full,
             displayName=strings.properties[ts, 'plotFullModelFit'],
             groupName='currentFEATSettings')
+
+        for i, pe in enumerate(pes):
+            peName = 'PE {}'.format(i + 1)
+            self.__widgets.AddWidget(
+                pe,
+                displayName=strings.properties[ts, 'plotPEFits'].format(
+                    peName),
+                groupName='currentFEATSettings') 
+
+
+        copeNames = overlay.contrastNames()
+        for i, (cope, name) in enumerate(zip(copes, copeNames)):
+            self.__widgets.AddWidget(
+                cope,
+                displayName=strings.properties[ts, 'plotCOPEFits'].format(
+                    i, name),
+                groupName='currentFEATSettings') 

fsl/fslview/views/timeseriespanel.py

@@ -38,9 +38,16 @@ class TimeSeries(plotpanel.DataSeries):
         self.coords  = map(int, coords)
         self.data    = overlay.data[coords[0], coords[1], coords[2], :]
 
+        
     def update(self, coords):
-        self.coords = map(int, coords)
+        
+        coords = map(int, coords)
+        if coords == self.coords:
+            return False
+        
+        self.coords = coords
         self.data   = self.overlay.data[coords[0], coords[1], coords[2], :]
+        return True
 
         
     def getData(self):
@@ -63,19 +70,159 @@ class FEATTimeSeries(TimeSeries):
     """
 
     plotFullModelFit = props.Boolean(default=False)
-    # plotResiduals          =            props.Boolean(default=False)
-    # plotParameterEstimates = props.List(props.Boolean(default=False))
-    # plotCopes              = props.List(props.Boolean(default=False))
-
-    # Reduced against what? It has to
-    # be w.r.t. a specific PE/COPE. 
-    # plotReducedData = props.Boolean(default=False)
+    plotPEFits       = props.List(props.Boolean(default=False))
+    plotCOPEFits     = props.List(props.Boolean(default=False))
 
 
     def __init__(self, *args, **kwargs):
         TimeSeries.__init__(self, *args, **kwargs)
         self.name = '{}_{}'.format(type(self).__name__, id(self))
 
+        numEVs   = self.overlay.numEVs()
+        numCOPEs = self.overlay.numContrasts()
+
+        for i in range(numEVs):
+            self.plotPEFits.append(False)
+
+        for i in range(numCOPEs):
+            self.plotCOPEFits.append(False) 
+
+        self.__fullModelTs =  None
+        self.__peTs        = [None] * numEVs
+        self.__copeTs      = [None] * numCOPEs
+        
+        self.addListener('plotFullModelFit',
+                         self.name,
+                         self.__plotFullModelFitChanged)
+        
+        for i, plotPEFit in enumerate(
+                self.plotPEFits.getPropertyValueList()):
+
+            def onChange(ctx, value, valid, name, pe=i):
+                self.__plotPEFitChanged(pe)
+
+            plotPEFit.addListener(self.name, onChange)
+
+        
+        for i, plotCOPEFit in enumerate(
+                self.plotCOPEFits.getPropertyValueList()):
+
+            def onChange(ctx, value, valid, name, cope=i):
+                self.__plotCOPEFitChanged(cope)
+
+            plotCOPEFit.addListener(self.name, onChange)
+ 
+            
+
+    def getModelTimeSeries(self):
+        modelts = []
+
+        if self.plotFullModelFit:
+            modelts.append(self.__fullModelTs)
+
+        for i in range(self.overlay.numEVs()):
+            if self.plotPEFits[i]:
+                modelts.append(self.__peTs[i])
+
+        for i in range(self.overlay.numContrasts()):
+            if self.plotCOPEFits[i]:
+                modelts.append(self.__copeTs[i]) 
+        
+        return modelts
+
+    
+    def __plotCOPEFitChanged(self, copenum):
+        if not self.plotCOPEFits[copenum]:
+            self.__copeTs[copenum] = None
+            return
+
+        con  = self.overlay.contrasts()[copenum]
+
+        copets = FEATModelFitTimeSeries(
+            con,
+            self.tsPanel,
+            self.overlay,
+            self.coords)
+        
+        copets.colour    = (0, 1, 0)
+        copets.alpha     = self.alpha
+        copets.label     = self.label
+        copets.lineWidth = self.lineWidth
+        copets.lineStyle = self.lineStyle
+
+        self.__copeTs[copenum] = copets 
+
+
+    def __plotPEFitChanged(self, evnum):
+        if not self.plotPEFits[evnum]:
+            self.__peTs[evnum] = None
+            return
+
+        con     = [0] * self.overlay.numEVs()
+        con[evnum] = 1
+
+        pets = FEATModelFitTimeSeries(
+            con,
+            self.tsPanel,
+            self.overlay,
+            self.coords)
+        
+        pets.colour    = (1, 0, 0)
+        pets.alpha     = self.alpha
+        pets.label     = self.label
+        pets.lineWidth = self.lineWidth
+        pets.lineStyle = self.lineStyle
+
+        self.__peTs[evnum] = pets
+
+
+    def __plotFullModelFitChanged(self, *a):
+        if not self.plotFullModelFit:
+            self.__fullModelTs = None
+            return
+
+        self.__fullModelTs = FEATModelFitTimeSeries(
+            [1] * self.overlay.numEVs(),
+            self.tsPanel,
+            self.overlay,
+            self.coords)
+        self.__fullModelTs.colour    = (0, 0, 1)
+        self.__fullModelTs.alpha     = self.alpha
+        self.__fullModelTs.label     = self.label
+        self.__fullModelTs.lineWidth = self.lineWidth
+        self.__fullModelTs.lineStyle = self.lineStyle
+
+        
+    def update(self, coords):
+        if not TimeSeries.update(self, coords):
+            return False
+            
+        if self.__fullModelTs is not None:
+            self.__fullModelTs.update(coords)
+
+        return True
+
+
+class FEATModelFitTimeSeries(TimeSeries):
+    
+    def __init__(self, contrast, *args, **kwargs):
+        TimeSeries.__init__(self, *args, **kwargs)
+        self.contrast = contrast
+        self.updateModelFit()
+
+        
+    def update(self, coords):
+        if not TimeSeries.update(self, coords):
+            return
+        self.updateModelFit()
+        
+
+    def updateModelFit(self):
+        x, y, z = self.coords
+
+        self.data = self.overlay.fit(self.contrast, (x, y, z))
+ 
+
 
 class TimeSeriesPanel(plotpanel.PlotPanel):
     """A panel with a :mod:`matplotlib` canvas embedded within.
@@ -154,6 +301,9 @@ class TimeSeriesPanel(plotpanel.PlotPanel):
         prevTs      = self.__currentTs
         prevOverlay = self.__currentOverlay
 
+        if prevTs is not None:
+            prevTs.removeGlobalListener(self._name)
+
         self.__currentTs      = None
         self.__currentOverlay = None
 
@@ -201,6 +351,8 @@ class TimeSeriesPanel(plotpanel.PlotPanel):
             self.__currentTs      = ts
             self.__currentOverlay = overlay
 
+        self.__currentTs.addGlobalListener(self._name, self.draw)
+
         
     def getCurrent(self):
         return self.__currentTs
@@ -209,16 +361,16 @@ class TimeSeriesPanel(plotpanel.PlotPanel):
     def draw(self, *a):
 
         self.__calcCurrent()
-        current = self.getCurrent()
+        current = self.__currentTs
 
         if self.showCurrent and \
            current is not None:
 
-            # Turn current into a list
-            # if it is not already one
-            try:    len(current)
-            except: current = [current]
+            extras = [current]
+
+            if isinstance(current, FEATTimeSeries):
+                extras += current.getModelTimeSeries()
 
-            self.drawDataSeries(current)
+            self.drawDataSeries(extras)
         else:
             self.drawDataSeries()