1 #include "mrilib.h"
2 #include "thd_niftiwrite.h"
3
4 /* prototypes */
5
6 nifti_image * populate_nifti_image(THD_3dim_dataset *dset, niftiwr_opts_t options) ;
7
8 void nifti_set_afni_extension(THD_3dim_dataset *dset,nifti_image *nim) ;
9
10 static int get_slice_timing_pattern( float * times, int len, float * delta );
11 static int needs_conversion_to_float(THD_3dim_dataset *dset, int warn);
12
13 /*******************************************************************/
14 /*! Write an AFNI dataset as a NIfTI file.
15 - dset = AFNI dataset
16 - options = structure with options to control the output
17
18 Return value is 1 if went OK, 0 if not.
19 -------------------------------------------------------------------*/
20
21 // [PT,DRG: Mar 8, 2019] Updated to include [qs]form_code = 5
22
THD_write_nifti(THD_3dim_dataset * din,niftiwr_opts_t options)23 int THD_write_nifti( THD_3dim_dataset *din, niftiwr_opts_t options )
24 {
25 THD_3dim_dataset * dset=NULL;
26 nifti_image *nim ;
27 nifti_brick_list nbl ;
28 int ii ;
29 char *fname , *cpt ;
30
31 ENTRY("THD_write_nifti") ;
32
33 nifti_set_debug_level(options.debug_level) ;
34
35 /*-- check inputs for goodness --*/
36
37 fname = nifti_strdup(options.infile_name);
38
39 if( !THD_filename_ok(fname) || fname[0] == '-' ){
40 ERROR_message("Illegal filename for NIfTI output: %s\n",
41 (fname != NULL) ? fname : "(null)" ) ;
42 RETURN(0) ;
43 }
44
45 /* if we need a float dataset, make one (insted of failing)
46 * (wasteful, but simple and effective) 6 Sep 2012 [rickr] */
47 if( needs_conversion_to_float(din, 1) ) {
48 dset = EDIT_full_copy(din, NULL);
49 EDIT_floatize_dataset(dset);
50 tross_Copy_History(din, dset); /* ZSS May 23 2013 */
51
52 if( ! ISVALID_DSET(dset) ) {
53 ERROR_message("failed to copy dset for NIfTI write\n");
54 RETURN(0);
55 }
56 } else dset = din;
57
58 if( !ISVALID_DSET(dset) ){
59 ERROR_message("Illegal input dataset for NIfTI output: %s\n", fname) ;
60 RETURN(0) ;
61 }
62
63 /*-- load dataset from disk, if need be --*/
64
65 DSET_load(dset) ;
66 if( !DSET_LOADED(dset) ){
67 ERROR_message(
68 "Can't write NIfTI file since dataset isn't loaded: %s\n", fname) ;
69 RETURN(0) ;
70 }
71
72 nim = populate_nifti_image(dset,options) ;
73 if( !nim ) RETURN(0) ; /* catch failure 6 Apr 2006 [rickr] */
74
75 /*-- construct filename --*/
76
77 nim->fname = malloc( strlen(fname)+16 ) ;
78 nim->iname = malloc( strlen(fname)+16 ) ;
79 strcpy(nim->fname,fname) ;
80 strcpy(nim->iname,fname) ;
81
82 /* 11 Oct 2005: Allow for .hdr/.img file outputs -- RWCox */
83
84 cpt = nifti_find_file_extension( nim->iname ) ;
85 if( cpt != NULL && strcmp(cpt,".hdr") == 0 ){
86 nim->nifti_type = 2 ; /* indicate 2 file output */
87 memcpy(cpt,".img",4) ; /* convert .hdr name to .img */
88 }
89
90 /*-- construct nifti_brick_list of pointers to data briks */
91
92 if( options.debug_level > 2 ) nifti_image_infodump(nim) ;
93 nbl.bricks = (void **) malloc ( DSET_NVALS(dset) * sizeof(void*) ) ;
94 nbl.nbricks = DSET_NVALS(dset) ;
95 nbl.bsize = DSET_BRICK_BYTES(dset,0) ;
96 for (ii = 0 ; ii < DSET_NVALS(dset) ; ii++ )
97 nbl.bricks[ii] = DSET_ARRAY(dset,ii) ;
98
99 /*-- 13 Mar 2006: check disk space --*/
100
101 { FILE *fp = fopen(nim->fname,"ab") ;
102 int mm = THD_freemegabytes(nim->fname) ;
103 int rr = (int)(dset->dblk->total_bytes/(1024*1024)) ;
104 if( fp != NULL ) fclose(fp) ;
105 if( mm >= 0 && mm <= rr )
106 WARNING_message("Disk space: writing dataset %s (%d MB),"
107 " but only %d free MB on disk" ,
108 nim->fname , rr , mm ) ;
109 }
110
111 /*-- use handy-dandy library function to write out data */
112
113 nifti_set_afni_extension( dset , nim ) ; /* 09 May 2005 - RWCox */
114
115 nifti_image_write_bricks (nim, &nbl ) ;
116
117 /* if we made a float copy, nuke it */
118 if( dset != din ) THD_delete_3dim_dataset(dset, True) ;
119
120 free(fname); /* free copied name */
121
122 RETURN(1) ;
123 }
124
125 /* if the dataset has inconsistent types or scale factors, then it needs
126 * to be converted fo floats in order to write 6 Sep 2012 [rickr] */
needs_conversion_to_float(THD_3dim_dataset * dset,int warn)127 static int needs_conversion_to_float(THD_3dim_dataset *dset, int warn)
128 {
129 float fac0; /* was int, ick! 11 Sep 2015 [rickr] */
130 int type0, ii;
131
132 ENTRY("needs_conversion_to_float");
133
134 if( ! ISVALID_DSET(dset) ) RETURN(0);
135
136 if( dset->dblk->nvals <= 1 ) RETURN(0); /* nothing can vary */
137
138 type0 = DSET_BRICK_TYPE(dset,0) ;
139 fac0 = DSET_BRICK_FACTOR(dset,0) ;
140
141 for( ii=1 ; ii < DSET_NVALS(dset) ; ii++ ){
142 if( DSET_BRICK_TYPE(dset,ii) != type0) {
143 if( warn )
144 WARNING_message("varying brick types, writing NIfTI as float");
145 RETURN(1);
146 }
147
148 if( DSET_BRICK_FACTOR(dset,ii) != fac0) {
149 if( warn )
150 WARNING_message("varying brick factors, writing NIfTI as float");
151 RETURN(1);
152 }
153 }
154
155 RETURN(0);
156 }
157
158
159 /*******************************************************************/
160
populate_nifti_image(THD_3dim_dataset * dset,niftiwr_opts_t options)161 nifti_image * populate_nifti_image(THD_3dim_dataset *dset, niftiwr_opts_t options)
162 {
163 int nparam, type0 , ii , jj;
164 int nif_x_axnum=0, nif_y_axnum=0, nif_z_axnum=0;
165 int slast, sfirst ;
166 int pattern, tlen ;
167 nifti_image *nim ;
168 char axcode[3], axsign[3] ;
169 float axstep[3] , axstart[3] ;
170 int axnum[3] ;
171 float fac0 ;
172 double dumqx, dumqy, dumqz, dumdx, dumdy, dumdz ;
173 float *tlist, fsdur;
174
175 ENTRY("populate_nifti_image") ;
176 /*-- create nifti_image structure --*/
177
178 nim = (nifti_image *) calloc( 1 , sizeof(nifti_image) ) ;
179 if( !nim ) {
180 fprintf(stderr, "** ERROR: failed to allocate nifti image\n");
181 RETURN(NULL) ;
182 }
183
184 /*-- calculate and set ndim and intents --*/
185 /* cases to allow for:
186 1. 3d+time dataset
187 2. 3d func bucket
188 (not going to handle this currently, may extend later)
189 -- RWC: modified so that 'u' dimension is for buckets
190 3. 3d single brick
191 4. 2d and 1d spatial
192 (if 2d or 1d spatial + time, treat as #1)
193 5. Don't know of any vectors in AFNI, so won't do those either
194 6. Single 3d (or 2d or 1d) functional brik */
195
196 if (dset->dblk->nvals > 1) {
197 STATUS("4D dataset") ;
198 /* RWC: bucket stored as 5th dimen (but require ntt > 1) */
199 if (dset->taxis && dset->taxis->ntt > 1) nim->ndim = 4;
200 else nim->ndim = 5;
201
202 /*-- check sub-bricks for uniformity in type and scale --*/
203
204 type0 = DSET_BRICK_TYPE(dset,0) ;
205 fac0 = DSET_BRICK_FACTOR(dset,0) ;
206
207 for( ii=1 ; ii < DSET_NVALS(dset) ; ii++ ){
208 if( DSET_BRICK_TYPE(dset,ii) != type0){
209 fprintf(stderr,
210 "** ERROR: CANNOT WRITE NIfTI FILE; BRICK DATA TYPES NOT CONSISTENT\n") ;
211 RETURN(NULL);
212 } else if( DSET_BRICK_FACTOR(dset,ii) != fac0) {
213 fprintf(stderr,
214 "** ERROR: CANNOT WRITE NIfTI FILE; BRICK FACTORS NOT CONSISTENT\n") ;
215 fprintf(stderr,
216 " (consider transforming to a float dataset before performing\n"
217 " this operation, or consider '3dAFNItoNIFTI -float')\n"
218 ) ;
219 RETURN(NULL);
220 }
221 }
222 } else { /* we only have one brick */
223 STATUS("3D dataset") ;
224 if( options.debug_level > 1 ) fprintf(stderr,"-- PNI: write one brick\n") ;
225 type0 = DSET_BRICK_TYPE(dset,0);
226 fac0 = DSET_BRICK_FACTOR(dset,0) ;
227 if (ISFUNC(dset)) {
228 STATUS("functional dataset") ;
229 if( options.debug_level > 1 )
230 fprintf(stderr,"-- PNI: functional brick\n") ;
231 nim->intent_code = DSET_BRICK_STATCODE(dset,0);
232 if (nim->intent_code < 0) nim->intent_code = dset->func_type ;
233 if (nim->intent_code < 0) nim->intent_code = NIFTI_INTENT_NONE ;
234 /* 3dbucket func_type=FUNC_BUCK_TYPE becomes NIFTI_INTENT_NORMAL, which
235 AFNI whines about... 3dD uses 'none' for stat type on betas, so
236 stick with that 6 Jul 2013 [rickr] */
237 if (nim->intent_code >= FUNC_BUCK_TYPE)
238 nim->intent_code = NIFTI_INTENT_NONE ;
239 if( options.debug_level > 1 )
240 fprintf(stderr,"-- PNI: stat code = %d !!!\n",nim->intent_code) ;
241 if(PRINT_TRACING){
242 char str[256]; sprintf(str,"intent_code = %d",nim->intent_code);STATUS(str);
243 }
244 if (nim->intent_code > -1) {
245 nparam = FUNC_need_stat_aux[nim->intent_code];
246 /* statpars should be 0-based 20 Oct 2010 [rickr] */
247 if (nparam >= 1) nim->intent_p1 = DSET_BRICK_STATPAR(dset,0,0);
248 if (nparam >= 2) nim->intent_p2 = DSET_BRICK_STATPAR(dset,0,1);
249 if (nparam == 3) nim->intent_p3 = DSET_BRICK_STATPAR(dset,0,2);
250 }
251 }
252 if (dset->daxes->nzz > 1) {
253 nim->ndim = 3 ;
254 } else if (dset->daxes->nyy > 1) {
255 nim->ndim = 2 ;
256 } else {
257 nim->ndim = 1;
258 }
259 }
260
261
262 /*-- set datatype, size, etc. --*/
263
264 STATUS("set datatype") ;
265 switch(type0) {
266 case MRI_byte:
267 nim->datatype = DT_UNSIGNED_CHAR;
268 nim->nbyper = 1 ;
269 break;
270 case MRI_short:
271 nim->datatype = DT_SIGNED_SHORT;
272 nim->nbyper = 2 ;
273 break;
274 case MRI_int:
275 nim->datatype = DT_SIGNED_INT;
276 nim->nbyper = 4 ;
277 break;
278 case MRI_float:
279 nim->datatype = DT_FLOAT;
280 nim->nbyper = 4 ;
281 break;
282 case MRI_double:
283 nim->datatype = DT_DOUBLE;
284 nim->nbyper = 8 ;
285 break;
286 case MRI_complex:
287 nim->datatype = DT_COMPLEX;
288 nim->nbyper = 8 ;
289 break;
290 case MRI_rgb:
291 nim->datatype = DT_RGB24;
292 nim->nbyper = 3 ;
293 break;
294 case MRI_rgba:
295 fprintf(stderr,
296 "** ERROR: Can't write NIfTI file since dataset is RGBA: %s\n",
297 options.infile_name) ;
298 RETURN(NULL) ;
299 break;
300 default:
301 fprintf(stderr,
302 "** ERROR: Can't write NIfTI file since datatype is unknown: %s\n",
303 options.infile_name) ;
304 RETURN(NULL) ;
305 break;
306 }
307
308 /*-- scaling --*/
309
310 nim->scl_slope = fac0 ;
311 nim->scl_inter = 0 ;
312
313 /*-- spatial transforms --*/
314
315 STATUS("set orientation") ;
316
317 axcode[0] = ORIENT_xyz[ dset->daxes->xxorient ] ; axnum[0] = dset->daxes->nxx ;
318 axcode[1] = ORIENT_xyz[ dset->daxes->yyorient ] ; axnum[1] = dset->daxes->nyy ;
319 axcode[2] = ORIENT_xyz[ dset->daxes->zzorient ] ; axnum[2] = dset->daxes->nzz ;
320
321 axsign[0] = ORIENT_sign[ dset->daxes->xxorient ] ;
322 axsign[1] = ORIENT_sign[ dset->daxes->yyorient ] ;
323 axsign[2] = ORIENT_sign[ dset->daxes->zzorient ] ;
324
325 axstep[0] = dset->daxes->xxdel ; axstart[0] = dset->daxes->xxorg ;
326 axstep[1] = dset->daxes->yydel ; axstart[1] = dset->daxes->yyorg ;
327 axstep[2] = dset->daxes->zzdel ; axstart[2] = dset->daxes->zzorg ;
328
329 for (ii = 0 ; ii < 3 ; ii++ ) {
330 if (axcode[ii] == 'x') {
331 nif_x_axnum = ii ;
332 } else if (axcode[ii] == 'y') {
333 nif_y_axnum = ii ;
334 } else nif_z_axnum = ii ;
335 }
336
337 nim->qto_xyz.m[0][0] = nim->qto_xyz.m[0][1] = nim->qto_xyz.m[0][2] =
338 nim->qto_xyz.m[1][0] = nim->qto_xyz.m[1][1] = nim->qto_xyz.m[1][2] =
339 nim->qto_xyz.m[2][0] = nim->qto_xyz.m[2][1] = nim->qto_xyz.m[2][2] = 0.0 ;
340
341 /*-- set voxel and time deltas and units --*/
342
343 nim->dx = nim->pixdim[1] = fabs ( axstep[0] ) ;
344 nim->dy = nim->pixdim[2] = fabs ( axstep[1] ) ;
345 nim->dz = nim->pixdim[3] = fabs ( axstep[2] ) ;
346
347 nim->du = nim->pixdim[5] = 0 ;
348 nim->dv = nim->pixdim[6] = 0 ;
349 nim->dw = nim->pixdim[7] = 0 ;
350
351 #if 0
352 val = (axsign[nif_x_axnum] == '+') ? -1 : 1 ;
353 nim->qto_xyz.m[0][nif_x_axnum] = val * nim->pixdim[nif_x_axnum + 1];
354 val = (axsign[nif_y_axnum] == '+') ? -1 : 1 ;
355 nim->qto_xyz.m[1][nif_y_axnum] = val * nim->pixdim[nif_y_axnum + 1];
356 val = (axsign[nif_x_axnum] == '-') ? -1 : 1 ;
357 nim->qto_xyz.m[2][nif_z_axnum] = val * nim->pixdim[nif_z_axnum + 1];
358 #else
359 nim->qto_xyz.m[0][nif_x_axnum] = - axstep[nif_x_axnum];
360 nim->qto_xyz.m[1][nif_y_axnum] = - axstep[nif_y_axnum];
361 nim->qto_xyz.m[2][nif_z_axnum] = axstep[nif_z_axnum];
362 #endif
363
364 /* nifti origin stuff */
365
366 #if 0
367 nim->qoffset_x = axstart[nif_x_axnum] ;
368 if (axsign[nif_x_axnum] == '+') nim->qoffset_x = - nim->qoffset_x ;
369 nim->qoffset_y = axstart[nif_y_axnum];
370 if (axsign[nif_y_axnum] == '+') nim->qoffset_y = - nim->qoffset_y ;
371 nim->qoffset_z = axstart[nif_z_axnum];
372 if (axsign[nif_z_axnum] == '-') nim->qoffset_z = - nim->qoffset_z ;
373 #endif
374
375 nim->qoffset_x = -axstart[nif_x_axnum] ;
376 nim->qoffset_y = -axstart[nif_y_axnum];
377 nim->qoffset_z = axstart[nif_z_axnum];
378
379 #if 0
380 nim->qoffset_x = -axstart[0] ;
381 nim->qoffset_y = -axstart[1];
382 nim->qoffset_z = axstart[2];
383 #endif
384
385 nim->qto_xyz.m[0][3] = nim->qoffset_x ;
386 nim->qto_xyz.m[1][3] = nim->qoffset_y ;
387 nim->qto_xyz.m[2][3] = nim->qoffset_z ;
388
389
390 /*-- from the same above info, set the sform matrix to equal the qform --*/
391 /* KRH 7/6/05 - using sform to duplicate qform for
392 interoperability with FSL */
393 /* update with oblique transformation if available DRG 24 May 2007 */
394 /* check for valid transformation matrix */
395 if(!ISVALID_MAT44(dset->daxes->ijk_to_dicom_real)) {
396 nim->sto_xyz = nim->qto_xyz; /* copy qform to sform */
397 }
398 else {
399 /* fill in sform with AFNI daxes transformation matrix */
400 /* n2 10 Jul, 2015 [rickr] */
401 nifti_mat44_to_dmat44(&dset->daxes->ijk_to_dicom_real, &nim->sto_xyz);
402 /* negate first two rows of sform for NIFTI - LPI standard versus
403 AFNI RAI "DICOM" standard */
404 for( ii = 0; ii < 2; ii++) {
405 for (jj = 0 ; jj < 4; jj++) {
406 nim->sto_xyz.m[ii][jj] = -(nim->sto_xyz.m[ii][jj]);
407 }
408 }
409 /* update qform too with struct copy from sform*/
410 nim->qto_xyz= nim->sto_xyz ;
411
412 }
413
414 /*-- from the above info, calculate the quaternion qform --*/
415
416 STATUS("set quaternion") ;
417
418 /* n2 10 Jul, 2015 [rickr] */
419 nifti_dmat44_to_quatern( nim->qto_xyz ,
420 &nim->quatern_b, &nim->quatern_c, &nim->quatern_d,
421 &dumqx, &dumqy, &dumqz, &dumdx, &dumdy, &dumdz,
422 &nim->qfac ) ;
423
424 /*-- verify dummy quaternion parameters --*/
425
426 if( options.debug_level > 2 )
427 /* n2 10 Jul, 2015 [rickr] */
428 fprintf(stderr,"++ Quaternion check:\n"
429 "%f , %f\n %f , %f\n %f , %f\n %f , %f\n %f , %f\n %f , %f\n; %f\n",
430 nim->qoffset_x, dumqx, nim->qoffset_y,dumqy, nim->qoffset_z, dumqz ,
431 nim->dx, dumdx , nim->dy, dumdy , nim->dz, dumdz, nim->qfac ) ;
432
433 /*-- calculate inverse qform --*/
434
435 /* n2 10 Jul, 2015 [rickr] */
436 nim->qto_ijk = nifti_dmat44_inverse( nim->qto_xyz ) ;
437
438 /*-- set dimensions of grid array --*/
439
440 nim->nt = nim->nu = nim->nv = nim->nw = 1 ;
441 nim->nx = axnum[0] ;
442 nim->ny = axnum[1] ;
443 nim->nz = axnum[2] ;
444
445 /* taxis might be set via EDIT_empty_copy, even if dset is not time series */
446 /* [1 Jun 2020 rickr] */
447 if (dset->taxis != NULL && (dset->taxis->ntt > 1 || DSET_NVALS(dset) == 1))
448 nim->nt = DSET_NUM_TIMES(dset) ; /* time is 4th dimension */
449 else
450 nim->nu = DSET_NVALS(dset) ; /* RWC: bucket is 5th dimension */
451
452 if ( nim->nt > 1){
453 float TR = dset->taxis->ttdel ;
454 if( DSET_TIMEUNITS(dset) == UNITS_MSEC_TYPE ) TR *= 0.001; /* 10 May 2005 */
455 nim->dt = nim->pixdim[4] = TR ;
456 }
457
458 nim->dim[0] = nim->ndim;
459 nim->dim[1] = nim->nx;
460 nim->dim[2] = nim->ny;
461 nim->dim[3] = nim->nz;
462 nim->dim[4] = nim->nt; /* RWC: at most one of nt and nu is > 1 */
463 nim->dim[5] = nim->nu;
464 nim->dim[6] = nim->nv;
465 nim->dim[7] = nim->nw;
466
467 nim->nvox = nim->nx * nim->ny * nim->nz * nim->nt
468 * nim->nu * nim->nv * nim->nw ;
469
470 /*-- slice timing --*/
471
472 nim->freq_dim = nim->phase_dim = 0 ;
473 if (dset->taxis != NULL) { /* if time axis exists */
474 nim->slice_dim = 3 ;
475 nim->slice_duration = 0 ;
476 nim->slice_start = 0 ;
477 nim->slice_end = nim->nz - 1;
478 nim->toffset = DSET_TIMEORIGIN(dset);
479
480 /*-- this bit assumes that afni slice timing offsets *
481 *-- are created starting from zero and including all *
482 *-- slices initially. They may later be modified by *
483 *-- zero padding at either end. No other *
484 *-- modifications are intentionally accepted right now. */
485
486 if (DSET_NUM_TTOFF(dset) > 0 ) { /* if time offset exists */
487
488 /*-- Find first and last non-zero element */
489 #define MYEPSILON 0.00001
490 #define MYFPEQ(a, b) (fabs((a) - (b)) < MYEPSILON)
491
492 tlist = dset->taxis->toff_sl;
493 for (ii = 0 ; ii < nim->nz ; ii++ ) {
494 if (!MYFPEQ(tlist[ii],0.0)) break ;
495 }
496 sfirst = ii ;
497 for (ii = nim->nz - 1 ; ii >= sfirst ; ii-- ) {
498 if (!MYFPEQ(tlist[ii],0.0)) break ;
499 }
500 slast = ii ;
501
502 /* pattern check re-written to deal with including zeros */
503 /* on either end 14 Jun 2006 [rickr] */
504
505 pattern = NIFTI_SLICE_UNKNOWN;
506
507 /* do we have all zeros? */
508 if( sfirst == slast && MYFPEQ(tlist[sfirst],0.0) ) {
509 nim->slice_duration = 0.0;
510 } else { /* see if there is a known pattern in the list */
511 tlen = slast-sfirst+2;
512
513 /* try including leading adjacent zero in the pattern, first */
514 if( sfirst > 0 ) {
515 /* n2 10 Jul, 2015 [rickr] */
516 pattern = get_slice_timing_pattern(tlist+sfirst-1, tlen, &fsdur);
517 nim->slice_duration = (double)fsdur;
518 if( pattern != NIFTI_SLICE_UNKNOWN ) sfirst--;
519 }
520
521 /* try including trailing adjacent zero in the pattern, next */
522 if( pattern == NIFTI_SLICE_UNKNOWN && slast < nim->nz-1 ) {
523 pattern = get_slice_timing_pattern(tlist+sfirst, tlen, &fsdur);
524 nim->slice_duration = (double)fsdur;
525 if( pattern != NIFTI_SLICE_UNKNOWN ) slast++;
526 }
527
528 /* if no pattern yet, try list without zeros */
529 if( pattern == NIFTI_SLICE_UNKNOWN ) {
530 pattern = get_slice_timing_pattern(tlist+sfirst, tlen-1, &fsdur);
531 nim->slice_duration = (double)fsdur;
532 }
533
534 if( pattern == NIFTI_SLICE_UNKNOWN ) {
535 nim->slice_code = pattern ;
536 nim->slice_start = 0 ;
537 nim->slice_end = 0 ;
538 nim->slice_duration = 0.0 ;
539 } else {
540 nim->slice_start = sfirst ;
541 nim->slice_end = slast ;
542 nim->slice_code = pattern;
543 }
544
545 if( options.debug_level > 1)
546 fprintf(stderr,"+d timing pattern '%s', slice %d to %d, stime %f\n",
547 nifti_slice_string(pattern), sfirst, slast, nim->slice_duration);
548 }
549
550 /* if toffset is 0 and the timing patter is known and the minimum
551 * slice offset is positive, the toffset to that minimum
552 * 12 Oct 2007 [rickr] */
553 if( nim->toffset == 0.0 && nim->slice_code != NIFTI_SLICE_UNKNOWN ){
554 float tmin = tlist[0];
555 for (ii = 1 ; ii < nim->nz ; ii++ )
556 if( tlist[ii] < tmin ) tmin = tlist[ii] ;
557 if( tmin > 0.0 ) nim->toffset = tmin ;
558 }
559 }
560
561 nim->time_units = NIFTI_UNITS_SEC ;
562
563 } else { /* if time axis not exists */
564 nim->slice_dim = 0 ;
565 nim->time_units = NIFTI_UNITS_UNKNOWN ;
566 }
567
568 /*-- byte order --*/
569
570 nim->byteorder = nifti_short_order() ;
571
572 /* KRH 7/25/05 modified to note talairach view into NIfTI file */
573
574 nim->qform_code = space_to_NIFTI_code(dset);
575
576 // NB: this is commented out. If ever brought back into the fold
577 // (unlikely!), take care with [qs]form changes. Probably delete.
578 /*#if 0
579 if ( dset->view_type == VIEW_TALAIRACH_TYPE ) {
580 nim->qform_code = NIFTI_XFORM_TALAIRACH ;
581 } else {
582 nim->qform_code = NIFTI_XFORM_SCANNER_ANAT ;
583 }
584 #endif*/
585
586 nim->sform_code = nim->qform_code ; /* KRH 7/6/05 - using */
587 /* sform to duplicate qform for interoperability with FSL */
588
589
590 /*-- odds and ends that are constant for AFNI files --*/
591 nim->cal_min = nim->cal_max = 0 ;
592 nim->nifti_type = 1 ;
593 nim->xyz_units = NIFTI_UNITS_MM ;
594 nim->num_ext = 0;
595 nim->ext_list = NULL ;
596 nim->iname_offset = 352 ; /* until extensions are added */
597 nim->data = NULL ;
598
599 RETURN(nim) ;
600 }
601
602 /*-------------------------------------------------------------------*/
603 /*! List of dataset attributes NOT to save in a NIfTI-1.1 file. -----*/
604
605 static char *badlist[] = {
606 "IDCODE_STRING" , /* this goes in the NI_group header */
607 "DATASET_RANK" ,
608 "DATASET_DIMENSIONS" ,
609 "TYPESTRING" ,
610 "SCENE_DATA" ,
611 "ORIENT_SPECIFIC" ,
612 "ORIGIN" ,
613 "DELTA" ,
614 "TAXIS_NUMS" ,
615 "TAXIS_FLOATS" ,
616 "TAXIS_OFFSETS" ,
617 "BYTEORDER_STRING" ,
618 "BRICK_TYPES" ,
619 "BRICK_FLOAT_FACS" ,
620 "STAT_AUX" ,
621 "LABEL_1" ,
622 "LABEL_2" ,
623 "DATASET_NAME" ,
624 NULL } ;
625
626 /*-------------------------------------------------------------------*/
627 /*! Create the AFNI extension string for a NIfTI-1.1 file, and insert
628 this metadata into the nifti_image struct for output to disk.
629 - If something bad happens, fails silently
630 - 09 May 2005 - RWCox
631 ---------------------------------------------------------------------*/
632
nifti_set_afni_extension(THD_3dim_dataset * dset,nifti_image * nim)633 void nifti_set_afni_extension( THD_3dim_dataset *dset , nifti_image *nim )
634 {
635 NI_group *ngr ;
636 NI_element *nel ;
637 NI_stream ns ;
638 char *rhs , buf[128] ;
639 int ii,bb , npart,*bpart ;
640
641 if( nim == NULL ) return ; /* stupid or evil caller */
642 if( AFNI_yesenv("AFNI_NIFTI_NOEXT") ) return ; /* not allowed */
643
644 /** write all dataset 'attributes' into a NIML group */
645
646 ngr = THD_nimlize_dsetatr( dset ) ;
647 if( ngr == NULL ) return ; /* bad */
648 NI_rename_group( ngr , "AFNI_attributes" ) ;
649
650 /* 12 May 2005: add a signature to check the file on input to AFNI */
651
652 /* n2 10 Jul, 2015 [rickr] */
653 sprintf(buf,"%d,%d,%d,%d,%d,%d" ,
654 (int)nim->nx, (int)nim->ny, (int)nim->nz, (int)nim->nt, (int)nim->nu, (int)nim->datatype ) ;
655 NI_set_attribute( ngr , "NIfTI_nums" , buf ) ;
656
657 /** now, scan attribute elements in the group, and mark some
658 of them as being useless or redundant in the NIfTI world **/
659
660 npart = ngr->part_num ;
661 bpart = (int *)calloc(sizeof(int),npart) ;
662 for( ii=0 ; ii < npart ; ii++ ){
663 if( ngr->part_typ[ii] != NI_ELEMENT_TYPE ) continue ;
664 nel = (NI_element *) ngr->part[ii] ;
665 if( strcmp(nel->name,"AFNI_atr") != 0 ) continue ;
666
667 /* to make this effective, change AFNI_name -> atr_name */
668 rhs = NI_get_attribute( nel , "AFNI_name" ) ;
669 if( rhs == NULL ) continue ;
670
671 for( bb=0 ; badlist[bb] != NULL ; bb++ )
672 if( strcmp(rhs,badlist[bb]) == 0 ){ bpart[ii] = 1; break; }
673 }
674
675 /** remove marked attributes from the NIML group **/
676
677 for( ii=npart-1 ; ii >= 0 ; ii-- ){
678 if( bpart[ii] )
679 NI_remove_from_group( ngr , ngr->part[ii] ) ;
680 }
681 free((void *)bpart) ; /* done with this */
682 if( ngr->part_num <= 0 ){ NI_free_element(ngr); return; }
683
684 /** format into a character string to put in the NIfTI-1.1 extension **/
685
686 ns = NI_stream_open( "str:" , "w" ) ;
687 NI_stream_writestring( ns , "<?xml version='1.0' ?>\n" ) ;
688 NI_write_element( ns , ngr , NI_TEXT_MODE ) ;
689 rhs = NI_stream_getbuf( ns ) ;
690
691 /** write contents of the string into the nifti_image struct **/
692
693 nifti_add_extension( nim , rhs , strlen(rhs)+1 , NIFTI_ECODE_AFNI ) ;
694
695 NI_stream_close(ns) ; /* frees the string buffer, too */
696 NI_free_element(ngr) ; /* done with this trashola */
697 return ;
698 }
699
700
701 /*! given a list of floats, detect any slice timing pattern
702 * 14 Jun 2006 [rickr]
703 *
704 * - if a pattern is found and delta is set, return the time delta
705 * - return one of:
706 * NIFTI_SLICE_UNKNOWN,
707 * NIFTI_SLICE_SEQ_INC, NIFTI_SLICE_SEQ_DEC,
708 * NIFTI_SLICE_ALT_INC, NIFTI_SLICE_ALT_DEC,
709 * NIFTI_SLICE_ALT_INC2, NIFTI_SLICE_ALT_DEC2,
710 */
get_slice_timing_pattern(float * times,int len,float * delta)711 static int get_slice_timing_pattern( float * times, int len, float * delta )
712 {
713 float * flist, diff;
714 int * ilist;
715 int c, index, pattern;
716
717 ENTRY("get_slice_timing_pattern");
718
719 if( delta ) *delta = 0.0; /* init, in case of early return */
720 if( ! times || len < 2 ) RETURN(NIFTI_SLICE_UNKNOWN);
721
722 /* if the length is very short, deal with it separately */
723 if( len == 2 ){
724 if( delta ) *delta = fabs(times[1]-times[0]);
725 if( times[1] > times[0] ) RETURN(NIFTI_SLICE_SEQ_INC);
726 else RETURN(NIFTI_SLICE_SEQ_DEC);
727 }
728
729 /*** sort the list, and look for a linear pattern ***/
730
731 /* duplicate list */
732 flist = (float *)malloc(len * sizeof(float));
733 ilist = (int *)malloc(len * sizeof(int));
734 if(!flist || !ilist) { /* yeah, lazy... */
735 ERROR_message(" GSTP: cannot dupe timing list\n");
736 RETURN(NIFTI_SLICE_UNKNOWN);
737 }
738 memcpy(flist, times, len*sizeof(float));
739 for(c = 0; c < len; c++) ilist[c] = c; /* init ilist with current indices */
740
741 /* sort flist, with ilist returning original indices */
742 qsort_floatint(len, flist, ilist);
743
744 /* and check for a fixed difference */
745 diff = flist[1] - flist[0];
746 pattern = 1;
747 for( c = 1; c < len-1; c++ )
748 if( !MYFPEQ(diff, (flist[c+1]-flist[c])) ) { pattern = 0; break; }
749
750 /* if no pattern, just return failure */
751 if( !pattern ) {
752 free(flist); free(ilist); RETURN(NIFTI_SLICE_UNKNOWN);
753 }
754
755 /* we have linear offsets, now see if the slices match a known pattern */
756 /* repeatedly: init to a pattern, and see if it fails */
757
758 /* SEQ_INC (0,1,2,3...,l-1) */
759 pattern = NIFTI_SLICE_SEQ_INC; index = 0;
760 for( c = 0; c < len; c++ ) {
761 if( ilist[c] != index ) { pattern = NIFTI_SLICE_UNKNOWN; break; }
762 index++;
763 }
764
765 if( pattern == NIFTI_SLICE_UNKNOWN ) { /* (l-1,l-2,...2,1,0) */
766 pattern = NIFTI_SLICE_SEQ_DEC; index = len-1;
767 for( c = 0; c < len; c++ ) {
768 if( ilist[c] != index ) { pattern = NIFTI_SLICE_UNKNOWN; break; }
769 index--;
770 }
771 }
772
773 if( pattern == NIFTI_SLICE_UNKNOWN ) { /* (0,2,4,6,...,1,3,5,...) */
774 pattern = NIFTI_SLICE_ALT_INC; index = 0;
775 for( c = 0; c < len; c++ ) {
776 if( ilist[c] != index ) { pattern = NIFTI_SLICE_UNKNOWN; break; }
777 index += 2; if( index >= len ) index = 1; /* so no parity issue */
778 }
779 }
780
781 if( pattern == NIFTI_SLICE_UNKNOWN ) { /* (l-1,l-3,...1/0,l-2,l-4,...,0/1) */
782 pattern = NIFTI_SLICE_ALT_DEC; index = len-1;
783 for( c = 0; c < len; c++ ) {
784 if( ilist[c] != index ) { pattern = NIFTI_SLICE_UNKNOWN; break; }
785 index -= 2; if( index < 0 ) index = len-2;
786 }
787 }
788
789 if( pattern == NIFTI_SLICE_UNKNOWN ) { /* (1,3,5,...,0,2,4...) */
790 pattern = NIFTI_SLICE_ALT_INC2; index = 1;
791 for( c = 0; c < len; c++ ) {
792 if( ilist[c] != index ) { pattern = NIFTI_SLICE_UNKNOWN; break; }
793 index += 2; if( index >= len ) index = 0;
794 }
795 }
796
797 if( pattern == NIFTI_SLICE_UNKNOWN ) { /* (l-2,l-4,...4,2,0,l-1,...,5,3,1) */
798 pattern = NIFTI_SLICE_ALT_DEC2; index = len-2;
799 for( c = 0; c < len; c++ ) {
800 if( ilist[c] != index ) { pattern = NIFTI_SLICE_UNKNOWN; break; }
801 index -= 2; if( index < 0 ) index = len-1;
802 }
803 }
804
805 if( delta && pattern != NIFTI_SLICE_UNKNOWN ) *delta = diff;
806
807 /** done, whatever the case may be **/
808 free(flist); free(ilist);
809 RETURN(pattern);
810 }
811
812 /* set NIFTI sform code based on atlas space */
space_to_NIFTI_code(THD_3dim_dataset * dset)813 int space_to_NIFTI_code(THD_3dim_dataset *dset)
814 {
815 char *genspc = NULL;
816 /* several changes for generic spaces and defaults 05/02/2012 -mod drg */
817 /* use generic space or space of dataset to choose xform code */
818 genspc = THD_get_generic_space(dset);
819 if(genspc == NULL)
820 return(NIFTI_XFORM_SCANNER_ANAT); // sform=1
821
822 if (strcmp(genspc,"TLRC") == 0) { // sform=3
823 return(NIFTI_XFORM_TALAIRACH);
824 }
825 if (strcmp(genspc,"MNI") == 0) { // sform=4
826 return(NIFTI_XFORM_MNI_152);
827 }
828 /* call MNI_ANAT as aligned to something else*/
829 /* [PT,DRG: Mar 8, 2019] CHANGE: MNI_ANAT will now go to sform=5.
830 NOTE that the SPM Anatomy Toolbox is distributed with sform=2.
831 This may cause issues with migrations across soft(be)wares. */
832 //if (strcmp(genspc,"MNI_ANAT") == 0) { // sform=2???!?!
833 // return(NIFTI_XFORM_ALIGNED_ANAT);
834 // }
835 if((strcmp(genspc,"ORIG") == 0) ||
836 (strcmp(genspc,"ACPC") == 0)){
837 return(NIFTI_XFORM_SCANNER_ANAT); // sform=1
838 }
839 /* make catch-all for other spaces as an aligned space
840 alternative is use SCANNER_ANAT again or UNKNOWN */
841 //return(NIFTI_XFORM_ALIGNED_ANAT);
842 // [PT,DRG: Mar 8, 2019]
843 return(NIFTI_XFORM_TEMPLATE_OTHER); // sform=5
844 }
845