Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
edt_calcmask.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "mrilib.h"
00012 #include "parser.h"
00013
00014
00015
00016 static int CALC_nvox ;
00017 static PARSER_code * CALC_code ;
00018
00019
00020
00021 #define DSHIFT_MODE_STOP 0
00022 #define DSHIFT_MODE_WRAP 1
00023 #define DSHIFT_MODE_ZERO 2
00024
00025 static int CALC_dshift [26] ;
00026 static int CALC_dshift_i [26] ;
00027 static int CALC_dshift_j [26] ;
00028 static int CALC_dshift_k [26] ;
00029 static int CALC_dshift_l [26] ;
00030 static int CALC_dshift_mode[26] ;
00031
00032 static int CALC_dshift_mode_current ;
00033
00034
00035
00036 static int CALC_has_sym[26] ;
00037 static char abet[] = "abcdefghijklmnopqrstuvwxyz" ;
00038
00039 #define HAS_I CALC_has_sym[ 8]
00040 #define HAS_J CALC_has_sym[ 9]
00041 #define HAS_K CALC_has_sym[10]
00042 #define HAS_X CALC_has_sym[23]
00043 #define HAS_Y CALC_has_sym[24]
00044 #define HAS_Z CALC_has_sym[25]
00045
00046 #define PREDEFINED_MASK ((1<< 8)|(1<< 9)|(1<<10)|(1<<23)|(1<<24)|(1<<25))
00047
00048 static int CALC_has_predefined ;
00049
00050 static THD_3dim_dataset * CALC_dset[26] ;
00051 static int CALC_type[26] ;
00052 static byte * CALC_byte[26] ;
00053 static short * CALC_short[26] ;
00054 static float * CALC_float[26] ;
00055 static float CALC_ffac[26] ;
00056
00057
00058
00059 #define VAR_DEFINED(kv) (CALC_dset[kv] != NULL || CALC_dshift[kv] >= 0)
00060
00061
00062
00063 static int CALC_read_opts( int argc , char * argv[] ) ;
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083 byte * EDT_calcmask( char * cmd , int * nxyz )
00084 {
00085 int Argc=0 ;
00086 char ** Argv=NULL ;
00087 byte * bmask ;
00088
00089 #define VSIZE 1024
00090
00091 double * atoz[26] ;
00092 int ii , ids , jj, ll, jbot, jtop ;
00093 THD_3dim_dataset * new_dset ;
00094 double temp[VSIZE];
00095
00096 int nx,nxy ;
00097 THD_dataxes * daxes ;
00098
00099 ENTRY("EDT_calcmask") ;
00100
00101
00102
00103 if( cmd == NULL ) RETURN( NULL );
00104 append_string_to_args( cmd , 0,NULL , &Argc , &Argv ) ;
00105 if( Argc == 0 || Argv == NULL ) RETURN( NULL );
00106
00107 jj = CALC_read_opts( Argc , Argv ) ;
00108
00109 for( ii=0 ; ii < Argc ; ii++ ) free(Argv[ii]) ;
00110 free(Argv) ;
00111
00112 if( jj != 0 ){
00113 if( CALC_code != NULL ) free(CALC_code) ;
00114 for( ids=0 ; ids < 26 ; ids++ ){
00115 if( CALC_dset[ids] != NULL ) DSET_delete( CALC_dset[ids] ) ;
00116 }
00117 RETURN( NULL );
00118 }
00119
00120
00121
00122 for( ids=0 ; ids < 26 ; ids++ ) if( CALC_dset[ids] != NULL ) break ;
00123
00124 new_dset = EDIT_empty_copy( CALC_dset[ids] ) ;
00125
00126 for (ids=0; ids<26; ids++)
00127 atoz[ids] = (double *) malloc(sizeof(double) * VSIZE ) ;
00128
00129 for( ids=0 ; ids < 26 ; ids++ )
00130 for (ii=0; ii<VSIZE; ii++)
00131 atoz[ids][ii] = 0.0 ;
00132
00133 nx = DSET_NX(new_dset) ;
00134 nxy = nx * DSET_NY(new_dset) ; daxes = new_dset->daxes ;
00135
00136 bmask = (byte *) malloc(sizeof(byte) * CALC_nvox) ;
00137
00138
00139
00140 for ( ii = 0 ; ii < CALC_nvox ; ii += VSIZE ) {
00141
00142 jbot = ii ;
00143 jtop = MIN( ii + VSIZE , CALC_nvox ) ;
00144
00145
00146
00147 for (ids = 0 ; ids < 26 ; ids ++ ) {
00148
00149
00150
00151 if( CALC_dshift[ids] >= 0 ){
00152 int jds = CALC_dshift[ids] ;
00153 int jjs , ix,jy,kz ;
00154 int id=CALC_dshift_i[ids] , jd=CALC_dshift_j[ids] ,
00155 kd=CALC_dshift_k[ids] , ld=CALC_dshift_l[ids] ;
00156 int ijkd = ((id!=0) || (jd!=0) || (kd!=0)) ;
00157 int dsx = DSET_NX(CALC_dset[jds]) - 1 ;
00158 int dsy = DSET_NY(CALC_dset[jds]) - 1 ;
00159 int dsz = DSET_NZ(CALC_dset[jds]) - 1 ;
00160 int mode = CALC_dshift_mode[ids] , dun ;
00161
00162 for( dun=0,jj=jbot ; jj < jtop ; jj++ ){
00163 jjs = jj ;
00164 if( ijkd ){
00165 ix = DSET_index_to_ix(CALC_dset[jds],jj) ;
00166 jy = DSET_index_to_jy(CALC_dset[jds],jj) ;
00167 kz = DSET_index_to_kz(CALC_dset[jds],jj) ;
00168
00169 ix += id ;
00170 if( ix < 0 || ix > dsx ){
00171 switch( mode ){
00172 case DSHIFT_MODE_ZERO:
00173 atoz[ids][jj-ii] = 0.0 ; dun = 1 ;
00174 break ;
00175 default:
00176 case DSHIFT_MODE_STOP:
00177 if( ix < 0 ) ix = 0 ;
00178 else if( ix > dsx ) ix = dsx ;
00179 break ;
00180 case DSHIFT_MODE_WRAP:
00181 while( ix < 0 ) ix += (dsx+1) ;
00182 while( ix > dsx ) ix -= (dsx+1) ;
00183 break ;
00184 }
00185 }
00186 if( dun ){ dun=0; continue; }
00187
00188 jy += jd ;
00189 if( jy < 0 || jy > dsy ){
00190 switch( mode ){
00191 case DSHIFT_MODE_ZERO:
00192 atoz[ids][jj-ii] = 0.0 ; dun = 1 ;
00193 break ;
00194 default:
00195 case DSHIFT_MODE_STOP:
00196 if( jy < 0 ) jy = 0 ;
00197 else if( jy > dsy ) jy = dsy ;
00198 break ;
00199 case DSHIFT_MODE_WRAP:
00200 while( jy < 0 ) jy += (dsy+1) ;
00201 while( jy > dsy ) jy -= (dsy+1) ;
00202 break ;
00203 }
00204 }
00205 if( dun ){ dun=0; continue; }
00206
00207 kz += kd ;
00208 if( kz < 0 || kz > dsz ){
00209 switch( mode ){
00210 case DSHIFT_MODE_ZERO:
00211 atoz[ids][jj-ii] = 0.0 ; dun = 1 ;
00212 break ;
00213 default:
00214 case DSHIFT_MODE_STOP:
00215 if( kz < 0 ) kz = 0 ;
00216 else if( kz > dsz ) kz = dsz ;
00217 break ;
00218 case DSHIFT_MODE_WRAP:
00219 while( kz < 0 ) kz += (dsz+1) ;
00220 while( kz > dsz ) kz -= (dsz+1) ;
00221 break ;
00222 }
00223 }
00224 if( dun ){ dun=0; continue; }
00225
00226 jjs = DSET_ixyz_to_index(CALC_dset[jds],ix,jy,kz) ;
00227 }
00228 switch( CALC_type[jds] ) {
00229 case MRI_short:
00230 atoz[ids][jj-ii] = CALC_short[jds][jjs]
00231 * CALC_ffac[jds];
00232 break ;
00233 case MRI_float:
00234 atoz[ids][jj-ii] = CALC_float[jds][jjs]
00235 * CALC_ffac[jds];
00236 break ;
00237 case MRI_byte:
00238 atoz[ids][jj-ii] = CALC_byte[jds][jjs]
00239 * CALC_ffac[jds];
00240 break ;
00241 }
00242 }
00243 }
00244
00245
00246
00247 else if ( CALC_type[ids] >= 0 ) {
00248 switch( CALC_type[ids] ) {
00249 case MRI_short:
00250 for (jj =jbot ; jj < jtop ; jj ++ ){
00251 atoz[ids][jj-ii] = CALC_short[ids][jj] * CALC_ffac[ids] ;
00252 }
00253 break;
00254
00255 case MRI_float:
00256 for (jj =jbot ; jj < jtop ; jj ++ ){
00257 atoz[ids][jj-ii] = CALC_float[ids][jj] * CALC_ffac[ids] ;
00258 }
00259 break;
00260
00261 case MRI_byte:
00262 for (jj =jbot ; jj < jtop ; jj ++ ){
00263 atoz[ids][jj-ii] = CALC_byte[ids][jj] * CALC_ffac[ids] ;
00264 }
00265 break;
00266 }
00267 }
00268
00269
00270
00271 else if( CALC_has_predefined ) {
00272
00273 switch( ids ){
00274 case 23:
00275 if( HAS_X )
00276 for( jj=jbot ; jj < jtop ; jj++ )
00277 atoz[ids][jj-ii] = daxes->xxorg +
00278 (jj%nx) * daxes->xxdel ;
00279 break ;
00280
00281 case 24:
00282 if( HAS_Y )
00283 for( jj=jbot ; jj < jtop ; jj++ )
00284 atoz[ids][jj-ii] = daxes->yyorg +
00285 ((jj%nxy)/nx) * daxes->yydel ;
00286 break ;
00287
00288 case 25:
00289 if( HAS_Z )
00290 for( jj=jbot ; jj < jtop ; jj++ )
00291 atoz[ids][jj-ii] = daxes->zzorg +
00292 (jj/nxy) * daxes->zzdel ;
00293 break ;
00294
00295 case 8:
00296 if( HAS_I )
00297 for( jj=jbot ; jj < jtop ; jj++ )
00298 atoz[ids][jj-ii] = (jj%nx) ;
00299 break ;
00300
00301 case 9:
00302 if( HAS_J )
00303 for( jj=jbot ; jj < jtop ; jj++ )
00304 atoz[ids][jj-ii] = ((jj%nxy)/nx) ;
00305 break ;
00306
00307 case 10:
00308 if( HAS_K )
00309 for( jj=jbot ; jj < jtop ; jj++ )
00310 atoz[ids][jj-ii] = (jj/nxy) ;
00311 break ;
00312
00313 }
00314
00315 }
00316 }
00317
00318
00319
00320 PARSER_evaluate_vector(CALC_code, atoz, jtop-jbot, temp);
00321 for ( jj = jbot ; jj < jtop ; jj ++ )
00322 bmask[jj] = (temp[jj-ii] != 0.0) ;
00323
00324 }
00325
00326
00327
00328 for( ids=0 ; ids < 26 ; ids++ ){
00329 free(atoz[ids]) ;
00330 if( CALC_dset[ids] != NULL ) DSET_delete( CALC_dset[ids] ) ;
00331 }
00332 DSET_delete(new_dset) ;
00333 free(CALC_code) ;
00334
00335 if( nxyz != NULL ) *nxyz = CALC_nvox ;
00336 RETURN( bmask );
00337 }
00338
00339
00340
00341
00342
00343 static int CALC_read_opts( int argc , char * argv[] )
00344 {
00345 int nopt = 0 ;
00346 int ids ;
00347 int ii ;
00348
00349 ENTRY("CALC_read_opts") ;
00350
00351 CALC_nvox = -1 ;
00352 CALC_code = NULL ;
00353 CALC_dshift_mode_current = DSHIFT_MODE_STOP ;
00354 CALC_has_predefined = 0 ;
00355
00356 for( ids=0 ; ids < 26 ; ids++ ){
00357 CALC_dset[ids] = NULL ;
00358 CALC_type[ids] = -1 ;
00359
00360 CALC_dshift[ids] = -1 ;
00361 CALC_dshift_mode[ids] = CALC_dshift_mode_current ;
00362 }
00363
00364 while( nopt < argc && argv[nopt][0] == '-' ){
00365
00366
00367
00368 if( strncmp(argv[nopt],"-expr",4) == 0 ){
00369 if( CALC_code != NULL ){
00370 fprintf(stderr,
00371 "** -cmask: cannot have 2 -expr options!\n") ; RETURN(1) ;
00372 }
00373 nopt++ ;
00374 if( nopt >= argc ){
00375 fprintf(stderr,
00376 "** -cmask: need argument after -expr!\n") ; RETURN(1) ;
00377 }
00378 CALC_code = PARSER_generate_code( argv[nopt++] ) ;
00379 if( CALC_code == NULL ){
00380 fprintf(stderr,
00381 "** -cmask: illegal expression!\n") ; RETURN(1) ;
00382 }
00383 PARSER_mark_symbols( CALC_code , CALC_has_sym ) ;
00384 continue ;
00385 }
00386
00387
00388
00389 if( strncmp(argv[nopt],"-dsSTOP",6) == 0 ){
00390 CALC_dshift_mode_current = DSHIFT_MODE_STOP ;
00391 nopt++ ; continue ;
00392 }
00393
00394
00395
00396 if( strncmp(argv[nopt],"-dsWRAP",6) == 0 ){
00397 CALC_dshift_mode_current = DSHIFT_MODE_WRAP ;
00398 nopt++ ; continue ;
00399 }
00400
00401
00402
00403 if( strncmp(argv[nopt],"-dsZERO",6) == 0 ){
00404 CALC_dshift_mode_current = DSHIFT_MODE_ZERO ;
00405 nopt++ ; continue ;
00406 }
00407
00408
00409
00410 ids = strlen( argv[nopt] ) ;
00411
00412 if( (argv[nopt][1] >= 'a' && argv[nopt][1] <= 'z') && ids == 2 ) {
00413
00414 int ival , nxyz , ll ;
00415 THD_3dim_dataset * dset ;
00416
00417 ival = argv[nopt][1] - 'a' ;
00418 if( VAR_DEFINED(ival) ){
00419 fprintf(stderr,
00420 "** -cmask: Can't define %c symbol twice\n",argv[nopt][1]);
00421 RETURN(1) ;
00422 }
00423
00424 nopt++ ;
00425 if( nopt >= argc ){
00426 fprintf(stderr,
00427 "** -cmask: need argument after %s\n",argv[nopt-1]);
00428 RETURN(1) ;
00429 }
00430
00431
00432
00433
00434 ll = strlen(argv[nopt]) ;
00435 if( (argv[nopt][0] >= 'a' && argv[nopt][0] <= 'z') &&
00436 ( (ll >= 3 && argv[nopt][1] == '[') ||
00437 (ll == 3 &&
00438 (argv[nopt][1] == '+' || argv[nopt][1] == '-'))
00439 ) ){
00440
00441 int jds = argv[nopt][0] - 'a' ;
00442 int * ijkl ;
00443
00444 if( CALC_dset[jds] == NULL ){
00445 fprintf(stderr,
00446 "** -cmask: Must define dataset %c before using it in %s\n",
00447 argv[nopt][0] , argv[nopt] ) ;
00448 RETURN(1) ;
00449 }
00450
00451
00452
00453 if( argv[nopt][1] == '[' ){
00454 MCW_intlist_allow_negative(1) ;
00455 ijkl = MCW_get_intlist( 9999 , argv[nopt]+1 ) ;
00456 MCW_intlist_allow_negative(0) ;
00457 if( ijkl == NULL || ijkl[0] != 4 ){
00458 fprintf(stderr,
00459 "** -cmask: Illegal differential subscripting %s\n",
00460 argv[nopt] ) ;
00461 RETURN(1) ;
00462 }
00463 } else {
00464 ijkl = (int *) malloc( sizeof(int) * 5 ) ;
00465 ijkl[1] = ijkl[2] = ijkl[3] = ijkl[4] = 0 ;
00466 switch( argv[nopt][2] ){
00467 default:
00468 fprintf(stderr,
00469 "** -cmask: Bad differential subscripting %s\n",
00470 argv[nopt] ) ;
00471 RETURN(1) ;
00472
00473 case 'i': ijkl[1] = (argv[nopt][1]=='+') ? 1 : -1 ; break ;
00474 case 'j': ijkl[2] = (argv[nopt][1]=='+') ? 1 : -1 ; break ;
00475 case 'k': ijkl[3] = (argv[nopt][1]=='+') ? 1 : -1 ; break ;
00476 case 'l': ijkl[4] = (argv[nopt][1]=='+') ? 1 : -1 ; break ;
00477 }
00478 }
00479
00480 if( ijkl[4] != 0 ){
00481 fprintf(stderr,
00482 "++ -cmask: Warning: differential time shifting %s not allowed\n",
00483 argv[nopt] ) ;
00484 ijkl[4] = 0 ;
00485 }
00486
00487
00488
00489 if( ijkl[1]==0 && ijkl[2]==0 && ijkl[3]==0 && ijkl[4]==0 ){
00490 fprintf(stderr,
00491 "++ -cmask: differential subscript %s is all zero\n",
00492 argv[nopt] ) ;
00493 }
00494
00495
00496
00497 CALC_dshift [ival] = jds ;
00498 CALC_dshift_i[ival] = ijkl[1] ;
00499 CALC_dshift_j[ival] = ijkl[2] ;
00500 CALC_dshift_k[ival] = ijkl[3] ;
00501 CALC_dshift_l[ival] = ijkl[4] ;
00502
00503 CALC_dshift_mode[ival] = CALC_dshift_mode_current ;
00504
00505
00506
00507 free(ijkl) ; nopt++ ; goto DSET_DONE ;
00508
00509 }
00510
00511
00512
00513 { char dname[512] ;
00514
00515 if( strstr(argv[nopt],"[") == NULL ){
00516 sprintf(dname,"%s[0]",argv[nopt++]) ;
00517 } else {
00518 strcpy(dname,argv[nopt++]) ;
00519 }
00520 dset = THD_open_dataset( dname ) ;
00521 if( dset == NULL ){
00522 fprintf(stderr,
00523 "** -cmask: can't open dataset %s\n",dname) ;
00524 RETURN(1) ;
00525 }
00526 }
00527 CALC_dset[ival] = dset ;
00528
00529
00530
00531 nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz ;
00532 if( CALC_nvox < 0 ){
00533 CALC_nvox = nxyz ;
00534 } else if( nxyz != CALC_nvox ){
00535 fprintf(stderr,
00536 "** -cmask: dataset %s differs in size from others\n",argv[nopt-1]);
00537 RETURN(1) ;
00538 }
00539
00540 CALC_type[ival] = DSET_BRICK_TYPE(dset,0) ;
00541
00542
00543
00544 CALC_ffac[ival] = DSET_BRICK_FACTOR(dset,0) ;
00545 if (CALC_ffac[ival] == 0.0 ) CALC_ffac[ival] = 1.0 ;
00546
00547
00548
00549 THD_load_datablock( dset->dblk ) ;
00550 if( ! DSET_LOADED(dset) ){
00551 fprintf(stderr,
00552 "** -cmask: Can't read data brick for dataset %s\n",argv[nopt-1]) ;
00553 RETURN(1) ;
00554 }
00555
00556
00557
00558 switch (CALC_type[ival]) {
00559 case MRI_short:
00560 CALC_short[ival] = (short *) DSET_ARRAY(dset,0) ;
00561 break;
00562
00563 case MRI_float:
00564 CALC_float[ival] = (float *) DSET_ARRAY(dset,0) ;
00565 break;
00566
00567 case MRI_byte:
00568 CALC_byte[ival] = (byte *) DSET_ARRAY(dset,0) ;
00569 break;
00570
00571 }
00572
00573 DSET_DONE: continue;
00574
00575 }
00576
00577 fprintf(stderr,"** -cmask: Unknown option: %s\n",argv[nopt]) ;
00578 RETURN(1) ;
00579
00580 }
00581
00582
00583
00584
00585 for( ids=0 ; ids < 26 ; ids++ ) if( CALC_dset[ids] != NULL ) break ;
00586 if( ids == 26 ){
00587 fprintf(stderr,
00588 "** -cmask: No actual input datasets given!\n") ;
00589 RETURN(1) ;
00590 }
00591
00592 if( CALC_code == NULL ){
00593 fprintf(stderr,"** -cmask: No expression given!\n") ;
00594 RETURN(1) ;
00595 }
00596
00597
00598
00599
00600 for (ids=0; ids < 26; ids ++){
00601 if( VAR_DEFINED(ids) && !CALC_has_sym[ids] )
00602 fprintf(stderr ,
00603 "++ -cmask: input '%c' is not used in the expression\n" ,
00604 abet[ids] ) ;
00605
00606 else if( !VAR_DEFINED(ids) && CALC_has_sym[ids] ){
00607
00608 if( ((1<<ids) & PREDEFINED_MASK) == 0 )
00609 fprintf(stderr ,
00610 "++ -cmask: symbol %c is used but not defined\n" ,
00611 abet[ids] ) ;
00612 else {
00613 CALC_has_predefined++ ;
00614 }
00615 }
00616 }
00617
00618 RETURN(0) ;
00619 }