Skip to content

Fixing 16-bits counters overflow correction #1

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
srapacchi wants to merge 5 commits into
base: main
Choose a base branch
from
Open

Fixing 16-bits counters overflow correction #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
053e5b9
Fixing 16-bits counters overflow correction
srapacchi May 1, 2024
aea69f5
Add PMU chunk length for XA6*
srapacchi Jan 28, 2025
26b9bea
Update mapVBVD.m
srapacchi Apr 23, 2025
f744956
Merge branch 'main' into patch-1
srapacchi Apr 23, 2025
34bab3d
Merge pull request #1 from srapacchi/patch-1
srapacchi Apr 23, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
169 changes: 115 additions & 54 deletions mapVBVD.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -199,6 +199,7 @@
arg.bReadRTfeedback = true;
arg.bReadPhaseStab = true;
arg.bReadHeader = true;
arg.bReadPrescan = true;

k=1;
while k <= numel(varargin)
Expand Down Expand Up @@ -295,7 +296,7 @@
% lazy software version check (VB or VD?)
if and(firstInt < 10000, secondInt <= 64)
version = 'vd';
disp('Software version: VD (!?)');
disp('Software version: VD/VE/XA (!?)');

% number of different scans in file stored in 2nd in
NScans = secondInt;
Expand Down Expand Up @@ -382,21 +383,30 @@
isInplacePATref = str2double(twix_obj{s}.hdr.MeasYaps.sPat.ucRefScanMode(end))==2;
end
end

isPrescan = false;
try
if ( (twix_obj{s}.hdr.MeasYaps.sPreScanNormalizeFilter.ucOn == 1) && (s==1) )
isPrescan=true;
end
catch
%donothing
end
% declare data objects:
twix_obj{s}.image = twix_map_obj(arg,'image',filename,version,rstraj);
twix_obj{s}.noise = twix_map_obj(arg,'noise',filename,version);
twix_obj{s}.phasecor = twix_map_obj(arg,'phasecor',filename,version,rstraj);
twix_obj{s}.phasestab = twix_map_obj(arg,'phasestab',filename,version,rstraj);
twix_obj{s}.phasestabRef0 = twix_map_obj(arg,'phasestab_ref0',filename,version,rstraj);
twix_obj{s}.phasestabRef1 = twix_map_obj(arg,'phasestab_ref1',filename,version,rstraj);
twix_obj{s}.refscan = twix_map_obj(arg,'refscan',filename,version,rstraj);
twix_obj{s}.refscanPC = twix_map_obj(arg,'refscan_phasecor',filename,version,rstraj);
twix_obj{s}.refscanPS = twix_map_obj(arg,'refscan_phasestab',filename,version,rstraj);
twix_obj{s}.refscanPSRef0 = twix_map_obj(arg,'refscan_phasestab_ref0',filename,version,rstraj);
twix_obj{s}.refscanPSRef1 = twix_map_obj(arg,'refscan_phasestab_ref1',filename,version,rstraj);
twix_obj{s}.RTfeedback = twix_map_obj(arg,'rtfeedback',filename,version,rstraj);
twix_obj{s}.vop = twix_map_obj(arg,'vop',filename,version); % tx-array rf pulses
twix_obj{s}.image = twix_map_obj_lee(arg,'image',filename,version,rstraj);
twix_obj{s}.noise = twix_map_obj_lee(arg,'noise',filename,version);
twix_obj{s}.phasecor = twix_map_obj_lee(arg,'phasecor',filename,version,rstraj);
twix_obj{s}.phasestab = twix_map_obj_lee(arg,'phasestab',filename,version,rstraj);
twix_obj{s}.phasestabRef0 = twix_map_obj_lee(arg,'phasestab_ref0',filename,version,rstraj);
twix_obj{s}.phasestabRef1 = twix_map_obj_lee(arg,'phasestab_ref1',filename,version,rstraj);
twix_obj{s}.refscan = twix_map_obj_lee(arg,'refscan',filename,version,rstraj);
twix_obj{s}.refscanPC = twix_map_obj_lee(arg,'refscan_phasecor',filename,version,rstraj);
twix_obj{s}.refscanPS = twix_map_obj_lee(arg,'refscan_phasestab',filename,version,rstraj);
twix_obj{s}.refscanPSRef0 = twix_map_obj_lee(arg,'refscan_phasestab_ref0',filename,version,rstraj);
twix_obj{s}.refscanPSRef1 = twix_map_obj_lee(arg,'refscan_phasestab_ref1',filename,version,rstraj);
twix_obj{s}.RTfeedback = twix_map_obj_lee(arg,'rtfeedback',filename,version,rstraj);
twix_obj{s}.vop = twix_map_obj_lee(arg,'vop',filename,version); % tx-array rf pulses
twix_obj{s}.prescanbody = twix_map_obj_lee(arg,'prescanbody',filename,version,rstraj); %add prescan with body coil
twix_obj{s}.prescansurface = twix_map_obj_lee(arg,'prescansurface',filename,version,rstraj); %add prescan with surface coil

