do{ int nn ;                              \
      if( (cl)->num_clu == (cl)->num_all ){ \
         (cl)->num_all += INC_CLUSTER ; nn = (cl)->num_all ;           \
         (cl)->clar = (MCW_cluster **) XtRealloc( (char *)(cl)->clar , \
                                                 sizeof(MCW_cluster *) * nn ) ;\
      } \
      (cl)->clar[((cl)->num_clu)++] = (cc) ; break ; } while(0)

do{ int nn ;                                              \
      if( (cc)->num_pt == (cc)->num_all ){                  \
         (cc)->num_all = 1.25*(cc)->num_all + INC_CLUSTER ; \
         nn = (cc)->num_all ;                               \
         (cc)->i=(short *)   XtRealloc((char *)(cc)->i,sizeof(short)*nn  );\
         (cc)->j=(short *)   XtRealloc((char *)(cc)->j,sizeof(short)*nn  );\
         (cc)->k=(short *)   XtRealloc((char *)(cc)->k,sizeof(short)*nn  );\
         (cc)->mag=(float *) XtRealloc((char *)(cc)->mag,sizeof(float)*nn);\
         DBMALL(nn) ; }                      \
      nn = (cc)->num_pt ; ((cc)->num_pt)++ ; \
      (cc)->i[nn] = (ii) ; (cc)->j[nn] = (jj) ; (cc)->k[nn] = (kk) ; \
      (cc)->mag[nn] = (m) ; break ; } while(0)

do{ int ii ; if( cl != NULL ){                    \
         for( ii=0 ; ii < (cl)->num_clu ; ii++ )    \
            KILL_CLUSTER( (cl)->clar[ii] ) ;        \
         myXtFree((cl)->clar) ; (cl) = NULL ; \
      } break ; } while(0)

do{ if( ISVALID_DSET(nds) && ISVALID_DSET(ods) )              \
         (nds)->anat_parent_idcode = (ods)->anat_parent_idcode ; \
    } while(0)

do{ float sqq[6] ;                                           \
     if( ISVALID_DSET(ds) &&                                    \
         (iv) >= 0 && (iv) < DSET_NVALS(ds) &&                    \
         (ft) >= 0 && (ft) <= LAST_FUNC_TYPE   ){                   \
        sqq[0] = (ft) ; sqq[1] = FUNC_need_stat_aux[ft] ;             \
        sqq[2] = (a) ; sqq[3] = (b) ; sqq[4] = (c) ; sqq[5] = (d) ;     \
        EDIT_dset_items( (ds),ADN_brick_stataux_one+(iv),sqq,ADN_none ) ; \
     } } while(0)

do{ if( ISVALID_DSET(nds) )                             \
         ZERO_IDCODE((nds)->anat_parent_idcode); } while(0)

( (cl) = XtNew(MCW_cluster_array) , \
    (cl)->num_clu = (cl)->num_all = 0 , (cl)->clar = NULL )

( (cc) = XtNew(MCW_cluster) ,        \
    (cc)->num_pt = (cc)->num_all = 0 , \
    (cc)->i = NULL , (cc)->j = NULL , (cc)->k = NULL ,(cc)->mag = NULL )

( (edopt)->thtoin        = 0   , \
        (edopt)->noneg         = 0   , \
        (edopt)->abss          = 0   , \
        (edopt)->clip_bot      = 0.0 , \
        (edopt)->clip_top      = 0.0 , \
        (edopt)->thresh        = 0.0 , \
        (edopt)->clust_rmm     = -1.0, \
        (edopt)->clust_vmul    = 0.0 , \
        (edopt)->edit_clust    = 0   , \
        (edopt)->erode_pv      = 0.0 , \
        (edopt)->dilate        = 0   , \
        (edopt)->filter_rmm    = 0.0 , \
        (edopt)->filter_opt    = 0   , \
        (edopt)->thrfilter_rmm = 0.0 , \
        (edopt)->thrfilter_opt = 0   , \
        (edopt)->blur          = 0.0 , \
        (edopt)->thrblur       = 0.0 , \
        (edopt)->scale         = 0   , \
        (edopt)->mult          = 0.0 , \
        (edopt)->do_zvol       = 0   , \
        (edopt)->clip_unscaled = 0   , \
        (edopt)->iv_fim        = -1  , \
        (edopt)->iv_thr        = -1  , \
        (edopt)->zscore        = 0   , \
        (edopt)->verbose       = 0   , \
        (edopt)->nfmask        = 0   , \
        (edopt)->fmask         = NULL, \
        (edopt)->fexpr         = NULL, \
        (edopt)->fake_dxyz     = 0   , \
       0 )

do{ if( cc != NULL ){        \
         myXtFree((cc)->i) ;   \
         myXtFree((cc)->j) ;   \
         myXtFree((cc)->k) ;   \
         myXtFree((cc)->mag) ; \
         myXtFree((cc)) ;      \
         (cc) = NULL ;         \
      } break ; } while(0)

if( (name) != NULL && (name)->num_clu > 1 ){                             \
      int iic , jjc , sss ; MCW_cluster * ct ;                              \
      for( iic=0 ; iic < (name)->num_clu ; iic++ ){                         \
         sss = 0 ;                                                          \
         for( jjc=1 ; jjc < (name)->num_clu ; jjc++ ){                      \
            if( !CLUST_ORDERED( (name)->clar[jjc-1] , (name)->clar[jjc] ) ) \
               CLUST_SWAP( (name)->clar[jjc-1] , (name)->clar[jjc] ) ;      \
   } if( !sss ) break ; }}

00116 {
00117    int    ttt = typ ;
00118    float  fff = fac ;
00119    void * bbb = br ;
00120 
00121    EDIT_add_bricklist( dset , 1 , &ttt , &fff , &bbb ) ;
00122    return ;
00123 }

00026 {
00027    int ibr , typ , nx,ny,nz , nvals,new_nvals ;
00028    THD_datablock * dblk ;
00029    MRI_IMAGE * qim ;
00030    char str[32] ;
00031 
00032 ENTRY("EDIT_add_bricklist") ;
00033 
00034    /**-- Sanity Checks --**/
00035 
00036    if( ! ISVALID_3DIM_DATASET(dset) || nbr <= 0 )       EXRETURN; /* error! */
00037    if( dset->dblk->brick == NULL )                      EXRETURN; /* error! */
00038    if( dset->dblk->malloc_type != DATABLOCK_MEM_MALLOC )EXRETURN; /* error! */
00039 
00040    dblk  = dset->dblk ;
00041    nvals = dblk->nvals ;
00042    nx    = dblk->diskptr->dimsizes[0] ;
00043    ny    = dblk->diskptr->dimsizes[1] ;
00044    nz    = dblk->diskptr->dimsizes[2] ;
00045 
00046    /**-- reallocate the brick control information --**/
00047 
00048    new_nvals = nvals + nbr ;
00049    dblk->brick_bytes = (int *) XtRealloc( (char *) dblk->brick_bytes ,
00050                                           sizeof(int) * new_nvals ) ;
00051 
00052    dblk->brick_fac = (float *) XtRealloc( (char *) dblk->brick_fac ,
00053                                           sizeof(float) * new_nvals ) ;
00054 
00055    dblk->nvals = dblk->diskptr->nvals = new_nvals ;
00056 
00057    /** allocate new sub-brick images **/
00058 
00059    for( ibr=0 ; ibr < nbr ; ibr++ ){
00060       typ = (tbr != NULL ) ? tbr[ibr] : MRI_short ;
00061       qim = mri_new_vol_empty( nx,ny,nz , typ ) ;  /* image with no data */
00062 
00063       if( sbr != NULL && sbr[ibr] != NULL )        /* attach data to image */
00064          mri_fix_data_pointer( sbr[ibr] , qim ) ;
00065 
00066       ADDTO_IMARR( dblk->brick , qim ) ;           /* attach image to dset */
00067 
00068       dblk->brick_fac[nvals+ibr]   = (fbr != NULL) ? fbr[ibr] : 0.0 ;
00069       dblk->brick_bytes[nvals+ibr] = qim->pixel_size * qim->nvox ;
00070       dblk->total_bytes           += dblk->brick_bytes[ibr] ;
00071    }
00072 
00073    /** allocate new sub-brick auxiliary data: labels **/
00074 
00075    if( dblk->brick_lab == NULL )
00076       THD_init_datablock_labels( dblk ) ;
00077    else
00078       dblk->brick_lab = (char **) XtRealloc( (char *) dblk->brick_lab ,
00079                                              sizeof(char *) * new_nvals ) ;
00080    for( ibr=0 ; ibr < nbr ; ibr++ ){
00081       sprintf( str , "#%d" , nvals+ibr ) ;
00082       dblk->brick_lab[nvals+ibr] = NULL ;
00083       THD_store_datablock_label( dblk , nvals+ibr , str ) ;
00084    }
00085 
00086    /** keywords **/
00087 
00088    if( dblk->brick_keywords == NULL )
00089       THD_init_datablock_keywords( dblk ) ;
00090    else
00091       dblk->brick_keywords = (char **) XtRealloc( (char *) dblk->brick_keywords ,
00092                                                   sizeof(char *) * new_nvals ) ;
00093    for( ibr=0 ; ibr < nbr ; ibr++ ){
00094       dblk->brick_keywords[nvals+ibr] = NULL ;
00095       THD_store_datablock_keywords( dblk , nvals+ibr , NULL ) ;
00096    }
00097 
00098    /** stataux **/
00099 
00100    if( dblk->brick_statcode != NULL ){
00101       dblk->brick_statcode = (int *) XtRealloc( (char *) dblk->brick_statcode ,
00102                                                 sizeof(int) * new_nvals        ) ;
00103       dblk->brick_stataux  = (float **) XtRealloc( (char *) dblk->brick_stataux ,
00104                                                    sizeof(float *) * new_nvals ) ;
00105 
00106       for( ibr=0 ; ibr < nbr ; ibr++ ){
00107          dblk->brick_statcode[nvals+ibr] = 0 ;
00108          dblk->brick_stataux[nvals+ibr]  = NULL ;
00109       }
00110    }
00111 
00112   EXRETURN;
00113 }

00393 {
00394    MCW_cluster * mask ;
00395    int i, j, k , ij , ii ;
00396    int jk,jkadd , nxadd,nyadd,nzadd , nxyz_add , mnum ;
00397    float * ffim ;
00398    int * madd ;
00399    float fac , sum ;
00400 
00401 ENTRY("EDIT_aver_fvol") ;
00402 
00403    /*--- Make a cluster that is a mask of points closer than rmm ---*/
00404 
00405    mask = MCW_build_mask(nx,ny,nz, dx,dy,dz, rmm) ;
00406    if( mask == NULL || mask->num_pt < 2 ){
00407       fprintf(stderr,"Warning: EDIT_aver_volume has no effect.\n") ;
00408       EXRETURN ;
00409    }
00410    mnum = mask->num_pt ;
00411 
00412    /*--- Allocate workspaces ---*/
00413 
00414 #if 1
00415    nxadd = nyadd = nzadd = 1 ;
00416    for( ii=0 ; ii < mnum ; ii++ ){
00417       i = abs((int)mask->i[ii]) ; nxadd = MAX(i,nxadd) ;
00418       j = abs((int)mask->j[ii]) ; nyadd = MAX(j,nyadd) ;
00419       k = abs((int)mask->k[ii]) ; nzadd = MAX(k,nzadd) ;
00420    }
00421 #else
00422    nxadd    = (int)(rmm/dx) ;
00423    nyadd    = (int)(rmm/dy) ;
00424    nzadd    = (int)(rmm/dz) ;
00425 #endif
00426 
00427    nxyz_add = (nx+2*nxadd) * (ny+2*nyadd) * (nz+2*nzadd) ;
00428 
00429    ffim = (float *) malloc(sizeof(float) * nxyz_add) ;
00430    if(ffim == NULL){
00431       fprintf(stderr,"*** EDIT_aver_volume can't malloc workspace!\n") ;
00432       fprintf(stderr,"nx=%d ny=%d nz=%d nxadd=%d nyadd=%d nzadd=%d\n",
00433               nx,ny,nz , nxadd,nyadd,nzadd ) ;
00434       EXIT(1) ;
00435    }
00436    for( i=0 ; i < nxyz_add ; i++ ) ffim[i] = 0.0 ;
00437 
00438    madd = (int *) malloc( sizeof(int) * (mnum+1) ) ;
00439    if( madd == NULL ){
00440       fprintf(stderr,"*** EDIT_aver_volume can't malloc workspace!\n") ;
00441       EXIT(1) ;
00442    }
00443    madd[0] = 0 ;
00444    for( ii=0 ; ii < mnum ; ii++ ){
00445       madd[ii+1] = mask->i[ii] +
00446                    mask->j[ii] * (nx+2*nxadd) +
00447                    mask->k[ii] * ((nx+2*nxadd)*(ny+2*nyadd)) ;
00448    }
00449    mnum++ ; fac = 1.0 / mnum ;
00450 
00451    KILL_CLUSTER(mask) ;
00452 
00453    /*-- copy data into workspace --*/
00454 
00455    for( k=0 ; k < nz ; k++ ){
00456       for( j=0 ; j < ny ; j++ ){
00457          jkadd = j * (nx+2*nxadd) + k * ((nx+2*nxadd)*(ny+2*nyadd)) ;
00458          jk    = j * nx + k * (nx * ny) ;
00459          for( i=0 ; i < nx ; i++ ) ffim[i+jkadd] = fim[i+jk] ;
00460       }
00461    }
00462 
00463    /*-- average data from workspace back into original array --*/
00464 
00465    for( k=0 ; k < nz ; k++ ){
00466       for( j=0 ; j < ny ; j++ ){
00467          jkadd = j * (nx+2*nxadd) + k * ((nx+2*nxadd)*(ny+2*nyadd)) ;
00468          jk    = j * nx + k * (nx * ny) ;
00469          for( i=0 ; i < nx ; i++ ){
00470             sum = 0.0 ; ij = i+jkadd ;
00471             for( ii=0 ; ii < mnum ; ii++ ) sum += ffim[ij+madd[ii]] ;
00472             fim[i+jk] = fac * sum ;
00473          }
00474       }
00475    }
00476 
00477    free(ffim); free(madd);
00478    EXRETURN;
00479 }

00027 {
00028   EDIT_blur_volume_3d (nx, ny, nz, dx, dy, dz, ftype, vfim,
00029                        sigma, sigma, sigma);
00030 }

00044 {
00045    int jj,kk , nxy , base,nby2 ;
00046    float  dk , aa , k , fac ;
00047    register int ii , nup ;
00048 
00049    static int cx_size  = 0 ;     /* workspaces */
00050    static int gg_size  = 0 ;
00051    static complex * cx = NULL ;
00052    static float   * gg = NULL ;
00053 
00054    byte *    bfim = NULL ;
00055    short *   sfim = NULL ;
00056    float *   ffim = NULL ;
00057    complex * cfim = NULL ;
00058 
00059    float fbot,ftop ;     /* 10 Jan 2003 */
00060    int nxyz ;
00061 
00062    /*** initialize ***/
00063 
00064 ENTRY("EDIT_blur_volume_3d") ;
00065 
00066    if( vfim == NULL ||
00067        sigmax <= 0.0 || sigmay <= 0.0 || sigmaz <= 0.0)  EXRETURN ;
00068 
00069    switch( ftype ){
00070       default: EXRETURN ;
00071       case MRI_short:   sfim = (short *)   vfim ; break ;
00072       case MRI_float:   ffim = (float *)   vfim ; break ;
00073       case MRI_byte:    bfim = (byte *)    vfim ; break ;
00074       case MRI_complex: cfim = (complex *) vfim ; break ;
00075    }
00076    nxy = nx * ny ; nxyz = nxy * nz ;
00077 
00078    /*** 10 Jan 2003: find bot and top of data input */
00079 
00080    switch( ftype ){
00081      default:
00082        fbot = ftop = 0.0 ;
00083      break ;
00084 
00085      case MRI_short:
00086        fbot = ftop = sfim[0] ;
00087        for( ii=1 ; ii < nxyz ; ii++ )
00088               if( sfim[ii] < fbot ) fbot = sfim[ii] ;
00089          else if( sfim[ii] > ftop ) ftop = sfim[ii] ;
00090      break ;
00091 
00092      case MRI_float:
00093        fbot = ftop = ffim[0] ;
00094        for( ii=1 ; ii < nxyz ; ii++ )
00095               if( ffim[ii] < fbot ) fbot = ffim[ii] ;
00096          else if( ffim[ii] > ftop ) ftop = ffim[ii] ;
00097      break ;
00098 
00099      case MRI_byte:
00100        fbot = ftop = bfim[0] ;
00101        for( ii=1 ; ii < nxyz ; ii++ )
00102               if( bfim[ii] < fbot ) fbot = bfim[ii] ;
00103          else if( bfim[ii] > ftop ) ftop = bfim[ii] ;
00104      break ;
00105    }
00106 
00107    /*** do x-direction ***/
00108 
00109 STATUS("start x FFTs") ;
00110 
00111    aa  = sigmax * sigmax * 0.5 ;
00112    nup = nx + (int)(3.0 * sigmax / dx) ;      /* min FFT length */
00113 #if 0
00114    ii  = 2 ; while( ii < nup ){ ii *= 2 ; }  /* next power of 2 larger */
00115 #else
00116    ii = csfft_nextup_one35(nup) ;
00117 #endif
00118    nup = ii ; nby2 = nup / 2 ;
00119 
00120    GET_AS_BIG(cx,complex,nup) ; GET_AS_BIG(gg,float,nup) ;
00121 
00122    dk    = (2.0*PI) / (nup * dx) ;
00123    fac   = 1.0 / nup ;
00124    gg[0] = fac ;
00125    for( ii=1 ; ii<=nby2 ; ii++ ){ k=ii*dk; gg[nup-ii]=gg[ii]=fac*exp(-aa*k*k); }
00126 
00127    /** July 20: double up on FFTs **/
00128    /** Feb  09: extend to other data types besides shorts;
00129                 doubling up does not apply to complex data! **/
00130 
00131    switch( ftype ){
00132       case MRI_short:{
00133          register short * qfim ;
00134          for( kk=0 ; kk < nz ; kk++ ){
00135             for( jj=0 ; jj < ny ; jj+=2 ){
00136                base = jj*nx + kk*nxy ;
00137                qfim = sfim + base ;
00138                if( jj == ny-1 )
00139                   for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = 0.0 ; }
00140                else
00141                   for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = qfim[ii+nx] ; }
00142                for( ii=nx; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00143                csfft_cox( -1 , nup , cx ) ;
00144                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00145                csfft_cox(  1 , nup , cx ) ;
00146                if( jj == ny-1 )
00147                   for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; }
00148                else
00149                   for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; qfim[ii+nx] = cx[ii].i ; }
00150             }
00151          }
00152       }
00153       break ;
00154 
00155       case MRI_float:{
00156          register float * qfim ;
00157          for( kk=0 ; kk < nz ; kk++ ){
00158             for( jj=0 ; jj < ny ; jj+=2 ){
00159                base = jj*nx + kk*nxy ;
00160                qfim = ffim + base ;
00161                if( jj == ny-1 )
00162                   for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = 0.0 ; }
00163                else
00164                   for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = qfim[ii+nx] ; }
00165                for( ii=nx; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00166                csfft_cox( -1 , nup , cx ) ;
00167                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00168                csfft_cox(  1 , nup , cx ) ;
00169                if( jj == ny-1 )
00170                   for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; }
00171                else
00172                   for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; qfim[ii+nx] = cx[ii].i ; }
00173             }
00174          }
00175       }
00176       break ;
00177 
00178       case MRI_byte:{
00179          register byte * qfim ;
00180          for( kk=0 ; kk < nz ; kk++ ){
00181             for( jj=0 ; jj < ny ; jj+=2 ){
00182                base = jj*nx + kk*nxy ;
00183                qfim = bfim + base ;
00184                if( jj == ny-1 )
00185                   for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = 0.0 ; }
00186                else
00187                   for( ii=0 ; ii<nx ; ii++){ cx[ii].r = qfim[ii] ; cx[ii].i = qfim[ii+nx] ; }
00188                for( ii=nx; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00189                csfft_cox( -1 , nup , cx ) ;
00190                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00191                csfft_cox(  1 , nup , cx ) ;
00192                if( jj == ny-1 )
00193                   for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; }
00194                else
00195                   for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii].r ; qfim[ii+nx] = cx[ii].i ; }
00196             }
00197          }
00198       }
00199       break ;
00200 
00201       case MRI_complex:{
00202          register complex * qfim ;
00203          for( kk=0 ; kk < nz ; kk++ ){
00204             for( jj=0 ; jj < ny ; jj++ ){
00205                base = jj*nx + kk*nxy ;
00206                qfim = cfim + base ;
00207                for( ii=0 ; ii<nx ; ii++) { cx[ii] = qfim[ii] ; }
00208                for( ii=nx; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00209                csfft_cox( -1 , nup , cx ) ;
00210                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00211                csfft_cox(  1 , nup , cx ) ;
00212                for( ii=0 ; ii<nx ; ii++){ qfim[ii] = cx[ii] ; }
00213             }
00214          }
00215       }
00216       break ;
00217    }
00218 
00219    /*** do y-direction ***/
00220 
00221 STATUS("start y FFTs") ;
00222 
00223    aa  = sigmay * sigmay * 0.5 ;
00224    nup = ny + (int)(3.0 * sigmay / dy) ;      /* min FFT length */
00225 #if 0
00226    ii  = 2 ; while( ii < nup ){ ii *= 2 ; }  /* next power of 2 larger */
00227 #else
00228    ii = csfft_nextup_one35(nup) ;
00229 #endif
00230    nup = ii ; nby2 = nup / 2 ;
00231 
00232    GET_AS_BIG(cx,complex,nup) ; GET_AS_BIG(gg,float,nup) ;
00233 
00234    dk    = (2.0*PI) / (nup * dy) ;
00235    fac   = 1.0 / nup ;
00236    gg[0] = fac ;
00237    for( ii=1 ; ii<=nby2 ; ii++ ){ k=ii*dk; gg[nup-ii]=gg[ii]=fac*exp(-aa*k*k); }
00238 
00239    switch( ftype ){
00240       case MRI_short:{
00241          register short * qfim ;
00242          for( kk=0 ; kk < nz ; kk++ ){
00243             for( jj=0 ; jj < nx ; jj+=2 ){
00244                base = jj + kk*nxy ;
00245                qfim = sfim + base ;
00246                if( jj == nx-1 )
00247                   for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = 0.0 ; }
00248                else
00249                   for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = qfim[ii*nx+1] ; }
00250                for( ii=ny; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00251                csfft_cox( -1 , nup , cx ) ;
00252                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00253                csfft_cox(  1 , nup , cx ) ;
00254                if( jj == nx-1 )
00255                   for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; }
00256                else
00257                   for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; qfim[ii*nx+1] = cx[ii].i ; }
00258             }
00259          }
00260       }
00261       break ;
00262 
00263       case MRI_byte:{
00264          register byte * qfim ;
00265          for( kk=0 ; kk < nz ; kk++ ){
00266             for( jj=0 ; jj < nx ; jj+=2 ){
00267                base = jj + kk*nxy ;
00268                qfim = bfim + base ;
00269                if( jj == nx-1 )
00270                   for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = 0.0 ; }
00271                else
00272                   for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = qfim[ii*nx+1] ; }
00273                for( ii=ny; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00274                csfft_cox( -1 , nup , cx ) ;
00275                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00276                csfft_cox(  1 , nup , cx ) ;
00277                if( jj == nx-1 )
00278                   for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; }
00279                else
00280                   for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; qfim[ii*nx+1] = cx[ii].i ; }
00281             }
00282          }
00283       }
00284       break ;
00285 
00286       case MRI_float:{
00287          register float * qfim ;
00288          for( kk=0 ; kk < nz ; kk++ ){
00289             for( jj=0 ; jj < nx ; jj+=2 ){
00290                base = jj + kk*nxy ;
00291                qfim = ffim + base ;
00292                if( jj == nx-1 )
00293                   for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = 0.0 ; }
00294                else
00295                   for( ii=0 ; ii<ny ; ii++){ cx[ii].r = qfim[ii*nx] ; cx[ii].i = qfim[ii*nx+1] ; }
00296                for( ii=ny; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00297                csfft_cox( -1 , nup , cx ) ;
00298                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00299                csfft_cox(  1 , nup , cx ) ;
00300                if( jj == nx-1 )
00301                   for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; }
00302                else
00303                   for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii].r ; qfim[ii*nx+1] = cx[ii].i ; }
00304             }
00305          }
00306       }
00307       break ;
00308 
00309       case MRI_complex:{
00310          register complex * qfim ;
00311          for( kk=0 ; kk < nz ; kk++ ){
00312             for( jj=0 ; jj < nx ; jj++ ){
00313                base = jj + kk*nxy ;
00314                qfim = cfim + base ;
00315                for( ii=0 ; ii<ny ; ii++){ cx[ii] = qfim[ii*nx] ; }
00316                for( ii=ny; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00317                csfft_cox( -1 , nup , cx ) ;
00318                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00319                csfft_cox(  1 , nup , cx ) ;
00320                for( ii=0 ; ii<ny ; ii++){ qfim[ii*nx] = cx[ii] ; }
00321             }
00322          }
00323       }
00324       break ;
00325    }
00326 
00327    /*** do z-direction ***/
00328 
00329 STATUS("start z FFTs") ;
00330 
00331    aa  = sigmaz * sigmaz * 0.5 ;
00332    nup = nz + (int)(3.0 * sigmaz / dz) ;      /* min FFT length */
00333 #if 0
00334    ii  = 2 ; while( ii < nup ){ ii *= 2 ; }  /* next power of 2 larger */
00335 #else
00336    ii = csfft_nextup_one35(nup) ;
00337 #endif
00338    nup = ii ; nby2 = nup / 2 ;
00339 
00340    GET_AS_BIG(cx,complex,nup) ; GET_AS_BIG(gg,float,nup) ;
00341 
00342    dk    = (2.0*PI) / (nup * dz) ;
00343    fac   = 1.0 / nup ;
00344    gg[0] = fac ;
00345    for( ii=1 ; ii<=nby2 ; ii++ ){ k=ii*dk; gg[nup-ii]=gg[ii]=fac*exp(-aa*k*k); }
00346 
00347    switch( ftype ){
00348       case MRI_short:{
00349          register short * qfim ;
00350          for( kk=0 ; kk < ny ; kk++ ){
00351             for( jj=0 ; jj < nx ; jj+=2 ){
00352                base = jj + kk*nx ;
00353                qfim = sfim + base ;
00354                if( jj == nx-1 )
00355                   for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = 0.0 ; }
00356                else
00357                   for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = qfim[ii*nxy+1] ; }
00358                for( ii=nz; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00359                csfft_cox( -1 , nup , cx ) ;
00360                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00361                csfft_cox(  1 , nup , cx ) ;
00362                if( jj == nx-1 )
00363                   for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; }
00364                else
00365                   for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; qfim[ii*nxy+1] = cx[ii].i ; }
00366             }
00367          }
00368       }
00369       break ;
00370 
00371       case MRI_float:{
00372          register float * qfim ;
00373          for( kk=0 ; kk < ny ; kk++ ){
00374             for( jj=0 ; jj < nx ; jj+=2 ){
00375                base = jj + kk*nx ;
00376                qfim = ffim + base ;
00377                if( jj == nx-1 )
00378                   for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = 0.0 ; }
00379                else
00380                   for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = qfim[ii*nxy+1] ; }
00381                for( ii=nz; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00382                csfft_cox( -1 , nup , cx ) ;
00383                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00384                csfft_cox(  1 , nup , cx ) ;
00385                if( jj == nx-1 )
00386                   for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; }
00387                else
00388                   for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; qfim[ii*nxy+1] = cx[ii].i ; }
00389             }
00390          }
00391       }
00392       break ;
00393 
00394       case MRI_byte:{
00395          register byte * qfim ;
00396          for( kk=0 ; kk < ny ; kk++ ){
00397             for( jj=0 ; jj < nx ; jj+=2 ){
00398                base = jj + kk*nx ;
00399                qfim = bfim + base ;
00400                if( jj == nx-1 )
00401                   for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = 0.0 ; }
00402                else
00403                   for( ii=0 ; ii<nz ; ii++){ cx[ii].r = qfim[ii*nxy] ; cx[ii].i = qfim[ii*nxy+1] ; }
00404                for( ii=nz; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00405                csfft_cox( -1 , nup , cx ) ;
00406                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00407                csfft_cox(  1 , nup , cx ) ;
00408                if( jj == nx-1 )
00409                   for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; }
00410                else
00411                   for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii].r ; qfim[ii*nxy+1] = cx[ii].i ; }
00412             }
00413          }
00414       }
00415       break ;
00416 
00417       case MRI_complex:{
00418          register complex * qfim ;
00419          for( kk=0 ; kk < ny ; kk++ ){
00420             for( jj=0 ; jj < nx ; jj++ ){
00421                base = jj + kk*nx ;
00422                qfim = cfim + base ;
00423                for( ii=0 ; ii<nz ; ii++){ cx[ii] = qfim[ii*nxy] ; }
00424                for( ii=nz; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00425                csfft_cox( -1 , nup , cx ) ;
00426                for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00427                csfft_cox(  1 , nup , cx ) ;
00428                for( ii=0 ; ii<nz ; ii++){ qfim[ii*nxy] = cx[ii] ; }
00429             }
00430          }
00431       }
00432       break ;
00433    }
00434 
00435    /*** 10 Jan 2003: clip data to bot and top found above ***/
00436    /***              to minimize Gibbs ringing artifacts  ***/
00437 
00438    switch( ftype ){
00439 
00440      case MRI_short:
00441        for( ii=0 ; ii < nxyz ; ii++ )
00442               if( sfim[ii] < fbot ) sfim[ii] = fbot ;
00443          else if( sfim[ii] > ftop ) sfim[ii] = ftop ;
00444      break ;
00445 
00446      case MRI_float:
00447        for( ii=0 ; ii < nxyz ; ii++ )
00448               if( ffim[ii] < fbot ) ffim[ii] = fbot ;
00449          else if( ffim[ii] > ftop ) ffim[ii] = ftop ;
00450      break ;
00451 
00452      case MRI_byte:
00453        for( ii=0 ; ii < nxyz ; ii++ )
00454               if( bfim[ii] < fbot ) bfim[ii] = fbot ;
00455          else if( bfim[ii] > ftop ) bfim[ii] = ftop ;
00456      break ;
00457    }
00458 
00459    /*** done! ***/
00460 
00461    EXRETURN ;
00462 }

