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 } |