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    Input:  far   = array of floats length nar
00011            shift = fractional shift
00012 
00013    Output: Newly malloc-ed array that is supposed to be like
00014            far[ii-shift] for ii=0..nar-1.
00015 
00016    Notes: * The shift is to the RIGHT (a good Republican, of course).
00017           * If nar==1, then the output is an array of length 1 that is
00018              just a copy of the input.
00019           * Otherwise, cubic interpolation is used.
00020 ------------------------------------------------------------------------*/
00021 
00022 #define GET_AS_BIG(name,type,dim)                                           \
00023    do{ if( (dim) > name ## _size ){                                         \
00024           if( name != NULL ) free(name) ;                                   \
00025           name = (type *) malloc( sizeof(type) * (dim) ) ;                  \
00026           if( name == NULL ){                                               \
00027              fprintf(stderr,"*** can't malloc shifter space\n"); EXIT(1); } \
00028           name ## _size = (dim) ; } } while(0)
00029 
00030 /* cubic interpolation polynomials */
00031 
00032 #define P_M1(x)  ((x)*(1.0-(x))*((x)-2.0))
00033 #define P_00(x)  (3.0*((x)+1.0)*((x)-1.0)*((x)-2.0))
00034 #define P_P1(x)  (3.0*(x)*((x)+1.0)*(2.0-(x)))
00035 #define P_P2(x)  ((x)*((x)+1.0)*((x)-1.0))
00036 #define SIXTH    0.1666667
00037 
00038 float * shifter( int nar , float * far , float shift )
00039 {
00040    int ii,jj , nup , nmid , ix ;
00041    float xx , wt_m1,wt_00,wt_p1,wt_p2 , fmin,fmax ;
00042    float * fnew ;
00043 
00044    static int fl_size  = 0 ;     /* workspace (will hang around between calls) */
00045    static float * fl = NULL ;
00046 
00047    /*-- sanity checks --*/
00048 
00049    if( nar <= 0 || far == NULL ) return NULL ;
00050 
00051    if( nar == 1 ){
00052       fnew = (float *) malloc( sizeof(float) ) ;
00053       if( fnew == NULL ){
00054           fprintf(stderr,"*** can't malloc shifter output\n"); EXIT(1);
00055       }
00056      *fnew = far[0] ;
00057       return fnew ;
00058    }
00059 
00060    /*-- get workspace --*/
00061 
00062    nup = nar + (int)( 2.0*fabs(shift) + 6.0 ) ;
00063    GET_AS_BIG(fl,float,nup) ;
00064 
00065    /*-- Insert data:
00066           far[0] --> fl[nmid], etc.;
00067           fl[] points before nmid are copies of far[0];
00068           fl[] points after nmid+nar-1 are copiex of far[nar-1]. --*/
00069 
00070    nmid  = (nup-nar) / 2 ;
00071    for( ii=0 ; ii < nup ; ii++ ){
00072       jj = ii - nmid ;
00073       if( jj < 0 ) jj = 0 ; else if( jj >= nar ) jj = nar-1 ;
00074       fl[ii] = far[jj] ;
00075    }
00076 
00077    /*-- put results into output array --*/
00078 
00079    fnew = (float *) malloc( sizeof(float) * nar ) ;
00080    if( fnew == NULL ){
00081        fprintf(stderr,"*** can't malloc shifter output\n"); EXIT(1);
00082    }
00083 
00084    fmax = fmin = far[0] ;          /* find min and max of input */
00085    for( ii=1 ; ii < nar ; ii++ ){  /* for "clipping" purposes   */
00086       fmax = MAX(fmax,far[ii]) ;
00087       fmin = MIN(fmin,far[ii]) ;
00088    }
00089 
00090    for( ii=0 ; ii < nar ; ii++ ){
00091       xx = ii+nmid - shift ;  /* "index" in fl we want */
00092       ix = (int) xx ; xx = xx - ix ;
00093       wt_m1 = P_M1(xx) ; wt_00 = P_00(xx) ;
00094       wt_p1 = P_P1(xx) ; wt_p2 = P_P2(xx) ;
00095       fnew[ii] = SIXTH * (  wt_m1 * fl[ix-1] + wt_00 * fl[ix]
00096                           + wt_p1 * fl[ix+1] + wt_p2 * fl[ix+2] ) ;
00097 
00098            if( fnew[ii] < fmin ) fnew[ii] = fmin ;
00099       else if( fnew[ii] > fmax ) fnew[ii] = fmax ;
00100    }
00101 
00102    return fnew ;
00103 }
00104 
00105 /*---------------------------------------------------------------------
00106    Shift a 1D image (timeseries).  Values that are WAY_BIG act as
00107    blocks to the shift.
00108 -----------------------------------------------------------------------*/
00109 
00110 MRI_IMAGE * mri_shift_1D( MRI_IMAGE * im , float shift )
00111 {
00112    MRI_IMAGE * newim , * flim ;
00113    float * newar , * flar , * shar ;
00114    int ii , ibot,itop , nx ;
00115 
00116    /*-- sanity check --*/
00117 
00118    if( im == NULL ) return NULL ;
00119 
00120    /*-- create output image --*/
00121 
00122    if( im->kind != MRI_float ) flim = mri_to_float( im ) ;
00123    else                        flim = im ;
00124    flar = MRI_FLOAT_PTR(flim) ;
00125 
00126    nx    = flim->nx ;
00127    newim = mri_new( nx , 1 , MRI_float ) ;
00128    newar = MRI_FLOAT_PTR(newim) ;
00129 
00130    /*-- scan for unbroken blocks to shift --*/
00131 
00132    ibot = 0 ;
00133    while( ibot < nx ){
00134 
00135       if( flar[ibot] >= WAY_BIG ){    /* just copy values */
00136          newar[ibot] = flar[ibot] ;   /* that are WAY_BIG */
00137          ibot++ ;
00138          continue ;
00139       }
00140 
00141       for( ii=ibot+1 ; ii < nx ; ii++ )    /* scan for next WAY_BIG */
00142          if( flar[ii] >= WAY_BIG ) break ;
00143 
00144       itop = ii ;  /* values from ibot to itop-1 are OK to shift */
00145 
00146       /* shift and copy output into new image */
00147 
00148       shar = shifter( itop-ibot , flar+ibot , shift ) ;
00149       for( ii=ibot ; ii < itop ; ii++ ) newar[ii] = shar[ii-ibot] ;
00150       free(shar) ; shar = NULL ;
00151 
00152       ibot = itop ;  /* start here next loop */
00153    }
00154 
00155    /*-- cleanup and exit --*/
00156 
00157    if( flim != im ) mri_free(flim) ;
00158    return newim ;
00159 }