00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008 #include "thd.h"
00009
00010 static int native_order = -1 ;
00011 static int no_mmap = -1 ;
00012 static int floatscan = -1 ;
00013
00014 #define PRINT_SIZE 66600000
00015 #define PRINT_STEP 10
00016
00017 static int verbose = 0 ;
00018
00019 void THD_load_datablock_verbose( int v ){ verbose = v; }
00020 int THD_alloc_datablock( THD_datablock *blk ) ;
00021
00022
00023
00024
00025
00026 int THD_datum_constant( THD_datablock *blk )
00027 {
00028 int ibr , dzero , nv=blk->nvals ;
00029 if( nv == 1 ) return 1 ;
00030 dzero = DBLK_BRICK_TYPE(blk,0) ;
00031 for( ibr=1 ; ibr < nv ; ibr++ )
00032 if( dzero != DBLK_BRICK_TYPE(blk,ibr) ) return 0 ;
00033 return 1 ;
00034 }
00035
00036
00037
00038
00039
00040 float THD_get_voxel( THD_3dim_dataset *dset , int ijk , int ival )
00041 {
00042 void *ar ;
00043 float val , fac ;
00044
00045 if( ival < 0 || ival >= DSET_NVALS(dset) ) return 0.0f ;
00046 if( ijk < 0 || ijk >= DSET_NVOX(dset) ) return 0.0f ;
00047
00048 ar = DSET_ARRAY(dset,ival) ;
00049 if( ar == NULL ){ DSET_load(dset) ; ar = DSET_ARRAY(dset,ival) ; }
00050 if( ar == NULL ) return 0.0f ;
00051
00052 switch( DSET_BRICK_TYPE(dset,ival) ){
00053
00054 default: return 0.0f ;
00055
00056 case MRI_byte:
00057 val = (float)(((byte *)ar)[ijk]) ; break ;
00058 case MRI_short:
00059 val = (float)(((short *)ar)[ijk]) ; break ;
00060 case MRI_int:
00061 val = (float)(((int *)ar)[ijk]) ; break ;
00062 case MRI_float:
00063 val = (float)(((float *)ar)[ijk]) ; break ;
00064 case MRI_double:
00065 val = (float)(((double *)ar)[ijk]) ; break ;
00066
00067 case MRI_complex:{
00068 complex c = (((complex *)ar)[ijk]) ;
00069 val = sqrt(c.r*c.r+c.i*c.i) ;
00070 break ;
00071 }
00072
00073 case MRI_rgb:{
00074 rgbyte c = (((rgbyte *)ar)[ijk]) ;
00075 val = 0.299f*(float)(c.r) + 0.587f*(float)(c.g) + 0.114f*(float)(c.b) ;
00076 break ;
00077 }
00078
00079 case MRI_rgba:{
00080 rgba c = (((rgba *)ar)[ijk]) ;
00081 val = 0.299f*(float)(c.r) + 0.587f*(float)(c.g) + 0.114f*(float)(c.b) ;
00082 val *= 0.00392157f*(float)(c.a) ;
00083 break ;
00084 }
00085 }
00086
00087 fac = DSET_BRICK_FACTOR(dset,ival) ;
00088 if( fac > 0.0f ) val *= fac ;
00089 return val ;
00090 }
00091
00092
00093
00094
00095
00096
00097
00098 static generic_func *freeup=NULL ;
00099
00100 void THD_set_freeup( generic_func *ff ){ freeup = ff; }
00101
00102
00103
00104 Boolean THD_load_datablock( THD_datablock *blk )
00105 {
00106 THD_diskptr *dkptr ;
00107 int id , offset ;
00108 int nx,ny,nz , nxy,nxyz , nv,vv , ii , ibr , nbad ;
00109 char *ptr ;
00110 MRI_IMAGE *im ;
00111 int verb=verbose ;
00112 int64_t idone , print_size=PRINT_SIZE ;
00113
00114 ENTRY("THD_load_datablock") ;
00115
00116 if( native_order < 0 ) native_order = mri_short_order() ;
00117
00118 floatscan = AFNI_yesenv("AFNI_FLOATSCAN") ;
00119
00120 if( floatscan ) no_mmap = 1 ;
00121 else no_mmap = AFNI_yesenv("AFNI_NOMMAP") ;
00122
00123
00124
00125 if( ! ISVALID_DATABLOCK(blk) || blk->brick == NULL ){
00126 STATUS("Illegal inputs"); RETURN( False );
00127 }
00128
00129 ii = THD_count_databricks( blk ) ;
00130 if( ii == blk->nvals ) RETURN( True );
00131
00132 if( blk->malloc_type == DATABLOCK_MEM_UNDEFINED ){
00133 STATUS("malloc_type forced to DATABLOCK_MEM_MALLOC") ;
00134 blk->malloc_type = DATABLOCK_MEM_MALLOC ;
00135 }
00136
00137
00138
00139 dkptr = blk->diskptr ;
00140 if( ! ISVALID_DISKPTR(dkptr) ){
00141 STATUS("invalid dkptr!!!"); RETURN( False );
00142 }
00143 if( dkptr->storage_mode == STORAGE_UNDEFINED ){
00144 if( PRINT_TRACING ){
00145 char str[512] ; THD_3dim_dataset *dset=(THD_3dim_dataset *)(blk->parent) ;
00146 sprintf(str,"dataset %s == STORAGE_UNDEFINED",DSET_BRIKNAME(dset)) ;
00147 STATUS(str) ;
00148 }
00149 RETURN(False) ;
00150 }
00151
00152 if( dkptr->rank != 3 ){
00153 fprintf(stderr,"\n*** Cannot read non 3D datablocks ***\n"); RETURN(False);
00154 }
00155
00156 if( dkptr->storage_mode == STORAGE_BY_VOLUMES ) no_mmap = 1 ;
00157
00158
00159
00160 if( dkptr->storage_mode == STORAGE_BY_MINC ){
00161 THD_load_minc( blk ) ;
00162 ii = THD_count_databricks( blk ) ;
00163 if( ii == blk->nvals ){
00164 THD_update_statistics( (THD_3dim_dataset *)blk->parent ) ;
00165 RETURN( True ) ;
00166 }
00167 STATUS("can't read MINC file?!") ;
00168 RETURN( False ) ;
00169 }
00170
00171 { THD_3dim_dataset *ds = (THD_3dim_dataset *)blk->parent ;
00172 if( ISVALID_DSET(ds) && DSET_IS_TCAT(ds) ){
00173 THD_load_tcat( blk ) ;
00174 ii = THD_count_databricks( blk ) ;
00175 if( ii == blk->nvals ){
00176 THD_update_statistics( (THD_3dim_dataset *)blk->parent ) ;
00177 RETURN( True ) ;
00178 }
00179 STATUS("can't read tcat-ed file?!") ;
00180 RETURN( False ) ;
00181 }
00182 }
00183
00184 if( dkptr->storage_mode == STORAGE_BY_ANALYZE ){
00185 THD_load_analyze( blk ) ;
00186 ii = THD_count_databricks( blk ) ;
00187 if( ii == blk->nvals ){
00188 THD_update_statistics( (THD_3dim_dataset *)blk->parent ) ;
00189 RETURN( True ) ;
00190 }
00191 STATUS("can't read ANALYZE file?!") ;
00192 RETURN( False ) ;
00193 }
00194
00195 if( dkptr->storage_mode == STORAGE_BY_CTFMRI ){
00196 THD_load_ctfmri( blk ) ;
00197 ii = THD_count_databricks( blk ) ;
00198 if( ii == blk->nvals ){
00199 THD_update_statistics( (THD_3dim_dataset *)blk->parent ) ;
00200 RETURN( True ) ;
00201 }
00202 STATUS("can't read CTF MRI file?!") ;
00203 RETURN( False ) ;
00204 }
00205
00206 if( dkptr->storage_mode == STORAGE_BY_CTFSAM ){
00207 THD_load_ctfsam( blk ) ;
00208 ii = THD_count_databricks( blk ) ;
00209 if( ii == blk->nvals ){
00210 THD_update_statistics( (THD_3dim_dataset *)blk->parent ) ;
00211 RETURN( True ) ;
00212 }
00213 STATUS("can't read CTF SAM file?!") ;
00214 RETURN( False ) ;
00215 }
00216
00217 if( dkptr->storage_mode == STORAGE_BY_1D ){
00218 THD_load_1D( blk ) ;
00219 ii = THD_count_databricks( blk ) ;
00220 if( ii == blk->nvals ){
00221 THD_update_statistics( (THD_3dim_dataset *)blk->parent ) ;
00222 RETURN( True ) ;
00223 }
00224 STATUS("can't read 1D dataset file?!") ;
00225 RETURN( False ) ;
00226 }
00227
00228 if( dkptr->storage_mode == STORAGE_BY_3D ){
00229 THD_load_3D( blk ) ;
00230 ii = THD_count_databricks( blk ) ;
00231 if( ii == blk->nvals ){
00232 THD_update_statistics( (THD_3dim_dataset *)blk->parent ) ;
00233 RETURN( True ) ;
00234 }
00235 STATUS("can't read 3D dataset file?!") ;
00236 RETURN( False ) ;
00237 }
00238
00239 if( dkptr->storage_mode == STORAGE_BY_NIFTI ){
00240 THD_load_nifti( blk ) ;
00241 ii = THD_count_databricks( blk ) ;
00242 if( ii == blk->nvals ){
00243 THD_update_statistics( (THD_3dim_dataset *)blk->parent ) ;
00244 RETURN( True ) ;
00245 }
00246 STATUS("can't read NIFTI dataset file?!") ;
00247 RETURN( False ) ;
00248 }
00249
00250 if( dkptr->storage_mode == STORAGE_BY_MPEG ){
00251 THD_load_mpeg( blk ) ;
00252 ii = THD_count_databricks( blk ) ;
00253 if( ii == blk->nvals ){
00254 THD_update_statistics( (THD_3dim_dataset *)blk->parent ) ;
00255 RETURN( True ) ;
00256 }
00257 STATUS("can't read MPEG dataset file?!") ;
00258 RETURN( False ) ;
00259 }
00260
00261
00262
00263
00264
00265 nx = dkptr->dimsizes[0] ;
00266 ny = dkptr->dimsizes[1] ; nxy = nx * ny ;
00267 nz = dkptr->dimsizes[2] ; nxyz = nxy * nz ;
00268 nv = dkptr->nvals ;
00269
00270 if( DBLK_IS_MASTERED(blk) && blk->malloc_type == DATABLOCK_MEM_MMAP )
00271 blk->malloc_type = DATABLOCK_MEM_MALLOC ;
00272
00273
00274
00275 if( !THD_datum_constant(blk) && blk->malloc_type != DATABLOCK_MEM_MALLOC ){
00276 fprintf(stderr,"++ WARNING: dataset %s: non-uniform sub-brick types\n",
00277 blk->diskptr->filecode) ;
00278 blk->malloc_type = DATABLOCK_MEM_MALLOC ;
00279 }
00280
00281
00282
00283 if( dkptr->byte_order <= 0 ) dkptr->byte_order = native_order ;
00284
00285
00286
00287
00288
00289 if( dkptr->byte_order != native_order ){
00290 for( ii=0 ; ii < nv ; ii++ )
00291 if( DBLK_BRICK_TYPE(blk,ii) != MRI_byte &&
00292 DBLK_BRICK_TYPE(blk,ii) != MRI_rgb ) break ;
00293 if( ii == nv ) dkptr->byte_order = native_order ;
00294 }
00295
00296
00297
00298 if( dkptr->byte_order != native_order || no_mmap ){
00299 if( blk->malloc_type == DATABLOCK_MEM_MMAP )
00300 blk->malloc_type = DATABLOCK_MEM_MALLOC ;
00301 }
00302
00303 DBLK_mmapfix(blk) ;
00304
00305
00306
00307 if( blk->malloc_type == DATABLOCK_MEM_MALLOC ||
00308 blk->malloc_type == DATABLOCK_MEM_SHARED ){
00309
00310 ii = THD_alloc_datablock( blk ) ;
00311 if( ii == 0 ) RETURN(False) ;
00312
00313
00314
00315
00316 } else if( blk->malloc_type == DATABLOCK_MEM_MMAP ){
00317
00318 int fd , fsize ;
00319 fd = open( dkptr->brick_name , O_RDONLY ) ;
00320 if( fd < 0 ){
00321 fprintf( stderr , "\n*** cannot open brick file %s for mmap\n"
00322 " - do you have permission? does file exist?\n" ,
00323 dkptr->brick_name ) ;
00324 perror("*** Unix error message") ;
00325 STATUS("open failed") ;
00326 RETURN( False );
00327 }
00328
00329
00330
00331 fsize = THD_filesize( dkptr->brick_name ) ;
00332 if( fsize < blk->total_bytes )
00333 fprintf(stderr ,
00334 "\n*** WARNING: file %s size is %d, but should be at least %lld!\n" ,
00335 dkptr->brick_name , fsize , blk->total_bytes ) ;
00336
00337
00338
00339 for( ibr=0 ; ibr < nv ; ibr++ )
00340 mri_clear_data_pointer( DBLK_BRICK(blk,ibr) ) ;
00341
00342
00343
00344 ptr = (char *) mmap( 0 , (size_t)blk->total_bytes ,
00345 PROT_READ , THD_MMAP_FLAG , fd , 0 ) ;
00346
00347
00348
00349 if( ptr == (char *)(-1) ){
00350 fprintf(stderr ,
00351 "\n*** cannot mmap brick file %s - maybe hit a system limit?\n" ,
00352 dkptr->brick_name ) ;
00353 perror("*** Unix error message") ;
00354 if( freeup != NULL ){
00355 fprintf(stderr,"*** trying to fix problem\n") ;
00356 freeup() ;
00357 ptr = (char *) mmap( 0 , (size_t)blk->total_bytes ,
00358 PROT_READ , THD_MMAP_FLAG , fd , 0 ) ;
00359 if( ptr == (char *)(-1) ){
00360 fprintf(stderr,"*** cannot fix problem!\n") ;
00361 close(fd) ;
00362 STATUS("freeup failed") ; RETURN( False );
00363 } else {
00364 fprintf(stderr,"*** was able to fix problem!\n") ;
00365 }
00366 } else {
00367 close(fd) ;
00368 STATUS("mmap failed") ; RETURN( False );
00369 }
00370 }
00371
00372 close(fd) ;
00373
00374
00375
00376 offset = 0 ;
00377 for( ibr=0 ; ibr < nv ; ibr++ ){
00378 mri_fix_data_pointer( ptr , DBLK_BRICK(blk,ibr) ) ;
00379 ptr += DBLK_BRICK_BYTES(blk,ibr) ;
00380 }
00381
00382 STATUS("mmap succeeded") ;
00383 RETURN( True );
00384
00385 } else {
00386 STATUS("unknown malloc_type code?!") ;
00387 RETURN( False ) ;
00388 }
00389
00390
00391
00392
00393 ptr = getenv("AFNI_LOAD_PRINTSIZE") ;
00394 if( verb && ptr != NULL ){
00395 char *ept ;
00396 idone = strtol( ptr , &ept , 10 ) ;
00397 if( idone > 0 ){
00398 if( *ept == 'K' || *ept == 'k' ) id *= 1024 ;
00399 else if( *ept == 'M' || *ept == 'm' ) id *= 1024*1024 ;
00400 print_size = id ;
00401 } else {
00402 print_size = 666000000 ;
00403 }
00404 }
00405
00406 if( verb ) verb = (blk->total_bytes > print_size) ;
00407 if( verb ) fprintf(stderr,"reading dataset %s",dkptr->filecode) ;
00408
00409 switch( dkptr->storage_mode ){
00410
00411
00412
00413 default:
00414 fprintf(stderr,"\n*** illegal storage mode in read ***\n") ;
00415 RETURN( False );
00416 break ;
00417
00418
00419
00420 case STORAGE_BY_BRICK:{
00421 FILE *far ;
00422
00423 STATUS("reading from BRIK file") ;
00424
00425 far = COMPRESS_fopen_read( dkptr->brick_name ) ;
00426
00427 if( far == NULL ){
00428 THD_purge_datablock( blk , blk->malloc_type ) ;
00429 fprintf(stderr,
00430 "\n*** failure while opening brick file %s "
00431 "- do you have permission?\n" ,
00432 dkptr->brick_name ) ;
00433 perror("*** Unix error message") ;
00434 RETURN( False );
00435 }
00436
00437
00438
00439 idone = 0 ;
00440 if( ! DBLK_IS_MASTERED(blk) ){
00441
00442 for( ibr=0 ; ibr < nv ; ibr++ ){
00443 idone += fread( DBLK_ARRAY(blk,ibr), 1,
00444 DBLK_BRICK_BYTES(blk,ibr), far ) ;
00445
00446 if( verb && ibr%PRINT_STEP == 0 ) fprintf(stderr,".") ;
00447 }
00448
00449 } else {
00450
00451 int nfilled = 0 , nbuf=0 , jbr, nbr ;
00452 char *buf=NULL ;
00453
00454
00455
00456
00457
00458 for( ibr=0 ; nfilled < nv && ibr < blk->master_nvals ; ibr++ ){
00459
00460 if( nbuf < blk->master_bytes[ibr] ){
00461 if( buf != NULL ) free(buf) ;
00462 nbuf = blk->master_bytes[ibr] ;
00463 buf = AFMALL(char, sizeof(char) * nbuf ) ;
00464 if( buf == NULL ) break ;
00465 }
00466
00467
00468
00469 nbr = fread( buf , 1 , blk->master_bytes[ibr] , far ) ;
00470 if( verb && ibr%PRINT_STEP == 0 ) fprintf(stderr,".") ;
00471 if( nbr < blk->master_bytes[ibr] ) break ;
00472
00473
00474
00475 for( jbr=0 ; jbr < nv ; jbr++ ){
00476 if( blk->master_ival[jbr] == ibr ){
00477 memcpy( DBLK_ARRAY(blk,jbr) , buf , blk->master_bytes[ibr] ) ;
00478 nfilled++ ;
00479 idone += nbr ;
00480 }
00481 }
00482 }
00483
00484 if( buf != NULL ) free(buf) ;
00485 }
00486
00487
00488
00489 COMPRESS_fclose(far) ;
00490
00491
00492
00493 if( idone != blk->total_bytes ){
00494 THD_purge_datablock( blk , blk->malloc_type ) ;
00495 fprintf(stderr ,
00496 "\n*** failure while reading from brick file %s\n"
00497 "*** desired %lld bytes but only got %lld\n" ,
00498 dkptr->brick_name , blk->total_bytes , idone ) ;
00499 perror("*** Unix error message") ;
00500 RETURN( False );
00501 }
00502 }
00503 break ;
00504
00505
00506
00507 case STORAGE_BY_VOLUMES:{
00508 ATR_string *atr ;
00509 char **fnam , *ptr , *flist ;
00510 FILE *far ;
00511
00512 STATUS("reading from volume files") ;
00513
00514
00515
00516 atr = THD_find_string_atr(blk,"VOLUME_FILENAMES") ;
00517
00518
00519
00520 if( atr == NULL || atr->nch < 2*nv ){
00521 THD_purge_datablock( blk , blk->malloc_type ) ;
00522 fprintf(stderr,
00523 "Dataset %s does not have legal VOLUME_FILENAMES attribute!\n",
00524 blk->diskptr->filecode) ;
00525 RETURN( False );
00526 }
00527
00528
00529
00530 flist = AFMALL(char, atr->nch+1) ;
00531 memcpy( flist , atr->ch , atr->nch ) ;
00532 flist[atr->nch] = '\0' ;
00533
00534 #if 0
00535 fprintf(stderr,"VOL: flist:\n%s\n",flist) ;
00536 #endif
00537
00538
00539
00540 fnam = (char **) malloc(sizeof(char *)*nv) ;
00541 for( ibr=0 ; ibr < nv ; ibr++ ){
00542
00543
00544
00545 ptr = strtok( (ibr==0)?flist:NULL , " \t\n\r\f\v" ) ;
00546
00547
00548
00549 if( ptr == NULL ){
00550 fprintf(stderr,
00551 "Dataset %s has illegal VOLUME_FILENAMES attribute!\n",
00552 blk->diskptr->filecode) ;
00553 for( ii=0 ; ii < nv ; ii++ ){
00554 free( DBLK_ARRAY(blk,ii) ) ;
00555 mri_clear_data_pointer( DBLK_BRICK(blk,ii) ) ;
00556 }
00557 for( ii=0 ; ii < ibr ; ii++ ) free(fnam[ii]) ;
00558 free(flist) ; free(fnam) ;
00559 RETURN( False );
00560 }
00561
00562
00563
00564 #if 0
00565 fprintf(stderr,"VOL[%d]: ptr=%s\n",ibr,ptr) ;
00566 #endif
00567
00568 if( *ptr == '/' ){
00569 fnam[ibr] = strdup(ptr) ;
00570 } else {
00571 ii = strlen(blk->diskptr->directory_name) + strlen(ptr) ;
00572 fnam[ibr] = AFMALL(char,ii+1) ; fnam[ibr][0] = '\0' ;
00573 strcat(fnam[ibr],blk->diskptr->directory_name) ;
00574 strcat(fnam[ibr],ptr) ;
00575 }
00576 #if 0
00577 fprintf(stderr,"VOL[%d]: fnam=%s\n",ibr,fnam[ibr]) ;
00578 #endif
00579 }
00580 free(flist) ;
00581
00582
00583
00584 if( ! DBLK_IS_MASTERED(blk) ){
00585
00586 for( ibr=0 ; ibr < nv ; ibr++ ){
00587
00588 #if 0
00589 fprintf(stderr,"VOL[%d]: opening %s\n",ibr,fnam[ibr]) ;
00590 #endif
00591
00592 far = COMPRESS_fopen_read( fnam[ibr] ) ;
00593
00594 if( far == NULL ){
00595 fprintf(stderr,
00596 "\n*** Failure while opening volume file %s "
00597 "- do you have permission?\n" ,
00598 fnam[ibr] ) ;
00599 perror("*** Unix error message") ;
00600 THD_purge_datablock( blk , blk->malloc_type ) ;
00601 for( ii=0 ; ii < nv ; ii++ ) free(fnam[ii]) ;
00602 free(fnam); RETURN( False );
00603 }
00604
00605
00606
00607 id = fread( DBLK_ARRAY(blk,ibr), 1,
00608 DBLK_BRICK_BYTES(blk,ibr), far ) ;
00609 if( verb && ibr%PRINT_STEP==0 ) fprintf(stderr,".") ;
00610
00611 #if 0
00612 fprintf(stderr,"VOL[%d]: id=%d\n",ibr,id) ;
00613 #endif
00614
00615 COMPRESS_fclose(far) ;
00616
00617 if( id < DBLK_BRICK_BYTES(blk,ibr) ){
00618 fprintf(stderr,
00619 "\n*** Volume file %s only gave %d out of %d bytes needed!\n",
00620 fnam[ibr] , id , DBLK_BRICK_BYTES(blk,ibr) ) ;
00621 THD_purge_datablock( blk , blk->malloc_type ) ;
00622 for( ii=0 ; ii < nv ; ii++ ) free(fnam[ii]) ;
00623 free(fnam); RETURN( False );
00624 }
00625
00626 }
00627
00628 } else {
00629
00630 int jbr ;
00631
00632
00633
00634 for( jbr=0 ; jbr < nv ; jbr++ ){
00635 ibr = blk->master_ival[jbr] ;
00636 far = COMPRESS_fopen_read( fnam[ibr] ) ;
00637
00638 if( far == NULL ){
00639 fprintf(stderr,
00640 "\n*** Failure while opening volume file %s "
00641 "- do you have permission?\n" ,
00642 fnam[ibr] ) ;
00643 perror("*** Unix error message") ;
00644 THD_purge_datablock( blk , blk->malloc_type ) ;
00645 for( ii=0 ; ii < nv ; ii++ ) free(fnam[ii]) ;
00646 free(fnam); RETURN( False );
00647 }
00648
00649
00650
00651 id = fread( DBLK_ARRAY(blk,jbr), 1,
00652 DBLK_BRICK_BYTES(blk,jbr), far ) ;
00653 if( verb && jbr%PRINT_STEP == 0 ) fprintf(stderr,".") ;
00654
00655 COMPRESS_fclose(far) ;
00656
00657 if( id < DBLK_BRICK_BYTES(blk,jbr) ){
00658 fprintf(stderr,
00659 "\n*** Volume file %s only gave %d out of %d bytes needed!\n",
00660 fnam[ibr] , id , DBLK_BRICK_BYTES(blk,jbr) ) ;
00661 THD_purge_datablock( blk , blk->malloc_type ) ;
00662 for( ii=0 ; ii < nv ; ii++ ) free(fnam[ii]) ;
00663 free(fnam); RETURN( False );
00664 }
00665
00666 }
00667
00668 }
00669
00670
00671
00672
00673
00674 for( ii=0 ; ii < nv ; ii++ ) free(fnam[ii]) ;
00675 free(fnam) ;
00676
00677 }
00678 break ;
00679
00680 }
00681
00682
00683
00684
00685 STATUS("data has been read in") ;
00686
00687
00688
00689 if( dkptr->byte_order != native_order ){
00690 STATUS("byte swapping") ;
00691 if( verb ) fprintf(stderr,".byte swap") ;
00692
00693 for( ibr=0 ; ibr < nv ; ibr++ ){
00694 switch( DBLK_BRICK_TYPE(blk,ibr) ){
00695 case MRI_short:
00696 mri_swap2( DBLK_BRICK_NVOX(blk,ibr) , DBLK_ARRAY(blk,ibr) ) ;
00697 break ;
00698
00699 case MRI_complex:
00700 mri_swap4( 2*DBLK_BRICK_NVOX(blk,ibr), DBLK_ARRAY(blk,ibr)) ;
00701 break ;
00702
00703 case MRI_float:
00704 case MRI_int:
00705 mri_swap4( DBLK_BRICK_NVOX(blk,ibr) , DBLK_ARRAY(blk,ibr) ) ;
00706 break ;
00707 }
00708 }
00709 }
00710
00711
00712
00713
00714 if( floatscan ){
00715 int nerr=0 ;
00716 STATUS("float scanning") ;
00717 if( verb ) fprintf(stderr,".float scan") ;
00718
00719 for( ibr=0 ; ibr < nv ; ibr++ ){
00720 if( DBLK_BRICK_TYPE(blk,ibr) == MRI_float ){
00721 nerr += thd_floatscan( DBLK_BRICK_NVOX(blk,ibr) ,
00722 DBLK_ARRAY(blk,ibr) ) ;
00723
00724 } else if( DBLK_BRICK_TYPE(blk,ibr) == MRI_complex ) {
00725 nerr += thd_complexscan( DBLK_BRICK_NVOX(blk,ibr) ,
00726 DBLK_ARRAY(blk,ibr) ) ;
00727 }
00728 }
00729 if( nerr > 0 )
00730 fprintf(stderr ,
00731 "*** %s: found %d float errors -- see program float_scan\n" ,
00732 dkptr->brick_name , nerr ) ;
00733 }
00734
00735
00736
00737 #if 0
00738 fprintf(stderr,"master_bot=%g master_top=%g\n",blk->master_bot,blk->master_top) ;
00739 #endif
00740 if( DBLK_IS_MASTERED(blk) && blk->master_bot <= blk->master_top ){
00741 float bot = blk->master_bot , top = blk->master_top , fac ;
00742 int jbr ;
00743
00744 STATUS("sub-ranging") ;
00745 if( verb ) fprintf(stderr,".sub-range") ;
00746
00747 for( jbr=0 ; jbr < nv ; jbr++ ){
00748 switch( DBLK_BRICK_TYPE(blk,jbr) ){
00749
00750 default:
00751 fprintf(stderr,
00752 "** Can't sub-range datum type %s!\n",
00753 MRI_TYPE_name[DBLK_BRICK_TYPE(blk,jbr)]) ;
00754 break ;
00755
00756 case MRI_short:{
00757 short mbot, mtop, *mar = (short *) DBLK_ARRAY(blk,jbr) ;
00758 float mfac = DBLK_BRICK_FACTOR(blk,jbr) ;
00759 if( mfac == 0.0 ) mfac = 1.0 ;
00760 mbot = SHORTIZE(bot/mfac) ; mtop = SHORTIZE(top/mfac) ;
00761 #if 0
00762 fprintf(stderr,"mbot=%d mtop=%d\n",(int)mbot,(int)mtop) ;
00763 #endif
00764 for( ii=0 ; ii < nxyz ; ii++ )
00765 if( mar[ii] < mbot || mar[ii] > mtop ) mar[ii] = 0 ;
00766 }
00767 break ;
00768
00769 case MRI_int:{
00770 int mbot, mtop, *mar = (int *) DBLK_ARRAY(blk,jbr) ;
00771 float mfac = DBLK_BRICK_FACTOR(blk,jbr) ;
00772 if( mfac == 0.0 ) mfac = 1.0 ;
00773 mbot = rint(bot/mfac) ; mtop = rint(top/mfac) ;
00774 for( ii=0 ; ii < nxyz ; ii++ )
00775 if( mar[ii] < mbot || mar[ii] > mtop ) mar[ii] = 0 ;
00776 }
00777 break ;
00778
00779 case MRI_byte:{
00780 byte mbot, mtop, *mar = (byte *) DBLK_ARRAY(blk,jbr) ;
00781 float mfac = DBLK_BRICK_FACTOR(blk,jbr) ;
00782 if( mfac == 0.0 ) mfac = 1.0 ;
00783 mbot = BYTEIZE(bot/mfac) ; mtop = BYTEIZE(top/mfac) ;
00784 for( ii=0 ; ii < nxyz ; ii++ )
00785 if( mar[ii] < mbot || mar[ii] > mtop ) mar[ii] = 0 ;
00786 }
00787 break ;
00788
00789 case MRI_float:{
00790 float mbot, mtop, *mar = (float *) DBLK_ARRAY(blk,jbr) ;
00791 float mfac = DBLK_BRICK_FACTOR(blk,jbr) ;
00792 if( mfac == 0.0 ) mfac = 1.0 ;
00793 mbot = (bot/mfac) ; mtop = (top/mfac) ;
00794 for( ii=0 ; ii < nxyz ; ii++ )
00795 if( mar[ii] < mbot || mar[ii] > mtop ) mar[ii] = 0 ;
00796 }
00797 break ;
00798
00799 case MRI_complex:{
00800 float mbot, mtop , val ;
00801 complex *mar = (complex *) DBLK_ARRAY(blk,jbr) ;
00802 float mfac = DBLK_BRICK_FACTOR(blk,jbr) ;
00803 if( mfac == 0.0 ) mfac = 1.0 ;
00804 mbot = (bot/mfac) ; mtop = (top/mfac) ;
00805 mbot = (mbot > 0) ? mbot*mbot : 0 ;
00806 mtop = (mtop > 0) ? mtop*mtop : 0 ;
00807 for( ii=0 ; ii < nxyz ; ii++ ){
00808 val = CSQR(mar[ii]) ;
00809 if( val < mbot || val > mtop ) mar[ii].r = mar[ii].i = 0 ;
00810 }
00811 }
00812 break ;
00813 }
00814 }
00815 }
00816
00817 if( verb ) fprintf(stderr,".done\n") ;
00818
00819 RETURN( True ) ;
00820 }
00821
00822
00823
00824
00825
00826
00827 int THD_alloc_datablock( THD_datablock *blk )
00828 {
00829 int ii,nbad,ibr, nv ;
00830 char *ptr ;
00831
00832 ENTRY("THD_alloc_datablock") ;
00833
00834
00835
00836 if( ! ISVALID_DATABLOCK(blk) || blk->brick == NULL ){
00837 STATUS("Illegal inputs"); RETURN(0);
00838 }
00839
00840 nv = blk->nvals ;
00841 ii = THD_count_databricks( blk ) ;
00842 if( ii == nv ) RETURN(1);
00843
00844 switch( blk->malloc_type ){
00845
00846 default: RETURN(0) ;
00847
00848
00849
00850
00851
00852 case DATABLOCK_MEM_MALLOC:{
00853
00854
00855
00856 STATUS("trying to malloc sub-bricks") ;
00857 for( nbad=ibr=0 ; ibr < nv ; ibr++ ){
00858 if( DBLK_ARRAY(blk,ibr) == NULL ){
00859 ptr = AFMALL(char, DBLK_BRICK_BYTES(blk,ibr) ) ;
00860 mri_fix_data_pointer( ptr , DBLK_BRICK(blk,ibr) ) ;
00861 if( ptr == NULL ) nbad++ ;
00862 }
00863 }
00864 if( nbad == 0 ) RETURN(1) ;
00865
00866
00867
00868 fprintf(stderr,
00869 "\n** failed to malloc %d dataset bricks out of %d - is memory exhausted?\n",
00870 nbad,nv ) ;
00871 #ifdef USING_MCW_MALLOC
00872 { char *str = MCW_MALLOC_status ;
00873 if( str != NULL ) fprintf(stderr,"*** MCW_malloc summary: %s\n",str); }
00874 #endif
00875
00876 if( freeup == NULL ){
00877 for( ibr=0 ; ibr < nv ; ibr++ ){
00878 if( DBLK_ARRAY(blk,ibr) != NULL ){
00879 free(DBLK_ARRAY(blk,ibr)) ;
00880 mri_fix_data_pointer( NULL , DBLK_BRICK(blk,ibr) ) ;
00881 }
00882 }
00883 STATUS("malloc failed, no freeup"); RETURN(0);
00884 }
00885
00886
00887
00888 fprintf(stderr,"** trying to free some memory\n") ;
00889 freeup() ;
00890
00891
00892
00893 for( ibr=0 ; ibr < nv ; ibr++ ){
00894 if( DBLK_ARRAY(blk,ibr) == NULL ){
00895 ptr = AFMALL(char, DBLK_BRICK_BYTES(blk,ibr) ) ;
00896 mri_fix_data_pointer( ptr , DBLK_BRICK(blk,ibr) ) ;
00897 }
00898 }
00899
00900
00901
00902 if( THD_count_databricks(blk) < nv ){
00903 fprintf(stderr,"** cannot free up enough memory\n") ;
00904 #ifdef USING_MCW_MALLOC
00905 { char *str = MCW_MALLOC_status ;
00906 if( str != NULL ) fprintf(stderr,"*** MCW_malloc summary: %s\n",str); }
00907 #endif
00908 for( ibr=0 ; ibr < nv ; ibr++ ){
00909 if( DBLK_ARRAY(blk,ibr) != NULL ){
00910 free(DBLK_ARRAY(blk,ibr)) ;
00911 mri_fix_data_pointer( NULL , DBLK_BRICK(blk,ibr) ) ;
00912 }
00913 }
00914 STATUS("malloc failed; freeup failed"); RETURN(0);
00915 }
00916
00917
00918
00919 STATUS("malloc failed; freeup worked") ;
00920 fprintf(stderr,"*** was able to free up enough memory\n") ;
00921 #ifdef USING_MCW_MALLOC
00922 { char *str = MCW_MALLOC_status ;
00923 if( str != NULL ) fprintf(stderr,"*** MCW_malloc summary: %s\n",str); }
00924 #endif
00925 }
00926 RETURN(1) ;
00927
00928
00929
00930
00931
00932 case DATABLOCK_MEM_SHARED:
00933 #if !defined(DONT_USE_SHM) && !defined(CYGWIN)
00934 { unsigned int offset ;
00935 if( blk->shm_idcode[0] == '\0' ){
00936 UNIQ_idcode_fill( blk->shm_idcode ) ;
00937 blk->shm_idint = shm_create( blk->shm_idcode , (int)blk->total_bytes ) ;
00938 if( blk->shm_idint < 0 ){ blk->shm_idcode[0] = '\0'; RETURN(0); }
00939 ptr = shm_attach( blk->shm_idint ) ;
00940 if( ptr == NULL ){
00941 shmctl( blk->shm_idint , IPC_RMID , NULL ) ;
00942 blk->shm_idint = -1; blk->shm_idcode[0] = '\0'; RETURN(0);
00943 }
00944 offset = 0 ;
00945 for( ibr=0 ; ibr < nv ; ibr++ ){
00946 mri_fix_data_pointer( ptr , DBLK_BRICK(blk,ibr) ) ;
00947 ptr += DBLK_BRICK_BYTES(blk,ibr) ;
00948 }
00949 }
00950 }
00951 RETURN(1) ;
00952 #else
00953 RETURN(0) ;
00954 #endif
00955
00956 }
00957
00958 RETURN(0) ;
00959 }
00960
00961
00962
00963
00964
00965 void THD_zerofill_dataset( THD_3dim_dataset *dset )
00966 {
00967 int ii ;
00968 void *vpt ;
00969
00970 ENTRY("THD_zerofill_dataset") ;
00971
00972 if( !ISVALID_DSET(dset) || !ISVALID_DATABLOCK(dset->dblk) ) EXRETURN ;
00973
00974 for( ii=0 ; ii < DSET_NVALS(dset) ; ii++ ){
00975 if( DSET_ARRAY(dset,ii) == NULL ){
00976 vpt = calloc(1,DSET_BRICK_BYTES(dset,ii)) ;
00977 mri_fix_data_pointer( vpt , DSET_BRICK(dset,ii) ) ;
00978 }
00979 }
00980 EXRETURN ;
00981 }