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 MRI_IMAGE * mri_quantize( int newcolors , MRI_IMAGE * im ) ;
00010 
00011 int main( int argc , char * argv[] )
00012 {
00013    MRI_IMAGE * fred ;
00014    fred = mri_read_ppm( "fred.pnm" ) ;
00015    mri_quantize( 100 , fred ) ;
00016    exit(0) ;
00017 }
00018 
00019 #define RGB_TO_INT(r,g,b) ((((int)(r))<<16) | (((int)(g))<< 8) | ((int)(b)))
00020 #define INT_TO_RRR(i)     (((i) & 0xff0000) >> 16)
00021 #define INT_TO_GGG(i)     (((i) & 0xff00)   >>  8)
00022 #define INT_TO_BBB(i)     ( (i) & 0xff)
00023 
00024 #define RGB_MASK          RGB_TO_INT( 0xff , 0xff , 0xff )
00025 #define RRR_MASK          RGB_TO_INT( 0xff , 0x00 , 0x00 )
00026 #define GGG_MASK          RGB_TO_INT( 0x00 , 0xff , 0x00 )
00027 #define BBB_MASK          RGB_TO_INT( 0x00 , 0x00 , 0xff )
00028 
00029 #define MAXCOL 32767
00030 
00031 void qsort_int( int , int * , int ) ;
00032 
00033 MRI_IMAGE * mri_quantize( int newcolors , MRI_IMAGE * im )
00034 {
00035    MRI_IMAGE * intim , * outim ;
00036    int * intar ;
00037    byte * inar , * outar ;
00038    byte mask ;
00039    int ii , ncol , imask , smask ;
00040    int * color , * count ;
00041 
00042    /** sanity checks **/
00043 
00044    if( im == NULL             || newcolors < 2      ||   
00045       MRI_RGB_PTR(im) == NULL || im->kind != MRI_rgb  ) return NULL ;
00046 
00047    /** make copy of input image as ints **/
00048 
00049    inar  = MRI_RGB_PTR(im) ;
00050    intim = mri_new_conforming( im , MRI_int ) ;
00051    intar = MRI_INT_PTR(intim) ;
00052 
00053    for( ii=0 ; ii < im->nvox ; ii++ )
00054       intar[ii] = RGB_TO_INT( inar[3*ii] , inar[3*ii+1] , inar[3*ii+2] ) ;
00055 
00056    /** sort copy, count unique colors **/
00057 
00058    qsort_int( im->nvox , intar , RGB_MASK ) ;
00059    ncol = 1 ;
00060    for( ii=1 ; ii < im->nvox ; ii++ )
00061       if( intar[ii] != intar[ii-1] ) ncol++ ;
00062 
00063    /** if too many colors, get rid of some by masking out low order bits **/
00064 
00065 fprintf(stderr,"%d colors in input image\n",ncol) ;
00066 
00067    imask = 0 ;
00068    while( ncol > MAXCOL ){
00069       imask++ ; mask = (byte)( 256 - (1<<imask) ) ;
00070       smask = RGB_TO_INT( mask , mask , mask ) ;
00071 
00072       for( ii=0 ; ii < im->nvox ; ii++ )
00073 #if 0
00074          intar[ii] &= smask ;
00075 #else
00076          intar[ii] = RGB_TO_INT( inar[3*ii]   & mask ,
00077                                  inar[3*ii+1] & mask , inar[3*ii+2] & mask ) ;
00078 #endif
00079 
00080       qsort_int( im->nvox , intar , RGB_MASK ) ;
00081       ncol = 1 ;
00082       for( ii=1 ; ii < im->nvox ; ii++ )
00083          if( intar[ii] != intar[ii-1] ) ncol++ ;
00084 
00085 fprintf(stderr,"%d colors in input image with mask=%d\n",ncol,(int)mask) ;
00086    }
00087 
00088    /** make color histogram **/
00089 
00090    color = (int *) malloc( sizeof(int) * ncol ) ;
00091    count = (int *) malloc( sizeof(int) * ncol ) ;
00092 
00093    color[0] = intar[0] ; count[0] = 1 ; ncol = 0 ;
00094    for( ii=1 ; ii < im->nvox ; ii++ ){
00095       if( intar[ii] != intar[ii-1] ){
00096          ncol++ ;
00097          color[ncol] = intar[ii] ;
00098          count[ncol] = 1 ;
00099       } else {
00100          count[ncol]++ ;
00101       }
00102    }
00103    ncol++ ;
00104    mri_free( intim ) ;
00105 }
00106 
00107 /******************************************************************************/
00108 
00109 #define QS_CUTOFF     40       /* cutoff to switch from qsort to isort */
00110 #define QS_STACK      4096     /* qsort stack size */
00111 #define QS_SWAP(x,y)  (temp=(x), (x)=(y),(y)=temp)
00112 
00113 static int stack[QS_STACK] ;  /* stack for qsort "recursion" */
00114 
00115 /*----------------------------------------------------------------------------*/
00116 /*------------- insertion sort : sort an array of int in-place ---------------*/
00117 
00118 #define ILESS(a,b)  ( ((a)&smask) <  ((b)&smask) )
00119 #define IMORE(a,b)  ( ((a)&smask) >  ((b)&smask) )
00120 #define IEQUAL(a,b) ( ((a)&smask) == ((b)&smask) )
00121 
00122 void isort_int( int n , int * ar , int smask )
00123 {
00124    register int  j , p ;  /* array indices */
00125    register int temp ;  /* a[j] holding place */
00126    register int * a = ar ;
00127 
00128    if( n < 2 || ar == NULL ) return ;
00129 
00130    for( j=1 ; j < n ; j++ ){
00131 
00132      if( ILESS(a[j],a[j-1]) ){  /* out of order */
00133         p    = j ;
00134         temp = a[j] ;
00135         do{
00136           a[p] = a[p-1] ;       /* at this point, a[p-1] > temp, so move it up */
00137           p-- ;
00138         } while( p > 0 && ILESS(temp,a[p-1]) ) ;
00139         a[p] = temp ;           /* finally, put temp in its place */
00140      }
00141    }
00142    return ;
00143 }
00144 
00145 /*--------- qsrec : recursive part of quicksort (stack implementation) ----------*/
00146 
00147 void qsrec_int( int n , int * ar , int cutoff , int smask )
00148 {
00149    register int i , j ;         /* scanning indices */
00150    register int temp , pivot ;  /* holding places */
00151    register int * a = ar ;
00152 
00153    int left , right , mst ;
00154 
00155    /* return if too short (insertion sort will clean up) */
00156 
00157    if( cutoff < 3 ) cutoff = 3 ;
00158    if( n < cutoff || ar == NULL ) return ;
00159 
00160    /* initialize stack to start with whole array */
00161 
00162    stack[0] = 0   ;
00163    stack[1] = n-1 ;
00164    mst      = 2   ;
00165 
00166    /* loop while the stack is nonempty */
00167 
00168    while( mst > 0 ){
00169       right = stack[--mst] ;  /* work on subarray from left -> right */
00170       left  = stack[--mst] ;
00171 
00172       i = ( left + right ) / 2 ;           /* middle of subarray */
00173 
00174       /*----- sort the left, middle, and right a[]'s -----*/
00175 
00176       if( IMORE(a[left],a[i]    ) ) QS_SWAP( a[left]  , a[i]     ) ;
00177       if( IMORE(a[left],a[right]) ) QS_SWAP( a[left]  , a[right] ) ;
00178       if( IMORE(a[i]   ,a[right]) ) QS_SWAP( a[right] , a[i]     ) ;
00179 
00180       pivot = a[i] ;                       /* a[i] is the median-of-3 pivot! */
00181       a[i]  = a[right] ;
00182 
00183       i = left ; j = right ;               /* initialize scanning */
00184 
00185       /*----- partition:  move elements bigger than pivot up and elements
00186                           smaller than pivot down, scanning in from ends -----*/
00187 
00188       do{
00189         for( ; ILESS(a[++i],pivot) ; ) ;  /* scan i up,   until a[i] >= pivot */
00190         for( ; IMORE(a[--j],pivot) ; ) ;  /* scan j down, until a[j] <= pivot */
00191 
00192         if( j <= i ) break ;         /* if j meets i, quit */
00193 
00194         QS_SWAP( a[i] , a[j] ) ;
00195       } while( 1 ) ;
00196 
00197       /*----- at this point, the array is partitioned -----*/
00198 
00199       a[right] = a[i] ; a[i] = pivot ;  /* restore the pivot */
00200 
00201       /*----- signal the subarrays that need further work -----*/
00202 
00203       if( (i-left)  > cutoff ){ stack[mst++] = left ; stack[mst++] = i-1   ; }
00204       if( (right-i) > cutoff ){ stack[mst++] = i+1  ; stack[mst++] = right ; }
00205 
00206    }  /* end of while stack is non-empty */
00207    return ;
00208 }
00209 
00210 /* quick_sort :  sort an array partially recursively, and partially insertion */
00211 
00212 void qsort_int( int n , int * a , int smask )
00213 {
00214    qsrec_int( n , a , QS_CUTOFF , smask ) ;
00215    isort_int( n , a , smask ) ;
00216    return ;
00217 }