00001
00002
00003
00004
00005
00006
00007 #undef MAIN
00008 #include "afni_pcor.h"
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 PCOR_references * new_PCOR_references(int numref)
00019 {
00020 PCOR_references *ref ;
00021 int ii,jj , nbad ;
00022
00023
00024
00025 if( numref < 1 ){
00026 fprintf( stderr , "new_PCOR_references called with numref=%d\n" , numref ) ;
00027 return NULL ;
00028 }
00029
00030
00031
00032 ref = (PCOR_references *) malloc( sizeof(PCOR_references) ) ;
00033 if( ref == NULL ){
00034 fprintf( stderr , "new_PCOR_references: malloc error for base\n" ) ;
00035 return NULL ;
00036 }
00037 ref->nref = numref ;
00038 ref->nupdate = 0 ;
00039
00040
00041
00042
00043 ref->chol = (float **) malloc( sizeof(float *) * numref ) ;
00044 if( ref->chol == NULL ){
00045 fprintf( stderr , "new_PCOR_references: malloc error for chol\n" ) ;
00046 free(ref) ; return NULL ;
00047 }
00048
00049 for( ii=0,nbad=0 ; ii < numref ; ii++ ){
00050 ref->chol[ii] = (float *) malloc( sizeof(float) * (ii+1) ) ;
00051 if( ref->chol[ii] == NULL ) nbad++ ;
00052 }
00053
00054 if( nbad > 0 ){
00055 fprintf( stderr , "new_PCOR_references: malloc error for chol[ii]\n" ) ;
00056 free_PCOR_references( ref ) ; return NULL ;
00057 }
00058
00059
00060
00061 ref->alp = (float *) malloc( sizeof(float) * numref ) ;
00062 ref->ff = (float *) malloc( sizeof(float) * numref ) ;
00063 ref->gg = (float *) malloc( sizeof(float) * numref ) ;
00064
00065 if( ref->alp == NULL || ref->ff == NULL || ref->gg == NULL ){
00066 fprintf( stderr , "new_PCOR_references: malloc error for data\n" ) ;
00067 free_PCOR_references( ref ) ; return NULL ;
00068 }
00069
00070
00071
00072 ref->rmin = (float *) malloc( sizeof(float) * numref ) ;
00073 ref->rmax = (float *) malloc( sizeof(float) * numref ) ;
00074 ref->rsum = (float *) malloc( sizeof(float) * numref ) ;
00075
00076 if( ref->rmin == NULL || ref->rmax == NULL || ref->rsum == NULL ){
00077 fprintf( stderr , "new_PCOR_references: malloc error for data\n" ) ;
00078 free_PCOR_references( ref ) ; return NULL ;
00079 }
00080
00081
00082
00083 for( ii=0 ; ii < numref ; ii++ ){
00084 for( jj=0 ; jj < ii ; jj++ ) RCH(ref,ii,jj) = 0.0 ;
00085 RCH(ref,ii,ii) = REF_EPS ;
00086 }
00087
00088 return ref ;
00089 }
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103 void update_PCOR_references(float * vec, PCOR_references * ref)
00104 {
00105 int nr = ref->nref , jj,kk ;
00106 float bold , bnew , aaa,fff,ggg ;
00107
00108 static float * zz = NULL ;
00109 static int zz_size = -1 ;
00110
00111
00112
00113 if( zz_size < nr ){
00114
00115 if( zz != NULL ) free( zz ) ;
00116 zz = (float *) malloc( sizeof(float) * nr ) ;
00117 zz_size = nr ;
00118 if( zz == NULL ){
00119 fprintf( stderr , "\nupdate_PCOR_references: can't malloc!\n" ) ;
00120 EXIT(1) ;
00121 }
00122 }
00123
00124 for( jj=0 ; jj < nr ; jj++) zz[jj] = vec[jj] ;
00125
00126
00127
00128 if( ref->nupdate == 0 ){
00129 for( jj=0 ; jj < nr ; jj++ ){
00130 ref->rmin[jj] = ref->rmax[jj] = ref->rsum[jj] = zz[jj] ;
00131 }
00132 } else {
00133 register float val ;
00134 for( jj=0 ; jj < nr ; jj++ ){
00135 val = zz[jj] ; ref->rsum[jj] += val ;
00136 if( val < ref->rmin[jj] ) ref->rmin[jj] = val ;
00137 else if( val > ref->rmax[jj] ) ref->rmax[jj] = val ;
00138 }
00139 }
00140
00141
00142
00143 bold = 1.0 ;
00144
00145 for( jj=0 ; jj < nr ; jj++ ){
00146
00147 aaa = zz[jj] / RCH(ref,jj,jj) ;
00148 bnew = sqrt( bold*bold + aaa*aaa ) ;
00149 fff = bnew / bold ;
00150 ggg = aaa / (bnew*bold) ;
00151 bold = bnew ;
00152
00153 ref->alp[jj] = aaa ;
00154 ref->ff[jj] = fff ;
00155 ref->gg[jj] = ggg ;
00156
00157 for( kk=jj ; kk < nr ; kk++ ){
00158 zz[kk] -= aaa * RCH(ref,kk,jj) ;
00159 RCH(ref,kk,jj) = fff * RCH(ref,kk,jj) + ggg * zz[kk] ;
00160 }
00161 }
00162
00163 ref->betasq = 1.0 / ( bold * bold ) ;
00164
00165 (ref->nupdate)++ ;
00166 return ;
00167 }
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178 PCOR_voxel_corr * new_PCOR_voxel_corr(int numvox, int numref)
00179 {
00180 int vox , jj ;
00181 PCOR_voxel_corr *vc ;
00182
00183
00184
00185 if( numvox < 1 ){
00186 fprintf( stderr , "new_PCOR_voxel_corr: numvox=%d\n" , numvox ) ;
00187 return NULL ;
00188 }
00189
00190
00191
00192 vc = (PCOR_voxel_corr *) malloc( sizeof(PCOR_voxel_corr) ) ;
00193 if( vc == NULL ){
00194 fprintf( stderr , "new_PCOR_voxel_corr: can't malloc base\n" ) ;
00195 return NULL ;
00196 }
00197
00198
00199
00200 vc->nvox = numvox ;
00201 vc->nref = numref ;
00202 vc->nupdate = 0 ;
00203
00204
00205
00206 vc->chrow = (float *)malloc( sizeof(float) * numvox*(numref+1) );
00207 if( vc->chrow == NULL ){
00208 fprintf( stderr , "new_PCOR_voxel_corr: can't malloc last rows\n" ) ;
00209 free( vc ) ; return NULL ;
00210 }
00211
00212
00213
00214 for( vox=0 ; vox < numvox ; vox++ ){
00215 for( jj=0 ; jj < numref ; jj++ ) VCH(vc,vox,jj) = 0 ;
00216 VCH(vc,vox,numref) = REF_EPS ;
00217 }
00218
00219 return vc ;
00220 }
00221
00222
00223
00224
00225
00226 void free_PCOR_references(PCOR_references * ref)
00227 {
00228 int ii , nr ;
00229
00230 if( ref == NULL ) return ;
00231
00232 nr = ref->nref ; if( nr <= 0 ) return ;
00233
00234 if( ref->alp != NULL ) free(ref->alp) ;
00235 if( ref->ff != NULL ) free(ref->ff) ;
00236 if( ref->gg != NULL ) free(ref->gg) ;
00237 if( ref->rmin != NULL ) free(ref->rmin);
00238 if( ref->rmax != NULL ) free(ref->rmax);
00239 if( ref->rsum != NULL ) free(ref->rsum);
00240
00241 if( ref->chol != NULL ){
00242 for( ii=0 ; ii < nr ; ii++ )
00243 if( ref->chol[ii] != NULL ) free( ref->chol[ii] ) ;
00244 free(ref->chol) ;
00245 }
00246
00247 free(ref) ;
00248 return ;
00249 }
00250
00251
00252
00253 void free_PCOR_voxel_corr(PCOR_voxel_corr * vc)
00254 {
00255 if( vc != NULL ){
00256 if( vc->chrow != NULL ) free( vc->chrow ) ;
00257 free( vc ) ;
00258 }
00259 return ;
00260 }
00261
00262
00263
00264
00265
00266
00267 #define DENEPS 1.e-5
00268
00269 void PCOR_get_pcor(PCOR_references * ref, PCOR_voxel_corr * vc, float * pcor)
00270 {
00271 int vox , nv = vc->nvox , nr = vc->nref ;
00272 float den ;
00273 static float deneps=-666.0 ;
00274
00275
00276
00277 if( vc->nref != ref->nref ){
00278 fprintf( stderr , "\nPCOR_get_pcor: reference size mismatch!\n" ) ;
00279 EXIT(1) ;
00280 }
00281
00282
00283
00284 if( deneps < 0.0 ){
00285 char * ccc = my_getenv("AFNI_PCOR_DENEPS") ;
00286 if( ccc != NULL ) deneps = strtod( ccc , NULL ) ;
00287 if( deneps < 0.0 ) deneps = DENEPS ;
00288 }
00289
00290
00291
00292 for( vox=0 ; vox < nv ; vox++ ){
00293
00294
00295 #if 0
00296 den = VCH(vc,vox,nr) ;
00297 if( den > DENEPS ){
00298 pcor[vox] = VCH(vc,vox,nr-1)
00299 / sqrt( den + SQR(VCH(vc,vox,nr-1)) ) ;
00300 } else {
00301 pcor[vox] = 0.0 ;
00302 }
00303 #else
00304
00305
00306
00307 den = VCH(vc,vox,nr) + SQR(VCH(vc,vox,nr-1)) ;
00308 if( den > deneps )
00309 pcor[vox] = VCH(vc,vox,nr-1) / sqrt(den) ;
00310 else
00311 pcor[vox] = 0.0 ;
00312 #endif
00313
00314 }
00315
00316 return ;
00317 }
00318
00319
00320
00321 void PCOR_get_mcor(PCOR_references * ref, PCOR_voxel_corr * vc, int m , float * pcor)
00322 {
00323 int vox , nv = vc->nvox , nr = vc->nref , ii ;
00324 double den ;
00325
00326
00327
00328 if( vc->nref != ref->nref ){
00329 fprintf( stderr , "\nPCOR_get_mcor: reference size mismatch!\n" ) ;
00330 EXIT(1) ;
00331 }
00332 if( m >= nr ){
00333 fprintf( stderr , "\nPCOR_get_mcor: m=%d but nref=%d\n",m,nr) ;
00334 EXIT(1) ;
00335 }
00336
00337
00338
00339 for( vox=0 ; vox < nv ; vox++ ){
00340 den = VCH(vc,vox,nr) ;
00341 switch(m){
00342 default:
00343 for( ii=1 ; ii <= m ; ii++ ) den += SQR(VCH(vc,vox,nr-ii)) ;
00344 break ;
00345
00346 case 1: den += SQR(VCH(vc,vox,nr-1)) ; break ;
00347
00348 case 2: den += SQR(VCH(vc,vox,nr-1))
00349 + SQR(VCH(vc,vox,nr-2)) ; break ;
00350
00351 case 3: den += SQR(VCH(vc,vox,nr-1))
00352 + SQR(VCH(vc,vox,nr-2))
00353 + SQR(VCH(vc,vox,nr-3)) ; break ;
00354
00355 case 4: den += SQR(VCH(vc,vox,nr-1))
00356 + SQR(VCH(vc,vox,nr-2))
00357 + SQR(VCH(vc,vox,nr-3))
00358 + SQR(VCH(vc,vox,nr-4)) ; break ;
00359 }
00360
00361 den = 1.0 - VCH(vc,vox,nr) / den ;
00362 pcor[vox] = (den > 0.0) ? sqrt(den) : 0.0 ;
00363 }
00364
00365 return ;
00366 }
00367
00368
00369
00370
00371
00372
00373 void PCOR_get_coef(PCOR_references * ref, PCOR_voxel_corr * vc, float * coef)
00374 {
00375 int vox , nv = vc->nvox , nr = vc->nref ;
00376 float scale ;
00377
00378
00379
00380 if( vc->nref != ref->nref ){
00381 fprintf( stderr , "\nPCOR_get_coef: reference size mismatch!\n" ) ;
00382 EXIT(1) ;
00383 }
00384
00385
00386
00387 scale = 1.0 / RCH(ref,nr-1,nr-1) ;
00388
00389 for( vox=0 ; vox < nv ; vox++ ){
00390 coef[vox] = scale * VCH(vc,vox,nr-1) ;
00391 }
00392
00393 return ;
00394 }
00395
00396
00397
00398
00399
00400 void PCOR_get_variance(PCOR_voxel_corr * vc, float * var)
00401 {
00402 int vox , nv = vc->nvox , nr = vc->nref , nup = vc->nupdate ;
00403 float scale ;
00404
00405
00406
00407 if( nup <= nr ){
00408 fprintf(stderr,"PCOR_get_variance: not enough data to compute!\n") ;
00409 for( vox=0 ; vox < nv ; vox++ ) var[vox] = 0.0 ;
00410 return ;
00411 }
00412
00413
00414
00415 scale = 1.0 / ( nup - nr ) ;
00416
00417 for( vox=0 ; vox < nv ; vox++ ){
00418 var[vox] = scale * VCH(vc,vox,nr) ;
00419 }
00420
00421 return ;
00422 }
00423
00424 void PCOR_get_stdev(PCOR_voxel_corr * vc, float * sig)
00425 {
00426 int vox , nv = vc->nvox ;
00427 PCOR_get_variance( vc , sig ) ;
00428 for( vox=0 ; vox < nv ; vox++ ) sig[vox] = sqrt(fabs(sig[vox])) ;
00429 return ;
00430 }
00431
00432
00433
00434
00435
00436
00437
00438
00439 void PCOR_get_lsqfit(PCOR_references * ref, PCOR_voxel_corr * vc, float *fit[] )
00440 {
00441 int vox,jj,kk , nv = vc->nvox , nr = vc->nref ;
00442 float sum ;
00443 float * ff ;
00444
00445
00446
00447 if( vc->nref != ref->nref ){
00448 fprintf( stderr , "\nPCOR_get_lsqfit: reference size mismatch!\n" ) ;
00449 EXIT(1) ;
00450 }
00451
00452 kk = 0 ;
00453 for( jj=0 ; jj < nr ; jj++ ) kk += (fit[jj] != NULL) ;
00454 if( kk == 0 ){
00455 fprintf(stderr,"PCOR_get_lsqfit: NO OUTPUT REQUESTED!\n") ;
00456 return ;
00457 }
00458
00459 ff = (float *) malloc( sizeof(float) * nr ) ;
00460 if( ff == NULL ){
00461 fprintf( stderr, "\nPCOR_get_lsqfit: can't malloc workspace!\n") ;
00462 EXIT(1) ;
00463 }
00464
00465
00466
00467 for( vox=0 ; vox < nv ; vox++ ){
00468
00469 for( jj=nr-1 ; jj >=0 ; jj-- ){
00470 sum = VCH(vc,vox,jj) ;
00471 for( kk=jj+1 ; kk < nr ; kk++ ) sum -= ff[kk] * RCH(ref,kk,jj) ;
00472 ff[jj] = sum / RCH(ref,jj,jj) ;
00473 }
00474
00475 for( jj=0 ; jj < nr ; jj++ )
00476 if( fit[jj] != NULL ) fit[jj][vox] = ff[jj] ;
00477 }
00478
00479 free( ff ) ;
00480 return ;
00481 }
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494 void PCOR_get_perc(PCOR_references * ref, PCOR_voxel_corr * vc,
00495 float * coef, float * bline, int basaver )
00496 {
00497 int vox,jj,kk , nv=vc->nvox , nr=vc->nref , nup=ref->nupdate ;
00498 float sum , base , rdif , rmin , rmax ;
00499 float * ff , * bb , * dd ;
00500
00501
00502
00503 if( vc->nref != ref->nref ){
00504 fprintf( stderr , "\nPCOR_get_perc: reference size mismatch!\n" ) ;
00505 EXIT(1) ;
00506 }
00507
00508 if( coef == NULL && bline == NULL ) return ;
00509
00510
00511
00512 ff = (float *) malloc( sizeof(float) * nr ) ;
00513 bb = (float *) malloc( sizeof(float) * nr ) ;
00514 dd = (float *) malloc( sizeof(float) * nr ) ;
00515 if( ff == NULL || bb == NULL || dd == NULL ){
00516 fprintf( stderr, "\nPCOR_get_perc: can't malloc workspace!\n") ;
00517 EXIT(1) ;
00518 }
00519
00520
00521
00522 rmin = ref->rmin[nr-1] ;
00523 rmax = ref->rmax[nr-1] ; rdif = 100.0 * (rmax-rmin) ;
00524 if( rdif == 0.0 ){
00525
00526 if( coef != NULL ) for( vox=0; vox < nv; vox++ ) coef[vox] = 0.0;
00527 if( bline != NULL ) for( vox=0; vox < nv; vox++ ) bline[vox] = 0.0;
00528 fprintf(stderr,"\nPCOR_get_perc: ref vector has no range!\n") ;
00529 return ;
00530 }
00531
00532 for( jj=0 ; jj < nr ; jj++ ){
00533 bb[jj] = ref->rsum[jj] / nup ;
00534 dd[jj] = 1.0 / RCH(ref,jj,jj) ;
00535 }
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546 for( vox=0 ; vox < nv ; vox++ ){
00547 for( jj=nr-1 ; jj >=0 ; jj-- ){
00548 sum = VCH(vc,vox,jj) ;
00549 for( kk=jj+1 ; kk < nr ; kk++ ) sum -= ff[kk] * RCH(ref,kk,jj) ;
00550 ff[jj] = sum * dd[jj] ;
00551 }
00552
00553 switch( basaver ){
00554 default:
00555 case 0: base = ff[nr-1] * rmin ; break;
00556 case 1: base = ff[nr-1] * bb[nr-1]; break;
00557 case 2: base = ff[nr-1] * rmax ; break;
00558 }
00559
00560 for( jj=0 ; jj < nr-1 ; jj++ )
00561 base += ff[jj] * bb[jj] ;
00562
00563 if( coef != NULL ) coef[vox] = (base > 0.0) ? ff[nr-1] * rdif / base
00564 : 0.0 ;
00565 if( bline != NULL ) bline[vox] = base ;
00566 }
00567
00568 free(ff) ; free(bb) ; free(dd) ;
00569 return ;
00570 }
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581 void PCOR_get_pcor_and_coef(PCOR_references * ref , PCOR_voxel_corr * vc ,
00582 float pcthresh , float * pcor , float * coef )
00583 {
00584 int vox , nv = vc->nvox , nr = vc->nref ;
00585 float den , num , scale ;
00586 float pc , co , thfac ;
00587
00588
00589
00590 if( vc->nref != ref->nref ){
00591 fprintf( stderr , "\nPCOR_get_pcor_and_coef: reference size mismatch!\n" ) ;
00592 EXIT(1) ;
00593 }
00594
00595 scale = 1.0 / RCH(ref,nr-1,nr-1) ;
00596 thfac = SQR(pcthresh)/(1.0-SQR(pcthresh)) ;
00597
00598
00599
00600 if( pcthresh <= 0.0 ){
00601 for( vox=0 ; vox < nv ; vox++ ){
00602 den = VCH(vc,vox,nr) ;
00603 num = VCH(vc,vox,nr-1) ;
00604 pcor[vox] = num / sqrt(den+SQR(num)) ;
00605 coef[vox] = scale * num ;
00606 }
00607 } else {
00608 thfac = SQR(pcthresh)/(1.0-SQR(pcthresh)) ;
00609 for( vox=0 ; vox < nv ; vox++ ){
00610 den = VCH(vc,vox,nr) ;
00611 num = VCH(vc,vox,nr-1) ;
00612 if( SQR(num) > thfac*den ){
00613 pcor[vox] = num / sqrt(den+SQR(num)) ;
00614 coef[vox] = scale * num ;
00615 } else {
00616 pcor[vox] = coef[vox] = 0.0 ;
00617 }
00618 }
00619 }
00620
00621 return ;
00622 }