00001 #include "afni.h"
00002 
00003 /*--------------- Sample 0D function: log10 of each input point ------------*/
00004 
00005 void log10_func( int num , float *vec )
00006 {
00007    int ii ;
00008    float vmax , vmin ;
00009 
00010    if( num <= 0 || vec == NULL ) return ;
00011 
00012    /* find largest element */
00013 
00014    vmax = vec[0] ;
00015    for( ii=1 ; ii < num ; ii++ ) vmax = MAX( vmax , vec[ii] ) ;
00016 
00017    /* if all are nonpositive, return all zeros */
00018 
00019    if( vmax <= 0.0 ){
00020       for( ii=0 ; ii < num ; ii++ ) vec[ii] = 0.0 ;
00021       return ;
00022    }
00023 
00024    /* find smallest positive element */
00025 
00026    vmin = vmax ;
00027    for( ii=0 ; ii < num ; ii++ )
00028       if( vec[ii] > 0.0 ) vmin = MIN( vmin , vec[ii] ) ;
00029 
00030    /* take log10 of each positive element;
00031       nonpositive elements get the log10 of vmin instead */
00032 
00033    vmin = log10(vmin) ;
00034    for( ii=0 ; ii < num ; ii++ )
00035       vec[ii] = (vec[ii] > 0.0) ? log10(vec[ii]) : vmin ;
00036 
00037    return ;
00038 }
00039 
00040 /*------------ Sample 0D function: Signed sqrt of each input point ---------*/
00041 
00042 void ssqrt_func( int num , float *vec )
00043 {
00044    int ii ;
00045    double val ;
00046 
00047    if( num <= 0 || vec == NULL ) return ;
00048 
00049    for( ii=0 ; ii < num ; ii++ ){
00050       val = sqrt(fabs(vec[ii])) ;                /* will be positive */
00051       vec[ii] = (vec[ii] >= 0.0) ? val : -val ;  /* output sign = input sign */
00052    }
00053 
00054    return ;
00055 }
00056 
00057 /*--------------- Sample 1D function: Order Statistics Filter -------------*/
00058 
00059 void osfilt3_func( int num , double to,double dt, float *vec )
00060 {
00061    int ii ;
00062    float aa,bb,cc ;
00063 
00064    bb = vec[0] ; cc = vec[1] ;
00065    for( ii=1 ; ii < num-1 ; ii++ ){
00066       aa = bb ; bb = cc ; cc = vec[ii+1] ;
00067       vec[ii] = OSFILT(aa,bb,cc) ;         /* see mrilib.h */
00068    }
00069 
00070    return ;
00071 }
00072 
00073 /*--------------- Sample 1D function: Median of 3 Filter ----------------*/
00074 
00075 void median3_func( int num , double to,double dt, float *vec )
00076 {
00077    int ii ;
00078    float aa,bb,cc ;
00079 
00080    bb = vec[0] ; cc = vec[1] ;
00081    for( ii=1 ; ii < num-1 ; ii++ ){
00082       aa = bb ; bb = cc ; cc = vec[ii+1] ;
00083       vec[ii] = MEDIAN(aa,bb,cc) ;         /* see mrilib.h */
00084    }
00085 
00086    return ;
00087 }
00088 
00089 /*---------------- Sample 1D function: abs(FFT) [30 Jun 2000] --------------*/
00090 
00091 void absfft_func( int num , double to,double dt, float *vec )
00092 {
00093    static complex *cx=NULL ;
00094    static int      ncx=0 , numold=0 ;
00095    float f0,f1 ;
00096    int ii ;
00097 
00098    if( num < 2 ) return ;
00099    if( num > numold ){
00100       numold = num ;
00101       ncx    = csfft_nextup(numold) ;
00102       if( cx != NULL ) free(cx) ;
00103       cx = (complex *) malloc(sizeof(complex)*ncx) ;
00104    }
00105 
00106    get_linear_trend( num , vec , &f0,&f1 ) ;  /* thd_detrend.c */
00107 
00108    for( ii=0 ; ii < num ; ii++ ){ cx[ii].r = vec[ii]-(f0+f1*ii); cx[ii].i = 0.0; }
00109    for(      ; ii < ncx ; ii++ ){ cx[ii].r = cx[ii].i = 0.0 ; }
00110 
00111    csfft_cox( -1 , ncx , cx ) ;               /* csfft.c */
00112 
00113    vec[0] = 0.0 ;
00114    for( ii=1 ; ii < num ; ii++ ) vec[ii] = CABS(cx[ii]) ;
00115 
00116    return ;
00117 }
00118 
00119 /*----------- Sample slice projection functions [31 Jan 2002] -----------*/
00120 
00121 float min_proj( int n , float *ar )
00122 {
00123    float v = ar[0] ;
00124    int ii ;
00125    for( ii=1 ; ii < n ; ii++ ) if( v > ar[ii] ) v = ar[ii] ;
00126    return v ;
00127 }
00128 
00129 float max_proj( int n , float *ar )
00130 {
00131    float v = ar[0] ;
00132    int ii ;
00133    for( ii=1 ; ii < n ; ii++ ) if( v < ar[ii] ) v = ar[ii] ;
00134    return v ;
00135 }
00136 
00137 float mean_proj( int n , float *ar )
00138 {
00139    float v=0.0 ;
00140    int ii ;
00141    for( ii=0 ; ii < n ; ii++ ) v += ar[ii] ;
00142    return (v/n) ;
00143 }
00144 
00145 float extreme_proj( int n , float *ar )  /* 02 Feb 2002 */
00146 {
00147    float vv,ww , med=qmed_float(n,ar) ;
00148    int ii , jj ;
00149 
00150    jj = 0 ; vv = fabs(ar[0]-med) ;       /* Find the value */
00151    for( ii=1 ; ii < n ; ii++ ){          /* furthest from */
00152       ww = fabs(ar[ii]-med) ;            /* the median.  */
00153       if( ww > vv ){ vv=ww; jj=ii; }
00154    }
00155    return ar[jj] ;
00156 }
00157 
00158 /*======================================================================*/
00159 /*----------------------- Sample 2D transformations --------------------*/
00160 /*======================================================================*/
00161 
00162 static float *atemp = NULL ;
00163 static int   natemp = -666 ;
00164 
00165 #define MAKE_ATEMP(nvox)                     \
00166   do{ if( natemp < (nvox) ){                 \
00167          if( atemp != NULL ) free(atemp) ;   \
00168          natemp = (nvox) ;                   \
00169          atemp  = (float *) malloc( sizeof(float) * natemp ) ; } } while(0)
00170 
00171 #define AT(i,j) atemp[(i)+(j)*nx]
00172 
00173 /*------------------------------------------------------------------------*/
00174 
00175 void median9_box_func( int nx , int ny , double dx, double dy, float *ar )
00176 {
00177    int ii , jj , nxy , joff ;
00178    float aa[9] ;
00179    float *ajj , *ajm , *ajp ;
00180 
00181    if( nx < 3 || ny < 3 ) return ;
00182 
00183    /** make space and copy input into it **/
00184 
00185    nxy = nx * ny ;
00186    MAKE_ATEMP(nxy) ; if( atemp == NULL ) return ;
00187    memcpy(atemp,ar,sizeof(float)*nxy) ;
00188 
00189    /** process copy of input back into the input array **/
00190 
00191    for( jj=0 ; jj < ny ; jj++ ){
00192 
00193       joff = jj * nx ;      /* offset into this row */
00194       ajj  = atemp + joff ; /* pointer to this row */
00195 
00196       ajm  = (jj==0   ) ? ajj : ajj-nx ;  /* pointer to last row */
00197       ajp  = (jj==ny-1) ? ajj : ajj+nx ;  /* pointer to next row */
00198 
00199       /* do interior points of this row */
00200 
00201       for( ii=1 ; ii < nx-1 ; ii++ ){
00202          aa[0] = ajm[ii-1] ; aa[1] = ajm[ii] ; aa[2] = ajm[ii+1] ;
00203          aa[3] = ajj[ii-1] ; aa[4] = ajj[ii] ; aa[5] = ajj[ii+1] ;
00204          aa[6] = ajp[ii-1] ; aa[7] = ajp[ii] ; aa[8] = ajp[ii+1] ;
00205 #if 0
00206          isort_float( 9 , aa ) ; ar[ii+joff] = aa[4] ;
00207 #else
00208          ar[ii+joff] = qmed_float(9,aa) ;  /* 25 Oct 2000 */
00209 #endif
00210       }
00211 
00212       /* do leading edge point (ii=0) */
00213 
00214       aa[0] = ajm[0] ; aa[1] = ajm[0] ; aa[2] = ajm[1] ;
00215       aa[3] = ajj[0] ; aa[4] = ajj[0] ; aa[5] = ajj[1] ;
00216       aa[6] = ajp[0] ; aa[7] = ajp[0] ; aa[8] = ajp[1] ;
00217 #if 0
00218       isort_float( 9 , aa ) ; ar[joff] = aa[4] ;
00219 #else
00220       ar[joff] = qmed_float(9,aa) ;  /* 25 Oct 2000 */
00221 #endif
00222 
00223       /* do trailing edge point (ii=nx-1) */
00224 
00225       aa[0] = ajm[nx-2] ; aa[1] = ajm[nx-1] ; aa[2] = ajm[nx-1] ;
00226       aa[3] = ajj[nx-2] ; aa[4] = ajj[nx-1] ; aa[5] = ajj[nx-1] ;
00227       aa[6] = ajp[nx-2] ; aa[7] = ajp[nx-1] ; aa[8] = ajp[nx-1] ;
00228 #if 0
00229       isort_float( 9 , aa ) ; ar[nx-1+joff] = aa[4] ;
00230 #else
00231       ar[nx-1+joff] = qmed_float(9,aa) ;  /* 25 Oct 2000 */
00232 #endif
00233    }
00234    return ;
00235 }
00236 
00237 /*------------------------------------------------------------------------*/
00238 
00239 void winsor9_box_func( int nx , int ny , double dx, double dy, float *ar )
00240 {
00241    int ii , jj , nxy , joff ;
00242    float aa[9] ;
00243    float *ajj , *ajm , *ajp ;
00244 
00245    if( nx < 3 || ny < 3 ) return ;
00246 
00247    /** make space and copy input into it **/
00248 
00249    nxy = nx * ny ;
00250    MAKE_ATEMP(nxy) ; if( atemp == NULL ) return ;
00251    memcpy(atemp,ar,sizeof(float)*nxy) ;
00252 
00253    /** process copy of input back into the input array **/
00254 
00255    for( jj=0 ; jj < ny ; jj++ ){
00256 
00257       joff = jj * nx ;      /* offset into this row */
00258       ajj  = atemp + joff ; /* pointer to this row */
00259 
00260       ajm  = (jj==0   ) ? ajj : ajj-nx ;  /* pointer to last row */
00261       ajp  = (jj==ny-1) ? ajj : ajj+nx ;  /* pointer to next row */
00262 
00263       /* do interior points of this row */
00264 
00265       for( ii=1 ; ii < nx-1 ; ii++ ){
00266          aa[0] = ajm[ii-1] ; aa[1] = ajm[ii] ; aa[2] = ajm[ii+1] ;
00267          aa[3] = ajj[ii-1] ; aa[4] = ajj[ii] ; aa[5] = ajj[ii+1] ;
00268          aa[6] = ajp[ii-1] ; aa[7] = ajp[ii] ; aa[8] = ajp[ii+1] ;
00269          isort_float( 9 , aa ) ;
00270               if( ar[ii+joff] < aa[2] ) ar[ii+joff] = aa[2] ;
00271          else if( ar[ii+joff] > aa[6] ) ar[ii+joff] = aa[6] ;
00272       }
00273 
00274       /* do leading edge point (ii=0) */
00275 
00276       aa[0] = ajm[0] ; aa[1] = ajm[0] ; aa[2] = ajm[1] ;
00277       aa[3] = ajj[0] ; aa[4] = ajj[0] ; aa[5] = ajj[1] ;
00278       aa[6] = ajp[0] ; aa[7] = ajp[0] ; aa[8] = ajp[1] ;
00279       isort_float( 9 , aa ) ;
00280            if( ar[joff] < aa[2] ) ar[joff] = aa[2] ;
00281       else if( ar[joff] > aa[6] ) ar[joff] = aa[6] ;
00282 
00283       /* do trailing edge point (ii=nx-1) */
00284 
00285       aa[0] = ajm[nx-2] ; aa[1] = ajm[nx-1] ; aa[2] = ajm[nx-1] ;
00286       aa[3] = ajj[nx-2] ; aa[4] = ajj[nx-1] ; aa[5] = ajj[nx-1] ;
00287       aa[6] = ajp[nx-2] ; aa[7] = ajp[nx-1] ; aa[8] = ajp[nx-1] ;
00288       isort_float( 9 , aa ) ;
00289            if( ar[nx-1+joff] < aa[2] ) ar[nx-1+joff] = aa[2] ;
00290       else if( ar[nx-1+joff] > aa[6] ) ar[nx-1+joff] = aa[6] ;
00291    }
00292    return ;
00293 }
00294 
00295 /*------------------------------------------------------------------------*/
00296 
00297 void osfilt9_box_func( int nx , int ny , double dx, double dy, float *ar )
00298 {
00299    int ii , jj , nxy , joff ;
00300    float aa[9] ;
00301    float *ajj , *ajm , *ajp ;
00302 
00303    if( nx < 3 || ny < 3 ) return ;
00304 
00305    /** make space and copy input into it **/
00306 
00307    nxy = nx * ny ;
00308    MAKE_ATEMP(nxy) ; if( atemp == NULL ) return ;
00309    memcpy(atemp,ar,sizeof(float)*nxy) ;
00310 
00311    /** process copy of input back into the input array **/
00312 
00313    for( jj=0 ; jj < ny ; jj++ ){
00314 
00315       joff = jj * nx ;      /* offset into this row */
00316       ajj  = atemp + joff ; /* pointer to this row */
00317 
00318       ajm  = (jj==0   ) ? ajj : ajj-nx ;  /* pointer to last row */
00319       ajp  = (jj==ny-1) ? ajj : ajj+nx ;  /* pointer to next row */
00320 
00321       /* do interior points of this row */
00322 
00323 #undef  OSUM
00324 #define OSUM(a,b,c,d,e) ( 0.1*((a)+(e)) + 0.2*((b)+(d)) + 0.4*(c) )
00325 
00326       for( ii=1 ; ii < nx-1 ; ii++ ){
00327          aa[0] = ajm[ii-1] ; aa[1] = ajm[ii] ; aa[2] = ajm[ii+1] ;
00328          aa[3] = ajj[ii-1] ; aa[4] = ajj[ii] ; aa[5] = ajj[ii+1] ;
00329          aa[6] = ajp[ii-1] ; aa[7] = ajp[ii] ; aa[8] = ajp[ii+1] ;
00330          isort_float( 9 , aa ) ;
00331          ar[ii+joff] = OSUM( aa[2],aa[3],aa[4],aa[5],aa[6] ) ;
00332       }
00333 
00334       /* do leading edge point (ii=0) */
00335 
00336       aa[0] = ajm[0] ; aa[1] = ajm[0] ; aa[2] = ajm[1] ;
00337       aa[3] = ajj[0] ; aa[4] = ajj[0] ; aa[5] = ajj[1] ;
00338       aa[6] = ajp[0] ; aa[7] = ajp[0] ; aa[8] = ajp[1] ;
00339       isort_float( 9 , aa ) ;
00340       ar[joff] = OSUM( aa[2],aa[3],aa[4],aa[5],aa[6] ) ;
00341 
00342       /* do trailing edge point (ii=nx-1) */
00343 
00344       aa[0] = ajm[nx-2] ; aa[1] = ajm[nx-1] ; aa[2] = ajm[nx-1] ;
00345       aa[3] = ajj[nx-2] ; aa[4] = ajj[nx-1] ; aa[5] = ajj[nx-1] ;
00346       aa[6] = ajp[nx-2] ; aa[7] = ajp[nx-1] ; aa[8] = ajp[nx-1] ;
00347       isort_float( 9 , aa ) ;
00348       ar[nx-1+joff] = OSUM( aa[2],aa[3],aa[4],aa[5],aa[6] ) ;
00349    }
00350    return ;
00351 }
00352 
00353 /*------------------------------------------------------------------------*/
00354 
00355 void median21_box_func( int nx , int ny , double dx, double dy, float *ar )
00356 {
00357    int ii , jj , nxy , joff ;
00358    float aa[21] ;
00359    float *ajj , *ajm , *ajp , *ajmm , *ajpp ;
00360 
00361    if( nx < 5 || ny < 5 ) return ;
00362 
00363    /** make space and copy input into it **/
00364 
00365    nxy = nx * ny ;
00366    MAKE_ATEMP(nxy) ; if( atemp == NULL ) return ;
00367    memcpy(atemp,ar,sizeof(float)*nxy) ;
00368 
00369    /** process copy of input back into the input array **/
00370 
00371    for( jj=1 ; jj < ny-1 ; jj++ ){
00372 
00373       joff = jj * nx ;      /* offset into this row */
00374       ajj  = atemp + joff ; /* pointer to this row */
00375 
00376       ajm  = ajj-nx ;  /* pointer to last row */
00377       ajp  = ajj+nx ;  /* pointer to next row */
00378 
00379       ajmm = (jj == 1  ) ? ajm : ajm-nx ;  /* to last last row */
00380       ajpp = (jj ==ny-2) ? ajp : ajp+nx ;  /* to next next row */
00381 
00382       /* do interior points of this row */
00383 
00384       for( ii=2 ; ii < nx-2 ; ii++ ){
00385          aa[0]=ajmm[ii-1]; aa[1]=ajmm[ii]; aa[2]=ajmm[ii+1];
00386 
00387          aa[ 3]=ajm[ii-2]; aa[ 4]=ajm[ii-1]; aa[ 5]=ajm[ii]; aa[ 6]=ajm[ii+1]; aa[ 7]=ajm[ii+2];
00388          aa[ 8]=ajj[ii-2]; aa[ 9]=ajj[ii-1]; aa[10]=ajj[ii]; aa[11]=ajj[ii+1]; aa[12]=ajj[ii+2];
00389          aa[13]=ajp[ii-2]; aa[14]=ajp[ii-1]; aa[15]=ajp[ii]; aa[16]=ajp[ii+1]; aa[17]=ajp[ii+2];
00390 
00391          aa[18]=ajpp[ii-1]; aa[19]=ajpp[ii]; aa[20]=ajpp[ii+1];
00392 
00393 #if 0
00394          isort_float( 21 , aa ) ; ar[ii+joff] = aa[10] ;
00395 #else
00396          ar[ii+joff] = qmed_float(21,aa) ; /* 25 Oct 2000 */
00397 #endif
00398       }
00399 
00400    }
00401    return ;
00402 }
00403 
00404 /*------------------------------------------------------------------------*/
00405 
00406 void winsor21_box_func( int nx , int ny , double dx, double dy, float *ar )
00407 {
00408    int ii , jj , nxy , joff ;
00409    float aa[21] ;
00410    float *ajj , *ajm , *ajp , *ajmm , *ajpp ;
00411 
00412    static int kbot=-1 , ktop ;
00413 
00414    if( nx < 5 || ny < 5 ) return ;
00415 
00416    /** initialize cutoffs [07 Dec 1999] **/
00417 
00418    if( kbot < 0 ){
00419       char *ee = my_getenv("AFNI_WINSOR21_CUTOFF") ;
00420       kbot = 6 ;   /* default */
00421       if( ee != NULL ){
00422          ii = strtol( ee , NULL , 10 ) ;
00423          if( ii > 0 && ii < 10 ) kbot = ii ;
00424       }
00425       ktop = 20 - kbot ;
00426    }
00427 
00428    /** make space and copy input into it **/
00429 
00430    nxy = nx * ny ;
00431    MAKE_ATEMP(nxy) ; if( atemp == NULL ) return ;
00432    memcpy(atemp,ar,sizeof(float)*nxy) ;
00433 
00434    /** process copy of input back into the input array **/
00435 
00436    for( jj=1 ; jj < ny-1 ; jj++ ){
00437 
00438       joff = jj * nx ;      /* offset into this row */
00439       ajj  = atemp + joff ; /* pointer to this row */
00440 
00441       ajm  = ajj-nx ;  /* pointer to last row */
00442       ajp  = ajj+nx ;  /* pointer to next row */
00443 
00444       ajmm = (jj == 1  ) ? ajm : ajm-nx ;  /* to last last row */
00445       ajpp = (jj ==ny-2) ? ajp : ajp+nx ;  /* to next next row */
00446 
00447       /* do interior points of this row */
00448 
00449       for( ii=2 ; ii < nx-2 ; ii++ ){
00450          aa[0]=ajmm[ii-1]; aa[1]=ajmm[ii]; aa[2]=ajmm[ii+1];
00451 
00452          aa[ 3]=ajm[ii-2]; aa[ 4]=ajm[ii-1]; aa[ 5]=ajm[ii]; aa[ 6]=ajm[ii+1]; aa[ 7]=ajm[ii+2];
00453          aa[ 8]=ajj[ii-2]; aa[ 9]=ajj[ii-1]; aa[10]=ajj[ii]; aa[11]=ajj[ii+1]; aa[12]=ajj[ii+2];
00454          aa[13]=ajp[ii-2]; aa[14]=ajp[ii-1]; aa[15]=ajp[ii]; aa[16]=ajp[ii+1]; aa[17]=ajp[ii+2];
00455 
00456          aa[18]=ajpp[ii-1]; aa[19]=ajpp[ii]; aa[20]=ajpp[ii+1];
00457 
00458          isort_float( 21 , aa ) ;
00459 
00460               if( ar[ii+joff] < aa[kbot] ) ar[ii+joff] = aa[kbot] ;
00461          else if( ar[ii+joff] > aa[ktop] ) ar[ii+joff] = aa[ktop] ;
00462       }
00463 
00464    }
00465    return ;
00466 }
00467 
00468 /*------- [30 Jun 2000: abs(2D FFT) function] ----------------------------*/
00469 
00470 void fft2D_func( int nx , int ny , double dx, double dy, float *ar )
00471 {
00472    complex *cxar , *cpt ;
00473    int nxup,nyup , ii,jj ;
00474    float fi,fj , *fpt ;
00475 
00476    if( nx < 5 || ny < 5 ) return ;
00477 
00478    nxup = csfft_nextup_one35(nx) ;  /* get FFT size */
00479    nyup = csfft_nextup_one35(ny) ;
00480 
00481    cxar = (complex *) malloc(sizeof(complex)*nxup*nyup) ;
00482 
00483    /* copy input to output, sign-alternating and zero-padding along the way */
00484 
00485    cpt = cxar ;
00486    fpt = ar   ;
00487    fj  = 1.0  ;
00488    for( jj=0 ; jj < ny ; jj++ ){
00489       fi = fj ; fj = -fj ;
00490       for(ii=0; ii<nx  ; ii++){cpt->r=*fpt*fi; cpt->i=0.0; cpt++;fpt++;fi=-fi;}
00491       for(    ; ii<nxup; ii++){cpt->r=cpt->i=0.0; cpt++;}
00492    }
00493    for( ; jj < nyup ; jj++ ){cpt->r=cpt->i=0.0; cpt++;}
00494 
00495    /* row FFTs */
00496 
00497    for( jj=0 ; jj < ny ; jj++ )
00498       csfft_cox( -1 , nxup , cxar+jj*nxup ) ;
00499 
00500    /* column FFTs */
00501 
00502    cpt = (complex *) malloc(sizeof(complex)*nyup) ;
00503 
00504    for( ii=0 ; ii < nxup ; ii++ ){
00505       for( jj=0 ; jj < nyup ; jj++ ) cpt[jj] = cxar[ii+jj*nxup] ;
00506       csfft_cox( -1 , nyup , cpt ) ;
00507       for( jj=0 ; jj < nyup ; jj++ ) cxar[ii+jj*nxup] = cpt[jj] ;
00508    }
00509 
00510    /* copy to output */
00511 
00512    for( jj=0 ; jj < ny ; jj++ )
00513       for( ii=0 ; ii < nx ; ii++ )
00514          ar[ii+jj*nx] = CABS(cxar[ii+jj*nxup]) ;
00515 
00516    free(cxar) ; free(cpt) ; return ;
00517 }