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