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 #include <math.h>
9
10 /*** 7D SAFE ***/
11
12 /*---------------------------------------------------------------------------*/
13
mri_max(MRI_IMAGE * im)14 double mri_max( MRI_IMAGE *im )
15 {
16 register int ii , npix ;
17 byte byte_max = 0 ;
18 short short_max = -32767 ; /* 23 Oct 1998: changed from 0 */
19 int int_max = -2147483647 ; /* ditto */
20 float float_max = -1.e+38 ; /* changed from -9999999.0 */
21 double double_max = -1.e+38 ; /* ditto */
22
23 ENTRY("mri_max") ;
24
25 npix = im->nvox ;
26
27 switch( im->kind ){
28
29 case MRI_byte:{
30 byte *qar = MRI_BYTE_PTR(im) ;
31 for( ii=0 ; ii < npix ; ii++ )
32 byte_max = MAX( byte_max , qar[ii] ) ;
33 RETURN( (double) byte_max );
34 }
35
36 case MRI_short:{
37 short *qar = MRI_SHORT_PTR(im) ;
38 for( ii=0 ; ii < npix ; ii++ )
39 short_max = MAX( short_max , qar[ii] ) ;
40 RETURN( (double) short_max );
41 }
42
43 case MRI_int:{
44 int *qar = MRI_INT_PTR(im) ;
45 for( ii=0 ; ii < npix ; ii++ )
46 int_max = MAX( int_max , qar[ii] ) ;
47 RETURN( (double) int_max );
48 }
49
50 case MRI_float:{
51 float *qar = MRI_FLOAT_PTR(im) ;
52 for( ii=0 ; ii < npix ; ii++ )
53 float_max = MAX( float_max , qar[ii] ) ;
54 RETURN( (double) float_max );
55 }
56
57 case MRI_double:{
58 double *qar = MRI_DOUBLE_PTR(im) ;
59 for( ii=0 ; ii < npix ; ii++ )
60 double_max = MAX( double_max , qar[ii] ) ;
61 RETURN( double_max );
62 }
63
64 case MRI_complex:{
65 complex *qar = MRI_COMPLEX_PTR(im) ;
66 for( ii=0 ; ii < npix ; ii++ )
67 float_max = MAX( float_max , CABS(qar[ii]) ) ;
68 RETURN( float_max );
69 }
70
71 case MRI_rgb:{
72 byte *rgb = MRI_RGB_PTR(im) ;
73 double val , top=0.0 ;
74 for( ii=0 ; ii < npix ; ii++ ){ /* scale to brightness */
75 val = 0.299 * rgb[3*ii] /* between 0 and 255 */
76 + 0.587 * rgb[3*ii+1]
77 + 0.114 * rgb[3*ii+2] ;
78 if( val > top ) top = val ;
79 }
80 RETURN( top );
81 }
82
83 default:
84 ERROR_message("mri_max: unknown image kind") ;
85 }
86 RETURN( 0.0 );
87 }
88
89 /*---------------------------------------------------------------------------*/
90
mri_maxabs(MRI_IMAGE * im)91 double mri_maxabs( MRI_IMAGE * im )
92 {
93 register int ii , npix ;
94 byte byte_max = 0 ;
95 int int_max = 0 ;
96 double double_max = 0.0 ;
97
98 ENTRY("mri_maxabs") ;
99
100 npix = im->nvox ;
101
102 switch( im->kind ){
103
104 case MRI_byte:{
105 byte *qar = MRI_BYTE_PTR(im) ;
106 for( ii=0 ; ii < npix ; ii++ )
107 byte_max = MAX( byte_max , qar[ii] ) ;
108 RETURN( (double) byte_max );
109 }
110
111 case MRI_short:{
112 short *qar = MRI_SHORT_PTR(im) ;
113 for( ii=0 ; ii < npix ; ii++ )
114 int_max = MAX( int_max , abs(qar[ii]) ) ;
115 RETURN( (double) int_max );
116 }
117
118 case MRI_int:{
119 int *qar = MRI_INT_PTR(im) ;
120 for( ii=0 ; ii < npix ; ii++ )
121 int_max = MAX( int_max , abs(qar[ii]) ) ;
122 RETURN( (double) int_max );
123 }
124
125 case MRI_float:{
126 float *qar = MRI_FLOAT_PTR(im) ;
127 for( ii=0 ; ii < npix ; ii++ )
128 double_max = MAX( double_max , fabs(qar[ii]) ) ;
129 RETURN( double_max );
130 }
131
132 case MRI_double:{
133 double *qar = MRI_DOUBLE_PTR(im) ;
134 for( ii=0 ; ii < npix ; ii++ )
135 double_max = MAX( double_max , fabs(qar[ii]) ) ;
136 RETURN( double_max );
137 }
138
139 case MRI_complex:{
140 complex *qar = MRI_COMPLEX_PTR(im) ;
141 for( ii=0 ; ii < npix ; ii++ )
142 double_max = MAX( double_max , CABS(qar[ii]) ) ;
143 RETURN( double_max );
144 }
145
146 case MRI_rgb:
147 RETURN( mri_max( im ) );
148
149 default:
150 ERROR_message("mri_max: unknown image kind") ;
151 }
152 RETURN( 0 );
153 }
154
155 /*---------------------------------------------------------------------------*/
156
mri_min(MRI_IMAGE * im)157 double mri_min( MRI_IMAGE *im )
158 {
159 register int ii , npix ;
160 byte byte_min = 255 ;
161 short short_min = 32767 ;
162 int int_min = 2147483647 ;
163 float float_min = 1.e+38 ;
164 double double_min = 1.e+38 ;
165
166 ENTRY("mri_min") ;
167
168 npix = im->nvox ;
169
170 switch( im->kind ){
171
172 case MRI_byte:{
173 byte *qar = MRI_BYTE_PTR(im) ;
174 for( ii=0 ; ii < npix ; ii++ )
175 byte_min = MIN( byte_min , qar[ii] ) ;
176 RETURN( (double) byte_min );
177 }
178
179 case MRI_short:{
180 short *qar = MRI_SHORT_PTR(im) ;
181 for( ii=0 ; ii < npix ; ii++ )
182 short_min = MIN( short_min , qar[ii] ) ;
183 RETURN( (double) short_min );
184 }
185
186 case MRI_int:{
187 int *qar = MRI_INT_PTR(im) ;
188 for( ii=0 ; ii < npix ; ii++ )
189 int_min = MIN( int_min , qar[ii] ) ;
190 RETURN( (double) int_min );
191 }
192
193 case MRI_float:{
194 float *qar = MRI_FLOAT_PTR(im) ;
195 for( ii=0 ; ii < npix ; ii++ )
196 float_min = MIN( float_min , qar[ii] ) ;
197 RETURN( (double) float_min );
198 }
199
200 case MRI_double:{
201 double *qar = MRI_DOUBLE_PTR(im) ;
202 for( ii=0 ; ii < npix ; ii++ )
203 double_min = MIN( double_min , qar[ii] ) ;
204 RETURN( double_min );
205 }
206
207 case MRI_complex:{
208 complex *qar = MRI_COMPLEX_PTR(im) ;
209 for( ii=0 ; ii < npix ; ii++ )
210 float_min = MIN( float_min , CABS(qar[ii]) ) ;
211 RETURN( float_min );
212 }
213
214 case MRI_rgb:{
215 byte *rgb = MRI_RGB_PTR(im) ;
216 double val , bot=255.9 ;
217 for( ii=0 ; ii < npix ; ii++ ){ /* scale to brightness */
218 val = 0.299 * rgb[3*ii] /* between 0 and 255 */
219 + 0.587 * rgb[3*ii+1]
220 + 0.114 * rgb[3*ii+2] ;
221 if( val < bot ) bot = val ;
222 }
223 RETURN( bot );
224 }
225
226 default:
227 ERROR_message("mri_min: unknown image kind") ;
228 }
229 RETURN( 0 );
230 }
231
232 /*---------------------------------------------------------------------------*/
233
mri_minmax(MRI_IMAGE * im)234 double_pair mri_minmax( MRI_IMAGE *im )
235 {
236 register int ii , npix ;
237 byte byte_min = 255 ;
238 short short_min = 32767 ;
239 int int_min = 2147483647 ;
240 float float_min = 1.e+38 ;
241 double double_min = 1.e+38 ;
242 byte byte_max = 0 ;
243 short short_max = -32767 ;
244 int int_max = -2147483647 ;
245 float float_max = -1.e+38 ;
246 double double_max = -1.e+38 ;
247 double_pair dp = {0.0,0.0} ;
248
249 ENTRY("mri_minmax") ;
250
251 npix = im->nvox ;
252
253 switch( im->kind ){
254
255 case MRI_byte:{
256 byte *qar = MRI_BYTE_PTR(im) ;
257 for( ii=0 ; ii < npix ; ii++ ){
258 byte_min = MIN( byte_min , qar[ii] ) ;
259 byte_max = MAX( byte_max , qar[ii] ) ;
260 }
261 dp.a = (double)byte_min ; dp.b = (double)byte_max ; RETURN(dp) ;
262 }
263
264 case MRI_short:{
265 short *qar = MRI_SHORT_PTR(im) ;
266 for( ii=0 ; ii < npix ; ii++ ){
267 short_min = MIN( short_min , qar[ii] ) ;
268 short_max = MAX( short_max , qar[ii] ) ;
269 }
270 dp.a = (double)short_min ; dp.b = (double)short_max ; RETURN(dp) ;
271 }
272
273 case MRI_float:{
274 float *qar = MRI_FLOAT_PTR(im) ;
275 for( ii=0 ; ii < npix ; ii++ ){
276 float_min = MIN( float_min , qar[ii] ) ;
277 float_max = MAX( float_max , qar[ii] ) ;
278 }
279 dp.a = (double)float_min ; dp.b = (double)float_max ; RETURN(dp) ;
280 }
281
282 default:
283 ERROR_message("mri_minmax: unknown image kind") ;
284 }
285 RETURN( dp );
286 }
287
288 /*---------------------------------------------------------------------------*/
289
mri_minmax_nz(MRI_IMAGE * im)290 double_pair mri_minmax_nz( MRI_IMAGE *im )
291 {
292 register int ii , npix ;
293 byte byte_min = 255 ;
294 short short_min = 32767 ;
295 int int_min = 2147483647 ;
296 float float_min = 1.e+38 ;
297 double double_min = 1.e+38 ;
298 byte byte_max = 0 ;
299 short short_max = -32767 ;
300 int int_max = -2147483647 ;
301 float float_max = -1.e+38 ;
302 double double_max = -1.e+38 ;
303 double_pair dp = {0.0,0.0} ;
304
305 ENTRY("mri_minmax_nz") ;
306
307 npix = im->nvox ;
308
309 switch( im->kind ){
310
311 case MRI_byte:{
312 byte *qar = MRI_BYTE_PTR(im) ;
313 for( ii=0 ; ii < npix ; ii++ ){
314 if( qar[ii] == 0 ) continue ;
315 byte_min = MIN( byte_min , qar[ii] ) ;
316 byte_max = MAX( byte_max , qar[ii] ) ;
317 }
318 if( byte_min <= byte_max ){
319 dp.a = (double)byte_min ; dp.b = (double)byte_max ;
320 }
321 RETURN(dp) ;
322 }
323
324 case MRI_short:{
325 short *qar = MRI_SHORT_PTR(im) ;
326 for( ii=0 ; ii < npix ; ii++ ){
327 if( qar[ii] == 0 ) continue ;
328 short_min = MIN( short_min , qar[ii] ) ;
329 short_max = MAX( short_max , qar[ii] ) ;
330 }
331 if( short_min <= short_max ){
332 dp.a = (double)short_min ; dp.b = (double)short_max ;
333 }
334 RETURN(dp) ;
335 }
336
337 case MRI_float:{
338 float *qar = MRI_FLOAT_PTR(im) ;
339 for( ii=0 ; ii < npix ; ii++ ){
340 if( qar[ii] == 0 ) continue ;
341 float_min = MIN( float_min , qar[ii] ) ;
342 float_max = MAX( float_max , qar[ii] ) ;
343 }
344 if( float_min <= float_max ){
345 dp.a = (double)float_min ; dp.b = (double)float_max ;
346 }
347 RETURN(dp) ;
348 }
349
350 case MRI_complex:{
351 complex *qar = MRI_COMPLEX_PTR(im) ; float aval ;
352 for( ii=0 ; ii < npix ; ii++ ){
353 aval = CABS(qar[ii]) ; if( aval == 0.0f ) continue ;
354 float_min = MIN( float_min , aval ) ;
355 float_max = MAX( float_max , aval ) ;
356 }
357 if( float_min <= float_max ){
358 dp.a = (double)float_min ; dp.b = (double)float_max ;
359 }
360 RETURN(dp) ;
361 }
362
363 case MRI_rgb:{
364 byte *rgb = MRI_RGB_PTR(im) ;
365 double val ;
366 for( ii=0 ; ii < npix ; ii++ ){ /* scale to brightness */
367 val = 0.299 * rgb[3*ii] /* between 0 and 255 */
368 + 0.587 * rgb[3*ii+1]
369 + 0.114 * rgb[3*ii+2] ; if( val == 0.0f ) continue ;
370 float_min = MIN( float_min , val ) ;
371 float_max = MAX( float_max , val ) ;
372 }
373 if( float_min <= float_max ){
374 dp.a = (double)float_min ; dp.b = (double)float_max ;
375 }
376 RETURN(dp) ;
377 }
378
379 default:
380 ERROR_message("mri_minmax_nz: unknown image kind") ;
381 }
382 RETURN( dp );
383 }
384
385 /*---------------------------------------------------------------------------*/
386
mri_indmax_nz(MRI_IMAGE * im)387 intfloat mri_indmax_nz( MRI_IMAGE *im )
388 {
389 register int ii , npix ;
390 byte byte_max = 0 ;
391 short short_max = -32767 ; /* 23 Oct 1998: changed from 0 */
392 int int_max = -2147483647 ; /* ditto */
393 float float_max = -1.e+38 ; /* changed from -9999999.0 */
394 double double_max = -1.e+38 ; /* ditto */
395 intfloat result = { -1 , 0.0f } ; int imax=-1 ;
396
397 ENTRY("mri_indmax_nz") ;
398
399 npix = im->nvox ;
400
401 switch( im->kind ){
402 default: break ;
403
404 case MRI_byte:{
405 byte *qar = MRI_BYTE_PTR(im) ;
406 for( ii=0 ; ii < npix ; ii++ ){
407 if( qar[ii] == 0 ) continue ;
408 if( qar[ii] > byte_max ){ byte_max = qar[ii]; imax = ii; }
409 }
410 result.i = imax; result.a = (float)byte_max; RETURN(result);
411 }
412
413 case MRI_short:{
414 short *qar = MRI_SHORT_PTR(im) ;
415 for( ii=0 ; ii < npix ; ii++ ){
416 if( qar[ii] == 0 ) continue ;
417 if( qar[ii] > short_max ){ short_max = qar[ii]; imax = ii; }
418 }
419 result.i = imax; result.a = (float)short_max; RETURN(result);
420 }
421
422 case MRI_int:{
423 int *qar = MRI_INT_PTR(im) ;
424 for( ii=0 ; ii < npix ; ii++ ){
425 if( qar[ii] == 0 ) continue ;
426 if( qar[ii] > int_max ){ int_max = qar[ii]; imax = ii; }
427 }
428 result.i = imax; result.a = (float)int_max; RETURN(result);
429 }
430
431 case MRI_float:{
432 float *qar = MRI_FLOAT_PTR(im) ;
433 for( ii=0 ; ii < npix ; ii++ ){
434 if( qar[ii] == 0.0f ) continue ;
435 if( qar[ii] > float_max ){ float_max = qar[ii]; imax = ii; }
436 }
437 result.i = imax; result.a = float_max; RETURN(result);
438 }
439
440 case MRI_double:{
441 double *qar = MRI_DOUBLE_PTR(im) ;
442 for( ii=0 ; ii < npix ; ii++ ){
443 if( qar[ii] == 0.0 ) continue ;
444 if( qar[ii] > double_max ){ double_max = qar[ii]; imax = ii; }
445 }
446 result.i = imax; result.a = (float)double_max; RETURN(result);
447 }
448
449 case MRI_complex:{
450 complex *qar = MRI_COMPLEX_PTR(im) ; float aval ;
451 for( ii=0 ; ii < npix ; ii++ ){
452 aval = CABS(qar[ii]) ;
453 if( aval == 0.0f ) continue ;
454 if( aval > float_max ){ float_max = aval; imax = ii; }
455 }
456 result.i = imax; result.a = float_max; RETURN(result);
457 }
458
459 case MRI_rgb:{
460 byte *rgb = MRI_RGB_PTR(im) ;
461 double val , top=0.0 ;
462 for( ii=0 ; ii < npix ; ii++ ){ /* scale to brightness */
463 val = 0.299 * rgb[3*ii] /* between 0 and 255 */
464 + 0.587 * rgb[3*ii+1]
465 + 0.114 * rgb[3*ii+2] ;
466 if( val != 0.0 && val > top ){ top = val; imax = ii; }
467 }
468 result.i = imax; result.a = (float)top; RETURN(result);
469 }
470
471 }
472 RETURN(result);
473 }
474
475 /*---------------------------------------------------------------------------*/
476
mri_indmin_nz(MRI_IMAGE * im)477 intfloat mri_indmin_nz( MRI_IMAGE *im )
478 {
479 register int ii , npix ;
480 byte byte_min = 255 ;
481 short short_min = 32767 ; /* 23 Oct 1998: changed from 0 */
482 int int_min = 2147483647 ; /* ditto */
483 float float_min = 1.e+38 ; /* changed from -9999999.0 */
484 double double_min = 1.e+38 ; /* ditto */
485 intfloat result = { -1 , 0.0f } ; int imin=-1 ;
486
487 ENTRY("mri_indmin_nz") ;
488
489 npix = im->nvox ;
490
491 switch( im->kind ){
492 default: break ;
493
494 case MRI_byte:{
495 byte *qar = MRI_BYTE_PTR(im) ;
496 for( ii=0 ; ii < npix ; ii++ ){
497 if( qar[ii] == 0 ) continue ;
498 if( qar[ii] < byte_min ){ byte_min = qar[ii]; imin = ii; }
499 }
500 result.i = imin; result.a = (float)byte_min; RETURN(result);
501 }
502
503 case MRI_short:{
504 short *qar = MRI_SHORT_PTR(im) ;
505 for( ii=0 ; ii < npix ; ii++ ){
506 if( qar[ii] == 0 ) continue ;
507 if( qar[ii] < short_min ){ short_min = qar[ii]; imin = ii; }
508 }
509 result.i = imin; result.a = (float)short_min; RETURN(result);
510 }
511
512 case MRI_int:{
513 int *qar = MRI_INT_PTR(im) ;
514 for( ii=0 ; ii < npix ; ii++ ){
515 if( qar[ii] == 0 ) continue ;
516 if( qar[ii] < int_min ){ int_min = qar[ii]; imin = ii; }
517 }
518 result.i = imin; result.a = (float)int_min; RETURN(result);
519 }
520
521 case MRI_float:{
522 float *qar = MRI_FLOAT_PTR(im) ;
523 for( ii=0 ; ii < npix ; ii++ ){
524 if( qar[ii] == 0.0f ) continue ;
525 if( qar[ii] < float_min ){ float_min = qar[ii]; imin = ii; }
526 }
527 result.i = imin; result.a = float_min; RETURN(result);
528 }
529
530 case MRI_double:{
531 double *qar = MRI_DOUBLE_PTR(im) ;
532 for( ii=0 ; ii < npix ; ii++ ){
533 if( qar[ii] == 0.0 ) continue ;
534 if( qar[ii] < double_min ){ double_min = qar[ii]; imin = ii; }
535 }
536 result.i = imin; result.a = (float)double_min; RETURN(result);
537 }
538
539 case MRI_complex:{
540 complex *qar = MRI_COMPLEX_PTR(im) ; float aval ;
541 for( ii=0 ; ii < npix ; ii++ ){
542 aval = CABS(qar[ii]) ;
543 if( aval == 0.0f ) continue ;
544 if( aval < float_min ){ float_min = aval; imin = ii; }
545 }
546 result.i = imin; result.a = float_min; RETURN(result);
547 }
548
549 case MRI_rgb:{
550 byte *rgb = MRI_RGB_PTR(im) ;
551 double val , bot=255.0 ;
552 for( ii=0 ; ii < npix ; ii++ ){ /* scale to brightness */
553 val = 0.299 * rgb[3*ii] /* between 0 and 255 */
554 + 0.587 * rgb[3*ii+1]
555 + 0.114 * rgb[3*ii+2] ;
556 if( val != 0.0 && val < bot ){ bot = val; imin = ii; }
557 }
558 result.i = imin; result.a = (float)bot; RETURN(result);
559 }
560
561 }
562 RETURN(result);
563 }
564