Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
3duuu.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007 #include <string.h>
00008 #include "mrilib.h"
00009 #include <stdlib.h>
00010 #include <ctype.h>
00011
00012
00013
00014
00015
00016 static THD_3dim_dataset * UC_dset = NULL ;
00017
00018 static char UC_prefix[THD_MAX_PREFIX] = "uuu" ;
00019
00020 static int UC_be_quiet = 1 ;
00021
00022 static byte * UC_mask = NULL ;
00023 static int UC_mask_nvox = 0 ;
00024 static int UC_mask_hits = 0 ;
00025
00026 static int UC_nvec = 0 ;
00027 static int UC_vdim = 0 ;
00028
00029 static float ** UC_vec = NULL ;
00030
00031
00032
00033 static int * UC_iv = NULL ;
00034
00035 static float UC_ptail = 0.0001 ;
00036
00037 void UC_syntax(char * msg) ;
00038
00039 #include "uuu.c"
00040
00041
00042
00043
00044
00045 void detrend( int n , float vec[] )
00046 {
00047 register int ii ;
00048 register float sum0 , sum1 , cf , lf ;
00049 float sum2 , det ;
00050
00051 static int nold = -1 ;
00052 static float cf0,cf1 , lf0,lf1 ;
00053
00054
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 ;
00063 cf1 = -sum1 / det ;
00064 lf0 = cf1 ;
00065 lf1 = sum0 / det ;
00066 }
00067
00068
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 }
00079
00080
00081
00082
00083
00084 void normalize( int n , float vec[] )
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 }
00101
00102
00103
00104 void UC_read_opts( int argc , char * argv[] )
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
00114
00115 if( strncmp(argv[nopt],"-verbose",5) == 0 ){
00116 UC_be_quiet = 0 ;
00117 nopt++ ; continue ;
00118 }
00119
00120
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
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
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
00163
00164 fprintf(stderr,"\n*** unrecognized option %s\n",argv[nopt]) ;
00165 exit(1) ;
00166
00167 }
00168
00169
00170
00171
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
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
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
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 }
00232
00233
00234
00235
00236
00237 float UC_unusuality( int ndim , float * ref , int nvec , float ** vec )
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
00248
00249 if( nvold != nvec ){
00250 if( zval != NULL ) free(zval) ;
00251 zval = (float *) malloc(sizeof(float)*nvec) ;
00252 nvold = nvec ;
00253 }
00254
00255
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 }
00267
00268
00269
00270 void UC_syntax(char * msg)
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 }
00289
00290
00291
00292 int main( int argc , char * argv[] )
00293 {
00294 int kk , nvox , ii ;
00295 THD_3dim_dataset * oset ;
00296 float * far ;
00297
00298
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 }