Skip to content
GitLab
Commit 542a3b17 authored by Mo Shahdloo's avatar Mo Shahdloo
Browse files

fullSize memory rewrite fixed

parent eddc85cf
Hide whitespace changes
Inline Side-by-side
mapVBVD/__init__.py → core/__init__.py
View file @ 542a3b17
from mapVBVD.mapVBVD import mapVBVD
from core.mapVBVD import mapVBVD
from ._version import get_versions
__version__ = get_versions()['version']
del get_versions
mapVBVD/_version.py → core/_version.py
View file @ 542a3b17
......@@ -43,7 +43,7 @@ def get_config():
cfg.style = "pep440"
cfg.tag_prefix = ""
cfg.parentdir_prefix = ""
cfg.versionfile_source = "mapVBVD/_version.py"
cfg.versionfile_source = "core/_version.py"
cfg.verbose = False
return cfg
......
mapVBVD/mapVBVD.py → core/mapVBVD.py
View file @ 542a3b17
import numpy as np
from dataclasses import dataclass, field
import mapVBVD as pkg
import core as pkg
from attrdict import AttrDict, AttrMap, AttrDefault
from tqdm import tqdm, trange
from mapVBVD.read_twix_hdr import read_twix_hdr
from mapVBVD.twix_map_obj import twix_map_obj
from .read_twix_hdr import read_twix_hdr, twix_hdr
from .twix_map_obj import twix_map_obj
def bitget(number, pos):
def get_bit(number, pos):
return (number >> pos) & 1
......@@ -43,7 +43,7 @@ def loop_mdh_read(fid, version, Nscans, scan, measOffset, measLength, print_prog
mdh_blob = np.zeros((byteMDH, 0), dtype=np.uint8)
szBlob = mdh_blob.shape[1] # pylint: disable=E1136 # pylint/issues/3139
filePos = np.zeros((0), dtype=float)
filePos = np.zeros(0, dtype=float)
fid.seek(cPos, 0)
......@@ -94,7 +94,7 @@ def loop_mdh_read(fid, version, Nscans, scan, measOffset, measLength, print_prog
if ((data_u8[0:3] == u8_000).all()) or (bitMask & bit_0):
# ok, look closer if really all *4* bytes are 0
data_u8[3] = bitget(data_u8[3], 0) # ubit24: keep only 1 bit from the 4th byte
data_u8[3] = get_bit(data_u8[3], 0) # ubit24: keep only 1 bit from the 4th byte
tmp = data_u8[0:4]
tmp.dtype = np.uint32
ulDMALength = float(tmp)
......@@ -107,8 +107,8 @@ def loop_mdh_read(fid, version, Nscans, scan, measOffset, measLength, print_prog
# isEOF = True
break
if (bitMask & bit_5): # MDH_SYNCDATA
data_u8[3] = bitget(data_u8[3], 0) # ubit24: keep only 1 bit from the 4th byte
if bitMask & bit_5: # MDH_SYNCDATA
data_u8[3] = get_bit(data_u8[3], 0) # ubit24: keep only 1 bit from the 4th byte
tmp = data_u8[0:4]
tmp.dtype = np.uint32
ulDMALength = float(tmp)
......@@ -169,10 +169,10 @@ def evalMDH(mdh_blob, version):
Nmeas = mdh_blob.shape[1]
ulPackBit = bitget(mdh_blob[3, :], 2)
ulPackBit = get_bit(mdh_blob[3, :], 2)
ulPCI_rx = set_bit(mdh_blob[3, :], 7, False) # keep 6 relevant bits
ulPCI_rx = set_bit(ulPCI_rx, 8, False)
mdh_blob[3, :] = bitget(mdh_blob[3, :], 1) # ubit24: keep only 1 bit from the 4th byte
mdh_blob[3, :] = get_bit(mdh_blob[3, :], 1) # ubit24: keep only 1 bit from the 4th byte
data_uint32 = np.ascontiguousarray(mdh_blob[0:76, :].transpose())
data_uint32.dtype = np.uint32
......@@ -241,6 +241,7 @@ def evalMDH(mdh_blob, version):
return mdh, mask
def mapVBVD(filename, quiet=False, **kwargs):
if not quiet:
print(f'pymapVBVD version {pkg.__version__}')
......@@ -309,7 +310,7 @@ def mapVBVD(filename, quiet=False, **kwargs):
hdr_len = np.fromfile(fid, dtype=np.uint32, count=1, offset=0)
currTwixObj = AttrDict()
currTwixObjHdr = AttrDict()
currTwixObjHdr = twix_hdr()
# rstraj = 0
# read header
currTwixObjHdr, rstraj = read_twix_hdr(fid, currTwixObjHdr)
......
mapVBVD/read_twix_hdr.py → core/read_twix_hdr.py
View file @ 542a3b17
......@@ -4,17 +4,18 @@ from scipy.integrate import cumtrapz
from itertools import chain
from attrdict import AttrDict
class twix_hdr(AttrDict):
def __init__(self,*args):
def __init__(self, *args):
super().__init__(*args)
@staticmethod
def search_using_tuple(s_terms,key,regex=True):
def search_using_tuple(s_terms, key, regex=True):
if regex:
def regex_tuple(pattern,key_tuple):
def regex_tuple(pattern, key_tuple):
inner_out = []
re_comp = re.compile(pattern,re.IGNORECASE)
re_comp = re.compile(pattern, re.IGNORECASE)
for k in key_tuple:
m = re_comp.search(k)
if m:
......@@ -25,12 +26,12 @@ class twix_hdr(AttrDict):
out = []
for st in s_terms:
out.append(regex_tuple(st,key))
return all(out)
out.append(regex_tuple(st, key))
return all(out)
else:
return all([st in key for st in s_terms])
def search_for_keys(self,search_terms,top_lvl=None,print_flag=True,regex=True):
def search_for_keys(self, search_terms, top_lvl=None, print_flag=True, regex=True):
"""Search header keys for terms.
Args:
......@@ -38,27 +39,27 @@ class twix_hdr(AttrDict):
recursive (optional): Search using regex or for exact strings.
top_lvl (optional) : Specify list of parameter sets to search (e.g. YAPS)
print_flag(optional): If False no output will be printed.
"""
"""
if top_lvl is None:
top_lvl = self.keys()
elif isinstance(top_lvl,str):
top_lvl = [top_lvl,]
elif isinstance(top_lvl, str):
top_lvl = [top_lvl, ]
out = {}
for key in top_lvl:
matching_keys = []
list_of_keys = self[key].keys()
for sub_key in list_of_keys:
if twix_hdr.search_using_tuple(search_terms,sub_key,regex=regex):
if twix_hdr.search_using_tuple(search_terms, sub_key, regex=regex):
matching_keys.append(sub_key)
if print_flag:
print(f'{key}:')
for mk in matching_keys:
print(f'\t{mk}: {self[key][mk]}')
out.update({key:matching_keys})
out.update({key: matching_keys})
return out
......@@ -92,9 +93,9 @@ def parse_ascconv(buffer):
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)
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))
......@@ -170,7 +171,7 @@ def read_twix_hdr(fid, prot):
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])
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:
......@@ -186,8 +187,10 @@ def read_twix_hdr(fid, prot):
# 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])
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')]
# TODO: these are examples of the keys not arranged correctly
kmax = prot.MeasYaps[('sKSpace', 'lBaseResolution')] / prot.MeasYaps[
('sSliceArray', 'asSlice', '0', 'dReadoutFOV')]
rstraj *= kmax
return prot, rstraj
mapVBVD/twix_map_obj.py → core/twix_map_obj.py
View file @ 542a3b17
......@@ -10,14 +10,19 @@ import logging
class twix_map_obj:
@property
def dataSize(self):
self.clean()
# print(self.fullSize)
out = self.fullSize
if out is None:
def fullSize(self):
if self.full_size is None:
self.clean()
out = self.fullSize
# Not yet implemented
return self.full_size
# @fullSize.setter
# def fullSize(self, val):
# self.full_size = val
@property
def dataSize(self):
out = self.fullSize.copy()
if self.removeOS:
ix = self.dataDims.index('Col')
out[ix] = self.NCol / 2
......@@ -106,8 +111,8 @@ class twix_map_obj:
self.skipLin = 0
self.skipPar = 0
self.fullSize[2] = np.maximum(1, self.NLin - self.skipLin)
self.fullSize[3] = np.maximum(1, self.NPar - self.skipPar)
self.full_size[2] = np.maximum(1, self.NLin - self.skipLin)
self.full_size[3] = np.maximum(1, self.NPar - self.skipPar)
@property
def flagRampSampRegrid(self):
......@@ -217,7 +222,7 @@ class twix_map_obj:
self.skipLin = None
self.skipPar = None
self.fullSize = None
self.full_size = None
# Flags
self.flagAverageDim = np.full(16, False, dtype=np.bool)
......@@ -318,7 +323,7 @@ class twix_map_obj:
while not isLastAcqGood and self.NAcq > 0 and cnt < 100:
try:
self.clean()
# self.clean()
self.unsorted(self.NAcq)
isLastAcqGood = True
except Exception as e:
......@@ -404,12 +409,12 @@ class twix_map_obj:
NLinAlloc = np.maximum(1, self.NLin - self.skipLin)
NParAlloc = np.maximum(1, self.NPar - self.skipPar)
self.fullSize = np.array(
self.full_size = np.array(
[self.NCol, self.NCha, NLinAlloc, NParAlloc,
self.NSli, self.NAve, self.NPhs, self.NEco,
self.NRep, self.NSet, self.NSeg, self.NIda,
self.NIdb, self.NIdc, self.NIdd, self.NIde]
)
, dtype=np.int)
nByte = self.NCha * (self.freadInfo.szChannelHeader + 8 * self.NCol)
......@@ -630,20 +635,21 @@ class twix_map_obj:
mem = mem.astype(int)
if outSize is None:
if selRange is None:
selRange = [np.arange(0, self.dataSize[0]).astype(int), np.arange(0, self.dataSize[1]).astype(
int)] # [slice(None,None,None),slice(None,None,None)]
selRange = [slice(None,None,None),slice(None,None,None)]
#[np.arange(0, self.dataSize[0]).astype(int), np.arange(0, self.dataSize[1]).astype(int)]
# [slice(None,None,None),slice(None,None,None)]
else:
selRange[0] = np.arange(0, self.dataSize[0]).astype(int) # slice(None,None,None)
selRange[1] = np.arange(0, self.dataSize[0]).astype(int) # slice(None,None,None)
selRange[0] = np.arange(self.dataSize[0]).astype(int) # slice(None,None,None)
selRange[1] = np.arange(self.dataSize[0]).astype(int) # slice(None,None,None)
outSize = np.append(self.dataSize[0:2], mem.size).astype(int)
selRangeSz = outSize
cIxToTarg = np.arange(0, selRangeSz[2])
cIxToRaw = cIxToTarg
else:
if np.array_equiv(selRange[0], np.arange(0, self.dataSize[0]).astype(int)):
if np.array_equiv(selRange[0], np.arange(self.dataSize[0]).astype(int)):
selRange[0] = slice(None, None, None)
if np.array_equiv(selRange[1], np.arange(0, self.dataSize[1]).astype(int)):
if np.array_equiv(selRange[1], np.arange(self.dataSize[1]).astype(int)):
selRange[1] = slice(None, None, None)
out = np.zeros(outSize, dtype=np.csingle)
......@@ -706,9 +712,9 @@ class twix_map_obj:
fid.seek(mem[k] + szScanHeader, 0)
raw = np.fromfile(fid, dtype=np.float32, count=readSize.prod()).reshape(
(readSize[1], readSize[0])) # do transpose by switching readSize order
# % MiVö: With incomplete files fread() returns less than readSize points. The subsequent reshape will therefore error out.
# % We could check if numel(raw) == prod(readSize), but people recommend exception handling for performance
# % reasons. Do it.
# With incomplete files fread() returns less than readSize points. The subsequent reshape will
# therefore error out. We could check if numel(raw) == prod(readSize), but people recommend
# exception handling for performance reasons. Do it.
try:
raw = (raw[:, 0] + 1j * raw[:, 1]).reshape(readShape, order='F')
except ValueError:
......@@ -722,8 +728,9 @@ class twix_map_obj:
raw = raw.reshape(readShape)
isBrokenRead = True # remember it and bail out later
block[:, :, blockCtr, None] = copy.deepcopy(raw).reshape(np.append(readShape,
1)) # fast serial storage in a cache array - this is probably all very dependent on whether I've got things contiguous in memory. I highly doubt that I have on this first pass. WTC
block[:, :, blockCtr, None] = copy.deepcopy(raw).reshape(np.append(readShape, 1))
# fast serial storage in a cache array - this is probably all very dependent on whether I've got things
# contiguous in memory. I highly doubt that I have on this first pass. WTC
blockCtr += 1
# Do expensive computations and reorderings on the gathered block.
......@@ -868,4 +875,4 @@ class twix_map_obj:
out = np.ascontiguousarray(out.reshape(outSize))
return out if not self.squeeze else np.squeeze(out)
\ No newline at end of file
return out if not self.squeeze else np.squeeze(out)
requirements.txt
View file @ 542a3b17
numpy
attrdict
tqdm
numpy~=1.18.3
attrdict~=2.0.1
tqdm~=4.45.0
dataclasses~=0.7
scipy~=1.4.1
matplotlib~=3.2.1
setuptools~=46.1.3
\ No newline at end of file
tests/test_svs.py
View file @ 542a3b17
import os.path as op
import numpy as np
from mapVBVD import mapVBVD
from core import mapVBVD
test_data_vb = op.join(op.dirname(__file__),'test_data','meas_MID311_STEAM_wref1_FID115674.dat')
test_data_ve = op.join(op.dirname(__file__),'test_data','meas_MID00305_FID74175_VOI_slaser_wref1.dat')
test_data_vb_broken = op.join(op.dirname(__file__),'test_data','meas_MID111_sLaser_broken_FID4873.dat')
test_data_vb = op.join(op.dirname(__file__), 'test_data', 'meas_MID311_STEAM_wref1_FID115674.dat')
test_data_ve = op.join(op.dirname(__file__), 'test_data', 'meas_MID00305_FID74175_VOI_slaser_wref1.dat')
test_data_vb_broken = op.join(op.dirname(__file__), 'test_data', 'meas_MID111_sLaser_broken_FID4873.dat')
def test_vb():
twixObj = mapVBVD(test_data_vb,quiet=True)
assert np.allclose(twixObj.image.fullSize,[4096,32,1,1,1,1,1,1,2,1,1,1,1,1,1,1])
assert np.allclose(twixObj.image.sqzSize(),[4096,32,2])
twixObj = mapVBVD(test_data_vb, quiet=False)
assert np.allclose(twixObj.image.fullSize, [4096, 32, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1])
assert np.allclose(twixObj.image.sqzSize, [2048, 32, 2])
keys = twixObj.search_header_for_keys(('sTXSPEC', 'asNucleusInfo'),top_lvl='MeasYaps')
keys = twixObj.search_header_for_keys(('sTXSPEC', 'asNucleusInfo'), top_lvl='MeasYaps')
assert ('sTXSPEC', 'asNucleusInfo', '0', 'tNucleus') in keys['MeasYaps']
key_value = twixObj.search_header_for_val('MeasYaps',('sTXSPEC', 'asNucleusInfo', '0', 'tNucleus'))
key_value = twixObj.search_header_for_val('MeasYaps', ('sTXSPEC', 'asNucleusInfo', '0', 'tNucleus'))
assert key_value[0] == '"1H"'
def test_ve():
twixObj = mapVBVD(test_data_ve,quiet=True)
assert len(twixObj)==2
twixObj = mapVBVD(test_data_ve, quiet=True)
assert len(twixObj) == 2
twixObj[1].image
assert np.allclose(twixObj[1].image.fullSize,[4096,32,1,1,1,1,1,1,2,1,1,1,1,1,1,1])
assert np.allclose(twixObj[1].image.sqzSize(),[4096,32,2])
assert np.allclose(twixObj[1].image.fullSize, [4096, 32, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1])
assert np.allclose(twixObj[1].image.sqzSize, [2048, 32, 2])
keys = twixObj[1].search_header_for_keys(('sTXSPEC', 'asNucleusInfo'),top_lvl='MeasYaps')
keys = twixObj[1].search_header_for_keys(('sTXSPEC', 'asNucleusInfo'), top_lvl='MeasYaps')
assert ('sTXSPEC', 'asNucleusInfo', '0', 'tNucleus') in keys['MeasYaps']
key_value = twixObj[1].search_header_for_val('MeasYaps',('sTXSPEC', 'asNucleusInfo', '0', 'tNucleus'))
key_value = twixObj[1].search_header_for_val('MeasYaps', ('sTXSPEC', 'asNucleusInfo', '0', 'tNucleus'))
assert key_value[0] == '"1H"'
def test_vb_broken():
twixObj = mapVBVD(test_data_vb_broken,quiet=True)
assert np.allclose(twixObj.image.fullSize,[4096,32,1,1,1,1,1,1,1,97,1,1,1,1,1,1])
assert np.allclose(twixObj.image.sqzSize(),[4096,32,97])
twixObj = mapVBVD(test_data_vb_broken, quiet=True)
assert np.allclose(twixObj.image.fullSize, [4096, 32, 1, 1, 1, 1, 1, 1, 1, 97, 1, 1, 1, 1, 1, 1])
assert np.allclose(twixObj.image.sqzSize, [2048, 32, 97])
keys = twixObj.search_header_for_keys(('sTXSPEC', 'asNucleusInfo'),top_lvl='MeasYaps')
keys = twixObj.search_header_for_keys(('sTXSPEC', 'asNucleusInfo'), top_lvl='MeasYaps')
assert ('sTXSPEC', 'asNucleusInfo', '0', 'tNucleus') in keys['MeasYaps']
key_value = twixObj.search_header_for_val('MeasYaps',('sTXSPEC', 'asNucleusInfo', '0', 'tNucleus'))
assert key_value[0] == '"1H"'
\ No newline at end of file
key_value = twixObj.search_header_for_val('MeasYaps', ('sTXSPEC', 'asNucleusInfo', '0', 'tNucleus'))
assert key_value[0] == '"1H"'
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment