00046 {
00047    register int ii ;
00048    register float sum0 , sum1 , cf , lf ;
00049    float sum2 , det ;
00050 
00051    static int nold = -1 ;             /* initialization flag */
00052    static float cf0,cf1 , lf0,lf1 ;   /* to be initialized */
00053 
00054  /*** initialize coefficients for detrending ***/
00055 
00056    if( n != nold ){
00057       nold = n ; sum0 = sum1 = sum2 = 0.0 ;
00058       for( ii=0 ; ii < n ; ii++ ){
00059          sum0 += 1.0 ; sum1 += ii ; sum2 += ii*ii ;
00060       }
00061       det = sum0 * sum2 - sum1 * sum1 ;
00062       cf0 =  sum2 / det ;     /* constant factor for sum0 */
00063       cf1 = -sum1 / det ;     /* constant factor for sum1 */
00064       lf0 = cf1 ;             /* linear factor for sum0 */
00065       lf1 =  sum0 / det ;     /* linear factor for sum1 */
00066    }
00067 
00068  /*** remove mean and linear trend ***/
00069 
00070    sum0 = sum1 = 0.0 ;
00071    for( ii=0 ; ii < n ; ii++ ){
00072       sum0 += vec[ii] ; sum1 += vec[ii] * ii ;
00073    }
00074 
00075    cf = cf0 * sum0 + cf1 * sum1 ;
00076    lf = lf0 * sum0 + lf1 * sum1 ;
00077    for( ii=0 ; ii < n ; ii++ ) vec[ii] -= cf + ii*lf ;
00078 }

00293 {
00294    int kk , nvox , ii ;
00295    THD_3dim_dataset * oset ;
00296    float * far ;
00297 
00298    /*-- read command line arguments --*/
00299 
00300    if( argc < 2 || strncmp(argv[1],"-help",5) == 0 ) UC_syntax(NULL) ;
00301 
00302    (void) my_getenv("junk") ;
00303 
00304    UC_read_opts( argc , argv ) ;
00305    set_unusuality_tail( UC_ptail ) ;
00306 
00307    oset = EDIT_empty_copy( UC_dset ) ;
00308    EDIT_dset_items( oset ,
00309                        ADN_prefix      , UC_prefix ,
00310                        ADN_ntt         , 0 ,
00311                        ADN_nvals       , 1 ,
00312                        ADN_datum_all   , MRI_float ,
00313                        ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
00314                     ADN_none ) ;
00315 
00316    nvox = DSET_NVOX(oset) ;
00317    far = (float *) malloc( sizeof(float) * nvox ) ;
00318    for( kk=0 ; kk < nvox ; kk++ ) far[kk] = 0.0 ;
00319    EDIT_substitute_brick( oset , 0 , MRI_float , far ) ;
00320 
00321    if( !UC_be_quiet ){ printf("--- computing u") ; fflush(stdout) ; }
00322 
00323    for( kk=0 ; kk < UC_nvec ; kk++ ){
00324       ii = (UC_iv == NULL) ? kk : UC_iv[kk] ;
00325       far[ii] = UC_unusuality( UC_vdim, UC_vec[kk] , UC_nvec, UC_vec ) ;
00326       if( !UC_be_quiet && kk%1000==999 ){ printf(".");fflush(stdout); }
00327    }
00328    if( !UC_be_quiet ) printf("\n--- writing output\n") ;
00329 
00330    DSET_write(oset) ;
00331    exit(0) ;
00332 }

00085 {
00086    register int ii ;
00087    register float sqsum ;
00088 
00089    detrend( n , vec ) ;
00090 
00091    sqsum = 0.0 ;
00092    for( ii=0 ; ii < n ; ii++ ) sqsum += vec[ii] * vec[ii] ;
00093 
00094    if( sqsum < 1.e-10 ){
00095       for( ii=0 ; ii < n ; ii++ ) vec[ii] = 0.0 ;
00096    } else {
00097       sqsum = 1.0 / sqrt(sqsum) ;
00098       for( ii=0 ; ii < n ; ii++ ) vec[ii] *= sqsum ;
00099    }
00100 }

