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 "mrilib.h"
00008 
00009 /*-------------------------------------------------------------------------*/
00010 
00011 void THD_const_detrend( int npt, float *xx, float *xx0 ) /* 24 Aug 2001 */
00012 {
00013    int ii ; float xbar ;
00014 
00015    if( npt < 2 || xx == NULL ) return ;
00016 
00017    xbar = 0.0 ;
00018    for( ii=0 ; ii < npt ; ii++ ) xbar += xx[ii] ;
00019    xbar /= npt ;
00020    for( ii=0 ; ii < npt ; ii++ ) xx[ii] -= xbar ;
00021    if( xx0 != NULL ) *xx0 = xbar ;
00022    return ;
00023 }
00024 
00025 /*-------------------------------------------------------------------------*/
00026 
00027 void get_linear_trend( int npt, float *xx, float *f0, float *f1 )
00028 {
00029    double t1,t3,t10 , x0,x1 ;
00030    int ii ;
00031 
00032    if( npt < 2 || xx == NULL || f0 == NULL || f1 == NULL ) return ;
00033 
00034    x0 = xx[0] ; x1 = 0.0 ;
00035    for( ii=1 ; ii < npt ; ii++ ){
00036       x0 += xx[ii] ;
00037       x1 += xx[ii] * ii ;
00038    }
00039 
00040    t1 = npt*x0; t3 = 1.0/npt; t10 = npt*npt;
00041 
00042    *f0 = (float)(2.0/(npt+1.0)*t3*(2.0*t1-3.0*x1-x0));
00043    *f1 = (float)(-6.0/(t10-1.0)*t3*(-x0-2.0*x1+t1));
00044    return ;
00045 }
00046 
00047 /*-------------------------------------------------------------------------
00048    Linear detrend a 1D float array in place.
00049    If xx0 != NULL and xx1 != NULL, then the trend removed is
00050       far[i] -= (*xx0) + (*xx1) * i, for i=0..npt-1
00051 ---------------------------------------------------------------------------*/
00052 
00053 void THD_linear_detrend( int npt, float *far, float *xx0, float *xx1 )
00054 {
00055    register int ii ;
00056    float f0 , f1 ;
00057 
00058    if( npt < 3 || far == NULL ) return ;
00059 
00060    get_linear_trend( npt , far , &f0 , &f1 ) ;
00061 
00062    far[0] -= f0 ;
00063    for( ii=1 ; ii < npt ; ii++ ) far[ii] -= (f0 + f1*ii) ;
00064 
00065    if( xx0 != NULL ) *xx0 = f0 ;
00066    if( xx1 != NULL ) *xx1 = f1 ;
00067 
00068    return ;
00069 }
00070 
00071 /*---------------------------------------------------------------------------
00072    Given x[0..npt-1], return f0,f1,f2 as the least squares coefficients to
00073      x[j] = f0 + f1*j + f2*j*j
00074 -----------------------------------------------------------------------------*/
00075 
00076 void get_quadratic_trend( int npt, float *xx, float *f0, float *f1, float *f2 )
00077 {
00078    double  x0,x1,x2 , N=npt ;
00079    int ii ;
00080 
00081    if( npt < 3 || xx == NULL || f0 == NULL || f1 == NULL || f2 == NULL ) return;
00082 
00083    x0 = xx[0] ; x1 = x2 = 0.0 ;
00084    for( ii=1 ; ii < npt ; ii++ ){
00085      x0 +=  xx[ii] ;
00086      x1 += (xx[ii] * ii) ;
00087      x2 += (xx[ii] * ii) * ii ;
00088    }
00089 
00090    *f0 = (  3.0*(3.0*N*N-3.0*N+2.0) * x0
00091           -18.0*(-1.0+2.0*N)        * x1
00092           +30.0                     * x2 ) / (N*(N+2.0)*(N+1.0)) ;
00093 
00094    *f1 = ( -18.0*(-1.0+2.0*N)              * x0
00095            +12.0*(-1.0+2.0*N)*(8.0*N-11.0) * x1 /((N-1.0)*(N-2.0))
00096            -180.0                          * x2 /(N-2.0)          )
00097         / (N*(N+2.0)*(N+1.0)) ;
00098 
00099    *f2 = ( 30.0  * x0
00100           -180.0 * x1 / (N-2.0)
00101           +180.0 * x2 / ((N-1.0)*(N-2.0)) ) / (N*(N+2.0)*(N+1.0))  ;
00102    return ;
00103 }
00104 
00105 /*-------------------------------------------------------------------------
00106    Quadratic detrend a 1D float array in place.
00107    If xx0 != NULL, xx1 != NULL, xx2 != NULL, then the trend removed is
00108       far[i] -= (*xx0) + (*xx1) * i + (*xx2)*(i*i) , for i=0..npt-1
00109 ---------------------------------------------------------------------------*/
00110 
00111 void THD_quadratic_detrend( int npt, float *far,
00112                             float *xx0, float *xx1, float *xx2 )
00113 {
00114    register int ii ;
00115    float f0 , f1 , f2 ;
00116 
00117    if( npt < 4 || far == NULL ) return ;
00118 
00119    get_quadratic_trend( npt , far , &f0 , &f1 , &f2 ) ;
00120 
00121    far[0] -= f0 ;
00122    for( ii=1 ; ii < npt ; ii++ ) far[ii] -= ( (f2*ii + f1)*ii + f0 ) ;
00123 
00124    if( xx0 != NULL ) *xx0 = f0 ;
00125    if( xx1 != NULL ) *xx1 = f1 ;
00126    if( xx2 != NULL ) *xx2 = f2 ;
00127 
00128    return ;
00129 }
00130 
00131 /*------------------------------------------------------------------------------*/
00132 /*! Cubic detrend a float array in place. */
00133 
00134 void THD_cubic_detrend( int npt , float *far )  /* 15 Nov 1999 */
00135 {
00136    register int ii ;
00137    float g0,g1,g2,g3 , f0,f1,f2,f3 , t1,t2,t5,t8 , t95,t56,t22,t25,txx ;
00138 
00139    if( npt < 5 || far == NULL ) return ;
00140 
00141    t8 = npt*npt ; t2 = npt-1.0 ; t5 = t2*(npt-2.0) ;
00142    t95 = 0.05*t5*(npt-3.0) ;
00143    t56 = 0.16666667*t5 ;
00144    t22 = 0.5*t2 ;
00145    t25 = 1.5*t2 ;
00146    txx = 0.6*t8-1.5*npt+1.1 ;
00147 
00148    g0=g1=g2=g3=0.0 ;
00149    for( ii=0 ; ii < npt ; ii++ ){
00150       t1 = ii*ii ;
00151       f1 = ii - t22 ;
00152       f2 = t1 - t2*ii + t56 ;
00153       f3 = t1*(ii - t25) + txx*ii - t95 ;
00154 
00155       g0 += far[ii] ;
00156       g1 += far[ii] * f1 ;
00157       g2 += far[ii] * f2 ;
00158       g3 += far[ii] * f3 ;
00159    }
00160    g0 *= (1.0/npt) ;
00161    g1 *= 12.0/(npt*(t8-1.0)) ;
00162    g2 *= 180.0/(npt*(t8-1.0)*(t8-4.0)) ;
00163    g3 *= 2800.0/(npt*(t8-1.0)*(t8-4.0)*(t8-9.0)) ;
00164 
00165    for( ii=0 ; ii < npt ; ii++ ){
00166       t1 = ii*ii ;
00167       f1 = ii- t22 ;
00168       f2 = t1 - t2*ii + t56 ;
00169       f3 = t1*(ii - t25) + txx*ii - t95 ;
00170 
00171       far[ii] -= ( g0 + g1*f1 + g2*f2 + g3*f3 ) ;
00172    }
00173 
00174    return ;
00175 }
00176 
00177 /*-------------------------------------------------------------------------
00178    Make a vector have L2 norm 1
00179 ---------------------------------------------------------------------------*/
00180 
00181 void THD_normalize( int npt , float *far )
00182 {
00183    register int ii ;
00184    register float fac ;
00185 
00186    if( npt <= 0 || far == NULL ) return ;
00187 
00188    fac = 0 ;
00189    for( ii=0 ; ii < npt ; ii++ ) fac += far[ii]*far[ii] ;
00190    if( fac == 0.0 ) return ;
00191    fac = 1.0 / sqrt(fac) ;
00192    for( ii=0 ; ii < npt ; ii++ ) far[ii] /= fac ;
00193    return ;
00194 }
00195 
00196 /*-----------------------------------------------------------------------*/
00197 /*! Detrend a vector with a given polort level, plus some others.
00198 -------------------------------------------------------------------------*/
00199 
00200 void THD_generic_detrend( int npt, float *far ,
00201                           int polort, int nort, float **ort )
00202 {
00203    int ii,jj , nref ;
00204    float **ref , *fit , xmid , xfac , val ;
00205 
00206    /* check inputs */
00207 
00208    if( npt <= 1 || far == NULL ) return ;
00209    if( nort > 0 ){
00210      if( ort == NULL ) return ;
00211      for( jj=0 ; jj < nort ; jj++ ) if( ort[jj] == NULL ) return ;
00212    }
00213    if( polort <  0 ) polort = -1 ;
00214    if( nort   <  0 ) nort   =  0 ;
00215 
00216    nref = polort+1+nort ;
00217    if( nref == 0 || nref >= npt-1 ) return ;
00218 
00219    ref  = (float **) malloc( sizeof(float *) * nref ) ;
00220    xmid = 0.5*(npt-1) ; xfac = 1.0 / xmid ;
00221    for( jj=0 ; jj <= polort ; jj++ ){
00222      ref[jj] = (float *) malloc( sizeof(float) * npt ) ;
00223      switch( jj ){
00224        case 0:
00225          for( ii=0 ; ii < npt ; ii++ ) ref[jj][ii] = 1.0 ;
00226        break ;
00227 
00228        case 1:
00229          for( ii=0 ; ii < npt ; ii++ ) ref[jj][ii] = xfac*(ii-xmid) ;
00230        break ;
00231 
00232        case 2:
00233          for( ii=0 ; ii < npt ; ii++ ){
00234            val = xfac*(ii-xmid) ; ref[jj][ii] = val*val ;
00235          }
00236        break ;
00237 
00238        case 3:
00239          for( ii=0 ; ii < npt ; ii++ ){
00240            val = xfac*(ii-xmid) ; ref[jj][ii] = val*val*val ;
00241          }
00242        break ;
00243 
00244        default:
00245          for( ii=0 ; ii < npt ; ii++ ){
00246            val = xfac*(ii-xmid) ; ref[jj][ii] = pow(val,(double)(jj)) ;
00247          }
00248        break ;
00249      }
00250    }
00251    for( jj=0 ; jj < nort ; jj++ )   /* user supplied refs */
00252      ref[polort+1+jj] = ort[jj] ;
00253 
00254    fit = lsqfit( npt , far , NULL , nref , ref ) ;
00255 
00256    if( fit != NULL ){                                   /* good */
00257      for( ii=0 ; ii < npt ; ii++ ){
00258        val = far[ii] ;
00259        for( jj=0 ; jj < nref ; jj++ )
00260          val -= fit[jj] * ref[jj][ii] ;
00261        far[ii] = val ;
00262      }
00263      free(fit) ;
00264    } else {
00265      fprintf(stderr,"THD_generic_detrend: lsqfit fails - no detrending!\n") ;
00266    }
00267 
00268    for( jj=0 ; jj <= polort ; jj++ ) free(ref[jj]) ;
00269    free(ref) ; return ;
00270 }