src/engine/pulseq/+mr/@Sequence/readBinary.m at main

concept-collection / seqlab

321 lines · 13.2 KBBlame History Raw

1function readBinary(obj,filename)

2%READBINARY Load sequence from binary file.

3% READBINARY(seqObj, filename) Read the given filename and load sequence

4% data stored in the binary version of the Pulseq open file format. The

5% binary format is described in the specficiation available at

6% http://pulseq.github.io

8% Examples:

9% Load the sequence defined in gre.bin in sequences directory

10%

11% readBinary(seqObj,'sequences/gre.bin')

12%

13% See also writeBinary

15binaryCodes = obj.getBinaryCodes();

16fid=fopen(filename);

17magicNum = fread(fid,1,'int64=>int64');

18assert(magicNum==binaryCodes.fileHeader,'Not a Pulseq binary file');

19version_major = fread(fid,1,'int64');

20version_minor = fread(fid,1,'int64');

21version_revision = fread(fid,1,'int64');

22assert(version_major==obj.version_major,'Unsupported version_major %d', version_major) % too strict - no reverse compatibility yet

23assert(version_minor==obj.version_minor,'Unsupported version_minor %d', version_minor) % too strict - no reverse compatibility yet

24assert(version_revision==obj.version_revision,'Unsupported version_revision %d', version_revision) % too strict - no reverse compatibility yet

26% set version

27obj.version_major = version_major;

28obj.version_minor = version_minor;

29obj.version_revision = version_revision;

31% Clear sequence data

32obj.blockEvents={};

33obj.blockDurations=[];

34obj.definitions=containers.Map();

35obj.gradLibrary=mr.EventLibrary();

36obj.shapeLibrary=mr.EventLibrary();

37obj.rfLibrary=mr.EventLibrary();

38obj.adcLibrary=mr.EventLibrary();

39obj.trigLibrary=mr.EventLibrary();

40obj.labelsetLibrary=mr.EventLibrary();

41obj.labelincLibrary=mr.EventLibrary();

42obj.extensionLibrary=mr.EventLibrary();

43obj.rfShimLibrary=mr.EventLibrary();

44obj.softDelayLibrary=mr.EventLibrary();

45obj.softDelayHints1=containers.Map();

46obj.softDelayHints2={};

47obj.rotationLibrary=mr.EventLibrary();

48obj.extensionStringIDs={};

49obj.extensionNumericIDs=[];

50obj.signatureType='';

51obj.signatureFile='';

52obj.signatureValue='';

54% Load data from file

55while true

56 section = int64(fread(fid,1,'int64'));

57 if isempty(section)

58 break

59 end

61 switch section

62 case binaryCodes.section.definitions

63 obj.definitions = readDefinitions(fid);

64 v=obj.getDefinition('GradientRasterTime');

65 if ~isempty(v), obj.gradRasterTime=v; end

66 v=obj.getDefinition('RadiofrequencyRasterTime');

67 if ~isempty(v), obj.rfRasterTime=v; end

68 v=obj.getDefinition('AdcRasterTime');

69 if ~isempty(v), obj.adcRasterTime=v; end

70 v=obj.getDefinition('BlockDurationRaster');

71 if ~isempty(v), obj.blockDurationRaster=v; end

73 case binaryCodes.section.blocks

74 [obj.blockEvents, obj.blockDurations] = readBlocks(fid, obj.blockDurationRaster);

76 case binaryCodes.section.rf

77 % array: [amp mag_id phase_id time_shape_id center delay freqPPM phasePPM freq phase]

78 % type (use) stored separately as a char

79 numEvents = double(fread(fid,1,'int64'));

80 for i=1:numEvents

81 id = double(fread(fid,1,'int32'));

82 amp = double(fread(fid,1,'float64'));

83 ids = double(fread(fid,3,'int32'))'; % mag_id, phase_id, time_shape_id

84 ctr = double(fread(fid,1,'int64')) * 1e-12; % center (ps -> s)

85 dly = double(fread(fid,1,'int64')) * 1e-12; % delay (ps -> s)

86 fpp = double(fread(fid,4,'float64'))'; % freqPPM, phasePPM, freq, phase

87 use = char(fread(fid,1,'char'));

88 obj.rfLibrary.insert(id, [amp ids(1) ids(2) ids(3) ctr dly fpp(1) fpp(2) fpp(3) fpp(4)], use);

89 end

91 case binaryCodes.section.gradients

92 % array: [amp first last amp_shape_id time_shape_id delay]

93 numEvents = double(fread(fid,1,'int64'));

94 for i=1:numEvents

95 id = double(fread(fid,1,'int32'));

96 afl = double(fread(fid,3,'float64'))'; % amp, first, last

97 ids2 = double(fread(fid,2,'int32'))'; % amp_shape_id, time_shape_id