00019 {
00020    float val ;
00021    int  ival , nopt_in=nopt ;
00022 
00023 ENTRY("EDIT_check_argv") ;
00024 
00025    /**** -1clip val ****/
00026 
00027    if( strncmp(argv[nopt],"-1clip" ,6) == 0 ||
00028        strncmp(argv[nopt],"-1uclip",6) == 0 ){
00029       nopt++ ;
00030       if( nopt >= argc ){
00031          fprintf(stderr,"no argument after %s?\n",argv[nopt-1]) ;
00032          EXIT(1);
00033       }
00034       edopt->clip_top = strtod( argv[nopt++] , NULL ) ;
00035       if( edopt->clip_top <= 0 ){
00036          fprintf(stderr,"illegal value after %s!\n",argv[nopt-2]) ;
00037          EXIT(1) ;
00038       }
00039       edopt->clip_bot = -edopt->clip_top ;
00040       edopt->clip_unscaled = (strncmp(argv[nopt-2],"-1uclip",6) == 0) ;
00041       CHECK_DONE ;
00042    }
00043 
00044    /**** -2clip val1 val2 ****/
00045 
00046    if( strncmp(argv[nopt],"-2clip" ,6) == 0 ||
00047        strncmp(argv[nopt],"-2uclip",6) == 0   ){
00048       nopt++ ;
00049       if( nopt+1 >= argc ){
00050          fprintf(stderr,"no arguments after %s?\n",argv[nopt-1]) ;
00051          EXIT(1) ;
00052       }
00053       edopt->clip_bot = strtod( argv[nopt++] , NULL ) ;  /* bot */
00054       edopt->clip_top = strtod( argv[nopt++] , NULL ) ;  /* top */
00055 
00056       if( edopt->clip_bot >= edopt->clip_top ){
00057          fprintf(stderr,
00058                  "*** %s %f %f is illegal:\n"
00059                  "*** first value must be less than second value!\n",
00060                  argv[nopt-3] , edopt->clip_bot , edopt->clip_top ) ;
00061          EXIT(1) ;
00062       }
00063       edopt->clip_unscaled = (strncmp(argv[nopt-3],"-2uclip",6) == 0) ;
00064       CHECK_DONE ;
00065    }
00066 
00067    /**** -1thtoin ****/
00068 
00069    if( strncmp(argv[nopt],"-1thtoin",6) == 0 ){
00070       edopt->thtoin = 1 ;
00071       nopt++ ; CHECK_DONE ;
00072    }
00073 
00074    /**** -2thtoin ****/
00075 
00076    if( strncmp(argv[nopt],"-2thtoin",6) == 0 ){
00077       edopt->thtoin = 2 ;
00078       nopt++ ; CHECK_DONE ;
00079    }
00080 
00081    /**** -1zscore (17 Sep 1998) ****/
00082 
00083    if( strncmp(argv[nopt],"-1zscore",6) == 0 ){
00084       edopt->zscore = 1 ;
00085       nopt++ ; CHECK_DONE ;
00086    }
00087 
00088    /**** -dxyz=1 (11 Sep 2000) ****/
00089 
00090    if( strcmp(argv[nopt],"-dxyz=1") == 0 ){
00091       edopt->fake_dxyz = 1 ;
00092       nopt++ ; CHECK_DONE ;
00093    }
00094 
00095    /**** -1noneg ****/
00096 
00097    if( strncmp(argv[nopt],"-1noneg",6) == 0 ){
00098       edopt->noneg = 1 ;
00099       nopt++ ; CHECK_DONE ;
00100    }
00101 
00102    /**** -1abs ****/
00103 
00104    if( strncmp(argv[nopt],"-1abs",6) == 0 ){
00105       edopt->abss = 1 ;
00106       nopt++ ; CHECK_DONE ;
00107    }
00108 
00109    /**** -1thresh thr ****/
00110 
00111    if( strncmp(argv[nopt],"-1thresh",6) == 0 ){
00112       nopt++ ;
00113       if( nopt >= argc ){
00114          fprintf(stderr,"no argument after -1thresh!\n") ;
00115          EXIT(1);
00116       }
00117       val = strtod( argv[nopt++] , NULL ) ;
00118       if( val < 0.0 ){
00119          fprintf(stderr,"illegal value after -1thresh!\n") ;
00120          EXIT(1) ;
00121       }
00122       edopt->thresh = val ;
00123       CHECK_DONE ;
00124    }
00125 
00126    /**** -1clust rmm vmul ****/
00127 
00128    if( strncmp(argv[nopt],"-1clust",12) == 0 ){
00129       nopt++ ;
00130       if( nopt+1 >= argc ){
00131          fprintf(stderr,"need 2 arguments after -1clust!\n") ;
00132          EXIT(1) ;
00133       }
00134       edopt->edit_clust = ECFLAG_SAME;
00135       edopt->clust_rmm  = strtod( argv[nopt++] , NULL ) ;
00136       edopt->clust_vmul = strtod( argv[nopt++] , NULL ) ;
00137       if( edopt->clust_rmm < 0 ){
00138          fprintf(stderr,"illegal value after -1clust\n") ;
00139          EXIT(1) ;
00140       }
00141       CHECK_DONE ;
00142    }
00143 
00144 
00145    /**** -1clust_mean rmm vmul ****/   /* 10 Sept 1996 */
00146 
00147    if( strncmp(argv[nopt],"-1clust_mean",12) == 0 ){
00148       nopt++ ;
00149       if( nopt+1 >= argc ){
00150          fprintf(stderr,"need 2 arguments after -1clust_mean!\n") ;
00151          EXIT(1) ;
00152       }
00153       edopt->edit_clust = ECFLAG_MEAN;
00154       edopt->clust_rmm  = strtod( argv[nopt++] , NULL ) ;
00155       edopt->clust_vmul = strtod( argv[nopt++] , NULL ) ;
00156       if( edopt->clust_rmm < 0 ){
00157          fprintf(stderr,"illegal value after -1clust_mean\n") ;
00158          EXIT(1) ;
00159       }
00160       CHECK_DONE ;
00161    }
00162 
00163    /**** -1clust_max rmm vmul ****/   /* 10 Sept 1996 */
00164 
00165    if( strncmp(argv[nopt],"-1clust_max",12) == 0 ){
00166       nopt++ ;
00167       if( nopt+1 >= argc ){
00168          fprintf(stderr,"need 2 arguments after -1clust_max!\n") ;
00169          EXIT(1) ;
00170       }
00171       edopt->edit_clust = ECFLAG_MAX;
00172       edopt->clust_rmm  = strtod( argv[nopt++] , NULL ) ;
00173       edopt->clust_vmul = strtod( argv[nopt++] , NULL ) ;
00174       if( edopt->clust_rmm < 0 ){
00175          fprintf(stderr,"illegal value after -1clust_max\n") ;
00176          EXIT(1) ;
00177       }
00178       CHECK_DONE ;
00179    }
00180 
00181    /**** -1clust_amax rmm vmul ****/   /* 10 Sept 1996 */
00182 
00183    if( strncmp(argv[nopt],"-1clust_amax",12) == 0 ){
00184       nopt++ ;
00185       if( nopt+1 >= argc ){
00186          fprintf(stderr,"need 2 arguments after -1clust_amax!\n") ;
00187          EXIT(1) ;
00188       }
00189       edopt->edit_clust = ECFLAG_AMAX;
00190       edopt->clust_rmm  = strtod( argv[nopt++] , NULL ) ;
00191       edopt->clust_vmul = strtod( argv[nopt++] , NULL ) ;
00192       if( edopt->clust_rmm < 0 ){
00193          fprintf(stderr,"illegal value after -1clust_amax\n") ;
00194          EXIT(1) ;
00195       }
00196       CHECK_DONE ;
00197    }
00198 
00199    /**** -1clust_smax rmm vmul ****/   /* 10 Sept 1996 */
00200 
00201    if( strncmp(argv[nopt],"-1clust_smax",12) == 0 ){
00202       nopt++ ;
00203       if( nopt+1 >= argc ){
00204          fprintf(stderr,"need 2 arguments after -1clust_smax!\n") ;
00205          EXIT(1) ;
00206       }
00207       edopt->edit_clust = ECFLAG_SMAX;
00208       edopt->clust_rmm  = strtod( argv[nopt++] , NULL ) ;
00209       edopt->clust_vmul = strtod( argv[nopt++] , NULL ) ;
00210       if( edopt->clust_rmm < 0 ){
00211          fprintf(stderr,"illegal value after -1clust_smax\n") ;
00212          EXIT(1) ;
00213       }
00214       CHECK_DONE ;
00215    }
00216 
00217    /**** -1clust_size rmm vmul ****/   /* 10 Sept 1996 */
00218 
00219    if( strncmp(argv[nopt],"-1clust_size",12) == 0 ){
00220       nopt++ ;
00221       if( nopt+1 >= argc ){
00222          fprintf(stderr,"need 2 arguments after -1clust_size!\n") ;
00223          EXIT(1) ;
00224       }
00225       edopt->edit_clust = ECFLAG_SIZE;
00226       edopt->clust_rmm  = strtod( argv[nopt++] , NULL ) ;
00227       edopt->clust_vmul = strtod( argv[nopt++] , NULL ) ;
00228       if( edopt->clust_rmm < 0 ){
00229          fprintf(stderr,"illegal value after -1clust_size\n") ;
00230          EXIT(1) ;
00231       }
00232       CHECK_DONE ;
00233    }
00234 
00235    /**** -1clust_order rmm vmul ****/   /* 09 June 1998 */
00236 
00237    if( strncmp(argv[nopt],"-1clust_order",12) == 0 ){
00238       nopt++ ;
00239       if( nopt+1 >= argc ){
00240          fprintf(stderr,"need 2 arguments after -1clust_order!\n") ;
00241          EXIT(1) ;
00242       }
00243       edopt->edit_clust = ECFLAG_ORDER;
00244       edopt->clust_rmm  = strtod( argv[nopt++] , NULL ) ;
00245       edopt->clust_vmul = strtod( argv[nopt++] , NULL ) ;
00246       if( edopt->clust_rmm < 0 ){
00247          fprintf(stderr,"illegal value after -1clust_order\n") ;
00248          EXIT(1) ;
00249       }
00250       CHECK_DONE ;
00251    }
00252 
00253    /**** -1erode pv ****/   /* 17 June 1998 */
00254 
00255    if( strncmp(argv[nopt],"-1erode",7) == 0 ){
00256       nopt++ ;
00257       if( nopt >= argc ){
00258          fprintf(stderr,"need 1 argument after -1erode!\n") ;
00259          EXIT(1) ;
00260       }
00261       edopt->erode_pv  = strtod( argv[nopt++] , NULL ) ;
00262       if (edopt->erode_pv > 1.0)  edopt->erode_pv /= 100.0;
00263       if( edopt->erode_pv < 0.0 || edopt->erode_pv > 1.0 ){
00264          fprintf(stderr,"illegal value after -1erode \n") ;
00265          EXIT(1) ;
00266       }
00267       CHECK_DONE ;
00268    }
00269 
00270    /**** -1dilate ****/   /* 17 June 1998 */
00271 
00272    if( strncmp(argv[nopt],"-1dilate",8) == 0 ){
00273       nopt++ ;
00274       edopt->dilate = 1;
00275       CHECK_DONE ;
00276    }
00277 
00278    /**** -1filter_mean rmm ****/   /* 11 Sept 1996 */
00279 
00280    if( strncmp(argv[nopt],"-1filter_mean",15) == 0 ){
00281       nopt++ ;
00282       if( nopt >= argc ){
00283          fprintf(stderr,"need 1 argument  after -1filter_mean \n") ;
00284          EXIT(1) ;
00285       }
00286       edopt->filter_opt = FCFLAG_MEAN;
00287       edopt->filter_rmm  = strtod( argv[nopt++] , NULL ) ;
00288       if( edopt->filter_rmm <= 0 ){
00289          fprintf(stderr,"illegal value after -1filter_mean \n") ;
00290          EXIT(1) ;
00291       }
00292       CHECK_DONE ;
00293    }
00294 
00295    /**** -1filter_nzmean rmm ****/   /* 11 Sept 1996 */
00296 
00297    if( strncmp(argv[nopt],"-1filter_nzmean",15) == 0 ){
00298       nopt++ ;
00299       if( nopt >= argc ){
00300          fprintf(stderr,"need 1 argument  after -1filter_nzmean \n") ;
00301          EXIT(1) ;
00302       }
00303       edopt->filter_opt = FCFLAG_NZMEAN;
00304       edopt->filter_rmm  = strtod( argv[nopt++] , NULL ) ;
00305       if( edopt->filter_rmm <= 0 ){
00306          fprintf(stderr,"illegal value after -1filter_nzmean \n") ;
00307          EXIT(1) ;
00308       }
00309       CHECK_DONE ;
00310    }
00311 
00312    /**** -1filter_max rmm ****/   /* 11 Sept 1996 */
00313 
00314    if( strncmp(argv[nopt],"-1filter_max",15) == 0 ){
00315       nopt++ ;
00316       if( nopt >= argc ){
00317          fprintf(stderr,"need 1 argument  after -1filter_max \n") ;
00318          EXIT(1) ;
00319       }
00320       edopt->filter_opt = FCFLAG_MAX;
00321       edopt->filter_rmm  = strtod( argv[nopt++] , NULL ) ;
00322       if( edopt->filter_rmm <= 0 ){
00323          fprintf(stderr,"illegal value after -1filter_max \n") ;
00324          EXIT(1) ;
00325       }
00326       CHECK_DONE ;
00327    }
00328 
00329    /**** -1filter_amax rmm ****/   /* 11 Sept 1996 */
00330 
00331    if( strncmp(argv[nopt],"-1filter_amax",15) == 0 ){
00332       nopt++ ;
00333       if( nopt >= argc ){
00334          fprintf(stderr,"need 1 argument  after -1filter_amax \n") ;
00335          EXIT(1) ;
00336       }
00337       edopt->filter_opt = FCFLAG_AMAX;
00338       edopt->filter_rmm  = strtod( argv[nopt++] , NULL ) ;
00339       if( edopt->filter_rmm <= 0 ){
00340          fprintf(stderr,"illegal value after -1filter_amax \n") ;
00341          EXIT(1) ;
00342       }
00343       CHECK_DONE ;
00344    }
00345 
00346    /**** -1filter_smax rmm ****/   /* 11 Sept 1996 */
00347 
00348    if( strncmp(argv[nopt],"-1filter_smax",15) == 0 ){
00349       nopt++ ;
00350       if( nopt >= argc ){
00351          fprintf(stderr,"need 1 argument  after -1filter_smax \n") ;
00352          EXIT(1) ;
00353       }
00354       edopt->filter_opt = FCFLAG_SMAX;
00355       edopt->filter_rmm  = strtod( argv[nopt++] , NULL ) ;
00356       if( edopt->filter_rmm <= 0 ){
00357          fprintf(stderr,"illegal value after -1filter_smax \n") ;
00358          EXIT(1) ;
00359       }
00360       CHECK_DONE ;
00361    }
00362 
00363    /**** -1filter_aver rmm ****/   /* 07 Jan 1998 */
00364 
00365    if( strncmp(argv[nopt],"-1filter_aver",15) == 0 ){
00366       nopt++ ;
00367       if( nopt >= argc ){
00368          fprintf(stderr,"need 1 argument  after -1filter_aver \n") ;
00369          EXIT(1) ;
00370       }
00371       edopt->filter_opt = FCFLAG_AVER ;
00372       edopt->filter_rmm  = strtod( argv[nopt++] , NULL ) ;
00373       if( edopt->filter_rmm <= 0 ){
00374          fprintf(stderr,"illegal value after -1filter_aver \n") ;
00375          EXIT(1) ;
00376       }
00377 
00378       if( edopt->nfmask > 0 ) edopt->filter_opt = FCFLAG_MEAN ;
00379       CHECK_DONE ;
00380    }
00381 
00382 
00383    /**** -t1filter_aver rmm ****/   /* 07 Jan 1998 */
00384 
00385    if( strncmp(argv[nopt],"-t1filter_aver",15) == 0 ){
00386       nopt++ ;
00387       if( nopt >= argc ){
00388          fprintf(stderr,"need 1 argument  after -t1filter_aver \n") ;
00389          EXIT(1) ;
00390       }
00391       edopt->thrfilter_opt = FCFLAG_AVER ;
00392       edopt->thrfilter_rmm  = strtod( argv[nopt++] , NULL ) ;
00393       if( edopt->thrfilter_rmm <= 0 ){
00394          fprintf(stderr,"illegal value after -t1filter_aver \n") ;
00395          EXIT(1) ;
00396       }
00397 
00398       if( edopt->nfmask > 0 ) edopt->thrfilter_opt = FCFLAG_MEAN ;
00399       CHECK_DONE ;
00400    }
00401 
00402 
00403    /**** -t1filter_mean rmm ****/   /* 1 Oct 1996 */
00404 
00405    if( strncmp(argv[nopt],"-t1filter_mean",15) == 0 ){
00406       nopt++ ;
00407       if( nopt >= argc ){
00408          fprintf(stderr,"need 1 argument  after -t1filter_mean \n") ;
00409          EXIT(1) ;
00410       }
00411       edopt->thrfilter_opt = FCFLAG_MEAN;
00412       edopt->thrfilter_rmm  = strtod( argv[nopt++] , NULL ) ;
00413       if( edopt->thrfilter_rmm <= 0 ){
00414          fprintf(stderr,"illegal value after -t1filter_mean \n") ;
00415          EXIT(1) ;
00416       }
00417       CHECK_DONE ;
00418    }
00419 
00420 
00421    /**** -t1filter_nzmean rmm ****/   /* 1 Oct 1996 */
00422 
00423    if( strncmp(argv[nopt],"-t1filter_nzmean",15) == 0 ){
00424       nopt++ ;
00425       if( nopt >= argc ){
00426          fprintf(stderr,"need 1 argument  after -t1filter_nzmean \n") ;
00427          EXIT(1) ;
00428       }
00429       edopt->thrfilter_opt = FCFLAG_NZMEAN;
00430       edopt->thrfilter_rmm  = strtod( argv[nopt++] , NULL ) ;
00431       if( edopt->thrfilter_rmm <= 0 ){
00432          fprintf(stderr,"illegal value after -t1filter_nzmean \n") ;
00433          EXIT(1) ;
00434       }
00435       CHECK_DONE ;
00436    }
00437 
00438 
00439    /**** -t1filter_max rmm ****/   /* 1 Oct 1996 */
00440 
00441    if( strncmp(argv[nopt],"-t1filter_max",15) == 0 ){
00442       nopt++ ;
00443       if( nopt >= argc ){
00444          fprintf(stderr,"need 1 argument  after -t1filter_max \n") ;
00445          EXIT(1) ;
00446       }
00447       edopt->thrfilter_opt = FCFLAG_MAX;
00448       edopt->thrfilter_rmm  = strtod( argv[nopt++] , NULL ) ;
00449       if( edopt->thrfilter_rmm <= 0 ){
00450          fprintf(stderr,"illegal value after -t1filter_max \n") ;
00451          EXIT(1) ;
00452       }
00453       CHECK_DONE ;
00454    }
00455 
00456 
00457    /**** -t1filter_amax rmm ****/   /* 1 Oct 1996 */
00458 
00459    if( strncmp(argv[nopt],"-t1filter_amax",15) == 0 ){
00460       nopt++ ;
00461       if( nopt >= argc ){
00462          fprintf(stderr,"need 1 argument  after -t1filter_amax \n") ;
00463          EXIT(1) ;
00464       }
00465       edopt->thrfilter_opt = FCFLAG_AMAX;
00466       edopt->thrfilter_rmm  = strtod( argv[nopt++] , NULL ) ;
00467       if( edopt->thrfilter_rmm <= 0 ){
00468          fprintf(stderr,"illegal value after -t1filter_amax \n") ;
00469          EXIT(1) ;
00470       }
00471       CHECK_DONE ;
00472    }
00473 
00474 
00475    /**** -t1filter_smax rmm ****/   /* 1 Oct 1996 */
00476 
00477    if( strncmp(argv[nopt],"-t1filter_smax",15) == 0 ){
00478       nopt++ ;
00479       if( nopt >= argc ){
00480          fprintf(stderr,"need 1 argument  after -t1filter_smax \n") ;
00481          EXIT(1) ;
00482       }
00483       edopt->thrfilter_opt = FCFLAG_SMAX;
00484       edopt->thrfilter_rmm  = strtod( argv[nopt++] , NULL ) ;
00485       if( edopt->thrfilter_rmm <= 0 ){
00486          fprintf(stderr,"illegal value after -t1filter_smax \n") ;
00487          EXIT(1) ;
00488       }
00489       CHECK_DONE ;
00490    }
00491 
00492 
00493 
00494 
00495    /**** -1blur_sigma size ****/
00496 
00497    if( strncmp(argv[nopt],"-1blur_sigma",12) == 0 ){
00498       nopt++ ;
00499       if( nopt >= argc ){
00500          fprintf(stderr,"need argument after -1blur_sigma!\n") ;
00501          EXIT(1) ;
00502       }
00503       edopt->blur = strtod( argv[nopt++] , NULL ) ;
00504       if( edopt->blur <= 0 ){
00505          fprintf(stderr,"illegal value after -1blur_sigma\n") ;
00506          EXIT(1) ;
00507       }
00508       CHECK_DONE ;
00509    }
00510 
00511    /**** -1blur_rms size ****/
00512 
00513    if( strncmp(argv[nopt],"-1blur_rms",12) == 0 ){
00514       nopt++ ;
00515       if( nopt >= argc ){
00516          fprintf(stderr,"need argument after -1blur_rms!\n") ;
00517          EXIT(1) ;
00518       }
00519       edopt->blur = strtod( argv[nopt++] , NULL ) ;
00520       if( edopt->blur <= 0 ){
00521          fprintf(stderr,"illegal value after -1blur_rms\n") ;
00522          EXIT(1) ;
00523       }
00524       edopt->blur = RMS_TO_SIGMA(edopt->blur) ;
00525       CHECK_DONE ;
00526    }
00527 
00528    /**** -1blur_fwhm size ****/
00529 
00530    if( strncmp(argv[nopt],"-1blur_fwhm",12) == 0 ){
00531       nopt++ ;
00532       if( nopt >= argc ){
00533          fprintf(stderr,"need argument after -1blur_fwhm!\n") ;
00534          EXIT(1) ;
00535       }
00536       edopt->blur = strtod( argv[nopt++] , NULL ) ;
00537       if( edopt->blur <= 0 ){
00538          fprintf(stderr,"illegal value after -1blur_fwhm\n") ;
00539          EXIT(1) ;
00540       }
00541       edopt->blur = FWHM_TO_SIGMA(edopt->blur) ;
00542       CHECK_DONE ;
00543    }
00544 
00545    /**** -1blur ****/
00546 
00547    if( strncmp(argv[nopt],"-1blur",6) == 0 ){
00548       fprintf(stderr,
00549               "*** the old -1blur option is no longer valid! ***\n") ;
00550       EXIT(1) ;
00551    }
00552 
00553    /**** -t1blur_sigma size ****/   /* 4 Oct 1996 */
00554 
00555    if( strncmp(argv[nopt],"-t1blur_sigma",12) == 0 ){
00556       nopt++ ;
00557       if( nopt >= argc ){
00558          fprintf(stderr,"need argument after -t1blur_sigma!\n") ;
00559          EXIT(1) ;
00560       }
00561       edopt->thrblur = strtod( argv[nopt++] , NULL ) ;
00562       if( edopt->thrblur <= 0 ){
00563          fprintf(stderr,"illegal value after -t1blur_sigma\n") ;
00564          EXIT(1) ;
00565       }
00566       CHECK_DONE ;
00567    }
00568 
00569    /**** -t1blur_rms size ****/   /* 4 Oct 1996 */
00570 
00571    if( strncmp(argv[nopt],"-t1blur_rms",12) == 0 ){
00572       nopt++ ;
00573       if( nopt >= argc ){
00574          fprintf(stderr,"need argument after -t1blur_rms!\n") ;
00575          EXIT(1) ;
00576       }
00577       edopt->thrblur = strtod( argv[nopt++] , NULL ) ;
00578       if( edopt->thrblur <= 0 ){
00579          fprintf(stderr,"illegal value after -t1blur_rms\n") ;
00580          EXIT(1) ;
00581       }
00582       edopt->thrblur = RMS_TO_SIGMA(edopt->thrblur) ;
00583       CHECK_DONE ;
00584    }
00585 
00586    /**** -t1blur_fwhm size ****/   /* 4 Oct 1996 */
00587 
00588    if( strncmp(argv[nopt],"-t1blur_fwhm",12) == 0 ){
00589       nopt++ ;
00590       if( nopt >= argc ){
00591          fprintf(stderr,"need argument after -t1blur_fwhm!\n") ;
00592          EXIT(1) ;
00593       }
00594       edopt->thrblur = strtod( argv[nopt++] , NULL ) ;
00595       if( edopt->thrblur <= 0 ){
00596          fprintf(stderr,"illegal value after -t1blur_fwhm\n") ;
00597          EXIT(1) ;
00598       }
00599       edopt->thrblur = FWHM_TO_SIGMA(edopt->thrblur) ;
00600       CHECK_DONE ;
00601    }
00602 
00603 
00604    /**** -1scale ****/
00605 
00606    if( strncmp(argv[nopt],"-1scale",6) == 0 ){
00607 #ifdef ALLOW_SCALE_TO_MAX
00608       edopt->scale = 1 ;
00609 #else
00610       fprintf(stderr,
00611               "*** the old -1scale option is no longer valid! ***\n") ;
00612 #endif
00613       nopt++ ; CHECK_DONE ;
00614    }
00615 
00616    /**** -1mult mult ****/
00617 
00618    if( strncmp(argv[nopt],"-1mult",6) == 0 ){
00619       nopt++ ;
00620       if( nopt >= argc ){
00621          fprintf(stderr,"no argument after -1mult!\n") ;
00622          EXIT(1);
00623       }
00624       val = strtod( argv[nopt++] , NULL ) ;
00625       if( val == 0.0 ){
00626          fprintf(stderr,"illegal value after -1mult!\n") ;
00627          EXIT(1) ;
00628       }
00629       edopt->mult = val ;
00630       CHECK_DONE ;
00631     }
00632 
00633     /**** -1zvol x1 x2 y1 y2 z1 z2 ***/
00634 
00635     if( strncmp(argv[nopt],"-1zvol",6) == 0 ){
00636       char * cerr ;
00637 
00638       if( nopt+6 >= argc ){
00639          fprintf(stderr,"need 6 arguments after -1zvol!\a\n") ;
00640          EXIT(1) ;
00641       }
00642 
00643       edopt->zv_x1 = strtod( argv[nopt+1] , &cerr ) ;
00644       if( cerr == argv[nopt+1] ){
00645          fprintf(stderr,"illegal 1st argument after -1zvol!\a\n") ;
00646          EXIT(1) ;
00647       }
00648 
00649       edopt->zv_x2 = strtod( argv[nopt+2] , &cerr ) ;
00650       if( cerr == argv[nopt+2] ){
00651          fprintf(stderr,"illegal 2nd argument after -1zvol!\a\n") ;
00652          EXIT(1) ;
00653       }
00654 
00655       edopt->zv_y1 = strtod( argv[nopt+3] , &cerr ) ;
00656       if( cerr == argv[nopt+3] ){
00657          fprintf(stderr,"illegal 3rd argument after -1zvol!\a\n") ;
00658          EXIT(1) ;
00659       }
00660 
00661       edopt->zv_y2 = strtod( argv[nopt+4] , &cerr ) ;
00662       if( cerr == argv[nopt+4] ){
00663          fprintf(stderr,"illegal 4th argument after -1zvol!\a\n") ;
00664          EXIT(1) ;
00665       }
00666 
00667       edopt->zv_z1 = strtod( argv[nopt+5] , &cerr ) ;
00668       if( cerr == argv[nopt+5] ){
00669          fprintf(stderr,"illegal 5th argument after -1zvol!\a\n") ;
00670          EXIT(1) ;
00671       }
00672 
00673       edopt->zv_z2 = strtod( argv[nopt+6] , &cerr ) ;
00674       if( cerr == argv[nopt+6] ){
00675          fprintf(stderr,"illegal 6th argument after -1zvol!\a\n") ;
00676          EXIT(1) ;
00677       }
00678       edopt->do_zvol = 1 ;
00679 
00680       nopt += 7 ; CHECK_DONE ;
00681    }
00682 
00683    RETURN( 0 );
00684 }

