Merge branch 'master' into dev

.gitignore

@@ -8,4 +8,6 @@ __pycache__
 .DS_Store
 .idea/*
 test.py
-.gitignore
\ No newline at end of file
+.gitignore
+.idea/*
+test.py

mapVBVD/read_twix_hdr.py

 import re
 import numpy as np
+from scipy.integrate import cumtrapz
 from itertools import chain
 from attrdict import AttrDict
 
@@ -63,94 +64,130 @@ class twix_hdr(AttrDict):
 
 
 def parse_ascconv(buffer):
-    #print(buffer)
-    vararray = re.finditer(r'(?P<name>\S*)\s*=\s*(?P<value>\S*)\n',buffer)
-    #print(vararray)
-    mrprot ={}
-    for v in vararray:        
+    # print(buffer)
+    vararray = re.finditer(r'(?P<name>\S*)\s*=\s*(?P<value>\S*)\n', buffer)
+    # print(vararray)
+    mrprot = AttrDict()
+    for v in vararray:
         try:
             value = float(v.group('value'))
         except ValueError:
             value = v.group('value')
-            
-        #now split array name and index (if present)
-        vvarray = re.finditer(r'(?P<name>\w+)(\[(?P<ix>[0-9]+)\])?',v.group('name'))
-        
+
+        # now split array name and index (if present)
+        vvarray = re.finditer(r'(?P<name>\w+)(\[(?P<ix>[0-9]+)\])?', v.group('name'))
+
         currKey = []
         for vv in vvarray:
-            #print(vv.group('name'))
+            # print(vv.group('name'))
             currKey.append(vv.group('name'))
             if vv.group('ix') is not None:
-                #print(vv.group('ix'))
+                # print(vv.group('ix'))
                 currKey.append(vv.group('ix'))
-        
-        mrprot.update({tuple(currKey):value})
-        
+
+        mrprot.update({tuple(currKey): value})
+
     return mrprot
 
+
 def parse_xprot(buffer):
     xprot = {}
     tokens = re.finditer(r'<Param(?:Bool|Long|String)\."(\w+)">\s*{([^}]*)',buffer)
     tokensDouble = re.finditer(r'<ParamDouble\."(\w+)">\s*{\s*(<Precision>\s*[0-9]*)?\s*([^}]*)',buffer)
     alltokens = chain(tokens,tokensDouble)
-    
+
     for t in alltokens:
-        #print(t.group(1))
-        #print(t.group(2))
+        # print(t.group(1))
+        # print(t.group(2))
         name = t.group(1)
-        
-        value = re.sub(r'("*)|( *<\w*> *[^\n]*)','',t.groups()[-1])
-        value = re.sub(r'[\t\n\r\f\v]*','',value.strip()) #value = re.sub(r'\s*',' ',value) for some bonkers reason this inserts whitespace between all the letters! Just look for other whitespace that \s usually does.        
+
+        value = re.sub(r'("*)|( *<\w*> *[^\n]*)', '', t.groups()[-1])
+        value = re.sub(r'[\t\n\r\f\v]*', '',
+                       value.strip())  # value = re.sub(r'\s*',' ',value) for some bonkers reason this inserts whitespace between all the letters! Just look for other whitespace that \s usually does.
         try:
             value = float(value)
         except ValueError:
             pass
-        
-        xprot.update({name : value})
-        
+
+        xprot.update({name: value})
+
     return xprot
 
+
 def parse_buffer(buffer):
-    reASCCONV = re.compile(r'### ASCCONV BEGIN[^\n]*\n(.*)\s### ASCCONV END ###',re.DOTALL)
-    #print(f'buffer = {buffer[0:10]}')
-    #import pdb; pdb.set_trace()
-    
+    reASCCONV = re.compile(r'### ASCCONV BEGIN[^\n]*\n(.*)\s### ASCCONV END ###', re.DOTALL)
+    # print(f'buffer = {buffer[0:10]}')
+    # import pdb; pdb.set_trace()
+
     ascconv = reASCCONV.search(buffer)
-    #print(f'ascconv = {ascconv}')
+    # print(f'ascconv = {ascconv}')
     if ascconv is not None:
-        prot =  parse_ascconv(ascconv.group(0)) 
+        prot = parse_ascconv(ascconv.group(0))
     else:
-        prot = {}
-        
-    
+        prot = AttrDict()
+
     xprot = reASCCONV.split(buffer)
-    #print(f'xprot = {xprot[0][0:10]}')
+    # print(f'xprot = {xprot[0][0:10]}')
     if xprot is not None:
         xprot = ''.join([found for found in xprot])
         prot2 = parse_xprot(xprot)
-        
+
         prot.update(prot2)
-        
+
     return prot
-    
-def read_twix_hdr(fid):
-    #function to read raw data header information from siemens MRI scanners 
-    #(currently VB and VD software versions are supported and tested).
-    
+
+
+def read_twix_hdr(fid, prot):
+    # function to read raw data header information from siemens MRI scanners
+    # (currently VB and VD software versions are supported and tested).
+
     nbuffers = np.fromfile(fid, dtype=np.uint32, count=1)
-    
-    prot = {}
+
     for _ in range(nbuffers[0]):
         tmpBuff = np.fromfile(fid, dtype=np.uint8, count=10)
         bufname = ''.join([chr(item) for item in tmpBuff])
         bufname = re.match(r'^\w*', bufname).group(0)
-        
-        fid.seek(len(bufname)-9,1)
+
+        fid.seek(len(bufname) - 9, 1)
         buflen = np.fromfile(fid, dtype=np.uint32, count=1)
         tmpBuff = np.fromfile(fid, dtype=np.uint8, count=buflen[0])
         buffer = ''.join([chr(item) for item in tmpBuff])
         buffer = re.sub(r'\n\s*\n', '', buffer)
-        
-        prot.update({bufname:parse_buffer(buffer)})
-    
-    return twix_hdr(prot)
\ No newline at end of file
+        prot.update({bufname: parse_buffer(buffer)})
+
+    rstraj = None
+    # read gridding info
+    if 'alRegridMode' in prot.Meas:
+        regrid_mode = int(prot.Meas.alRegridMode.split(' ')[0])
+        if regrid_mode > 1:
+            ncol = int(prot.Meas.alRegridDestSamples.split(' ')[0])
+            dwelltime = float(prot.Meas.aflRegridADCDuration.split(' ')[0]) / ncol
+            gr_adc = np.zeros(ncol, dtype=np.single)
+            start = float(prot.Meas.alRegridDelaySamplesTime.split(' ')[0])
+            time_adc = start + dwelltime * (np.array(range(ncol)) + 0.5)
+            rampup_time = float(prot.Meas.alRegridRampupTime.split(' ')[0])
+            flattop_time = float(prot.Meas.alRegridFlattopTime.split(' ')[0])
+            rampdown_time = float(prot.Meas.alRegridRampdownTime.split(' ')[0])
+            ixUp = np.where(time_adc < rampup_time)[0]
+            ixFlat = np.setdiff1d(np.where(time_adc <= rampup_time + flattop_time)[0],
+                                 np.where(time_adc < rampup_time)[0])
+            ixDn = np.setdiff1d(np.setdiff1d(list(range(ncol)), ixFlat), ixUp)
+            gr_adc[ixFlat] = 1
+            if regrid_mode == 2:
+                # trapezoidal gradient
+                gr_adc[ixUp] = time_adc[ixUp] / float(prot.Meas.alRegridRampupTime.split(' ')[0])
+                gr_adc[ixDn] = 1 - (time_adc[ixDn] - rampup_time - flattop_time) / rampdown_time
+            elif regrid_mode == 4:
+                gr_adc[ixUp] = np.sin(np.pi / 2 * time_adc[ixUp] / rampup_time)
+                gr_adc[ixDn] = np.sin(np.pi / 2 * (1 + (time_adc[ixDn] - rampup_time - flattop_time) / rampdown_time))
+            else:
+                raise Exception('regridding mode unknown')
+
+            # make sure that gr_adc is always positive (rstraj needs to be strictly monotonic)
+            gr_adc = np.maximum(gr_adc, 1e-4)
+            rstraj = (np.append(0, cumtrapz(gr_adc)) - ncol / 2) / np.sum(gr_adc)
+            rstraj -= np.mean(rstraj[int(ncol/2)-1:int(ncol/2)+1])
+            # scale rstraj by kmax (only works if all slices have same FoV!!!)
+            kmax = prot.MeasYaps[('sKSpace', 'lBaseResolution')] / prot.MeasYaps[('sSliceArray', 'asSlice', '0', 'dReadoutFOV')]
+            rstraj *= kmax
+    return prot, rstraj

requirements.txt

 numpy
 attrdict
+tqdm