% print reader version information
if s==1
Expand Down Expand Up @@ -427,7 +437,7 @@
filePos( end ) = [];

% get mdhs and masks for each scan, no matter if noise, image, RTfeedback etc:
[mdh, mask] = evalMDH( mdh_blob, version, isInplacePATref); % this is quasi-instant (< 1s) :-)
[mdh, mask] = evalMDH( mdh_blob, version, isInplacePATref,isPrescan); % this is quasi-instant (< 1s) :-)
clear mdh_blob

% Assign mdhs to their respective scans and parse it in the correct twix objects.
Expand Down Expand Up @@ -540,6 +550,26 @@
end
twix_obj{s}.refscanPSRef1.readMDH( tmpMdh, filePos(isCurrScan) );
end
% prescan
if isPrescan && arg.bReadPrescan
if s==1
clear tmpMdh

isCurrScan = logical( mask.MDH_prescanbody );
for f = fieldnames( mdh ).'
tmpMdh.(f{1}) = mdh.(f{1})( isCurrScan, : );
end
twix_obj{s}.prescanbody.readMDH( tmpMdh, filePos(isCurrScan) );

clear tmpMdh
isCurrScan = logical( mask.MDH_prescansurface );
for f = fieldnames( mdh ).'
tmpMdh.(f{1}) = mdh.(f{1})( isCurrScan, : );
end
twix_obj{s}.prescansurface.readMDH( tmpMdh, filePos(isCurrScan) );
end
end

clear mdh tmpMdh filePos isCurrScan

% H.-L. Lee, extract PMU data from mdh and interpolate the waveforms to
Expand All @@ -555,7 +585,8 @@
for scan = { 'image', 'noise', 'phasecor', 'phasestab', ...
'phasestabRef0', 'phasestabRef1', 'refscan', ...
'refscanPC', 'refscanPS', 'refscanPSRef0', ...
'refscanPSRef1', 'RTfeedback', 'vop' }
'refscanPSRef1', 'RTfeedback', 'vop' ...
,'prescanbody','prescansurface'}
f = scan{1};

% remove unused fields
Expand Down Expand Up @@ -627,10 +658,14 @@
filePos = zeros(0, 1, class(cPos)); % avoid bug in Matlab 2013b: https://scivision.co/matlab-fseek-bug-with-uint64-offset/

% H.-L. Lee, setup for PMU data
if strncmp(VerString,'XA3',3) || strncmp(VerString,'XA5',3) || strncmp(VerString,'XA6',3)
if strncmp(VerString,'XA3',3)
syncdata_length = 1632;
elseif strncmp(VerString,'XA2',3)
syncdata_length = 1120;
elseif strncmp(VerString,'XA60',4)
syncdata_length = 960;
elseif strncmp(VerString,'XA61',4)
syncdata_length = 1204;%9632;
else
syncdata_length = 0;
end
Expand Down Expand Up @@ -783,7 +818,7 @@



function [mdh,mask] = evalMDH( mdh_blob, version, isInplacePATref )
function [mdh,mask] = evalMDH( mdh_blob, version, isInplacePATref, isPrescan )
% see pkg/MrServers/MrMeasSrv/SeqIF/MDH/mdh.h
% and pkg/MrServers/MrMeasSrv/SeqIF/MDH/MdhProxy.h

Expand Down Expand Up @@ -865,13 +900,23 @@
if ~isInplacePATref
noImaScan = (noImaScan | (mask.MDH_PATREFSCAN & ~mask.MDH_PATREFANDIMASCAN));
end

