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 #ifndef _AFNI_PCOR_HEADER_
00008 #define _AFNI_PCOR_HEADER_
00009 
00010 /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00011      header file for recursive partial correlation calculations
00012      -- RWCox, Feb 1994
00013      -- March 1995: this work is the basis for the paper
00014          "Real-Time Functional Magnetic Resonance Imaging"
00015          by RW Cox, A Jesmanowicz, and JS Hyde,
00016          Magnetic Resonance in Medicine, 33:230-236.
00017      -- June 1995: added get_variance and some comments
00018      -- June 1996: modified to fit into AFNI
00019      -- Jan  2000: added get_stdev
00020 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00021 
00022 #include <stdlib.h>
00023 #include <stdio.h>
00024 #include <math.h>
00025 #include "mrilib.h"
00026 
00027 /*** some macros ***/
00028 
00029 /** combine two interpreted tokens into one using TWO_TWO **/
00030 
00031 #ifndef TWO_TWO
00032 #  define TWO_ONE(x,y) x ## y
00033 #  define TWO_TWO(x,y) TWO_ONE(x,y)
00034 #endif
00035 
00036 #ifndef  SQR
00037 #  define SQR(x)  ((x)*(x))
00038 #endif
00039 
00040 #ifndef MAX
00041 #  define MAX(x,y) (((x)>(y)) ? (x) : (y))
00042 #endif
00043 
00044 #ifndef MIN
00045 #  define MIN(x,y) (((x)<(y)) ? (x) : (y))
00046 #endif
00047 
00048 /*** default is to expand the inner loop in update_voxel_corr for small
00049      values of nref;  if this is not desired, comment the next line out ***/
00050 
00051 #define EXPAND_UPDATE
00052 
00053 #define REF_EPS  1.0e-7
00054 
00055 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00056 
00057 /*** this type holds the independent-of-voxel references information ***/
00058 
00059    /* note that the partial correlation coefficient of the data vectors
00060       is calculated with respect to the LAST reference vector,
00061       removing the effects of all the previous reference vectors */
00062 
00063 typedef struct {
00064           int nref ;                /* number of ref vectors */
00065           int nupdate ;             /* number of times has been updated */
00066           float **chol ,            /* nref X nref Cholesky factor */
00067                 *alp , *ff , *gg ;  /* alpha, f, and g factors */
00068           float *rmin , *rmax ,     /* 14 Jan 1998 additions: */
00069                 *rsum ;             /* statistics of references */
00070           float betasq ;            /* 1/beta^2 factor */
00071    } PCOR_references ;
00072 
00073  /*** RCH is for access to the (ii,jj) element
00074       of the Cholesky factor of a given set of references;
00075       ii >= jj, since this is the lower triangular factor
00076 
00077       N.B.: compared to the paper's notation,
00078             L+1     = nref
00079             c_{i,j} = RCH(rr,i-1,j-1)  for i=1..L+1 , j=1..i,
00080         and c_{i,j} = VCH(vv,vox,j-1)  for i=L+2, j=1..L+2    ***/
00081 
00082 #define RCH(rr,ii,jj)  (rr->chol[(ii)][(jj)])
00083 
00084 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00085 
00086 /*** this type holds all the voxel-specific data ***/
00087 
00088 typedef struct {
00089           int nvox ;      /* number of voxels */
00090           int nref ;      /* number of references to allow for */
00091           int nupdate ;   /* number of times it has been updated */
00092           float *chrow ;  /* last row (length nref+1)
00093                              of Cholesky factor for each voxel
00094                              (N.B.:  last element is actually squared) */
00095    } PCOR_voxel_corr ;
00096 
00097    /* VCH is for access to the jj-th element of the last row of
00098       the Cholesky factor for the vox-th voxel in a given voxel_corr struct;
00099 
00100       N.B.:  the last element [VCH(vv,vox,nref)] is the SQUARE of
00101              the actual Cholesky diagonal element, since that is all that the
00102              partial correlation coefficient algorithm actually needs;
00103              this prevents taking a square root at each time step for each voxel */
00104 
00105 #define VCH(vv,vox,jj) (vv->chrow[(jj)+(vox)*(vv->nref+1)])
00106 
00107 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00108 
00109 /*** prototypes ***/
00110 
00111 extern PCOR_references * new_PCOR_references(int) ;
00112 extern PCOR_voxel_corr * new_PCOR_voxel_corr(int,int) ;
00113 
00114 extern void free_PCOR_voxel_corr(PCOR_voxel_corr *) ;
00115 extern void free_PCOR_references(PCOR_references *) ;
00116 
00117 extern void PCOR_get_lsqfit(PCOR_references *, PCOR_voxel_corr *, float *fit[]) ;
00118 extern void PCOR_get_coef  (PCOR_references *, PCOR_voxel_corr *, float *) ;
00119 extern void PCOR_get_pcor  (PCOR_references *, PCOR_voxel_corr *, float *) ;
00120 extern void PCOR_get_mcor  (PCOR_references *, PCOR_voxel_corr *, int , float *) ;
00121 extern void PCOR_get_perc  (PCOR_references *, PCOR_voxel_corr *, float *,float *, int );
00122 
00123 extern void PCOR_get_pcor_and_coef(PCOR_references *, PCOR_voxel_corr *,
00124                                    float, float *,float * ) ;
00125 
00126 extern void PCOR_get_variance(PCOR_voxel_corr *, float *) ;
00127 extern void PCOR_get_stdev   (PCOR_voxel_corr *, float *) ;  /* 03 Jan 2000 */
00128 
00129 extern void update_PCOR_references(float *,PCOR_references *) ;
00130 
00131 extern void PCOR_update_short( short * , PCOR_references * , PCOR_voxel_corr * ) ;
00132 extern void PCOR_update_float( float * , PCOR_references * , PCOR_voxel_corr * ) ;
00133 extern void PCOR_update_byte ( byte  * , PCOR_references * , PCOR_voxel_corr * ) ;
00134 
00135 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00136 #endif