00402 {
00403    int ii ;
00404    float bot ;
00405 
00406 ENTRY("EDIT_clip_float") ;
00407 
00408    if( top <= 0.0 || nxyz <= 0 || vol == NULL ) EXRETURN ;
00409 
00410    bot = -top ;
00411 
00412    for( ii=0 ; ii < nxyz ; ii++ )
00413            if( vol[ii] > top ) vol[ii] = top ;
00414       else if( vol[ii] < bot ) vol[ii] = bot ;
00415 
00416    EXRETURN ;
00417 }

00027 {
00028    int iclu;       /* cluster index */
00029    int nclu;       /* non-empty cluster index */
00030    int ii;         /* voxel index */
00031    float
00032       mag,         /* voxel intensity */
00033       sum,         /* sum of voxel intensities */
00034       max,         /* maximum of voxel intensities */
00035       amax,        /* maximum of absolute voxel intensities */
00036       smax,        /* signed maximum of absolute voxel intensities */
00037       mean,        /* mean of voxel intensities */
00038       size;        /* size of cluster (multiples of vmul) */
00039 
00040 ENTRY("EDIT_cluster_array") ;
00041 
00042    if( edit_clust == ECFLAG_ORDER ){
00043       SORT_CLARR(clar) ;
00044    }
00045 
00046    nclu = 0;
00047    for (iclu = 0; iclu < clar->num_clu; iclu++)
00048    {
00049       if ((clar->clar[iclu] != NULL) && (clar->clar[iclu]->num_pt > 0))
00050       {
00051          nclu++;
00052 
00053          /* initialization of basic statistics for this cluster */
00054          sum = max = smax = clar->clar[iclu]->mag[0];
00055          amax = fabs(smax);
00056 
00057          /* calculate basic statistics for this cluster */
00058          for (ii = 1; ii < clar->clar[iclu]->num_pt; ii++)
00059          {
00060             mag = clar->clar[iclu]->mag[ii];
00061             switch (edit_clust)
00062             {
00063                case ECFLAG_MEAN :
00064                   sum += mag;  break;
00065                case ECFLAG_MAX  :
00066                   if (mag > max)  max = mag;   break;
00067                case ECFLAG_AMAX :
00068                   if (fabs(mag) > amax)  amax = fabs(mag);  break;
00069                case ECFLAG_SMAX :
00070                   if (fabs(mag) > fabs(smax))  smax = mag;  break;
00071                case ECFLAG_SIZE : break;
00072                default          : break;
00073             }
00074 
00075          }
00076 
00077          /* additional calculations */
00078          if (edit_clust == ECFLAG_MEAN)
00079             mean = sum / clar->clar[iclu]->num_pt;
00080          if (edit_clust == ECFLAG_SIZE)
00081             size = clar->clar[iclu]->num_pt * dxyz / vmul;
00082 
00083          /* set all voxel intensities in this cluster to the same value */
00084          for (ii = 0; ii < clar->clar[iclu]->num_pt; ii++)
00085          {
00086             switch (edit_clust)
00087             {
00088                case ECFLAG_MEAN :  clar->clar[iclu]->mag[ii] = mean;
00089                case ECFLAG_MAX  :  clar->clar[iclu]->mag[ii] = max;   break;
00090                case ECFLAG_AMAX :  clar->clar[iclu]->mag[ii] = amax;  break;
00091                case ECFLAG_SMAX :  clar->clar[iclu]->mag[ii] = smax;  break;
00092                case ECFLAG_SIZE :  clar->clar[iclu]->mag[ii] = size;  break;
00093                case ECFLAG_ORDER:  clar->clar[iclu]->mag[ii] = nclu;  break;
00094                default          :                                     break;
00095             }
00096          }
00097       }
00098    }  /* iclu */
00099 
00100    EXRETURN ;
00101 }