00105 {
00106    int nopt = 1 ;
00107    float val ;
00108    int  kk, nxyz, mm,nn ;
00109    float * vv , * bb ;
00110 
00111    while( nopt < argc && argv[nopt][0] == '-' ){
00112 
00113       /**** -verbose ****/
00114 
00115       if( strncmp(argv[nopt],"-verbose",5) == 0 ){
00116          UC_be_quiet = 0 ;
00117          nopt++ ; continue ;
00118       }
00119 
00120       /**** -prefix prefix ****/
00121 
00122       if( strncmp(argv[nopt],"-prefix",6) == 0 ){
00123          nopt++ ;
00124          if( nopt >= argc ) UC_syntax("-prefix needs an argument!") ;
00125          MCW_strncpy( UC_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
00126          continue ;
00127       }
00128 
00129       /**** -mask mset ****/
00130 
00131       if( strncmp(argv[nopt],"-mask",5) == 0 ){
00132          THD_3dim_dataset * mset ; int ii,nn ;
00133          nopt++ ;
00134          if( nopt >= argc ) UC_syntax("need arguments after -mask!") ;
00135          mset = THD_open_dataset( argv[nopt] ) ;
00136          if( mset == NULL ) UC_syntax("can't open -mask dataset!") ;
00137          UC_mask = THD_makemask( mset , 0 , 1.0,0.0 ) ;
00138          UC_mask_nvox = DSET_NVOX(mset) ;
00139          DSET_delete(mset) ;
00140          if( UC_mask == NULL ) UC_syntax("can't use -mask dataset!") ;
00141          UC_mask_hits = THD_countmask( UC_mask_nvox , UC_mask ) ;
00142          if( UC_mask_hits == 0 ) UC_syntax("mask is all zeros!") ;
00143          if( !UC_be_quiet ) printf("--- %d voxels in mask\n",UC_mask_hits) ;
00144 
00145          UC_iv = (int *) malloc( sizeof(int) * UC_mask_hits ) ;
00146          for( nn=ii=0 ; ii < UC_mask_nvox ; ii++ )
00147             if( UC_mask[ii] ) UC_iv[nn++] = ii ;
00148 
00149          nopt++ ; continue ;
00150       }
00151 
00152       /**** -ptail p ****/
00153 
00154       if( strcmp(argv[nopt],"-ptail") == 0 ){
00155          if( ++nopt >= argc ) UC_syntax("-ptail needs an argument!") ;
00156          UC_ptail = strtod( argv[nopt] , NULL ) ;
00157          if( UC_ptail <= 0.0 || UC_ptail >= 0.499 )
00158             UC_syntax("value after -ptail is illegal!") ;
00159          nopt++ ; continue ;
00160       }
00161 
00162       /**** unknown switch ****/
00163 
00164       fprintf(stderr,"\n*** unrecognized option %s\n",argv[nopt]) ;
00165       exit(1) ;
00166 
00167    }  /* end of loop over options */
00168 
00169    /*--- a simple consistency check ---*/
00170 
00171    /*--- last input is dataset name ---*/
00172 
00173    if( nopt >= argc ) UC_syntax("no input dataset name?") ;
00174 
00175    UC_dset = THD_open_dataset( argv[nopt] ) ;
00176    if( !ISVALID_3DIM_DATASET(UC_dset) ){
00177       fprintf(stderr,"\n*** can't open dataset file %s\n",argv[nopt]) ;
00178       exit(1) ;
00179    }
00180 
00181    nxyz = DSET_NVOX(UC_dset) ;
00182    if( UC_mask != NULL && nxyz != UC_mask_nvox )
00183       UC_syntax("mask and input dataset size mismatch!") ;
00184 
00185    /*--- load vectors ---*/
00186 
00187    UC_nvec = (UC_mask_hits > 0) ? UC_mask_hits : nxyz ;
00188    UC_vdim = DSET_NVALS(UC_dset) ;
00189    if( UC_vdim < 4 )
00190       UC_syntax("input dataset needs at least 4 sub-bricks!") ;
00191 
00192    vv     = (float *) malloc( sizeof(float) * UC_nvec * UC_vdim ) ;
00193    UC_vec = (float **) malloc( sizeof(float *) * UC_nvec ) ;
00194    for( kk=0 ; kk < UC_nvec ; kk++ ) UC_vec[kk] = vv + (kk*UC_vdim) ;
00195 
00196    if( !UC_be_quiet ) printf("--- reading input dataset\n") ;
00197    DSET_load(UC_dset) ;
00198    if( ! DSET_LOADED(UC_dset) )
00199       UC_syntax("Can't load input dataset bricks!") ;
00200 
00201    /* copy brick data into float storage */
00202 
00203    if( !UC_be_quiet ) printf("--- loading vectors\n") ;
00204 
00205    bb = (float *) malloc( sizeof(float) * nxyz ) ;
00206    for( mm=0 ; mm < UC_vdim ; mm++ ){
00207 
00208       EDIT_coerce_type( nxyz ,
00209                         DSET_BRICK_TYPE(UC_dset,mm) , DSET_ARRAY(UC_dset,mm) ,
00210                         MRI_float , bb ) ;
00211 
00212       DSET_unload_one( UC_dset , mm ) ;
00213 
00214       if( UC_mask == NULL ){
00215          for( kk=0 ; kk < nxyz ; kk++ ) UC_vec[kk][mm] = bb[kk] ;
00216       } else {
00217          for( nn=kk=0 ; kk < nxyz ; kk++ )
00218             if( UC_mask[kk] ) UC_vec[nn++][mm] = bb[kk] ;
00219       }
00220    }
00221    free(bb) ; DSET_unload( UC_dset ) ;
00222 
00223    /* detrend and normalize vectors */
00224 
00225    if( !UC_be_quiet ) printf("--- normalizing vectors\n") ;
00226 
00227    for( kk=0 ; kk < UC_nvec ; kk++ )
00228       normalize( UC_vdim , UC_vec[kk] ) ;
00229 
00230    return ;
00231 }

00271 {
00272    if( msg != NULL ){ fprintf(stderr,"\n*** %s\n",msg) ; exit(1) ; }
00273 
00274    printf(
00275     "Usage: 3duuu [options] dataset ...\n"
00276     "\n"
00277     "The input dataset may have a sub-brick selector list.\n"
00278     "Otherwise, all sub-bricks from a dataset will be used.\n"
00279     "\n"
00280     "OPTIONS:\n"
00281     "  -prefix pname \n"
00282     "  -verbose\n"
00283     "  -mask mset\n"
00284     "  -ptail p\n"
00285    ) ;
00286 
00287    exit(0) ;
00288 }

00238 {
00239    register int ii , kk ;
00240    register float psum , * vv ;
00241 
00242    static int     nvold=-1   ;
00243    static float * zval =NULL ;
00244 
00245    if( ndim < 4 || nvec < 4 || ref == NULL || vec == NULL ) return 0.0 ;
00246 
00247    /* initialize if number of vectors has changed */
00248 
00249    if( nvold != nvec ){
00250       if( zval != NULL ) free(zval) ;
00251       zval = (float *) malloc(sizeof(float)*nvec) ;
00252       nvold = nvec ;
00253    }
00254 
00255    /* compute dot products */
00256 
00257    for( kk=0 ; kk < nvec ; kk++ ){
00258       psum = 0.0 ; vv = vec[kk] ;
00259       for( ii=0 ; ii < ndim ; ii++ ) psum += ref[ii] * vv[ii] ;
00260       zval[kk] = psum ;
00261    }
00262 
00263    psum = unusuality( nvec, zval ) ;
00264 
00265    return psum ;
00266 }