1function write_v141(obj,filename,create_signature)
2%WRITE Write sequence to file.
3% WRITE(seqObj, filename) Write the sequence data to the given
4% filename using the open file format for MR sequences.
5%
6% Examples:
7% Write the sequence file to the my_sequences directory, relying on the
8% previous stable file format (v143).
9%
10% write_v143(seqObj,'my_sequences/gre.seq')
11%
12% See also write, read
14if (nargin<3)
15 create_signature=true;
16end
18fid=fopen(filename, 'w');
19assert(fid ~= -1, 'Cannot open file: %s', filename);
20fprintf(fid, '# Pulseq sequence file\n');
21fprintf(fid, '# Created by MATLAB mr toolbox\n\n');
23% we always write files in the default current version, which may be
24% differen to one, loaded (and stored in the seq object)
25[version_major, version_minor, version_revision]=mr.aux.version('output');
26fprintf(fid, '[VERSION]\n');
27fprintf(fid, 'major %s\n', num2str(1));
28fprintf(fid, 'minor %s\n', num2str(4));
29fprintf(fid, 'revision %s\n', num2str(1));
30fprintf(fid, '\n');
32if ~isempty(obj.definitions)
33 fprintf(fid, '[DEFINITIONS]\n');
34 keys = obj.definitions.keys;
35 values = obj.definitions.values;
36 for i=1:length(keys)
37 fprintf(fid, '%s ', keys{i});
38 if (ischar(values{i}))
39 fprintf(fid, '%s ', values{i});
40 else
41 fprintf(fid, '%.9g ', values{i});
42 end
43 fprintf(fid, '\n');
44 end
45 fprintf(fid, '\n');
46end
48fprintf(fid, '# Format of blocks:\n');
49fprintf(fid, '# NUM DUR RF GX GY GZ ADC EXT\n');
50fprintf(fid, '[BLOCKS]\n');
51idFormatWidth = length(num2str(length(obj.blockEvents)));
52idFormatStr = ['%' num2str(idFormatWidth) 'd'];
53for i = 1:length(obj.blockEvents)
54 %fprintf(fid,[idFormatStr ' %2d %2d %3d %3d %3d %2d 0\n'],[i obj.blockEvents(i,:)]);
55 %fprintf(fid,[idFormatStr ' %2d %2d %3d %3d %3d %2d 0\n'],[i obj.blockEvents{i}]);
56 bd=obj.blockDurations(i)/obj.blockDurationRaster;
57 bdr=round(bd);
58 assert(abs(bdr-bd)<1e-6); % this may still trigger false alarms for very long delays due to the limited accuracy of the double
59 fprintf(fid,[idFormatStr ' %3d %3d %3d %3d %3d %2d %2d\n'], ...
60 [i bdr obj.blockEvents{i}(2:end)]);
61end
62fprintf(fid, '\n');
64if ~isempty(obj.rfLibrary.keys)
65 fprintf(fid, '# Format of RF events:\n');
66 fprintf(fid, '# id amplitude mag_id phase_id time_shape_id delay freq phase\n');
67 fprintf(fid, '# .. Hz .... .... .... us Hz rad\n');
68 fprintf(fid, '[RF]\n');
69 keys = obj.rfLibrary.keys;
70 for k = keys
71 libData1 = obj.rfLibrary.data(k).array(1:4);
72 full_freqPhaseOffset=obj.rfLibrary.data(k).array(9:10)+obj.rfLibrary.data(k).array(7:8)*1e-6*obj.sys.gamma*obj.sys.B0; % convert freqPPM and phasePPM to fixed offsets
73 delay = round(obj.rfLibrary.data(k).array(6)/obj.rfRasterTime)*obj.rfRasterTime*1e6; % a bit of a hack: round the delay
74 fprintf(fid, '%d %12g %d %d %d %g %g %g\n', [k libData1 delay full_freqPhaseOffset]);
75 end
76 fprintf(fid, '\n');
77end
79arbGradMask = obj.gradLibrary.type == 'g';
80trapGradMask = obj.gradLibrary.type == 't';
82if any(arbGradMask)
83 fprintf(fid, '# Format of arbitrary gradients:\n');
84 fprintf(fid, '# time_shape_id of 0 means default timing (stepping with grad_raster starting at 1/2 of grad_raster)\n');
85 fprintf(fid, '# id amplitude amp_shape_id time_shape_id delay\n'); % do we need delay ???
86 fprintf(fid, '# .. Hz/m .. .. us\n');
87 fprintf(fid, '[GRADIENTS]\n');
88 keys = obj.gradLibrary.keys;
89 for k = keys(arbGradMask)
90 if obj.gradLibrary.data(k).array(5)<0
91 error("Current sequence contains arbitrary gradients with oversampling, which are not compatible with the format 1.4.1");
92 end
93 fprintf(fid, '%d %12g %d %d %d\n', ...
94 [k obj.gradLibrary.data(k).array([1 4 5]) ...
95 round(obj.gradLibrary.data(k).array(6)*1e6)]);
96 end
97 fprintf(fid, '\n');
98end
100if any(trapGradMask)
101 fprintf(fid, '# Format of trapezoid gradients:\n');
102 fprintf(fid, '# id amplitude rise flat fall delay\n');
103 fprintf(fid, '# .. Hz/m us us us us\n');
104 fprintf(fid, '[TRAP]\n');
105 keys = obj.gradLibrary.keys;
106 for k = keys(trapGradMask)
107 data = obj.gradLibrary.data(k).array;
108 data(2:end) = round(1e6*data(2:end));
109 fprintf(fid, '%2d %12g %3d %4d %3d %3d\n', [k data]);
110 end
111 fprintf(fid, '\n');
112end
114if ~isempty(obj.adcLibrary.keys)
115 fprintf(fid, '# Format of ADC events:\n');
116 fprintf(fid, '# id num dwell delay freq phase\n');
117 fprintf(fid, '# .. .. ns us Hz rad\n');
118 fprintf(fid, '[ADC]\n');
119 keys = obj.adcLibrary.keys;
120 for k = keys
121 data = [obj.adcLibrary.data(k).array(1:3) obj.adcLibrary.data(k).array(6:7)+obj.adcLibrary.data(k).array(4:5)*1e-6*obj.sys.gamma*obj.sys.B0].*[1 1e9 1e6 1 1]; % convert freqPPM and phasePPM to fixed offsets
122 fprintf(fid, '%d %d %.0f %.0f %g %g\n', [k data]);
123 end
124 fprintf(fid, '\n');
125end
127%if ~isempty(obj.delayLibrary.keys)
128% fprintf(fid, '# Format of delays:\n');
129% fprintf(fid, '# id delay (us)\n');
130% fprintf(fid, '[DELAYS]\n');
131% keys = obj.delayLibrary.keys;
132% for k = keys
133% fprintf(fid, '%d %d\n', ...
134% [k round(1e6*obj.delayLibrary.data(k).array)]);
135% end
136% fprintf(fid, '\n');
137%end
139if ~isempty(obj.extensionLibrary.keys)
140 fprintf(fid, '# Format of extension lists:\n');
141 fprintf(fid, '# id type ref next_id\n');
142 fprintf(fid, '# next_id of 0 terminates the list\n');
143 fprintf(fid, '# Extension list is followed by extension specifications\n');
144 fprintf(fid, '[EXTENSIONS]\n');
145 keys = obj.extensionLibrary.keys;
146 for k = keys
147 fprintf(fid, '%d %d %d %d\n', ...
148 [k round(obj.extensionLibrary.data(k).array)]);
149 end
150 fprintf(fid, '\n');
151end
153if ~isempty(obj.trigLibrary.keys)
154 fprintf(fid, '# Extension specification for digital output and input triggers:\n');
155 fprintf(fid, '# id type channel delay (us) duration (us)\n');
156% fprintf(fid, 'extension TRIGGERS 1\n'); % fixme: extension ID 1 is hardcoded here for triggers
157 fprintf(fid, ['extension TRIGGERS ',num2str(obj.getExtensionTypeID('TRIGGERS')),'\n']);
159 keys = obj.trigLibrary.keys;
160 for k = keys
161 fprintf(fid, '%d %d %d %d %d\n', ...
162 [k round(obj.trigLibrary.data(k).array.*[1 1 1e6 1e6])]);
163 end
164 fprintf(fid, '\n');
165end
167if ~isempty(obj.labelsetLibrary.keys) || ~isempty(obj.labelincLibrary.keys)
168 lbls=mr.getSupportedLabels();
170 if ~isempty(obj.labelsetLibrary.keys)
171 fprintf(fid, '# Extension specification for setting labels:\n');
172 fprintf(fid, '# id set labelstring\n');
173 tid=obj.getExtensionTypeID('LABELSET');
174 fprintf(fid, ['extension LABELSET ',num2str(tid),'\n']);
175 keys = obj.labelsetLibrary.keys;
176 for k = keys
177 fprintf(fid, '%d %d %s\n', ...
178 k, obj.labelsetLibrary.data(k).array(1),lbls{obj.labelsetLibrary.data(k).array(2)});
179 end
180 fprintf(fid, '\n');
181 end
182 if ~isempty(obj.labelincLibrary.keys)
183 fprintf(fid, '# Extension specification for increasing labels:\n');
184 fprintf(fid, '# id inc labelstring\n');
185 tid=obj.getExtensionTypeID('LABELINC');
186 fprintf(fid, ['extension LABELINC ',num2str(tid),'\n']);
187 lbls=mr.getSupportedLabels();
188 keys = obj.labelincLibrary.keys;
189 for k = keys
190 fprintf(fid, '%d %d %s\n', ...
191 k, obj.labelincLibrary.data(k).array(1),lbls{obj.labelincLibrary.data(k).array(2)});
192 end
193 fprintf(fid, '\n');
194 end
195end
197if ~isempty(obj.softDelayLibrary.keys)
198 warning('WARNING! The sequence in memory uses ''soft delay'' extension, which is incompatible with the file format v1.4.1. The produced Pulseq file is only partially valid and may fail to load or operate in some cases');
199end
201if ~isempty(obj.rotationLibrary.keys)
202 error('WARNING! The sequence in memory uses the ''rotations'' extension, which is incompatible with the file format v1.4.1. The produced Pulseq file is likely to be invalid and would probably fail to operate');
203end
205if ~isempty(obj.shapeLibrary.keys)
206 fprintf(fid, '# Sequence Shapes\n');
207 fprintf(fid, '[SHAPES]\n\n');
208 keys = obj.shapeLibrary.keys;
209 for k = keys
210 shape_dat = obj.shapeLibrary.data(k).array;
211 fprintf(fid, 'shape_id %d\n', k);
212 fprintf(fid, 'num_samples %d\n', shape_dat(1));
213 fprintf(fid, '%.9g\n', shape_dat(2:end));
214 fprintf(fid, '\n');
215 end
216end
218fclose(fid);
220if create_signature
221 % sign the file (this version with re/loading the file is a factor 2 faster than the sprintf() based one that kept a memory-copy of the data written)
223 % re-open and read in the file
224 fid=fopen(filename, 'r');
225 buf=fread(fid);
226 fclose(fid);
228 % calculate the digest
229 if mr.aux.isOctave()
230 md5hash=hash('MD5',char(buf(:)')); % Octave-specific function
231 else
232 md5hash=md5_java(buf); % Matlab Java hack
233 end
234 %fprintf('%s\n',md5hash);
236 % store the signature in the object
237 obj.signatureType='md5';
238 obj.signatureFile='text';
239 obj.signatureValue=md5hash;
241 % re-open the file for appending
242 fid=fopen(filename, 'a');
243 fprintf(fid, '\n[SIGNATURE]\n'); % the preceding new line BELONGS to the signature (and needs to be sripped away to recalculate the signature)
244 fprintf(fid, '# This is the hash of the Pulseq file, calculated right before the [SIGNATURE] section was added\n');
245 fprintf(fid, '# It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE]\n');
246 fprintf(fid, '# The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be sripped away for recalculating/verification)\n');
247 fprintf(fid, 'Type md5\n');
248 fprintf(fid, 'Hash %s\n', md5hash);
249 fclose(fid);
250end
252end
254function out=md5_java(buf)
255 import java.security.*;
256 import java.math.*;
257 import java.lang.String;
259 md = MessageDigest.getInstance('MD5');
260 hash = md.digest(double(buf));
261 bi = BigInteger(1, hash);
263 out=char(String.format('%032x', bi));
264end