Doxygen Source Code Documentation
afni_pcor.h File Reference
#include <stdlib.h>#include <stdio.h>#include <math.h>#include "mrilib.h"Go to the source code of this file.
Define Documentation
|
|
combine two interpreted tokens into one using TWO_TWO * Definition at line 51 of file afni_pcor.h. |
|
|
Definition at line 82 of file afni_pcor.h. Referenced by new_PCOR_references(), PCOR_get_coef(), PCOR_get_lsqfit(), PCOR_get_pcor_and_coef(), PCOR_get_perc(), and update_PCOR_references(). |
|
|
Definition at line 53 of file afni_pcor.h. |
|
|
Definition at line 105 of file afni_pcor.h. Referenced by new_PCOR_voxel_corr(), PCOR_get_coef(), PCOR_get_lsqfit(), PCOR_get_mcor(), PCOR_get_pcor(), PCOR_get_pcor_and_coef(), PCOR_get_perc(), PCOR_get_variance(), and PCOR_UPDATE(). |
Function Documentation
|
|
Definition at line 226 of file afni_pcor.c. References PCOR_references::alp, PCOR_references::chol, PCOR_references::ff, free, PCOR_references::gg, nr, PCOR_references::nref, ref, PCOR_references::rmax, PCOR_references::rmin, and PCOR_references::rsum. Referenced by AFNI_fimmer_compute(), CORREL_main(), fim3d_fimmer_compute(), new_PCOR_references(), and RT_fim_recurse().
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); /* 14 Jan 1998 */
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 }
|
|
|
Definition at line 253 of file afni_pcor.c. References PCOR_voxel_corr::chrow, free, and vc. Referenced by AFNI_fimmer_compute(), CORREL_main(), fim3d_fimmer_compute(), and RT_fim_recurse().
|
|
|
Definition at line 18 of file afni_pcor.c. References PCOR_references::alp, PCOR_references::chol, PCOR_references::ff, free, free_PCOR_references(), PCOR_references::gg, malloc, PCOR_references::nref, numref, PCOR_references::nupdate, RCH, ref, PCOR_references::rmax, PCOR_references::rmin, and PCOR_references::rsum. Referenced by AFNI_fimmer_compute(), CORREL_main(), fim3d_fimmer_compute(), and RT_fim_recurse().
00019 {
00020 PCOR_references *ref ;
00021 int ii,jj , nbad ;
00022
00023 /*** check input for reasonableness ***/
00024
00025 if( numref < 1 ){
00026 fprintf( stderr , "new_PCOR_references called with numref=%d\n" , numref ) ;
00027 return NULL ;
00028 }
00029
00030 /*** allocate storage for top level data ***/
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 ; /* June 1995: not updated at all yet */
00039
00040 /*** allocate storage for Cholesky factor
00041 (an array of rows, row #ii is length ii+1, for ii=0..numref-1) ***/
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 /*** allocate storage for vectors of alpha, f, g ***/
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 /*** 14 Jan 1998: space for ref vector statistics ***/
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 /*** initialize Cholesky factor ***/
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 }
|
|
||||||||||||
|
Definition at line 178 of file afni_pcor.c. References PCOR_voxel_corr::chrow, free, malloc, PCOR_voxel_corr::nref, numref, PCOR_voxel_corr::nupdate, PCOR_voxel_corr::nvox, vc, and VCH. Referenced by AFNI_fimmer_compute(), CORREL_main(), fim3d_fimmer_compute(), and RT_fim_recurse().
00179 {
00180 int vox , jj ;
00181 PCOR_voxel_corr *vc ;
00182
00183 /*** check input for OK-osity ***/
00184
00185 if( numvox < 1 ){
00186 fprintf( stderr , "new_PCOR_voxel_corr: numvox=%d\n" , numvox ) ;
00187 return NULL ;
00188 }
00189
00190 /*** get the base storage ***/
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 /*** setup the references common to all voxels ***/
00199
00200 vc->nvox = numvox ;
00201 vc->nref = numref ;
00202 vc->nupdate = 0 ; /* June 1995: not updated at all yet */
00203
00204 /*** setup the storage of the last row for each voxel ***/
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 /*** initialize each voxel ***/
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 }
|
|
||||||||||||||||
|
Definition at line 373 of file afni_pcor.c. References EXIT, nr, PCOR_references::nref, PCOR_voxel_corr::nref, PCOR_voxel_corr::nvox, RCH, ref, scale, vc, and VCH. Referenced by AFNI_fimmer_compute(), fim3d_fimmer_compute(), and RT_fim_recurse().
00374 {
00375 int vox , nv = vc->nvox , nr = vc->nref ;
00376 float scale ;
00377
00378 /*** check inputs for OK-ness ***/
00379
00380 if( vc->nref != ref->nref ){
00381 fprintf( stderr , "\nPCOR_get_coef: reference size mismatch!\n" ) ;
00382 EXIT(1) ;
00383 }
00384
00385 /*** Work ***/
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 }
|
|
||||||||||||||||
|
Definition at line 439 of file afni_pcor.c. References EXIT, fit, free, malloc, nr, PCOR_references::nref, PCOR_voxel_corr::nref, PCOR_voxel_corr::nvox, RCH, ref, vc, and VCH.
00440 {
00441 int vox,jj,kk , nv = vc->nvox , nr = vc->nref ;
00442 float sum ;
00443 float * ff ;
00444
00445 /*** check inputs for OK-ness ***/
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 /*** for each voxel, compute the nr fit coefficients (backwards) ***/
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 }
|
|
||||||||||||||||||||
|
Definition at line 321 of file afni_pcor.c. References EXIT, nr, PCOR_references::nref, PCOR_voxel_corr::nref, PCOR_voxel_corr::nvox, pcor, ref, SQR, vc, and VCH.
00322 {
00323 int vox , nv = vc->nvox , nr = vc->nref , ii ;
00324 double den ;
00325
00326 /*** check inputs for OK-ness ***/
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 /*** Work ***/
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 }
|
|
||||||||||||||||
|
Definition at line 269 of file afni_pcor.c. References DENEPS, EXIT, my_getenv(), nr, PCOR_references::nref, PCOR_voxel_corr::nref, PCOR_voxel_corr::nvox, pcor, ref, SQR, strtod(), vc, and VCH. Referenced by AFNI_fimmer_compute(), CORREL_main(), fim3d_fimmer_compute(), and RT_fim_recurse().
00270 {
00271 int vox , nv = vc->nvox , nr = vc->nref ;
00272 float den ;
00273 static float deneps=-666.0 ; /* 28 Sep 1999 */
00274
00275 /*** check inputs for OK-ness ***/
00276
00277 if( vc->nref != ref->nref ){
00278 fprintf( stderr , "\nPCOR_get_pcor: reference size mismatch!\n" ) ;
00279 EXIT(1) ;
00280 }
00281
00282 /* 28 Sep 1999: load user option for denominator epsilon */
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 /*** Work ***/
00291
00292 for( vox=0 ; vox < nv ; vox++ ){
00293
00294 /* change below made 15 July 1998 */
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 /*----- Allow for numerical roundoff error, 26 August 1998 BDW -----*/
00305 /*----- Replace DENEPS with deneps: 28 September 1999 RWC ---*/
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 }
|
|
||||||||||||||||||||||||
|
06 Dec 2000: zero out both coef and bline in this case * Definition at line 581 of file afni_pcor.c. References EXIT, nr, PCOR_references::nref, PCOR_voxel_corr::nref, PCOR_voxel_corr::nvox, pcor, pcthresh, RCH, ref, scale, SQR, vc, and VCH.
00583 {
00584 int vox , nv = vc->nvox , nr = vc->nref ;
00585 float den , num , scale ;
00586 float pc , co , thfac ;
00587
00588 /*** check inputs for OK-ness ***/
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) ; /* for coef calculation */
00596 thfac = SQR(pcthresh)/(1.0-SQR(pcthresh)) ;
00597
00598 /*** Compute pcor and coef, thresholded on pcthresh ***/
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 ){ /* fancy threshold test */
00613 pcor[vox] = num / sqrt(den+SQR(num)) ;
00614 coef[vox] = scale * num ;
00615 } else { /* fails pcor thresh */
00616 pcor[vox] = coef[vox] = 0.0 ;
00617 }
00618 }
00619 }
00620
00621 return ;
00622 }
|
|
||||||||||||||||||||||||
|
Definition at line 494 of file afni_pcor.c. References base, EXIT, free, malloc, nr, PCOR_references::nref, PCOR_voxel_corr::nref, PCOR_references::nupdate, PCOR_voxel_corr::nvox, RCH, ref, PCOR_references::rmax, PCOR_references::rmin, PCOR_references::rsum, vc, and VCH. Referenced by AFNI_fimmer_compute().
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 /*** check inputs for OK-ness ***/
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 ; /* nothing to do */
00509
00510 /*** Setup ***/
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 /* range of last reference */
00521
00522 rmin = ref->rmin[nr-1] ;
00523 rmax = ref->rmax[nr-1] ; rdif = 100.0 * (rmax-rmin) ;
00524 if( rdif == 0.0 ){
00525 /** 06 Dec 2000: zero out both coef and bline in this case **/
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 ; /* average of each ref */
00534 dd[jj] = 1.0 / RCH(ref,jj,jj) ; /* factor for loop below */
00535 }
00536
00537 /*** Work: jj=nr-2
00538 x(t) = fit[nr-1] * ref(t,nr-1) + SUM fit[jj] * ref(t,jj)
00539 jj=0
00540
00541 The percent change due to ref(t,nr-1) is computed by scaling
00542 this ref to run from 0 to 1 (that's why rmin and rmax are used),
00543 and by computing the baseline as fit[nr-1]*rmin + the sum
00544 of the averages of the other refs times their fit coefficients. ***/
00545
00546 for( vox=0 ; vox < nv ; vox++ ){
00547 for( jj=nr-1 ; jj >=0 ; jj-- ){ /* compute fit coefficients */
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; /* baseline = bottom */
00556 case 1: base = ff[nr-1] * bb[nr-1]; break; /* baseline = average */
00557 case 2: base = ff[nr-1] * rmax ; break; /* baseline = top */
00558 }
00559
00560 for( jj=0 ; jj < nr-1 ; jj++ ) /* 30 May 1999: used to have nr-2 */
00561 base += ff[jj] * bb[jj] ; /* which was wrong! */
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 }
|
|
||||||||||||
|
Definition at line 424 of file afni_pcor.c. References PCOR_voxel_corr::nvox, PCOR_get_variance(), and vc. Referenced by AFNI_fimmer_compute().
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 }
|
|
||||||||||||
|
Definition at line 400 of file afni_pcor.c. References nr, PCOR_voxel_corr::nref, PCOR_voxel_corr::nupdate, PCOR_voxel_corr::nvox, scale, var, vc, and VCH. Referenced by PCOR_get_stdev().
00401 {
00402 int vox , nv = vc->nvox , nr = vc->nref , nup = vc->nupdate ;
00403 float scale ;
00404
00405 /*** check inputs for OK-ness ***/
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 /*** Work ***/
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 }
|
|
||||||||||||||||
|
|
|
||||||||||||||||
|
Referenced by AFNI_fimmer_compute(), CORREL_main(), fim3d_fimmer_compute(), and RT_fim_recurse(). |
|
||||||||||||||||
|
|
|
||||||||||||
|
Definition at line 103 of file afni_pcor.c. References PCOR_references::alp, PCOR_references::betasq, EXIT, PCOR_references::ff, free, PCOR_references::gg, malloc, nr, PCOR_references::nref, PCOR_references::nupdate, RCH, ref, PCOR_references::rmax, PCOR_references::rmin, PCOR_references::rsum, and vec. Referenced by AFNI_fimmer_compute(), CORREL_main(), fim3d_fimmer_compute(), and RT_fim_recurse().
00104 {
00105 int nr = ref->nref , jj,kk ;
00106 float bold , bnew , aaa,fff,ggg ;
00107
00108 static float * zz = NULL ; /* local storage for a copy of vec */
00109 static int zz_size = -1 ;
00110
00111 /*** copy vector data into local storage (will be altered below) ***/
00112
00113 if( zz_size < nr ){ /* get new space, if not enough is present */
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 /*** 14 Jan 1998: collect ref vector stats ***/
00127
00128 if( ref->nupdate == 0 ){ /* initialize */
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 /*** Carlson algorithm ***/
00142
00143 bold = 1.0 ;
00144
00145 for( jj=0 ; jj < nr ; jj++ ){
00146
00147 aaa = zz[jj] / RCH(ref,jj,jj) ; /* alpha */
00148 bnew = sqrt( bold*bold + aaa*aaa ) ; /* new beta */
00149 fff = bnew / bold ; /* f factor */
00150 ggg = aaa / (bnew*bold) ; /* g factor */
00151 bold = bnew ; /* new beta becomes old beta */
00152
00153 ref->alp[jj] = aaa ; /* save these for later use */
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 ) ; /* and save this too! */
00164
00165 (ref->nupdate)++ ; /* June 1995: another update! */
00166 return ;
00167 }
|
