00014 {
00015    int iarg , pos = 0 ;
00016    float thresh=0.0 ;
00017    MRI_IMAGE * maskim=NULL , *imin , *imout ;
00018    float * maskar ;
00019    int nxim , nyim , ii , npix ;
00020 
00021    if( argc < 3 || strncmp(argv[1],"-help",4) == 0 ){
00022       printf("Usage: immask [-thresh #] [-mask mask_image] [-pos] input_image output_image\n"
00023              "* Masks the input_image and produces the output_image;\n"
00024              "* Use of -thresh # means all pixels with absolute value below # in\n"
00025              "   input_image will be set to zero in the output_image\n"
00026              "* Use of -mask mask_image means that only locations that are nonzero\n"
00027              "   in the mask_image will be nonzero in the output_image\n"
00028              "* Use of -pos means only positive pixels from input_image will be used\n"
00029              "* At least one of -thresh, -mask, -pos must be used; more than one is OK.\n"
00030             ) ;
00031      exit(0) ;
00032    }
00033 
00034    machdep() ;
00035 
00036    iarg = 1 ;
00037    while( iarg < argc && argv[iarg][0] == '-' ){
00038 
00039       /*** -pos ***/
00040 
00041       if( strncmp(argv[iarg],"-pos",4) == 0 ){
00042          pos = 1 ;
00043          iarg++ ; continue ;
00044       }
00045 
00046       /*** -thresh # ***/
00047 
00048       if( strncmp(argv[iarg],"-thresh",5) == 0 ){
00049          thresh = strtod( argv[++iarg] , NULL ) ;
00050          if( iarg >= argc || thresh <= 0.0 ){
00051             fprintf(stderr,"Illegal -thresh!\a\n") ; exit(1) ;
00052          }
00053          iarg++ ; continue ;
00054       }
00055 
00056       if( strncmp(argv[iarg],"-mask",5) == 0 ){
00057          maskim = mri_read_just_one( argv[++iarg] ) ;
00058          if( maskim == NULL || iarg >= argc || ! MRI_IS_2D(maskim) ){
00059             fprintf(stderr,"Illegal -mask!\a\n") ; exit(1) ;
00060          }
00061          if( maskim->kind != MRI_float ){
00062             imin = mri_to_float( maskim ) ;
00063             mri_free( maskim ) ;
00064             maskim = imin ;
00065          }
00066          iarg++ ; continue ;
00067       }
00068 
00069       fprintf(stderr,"** Illegal option: %s\a\n",argv[iarg]) ;
00070       exit(1) ;
00071    }
00072    if( thresh <= 0.0 && maskim == NULL && pos == 0 ){
00073       fprintf(stderr,"No -thresh, -mask, -pos ==> can't go on!\a\n") ; exit(1) ;
00074    }
00075    if( iarg+1 >= argc ){
00076       fprintf(stderr,"Must have input_image and output_image on command line!\a\n") ;
00077       exit(1) ;
00078    }
00079 
00080    imin = mri_read_just_one( argv[iarg++] ) ;
00081    if( imin == NULL ) exit(1) ;
00082    if( ! MRI_IS_2D(imin) ){
00083       fprintf(stderr,"can only process 2D images!\a\n") ;
00084       exit(1) ;
00085    }
00086 
00087    nxim = imin->nx ;
00088    nyim = imin->ny ;
00089    npix = nxim * nyim ;
00090 
00091    if( maskim == NULL ){
00092       maskim = mri_new( nxim , nyim , MRI_float ) ;
00093       maskar = MRI_FLOAT_PTR(maskim) ;
00094       for( ii=0 ; ii < npix ; ii++ ) maskar[ii] = 1.0 ;
00095    } else if( maskim->nx != nxim || maskim->ny != nyim ){
00096       fprintf(stderr,"Mask and input image not same size!\a\n") ;
00097       exit(1) ;
00098    } else {
00099       maskar = MRI_FLOAT_PTR(maskim) ;
00100    }
00101    imout = mri_new( nxim , nyim , imin->kind ) ;
00102 
00103    switch( imin->kind ){
00104 
00105       default:
00106          fprintf(stderr,"Unrecognized input image type!\a\n") ;
00107          exit(1) ;
00108 
00109       case MRI_byte:{
00110          byte * arin , * arout , val ;
00111          arin  = mri_data_pointer(imin) ;
00112          arout = mri_data_pointer(imout) ;
00113          for( ii=0 ; ii < npix ; ii++ ){
00114             val = arin[ii] ;
00115             if( maskar[ii] != 0.0 && ABS(val) >= thresh ) arout[ii] = val ;
00116             else                                          arout[ii] = 0 ;
00117          }
00118       } break ;
00119 
00120       case MRI_short:{
00121          short * arin , * arout , val ;
00122          arin  = mri_data_pointer(imin) ;
00123          arout = mri_data_pointer(imout) ;
00124          for( ii=0 ; ii < npix ; ii++ ){
00125             val = arin[ii] ;
00126             if( maskar[ii] != 0.0 && ABS(val) >= thresh ) arout[ii] = val ;
00127             else                                          arout[ii] = 0 ;
00128          }
00129          if( pos ) for( ii=0 ; ii < npix ; ii++ ) if( arout[ii] < 0 ) arout[ii] = 0 ;
00130       } break ;
00131 
00132       case MRI_float:{
00133          float * arin , * arout , val ;
00134          arin  = mri_data_pointer(imin) ;
00135          arout = mri_data_pointer(imout) ;
00136          for( ii=0 ; ii < npix ; ii++ ){
00137             val = arin[ii] ;
00138             if( maskar[ii] != 0.0 && ABS(val) >= thresh ) arout[ii] = val ;
00139             else                                          arout[ii] = 0 ;
00140          }
00141          if( pos ) for( ii=0 ; ii < npix ; ii++ ) if( arout[ii] < 0 ) arout[ii] = 0 ;
00142       } break ;
00143 
00144       case MRI_int:{
00145          int * arin , * arout , val ;
00146          arin  = mri_data_pointer(imin) ;
00147          arout = mri_data_pointer(imout) ;
00148          for( ii=0 ; ii < npix ; ii++ ){
00149             val = arin[ii] ;
00150             if( maskar[ii] != 0.0 && ABS(val) >= thresh ) arout[ii] = val ;
00151             else                                          arout[ii] = 0 ;
00152          }
00153          if( pos ) for( ii=0 ; ii < npix ; ii++ ) if( arout[ii] < 0 ) arout[ii] = 0 ;
00154       } break ;
00155 
00156       case MRI_double:{
00157          double * arin , * arout , val ;
00158          arin  = mri_data_pointer(imin) ;
00159          arout = mri_data_pointer(imout) ;
00160          for( ii=0 ; ii < npix ; ii++ ){
00161             val = arin[ii] ;
00162             if( maskar[ii] != 0.0 && ABS(val) >= thresh ) arout[ii] = val ;
00163             else                                          arout[ii] = 0 ;
00164          }
00165          if( pos ) for( ii=0 ; ii < npix ; ii++ ) if( arout[ii] < 0 ) arout[ii] = 0 ;
00166       } break ;
00167 
00168       case MRI_complex:{
00169          complex * arin , * arout , val ;
00170          arin  = mri_data_pointer(imin) ;
00171          arout = mri_data_pointer(imout) ;
00172          for( ii=0 ; ii < npix ; ii++ ){
00173             val = arin[ii] ;
00174             if( maskar[ii] != 0.0 && CABS(val) >= thresh ) arout[ii] = val ;
00175             else                                           arout[ii] = CMPLX(0,0) ;
00176          }
00177       } break ;
00178    }
00179 
00180    mri_write( argv[iarg] , imout ) ;
00181    exit(0) ;
00182 }