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 /*******************************************************************
8    Adapted from 3dcalc.c - RWCox - 16 Mar 2000
9 ********************************************************************/
10 
11 #include "mrilib.h"
12 #include "parser.h"
13 
14 /*-------------------------- global data --------------------------*/
15 
16 static int                CALC_nvox ;
17 static PARSER_code *      CALC_code ;
18 
19 /*---------- dshift stuff [22 Nov 1999] ----------*/
20 
21 #define DSHIFT_MODE_STOP  0
22 #define DSHIFT_MODE_WRAP  1
23 #define DSHIFT_MODE_ZERO  2
24 
25 static int                CALC_dshift     [26] ; /* 22 Nov 1999 */
26 static int                CALC_dshift_i   [26] ;
27 static int                CALC_dshift_j   [26] ;
28 static int                CALC_dshift_k   [26] ;
29 static int                CALC_dshift_l   [26] ;
30 static int                CALC_dshift_mode[26] ;
31 
32 static int                CALC_dshift_mode_current ;
33 
34 /*------------------------------------------------*/
35 
36 static int   CALC_has_sym[26] ;                      /* 15 Sep 1999 */
37 static char  abet[] = "abcdefghijklmnopqrstuvwxyz" ;
38 
39 #define HAS_I  CALC_has_sym[ 8]
40 #define HAS_J  CALC_has_sym[ 9]
41 #define HAS_K  CALC_has_sym[10]
42 #define HAS_X  CALC_has_sym[23]
43 #define HAS_Y  CALC_has_sym[24]
44 #define HAS_Z  CALC_has_sym[25]
45 
46 #define PREDEFINED_MASK ((1<< 8)|(1<< 9)|(1<<10)|(1<<23)|(1<<24)|(1<<25))
47 
48 static int     CALC_has_predefined ;  /* 19 Nov 1999 */
49 
50 static THD_3dim_dataset *  CALC_dset[26] ;
51 static int                 CALC_type[26] ;
52 static byte *              CALC_byte[26] ;
53 static short *             CALC_short[26] ;
54 static float *             CALC_float[26] ;
55 static float               CALC_ffac[26] ;
56 
57 /* this macro tells if a variable (index 0..25) is defined  */
58 
59 #define VAR_DEFINED(kv) (CALC_dset[kv] != NULL || CALC_dshift[kv] >= 0)
60 
61 /* prototype */
62 
63 static int CALC_read_opts( int argc , char * argv[] ) ;
64 
65 /*------------------------------------------------------------------
66   Input: cmd  = a command string, like the options for 3dcalc
67          nxyz = pointer to integer
68          forced_mask_length = force length of mask to be
69                               forced_mask_length. This is used
70                               for sparse surface dsets.
71                               Use 0 when not dealing with surface
72                               related dsets.
73                            This parameter is not used if
74                            CALC_dset[ids]->blk->node_list = NULL
75   Output: return value is a byte mask (array of 0 or 1)
76           *nxyz = number of voxels in output array
77   The returned array should be free()-ed when its usefulness ends.
78 
79   Example:
80     byte * bm ; int ibm ;
81     bm = EDT_calcmask( "-a fred+orig[7] -b ethel+orig[0]"
82                        "-expr AND(step(a-99),b)" , &ibm  ,
83                        0 ) ;
84 
85     Here, bm[i] is 1 if the 7th sub-brick of fred+orig is
86     greater than 99 at the i-th voxel, and at the same time
87     the 0th sub-brick of ethel+orig is nonzero at the i-th voxel.
88 --------------------------------------------------------------------*/
89 
EDT_calcmask(char * cmd,int * nxyz,int forced_mask_length)90 byte * EDT_calcmask( char * cmd , int * nxyz, int forced_mask_length )
91 {
92    int Argc=0 ;
93    char ** Argv=NULL ;
94    byte * bmask=NULL ;
95 
96 #define VSIZE 1024
97 
98    double * atoz[26] ;
99    int zii, ii , ids , jj, ll, jbot, jtop, nodemax=-1 ;
100    THD_3dim_dataset * new_dset, *node_dset=NULL;
101    double   temp[VSIZE];
102 
103    int   nx,nxy ;
104    THD_dataxes * daxes ;
105 
106 ENTRY("EDT_calcmask") ;
107 
108    /*** parse input options ***/
109 
110    if( cmd == NULL ) RETURN( NULL );
111    append_string_to_args( cmd , 0,NULL , &Argc , &Argv ) ;
112    if( Argc == 0 || Argv == NULL ) RETURN( NULL );
113 
114    jj = CALC_read_opts( Argc , Argv ) ;
115 
116    for( ii=0 ; ii < Argc ; ii++ ) free(Argv[ii]) ;
117    free(Argv) ;
118 
119    if( jj != 0 ){
120       if( CALC_code != NULL ) free(CALC_code) ;
121       for( ids=0 ; ids < 26 ; ids++ ){
122          if( CALC_dset[ids] != NULL ) DSET_delete( CALC_dset[ids] ) ;
123       }
124       RETURN( NULL );
125    }
126 
127    /*** make output dataset ***/
128 
129    for( ids=0 ; ids < 26 ; ids++ ) if( CALC_dset[ids] != NULL ) break ;
130 
131    new_dset = EDIT_empty_copy( CALC_dset[ids] ) ;
132    if (nodemax < 0 && CALC_dset[ids]->dblk->node_list) {
133       /*keep track of the maximum node index */
134       node_dset = CALC_dset[ids];
135       nodemax = CALC_dset[ids]->dblk->node_list[0];
136       for (zii=0; zii<CALC_dset[ids]->dblk->nnodes; ++zii)
137          if (nodemax < CALC_dset[ids]->dblk->node_list[zii])  nodemax = CALC_dset[ids]->dblk->node_list[zii];
138       if (forced_mask_length > 0) {
139          if (nodemax+1 > forced_mask_length) {
140             fprintf(stderr,"** nodemax(+1) of %d(+1) is larger than the forced_mask_length of %d\n"
141                         "   -cmask datasets may be inappropriate for surface used\n",
142                         nodemax, forced_mask_length);
143             CALC_nvox = 0;
144             goto CLEANUP;
145          }
146          nodemax = forced_mask_length-1;
147       }
148    }
149 
150    for (ids=0; ids<26; ids++)
151       atoz[ids] = (double *) malloc(sizeof(double) * VSIZE ) ;
152 
153    for( ids=0 ; ids < 26 ; ids++ )  /* initialize to all zeros */
154       for (ii=0; ii<VSIZE; ii++)
155           atoz[ids][ii] = 0.0 ;
156 
157    nx  =      DSET_NX(new_dset) ;
158    nxy = nx * DSET_NY(new_dset) ; daxes = new_dset->daxes ;
159 
160    bmask = (byte *) malloc(sizeof(byte) * CALC_nvox) ;
161 
162       /*** loop over voxels ***/
163 
164       for ( ii = 0 ; ii < CALC_nvox ; ii += VSIZE ) {
165 
166           jbot = ii ;
167           jtop = MIN( ii + VSIZE , CALC_nvox ) ;
168 
169          /* loop over datasets or other symbol definitions */
170 
171           for (ids = 0 ; ids < 26 ; ids ++ ) {
172 
173             /* 22 Nov 1999: if a differentially subscripted dataset is here */
174 
175             if( CALC_dshift[ids] >= 0 ){
176                int jds = CALC_dshift[ids] ;     /* actual dataset index */
177                int jjs , ix,jy,kz ;
178                int id=CALC_dshift_i[ids] , jd=CALC_dshift_j[ids] ,
179                    kd=CALC_dshift_k[ids] , ld=CALC_dshift_l[ids] ;
180                int ijkd = ((id!=0) || (jd!=0) || (kd!=0)) ;
181                int dsx = DSET_NX(CALC_dset[jds]) - 1 ;
182                int dsy = DSET_NY(CALC_dset[jds]) - 1 ;
183                int dsz = DSET_NZ(CALC_dset[jds]) - 1 ;
184                int mode = CALC_dshift_mode[ids] , dun ;
185 
186                   for( dun=0,jj=jbot ; jj < jtop ; jj++ ){
187                      jjs = jj ;
188                      if( ijkd ){
189                         ix = DSET_index_to_ix(CALC_dset[jds],jj) ;
190                         jy = DSET_index_to_jy(CALC_dset[jds],jj) ;
191                         kz = DSET_index_to_kz(CALC_dset[jds],jj) ;
192 
193                         ix += id ;                  /* x shift */
194                         if( ix < 0 || ix > dsx ){
195                            switch( mode ){
196                               case DSHIFT_MODE_ZERO:
197                                  atoz[ids][jj-ii] = 0.0 ; dun = 1 ;
198                               break ;
199                               default:
200                               case DSHIFT_MODE_STOP:
201                                      if( ix <  0  ) ix = 0   ;
202                                 else if( ix > dsx ) ix = dsx ;
203                               break ;
204                               case DSHIFT_MODE_WRAP:
205                                  while( ix <  0  ) ix += (dsx+1) ;
206                                  while( ix > dsx ) ix -= (dsx+1) ;
207                               break ;
208                            }
209                         }
210                         if( dun ){ dun=0; continue; } /* go to next jj */
211 
212                         jy += jd ;                  /* y shift */
213                         if( jy < 0 || jy > dsy ){
214                            switch( mode ){
215                               case DSHIFT_MODE_ZERO:
216                                  atoz[ids][jj-ii] = 0.0 ; dun = 1 ;
217                               break ;
218                               default:
219                               case DSHIFT_MODE_STOP:
220                                      if( jy <  0  ) jy = 0   ;
221                                 else if( jy > dsy ) jy = dsy ;
222                               break ;
223                               case DSHIFT_MODE_WRAP:
224                                  while( jy <  0  ) jy += (dsy+1) ;
225                                  while( jy > dsy ) jy -= (dsy+1) ;
226                               break ;
227                            }
228                         }
229                         if( dun ){ dun=0; continue; } /* go to next jj */
230 
231                         kz += kd ;                  /* z shift */
232                         if( kz < 0 || kz > dsz ){
233                            switch( mode ){
234                               case DSHIFT_MODE_ZERO:
235                                  atoz[ids][jj-ii] = 0.0 ; dun = 1 ;
236                               break ;
237                               default:
238                               case DSHIFT_MODE_STOP:
239                                      if( kz <  0  ) kz = 0   ;
240                                 else if( kz > dsz ) kz = dsz ;
241                               break ;
242                               case DSHIFT_MODE_WRAP:
243                                  while( kz <  0  ) kz += (dsz+1) ;
244                                  while( kz > dsz ) kz -= (dsz+1) ;
245                               break ;
246                            }
247                         }
248                         if( dun ){ dun=0; continue; } /* go to next jj */
249 
250                         jjs = DSET_ixyz_to_index(CALC_dset[jds],ix,jy,kz) ;
251                      }
252                      switch( CALC_type[jds] ) {
253                         case MRI_short:
254                            atoz[ids][jj-ii] =  CALC_short[jds][jjs]
255                                              * CALC_ffac[jds];
256                         break ;
257                         case MRI_float:
258                            atoz[ids][jj-ii] =  CALC_float[jds][jjs]
259                                              * CALC_ffac[jds];
260                         break ;
261                         case MRI_byte:
262                            atoz[ids][jj-ii] =  CALC_byte[jds][jjs]
263                                              * CALC_ffac[jds];
264                         break ;
265                      }
266                   }
267             }
268 
269             /* the case of a 3D dataset (i.e., only 1 sub-brick) */
270 
271             else if ( CALC_type[ids] >= 0 ) {
272                switch( CALC_type[ids] ) {
273                     case MRI_short:
274                        for (jj =jbot ; jj < jtop ; jj ++ ){
275                            atoz[ids][jj-ii] = CALC_short[ids][jj] * CALC_ffac[ids] ;
276                      }
277                     break;
278 
279                   case MRI_float:
280                      for (jj =jbot ; jj < jtop ; jj ++ ){
281                         atoz[ids][jj-ii] = CALC_float[ids][jj] * CALC_ffac[ids] ;
282                      }
283                   break;
284 
285                   case MRI_byte:
286                      for (jj =jbot ; jj < jtop ; jj ++ ){
287                         atoz[ids][jj-ii] = CALC_byte[ids][jj] * CALC_ffac[ids] ;
288                      }
289                   break;
290                }
291             }
292 
293            /* the case of a voxel (x,y,z) or (i,j,k) coordinate */
294 
295            else if( CALC_has_predefined ) {
296 
297               switch( ids ){
298                  case 23:     /* x */
299                     if( HAS_X )
300                      for( jj=jbot ; jj < jtop ; jj++ )
301                        atoz[ids][jj-ii] = daxes->xxorg +
302                                           (jj%nx) * daxes->xxdel ;
303                  break ;
304 
305                  case 24:     /* y */
306                     if( HAS_Y )
307                      for( jj=jbot ; jj < jtop ; jj++ )
308                        atoz[ids][jj-ii] = daxes->yyorg +
309                                           ((jj%nxy)/nx) * daxes->yydel ;
310                  break ;
311 
312                  case 25:     /* z */
313                     if( HAS_Z )
314                      for( jj=jbot ; jj < jtop ; jj++ )
315                        atoz[ids][jj-ii] = daxes->zzorg +
316                                           (jj/nxy) * daxes->zzdel ;
317                  break ;
318 
319                  case 8:     /* i */
320                     if( HAS_I )
321                      for( jj=jbot ; jj < jtop ; jj++ )
322                        atoz[ids][jj-ii] = (jj%nx) ;
323                  break ;
324 
325                  case 9:     /* j */
326                     if( HAS_J )
327                      for( jj=jbot ; jj < jtop ; jj++ )
328                        atoz[ids][jj-ii] = ((jj%nxy)/nx) ;
329                  break ;
330 
331                  case 10:    /* k */
332                     if( HAS_K )
333                      for( jj=jbot ; jj < jtop ; jj++ )
334                        atoz[ids][jj-ii] = (jj/nxy) ;
335                  break ;
336 
337                } /* end of switch on symbol subscript */
338 
339               } /* end of choice over data type (if-else cascade) */
340              } /* end of loop over datasets/symbols */
341 
342             /**** actually do the work! ****/
343 
344             PARSER_evaluate_vector(CALC_code, atoz, jtop-jbot, temp);
345              for ( jj = jbot ; jj < jtop ; jj ++ )
346                 bmask[jj] = (temp[jj-ii] != 0.0) ;
347 
348          } /* end of loop over space (voxels) */
349 
350    if (nodemax >= 0) {
351       byte *m2 = NULL;
352       if (nodemax+1 < CALC_nvox) {
353          fprintf(stderr,   "** Gosh darn it Cleatus! Now how on earth can that happen?\n"
354                            "   nodemax(+1) is %d(+1) and CALC_nvox is %d\n"
355                            "   Al Gore was right now, all along!\n", nodemax, CALC_nvox);
356          CALC_nvox = 0;
357          free(bmask); bmask=NULL;
358          goto CLEANUP;
359       }
360       m2 = (byte *) calloc(nodemax+1, sizeof(byte)) ;
361       for (zii=0; zii<CALC_nvox; ++zii) {
362          if (bmask[zii]) m2[node_dset->dblk->node_list[zii]] = 1;
363       }
364       free(bmask);
365       bmask = m2; m2 = NULL;
366       CALC_nvox = nodemax+1;
367    }
368 
369    /* cleanup and go home */
370    CLEANUP:
371    for( ids=0 ; ids < 26 ; ids++ ){
372       free(atoz[ids]) ;
373       if( CALC_dset[ids] != NULL ) DSET_delete( CALC_dset[ids] ) ;
374    }
375    DSET_delete(new_dset) ;
376    free(CALC_code) ;
377    if( nxyz != NULL ) *nxyz = CALC_nvox ;
378    RETURN( bmask );
379 }
380 
381 /*--------------------------------------------------------------------
382    read the arguments, load the global variables
383 ----------------------------------------------------------------------*/
384 
CALC_read_opts(int argc,char * argv[])385 static int CALC_read_opts( int argc , char * argv[] )
386 {
387    int nopt = 0 ;
388    int ids ;
389    int ii ;
390 
391 ENTRY("CALC_read_opts") ;
392 
393    CALC_nvox  = -1 ;
394    CALC_code  = NULL ;
395    CALC_dshift_mode_current = DSHIFT_MODE_STOP ;
396    CALC_has_predefined = 0 ;
397 
398    for( ids=0 ; ids < 26 ; ids++ ){
399       CALC_dset[ids]   = NULL ;
400       CALC_type[ids]   = -1 ;
401 
402       CALC_dshift[ids]      = -1 ;                        /* 22 Nov 1999 */
403       CALC_dshift_mode[ids] = CALC_dshift_mode_current ;
404    }
405 
406    while( nopt < argc && argv[nopt][0] == '-' ){
407 
408       /**** -expr expression ****/
409 
410       if( strncmp(argv[nopt],"-expr",4) == 0 ){
411          if( CALC_code != NULL ){
412             fprintf(stderr,
413              "** -cmask: cannot have 2 -expr options!\n") ; RETURN(1) ;
414          }
415          nopt++ ;
416          if( nopt >= argc ){
417             fprintf(stderr,
418              "** -cmask: need argument after -expr!\n") ; RETURN(1) ;
419          }
420          CALC_code = PARSER_generate_code( argv[nopt++] ) ;
421          if( CALC_code == NULL ){
422             fprintf(stderr,
423              "** -cmask: illegal expression!\n") ; RETURN(1) ;
424          }
425          PARSER_mark_symbols( CALC_code , CALC_has_sym ) ; /* 15 Sep 1999 */
426          continue ;
427       }
428 
429       /**** -dsSTOP [22 Nov 1999] ****/
430 
431       if( strncmp(argv[nopt],"-dsSTOP",6) == 0 ){
432          CALC_dshift_mode_current = DSHIFT_MODE_STOP ;
433          nopt++ ; continue ;
434       }
435 
436       /**** -dsWRAP [22 Nov 1999] ****/
437 
438       if( strncmp(argv[nopt],"-dsWRAP",6) == 0 ){
439          CALC_dshift_mode_current = DSHIFT_MODE_WRAP ;
440          nopt++ ; continue ;
441       }
442 
443       /**** -dsZERO [22 Nov 1999] ****/
444 
445       if( strncmp(argv[nopt],"-dsZERO",6) == 0 ){
446          CALC_dshift_mode_current = DSHIFT_MODE_ZERO ;
447          nopt++ ; continue ;
448       }
449 
450       /**** -<letter> dataset ****/
451 
452       ids = strlen( argv[nopt] ) ;
453 
454       if( (argv[nopt][1] >= 'a' && argv[nopt][1] <= 'z') && ids == 2 ) {
455 
456          int ival , nxyz , ll ;
457          THD_3dim_dataset * dset ;
458 
459          ival = argv[nopt][1] - 'a' ;
460          if( VAR_DEFINED(ival) ){
461             fprintf(stderr,
462              "** -cmask: Can't define %c symbol twice\n",argv[nopt][1]);
463             RETURN(1) ;
464          }
465 
466          nopt++ ;
467          if( nopt >= argc ){
468             fprintf(stderr,
469              "** -cmask: need argument after %s\n",argv[nopt-1]);
470             RETURN(1) ;
471          }
472 
473          /*-- 22 Nov 1999: allow for a differentially
474                            subscripted name, as in "-b a[1,0,0,0]" --*/
475 
476          ll = strlen(argv[nopt]) ;
477          if( (argv[nopt][0] >= 'a' && argv[nopt][0] <= 'z') &&  /* legal name */
478              ( (ll >= 3 && argv[nopt][1] == '[') ||             /* subscript */
479                (ll == 3 &&                                      /*    OR    */
480                 (argv[nopt][1] == '+' || argv[nopt][1] == '-')) /* +- ijkl */
481              ) ){
482 
483             int jds = argv[nopt][0] - 'a' ;  /* actual dataset index */
484             int * ijkl ;                     /* array of subscripts */
485 
486             if( CALC_dset[jds] == NULL ){
487                fprintf(stderr,
488                 "** -cmask: Must define dataset %c before using it in %s\n",
489                        argv[nopt][0] , argv[nopt] ) ;
490                RETURN(1) ;
491             }
492 
493             /*- get subscripts -*/
494 
495             if( argv[nopt][1] == '[' ){            /* format is [i,j,k,l] */
496                MCW_intlist_allow_negative(1) ;
497                ijkl = MCW_get_intlist( 9999 , argv[nopt]+1 ) ;
498                MCW_intlist_allow_negative(0) ;
499                if( ijkl == NULL || ijkl[0] != 4 ){
500                   fprintf(stderr,
501                    "** -cmask: Illegal differential subscripting %s\n",
502                                  argv[nopt] ) ;
503                   RETURN(1) ;
504                }
505             } else {                               /* format is +i, -j, etc */
506                 ijkl = (int *) malloc( sizeof(int) * 5 ) ;
507                 ijkl[1] = ijkl[2] = ijkl[3] = ijkl[4] = 0 ;  /* initialize */
508                 switch( argv[nopt][2] ){
509                    default:
510                       fprintf(stderr,
511                        "** -cmask: Bad differential subscripting %s\n",
512                                  argv[nopt] ) ;
513                    RETURN(1) ;
514 
515                    case 'i': ijkl[1] = (argv[nopt][1]=='+') ? 1 : -1 ; break ;
516                    case 'j': ijkl[2] = (argv[nopt][1]=='+') ? 1 : -1 ; break ;
517                    case 'k': ijkl[3] = (argv[nopt][1]=='+') ? 1 : -1 ; break ;
518                    case 'l': ijkl[4] = (argv[nopt][1]=='+') ? 1 : -1 ; break ;
519                 }
520             }
521 
522             if( ijkl[4] != 0 ){  /* disallow time subscripting */
523                fprintf(stderr,
524                 "++ -cmask: Warning: differential time shifting %s not allowed\n",
525                        argv[nopt] ) ;
526                ijkl[4] = 0 ;
527             }
528 
529             /*- more sanity checks -*/
530 
531             if( ijkl[1]==0 && ijkl[2]==0 && ijkl[3]==0 && ijkl[4]==0 ){
532                fprintf(stderr,
533                 "++ -cmask: differential subscript %s is all zero\n",
534                        argv[nopt] ) ;
535             }
536 
537             /*- set values for later use -*/
538 
539             CALC_dshift  [ival] = jds ;
540             CALC_dshift_i[ival] = ijkl[1] ;
541             CALC_dshift_j[ival] = ijkl[2] ;
542             CALC_dshift_k[ival] = ijkl[3] ;
543             CALC_dshift_l[ival] = ijkl[4] ;
544 
545             CALC_dshift_mode[ival] = CALC_dshift_mode_current ;
546 
547             /*- time to trot, Bwana -*/
548 
549             free(ijkl) ; nopt++ ; goto DSET_DONE ;
550 
551          } /* end of _dshift */
552 
553          /*-- meanwhile, back at the "normal" dataset opening ranch --*/
554 
555          { char dname[512] ;                               /* 02 Nov 1999 */
556 
557            if( strstr(argv[nopt],"[") == NULL ){
558               sprintf(dname,"%s[0]",argv[nopt++]) ;        /* add [0] */
559            } else {
560               strcpy(dname,argv[nopt++]) ;                 /* don't mangle */
561            }
562            dset = THD_open_dataset( dname ) ;              /* open it */
563            if( dset == NULL ){
564               fprintf(stderr,
565                "** -cmask: can't open dataset %s\n",dname) ;
566               RETURN(1) ;
567            }
568          }
569          CALC_dset[ival] = dset ;
570 
571          /* set some parameters based on the dataset */
572 
573          nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz ;
574          if( CALC_nvox < 0 ){
575             CALC_nvox = nxyz ;
576          } else if( nxyz != CALC_nvox ){
577             fprintf(stderr,
578              "** -cmask: dataset %s differs in size from others\n",argv[nopt-1]);
579             RETURN(1) ;
580          }
581 
582          CALC_type[ival] = DSET_BRICK_TYPE(dset,0) ;
583 
584          /* load floating scale factors */
585 
586          CALC_ffac[ival] = DSET_BRICK_FACTOR(dset,0) ;
587          if (CALC_ffac[ival] == 0.0 ) CALC_ffac[ival] = 1.0 ;
588 
589          /* read data from disk */
590 
591          THD_load_datablock( dset->dblk ) ;
592          if( ! DSET_LOADED(dset) ){
593             fprintf(stderr,
594              "** -cmask: Can't read data brick for dataset %s\n",argv[nopt-1]) ;
595             RETURN(1) ;
596          }
597 
598          /* set pointers for actual dataset arrays */
599 
600          switch (CALC_type[ival]) {
601              case MRI_short:
602                 CALC_short[ival] = (short *) DSET_ARRAY(dset,0) ;
603              break;
604 
605              case MRI_float:
606                 CALC_float[ival] = (float *) DSET_ARRAY(dset,0) ;
607              break;
608 
609              case MRI_byte:
610                 CALC_byte[ival] = (byte *) DSET_ARRAY(dset,0) ;
611              break;
612 
613           } /* end of switch over type */
614 
615 DSET_DONE: continue;
616 
617       } /* end of dataset input */
618 
619       fprintf(stderr,"** -cmask: Unknown option: %s\n",argv[nopt]) ;
620       RETURN(1) ;
621 
622    }  /* end of loop over options */
623 
624    /*---------------------------------------*/
625    /*** cleanup: check for various errors ***/
626 
627    for( ids=0 ; ids < 26 ; ids++ ) if( CALC_dset[ids] != NULL ) break ;
628    if( ids == 26 ){
629       fprintf(stderr,
630        "** -cmask: No actual input datasets given!\n") ;
631       RETURN(1) ;
632    }
633 
634    if( CALC_code == NULL ){
635       fprintf(stderr,"** -cmask: No expression given!\n") ;
636       RETURN(1) ;
637    }
638 
639    /* 15 Apr 1999: check if each input dataset is used,
640                    or if an undefined symbol is used.   */
641 
642    for (ids=0; ids < 26; ids ++){
643       if( VAR_DEFINED(ids) && !CALC_has_sym[ids] )
644          fprintf(stderr ,
645           "++ -cmask: input '%c' is not used in the expression\n" ,
646               abet[ids] ) ;
647 
648       else if( !VAR_DEFINED(ids) && CALC_has_sym[ids] ){
649 
650          if( ((1<<ids) & PREDEFINED_MASK) == 0 )
651             fprintf(stderr ,
652              "++ -cmask: symbol %c is used but not defined\n" ,
653                     abet[ids] ) ;
654          else {
655             CALC_has_predefined++ ;
656          }
657       }
658    }
659 
660    RETURN(0) ;
661 }
662