00385 {
00386    float fac=0.0 , top ;
00387 
00388 ENTRY("EDIT_coerce_autoscale") ;
00389 
00390    if( MRI_IS_INT_TYPE(otype) ){
00391       top = MCW_vol_amax( nxyz,1,1 , itype,ivol ) ;
00392       fac = (top > MRI_TYPE_maxval[otype]) ? MRI_TYPE_maxval[otype]/top : 0.0 ;
00393    }
00394 
00395    EDIT_coerce_scale_type( nxyz , fac , itype,ivol , otype,ovol ) ;
00396    RETURN( fac );
00397 }

00196 {
00197    register int ii ;
00198    register float fac = scl , val ;
00199 
00200    complex * cin , * cout ;
00201    short   * sin , * sout ;
00202    float   * fin , * fout ;
00203    byte    * bin , * bout ;
00204    double  * din , * dout ;   /* 11 Jan 1999 */
00205 
00206 ENTRY("EDIT_coerce_scale_type") ;
00207 #ifdef AFNI_DEBUG
00208 { char str[256] ;
00209   sprintf(str,"voxels=%d scale=%g input type=%s output type=%s",
00210           nxyz,scl , MRI_TYPE_name[itype],MRI_TYPE_name[otype]) ;
00211   STATUS(str) ; }
00212 #endif
00213 
00214    if( nxyz <= 0 || ivol == NULL || ovol == NULL ) EXRETURN ;
00215 
00216    if( fac == 0.0 || fac == 1.0 ){
00217       EDIT_coerce_type( nxyz , itype,ivol , otype,ovol ) ;
00218       EXRETURN ;
00219    }
00220 
00221    switch( itype ){
00222       default:
00223         fprintf(stderr,"** Unknown itype=%d in EDIT_coerce_scale_type\n",itype);
00224       EXRETURN ;
00225       case MRI_complex:  cin = (complex *) ivol ; break ;
00226       case MRI_short  :  sin = (short   *) ivol ; break ;
00227       case MRI_float  :  fin = (float   *) ivol ; break ;
00228       case MRI_byte   :  bin = (byte    *) ivol ; break ;
00229       case MRI_double :  din = (double  *) ivol ; break ;
00230    }
00231    switch( otype ){
00232       default:
00233         fprintf(stderr,"** Unknown otype=%d in EDIT_coerce_scale_type\n",otype);
00234       EXRETURN ;
00235       case MRI_complex:  cout = (complex *) ovol ; break ;
00236       case MRI_short  :  sout = (short   *) ovol ; break ;
00237       case MRI_float  :  fout = (float   *) ovol ; break ;
00238       case MRI_byte   :  bout = (byte    *) ovol ; break ;
00239       case MRI_double :  dout = (double  *) ovol ; break ;
00240    }
00241 
00242    switch( otype ){
00243 
00244       /*** outputs are shorts ***/
00245 
00246       case MRI_short:
00247          switch( itype ){
00248             case MRI_short:   /* inputs are shorts */
00249                for( ii=0 ; ii < nxyz ; ii++ ) sout[ii] = ROUND(fac*sin[ii]) ;
00250                EXRETURN ;
00251             case MRI_float:   /* inputs are floats */
00252                for( ii=0 ; ii < nxyz ; ii++ ) sout[ii] = ROUND(fac*fin[ii]) ;
00253                EXRETURN ;
00254             case MRI_double:   /* inputs are double */
00255                for( ii=0 ; ii < nxyz ; ii++ ) sout[ii] = ROUND(fac*din[ii]) ;
00256                EXRETURN ;
00257             case MRI_byte:    /* inputs are bytes */
00258                for( ii=0 ; ii < nxyz ; ii++ ) sout[ii] = ROUND(fac*bin[ii]) ;
00259                EXRETURN ;
00260             case MRI_complex:    /* inputs are complex */
00261                for( ii=0 ; ii < nxyz ; ii++ ) sout[ii] = ROUND(fac*CABS(cin[ii])) ;
00262                EXRETURN ;
00263          }
00264          EXRETURN ;
00265 
00266       /*** outputs are floats ***/
00267 
00268       case MRI_float:
00269          switch( itype ){
00270             case MRI_short:   /* inputs are shorts */
00271                for( ii=0 ; ii < nxyz ; ii++ ) fout[ii] = fac*sin[ii] ;
00272                EXRETURN ;
00273             case MRI_float:   /* inputs are floats */
00274                for( ii=0 ; ii < nxyz ; ii++ ) fout[ii] = fac*fin[ii] ;
00275                EXRETURN ;
00276             case MRI_double:   /* inputs are doubles */
00277                for( ii=0 ; ii < nxyz ; ii++ ) fout[ii] = fac*din[ii] ;
00278                EXRETURN ;
00279             case MRI_byte:    /* inputs are bytes */
00280                for( ii=0 ; ii < nxyz ; ii++ ) fout[ii] = fac*bin[ii] ;
00281                EXRETURN ;
00282             case MRI_complex:    /* inputs are complex */
00283                for( ii=0 ; ii < nxyz ; ii++ ) fout[ii] = fac*CABS(cin[ii]) ;
00284                EXRETURN ;
00285          }
00286          EXRETURN ;
00287 
00288       /*** outputs are doubles ***/
00289 
00290       case MRI_double:
00291          switch( itype ){
00292             case MRI_short:   /* inputs are shorts */
00293                for( ii=0 ; ii < nxyz ; ii++ ) dout[ii] = fac*sin[ii] ;
00294                EXRETURN ;
00295             case MRI_float:   /* inputs are floats */
00296                for( ii=0 ; ii < nxyz ; ii++ ) dout[ii] = fac*fin[ii] ;
00297                EXRETURN ;
00298             case MRI_double:   /* inputs are doubles */
00299                for( ii=0 ; ii < nxyz ; ii++ ) dout[ii] = fac*din[ii] ;
00300                EXRETURN ;
00301             case MRI_byte:    /* inputs are bytes */
00302                for( ii=0 ; ii < nxyz ; ii++ ) dout[ii] = fac*bin[ii] ;
00303                EXRETURN ;
00304             case MRI_complex:    /* inputs are complex */
00305                for( ii=0 ; ii < nxyz ; ii++ ) dout[ii] = fac*CABS(cin[ii]) ;
00306                EXRETURN ;
00307          }
00308          EXRETURN ;
00309 
00310       /*** outputs are bytes ***/
00311 
00312       case MRI_byte:
00313          switch( itype ){
00314             case MRI_short:   /* inputs are shorts */
00315                for( ii=0 ; ii < nxyz ; ii++ ){
00316                   val = fac*sin[ii] ; bout[ii] = FLOAT_TO_BYTE(val) ;
00317                }
00318                EXRETURN ;
00319             case MRI_float:   /* inputs are floats */
00320                for( ii=0 ; ii < nxyz ; ii++ ){
00321                   val = fac*fin[ii] ; bout[ii] = FLOAT_TO_BYTE(val) ;
00322                }
00323                EXRETURN ;
00324             case MRI_double:   /* inputs are doubles */
00325                for( ii=0 ; ii < nxyz ; ii++ ){
00326                   val = fac*din[ii] ; bout[ii] = FLOAT_TO_BYTE(val) ;
00327                }
00328                EXRETURN ;
00329             case MRI_byte:    /* inputs are bytes */
00330                for( ii=0 ; ii < nxyz ; ii++ ){
00331                   val = fac*bin[ii] ; bout[ii] = FLOAT_TO_BYTE(val) ;
00332                }
00333                EXRETURN ;
00334             case MRI_complex:{    /* inputs are complex */
00335                for( ii=0 ; ii < nxyz ; ii++ ){
00336                   val = fac*CABS(cin[ii]) ; bout[ii] = FLOAT_TO_BYTE(val) ;
00337                }
00338             }
00339             EXRETURN ;
00340          }
00341          EXRETURN ;
00342 
00343       /*** outputs are complex ***/
00344 
00345       case MRI_complex:
00346          switch( itype ){
00347             case MRI_short:   /* inputs are shorts */
00348                for( ii=0 ; ii < nxyz ; ii++ )
00349                   cout[ii].r = fac*sin[ii] , cout[ii].i = 0.0 ;
00350                EXRETURN ;
00351             case MRI_float:   /* inputs are floats */
00352                for( ii=0 ; ii < nxyz ; ii++ )
00353                   cout[ii].r = fac*fin[ii] , cout[ii].i = 0.0 ;
00354                EXRETURN ;
00355             case MRI_double:   /* inputs are doubles */
00356                for( ii=0 ; ii < nxyz ; ii++ )
00357                   cout[ii].r = fac*din[ii] , cout[ii].i = 0.0 ;
00358                EXRETURN ;
00359             case MRI_byte:    /* inputs are bytes */
00360                for( ii=0 ; ii < nxyz ; ii++ )
00361                   cout[ii].r = fac*bin[ii] , cout[ii].i = 0.0 ;
00362                EXRETURN ;
00363             case MRI_complex:    /* inputs are complex */
00364                for( ii=0 ; ii < nxyz ; ii++ )
00365                   cout[ii].r = fac*cin[ii].r , cout[ii].i = fac*cin[ii].i ;
00366                EXRETURN ;
00367          }
00368          EXRETURN ;
00369    }
00370    EXRETURN ;
00371 }

