1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8 #include "thd.h"
9
10 static int native_order = -1 ;
11 static int no_mmap = -1 ;
12 static int floatscan = -1 ; /* 30 Jul 1999 */
13
14 #define PRINT_SIZE 123456789
15 #define PRINT_STEP 10
16
17 static int verbose = 0 ;
18
THD_load_datablock_verbose(int v)19 void THD_load_datablock_verbose( int v ){ verbose = v; }
20
THD_load_no_mmap(void)21 void THD_load_no_mmap(void){ no_mmap = 2 ; }
22
23 int THD_alloc_datablock( THD_datablock *blk ) ; /* prototype */
24
25 /*-----------------------------------------------------------------*/
26 /*! Check if all sub-bricks have the same datum type. [14 Mar 2002]
27 -------------------------------------------------------------------*/
28
THD_datum_constant(THD_datablock * blk)29 int THD_datum_constant( THD_datablock *blk )
30 {
31 int ibr , dzero , nv=blk->nvals ;
32 if( nv == 1 ) return 1 ; /* of course */
33 dzero = DBLK_BRICK_TYPE(blk,0) ; /* #0 type */
34 for( ibr=1 ; ibr < nv ; ibr++ )
35 if( dzero != DBLK_BRICK_TYPE(blk,ibr) ) return 0 ;
36 return 1 ;
37 }
38
39 /*---------------------------------------------------------------*/
40 /*! Return a float representation of a given voxel. [15 Sep 2004]
41 *//*-------------------------------------------------------------*/
42
THD_get_voxel(THD_3dim_dataset * dset,int ijk,int ival)43 float THD_get_voxel( THD_3dim_dataset *dset , int ijk , int ival )
44 {
45 void *ar ;
46 float val , fac ;
47
48 if( !ISVALID_DSET(dset) ) return 0.0f ;
49 if( ival < 0 || ival >= DSET_NVALS(dset) ) return 0.0f ;
50 if( ijk < 0 || ijk >= DSET_NVOX(dset) ) return 0.0f ;
51
52 ar = DSET_ARRAY(dset,ival) ;
53 if( ar == NULL ){ DSET_load(dset) ; ar = DSET_ARRAY(dset,ival) ; }
54 if( ar == NULL ) return 0.0f ;
55
56 switch( DSET_BRICK_TYPE(dset,ival) ){
57
58 default: return 0.0f ;
59
60 case MRI_byte:
61 val = (float)(((byte *)ar)[ijk]) ; break ;
62 case MRI_short:
63 val = (float)(((short *)ar)[ijk]) ; break ;
64 case MRI_int:
65 val = (float)(((int *)ar)[ijk]) ; break ;
66 case MRI_float:
67 val = (float)(((float *)ar)[ijk]) ; break ;
68 case MRI_double:
69 val = (float)(((double *)ar)[ijk]) ; break ;
70
71 case MRI_complex:{
72 complex c = (((complex *)ar)[ijk]) ;
73 val = sqrtf(c.r*c.r+c.i*c.i) ;
74 break ;
75 }
76
77 case MRI_rgb:{
78 rgbyte c = (((rgbyte *)ar)[ijk]) ;
79 val = 0.299f*(float)(c.r) + 0.587f*(float)(c.g) + 0.114f*(float)(c.b) ;
80 break ;
81 }
82
83 case MRI_rgba:{
84 rgba c = (((rgba *)ar)[ijk]) ;
85 val = 0.299f*(float)(c.r) + 0.587f*(float)(c.g) + 0.114f*(float)(c.b) ;
86 val *= 0.00392157f*(float)(c.a) ;
87 break ;
88 }
89 }
90
91 fac = DSET_BRICK_FACTOR(dset,ival) ;
92 if( fac > 0.0f ) val *= fac ;
93 return val ;
94 }
95
96 /*----------------------------------------------------------------------------*/
97
THD_get_voxel_dicom(THD_3dim_dataset * dset,float x,float y,float z,int ival)98 float THD_get_voxel_dicom( THD_3dim_dataset *dset, float x,float y,float z, int ival )
99 {
100 THD_fvec3 vv , uu ; THD_ivec3 ii ; int ijk ;
101
102 if( !ISVALID_DSET(dset) ) return 0.0f ;
103 if( ival < 0 || ival >= DSET_NVALS(dset) ) return 0.0f ;
104
105 LOAD_FVEC3( vv , x,y,z ) ;
106 uu = THD_dicomm_to_3dmm( dset , vv ) ;
107 ii = THD_3dmm_to_3dind_no_wod( dset , uu ) ;
108 ijk = DSET_ixyz_to_index(dset,ii.ijk[0],ii.ijk[1],ii.ijk[2]) ;
109
110 return THD_get_voxel( dset , ijk , ival ) ;
111 }
112
113 /*---------------------------------------------------------------
114 18 Oct 2001:
115 Put freeup function here, and set it by a function, rather
116 than have it provided on every call to THD_load_datablock()
117 ----------------------------------------------------------------*/
118
119 static generic_func *freeup=NULL ; /* 18 Oct 2001 */
120
THD_set_freeup(generic_func * ff)121 void THD_set_freeup( generic_func *ff ){ freeup = ff; }
122
123 /*---------------------------------------------------------------*/
124
THD_load_datablock(THD_datablock * blk)125 RwcBoolean THD_load_datablock( THD_datablock *blk )
126 {
127 THD_diskptr *dkptr ;
128 int nx,ny,nz , nxy,nxyz , nv , ii , ibr ;
129 char *ptr ;
130 int verb=verbose ;
131 int64_t idone , print_size=PRINT_SIZE , id ;
132
133 ENTRY("THD_load_datablock") ; /* 29 Aug 2001 */
134
135 if( native_order < 0 ) native_order = mri_short_order() ; /* 1st time in */
136
137 floatscan = AFNI_yesenv("AFNI_FLOATSCAN") ; /* check float bricks? */
138
139 if( no_mmap < 2 ){
140 if( floatscan ) no_mmap = 1 ; /* perhaps disable */
141 else no_mmap = AFNI_yesenv("AFNI_NOMMAP") ; /* use of mmap() */
142 }
143
144 /*-- sanity checks --*/
145
146 if( ! ISVALID_DATABLOCK(blk) || blk->brick == NULL ){
147 STATUS("Illegal inputs"); RETURN( False );
148 }
149 ii = THD_count_databricks( blk ) ;
150 if( ii == blk->nvals ) RETURN( True ); /* already loaded! */
151
152 if( blk->malloc_type == DATABLOCK_MEM_UNDEFINED ){
153 STATUS("malloc_type forced to DATABLOCK_MEM_MALLOC") ;
154 blk->malloc_type = DATABLOCK_MEM_MALLOC ;
155 }
156
157 /* these next 2 conditions should never happen */
158
159 dkptr = blk->diskptr ;
160 if( ! ISVALID_DISKPTR(dkptr) ){
161 STATUS("invalid dkptr!!!"); RETURN( False );
162 }
163 if( dkptr->storage_mode == STORAGE_UNDEFINED ){
164 if( PRINT_TRACING ){
165 char str[512] ; THD_3dim_dataset *dset=(THD_3dim_dataset *)(blk->parent) ;
166 sprintf(str,"dataset %s == STORAGE_UNDEFINED",DSET_BRIKNAME(dset)) ;
167 STATUS(str) ;
168 }
169 RETURN(False) ;
170 }
171
172 if( dkptr->rank != 3 ){
173 fprintf(stderr,"\n*** Cannot read non 3D datablocks ***\n"); RETURN(False);
174 }
175
176 if( !no_mmap && dkptr->storage_mode == STORAGE_BY_VOLUMES ) no_mmap = 1 ; /* 20 Jun 2002 */
177
178 /*-- 29 Oct 2001: MINC input (etc.) --*/
179
180 if( dkptr->storage_mode == STORAGE_BY_MINC ){
181 ERROR_message("MINC-1 dataset input support is disabled") ;
182 }
183
184 { THD_3dim_dataset *ds = (THD_3dim_dataset *)blk->parent ; /* 04 Aug 2004 */
185 if( ISVALID_DSET(ds) && DSET_IS_TCAT(ds) ){
186 THD_load_tcat( blk ) ;
187 ii = THD_count_databricks( blk ) ;
188 if( ii == blk->nvals ){
189 THD_update_statistics( (THD_3dim_dataset *)blk->parent ) ;
190 RETURN( True ) ;
191 }
192 STATUS("can't read tcat-ed file?!") ;
193 RETURN( False ) ;
194 }
195 }
196
197 if( dkptr->storage_mode == STORAGE_BY_ANALYZE ){
198 THD_load_analyze( blk ) ;
199 ii = THD_count_databricks( blk ) ;
200 if( ii == blk->nvals ){
201 THD_update_statistics( (THD_3dim_dataset *)blk->parent ) ;
202 ii = dblk_floatscan(blk) ; /* 22 Feb 2007 */
203 if(ii>0) WARNING_message("fixed %d bad floats in %s",ii,dkptr->brick_name);
204 RETURN( True ) ;
205 }
206 STATUS("can't read ANALYZE file?!") ;
207 RETURN( False ) ;
208 }
209
210 if( dkptr->storage_mode == STORAGE_BY_CTFMRI ){ /* 04 Dec 2002 */
211 THD_load_ctfmri( blk ) ;
212 ii = THD_count_databricks( blk ) ;
213 if( ii == blk->nvals ){
214 THD_update_statistics( (THD_3dim_dataset *)blk->parent ) ;
215 ii = dblk_floatscan(blk) ; /* 22 Feb 2007 */
216 if(ii>0) WARNING_message("fixed %d bad floats in %s",ii,dkptr->brick_name);
217 RETURN( True ) ;
218 }
219 STATUS("can't read CTF MRI file?!") ;
220 RETURN( False ) ;
221 }
222
223 if( dkptr->storage_mode == STORAGE_BY_CTFSAM ){ /* 04 Dec 2002 */
224 THD_load_ctfsam( blk ) ;
225 ii = THD_count_databricks( blk ) ;
226 if( ii == blk->nvals ){
227 THD_update_statistics( (THD_3dim_dataset *)blk->parent ) ;
228 ii = dblk_floatscan(blk) ; /* 22 Feb 2007 */
229 if(ii>0) WARNING_message("fixed %d bad floats in %s",ii,dkptr->brick_name);
230 RETURN( True ) ;
231 }
232 STATUS("can't read CTF SAM file?!") ;
233 RETURN( False ) ;
234 }
235
236 if( dkptr->storage_mode == STORAGE_BY_1D ){ /* 04 Mar 2003 */
237 THD_load_1D( blk ) ;
238 ii = THD_count_databricks( blk ) ;
239 if( ii == blk->nvals ){
240 THD_update_statistics( (THD_3dim_dataset *)blk->parent ) ;
241 RETURN( True ) ;
242 }
243 STATUS("can't read 1D dataset file?!") ;
244 RETURN( False ) ;
245 }
246
247 if( dkptr->storage_mode == STORAGE_BY_3D ){ /* 21 Mar 2003 */
248 THD_load_3D( blk ) ;
249 ii = THD_count_databricks( blk ) ;
250 if( ii == blk->nvals ){
251 THD_update_statistics( (THD_3dim_dataset *)blk->parent ) ;
252 RETURN( True ) ;
253 }
254 STATUS("can't read 3D dataset file?!") ;
255 RETURN( False ) ;
256 }
257
258 if( dkptr->storage_mode == STORAGE_BY_NIFTI ){ /* 28 Aug 2003 */
259 THD_load_nifti( blk ) ;
260 ii = THD_count_databricks( blk ) ;
261 if( ii == blk->nvals ){
262 THD_update_statistics( (THD_3dim_dataset *)blk->parent ) ;
263 RETURN( True ) ;
264 }
265 STATUS("can't read NIFTI dataset file?!") ;
266 RETURN( False ) ;
267 }
268
269 if( dkptr->storage_mode == STORAGE_BY_MPEG ){ /* 03 Dec 2003 */
270 THD_load_mpeg( blk ) ;
271 ii = THD_count_databricks( blk ) ;
272 if( ii == blk->nvals ){
273 THD_update_statistics( (THD_3dim_dataset *)blk->parent ) ;
274 RETURN( True ) ;
275 }
276 STATUS("can't read MPEG dataset file?!") ;
277 RETURN( False ) ;
278 }
279
280 if( dkptr->storage_mode == STORAGE_BY_NIML ){ /* 26 May 2006 [rickr] */
281 if( THD_load_niml( blk ) ){
282 STATUS("failed to load NIML file"); RETURN( False );
283 }
284 if( THD_count_databricks( blk ) == blk->nvals ){
285 THD_update_statistics( (THD_3dim_dataset *)blk->parent ) ;
286 ii = dblk_floatscan(blk) ; /* 22 Feb 2007 */
287 if(ii>0) WARNING_message("fixed %d bad floats in %s",ii,dkptr->brick_name);
288 RETURN( True ) ;
289 }
290 STATUS("can't read NIML dataset file?!") ;
291 RETURN( False ) ;
292 }
293
294 if( dkptr->storage_mode == STORAGE_BY_NI_SURF_DSET ) {
295 THD_load_niml( blk ) ;
296 ii = THD_count_databricks( blk ) ;
297 if( ii == blk->nvals ){
298 THD_update_statistics( (THD_3dim_dataset *)blk->parent ) ;
299 RETURN( True ) ;
300 }
301 STATUS("can't read NI_SURF_DSET dataset file?!") ;
302 RETURN( False ) ;
303 }
304
305 if( dkptr->storage_mode == STORAGE_BY_NI_TRACT ) {
306 ERROR_message("Cannot handle tracts here");
307 RETURN( False ) ;
308 }
309
310 if( dkptr->storage_mode == STORAGE_BY_GIFTI ) { /* 13 Feb 2008 */
311 THD_load_gifti( blk ) ;
312 ii = THD_count_databricks( blk ) ;
313 if( ii == blk->nvals ){
314 THD_update_statistics( (THD_3dim_dataset *)blk->parent ) ;
315 RETURN( True ) ;
316 }
317 STATUS("can't read GIFTI dataset file?!") ;
318 RETURN( False ) ;
319 }
320
321 if( dkptr->storage_mode == STORAGE_BY_IMAGE_FILE ){ /* 06 Jul 2016 */
322 WARNING_message("attempt to reload an image file 'dataset' -- not supported") ;
323 STATUS("can't reload image file dataset") ;
324 RETURN( False ) ;
325 }
326
327 /*** END OF SPECIAL CASES ABOVE; NOW FOR THE AFNI FORMAT! ***/
328
329 /*-- allocate data space --*/
330
331 nx = dkptr->dimsizes[0] ;
332 ny = dkptr->dimsizes[1] ; nxy = nx * ny ;
333 nz = dkptr->dimsizes[2] ; nxyz = nxy * nz ;
334 nv = dkptr->nvals ;
335
336 if( DBLK_IS_MASTERED(blk) && blk->malloc_type == DATABLOCK_MEM_MMAP ) /* 11 Jan 1999 */
337 blk->malloc_type = DATABLOCK_MEM_MALLOC ;
338
339 /* the following code is due to Mike Beauchamp's idiocy */
340
341 if( !THD_datum_constant(blk) && blk->malloc_type != DATABLOCK_MEM_MALLOC ){
342 fprintf(stderr,"++ WARNING: dataset %s: non-uniform sub-brick types\n",
343 blk->diskptr->filecode) ;
344 blk->malloc_type = DATABLOCK_MEM_MALLOC ;
345 }
346
347 /* 25 April 1998: byte order issues */
348
349 if( dkptr->byte_order <= 0 ) dkptr->byte_order = native_order ;
350
351 /* 05 Jul 2001: if all sub-bricks are bytes,
352 mark dataset as being in native order
353 20 Jun 2002: also don't have to swap RGB datasets */
354
355 if( dkptr->byte_order != native_order ){
356 for( ii=0 ; ii < nv ; ii++ )
357 if( DBLK_BRICK_TYPE(blk,ii) != MRI_byte &&
358 DBLK_BRICK_TYPE(blk,ii) != MRI_rgb ) break ;
359 if( ii == nv ) dkptr->byte_order = native_order ;
360 }
361
362 /* under some circumstances, we must force use of malloc() */
363
364 if( dkptr->byte_order != native_order || no_mmap ){
365 if( blk->malloc_type == DATABLOCK_MEM_MMAP )
366 blk->malloc_type = DATABLOCK_MEM_MALLOC ;
367 }
368
369 DBLK_mmapfix(blk) ; /* 18 Mar 2005 */
370
371 /** set up space for bricks via malloc, if so ordered **/
372
373 if( blk->malloc_type == DATABLOCK_MEM_MALLOC ||
374 blk->malloc_type == DATABLOCK_MEM_SHARED ){
375
376 ii = THD_alloc_datablock( blk ) ; /* 02 May 2003 */
377 if( ii == 0 ) RETURN(False) ;
378
379 /** mmap whole file (makes space and reads it all at once) **/
380 /** [if this route is followed, then we will finish here] **/
381
382 } else if( blk->malloc_type == DATABLOCK_MEM_MMAP ){
383
384 int fd ; long long fsize ;
385 fd = open( dkptr->brick_name , O_RDONLY ) ; /* N.B.: readonly mode */
386 if( fd < 0 ){
387 fprintf( stderr , "\n*** cannot open brick file %s for mmap\n"
388 " - do you have permission? does file exist?\n" ,
389 dkptr->brick_name ) ;
390 perror("*** Unix error message") ; TRACEBACK ;
391 STATUS("open failed") ;
392 RETURN( False );
393 }
394
395 /* 04 May 2001: check file size (the Katie Lee bug) */
396
397 fsize = THD_filesize( dkptr->brick_name ) ;
398 if( fsize < blk->total_bytes )
399 fprintf(stderr ,
400 "\n*** WARNING: file %s size is %lld, but should be at least %lld!\n" ,
401 dkptr->brick_name , fsize , (long long)blk->total_bytes ) ;
402
403 /* clear the sub-brick pointers */
404
405 for( ibr=0 ; ibr < nv ; ibr++ )
406 mri_clear_data_pointer( DBLK_BRICK(blk,ibr) ) ;
407
408 /* map the file into memory */
409
410 ptr = (char *) mmap( 0 , (size_t)blk->total_bytes ,
411 PROT_READ , THD_MMAP_FLAG , fd , 0 ) ;
412
413 if( verbose && blk->total_bytes >= 1000000000ll )
414 fprintf(stderr , "%s: mmap(%s bytes) -- %s\n" ,
415 dkptr->brick_name ,
416 approximate_number_string((double)blk->total_bytes) ,
417 (ptr == (char*)(-1)) ? "FAILED" : "SUCCEEDED" ) ;
418
419 /* if that fails, maybe try again (via freeup) */
420
421 if( ptr == (char *)(-1) ){
422 fprintf(stderr ,
423 "\n*** cannot mmap brick file %s - maybe hit a system limit?\n" ,
424 dkptr->brick_name ) ;
425 perror("*** Unix error message") ; TRACEBACK ;
426 if( freeup != NULL ){
427 freeup() ; /* AFNI_purge_unused_dsets */
428 fprintf(stderr,"*** purge and retry mmap(%s bytes)",
429 approximate_number_string((double)blk->total_bytes) ) ;
430 ptr = (char *) mmap( 0 , (size_t)blk->total_bytes ,
431 PROT_READ , THD_MMAP_FLAG , fd , 0 ) ;
432 if( ptr == (char *)(-1) ){
433 fprintf(stderr,"*** FAILED: cannot fix problem :-(\n") ;
434 close(fd) ;
435 STATUS("freeup failed") ; RETURN( False );
436 } else {
437 fprintf(stderr,"*** SUCCEEDED: was able to fix problem :-)\n") ;
438 }
439 } else { /* no freeup function to try */
440 close(fd) ;
441 STATUS("mmap failed") ; RETURN( False );
442 }
443 }
444
445 close(fd) ; /* can close file after mmap-ing it */
446
447 /* (re)create pointers to all sub-bricks */
448
449 for( ibr=0 ; ibr < nv ; ibr++ ){
450 mri_fix_data_pointer( ptr , DBLK_BRICK(blk,ibr) ) ;
451 ptr += DBLK_BRICK_BYTES(blk,ibr) ;
452 }
453
454 STATUS("mmap succeeded") ;
455 RETURN( True ); /* finito */
456
457 } else {
458 STATUS("unknown malloc_type code?!") ;
459 RETURN( False ) ; /* should never happen */
460 }
461
462 /*** Below here, space for brick images was malloc()-ed,
463 and now we have to read data into them, somehow ***/
464
465 ptr = getenv("AFNI_LOAD_PRINTSIZE") ; /* 23 Aug 2002 */
466 if( verb && ptr != NULL ){
467 char *ept ;
468 idone = strtol( ptr , &ept , 10 ) ;
469 if( idone > 0 ){
470 if( *ept == 'K' || *ept == 'k' ) idone *= 1024 ;
471 else if( *ept == 'M' || *ept == 'm' ) idone *= 1024*1024 ;
472 print_size = idone ;
473 } else {
474 print_size = PRINT_SIZE ;
475 }
476 }
477
478 if( verb ) verb = (blk->total_bytes > print_size) ;
479 if( verb )
480 fprintf(stderr,"reading %s(%s bytes)",
481 dkptr->filecode ,
482 approximate_number_string((double)blk->total_bytes) ) ;
483
484 switch( dkptr->storage_mode ){
485
486 /*-- should never ever happen --*/
487
488 default:
489 fprintf(stderr,"\n*** illegal storage mode in read ***\n") ;
490 RETURN( False );
491 break ;
492
493 /*-- read everything from .BRIK file --*/
494
495 case STORAGE_BY_BRICK:{
496 FILE *far ;
497
498 STATUS("reading from BRIK file") ;
499
500 far = COMPRESS_fopen_read( dkptr->brick_name ) ;
501
502 if( far == NULL ){
503 THD_purge_datablock( blk , blk->malloc_type ) ;
504 fprintf(stderr,
505 "\n*** failure while opening brick file %s "
506 "- do you have permission?\n" ,
507 dkptr->brick_name ) ;
508 perror("*** Unix error message") ; TRACEBACK ;
509 RETURN( False );
510 }
511
512 /* read each sub-brick all at once */
513
514 idone = 0 ;
515 if( ! DBLK_IS_MASTERED(blk) ){ /* read each brick */
516
517 for( ibr=0 ; ibr < nv ; ibr++ ){
518 idone += fread( DBLK_ARRAY(blk,ibr), 1,
519 DBLK_BRICK_BYTES(blk,ibr), far ) ;
520
521 if( verb && ibr%PRINT_STEP == 0 ) fprintf(stderr,".") ;
522 }
523
524 } else { /* 11 Jan 1999: read brick from master, put into place(s) */
525
526 int nfilled = 0 , jbr ; int64_t nbuf=0 , nbr ;
527 char *buf=NULL ; /* temp buffer for master sub-brick */
528
529 /* loop over master sub-bricks until dataset is filled */
530 /* [because dataset might be compressed, must read them in
531 sequence, even if some intermediate ones aren't used. ] */
532
533 for( ibr=0 ; nfilled < nv && ibr < blk->master_nvals ; ibr++ ){
534
535 if( nbuf < blk->master_bytes[ibr] ){ /* make more space for it */
536 if( buf != NULL ) free(buf) ;
537 nbuf = blk->master_bytes[ibr] ;
538 buf = AFMALL(char, sizeof(char) * nbuf ) ;
539 if( buf == NULL ) break ;
540 }
541
542 /* read the master sub-brick #ibr */
543
544 nbr = fread( buf , 1 , blk->master_bytes[ibr] , far ) ;
545 if( verb && ibr%PRINT_STEP == 0 ) fprintf(stderr,".") ;
546 if( nbr < blk->master_bytes[ibr] ) break ;
547
548 /* find all the dataset sub-bricks that are copies of this */
549
550 for( jbr=0 ; jbr < nv ; jbr++ ){
551 if( blk->master_ival[jbr] == ibr ){ /* copy it in */
552 memcpy( DBLK_ARRAY(blk,jbr) , buf , blk->master_bytes[ibr] ) ;
553 nfilled++ ; /* number of bricks filled */
554 idone += nbr ; /* number of bytes read into dataset */
555 }
556 }
557 } /* end of loop over master sub-bricks */
558
559 if( buf != NULL ) free(buf) ; /* free temp sub-brick */
560 }
561
562 /* close input file */
563
564 COMPRESS_fclose(far) ;
565
566 /* check if total amount of data read is correct */
567
568 if( idone != blk->total_bytes ){
569 THD_purge_datablock( blk , blk->malloc_type ) ;
570 fprintf(stderr ,
571 "\n*** failure while reading from brick file %s\n"
572 "*** desired %lld bytes but only got %lld\n" ,
573 dkptr->brick_name ,
574 (long long)blk->total_bytes , (long long)idone ) ;
575 perror("*** Unix error message") ; TRACEBACK ;
576 RETURN( False );
577 }
578 }
579 break ; /* end of STORAGE_BY_BRICK */
580
581 /*** Read from a sequence of volume files (1 per sub-brick) [20 Jun 2002] ***/
582
583 case STORAGE_BY_VOLUMES:{
584 ATR_string *atr ;
585 char **fnam , *ptr , *flist ;
586 FILE *far ;
587
588 STATUS("reading from volume files") ;
589
590 /* get list of filenames */
591
592 atr = THD_find_string_atr(blk,"VOLUME_FILENAMES") ;
593
594 /* the following should never happen */
595
596 if( atr == NULL || atr->nch < 2*nv ){
597 THD_purge_datablock( blk , blk->malloc_type ) ;
598 fprintf(stderr,
599 "Dataset %s does not have legal VOLUME_FILENAMES attribute!\n",
600 blk->diskptr->filecode) ;
601 RETURN( False );
602 }
603
604 /* copy filename list into NUL terminated local string */
605
606 flist = AFMALL(char, atr->nch+1) ;
607 memcpy( flist , atr->ch , atr->nch ) ;
608 flist[atr->nch] = '\0' ;
609
610 #if 0
611 fprintf(stderr,"VOL: flist:\n%s\n",flist) ;
612 #endif
613
614 /* break filename list into component filenames (1 per sub-brick) */
615
616 fnam = (char **) malloc(sizeof(char *)*nv) ;
617 for( ibr=0 ; ibr < nv ; ibr++ ){
618
619 /* find sub-string with next filename */
620
621 ptr = strtok( (ibr==0)?flist:NULL , " \t\n\r\f\v" ) ;
622
623 /* didn't get one => bad news */
624
625 if( ptr == NULL ){
626 fprintf(stderr,
627 "Dataset %s has illegal VOLUME_FILENAMES attribute!\n",
628 blk->diskptr->filecode) ;
629 for( ii=0 ; ii < nv ; ii++ ){
630 free( DBLK_ARRAY(blk,ii) ) ;
631 mri_clear_data_pointer( DBLK_BRICK(blk,ii) ) ;
632 }
633 for( ii=0 ; ii < ibr ; ii++ ) free(fnam[ii]) ;
634 free(flist) ; free(fnam) ;
635 RETURN( False );
636 }
637
638 /* make fnam[ibr] = name ibr-th volume file */
639
640 #if 0
641 fprintf(stderr,"VOL[%d]: ptr=%s\n",ibr,ptr) ;
642 #endif
643
644 if( *ptr == '/' ){ /* if filename is absolute, duplicate it */
645 fnam[ibr] = strdup(ptr) ;
646 } else { /* otherwise, put directory name in front */
647 ii = strlen(blk->diskptr->directory_name) + strlen(ptr) ;
648 fnam[ibr] = AFMALL(char,ii+1) ; fnam[ibr][0] = '\0' ;
649 strcat(fnam[ibr],blk->diskptr->directory_name) ;
650 strcat(fnam[ibr],ptr) ;
651 }
652 #if 0
653 fprintf(stderr,"VOL[%d]: fnam=%s\n",ibr,fnam[ibr]) ;
654 #endif
655 }
656 free(flist) ; /* done with unparsed filename list copy */
657
658 /*** loop over sub-bricks and read them in ***/
659
660 if( ! DBLK_IS_MASTERED(blk) ){ /* read each brick */
661
662 for( ibr=0 ; ibr < nv ; ibr++ ){
663
664 #if 0
665 fprintf(stderr,"VOL[%d]: opening %s\n",ibr,fnam[ibr]) ;
666 #endif
667
668 far = COMPRESS_fopen_read( fnam[ibr] ) ; /* open file */
669
670 if( far == NULL ){ /* can't open it? */
671 fprintf(stderr,
672 "\n*** Failure while opening volume file %s "
673 "- do you have permission?\n" ,
674 fnam[ibr] ) ;
675 perror("*** Unix error message") ; TRACEBACK ;
676 THD_purge_datablock( blk , blk->malloc_type ) ;
677 for( ii=0 ; ii < nv ; ii++ ) free(fnam[ii]) ;
678 free(fnam); RETURN( False );
679 }
680
681 /* read all of file at once */
682
683 id = fread( DBLK_ARRAY(blk,ibr), 1,
684 DBLK_BRICK_BYTES(blk,ibr), far ) ;
685 if( verb && ibr%PRINT_STEP==0 ) fprintf(stderr,".") ;
686
687 #if 0
688 fprintf(stderr,"VOL[%d]: id=%d\n",ibr,id) ;
689 #endif
690
691 COMPRESS_fclose(far) ;
692
693 if( id < DBLK_BRICK_BYTES(blk,ibr) ){
694 fprintf(stderr,
695 "\n*** Volume file %s only gave %lld out of %lld bytes needed!\n",
696 fnam[ibr] , id , DBLK_BRICK_BYTES(blk,ibr) ) ;
697 THD_purge_datablock( blk , blk->malloc_type ) ;
698 for( ii=0 ; ii < nv ; ii++ ) free(fnam[ii]) ;
699 free(fnam); RETURN( False );
700 }
701
702 } /* end of loop over sub-bricks */
703
704 } else { /*** Mastered dataset ==> must read only selected sub-bricks ***/
705
706 int jbr ;
707
708 /** loop over output sub-bricks **/
709
710 for( jbr=0 ; jbr < nv ; jbr++ ){
711 ibr = blk->master_ival[jbr] ; /* master sub-brick index */
712 far = COMPRESS_fopen_read( fnam[ibr] ) ; /* open file */
713
714 if( far == NULL ){ /* can't open it? */
715 fprintf(stderr,
716 "\n*** Failure while opening volume file %s "
717 "- do you have permission?\n" ,
718 fnam[ibr] ) ;
719 perror("*** Unix error message") ; TRACEBACK ;
720 THD_purge_datablock( blk , blk->malloc_type ) ;
721 for( ii=0 ; ii < nv ; ii++ ) free(fnam[ii]) ;
722 free(fnam); RETURN( False );
723 }
724
725 /* read all of file at once */
726
727 id = fread( DBLK_ARRAY(blk,jbr), 1,
728 DBLK_BRICK_BYTES(blk,jbr), far ) ;
729 if( verb && jbr%PRINT_STEP == 0 ) fprintf(stderr,".") ;
730
731 COMPRESS_fclose(far) ;
732
733 if( id < DBLK_BRICK_BYTES(blk,jbr) ){
734 fprintf(stderr,
735 "\n*** Volume file %s only gave %lld out of %lld bytes needed!\n",
736 fnam[ibr] , id , DBLK_BRICK_BYTES(blk,jbr) ) ;
737 THD_purge_datablock( blk , blk->malloc_type ) ;
738 for( ii=0 ; ii < nv ; ii++ ) free(fnam[ii]) ;
739 free(fnam); RETURN( False );
740 }
741
742 } /* end of loop over sub-bricks */
743
744 } /* end of mastered vs. unmastered dataset */
745
746 /*** at this point, all the data has been read in ***/
747
748 /* volume filenames no longer needed */
749
750 for( ii=0 ; ii < nv ; ii++ ) free(fnam[ii]) ;
751 free(fnam) ;
752
753 }
754 break ; /* end of STORAGE_BY_VOLUMES */
755
756 } /* end of switch on storage modes */
757
758 /************* At this point, data is all in bricks;
759 now do any post-processing before exiting *************/
760
761 STATUS("data has been read in") ;
762
763 /* 25 April 1998: check and fix byte ordering */
764
765 if( dkptr->byte_order != native_order ){
766 STATUS("byte swapping") ;
767 if( verb ) fprintf(stderr,".byte swap") ;
768
769 for( ibr=0 ; ibr < nv ; ibr++ ){
770 switch( DBLK_BRICK_TYPE(blk,ibr) ){
771 default: break ;
772 case MRI_short:
773 mri_swap2( DBLK_BRICK_NVOX(blk,ibr) , DBLK_ARRAY(blk,ibr) ) ;
774 break ;
775
776 case MRI_complex: /* 14 Sep 1999: swap complex also! */
777 mri_swap4( 2*DBLK_BRICK_NVOX(blk,ibr), DBLK_ARRAY(blk,ibr)) ;
778 break ;
779
780 case MRI_float:
781 case MRI_int:
782 mri_swap4( DBLK_BRICK_NVOX(blk,ibr) , DBLK_ARRAY(blk,ibr) ) ;
783 break ;
784 }
785 }
786 }
787
788 /* 30 July 1999: check float sub-brick for errors? */
789 /* 14 Sep 1999: also check complex sub-bricks! */
790
791 if( floatscan &&
792 (( DBLK_BRICK_TYPE(blk,0) == MRI_float ) ||
793 ( DBLK_BRICK_TYPE(blk,0) == MRI_complex ))){
794 int nerr=0 ;
795 STATUS("float scanning") ;
796 if( verb ) fprintf(stderr,".float scan") ;
797
798 for( ibr=0 ; ibr < nv ; ibr++ ){
799 if( DBLK_BRICK_TYPE(blk,ibr) == MRI_float ){
800 nerr += thd_floatscan( DBLK_BRICK_NVOX(blk,ibr) ,
801 DBLK_ARRAY(blk,ibr) ) ;
802
803 } else if( DBLK_BRICK_TYPE(blk,ibr) == MRI_complex ) { /* 14 Sep 1999 */
804 nerr += thd_complexscan( DBLK_BRICK_NVOX(blk,ibr) ,
805 DBLK_ARRAY(blk,ibr) ) ;
806 }
807 }
808 if( nerr > 0 )
809 WARNING_message("file %s: fixed %d float errors" ,
810 dkptr->brick_name , nerr ) ;
811 }
812
813 /* 21 Feb 2001: if sub-ranging also used, clip values in brick */
814
815 #if 0
816 fprintf(stderr,"master_bot=%g master_top=%g\n",blk->master_bot,blk->master_top) ;
817 #endif
818 /* rcr - replace this with DBLK_IS_RANGE_MASTERED to also check for
819 * csv list
820 * - then call a new parent function to THD_apply_master_subrange
821 *
822 * *** same as in thd_load_nifti and thd_niml.c
823 */
824
825 if( DBLK_IS_MASTER_SUBRANGED(blk) )
826 THD_apply_master_subrange(blk) ;
827
828 if( verb ) fprintf(stderr,".done\n") ;
829
830 RETURN( True ) ; /* things are now cool */
831 }
832
833 /*----------------------------------------------------------------------------*/
834 /*! Allocate memory for the volumes in a datablock -- 02 May 2003.
835 Return value is 1 if OK, 0 if not.
836 *//*--------------------------------------------------------------------------*/
837
THD_alloc_datablock(THD_datablock * blk)838 int THD_alloc_datablock( THD_datablock *blk )
839 {
840 int ii,nbad,ibr, nv ;
841 char *ptr ;
842
843 ENTRY("THD_alloc_datablock") ;
844
845 /*-- sanity checks --*/
846
847 if( ! ISVALID_DATABLOCK(blk) || blk->brick == NULL ){
848 STATUS("Illegal inputs"); RETURN(0);
849 }
850
851 nv = blk->nvals ;
852 ii = THD_count_databricks( blk ) ;
853 if( ii == nv ) RETURN(1); /* already loaded! */
854
855 switch( blk->malloc_type ){
856
857 default: RETURN(0) ; /* stupid datablock */
858
859 /*-------------------------------------------*/
860 /* malloc separate arrays for each sub-brick */
861 /*-------------------------------------------*/
862
863 case DATABLOCK_MEM_MALLOC:{
864
865 /** malloc space for each brick separately **/
866
867 STATUS("trying to malloc sub-bricks") ;
868
869 if( verbose && blk->total_bytes >= 1000000000ll )
870 fprintf(stderr,"malloc(%s bytes",
871 approximate_number_string((double)blk->total_bytes) ) ;
872
873 for( nbad=ibr=0 ; ibr < nv ; ibr++ ){
874 if( DBLK_ARRAY(blk,ibr) == NULL ){
875 ptr = AFMALL(char, DBLK_BRICK_BYTES(blk,ibr) ) ;
876 mri_fix_data_pointer( ptr , DBLK_BRICK(blk,ibr) ) ;
877 if( ptr == NULL ) nbad++ ;
878 }
879 }
880 if( verbose && blk->total_bytes >= 1000000000ll ) fprintf(stderr,")") ;
881 if( nbad == 0 ) RETURN(1) ; /* things are cool */
882
883 /* at least one malloc() failed, so possibly try to free some space */
884
885 fprintf(stderr,
886 "\n** failed to malloc %d dataset bricks out of %d - is memory exhausted?\n",
887 nbad,nv ) ;
888 #ifdef USING_MCW_MALLOC
889 { char *str = MCW_MALLOC_status ;
890 if( str != NULL ) fprintf(stderr,"*** MCW_malloc summary: %s\n",str); }
891 #endif
892
893 if( freeup == NULL ){ /* don't have a freeup() function? */
894 for( ibr=0 ; ibr < nv ; ibr++ ){ /* then toss all data bricks and scram */
895 if( DBLK_ARRAY(blk,ibr) != NULL ){
896 free(DBLK_ARRAY(blk,ibr)) ;
897 mri_fix_data_pointer( NULL , DBLK_BRICK(blk,ibr) ) ;
898 }
899 }
900 STATUS("malloc failed, no freeup"); RETURN(0); /* go away */
901 }
902
903 /* space freeing is done by caller-supplied function freeup() */
904
905 fprintf(stderr,"** trying to free some memory\n") ; /* 18 Oct 2001 */
906 freeup() ; /* cf. AFNI_purge_unused_dsets() */
907
908 /* now try to malloc() those that failed before */
909
910 for( ibr=0 ; ibr < nv ; ibr++ ){
911 if( DBLK_ARRAY(blk,ibr) == NULL ){
912 ptr = AFMALL(char, DBLK_BRICK_BYTES(blk,ibr) ) ;
913 mri_fix_data_pointer( ptr , DBLK_BRICK(blk,ibr) ) ;
914 }
915 }
916
917 /* if it still failed, then free everything and go away */
918
919 if( THD_count_databricks(blk) < nv ){
920 fprintf(stderr,"** cannot free up enough memory\n") ;
921 #ifdef USING_MCW_MALLOC
922 { char *str = MCW_MALLOC_status ;
923 if( str != NULL ) fprintf(stderr,"*** MCW_malloc summary: %s\n",str); }
924 #endif
925 for( ibr=0 ; ibr < nv ; ibr++ ){ /* 18 Oct 2001 */
926 if( DBLK_ARRAY(blk,ibr) != NULL ){
927 free(DBLK_ARRAY(blk,ibr)) ;
928 mri_fix_data_pointer( NULL , DBLK_BRICK(blk,ibr) ) ;
929 }
930 }
931 STATUS("malloc failed; freeup failed"); RETURN(0); /* go away */
932 }
933
934 /* malloc() didn't fail this time! */
935
936 STATUS("malloc failed; freeup worked") ;
937 fprintf(stderr,"*** was able to free up enough memory\n") ;
938 #ifdef USING_MCW_MALLOC
939 { char *str = MCW_MALLOC_status ;
940 if( str != NULL ) fprintf(stderr,"*** MCW_malloc summary: %s\n",str); }
941 #endif
942 }
943 RETURN(1) ; /* get to here is good */
944
945 /*-------------------------------------------------------------------*/
946 /* create a shared memory segment, then setup sub-bricks inside that */
947 /*-------------------------------------------------------------------*/
948
949 case DATABLOCK_MEM_SHARED:
950 #if !defined(DONT_USE_SHM) && !defined(CYGWIN)
951 { if( blk->shm_idcode[0] == '\0' ){ /* new segment */
952 UNIQ_idcode_fill( blk->shm_idcode ) ;
953 blk->shm_idint = shm_create( blk->shm_idcode , (int)blk->total_bytes ) ;
954 if( blk->shm_idint < 0 ){ blk->shm_idcode[0] = '\0'; RETURN(0); }
955 ptr = shm_attach( blk->shm_idint ) ;
956 if( ptr == NULL ){
957 shmctl( blk->shm_idint , IPC_RMID , NULL ) ;
958 blk->shm_idint = -1; blk->shm_idcode[0] = '\0'; RETURN(0);
959 }
960 for( ibr=0 ; ibr < nv ; ibr++ ){
961 mri_fix_data_pointer( ptr , DBLK_BRICK(blk,ibr) ) ;
962 ptr += DBLK_BRICK_BYTES(blk,ibr) ;
963 }
964 }
965 }
966 RETURN(1) ;
967 #else
968 RETURN(0) ;
969 #endif
970
971 }
972
973 RETURN(0) ; /* should never get to here */
974 }
975
976 /*----------------------------------------------------------------------------*/
977 /*! Fill up a dataset with zero-ed out sub-bricks, if needed.
978 *//*--------------------------------------------------------------------------*/
979
THD_zerofill_dataset(THD_3dim_dataset * dset)980 void THD_zerofill_dataset( THD_3dim_dataset *dset )
981 {
982 int ii ;
983 void *vpt ;
984
985 ENTRY("THD_zerofill_dataset") ;
986
987 if( !ISVALID_DSET(dset) || !ISVALID_DATABLOCK(dset->dblk) ) EXRETURN ;
988
989 for( ii=0 ; ii < DSET_NVALS(dset) ; ii++ ){
990 if( DSET_ARRAY(dset,ii) == NULL ){
991 vpt = calloc(1,DSET_BRICK_BYTES(dset,ii)) ;
992 mri_fix_data_pointer( vpt , DSET_BRICK(dset,ii) ) ;
993 }
994 }
995 EXRETURN ;
996 }
997
998 /*----------------------------------------------------------------------------*/
999 /*! Apply master limits to data sub-bricks. 14 Apr 2006 [rickr]
1000 (from THD_load_datablock())
1001
1002 return 0 on success
1003 *//*--------------------------------------------------------------------------*/
THD_apply_master_subrange(THD_datablock * blk)1004 int THD_apply_master_subrange( THD_datablock * blk )
1005 {
1006 THD_diskptr * dkptr ;
1007 float bot = blk->master_bot , top = blk->master_top ;
1008 int jbr, nv, ii, nxyz ;
1009
1010 ENTRY("THD_apply_master_limits") ;
1011
1012 if( ! DBLK_IS_MASTERED(blk) )
1013 RETURN(0);
1014
1015 /* if set, process csv list instead of bot..top 30 Nov 2016 [rickr] */
1016 if( (blk->master_ncsv > 0) && (blk->master_csv != NULL) )
1017 RETURN(THD_apply_master_subrange_list(blk));
1018
1019 /* if not applicable, skip bot..top range */
1020 if( blk->master_bot > blk->master_top )
1021 RETURN(0);
1022
1023 dkptr = blk->diskptr ;
1024 nv = dkptr->nvals ;
1025 nxyz = dkptr->dimsizes[0] * dkptr->dimsizes[1] * dkptr->dimsizes[2];
1026
1027 for( jbr=0 ; jbr < nv ; jbr++ ){
1028 switch( DBLK_BRICK_TYPE(blk,jbr) ){
1029
1030 default:
1031 fprintf(stderr,
1032 "** Can't sub-range datum type %s!\n",
1033 MRI_TYPE_name[DBLK_BRICK_TYPE(blk,jbr)]) ;
1034 RETURN(1);
1035 break ;
1036
1037 case MRI_short:{
1038 short mbot, mtop, *mar = (short *) DBLK_ARRAY(blk,jbr) ;
1039 float mfac = DBLK_BRICK_FACTOR(blk,jbr) ;
1040 float fval, fvmax;
1041 if( mfac == 0.0 ) mfac = 1.0 ;
1042 /* - do not use SHORTIZE, rounding can include unwanted values */
1043 /* - use ceil() for bot and floor() for top 21 Nov 2016 [rickr] */
1044 fval = bot/mfac;
1045 fvmax = MRI_TYPE_maxval[MRI_short];
1046 if ( fval < -fvmax ) mbot = (short)-fvmax;
1047 else if( fval > fvmax ) mbot = (short)fvmax;
1048 else mbot = (short)ceilf(fval);
1049
1050 fval = top/mfac;
1051 if ( fval < -fvmax ) mtop = (short)-fvmax;
1052 else if( fval > fvmax ) mtop = (short)fvmax;
1053 else mtop = (short)floorf(fval);
1054 /* mbot = SHORTIZE(bot/mfac) ; mtop = SHORTIZE(top/mfac) ; */
1055 #if 0
1056 fprintf(stderr,"bot=%f top=%f\n",bot,top) ;
1057 fprintf(stderr,"mbot=%d mtop=%d\n",(int)mbot,(int)mtop) ;
1058 #endif
1059 for( ii=0 ; ii < nxyz ; ii++ )
1060 if( mar[ii] < mbot || mar[ii] > mtop ) mar[ii] = 0 ;
1061 }
1062 break ;
1063
1064 case MRI_int:{
1065 int mbot, mtop, *mar = (int *) DBLK_ARRAY(blk,jbr) ;
1066 float mfac = DBLK_BRICK_FACTOR(blk,jbr) ;
1067 if( mfac == 0.0 ) mfac = 1.0 ;
1068 /* do not include unrequested values 21 Nov 2016 [rickr] */
1069 mbot = ceilf(bot/mfac) ; mtop = floor(top/mfac) ;
1070 for( ii=0 ; ii < nxyz ; ii++ )
1071 if( mar[ii] < mbot || mar[ii] > mtop ) mar[ii] = 0 ;
1072 }
1073 break ;
1074
1075 case MRI_byte:{
1076 byte mbot, mtop, *mar = (byte *) DBLK_ARRAY(blk,jbr) ;
1077 float fval, fvmax, mfac = DBLK_BRICK_FACTOR(blk,jbr) ;
1078 if( mfac == 0.0 ) mfac = 1.0 ;
1079 fval = bot/mfac;
1080 fvmax = MRI_TYPE_maxval[MRI_byte];
1081 if ( fval < 0 ) mbot = (byte)0;
1082 else if( fval > fvmax ) mbot = (byte)fvmax;
1083 else mbot = (byte)ceilf(fval);
1084
1085 fval = top/mfac;
1086 if ( fval < 0 ) mtop = (byte)0;
1087 else if( fval > fvmax ) mtop = (byte)fvmax;
1088 else mtop = (byte)floorf(fval);
1089
1090 /* old way: mbot = BYTEIZE(bot/mfac) ; mtop = BYTEIZE(top/mfac) ; */
1091 for( ii=0 ; ii < nxyz ; ii++ )
1092 if( mar[ii] < mbot || mar[ii] > mtop ) mar[ii] = 0 ;
1093 }
1094 break ;
1095
1096 case MRI_float:{
1097 float mbot, mtop, *mar = (float *) DBLK_ARRAY(blk,jbr) ;
1098 float mfac = DBLK_BRICK_FACTOR(blk,jbr) ;
1099 if( mfac == 0.0 ) mfac = 1.0 ;
1100 mbot = (bot/mfac) ; mtop = (top/mfac) ;
1101 for( ii=0 ; ii < nxyz ; ii++ )
1102 if( mar[ii] < mbot || mar[ii] > mtop ) mar[ii] = 0 ;
1103 }
1104 break ;
1105
1106 case MRI_complex:{
1107 float mbot, mtop , val ;
1108 complex *mar = (complex *) DBLK_ARRAY(blk,jbr) ;
1109 float mfac = DBLK_BRICK_FACTOR(blk,jbr) ;
1110 if( mfac == 0.0 ) mfac = 1.0 ;
1111 mbot = (bot/mfac) ; mtop = (top/mfac) ;
1112 for( ii=0 ; ii < nxyz ; ii++ ){
1113 val = CABS(mar[ii]) ;
1114 if( val < mbot || val > mtop ) mar[ii].r = mar[ii].i = 0 ;
1115 }
1116 }
1117 break ;
1118 }
1119 }
1120
1121 RETURN(0) ;
1122 }
1123
1124 /*----------------------------------------------------------------------------*/
1125 /*! Apply master limit lists to data sub-bricks. 30 Nov 2016 [rickr]
1126
1127 Like THD_apply_master_subrange(), but require all values to have
1128 corresponding entries in master_csv list.
1129
1130 This will be a slower function.
1131
1132 Perhaps this should all go directly through floats? Perhaps generalize
1133 to list of float ranges? It would not be difficult.
1134
1135 return 0 on success
1136 *//*--------------------------------------------------------------------------*/
THD_apply_master_subrange_list(THD_datablock * blk)1137 int THD_apply_master_subrange_list( THD_datablock * blk )
1138 {
1139 THD_diskptr * dkptr ;
1140 float * fcsv=NULL;
1141 int jbr, ii, ic, ncsv, nxyz ;
1142
1143 ENTRY("THD_apply_master_limits") ;
1144
1145 if( ! DBLK_IS_MASTERED(blk) || blk->master_ncsv <= 0 ||
1146 blk->master_csv == NULL )
1147 RETURN(0);
1148
1149 /* should we warn if DBLK_BRICK_FACTOR() != 0.0 or 1.0? */
1150
1151 /* do all work as floats, to ease issues with brick factors */
1152 /* (unless we do not allow factors, and then re-write this) */
1153 ncsv = blk->master_ncsv;
1154 fcsv = (float *)malloc(ncsv * sizeof(float));
1155 for( ii=0; ii < ncsv; ii++)
1156 fcsv[ii] = (float)blk->master_csv[ii];
1157
1158 dkptr = blk->diskptr ;
1159 nxyz = dkptr->dimsizes[0] * dkptr->dimsizes[1] * dkptr->dimsizes[2];
1160
1161 for( jbr=0 ; jbr < dkptr->nvals ; jbr++ ){
1162 switch( DBLK_BRICK_TYPE(blk,jbr) ){
1163
1164 default:
1165 fprintf(stderr, "** TAMSL: cannot sub-range datum type %s\n",
1166 MRI_TYPE_name[DBLK_BRICK_TYPE(blk,jbr)]) ;
1167 RETURN(1);
1168 break ;
1169
1170 case MRI_short:{
1171 short *mar = (short *) DBLK_ARRAY(blk,jbr) ;
1172 float dval, mfac = DBLK_BRICK_FACTOR(blk,jbr) ;
1173 if( mfac == 0.0 ) mfac = 1.0 ;
1174 for( ii=0 ; ii < nxyz ; ii++ ) {
1175 /* write as function? overhead per voxel, hmmmm... */
1176 dval = mar[ii] * mfac;
1177 for( ic=0; ic < ncsv; ic++ )
1178 if( dval == fcsv[ic] ) break;
1179 /* if not found, clear */
1180 if( ic == ncsv ) mar[ii] = 0 ;
1181 }
1182 }
1183 break ;
1184
1185 /* going through float limits to 24 of 32 bits, ~16 million */
1186 case MRI_int:{
1187 int *mar = (int *) DBLK_ARRAY(blk,jbr) ;
1188 float dval, mfac = DBLK_BRICK_FACTOR(blk,jbr) ;
1189 if( mfac == 0.0 ) mfac = 1.0 ;
1190 for( ii=0 ; ii < nxyz ; ii++ ) {
1191 dval = mar[ii] * mfac;
1192 for( ic=0; ic < ncsv; ic++ )
1193 if( dval == fcsv[ic] ) break;
1194 if( ic == ncsv ) mar[ii] = 0 ; /* if not found, clear */
1195 }
1196 }
1197 break ;
1198
1199 case MRI_byte:{
1200 byte *mar = (byte *) DBLK_ARRAY(blk,jbr) ;
1201 float dval, mfac = DBLK_BRICK_FACTOR(blk,jbr) ;
1202 if( mfac == 0.0 ) mfac = 1.0 ;
1203 for( ii=0 ; ii < nxyz ; ii++ ) {
1204 dval = mar[ii] * mfac;
1205 for( ic=0; ic < ncsv; ic++ )
1206 if( dval == fcsv[ic] ) break;
1207 if( ic == ncsv ) mar[ii] = 0 ; /* if not found, clear */
1208 }
1209 }
1210 break ;
1211
1212 case MRI_float:{
1213 float *mar = (float *) DBLK_ARRAY(blk,jbr) ;
1214 float dval, mfac = DBLK_BRICK_FACTOR(blk,jbr) ;
1215 if( mfac == 0.0 ) mfac = 1.0 ;
1216 for( ii=0 ; ii < nxyz ; ii++ ) {
1217 dval = mar[ii] * mfac;
1218 for( ic=0; ic < ncsv; ic++ )
1219 if( dval == fcsv[ic] ) break;
1220 if( ic == ncsv ) mar[ii] = 0 ; /* if not found, clear */
1221 }
1222 }
1223 break ;
1224
1225 case MRI_complex:{
1226 complex *mar = (complex *) DBLK_ARRAY(blk,jbr) ;
1227 float dval, mfac = DBLK_BRICK_FACTOR(blk,jbr) ;
1228 if( mfac == 0.0 ) mfac = 1.0 ;
1229 for( ii=0 ; ii < nxyz ; ii++ ){
1230 dval = CABS(mar[ii]) ;
1231 for( ic=0; ic < ncsv; ic++ )
1232 if( dval == fcsv[ic] ) break;
1233 /* if not found, clear */
1234 if( ic == ncsv ) mar[ii].r = mar[ii].i = 0 ;
1235 }
1236 }
1237 break ;
1238 }
1239 }
1240
1241 free(fcsv);
1242
1243 RETURN(0) ;
1244 }
1245
1246 /*----------------------------------------------------------------------*/
1247 /*! Set dx,dy,dz fields for MRI_IMAGEs in a dataset.
1248 [19 Dec 2011] Also patch the dataset grid spacing if it is 0.
1249 *//*--------------------------------------------------------------------*/
1250
THD_patch_brickim(THD_3dim_dataset * dset)1251 void THD_patch_brickim( THD_3dim_dataset *dset )
1252 {
1253 float dx,dy,dz , dm ;
1254 int iv , nvals , nfix=0 ;
1255 static char *cfix[8] = { NULL ,
1256 "x-axis" , "y-axis" , "x- & y-axes" ,
1257 "z-axis" , "x- & z-axes" , "y- & z-axes" ,
1258 "x- & y- & z-axes" } ;
1259 static char **idstr=NULL ; static int nidstr=0 ; /* no repeating */
1260
1261 ENTRY("THD_patch_brickim") ;
1262
1263 if( !ISVALID_DSET(dset) ) EXRETURN ;
1264
1265 dx = fabsf(DSET_DX(dset)) ; /* unsigned grid spacings */
1266 dy = fabsf(DSET_DY(dset)) ;
1267 dz = fabsf(DSET_DZ(dset)) ;
1268
1269 /* check to see if any axis spacing was 0; if so, spackle over it */
1270
1271 dm = dx + dy + dz ;
1272 dm = (dm == 0.0f) ? 1.0f : dm*0.5f ;
1273
1274 if( dx == 0.0f ){ dx = dset->daxes->xxdel = dm ; nfix += 1 ; }
1275 if( dy == 0.0f ){ dy = dset->daxes->yydel = dm ; nfix += 2 ; }
1276 if( dz == 0.0f ){ dz = dset->daxes->zzdel = dm ; nfix += 4 ; }
1277
1278 if( nfix > 0 ){ /* have to patch at least one axis */
1279 int ii ;
1280 for( ii=0 ; ii < nidstr ; ii++ ) /* see if we did this dataset before */
1281 if( strcmp(DSET_IDCODE_STR(dset),idstr[ii]) == 0 ) break ;
1282 if( ii == nidstr ){ /* this is a new dataset (to us) */
1283 if( nidstr == 0 ) fprintf(stderr,"\n") ;
1284 WARNING_message("Dataset %s : patched zero grid spacing along %s to %g",
1285 THD_trailname(DSET_HEADNAME(dset),0) , cfix[nfix] , dm ) ;
1286 idstr = (char **)realloc(idstr,sizeof(char *)*(nidstr+1)) ;
1287 idstr[nidstr++] = strdup(DSET_IDCODE_STR(dset)) ; /* remember the past */
1288 }
1289 }
1290
1291 /* fix the sub-brick (MRI_IMAGE) headers */
1292
1293 nvals = DSET_NVALS(dset) ;
1294 for( iv=0 ; iv < nvals ; iv++ ){
1295 DSET_BRICK(dset,iv)->dx = dx ;
1296 DSET_BRICK(dset,iv)->dy = dy ;
1297 DSET_BRICK(dset,iv)->dz = dz ;
1298 }
1299
1300 EXRETURN ;
1301 }
1302