~registry/ssacontrolled/saxescriptsaes : revision 33

1

function save = write_to_spm(subj,objtype,objin,varargin)

2

3

% Writes an MVPA data object or group to a SPM file.

4

%

5

% [] = WRITE_TO_SPM(SUBJ,OBJTYPE,OBJNAME,...)

6

%

7

%% Required Fields

8

%

9

% Writes an MVPA object in the subject structure SUBJ to an SPM file. This

10

% structure is to be of type OBJTYPE and must corespond to the specified

11

% OBJNAME. This OBJNAME can be either and object name, a group name or a

12

% cell arry of object names.

13

%

14

% SUBJ: The subject file that contains the data to be saved. This object

15

% is required.

16

%

17

% OBJTYPE: The type of object to be saved by the function. Currently it

18

% supports patterns and masks. This is simply a string indicating the

19

% type.

20

%

21

% OBJNAME: The name of the object to be saved. This is to be a string

22

% defining the name of the object.

23

%

24

%% Options:

25

%

26

% OUTPUT_FILENAME(string/*OBJNAME): Subverts OBJNAME as the string to be used when

27

% generating the output file name. It must be of the same type and form of

28

% the OBJNAME provided earlier. If this option is not used then the prefix

29

% for all output files will be the value of OBJNAME.

30

%

31

% PATHNAME(string/*''): This alters the distination of the final files to be

32

% created This path is relative to the current working directory. This

33

% directory must also already exsist.

34

%

35

% ONEMINUS(True/*False): If true, this will write out 1-M (where M is your

36

% matrix). This is useful if you're writing out an anova statmap of p

37

% values, where low is better. 1-p will reverse things so that higher is

38

% better when you write it out for convenience. Note: this has been tweaked

39

% so that zero values stay as zero

40

%

41

% NEW_HEADER(True/*False): This command option tells the script to generate

42

% a new header.

43

%

44

% VOXEL_DIMS([double double double]/*[2 2 2]): This function takes an array

45

% of 3 numbers that indicate the mm resolution of the voxels being saved.

46

% This value set can be retrieved using the function

47

% 'spm_extrac_voxel_dims' on a pre existing pattern

48

%

49

%

50

% PADDING(integer): This option is to be included if you wish to specify

51

% the number of numbers used in a naming scheme so if this is a five and

52

% you had 250 entries the filenames would include a number ranging from

53

% 00001 to 00250. (counting starts from zero.) This will however be

54

% ignored if the number is smaller than the maximum width. So if you

55

% specify a 2 with 300 entries, the numbers will display 001 - 300 not 1 -

56

% 300. This function will only have an effect if you are saving a group of

57

% files. Single patterns have no innate numeric association so this is not

58

% used in that context, files can still be tagged with arbitrary numbers by

59

% using the OUTPUT_FILENAME variable.

60

%

61

% FEXTENSION(*'.nii'/'.img')

62

% If unset, this will assume you wish to run the tutorial against nifti

63

% data. If set to '.img' it will change to using the analyze data set.

64

%

65

% Notes: This function currently requires the use of the SPM5 library it

66

% will save files in 4D format, so all time series of a given pattern are

67

% in a single file. Files in a group will each have a separate number

68

% starting at 1 and extending to the last object in the group.

69

%

70

% Other Potholes: There is a special case of patterns in which this function

71

% will not work (currently). This is the instance where multiple files

72

% worth of data are stores in one pattern area. This will happen if you

73

% load several files into a single pattern object (eg. the 'epi' pattern

74

% created when running tutorial_easy_spm). You can however save these

75

% patterns to disk but you will have to generate a new header and all of

76

% the data will be stored in a single spm file instead of multiple spm

77

% files.

78

79

80

%% License:

81

%=====================================================================

82

%

83

% This is part of the Princeton MVPA toolbox, released under

84

% the GPL. See http://www.csbmb.princeton.edu/mvpa for more

85

% information.

86

%

87

% The Princeton MVPA toolbox is available free and

88

% unsupported to those who might find it useful. We do not

89

% take any responsibility whatsoever for any problems that

90

% you have related to the use of the MVPA toolbox.

91

%

92

% ======================================================================

93

94

%% Test the basic validity of the object and setup some basics reqs.

95

% Test that the requested save is both a valid name and only a valid name

96

% for either a group or an object not both. Also it makes sure there is

97

% only one instance of the name.

98

[obj_name isgroup] = find_group_single(subj,objtype,objin);

99

100

101

102

%% Handle the setup of all defaults for the function.

103

% assume the output filenames can be based off the object names

104

% assume there is no change to the current path by default

105

defaults.output_filename = '';

106

defaults.pathname = '';

107

% Turn off Oneminus effect unless it's called for in optional arguements.

108

defaults.oneminus = false;

109

% set the initial file name padding to zero unless overwritten:

110

defaults.padding = 0;

111

%preload new headers option to false.

112

defaults.new_header = false;

113

%setup the default voxel dimensions

114

defaults.voxel_dims = [2 2 2];

115

defaults.fextension = '.nii';

116

117

118

% collect that other arguements made into a variable for analysis. and load

119

% them in as replacements for the defaults as applicable.

120

121

args = propval(varargin,defaults);

122

123

%populate the requisit values based on defaults or args.

124

%output_filename

125

if (isfield(args,'output_filename'))

126

output_filename=args.output_filename;

127

else

128

output_filename=defaults.output_filename;

129

end

130

131

%pathname

132

133

if (isfield(args,'pathname'))

134

pathname=args.pathname;

135

else

136

pathname=defaults.pathname;

137

end

138

%oneminus

139

if (isfield(args,'oneminus'))

140

oneminus=args.oneminus;

141

else

142

oneminus=defaults.oneminus;

143

end

144

%padding

145

if (isfield(args,'padding'))

146

padding=args.padding;

147

else

148

padding=defaults.padding;

149

end

150

%new_header

151

if (isfield(args,'new_header'))

152

new_header=args.new_header;

153

else

154

new_header=defaults.new_header;

155

end

156

%default voxel size

157

if (isfield(args,'voxel_dims'))

158

voxel_dims=args.voxel_dims;

159

else

160

voxel_dims=defaults.voxel_dims;

161

end

162

%setup the file extension to be used by the function.

163

if (isfield(args,'fextension'))

164

165

%if STRCMP(args.type,'analyze')

166

fextension=args.fextension;

167

168

else

169

170

fextension = defaults.fextension;

171

172

173

end

174

175

176

%% Start loop based on the number of objects to be processed.

177

% The code should be run for each object in the group, or for the single

178

% object. So cur_object is the single object or the current object in the group.

179

num_in_group=length(obj_name);

180

181

for cur_object_index=1:num_in_group

182

183

184

%% Path Handling

185

% generate the padded index string for this index

186

s_cur_object_index=num2str(cur_object_index);

187

188

if (padding > length(s_cur_object_index))

189

tmp(1:padding)='0';

190

191

tmp(end-(length(s_cur_object_index)-1):end)=s_cur_object_index(:);

192

s_cur_object_index=tmp;

193

194

end

195

196

197

% Generate the final path if there is one to generate at all.

198

final_path='';

199

%first generate it based on output_filename

200

if ~(strcmp(output_filename,''))

201

202

203

if (isgroup)

204

final_path= [output_filename '_' s_cur_object_index];

205

else

206

final_path = output_filename;

207

end

208

209

final_path=[final_path fextension]; %#ok<AGROW>

210

211

end

212

%if you have to add a path, overwrite the output_filename (since you

213

%don't have internal knowledge of whether it was done or not, this

214

%should be as fast as testing for final_path having changed).

215

if ~(strcmp(pathname,''))

216

217

if (strcmp(output_filename,''))

218

final_path = [pathname '/' objname];

219

else

220

final_path = [pathname '/' output_filename];

221

end

222

223

if (isgroup)

224

final_path= [final_path '_' s_cur_object_index]; %#ok<AGROW>

225

end

226

227

final_path = [final_path fextension]; %#ok<AGROW>

228

229

end

230

231

%% Capture the current object

232

% capture the object name for future use.

233

cur_obj_name = obj_name{cur_object_index};

234

% Load the data into a local object to be processed.

235

cur_pattern_matrix = get_mat(subj, objtype, cur_obj_name);

236

% find out how many specific voxels you will be saving.

237

cur_voxel_count = size(cur_pattern_matrix);

238

239

240

%% Adjust for One Minus

241

% Test for the oneminus flag and act on it now that you have the

242

% cur_pattern_matrix.

243

if (oneminus)

244

cur_pattern_matrix=1-cur_pattern_matrix;

245

end

246

247

248

% switch statement will handle this nicely, allowing for setup before

249

% saving each file. The object type will also control how the loop is

250

% executed and some early setup but otherwise the function is the same and

251

% will simply be called once the setup is done.

252

253

% Calculate the size of the pattern your working with.

254

255

256

%% Selector for either mask or pattern saving.

257

switch lower(objtype)

258

case {'mask'}

259

%% Switch/MASK: capture pattern

260

masked_pattern_matrix=cur_pattern_matrix; %fix so that the mask and the pattern saves are inline with each other.

261

% if you are working with a mask the 'mask' of the item being written is

262

% simply an array of ones the size of the data set you are writing. This

263

% is due to the fact that we want to save a mask which means saving all of

264

% it's ones and zero's not just it's non zero values (which are the

265

% relevant values in other patterns). The dimensions of the mask are

266

% simply the size of the pattern_matrix (mask) being saved. And the

267

% pattern projection is simply the pattern matrix. The loop must be based

268

% on the number of different data sets that must be saved, these data sets

269

% are 3D so the 4th value in the masks voxel count is in fact the 'time

270

% index' which realistically should be one for all masks. Unless you've

271

% somehow defined a dynamic mask for relevant data.

272

%% Switch/MASK: Adjust pattern to be treated as 4D

273

pattern_size=size(cur_pattern_matrix);

274

%if (length(pattern_size) ==3)

275

if isequal(length(pattern_size),3)

276

pattern_size(4)=1;

277

end

278

279

% capture the volume information from the previous spm load command. (this

280

% should be modified to be more dynamic by say calling a function to parse

281

% this information out of a centrally stored set of data. This would allow

282

% for on the fly format changes.)

283

%% Switch/MASK: Setup volume

284

% as a cell structure, generate new header if needed, also

285

% set filenames if required

286

cur_vol=cell(1);

287

if (~new_header)

288

cur_vol{:} = get_objsubfield(subj, objtype, cur_obj_name,'header','vol');

289

290

if ~(strcmp(final_path,''))

291

for fname_index=1:pattern_size(4)

292

cur_vol{1}(fname_index).fname=final_path;

293

if (isfield(cur_vol,'private'))

294

cur_vol.private.dat.fname=final_path;

295

end

296

end

297

298

end

299

else

300

%if you ask for a new header then it will be generated here instead of being retrieved from the obj sub field.

301

% Current Work.

302

303

if (strcmp(final_path,''))

304

if (isgroup)

305

cur_vol=gen_vol_info(pattern_size(1:3),pattern_size(4),[cur_obj_name '_' s_cur_object_index fextension],voxel_dims);

306

else

307

cur_vol=gen_vol_info(pattern_size(1:3),pattern_size(4),[cur_obj_name fextension],voxel_dims);

308

end

309

else

310

cur_vol=gen_vol_info(pattern_size(1:3),pattern_size(4),final_path,voxel_dims);

311

end

312

end

313

314

315

%% Switch/MASK: Legacy cell insurance

316

% This being matlab a variety of the objects must be cell arrays for

317

% things to behave properly. As such some corrections must be made to the

318

% volume that was extracted if it is not a cell array. Rebuild it as a

319

% cell array to fix this problem.

320

if ~iscell(cur_vol)

321

temp = cur_vol;

322

cur_vol = cell(1);

323

cur_vol{1} = temp;

324

clear temp;

325

end

326

%% Switch/MASK: Loop to write out the file, masks should only loop once.

327

for cur_vol_index=1:pattern_size(4)%i need the 1xn value, so this gets it.

328

329

%This 'magic number' 1 in the cell call to vol is a necessary

330

%evil. Vol is populated as a 1x1 cell, this 1x1 cell contains

331

%an array of structures that control the actual layout of each

332

%section of the volume.

333

save = spm_write_vol(cur_vol{1}(cur_vol_index), masked_pattern_matrix(:,:,:,cur_vol_index));

334

end

335

336

337

338

339

case {'pattern'}

340

341

%% Switch/PATTERN:

342

343

344

345

% first. Establish the mask that you will be working with for a given

346

% pattern matrix.

347

%% Switch/PATTERN: Capture data and mask.

348

masked_by = get_objfield(subj,'pattern',cur_obj_name,'masked_by');

349

mask_matrix = get_mat(subj,'mask',masked_by);

350

% second. Calculate the size of the mask so you know what's going to be

351

% saved.

352

mask_dims = size(mask_matrix);

353

% Masked pattern projection creation: basically you take a mask and a

354

% pattern. For each time point in the pattern you apply it to the next 1

355

% in the mask you have associated with it. This builds a 3D image of the

356

% particular time relevant pattern you are working with for saving.

357

masked_pattern_matrix=zeros(mask_dims(1),mask_dims(2),mask_dims(3),cur_voxel_count(2));

358

% Generate a map of the relative mask information, this is

359

% basically an array of the relevant indexes in 1D space

360

% instead of 3D space.

361

relative_map_of_mask = find(mask_matrix);

362

% calculate the mask volume.(LxWxH)

363

mask_size=mask_dims(1)*mask_dims(2)*mask_dims(3);

364

% Populate the masked_pattern_matrix based on the relative

365

% maps list of indexes offset by so that you populate the

366

% mask in 4D.

367

%% Switch/PATTERN: Set to 4D if it is not already

368

pattern_size=size(masked_pattern_matrix);

369

%if (length(pattern_size) ==3)

370

if isequal(length(pattern_size),3)

371

pattern_size(4)=1;

372

end

373

374

375

%% Switch/PATTERN:Setup volume, header and filename info

376

% capture the volume information from the previous spm load command. (this

377

% should be modified to be more dynamic by say calling a function to parse

378

% this information out of a centrally stored set of data. This would allow

379

% for on the fly format changes.)

380

if (~new_header)

381

cur_vol = get_objsubfield(subj, objtype, cur_obj_name,'header','vol');

382

383

if ~(strcmp(final_path,''))

384

for fname_index=1:pattern_size(4)

385

cur_vol{1}(fname_index).fname=final_path;

386

if (isfield(cur_vol,'private'))

387

cur_vol.private.dat.fname=final_path;

388

end

389

end

390

391

end

392

else

393

394

if (strcmp(final_path,''))

395

if (isgroup)

396

cur_vol=gen_vol_info(pattern_size(1:3),pattern_size(4),[cur_obj_name '_' s_cur_object_index fextension],voxel_dims);

397

else

398

cur_vol=gen_vol_info(pattern_size(1:3),pattern_size(4),[cur_obj_name fextension],voxel_dims);

399

end

400

else

401

cur_vol=gen_vol_info(pattern_size(1:3),pattern_size(4),final_path,voxel_dims);

402

end

403

end

404

405

406

%% Switch/PATTERN: Legacy cell insurance. (may be able to deprecate)

407

% This being matlab a variety of the objects must be cell arrays for

408

% things to behave properly. As such some corrections must be made to the

409

% volume that was extracted if it is not a cell array. Rebuild it as a

410

% cell array to fix this problem.

411

if ~iscell(cur_vol)

412

temp = cur_vol;

413

cur_vol = cell(1);

414

cur_vol{1} = temp;

415

clear temp;

416

end

417

418

%% Switch/PATTERN: Load 4D pattern using mask

419

for index=1:cur_voxel_count(2)

420

masked_pattern_matrix(relative_map_of_mask(:)+((index-1)*mask_size))=cur_pattern_matrix(:,index);

421

end

422

423

%% Switch/PATTERN: Loop to save file in 4D

424

for cur_vol_index=1:pattern_size(4)%i need the 1xn value, so this gets it.

425

426

%This 'magic number' 1 in the cell call to vol is a necessary

427

%evil. Vol is populated as a 1x1 cell, this 1x1 cell contains

428

%an array of structures that control the actual layout of each

429

%section of the volume.

430

save = spm_write_vol(cur_vol{1}(cur_vol_index), masked_pattern_matrix(:,:,:,cur_vol_index));

431

end

432

433

end

434

end

435

436

end

437

438

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

439

%% Sub Function for Making fake headers.

440

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

441

442

function [header] = gen_vol_info(dims, time_slices, file_name,voxel_dims)

443

444

%% code recommended by Thomas Nichols, SPM mailing list

445

446

unmasked_size = dims(1)*dims(2)*dims(3);

447

448

% blind default for sizing of voxels as a 'happy medium' value.

449

Origin = dims/2;

450

mat = diag([voxel_dims 1]);

451

mat(1:3,4)= -Origin.*voxel_dims;

452

%% end recommended code.

453

454

455

456

457

458

for index=1:time_slices

459

header(index).dt = [spm_type('float64') spm_platform('bigend')]; %#ok<AGROW>

460

header(index).dim = dims; %#ok<AGROW>

461

header(index).mat=mat; %#ok<AGROW>

462

header(index).n=[index 1]; %#ok<AGROW> %guarantees that the indexing is correct

463

header(index).descrip = 'Generated by MVPA write_to analyze script'; %#ok<AGROW>

464

header(index).fname = file_name; %#ok<AGROW>

465

header(index).pinfo = [1; 0; (unmasked_size*(index-1))]; %#ok<AGROW>

466

467

end

468

469

470

471

end

472

473

474

475

476

477

478

479

480

481

482

483

484

485

486

487

488

489