00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
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      recursive calculation of partial correlation coefficients for
00018      a lot of voxels at once. -- RWCox, Feb 1994
00019 ***/
00020 
00021 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00022 
00023 /***
00024    create a new references data structure:
00025       input:  numref = number of reference vectors to allow for
00026       output: pointer to the data structure
00027 ***/
00028 
00029 references * new_references(numref)
00030      int numref;
00031 {
00032    references *ref ;
00033    int ii,jj ;
00034 
00035    /*** check input for reasonableness ***/
00036 
00037    if( numref < 1 ){
00038       fprintf( stderr , "new_references called with numref=%d\n" , numref ) ;
00039       exit(-1) ;
00040    }
00041 
00042    /*** allocate storage for top level data ***/
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 ;      /* June 1995: not updated at all yet */
00051 
00052    /*** allocate storage for Cholesky factor
00053         (an array of rows, row #ii is length ii+1, for ii=0..numref-1) ***/
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    /*** allocate storage for vectors of alpha, f, g ***/
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    /*** initialize Cholesky factor ***/
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    update a references structure:
00101       input:  vec = pointer to nref values of the reference vectors
00102                     at the new time point
00103               ref = references data structure
00104       output: ref is updated;
00105                 the Cholesky factor is modified via the Carlson algorithm,
00106                 the alpha, f, and g factors are saved for later use
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] ;  /* to hold the sums of squares of inputs */
00120 #endif
00121 
00122    /*** copy vector data into local storage ***/
00123 
00124    if( zz_size < nr ){   /* get new space, if not enough is present */
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]) ;  /* for later */
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    /*** Carlson algorithm ***/
00147 
00148    bold = 1.0 ;
00149 
00150    for( jj=0 ; jj < nr ; jj++ ){
00151 
00152       aaa  = zz[jj] / RCH(ref,jj,jj) ;        /* alpha */
00153       bnew = sqrt( bold*bold + aaa*aaa ) ;    /* new beta */
00154       fff  = bnew / bold ;                    /* f factor */
00155       ggg  = aaa  / (bnew*bold) ;             /* g factor */
00156       bold = bnew ;                           /* new beta becomes old beta */
00157 
00158       ref->alp[jj] = aaa ;   /* save these for later use */
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 ) ;  /* and save this too! */
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)++ ;  /* June 1995: another update! */
00178    return ;
00179 }
00180 
00181 /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00182 
00183 /***
00184    create a new voxel partial correlation data structure
00185       inputs:  numvox = number of voxels in image
00186                numref = number of reference vectors to allow for
00187       output:  pointer to voxel_corr data structure
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    /*** check input for OK-osity ***/
00198 
00199    if( numvox < 1 ){
00200       fprintf( stderr , "new_voxel_corr: numvox=%d\n" , numvox ) ;
00201       exit(-1) ;
00202    }
00203 
00204    /*** get the base storage ***/
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    /*** setup the references common to all voxels ***/
00213 
00214    vc->nvox    = numvox ;
00215    vc->nref    = numref ;
00216    vc->nupdate = 0 ;      /* June 1995: not updated at all yet */
00217 
00218    /*** setup the storage of the last row for each voxel ***/
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    /*** initialize each voxel ***/
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 /*** de-allocate a references data structure ***/
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 /*** update all voxels with a new array of data
00273       inputs:  vdata = array[nvox] of new data for each voxel
00274                ref   = pointer to references structure to use
00275                vc    = pointer to correlation data structure
00276       output:  updated vc
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 ,       /* loop indices */
00285        nv = vc->nvox ,  /* number of voxels */
00286        nr = vc->nref ;  /* number of references */
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    /*** check inputs for OK-ness ***/
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 /** innermost loop expansion is for speedup if nref is small, if enabled **/
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    /*** for each voxel ***/
00328 
00329    for( vox=0 ; vox < nv ; vox++ ){
00330 
00331       /*** update last row of each Cholesky factor ***/
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 ; /* square of true Cholesky diagonal */
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 /* EXPAND_UPDATE */
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)++ ;  /* June 1995: another update */
00452    return ;
00453 }
00454 
00455 /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00456 
00457 /*** compute the partial correlation coefficients ***/
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    /*** check inputs for OK-ness ***/
00469 
00470    if( vc->nref != ref->nref ){
00471       fprintf( stderr , "get_pcor: reference size mismatch!\n" ) ;
00472       exit(-1) ;
00473    }
00474 
00475    /*** Work ***/
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 /*** get activation coefficient ***/
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    /*** check inputs for OK-ness ***/
00509 
00510    if( vc->nref != ref->nref ){
00511       fprintf( stderr , "get_coef: reference size mismatch!\n" ) ;
00512       exit(-1) ;
00513    }
00514 
00515    /*** Work ***/
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 /*** get variance estimate (June 1995) ***/
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    /*** check inputs for OK-ness ***/
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    /*** Work ***/
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 /*** Get least squares fit coefficients (all of them, Frank).
00567      "fit" is an array of pointers to floats, of length nref.
00568      If fit[j] != NULL, then it points to an array of size nvox that
00569      will get the coefficient for reference #j (j=0..nref-1).  [June 1995] ***/
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    /*** check inputs for OK-ness ***/
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    /*** for each voxel, compute the nr fit coefficients (backwards) ***/
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 /*** get correlation and thresholded alpha:
00622 
00623      only |pcor| >= pcthresh will be computed;
00624      only those voxels will have coef computed;
00625      if cothresh > 0, then voxels whose |coef| is less than
00626        cothresh * max|coef| will be also be set to zero
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    /*** check inputs for OK-ness ***/
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) ;      /* for coef calculation */
00656    thfac   = SQR(pcthresh)/(1.0-SQR(pcthresh)) ;
00657    do_pcth = pcthresh <= 0.0 ;               /* whether to do these tests */
00658    do_coth = cothresh  > 0.0 ;
00659 
00660    /*** Compute pcor and coef, thresholded on pcthresh ***/
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 ){   /* fancy threshold test */
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 {                                    /* fails pcor thresh */
00682          pcor[vox] = coef[vox] = 0.0 ;
00683       }
00684 
00685    }
00686 
00687    nco_pos = npc_pos ;
00688    nco_neg = npc_neg ;
00689 
00690 /*** threshold coef on cothresh as well ***/
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 /*** load threshold output report ***/
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 }