Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
mri_shifter.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
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
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 ;
00045 static float * fl = NULL ;
00046
00047
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
00061
00062 nup = nar + (int)( 2.0*fabs(shift) + 6.0 ) ;
00063 GET_AS_BIG(fl,float,nup) ;
00064
00065
00066
00067
00068
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
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] ;
00085 for( ii=1 ; ii < nar ; ii++ ){
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 ;
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
00107
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
00117
00118 if( im == NULL ) return NULL ;
00119
00120
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
00131
00132 ibot = 0 ;
00133 while( ibot < nx ){
00134
00135 if( flar[ibot] >= WAY_BIG ){
00136 newar[ibot] = flar[ibot] ;
00137 ibot++ ;
00138 continue ;
00139 }
00140
00141 for( ii=ibot+1 ; ii < nx ; ii++ )
00142 if( flar[ii] >= WAY_BIG ) break ;
00143
00144 itop = ii ;
00145
00146
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 ;
00153 }
00154
00155
00156
00157 if( flim != im ) mri_free(flim) ;
00158 return newim ;
00159 }