00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #include "mrilib.h"
00008 
00009 /*-- 14 Oct 1999: modified to allow for inverse warping --*/
00010 
00011 /*-- 18 Oct 1999: modified to implement -preserve option --*/
00012 
00013 THD_fvec3 AFNI_backward_warp_vector( THD_warp * warp , THD_fvec3 old_fv ) ;
00014 
00015 int main( int argc , char * argv[] )
00016 {
00017    char *prefix = "fractionize" ;
00018    THD_3dim_dataset *tset=NULL , *iset=NULL , *dset=NULL , *wset=NULL ;
00019    int iarg=1 ;
00020    int nxin,nyin,nzin , nxyin ;
00021    float dxin,dyin,dzin , xorgin,yorgin,zorgin , clip=0.0 ;
00022    int nxout,nyout,nzout , nxyout , nvoxout ;
00023    float dxout,dyout,dzout , xorgout,yorgout,zorgout ;
00024    float f1,f2,f , g1,g2,g , h1,h2,h , sx,sy,sz , tx,ty,tz ;
00025    float x1,x2 , y1,y2 , z1,z2 , xx1,xx2 , yy1,yy2 , zz1,zz2 ;
00026    int i,ip , j,jp , k,kp , iv,jv,kv , vtype , ijk ;
00027    float *voxout ;
00028    byte  *bin   = NULL ;
00029    short *sin   = NULL , sclip=0   ;
00030    float *fin   = NULL , fclip=0.00001 ;
00031    void  *voxin = NULL ;
00032    THD_fvec3 vv ;
00033 
00034    THD_warp *warp=NULL ;  /* 14 Oct 1999 */
00035 
00036    int do_vote=0 ;          /* 18 Oct 1999 */
00037    int *vote_val = NULL ;
00038    int  nvote_val , ivote , voter , vote_print=0 ;
00039    byte  *vote_bout = NULL ;
00040    short *vote_sout = NULL ;
00041    float *vote_best = NULL ;
00042 
00043    /*-- help? --*/
00044 
00045    if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00046       printf("Usage: 3dfractionize [options]\n"
00047              "\n"
00048              "* For each voxel in the output dataset, computes the fraction\n"
00049              "    of it that is occupied by nonzero voxels from the input.\n"
00050              "* The fraction is stored as a short in the range 0..10000,\n"
00051              "    indicating fractions running from 0..1.\n"
00052              "* The template dataset is used only to define the output grid;\n"
00053              "    its brick(s) will not be read into memory.  (The same is\n"
00054              "    true of the warp dataset, if it is used.)\n"
00055              "* The actual values stored in the input dataset are irrelevant,\n"
00056              "    except in that they are zero or nonzero (UNLESS the -preserve\n"
00057              "    option is used).\n"
00058              "\n"
00059              "The purpose of this program is to allow the resampling of a mask\n"
00060              "dataset (the input) from a fine grid to a coarse grid (defined by\n"
00061              "the template).  When you are using the output, you will probably\n"
00062              "want to threshold the mask so that voxels with a tiny occupancy\n"
00063              "fraction aren't used.  This can be done in 3dmaskave, by using\n"
00064              "3calc, or with the '-clip' option below.\n"
00065              "\n"
00066              "Options are [the first 2 are 'mandatory options']:\n"
00067              "  -template tset  = Use dataset 'tset' as a template for the output.\n"
00068              "                      The output dataset will be on the same grid as\n"
00069              "                      this dataset.\n"
00070              "\n"
00071              "  -input iset     = Use dataset 'iset' for the input.\n"
00072              "                      Only the sub-brick #0 of the input is used.\n"
00073              "                      You can use the sub-brick selection technique\n"
00074              "                      described in '3dcalc -help' to choose the\n"
00075              "                      desired sub-brick from a multi-brick dataset.\n"
00076              "\n"
00077              "  -prefix ppp     = Use 'ppp' for the prefix of the output.\n"
00078              "                      [default prefix = 'fractionize']\n"
00079              "\n"
00080              "  -clip fff       = Clip off voxels that are less than 'fff' occupied.\n"
00081              "                      'fff' can be a number between 0.0 and 1.0, meaning\n"
00082              "                      the fraction occupied, can be a number between 1.0\n"
00083              "                      and 100.0, meaning the percent occupied, or can be\n"
00084              "                      a number between 100.0 and 10000.0, meaning the\n"
00085              "                      direct output value to use as a clip level.\n"
00086              "                   ** Some sort of clipping is desirable; otherwise,\n"
00087              "                        an output voxel that is barely overlapped by a\n"
00088              "                        single nonzero input voxel will enter the mask.\n"
00089              "                      [default clip = 0.0]\n"
00090              "\n"
00091              "  -warp wset      = If this option is used, 'wset' is a dataset that\n"
00092              "                      provides a transformation (warp) from +orig\n"
00093              "                      coordinates to the coordinates of 'iset'.\n"
00094              "                      In this case, the output dataset will be in\n"
00095              "                      +orig coordinates rather than the coordinates\n"
00096              "                      of 'iset'.  With this option:\n"
00097              "                   ** 'tset' must be in +orig coordinates\n"
00098              "                   ** 'iset' must be in +acpc or +tlrc coordinates\n"
00099              "                   ** 'wset' must be in the same coordinates as 'iset'\n"
00100              "\n"
00101              "  -preserve       = When this option is used, the program will copy\n"
00102              "     or               the nonzero values of input voxels to the output\n"
00103              "  -vote               dataset, rather than create a fractional mask.\n"
00104              "                      Since each output voxel might be overlapped\n"
00105              "                      by more than one input voxel, the program 'votes'\n"
00106              "                      for which input value to preserve.  For example,\n"
00107              "                      if input voxels with value=1 occupy 10%% of an\n"
00108              "                      output voxel, and inputs with value=2 occupy 20%%\n"
00109              "                      of the same voxel, then the output value in that\n"
00110              "                      voxel will be set to 2 (provided that 20%% is >=\n"
00111              "                      to the clip fraction).\n"
00112              "                   ** Voting can only be done on short-valued datasets,\n"
00113              "                        or on byte-valued datasets.\n"
00114              "                   ** Voting is a relatively time-consuming option,\n"
00115              "                        since a separate loop is made through the\n"
00116              "                        input dataset for each distinct value found.\n"
00117              "                   ** Combining this with the -warp option does NOT\n"
00118              "                        make a general +trlc to +orig transformer!\n"
00119              "                        This is because for any value to survive the\n"
00120              "                        vote, its fraction in the output voxel must be\n"
00121              "                        >= clip fraction, regardless of other values\n"
00122              "                        present in the output voxel.\n"
00123              "\n"
00124              "Example usage:\n"
00125              " 3dfractionize -template a+orig -input b+tlrc -warp anat+tlrc -clip 0.2\n"
00126              "\n"
00127              "This program will also work in going from a coarse grid to a fine grid,\n"
00128              "but it isn't clear that this capability has any purpose.\n"
00129              "-- RWCox - February 1999\n"
00130              "         - October 1999: added -warp and -preserve options\n"
00131             ) ;
00132       exit(0) ;
00133    }
00134 
00135    mainENTRY("3dfractionize main"); machdep();
00136    AFNI_logger("3dfractionize",argc,argv); PRINT_VERSION("3dfractionize");
00137 
00138    /*-- read command line args --*/
00139 
00140    while( iarg < argc && argv[iarg][0] == '-' ){
00141 
00142       if( strcmp(argv[iarg],"-clip") == 0 ){
00143          clip = strtod( argv[++iarg] , NULL ) ;
00144 
00145               if( clip <= 1.0     ) sclip = (short)(10000.0 * clip + 0.49999) ;
00146          else if( clip <= 100.0   ) sclip = (short)(  100.0 * clip + 0.49999) ;
00147          else if( clip <= 10000.0 ) sclip = (short)(clip) ;
00148          else                       sclip = -1 ;
00149 
00150          if( sclip < 0 || sclip > 10000 ){
00151             fprintf(stderr,"** Illegal clip value!\n") ; exit(1) ;
00152          }
00153 
00154          fclip = 0.0001 * sclip ;
00155          iarg++ ; continue ;
00156       }
00157 
00158       if( strcmp(argv[iarg],"-template") == 0 || strcmp(argv[iarg],"-tset") == 0 ){
00159          if( tset != NULL ){
00160             fprintf(stderr,"** Can't have more than one -template argument!\n") ;
00161             exit(1) ;
00162          }
00163          tset = THD_open_one_dataset( argv[++iarg] ) ;
00164          if( tset == NULL ){
00165             fprintf(stderr,"** Can't open template %s\n",argv[iarg]) ;
00166             exit(1) ;
00167          }
00168          iarg++ ; continue ;
00169       }
00170 
00171       if( strcmp(argv[iarg],"-input") == 0 || strcmp(argv[iarg],"-iset") == 0 ){
00172          if( iset != NULL ){
00173             fprintf(stderr,"** Can't have more than one -input argument!\n") ;
00174             exit(1) ;
00175          }
00176          iset = THD_open_dataset( argv[++iarg] ) ;
00177          if( iset == NULL ){
00178             fprintf(stderr,"** Can't open input %s\n",argv[iarg]) ;
00179             exit(1) ;
00180          }
00181          iarg++ ; continue ;
00182       }
00183 
00184       if( strcmp(argv[iarg],"-warp") == 0 || strcmp(argv[iarg],"-wset") == 0 ){
00185          if( wset != NULL ){
00186             fprintf(stderr,"** Can't have more than one -warp argument!\n") ;
00187             exit(1) ;
00188          }
00189          wset = THD_open_dataset( argv[++iarg] ) ;
00190          if( wset == NULL ){
00191             fprintf(stderr,"** Can't open warp %s\n",argv[iarg]) ;
00192             exit(1) ;
00193          }
00194          iarg++ ; continue ;
00195       }
00196 
00197       if( strcmp(argv[iarg],"-prefix") == 0 ){
00198          prefix = argv[++iarg] ;
00199          iarg++ ; continue ;
00200       }
00201 
00202       if( strcmp(argv[iarg],"-preserve") == 0 || strcmp(argv[iarg],"-vote") == 0 ){
00203          do_vote = 1 ;
00204          iarg++ ; continue ;
00205       }
00206 
00207       fprintf(stderr,"** Illegal option: %s\n",argv[iarg]) ; exit(1) ;
00208    }
00209 
00210    /*-- check inputs for sanity --*/
00211 
00212    if( tset == NULL ){
00213       fprintf(stderr,"** No template dataset?\n") ; exit(1) ;
00214    }
00215    if( iset == NULL ){
00216       fprintf(stderr,"** No input dataset?\n") ; exit(1) ;
00217    }
00218    if( !THD_filename_ok(prefix) ){
00219       fprintf(stderr,"** Illegal prefix?\n") ; exit(1) ;
00220    }
00221 
00222    /*- 14 Oct 1999: additional checks with the -warp option -*/
00223 
00224    if( wset != NULL ){
00225 
00226       if( tset->view_type != VIEW_ORIGINAL_TYPE ){
00227          fprintf(stderr,"** Template is not in +orig view - this is illegal!\n");
00228          exit(1);
00229       }
00230 
00231       if( iset->view_type == VIEW_ORIGINAL_TYPE ){
00232          fprintf(stderr,"** Input is in +orig view - this is illegal!\n"); exit(1);
00233       }
00234 
00235       if( wset != NULL && wset->view_type != iset->view_type ){
00236          fprintf(stderr,"** Warp and Input are not in same view - this is illegal!\n");
00237          exit(1);
00238       }
00239 
00240       warp = wset->warp ;
00241       if( warp == NULL ){
00242          fprintf(stderr,"** Warp dataset doesn't contain a warp transformation!\n");
00243          exit(1) ;
00244       }
00245    }
00246 
00247    /*- 18 Oct 1999: check input for legality if we are voting --*/
00248 
00249    if( do_vote ){
00250      vtype = DSET_BRICK_TYPE(iset,0) ;
00251      if( vtype != MRI_short && vtype != MRI_byte ){
00252        fprintf(stderr,"** -preserve option requires short- or byte-valued input dataset!\n");
00253        exit(1);
00254      }
00255    }
00256 
00257    /*-- start to create output dataset --*/
00258 
00259    dset = EDIT_empty_copy( tset ) ;
00260 
00261    tross_Copy_History( iset , dset ) ;                          /* 14 Oct 1999 */
00262    tross_Make_History( "3dfractionize" , argc,argv , dset ) ;
00263 
00264    EDIT_dset_items( dset ,
00265                        ADN_prefix    , prefix ,
00266                        ADN_nvals     , 1 ,
00267                        ADN_ntt       , 0 ,
00268                        ADN_brick_fac , NULL ,
00269                        ADN_datum_all , MRI_short ,
00270                     ADN_none ) ;
00271 
00272    if( ISFUNC(dset) )
00273       EDIT_dset_items( dset , ADN_func_type,FUNC_FIM_TYPE , ADN_none ) ;
00274 
00275    if( THD_is_file(dset->dblk->diskptr->header_name) ){
00276       fprintf(stderr,
00277               "** Output file %s already exists -- cannot continue!\n",
00278               dset->dblk->diskptr->header_name ) ;
00279       exit(1) ;
00280    }
00281 
00282    DSET_delete( tset ) ;  /* don't need template any more */
00283 
00284    /* load the input dataset */
00285 
00286    DSET_load( iset ) ;
00287    if( ! DSET_LOADED(iset) ){
00288       fprintf(stderr,"** Can't read input dataset brick!\n") ; exit(1) ;
00289    }
00290    voxin = DSET_ARRAY(iset,0) ; vtype = DSET_BRICK_TYPE(iset,0) ;
00291    switch( vtype ){
00292       default: fprintf(stderr,"** Illegal brick type in input dataset!\n"); exit(1);
00293 
00294       case MRI_byte:  bin = (byte * ) voxin ; break ;
00295       case MRI_short: sin = (short *) voxin ; break ;
00296       case MRI_float: fin = (float *) voxin ; break ;
00297    }
00298 
00299    /*-- setup voxel index (iv,jv,kv) to coordinate (x,y,z) calculations --*/
00300 
00301    nxin   =iset->daxes->nxx  ; nyin   =iset->daxes->nyy  ; nzin   =iset->daxes->nzz  ;
00302    dxin   =iset->daxes->xxdel; dyin   =iset->daxes->yydel; dzin   =iset->daxes->zzdel;
00303    xorgin =iset->daxes->xxorg; yorgin =iset->daxes->yyorg; zorgin =iset->daxes->zzorg;
00304 
00305    nxout  =dset->daxes->nxx  ; nyout  =dset->daxes->nyy  ; nzout  =dset->daxes->nzz  ;
00306    dxout  =dset->daxes->xxdel; dyout  =dset->daxes->yydel; dzout  =dset->daxes->zzdel;
00307    xorgout=dset->daxes->xxorg; yorgout=dset->daxes->yyorg; zorgout=dset->daxes->zzorg;
00308 
00309    /* voxel (i,j,k) center is at (xorg+i*dx,yorg+j*dy,zorg+k*dz) */
00310 
00311    nxyout = nxout * nyout ; nvoxout = nxyout * nzout ;
00312 
00313    voxout = (float *) malloc( sizeof(float) * nvoxout ) ;
00314 
00315    nxyin = nxin * nyin ;
00316 
00317    /*-- if voting, do some setup --*/
00318 
00319    if( do_vote ){
00320       int nvoxin = nxyin * nzin ;
00321 
00322       /* 1: find all distinct nonzero values in dataset */
00323 
00324       vote_val  = (int *) malloc( sizeof(int) ) ;
00325       nvote_val = 0 ;
00326       for( iv=0 ; iv < nvoxin ; iv++ ){
00327          kv = (vtype == MRI_short) ? sin[iv] : bin[iv] ;   /* voxel value */
00328          if( kv == 0 ) continue ;                          /* skip zeroes */
00329          for( jv=0 ; jv < nvote_val ; jv++ )              /* find in list */
00330             if( vote_val[jv] == kv ) break ;
00331          if( jv < nvote_val ) continue ;       /* skip if already in list */
00332 
00333          vote_val = (int *) realloc( vote_val , sizeof(int)*(nvote_val+1) ) ;
00334          vote_val[nvote_val++] = kv ;
00335       }
00336       if( nvote_val == 0 ){
00337          fprintf(stderr,"** Input dataset is all zero!\n") ; exit(1) ;
00338       }
00339       fprintf(stderr,"++ Found %d distinct nonzero values in input.\n",nvote_val) ;
00340 
00341       if( nvote_val > 1 )                    /* arrange into ascending value */
00342          qsort_int( nvote_val , vote_val ) ;
00343 
00344       /* 2: make a receptacle for the voting results */
00345 
00346       if( vtype == MRI_byte ){
00347           vote_bout = (byte *) malloc( sizeof(byte) * nvoxout ) ;
00348           memset( vote_bout , 0 , sizeof(byte) * nvoxout ) ;
00349       } else {
00350           vote_sout = (short *) malloc( sizeof(short) * nvoxout ) ;
00351           memset( vote_sout , 0 , sizeof(short) * nvoxout ) ;
00352       }
00353 
00354       /* 3: make a place to hold the best fraction yet */
00355 
00356       vote_best = (float *) malloc( sizeof(float) * nvoxout ) ;
00357       for( i=0 ; i < nvoxout ; i++ ) vote_best[i] = 0.0 ;
00358 
00359       /* 4: if needed, attach the correct brick factor to the output */
00360 
00361       f = DSET_BRICK_FACTOR(iset,0) ;
00362       if( f != 0.0 && f != 1.0 ) EDIT_BRICK_FACTOR(dset,0,f) ;
00363 
00364       /* 5: if input is bytes, make the output be that too */
00365 
00366       if( vtype == MRI_byte )
00367          EDIT_dset_items( dset , ADN_datum_all , MRI_byte , ADN_none ) ;
00368    }
00369 
00370    /*-- loop over values to vote on --*/
00371 
00372    ivote = 0 ;
00373    do{
00374 
00375       if( do_vote ){
00376          voter = vote_val[ivote] ;  /* who gets to vote this time */
00377          if( ivote%10 == 9 ){
00378             if( !vote_print ) fprintf(stderr,"++ Voting:") ;
00379             fprintf(stderr,"%d",(ivote/10)%10) ;
00380             vote_print++ ;
00381          }
00382       }
00383 
00384       for( i=0 ; i < nvoxout ; i++ ) voxout[i] = 0.0 ;  /* fractions */
00385 
00386       /*-- loop over input voxels --*/
00387 
00388       for( kv=0 ; kv < nzin ; kv++ ){
00389 
00390        z1 = zorgin + dzin * (kv-0.5) ; z2 = zorgin + dzin * (kv+0.49999) ;
00391 
00392        for( jv=0 ; jv < nyin ; jv++ ){
00393 
00394         y1 = yorgin + dyin * (jv-0.5) ; y2 = yorgin + dyin * (jv+0.49999) ;
00395 
00396         for( iv=0 ; iv < nxin ; iv++ ){
00397 
00398          ijk = iv + jv*nxin + kv*nxyin ;  /* 1D index of voxel (iv,jv,kv) */
00399 
00400          /* decide if we use this voxel, based on its value */
00401 
00402          if( do_vote ){
00403                  if( vtype == MRI_short && sin[ijk] != voter ) continue ;
00404             else if( vtype == MRI_byte  && bin[ijk] != voter ) continue ;
00405          } else {
00406                  if( vtype == MRI_short && sin[ijk] == 0 ) continue ;
00407             else if( vtype == MRI_byte  && bin[ijk] == 0 ) continue ;
00408             else if( vtype == MRI_float && fin[ijk] == 0 ) continue ;
00409          }
00410 
00411          x1 = xorgin + dxin * (iv-0.5) ; x2 = xorgin + dxin * (iv+0.49999) ;
00412 
00413          /* input voxel (iv,jv,kv) spans coordinates [x1,x2] X [y1,y2] X [z1,z2] */
00414 
00415          /* transform these corner coordinates to output dataset grid coordinates */
00416 
00417          LOAD_FVEC3(vv , x1,y1,z1) ;
00418          vv = THD_3dmm_to_dicomm( iset , vv );        /* transpose from iset to Dicom */
00419          vv = AFNI_backward_warp_vector( warp , vv ); /* transform to +orig ???*/
00420          vv = THD_dicomm_to_3dmm( dset , vv );        /* and back from Dicom to dset  */
00421          UNLOAD_FVEC3(vv , xx1,yy1,zz1) ;
00422 
00423          LOAD_FVEC3(vv , x2,y2,z2) ;
00424          vv = THD_3dmm_to_dicomm( iset , vv ) ;
00425          vv = AFNI_backward_warp_vector( warp , vv ) ;
00426          vv = THD_dicomm_to_3dmm( dset , vv ) ;
00427          UNLOAD_FVEC3(vv , xx2,yy2,zz2) ;
00428 
00429          /* [xx1,xx2] X [yy1,yy2] X [zz1,zz2] is now in coordinates of output dataset */
00430 
00431          /* compute indices into output dataset voxel (keeping fractions) */
00432 
00433          f1 = (xx1-xorgout)/dxout + 0.49999 ; f2 = (xx2-xorgout)/dxout + 0.49999 ;
00434          if( f1 > f2 ){ tx = f1 ; f1 = f2 ; f2 = tx ; }
00435          if( f1 >= nxout || f2 <= 0.0 ) continue ;
00436          if( f1 < 0.0 ) f1 = 0.0 ;  if( f2 >= nxout ) f2 = nxout - 0.001 ;
00437 
00438          g1 = (yy1-yorgout)/dyout + 0.49999 ; g2 = (yy2-yorgout)/dyout + 0.49999 ;
00439          if( g1 > g2 ){ ty = g1 ; g1 = g2 ; g2 = ty ; }
00440          if( g1 >= nyout || g2 <= 0.0 ) continue ;
00441          if( g1 < 0.0 ) g1 = 0.0 ;  if( g2 >= nyout ) g2 = nyout - 0.001 ;
00442 
00443          h1 = (zz1-zorgout)/dzout + 0.49999 ; h2 = (zz2-zorgout)/dzout + 0.49999 ;
00444          if( h1 > h2 ){ tz = h1 ; h1 = h2 ; h2 = tz ; }
00445          if( h1 >= nzout || h2 <= 0.0 ) continue ;
00446          if( h1 < 0.0 ) h1 = 0.0 ;  if( h2 >= nzout ) h2 = nzout - 0.001 ;
00447 
00448          /* input voxel covers voxels [f1,f2] X [g1,g2] X [h1,h2] in the output */
00449 
00450          /* For example, [6.3,7.2] X [9.3,9.6] X [11.7,13.4], which must be     */
00451          /* distributed into these voxels:                                      */
00452          /*  (6,9,11), (7,9,11), (6,9,12), (7,9,12), (6,9,13), and (7,9,13)     */
00453 
00454          for( f=f1 ; f < f2 ; f = ip ){
00455             i = (int) f ; ip = i+1 ; tx = MIN(ip,f2) ; sx = tx - f ;
00456             for( g=g1 ; g < g2 ; g = jp ){
00457                j = (int) g ; jp = j+1 ; ty = MIN(jp,g2) ; sy = ty - g ;
00458                for( h=h1 ; h < h2 ; h = kp ){
00459                   k = (int) h ; kp = k+1 ; tz = MIN(kp,h2) ; sz = tz - h ;
00460                   voxout[ i + j*nxout + k * nxyout ] += sx * sy * sz ;
00461 
00462                }
00463             }
00464          }
00465 
00466       }}} /* end of loops over voxels */
00467 
00468       /* at this point, voxout[i] contains the fraction that output voxel [i]
00469          is occupied by input voxels that were allowed to vote (or were nonzero) */
00470 
00471       /* if not voting, we are done; otherwise, we must tally the count so far */
00472 
00473       if( do_vote ){
00474          for( i=0 ; i < nvoxout ; i++ ){
00475             if( voxout[i] > vote_best[i] && voxout[i] >= fclip ){  /* wins */
00476                vote_best[i] = voxout[i] ;
00477                     if( vtype == MRI_byte  ) vote_bout[i] = (byte)  voter ;
00478                else if( vtype == MRI_short ) vote_sout[i] = (short) voter ;
00479             }
00480          }
00481       }
00482 
00483       /* if we are voting, loop back for the next ballot */
00484 
00485    } while( do_vote && ++ivote < nvote_val ) ; /* end of voting loop */
00486 
00487    /*- throw out some trash -*/
00488 
00489    DSET_delete(iset) ;
00490    if( wset != NULL ) DSET_delete(wset) ;
00491 
00492    if( do_vote ){
00493       free(voxout) ; free(vote_best) ; free(vote_val) ;
00494       if( vote_print ) fprintf(stderr,"\n") ;
00495    }
00496 
00497    /*-- if not voting, compute the output brick --*/
00498 
00499    if( ! do_vote ){
00500       sin = (short *) malloc( sizeof(short) * nvoxout ) ;
00501       if( sin == NULL ){
00502          fprintf(stderr,"** Can't malloc output brick!\n") ; exit(1) ;
00503       }
00504 
00505       for( i=ijk=0 ; i < nvoxout ; i++ ){
00506          sin[i] = (short) (10000.0 * voxout[i] + 0.49999) ;
00507               if( sin[i] < sclip ) sin[i] = 0 ;
00508          else if( sin[i] > 10000 ) sin[i] = 10000 ;
00509 
00510          if( sin[i] != 0 ) ijk++ ;
00511       }
00512 
00513       free(voxout) ;
00514       EDIT_substitute_brick( dset , 0 , MRI_short , sin ) ;
00515 
00516       fprintf(stderr,"-- Writing %d nonzero mask voxels to dataset %s\n",
00517              ijk , DSET_BRIKNAME(dset) ) ;
00518 
00519    /*-- if voting, output brick will be in vote_bout or vote_sout --*/
00520 
00521    } else {
00522       if( vtype == MRI_byte ){
00523          EDIT_substitute_brick( dset , 0 , MRI_byte , vote_bout ) ;
00524          for( i=ijk=0 ; i < nvoxout ; i++ )
00525             if( vote_bout[i] != 0 ) ijk++ ;
00526       } else {
00527          EDIT_substitute_brick( dset , 0 , MRI_short , vote_sout ) ;
00528          for( i=ijk=0 ; i < nvoxout ; i++ )
00529             if( vote_sout[i] != 0 ) ijk++ ;
00530       }
00531 
00532       fprintf(stderr,"-- Writing %d nonzero voted voxels to dataset %s\n",
00533              ijk , DSET_BRIKNAME(dset) ) ;
00534    }
00535 
00536    DSET_write(dset) ;
00537    exit(0) ;
00538 }
00539 
00540 /*------------------------------------------------------------------------
00541    Backward transform a vector following a warp
00542    - lifted from afni.c
00543    - note that a NULL warp is equivalent to the identity
00544 --------------------------------------------------------------------------*/
00545 
00546 THD_fvec3 AFNI_backward_warp_vector( THD_warp * warp , THD_fvec3 old_fv )
00547 {
00548    THD_fvec3 new_fv ;
00549 
00550    if( warp == NULL ) return old_fv ;
00551 
00552    switch( warp->type ){
00553 
00554       default: new_fv = old_fv ; break ;
00555 
00556       case WARP_TALAIRACH_12_TYPE:{
00557          THD_linear_mapping map ;
00558          int iw ;
00559 
00560          /* test if input is in bot..top of each defined map */
00561 
00562          for( iw=0 ; iw < 12 ; iw++ ){
00563             map = warp->tal_12.warp[iw] ;
00564 
00565             if( old_fv.xyz[0] >= map.bot.xyz[0] &&
00566                 old_fv.xyz[1] >= map.bot.xyz[1] &&
00567                 old_fv.xyz[2] >= map.bot.xyz[2] &&
00568                 old_fv.xyz[0] <= map.top.xyz[0] &&
00569                 old_fv.xyz[1] <= map.top.xyz[1] &&
00570                 old_fv.xyz[2] <= map.top.xyz[2]   ) break ;  /* leave loop */
00571          }
00572          new_fv = MATVEC_SUB(map.mbac,old_fv,map.svec) ;
00573       }
00574       break ;
00575 
00576       case WARP_AFFINE_TYPE:{
00577          THD_linear_mapping map = warp->rig_bod.warp ;
00578          new_fv = MATVEC_SUB(map.mbac,old_fv,map.svec) ;
00579       }
00580       break ;
00581 
00582    }
00583    return new_fv ;
00584 }