if isPrescan
mask.MDH_prescanbody = mdh.ushUsedChannels==2;

tempmask=false(Nmeas,1);
%sssr: does that mean we only use the 2nd body coil?
cind=find(mdh.ushUsedChannels==2);
tempmask(cind-1)=true;
mask.MDH_prescansurface = tempmask;
end

mask.MDH_IMASCAN( noImaScan ) = 0;

end % of evalMDH()


% H.-L. Lee, extract PMU data and timestamps and store them in the mdh struct
% H.-L. Lee, read PMU data and timestamps and store them in the mdh struct
function PMUdata = evalMDH_syncData( mdh_blob, version ,timestamps)

vec = @(x) x(:);
Expand Down Expand Up @@ -903,42 +948,58 @@
mdh.ulDMALength = data_uint32(1,:); % 1 : 4
mdh.ulTimeStamp = data_uint32(4,:).'; % 13 : 16

mdh.EKG.data = [vec(data_uint16(1:2:80,:)) vec(data_uint16(85:2:164,:)) vec(data_uint16(169:2:248,:)) vec(data_uint16(253:2:332,:))];
mdh.PULS.data = vec(data_uint16(337:2:376,:));
mdh.RESP.data = vec(data_uint16(381:2:390,:));
mdh.EXT.data = [vec(data_uint16(395:2:404,:)) vec(data_uint16(409:2:418,:))];

timeStamp40 = interp1(0:40:40*Nmeas-1,double(mdh.ulTimeStamp),1:40*Nmeas,'linear','extrap').';
timeStamp20 = interp1(0:20:20*Nmeas-1,double(mdh.ulTimeStamp),1:20*Nmeas,'linear','extrap').';
timeStamp05 = interp1(0:5:5*Nmeas-1, double(mdh.ulTimeStamp), 1:5*Nmeas,'linear','extrap').';

mdh.EKG.TimeStamp = timeStamp40;
mdh.PULS.TimeStamp = timeStamp20;
mdh.RESP.TimeStamp = timeStamp05;
mdh.EXT.TimeStamp = timeStamp05;
try
mdh.EVNT.data = vec(data_uint16(423:2:502,:));
mdh.EVNT.TimeStamp = timeStamp40;
catch
mdh.EVNT.data = zeros(40,1);
mdh.EVNT.TimeStamp = timeStamp40;
end
if(size(data_uint16,1)>342)

mdh.EKG.data = [vec(data_uint16(1:2:80,:)) vec(data_uint16(85:2:164,:)) vec(data_uint16(169:2:248,:)) vec(data_uint16(253:2:332,:))];
mdh.PULS.data = vec(data_uint16(337:2:376,:));
mdh.RESP.data = vec(data_uint16(381:2:390,:));
mdh.EXT.data = [vec(data_uint16(395:2:404,:)) vec(data_uint16(409:2:418,:))];

timeStamp40 = interp1(0:40:40*Nmeas-1,double(mdh.ulTimeStamp),1:40*Nmeas,'linear','extrap').';
timeStamp20 = interp1(0:20:20*Nmeas-1,double(mdh.ulTimeStamp),1:20*Nmeas,'linear','extrap').';
timeStamp05 = interp1(0:5:5*Nmeas-1, double(mdh.ulTimeStamp), 1:5*Nmeas,'linear','extrap').';

mdh.EKG.TimeStamp = timeStamp40;
mdh.PULS.TimeStamp = timeStamp20;
mdh.RESP.TimeStamp = timeStamp05;
mdh.EXT.TimeStamp = timeStamp05;
try
mdh.EVNT.data = vec(data_uint16(423:2:502,:));
mdh.EVNT.TimeStamp = timeStamp40;
catch
mdh.EVNT.data = zeros(40,1);
mdh.EVNT.TimeStamp = timeStamp40;
end

