Doxygen Source Code Documentation
3duca.c File Reference
#include <string.h>
#include "mrilib.h"
#include <stdlib.h>
#include <ctype.h>
#include "uuu.c"
Go to the source code of this file.
Functions | |
void | UC_syntax (char *msg) |
void | detrend (int n, float vec[]) |
void | normalize (int n, float vec[]) |
void | UC_read_opts (int argc, char *argv[]) |
float | UC_unusuality (int ndim, float *ref, int nvec, float **vec) |
int | main (int argc, char *argv[]) |
Variables | |
THD_3dim_dataset * | UC_dset = NULL |
MRI_IMAGE * | UC_ref = NULL |
char | UC_prefix [THD_MAX_PREFIX] = "uc" |
int | UC_be_quiet = 1 |
byte * | UC_mask = NULL |
int | UC_mask_nvox = 0 |
int | UC_mask_hits = 0 |
int | UC_nvec = 0 |
int | UC_vdim = 0 |
float ** | UC_vec = NULL |
Function Documentation
|
Definition at line 40 of file 3duca.c. References vec.
00041 { 00042 register int ii ; 00043 register float sum0 , sum1 , cf , lf ; 00044 float sum2 , det ; 00045 00046 static int nold = -1 ; /* initialization flag */ 00047 static float cf0,cf1 , lf0,lf1 ; /* to be initialized */ 00048 00049 /*** initialize coefficients for detrending ***/ 00050 00051 if( n != nold ){ 00052 nold = n ; sum0 = sum1 = sum2 = 0.0 ; 00053 for( ii=0 ; ii < n ; ii++ ){ 00054 sum0 += 1.0 ; sum1 += ii ; sum2 += ii*ii ; 00055 } 00056 det = sum0 * sum2 - sum1 * sum1 ; 00057 cf0 = sum2 / det ; /* constant factor for sum0 */ 00058 cf1 = -sum1 / det ; /* constant factor for sum1 */ 00059 lf0 = cf1 ; /* linear factor for sum0 */ 00060 lf1 = sum0 / det ; /* linear factor for sum1 */ 00061 } 00062 00063 /*** remove mean and linear trend ***/ 00064 00065 sum0 = sum1 = 0.0 ; 00066 for( ii=0 ; ii < n ; ii++ ){ 00067 sum0 += vec[ii] ; sum1 += vec[ii] * ii ; 00068 } 00069 00070 cf = cf0 * sum0 + cf1 * sum1 ; 00071 lf = lf0 * sum0 + lf1 * sum1 ; 00072 for( ii=0 ; ii < n ; ii++ ) vec[ii] -= cf + ii*lf ; 00073 } |
|
---------- Adapted from 3dZeropad.c by RWCox - 08 Aug 2001 ----------* Definition at line 319 of file 3duca.c. References argc, MRI_FLOAT_PTR, my_getenv(), MRI_IMAGE::nx, MRI_IMAGE::ny, UC_nvec, UC_read_opts(), UC_syntax(), UC_unusuality(), UC_vdim, and UC_vec.
00320 { 00321 int kk ; 00322 00323 /*-- read command line arguments --*/ 00324 00325 if( argc < 2 || strncmp(argv[1],"-help",5) == 0 ) UC_syntax(NULL) ; 00326 00327 (void) my_getenv("junk") ; 00328 00329 UC_read_opts( argc , argv ) ; 00330 00331 for( kk=0 ; kk < UC_ref->ny ; kk++ ) 00332 (void) UC_unusuality( UC_vdim, 00333 MRI_FLOAT_PTR(UC_ref) + kk*UC_ref->nx, 00334 UC_nvec, UC_vec ) ; 00335 exit(0) ; 00336 } |
|
Definition at line 79 of file 3duca.c. References detrend(), and vec.
00080 { 00081 register int ii ; 00082 register float sqsum ; 00083 00084 detrend( n , vec ) ; 00085 00086 sqsum = 0.0 ; 00087 for( ii=0 ; ii < n ; ii++ ) sqsum += vec[ii] * vec[ii] ; 00088 00089 if( sqsum < 1.e-10 ){ 00090 for( ii=0 ; ii < n ; ii++ ) vec[ii] = 0.0 ; 00091 } else { 00092 sqsum = 1.0 / sqrt(sqsum) ; 00093 for( ii=0 ; ii < n ; ii++ ) vec[ii] *= sqsum ; 00094 } 00095 } |
|
Definition at line 99 of file 3duca.c. References argc, DSET_ARRAY, DSET_BRICK_TYPE, DSET_delete, DSET_load, DSET_LOADED, DSET_NVALS, DSET_NVOX, DSET_unload, DSET_unload_one, EDIT_coerce_type(), free, ISVALID_3DIM_DATASET, MRI_IMAGE::kind, malloc, MCW_strncpy, MRI_FLOAT_PTR, mri_free(), mri_read(), mri_to_float(), mri_transpose(), normalize(), MRI_IMAGE::nx, MRI_IMAGE::ny, THD_countmask(), THD_makemask(), THD_MAX_PREFIX, THD_open_dataset(), UC_be_quiet, UC_mask, UC_mask_hits, UC_mask_nvox, UC_nvec, UC_prefix, UC_syntax(), UC_vdim, and UC_vec.
00100 { 00101 int nopt = 1 ; 00102 float val ; 00103 int kk, nxyz, mm,nn ; 00104 float * vv , * bb ; 00105 00106 while( nopt < argc && argv[nopt][0] == '-' ){ 00107 00108 /**** -verbose ****/ 00109 00110 if( strncmp(argv[nopt],"-verbose",5) == 0 ){ 00111 UC_be_quiet = 0 ; 00112 nopt++ ; continue ; 00113 } 00114 00115 /**** -ref file.1D ****/ 00116 00117 if( strncmp(argv[nopt],"-ref",4) == 0 ){ 00118 MRI_IMAGE * im ; 00119 nopt++ ; 00120 if( nopt >= argc ) UC_syntax("-ref needs an argument!") ; 00121 im = mri_read( argv[nopt] ) ; 00122 if( im == NULL ) UC_syntax("Can't read -ref file!") ; 00123 if( im->kind == MRI_float ){ 00124 UC_ref = im ; 00125 } else { 00126 UC_ref = mri_to_float(im) ; mri_free(im) ; 00127 } 00128 im = mri_transpose(UC_ref) ; mri_free(UC_ref) ; UC_ref = im ; 00129 nopt++ ; continue ; 00130 } 00131 00132 /**** -prefix prefix ****/ 00133 00134 if( strncmp(argv[nopt],"-prefix",6) == 0 ){ 00135 nopt++ ; 00136 if( nopt >= argc ) UC_syntax("-prefix needs an argument!") ; 00137 MCW_strncpy( UC_prefix , argv[nopt++] , THD_MAX_PREFIX ) ; 00138 continue ; 00139 } 00140 00141 /**** -mask mset ****/ 00142 00143 if( strncmp(argv[nopt],"-mask",5) == 0 ){ 00144 THD_3dim_dataset * mset ; int ii ; 00145 nopt++ ; 00146 if( nopt >= argc ) UC_syntax("need arguments after -mask!") ; 00147 mset = THD_open_dataset( argv[nopt] ) ; 00148 if( mset == NULL ) UC_syntax("can't open -mask dataset!") ; 00149 UC_mask = THD_makemask( mset , 0 , 1.0,0.0 ) ; 00150 UC_mask_nvox = DSET_NVOX(mset) ; 00151 DSET_delete(mset) ; 00152 if( UC_mask == NULL ) UC_syntax("can't use -mask dataset!") ; 00153 UC_mask_hits = THD_countmask( UC_mask_nvox , UC_mask ) ; 00154 if( UC_mask_hits == 0 ) UC_syntax("mask is all zeros!") ; 00155 if( !UC_be_quiet ) printf("--- %d voxels in mask\n",UC_mask_hits) ; 00156 nopt++ ; continue ; 00157 } 00158 00159 /**** unknown switch ****/ 00160 00161 fprintf(stderr,"\n*** unrecognized option %s\n",argv[nopt]) ; 00162 exit(1) ; 00163 00164 } /* end of loop over options */ 00165 00166 /*--- a simple consistency check ---*/ 00167 00168 /*--- last input is dataset name ---*/ 00169 00170 if( nopt >= argc ) UC_syntax("no input dataset name?") ; 00171 00172 UC_dset = THD_open_dataset( argv[nopt] ) ; 00173 if( !ISVALID_3DIM_DATASET(UC_dset) ){ 00174 fprintf(stderr,"\n*** can't open dataset file %s\n",argv[nopt]) ; 00175 exit(1) ; 00176 } 00177 00178 nxyz = DSET_NVOX(UC_dset) ; 00179 if( UC_mask != NULL && nxyz != UC_mask_nvox ) 00180 UC_syntax("mask and input dataset size mismatch!") ; 00181 00182 /*--- load vectors ---*/ 00183 00184 UC_nvec = (UC_mask_hits > 0) ? UC_mask_hits : nxyz ; 00185 UC_vdim = DSET_NVALS(UC_dset) ; 00186 if( UC_vdim < 4 ) 00187 UC_syntax("input dataset needs at least 4 sub-bricks!") ; 00188 00189 if( UC_ref == NULL || UC_ref->nx < UC_vdim ) 00190 UC_syntax("input ref not long enough for input dataset!") ; 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 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 , DSET_BRICK_TYPE(UC_dset,mm) , 00209 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 for( kk=0 ; kk < UC_ref->ny ; kk++ ) 00231 normalize( UC_vdim , MRI_FLOAT_PTR(UC_ref) + kk*UC_ref->nx ) ; 00232 00233 return ; 00234 } |
|
Definition at line 295 of file 3duca.c. References MASTER_SHORTHELP_STRING.
00296 { 00297 if( msg != NULL ){ fprintf(stderr,"\n*** %s\n",msg) ; exit(1) ; } 00298 00299 printf( 00300 "Unusual Component Analysis of 3D Datasets\n" 00301 "Usage: 3duca [options] dataset ...\n" 00302 "\n" 00303 "The input dataset may have a sub-brick selector list.\n" 00304 "Otherwise, all sub-bricks from a dataset will be used.\n" 00305 "\n" 00306 "OPTIONS:\n" 00307 " -prefix pname \n" 00308 " -verbose\n" 00309 " -mask mset \n" 00310 " -ref file.1D\n" 00311 printf("\n" MASTER_SHORTHELP_STRING ) ; 00312 ) ; 00313 00314 exit(0) ; 00315 } |
|
Definition at line 242 of file 3duca.c. References free, getenv(), malloc, qginv(), ref, set_unusuality_tail(), strtod(), unusuality(), and vec.
00243 { 00244 register int ii , kk ; 00245 float psum,msum , * vv , val ; 00246 float zmid , zmed , zsig , zplus,zminus , uval ; 00247 00248 static int nvold=-1 ; 00249 static float * zval =NULL , *aval=NULL ; 00250 static float pstar , zstar ; 00251 00252 if( ndim < 4 || nvec < 4 || ref == NULL || vec == NULL ) return 0.0 ; 00253 00254 /* initialize if number of vectors has changed */ 00255 00256 if( nvold != nvec ){ 00257 if( zval != NULL ) free(zval) ; 00258 if( aval != NULL ) free(aval) ; 00259 zval = (float *) malloc(sizeof(float)*nvec) ; 00260 aval = (float *) malloc(sizeof(float)*nvec) ; 00261 nvold = nvec ; 00262 pstar = 10.0 / nvec ; 00263 zstar = qginv(0.5*pstar) ; 00264 } 00265 00266 /* compute dot products */ 00267 00268 for( kk=0 ; kk < nvec ; kk++ ){ 00269 psum = 0.0 ; vv = vec[kk] ; 00270 for( ii=0 ; ii < ndim ; ii++ ) psum += ref[ii] * vv[ii] ; 00271 zval[kk] = psum ; 00272 } 00273 00274 { char * ps = getenv("PTAIL") ; 00275 float pp=0.0 ; 00276 if( ps != NULL ) pp = strtod(ps,NULL) ; 00277 set_unusuality_tail(pp) ; 00278 } 00279 00280 psum = unusuality( nvec, zval ) ; 00281 00282 for( kk=0 ; kk < nvec ; kk++ ) zval[kk] = -zval[kk] ; 00283 00284 msum = unusuality( nvec, zval ) ; 00285 00286 uval = psum - msum ; 00287 00288 printf("psum=%.1f msum=%.1f total=%.1f uval=%.1f\n", 00289 psum,msum,psum+msum,uval) ; 00290 return (psum+msum) ; 00291 } |
Variable Documentation
|
Definition at line 22 of file 3duca.c. Referenced by UC_read_opts(). |
|
inputs * |
|
Definition at line 24 of file 3duca.c. Referenced by UC_read_opts(). |
|
Definition at line 26 of file 3duca.c. Referenced by UC_read_opts(). |
|
Definition at line 25 of file 3duca.c. Referenced by UC_read_opts(). |
|
Definition at line 28 of file 3duca.c. Referenced by main(), and UC_read_opts(). |
|
Definition at line 20 of file 3duca.c. Referenced by UC_read_opts(). |
|
|
|
Definition at line 29 of file 3duca.c. Referenced by main(), and UC_read_opts(). |
|
Definition at line 31 of file 3duca.c. Referenced by main(), and UC_read_opts(). |