1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College of
3 Wisconsin, 1994-2000, and are released under the Gnu General Public License,
4 Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8 #include "thd.h"
9 /*----------------------------------------------------------------
10 Given a directory name, and a header filename, create a
11 datablock that corresponds (return NULL if impossible).
12 ------------------------------------------------------------------*/
13
14 static int native_order = -1 ;
15 static int no_mmap = -1 ;
16 static int no_ordwarn = -1 ;
17
THD_init_one_datablock(char * dirname,char * headname)18 THD_datablock * THD_init_one_datablock( char *dirname , char *headname )
19 {
20 THD_datablock *dblk ;
21 THD_diskptr *dkptr ;
22 int ii ;
23 char prefix[THD_MAX_NAME] = "\0" ;
24 int default_order ; /* 21 Jun 2000 */
25
26 ENTRY("THD_init_one_datablock") ;
27
28 /*-- sanity checks --*/
29
30 if( dirname == NULL || strlen(dirname) == 0 ||
31 headname == NULL || strlen(headname) == 0 ) RETURN( NULL ) ;
32
33 FILENAME_TO_PREFIX(headname,prefix) ;
34 if( strlen(prefix) == 0 ||
35 strstr(headname,DATASET_HEADER_SUFFIX) == NULL ) RETURN( NULL ) ;
36
37 /*-- byte ordering stuff --*/
38
39 if( native_order < 0 ) native_order = mri_short_order() ;
40
41 no_mmap = AFNI_yesenv("AFNI_NOMMAP") ;
42 no_ordwarn = AFNI_yesenv("AFNI_NO_BYTEORDER_WARNING") ;
43
44 { char *hh = getenv("AFNI_BYTEORDER_INPUT") ; /* 21 Jun 2000 */
45 default_order = native_order ;
46 if( hh != NULL ){
47 if( strncmp(hh,LSB_FIRST_STRING,ORDER_LEN) == 0 )
48 default_order = LSB_FIRST ;
49 else if( strncmp(hh,MSB_FIRST_STRING,ORDER_LEN) == 0 )
50 default_order = MSB_FIRST ;
51 }
52 }
53
54 #if 1
55
56 dblk = EDIT_empty_datablock() ; /* 11 Mar 2005 -- the new way */
57 dkptr = dblk->diskptr ;
58
59 #else
60 /*-- create output datablock (the old way) --*/
61
62 dblk = myRwcNew( THD_datablock ) ;
63 dblk->type = DATABLOCK_TYPE ;
64 dblk->brick = NULL ; /* will be filled in below */
65 dblk->brick_bytes = NULL ; /* ditto */
66 dblk->brick_fac = NULL ; /* ditto */
67 dblk->total_bytes = 0 ; /* ditto */
68 dblk->malloc_type = DATABLOCK_MEM_UNDEFINED ;
69 dblk->parent = NULL ;
70
71 dblk->brick_lab = NULL ; /* 30 Nov 1997 */
72 dblk->brick_keywords = NULL ;
73 dblk->brick_statcode = NULL ;
74 dblk->brick_stataux = NULL ;
75
76 dblk->master_nvals = 0 ; /* 11 Jan 1999 */
77 dblk->master_ival = NULL ;
78 dblk->master_bytes = NULL ;
79
80 dblk->vedim = NULL ; /* 05 Sep 2006 */
81
82 dblk->brick_fdrcurve = NULL ; /* 23 Jan 2008 */
83 dblk->brick_mdfcurve = NULL ; /* 22 Oct 2008 */
84
85 dblk->master_bot = 1.0 ; /* 21 Feb 2001 */
86 dblk->master_top = 0.0 ;
87
88 DBLK_unlock(dblk) ; /* Feb 1998 */
89
90 dblk->shm_idcode[0] = '\0' ; /* 02 May 2003 */
91
92 INIT_KILL(dblk->kl) ;
93
94 dblk->diskptr = dkptr = myRwcNew( THD_diskptr ) ;
95 dkptr->type = DISKPTR_TYPE ;
96 dkptr->storage_mode = STORAGE_UNDEFINED ;
97 #if 0
98 dkptr->byte_order = native_order ; /* 25 April 1998 */
99 #else
100 dkptr->byte_order = default_order; /* 21 June 2000 */
101 #endif
102
103 ADDTO_KILL(dblk->kl,dkptr) ;
104
105 #endif /* end of initializing empty datablock and diskptr */
106
107 /*-- read attributes from disk, store in the datablock --*/
108
109 THD_read_all_atr( headname , dblk ) ;
110
111 /*-- 09 Mar 2005: all the attribute processing is moved away --*/
112
113 ii = THD_datablock_from_atr( dblk, dirname, headname ) ;
114 if( ii == 0 ){
115 THD_delete_datablock( dblk ) ;
116 myRwcFree(dblk) ;
117 RETURN( NULL ) ;
118 }
119
120 #if 0
121 if( PRINT_TRACING ){
122 char str[256] ;
123 sprintf(str,"rank=%d nvals=%d dim[0]=%d dim[1]=%d dim[2]=%d",
124 dkptr->rank , dkptr->nvals ,
125 dkptr->dimsizes[0] , dkptr->dimsizes[1] , dkptr->dimsizes[2] ) ;
126 STATUS(str) ;
127 }
128 #endif
129
130 RETURN( dblk ) ;
131 }
132
133 /*-----*/
134
135 #undef MYHEAD
136 #define MYHEAD ((headname==NULL) ? "UNKNOWN" : headname)
137
138 /*---------------------------------------------------------------------------*/
139 /*! Take the internal attributes and load the datablock struct up.
140 -----------------------------------------------------------------------------*/
141
THD_datablock_from_atr(THD_datablock * dblk,char * dirname,char * headname)142 int THD_datablock_from_atr( THD_datablock *dblk, char *dirname, char *headname )
143 {
144 THD_diskptr *dkptr ;
145 ATR_int *atr_rank , *atr_dimen , *atr_scene , *atr_btype ;
146 ATR_float *atr_flt ;
147 ATR_string *atr_labs ;
148 int ii , view_type , func_type , dset_type ,
149 nx,ny,nz,nvox , nvals , ibr,typ ;
150 RwcBoolean ok ;
151 char prefix[THD_MAX_NAME]="Unknown" ;
152 MRI_IMAGE *qim ;
153 int brick_ccode ;
154 char name[666] ;
155
156 ENTRY("THD_datablock_from_atr") ;
157
158 if( dblk == NULL || dblk->natr <= 0 ) RETURN(0) ; /* bad input */
159
160 dkptr = dblk->diskptr ;
161
162 /*-- get relevant attributes: rank, dimensions, view_type & func_type --*/
163
164 atr_rank = THD_find_int_atr( dblk , ATRNAME_DATASET_RANK ) ;
165 atr_dimen = THD_find_int_atr( dblk , ATRNAME_DATASET_DIMENSIONS ) ;
166 atr_scene = THD_find_int_atr( dblk , ATRNAME_SCENE_TYPE ) ;
167
168 /*-- missing an attribute ==> quit now --*/
169
170 if( atr_rank == NULL || atr_dimen == NULL || atr_scene == NULL ) RETURN(0) ;
171
172 /*-- load type codes from SCENE attribute --*/
173
174 STATUS("loading *_type from SCENE") ;
175
176 view_type = atr_scene->in[0] ;
177 func_type = atr_scene->in[1] ;
178 dset_type = atr_scene->in[2] ;
179
180 /*-- load other values from attributes into relevant places --*/
181
182 ok = True ;
183 nvox = 1 ;
184
185 STATUS("loading from RANK") ;
186
187 dkptr->rank = atr_rank->in[0] ; /* N.B.: rank isn't used much */
188 dkptr->nvals = dblk->nvals = nvals = atr_rank->in[1] ; /* but nvals is used */
189
190 STATUS("loading from DIMENSIONS") ;
191
192 for( ii=0 ; ii < dkptr->rank ; ii++ ){
193 dkptr->dimsizes[ii] = atr_dimen->in[ii] ;
194 ok = ( ok && dkptr->dimsizes[ii] >= 1 ) ;
195 nvox *= dkptr->dimsizes[ii] ;
196 }
197
198 #if 0
199 if( PRINT_TRACING ){
200 char str[256] ;
201 sprintf(str,"rank=%d nvals=%d dim[0]=%d dim[1]=%d dim[2]=%d nvox=%d",
202 dkptr->rank , dkptr->nvals ,
203 dkptr->dimsizes[0] , dkptr->dimsizes[1] , dkptr->dimsizes[2] , nvox ) ;
204 STATUS(str) ;
205 }
206 #endif
207
208 if( !ok || nvals < 1 ||
209 dkptr->rank < THD_MIN_RANK || dkptr->rank > THD_MAX_RANK ){
210 STATUS("bad rank!!??") ;
211 RETURN(0) ;
212 }
213
214 /*-- create the storage filenames --*/
215
216 STATUS("creating storage filenames") ;
217
218 if( headname != NULL && strchr(headname,'+') != NULL ){
219 FILENAME_TO_PREFIX(headname,prefix) ;
220 THD_init_diskptr_names( dkptr, dirname,NULL,prefix , view_type , True ) ;
221 } else {
222 if( headname != NULL ) MCW_strncpy(prefix,headname,THD_MAX_NAME) ;
223 THD_init_diskptr_names( dkptr, dirname,NULL,prefix , view_type , True ) ;
224 }
225
226 /*-- determine if the BRIK file exists --*/
227
228 STATUS("checking if .BRIK file exists") ;
229
230 brick_ccode = COMPRESS_filecode(dkptr->brick_name) ;
231 if (dkptr->storage_mode == STORAGE_UNDEFINED) { /* ZSS: Oct. 2011
232 the next line was being called all the time before */
233 if( brick_ccode != COMPRESS_NOFILE )
234 dkptr->storage_mode = STORAGE_BY_BRICK ; /* a .BRIK file */
235 }
236
237 /*-- if VOLUME_FILENAMES attribute exists, make it so [20 Jun 2002] --*/
238
239 if( headname != NULL && dkptr->storage_mode == STORAGE_UNDEFINED ){
240 atr_labs = THD_find_string_atr(dblk,"VOLUME_FILENAMES") ;
241 if( atr_labs != NULL ){
242 dkptr->storage_mode = STORAGE_BY_VOLUMES ;
243 dblk->malloc_type = DATABLOCK_MEM_MALLOC ;
244 }
245 }
246
247 /*-- now set the memory allocation codes, etc. --*/
248
249 dblk->brick_fac = (float *) RwcMalloc( sizeof(float) * nvals ) ;
250 for( ibr=0 ; ibr < nvals ; ibr++ ) dblk->brick_fac[ibr] = 0.0 ;
251
252 /* scaling factors from short type to float type, if nonzero */
253
254 if( !AFNI_yesenv("AFNI_IGNORE_BRICK_FLTFAC") ){
255 atr_flt = THD_find_float_atr( dblk , ATRNAME_BRICK_FLTFAC ) ;
256 if( atr_flt != NULL ){
257 for( ibr=0 ; ibr < nvals && ibr < atr_flt->nfl ; ibr++ )
258 dblk->brick_fac[ibr] = atr_flt->fl[ibr] ;
259 }
260 }
261
262 /** Now create an empty shell of the "brick" == the data structure
263 that will hold all the voxel data. Note that all datablocks
264 will have a brick, even if they never actually contain data
265 themselves (are only warp-on-demand).
266
267 If the BRICK_TYPES input attribute doesn't exist, then all
268 sub-bricks are shorts. This makes the code work with old-style
269 datasets, which were always made up of shorts.
270 **/
271
272 atr_btype = THD_find_int_atr( dblk , ATRNAME_BRICK_TYPES ) ;
273
274 if( atr_btype == NULL ){
275 THD_init_datablock_brick( dblk , MRI_short , NULL ) ;
276 } else {
277 THD_init_datablock_brick( dblk , atr_btype->nin , atr_btype->in ) ;
278 }
279
280 if( !THD_datum_constant(dblk) ){ /* 15 Sep 2004 */
281 fprintf(stderr,
282 "\n** WARNING: File %s has mixed-type sub-bricks. ", MYHEAD ) ;
283 }
284
285 /* 25 April 1998: check if the byte order is stored inside */
286
287 atr_labs = THD_find_string_atr( dblk , ATRNAME_BYTEORDER ) ;
288 if( atr_labs != NULL && atr_labs->nch > 0 ){
289
290 if( strncmp(atr_labs->ch,LSB_FIRST_STRING,ORDER_LEN) == 0 )
291 dkptr->byte_order = LSB_FIRST ;
292 else if( strncmp(atr_labs->ch,MSB_FIRST_STRING,ORDER_LEN) == 0 )
293 dkptr->byte_order = MSB_FIRST ;
294 else
295 fprintf(stderr,"*** Unknown %s found in dataset %s\n",
296 ATRNAME_BYTEORDER , MYHEAD ) ;
297
298 } else if( !no_ordwarn &&
299 DBLK_BRICK_TYPE(dblk,0) != MRI_byte &&
300 dblk->diskptr->storage_mode == STORAGE_BY_BRICK ){ /* 20 Sep 1999 */
301
302 static int first=1 ;
303 if( first ){
304 fprintf(stderr,
305 "\n*** The situation below can be rectified with program '3drefit -byteorder':\n");
306 first = 0 ;
307 }
308 fprintf(stderr," ** Dataset %s: assuming byteorder %s\n",
309 MYHEAD , BYTE_ORDER_STRING(dkptr->byte_order) ) ;
310 }
311
312 /* if the data is not on disk, the flag remains at DATABLOCK_MEM_UNDEFINED,
313 otherwise the flag says how the memory for the bricks is to be created. */
314
315 if( dkptr->storage_mode == STORAGE_BY_BRICK ){
316 #if MMAP_THRESHOLD > 0
317 dblk->malloc_type = (dblk->total_bytes > MMAP_THRESHOLD)
318 ? DATABLOCK_MEM_MMAP : DATABLOCK_MEM_MALLOC ;
319 DBLK_mmapfix(dblk) ; /* 18 Mar 2005 */
320 #else
321 dblk->malloc_type = DATABLOCK_MEM_MALLOC ;
322 #endif
323
324 /* must be malloc-ed if:
325 data is compressed,
326 data is not in native byte order, or
327 user explicity forbids use of mmap */
328
329 if( brick_ccode >= 0 || dkptr->byte_order != native_order || no_mmap )
330 dblk->malloc_type = DATABLOCK_MEM_MALLOC ;
331 }
332
333 /* 30 Nov 1997: create the labels for sub-bricks */
334
335 THD_init_datablock_labels( dblk ) ;
336
337 atr_labs = THD_find_string_atr( dblk , ATRNAME_BRICK_LABS ) ;
338 if( atr_labs != NULL && atr_labs->nch > 0 ){ /* create labels from attribute */
339 int ipos = -1 , ipold , ngood ;
340
341 for( ibr=0 ; ibr < nvals ; ibr++ ){ /* loop over bricks */
342
343 for( ipold = ipos++ ; /* skip to */
344 ipos < atr_labs->nch && atr_labs->ch[ipos] != '\0' ; /* next \0 */
345 ipos++ ) /* nada */ ; /* or end. */
346
347 ngood = ipos - ipold - 1 ; /* number of good chars */
348 if( ngood > 0 ){
349 RwcFree(dblk->brick_lab[ibr]) ;
350 /* 27 Oct 2011 - increase to 64 */
351 if( ngood > THD_MAX_SBLABEL ) ngood = THD_MAX_SBLABEL;
352 dblk->brick_lab[ibr] = (char *) RwcMalloc(sizeof(char)*(ngood+2)) ;
353 memcpy( dblk->brick_lab[ibr] , atr_labs->ch+(ipold+1) , ngood ) ;
354 dblk->brick_lab[ibr][ngood] = '\0' ;
355
356 }
357
358 if( ipos >= atr_labs->nch ) break ; /* nothing more to do */
359 } /* end of loop over sub-bricks */
360 }
361
362 /* create the keywords for sub-bricks */
363
364 THD_init_datablock_keywords( dblk ) ;
365
366 atr_labs = THD_find_string_atr( dblk , ATRNAME_BRICK_KEYWORDS ) ;
367
368 if( atr_labs != NULL && atr_labs->nch > 0 ){ /* create keywords from attribute */
369 int ipos = -1 , ipold , ngood ;
370
371 for( ibr=0 ; ibr < nvals ; ibr++ ){ /* loop over bricks */
372
373 for( ipold = ipos++ ; /* skip to */
374 ipos < atr_labs->nch && atr_labs->ch[ipos] != '\0' ; /* next \0 */
375 ipos++ ) /* nada */ ; /* or end. */
376
377 ngood = ipos - ipold - 1 ; /* number of good chars */
378 if( ngood > 0 ){
379 RwcFree(dblk->brick_keywords[ibr]) ;
380 dblk->brick_keywords[ibr] = (char *) RwcMalloc(sizeof(char)*(ngood+2)) ;
381 memcpy( dblk->brick_keywords[ibr] , atr_labs->ch+(ipold+1) , ngood ) ;
382 dblk->brick_keywords[ibr][ngood] = '\0' ;
383 }
384
385 if( ipos >= atr_labs->nch ) break ; /* nothing more to do */
386 } /* end of loop over sub-bricks */
387 }
388
389 /* create the auxiliary statistics stuff for each brick, if present */
390
391 atr_labs = THD_find_string_atr( dblk , "BRICK_STATSYM" ) ; /* 01 Jun 2005 */
392 if( atr_labs != NULL && atr_labs->nch > 0 ){
393 NI_str_array *sar ; int scode,np ; float parm[3] ;
394 sar = NI_decode_string_list( atr_labs->ch , ";" ) ;
395 if( sar != NULL && sar->num > 0 ){
396 for( ibr=0 ; ibr < nvals && ibr < sar->num ; ibr++ ){
397 NI_stat_decode( sar->str[ibr] , &scode , parm,parm+1,parm+2 ) ;
398 if( scode >= AFNI_FIRST_STATCODE && scode <= AFNI_LAST_STATCODE ){
399 np = NI_stat_numparam(scode) ;
400 THD_store_datablock_stataux( dblk , ibr,scode,np,parm ) ;
401 }
402 }
403 NI_delete_str_array(sar) ;
404 }
405 } else { /*--- the olde way to get ye brick stataux parameters ---*/
406 atr_flt = THD_find_float_atr( dblk , ATRNAME_BRICK_STATAUX ) ;
407 if( atr_flt != NULL && atr_flt->nfl >= 3 ){
408 int ipos=0 , iv,nv,jv ;
409
410 /* attribute stores all stataux stuff as follows:
411 sub-brick-index statcode no.-of-values value ... value
412 sub-brick-index statcode no.-of-values value ... value, etc. */
413
414 while( ipos <= atr_flt->nfl - 3 ){
415 iv = (int) ( atr_flt->fl[ipos++] ) ; /* which sub-brick */
416 jv = (int) ( atr_flt->fl[ipos++] ) ; /* statcode */
417 nv = (int) ( atr_flt->fl[ipos++] ) ; /* # of values that follow */
418
419 if( nv > atr_flt->nfl - ipos ) nv = atr_flt->nfl - ipos ;
420
421 THD_store_datablock_stataux( dblk , iv , jv , nv , atr_flt->fl + ipos ) ;
422 ipos += nv ;
423 }
424 }
425 }
426 #if 0
427 if( PRINT_TRACING ){
428 char str[256] ;
429 sprintf(str,"rank=%d nvals=%d dim[0]=%d dim[1]=%d dim[2]=%d",
430 dkptr->rank , dkptr->nvals ,
431 dkptr->dimsizes[0] , dkptr->dimsizes[1] , dkptr->dimsizes[2] ) ;
432 STATUS(str) ;
433 }
434 #endif
435
436 /*-- FDR curves [23 Jan 2008] --*/
437
438 for( ibr=0 ; ibr < dblk->nvals ; ibr++ ){
439 sprintf(name,"FDRCURVE_%06d",ibr) ;
440 atr_flt = THD_find_float_atr( dblk , name ) ;
441 if( atr_flt != NULL && atr_flt->nfl > 3 ){
442 int nv = atr_flt->nfl - 2 ; floatvec *fv ;
443 MAKE_floatvec(fv,nv) ;
444 fv->x0 = atr_flt->fl[0] ; fv->dx = atr_flt->fl[1] ;
445 memcpy( fv->ar , atr_flt->fl + 2 , sizeof(float)*nv ) ;
446 if( dblk->brick_fdrcurve == NULL )
447 dblk->brick_fdrcurve = (floatvec **)calloc(sizeof(floatvec *),dblk->nvals);
448 dblk->brick_fdrcurve[ibr] = fv ;
449 }
450 }
451
452 for( ibr=0 ; ibr < dblk->nvals ; ibr++ ){
453 sprintf(name,"MDFCURVE_%06d",ibr) ;
454 atr_flt = THD_find_float_atr( dblk , name ) ;
455 if( atr_flt != NULL && atr_flt->nfl > 3 ){
456 int nv = atr_flt->nfl - 2 ; floatvec *fv ;
457 MAKE_floatvec(fv,nv) ;
458 fv->x0 = atr_flt->fl[0] ; fv->dx = atr_flt->fl[1] ;
459 memcpy( fv->ar , atr_flt->fl + 2 , sizeof(float)*nv ) ;
460 if( dblk->brick_mdfcurve == NULL )
461 dblk->brick_mdfcurve = (floatvec **)calloc(sizeof(floatvec *),dblk->nvals);
462 dblk->brick_mdfcurve[ibr] = fv ;
463 }
464 }
465
466 RETURN(1) ;
467 }
468
469 #if 0
470 /* update daxes structure in dataset header from datablock attributes */
471 int THD_daxes_from_atr( THD_datablock *dblk, THD_dataxes *daxes)
472 {
473 ATR_int *atr_rank , *atr_dimen , *atr_scene , *atr_btype ;
474 ATR_float *atr_flt ;
475 ATR_string *atr_labs ;
476 int ii , view_type , func_type , dset_type ,
477 nx,ny,nz,nvox , nvals , ibr,typ ;
478 RwcBoolean ok ;
479 char prefix[THD_MAX_NAME]="Unknown" ;
480 MRI_IMAGE *qim ;
481 int brick_ccode ;
482
483 ENTRY("THD_daxes_from_atr") ;
484
485 if( dblk == NULL || dblk->natr <= 0 ) RETURN(0) ; /* bad input */
486
487 dkptr = dblk->diskptr ;
488
489 /*-- get relevant attributes: rank, dimensions, view_type & func_type --*/
490
491 atr_rank = THD_find_int_atr( dblk , ATRNAME_DATASET_RANK ) ;
492 atr_dimen = THD_find_int_atr( dblk , ATRNAME_DATASET_DIMENSIONS ) ;
493 atr_scene = THD_find_int_atr( dblk , ATRNAME_SCENE_TYPE ) ;
494
495 /*-- missing an attribute ==> quit now --*/
496
497 if( atr_rank == NULL || atr_dimen == NULL || atr_scene == NULL ) RETURN(0) ;
498
499 /*-- load type codes from SCENE attribute --*/
500
501 STATUS("loading *_type from SCENE") ;
502
503 view_type = atr_scene->in[0] ;
504 func_type = atr_scene->in[1] ;
505 dset_type = atr_scene->in[2] ;
506
507 /*-- load other values from attributes into relevant places --*/
508
509 ok = True ;
510 nvox = 1 ;
511
512 RETURN(1) ;
513 }
514
515 #endif
516
517 /*-------------------------------------------------------------------*/
518 /*! Given a 12 parameter affine transform created a la 3dWarpDrive's
519 WARPDRIVE_MATVEC_INV_000000, turn it into a WARP_DATA type array of 30
520 parameters. See matlab function MATVEC_to_WARP.m for more inspiration.
521
522 ZSS: June 06
523 ---------------------------------------------------------------------*/
524
Matvec_2_WarpData(ATR_float * atr_flo,THD_affine_warp * ww,float * wdv)525 int Matvec_2_WarpData(ATR_float *atr_flo, THD_affine_warp *ww, float *wdv)
526 {
527 int ans=0;
528 mat44 Mfor, Mbac;
529 int k;
530 int dbg = 0;
531
532 ENTRY("Matvec_2_WarpData") ;
533 if (!atr_flo) {
534 fprintf(stderr,"NULL atr_flo!\n");
535 RETURN(ans);
536 }
537 if (atr_flo->nfl != 12) {
538 fprintf(stderr,"atr_flo->nfl != 12\n");
539 RETURN(ans);
540 }
541
542 if (!ww) {
543 fprintf(stderr,"NULL ww\n");
544 RETURN(ans);
545 }
546
547 ww->type = WARP_AFFINE_TYPE;
548 ww->resam_type = 0 ; /* not used */
549 ww->warp.type = MAPPING_LINEAR_TYPE ;
550
551 k = 0;
552 Mfor.m[0][0] = atr_flo->fl[k]; ++k;
553 Mfor.m[0][1] = atr_flo->fl[k]; ++k;
554 Mfor.m[0][2] = atr_flo->fl[k]; ++k;
555 Mfor.m[0][3] = atr_flo->fl[k]; ++k;
556 Mfor.m[1][0] = atr_flo->fl[k]; ++k;
557 Mfor.m[1][1] = atr_flo->fl[k]; ++k;
558 Mfor.m[1][2] = atr_flo->fl[k]; ++k;
559 Mfor.m[1][3] = atr_flo->fl[k]; ++k;
560 Mfor.m[2][0] = atr_flo->fl[k]; ++k;
561 Mfor.m[2][1] = atr_flo->fl[k]; ++k;
562 Mfor.m[2][2] = atr_flo->fl[k]; ++k;
563 Mfor.m[2][3] = atr_flo->fl[k]; ++k;
564 Mfor.m[3][0] = 0.0;
565 Mfor.m[3][1] = 0.0;
566 Mfor.m[3][2] = 0.0;
567 Mfor.m[3][3] = 0.0;
568 if (dbg) {
569 DUMP_MAT44("Mfor:\n",Mfor);
570 }
571
572 /* calculate the backward transform */
573 Mbac = nifti_mat44_inverse(Mfor);
574 if (dbg) {
575 DUMP_MAT44("Mbac:\n",Mbac);
576 }
577
578 if (wdv) {
579 /* Load the forward transform */
580 /* Now load the 30 values of Wd */
581 k=0;
582 wdv[k] = Mfor.m[0][0]; ++k;
583 wdv[k] = Mfor.m[0][1]; ++k;
584 wdv[k] = Mfor.m[0][2]; ++k;
585 wdv[k] = Mfor.m[1][0]; ++k;
586 wdv[k] = Mfor.m[1][1]; ++k;
587 wdv[k] = Mfor.m[1][2]; ++k;
588 wdv[k] = Mfor.m[2][0]; ++k;
589 wdv[k] = Mfor.m[2][1]; ++k;
590 wdv[k] = Mfor.m[2][2]; ++k;
591
592 wdv[k] = Mbac.m[0][0]; ++k;
593 wdv[k] = Mbac.m[0][1]; ++k;
594 wdv[k] = Mbac.m[0][2]; ++k;
595 wdv[k] = Mbac.m[1][0]; ++k;
596 wdv[k] = Mbac.m[1][1]; ++k;
597 wdv[k] = Mbac.m[1][2]; ++k;
598 wdv[k] = Mbac.m[2][0]; ++k;
599 wdv[k] = Mbac.m[2][1]; ++k;
600 wdv[k] = Mbac.m[2][2]; ++k;
601
602 wdv[k] = -Mfor.m[0][3]; ++k;
603 wdv[k] = -Mfor.m[1][3]; ++k;
604 wdv[k] = -Mfor.m[2][3]; ++k;
605
606 wdv[k] = -Mbac.m[0][3]; ++k;
607 wdv[k] = -Mbac.m[1][3]; ++k;
608 wdv[k] = -Mbac.m[2][3]; ++k;
609
610 /* bot and top are filled as to not cause trouble for Talairach bounding box at a minimum from:
611 [-80 ; -80 ; -65] to [80 ; 110 ; 85] */
612 wdv[k] = -80 * 2; ++k; /* x 2 ? be generous, it is free! */
613 wdv[k] = -80 * 2; ++k;
614 wdv[k] = -65 * 2; ++k;
615
616 wdv[k] = 80 * 2; ++k;
617 wdv[k] = 110 * 2; ++k;
618 wdv[k] = 85 * 2; ++k;
619 }
620
621 /* Now load the 30 values into warp */
622 ww->warp.mfor.mat[0][0] = Mfor.m[0][0];
623 ww->warp.mfor.mat[0][1] = Mfor.m[0][1];
624 ww->warp.mfor.mat[0][2] = Mfor.m[0][2];
625 ww->warp.mfor.mat[1][0] = Mfor.m[1][0];
626 ww->warp.mfor.mat[1][1] = Mfor.m[1][1];
627 ww->warp.mfor.mat[1][2] = Mfor.m[1][2];
628 ww->warp.mfor.mat[2][0] = Mfor.m[2][0];
629 ww->warp.mfor.mat[2][1] = Mfor.m[2][1];
630 ww->warp.mfor.mat[2][2] = Mfor.m[2][2];
631
632 ww->warp.mbac.mat[0][0] = Mbac.m[0][0];
633 ww->warp.mbac.mat[0][1] = Mbac.m[0][1];
634 ww->warp.mbac.mat[0][2] = Mbac.m[0][2];
635 ww->warp.mbac.mat[1][0] = Mbac.m[1][0];
636 ww->warp.mbac.mat[1][1] = Mbac.m[1][1];
637 ww->warp.mbac.mat[1][2] = Mbac.m[1][2];
638 ww->warp.mbac.mat[2][0] = Mbac.m[2][0];
639 ww->warp.mbac.mat[2][1] = Mbac.m[2][1];
640 ww->warp.mbac.mat[2][2] = Mbac.m[2][2];
641
642 ww->warp.bvec.xyz[0] = -Mfor.m[0][3];
643 ww->warp.bvec.xyz[1] = -Mfor.m[1][3];
644 ww->warp.bvec.xyz[2] = -Mfor.m[2][3];
645
646 ww->warp.svec.xyz[0] = -Mbac.m[0][3];
647 ww->warp.svec.xyz[1] = -Mbac.m[1][3];
648 ww->warp.svec.xyz[2] = -Mbac.m[2][3];
649
650 /* bot and top are filled as to not cause trouble for Talairach bounding box at a minimum from:
651 [-80 ; -80 ; -65] to [80 ; 110 ; 85] */
652 ww->warp.bot.xyz[0] = -80 * 2; /* x 2 ? be generous, it is free! */
653 ww->warp.bot.xyz[1] = -80 * 2;
654 ww->warp.bot.xyz[2] = -65 * 2;
655
656 ww->warp.top.xyz[0] = 80 * 2;
657 ww->warp.top.xyz[1] = 110 * 2;
658 ww->warp.top.xyz[2] = 85 * 2;
659
660 RETURN(1);
661 }
662
THD_WarpData_From_3dWarpDrive(THD_3dim_dataset * dset,ATR_float * atr_flt)663 int THD_WarpData_From_3dWarpDrive(THD_3dim_dataset *dset, ATR_float *atr_flt)
664 {
665 int dbg = 0;
666
667 ENTRY("THD_WarpData_From_3dWarpDrive");
668
669 if (!dset) {
670 fprintf(stderr,"NULL dset!");
671 RETURN(0);
672 }
673 if (dset->warp) {
674 /* fprintf(stderr,"Warp already there!");
675 Used to return(0) DRG,ZSS:01/03/13*/
676 SINGLE_KILL( dset->kl , dset->warp ) ;
677 dset->warp = NULL;
678 }
679 if (!atr_flt) {
680 fprintf(stderr,"No attribute!");
681 RETURN(0);
682 }
683 if (atr_flt->nfl != 12) {
684 fprintf( stderr,
685 "Number of parameters in TLRC transform is not 12.\n"
686 "I won't float your boat.\n");
687 RETURN(0);
688 }
689 dset->warp = myRwcNew( THD_warp ) ;
690 ADDTO_KILL( dset->kl , dset->warp ) ;
691 {
692 THD_affine_warp *ww = (THD_affine_warp *) dset->warp ;
693 if (!Matvec_2_WarpData(atr_flt, ww, NULL)) {
694 fprintf(stderr,"Failed to create warp!");
695 RETURN(0);
696 }
697 }
698 /* If you have a warp, you must have a warp_parent
699 However, previous versions of @auto_tlrc did not set
700 that, so use some defaults */
701 if( strlen(dset->warp_parent_name) <= 0
702 && ISZERO_IDCODE(dset->warp_parent_idcode)) {
703 if (dbg) fprintf(stderr,"Assigning a dummy warp parent name\n");
704 sprintf(dset->warp_parent_name,"Not_Set");
705 }
706
707 RETURN(1);
708 }
709 /*---------------------------------------------------------------------------*/
710 /* Macros to fetch an attribute named nnn and test if it exists. */
711
712 #define ATR_IS_STR(nnn) ( (atr_str = THD_find_string_atr(blk,nnn)) != NULL )
713 #define ATR_IS_FLT(nnn) ( (atr_flt = THD_find_float_atr (blk,nnn)) != NULL )
714 #define ATR_IS_INT(nnn) ( (atr_int = THD_find_int_atr (blk,nnn)) != NULL )
715
716 /*---------------------------------------------------------------------------*/
717 /*! Apply attributes to modify an existing datablock.
718 Only some attributes have an effect.
719 09 May 2005 -- written to support NIfTI-ization, by allowing
720 attributes to be applied AFTER a dataset is created.
721 -----------------------------------------------------------------------------*/
722
THD_datablock_apply_atr(THD_3dim_dataset * dset)723 void THD_datablock_apply_atr( THD_3dim_dataset *dset )
724 {
725 THD_datablock *blk ;
726 THD_diskptr *dkptr ;
727 THD_dataxes *daxes ;
728 ATR_int *atr_int = NULL ;
729 ATR_float *atr_flt = NULL;
730 ATR_string *atr_str = NULL ;
731 int ii , view_type , func_type , dset_type , nx,ny,nz,nvox , nvals , ibr,typ ;
732 RwcBoolean ok ;
733 char prefix[THD_MAX_NAME]="Unknown" ;
734 MRI_IMAGE *qim ;
735 int brick_ccode ;
736 char name[666] ;
737
738 ENTRY("THD_datablock_apply_atr") ;
739
740 if( !ISVALID_DSET(dset) ) EXRETURN ; /* bad input */
741
742 blk = dset->dblk ; if( blk == NULL ) EXRETURN ;
743 nvals = blk->nvals ; if( nvals <= 0 ) EXRETURN ;
744 daxes = dset->daxes ; if( daxes == NULL ) EXRETURN ;
745 dkptr = blk->diskptr ;
746
747 /*-- brick labels --*/
748
749 if( ATR_IS_STR(ATRNAME_BRICK_LABS) ){
750 int ipos = -1 , ipold , ngood ;
751
752 STATUS("brick labels") ;
753 if( blk->brick_lab == NULL ) THD_init_datablock_labels( blk ) ;
754
755 for( ibr=0 ; ibr < nvals ; ibr++ ){ /* loop over bricks */
756
757 for( ipold = ipos++ ; /* skip to */
758 ipos < atr_str->nch && atr_str->ch[ipos] != '\0' ; /* next \0 */
759 ipos++ ) /* nada */ ; /* or end. */
760
761 ngood = ipos - ipold - 1 ; /* number of good chars */
762 if( ngood > 0 ){
763 RwcFree(blk->brick_lab[ibr]) ;
764 if( ngood > THD_MAX_SBLABEL ) ngood = THD_MAX_SBLABEL ;
765 blk->brick_lab[ibr] = (char *) RwcMalloc(sizeof(char)*(ngood+2)) ;
766 memcpy( blk->brick_lab[ibr] , atr_str->ch+(ipold+1) , ngood ) ;
767 blk->brick_lab[ibr][ngood] = '\0' ;
768 }
769
770 if( ipos >= atr_str->nch ) break ; /* nothing more to do */
771 } /* end of loop over sub-bricks */
772 }
773
774 /*-- keywords for the dataset itself --*/
775
776 if( ATR_IS_STR(ATRNAME_KEYWORDS) ){
777 STATUS("dataset keywords") ;
778 dset->keywords = RwcNewString( atr_str->ch ) ;
779 }
780
781 /*-- keywords for sub-bricks --*/
782
783 if( ATR_IS_STR(ATRNAME_BRICK_KEYWORDS) ){
784 int ipos = -1 , ipold , ngood ;
785
786 STATUS("brick keywords") ;
787 if( blk->brick_keywords == NULL ) THD_init_datablock_keywords( blk ) ;
788
789 for( ibr=0 ; ibr < nvals ; ibr++ ){ /* loop over bricks */
790
791 for( ipold = ipos++ ; /* skip to */
792 ipos < atr_str->nch && atr_str->ch[ipos] != '\0' ; /* next \0 */
793 ipos++ ) /* nada */ ; /* or end. */
794
795 ngood = ipos - ipold - 1 ; /* number of good chars */
796 if( ngood > 0 ){
797 RwcFree(blk->brick_keywords[ibr]) ;
798 blk->brick_keywords[ibr] = (char *) RwcMalloc(sizeof(char)*(ngood+2)) ;
799 memcpy( blk->brick_keywords[ibr] , atr_str->ch+(ipold+1) , ngood ) ;
800 blk->brick_keywords[ibr][ngood] = '\0' ;
801 }
802
803 if( ipos >= atr_str->nch ) break ; /* nothing more to do */
804 } /* end of loop over sub-bricks */
805 }
806
807 /*-- auxiliary statistics stuff for each brick --*/
808
809 if( ATR_IS_FLT(ATRNAME_BRICK_STATAUX) ){
810 int ipos=0 , iv,nv,jv ;
811
812 STATUS("brick stataux") ;
813
814 /* attribute stores all stataux stuff as follows:
815 sub-brick-index statcode no.-of-values value ... value
816 sub-brick-index statcode no.-of-values value ... value, etc. */
817
818 while( ipos <= atr_flt->nfl - 3 ){
819 iv = (int) ( atr_flt->fl[ipos++] ) ; /* which sub-brick */
820 jv = (int) ( atr_flt->fl[ipos++] ) ; /* statcode */
821 nv = (int) ( atr_flt->fl[ipos++] ) ; /* # of values that follow */
822
823 if( nv > atr_flt->nfl - ipos ) nv = atr_flt->nfl - ipos ;
824
825 THD_store_datablock_stataux( blk , iv , jv , nv , atr_flt->fl + ipos ) ;
826 ipos += nv ;
827 }
828 }
829
830 /*-- FDR curves [23 Jan 2008] --*/
831
832 for( ibr=0 ; ibr < blk->nvals ; ibr++ ){
833 sprintf(name,"FDRCURVE_%06d",ibr) ;
834 atr_flt = THD_find_float_atr( blk , name ) ;
835 if( atr_flt != NULL && atr_flt->nfl > 3 ){
836 int nv = atr_flt->nfl - 2 ; floatvec *fv ;
837 MAKE_floatvec(fv,nv) ;
838 fv->x0 = atr_flt->fl[0] ; fv->dx = atr_flt->fl[1] ;
839 memcpy( fv->ar , atr_flt->fl + 2 , sizeof(float)*nv ) ;
840 if( blk->brick_fdrcurve == NULL )
841 blk->brick_fdrcurve = (floatvec **)calloc(sizeof(floatvec *),blk->nvals);
842 blk->brick_fdrcurve[ibr] = fv ;
843 }
844 }
845
846 for( ibr=0 ; ibr < blk->nvals ; ibr++ ){
847 sprintf(name,"MDFCURVE_%06d",ibr) ;
848 atr_flt = THD_find_float_atr( blk , name ) ;
849 if( atr_flt != NULL && atr_flt->nfl > 3 ){
850 int nv = atr_flt->nfl - 2 ; floatvec *fv ;
851 MAKE_floatvec(fv,nv) ;
852 fv->x0 = atr_flt->fl[0] ; fv->dx = atr_flt->fl[1] ;
853 memcpy( fv->ar , atr_flt->fl + 2 , sizeof(float)*nv ) ;
854 if( blk->brick_mdfcurve == NULL )
855 blk->brick_mdfcurve = (floatvec **)calloc(sizeof(floatvec *),blk->nvals);
856 blk->brick_mdfcurve[ibr] = fv ;
857 }
858 }
859
860 /*-- ID codes --*/
861
862 if( ATR_IS_STR(ATRNAME_IDSTRING) )
863 MCW_strncpy( dset->idcode.str , atr_str->ch , MCW_IDSIZE ) ;
864
865 if( ATR_IS_STR(ATRNAME_IDDATE) )
866 MCW_strncpy( dset->idcode.date , atr_str->ch , MCW_IDDATE ) ;
867
868 if( ATR_IS_STR(ATRNAME_IDANATPAR) )
869 MCW_strncpy( dset->anat_parent_idcode.str , atr_str->ch , MCW_IDSIZE ) ;
870
871 if( ATR_IS_STR(ATRNAME_IDWARPPAR) )
872 MCW_strncpy( dset->warp_parent_idcode.str , atr_str->ch , MCW_IDSIZE ) ;
873
874 /*-- parent names --*/
875
876 if( ATR_IS_STR(ATRNAME_ANATOMY_PARENT) &&
877 ISZERO_IDCODE(dset->anat_parent_idcode) )
878 MCW_strncpy( dset->anat_parent_name , atr_str->ch , THD_MAX_NAME ) ;
879
880 if( ATR_IS_STR(ATRNAME_WARP_PARENT) &&
881 ISZERO_IDCODE(dset->warp_parent_idcode) )
882 MCW_strncpy( dset->warp_parent_name , atr_str->ch , THD_MAX_NAME ) ;
883
884 if( ATR_IS_STR(ATRNAME_DATANAME) )
885 MCW_strncpy( dset->self_name , atr_str->ch , THD_MAX_NAME ) ;
886
887 /*-- markers --*/
888
889 if( ATR_IS_FLT(ATRNAME_MARKSXYZ) && ATR_IS_STR(ATRNAME_MARKSLAB) ){
890 int im , llen ;
891 THD_ivec3 iv ;
892 float xxdown,xxup , yydown,yyup , zzdown,zzup ;
893
894 STATUS("markers") ;
895
896 if( dset->markers == NULL ){
897 dset->markers = myRwcNew( THD_marker_set ) ; /* new set */
898 ADDTO_KILL(dset->kl , dset->markers) ;
899 }
900
901 COPY_INTO_STRUCT( *(dset->markers) , /* actual struct */
902 MARKS_FSTART , /* byte offset */
903 float , /* type being copied */
904 atr_flt->fl , /* start of source */
905 MARKS_FSIZE ) ; /* number of floats */
906
907 COPY_INTO_STRUCT( *(dset->markers) ,
908 MARKS_LSTART ,
909 char ,
910 atr_str->ch ,
911 MARKS_LSIZE ) ;
912
913 xxdown = daxes->xxmin ; xxup = daxes->xxmax ;
914 yydown = daxes->yymin ; yyup = daxes->yymax ;
915 zzdown = daxes->zzmin ; zzup = daxes->zzmax ;
916
917 dset->markers->numdef = dset->markers->numset = 0 ;
918
919 for( im=0 ; im < MARKS_MAXNUM ; im++ ){
920 llen = strlen( &(dset->markers->label[im][0]) ) ;
921 dset->markers->valid[im] =
922 (llen > 0) &&
923 ( dset->markers->xyz[im][0] >= xxdown ) &&
924 ( dset->markers->xyz[im][0] <= xxup ) &&
925 ( dset->markers->xyz[im][1] >= yydown ) &&
926 ( dset->markers->xyz[im][1] <= yyup ) &&
927 ( dset->markers->xyz[im][2] >= zzdown ) &&
928 ( dset->markers->xyz[im][2] <= zzup ) ;
929
930 if( dset->markers->valid[im] ) (dset->markers->numset)++ ;
931
932 if( llen > 0 ) (dset->markers->numdef)++ ;
933
934 dset->markers->ovcolor[im] = -1 ; /* default color */
935 }
936
937 if( ATR_IS_STR(ATRNAME_MARKSHELP) ){
938 COPY_INTO_STRUCT( *(dset->markers) ,
939 MARKS_HSTART ,
940 char ,
941 atr_str->ch ,
942 MARKS_HSIZE ) ;
943 } else {
944 for( im=0 ; im < MARKS_MAXNUM ; im++ ) /* no help */
945 dset->markers->help[im][0] = '\0' ;
946 }
947
948 if( ATR_IS_INT(ATRNAME_MARKSFLAG) ){
949 COPY_INTO_STRUCT( *(dset->markers) ,
950 MARKS_ASTART ,
951 int ,
952 atr_int->in ,
953 MARKS_ASIZE ) ;
954 dset->markers->type = dset->markers->aflags[0] ;
955 } else {
956 for( im=0 ; im < MARKS_MAXFLAG ; im++ )
957 dset->markers->aflags[im] = -1 ;
958 }
959 }
960
961 /*-- warp --*/
962
963 if( ATR_IS_INT(ATRNAME_WARP_TYPE) && ATR_IS_FLT(ATRNAME_WARP_DATA) ){
964 int wtype = atr_int->in[0] , rtype = atr_int->in[1] ;
965
966 STATUS("warp") ;
967
968 dset->warp = myRwcNew( THD_warp ) ;
969 ADDTO_KILL( dset->kl , dset->warp ) ;
970 switch( wtype ){
971 case WARP_AFFINE_TYPE:{
972 THD_affine_warp *ww = (THD_affine_warp *) dset->warp ;
973 ww->type = wtype ;
974 ww->resam_type = rtype ;
975 ww->warp.type = MAPPING_LINEAR_TYPE ;
976
977 COPY_INTO_STRUCT( ww->warp ,
978 MAPPING_LINEAR_FSTART ,
979 float ,
980 atr_flt->fl ,
981 MAPPING_LINEAR_FSIZE ) ;
982 }
983 break ; /* end affine warp */
984
985 case WARP_TALAIRACH_12_TYPE:{
986 THD_talairach_12_warp *ww =
987 (THD_talairach_12_warp *) dset->warp ;
988 int iw , ioff ;
989 ww->type = wtype ;
990 ww->resam_type = rtype ;
991 for( iw=0 ; iw < 12 ; iw++ ){
992 ww->warp[iw].type = MAPPING_LINEAR_TYPE ;
993
994 ioff = iw * MAPPING_LINEAR_FSIZE ;
995
996 COPY_INTO_STRUCT( ww->warp[iw] ,
997 MAPPING_LINEAR_FSTART ,
998 float ,
999 &(atr_flt->fl[ioff]) ,
1000 MAPPING_LINEAR_FSIZE ) ;
1001
1002 } /* end loop over 12 warps */
1003 }
1004 break ; /* end talairach_12 warp */
1005
1006 } /* end of switch on warp type */
1007 } /* end of if on warp existing */ else { /* But perhaps there is a little
1008 something from auto talairaching ZSS, June 06 */
1009 if (dset->view_type == VIEW_TALAIRACH_TYPE) { /* something to do */
1010 int dbg = 0;
1011 atr_flt = THD_find_float_atr( blk , ATRNAME_WARP_DATA_3DWD_AF ) ;
1012 if ( atr_flt == NULL ){
1013 /* A tlrc set with no transform. No problem */
1014 /* fprintf(stderr,"Dude, where's my transform?\n"); */
1015 } else {
1016 STATUS("AutoTlrc Warp") ;
1017 if (dbg)
1018 fprintf(stderr,
1019 "++ Will be using %s attribute for talairach warp in"
1020 " dset %s\n",
1021 ATRNAME_WARP_DATA_3DWD_AF, dset->self_name) ;
1022 if (!THD_WarpData_From_3dWarpDrive(dset, atr_flt)) {
1023 fprintf(stderr,"Error: Failed to create WarpData!\n");
1024 }
1025 }
1026 } else {
1027 /* fprintf(stderr,"Not in TLRC space, bother not.\n"); */
1028 }
1029 } /* the very end of if on warp existing */
1030
1031 /*-- brick stats --*/
1032
1033 if( ATR_IS_FLT(ATRNAME_BRICK_STATS) ){
1034 int qq ;
1035 STATUS("brick statistics") ;
1036 dset->stats = myRwcNew( THD_statistics ) ;
1037 dset->stats->type = STATISTICS_TYPE ;
1038 dset->stats->parent = (RwcPointer) dset ;
1039 dset->stats->nbstat = blk->nvals ;
1040 dset->stats->bstat = (THD_brick_stats *)
1041 RwcMalloc( sizeof(THD_brick_stats) * blk->nvals ) ;
1042 for( qq=0 ; qq < blk->nvals ; qq++ ){
1043 if( 2*qq+1 < atr_flt->nfl ){
1044 dset->stats->bstat[qq].min = atr_flt->fl[2*qq] ;
1045 dset->stats->bstat[qq].max = atr_flt->fl[2*qq+1] ;
1046 } else {
1047 INVALIDATE_BSTAT( dset->stats->bstat[qq] ) ;
1048 }
1049 }
1050 ADDTO_KILL( dset->kl , dset->stats->bstat ) ;
1051 ADDTO_KILL( dset->kl , dset->stats ) ;
1052 }
1053
1054 /*-- tagset --*/
1055
1056 if( ATR_IS_INT(ATRNAME_TAGSET_NUM) &&
1057 ATR_IS_FLT(ATRNAME_TAGSET_FLOATS) &&
1058 ATR_IS_STR(ATRNAME_TAGSET_LABELS) ){
1059
1060 int nin=atr_int->nin , nfl=atr_flt->nfl , nch=atr_str->nch ;
1061 int ii , ntag , nfper , jj , kk ;
1062
1063 STATUS("tagset") ;
1064 ntag = atr_int->in[0] ; /* number of tags */
1065 nfper = atr_int->in[1] ; /* number of floats per tag */
1066
1067 if( ntag > MAX_TAG_NUM ) ntag = MAX_TAG_NUM ;
1068
1069 dset->tagset = myRwcNew( THD_usertaglist ) ; /* create tagset */
1070 ADDTO_KILL( dset->kl , dset->tagset ) ;
1071
1072 dset->tagset->num = ntag ;
1073 strcpy( dset->tagset->label , "Bebe Rebozo" ) ; /* not used */
1074
1075 /* read out tag values; allow for chance there isn't enough data */
1076
1077 #undef TF
1078 #define TF(i,j) \
1079 ( ((j)<nfper && (i)*nfper+(j)<nfl) ? atr_flt->fl[(i)*nfper+(j)] : -666.0 )
1080
1081 for( ii=0 ; ii < ntag ; ii++ ){
1082 dset->tagset->tag[ii].x = TF(ii,0) ; /* coords */
1083 dset->tagset->tag[ii].y = TF(ii,1) ;
1084 dset->tagset->tag[ii].z = TF(ii,2) ;
1085 dset->tagset->tag[ii].val = TF(ii,3) ; /* value */
1086 dset->tagset->tag[ii].ti = TF(ii,4) ; /* time index: if < 0, not set */
1087 if( dset->tagset->tag[ii].ti >= 0 ){
1088 dset->tagset->tag[ii].set = 1 ;
1089 } else {
1090 dset->tagset->tag[ii].set = 0 ; dset->tagset->tag[ii].ti = 0 ;
1091 }
1092 }
1093 #undef TF
1094
1095 /* read out tag labels; allow for empty labels */
1096
1097 jj = 0 ;
1098 for( ii=0 ; ii < ntag ; ii++ ){
1099 if( jj < nch ){
1100 kk = strlen( atr_str->ch + jj ) ;
1101 if( kk > 0 ) TAG_SETLABEL( dset->tagset->tag[ii] , atr_str->ch + jj );
1102 else sprintf( dset->tagset->tag[ii].label , "Tag %d" , ii+1 );
1103 jj += kk+1 ;
1104 } else {
1105 sprintf( dset->tagset->tag[ii].label , "Tag %d" , ii+1 ) ;
1106 }
1107 }
1108 }
1109
1110 if( ( atr_flt = THD_find_float_atr(blk,"IJK_TO_DICOM_REAL") ) ){
1111 /* load oblique transformation matrix */
1112 if(atr_flt) {
1113 LOAD_MAT44(dset->daxes->ijk_to_dicom_real, \
1114 atr_flt->fl[0], atr_flt->fl[1], atr_flt->fl[2], atr_flt->fl[3], \
1115 atr_flt->fl[4], atr_flt->fl[5], atr_flt->fl[6], atr_flt->fl[7], \
1116 atr_flt->fl[8], atr_flt->fl[9], atr_flt->fl[10], atr_flt->fl[11]);
1117 }
1118 }
1119
1120 /* update attributes for time axes - copied from thd_dsetdblk.c */
1121 /*--- read time-dependent information, if any ---*/
1122 atr_int = THD_find_int_atr ( blk , ATRNAME_TAXIS_NUMS ) ;
1123 atr_flt = THD_find_float_atr( blk , ATRNAME_TAXIS_FLOATS ) ;
1124 if( atr_int != NULL && atr_flt != NULL ){
1125 int isfunc , nvals ;
1126
1127 dset->taxis = myRwcNew( THD_timeaxis ) ;
1128
1129 dset->taxis->type = TIMEAXIS_TYPE ;
1130 dset->taxis->ntt = atr_int->in[0] ;
1131 dset->taxis->nsl = atr_int->in[1] ;
1132 dset->taxis->ttorg = atr_flt->fl[0] ;
1133 dset->taxis->ttdel = atr_flt->fl[1] ;
1134 dset->taxis->ttdur = atr_flt->fl[2] ;
1135 dset->taxis->zorg_sl = atr_flt->fl[3] ;
1136 dset->taxis->dz_sl = atr_flt->fl[4] ;
1137
1138 dset->taxis->units_type = atr_int->in[2] ; /* 21 Oct 1996 */
1139 if( dset->taxis->units_type < 0 ) /* assign units */
1140 dset->taxis->units_type = UNITS_SEC_TYPE ; /* to the time axis */
1141
1142 if( dset->taxis->nsl > 0 ){
1143 atr_flt = THD_find_float_atr( blk , ATRNAME_TAXIS_OFFSETS ) ;
1144 if( atr_flt == NULL || atr_flt->nfl < dset->taxis->nsl ){
1145 dset->taxis->nsl = 0 ;
1146 dset->taxis->toff_sl = NULL ;
1147 dset->taxis->zorg_sl = 0.0 ;
1148 dset->taxis->dz_sl = 0.0 ;
1149 } else {
1150 int ii ;
1151 dset->taxis->toff_sl = (float *) RwcMalloc(sizeof(float)*dset->taxis->nsl) ;
1152 for( ii=0 ; ii < dset->taxis->nsl ; ii++ )
1153 dset->taxis->toff_sl[ii] = atr_flt->fl[ii] ;
1154 }
1155 } else {
1156 dset->taxis->nsl = 0 ;
1157 dset->taxis->toff_sl = NULL ;
1158 dset->taxis->zorg_sl = 0.0 ;
1159 dset->taxis->dz_sl = 0.0 ;
1160 }
1161
1162 }
1163
1164 /* get template space for dataset */
1165 if(ATR_IS_STR("TEMPLATE_SPACE"))
1166 MCW_strncpy( dset->atlas_space , atr_str->ch , THD_MAX_NAME ) ;
1167 /* otherwise guess from view and environment */
1168 else
1169 THD_get_space(dset);
1170
1171 /* flag for whether dataset should be displayed with integer colormap */
1172 if(ATR_IS_INT("INT_CMAP"))
1173 dset->int_cmap = atr_int->in[0];
1174 /* otherwise guess it's just regular old data */
1175 else
1176 dset->int_cmap = CONT_CMAP;
1177
1178 EXRETURN ;
1179 }
1180
1181 /*----------------------------------------------------------------
1182 Initialize the brick structure of a datablock. This will
1183 contain no data for now, but provides a place for it to go.
1184 This routine presumes that the datablock is already set up
1185 with the dimension information, etc.
1186
1187 Inputs ntype and btype determine the types of data in the
1188 sub-bricks.
1189 If btype == NULL, then the type code for all sub-bricks
1190 is given by ntype;
1191 If ntype < 0, then the type code for all sub-bricks
1192 is taken from the datablock pointed
1193 to by (THD_datablock *) btype;
1194 If btype != NULL then the type code for each sub-brick
1195 && ntype > 0 is taken from the array pointed to
1196 by (int *) btype.
1197 ------------------------------------------------------------------*/
1198
THD_init_datablock_brick(THD_datablock * dblk,int ntype,void * btype)1199 void THD_init_datablock_brick( THD_datablock *dblk, int ntype, void *btype )
1200 {
1201 int ibr , nx,ny,nz , typ , nvals ;
1202 MRI_IMAGE *qim ;
1203 THD_datablock *pblk = NULL ;
1204 int *itype = NULL ;
1205
1206 ENTRY("THD_init_datablock_brick") ;
1207
1208 if( ! ISVALID_DATABLOCK(dblk) ) EXRETURN ; /* bad inputs */
1209 if( ntype < 0 && btype == NULL ) EXRETURN ;
1210 if( ntype == 0 && btype != NULL ) EXRETURN ;
1211
1212 if( ntype < 0 ){ /* copy types from */
1213 pblk = (THD_datablock *) btype ; /* datablock pblk */
1214 if( ! ISVALID_DATABLOCK(pblk) ) EXRETURN ;
1215 } else {
1216 itype = (int *) btype ;
1217 }
1218
1219 nx = dblk->diskptr->dimsizes[0] ;
1220 ny = dblk->diskptr->dimsizes[1] ;
1221 nz = dblk->diskptr->dimsizes[2] ;
1222 nvals = dblk->nvals ; if( nvals < 1 ) EXRETURN ; /* something wrong */
1223
1224 /** make brick information arrays, if not pre-existing **/
1225
1226 if( dblk->brick_bytes == NULL ){
1227 STATUS("making dblk->brick_bytes") ;
1228 dblk->brick_bytes = (int64_t *) RwcMalloc( sizeof(int64_t) * nvals ) ;
1229 }
1230
1231 if( dblk->brick_fac == NULL ){
1232 STATUS("making dblk->brick_fac") ;
1233 dblk->brick_fac = (float *) RwcMalloc( sizeof(float) * nvals ) ;
1234 for( ibr=0 ; ibr < nvals ; ibr++ )
1235 dblk->brick_fac[ibr] = (ntype < 0) ? pblk->brick_fac[ibr] : 0.0 ;
1236 }
1237
1238 dblk->total_bytes = 0 ;
1239
1240 if( dblk->brick != NULL ){
1241 STATUS("destroying old dblk->brick") ;
1242 DESTROY_IMARR( dblk->brick ) ;
1243 }
1244
1245 STATUS("making new dblk->brick") ;
1246 INIT_IMARR( dblk->brick ) ; /* make the new brick */
1247
1248 /** set up each sub-brick **/
1249
1250 STATUS("starting sub-brick creations") ;
1251 for( ibr=0 ; ibr < nvals ; ibr++ ){
1252 if( ntype < 0 ){ typ = DBLK_BRICK_TYPE(pblk,ibr) ;
1253 } else if( itype == NULL ){
1254 typ = ntype ; /* all types are the same */
1255 } else {
1256 if( ibr < ntype ) typ = itype[ibr] ; /* types may vary */
1257 else typ = itype[ntype-1] ;
1258 }
1259
1260 qim = mri_new_vol_empty( nx,ny,nz , typ ) ; /* image with no data */
1261 ADDTO_IMARR( dblk->brick , qim ) ;
1262
1263 dblk->brick_bytes[ibr] = (int64_t)qim->pixel_size * (int64_t)qim->nvox ;
1264 dblk->total_bytes += dblk->brick_bytes[ibr] ;
1265 }
1266
1267 STATUS("exiting") ;
1268 EXRETURN ;
1269 }
1270
1271 /*----------------------------------------------------------------
1272 Determine if the brick factors are needed.
1273 ------------------------------------------------------------------*/
1274
THD_need_brick_factor(THD_3dim_dataset * dset)1275 int THD_need_brick_factor( THD_3dim_dataset * dset )
1276 {
1277 int ii , nval ;
1278
1279 /** ENTRY("THD_need_brick_factor") ; **/
1280
1281 if( ! ISVALID_DSET(dset) ) return( 0 ) ;
1282 if( ! ISVALID_DATABLOCK(dset->dblk) ) return( 0 ) ;
1283 if( dset->dblk->brick_fac == NULL ) return( 0 ) ;
1284
1285 nval = DSET_NVALS(dset) ;
1286 for( ii=0 ; ii < nval ; ii++ )
1287 if( DSET_BRICK_FACTOR(dset,ii) != 0.0 &&
1288 DSET_BRICK_FACTOR(dset,ii) != 1.0 ) return( 1 ) ;
1289
1290 return( 0 ) ;
1291 }
1292