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 ;                        /* 22 Nov 1999 */
00361       CALC_dshift_mode[ids] = CALC_dshift_mode_current ;
00362    }
00363 
00364    while( nopt < argc && argv[nopt][0] == '-' ){
00365 
00366       /**** -expr expression ****/
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 ) ; /* 15 Sep 1999 */
00384          continue ;
00385       }
00386 
00387       /**** -dsSTOP [22 Nov 1999] ****/
00388 
00389       if( strncmp(argv[nopt],"-dsSTOP",6) == 0 ){
00390          CALC_dshift_mode_current = DSHIFT_MODE_STOP ;
00391          nopt++ ; continue ;
00392       }
00393 
00394       /**** -dsWRAP [22 Nov 1999] ****/
00395 
00396       if( strncmp(argv[nopt],"-dsWRAP",6) == 0 ){
00397          CALC_dshift_mode_current = DSHIFT_MODE_WRAP ;
00398          nopt++ ; continue ;
00399       }
00400 
00401       /**** -dsZERO [22 Nov 1999] ****/
00402 
00403       if( strncmp(argv[nopt],"-dsZERO",6) == 0 ){
00404          CALC_dshift_mode_current = DSHIFT_MODE_ZERO ;
00405          nopt++ ; continue ;
00406       }
00407 
00408       /**** -<letter> dataset ****/
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          /*-- 22 Nov 1999: allow for a differentially
00432                            subscripted name, as in "-b a[1,0,0,0]" --*/
00433 
00434          ll = strlen(argv[nopt]) ;
00435          if( (argv[nopt][0] >= 'a' && argv[nopt][0] <= 'z') &&  /* legal name */
00436              ( (ll >= 3 && argv[nopt][1] == '[') ||             /* subscript */
00437                (ll == 3 &&                                      /*    OR    */
00438                 (argv[nopt][1] == '+' || argv[nopt][1] == '-')) /* +- ijkl */
00439              ) ){
00440 
00441             int jds = argv[nopt][0] - 'a' ;  /* actual dataset index */
00442             int * ijkl ;                     /* array of subscripts */
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             /*- get subscripts -*/
00452 
00453             if( argv[nopt][1] == '[' ){            /* format is [i,j,k,l] */
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 {                               /* format is +i, -j, etc */
00464                 ijkl = (int *) malloc( sizeof(int) * 5 ) ;
00465                 ijkl[1] = ijkl[2] = ijkl[3] = ijkl[4] = 0 ;  /* initialize */
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 ){  /* disallow time subscripting */
00481                fprintf(stderr,
00482                 "++ -cmask: Warning: differential time shifting %s not allowed\n",
00483                        argv[nopt] ) ;
00484                ijkl[4] = 0 ;
00485             }
00486 
00487             /*- more sanity checks -*/
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             /*- set values for later use -*/
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             /*- time to trot, Bwana -*/
00506 
00507             free(ijkl) ; nopt++ ; goto DSET_DONE ;
00508 
00509          } /* end of _dshift */
00510 
00511          /*-- meanwhile, back at the "normal" dataset opening ranch --*/
00512 
00513          { char dname[512] ;                               /* 02 Nov 1999 */
00514 
00515            if( strstr(argv[nopt],"[") == NULL ){
00516               sprintf(dname,"%s[0]",argv[nopt++]) ;        /* add [0] */
00517            } else {
00518               strcpy(dname,argv[nopt++]) ;                 /* don't mangle */
00519            }
00520            dset = THD_open_dataset( dname ) ;              /* open it */
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          /* set some parameters based on the dataset */
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          /* load floating scale factors */
00543 
00544          CALC_ffac[ival] = DSET_BRICK_FACTOR(dset,0) ;
00545          if (CALC_ffac[ival] == 0.0 ) CALC_ffac[ival] = 1.0 ;
00546 
00547          /* read data from disk */
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          /* set pointers for actual dataset arrays */
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           } /* end of switch over type */
00572 
00573 DSET_DONE: continue;
00574 
00575       } /* end of dataset input */
00576 
00577       fprintf(stderr,"** -cmask: Unknown option: %s\n",argv[nopt]) ;
00578       RETURN(1) ;
00579 
00580    }  /* end of loop over options */
00581 
00582    /*---------------------------------------*/
00583    /*** cleanup: check for various errors ***/
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    /* 15 Apr 1999: check if each input dataset is used,
00598                    or if an undefined symbol is used.   */
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 }

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    /*** parse input options ***/
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    /*** make output dataset ***/
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++ )  /* initialize to all zeros */
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       /*** loop over voxels ***/
00139 
00140       for ( ii = 0 ; ii < CALC_nvox ; ii += VSIZE ) {
00141 
00142           jbot = ii ;
00143           jtop = MIN( ii + VSIZE , CALC_nvox ) ;
00144 
00145          /* loop over datasets or other symbol definitions */
00146 
00147           for (ids = 0 ; ids < 26 ; ids ++ ) {
00148 
00149             /* 22 Nov 1999: if a differentially subscripted dataset is here */
00150 
00151             if( CALC_dshift[ids] >= 0 ){
00152                int jds = CALC_dshift[ids] ;     /* actual dataset index */
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 ;                  /* x shift */
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; } /* go to next jj */
00187 
00188                         jy += jd ;                  /* y shift */
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; } /* go to next jj */
00206 
00207                         kz += kd ;                  /* z shift */
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; } /* go to next jj */
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             /* the case of a 3D dataset (i.e., only 1 sub-brick) */
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            /* the case of a voxel (x,y,z) or (i,j,k) coordinate */
00270 
00271            else if( CALC_has_predefined ) {
00272 
00273               switch( ids ){
00274                  case 23:     /* x */
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:     /* y */
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:     /* z */
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:     /* i */
00296                     if( HAS_I )
00297                      for( jj=jbot ; jj < jtop ; jj++ )
00298                        atoz[ids][jj-ii] = (jj%nx) ;
00299                  break ;
00300 
00301                  case 9:     /* j */
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:    /* k */
00308                     if( HAS_K )
00309                      for( jj=jbot ; jj < jtop ; jj++ )
00310                        atoz[ids][jj-ii] = (jj/nxy) ;
00311                  break ;
00312 
00313                } /* end of switch on symbol subscript */
00314 
00315               } /* end of choice over data type (if-else cascade) */
00316              } /* end of loop over datasets/symbols */
00317 
00318             /**** actually do the work! ****/
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          } /* end of loop over space (voxels) */
00325 
00326    /* cleanup and go home */
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 }