1 /***************************************************************************** 2 Major portions of this software are copyrighted by the Medical College 3 of Wisconsin, 1994-2000, and are released under the Gnu General Public 4 License, Version 2. See the file README.Copyright for details. 5 ******************************************************************************/ 6 7 /*** 8 header file for recursive partial correlation calculations 9 -- RWCox, Feb 1994 10 -- March 1995: this work is the basis for the paper 11 "Real-Time Functional Magnetic Resonance Imaging" 12 by RW Cox, A Jesmanowicz, and JS Hyde, 13 Magnetic Resonance in Medicine, 33:230-236. 14 -- June 1995: added get_variance and some comments 15 ***/ 16 17 #ifndef PCOR_HEADER 18 #define PCOR_HEADER 19 20 #include <stdlib.h> 21 #include <stdio.h> 22 #include <math.h> 23 24 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ 25 26 /*** some macros ***/ 27 28 #ifndef SQR 29 #define SQR(x) ((x)*(x)) 30 #endif 31 32 #ifndef MAX 33 # define MAX(x,y) (((x)>(y)) ? (x) : (y)) 34 #endif 35 36 #ifndef MIN 37 # define MIN(x,y) (((x)<(y)) ? (x) : (y)) 38 #endif 39 40 /*** default is to expand the inner loop in update_voxel_corr for 41 small values of nref; if this is not desired, compile with 42 the switch -DNO_EXPAND_UPDATE, or comment the next line out ***/ 43 44 #define EXPAND_UPDATE 45 46 #ifdef NO_EXPAND_UPDATE 47 #undef EXPAND_UPDATE 48 #endif 49 50 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ 51 52 /*** Define atomic data types ***/ 53 54 #ifndef REF_FLOAT_SINGLE /* compile time option to save space */ 55 typedef double ref_float ; /* internal storage of reference data */ 56 # define REF_EPS 1.0e-13 /* a small number of this type */ 57 #else 58 typedef float ref_float ; 59 # define REF_EPS 1.0e-7 60 #endif /* REF_FLOAT_SINGLE */ 61 62 /*-------------------------------*/ 63 64 #ifndef VOX_SHORT /* define voxel data as shorts or floats */ 65 typedef float vox_data ; 66 #else 67 typedef short vox_data ; 68 #endif /* VOX_SHORT */ 69 70 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ 71 72 /*** this type holds the independent-of-voxel references information ***/ 73 74 /* note that the partial correlation coefficient of the data vectors 75 is calculated with respect to the LAST reference vector, 76 removing the effects of all the previous reference vectors */ 77 78 typedef struct { 79 int nref ; /* number of ref vectors */ 80 int nupdate ; /* number of times has been updated */ 81 ref_float **chol , /* nref X nref Cholesky factor */ 82 *alp , *ff , *gg ; /* alpha, f, and g factors */ 83 ref_float betasq ; /* 1/beta^2 factor */ 84 } references ; 85 86 /*** RCH is for access to the (ii,jj) element 87 of the Cholesky factor of a given set of references; 88 ii >= jj, since this is the lower triangular factor 89 90 N.B.: compared to the paper's notation, 91 L+1 = nref 92 c_{i,j} = RCH(rr,i-1,j-1) for i=1..L+1 , j=1..i, 93 and c_{i,j} = VCH(vv,vox,j-1) for i=L+2, j=1..L+2 ***/ 94 95 #define RCH(rr,ii,jj) (rr->chol[(ii)][(jj)]) 96 97 #ifdef OV_DEBUG1 98 #define REF_DUMP(rr,str) \ 99 { int iq,jq ; ref_float qsum ; \ 100 fprintf(stderr,"%s: reference dump, nref=%d betasq=%11.4e\n", \ 101 str , rr->nref , rr->betasq ) ; \ 102 for( iq=0 ; iq < rr->nref ; iq++){ \ 103 fprintf(stderr," ROW %d: ",iq) ; \ 104 qsum = 0.0 ; \ 105 for( jq=0 ; jq <= iq ; jq++ ){ \ 106 fprintf(stderr,"%11.4e ",RCH(rr,iq,jq)); \ 107 qsum += SQR(RCH(rr,iq,jq)) ; \ 108 } \ 109 fprintf(stderr,": qsum=%11.4e\n",qsum) ; \ 110 fprintf(stderr," alpha=%11.4e ff=%11.4e gg=%11.4e\n", \ 111 rr->alp[iq] , rr->ff[iq] , rr->gg[iq] ) ; \ 112 } \ 113 } 114 #endif 115 116 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ 117 118 /*** this type holds all the voxel-specific data ***/ 119 120 typedef struct { 121 int nvox ; /* number of voxels */ 122 int nref ; /* number of references to allow for */ 123 int nupdate ; /* number of times it has been updated */ 124 ref_float *chrow ; /* last row (length nref+1) 125 of Cholesky factor for each voxel 126 (N.B.: last element is actually squared) */ 127 } voxel_corr ; 128 129 /* VCH is for access to the jj-th element of the last row of 130 the Cholesky factor for the vox-th voxel in a given voxel_corr struct; 131 132 N.B.: the last element [VCH(vv,vox,nref)] is the SQUARE of 133 the actual Cholesky diagonal element, since that is all that the 134 partial correlation coefficient algorithm actually needs; 135 this prevents taking a square root at each time step for each voxel */ 136 137 #define VCH(vv,vox,jj) (vv->chrow[(jj)+(vox)*(vv->nref+1)]) 138 139 #ifdef OV_DEBUG1 140 #define VD 2001 141 #define VOX_DUMP(vv,vox,str) \ 142 { int jq ; ref_float qsum = 0.0 ; \ 143 fprintf(stderr,"%s: voxel_corr dump #%d\n ",str,vox) ; \ 144 for( jq=0 ; jq < vv->nref ; jq++ ){ \ 145 fprintf(stderr,"%11.4e ",VCH(vv,vox,jq)) ; \ 146 qsum += SQR(VCH(vv,vox,jq)) ; \ 147 } \ 148 qsum += VCH(vv,vox,vv->nref) ; \ 149 fprintf(stderr,"%11.4e : qsum=%11.4e\n",VCH(vv,vox,vv->nref),qsum); \ 150 } 151 #endif 152 153 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ 154 155 /*** this type holds the return data for thresholding results ***/ 156 157 typedef struct { 158 int num_pcor_pos , num_pcor_neg , num_coef_pos , num_coef_neg ; 159 float max_pcor_pos , max_pcor_neg , max_coef_pos , max_coef_neg ; 160 } thresh_result ; 161 162 #ifdef OV_DEBUG1 163 #define THR_DUMP(thr,str) \ 164 fprintf(stderr,"thresh_results dump for %s\n:") ; \ 165 fprintf(stderr," num_pcor_pos=%d neg=%d num_coef_pos=%d neg=%d\n", \ 166 (thr).num_pcor_pos,(thr).num_pcor_neg,(thr).num_coef_pos,(thr).num_coef_neg ) ; \ 167 fprintf(stderr," max_pcor_pos=%11.4g neg=%11.4g max_coef_pos=%11.4g neg=%11.4g\n", \ 168 (thr).max_pcor_pos,(thr).max_pcor_neg,(thr).max_coef_pos,(thr).max_coef_neg ) ; 169 #endif 170 171 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ 172 173 /*** prototypes ***/ 174 175 extern references * new_references() ; 176 177 extern void update_references() ; 178 179 extern voxel_corr * new_voxel_corr() ; 180 181 extern void free_voxel_corr() ; 182 183 extern void free_references() ; 184 185 extern void update_voxel_corr() ; 186 187 extern void get_pcor() ; 188 189 extern void get_coef() ; 190 191 extern void get_pcor_thresh_coef() ; 192 193 extern void get_variances() ; 194 195 extern void get_lsqfit() ; 196 197 extern void get_variance() ; 198 199 #endif 200