98 dly = double(fread(fid,1,'int64')) * 1e-12; % delay (ps -> s)

99 obj.gradLibrary.insert(id, [afl ids2 dly], 'g');

100 end

101

102 case binaryCodes.section.trapezoids

103 % array: [amp rise flat fall delay]

104 numEvents = double(fread(fid,1,'int64'));

105 for i=1:numEvents

106 id = double(fread(fid,1,'int32'));

107 amp = double(fread(fid,1,'float64'));

108 t4 = double(fread(fid,4,'int64'))' * 1e-12; % rise,flat,fall,delay (ps->s)

109 obj.gradLibrary.insert(id, [amp t4], 't');

110 end

111

112 case binaryCodes.section.adc

113 % array: [num dwell delay freqPPM phasePPM freq phase phase_id]

114 numEvents = double(fread(fid,1,'int64'));

115 for i=1:numEvents

116 id = double(fread(fid,1,'int32'));

117 num = double(fread(fid,1,'int64'));

118 dwell = double(fread(fid,1,'int64')) * 1e-12; % ps -> s

119 delay = double(fread(fid,1,'int64')) * 1e-12; % ps -> s

120 f4 = double(fread(fid,4,'float64'))'; % freqPPM, phasePPM, freq, phase

121 phid = double(fread(fid,1,'int32'));

122 obj.adcLibrary.insert(id, [num dwell delay f4(1) f4(2) f4(3) f4(4) phid]);

123 end

124

125 case binaryCodes.section.delays

126 readLegacyDelays(fid);

127

128 case binaryCodes.section.shapes

129 obj.shapeLibrary = readShapes(fid);

130

131 case binaryCodes.section.extensions

132 % array per entry: [type ref next_id]

133 numEvents = double(fread(fid,1,'int64'));

134 for i=1:numEvents

135 id = double(fread(fid,1,'int32'));

136 data = double(fread(fid,3,'int32'))';

137 obj.extensionLibrary.insert(id, data);

138 end

139

140 case binaryCodes.section.triggers

141 % type id(i32) then events: id(i32) type(i32) channel(i32) delay(i32,us) duration(i32,us)

142 ext_id = double(fread(fid,1,'int32'));

143 obj.setExtensionStringAndID('TRIGGERS', ext_id);

144 numEvents = double(fread(fid,1,'int64'));

145 for i=1:numEvents

146 id = double(fread(fid,1,'int32'));

147 tc = double(fread(fid,2,'int32'))'; % type, channel

148 dd = double(fread(fid,2,'int64'))' * 1e-12; % delay, duration (ps->s)

149 obj.trigLibrary.insert(id, [tc dd]);

150 end

151

152 case binaryCodes.section.labelset

153 ext_id = double(fread(fid,1,'int32'));

154 obj.setExtensionStringAndID('LABELSET', ext_id);

155 numEvents = double(fread(fid,1,'int64'));

156 for i=1:numEvents

157 id = double(fread(fid,1,'int32'));

158 data = double(fread(fid,2,'int32'))'; % value, label_index

159 obj.labelsetLibrary.insert(id, data);

160 end

161

162 case binaryCodes.section.labelinc

163 ext_id = double(fread(fid,1,'int32'));

164 obj.setExtensionStringAndID('LABELINC', ext_id);

165 numEvents = double(fread(fid,1,'int64'));

166 for i=1:numEvents

167 id = double(fread(fid,1,'int32'));

168 data = double(fread(fid,2,'int32'))'; % value, label_index

169 obj.labelincLibrary.insert(id, data);

170 end

171

172 case binaryCodes.section.softdelays

173 ext_id = double(fread(fid,1,'int32'));

174 obj.setExtensionStringAndID('DELAYS', ext_id);

175 numEvents = double(fread(fid,1,'int64'));

176 for i=1:numEvents

177 id = double(fread(fid,1,'int32'));

178 num = double(fread(fid,1,'int32'));

179 offset = double(fread(fid,1,'int64')) * 1e-12; % ps -> s

180 factor = double(fread(fid,1,'float64'));

181 hlen = double(fread(fid,1,'int32'));

