1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8
9 /*** 7D SAFE ***/
10
11 /* prototypes for internal sorting routines:
12 array of shorts
13 array of ints
14 array of floats
15 array of floats, with an integer array carried along in the swaps */
16
17 void qsrec_short( int , short * , int ) ;
18 void qsrec_float( int , float * , int ) ;
19 void qsrec_pair ( int , float * , int * , int ) ;
20
21 static void qsrec_int ( int , int * , int ) ;
22
23 #undef QS_CUTOFF
24 #undef QS_SWAP
25 #undef QS_SWAPI
26 #undef QS_SMALL
27 #define QS_CUTOFF 40 /* cutoff to switch from qsort to isort */
28 #define QS_SMALL 21
29 #define QS_SWAP(x,y) (temp=(x), (x)=(y),(y)=temp)
30 #ifndef QS_STACK
31 # define QS_STACK 9999
32 #endif
33
34 /***************************************************************************
35 Each qsort_TYPE routine (TYPE=short, int, float, or pair) has two
36 pieces. The isort_TYPE routine does an insertion sort routine,
37 which is only fast for nearly sorted routines. The qsrec_TYPE
38 routine carries out the quicksort algorithm down to the level
39 QS_CUTOFF -- at that point, the array is nearly sorted.
40
41 In the qsrec_TYPE routines, the fundamental operation is
42 the SWAP of two items. In the isort_TYPE routines, the
43 value at the j index is determined to be out of place, so
44 it is stored in a temporary variable and the values below
45 it are moved up in the array until the proper place for
46 the former a[j] is found. Then it is stored where it belongs.
47 Compare the isort_short and isort_pair routines to see how
48 this should be generalized to more complex objects.
49
50 -- Robert W. Cox
51 ***************************************************************************/
52
53 /*------------------------------------------------------------------------------*/
54 /*------------- insertion sort : sort an array of short in-place ---------------*/
55
isort_short(int n,short * ar)56 void isort_short( int n , short *ar )
57 {
58 register int j , p ; /* array indices */
59 register short temp ; /* a[j] holding place */
60 register short *a = ar ;
61
62 if( n < 2 || ar == NULL ) return ;
63
64 for( j=1 ; j < n ; j++ ){
65
66 if( a[j] < a[j-1] ){ /* out of order */
67 p = j ;
68 temp = a[j] ;
69 do{
70 a[p] = a[p-1] ; /* at this point, a[p-1] > temp, so move it up */
71 p-- ;
72 } while( p > 0 && temp < a[p-1] ) ;
73 a[p] = temp ; /* finally, put temp in its place */
74 }
75 }
76 return ;
77 }
78
79 /*--------- qsrec : recursive part of quicksort (stack implementation) ----------*/
80
qsrec_short(int n,short * ar,int cutoff)81 void qsrec_short( int n , short *ar , int cutoff )
82 {
83 register int i , j ; /* scanning indices */
84 register short temp , pivot ; /* holding places */
85 register short *a = ar ;
86
87 int left , right , mst ;
88 int stack[QS_STACK] ; /* stack for qsort "recursion" */
89
90 /* return if too short (insertion sort will clean up) */
91
92 if( cutoff < 3 ) cutoff = 3 ;
93 if( n < cutoff || ar == NULL ) return ;
94
95 /* initialize stack to start with whole array */
96
97 stack[0] = 0 ;
98 stack[1] = n-1 ;
99 mst = 2 ;
100
101 /* loop while the stack is nonempty */
102
103 while( mst > 0 ){
104 right = stack[--mst] ; /* work on subarray from left -> right */
105 left = stack[--mst] ;
106
107 i = ( left + right ) / 2 ; /* middle of subarray */
108
109 /*----- sort the left, middle, and right a[]'s -----*/
110
111 if( a[left] > a[i] ) QS_SWAP( a[left] , a[i] ) ;
112 if( a[left] > a[right] ) QS_SWAP( a[left] , a[right] ) ;
113 if( a[i] > a[right] ) QS_SWAP( a[right] , a[i] ) ;
114
115 pivot = a[i] ; /* a[i] is the median-of-3 pivot! */
116 a[i] = a[right] ;
117
118 i = left ; j = right ; /* initialize scanning */
119
120 /*----- partition: move elements bigger than pivot up and elements
121 smaller than pivot down, scanning in from ends -----*/
122
123 do{
124 for( ; a[++i] < pivot ; ) ; /* scan i up, until a[i] >= pivot */
125 for( ; a[--j] > pivot ; ) ; /* scan j down, until a[j] <= pivot */
126
127 if( j <= i ) break ; /* if j meets i, quit */
128
129 QS_SWAP( a[i] , a[j] ) ;
130 } while( 1 ) ;
131
132 /*----- at this point, the array is partitioned -----*/
133
134 a[right] = a[i] ; a[i] = pivot ; /* restore the pivot */
135
136 /*----- signal the subarrays that need further work -----*/
137
138 if( (i-left) > cutoff ){ stack[mst++] = left ; stack[mst++] = i-1 ; }
139 if( (right-i) > cutoff ){ stack[mst++] = i+1 ; stack[mst++] = right ; }
140
141 } /* end of while stack is non-empty */
142 return ;
143 }
144
145 /* quick_sort : sort an array partially recursively, and partially insertion */
146
147 #undef DTYPE
148 #define DTYPE short
149 #include "cs_qsort_small.h"
150
qsort_short(int n,short * a)151 void qsort_short( int n , short *a )
152 {
153 if( n <= 1 ) return ;
154 switch(n){
155 default: break ; /* handled below */
156 case 2: qsort2_short(a) ; return ;
157 case 3: qsort3_short(a) ; return ;
158 case 4: qsort4_short(a) ; return ;
159 case 5: qsort5_short(a) ; return ;
160 case 6: qsort6_short(a) ; return ;
161 case 7: qsort7_short(a) ; return ;
162 case 8: qsort8_short(a) ; return ;
163 case 9: qsort9_short(a) ; return ;
164 case 10: qsort10_short(a); return ;
165 case 11: qsort11_short(a); return ;
166 case 12: qsort12_short(a); return ;
167 case 13: qsort13_short(a); return ;
168 case 14: qsort14_short(a); return ;
169 case 15: qsort15_short(a); return ;
170 case 16: qsort16_short(a); return ;
171 case 17: qsort17_short(a); return ;
172 case 18: qsort18_short(a); return ;
173 case 19: qsort19_short(a); return ;
174 case 20: qsort20_short(a); return ;
175 case 21: qsort21_short(a); return ;
176 case 25: qsort25_short(a); return ;
177 case 27: qsort27_short(a); return ;
178 }
179 qsrec_short( n , a , QS_CUTOFF ) ;
180 isort_short( n , a ) ;
181 return ;
182 }
183
184 /*----------------------------------------------------------------------------*/
185 /*------------- insertion sort : sort an array of int in-place ---------------*/
186
isort_int(int n,int * ar)187 static void isort_int( int n , int *ar )
188 {
189 register int j , p ; /* array indices */
190 register int temp ; /* a[j] holding place */
191 register int *a = ar ;
192
193 if( n < 2 || ar == NULL ) return ;
194
195 for( j=1 ; j < n ; j++ ){
196
197 if( a[j] < a[j-1] ){ /* out of order */
198 p = j ;
199 temp = a[j] ;
200 do{
201 a[p] = a[p-1] ; /* at this point, a[p-1] > temp, so move it up */
202 p-- ;
203 } while( p > 0 && temp < a[p-1] ) ;
204 a[p] = temp ; /* finally, put temp in its place */
205 }
206 }
207 return ;
208 }
209
210 /*--------- qsrec : recursive part of quicksort (stack implementation) ----------*/
211
qsrec_int(int n,int * ar,int cutoff)212 static void qsrec_int( int n , int *ar , int cutoff )
213 {
214 register int i , j ; /* scanning indices */
215 register int temp , pivot ; /* holding places */
216 register int *a = ar ;
217
218 int left , right , mst ;
219 int stack[QS_STACK] ; /* stack for qsort "recursion" */
220
221 /* return if too short (insertion sort will clean up) */
222
223 if( cutoff < 3 ) cutoff = 3 ;
224 if( n < cutoff || ar == NULL ) return ;
225
226 /* initialize stack to start with whole array */
227
228 stack[0] = 0 ;
229 stack[1] = n-1 ;
230 mst = 2 ;
231
232 /* loop while the stack is nonempty */
233
234 while( mst > 0 ){
235 right = stack[--mst] ; /* work on subarray from left -> right */
236 left = stack[--mst] ;
237
238 i = ( left + right ) / 2 ; /* middle of subarray */
239
240 /*----- sort the left, middle, and right a[]'s -----*/
241
242 if( a[left] > a[i] ) QS_SWAP( a[left] , a[i] ) ;
243 if( a[left] > a[right] ) QS_SWAP( a[left] , a[right] ) ;
244 if( a[i] > a[right] ) QS_SWAP( a[right] , a[i] ) ;
245
246 pivot = a[i] ; /* a[i] is the median-of-3 pivot! */
247 a[i] = a[right] ;
248
249 i = left ; j = right ; /* initialize scanning */
250
251 /*----- partition: move elements bigger than pivot up and elements
252 smaller than pivot down, scanning in from ends -----*/
253
254 do{
255 for( ; a[++i] < pivot ; ) ; /* scan i up, until a[i] >= pivot */
256 for( ; a[--j] > pivot ; ) ; /* scan j down, until a[j] <= pivot */
257
258 if( j <= i ) break ; /* if j meets i, quit */
259
260 QS_SWAP( a[i] , a[j] ) ;
261 } while( 1 ) ;
262
263 /*----- at this point, the array is partitioned -----*/
264
265 a[right] = a[i] ; a[i] = pivot ; /* restore the pivot */
266
267 /*----- signal the subarrays that need further work -----*/
268
269 if( (i-left) > cutoff ){ stack[mst++] = left ; stack[mst++] = i-1 ; }
270 if( (right-i) > cutoff ){ stack[mst++] = i+1 ; stack[mst++] = right ; }
271
272 } /* end of while stack is non-empty */
273 return ;
274 }
275
276 /* quick_sort: sort an array partially recursively, and partially insertion */
277
278 #undef DTYPE
279 #define DTYPE int
280 #include "cs_qsort_small.h"
281
qsort_int(int n,int * a)282 void qsort_int( int n , int *a )
283 {
284 if( a == NULL || n <= 1 ) return ;
285 switch(n){
286 default: break ; /* handled below */
287 case 2: qsort2_int(a) ; return ;
288 case 3: qsort3_int(a) ; return ;
289 case 4: qsort4_int(a) ; return ;
290 case 5: qsort5_int(a) ; return ;
291 case 6: qsort6_int(a) ; return ;
292 case 7: qsort7_int(a) ; return ;
293 case 8: qsort8_int(a) ; return ;
294 case 9: qsort9_int(a) ; return ;
295 case 10: qsort10_int(a); return ;
296 case 11: qsort11_int(a); return ;
297 case 12: qsort12_int(a); return ;
298 case 13: qsort13_int(a); return ;
299 case 14: qsort14_int(a); return ;
300 case 15: qsort15_int(a); return ;
301 case 16: qsort16_int(a); return ;
302 case 17: qsort17_int(a); return ;
303 case 18: qsort18_int(a); return ;
304 case 19: qsort19_int(a); return ;
305 case 20: qsort20_int(a); return ;
306 case 21: qsort21_int(a); return ;
307 case 25: qsort25_int(a); return ;
308 case 27: qsort27_int(a); return ;
309 }
310 qsrec_int( n , a , QS_CUTOFF ) ;
311 isort_int( n , a ) ;
312 return ;
313 }
314
315 /*------------------------------------------------------------------------------*/
316
qsort_int_mostly(int n,int * a,int cut)317 void qsort_int_mostly( int n , int *a , int cut ) /* 12 Aug 2017 */
318 {
319 if( cut < 3 || cut >= n || n < 28 ){
320 qsort_int(n,a) ; return ;
321 }
322 qsrec_int( n , a , cut ) ;
323 return ;
324 }
325
326 /*------------------------------------------------------------------------------*/
327 /*------------- insertion sort : sort an array of float in-place ---------------*/
328
isort_float(int n,float * ar)329 void isort_float( int n , float *ar )
330 {
331 register int j , p ; /* array indices */
332 register float temp ; /* a[j] holding place */
333 register float *a = ar ;
334
335 if( n < 2 || ar == NULL ) return ;
336
337 for( j=1 ; j < n ; j++ ){
338
339 if( a[j] < a[j-1] ){ /* out of order */
340 p = j ;
341 temp = a[j] ;
342 do{
343 a[p] = a[p-1] ; /* at this point, a[p-1] > temp, so move it up */
344 p-- ;
345 } while( p > 0 && temp < a[p-1] ) ;
346 a[p] = temp ; /* finally, put temp in its place */
347 }
348 }
349 return ;
350 }
351
352 /*--------- qsrec : recursive part of quicksort (stack implementation) ----------*/
353
qsrec_float(int n,float * ar,int cutoff)354 void qsrec_float( int n , float *ar , int cutoff )
355 {
356 register int i , j ; /* scanning indices */
357 register float temp , pivot ; /* holding places */
358 register float *a = ar ;
359
360 int left , right , mst ;
361 int stack[QS_STACK] ; /* stack for qsort "recursion" */
362
363 /* return if too short (insertion sort will clean up) */
364
365 if( cutoff < 3 ) cutoff = 3 ;
366 if( n < cutoff || ar == NULL ) return ;
367
368 /* initialize stack to start with whole array */
369
370 stack[0] = 0 ;
371 stack[1] = n-1 ;
372 mst = 2 ;
373
374 /* loop while the stack is nonempty */
375
376 while( mst > 0 ){
377 right = stack[--mst] ; /* work on subarray from left -> right */
378 left = stack[--mst] ;
379
380 i = ( left + right ) / 2 ; /* middle of subarray */
381
382 /*----- sort the left, middle, and right a[]'s -----*/
383
384 if( a[left] > a[i] ) QS_SWAP( a[left] , a[i] ) ;
385 if( a[left] > a[right] ) QS_SWAP( a[left] , a[right] ) ;
386 if( a[i] > a[right] ) QS_SWAP( a[right] , a[i] ) ;
387
388 pivot = a[i] ; /* a[i] is the median-of-3 pivot! */
389 a[i] = a[right] ;
390
391 i = left ; j = right ; /* initialize scanning */
392
393 /*----- partition: move elements bigger than pivot up and elements
394 smaller than pivot down, scanning in from ends -----*/
395
396 do{
397 for( ; a[++i] < pivot ; ) ; /* scan i up, until a[i] >= pivot */
398 for( ; a[--j] > pivot ; ) ; /* scan j down, until a[j] <= pivot */
399
400 if( j <= i ) break ; /* if j meets i, quit */
401
402 QS_SWAP( a[i] , a[j] ) ;
403 } while( 1 ) ;
404
405 /*----- at this point, the array is partitioned -----*/
406
407 a[right] = a[i] ; a[i] = pivot ; /* restore the pivot */
408
409 /*----- signal the subarrays that need further work -----*/
410
411 if( (i-left) > cutoff ){ stack[mst++] = left; stack[mst++] = i-1; }
412 else if( (i-left) <= QS_SMALL ){ qsort_float( (i-left) , a+left ) ; }
413
414 if( (right-i) > cutoff ){ stack[mst++] = i+1; stack[mst++] = right; }
415 else if( (right-i) <= QS_SMALL ){ qsort_float( (right-i) , a+(i+1) ) ; }
416
417 } /* end of while stack is non-empty */
418 return ;
419 }
420
421 /* quick_sort : sort an array partially recursively, and partially insertion */
422
423 #undef DTYPE
424 #define DTYPE float
425 #include "cs_qsort_small.h"
426
qsort_float(int n,float * a)427 void qsort_float( int n , float *a )
428 {
429 if( n <= 1 ) return ;
430 switch(n){
431 default: break ; /* handled below */
432 case 2: qsort2_float(a) ; return ;
433 case 3: qsort3_float(a) ; return ;
434 case 4: qsort4_float(a) ; return ;
435 case 5: qsort5_float(a) ; return ;
436 case 6: qsort6_float(a) ; return ;
437 case 7: qsort7_float(a) ; return ;
438 case 8: qsort8_float(a) ; return ;
439 case 9: qsort9_float(a) ; return ;
440 case 10: qsort10_float(a); return ;
441 case 11: qsort11_float(a); return ;
442 case 12: qsort12_float(a); return ;
443 case 13: qsort13_float(a); return ;
444 case 14: qsort14_float(a); return ;
445 case 15: qsort15_float(a); return ;
446 case 16: qsort16_float(a); return ;
447 case 17: qsort17_float(a); return ;
448 case 18: qsort18_float(a); return ;
449 case 19: qsort19_float(a); return ;
450 case 20: qsort20_float(a); return ;
451 case 21: qsort21_float(a); return ;
452 case 25: qsort25_float(a); return ;
453 case 27: qsort27_float(a); return ;
454 }
455 qsrec_float( n , a , QS_CUTOFF ) ;
456 isort_float( n , a ) ;
457 return ;
458 }
459
qsort_float_rev(int n,float * a)460 void qsort_float_rev( int n , float *a )
461 {
462 register int ii ;
463 if( n < 2 || a == NULL ) return ;
464 for( ii=0 ; ii < n ; ii++ ) a[ii] = -a[ii] ;
465 qsort_float(n,a) ;
466 for( ii=0 ; ii < n ; ii++ ) a[ii] = -a[ii] ;
467 return ;
468 }
469
470 /*------------------------------------------------------------------------------*/
471 /*--------------- insertion sort of a float-int pair of arrays -----------------*/
472 /* IN THIS CASE: the floats determine the order */
473
isort_pair(int n,float * ar,int * iar)474 void isort_pair( int n , float *ar , int *iar )
475 {
476 register int j , p ; /* array indices */
477 register float temp ; /* a[j] holding place */
478 register int itemp ;
479 register float * a = ar ;
480 register int *ia = iar ;
481
482 if( n < 2 || ar == NULL || iar == NULL ) return ;
483
484 for( j=1 ; j < n ; j++ ){
485
486 if( a[j] < a[j-1] ){ /* out of order */
487 p = j ;
488 temp = a[j] ; itemp = ia[j] ;
489 do{
490 a[p] = a[p-1] ; ia[p] = ia[p-1] ;
491 p-- ;
492 } while( p > 0 && temp < a[p-1] ) ;
493 a[p] = temp ; ia[p] = itemp ;
494 }
495 }
496 return ;
497 }
498
499 /*--------- qsrec : recursive part of quicksort (stack implementation) ----------*/
500
501 #define QS_ISWAP(x,y) (itemp=(x),(x)=(y),(y)=itemp)
502
qsrec_pair(int n,float * ar,int * iar,int cutoff)503 void qsrec_pair( int n , float *ar , int *iar , int cutoff )
504 {
505 register int i , j ; /* scanning indices */
506 register float temp , pivot ; /* holding places */
507 register int itemp ,ipivot ;
508 register float *a = ar ;
509 register int *ia = iar ;
510
511 int left , right , mst ;
512 int stack[QS_STACK] ; /* stack for qsort "recursion" */
513
514 /* return if too short (insertion sort will clean up) */
515
516 if( cutoff < 3 ) cutoff = 3 ;
517 if( n < cutoff || ar == NULL || iar == NULL ) return ;
518
519 /* initialize stack to start with whole array */
520
521 stack[0] = 0 ;
522 stack[1] = n-1 ;
523 mst = 2 ;
524
525 /* loop while the stack is nonempty */
526
527 while( mst > 0 ){
528 right = stack[--mst] ; /* work on subarray from left -> right */
529 left = stack[--mst] ;
530
531 i = ( left + right ) / 2 ; /* middle of subarray */
532
533 /*----- sort the left, middle, and right a[]'s -----*/
534
535 if( a[left] > a[i] ){ QS_SWAP ( a[left] , a[i] ) ;
536 QS_ISWAP(ia[left] ,ia[i] ) ; }
537 if( a[left] > a[right] ){ QS_SWAP ( a[left] , a[right] ) ;
538 QS_ISWAP(ia[left] ,ia[right] ) ; }
539 if( a[i] > a[right] ){ QS_SWAP ( a[right] , a[i] ) ;
540 QS_ISWAP(ia[right] ,ia[i] ) ; }
541
542 pivot = a[i] ; a[i] = a[right] ;
543 ipivot =ia[i] ; ia[i] =ia[right] ;
544
545 i = left ; j = right ; /* initialize scanning */
546
547 /*----- partition: move elements bigger than pivot up and elements
548 smaller than pivot down, scanning in from ends -----*/
549
550 do{
551 for( ; a[++i] < pivot ; ) ; /* scan i up, until a[i] >= pivot */
552 for( ; a[--j] > pivot ; ) ; /* scan j down, until a[j] <= pivot */
553
554 if( j <= i ) break ; /* if j meets i, quit */
555
556 QS_SWAP ( a[i] , a[j] ) ;
557 QS_ISWAP(ia[i] ,ia[j] ) ;
558 } while( 1 ) ;
559
560 /*----- at this point, the array is partitioned -----*/
561
562 a[right] = a[i] ; a[i] = pivot ; /* restore the pivot */
563 ia[right]=ia[i] ;ia[i] =ipivot ;
564
565 /*----- signal the subarrays that need further work -----*/
566
567 if( (i-left) > cutoff ){ stack[mst++] = left ; stack[mst++] = i-1 ; }
568 if( (right-i) > cutoff ){ stack[mst++] = i+1 ; stack[mst++] = right ; }
569
570 } /* end of while stack is non-empty */
571 return ;
572 }
573
574 /* quick_sort : sort an array partially recursively, and partially insertion */
575
qsort_pair(int n,float * a,int * ia)576 void qsort_pair( int n , float *a , int *ia )
577 {
578 qsrec_pair( n , a , ia , QS_CUTOFF ) ;
579 isort_pair( n , a , ia ) ;
580 return ;
581 }
582
qsort_pair_rev(int n,float * a,int * ia)583 void qsort_pair_rev( int n , float *a , int *ia ) /* 08 Mar 2019 */
584 {
585 register int ii ;
586 if( n < 2 || a == NULL ) return ;
587 for( ii=0 ; ii < n ; ii++ ) a[ii] = -a[ii] ;
588 qsort_pair(n,a,ia) ;
589 for( ii=0 ; ii < n ; ii++ ) a[ii] = -a[ii] ;
590 return ;
591 }
592
593 /*------------------------------------------------------------------------------*/
594 /*--------------- insertion sort of a float-int pair of arrays -----------------*/
595 /* IN THIS CASE: the ints determine the order [08 Mar 2019] */
596
isort_pairX(int n,float * ar,int * iar)597 void isort_pairX( int n , float *ar , int *iar )
598 {
599 register int j , p ; /* array indices */
600 register float temp ; /* a[j] holding place */
601 register int itemp ;
602 register float * a = ar ;
603 register int *ia = iar ;
604
605 if( n < 2 || ar == NULL || iar == NULL ) return ;
606
607 for( j=1 ; j < n ; j++ ){
608
609 if( ia[j] < ia[j-1] ){ /* out of order */
610 p = j ;
611 temp = a[j] ; itemp = ia[j] ;
612 do{
613 a[p] = a[p-1] ; ia[p] = ia[p-1] ;
614 p-- ;
615 } while( p > 0 && itemp < ia[p-1] ) ;
616 a[p] = temp ; ia[p] = itemp ;
617 }
618 }
619 return ;
620 }
621
622 /*--------- qsrec : recursive part of quicksort (stack implementation) ----------*/
623
qsrec_pairX(int n,float * ar,int * iar,int cutoff)624 void qsrec_pairX( int n , float *ar , int *iar , int cutoff )
625 {
626 register int i , j ; /* scanning indices */
627 register float temp , pivot ; /* holding places */
628 register int itemp ,ipivot ;
629 register float *a = ar ;
630 register int *ia = iar ;
631
632 int left , right , mst ;
633 int stack[QS_STACK] ; /* stack for qsort "recursion" */
634
635 /* return if too short (insertion sort will clean up) */
636
637 if( cutoff < 3 ) cutoff = 3 ;
638 if( n < cutoff || ar == NULL || iar == NULL ) return ;
639
640 /* initialize stack to start with whole array */
641
642 stack[0] = 0 ;
643 stack[1] = n-1 ;
644 mst = 2 ;
645
646 /* loop while the stack is nonempty */
647
648 while( mst > 0 ){
649 right = stack[--mst] ; /* work on subarray from left -> right */
650 left = stack[--mst] ;
651
652 i = ( left + right ) / 2 ; /* middle of subarray */
653
654 /*----- sort the left, middle, and right a[]'s -----*/
655
656 if( ia[left] > ia[i] ){ QS_SWAP ( a[left] , a[i] ) ;
657 QS_ISWAP(ia[left] ,ia[i] ) ; }
658 if( ia[left] > ia[right] ){ QS_SWAP ( a[left] , a[right] ) ;
659 QS_ISWAP(ia[left] ,ia[right] ) ; }
660 if( ia[i] > ia[right] ){ QS_SWAP ( a[right] , a[i] ) ;
661 QS_ISWAP(ia[right] ,ia[i] ) ; }
662
663 pivot = a[i] ; a[i] = a[right] ;
664 ipivot =ia[i] ; ia[i] =ia[right] ;
665
666 i = left ; j = right ; /* initialize scanning */
667
668 /*----- partition: move elements bigger than pivot up and elements
669 smaller than pivot down, scanning in from ends -----*/
670
671 do{
672 for( ; ia[++i] < ipivot ; ) ; /* scan i up, until a[i] >= pivot */
673 for( ; ia[--j] > ipivot ; ) ; /* scan j down, until a[j] <= pivot */
674
675 if( j <= i ) break ; /* if j meets i, quit */
676
677 QS_SWAP ( a[i] , a[j] ) ;
678 QS_ISWAP(ia[i] ,ia[j] ) ;
679 } while( 1 ) ;
680
681 /*----- at this point, the array is partitioned -----*/
682
683 a[right] = a[i] ; a[i] = pivot ; /* restore the pivot */
684 ia[right]=ia[i] ;ia[i] =ipivot ;
685
686 /*----- signal the subarrays that need further work -----*/
687
688 if( (i-left) > cutoff ){ stack[mst++] = left ; stack[mst++] = i-1 ; }
689 if( (right-i) > cutoff ){ stack[mst++] = i+1 ; stack[mst++] = right ; }
690
691 } /* end of while stack is non-empty */
692 return ;
693 }
694
695 /* quick_sort : sort an array partially recursively, and partially insertion */
696
qsort_pairX(int n,float * a,int * ia)697 void qsort_pairX( int n , float *a , int *ia )
698 {
699 qsrec_pairX( n , a , ia , QS_CUTOFF ) ;
700 isort_pairX( n , a , ia ) ;
701 return ;
702 }
703
704 /*******************************************************************************
705 Compute the "percentage points" of the histogram of an image.
706 "per" should be a pre-allocated array of "nper+1" floats; per[i] will be set
707 to the image intensity below which fraction i/nper of the pixels fall,
708 for i=0..nper.
709 N.B.: per[0] = image minimum and per[nper] = maximum.
710 "per" is float, no matter what the input image type is.
711 ********************************************************************************/
712
mri_percents(MRI_IMAGE * im,int nper,float per[])713 void mri_percents( MRI_IMAGE *im , int nper , float per[] )
714 {
715 register int pp , ii , nvox ;
716 register float fi , frac ;
717
718 /*** sanity checks ***/
719
720 if( im == NULL || per == NULL || nper < 2 ) return ;
721
722 nvox = im->nvox ;
723 frac = nvox / ((float) nper) ;
724
725 switch( im->kind ){
726
727 /*** create a float image copy of the data,
728 sort it, then interpolate the percentage points ***/
729
730 default:{
731 MRI_IMAGE *inim ;
732 float *far ;
733
734 inim = mri_to_float( im ) ;
735 far = MRI_FLOAT_PTR(inim) ;
736 qsort_float( nvox , far ) ;
737
738 per[0] = far[0] ;
739 for( pp=1 ; pp < nper ; pp++ ){
740 fi = frac *pp ; ii = fi ; fi = fi - ii ;
741 per[pp] = (1.0-fi) * far[ii] + fi * far[ii+1] ;
742 }
743 per[nper] = far[nvox-1] ;
744 mri_free( inim ) ;
745 }
746 break ;
747
748 /*** create a short image copy of the data,
749 sort it, then interpolate the percentage points ***/
750
751 case MRI_short:
752 case MRI_byte:{
753 MRI_IMAGE *inim ;
754 short *sar ;
755
756 inim = mri_to_short( 1.0 , im ) ;
757 sar = MRI_SHORT_PTR(inim) ;
758 qsort_short( nvox , sar ) ;
759
760 per[0] = sar[0] ;
761 for( pp=1 ; pp < nper ; pp++ ){
762 fi = frac *pp ; ii = fi ; fi = fi - ii ;
763 per[pp] = (1.0-fi) * sar[ii] + fi * sar[ii+1] ;
764 }
765 per[nper] = sar[nvox-1] ;
766 mri_free( inim ) ;
767 }
768 }
769
770 return ;
771 }
772
773 /******************************************************************************/
774 /* Get one percentile point of nonzero absolute values only. [24 May 2019] */
775 /******************************************************************************/
776
percentile_nzabs(int nvox,float * far,float per)777 float percentile_nzabs( int nvox , float *far , float per )
778 {
779 int pp , ii ;
780 float fi , frac , val ;
781 int nqq ; float *qar ;
782
783 if( nvox < 9 || far == NULL ) return 0.0f ;
784
785 if( per < 0.0f ) per = 0.0f ; /* min */
786 else if( per > 1.0f ) per *= 0.01f ; /* 37% -> 0.37 */
787 if( per > 1.0f ) per = 1.0f ; /* max */
788
789 qar = (float *)malloc(sizeof(float)*(nvox+1)) ;
790 for( nqq=ii=0 ; ii < nvox ; ii++ ){
791 val = fabsf(far[ii]) ; if( val == 0.0f ) continue ;
792 qar[nqq++] = val ;
793 }
794
795 /* special cases of very small values of nonzero voxels */
796
797 if( nqq == 0 ){ free(qar) ; return 0.0f ; }
798 if( nqq == 1 ){ val = qar[0] ; free(qar) ; return val ; }
799 if( nqq == 2 ){
800 val = qar[ (per < 0.5f) ? 0 : 1 ] ;
801 free(qar) ; return val ;
802 }
803
804 /* general case */
805
806 qsort_float( nqq , qar ) ;
807
808 if( per <= 0.0f ){ val = qar[0] ; free(qar) ; return val ; }
809
810 qar[nqq] = qar[nqq-1] ;
811 fi = nqq * per ; ii = (int)fi ; fi = fi - ii ;
812 val = (1.0f-fi) * qar[ii] + fi * qar[ii+1] ;
813
814 #if 0
815 INFO_message("percentile_nzabs: nvox=%d nqq=%d per=%.4f ii=%d fi=%.4f val=%g",nvox,nqq,per,ii,fi,val ) ;
816 #endif
817
818 free(qar) ; return val ;
819 }
820
821 /****************************************************************************
822 Produce a "histogram flattened" version of the input image. That is, the
823 output value at each pixel will be proportional to its place in the
824 histogram of intensities. That is, if the value at a pixel is the 173rd
825 in intensity (0th is darkest) out of 1000 pixels total, then the flattened
826 image value will be 0.173. The output image is in MRI_float format.
827 *****************************************************************************/
828
829 /* 'softening' of the flattening thru a parameter 'bfac' between 0 and 1.
830 bfac = amount of histogram that is pure flat (1 = all of it);
831 afac = quadratic part of histogram = 6*(1-bfac) */
832
833 #define USE_BFAC
834
835 #ifdef USE_BFAC
836 static float afac=6.0f ;
837 static float bfac=0.0f ;
mri_flatten_set_bfac(float b)838 void mri_flatten_set_bfac(float b)
839 {
840 bfac = ( b >= 0.0f && b <= 1.0f ) ? b : 1.0f ;
841 afac = 6.0f * (1.0f-bfac) ;
842 }
843 # define FLATV(x) ( (x)*( bfac + afac*(x)*( 0.5f - 0.3333333f*(x) ) ) )
844 #else
845 # define FLATV(x) (x)
846 #endif
847
mri_flatten(MRI_IMAGE * im)848 MRI_IMAGE * mri_flatten( MRI_IMAGE *im )
849 {
850 MRI_IMAGE *flim , *intim , *outim ;
851 float *far , *outar ;
852 int *iar ;
853 int ii , nvox , ibot,itop , nvox1,isub ;
854 float fac , val , xval ;
855
856 ENTRY("mri_flatten") ;
857
858 if( im == NULL ) RETURN(NULL) ;
859
860 /*** make an image that is just the voxel index in its array ***/
861 /*** also, make the output image while we are at it ***/
862
863 nvox = im->nvox ;
864 intim = mri_new_conforming( im , MRI_int ) ;
865 outim = mri_new_conforming( im , MRI_float ) ;
866
867 iar = MRI_INT_PTR(intim) ; outar = MRI_FLOAT_PTR(outim) ;
868
869 for( ii=0 ; ii < nvox ; ii++ ) iar[ii] = ii ;
870
871 /*** copy the input data to a floating point image ***/
872
873 flim = mri_to_float( im ) ; far = MRI_FLOAT_PTR(flim) ;
874
875 /*** sort this image, with the index array being carried along
876 so that we know where every pixel came from originally ***/
877
878 qsort_pair( nvox , far , iar ) ;
879
880 /*** The "far" array is now sorted. Thus, if the pixel that was in
881 voxel i is now in voxel j, then its place in the histogram is
882 j/nvox. The only difficulty is that there may be ties. We need
883 to resolve these ties so that pixels with the same intensity
884 don't get different output values. We do this by scanning
885 through far, finding blocks of equal values, and replacing
886 them by their average position in the histogram.
887 ***/
888
889 nvox1 = nvox - 1 ;
890 if( far[0] != 0.0f ){
891 fac = 1.0f/nvox ; ibot = 0 ; isub = 0 ;
892 } else {
893 for( ibot=1 ; ibot < nvox1 && far[ibot]==0.0f ; ibot++ ) ; /*nada*/
894 if( ibot == nvox1 ){
895 mri_free(flim) ; mri_free(intim) ; RETURN(NULL) ;
896 }
897 fac = 1.0f/(nvox-ibot-1) ;
898 isub = ibot-1 ; fac = 1.0f/(nvox-ibot) ;
899 }
900
901 for( ; ibot < nvox1 ; ){
902
903 /** if this value is unique, just set the value and move on **/
904
905 val = far[ibot] ; itop = ibot+1 ;
906 if( val != far[itop] ){
907 xval = fac*(ibot-isub) ; far[ibot] = FLATV(xval) ;
908 ibot++ ; continue ;
909 }
910
911 /** scan itop up until value is distinct **/
912
913 for( ; itop < nvox1 && val == far[itop] ; itop++ ) ; /* nada */
914
915 val = 0.5f*fac*(ibot+itop-1-2*isub) ; xval = FLATV(val) ;
916 for( ii=ibot ; ii < itop ; ii++ ) far[ii] = xval ;
917 ibot = itop ;
918 }
919 far[nvox1] = 1.0f ;
920
921 /*** now propagate these values back to the output image ***/
922
923 for( ii=0 ; ii < nvox ; ii++ ) outar[iar[ii]] = far[ii] ;
924
925 mri_free(flim) ; mri_free(intim) ;
926
927 MRI_COPY_AUX( outim , im ) ;
928 RETURN(outim) ;
929 }
930
931 /*-------------------------------------------------------------------*/
932 /*! Find the intensity in an image that is at the alpha-th
933 quantile of the distribution. That is, for 0 <= alpha <= 1,
934 alpha*npix of the image values are below, and (1-alpha)*npix
935 are above. If alpha is 0, this is the minimum. If alpha
936 is 1, this is the maximum.
937 ---------------------------------------------------------------------*/
938
mri_quantile(MRI_IMAGE * im,float alpha)939 float mri_quantile( MRI_IMAGE *im , float alpha )
940 {
941 int ii , nvox ;
942 float fi , quan ;
943
944 ENTRY("mri_quantile") ;
945
946 /*** sanity checks ***/
947
948 if( im == NULL ) RETURN( 0.0 );
949
950 if( alpha <= 0.0 ) RETURN( (float) mri_min(im) );
951 if( alpha >= 1.0 ) RETURN( (float) mri_max(im) );
952
953 nvox = im->nvox ;
954
955 switch( im->kind ){
956
957 /*** create a float image copy of the data,
958 sort it, then interpolate the percentage points ***/
959
960 default:{
961 MRI_IMAGE *inim ;
962 float *far ;
963
964 inim = mri_to_float( im ) ;
965 far = MRI_FLOAT_PTR(inim) ;
966 qsort_float( nvox , far ) ;
967
968 fi = alpha * nvox ;
969 ii = (int) fi ; if( ii >= nvox ) ii = nvox-1 ;
970 fi = fi - ii ;
971 quan = (1.0-fi) * far[ii] + fi * far[ii+1] ;
972 mri_free( inim ) ;
973 }
974 break ;
975
976 /*** create a short image copy of the data,
977 sort it, then interpolate the percentage points ***/
978
979 case MRI_short:
980 case MRI_byte:{
981 MRI_IMAGE *inim ;
982 short *sar ;
983
984 inim = mri_to_short( 1.0 , im ) ;
985 sar = MRI_SHORT_PTR(inim) ;
986 qsort_short( nvox , sar ) ;
987
988 fi = alpha * nvox ;
989 ii = (int) fi ; if( ii >= nvox ) ii = nvox-1 ;
990 fi = fi - ii ;
991 quan = (1.0-fi) * sar[ii] + fi * sar[ii+1] ;
992 mri_free( inim ) ;
993 }
994 break ;
995 }
996
997 RETURN( quan );
998 }
999
1000 /*-------------------------------------------------------------------*/
1001 /*! Return TWO quantile levels at once; cf. mri_quantile().
1002 ---------------------------------------------------------------------*/
1003
mri_twoquantiles(MRI_IMAGE * im,float alpha,float beta)1004 float_pair mri_twoquantiles( MRI_IMAGE *im, float alpha, float beta )
1005 {
1006 int ii , nvox ;
1007 float fi ;
1008 float_pair qt = {0.0f,0.0f} ;
1009 float qalph=WAY_BIG,qbeta=WAY_BIG ;
1010
1011 ENTRY("mri_twoquantiles") ;
1012
1013 /*** sanity checks ***/
1014
1015 if( im == NULL ) RETURN( qt );
1016
1017 if( alpha == beta ){
1018 qt.a = qt.b = mri_quantile(im,alpha) ; RETURN( qt );
1019 }
1020
1021 if( alpha <= 0.0f ) qalph = (float) mri_min(im) ;
1022 else if( alpha >= 1.0f ) qalph = (float) mri_max(im) ;
1023 if( beta <= 0.0f ) qbeta = (float) mri_min(im) ;
1024 else if( beta >= 1.0f ) qbeta = (float) mri_max(im) ;
1025
1026 if( qalph != WAY_BIG && qbeta != WAY_BIG ){
1027 qt.a = qalph; qt.b = qbeta; RETURN(qt);
1028 }
1029
1030 nvox = im->nvox ;
1031
1032 switch( im->kind ){
1033
1034 /*** create a float image copy of the data,
1035 sort it, then interpolate the percentage points ***/
1036
1037 default:{
1038 MRI_IMAGE *inim ;
1039 float *far ;
1040
1041 inim = mri_to_float( im ) ;
1042 far = MRI_FLOAT_PTR(inim) ;
1043 qsort_float( nvox , far ) ;
1044
1045 if( alpha > 0.0f && alpha < 1.0f ){
1046 fi = alpha * nvox ;
1047 ii = (int) fi ; if( ii >= nvox ) ii = nvox-1 ;
1048 fi = fi - ii ;
1049 qalph = (1.0-fi) * far[ii] + fi * far[ii+1] ;
1050 }
1051 if( beta > 0.0f && beta < 1.0f ){
1052 fi = beta * nvox ;
1053 ii = (int) fi ; if( ii >= nvox ) ii = nvox-1 ;
1054 fi = fi - ii ;
1055 qbeta = (1.0-fi) * far[ii] + fi * far[ii+1] ;
1056 }
1057 mri_free( inim ) ;
1058 }
1059 break ;
1060
1061 /*** create a short image copy of the data,
1062 sort it, then interpolate the percentage points ***/
1063
1064 case MRI_short:
1065 case MRI_byte:{
1066 MRI_IMAGE *inim ;
1067 short *sar ;
1068
1069 inim = mri_to_short( 1.0 , im ) ;
1070 sar = MRI_SHORT_PTR(inim) ;
1071 qsort_short( nvox , sar ) ;
1072
1073 if( alpha > 0.0f && alpha < 1.0f ){
1074 fi = alpha * nvox ;
1075 ii = (int) fi ; if( ii >= nvox ) ii = nvox-1 ;
1076 fi = fi - ii ;
1077 qalph = (1.0-fi) * sar[ii] + fi * sar[ii+1] ;
1078 }
1079 if( beta > 0.0f && beta < 1.0f ){
1080 fi = beta * nvox ;
1081 ii = (int) fi ; if( ii >= nvox ) ii = nvox-1 ;
1082 fi = fi - ii ;
1083 qbeta = (1.0-fi) * sar[ii] + fi * sar[ii+1] ;
1084 }
1085 mri_free( inim ) ;
1086 }
1087 break ;
1088 }
1089
1090 qt.a = qalph; qt.b = qbeta; RETURN(qt);
1091 }
1092