PMUdata.raw = mdh;
try
% sync PMU data to ADC timestamps
PMUdata.EKG = interp1(mdh.EKG.TimeStamp, double(mdh.EKG.data), double(timestamps),'linear','extrap').';
PMUdata.PULS = interp1(mdh.PULS.TimeStamp,double(mdh.PULS.data),double(timestamps),'linear','extrap');
PMUdata.RESP = interp1(mdh.RESP.TimeStamp,double(mdh.RESP.data),double(timestamps),'linear','extrap');
PMUdata.EXT = interp1(mdh.EXT.TimeStamp, double(mdh.EXT.data), double(timestamps),'linear','extrap').';
PMUdata.EVNT = interp1(mdh.EVNT.TimeStamp,double(mdh.EVNT.data),double(timestamps),'linear','extrap');
catch
PMUdata.EKG = zeros(4,numel(timestamps));
PMUdata.PULS = zeros(1,numel(timestamps));
PMUdata.RESP = zeros(1,numel(timestamps));
PMUdata.EXT = zeros(2,numel(timestamps));
PMUdata.EVNT = zeros(1,numel(timestamps));
fprintf('PMU data interpolation failed.\n');
end


else
mdh.EKG.data = [vec(data_uint16(1:2:80,:)) vec(data_uint16(85:2:164,:)) vec(data_uint16(169:2:248,:)) vec(data_uint16(253:2:332,:))];
mdh.EXT.data = vec(data_uint16(335:342,:));

timeStamp40 = interp1(0:40:40*Nmeas-1,double(mdh.ulTimeStamp),1:40*Nmeas,'linear','extrap').';
timeStamp08 = interp1(0:8:8*Nmeas-1, double(mdh.ulTimeStamp), 1:8*Nmeas,'linear','extrap').';

mdh.EKG.TimeStamp = timeStamp40;
mdh.EXT.TimeStamp = timeStamp08;

PMUdata.raw = mdh;
try
% sync PMU data to ADC timestamps
PMUdata.EKG = interp1(mdh.EKG.TimeStamp, double(mdh.EKG.data), double(timestamps),'linear','extrap').';
PMUdata.PULS = interp1(mdh.PULS.TimeStamp,double(mdh.PULS.data),double(timestamps),'linear','extrap');
PMUdata.RESP = interp1(mdh.RESP.TimeStamp,double(mdh.RESP.data),double(timestamps),'linear','extrap');
PMUdata.EXT = interp1(mdh.EXT.TimeStamp, double(mdh.EXT.data), double(timestamps),'linear','extrap').';
PMUdata.EVNT = interp1(mdh.EVNT.TimeStamp,double(mdh.EVNT.data),double(timestamps),'linear','extrap');
catch
PMUdata.EKG = zeros(4,numel(timestamps));
PMUdata.PULS = zeros(1,numel(timestamps));
PMUdata.RESP = zeros(1,numel(timestamps));
PMUdata.EXT = zeros(2,numel(timestamps));
PMUdata.EVNT = zeros(1,numel(timestamps));
fprintf('PMU data interpolation failed.\n');
end

end % of evalMDH_syncData()
22 changes: 14 additions & 8 deletions twix_map_obj.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -353,17 +353,23 @@
this.NCol = max(this.NCol);
this.NCha = max(this.NCha);

if strcmp(this.dataType,'refscan')
if strcmp(this.dataType,'image') || strcmp(this.dataType,'refscan')
%pehses: check for lines with 'negative' line/partition numbers
%this can happen when the reference scan line/partition range
%exceeds the one of the actual imaging scan
if this.NLin>65500 %uint overflow check
this.Lin = mod(this.Lin + (65536 - min(this.Lin(this.Lin>65500))),65536)+1;
this.NLin = max(this.Lin);
end
if this.NPar>65500 %uint overflow check
this.Par = mod(this.Par + (65536 - min(this.Par(this.Par>65500))),65536)+1;
this.NPar = max(this.Par);
Nacqdims = this.NPar*this.NSli*this.NAve*this.NPhs...
*this.NEco*this.NRep*this.NSet;

Noverflow= 1+floor((this.NAcq-(Nacqdims*65536))/65536);
if(Noverflow>0)
for of=1:(Noverflow)
this.Lin((65536*of*Nacqdims)+1:end) = this.Lin((65536*of*Nacqdims)+1:end) + 65536;
if(this.NPar==(65536*of+1))
this.Par((65536*of*Nacqdims)+1:end) = this.Par((65536*of*Nacqdims)+1:end) + 65536;
end
this.NLin = max(this.Lin);
this.NPar = max(this.Par);
end
end
end

Expand Down