Skip to content

AFNI/NIfTI Server

Sections
Personal tools
You are here: Home » AFNI » Documentation

Doxygen Source Code Documentation


Main Page   Alphabetical List   Data Structures   File List   Data Fields   Globals   Search  

cs_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)
#define SORT2(a, b)   if(a>b) SWAP(a,b);

Functions

float median_float4 (float, float, float, float)
float median_float5 (float *)
float median_float7 (float *)
float median_float9 (float *)
float qmed_float (int n, float *ar)
void qmedmad_float (int n, float *ar, float *med, float *mad)

Define Documentation

#define MED3 a,
b,
c   
 

Value:

( (a < b) ? ( (a > c) ? a : (b > c) ? c : b )  \
                              : ( (b > c) ? b : (a > c) ? c : a ) )

Definition at line 24 of file cs_qmed.c.

Referenced by qmed_float().

#define SORT2 a,
b       if(a>b) SWAP(a,b);
 

Definition at line 172 of file cs_qmed.c.

Referenced by median_float5(), median_float7(), and median_float9().

#define SWAP x,
y       (temp=x,x=y,y=temp)
 

Definition at line 30 of file cs_qmed.c.

Referenced by qmed_float().


Function Documentation

float median_float4 float   ,
float   ,
float   ,
float   
[static]
 

Definition at line 130 of file cs_qmed.c.

References a, and c.

Referenced by qmed_float().

00131 {
00132   register float t1,t2,t3;
00133 
00134   if (a > b){t1 = a; a = b;} else t1 = b;  /* t1 = max(a,b) */
00135   if (c > d){t2 = d; d = c;} else t2 = c;  /* t2 = min(c,d) */
00136 
00137   /* now have: a <= t1 and t2 <= d */
00138 
00139   t3 = (a > t2) ? a : t2 ;
00140   t2 = (d < t1) ? d : t1 ;
00141 
00142   return 0.5*(t2+t3) ;
00143 }

float median_float5 float *    [static]
 

Definition at line 174 of file cs_qmed.c.

References p, and SORT2.

Referenced by qmed_float().

00175 {
00176     register float temp ;
00177     SORT2(p[0],p[1]) ; SORT2(p[3],p[4]) ; SORT2(p[0],p[3]) ;
00178     SORT2(p[1],p[4]) ; SORT2(p[1],p[2]) ; SORT2(p[2],p[3]) ;
00179     SORT2(p[1],p[2]) ; return(p[2]) ;
00180 }

float median_float7 float *    [static]
 

Definition at line 182 of file cs_qmed.c.

References p, and SORT2.

Referenced by qmed_float().

00183 {
00184     register float temp ;
00185     SORT2(p[0],p[1]) ; SORT2(p[4],p[5]) ; SORT2(p[1],p[2]) ;
00186     SORT2(p[5],p[6]) ; SORT2(p[0],p[1]) ; SORT2(p[4],p[5]) ;
00187     SORT2(p[0],p[4]) ; SORT2(p[2],p[6]) ; SORT2(p[1],p[3]) ;
00188     SORT2(p[3],p[5]) ; SORT2(p[1],p[3]) ; SORT2(p[2],p[3]) ;
00189     SORT2(p[3],p[4]) ; SORT2(p[2],p[3]) ; return(p[3]) ;
00190 }

float median_float9 float *    [static]
 

Definition at line 192 of file cs_qmed.c.

References p, and SORT2.

Referenced by qmed_float().

00193 {
00194     register float temp ;
00195     SORT2(p[1],p[2]) ; SORT2(p[4],p[5]) ; SORT2(p[7],p[8]) ;
00196     SORT2(p[0],p[1]) ; SORT2(p[3],p[4]) ; SORT2(p[6],p[7]) ;
00197     SORT2(p[1],p[2]) ; SORT2(p[4],p[5]) ; SORT2(p[7],p[8]) ;
00198     SORT2(p[0],p[3]) ; SORT2(p[5],p[8]) ; SORT2(p[4],p[7]) ;
00199     SORT2(p[3],p[6]) ; SORT2(p[1],p[4]) ; SORT2(p[2],p[5]) ;
00200     SORT2(p[4],p[7]) ; SORT2(p[4],p[2]) ; SORT2(p[6],p[4]) ;
00201     SORT2(p[4],p[2]) ; return(p[4]) ;
00202 }

