Doxygen Source Code Documentation
qmed.c File Reference
#include "mrilib.h"Go to the source code of this file.
| Defines | |
| #define | MED3(a, b, c) | 
| #define | SWAP(x, y) (temp=(x),(x)=(y),(y)=temp) | 
| Functions | |
| float | qmed_float (int n, float *ar) | 
| int | main (int argc, char *argv[]) | 
Define Documentation
| 
 | 
| Value: | 
| 
 | 
| 
 | 
Function Documentation
| 
 | ||||||||||||
| \** File : SUMA.c 
 Input paramters : 
 
 
 Definition at line 90 of file qmed.c. References argc, COX_cpu_time(), malloc, qmed_float(), and qsort_float(). 
 00091 {
00092    int nn,ii , med , mrep=1 , jj ;
00093    float * ar , * br ;
00094    float mm ;
00095    double ct , ctsum=0.0 ;
00096 
00097    if( argc < 2 ){printf("Usage: qmed N [m]\n");exit(0);}
00098    nn = strtol(argv[1],NULL,10) ; if( nn == 0 ) exit(1);
00099 
00100    if( argc == 3 ){
00101       mrep = strtol(argv[2],NULL,10) ; if( mrep < 1 ) mrep = 1 ;
00102    }
00103 
00104    med = ( nn > 0 ) ; if( !med ) nn = -nn ;
00105 
00106    ar = (float *) malloc( sizeof(float)*nn ) ;
00107    br = (float *) malloc( sizeof(float)*nn ) ;
00108    srand48((long)(237+jj)) ;
00109    for( ii=0 ; ii < nn ; ii++ ) br[ii] = (float)(lrand48() % nn)*0.1 ;
00110 
00111    ct = COX_cpu_time() ;
00112    for( jj=0 ; jj < mrep ; jj++ ){
00113 
00114       memcpy( ar , br , sizeof(float)*nn ) ;
00115 
00116       if( med ){
00117          mm = qmed_float( nn , ar ) ;
00118       } else {
00119          qsort_float( nn , ar ) ;
00120          if( nn%2 == 1 ) mm = ar[nn/2] ;
00121          else            mm = 0.5*(ar[nn/2]+ar[nn/2-1]) ;
00122       }
00123    }
00124    ct = COX_cpu_time() - ct ;
00125    printf("Median = %g  CPU = %g\n",mm,ct) ;
00126    exit(0) ;
00127 }
 | 
| 
 | ||||||||||||
| 
 Definition at line 21 of file qmed.c. References a, i, left, MED3, right, and SWAP. Referenced by extreme_proj(), find_unusual_correlations(), main(), median21_box_func(), median9_box_func(), mri_medianfilter(), qmedmad_float(), STATS_tsfunc(), THD_median_brick(), THD_outlier_count(), and unusuality(). 
 00022 {
00023    register int i , j ;           /* scanning indices */
00024    register float temp , pivot ;  /* holding places */
00025    register float * a = ar ;
00026 
00027    int left , right , mid , nodd ;
00028 
00029    switch( n ){
00030       case 1: return ar[0] ;
00031       case 2: return 0.5*(ar[0]+ar[1]) ;
00032       case 3: return MED3( ar[0] , ar[1] , ar[2] ) ;
00033    }
00034 
00035    left = 0 ; right = n-1 ; mid = n/2 ; nodd = ((n & 1) != 0) ;
00036 
00037    /* loop while the subarray is at least 3 long */
00038 
00039    while( right-left > 1  ){  /* work on subarray from left -> right */
00040 
00041       i = ( left + right ) / 2 ;   /* middle of subarray */
00042 
00043       /* sort the left, middle, and right a[]'s */
00044 
00045       if( a[left] > a[i]     ) SWAP( a[left]  , a[i]     ) ;
00046       if( a[left] > a[right] ) SWAP( a[left]  , a[right] ) ;
00047       if( a[i] > a[right]    ) SWAP( a[right] , a[i]     ) ;
00048 
00049       pivot = a[i] ;                 /* a[i] is the median-of-3 pivot! */
00050       a[i]  = a[right] ;
00051 
00052       i = left ;                     /* initialize scanning */
00053       j = right ;
00054 
00055       /*----- partition:  move elements bigger than pivot up and elements
00056                           smaller than pivot down, scanning in from ends -----*/
00057 
00058       do{
00059         for( ; a[++i] < pivot ; ) ;  /* scan i up,   until a[i] >= pivot */
00060         for( ; a[--j] > pivot ; ) ;  /* scan j down, until a[j] <= pivot */
00061 
00062         if( j <= i ) break ;         /* if j meets i, quit */
00063 
00064         SWAP( a[i] , a[j] ) ;
00065       } while( 1 ) ;
00066 
00067       /*----- at this point, the array is partitioned -----*/
00068 
00069       a[right] = a[i] ;           /* restore the pivot */
00070       a[i]     = pivot ;
00071 
00072       if( i == mid ){             /* good luck */
00073          if( nodd ) return pivot ;   /* exact middle of array */
00074 
00075          temp = a[left] ;            /* must find next largest element below */
00076          for( j=left+1 ; j < i ; j++ )
00077             if( a[j] > temp ) temp = a[j] ;
00078 
00079          return 0.5*(pivot+temp) ;   /* return average */
00080       }
00081 
00082       if( i <  mid ) left  = i ; /* throw away bottom partition */
00083       else           right = i ; /* throw away top partition    */
00084 
00085    }  /* end of while sub-array is long */
00086 
00087    return (nodd) ? a[mid] : 0.5*(a[mid]+a[mid-1]) ;
00088 }
 | 
 
                             
                             
                             
                             
                             
                             
                             
                             
                             
                             
                             
                             
 
 
 
 
       
	   
	   
	   
	  