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 /***
00008      header file for recursive partial correlation calculations
00009      -- RWCox, Feb 1994
00010      -- March 1995: this work is the basis for the paper
00011          "Real-Time Functional Magnetic Resonance Imaging"
00012          by RW Cox, A Jesmanowicz, and JS Hyde,
00013          Magnetic Resonance in Medicine, 33:230-236.
00014      -- June 1995: added get_variance and some comments
00015 ***/
00016 
00017 #ifndef PCOR_HEADER
00018 #define PCOR_HEADER
00019 
00020 #include <stdlib.h>
00021 #include <stdio.h>
00022 #include <math.h>
00023 
00024 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00025 
00026 /*** some macros ***/
00027 
00028 #ifndef  SQR
00029 #define SQR(x)  ((x)*(x))
00030 #endif
00031 
00032 #ifndef MAX
00033 #  define MAX(x,y) (((x)>(y)) ? (x) : (y))
00034 #endif
00035 
00036 #ifndef MIN
00037 #  define MIN(x,y) (((x)<(y)) ? (x) : (y))
00038 #endif
00039 
00040 /*** default is to expand the inner loop in update_voxel_corr for
00041      small values of nref;  if this is not desired, compile with
00042      the switch -DNO_EXPAND_UPDATE, or comment the next line out ***/
00043 
00044 #define EXPAND_UPDATE
00045 
00046 #ifdef NO_EXPAND_UPDATE
00047 #undef EXPAND_UPDATE
00048 #endif
00049 
00050 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00051 
00052 /*** Define atomic data types ***/
00053 
00054 #ifndef REF_FLOAT_SINGLE           /* compile time option to save space */
00055    typedef double ref_float ;      /* internal storage of reference data */
00056 #  define REF_EPS  1.0e-13         /* a small number of this type */
00057 #else
00058    typedef float ref_float ;
00059 #  define REF_EPS  1.0e-7
00060 #endif  /* REF_FLOAT_SINGLE */
00061 
00062 /*-------------------------------*/
00063 
00064 #ifndef VOX_SHORT                  /* define voxel data as shorts or floats */
00065    typedef float vox_data ;
00066 #else
00067    typedef short vox_data ;
00068 #endif  /* VOX_SHORT */
00069 
00070 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00071 
00072 /*** this type holds the independent-of-voxel references information ***/
00073 
00074    /* note that the partial correlation coefficient of the data vectors
00075       is calculated with respect to the LAST reference vector,
00076       removing the effects of all the previous reference vectors */
00077 
00078 typedef struct {
00079           int nref ;                    /* number of ref vectors */
00080           int nupdate ;                 /* number of times has been updated */
00081           ref_float **chol ,            /* nref X nref Cholesky factor */
00082                     *alp , *ff , *gg ;  /* alpha, f, and g factors */
00083           ref_float betasq ;            /* 1/beta^2 factor */
00084    } references ;
00085 
00086  /*** RCH is for access to the (ii,jj) element
00087       of the Cholesky factor of a given set of references;
00088       ii >= jj, since this is the lower triangular factor
00089 
00090       N.B.: compared to the paper's notation,
00091             L+1     = nref
00092             c_{i,j} = RCH(rr,i-1,j-1)  for i=1..L+1 , j=1..i,
00093         and c_{i,j} = VCH(vv,vox,j-1)  for i=L+2, j=1..L+2    ***/
00094 
00095 #define RCH(rr,ii,jj)  (rr->chol[(ii)][(jj)])
00096 
00097 #ifdef OV_DEBUG1
00098 #define REF_DUMP(rr,str) \
00099   {  int iq,jq ; ref_float qsum ; \
00100      fprintf(stderr,"%s: reference dump, nref=%d betasq=%11.4e\n", \
00101              str , rr->nref , rr->betasq ) ; \
00102      for( iq=0 ; iq < rr->nref ; iq++){ \
00103        fprintf(stderr," ROW %d: ",iq) ; \
00104        qsum = 0.0 ; \
00105        for( jq=0 ; jq <= iq ; jq++ ){ \
00106           fprintf(stderr,"%11.4e ",RCH(rr,iq,jq)); \
00107           qsum += SQR(RCH(rr,iq,jq)) ; \
00108        } \
00109        fprintf(stderr,": qsum=%11.4e\n",qsum) ; \
00110        fprintf(stderr,"      alpha=%11.4e  ff=%11.4e  gg=%11.4e\n", \
00111                rr->alp[iq] , rr->ff[iq] , rr->gg[iq] ) ; \
00112      } \
00113    }
00114 #endif
00115 
00116 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00117 
00118 /*** this type holds all the voxel-specific data ***/
00119 
00120 typedef struct {
00121           int nvox ;          /* number of voxels */
00122           int nref ;          /* number of references to allow for */
00123           int nupdate ;       /* number of times it has been updated */
00124           ref_float *chrow ;  /* last row (length nref+1)
00125                                  of Cholesky factor for each voxel
00126                                  (N.B.:  last element is actually squared) */
00127    } voxel_corr ;
00128 
00129    /* VCH is for access to the jj-th element of the last row of
00130       the Cholesky factor for the vox-th voxel in a given voxel_corr struct;
00131 
00132       N.B.:  the last element [VCH(vv,vox,nref)] is the SQUARE of
00133              the actual Cholesky diagonal element, since that is all that the
00134              partial correlation coefficient algorithm actually needs;
00135              this prevents taking a square root at each time step for each voxel */
00136 
00137 #define VCH(vv,vox,jj) (vv->chrow[(jj)+(vox)*(vv->nref+1)])
00138 
00139 #ifdef OV_DEBUG1
00140 #define VD 2001
00141 #define VOX_DUMP(vv,vox,str) \
00142    {  int jq ; ref_float qsum = 0.0 ; \
00143       fprintf(stderr,"%s: voxel_corr dump #%d\n  ",str,vox) ; \
00144       for( jq=0 ; jq < vv->nref ; jq++ ){ \
00145         fprintf(stderr,"%11.4e ",VCH(vv,vox,jq)) ; \
00146         qsum += SQR(VCH(vv,vox,jq)) ; \
00147       } \
00148       qsum += VCH(vv,vox,vv->nref) ; \
00149       fprintf(stderr,"%11.4e : qsum=%11.4e\n",VCH(vv,vox,vv->nref),qsum); \
00150    }
00151 #endif
00152 
00153 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00154 
00155 /*** this type holds the return data for thresholding results ***/
00156 
00157 typedef struct {
00158           int   num_pcor_pos , num_pcor_neg , num_coef_pos , num_coef_neg ;
00159           float max_pcor_pos , max_pcor_neg , max_coef_pos , max_coef_neg ;
00160    } thresh_result ;
00161 
00162 #ifdef OV_DEBUG1
00163 #define THR_DUMP(thr,str) \
00164 fprintf(stderr,"thresh_results dump for %s\n:") ; \
00165 fprintf(stderr," num_pcor_pos=%d neg=%d  num_coef_pos=%d neg=%d\n", \
00166  (thr).num_pcor_pos,(thr).num_pcor_neg,(thr).num_coef_pos,(thr).num_coef_neg ) ; \
00167 fprintf(stderr," max_pcor_pos=%11.4g neg=%11.4g  max_coef_pos=%11.4g neg=%11.4g\n", \
00168  (thr).max_pcor_pos,(thr).max_pcor_neg,(thr).max_coef_pos,(thr).max_coef_neg ) ;
00169 #endif
00170 
00171 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00172 
00173 /*** prototypes ***/
00174 
00175 extern references * new_references() ;
00176 
00177 extern void update_references() ;
00178 
00179 extern voxel_corr * new_voxel_corr() ;
00180 
00181 extern void free_voxel_corr() ;
00182 
00183 extern void free_references() ;
00184 
00185 extern void update_voxel_corr() ;
00186 
00187 extern void get_pcor() ;
00188 
00189 extern void get_coef() ;
00190 
00191 extern void get_pcor_thresh_coef() ;
00192 
00193 extern void get_variances() ;
00194 
00195 extern void get_lsqfit() ;
00196 
00197 #endif