Blame src/engine/pulseq/+mr/@Sequence/calcPNS.m at 78d04f1a76b841cf1ae39c4675d7366ef2670742

78d04f1 seqlab: write and view pulseq MRI sequences in the browserJeremy Magland 26 days ago1function [ok, pns_norm, pns_comp, t_axis]=calcPNS(obj,hardware,doPlots,calcCNS)

2% calculate PNS using safe model implementation by Szczepankiewicz and Witzel

3% assumes safe_pns_prediction package has been downloaded and installed in

4% Matlab path. See http://github.com/filip-szczepankiewicz/safe_pns_prediction

6% returns pns levels due to respective axes (normalized to 1 and not to 100%)

8% inputs:

9% hardware - hardware specifications. see safe_example_hw() from

10% the safe_pns_prediction package. Alternatively a text file

11% in the .asc format (Siemens) can be passed, e.g. for Prisma

12% it is MP_GPA_K2309_2250V_951A_AS82.asc (we leave it as an

13% exercise to the interested user to find were these files

14% can be acquired from);

15% doPlots - optional parameter (defaluts to true)

16% calcCNS - optional parameter (defaults to false): calculate cardiac

17% stimulation (CNS) instead of the PNS if the modela data

18% are provided by the given platform

20if nargin < 3

21 doPlots=true;

22 %calcCNS=false;

23end

24%if nargin < 4

25% calcCNS=false;

26%end

28% acquire the entire gradient wave form

29gw=obj.waveforms_and_times();

30if doPlots

31 figure;

32 plot(gw{1}(1,:),gw{1}(2,:),gw{2}(1,:),gw{2}(2,:),gw{3}(1,:),gw{3}(2,:)); % plot the entire gradient shape

33 title('gradient wave form, in T/m');

34end

36% find beginning and end times and resample GWs to a regular sampling raster

37tf=[];

38tl=[];

39for i=1:3

40 if size(gw{i},2)>0

41 tf(end+1)=gw{i}(1,1);

42 tl(end+1)=gw{i}(1,end);

43 end

44end

45nt_min=floor(min(tf)/obj.gradRasterTime+eps);

46nt_max=ceil(max(tl)/obj.gradRasterTime-eps);

47% shift raster positions to the centers of the raster periods

48nt_min = nt_min + 0.5;

49nt_max = nt_max - 0.5;

50if (nt_min<0.5)

51 nt_min=0.5

52end

53t_axis=(nt_min:nt_max)*obj.gradRasterTime;

54gwr=zeros(length(t_axis),3);

55for i=1:3

56 if size(gw{i},2)>0

57 gwr(:,i)=interp1(gw{i}(1,:),gw{i}(2,:),t_axis,'linear',0);

58 end

59end

61asc=[];

62ascHasCNS=false;

63if ischar(hardware)

64 % this loads the parameters from the provided text file

65 asc=mr.Siemens.readasc(hardware);

66 ascHasCNS=isfield(asc, 'GradPatSup') && isfield(asc.GradPatSup, 'Phys') && isfield(asc.GradPatSup.Phys, 'CarNS');

67end

69% for Cima.X and later we may want to calculate both CNS and PNS

70if ~exist('calcCNS','var')

71 if ascHasCNS

72 calcCNS=[false,true];

73 else

74 calcCNS=false;

75 end

76end

78pns_comp=[];

79pns_norm=[];

80ok=[];

81for c=calcCNS

82 hardware=asc_to_hw(asc,c);

84 % use the Szczepankiewicz' and Witzel's implementation

85 [pns_comp0,res]=safe_gwf_to_pns(gwr/obj.sys.gamma, NaN*ones(length(t_axis),1), obj.gradRasterTime, hardware); % the RF vector is unused in the code inside but it is zeropaded and exported ...

86 % use the exported RF vector to detect and undo zerpopadding

87 pns_comp=vertcat(pns_comp,0.01*pns_comp0(~isfinite(res.rf),:)');

88 % calc pns_norm and the final ok/not_ok

89 pns_norm=vertcat(pns_norm,vecnorm(pns_comp(end-2:end,:)));

90 ok=vertcat(ok,all(pns_norm(end,:)<1));

91 % ready

92 if doPlots

93 % plot results

94 figure;

95 h=safe_plot(pns_comp0, obj.gradRasterTime);

96 if c

97 p=[h.Parent];