00022 {
00023    register int ii ;
00024 
00025    complex * cin , * cout ;
00026    short   * sin , * sout ;
00027    float   * fin , * fout ;
00028    byte    * bin , * bout ;
00029    double  * din , * dout ;  /* 10 Jan 1999 */
00030 
00031 ENTRY("EDIT_coerce_type") ;
00032 #ifdef AFNI_DEBUG
00033 { char str[256] ;
00034   sprintf(str,"voxels=%d input type=%s output type=%s",
00035           nxyz, MRI_TYPE_name[itype],MRI_TYPE_name[otype]) ;
00036   STATUS(str) ; }
00037 #endif
00038 
00039    if( nxyz <= 0 || ivol == NULL || ovol == NULL ) EXRETURN ;
00040 
00041    switch( itype ){
00042       default:
00043         fprintf(stderr,"** Unknown itype=%d in EDIT_coerce_type\n",itype);
00044       EXRETURN ;
00045       case MRI_complex:  cin = (complex *) ivol ; break ;
00046       case MRI_short  :  sin = (short   *) ivol ; break ;
00047       case MRI_float  :  fin = (float   *) ivol ; break ;
00048       case MRI_byte   :  bin = (byte    *) ivol ; break ;
00049       case MRI_double :  din = (double  *) ivol ; break ;
00050    }
00051    switch( otype ){
00052       default:
00053         fprintf(stderr,"** Unknown otype=%d in EDIT_coerce_type\n",otype);
00054       EXRETURN ;
00055       case MRI_complex:  cout = (complex *) ovol ; break ;
00056       case MRI_short  :  sout = (short   *) ovol ; break ;
00057       case MRI_float  :  fout = (float   *) ovol ; break ;
00058       case MRI_byte   :  bout = (byte    *) ovol ; break ;
00059       case MRI_double :  dout = (double  *) ovol ; break ;
00060    }
00061 
00062    switch( otype ){
00063 
00064       /*** outputs are shorts ***/
00065 
00066       case MRI_short:
00067          switch( itype ){
00068             case MRI_short:   /* inputs are shorts */
00069                memcpy( sout , sin , sizeof(short)*nxyz ) ;
00070                EXRETURN ;
00071             case MRI_float:   /* inputs are floats */
00072                for( ii=0 ; ii < nxyz ; ii++ ) sout[ii] = ROUND(fin[ii]) ;
00073                EXRETURN ;
00074             case MRI_double:   /* inputs are doubles */
00075                for( ii=0 ; ii < nxyz ; ii++ ) sout[ii] = ROUND(din[ii]) ;
00076                EXRETURN ;
00077             case MRI_byte:    /* inputs are bytes */
00078                for( ii=0 ; ii < nxyz ; ii++ ) sout[ii] = bin[ii] ;
00079                EXRETURN ;
00080             case MRI_complex:    /* inputs are complex */
00081                for( ii=0 ; ii < nxyz ; ii++ ) sout[ii] = ROUND(CABS(cin[ii])) ;
00082                EXRETURN ;
00083          }
00084          EXRETURN ;
00085 
00086       /*** outputs are floats ***/
00087 
00088       case MRI_float:
00089          switch( itype ){
00090             case MRI_short:   /* inputs are shorts */
00091                for( ii=0 ; ii < nxyz ; ii++ ) fout[ii] = sin[ii] ;
00092                EXRETURN ;
00093             case MRI_float:   /* inputs are floats */
00094                memcpy( fout , fin , sizeof(float)*nxyz ) ;
00095                EXRETURN ;
00096             case MRI_double:    /* inputs are doubles */
00097                for( ii=0 ; ii < nxyz ; ii++ ) fout[ii] = din[ii] ;
00098                EXRETURN ;
00099             case MRI_byte:    /* inputs are bytes */
00100                for( ii=0 ; ii < nxyz ; ii++ ) fout[ii] = bin[ii] ;
00101                EXRETURN ;
00102             case MRI_complex:    /* inputs are complex */
00103                for( ii=0 ; ii < nxyz ; ii++ ) fout[ii] = CABS(cin[ii]) ;
00104                EXRETURN ;
00105          }
00106          EXRETURN ;
00107 
00108       /*** outputs are doubles ***/
00109 
00110       case MRI_double:
00111          switch( itype ){
00112             case MRI_short:   /* inputs are shorts */
00113                for( ii=0 ; ii < nxyz ; ii++ ) dout[ii] = sin[ii] ;
00114                EXRETURN ;
00115             case MRI_float:   /* inputs are floats */
00116                for( ii=0 ; ii < nxyz ; ii++ ) dout[ii] = fin[ii] ;
00117                EXRETURN ;
00118             case MRI_double:   /* inputs are doubles */
00119                memcpy( dout , din , sizeof(double)*nxyz ) ;
00120                EXRETURN ;
00121             case MRI_byte:    /* inputs are bytes */
00122                for( ii=0 ; ii < nxyz ; ii++ ) dout[ii] = bin[ii] ;
00123                EXRETURN ;
00124             case MRI_complex:    /* inputs are complex */
00125                for( ii=0 ; ii < nxyz ; ii++ ) dout[ii] = CABS(cin[ii]) ;
00126                EXRETURN ;
00127          }
00128          EXRETURN ;
00129 
00130       /*** outputs are bytes ***/
00131 
00132       case MRI_byte:
00133          switch( itype ){
00134             case MRI_short:   /* inputs are shorts */
00135                for( ii=0 ; ii < nxyz ; ii++ ) bout[ii] = SHORT_TO_BYTE(sin[ii]) ;
00136                EXRETURN ;
00137             case MRI_float:   /* inputs are floats */
00138                for( ii=0 ; ii < nxyz ; ii++ ) bout[ii] = FLOAT_TO_BYTE(fin[ii]) ;
00139                EXRETURN ;
00140             case MRI_double:   /* inputs are doubles */
00141                for( ii=0 ; ii < nxyz ; ii++ ) bout[ii] = FLOAT_TO_BYTE(din[ii]) ;
00142                EXRETURN ;
00143             case MRI_byte:    /* inputs are bytes */
00144                memcpy( bout , bin , sizeof(byte)*nxyz ) ;
00145                EXRETURN ;
00146             case MRI_complex:{    /* inputs are complex */
00147                float val ;
00148                for( ii=0 ; ii < nxyz ; ii++ ){ val = CABS(cin[ii]) ;
00149                                                bout[ii] = FLOAT_TO_BYTE(val) ; }
00150             }
00151             EXRETURN ;
00152          }
00153          EXRETURN ;
00154 
00155       /*** outputs are complex ***/
00156 
00157       case MRI_complex:
00158          switch( itype ){
00159             case MRI_short:   /* inputs are shorts */
00160                for( ii=0 ; ii < nxyz ; ii++ )
00161                   cout[ii].r = sin[ii] , cout[ii].i = 0.0 ;
00162                EXRETURN ;
00163             case MRI_float:   /* inputs are floats */
00164                for( ii=0 ; ii < nxyz ; ii++ )
00165                   cout[ii].r = fin[ii] , cout[ii].i = 0.0 ;
00166                EXRETURN ;
00167             case MRI_double:   /* inputs are doubles */
00168                for( ii=0 ; ii < nxyz ; ii++ )
00169                   cout[ii].r = din[ii] , cout[ii].i = 0.0 ;
00170                EXRETURN ;
00171             case MRI_byte:    /* inputs are bytes */
00172                for( ii=0 ; ii < nxyz ; ii++ )
00173                   cout[ii].r = bin[ii] , cout[ii].i = 0.0 ;
00174                EXRETURN ;
00175             case MRI_complex:    /* inputs are complex */
00176                memcpy( cout , cin , sizeof(complex)*nxyz ) ;
00177                EXRETURN ;
00178          }
00179          EXRETURN ;
00180    }
00181    EXRETURN ;
00182 }

