1 #include "mrilib.h"
2
3 /*---------------------------------------------------------------------------*/
4 /* Invert the contrast of an image, overwriting itself.
5 uperc = upper percentage point at which to invert
6 mask = what do you think?
7 Works on 2D or 3D images.
8 *//*-------------------------------------------------------------------------*/
9
mri_invertcontrast_inplace(MRI_IMAGE * im,float uperc,byte * mask)10 void mri_invertcontrast_inplace( MRI_IMAGE *im , float uperc , byte *mask )
11 {
12 byte *mmm=NULL ; MRI_IMAGE *qim ;
13 int nvox , nhist , ii ; float *hist=NULL , *imar , ucut ;
14
15 if( im == NULL || im->ny == 1 ) return ;
16
17 if( uperc < 90.0f ) uperc = 90.0f ;
18 else if( uperc > 100.0f ) uperc = 100.0f ;
19
20 /* make a float copy */
21
22 qim = mri_to_float(im) ;
23
24 /* find the mask */
25
26 mmm = mask ;
27 if( mmm == NULL ) mmm = mri_automask_image(qim) ;
28
29 /* extract nonzero voxels */
30
31 nvox = im->nvox ;
32 hist = (float *)malloc(sizeof(float)*nvox) ;
33 imar = MRI_FLOAT_PTR(qim) ;
34 for( nhist=ii=0 ; ii < nvox ; ii++ ){ if( mmm[ii] ) hist[nhist++] = imar[ii]; }
35
36 if( nhist < 100 ){
37 if( mmm != mask ) free(mmm) ;
38 free(hist) ; mri_free(qim) ; return ;
39 }
40
41 /* sort them to find the correct value at which to invert */
42
43 qsort_float(nhist,hist) ;
44 ii = (int)rintf(nhist*uperc*0.01f) ; ucut = hist[ii] ; free(hist) ;
45
46 /* invert contrast */
47
48 for( ii=0 ; ii < nvox ; ii++ ){
49 if( mmm[ii] && imar[ii] > 0.0f ) imar[ii] = ucut - imar[ii] ;
50 if( !mmm[ii] || imar[ii] < 0.0f ) imar[ii] = 0.0f ;
51 }
52 if( mmm != mask ) free(mmm) ;
53
54 /* apply the inverted contrast to the original image */
55
56 switch( im->kind ){
57
58 case MRI_float:{
59 AAmemcpy( MRI_FLOAT_PTR(im) , imar , sizeof(float)*nvox ) ;
60 }
61 break ;
62
63 case MRI_short:{
64 short *shar = MRI_SHORT_PTR(im) ;
65 for( ii=0 ; ii < nvox ; ii++ ) shar[ii] = SHORTIZE(imar[ii]) ;
66 }
67 break ;
68
69 #undef INTENSITY
70 #define INTENSITY(r,g,b) (0.299*(r)+0.587*(g)+0.114*(b))
71 case MRI_rgb:{
72 byte *bar = MRI_RGB_PTR(im) ; float r,g,b, irgb,iinv ;
73 for( ii=0 ; ii < nvox ; ii++ ){
74 iinv = imar[ii] ;
75 if( iinv == 0.0f ){
76 bar[3*ii+0] = bar[3*ii+1] = bar[3*ii+2] = 0 ;
77 } else {
78 r = bar[3*ii+0] ;
79 g = bar[3*ii+1] ;
80 b = bar[3*ii+2] ; irgb = 0.299*r + 0.587*g + 0.144*b ;
81 if( irgb > 0.0f ){ /* should always be true */
82 irgb = iinv / irgb ;
83 r *= irgb ; bar[3*ii+0] = BYTEIZE(r) ;
84 g *= irgb ; bar[3*ii+1] = BYTEIZE(g) ;
85 b *= irgb ; bar[3*ii+2] = BYTEIZE(b) ;
86 }
87 }
88 }
89 }
90 break ;
91
92 default:
93 WARNING_message("unknown image type %s in mri_invertcontrast_inplace()",
94 MRI_type_string(im->kind) ) ;
95 break ;
96
97 }
98
99 mri_free(qim) ;
100 return ;
101 }
102