00001
00002
00003
00004
00005
00006
00007 #include <stdio.h>
00008 #include <math.h>
00009 #include <stdlib.h>
00010
00011 #include "pcor.h"
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 references * new_references(numref)
00030 int numref;
00031 {
00032 references *ref ;
00033 int ii,jj ;
00034
00035
00036
00037 if( numref < 1 ){
00038 fprintf( stderr , "new_references called with numref=%d\n" , numref ) ;
00039 exit(-1) ;
00040 }
00041
00042
00043
00044 ref = (references *) malloc( sizeof(references) ) ;
00045 if( ref == NULL ){
00046 fprintf( stderr , "new_references: malloc error for base\n" ) ;
00047 exit(-1) ;
00048 }
00049 ref->nref = numref ;
00050 ref->nupdate = 0 ;
00051
00052
00053
00054
00055 ref->chol = (ref_float **) malloc( sizeof(ref_float *) * numref ) ;
00056 if( ref->chol == NULL ){
00057 fprintf( stderr , "new_references: malloc error for chol\n" ) ;
00058 exit(-1) ;
00059 }
00060
00061 for( ii=0 ; ii < numref ; ii++ ){
00062 ref->chol[ii] = (ref_float *) malloc( sizeof(ref_float) * (ii+1) ) ;
00063 if( ref->chol[ii] == NULL ){
00064 fprintf( stderr , "new_references: malloc error for chol[ii]\n" ) ;
00065 exit(-1) ;
00066 }
00067 }
00068
00069
00070
00071 ref->alp = (ref_float *) malloc( sizeof(ref_float) * numref ) ;
00072 ref->ff = (ref_float *) malloc( sizeof(ref_float) * numref ) ;
00073 ref->gg = (ref_float *) malloc( sizeof(ref_float) * numref ) ;
00074
00075 if( ref->alp == NULL || ref->ff == NULL || ref->gg == NULL ){
00076 fprintf( stderr , "new_references: malloc error for data\n" ) ;
00077 exit(-1) ;
00078 }
00079
00080
00081
00082 for( ii=0 ; ii < numref ; ii++ ){
00083 for( jj=0 ; jj < ii ; jj++ ) RCH(ref,ii,jj) = 0.0 ;
00084 RCH(ref,ii,ii) = REF_EPS ;
00085 #ifdef OV_DEBUG2
00086 ref->alp[ii] = ref->ff[ii] = ref->gg[ii] = 0.0 ;
00087 #endif
00088 }
00089
00090 #ifdef OV_DEBUG2
00091 ref->betasq = 0.0 ;
00092 #endif
00093
00094 return ref ;
00095 }
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109 void update_references(vec, ref)
00110 float *vec;
00111 references *ref;
00112 {
00113 int nr = ref->nref , jj,kk ;
00114 ref_float bold , bnew , aaa,fff,ggg ;
00115 static ref_float * zz = NULL ;
00116 static int zz_size = -1 ;
00117
00118 #ifdef OV_DEBUG2
00119 static ref_float qinput[50] ;
00120 #endif
00121
00122
00123
00124 if( zz_size < nr ){
00125
00126 if( zz != NULL ) free( zz ) ;
00127 zz = (ref_float *) malloc( sizeof(ref_float) * nr ) ;
00128 zz_size = nr ;
00129 if( zz == NULL ){
00130 fprintf( stderr , "update_references: cannot malloc!\n" ) ;
00131 exit(-1) ;
00132 }
00133 }
00134
00135 for( jj=0 ; jj < nr ; jj++) zz[jj] = (ref_float) vec[jj] ;
00136
00137 #ifdef OV_DEBUG2
00138 for( jj=0 ; jj < nr ; jj++) qinput[jj] += SQR(zz[jj]) ;
00139
00140 REF_DUMP(ref,"before update") ;
00141 fprintf(stderr," input vec= ") ;
00142 for( jj=0 ; jj < nr ; jj++ ) fprintf(stderr,"%11.4e ",zz[jj]) ;
00143 fprintf(stderr,"\n") ;
00144 #endif
00145
00146
00147
00148 bold = 1.0 ;
00149
00150 for( jj=0 ; jj < nr ; jj++ ){
00151
00152 aaa = zz[jj] / RCH(ref,jj,jj) ;
00153 bnew = sqrt( bold*bold + aaa*aaa ) ;
00154 fff = bnew / bold ;
00155 ggg = aaa / (bnew*bold) ;
00156 bold = bnew ;
00157
00158 ref->alp[jj] = aaa ;
00159 ref->ff[jj] = fff ;
00160 ref->gg[jj] = ggg ;
00161
00162 for( kk=jj ; kk < nr ; kk++ ){
00163 zz[kk] -= aaa * RCH(ref,kk,jj) ;
00164 RCH(ref,kk,jj) = fff * RCH(ref,kk,jj) + ggg * zz[kk] ;
00165 }
00166 }
00167
00168 ref->betasq = 1.0 / ( bold * bold ) ;
00169
00170 #ifdef OV_DEBUG2
00171 REF_DUMP(ref,"after update") ;
00172 fprintf(stderr," qsum of input vecs= ") ;
00173 for( jj=0 ; jj < nr ; jj++ ) fprintf(stderr,"%11.4e ",qinput[jj]) ;
00174 fprintf(stderr,"\n") ;
00175 #endif
00176
00177 (ref->nupdate)++ ;
00178 return ;
00179 }
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190 voxel_corr * new_voxel_corr(numvox, numref)
00191 int numvox;
00192 int numref;
00193 {
00194 int vox , jj ;
00195 voxel_corr *vc ;
00196
00197
00198
00199 if( numvox < 1 ){
00200 fprintf( stderr , "new_voxel_corr: numvox=%d\n" , numvox ) ;
00201 exit(-1) ;
00202 }
00203
00204
00205
00206 vc = (voxel_corr *) malloc( sizeof(voxel_corr) ) ;
00207 if( vc == NULL ){
00208 fprintf( stderr , "new_voxel_corr: cannot malloc base\n" ) ;
00209 exit(-1) ;
00210 }
00211
00212
00213
00214 vc->nvox = numvox ;
00215 vc->nref = numref ;
00216 vc->nupdate = 0 ;
00217
00218
00219
00220 vc->chrow = (ref_float *)malloc( sizeof(ref_float) * numvox*(numref+1) );
00221 if( vc->chrow == NULL ){
00222 fprintf( stderr , "new_voxel_corr: cannot malloc last rows\n" ) ;
00223 exit(-1) ;
00224 }
00225
00226
00227
00228 for( vox=0 ; vox < numvox ; vox++ ){
00229 for( jj=0 ; jj < numref ; jj++ ) VCH(vc,vox,jj) = 0 ;
00230 VCH(vc,vox,numref) = REF_EPS ;
00231 }
00232
00233 return vc ;
00234 }
00235
00236
00237
00238
00239
00240 void free_references(ref)
00241 references *ref;
00242 {
00243 int ii , nr ;
00244
00245 if( ref == NULL ) return ;
00246
00247 nr = ref->nref ; if( nr <= 0 ) return ;
00248
00249 free(ref->alp) ; free(ref->ff) ; free(ref->gg) ;
00250
00251 for( ii=0 ; ii < nr ; ii++ ) free( ref->chol[ii] ) ;
00252
00253 free(ref->chol) ; free(ref) ;
00254
00255 return ;
00256 }
00257
00258
00259
00260 void free_voxel_corr(vc)
00261 voxel_corr *vc;
00262 {
00263 if( vc != NULL ){
00264 free( vc->chrow ) ;
00265 free( vc ) ;
00266 }
00267 return ;
00268 }
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279 void update_voxel_corr(vdata, ref, vc)
00280 vox_data *vdata;
00281 references *ref;
00282 voxel_corr *vc;
00283 {
00284 int vox , jj ,
00285 nv = vc->nvox ,
00286 nr = vc->nref ;
00287
00288 ref_float *aaa = ref->alp ,
00289 *fff = ref->ff ,
00290 *ggg = ref->gg ;
00291
00292 ref_float zz ,
00293 bq = ref->betasq ;
00294
00295 #ifdef OV_DEBUG2
00296 static ref_float qvox = 0.0 ;
00297 #endif
00298
00299
00300
00301 if( vc->nref != ref->nref ){
00302 fprintf( stderr , "update_voxel_corr: reference size mismatch!\n" ) ;
00303 exit(-1) ;
00304 }
00305
00306 #ifdef OV_DEBUG2
00307 VOX_DUMP(vc,VD,"before update") ;
00308 #ifdef VOX_SHORT
00309 fprintf(stderr," integer input data = %d\n" , vdata[VD] ) ;
00310 #else
00311 fprintf(stderr," float input data = %11.4e\n" , vdata[VD] ) ;
00312 #endif
00313 qvox += SQR(vdata[VD]) ;
00314 #endif
00315
00316
00317
00318 #ifdef EXPAND_UPDATE
00319 #define UPZZ(j) zz -= aaa[j] * VCH(vc,vox,j)
00320 #define UPCH(j) VCH(vc,vox,j) = fff[j] * VCH(vc,vox,j) + ggg[j] * zz
00321 #define UPLL(j) VCH(vc,vox,j) += bq * zz * zz
00322
00323 switch( nr ){
00324 default:
00325 #endif
00326
00327
00328
00329 for( vox=0 ; vox < nv ; vox++ ){
00330
00331
00332
00333 zz = (ref_float) vdata[vox] ;
00334 for( jj=0 ; jj < nr ; jj++ ){
00335 zz -= aaa[jj] * VCH(vc,vox,jj) ;
00336 VCH(vc,vox,jj) = fff[jj] * VCH(vc,vox,jj) + ggg[jj] * zz ;
00337 }
00338 VCH(vc,vox,nr) += bq * zz * zz ;
00339 }
00340
00341 #ifdef EXPAND_UPDATE
00342 break ;
00343
00344 case 1:
00345 for( vox=0 ; vox < nv ; vox++ ){
00346 zz = (ref_float) vdata[vox] ;
00347 UPZZ(0) ; UPCH(0) ; UPLL(1) ;
00348 }
00349 break ;
00350
00351 case 2:
00352 for( vox=0 ; vox < nv ; vox++ ){
00353 zz = (ref_float) vdata[vox] ;
00354 UPZZ(0) ; UPCH(0) ;
00355 UPZZ(1) ; UPCH(1) ; UPLL(2) ;
00356 }
00357 break ;
00358
00359 case 3:
00360 for( vox=0 ; vox < nv ; vox++ ){
00361 zz = (ref_float) vdata[vox] ;
00362 UPZZ(0) ; UPCH(0) ;
00363 UPZZ(1) ; UPCH(1) ;
00364 UPZZ(2) ; UPCH(2) ; UPLL(3) ;
00365 }
00366 break ;
00367
00368 case 4:
00369 for( vox=0 ; vox < nv ; vox++ ){
00370 zz = (ref_float) vdata[vox] ;
00371 UPZZ(0) ; UPCH(0) ;
00372 UPZZ(1) ; UPCH(1) ;
00373 UPZZ(2) ; UPCH(2) ;
00374 UPZZ(3) ; UPCH(3) ; UPLL(4) ;
00375 }
00376 break ;
00377
00378 case 5:
00379 for( vox=0 ; vox < nv ; vox++ ){
00380 zz = (ref_float) vdata[vox] ;
00381 UPZZ(0) ; UPCH(0) ;
00382 UPZZ(1) ; UPCH(1) ;
00383 UPZZ(2) ; UPCH(2) ;
00384 UPZZ(3) ; UPCH(3) ;
00385 UPZZ(4) ; UPCH(4) ; UPLL(5) ;
00386 }
00387 break ;
00388
00389 case 6:
00390 for( vox=0 ; vox < nv ; vox++ ){
00391 zz = (ref_float) vdata[vox] ;
00392 UPZZ(0) ; UPCH(0) ;
00393 UPZZ(1) ; UPCH(1) ;
00394 UPZZ(2) ; UPCH(2) ;
00395 UPZZ(3) ; UPCH(3) ;
00396 UPZZ(4) ; UPCH(4) ;
00397 UPZZ(5) ; UPCH(5) ; UPLL(6) ;
00398 }
00399 break ;
00400
00401 case 7:
00402 for( vox=0 ; vox < nv ; vox++ ){
00403 zz = (ref_float) vdata[vox] ;
00404 UPZZ(0) ; UPCH(0) ;
00405 UPZZ(1) ; UPCH(1) ;
00406 UPZZ(2) ; UPCH(2) ;
00407 UPZZ(3) ; UPCH(3) ;
00408 UPZZ(4) ; UPCH(4) ;
00409 UPZZ(5) ; UPCH(5) ;
00410 UPZZ(6) ; UPCH(6) ; UPLL(7) ;
00411 }
00412 break ;
00413
00414 case 8:
00415 for( vox=0 ; vox < nv ; vox++ ){
00416 zz = (ref_float) vdata[vox] ;
00417 UPZZ(0) ; UPCH(0) ;
00418 UPZZ(1) ; UPCH(1) ;
00419 UPZZ(2) ; UPCH(2) ;
00420 UPZZ(3) ; UPCH(3) ;
00421 UPZZ(4) ; UPCH(4) ;
00422 UPZZ(5) ; UPCH(5) ;
00423 UPZZ(6) ; UPCH(6) ;
00424 UPZZ(7) ; UPCH(7) ; UPLL(8) ;
00425 }
00426 break ;
00427
00428 case 9:
00429 for( vox=0 ; vox < nv ; vox++ ){
00430 zz = (ref_float) vdata[vox] ;
00431 UPZZ(0) ; UPCH(0) ;
00432 UPZZ(1) ; UPCH(1) ;
00433 UPZZ(2) ; UPCH(2) ;
00434 UPZZ(3) ; UPCH(3) ;
00435 UPZZ(4) ; UPCH(4) ;
00436 UPZZ(5) ; UPCH(5) ;
00437 UPZZ(6) ; UPCH(6) ;
00438 UPZZ(7) ; UPCH(7) ;
00439 UPZZ(8) ; UPCH(8) ; UPLL(9) ;
00440 }
00441 break ;
00442
00443 }
00444 #endif
00445
00446 #ifdef OV_DEBUG2
00447 VOX_DUMP(vc,VD,"after update") ;
00448 fprintf(stderr," qsum of vox[VD]=%11.4e\n",qvox) ;
00449 #endif
00450
00451 (vc->nupdate)++ ;
00452 return ;
00453 }
00454
00455
00456
00457
00458
00459 void get_pcor(ref, vc, pcor)
00460 references *ref;
00461 voxel_corr *vc;
00462 float *pcor;
00463 {
00464 int vox , nv = vc->nvox , nr = vc->nref ;
00465 ref_float den ;
00466 #define DENEPS 1.e-5
00467
00468
00469
00470 if( vc->nref != ref->nref ){
00471 fprintf( stderr , "get_pcor: reference size mismatch!\n" ) ;
00472 exit(-1) ;
00473 }
00474
00475
00476
00477 for( vox=0 ; vox < nv ; vox++ ){
00478
00479 den = VCH(vc,vox,nr) ;
00480 if( den > DENEPS ){
00481 pcor[vox] = VCH(vc,vox,nr-1)
00482 / sqrt( den + SQR(VCH(vc,vox,nr-1)) ) ;
00483 } else {
00484 pcor[vox] = 0.0 ;
00485 }
00486
00487 }
00488
00489 #ifdef OV_DEBUG2
00490 fprintf(stderr,"get_pcor: pcor[VD]=%11.4e\n",pcor[VD]) ;
00491 #endif
00492
00493 return ;
00494 }
00495
00496
00497
00498
00499
00500 void get_coef(ref, vc, coef)
00501 references *ref;
00502 voxel_corr *vc;
00503 float *coef;
00504 {
00505 int vox , nv = vc->nvox , nr = vc->nref ;
00506 ref_float scale ;
00507
00508
00509
00510 if( vc->nref != ref->nref ){
00511 fprintf( stderr , "get_coef: reference size mismatch!\n" ) ;
00512 exit(-1) ;
00513 }
00514
00515
00516
00517 scale = 1.0 / RCH(ref,nr-1,nr-1) ;
00518
00519 for( vox=0 ; vox < nv ; vox++ ){
00520 coef[vox] = scale * VCH(vc,vox,nr-1) ;
00521 }
00522
00523 #ifdef OV_DEBUG2
00524 fprintf(stderr,"get_coef: coef[VD]=%11.4e\n",coef[VD]) ;
00525 #endif
00526
00527 return ;
00528 }
00529
00530
00531
00532
00533
00534 void get_variance(vc, var)
00535 voxel_corr *vc;
00536 float *var;
00537 {
00538 int vox , nv = vc->nvox , nr = vc->nref , nup = vc->nupdate ;
00539 ref_float scale ;
00540
00541
00542
00543 if( nup <= nr ){
00544 fprintf(stderr,"get_variance: not enough data to compute!\n") ;
00545 for( vox=0 ; vox < nv ; vox++ ) var[vox] = 0.0 ;
00546 return ;
00547 }
00548
00549
00550
00551 scale = 1.0 / ( nup - nr ) ;
00552
00553 for( vox=0 ; vox < nv ; vox++ ){
00554 var[vox] = scale * VCH(vc,vox,nr) ;
00555 }
00556
00557 #ifdef OV_DEBUG2
00558 fprintf(stderr,"get_variance: var[VD]=%11.4e\n",var[VD]) ;
00559 #endif
00560
00561 return ;
00562 }
00563
00564
00565
00566
00567
00568
00569
00570
00571 void get_lsqfit(ref, vc, fit)
00572 references *ref;
00573 voxel_corr *vc;
00574 float *fit[] ;
00575 {
00576 int vox,jj,kk , nv = vc->nvox , nr = vc->nref ;
00577 ref_float sum ;
00578 ref_float * ff ;
00579
00580
00581
00582 if( vc->nref != ref->nref ){
00583 fprintf( stderr , "get_lsqfit: reference size mismatch!\n" ) ;
00584 exit(-1) ;
00585 }
00586
00587 kk = 0 ;
00588 for( jj=0 ; jj < nr ; jj++ ) kk += (fit[jj] != NULL) ;
00589 if( kk == 0 ){
00590 fprintf(stderr,"get_lsqfit: NO OUTPUT REQUESTED!\n") ;
00591 return ;
00592 }
00593
00594 ff = (ref_float *) malloc( sizeof(ref_float) * nr ) ;
00595 if( ff == NULL ){
00596 fprintf( stderr, "get_lsqfit: cannot malloc workspace!\n") ;
00597 exit(-1) ;
00598 }
00599
00600
00601
00602 for( vox=0 ; vox < nv ; vox++ ){
00603
00604 for( jj=nr-1 ; jj >=0 ; jj-- ){
00605 sum = VCH(vc,vox,jj) ;
00606 for( kk=jj+1 ; kk < nr ; kk++ ) sum -= ff[kk] * RCH(ref,kk,jj) ;
00607 ff[jj] = sum / RCH(ref,jj,jj) ;
00608 }
00609
00610 for( jj=0 ; jj < nr ; jj++ )
00611 if( fit[jj] != NULL ) fit[jj][vox] = ff[jj] ;
00612 }
00613
00614 free( ff ) ;
00615 return ;
00616 }
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629 void get_pcor_thresh_coef(ref, vc, pcthresh, cothresh, pcor, coef, thr)
00630 references *ref;
00631 voxel_corr *vc;
00632 float pcthresh;
00633 float cothresh;
00634 float *pcor;
00635 float *coef;
00636 thresh_result *thr;
00637 {
00638 int vox , nv = vc->nvox , nr = vc->nref ;
00639 ref_float den , num , scale ;
00640 #define DENEPS 1.e-5
00641
00642 float pc , co , thfac ;
00643 int npc_pos=0 , npc_neg=0 , nco_pos=0 , nco_neg=0 ;
00644 float mpc_pos=0., mpc_neg=0., mco_pos=0., mco_neg=0.;
00645
00646 int do_pcth , do_coth ;
00647
00648
00649
00650 if( vc->nref != ref->nref ){
00651 fprintf( stderr , "get_pcor: reference size mismatch!\n" ) ;
00652 exit(-1) ;
00653 }
00654
00655 scale = 1.0 / RCH(ref,nr-1,nr-1) ;
00656 thfac = SQR(pcthresh)/(1.0-SQR(pcthresh)) ;
00657 do_pcth = pcthresh <= 0.0 ;
00658 do_coth = cothresh > 0.0 ;
00659
00660
00661
00662 for( vox=0 ; vox < nv ; vox++ ){
00663
00664 den = VCH(vc,vox,nr) ;
00665 num = VCH(vc,vox,nr-1) ;
00666 if( do_pcth || SQR(num) > thfac*den ){
00667
00668 pc = pcor[vox] = num / sqrt(den+SQR(num)) ;
00669 co = coef[vox] = scale * num ;
00670
00671 if( pc > 0 ){
00672 npc_pos++ ;
00673 if( mpc_pos < pc ) mpc_pos = pc ;
00674 if( mco_pos < co ) mco_pos = co ;
00675 } else {
00676 npc_neg++ ;
00677 if( mpc_neg > pc ) mpc_neg = pc ;
00678 if( mco_neg > co ) mco_neg = co ;
00679 }
00680
00681 } else {
00682 pcor[vox] = coef[vox] = 0.0 ;
00683 }
00684
00685 }
00686
00687 nco_pos = npc_pos ;
00688 nco_neg = npc_neg ;
00689
00690
00691
00692 if( do_coth && nco_pos+nco_neg > 0 ){
00693
00694 thfac = cothresh * MAX(mco_pos,-mco_neg) ;
00695
00696 for( vox=0 ; vox < nv ; vox++ ){
00697 if( coef[vox] > 0.0 && coef[vox] < thfac ){
00698 coef[vox] = 0.0 ;
00699 nco_pos-- ;
00700 } else if( coef[vox] < 0.0 && coef[vox] > -thfac ){
00701 coef[vox] = 0.0 ;
00702 nco_neg-- ;
00703 }
00704 }
00705 }
00706
00707
00708
00709 thr->num_pcor_pos = npc_pos ;
00710 thr->num_pcor_neg = npc_neg ;
00711 thr->max_pcor_pos = mpc_pos ;
00712 thr->max_pcor_neg = mpc_neg ;
00713
00714 thr->num_coef_pos = nco_pos ;
00715 thr->num_coef_neg = nco_neg ;
00716 thr->max_coef_pos = mco_pos ;
00717 thr->max_coef_neg = mco_neg ;
00718
00719 return ;
00720 }