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 /*** 7D SAFE ***/
00010 
00011 /* prototypes for internal sorting routines:
00012      array of shorts
00013      array of ints
00014      array of floats
00015      array of floats, with an integer array carried along in the swaps */
00016 
00017 void qsrec_short( int , short * , int ) ;
00018 void qsrec_int  ( int , int   * , int ) ;
00019 void qsrec_float( int , float * , int ) ;
00020 void qsrec_pair ( int , float * , int * , int ) ;
00021 
00022 #define QS_CUTOFF     40       /* cutoff to switch from qsort to isort */
00023 #define QS_STACK      4096     /* qsort stack size */
00024 #define QS_SWAP(x,y)  (temp=(x), (x)=(y),(y)=temp)
00025 
00026 static int stack[QS_STACK] ;  /* stack for qsort "recursion" */
00027 
00028 /***************************************************************************
00029      Each qsort_TYPE routine (TYPE=short, int, float, or pair) has two
00030      pieces.  The isort_TYPE routine does an insertion sort routine,
00031      which is only fast for nearly sorted routines.  The qsrec_TYPE
00032      routine carries out the quicksort algorithm down to the level
00033      QS_CUTOFF -- at that point, the array is nearly sorted.
00034 
00035      In the qsrec_TYPE routines, the fundamental operation is
00036      the SWAP of two items.  In the isort_TYPE routines, the
00037      value at the j index is determined to be out of place, so
00038      it is stored in a temporary variable and the values below
00039      it are moved up in the array until the proper place for
00040      the former a[j] is found.  Then it is stored where it belongs.
00041      Compare the isort_short and isort_pair routines to see how
00042      this should be generalized to more complex objects.
00043 
00044                                                       -- Robert W. Cox
00045 ***************************************************************************/
00046 
00047 /*------------------------------------------------------------------------------*/
00048 /*------------- insertion sort : sort an array of short in-place ---------------*/
00049 
00050 void isort_short( int n , short * ar )
00051 {
00052    register int  j , p ;  /* array indices */
00053    register short temp ;  /* a[j] holding place */
00054    register short * a = ar ;
00055 
00056    if( n < 2 || ar == NULL ) return ;
00057 
00058    for( j=1 ; j < n ; j++ ){
00059 
00060      if( a[j] < a[j-1] ){  /* out of order */
00061         p    = j ;
00062         temp = a[j] ;
00063         do{
00064           a[p] = a[p-1] ;       /* at this point, a[p-1] > temp, so move it up */
00065           p-- ;
00066         } while( p > 0 && temp < a[p-1] ) ;
00067         a[p] = temp ;           /* finally, put temp in its place */
00068      }
00069    }
00070    return ;
00071 }
00072 
00073 /*--------- qsrec : recursive part of quicksort (stack implementation) ----------*/
00074 
00075 void qsrec_short( int n , short * ar , int cutoff )
00076 {
00077    register int i , j ;           /* scanning indices */
00078    register short temp , pivot ;  /* holding places */
00079    register short * a = ar ;
00080 
00081    int left , right , mst ;
00082 
00083    /* return if too short (insertion sort will clean up) */
00084 
00085    if( cutoff < 3 ) cutoff = 3 ;
00086    if( n < cutoff || ar == NULL ) return ;
00087 
00088    /* initialize stack to start with whole array */
00089 
00090    stack[0] = 0   ;
00091    stack[1] = n-1 ;
00092    mst      = 2   ;
00093 
00094    /* loop while the stack is nonempty */
00095 
00096    while( mst > 0 ){
00097       right = stack[--mst] ;  /* work on subarray from left -> right */
00098       left  = stack[--mst] ;
00099 
00100       i = ( left + right ) / 2 ;           /* middle of subarray */
00101 
00102       /*----- sort the left, middle, and right a[]'s -----*/
00103 
00104       if( a[left] > a[i]     ) QS_SWAP( a[left]  , a[i]     ) ;
00105       if( a[left] > a[right] ) QS_SWAP( a[left]  , a[right] ) ;
00106       if( a[i] > a[right]    ) QS_SWAP( a[right] , a[i]     ) ;
00107 
00108       pivot = a[i] ;                       /* a[i] is the median-of-3 pivot! */
00109       a[i]  = a[right] ;
00110 
00111       i = left ; j = right ;               /* initialize scanning */
00112 
00113       /*----- partition:  move elements bigger than pivot up and elements
00114                           smaller than pivot down, scanning in from ends -----*/
00115 
00116       do{
00117         for( ; a[++i] < pivot ; ) ;  /* scan i up,   until a[i] >= pivot */
00118         for( ; a[--j] > pivot ; ) ;  /* scan j down, until a[j] <= pivot */
00119 
00120         if( j <= i ) break ;         /* if j meets i, quit */
00121 
00122         QS_SWAP( a[i] , a[j] ) ;
00123       } while( 1 ) ;
00124 
00125       /*----- at this point, the array is partitioned -----*/
00126 
00127       a[right] = a[i] ; a[i] = pivot ;  /* restore the pivot */
00128 
00129       /*----- signal the subarrays that need further work -----*/
00130 
00131       if( (i-left)  > cutoff ){ stack[mst++] = left ; stack[mst++] = i-1   ; }
00132       if( (right-i) > cutoff ){ stack[mst++] = i+1  ; stack[mst++] = right ; }
00133 
00134    }  /* end of while stack is non-empty */
00135    return ;
00136 }
00137 
00138 /* quick_sort :  sort an array partially recursively, and partially insertion */
00139 
00140 void qsort_short( int n , short * a )
00141 {
00142    qsrec_short( n , a , QS_CUTOFF ) ;
00143    isort_short( n , a ) ;
00144    return ;
00145 }
00146 
00147 /*----------------------------------------------------------------------------*/
00148 /*------------- insertion sort : sort an array of int in-place ---------------*/
00149 
00150 void isort_int( int n , int * ar )
00151 {
00152    register int  j , p ;  /* array indices */
00153    register int temp ;    /* a[j] holding place */
00154    register int * a = ar ;
00155 
00156    if( n < 2 || ar == NULL ) return ;
00157 
00158    for( j=1 ; j < n ; j++ ){
00159 
00160      if( a[j] < a[j-1] ){  /* out of order */
00161         p    = j ;
00162         temp = a[j] ;
00163         do{
00164           a[p] = a[p-1] ;       /* at this point, a[p-1] > temp, so move it up */
00165           p-- ;
00166         } while( p > 0 && temp < a[p-1] ) ;
00167         a[p] = temp ;           /* finally, put temp in its place */
00168      }
00169    }
00170    return ;
00171 }
00172 
00173 /*--------- qsrec : recursive part of quicksort (stack implementation) ----------*/
00174 
00175 void qsrec_int( int n , int * ar , int cutoff )
00176 {
00177    register int i , j ;         /* scanning indices */
00178    register int temp , pivot ;  /* holding places */
00179    register int * a = ar ;
00180 
00181    int left , right , mst ;
00182 
00183    /* return if too short (insertion sort will clean up) */
00184 
00185    if( cutoff < 3 ) cutoff = 3 ;
00186    if( n < cutoff || ar == NULL ) return ;
00187 
00188    /* initialize stack to start with whole array */
00189 
00190    stack[0] = 0   ;
00191    stack[1] = n-1 ;
00192    mst      = 2   ;
00193 
00194    /* loop while the stack is nonempty */
00195 
00196    while( mst > 0 ){
00197       right = stack[--mst] ;  /* work on subarray from left -> right */
00198       left  = stack[--mst] ;
00199 
00200       i = ( left + right ) / 2 ;           /* middle of subarray */
00201 
00202       /*----- sort the left, middle, and right a[]'s -----*/
00203 
00204       if( a[left] > a[i]     ) QS_SWAP( a[left]  , a[i]     ) ;
00205       if( a[left] > a[right] ) QS_SWAP( a[left]  , a[right] ) ;
00206       if( a[i] > a[right]    ) QS_SWAP( a[right] , a[i]     ) ;
00207 
00208       pivot = a[i] ;                       /* a[i] is the median-of-3 pivot! */
00209       a[i]  = a[right] ;
00210 
00211       i = left ; j = right ;               /* initialize scanning */
00212 
00213       /*----- partition:  move elements bigger than pivot up and elements
00214                           smaller than pivot down, scanning in from ends -----*/
00215 
00216       do{
00217         for( ; a[++i] < pivot ; ) ;  /* scan i up,   until a[i] >= pivot */
00218         for( ; a[--j] > pivot ; ) ;  /* scan j down, until a[j] <= pivot */
00219 
00220         if( j <= i ) break ;         /* if j meets i, quit */
00221 
00222         QS_SWAP( a[i] , a[j] ) ;
00223       } while( 1 ) ;
00224 
00225       /*----- at this point, the array is partitioned -----*/
00226 
00227       a[right] = a[i] ; a[i] = pivot ;  /* restore the pivot */
00228 
00229       /*----- signal the subarrays that need further work -----*/
00230 
00231       if( (i-left)  > cutoff ){ stack[mst++] = left ; stack[mst++] = i-1   ; }
00232       if( (right-i) > cutoff ){ stack[mst++] = i+1  ; stack[mst++] = right ; }
00233 
00234    }  /* end of while stack is non-empty */
00235    return ;
00236 }
00237 
00238 /* quick_sort:  sort an array partially recursively, and partially insertion */
00239 
00240 void qsort_int( int n , int * a )
00241 {
00242    qsrec_int( n , a , QS_CUTOFF ) ;
00243    isort_int( n , a ) ;
00244    return ;
00245 }
00246 
00247 /*------------------------------------------------------------------------------*/
00248 /*------------- insertion sort : sort an array of float in-place ---------------*/
00249 
00250 void isort_float( int n , float * ar )
00251 {
00252    register int  j , p ;  /* array indices */
00253    register float temp ;  /* a[j] holding place */
00254    register float * a = ar ;
00255 
00256    if( n < 2 || ar == NULL ) return ;
00257 
00258    for( j=1 ; j < n ; j++ ){
00259 
00260      if( a[j] < a[j-1] ){  /* out of order */
00261         p    = j ;
00262         temp = a[j] ;
00263         do{
00264           a[p] = a[p-1] ;       /* at this point, a[p-1] > temp, so move it up */
00265           p-- ;
00266         } while( p > 0 && temp < a[p-1] ) ;
00267         a[p] = temp ;           /* finally, put temp in its place */
00268      }
00269    }
00270    return ;
00271 }
00272 
00273 /*--------- qsrec : recursive part of quicksort (stack implementation) ----------*/
00274 
00275 void qsrec_float( int n , float * ar , int cutoff )
00276 {
00277    register int i , j ;           /* scanning indices */
00278    register float temp , pivot ;  /* holding places */
00279    register float * a = ar ;
00280 
00281    int left , right , mst ;
00282 
00283    /* return if too short (insertion sort will clean up) */
00284 
00285    if( cutoff < 3 ) cutoff = 3 ;
00286    if( n < cutoff || ar == NULL ) return ;
00287 
00288    /* initialize stack to start with whole array */
00289 
00290    stack[0] = 0   ;
00291    stack[1] = n-1 ;
00292    mst      = 2   ;
00293 
00294    /* loop while the stack is nonempty */
00295 
00296    while( mst > 0 ){
00297       right = stack[--mst] ;  /* work on subarray from left -> right */
00298       left  = stack[--mst] ;
00299 
00300       i = ( left + right ) / 2 ;           /* middle of subarray */
00301 
00302       /*----- sort the left, middle, and right a[]'s -----*/
00303 
00304       if( a[left] > a[i]     ) QS_SWAP( a[left]  , a[i]     ) ;
00305       if( a[left] > a[right] ) QS_SWAP( a[left]  , a[right] ) ;
00306       if( a[i] > a[right]    ) QS_SWAP( a[right] , a[i]     ) ;
00307 
00308       pivot = a[i] ;                       /* a[i] is the median-of-3 pivot! */
00309       a[i]  = a[right] ;
00310 
00311       i = left ; j = right ;               /* initialize scanning */
00312 
00313       /*----- partition:  move elements bigger than pivot up and elements
00314                           smaller than pivot down, scanning in from ends -----*/
00315 
00316       do{
00317         for( ; a[++i] < pivot ; ) ;  /* scan i up,   until a[i] >= pivot */
00318         for( ; a[--j] > pivot ; ) ;  /* scan j down, until a[j] <= pivot */
00319 
00320         if( j <= i ) break ;         /* if j meets i, quit */
00321 
00322         QS_SWAP( a[i] , a[j] ) ;
00323       } while( 1 ) ;
00324 
00325       /*----- at this point, the array is partitioned -----*/
00326 
00327       a[right] = a[i] ; a[i] = pivot ;  /* restore the pivot */
00328 
00329       /*----- signal the subarrays that need further work -----*/
00330 
00331       if( (i-left)  > cutoff ){ stack[mst++] = left ; stack[mst++] = i-1   ; }
00332       if( (right-i) > cutoff ){ stack[mst++] = i+1  ; stack[mst++] = right ; }
00333 
00334    }  /* end of while stack is non-empty */
00335    return ;
00336 }
00337 
00338 /* quick_sort :  sort an array partially recursively, and partially insertion */
00339 
00340 void qsort_float( int n , float * a )
00341 {
00342    qsrec_float( n , a , QS_CUTOFF ) ;
00343    isort_float( n , a ) ;
00344    return ;
00345 }
00346 
00347 /*------------------------------------------------------------------------------*/
00348 /*--------------- insertion sort of a float-int pair of arrays -----------------*/
00349 
00350 void isort_pair( int n , float * ar , int * iar )
00351 {
00352    register int  j , p ;  /* array indices */
00353    register float temp ;  /* a[j] holding place */
00354    register int  itemp ;
00355    register float * a = ar ;
00356    register int   *ia = iar ;
00357 
00358    if( n < 2 || ar == NULL || iar == NULL ) return ;
00359 
00360    for( j=1 ; j < n ; j++ ){
00361 
00362      if( a[j] < a[j-1] ){  /* out of order */
00363         p    = j ;
00364         temp = a[j] ; itemp = ia[j] ;
00365         do{
00366           a[p] = a[p-1] ; ia[p] = ia[p-1] ;
00367           p-- ;
00368         } while( p > 0 && temp < a[p-1] ) ;
00369         a[p] = temp ; ia[p] = itemp ;
00370      }
00371    }
00372    return ;
00373 }
00374 
00375 /*--------- qsrec : recursive part of quicksort (stack implementation) ----------*/
00376 
00377 #define QS_ISWAP(x,y) (itemp=(x),(x)=(y),(y)=itemp)
00378 
00379 void qsrec_pair( int n , float * ar , int * iar , int cutoff )
00380 {
00381    register int i , j ;           /* scanning indices */
00382    register float temp , pivot ;  /* holding places */
00383    register int  itemp ,ipivot ;
00384    register float * a = ar ;
00385    register int   *ia = iar ;
00386 
00387    int left , right , mst ;
00388 
00389    /* return if too short (insertion sort will clean up) */
00390 
00391    if( cutoff < 3 ) cutoff = 3 ;
00392    if( n < cutoff || ar == NULL || iar == NULL ) return ;
00393 
00394    /* initialize stack to start with whole array */
00395 
00396    stack[0] = 0   ;
00397    stack[1] = n-1 ;
00398    mst      = 2   ;
00399 
00400    /* loop while the stack is nonempty */
00401 
00402    while( mst > 0 ){
00403       right = stack[--mst] ;  /* work on subarray from left -> right */
00404       left  = stack[--mst] ;
00405 
00406       i = ( left + right ) / 2 ;           /* middle of subarray */
00407 
00408       /*----- sort the left, middle, and right a[]'s -----*/
00409 
00410       if( a[left] > a[i]     ){ QS_SWAP ( a[left]  , a[i]     ) ;
00411                                 QS_ISWAP(ia[left]  ,ia[i]     ) ; }
00412       if( a[left] > a[right] ){ QS_SWAP ( a[left]  , a[right] ) ;
00413                                 QS_ISWAP(ia[left]  ,ia[right] ) ; }
00414       if( a[i] > a[right]    ){ QS_SWAP ( a[right] , a[i]     ) ;
00415                                 QS_ISWAP(ia[right] ,ia[i]     ) ; }
00416 
00417        pivot = a[i] ;  a[i] = a[right] ;
00418       ipivot =ia[i] ; ia[i] =ia[right] ;
00419 
00420       i = left ; j = right ;               /* initialize scanning */
00421 
00422       /*----- partition:  move elements bigger than pivot up and elements
00423                           smaller than pivot down, scanning in from ends -----*/
00424 
00425       do{
00426         for( ; a[++i] < pivot ; ) ;  /* scan i up,   until a[i] >= pivot */
00427         for( ; a[--j] > pivot ; ) ;  /* scan j down, until a[j] <= pivot */
00428 
00429         if( j <= i ) break ;         /* if j meets i, quit */
00430 
00431         QS_SWAP ( a[i] , a[j] ) ;
00432         QS_ISWAP(ia[i] ,ia[j] ) ;
00433       } while( 1 ) ;
00434 
00435       /*----- at this point, the array is partitioned -----*/
00436 
00437       a[right] = a[i] ; a[i] = pivot ;  /* restore the pivot */
00438       ia[right]=ia[i] ;ia[i] =ipivot ;
00439 
00440       /*----- signal the subarrays that need further work -----*/
00441 
00442       if( (i-left)  > cutoff ){ stack[mst++] = left ; stack[mst++] = i-1   ; }
00443       if( (right-i) > cutoff ){ stack[mst++] = i+1  ; stack[mst++] = right ; }
00444 
00445    }  /* end of while stack is non-empty */
00446    return ;
00447 }
00448 
00449 /* quick_sort :  sort an array partially recursively, and partially insertion */
00450 
00451 void qsort_pair( int n , float * a , int * ia )
00452 {
00453    qsrec_pair( n , a , ia , QS_CUTOFF ) ;
00454    isort_pair( n , a , ia ) ;
00455    return ;
00456 }
00457 
00458 /*******************************************************************************
00459   Compute the "percentage points" of the histogram of an image.
00460   "per" should be a pre-allocated array of "nper+1" floats; per[i] will be set
00461   to the image intensity below which fraction i/nper of the pixels fall,
00462   for i=0..nper.
00463   N.B.:  per[0] = image minimum and per[nper] = maximum.
00464         "per" is float, no matter what the input image type is.
00465 ********************************************************************************/
00466 
00467 void mri_percents( MRI_IMAGE * im , int nper , float per[] )
00468 {
00469    register int pp , ii , nvox ;
00470    register float fi , frac ;
00471 
00472    /*** sanity checks ***/
00473 
00474    if( im == NULL || per == NULL || nper < 2 ) return ;
00475 
00476 #ifdef MRILIB_7D
00477    nvox = im->nvox ;
00478 #else
00479    nvox = im->nx * im->ny ;
00480 #endif
00481    frac = nvox / ((float) nper) ;
00482 
00483    switch( im->kind ){
00484 
00485       /*** create a float image copy of the data,
00486            sort it, then interpolate the percentage points ***/
00487 
00488       default:{
00489          MRI_IMAGE * inim ;
00490          float * far ;
00491 
00492          inim = mri_to_float( im ) ;
00493          far  = MRI_FLOAT_PTR(inim) ;
00494          qsort_float( nvox , far ) ;
00495 
00496          per[0] = far[0] ;
00497          for( pp=1 ; pp < nper ; pp++ ){
00498             fi = frac * pp ; ii = fi ; fi = fi - ii ;
00499             per[pp] = (1.0-fi) * far[ii] + fi * far[ii+1] ;
00500          }
00501          per[nper] = far[nvox-1] ;
00502          mri_free( inim ) ;
00503       }
00504       break ;
00505 
00506       /*** create a short image copy of the data,
00507            sort it, then interpolate the percentage points ***/
00508 
00509       case MRI_short:
00510       case MRI_byte:{
00511          MRI_IMAGE * inim ;
00512          short * sar ;
00513 
00514          inim = mri_to_short( 1.0 , im ) ;
00515          sar  = MRI_SHORT_PTR(inim) ;
00516          qsort_short( nvox , sar ) ;
00517 
00518          per[0] = sar[0] ;
00519          for( pp=1 ; pp < nper ; pp++ ){
00520             fi = frac * pp ; ii = fi ; fi = fi - ii ;
00521             per[pp] = (1.0-fi) * sar[ii] + fi * sar[ii+1] ;
00522          }
00523          per[nper] = sar[nvox-1] ;
00524          mri_free( inim ) ;
00525       }
00526    }
00527 
00528    return ;
00529 }
00530 
00531 /****************************************************************************
00532    Produce a "histogram flattened" version of the input image.  That is, the
00533    output value at each pixel will be proportional to its place in the
00534    histogram of intensities.  That is, if the value at a pixel is the 173rd
00535    in intensity (0th is darkest) out of 1000 pixels total, then the flattened
00536    image value will be 0.173.  The output image is in MRI_float format.
00537 *****************************************************************************/
00538 
00539 MRI_IMAGE * mri_flatten( MRI_IMAGE * im )
00540 {
00541    MRI_IMAGE * flim , * intim , * outim ;
00542    float * far , * outar ;
00543    int * iar ;
00544    int ii , nvox , ibot,itop , nvox1 ;
00545    float fac , val ;
00546 
00547 #ifdef DEBUG
00548 printf("Entry: mri_flatten\n") ;
00549 #endif
00550 
00551    if( im == NULL ) return NULL ;
00552 
00553    /*** make an image that is just the voxel index in its array ***/
00554    /*** also, make the output image while we are at it          ***/
00555 
00556 #ifdef MRILIB_7D
00557    nvox  = im->nvox ;
00558    intim = mri_new_conforming( im , MRI_int ) ;
00559    outim = mri_new_conforming( im , MRI_float ) ;
00560 #else
00561    nvox  = im->nx * im->ny ;
00562    intim = mri_new( im->nx , im->ny , MRI_int ) ;
00563    outim = mri_new( im->nx , im->ny , MRI_float ) ;
00564 #endif
00565 
00566    iar = MRI_INT_PTR(intim) ; outar = MRI_FLOAT_PTR(outim) ;
00567 
00568    for( ii=0 ; ii < nvox ; ii++ ) iar[ii] = ii ;
00569 
00570    /*** copy the input data to a floating point image ***/
00571 
00572    flim = mri_to_float( im ) ; far = MRI_FLOAT_PTR(flim) ;
00573 
00574    /*** sort this image, with the index array being carried along
00575         so that we know where every pixel came from originally  ***/
00576 
00577    qsort_pair( nvox , far , iar ) ;
00578 
00579    /*** The "far" array is now sorted.  Thus, if the pixel that was in
00580         voxel i is now in voxel j, then its place in the histogram is
00581         j/nvox.  The only difficulty is that there may be ties.  We need
00582         to resolve these ties so that pixels with the same intensity
00583         don't get different output values.  We do this by scanning
00584         through far, finding blocks of equal values, and replacing
00585         them by their average position in the histogram.
00586    ***/
00587 
00588    fac = 1.0 / nvox ; nvox1 = nvox - 1 ;
00589 
00590    for( ibot=0 ; ibot < nvox1 ; ){
00591 
00592       /** if this value is unique, just set the value and move on **/
00593 
00594       val = far[ibot] ; itop = ibot+1 ;
00595       if( val != far[itop] ){
00596          far[ibot] = fac * ibot ;
00597          ibot++ ; continue ;
00598       }
00599 
00600       /** scan itop up until value is distinct **/
00601 
00602       for( ; itop < nvox1 && val == far[itop] ; itop++ ) ; /* nada */
00603 
00604       val = 0.5*fac * (ibot+itop-1) ;
00605       for( ii=ibot ; ii < itop ; ii++ ) far[ii] = val ;
00606       ibot = itop ;
00607    }
00608    far[nvox1] = 1.0 ;
00609 
00610    /*** now propagate these values back to the output image ***/
00611 
00612    for( ii=0 ; ii < nvox ; ii++ ) outar[iar[ii]] = far[ii] ;
00613 
00614    mri_free( flim ) ; mri_free( intim ) ;
00615 
00616    MRI_COPY_AUX( outim , im ) ;
00617    return outim ;
00618 }
00619 
00620 /*-------------------------------------------------------------
00621    Find the intensity in an image that is at the alpha-th
00622    quantile of the distribution.  That is, for 0 <= alpha <= 1,
00623    alpha*npix of the image values are below, and (1-alpha)*npix
00624    are above.  If alpha is 0, this is the minimum.  If alpha
00625    is 1, this is the maximum.
00626 ---------------------------------------------------------------*/
00627 
00628 float mri_quantile( MRI_IMAGE * im , float alpha )
00629 {
00630    int ii , nvox ;
00631    float fi , quan ;
00632 
00633    /*** sanity checks ***/
00634 
00635    if( im == NULL ) return 0.0 ;
00636 
00637    if( alpha <= 0.0 ) return (float) mri_min(im) ;
00638    if( alpha >= 1.0 ) return (float) mri_max(im) ;
00639 
00640 #ifdef MRILIB_7D
00641    nvox = im->nvox ;
00642 #else
00643    nvox = im->nx * im->ny ;
00644 #endif
00645 
00646    switch( im->kind ){
00647 
00648       /*** create a float image copy of the data,
00649            sort it, then interpolate the percentage points ***/
00650 
00651       default:{
00652          MRI_IMAGE * inim ;
00653          float * far ;
00654 
00655          inim = mri_to_float( im ) ;
00656          far  = MRI_FLOAT_PTR(inim) ;
00657          qsort_float( nvox , far ) ;
00658 
00659          fi   = alpha * nvox ;
00660          ii   = (int) fi ; if( ii >= nvox ) ii = nvox-1 ;
00661          fi   = fi - ii ;
00662          quan = (1.0-fi) * far[ii] + fi * far[ii+1] ;
00663          mri_free( inim ) ;
00664       }
00665       break ;
00666 
00667       /*** create a short image copy of the data,
00668            sort it, then interpolate the percentage points ***/
00669 
00670       case MRI_short:
00671       case MRI_byte:{
00672          MRI_IMAGE * inim ;
00673          short * sar ;
00674 
00675          inim = mri_to_short( 1.0 , im ) ;
00676          sar  = MRI_SHORT_PTR(inim) ;
00677          qsort_short( nvox , sar ) ;
00678 
00679          fi   = alpha * nvox ;
00680          ii   = (int) fi ; if( ii >= nvox ) ii = nvox-1 ;
00681          fi   = fi - ii ;
00682          quan = (1.0-fi) * sar[ii] + fi * sar[ii+1] ;
00683          mri_free( inim ) ;
00684       }
00685       break ;
00686    }
00687 
00688    return quan ;
00689 }