00001 #include "qsort.h"
00002 
00003 #ifdef DEBUG
00004 #include <time.h>
00005 #define FCLOCK ((float)clock()/((float)CLOCKS_PER_SEC))
00006 #endif
00007 
00008 /********************************************************************************/
00009 /* insertion_sort : sort an array of short */
00010 
00011 void isort_sh( int n , short * ar )
00012 {
00013    register int  j , p ;  /* array indices */
00014    register short temp ;  /* a[j] holding place */
00015    register short * a = ar ;
00016 
00017    if( n < 2 ) return ;
00018 
00019    for( j=1 ; j < n ; j++ ){
00020 
00021      if( a[j] < a[j-1] ){  /* out of order */
00022         p    = j ;
00023         temp = a[j] ;
00024 
00025        do{
00026           a[p] = a[p-1] ;       /*at this point, a[p-1] > temp, so move it up*/
00027          p-- ;
00028         } while( p > 0 && temp < a[p-1] ) ;
00029 
00030         a[p] = temp ;           /*finally, put temp in its place*/
00031      }
00032    }
00033 }
00034 
00035 
00036 /********************************************************************************/
00037 /* qsrec : recursive part of quicksort (stack implementation) */
00038 
00039 #define QS_STACK  1024  /* stack size */
00040 #define QS_SWAP(x,y) (temp=(x),(x)=(y),(y)=temp)
00041 
00042 void qsrec_sh( int n , short * ar , int cutoff )
00043 {
00044    register int i , j ;           /* scanning indices */
00045    register short temp , pivot ;  /* holding places */
00046    register short * a = ar ;
00047 
00048    int left , right , mst , stack[QS_STACK] ;
00049 
00050 #ifdef DEBUG
00051    int max_st = 2 ;
00052 #endif
00053 
00054    /* return if too short (insertion sort will clean up) */
00055 
00056    if( cutoff < 3 ) cutoff = 3 ;
00057    if( n < cutoff ) return ;
00058 
00059    /* initialize stack to start with whole array */
00060 
00061    stack[0] = 0   ;
00062    stack[1] = n-1 ;
00063    mst      = 2   ;
00064 
00065    /* loop while the stack is nonempty */
00066 
00067    while( mst > 0 ){
00068       right = stack[--mst] ;  /* work on subarray from left -> right */
00069       left  = stack[--mst] ;
00070 
00071       i = ( left + right ) / 2 ;           /* middle of subarray */
00072 
00073       /* sort the left, middle, and right a[]'s */
00074 
00075       if( a[left] > a[i]     ) QS_SWAP( a[left]  , a[i]     ) ;
00076       if( a[left] > a[right] ) QS_SWAP( a[left]  , a[right] ) ;
00077       if( a[i] > a[right]    ) QS_SWAP( a[right] , a[i]     ) ;
00078 
00079       pivot = a[i] ;                        /* a[i] is the median-of-3 pivot! */
00080       a[i]  = a[right] ;
00081 
00082       i = left ;                            /* initialize scanning */
00083       j = right ;
00084 
00085       /*----- partition:  move elements bigger than pivot up and elements
00086                           smaller than pivot down, scanning in from ends -----*/
00087 
00088       do{
00089         for( ; a[++i] < pivot ; ) ;  /* scan i up,   until a[i] >= pivot */
00090         for( ; a[--j] > pivot ; ) ;  /* scan j down, until a[j] <= pivot */
00091 
00092         if( j <= i ) break ;         /* if j meets i, quit */
00093 
00094         QS_SWAP( a[i] , a[j] ) ;
00095       } while( 1 ) ;
00096 
00097       /*----- at this point, the array is partitioned -----*/
00098 
00099       a[right] = a[i] ;           /*restore the pivot*/
00100       a[i]     = pivot ;
00101 
00102       if( (i-left)  > cutoff ){ stack[mst++] = left ; stack[mst++] = i-1   ; }
00103       if( (right-i) > cutoff ){ stack[mst++] = i+1  ; stack[mst++] = right ; }
00104 
00105 #ifdef DEBUG
00106    if( mst > max_st ) max_st = mst ;
00107 #endif
00108 
00109    }  /* end of while stack is non-empty */
00110 
00111 #ifdef DEBUG
00112    printf("qsort: for n=%d max stack depth=%d\n",n,max_st) ;
00113 #endif
00114 
00115 } /*qsrec*/
00116 
00117 /********************************************************************************/
00118 /* quick_sort :  sort an array partially recursively, and partially insertion */
00119 
00120 #ifndef QS_CUTOFF
00121 #define QS_CUTOFF 40
00122 #endif
00123 
00124 void qsort_sh( int n , short * a )
00125 {
00126 #ifdef DEBUG
00127    float clk1 , clk2 , clk3 ;
00128    clk1 = FCLOCK ;
00129 #endif
00130 
00131    qsrec_sh( n , a , QS_CUTOFF ) ;
00132 
00133 #ifdef DEBUG
00134    clk2 = FCLOCK ;
00135 #endif
00136 
00137    isort_sh( n , a ) ;
00138 
00139 #ifdef DEBUG
00140    clk3 = FCLOCK ;
00141    printf("qsort_sh cpu: qsrec=%f isort=%f\n",clk2-clk1,clk3-clk2) ;
00142 #endif
00143 
00144 }
00145 
00146 void qsort_partial_sh( int n , short * a , int cutoff )
00147 {
00148    qsrec_sh( n , a , cutoff ) ;
00149    if( cutoff < 2 ) isort_sh( n , a ) ;
00150 }
00151 
00152 /*-----------------------------------------------------------------------
00153    compute an approximation to the inverse histogram of an array of shorts:
00154 
00155      n    = number of input shorts
00156      ar   = array of input shorts
00157      nbin = number of output bins (e.g., 101 for percentage points)
00158      ih   = array [nbin] of shorts:
00159             ih[k] = value that has k/(nbin-1) fraction of the
00160                     points in ar below it (approximately) for k=0..nbin-1;
00161             ih[0] = min of ar, ih[nbin-1] = max of ar
00162 -------------------------------------------------------------------------*/
00163 
00164 int MCW_inverse_histogram_sh( int n , short * ar , int nbin , short * ih )
00165 {
00166    int k , i , nwin,nhalf , sum,ibot,itop ;
00167    float fwin ;
00168    short * at ;
00169 
00170    if( n < 2 || nbin < 2 || ar == NULL || ih == NULL ) return 0 ;  /* bad inputs */
00171 
00172    fwin = ((float) n) / (float) nbin ;
00173    nwin = (int) fwin ; nhalf = (int) (0.5*fwin) ;
00174 
00175    at = (short *) malloc( n * sizeof(short) ) ; if( at == NULL ) return 0 ;
00176    memcpy( at , ar , n * sizeof(short) ) ;
00177 
00178    qsort_partial_sh( n , at , nhalf ) ;
00179 
00180    ih[0] = at[0] ;
00181    for( i=1 ; i < nwin ; i++ ) if( ih[0] > at[i] ) ih[0] = at[i] ;
00182 
00183    ih[nbin-1] = at[n-1] ;
00184    for( i=1 ; i < nwin ; i++ ) if( ih[nbin-1] < at[n-1-i] ) ih[nbin-1] = at[n-1-i] ;
00185 
00186    for( k=1 ; k < nbin-1 ; k++ ){
00187 
00188       sum  = 0 ;
00189       ibot = (k-0.5) * fwin ;
00190       itop = (k+0.5) * fwin + 1.49 ;
00191       if( itop > nbin ) itop = nbin ; if( itop <= ibot ) itop = ibot+1 ;
00192 
00193       for( i=ibot ; i < itop ; i++ ) sum += at[i] ;
00194 
00195       ih[k] = (short)( sum / (itop-ibot) ) ;
00196 #ifdef DEBUG
00197    printf("MCW_inverse_historgram: k=%d ibot=%d itop=%d sum=%d ih[k]=%d\n",
00198           k,ibot,itop,ih[k]) ;
00199 #endif
00200    }
00201 
00202    free( at ) ;
00203    return 1 ;
00204 }