1 #include "mrilib.h"
2
3 #ifndef DONT_INCLUDE_ANALYZE_STRUCT
4 #define DONT_INCLUDE_ANALYZE_STRUCT
5 #endif
6 #include "nifti2_io.h" /** will include nifti1.h **/
7 static void NIFTI_code_to_space(int code,THD_3dim_dataset *dset);
8 static int NIFTI_code_to_view(int code, char *atlas_space);
9 static int NIFTI_default_view();
10 extern char *THD_get_space(THD_3dim_dataset *dset);
11 static int THD_nifti_process_afni_ext(THD_3dim_dataset *dset,
12 nifti_image *nim, char *pathnew);
13
14 #undef KILL_pathnew
15 #define KILL_pathnew if( pathnew != NULL && pathnew != pathname ) free(pathnew)
16
17 #undef NULLRET
18 #define NULLRET do{ KILL_pathnew ; RETURN(NULL) ; } while(0)
19
20 /*******************************************************************/
21 /********** 26 Aug 2003: read a NIFTI-1 file as a dataset **********/
22 /*******************************************************************/
23
24 /*
25 [PT,DRG: Mar 8, 2019] Updated to include [qs]form_code = 5
26
27 [PT,DRG: Mar 12, 2019] some fixes for [qs]form behavior. Still need
28 to look at [qs]form=2 case, changing space/views.
29
30 */
31
THD_open_nifti(char * pathname)32 THD_3dim_dataset * THD_open_nifti( char *pathname )
33 {
34 THD_3dim_dataset *dset=NULL ;
35 nifti_image *nim ;
36 int ntt , nbuc , nvals ;
37 int use_qform = 0 , use_sform = 0, form_code = 0 ;
38 int statcode = 0 , datum , iview , ibr ;
39 int scale_data = 0 ; /* flag based on scl_slope and inter 20 Jun 2008 */
40 int xform_data = 0;
41 THD_ivec3 orixyz , nxyz ;
42 THD_fvec3 dxyz , orgxyz ;
43 THD_mat33 R ;
44 mat44 ijk_to_dicom44 ;
45 char *ppp , prefix[THD_MAX_PREFIX] ;
46 char form_priority = 'S' ; /* 23 Mar 2006 */
47 static int n_xform_warn=0;
48 char *pathnew = NULL ;
49 char *rpath = NULL , *rp=NULL ;;
50
51 ENTRY("THD_open_nifti") ;
52
53 /*-- open input file --*/
54
55 { /* set the nifti_io debug level 8 Apr 2005 [rickr] */
56 char * ept = my_getenv("AFNI_NIFTI_DEBUG");
57 if( ept != NULL ) nifti_set_debug_level(atoi(ept));
58 }
59
60 nifti_set_alter_cifti(1) ; /* if CIFTI, shift dims 23 Jul 2015 [rickr] */
61
62 /* pathname munge by RWC [07 Feb 2019] */
63
64 if( pathname == NULL || *pathname == '\0' ) NULLRET ;
65 if( THD_is_file(pathname) ){
66 pathnew = pathname ;
67 } else if( strstr(pathname,".gz") == NULL ){
68 pathnew = (char *)malloc(sizeof(char)*(strlen(pathname)+16)) ;
69 sprintf(pathnew,"%s.gz",pathname) ;
70 if( ! THD_is_file(pathnew) ) NULLRET ;
71 } else {
72 NULLRET ;
73 }
74 if( !AFNI_yesenv("AFNI_NOREALPATH") && !THD_is_symlink(pathnew) ){ /* 21 Mar 2020 */
75 rpath = (char *)malloc(sizeof(char)*RPMAX) ;
76 rp = realpath( pathnew , rpath ) ;
77 if( rp == NULL ){ free(rpath) ; NULLRET ; }
78 KILL_pathnew ; pathnew = rpath ;
79 }
80 #if 0
81 INFO_message("NIFTI - opening dataset %s",pathnew) ;
82 #endif
83
84 /*-- Read read read --*/
85
86 nim = nifti_image_read( pathnew, 0 ) ;
87
88 if( nim == NULL || nim->nifti_type == 0 ) NULLRET ;
89
90 if( pathname != pathnew && pathnew != rpath )
91 ININFO_message("reading %s instead of non-existing %s",pathnew,pathname) ;
92
93 /*-- extract some useful AFNI-ish information from the nim struct --*/
94
95 /* we must have at least 2 spatial dimensions */
96
97 /* this should be okay 11 Jun 2007 */
98 /* if( nim->nx < 2 || nim->ny < 2 ) NULLRET ; */
99
100 /* 4th dimension = time; 5th dimension = bucket:
101 these are mutually exclusive in AFNI at present */
102
103 ntt = nim->nt ; nbuc = nim->nu ;
104
105 /* nt and nu might be 0 now (see irritating niftilib 1.17 update) */
106 /* so ensure that ntt and nbuc are positive 02 Mar 2006 [rickr] */
107
108 if( ntt <= 0 ) ntt = 1;
109 if( nbuc <= 0 ) nbuc = 1;
110
111 if( nim->nz <= 0 ) nim->nz = 1 ; /* 03 Mar 2006: RWCox */
112
113 if( ntt > 1 && nbuc > 1 ){
114 fprintf(stderr,
115 "** AFNI can't deal with 5 dimensional NIfTI(%s)\n",
116 pathnew ) ;
117 NULLRET ;
118 }
119
120 nvals = MAX(ntt,nbuc) ;
121
122 /* collapse higher-dimensional datasets 23 Jul 2015 [rickr] */
123 /* (this includes CIFTI) */
124 if( nim->nv > 1 ) nvals *= nim->nv;
125 if( nim->nw > 1 ) nvals *= nim->nw;
126 if( ntt > 1 ) ntt = nvals;
127 else nbuc = nvals;
128
129 /* determine type of dataset values:
130 if we are scaling, or if the data type in the NIfTI file
131 is something AFNI can't handle, then the result will be floats */
132
133 /* do not scale if slope is 0 or if slope is 1 and inter is 0 */
134 if( !isfinite(nim->scl_slope) || !isfinite(nim->scl_inter) ){
135 fprintf(stderr,"** bad scl_slope and inter = %f, %f, ignoring...\n",
136 nim->scl_slope, nim->scl_inter);
137 } else {
138 /* do not scale if either slope or inter is zero
139 * slope==0 : treat as unset (rather than constant data)
140 * inter==0 : use slope as brick_fac */
141
142 /* Any NIFTI scalar previously implied conversion to float.
143 * No longer scale if inter==0 && slope!=0 (set brick_fac).
144 * Thanks to C Caballero and S Moia for reporting this.
145 * 26 Jan 2021 [rickr] */
146 /* still scale if slope == 1 && inter != 0 */
147 scale_data = nim->scl_slope != 0.0 && nim->scl_inter != 0.0 ;
148 }
149 { char *eee = getenv("AFNI_NIFTI_SCALE") ;
150 if( eee != NULL && toupper(*eee) == 'N' ) scale_data = 0 ;
151 }
152
153 switch( nim->datatype ){
154 default:
155 fprintf(stderr,
156 "** AFNI can't handle NIFTI datatype=%d (%s) in file %s\n",
157 nim->datatype, nifti_datatype_string(nim->datatype), pathnew );
158 NULLRET ;
159 break ;
160
161 case DT_UINT8: datum = scale_data ? MRI_float : MRI_byte ;
162 xform_data = scale_data;
163 break ;
164 case DT_INT16: datum = scale_data ? MRI_float : MRI_short ;
165 xform_data = scale_data;
166 break ;
167 case DT_FLOAT32: datum = MRI_float ; break ;
168 case DT_COMPLEX64: datum = MRI_complex ; break ;
169 case DT_RGB24: datum = MRI_rgb ; break ;
170
171 case DT_INT8: /* NIfTI-1 data types that AFNI can't handle directly */
172 case DT_UINT16:
173 case DT_INT32:
174 case DT_UINT32:
175 case DT_FLOAT64: datum = MRI_float ; xform_data = 1 ; break ;
176
177 #if 0
178 case DT_COMPLEX128: /* this case would be too much like real work */
179 fprintf(stderr,
180 "** AFNI convert NIFTI_datatype=%d (%s) in file %s to COMPLEX64\n",
181 nim->datatype, nifti_datatype_string(nim->datatype), pathnew );
182 datum = MRI_complex ;
183 break ;
184 #endif
185 }
186
187 if( xform_data && !AFNI_noenv("AFNI_NIFTI_TYPE_WARN")) {
188 if (!n_xform_warn || AFNI_yesenv("AFNI_NIFTI_TYPE_WARN")) {/* ZSS 04/11 */
189 fprintf(stderr,
190 "** AFNI converts NIFTI_datatype=%d (%s) in file %s to FLOAT32\n",
191 nim->datatype, nifti_datatype_string(nim->datatype), pathnew );
192 if (!AFNI_yesenv("AFNI_NIFTI_TYPE_WARN")) {
193 fprintf(stderr,
194 " Warnings of this type will be muted for this session.\n"
195 " Set AFNI_NIFTI_TYPE_WARN to YES to see them all, NO to see none.\n");
196 }
197 }
198 ++n_xform_warn;
199 }
200 /* check for statistics code */
201
202 if( nim->intent_code >= NIFTI_FIRST_STATCODE &&
203 nim->intent_code <= NIFTI_LAST_STATCODE ){
204
205 if( nim->intent_code > FUNC_PT_TYPE ){
206 fprintf(stderr,
207 "** AFNI doesn't understand NIFTI statistic type %d (%s) in file %s\n",
208 nim->intent_code , nifti_intent_string(nim->intent_code) , pathnew ) ;
209 } else {
210 statcode = nim->intent_code ;
211 if( nbuc > 1 ){
212 fprintf(stderr,
213 "** AFNI doesn't support NIFTI voxel-dependent statistic parameters"
214 " in file %s\n" , pathnew ) ;
215 statcode = 0 ;
216 }
217 }
218 }
219
220 /* 23 Mar 2006: set qform or sform as having priority -- RWCox */
221
222 ppp = my_getenv("AFNI_NIFTI_PRIORITY") ;
223 if( ppp != NULL ){
224 char fp = toupper(*ppp) ;
225 if( fp == 'S' || fp == 'Q' ) form_priority = fp ;
226 else WARNING_message("Illegal AFNI_NIFTI_PRIORITY='%s'",ppp) ;
227 }
228
229 /** 24 Mar 2006: check determs of qform and sform, if have both **/
230
231 if( nim->qform_code > 0 && nim->sform_code > 0 ){
232 float qdet , sdet ;
233 LOAD_MAT(R, nim->qto_xyz.m[0][0] ,
234 nim->qto_xyz.m[0][1] ,
235 nim->qto_xyz.m[0][2] ,
236 nim->qto_xyz.m[1][0] ,
237 nim->qto_xyz.m[1][1] ,
238 nim->qto_xyz.m[1][2] ,
239 nim->qto_xyz.m[2][0] ,
240 nim->qto_xyz.m[2][1] ,
241 nim->qto_xyz.m[2][2] ) ; qdet = MAT_DET(R) ;
242
243 LOAD_MAT(R, nim->sto_xyz.m[0][0] ,
244 nim->sto_xyz.m[0][1] ,
245 nim->sto_xyz.m[0][2] ,
246 nim->sto_xyz.m[1][0] ,
247 nim->sto_xyz.m[1][1] ,
248 nim->sto_xyz.m[1][2] ,
249 nim->sto_xyz.m[2][0] ,
250 nim->sto_xyz.m[2][1] ,
251 nim->sto_xyz.m[2][2] ) ; sdet = MAT_DET(R) ;
252
253 if( qdet*sdet < 0.0f )
254 WARNING_message("NIfTI('%s'): Qform/Sform handedness differ; %c wins!",
255 pathnew , form_priority ) ;
256 }
257
258 /* KRH 07/11/05 -- adding ability to choose spatial transform
259 from the options of qform, sform, bothform, or noform.
260
261 If sform is present, it will be used.
262
263 If qform is absent, but sform present, then the sform
264 will be modified to be an orthogonal rotation and used.
265
266 If both qform and sform are absent, then we will have
267 an error.
268
269 Previously assumed qform present. */
270
271 /* 23 Mar 2006: use form_priority to choose between them */
272
273 if ((nim->qform_code > 0) && (nim->sform_code > 0) ) {
274 if( form_priority == 'Q' ) { use_qform = 1 ; use_sform = 0 ; }
275 else { use_qform = 0 ; use_sform = 1 ; }
276 } else if (nim->qform_code > 0){ use_qform = 1 ; use_sform = 0 ; }
277 else if (nim->sform_code > 0){ use_qform = 0 ; use_sform = 1 ; }
278 else {
279 use_qform = 0 ; use_sform = 0 ;
280 WARNING_message(
281 "NO spatial transform (neither qform nor sform), in NIfTI file '%s'" ,
282 pathnew ) ;
283 }
284
285 /** now take NIfTI-1.1 coords and transform to AFNI codes **/
286
287 if (use_qform) {
288
289 float orgx, orgy, orgz ;
290
291 form_code = nim->qform_code;
292
293 /* determine orientation from the qto_xyz matrix,
294 which transforms (i,j,k) voxel indexes to (x,y,z) LPI coordinates */
295
296 LOAD_MAT(R, -nim->qto_xyz.m[0][0] , /* negate x and y */
297 -nim->qto_xyz.m[0][1] , /* coefficients, */
298 -nim->qto_xyz.m[0][2] , /* since AFNI works */
299 -nim->qto_xyz.m[1][0] , /* with RAI coords, */
300 -nim->qto_xyz.m[1][1] , /* but NIFTI uses */
301 -nim->qto_xyz.m[1][2] , /* LPI coordinates. */
302 nim->qto_xyz.m[2][0] , /* [Which is my own] */
303 nim->qto_xyz.m[2][1] , /* [damn fault!!!!!] */
304 nim->qto_xyz.m[2][2] ) ;
305
306 LOAD_MAT44(ijk_to_dicom44,
307 -nim->qto_xyz.m[0][0] , /* negate x and y */
308 -nim->qto_xyz.m[0][1] , /* coefficients, */
309 -nim->qto_xyz.m[0][2] , /* since AFNI works */
310 -nim->qto_xyz.m[0][3] ,
311 -nim->qto_xyz.m[1][0] , /* with RAI coords, */
312 -nim->qto_xyz.m[1][1] , /* but NIFTI uses */
313 -nim->qto_xyz.m[1][2] , /* LPI coordinates. */
314 -nim->qto_xyz.m[1][3] ,
315 nim->qto_xyz.m[2][0] , /* [Which is my own] */
316 nim->qto_xyz.m[2][1] , /* [damn fault!!!!!] */
317 nim->qto_xyz.m[2][2] ,
318 nim->qto_xyz.m[2][3] ) ;
319
320 orixyz = THD_matrix_to_orientation( R ) ; /* compute orientation codes */
321
322 // [PT,DRG: Mar 8, 2019] Done below.
323 //iview = NIFTI_code_to_view(nim->qform_code, NULL);
324
325 /* load the offsets and the grid spacings */
326
327 if (ORIENT_xyz[orixyz.ijk[0]] == 'z' ) {
328 orgx = nim->qto_xyz.m[ORIENT_xyzint[orixyz.ijk[0]] - 1][3] ;
329 } else {
330 orgx = - nim->qto_xyz.m[ORIENT_xyzint[orixyz.ijk[0]] - 1][3] ;
331 }
332
333 if (ORIENT_xyz[orixyz.ijk[1]] == 'z' ) {
334 orgy = nim->qto_xyz.m[ORIENT_xyzint[orixyz.ijk[1]] - 1][3] ;
335 } else {
336 orgy = - nim->qto_xyz.m[ORIENT_xyzint[orixyz.ijk[1]] - 1][3] ;
337 }
338
339 if (ORIENT_xyz[orixyz.ijk[2]] == 'z' ) {
340 orgz = nim->qto_xyz.m[ORIENT_xyzint[orixyz.ijk[2]] - 1][3] ;
341 } else {
342 orgz = - nim->qto_xyz.m[ORIENT_xyzint[orixyz.ijk[2]] - 1][3] ;
343 }
344
345
346 LOAD_FVEC3( orgxyz , orgx ,
347 orgy ,
348 orgz ) ;
349 #if 0
350 LOAD_FVEC3( orgxyz , -nim->qto_xyz.m[0][3] , /* again, negate */
351 -nim->qto_xyz.m[1][3] , /* x and y coords */
352 nim->qto_xyz.m[2][3] ) ;
353 #endif
354
355 /* AFNI space units are always mm */
356
357 if( nim->xyz_units == NIFTI_UNITS_METER ){
358 nim->dx *= 1000.0 ; nim->dy *= 1000.0 ; nim->dz *= 1000.0 ;
359 } else if( nim->xyz_units == NIFTI_UNITS_MICRON ){
360 nim->dx *= 0.001 ; nim->dy *= 0.001 ; nim->dz *= 0.001 ;
361 }
362
363 LOAD_FVEC3( dxyz , (ORIENT_sign[orixyz.ijk[0]]=='+') ? nim->dx : -nim->dx ,
364 (ORIENT_sign[orixyz.ijk[1]]=='+') ? nim->dy : -nim->dy ,
365 (ORIENT_sign[orixyz.ijk[2]]=='+') ? nim->dz : -nim->dz ) ;
366
367 } else if (use_sform) {
368
369 int orimap[7] = { 6 , 1 , 0 , 2 , 3 , 4 , 5 } ;
370 int oritmp[3] ;
371 float dxtmp, dytmp, dztmp ;
372 float xmax, ymax, zmax ;
373 float orgx, orgy, orgz ;
374 float fig_merit, ang_merit ;
375
376 form_code = nim->sform_code;
377
378 /* convert sform to nifti orientation codes */
379
380 /* n2 10 Jul, 2015 [rickr] */
381 nifti_dmat44_to_orientation(nim->sto_xyz,
382 &oritmp[0], &oritmp[1], &oritmp[2] ) ;
383
384 /* convert nifti orientation codes to AFNI codes and store in vector */
385
386 LOAD_IVEC3( orixyz , orimap[oritmp[0]] ,
387 orimap[oritmp[1]] ,
388 orimap[oritmp[2]] ) ;
389
390 /* assume original view if there's no talairach id present */
391 // [PT,DRG: Mar 8, 2019] Done below.
392 //iview = NIFTI_code_to_view(nim->sform_code, NULL);
393
394 /* load the offsets and the grid spacings */
395
396 if (ORIENT_xyz[orixyz.ijk[0]] == 'z' ) {
397 orgx = nim->sto_xyz.m[ORIENT_xyzint[orixyz.ijk[0]] - 1][3] ;
398 } else {
399 orgx = - nim->sto_xyz.m[ORIENT_xyzint[orixyz.ijk[0]] - 1][3] ;
400 }
401
402 if (ORIENT_xyz[orixyz.ijk[1]] == 'z' ) {
403 orgy = nim->sto_xyz.m[ORIENT_xyzint[orixyz.ijk[1]] - 1][3] ;
404 } else {
405 orgy = - nim->sto_xyz.m[ORIENT_xyzint[orixyz.ijk[1]] - 1][3] ;
406 }
407
408 if (ORIENT_xyz[orixyz.ijk[2]] == 'z' ) {
409 orgz = nim->sto_xyz.m[ORIENT_xyzint[orixyz.ijk[2]] - 1][3] ;
410 } else {
411 orgz = - nim->sto_xyz.m[ORIENT_xyzint[orixyz.ijk[2]] - 1][3] ;
412 }
413
414
415 LOAD_FVEC3( orgxyz , orgx ,
416 orgy ,
417 orgz ) ;
418
419 #if 0
420 LOAD_FVEC3( orgxyz , -nim->sto_xyz.m[0][3] , /* again, negate */
421 -nim->sto_xyz.m[1][3] , /* x and y coords */
422 nim->sto_xyz.m[2][3] ) ;
423 #endif
424
425 #define MAXNUM(a,b) ( (a) > (b) ? (a):(b))
426 #define MAX3(a,b,c) ( (MAXNUM(a,b)) > (MAXNUM(a,c)) ? (MAXNUM(a,b)):(MAXNUM(a,c)))
427 #define MINNUM(a,b) ( (a) < (b) ? (a):(b))
428 #define MIN3(a,b,c) ( (MINNUM(a,b)) < (MINNUM(a,c)) ? (MINNUM(a,b)):(MINNUM(a,c)))
429
430 dxtmp = sqrt ( nim->sto_xyz.m[0][0] * nim->sto_xyz.m[0][0] +
431 nim->sto_xyz.m[1][0] * nim->sto_xyz.m[1][0] +
432 nim->sto_xyz.m[2][0] * nim->sto_xyz.m[2][0] ) ;
433
434 xmax = MAX3(fabs(nim->sto_xyz.m[0][0]),fabs(nim->sto_xyz.m[1][0]),fabs(nim->sto_xyz.m[2][0])) / dxtmp ;
435
436 dytmp = sqrt ( nim->sto_xyz.m[0][1] * nim->sto_xyz.m[0][1] +
437 nim->sto_xyz.m[1][1] * nim->sto_xyz.m[1][1] +
438 nim->sto_xyz.m[2][1] * nim->sto_xyz.m[2][1] ) ;
439
440 ymax = MAX3(fabs(nim->sto_xyz.m[0][1]),fabs(nim->sto_xyz.m[1][1]),fabs(nim->sto_xyz.m[2][1])) / dytmp ;
441
442 dztmp = sqrt ( nim->sto_xyz.m[0][2] * nim->sto_xyz.m[0][2] +
443 nim->sto_xyz.m[1][2] * nim->sto_xyz.m[1][2] +
444 nim->sto_xyz.m[2][2] * nim->sto_xyz.m[2][2] ) ;
445
446 zmax = MAX3(fabs(nim->sto_xyz.m[0][2]),fabs(nim->sto_xyz.m[1][2]),fabs(nim->sto_xyz.m[2][2])) / dztmp ;
447
448 fig_merit = MIN3(xmax,ymax,zmax) ;
449 ang_merit = acos (fig_merit) * 180.0 / 3.141592653 ;
450 #if 0
451 if (fabs(ang_merit) > .01) {
452 WARNING_message (
453 "qform not present in:\n"
454 " '%s'\n"
455 " oblique sform used, and the worst axis is\n"
456 " %f degrees from plumb.\n"
457 " If you are performing spatial transformations on this dset, \n"
458 " or viewing/combining it with volumes of differing obliquity,\n"
459 " you should consider running: \n"
460 " 3dWarp -deoblique \n"
461 " on this and other oblique datasets in the same session.\n"
462 ,pathnew, ang_merit ) ;
463 }
464 #endif
465
466 if( nim->xyz_units == NIFTI_UNITS_METER ){
467 dxtmp *= 1000.0 ; dytmp *= 1000.0 ; dztmp *= 1000.0 ;
468 } else if( nim->xyz_units == NIFTI_UNITS_MICRON ){
469 dxtmp *= 0.001 ; dytmp *= 0.001 ; dztmp *= 0.001 ;
470 }
471
472 LOAD_FVEC3( dxyz , (ORIENT_sign[orixyz.ijk[0]]=='+') ? dxtmp : -dxtmp ,
473 (ORIENT_sign[orixyz.ijk[1]]=='+') ? dytmp : -dytmp ,
474 (ORIENT_sign[orixyz.ijk[2]]=='+') ? dztmp : -dztmp ) ;
475
476 LOAD_MAT44(ijk_to_dicom44, -nim->sto_xyz.m[0][0] , /* negate x and y */
477 -nim->sto_xyz.m[0][1] , /* coefficients, */
478 -nim->sto_xyz.m[0][2] , /* since AFNI works */
479 -nim->sto_xyz.m[0][3] ,
480 -nim->sto_xyz.m[1][0] , /* with RAI coords, */
481 -nim->sto_xyz.m[1][1] , /* but NIFTI uses */
482 -nim->sto_xyz.m[1][2] , /* LPI coordinates. */
483 -nim->sto_xyz.m[1][3] ,
484 nim->sto_xyz.m[2][0] , /* [Which is my own] */
485 nim->sto_xyz.m[2][1] , /* [damn fault!!!!!] */
486 nim->sto_xyz.m[2][2] ,
487 nim->sto_xyz.m[2][3] ) ;
488
489 } else { /* NO SPATIAL XFORM. BAD BAD BAD BAD BAD BAD. */
490
491 float dxtmp, dytmp, dztmp ;
492
493 /* if pixdim data are present, use them in order to set pixel
494 dimensions. otherwise, set the dimensions to 1 unit. */
495
496 dxtmp = ((nim->pixdim[1] > 0) ? nim->pixdim[1] : 1) ;
497 dytmp = ((nim->pixdim[2] > 0) ? nim->pixdim[2] : 1) ;
498 dztmp = ((nim->pixdim[3] > 0) ? nim->pixdim[3] : 1) ;
499
500 if( nim->xyz_units == NIFTI_UNITS_METER ){
501 dxtmp *= 1000.0 ; dytmp *= 1000.0 ; dztmp *= 1000.0 ;
502 } else if( nim->xyz_units == NIFTI_UNITS_MICRON ){
503 dxtmp *= 0.001 ; dytmp *= 0.001 ; dztmp *= 0.001 ;
504 }
505
506 /* set orientation to LPI by default */
507
508 LOAD_IVEC3( orixyz , 1 ,
509 2 ,
510 4 ) ;
511
512 LOAD_FVEC3( dxyz , (ORIENT_sign[orixyz.ijk[0]]=='+') ? dxtmp : -dxtmp ,
513 (ORIENT_sign[orixyz.ijk[1]]=='+') ? dytmp : -dytmp ,
514 (ORIENT_sign[orixyz.ijk[2]]=='+') ? dztmp : -dztmp ) ;
515
516 /* no qform/sform should default to orig */
517 /* use VIEW_ORIGINAL_TYPE, not NIFTI_default_view() 5 Feb 2016 [rickr] */
518 // [PT,DRG: Mar 8, 2019] Done below.
519 //iview = VIEW_ORIGINAL_TYPE;
520
521 /* set origin to 0,0,0 */
522
523 LOAD_FVEC3( orgxyz , 0 ,
524 0 ,
525 0 ) ;
526 /* put scaled identity matrix by default */
527 LOAD_MAT44(ijk_to_dicom44, dxtmp, 0.0, 0.0, 0.0,
528 0.0, dytmp, 0.0, 0.0,
529 0.0, 0.0, dztmp, 0.0 );
530 }
531
532
533
534 /*-- make an AFNI dataset! --*/
535
536 dset = EDIT_empty_copy(NULL) ;
537
538 ppp = THD_trailname(pathnew,0) ; /* strip directory */
539 MCW_strncpy( prefix , ppp , THD_MAX_PREFIX ) ; /* to make prefix */
540
541 /* You need to set the path too
542 before, if you loaded ~/tmp/joe.nii the path appeared
543 to be ./joe.nii, troubling in multiple instances.
544 ZSS Dec 2011 */
545 THD_init_diskptr_names( dset->dblk->diskptr ,
546 THD_filepath(pathnew) ,
547 NULL , prefix ,
548 dset->view_type , True );
549
550 nxyz.ijk[0] = nim->nx ; /* grid dimensions */
551 nxyz.ijk[1] = nim->ny ;
552 nxyz.ijk[2] = nim->nz ;
553
554 dset->idcode.str[0] = 'N' ; /* overwrite 1st 3 bytes with something special */
555 dset->idcode.str[1] = 'I' ;
556 dset->idcode.str[2] = 'I' ;
557
558 MCW_hash_idcode( pathnew , dset ) ; /* 06 May 2005 */
559
560 // [PT,DRG: Mar 8, 2019] We assume that the SPACE, whether it
561 // existed upon dset read or not, has been set correctly by this
562 // point.
563 // Also note, using the *generic space*, so TT_N27 and TLRC are
564 // both considered TLRC.
565 if (use_qform) {
566 iview = NIFTI_code_to_view(nim->qform_code, NULL);
567 }
568 else if (use_sform) {
569 iview = NIFTI_code_to_view(nim->sform_code, NULL);
570 }
571 else{ /* NO SPATIAL XFORM. BAD BAD BAD BAD BAD BAD. */
572 iview = VIEW_ORIGINAL_TYPE;
573 }
574
575 // And set the dset header
576 EDIT_dset_items( dset ,
577 ADN_prefix , prefix ,
578 ADN_datum_all , datum ,
579 ADN_nxyz , nxyz ,
580 ADN_xyzdel , dxyz ,
581 ADN_xyzorg , orgxyz ,
582 ADN_view_type , iview ,
583 ADN_xyzorient , orixyz ,
584 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
585 ADN_type , (statcode != 0) ? HEAD_FUNC_TYPE
586 : HEAD_ANAT_TYPE ,
587 ADN_none ) ;
588 // [PT,DRG: Mar 8, 2019] Done below.
589 //ADN_view_type , iview ,
590
591 /* copy transformation matrix to dataset structure */
592 /* moved after setting grid 4 Apr 2014 [rickr,drg] */
593 dset->daxes->ijk_to_dicom_real = ijk_to_dicom44;
594
595 /* not a time dependent dataset */
596
597 if( ntt < 2 ){
598 EDIT_dset_items( dset ,
599 ADN_nvals , nbuc ,
600 ADN_datum_all , datum ,
601 ADN_func_type , (statcode != 0) ? FUNC_BUCK_TYPE
602 : ANAT_BUCK_TYPE ,
603 ADN_none ) ;
604
605 } else { /* is a time dependent dataset */
606
607 if( nim->time_units == NIFTI_UNITS_MSEC ){
608 nim->dt *= 0.001 ;
609 nim->toffset *= 0.001 ;
610 } else if( nim->time_units == NIFTI_UNITS_USEC ){
611 nim->dt *= 1.e-6 ;
612 nim->toffset *= 1.e-6 ;
613 }
614 EDIT_dset_items( dset ,
615 ADN_nvals , ntt ,
616 ADN_ntt , ntt ,
617 ADN_datum_all , datum ,
618 ADN_ttorg , nim->toffset , /* 12 Oct 2007 [rickr] */
619 ADN_ttdel , nim->dt ,
620 ADN_ttdur , 0.0 ,
621 ADN_tunits , UNITS_SEC_TYPE ,
622 ADN_func_type , (statcode != 0) ? FUNC_FIM_TYPE
623 : ANAT_EPI_TYPE ,
624 ADN_none ) ;
625
626 /* if present, add stuff about the slice-timing offsets */
627
628 if( nim->slice_dim == 3 && /* AFNI can only deal with */
629 nim->slice_code > 0 && /* slice timing offsets */
630 nim->slice_duration > 0.0 && /* along the k-axis of */
631 nim->slice_start >= 0 && /* the dataset volume */
632 nim->slice_start < nim->nz &&
633 nim->slice_end > nim->slice_start &&
634 nim->slice_end < nim->nz ){
635
636 float *toff=(float *)calloc(sizeof(float),nim->nz) , tsl ;
637 int kk ;
638
639 if( nim->time_units == NIFTI_UNITS_MSEC ) nim->slice_duration *= 0.001;
640 else if( nim->time_units == NIFTI_UNITS_USEC ) nim->slice_duration *= 1.e-6;
641
642 /* set up slice time offsets in the divers orders */
643
644 switch( nim->slice_code ){
645 case NIFTI_SLICE_SEQ_INC:
646 tsl = 0.0 ;
647 for( kk=nim->slice_start ; kk <= nim->slice_end ; kk++ ){
648 toff[kk] = tsl ; tsl += nim->slice_duration ;
649 }
650 break ;
651 case NIFTI_SLICE_SEQ_DEC:
652 tsl = 0.0 ;
653 for( kk=nim->slice_end ; kk >= nim->slice_end ; kk-- ){
654 toff[kk] = tsl ; tsl += nim->slice_duration ;
655 }
656 break ;
657 case NIFTI_SLICE_ALT_INC:
658 tsl = 0.0 ;
659 for( kk=nim->slice_start ; kk <= nim->slice_end ; kk+=2 ){
660 toff[kk] = tsl ; tsl += nim->slice_duration ;
661 }
662 for( kk=nim->slice_start+1 ; kk <= nim->slice_end ; kk+=2 ){
663 toff[kk] = tsl ; tsl += nim->slice_duration ;
664 }
665 break ;
666 case NIFTI_SLICE_ALT_INC2:
667 tsl = 0.0 ;
668 for( kk=nim->slice_start+1 ; kk <= nim->slice_end ; kk+=2 ){
669 toff[kk] = tsl ; tsl += nim->slice_duration ;
670 }
671 for( kk=nim->slice_start ; kk <= nim->slice_end ; kk+=2 ){
672 toff[kk] = tsl ; tsl += nim->slice_duration ;
673 }
674 break ;
675 case NIFTI_SLICE_ALT_DEC:
676 tsl = 0.0 ;
677 for( kk=nim->slice_end ; kk >= nim->slice_start ; kk-=2 ){
678 toff[kk] = tsl ; tsl += nim->slice_duration ;
679 }
680 for( kk=nim->slice_end-1 ; kk >= nim->slice_start ; kk-=2 ){
681 toff[kk] = tsl ; tsl += nim->slice_duration ;
682 }
683 break ;
684 case NIFTI_SLICE_ALT_DEC2:
685 tsl = 0.0 ;
686 for( kk=nim->slice_end-1 ; kk >= nim->slice_start ; kk-=2 ){
687 toff[kk] = tsl ; tsl += nim->slice_duration ;
688 }
689 for( kk=nim->slice_end ; kk >= nim->slice_start ; kk-=2 ){
690 toff[kk] = tsl ; tsl += nim->slice_duration ;
691 }
692 break ;
693 }
694
695 /* not ready for 64-bits, cast nz as (int) [17 Jul 2019 rickr] */
696 if( nim->nz >> 31 ) {
697 static int nnnz=0;
698 if( nnnz == 0 ) {
699 WARNING_message("have NIFTI nz > 2^31, chaos to ensue, enjoy ...");
700 nnnz = 1;
701 }
702 }
703
704 EDIT_dset_items( dset ,
705 ADN_nsl , (int)nim->nz ,
706 ADN_zorg_sl , orgxyz.xyz[2] ,
707 ADN_dz_sl , dxyz.xyz[2] ,
708 ADN_toff_sl , toff ,
709 ADN_none ) ;
710
711 free(toff) ;
712
713 } /* end of slice timing stuff */
714
715 } /* end of 3D+time dataset stuff */
716
717 /* if scalars are attached (and we did not xform data), set them */
718 /* note: this MUST be after EDIT_dset_items(ADN_datum_all), as it
719 will clear any existing values 8 Mar 2021 [rickr] */
720 if( ! xform_data && nim->scl_slope != 0.0 && nim->scl_slope != 1.0) {
721 for( ibr=0 ; ibr < nvals ; ibr++ )
722 DBLK_BRICK_FACTOR(dset->dblk, ibr) = nim->scl_slope;
723 }
724
725
726 /* set atlas space based on NIFTI s/qform code */
727 NIFTI_code_to_space(form_code,dset);
728
729 /* add statistics, if present */
730
731 if( statcode != 0 ){
732 for( ibr=0 ; ibr < nvals ; ibr++ )
733 EDIT_STATAUX4(dset,ibr,statcode,nim->intent_p1,nim->intent_p2,nim->intent_p3,0) ;
734 }
735
736 /*-- flag to read data from disk using NIFTI functions --*/
737
738 dset->dblk->diskptr->storage_mode = STORAGE_BY_NIFTI ;
739 strcpy( dset->dblk->diskptr->brick_name , pathnew ) ;
740 dset->dblk->diskptr->byte_order = nim->byteorder ;
741
742 #if 0
743 for( ibr=0 ; ibr < nvals ; ibr++ ){ /* make sub-brick labels */
744 sprintf(prefix,"%s[%d]",tname,ibr) ;
745 EDIT_BRICK_LABEL( dset , ibr , prefix ) ;
746 }
747 #endif
748
749 /** 10 May 2005: see if there is an AFNI extension;
750 if so, load attributes from it and
751 then edit the dataset appropriately **/
752 /** 12 Nov 2021: moved to new thd_nifti_process_afni_ext [rickr] **/
753 (void)THD_nifti_process_afni_ext(dset, nim, pathnew);
754
755 /* return unpopulated dataset */
756
757
758 // [PT,DRG: Mar 8, 2019] ... go back and fix case of sform_code==2
759 // from above: The fix is in!
760 if (form_code==2) {
761 iview = NIFTI_code_to_view(form_code,
762 THD_get_generic_space(dset));
763
764 /* duplicating here how Rick did this in 3drefit */
765 dset->view_type = iview ;
766 THD_init_diskptr_names( dset->dblk->diskptr ,
767 NULL , NULL , NULL , iview , True ) ;
768 }
769
770 nifti_image_free(nim) ; KILL_pathnew ; RETURN(dset) ;
771 }
772
773 /* If there is an AFNI extension in the NIFTI dataset, process it.
774 * (moved from THD_open_nifti)
775 *
776 * Modify logic to return on the "failed" tests, rather than having a
777 * set of nested if's to proceed. [15 Nov 2021 rickr]
778 */
THD_nifti_process_afni_ext(THD_3dim_dataset * dset,nifti_image * nim,char * pathnew)779 static int THD_nifti_process_afni_ext(THD_3dim_dataset * dset,
780 nifti_image * nim, char * pathnew)
781 {
782 NI_stream ns ;
783 NI_group *ngr , *nngr ;
784 void *nini ;
785 char *buf , *rhs , *cpt ;
786 int nbuf, ee ; /* buf size, extension index */
787
788 ENTRY("THD_nifti_process_afni_ext") ;
789
790 /* scan extension list to find the first AFNI extension */
791
792 for( ee=0 ; ee < nim->num_ext ; ee++ )
793 if( nim->ext_list[ee].ecode == NIFTI_ECODE_AFNI &&
794 nim->ext_list[ee].esize > 32 &&
795 nim->ext_list[ee].edata != NULL ) break ;
796
797 /* if no AFNI extension, we are done */
798 if( ee >= nim->num_ext ) RETURN(0);
799
800 /* point to extension buffer */
801 buf = nim->ext_list[ee].edata ;
802 nbuf = nim->ext_list[ee].esize - 8 ;
803
804 /* if have data, it's long enough, and starts properly, then ... */
805
806 /* if insufficent data, fail */
807 if( buf == NULL || nbuf <= 32 || strncmp(buf,"<?xml",5) != 0 )
808 RETURN(1);
809
810 /* be sure buf is terminated, and find XML prolog close */
811 if( buf[nbuf-1] != '\0' ) buf[nbuf-1] = '\0' ;
812 cpt = strstr(buf,"?>") ;
813
814 /* if no XML prolog close, fail */
815 if( cpt == NULL )
816 RETURN(1);
817
818 /* read contents of NIML stream */
819 ns = NI_stream_open( "str:" , "r" ) ;
820 NI_stream_setbuf( ns , cpt+2 ) ; /* start just after prolog */
821 nini = NI_read_element(ns,1) ; /* get root element */
822 NI_stream_close(ns) ;
823
824 /* if not a group, fail */
825 if( NI_element_type(nini) != NI_GROUP_TYPE )
826 RETURN(0);
827
828 /* have nngr point to AFNI_attributes */
829 ngr = (NI_group *)nini ;
830 if( strcmp(ngr->name,"AFNI_attributes") == 0 ){ /* root is OK */
831 nngr = ngr ;
832 } else { /* search in group for proper element */
833 int nn ; void **nnini ;
834 nn = NI_search_group_deep( ngr , "AFNI_attributes" , &nnini ) ;
835 if( nn <= 0 ) nngr = NULL ;
836 else { nngr = (NI_group *)nnini[0]; NI_free(nnini); }
837 }
838
839 /* if such a group is found, process it */
840 if( NI_element_type(nngr) == NI_GROUP_TYPE ) {
841
842 rhs = NI_get_attribute( nngr , "self_idcode" ) ;
843 if( rhs == NULL )
844 rhs = NI_get_attribute( nngr , "AFNI_idcode" ) ;
845 if( rhs != NULL ) /* set dataset ID code from XML attribute */
846 MCW_strncpy( dset->idcode.str , rhs , MCW_IDSIZE ) ;
847
848
849 /***** This is an example of what might go into consistency check. */
850
851 rhs = NI_get_attribute( nngr , "NIfTI_nums" ) ; /* check if */
852 if( rhs != NULL ){ /* dataset dimensions */
853 char buf[128] ; /* were altered */
854 /* %ld fails on older systems [17 Jul 2019 rickr] */
855 sprintf(buf,"%" PRId64 ",%" PRId64 ",%" PRId64 ",%" PRId64
856 ",%" PRId64 ",%d" , /* 12 May 2005 */
857 nim->nx, nim->ny, nim->nz,
858 nim->nt, nim->nu, nim->datatype );
859 if( strcmp(buf,rhs) != 0 ){
860 static int nnn=0 ;
861 if(nnn==0){fprintf(stderr,"\n"); nnn=1;}
862 fprintf(stderr,
863 "** WARNING: NIfTI file %s dimensions altered since "
864 "AFNI extension was added\n",pathnew ) ;
865 /* fprintf(stderr, " buf=%s\n", buf); */
866 }
867 }
868
869 /* the main point: apply AFNI extension attributes to dset */
870 THD_dblkatr_from_niml( nngr , dset->dblk ); /* load attributes */
871 THD_datablock_apply_atr( dset ) ; /* apply to dataset struct */
872 }
873
874 NI_free_element( ngr ) ; /* get rid of the root element */
875
876 RETURN(0);
877 }
878
879 /* n2 10 Jul, 2015 [rickr] */
copy_ints_as_i64(int * ivals,int nvals)880 int64_t * copy_ints_as_i64(int * ivals, int nvals)
881 {
882 int64_t * i64;
883 int c;
884
885 i64 = (int64_t *)malloc(nvals * sizeof(int64_t));
886 if( ! i64 ) {
887 fprintf(stderr,"** CIA64: failed to alloc %d int64_t's\n", nvals);
888 return NULL;
889 }
890 for( c=0; c<nvals; c++ ) i64[c] = ivals[c];
891
892 return i64;
893 }
894
895
896 /*-----------------------------------------------------------------
897 Load a NIFTI dataset's data into memory
898 (called from THD_load_datablock in thd_loaddblk.c)
899 - RWC: modified 07 Apr 2005 to read data bricks via nifti_io.c
900 'NBL' functions, rather than directly from disk
901 -------------------------------------------------------------------*/
902
THD_load_nifti(THD_datablock * dblk)903 void THD_load_nifti( THD_datablock *dblk )
904 {
905 THD_diskptr *dkptr ;
906 int nx,ny,nz,nxy,nxyz,nxyzv , nerr=0,ibr,nv, nslice ;
907 int datum, need_copy=0 ;
908 int scale_data=0 ;
909 void *ptr ;
910 nifti_image *nim ;
911 nifti_brick_list NBL ; /* holds the data read from disk */
912
913 ENTRY("THD_load_nifti") ;
914
915 /*-- open and read input [these errors should never occur] --*/
916
917 if( !ISVALID_DATABLOCK(dblk) ||
918 dblk->diskptr->storage_mode != STORAGE_BY_NIFTI ||
919 dblk->brick == NULL ) EXRETURN ;
920
921 dkptr = dblk->diskptr ;
922
923 /* purge any existing bricks [10 Mar 2014] */
924
925 STATUS("purging existing data bricks (if any)") ;
926 THD_purge_datablock(dblk,DATABLOCK_MEM_ANY) ;
927
928 STATUS("calling nifti_image_read_bricks") ;
929 NBL.nbricks = 0 ;
930
931 if( ! DBLK_IS_MASTERED(dblk) ) /* allow mastering 14 Apr 2006 [rickr] */
932 nim = nifti_image_read_bricks( dkptr->brick_name, 0,NULL , &NBL ) ;
933 else {
934 /* n2 10 Jul, 2015 [rickr] */
935 /* convert master_ival to an array of int64_t */
936 int64_t * i64_vals = copy_ints_as_i64(dblk->master_ival, dblk->nvals);
937 nim = nifti_image_read_bricks( dkptr->brick_name, dblk->nvals,
938 i64_vals, &NBL ) ;
939 }
940
941 if( nim == NULL || NBL.nbricks <= 0 ) EXRETURN ;
942
943 datum = DBLK_BRICK_TYPE(dblk,0) ; /* destination data type */
944
945 /*-- determine if we need to copy the data from the
946 bricks as loaded above because of a type conversion --*/
947
948 switch( nim->datatype ){
949 case DT_INT16:
950 case DT_UINT8: need_copy = (datum == MRI_float) ; break ;
951
952 case DT_FLOAT32:
953 case DT_COMPLEX64:
954 case DT_RGB24: need_copy = 0 ; break ;
955
956 case DT_INT8: /* these are the cases where AFNI can't */
957 case DT_UINT16: /* directly handle the NIFTI datatype, */
958 case DT_INT32: /* so we'll convert them to floats. */
959 case DT_UINT32:
960 case DT_FLOAT64: need_copy = 1 ; break ;
961 #if 0
962 case DT_COMPLEX128: need_copy = 1 ; break ;
963 #endif
964 }
965
966 /*-- various dimensions --*/
967
968 nx = dkptr->dimsizes[0] ;
969 ny = dkptr->dimsizes[1] ; nxy = nx * ny ;
970 nz = dkptr->dimsizes[2] ; nxyz = nxy * nz ;
971 nv = dkptr->nvals ; if( nv > NBL.nbricks ) nv = NBL.nbricks ;
972 nxyzv = nxyz * nv ; nslice = nz*nv ;
973
974 dblk->malloc_type = DATABLOCK_MEM_MALLOC ;
975
976 /*------ don't need to copy data ==> just copy pointers from NBL ------*/
977
978 if( !need_copy ){
979
980 STATUS("copying brick pointers directly") ;
981 for( ibr=0 ; ibr < nv ; ibr++ ){
982 mri_fix_data_pointer( NBL.bricks[ibr] ,DBLK_BRICK(dblk,ibr) ) ;
983 NBL.bricks[ibr] = NULL ; /* so it won't be deleted later */
984
985 if( DBLK_BRICK_TYPE(dblk,ibr) == MRI_float ){
986 STATUS("doing floatscan") ;
987 nerr += thd_floatscan( DBLK_BRICK_NVOX(dblk,ibr) ,
988 DBLK_ARRAY(dblk,ibr) ) ;
989 } else if( DBLK_BRICK_TYPE(dblk,ibr) == MRI_complex ){
990 STATUS("doing complexscan") ;
991 nerr += thd_complexscan( DBLK_BRICK_NVOX(dblk,ibr) ,
992 DBLK_ARRAY(dblk,ibr) ) ;
993 }
994 }
995 if( nerr > 0 ) WARNING_message("file %s: corrected %d float errors\n",
996 dkptr->brick_name , nerr ) ;
997
998 } else { /*---------- need to copy data ==> do some more work -----------*/
999
1000 register int ii ; void *nbuf ;
1001
1002 STATUS("converting input bricks to floats") ;
1003 for( ibr=0 ; ibr < nv ; ibr++ ){
1004
1005 if( DBLK_ARRAY(dblk,ibr) == NULL ){ /* make space */
1006 ptr = AFMALL(void, DBLK_BRICK_BYTES(dblk,ibr) ) ; /* for this */
1007 if( ptr == NULL ) ERROR_message("malloc fails for NIfTI sub-brick #%d",ibr) ;
1008 mri_fix_data_pointer( ptr , DBLK_BRICK(dblk,ibr) ) ; /* sub-brick! */
1009 }
1010 ptr = DBLK_ARRAY(dblk,ibr) ; if( ptr == NULL ) break ; /* bad news!! */
1011
1012 nbuf = NBL.bricks[ibr] ; /* data as read from NIfTI file */
1013
1014 /* macro to convert data from type "ityp" in nbuf to float in dataset */
1015
1016 #undef CPF
1017 #define CPF(ityp) do{ ityp *sar = (ityp *)nbuf ; float *far = (float *)ptr ; \
1018 for( ii=0 ; ii < nxyz ; ii++ ) far[ii] = (float)sar[ii]; \
1019 } while(0)
1020
1021 /* load from nbuf into brick array (will be float or complex) */
1022
1023 STATUS(" converting sub-brick") ;
1024
1025 switch( nim->datatype ){
1026 case DT_UINT8: CPF(unsigned char) ; break ;
1027 case DT_INT8: CPF(signed char) ; break ;
1028 case DT_INT16: CPF(signed short) ; break ;
1029 case DT_UINT16: CPF(unsigned short) ; break ;
1030 case DT_INT32: CPF(signed int) ; break ;
1031 case DT_UINT32: CPF(unsigned int) ; break ;
1032 case DT_FLOAT64: /* added floatscan 2 Dec, 2014 [rickr] */
1033 { CPF(double) ; thd_floatscan(nxyz, (float *)ptr) ; break ; }
1034 #if 0
1035 case DT_COMPLEX128: break ;
1036 #endif
1037 }
1038
1039 STATUS(" free-ing NIfTI volume") ;
1040
1041 free(NBL.bricks[ibr]) ; NBL.bricks[ibr] = NULL ;
1042 }
1043 }
1044
1045 STATUS("free-ing NBL") ;
1046
1047 nifti_free_NBL( &NBL ) ; /* done with this */
1048
1049 /*-- scale results? ---*/
1050
1051 /* errors for !isfinite() have been printed */
1052
1053 /* do not scale if either slope or inter is zero
1054 * slope==0 : treat as unset (rather than constant data)
1055 * inter==0 : use slope as brick_fac */
1056
1057 /* Any NIFTI scalar previously implied conversion to float.
1058 * No longer scale if inter==0 && slope!=0 (set brick_fac).
1059 * Thanks to C Caballero and S Moia for reporting this.
1060 * 26 Jan 2021 [rickr] */
1061 scale_data = isfinite(nim->scl_slope) && isfinite(nim->scl_inter)
1062 && (nim->scl_slope != 0.0) && (nim->scl_inter != 0.0) ;
1063
1064 if( scale_data ){
1065 STATUS("scaling sub-bricks") ;
1066 for( ibr=0 ; ibr < nv ; ibr++ ){
1067 if( DBLK_BRICK_TYPE(dblk,ibr) == MRI_float ){
1068 float *far = (float *) DBLK_ARRAY(dblk,ibr) ; int ii ;
1069 for( ii=0 ; ii < nxyz ; ii++ ){
1070 far[ii] = nim->scl_slope * far[ii] + nim->scl_inter ;
1071 if( !isfinite(far[ii]) ) far[ii] = 0.0f ;
1072 }
1073 } else if( DBLK_BRICK_TYPE(dblk,ibr) == MRI_complex ){
1074 complex *car = (complex *) DBLK_ARRAY(dblk,ibr) ; int ii ;
1075 for( ii=0 ; ii < nxyz ; ii++ ){
1076 car[ii].r = nim->scl_slope * car[ii].r + nim->scl_inter ;
1077 car[ii].i = nim->scl_slope * car[ii].i + nim->scl_inter ;
1078 }
1079 }
1080 }
1081 } else if ( isfinite(nim->scl_slope) && nim->scl_slope != 0.0 ) {
1082
1083 /* actually pass slope as brick_fac 26 Jan 2021 [rickr] */
1084 for( ibr=0 ; ibr < nv ; ibr++ )
1085 DBLK_BRICK_FACTOR(dblk, ibr) = nim->scl_slope;
1086 }
1087
1088 /* rcr - replace this with DBLK_IS_RANGE_MASTERED to also check for
1089 * csv list
1090 * - then call a new parent function to THD_apply_master_subrange
1091 */
1092 if( DBLK_IS_MASTER_SUBRANGED(dblk) )
1093 THD_apply_master_subrange(dblk) ;
1094
1095 /*-- throw away the trash and return --*/
1096
1097 nifti_image_free(nim) ; EXRETURN ;
1098 }
1099
1100
1101 /* set atlas space based on NIFTI sform code */
NIFTI_code_to_space(int code,THD_3dim_dataset * dset)1102 static void NIFTI_code_to_space(int code,THD_3dim_dataset *dset)
1103 {
1104 switch(code) {
1105 case NIFTI_XFORM_TALAIRACH:
1106 MCW_strncpy(dset->atlas_space, "TLRC", THD_MAX_NAME);
1107 break;
1108 case NIFTI_XFORM_MNI_152:
1109 MCW_strncpy(dset->atlas_space, "MNI", THD_MAX_NAME);
1110 break;
1111 default:
1112 THD_get_space(dset);
1113 }
1114 }
1115
1116 /* return dataset view code based on NIFTI sform/qform code */
1117 /* code for +tlrc, +orig */
NIFTI_code_to_view(int code,char * atlas_space)1118 static int NIFTI_code_to_view(int code, char *atlas_space)
1119 {
1120 int iview;
1121
1122 /*
1123 SFORM SPACE VIEW
1124 ----- ----- ----
1125 1 (or 0) ORIG orig
1126 4 MNI tlrc
1127 3 TLRC tlrc
1128 ?2 (->1) ACPC acpc
1129 5 A_TEMP tlrc
1130
1131 */
1132
1133 ENTRY("NIFTI_code_to_view");
1134
1135 //INFO_message("atlas_space = %s", atlas_space);
1136 //INFO_message("code = %d", code);
1137
1138 /* only two standard templates now defined */
1139 switch(code) {
1140 case NIFTI_XFORM_TALAIRACH: /* TLRC space -> tlrc view */
1141 // sform_code = 3
1142 iview = VIEW_TALAIRACH_TYPE;
1143 break;
1144 case NIFTI_XFORM_MNI_152: /* MNI space -> tlrc 'view' */
1145 // sform_code = 4
1146 iview = VIEW_TALAIRACH_TYPE;
1147 break;
1148 case NIFTI_XFORM_SCANNER_ANAT: /* no code set or scanner -> orig 'view' */
1149 // sform_code = 1
1150 case NIFTI_XFORM_UNKNOWN:
1151 // sform_code = 0
1152 iview = VIEW_ORIGINAL_TYPE;
1153 break;
1154 case NIFTI_XFORM_TEMPLATE_OTHER: // [PT,DRG: Mar 8, 2019] For sform_code=5
1155 // sform_code = 5 !!
1156 iview = VIEW_TALAIRACH_TYPE;
1157 break;
1158 case NIFTI_XFORM_ALIGNED_ANAT: /* aligned to something... */
1159 // sform_code = 2 (booo)
1160 if( !atlas_space ) {
1161 // *this* special scenario is to allow use of this function
1162 // before checking for AFNI extension; it can be called with
1163 // a 'NULL' atlas_space, which won't affect any
1164 // sform_code!=2, but will allow a dummy value for
1165 // sform_code=2 to exist until it is set "properly"
1166 iview = NIFTI_default_view(); // should just be a dummy case
1167 }
1168 else if ( !strcmp(atlas_space, "ORIG") )
1169 iview = VIEW_ORIGINAL_TYPE;
1170 else if ( !strcmp(atlas_space, "ACPC") )
1171 iview = VIEW_ACPCALIGNED_TYPE;
1172 else if ( strlen(atlas_space) )
1173 iview = VIEW_TALAIRACH_TYPE;
1174 else { // no space!
1175 WARNING_message("Hmm, {s,q}form_code = 2, but no space set??\n"
1176 "-> rollin' the dice here!");
1177 iview = NIFTI_default_view();
1178 }
1179 break;
1180 default: /* or something else we don't know
1181 about yet (higher form codes)*/
1182 // sform_code > 5 or sform_code<0
1183 WARNING_message("I don't think this sform code '%d' is allowed!", code);
1184 iview = NIFTI_default_view();
1185 }
1186 RETURN(iview);
1187 }
1188
1189 /* get default view from environment variable */
NIFTI_default_view()1190 static int NIFTI_default_view()
1191 {
1192 char *ppp;
1193 int iview = VIEW_TALAIRACH_TYPE; /* default view if not otherwise set */
1194
1195 ENTRY("NIFTI_default_view");
1196 ppp = my_getenv("AFNI_NIFTI_VIEW");
1197 if(ppp){
1198 if(strcasecmp(ppp, "TLRC")==0)
1199 iview = VIEW_TALAIRACH_TYPE;
1200 else if (strcasecmp(ppp,"ORIG")==0) {
1201 iview = VIEW_ORIGINAL_TYPE;
1202 }
1203 else if (strcasecmp(ppp,"ACPC")==0) {
1204 iview = VIEW_ACPCALIGNED_TYPE;
1205 }
1206 }
1207 RETURN(iview);
1208 }
1209