182 hint = char(fread(fid,hlen,'char')');

183 % register hint string and get its index

184 if obj.softDelayHints1.isKey(hint)

185 hint_idx = obj.softDelayHints1(hint);

186 else

187 hint_idx = length(obj.softDelayHints2) + 1;

188 obj.softDelayHints1(hint) = hint_idx;

189 obj.softDelayHints2{hint_idx} = hint;

190 end

191 obj.softDelayLibrary.insert(id, [num offset factor hint_idx]);

192 end

193

194 case binaryCodes.section.rfshims

195 ext_id = double(fread(fid,1,'int32'));

196 obj.setExtensionStringAndID('RF_SHIMS', ext_id);

197 numEvents = double(fread(fid,1,'int64'));

198 for i=1:numEvents

199 id = double(fread(fid,1,'int32'));

200 num_chan = double(fread(fid,1,'int32'));

201 chan_data = double(fread(fid,2*num_chan,'float64'))';

202 obj.rfShimLibrary.insert(id, chan_data);

203 end

204

205 case binaryCodes.section.rotations

206 ext_id = double(fread(fid,1,'int32'));

207 obj.setExtensionStringAndID('ROTATIONS', ext_id);

208 numEvents = double(fread(fid,1,'int64'));

209 for i=1:numEvents

210 id = double(fread(fid,1,'int32'));

211 quat = double(fread(fid,4,'float64'))';

212 obj.rotationLibrary.insert(id, mr.aux.quat.normalize(quat));

213 end

214

215 case binaryCodes.section.signature

216 type_len = double(fread(fid,1,'int32'));

217 sig_type = char(fread(fid,type_len,'char')');

218 hash_len = double(fread(fid,1,'int32'));

219 hash_raw = uint8(fread(fid,hash_len,'uint8'));

220 fread(fid,1,'int64'); % original file length prior to signature append

221

222 obj.signatureType = sig_type;

223 obj.signatureFile = 'bin';

224 if isempty(hash_raw)

225 obj.signatureValue = '';

226 else

227 obj.signatureValue = lower(reshape(dec2hex(hash_raw,2)',1,[]));

228 end

229

230 otherwise

231 error('Unknown section code: %s',dec2hex(section));

232 end

233end

234fclose(fid);

235

236

237return

238

239%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

240%%%%%%%%%%%%%%%%%%%%%%% Helper functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

241

242 function def = readDefinitions(fid)

243 %readDefinitions Read the [DEFINITIONS] section of a sequence file.

244 % defs=readDefinitions(fid) Read user definitions from file

245 % identifier of an open MR sequence file and return a map of

246 % key/value entries.

247

248 def = containers.Map();

249 numDefs = fread(fid,1,'int64');

250 for iDef=1:numDefs

251 count = fread(fid,1,'int32');

252 key = char(fread(fid,count,'char'));

253 if size(key,1)~=1

254 key=key.';

255 end

256 count = fread(fid,1,'int32');

257 type = char(fread(fid,1,'char'));

258 switch type

259 case 'f'

260 values = double(fread(fid,count,'float64'));

261 case 'i'

262 values = int32(fread(fid,count,'int32'));

263 case 'c'

264 values = char(fread(fid,count,'char')').';

265 if size(values,1)~=1

266 values=values.';

267 end

268 if ~isempty(values) && values(end)==0

269 values=values(1:(end-1));

270 end

271 otherwise

272 error('Unknown definition type: %s', type);

273 end

274 def(key) = values;

275 end

276 end

277

278 function [eventTable, blockDurations] = readBlocks(fid, blockDurationRaster)

279 %readBlocks Read the [BLOCKS] section of a binary sequence file.

280 % Each block stores: duration(int64) + 6 event IDs (int32):

281 % rf, gx, gy, gz, adc, ext

282 % Returns eventTable (cell array) with 7 elements per block

283 % [0 rf gx gy gz adc ext] and blockDurations in seconds.

284

285 numBlocks = double(fread(fid,1,'int64'));

286 eventTable = cell(1, numBlocks);

287 blockDurations = zeros(1, numBlocks);

288 for ii = 1:numBlocks

289 dur_raster = double(fread(fid,1,'int64'));

290 event_ids = double(fread(fid,6,'int32'))';

291 blockDurations(ii) = dur_raster * blockDurationRaster;

292 eventTable{ii} = [0 event_ids]; % prepend 0 for legacy delay placeholder

293 end

294 end

295

296 function shapeLibrary = readShapes(fid)

297 %readShapes Read the [SHAPES] section of a binary sequence file.

298

299 shapeLibrary=mr.EventLibrary();

300 numShapes = fread(fid,1,'int64');

301 for iShape=1:numShapes

302 id = double(fread(fid,1,'int32'));

303 numUncompressed = double(fread(fid,1,'int64'));

304 numCompressed = double(fread(fid,1,'int64'));

305 data = double(fread(fid,numCompressed,'float32'))';

306 shapeData = [numUncompressed data];

307 shapeLibrary.insert(id,shapeData);

308 end

309 end

310

311 function readLegacyDelays(fid)

312 % readLegacyDelays Read and ignore legacy [DELAYS] binary section.

313 % Delay events are no longer represented via obj.delayLibrary.

314 numEvents = double(fread(fid,1,'int64'));

315 for ii=1:numEvents

316 fread(fid,1,'int32'); % id

317 fread(fid,1,'int64'); % delay value

318 end

319 end

320

321end