00035 {
00036    va_list vararg_ptr ;
00037    int     flag_arg , errnum = 0 ;
00038    int     redo_bricks , redo_daxes , redo_taxis , ii ;
00039    void *  dummy ;
00040    int     iarg ;
00041 
00042    /**----- variables to flag and store presence of arguments -----**/
00043 
00044    int new_prefix         = 0 ; char *              prefix         = NULL ;
00045    int new_directory_name = 0 ; char *              directory_name = NULL ;
00046    int new_brick_fac      = 0 ; float *             brick_fac      = NULL ;
00047    int new_malloc_type    = 0 ; int                 malloc_type    = ILLEGAL_TYPE ;
00048    int new_datum_all      = 0 ; int                 datum_all      = ILLEGAL_TYPE ;
00049    int new_datum_array    = 0 ; int *               datum_array    = NULL ;
00050    int new_nvals          = 0 ; int                 nvals          = 0 ;
00051    int new_nxyz           = 0 ; THD_ivec3           nxyz           ;
00052    int new_xyzdel         = 0 ; THD_fvec3           xyzdel         ;
00053    int new_xyzorg         = 0 ; THD_fvec3           xyzorg         ;
00054    int new_xyzorient      = 0 ; THD_ivec3           xyzorient      ;
00055    int new_to_dicomm      = 0 ; THD_mat33           to_dicomm      ;
00056    int new_ntt            = 0 ; int                 ntt            = 0 ;
00057    int new_ttorg          = 0 ; float               ttorg          = 0.0 ;
00058    int new_ttdel          = 0 ; float               ttdel          = 0.0 ;
00059    int new_ttdur          = 0 ; float               ttdur          = 0.0 ;
00060    int new_nsl            = 0 ; int                 nsl            = 0 ;
00061    int new_zorg_sl        = 0 ; float               zorg_sl        = 0.0 ;
00062    int new_dz_sl          = 0 ; float               dz_sl          = 0.0 ;
00063    int new_toff_sl        = 0 ; float *             toff_sl        = NULL ;
00064    int new_type           = 0 ; int                 type           = ILLEGAL_TYPE ;
00065    int new_view_type      = 0 ; int                 view_type      = ILLEGAL_TYPE ;
00066    int new_func_type      = 0 ; int                 func_type      = ILLEGAL_TYPE ;
00067    int new_label1         = 0 ; char *              label1         = NULL ;
00068    int new_label2         = 0 ; char *              label2         = NULL ;
00069    int new_self_name      = 0 ; char *              self_name      = NULL ;
00070    int new_warp_parent    = 0 ; THD_3dim_dataset *  warp_parent    = NULL ;
00071    int new_anat_parent    = 0 ; THD_3dim_dataset *  anat_parent    = NULL ;
00072    int new_stat_aux       = 0 ; float *             stat_aux       = NULL ;
00073    int new_warp           = 0 ; THD_warp *          warp           = NULL ;
00074    int new_tunits         = 0 ; int                 tunits         = ILLEGAL_TYPE ;
00075    int new_keywords       = 0 ; char *              keywords       = NULL ;
00076 
00077    /* 30 Nov 1997 */
00078 
00079    int new_brick_label_one = 0 ; char * brick_label_one    = NULL ;
00080                                  int    brick_label_one_iv = -1 ;
00081 
00082    int new_brick_fac_one = 0 ; float brick_fac_one    = 0.0 ;
00083                                int   brick_fac_one_iv = -1 ;
00084 
00085    int new_brick_stataux_one = 0 ; float * brick_stataux_one    = NULL ;
00086                                    int     brick_stataux_one_iv = -1 ;
00087 
00088    int new_brick_keywords_one = 0 ; char * brick_keywords_one    = NULL ;
00089                                     int    brick_keywords_one_iv = -1 ;
00090 
00091    /****---------------------- Sanity Check ----------------------****/
00092 
00093 ENTRY("EDIT_dset_items") ;
00094 
00095    if( ! ISVALID_3DIM_DATASET(dset) ) RETURN(1) ;  /* bad data */
00096 
00097    /****----------- Scan input argument list;
00098                   - Load data into locals (va_arg);
00099                   - Check for legal values;
00100                   - Flag its presence (the new_ variables);
00101                   - Carry out simple processing that doesn't
00102                       depend on the presence of other arguments. ---------****/
00103 
00104    va_start( vararg_ptr , dset ) ;              /** Initialize arg reading **/
00105    iarg = 1 ;
00106 
00107    do{
00108       flag_arg = va_arg( vararg_ptr , int ) ;   /** Get next arg  **/
00109       if( flag_arg == ADN_none ) break ;        /** No more args! **/
00110 #if 0
00111 fprintf(stderr,"EDIT_dset_items: iarg=%d flag_arg=%d\n",iarg,flag_arg) ;
00112 #endif
00113       iarg++ ;
00114 
00115       switch( flag_arg ){
00116 
00117          default:{
00118             int iv ;
00119             char str[128] ;
00120 
00121             /** 30 Nov 1997: check for special cases **/
00122 
00123             iv = flag_arg - ADN_brick_label_one ;
00124             if( iv >= 0 && iv < ADN_ONE_STEP ){
00125                brick_label_one_iv  = iv ;
00126                brick_label_one     = va_arg( vararg_ptr , char * ) ;
00127                if( brick_label_one != NULL ) new_brick_label_one = 1 ;
00128                break ; /* exit switch */
00129             }
00130 
00131             iv = flag_arg - ADN_brick_fac_one ;
00132             if( iv >= 0 && iv < ADN_ONE_STEP ){
00133                brick_fac_one_iv  = iv ;
00134                brick_fac_one     = va_arg( vararg_ptr , double ) ;
00135                new_brick_fac_one = 1 ;
00136                break ; /* exit switch */
00137             }
00138 
00139             iv = flag_arg - ADN_brick_stataux_one ;
00140             if( iv >= 0 && iv < ADN_ONE_STEP ){
00141                brick_stataux_one_iv = iv ;
00142                brick_stataux_one    = va_arg( vararg_ptr , float * ) ;
00143                if( brick_stataux_one != NULL ) new_brick_stataux_one = 1 ;
00144                break ; /* exit switch */
00145             }
00146 
00147             iv = flag_arg - ADN_brick_keywords_replace_one ;
00148             if( iv >= 0 && iv < ADN_ONE_STEP ){
00149                brick_keywords_one_iv  = iv ;
00150                brick_keywords_one     = va_arg( vararg_ptr , char * ) ;
00151                new_brick_keywords_one = 1 ;
00152                break ; /* exit switch */
00153             }
00154 
00155             iv = flag_arg - ADN_brick_keywords_append_one ;
00156             if( iv >= 0 && iv < ADN_ONE_STEP ){
00157                brick_keywords_one_iv  = iv ;
00158                brick_keywords_one     = va_arg( vararg_ptr , char * ) ;
00159                new_brick_keywords_one = 2 ;
00160                break ; /* exit switch */
00161             }
00162 
00163             /** not a special case? error! **/
00164 
00165             sprintf(str,"illegal opcode = %d at arg #%d",flag_arg,iarg) ;
00166             EDERR(str) ; if( errnum > 9 ) RETURN(errnum) ;
00167             dummy = va_arg( vararg_ptr , void * ) ;  /* skip next arg */
00168          }
00169          break ;
00170 
00171          /** these two commands affect the disk file names **/
00172 
00173          case ADN_prefix:  /* processed later */
00174             prefix = va_arg( vararg_ptr , char * ) ;
00175             if( prefix != NULL ) new_prefix = 1 ;
00176             else EDERR("illegal new prefix") ;
00177          break ;
00178 
00179          case ADN_directory_name:  /* processed later */
00180             directory_name = va_arg( vararg_ptr , char * ) ;
00181             if( directory_name != NULL ) new_directory_name = 1 ;
00182             else EDERR("illegal new directory_name") ;
00183          break ;
00184 
00185          /** change the memory allocation type (mmap or malloc) **/
00186 
00187          case ADN_malloc_type:  /* processed now */
00188             malloc_type = va_arg( vararg_ptr , int ) ;
00189             if( ISVALID_MEM_CODE(malloc_type) ){
00190                new_malloc_type = 1 ;
00191                THD_force_malloc_type( dset->dblk , malloc_type ) ;
00192             }
00193             else EDERR("illegal new malloc_type") ;
00194          break ;
00195 
00196          /** these commands affect the data in the bricks **/
00197 
00198          case ADN_brick_fac:  /* processed later */
00199             brick_fac = va_arg( vararg_ptr , float * ) ;
00200             new_brick_fac = 1 ;
00201          break ;
00202 
00203          case ADN_datum_all:  /* processed later */
00204             datum_all = va_arg( vararg_ptr , int ) ;
00205             if( AFNI_GOOD_DTYPE(datum_all) ) new_datum_all = 1 ;
00206             else EDERR("illegal new datum_all") ;
00207          break ;
00208 
00209          case ADN_datum_array:  /* processed later */
00210             datum_array = va_arg( vararg_ptr , int * ) ;
00211             if( datum_array != NULL ) new_datum_array = 1 ;
00212             else EDERR("illegal new datum_array") ;
00213          break ;
00214 
00215          case ADN_nvals:  /* processed later */
00216             nvals = va_arg( vararg_ptr , int ) ;
00217             if( nvals > 0 ) new_nvals = 1 ;
00218             else EDERR("illegal new nvals") ;
00219          break ;
00220 
00221          /** these commands affect the spatial axes,
00222              which may also influence size of the data bricks **/
00223 
00224          case ADN_nxyz:  /* processed later */
00225             nxyz = va_arg( vararg_ptr , THD_ivec3 ) ;
00226             if( nxyz.ijk[0] >= 1 && nxyz.ijk[1] >= 1 && nxyz.ijk[2] >= 1 )
00227                new_nxyz = 1 ;
00228                else EDERR("illegal new nxyz") ;
00229          break ;
00230 
00231          case ADN_xyzdel:  /* processed later */
00232             xyzdel = va_arg( vararg_ptr , THD_fvec3 ) ;
00233             if( xyzdel.xyz[0]!=0.0 && xyzdel.xyz[1]!=0.0 && xyzdel.xyz[2]!=0.0 )
00234                new_xyzdel = 1 ;
00235                else EDERR("illegal new xyzdel") ;
00236          break ;
00237 
00238          case ADN_xyzorg:  /* processed later */
00239             xyzorg = va_arg( vararg_ptr , THD_fvec3 ) ;
00240             new_xyzorg = 1 ;
00241          break ;
00242 
00243          case ADN_xyzorient:  /* processed later */
00244             xyzorient = va_arg( vararg_ptr , THD_ivec3 ) ;
00245             if( xyzorient.ijk[0] >= FIRST_ORIENT_TYPE &&
00246                 xyzorient.ijk[0] <= LAST_ORIENT_TYPE  &&
00247                 xyzorient.ijk[1] >= FIRST_ORIENT_TYPE &&
00248                 xyzorient.ijk[1] <= LAST_ORIENT_TYPE  &&
00249                 xyzorient.ijk[2] >= FIRST_ORIENT_TYPE &&
00250                 xyzorient.ijk[2] <= LAST_ORIENT_TYPE  &&
00251                   ORIENT_xyzint[xyzorient.ijk[0]]
00252                 + ORIENT_xyzint[xyzorient.ijk[1]]
00253                 + ORIENT_xyzint[xyzorient.ijk[2]] == 6  )
00254 
00255                new_xyzorient = 1 ;
00256             else EDERR("illegal new xyzorient") ;
00257          break ;
00258 
00259          case ADN_to_dicomm:   /* illegal at this time */
00260             EDERR("to_dicomm not implemented") ;
00261          break ;
00262 
00263          /** these commands affect the time axis of the dataset **/
00264 
00265          case ADN_ntt:  /* processed later */
00266             ntt = va_arg( vararg_ptr , int ) ;
00267             if( ntt >= 0 ) new_ntt = 1 ;
00268             else EDERR("illegal new taxis ntt") ;
00269          break ;
00270 
00271          case ADN_tunits:  /* processed later */
00272             tunits = va_arg( vararg_ptr , int ) ;
00273             if( tunits >= 0 ) new_tunits = 1 ;
00274             else EDERR("illegal new taxis tunits") ;
00275          break ;
00276 
00277          case ADN_nsl:  /* processed later */
00278             nsl = va_arg( vararg_ptr , int ) ;
00279             if( nsl >= 0 ) new_nsl = 1 ;
00280             else EDERR("illegal new taxis nsl") ;
00281          break ;
00282 
00283          case ADN_ttorg:  /* processed later */
00284             ttorg = (float) va_arg( vararg_ptr , double ) ;
00285             new_ttorg = 1 ;
00286          break ;
00287 
00288          case ADN_ttdel:  /* processed later */
00289             ttdel = (float) va_arg( vararg_ptr , double ) ;
00290             new_ttdel = 1 ;
00291          break ;
00292 
00293          case ADN_ttdur:  /* processed later */
00294             ttdur = (float) va_arg( vararg_ptr , double ) ;
00295             new_ttdur = 1 ;
00296          break ;
00297 
00298          case ADN_zorg_sl:  /* processed later */
00299             zorg_sl = (float) va_arg( vararg_ptr , double ) ;
00300             new_zorg_sl = 1 ;
00301          break ;
00302 
00303          case ADN_dz_sl:  /* processed later */
00304             dz_sl = (float) va_arg( vararg_ptr , double ) ;
00305             if( dz_sl != 0.0 ) new_dz_sl = 1 ;
00306             else EDERR("illegal new taxis dz_sl") ;
00307          break ;
00308 
00309          case ADN_toff_sl:  /* processed later */
00310             toff_sl = va_arg( vararg_ptr , float * ) ;
00311             new_toff_sl = 1 ;
00312          break ;
00313 
00314          /** these commands affect the interpretation of the dataset
00315              (e.g., is it functional or anatomical, which view type, ...) **/
00316 
00317          case ADN_type:  /* processed later */
00318             type = va_arg( vararg_ptr , int ) ;
00319             if( type >= FIRST_3DIM_TYPE && type <= LAST_3DIM_TYPE )
00320                new_type = 1 ;
00321             else EDERR("illegal new type") ;
00322          break ;
00323 
00324          case ADN_view_type:  /* processed later */
00325             view_type = va_arg( vararg_ptr , int ) ;
00326             if( view_type >= FIRST_VIEW_TYPE && view_type <= LAST_VIEW_TYPE )
00327                new_view_type = 1 ;
00328             else EDERR("illegal new view_type") ;
00329          break ;
00330 
00331          case ADN_func_type:  /* processed later */
00332             func_type = va_arg( vararg_ptr , int ) ;
00333             if( func_type >= 0 ) new_func_type = 1 ;
00334             else EDERR("illegal new func_type") ;
00335          break ;
00336 
00337          /** auxiliary statistical data, for interpretation of functions **/
00338 
00339          case ADN_stat_aux:  /* processed now */
00340             stat_aux = va_arg( vararg_ptr , float * ) ;
00341             if( stat_aux != NULL ){
00342                new_stat_aux = 1 ;
00343                for( ii=0 ; ii < MAX_STAT_AUX ; ii++ )
00344                   dset->stat_aux[ii] = stat_aux[ii] ;
00345             } else EDERR("illegal new stat_aux") ;
00346          break ;
00347 
00348          /** dataset keywords **/
00349 
00350          case ADN_keywords_replace: /* processed now */
00351             keywords = va_arg( vararg_ptr , char * ) ;
00352             new_keywords = 1 ;
00353             THD_store_dataset_keywords( dset , keywords ) ;
00354          break ;
00355 
00356          case ADN_keywords_append: /* processed now */
00357             keywords = va_arg( vararg_ptr , char * ) ;
00358             new_keywords = 1 ;
00359             THD_append_dataset_keywords( dset , keywords ) ;
00360          break ;
00361 
00362          /** various labeling options **/
00363 
00364          case ADN_label1:  /* processed now */
00365             label1 = va_arg( vararg_ptr , char * ) ;
00366             if( label1 != NULL ){
00367                MCW_strncpy( dset->label1 , label1 , THD_MAX_LABEL ) ;
00368                new_label1 = 1 ;
00369             }
00370             else EDERR("illegal new label1") ;
00371          break ;
00372 
00373          case ADN_label2:  /* processed now */
00374             label2 = va_arg( vararg_ptr , char * ) ;
00375             if( label2 != NULL ){
00376                MCW_strncpy( dset->label2 , label2 , THD_MAX_LABEL ) ;
00377                new_label2 = 1 ;
00378             }
00379             else EDERR("illegal new label2") ;
00380          break ;
00381 
00382          case ADN_self_name:  /* processed now */
00383             self_name = va_arg( vararg_ptr , char * ) ;
00384             if( self_name != NULL ){
00385                MCW_strncpy( dset->self_name , self_name , THD_MAX_NAME ) ;
00386                new_self_name = 1 ;
00387             }
00388             else EDERR("illegal new self_name") ;
00389          break ;
00390 
00391          /** relationships to other datasets **/
00392 
00393          case ADN_warp_parent:  /* processed now */
00394             warp_parent = va_arg( vararg_ptr , THD_3dim_dataset * ) ;
00395             if( ISVALID_3DIM_DATASET(warp_parent) ){
00396                new_warp_parent = 1 ;
00397                dset->warp_parent = warp_parent ;
00398                MCW_strncpy(dset->warp_parent_name,warp_parent->self_name,THD_MAX_NAME) ;
00399                dset->warp_parent_idcode = warp_parent->idcode ;
00400             }
00401             else EDERR("illegal new warp_parent") ;
00402          break ;
00403 
00404          case ADN_warp:  /* processed now */
00405             warp = va_arg( vararg_ptr , THD_warp * ) ;
00406             if( ISVALID_WARP(warp) ){
00407                new_warp = 1 ;
00408                if( dset->warp == NULL ) dset->warp = myXtNew(THD_warp) ;
00409                *(dset->warp) =* warp ;
00410             } else EDERR("illegal new warp") ;
00411          break ;
00412 
00413          case ADN_anat_parent:  /* processed now */
00414             anat_parent = va_arg( vararg_ptr , THD_3dim_dataset * ) ;
00415             if( ISVALID_3DIM_DATASET(anat_parent) ){
00416                new_anat_parent = 1 ;
00417                dset->anat_parent = anat_parent ;
00418                MCW_strncpy(dset->anat_parent_name,anat_parent->self_name,THD_MAX_NAME) ;
00419                dset->anat_parent_idcode = anat_parent->idcode ;
00420             }
00421             else EDERR("illegal new anat_parent") ;
00422          break ;
00423 
00424          case ADN_anatpar_idcode:{ /* processed now [13 Dec 1999] */
00425             MCW_idcode * idc ;
00426             idc = va_arg( vararg_ptr , MCW_idcode * ) ;
00427             if( idc != NULL )
00428                dset->anat_parent_idcode = *idc ;
00429             else
00430                EDERR("illegal new anatpar_idcode") ;
00431          }
00432          break ;
00433 
00434       }  /*- end of switch on flag_arg -*/
00435 
00436       iarg++ ;
00437    } while( 1 ) ;  /* end of loop over arguments */
00438    va_end( vararg_ptr ) ;
00439    if( errnum > 0 ) RETURN(errnum) ;
00440 
00441    /**** carry out edits that were flagged above ****/
00442 
00443    /**---------- Need to reset the disk filename? ------------**/
00444    /** 22 Nov 2002: remove +orig etc. from prefix, if present **/
00445 
00446    if( new_prefix || new_directory_name || new_view_type ){
00447       char *nprefix = THD_deplus_prefix( prefix ) ;
00448 
00449       THD_init_diskptr_names( dset->dblk->diskptr ,
00450                               directory_name , NULL ,
00451                               nprefix , view_type , True ) ;
00452 
00453       if( DSET_IS_1D(dset) || DSET_IS_3D(dset) ){         /* 21 Mar 2003 */
00454         char *fname = dset->dblk->diskptr->brick_name ;
00455         int  ll = strlen(fname) ;
00456         fname[ll-10] = '\0' ;
00457         if( DSET_IS_1D(dset) || (DSET_NY(dset)==1 && DSET_NZ(dset)==1) )
00458           strcat(fname,".1D");
00459         else
00460           strcat(fname,".3D");
00461       }
00462 
00463       /** output of NIfTI-1.1 dataset: 06 May 2005 **/
00464 
00465       if( nprefix != NULL && ( STRING_HAS_SUFFIX(nprefix,".nii") ||
00466                                STRING_HAS_SUFFIX(nprefix,".nii.gz") ) ){
00467         char *fname = dset->dblk->diskptr->brick_name ;
00468         int  ll = strlen(fname) ;
00469         fname[ll-10] = '\0' ;
00470         if( STRING_HAS_SUFFIX(nprefix,".nii") ) strcat(fname,".nii") ;
00471         else                                    strcat(fname,".nii.gz") ;
00472       }
00473 
00474       if( nprefix != NULL ) free(nprefix) ;
00475    }
00476 
00477    /**----------- Need to reconfigure the spatial axes? -----------**/
00478    /**    Most of this code is from routine THD_3dim_from_block    **/
00479 
00480    redo_daxes = ( new_nxyz || new_xyzorg || new_xyzdel || new_xyzorient ) ;
00481 
00482    if( redo_daxes ){
00483       THD_dataxes * daxes = dset->daxes ;
00484       THD_diskptr * dkptr = dset->dblk->diskptr ;
00485 
00486       /** copy new stuff into the daxes structure **/
00487 
00488       if( new_nxyz ){
00489          daxes->nxx  = dkptr->dimsizes[0] = nxyz.ijk[0] ;
00490          daxes->nyy  = dkptr->dimsizes[1] = nxyz.ijk[1] ;
00491          daxes->nzz  = dkptr->dimsizes[2] = nxyz.ijk[2] ;
00492       }
00493 
00494       if( new_xyzorg ){
00495          daxes->xxorg = xyzorg.xyz[0] ;
00496          daxes->yyorg = xyzorg.xyz[1] ;
00497          daxes->zzorg = xyzorg.xyz[2] ;
00498       }
00499 
00500       if( new_xyzdel ){
00501          daxes->xxdel = xyzdel.xyz[0] ;
00502          daxes->yydel = xyzdel.xyz[1] ;
00503          daxes->zzdel = xyzdel.xyz[2] ;
00504       }
00505 
00506       if( new_xyzorient ){
00507          daxes->xxorient = xyzorient.ijk[0] ;
00508          daxes->yyorient = xyzorient.ijk[1] ;
00509          daxes->zzorient = xyzorient.ijk[2] ;
00510       }
00511 
00512       /*---------------------------------------*/
00513       /*-- set bounding box for this dataset --*/
00514       /*---------------------------------------*/
00515 
00516       daxes->xxmin = daxes->xxorg ;
00517       daxes->xxmax = daxes->xxorg + (daxes->nxx-1) * daxes->xxdel ;
00518       if( daxes->xxmin > daxes->xxmax ){
00519          float temp   = daxes->xxmin ;
00520          daxes->xxmin = daxes->xxmax ; daxes->xxmax = temp ;
00521       }
00522 
00523       daxes->yymin = daxes->yyorg ;
00524       daxes->yymax = daxes->yyorg + (daxes->nyy-1) * daxes->yydel ;
00525       if( daxes->yymin > daxes->yymax ){
00526          float temp   = daxes->yymin ;
00527          daxes->yymin = daxes->yymax ; daxes->yymax = temp ;
00528       }
00529 
00530       daxes->zzmin = daxes->zzorg ;
00531       daxes->zzmax = daxes->zzorg + (daxes->nzz-1) * daxes->zzdel ;
00532       if( daxes->zzmin > daxes->zzmax ){
00533          float temp   = daxes->zzmin ;
00534          daxes->zzmin = daxes->zzmax ; daxes->zzmax = temp ;
00535       }
00536 
00537 #ifdef EXTEND_BBOX
00538       daxes->xxmin -= 0.5 * daxes->xxdel ;  /* pushes edges back by 1/2  */
00539       daxes->xxmax += 0.5 * daxes->xxdel ;  /* voxel dimensions (the box */
00540       daxes->yymin -= 0.5 * daxes->yydel ;  /* defined above is based on */
00541       daxes->yymax += 0.5 * daxes->yydel ;  /* voxel centers, not edges) */
00542       daxes->zzmin -= 0.5 * daxes->zzdel ;
00543       daxes->zzmax += 0.5 * daxes->zzdel ;
00544 #endif
00545 
00546       /*----------------------------------------------------------------*/
00547       /*--  matrix that transforms to Dicom (left-posterior-superior) --*/
00548       /*----------------------------------------------------------------*/
00549 
00550       LOAD_ZERO_MAT(daxes->to_dicomm) ;
00551 
00552       switch( daxes->xxorient ){
00553          case ORI_R2L_TYPE:
00554          case ORI_L2R_TYPE: daxes->to_dicomm.mat[0][0] = 1.0 ; break ;
00555          case ORI_P2A_TYPE:
00556          case ORI_A2P_TYPE: daxes->to_dicomm.mat[1][0] = 1.0 ; break ;
00557          case ORI_I2S_TYPE:
00558          case ORI_S2I_TYPE: daxes->to_dicomm.mat[2][0] = 1.0 ; break ;
00559       }
00560 
00561       switch( daxes->yyorient ){
00562          case ORI_R2L_TYPE:
00563          case ORI_L2R_TYPE: daxes->to_dicomm.mat[0][1] = 1.0 ; break ;
00564          case ORI_P2A_TYPE:
00565          case ORI_A2P_TYPE: daxes->to_dicomm.mat[1][1] = 1.0 ; break ;
00566          case ORI_I2S_TYPE:
00567          case ORI_S2I_TYPE: daxes->to_dicomm.mat[2][1] = 1.0 ; break ;
00568       }
00569 
00570       switch( daxes->zzorient ){
00571          case ORI_R2L_TYPE:
00572          case ORI_L2R_TYPE: daxes->to_dicomm.mat[0][2] = 1.0 ; break ;
00573          case ORI_P2A_TYPE:
00574          case ORI_A2P_TYPE: daxes->to_dicomm.mat[1][2] = 1.0 ; break ;
00575          case ORI_I2S_TYPE:
00576          case ORI_S2I_TYPE: daxes->to_dicomm.mat[2][2] = 1.0 ; break ;
00577       }
00578    }
00579 
00580    /**---------- Need to reconfigure the sub-bricks? ----------**/
00581 
00582    if( new_datum_all && new_datum_array ){
00583        EDERR("datum_all and datum_array can't be used together") ;
00584        RETURN(errnum) ;
00585    }
00586 
00587    redo_bricks = ( new_datum_all || new_datum_array ||
00588                    new_nvals     || new_nxyz          ) ;
00589 
00590    if( redo_bricks && THD_count_databricks(dset->dblk) > 0 ){
00591       EDERR("cannot reconfigure bricks that already are full") ;
00592       RETURN(errnum) ;
00593    }
00594 
00595    if( redo_bricks ){
00596       int old_nvals = dset->dblk->nvals ;
00597 #if 0
00598 fprintf(stderr,"EDIT_dset_items: about to redo_bricks\n") ;
00599 #endif
00600       if( ! new_nvals ) nvals = old_nvals ;
00601 
00602       /** make an array of data types, if one not provided **/
00603 
00604       if( ! new_datum_array ){
00605          datum_array = (int *) XtMalloc( sizeof(int) * nvals ) ;
00606 
00607 #if 0
00608 fprintf(stderr,"EDIT_dset_items: about to make datum_array\n") ;
00609 #endif
00610          for( ii=0 ; ii < nvals ; ii++ )
00611             datum_array[ii] =  (new_datum_all)  ? datum_all
00612                              : (ii < old_nvals) ? DSET_BRICK_TYPE(dset,ii)
00613                                                 : DSET_BRICK_TYPE(dset,0) ;
00614       }                                           /* 06 Apr 2005 [rickr] */
00615 
00616       if( new_nvals ){
00617          if( dset->dblk->nvals != nvals )
00618             THD_copy_datablock_auxdata( NULL , dset->dblk ) ; /* 30 Nov 1997 */
00619 
00620          myXtFree( dset->dblk->brick_bytes ) ;
00621          myXtFree( dset->dblk->brick_fac   ) ;
00622 
00623          dset->dblk->nvals = dset->dblk->diskptr->nvals = nvals ;
00624       }
00625 
00626       THD_init_datablock_brick( dset->dblk , nvals , datum_array ) ;
00627 
00628       if( ! new_datum_array ) myXtFree(datum_array) ;
00629    }
00630 
00631    /**---------- Need to add new brick_fac values? ----------**/
00632 
00633    if( new_brick_fac ){
00634       if( brick_fac != NULL ){
00635          for( ii=0 ; ii < dset->dblk->nvals ; ii++ )
00636             dset->dblk->brick_fac[ii] = brick_fac[ii] ;
00637       } else {
00638          for( ii=0 ; ii < dset->dblk->nvals ; ii++ )
00639             dset->dblk->brick_fac[ii] = 0.0 ;
00640       }
00641    }
00642 
00643    /** 30 Nov 1997: do just one brick_fac value **/
00644 
00645    if( new_brick_fac_one ){
00646       if( brick_fac_one_iv < 0 || brick_fac_one_iv >= dset->dblk->nvals ){
00647          EDERR("illegal index for ADN_brick_fac_one") ;
00648          RETURN(errnum) ;
00649       }
00650       dset->dblk->brick_fac[ brick_fac_one_iv ] = brick_fac_one ;
00651    }
00652 
00653    /**--------- 30 Nov 1997: add a single brick label value --------**/
00654 
00655    if( new_brick_label_one ){
00656       if( brick_label_one_iv < 0 || brick_label_one_iv >= dset->dblk->nvals ){
00657          EDERR("illegal index for ADN_brick_label_one") ;
00658          RETURN(errnum) ;
00659       }
00660 
00661       THD_store_datablock_label( dset->dblk, brick_label_one_iv, brick_label_one ) ;
00662    }
00663 
00664    /*---- add a single brick keywords value ----*/
00665 
00666    if( new_brick_keywords_one ){
00667       if( brick_keywords_one_iv < 0 || brick_keywords_one_iv >= dset->dblk->nvals ){
00668          EDERR("illegal index for ADN_brick_keywords_one") ;
00669          RETURN(errnum) ;
00670       }
00671 
00672       if( new_brick_keywords_one == 1 )
00673          THD_store_datablock_keywords( dset->dblk, brick_keywords_one_iv,
00674                                                    brick_keywords_one    );
00675       else if( new_brick_keywords_one == 2 )
00676          THD_append_datablock_keywords( dset->dblk, brick_keywords_one_iv,
00677                                                     brick_keywords_one    );
00678    }
00679 
00680    /*---- Add a single brick stataux value.
00681           The input is a float array formatted like so:
00682             <statcode> <npar> <value> ... <value>
00683           where <statcode> is a FUNC_*_TYPE code
00684                 <npar>     is the number of values to follow (may be 0);
00685                              normally is FUNC_need_stat_aux[<statcode>]
00686                 <value>    is an auxiliary statistical parameter needed
00687                              for data of type <statcode>                ----*/
00688 
00689    if( new_brick_stataux_one ){
00690       int jv , npar , kv , iv ;
00691 
00692       iv = brick_stataux_one_iv ;
00693 
00694       if( iv < 0 || iv >= dset->dblk->nvals ){
00695          EDERR("illegal index for ADN_brick_stataux_one") ;
00696          RETURN(errnum) ;
00697       }
00698 
00699       jv = brick_stataux_one[0] ;  /* statcode */
00700 
00701       npar = brick_stataux_one[1] ;  /* # of values present */
00702       if( npar < 0 ){
00703          EDERR("illegal npar for ADN_brick_stataux_one") ;
00704          RETURN(errnum) ;
00705       }
00706 
00707       kv = FUNC_need_stat_aux[jv] ;  /* # of values needed */
00708       if( npar > kv ) npar = kv ;
00709 
00710       THD_store_datablock_stataux( dset->dblk ,
00711                                    iv , jv , npar , brick_stataux_one + 2 ) ;
00712    }
00713 
00714    /**---------- Need to reconfigure the time axis? ----------**/
00715 
00716    redo_taxis = ( new_ntt   || new_nsl     || new_ttorg || new_ttdel ||
00717                   new_ttdur || new_zorg_sl || new_dz_sl || new_toff_sl ) ;
00718 
00719    if( ! new_ntt ) ntt = ISVALID_TIMEAXIS(dset->taxis) ? dset->taxis->ntt : 0 ;
00720 
00721    if( ntt == 0 && dset->taxis != NULL ){
00722       myXtFree( dset->taxis->toff_sl ) ;
00723       myXtFree( dset->taxis ) ;
00724       dset->taxis = NULL ;
00725    }
00726 
00727    redo_taxis = ( redo_taxis && ntt > 0 ) ;
00728 
00729    if( (new_nsl && nsl > 0) && !new_toff_sl ){    /* if we have new slice count */
00730       EDERR("have new_nsl but not new_toff_sl") ; /* but no new slice offsets */
00731       RETURN(errnum) ;
00732    }
00733 
00734    if( redo_taxis ){
00735       THD_timeaxis * taxis = dset->taxis ;
00736 
00737       if( taxis == NULL ){
00738          taxis          = dset->taxis     = myXtNew( THD_timeaxis ) ;
00739          taxis->type    = TIMEAXIS_TYPE ;
00740          taxis->toff_sl = NULL ;
00741          taxis->nsl     = 0 ;
00742          taxis->ttorg   = taxis->ttdel = taxis->ttdur = 0.0 ;
00743          taxis->ntt     = ntt ;
00744       }
00745 
00746       if( new_ntt     ) taxis->ntt     = ntt ;
00747       if( new_ttorg   ) taxis->ttorg   = ttorg ;
00748       if( new_ttdel   ) taxis->ttdel   = ttdel ;
00749       if( new_ttdur   ) taxis->ttdur   = ttdur ;
00750       if( new_zorg_sl ) taxis->zorg_sl = zorg_sl ;
00751       if( new_dz_sl   ) taxis->dz_sl   = dz_sl ;
00752 
00753       if( new_nsl ){
00754          taxis->nsl = nsl ;
00755          if( nsl > 0 )
00756             taxis->toff_sl = (float *) XtRealloc( (char *) taxis->toff_sl ,
00757                                                   sizeof(float) * nsl      ) ;
00758          else
00759             myXtFree(taxis->toff_sl) ;
00760       }
00761 
00762       if( new_toff_sl )
00763          for( ii=0 ; ii < taxis->nsl ; ii++ ) taxis->toff_sl[ii] = toff_sl[ii] ;
00764    }
00765 
00766    if( new_tunits ){
00767       THD_timeaxis * taxis = dset->taxis ;
00768 
00769       if( taxis == NULL ){
00770          EDERR("have new_tunits but have no time axis") ;
00771          RETURN(errnum) ;
00772       }
00773 
00774       taxis->units_type = tunits ;
00775    }
00776 
00777    /**--------------- Need to redo dataset type codes? ------------**/
00778    /**  Note that changing the type codes by themselves won't fix  **/
00779    /**  nvals or other such stuff -- that must be done separately. **/
00780 
00781    if( new_type      ) dset->type      = type ;
00782    if( new_view_type ) dset->view_type = view_type ;
00783 
00784    if( new_func_type ){
00785       if( (ISANAT(dset) && func_type <= LAST_ANAT_TYPE) ||
00786           (ISFUNC(dset) && func_type <= LAST_FUNC_TYPE)   ){
00787 
00788          dset->func_type = func_type ;
00789 
00790       } else{
00791          EDERR("illegal new_func type combination") ; RETURN(errnum) ;
00792       }
00793    }
00794 
00795    /****--------------- hopefully, we are done! ---------------****/
00796 
00797    RETURN(errnum) ;
00798 }