00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #include "mrilib.h"
00008 
00009 /*------------------------------------------------------------------------
00010    Compute the median of an array of floats.  Will rearrange (partially
00011    sort) the array in the process.  The algorithm is based on Quicksort,
00012    where we only keep the partition that has the middle element.
00013    For large n, this is faster than finishing the whole sort.
00014 --------------------------------------------------------------------------*/
00015 
00016 static float median_float4(float,float,float,float) ;
00017 static float median_float5(float *) ;  /* 5,7,9 added 20 Oct 2000 */
00018 static float median_float7(float *) ;
00019 static float median_float9(float *) ;
00020 
00021 /* macro for median-of-3 */
00022 
00023 #undef  MED3
00024 #define MED3(a,b,c) ( (a < b) ? ( (a > c) ? a : (b > c) ? c : b )  \
00025                               : ( (b > c) ? b : (a > c) ? c : a ) )
00026 
00027 /* macro for swap (duh) */
00028 
00029 #undef  SWAP
00030 #define SWAP(x,y) (temp=x,x=y,y=temp)
00031 
00032 float qmed_float( int n , float * ar )
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 }
00104 
00105 /*---------------------------------------------------------------
00106   Return median and MAD, nondestructively -- 08 Mar 2001 - RWCox
00107 -----------------------------------------------------------------*/
00108 
00109 void qmedmad_float( int n, float *ar, float *med, float *mad )
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 }
00127 
00128 /*---------------------------------------------------------------*/
00129 
00130 static float median_float4(float a, float b, float c, float d)
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 }
00144 
00145 /*=======================================================================
00146  * Median sorting:
00147  * These highly optimized functions are used to find the median
00148  * value out of an array of 5, 7, or 9 pixel values.
00149  * They are a simple implementation of the minimum number of comparisons
00150  * required to find the median.
00151  *
00152  * Reference used: Implementing median filters in XC4000E FPGA's
00153  * John L. Smith <jsmith@univision.com>
00154  * http://www.xilinx.com/xcell/xcell23.htm/xl23_16.pdf
00155  *
00156  * N. Devillard <nDevil@eso.org> - September 1997
00157  * Please distribute this code freely -- no warranty
00158  *=======================================================================*/
00159 
00160 /*-----------------------------------------------------------------------
00161  * Function: median_float5(), median_float7(), median_float9()
00162  * In      : A pointer to an array containing resp. 5, 7, or 9 pixel values
00163  * Out     : The median pixel value of the input array
00164  * Job     : Fast computation of the median of 5, 7, or 9 values
00165  * Notice  : The input array is modified: partly sorted so that the
00166  *            middle element is the median
00167  *            No bound checking to gain time, it is up to the caller
00168  *            to make sure arrays are allocated
00169  *---------------------------------------------------------------------*/
00170 
00171 #undef  SORT2
00172 #define SORT2(a,b) if(a>b) SWAP(a,b);
00173 
00174 static float median_float5(float *p)
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 }
00181 
00182 static float median_float7(float *p)
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 }
00191 
00192 static float median_float9(float *p)
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 }