Merge pull request #79 from cindeem/master

ecat fix for mis-ordered mlist

Author: matthew-brett

nibabel/ecat.py

@@ -442,9 +442,9 @@ def get_mlist(self, fileobj):
         return mlist
 
     def get_frame_order(self):
-        """Returns the order of the frames in the file
-        Useful as sometimes Frames are not stored in the file in
-        chronological order, and can be used to extract frames
+        """Returns the order of the frames stored in the file
+        Sometimes Frames are not stored in the file in
+        chronological order, this can be used to extract frames
         in correct order
 
         Returns
@@ -465,20 +465,69 @@ def get_frame_order(self):
         >>> mlist.get_frame_order()
         {0: [0, 16842758.0]}
 
+        
         """
         mlist  = self._mlist
-        ind = mlist[:,0] > 0
-        nframes = ind.shape[0]
-        tmp = mlist[ind,0]
-        tmp.sort()
-        frame_dict = {}
+        ids = mlist[:, 0].copy()
+        n_valid = np.sum(ids > 0)
+        ids[ids <=0] = ids.max() + 1 # put invalid frames at end after sort
+        valid_order = np.argsort(ids)
+        if not all(valid_order == sorted(valid_order)):
+            #raise UserWarning if Frames stored out of order
+            warnings.warn_explicit('Frames stored out of order;'\
+                          'true order = %s\n'\
+                          'frames will be accessed in order '\
+                          'STORED, NOT true order'%(valid_order),
+                          UserWarning,'ecat', 0)
         id_dict = {}
-        for fn, id in enumerate(mlist[ind,0]):
-            mlist_n = np.where(mlist[:,0] == id)[0][0]
-            id_dict.update({fn:[mlist_n,id]})
-        
+        for i in range(n_valid):
+            id_dict[i] = [valid_order[i], ids[valid_order[i]]]
+                
         return id_dict
-    
+
+    def get_series_framenumbers(self):
+        """ Returns framenumber of data as it was collected,
+        as part of a series; not just the order of how it was
+        stored in this or across other files
+
+        For example, if the data is split between multiple files
+        this should give you the true location of this frame as
+        collected in the series
+        (Frames are numbered starting at ONE (1) not Zero)
+
+        Returns
+        -------
+        frame_dict: dict mapping order_stored -> frame in series
+               where frame in series counts from 1; [1,2,3,4...]
+
+        Examples
+        --------
+        >>> import os
+        >>> import nibabel as nib
+        >>> nibabel_dir = os.path.dirname(nib.__file__)
+        >>> from nibabel import ecat
+        >>> ecat_file = os.path.join(nibabel_dir,'tests','data','tinypet.v')
+        >>> img = ecat.load(ecat_file)
+        >>> mlist = img.get_mlist()
+        >>> mlist.get_series_framenumbers()
+        {0: 1}
+
+        
+        
+        """
+        frames_order = self.get_frame_order()
+        nframes = self.hdr['num_frames']
+        mlist_nframes = len(frames_order)
+        trueframenumbers = np.arange(nframes - mlist_nframes, nframes)
+        frame_dict = {}
+        try:
+            for frame_stored, (true_order, _) in frames_order.items():
+                #frame as stored in file -> true number in series
+                frame_dict[frame_stored] = trueframenumbers[true_order]+1
+            return frame_dict
+        except:
+            raise IOError('Error in header or mlist order unknown')
+            
 class EcatSubHeader(object):
 
     _subhdrdtype = subhdr_dtype
@@ -625,7 +674,8 @@ def data_from_fileobj(self, frame=0):
 
 
 class EcatImage(SpatialImage):
-    """This class returns a list of Ecat images, with one image(hdr/data) per frame
+    """This class returns a list of Ecat images,
+    with one image(hdr/data) per frame
     """
     _header = EcatHeader
     header_class = _header
@@ -650,7 +700,8 @@ def __init__(self, subheader):
         def __array__(self):
             ''' Cached read of data from file
             This reads ALL FRAMES into one array, can be memory expensive
-            use subheader.data_from_fileobj(frame) for less memeory intensive reads
+            use subheader.data_from_fileobj(frame) for less memory intensive
+            reads
             '''
             if self._data is None:
                 self._data = np.empty(self.shape)
@@ -664,8 +715,9 @@ def __init__(self, data, affine, header,
                  extra = None, file_map = None):
         """ Initialize Image
 
-        The image is a combination of (array, affine matrix, header, subheader,
-        mlist) with optional meta data in `extra`, and filename / file-like objects
+        The image is a combination of
+        (array, affine matrix, header, subheader, mlist)
+        with optional meta data in `extra`, and filename / file-like objects
         contained in the `file_map`.
 
         Parameters
@@ -739,7 +791,7 @@ def get_affine(self):
 
     def get_frame_affine(self, frame):
         """returns 4X4 affine"""
-        return self._subheader.get_affine(frame=frame)
+        return self._subheader.get_frame_affine(frame=frame)
 
     def get_frame(self,frame):
         return self._subheader.data_from_fileobj(frame)

nibabel/tests/test_ecat.py

@@ -100,10 +100,57 @@ def test_mlist(self):
         dt = dt.newbyteorder('>')
         mats = np.recarray(shape=(32,4), dtype=dt,  buf=dat)
         fid.close()
+        #tests
         assert_true(mats['matlist'][0,0] +  mats['matlist'][0,3] == 31)
         assert_true(mlist.get_frame_order()[0][0] == 0)
         assert_true(mlist.get_frame_order()[0][1] == 16842758.0)
-
+        # test badly ordered mlist
+        badordermlist = mlist
+        badordermlist._mlist = np.array([[  1.68427540e+07,   3.00000000e+00,
+                                            1.20350000e+04,   1.00000000e+00],
+                                         [  1.68427530e+07,   1.20360000e+04,
+                                            2.40680000e+04,   1.00000000e+00],
+                                         [  1.68427550e+07,   2.40690000e+04,
+                                            3.61010000e+04,   1.00000000e+00],
+                                         [  1.68427560e+07,   3.61020000e+04,
+                                            4.81340000e+04,   1.00000000e+00],
+                                         [  1.68427570e+07,   4.81350000e+04,
+                                            6.01670000e+04,   1.00000000e+00],
+                                         [  1.68427580e+07,   6.01680000e+04,
+                                            7.22000000e+04,   1.00000000e+00]])
+        assert_true(badordermlist.get_frame_order()[0][0] == 1)
+
+    def test_mlist_errors(self):
+        fid = open(self.example_file, 'rb')
+        hdr = self.header_class.from_fileobj(fid)
+        hdr['num_frames'] = 6
+        mlist =  self.mlist_class(fid, hdr)    
+        mlist._mlist = np.array([[  1.68427540e+07,   3.00000000e+00,
+                                    1.20350000e+04,   1.00000000e+00],
+                                 [  1.68427530e+07,   1.20360000e+04,
+                                    2.40680000e+04,   1.00000000e+00],
+                                 [  1.68427550e+07,   2.40690000e+04,
+                                    3.61010000e+04,   1.00000000e+00],
+                                 [  1.68427560e+07,   3.61020000e+04,
+                                    4.81340000e+04,   1.00000000e+00],
+                                 [  1.68427570e+07,   4.81350000e+04,
+                                    6.01670000e+04,   1.00000000e+00],
+                                 [  1.68427580e+07,   6.01680000e+04,
+                                    7.22000000e+04,   1.00000000e+00]])        
+        series_framenumbers = mlist.get_series_framenumbers()
+        # first frame stored was actually 2nd frame acquired
+        assert_true(series_framenumbers[0] == 2)
+        order = [series_framenumbers[x] for x in sorted(series_framenumbers)]
+        # true series order is [2,1,3,4,5,6], note counting starts at 1
+        assert_true(order == [2, 1, 3, 4, 5, 6])
+        mlist._mlist[0,0] = 0
+        frames_order = mlist.get_frame_order()
+        neworder =[frames_order[x][0] for x in sorted(frames_order)] 
+        assert_true(neworder == [1, 2, 3, 4, 5])
+        assert_raises(IOError,
+                      mlist.get_series_framenumbers)
+        
+        
 
 class TestEcatSubHeader(TestCase):
     header_class = EcatHeader