float qmed_float int    n,
float *    ar
 

Definition at line 32 of file cs_qmed.c.

References a, i, left, MED3, median_float4(), median_float5(), median_float7(), median_float9(), right, and SWAP.

00033 {
00034    register int i , j ;           /* scanning indices */
00035    register float temp , pivot ;  /* holding places */
00036    register float * a = ar ;
00037 
00038    int left , right , mid , nodd ;
00039 
00040    switch( n ){
00041       case 1: return ar[0] ;
00042       case 2: return 0.5*(ar[0]+ar[1]) ;
00043       case 3: return MED3( ar[0] , ar[1] , ar[2] ) ;
00044       case 4: return median_float4( ar[0],ar[1],ar[2],ar[3] ) ;
00045       case 5: return median_float5( ar ) ;
00046       case 7: return median_float7( ar ) ;
00047       case 9: return median_float9( ar ) ;
00048    }
00049 
00050    left = 0 ; right = n-1 ; mid = n/2 ; nodd = ((n & 1) != 0) ;
00051 
00052    /* loop while the subarray is at least 3 long */
00053 
00054    while( right-left > 1  ){  /* work on subarray from left -> right */
00055 
00056       i = ( left + right ) / 2 ;   /* middle of subarray */
00057 
00058       /* sort the left, middle, and right a[]'s */
00059 
00060       if( a[left] > a[i]     ) SWAP( a[left]  , a[i]     ) ;
00061       if( a[left] > a[right] ) SWAP( a[left]  , a[right] ) ;
00062       if( a[i] > a[right]    ) SWAP( a[right] , a[i]     ) ;
00063 
00064       pivot = a[i] ;                 /* a[i] is the median-of-3 pivot! */
00065       a[i]  = a[right] ;
00066 
00067       i = left ;                     /* initialize scanning */
00068       j = right ;
00069 
00070       /*----- partition:  move elements bigger than pivot up and elements
00071                           smaller than pivot down, scanning in from ends -----*/
00072 
00073       do{
00074         for( ; a[++i] < pivot ; ) ;  /* scan i up,   until a[i] >= pivot */
00075         for( ; a[--j] > pivot ; ) ;  /* scan j down, until a[j] <= pivot */
00076 
00077         if( j <= i ) break ;         /* if j meets i, quit */
00078 
00079         SWAP( a[i] , a[j] ) ;
00080       } while( 1 ) ;
00081 
00082       /*----- at this point, the array is partitioned -----*/
00083 
00084       a[right] = a[i] ;           /* restore the pivot */
00085       a[i]     = pivot ;
00086 
00087       if( i == mid ){             /* good luck */
00088          if( nodd ) return pivot ;   /* exact middle of array */
00089 
00090          temp = a[left] ;            /* must find next largest element below */
00091          for( j=left+1 ; j < i ; j++ )
00092             if( a[j] > temp ) temp = a[j] ;
00093 
00094          return 0.5*(pivot+temp) ;   /* return average */
00095       }
00096 
00097       if( i <  mid ) left  = i ; /* throw away bottom partition */
00098       else           right = i ; /* throw away top partition    */
00099 
00100    }  /* end of while sub-array is long */
00101 
00102    return (nodd) ? a[mid] : 0.5*(a[mid]+a[mid-1]) ;
00103 }

void qmedmad_float int    n,
float *    ar,
float *    med,
float *    mad
 

Definition at line 109 of file cs_qmed.c.

References free, malloc, q, and qmed_float().

Referenced by main(), plot_graphs(), and THD_outlier_count().

00110 {
00111    float * q = (float *) malloc(sizeof(float)*n) ;
00112    float me,ma ;
00113    register int ii ;
00114 
00115    memcpy(q,ar,sizeof(float)*n) ;  /* duplicate input array */
00116    me = qmed_float( n , q ) ;      /* compute median */
00117 
00118    for( ii=0 ; ii < n ; ii++ )     /* subtract off median */
00119       q[ii] = fabs(q[ii]-me) ;     /* (absolute deviation) */
00120 
00121    ma = qmed_float( n , q ) ;      /* MAD = median absolute deviation */
00122 
00123    free(q) ;                       /* 05 Nov 2001 */
00124 
00125    *med = me ; *mad = ma ; return ;
00126 }
 

Powered by Plone

This site conforms to the following standards: