00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008
00009
00010
00011
00012
00013
00014
00015
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
00023 #define QS_STACK 4096
00024 #define QS_SWAP(x,y) (temp=(x), (x)=(y),(y)=temp)
00025
00026 static int stack[QS_STACK] ;
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050 void isort_short( int n , short * ar )
00051 {
00052 register int j , p ;
00053 register short temp ;
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] ){
00061 p = j ;
00062 temp = a[j] ;
00063 do{
00064 a[p] = a[p-1] ;
00065 p-- ;
00066 } while( p > 0 && temp < a[p-1] ) ;
00067 a[p] = temp ;
00068 }
00069 }
00070 return ;
00071 }
00072
00073
00074
00075 void qsrec_short( int n , short * ar , int cutoff )
00076 {
00077 register int i , j ;
00078 register short temp , pivot ;
00079 register short * a = ar ;
00080
00081 int left , right , mst ;
00082
00083
00084
00085 if( cutoff < 3 ) cutoff = 3 ;
00086 if( n < cutoff || ar == NULL ) return ;
00087
00088
00089
00090 stack[0] = 0 ;
00091 stack[1] = n-1 ;
00092 mst = 2 ;
00093
00094
00095
00096 while( mst > 0 ){
00097 right = stack[--mst] ;
00098 left = stack[--mst] ;
00099
00100 i = ( left + right ) / 2 ;
00101
00102
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] ;
00109 a[i] = a[right] ;
00110
00111 i = left ; j = right ;
00112
00113
00114
00115
00116 do{
00117 for( ; a[++i] < pivot ; ) ;
00118 for( ; a[--j] > pivot ; ) ;
00119
00120 if( j <= i ) break ;
00121
00122 QS_SWAP( a[i] , a[j] ) ;
00123 } while( 1 ) ;
00124
00125
00126
00127 a[right] = a[i] ; a[i] = pivot ;
00128
00129
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 }
00135 return ;
00136 }
00137
00138
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
00149
00150 void isort_int( int n , int * ar )
00151 {
00152 register int j , p ;
00153 register int temp ;
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] ){
00161 p = j ;
00162 temp = a[j] ;
00163 do{
00164 a[p] = a[p-1] ;
00165 p-- ;
00166 } while( p > 0 && temp < a[p-1] ) ;
00167 a[p] = temp ;
00168 }
00169 }
00170 return ;
00171 }
00172
00173
00174
00175 void qsrec_int( int n , int * ar , int cutoff )
00176 {
00177 register int i , j ;
00178 register int temp , pivot ;
00179 register int * a = ar ;
00180
00181 int left , right , mst ;
00182
00183
00184
00185 if( cutoff < 3 ) cutoff = 3 ;
00186 if( n < cutoff || ar == NULL ) return ;
00187
00188
00189
00190 stack[0] = 0 ;
00191 stack[1] = n-1 ;
00192 mst = 2 ;
00193
00194
00195
00196 while( mst > 0 ){
00197 right = stack[--mst] ;
00198 left = stack[--mst] ;
00199
00200 i = ( left + right ) / 2 ;
00201
00202
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] ;
00209 a[i] = a[right] ;
00210
00211 i = left ; j = right ;
00212
00213
00214
00215
00216 do{
00217 for( ; a[++i] < pivot ; ) ;
00218 for( ; a[--j] > pivot ; ) ;
00219
00220 if( j <= i ) break ;
00221
00222 QS_SWAP( a[i] , a[j] ) ;
00223 } while( 1 ) ;
00224
00225
00226
00227 a[right] = a[i] ; a[i] = pivot ;
00228
00229
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 }
00235 return ;
00236 }
00237
00238
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
00249
00250 void isort_float( int n , float * ar )
00251 {
00252 register int j , p ;
00253 register float temp ;
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] ){
00261 p = j ;
00262 temp = a[j] ;
00263 do{
00264 a[p] = a[p-1] ;
00265 p-- ;
00266 } while( p > 0 && temp < a[p-1] ) ;
00267 a[p] = temp ;
00268 }
00269 }
00270 return ;
00271 }
00272
00273
00274
00275 void qsrec_float( int n , float * ar , int cutoff )
00276 {
00277 register int i , j ;
00278 register float temp , pivot ;
00279 register float * a = ar ;
00280
00281 int left , right , mst ;
00282
00283
00284
00285 if( cutoff < 3 ) cutoff = 3 ;
00286 if( n < cutoff || ar == NULL ) return ;
00287
00288
00289
00290 stack[0] = 0 ;
00291 stack[1] = n-1 ;
00292 mst = 2 ;
00293
00294
00295
00296 while( mst > 0 ){
00297 right = stack[--mst] ;
00298 left = stack[--mst] ;
00299
00300 i = ( left + right ) / 2 ;
00301
00302
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] ;
00309 a[i] = a[right] ;
00310
00311 i = left ; j = right ;
00312
00313
00314
00315
00316 do{
00317 for( ; a[++i] < pivot ; ) ;
00318 for( ; a[--j] > pivot ; ) ;
00319
00320 if( j <= i ) break ;
00321
00322 QS_SWAP( a[i] , a[j] ) ;
00323 } while( 1 ) ;
00324
00325
00326
00327 a[right] = a[i] ; a[i] = pivot ;
00328
00329
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 }
00335 return ;
00336 }
00337
00338
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
00349
00350 void isort_pair( int n , float * ar , int * iar )
00351 {
00352 register int j , p ;
00353 register float temp ;
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] ){
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
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 ;
00382 register float temp , pivot ;
00383 register int itemp ,ipivot ;
00384 register float * a = ar ;
00385 register int *ia = iar ;
00386
00387 int left , right , mst ;
00388
00389
00390
00391 if( cutoff < 3 ) cutoff = 3 ;
00392 if( n < cutoff || ar == NULL || iar == NULL ) return ;
00393
00394
00395
00396 stack[0] = 0 ;
00397 stack[1] = n-1 ;
00398 mst = 2 ;
00399
00400
00401
00402 while( mst > 0 ){
00403 right = stack[--mst] ;
00404 left = stack[--mst] ;
00405
00406 i = ( left + right ) / 2 ;
00407
00408
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 ;
00421
00422
00423
00424
00425 do{
00426 for( ; a[++i] < pivot ; ) ;
00427 for( ; a[--j] > pivot ; ) ;
00428
00429 if( j <= i ) break ;
00430
00431 QS_SWAP ( a[i] , a[j] ) ;
00432 QS_ISWAP(ia[i] ,ia[j] ) ;
00433 } while( 1 ) ;
00434
00435
00436
00437 a[right] = a[i] ; a[i] = pivot ;
00438 ia[right]=ia[i] ;ia[i] =ipivot ;
00439
00440
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 }
00446 return ;
00447 }
00448
00449
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
00460
00461
00462
00463
00464
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
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
00486
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
00507
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
00533
00534
00535
00536
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
00554
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
00571
00572 flim = mri_to_float( im ) ; far = MRI_FLOAT_PTR(flim) ;
00573
00574
00575
00576
00577 qsort_pair( nvox , far , iar ) ;
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588 fac = 1.0 / nvox ; nvox1 = nvox - 1 ;
00589
00590 for( ibot=0 ; ibot < nvox1 ; ){
00591
00592
00593
00594 val = far[ibot] ; itop = ibot+1 ;
00595 if( val != far[itop] ){
00596 far[ibot] = fac * ibot ;
00597 ibot++ ; continue ;
00598 }
00599
00600
00601
00602 for( ; itop < nvox1 && val == far[itop] ; itop++ ) ;
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
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
00622
00623
00624
00625
00626
00627
00628 float mri_quantile( MRI_IMAGE * im , float alpha )
00629 {
00630 int ii , nvox ;
00631 float fi , quan ;
00632
00633
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
00649
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
00668
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 }