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 #include "mrilib.h"
8 
9 /*----------------------------------------------------------------------
10    Input:  far   = array of floats length nar
11            shift = fractional shift
12 
13    Output: Newly malloc-ed array that is supposed to be like
14            far[ii-shift] for ii=0..nar-1.
15 
16    Notes: * The shift is to the RIGHT (a good Republican, of course).
17           * If nar==1, then the output is an array of length 1 that is
18              just a copy of the input.
19           * Otherwise, cubic interpolation is used.
20 ------------------------------------------------------------------------*/
21 
22 #define GET_AS_BIG(name,type,dim)                                           \
23    do{ if( (dim) > name ## _size ){                                         \
24           if( name != NULL ) free(name) ;                                   \
25           name = (type *) malloc( sizeof(type) * (dim) ) ;                  \
26           if( name == NULL ){                                               \
27              fprintf(stderr,"*** can't malloc shifter space\n"); EXIT(1); } \
28           name ## _size = (dim) ; } } while(0)
29 
30 /* cubic interpolation polynomials */
31 
32 #define P_M1(x)  ((x)*(1.0-(x))*((x)-2.0))
33 #define P_00(x)  (3.0*((x)+1.0)*((x)-1.0)*((x)-2.0))
34 #define P_P1(x)  (3.0*(x)*((x)+1.0)*(2.0-(x)))
35 #define P_P2(x)  ((x)*((x)+1.0)*((x)-1.0))
36 #define SIXTH    0.1666667
37 
shifter(int nar,float * far,float shift)38 float * shifter( int nar , float * far , float shift )
39 {
40    int ii,jj , nup , nmid , ix ;
41    float xx , wt_m1,wt_00,wt_p1,wt_p2 , fmin,fmax ;
42    float * fnew ;
43 
44    static int fl_size  = 0 ;     /* workspace (will hang around between calls) */
45    static float * fl = NULL ;
46 
47    /*-- sanity checks --*/
48 
49    if( nar <= 0 || far == NULL ) return NULL ;
50 
51    if( nar == 1 ){
52       fnew = (float *) malloc( sizeof(float) ) ;
53       if( fnew == NULL ){
54           fprintf(stderr,"*** can't malloc shifter output\n"); EXIT(1);
55       }
56      *fnew = far[0] ;
57       return fnew ;
58    }
59 
60    /*-- get workspace --*/
61 
62    nup = nar + (int)( 2.0*fabs(shift) + 6.0 ) ;
63    GET_AS_BIG(fl,float,nup) ;
64 
65    /*-- Insert data:
66           far[0] --> fl[nmid], etc.;
67           fl[] points before nmid are copies of far[0];
68           fl[] points after nmid+nar-1 are copiex of far[nar-1]. --*/
69 
70    nmid  = (nup-nar) / 2 ;
71    for( ii=0 ; ii < nup ; ii++ ){
72       jj = ii - nmid ;
73       if( jj < 0 ) jj = 0 ; else if( jj >= nar ) jj = nar-1 ;
74       fl[ii] = far[jj] ;
75    }
76 
77    /*-- put results into output array --*/
78 
79    fnew = (float *) malloc( sizeof(float) * nar ) ;
80    if( fnew == NULL ){
81        fprintf(stderr,"*** can't malloc shifter output\n"); EXIT(1);
82    }
83 
84    fmax = fmin = far[0] ;          /* find min and max of input */
85    for( ii=1 ; ii < nar ; ii++ ){  /* for "clipping" purposes   */
86       fmax = MAX(fmax,far[ii]) ;
87       fmin = MIN(fmin,far[ii]) ;
88    }
89 
90    for( ii=0 ; ii < nar ; ii++ ){
91       xx = ii+nmid - shift ;  /* "index" in fl we want */
92       ix = (int) xx ; xx = xx - ix ;
93       wt_m1 = P_M1(xx) ; wt_00 = P_00(xx) ;
94       wt_p1 = P_P1(xx) ; wt_p2 = P_P2(xx) ;
95       fnew[ii] = SIXTH * (  wt_m1 * fl[ix-1] + wt_00 * fl[ix]
96                           + wt_p1 * fl[ix+1] + wt_p2 * fl[ix+2] ) ;
97 
98            if( fnew[ii] < fmin ) fnew[ii] = fmin ;
99       else if( fnew[ii] > fmax ) fnew[ii] = fmax ;
100    }
101 
102    return fnew ;
103 }
104 
105 /*---------------------------------------------------------------------
106    Shift a 1D image (timeseries).  Values that are WAY_BIG act as
107    blocks to the shift.
108 -----------------------------------------------------------------------*/
109 
mri_shift_1D(MRI_IMAGE * im,float shift)110 MRI_IMAGE * mri_shift_1D( MRI_IMAGE * im , float shift )
111 {
112    MRI_IMAGE * newim , * flim ;
113    float * newar , * flar , * shar ;
114    int ii , ibot,itop , nx ;
115 
116    /*-- sanity check --*/
117 
118    if( im == NULL ) return NULL ;
119 
120    /*-- create output image --*/
121 
122    if( im->kind != MRI_float ) flim = mri_to_float( im ) ;
123    else                        flim = im ;
124    flar = MRI_FLOAT_PTR(flim) ;
125 
126    nx    = flim->nx ;
127    newim = mri_new( nx , 1 , MRI_float ) ;
128    newar = MRI_FLOAT_PTR(newim) ;
129 
130    /*-- scan for unbroken blocks to shift --*/
131 
132    ibot = 0 ;
133    while( ibot < nx ){
134 
135       if( flar[ibot] >= WAY_BIG ){    /* just copy values */
136          newar[ibot] = flar[ibot] ;   /* that are WAY_BIG */
137          ibot++ ;
138          continue ;
139       }
140 
141       for( ii=ibot+1 ; ii < nx ; ii++ )    /* scan for next WAY_BIG */
142          if( flar[ii] >= WAY_BIG ) break ;
143 
144       itop = ii ;  /* values from ibot to itop-1 are OK to shift */
145 
146       /* shift and copy output into new image */
147 
148       shar = shifter( itop-ibot , flar+ibot , shift ) ;
149       for( ii=ibot ; ii < itop ; ii++ ) newar[ii] = shar[ii-ibot] ;
150       free(shar) ; shar = NULL ;
151 
152       ibot = itop ;  /* start here next loop */
153    }
154 
155    /*-- cleanup and exit --*/
156 
157    if( flim != im ) mri_free(flim) ;
158    return newim ;
159 }
160