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