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 /*-------------------------------------
12 mode = 0 for normal multiplication
13 = 1 for conjugation of g
14 ---------------------------------------*/
15
mri_multiply_complex(int mode,MRI_IMAGE * f,MRI_IMAGE * g)16 MRI_IMAGE *mri_multiply_complex( int mode , MRI_IMAGE *f , MRI_IMAGE* g )
17 {
18 register int ii , npix ;
19 MRI_IMAGE *newImg ;
20 complex *nar , *gar , *far ;
21
22 if( f->nvox != g->nvox ){
23 fprintf( stderr , "mri_multiply_complex shapes imcompatible!\n" ) ;
24 MRI_FATAL_ERROR ;
25 }
26
27 if( f->kind != MRI_complex || g->kind != MRI_complex ){
28 fprintf( stderr , "mri_multiply_complex illegal image type!\n" ) ;
29 MRI_FATAL_ERROR ;
30 }
31
32 newImg = mri_new_conforming( f , MRI_complex ) ;
33 npix = f->nvox ;
34 MRI_COPY_AUX( newImg , f ) ;
35 far = MRI_COMPLEX_PTR(f); gar = MRI_COMPLEX_PTR(g); nar = MRI_COMPLEX_PTR(newImg);
36
37 switch( mode ){
38 case 0:
39 for( ii=0 ; ii < npix ; ii++ ){
40 nar[ii] = CMULT( far[ii] , gar[ii] ) ;
41 }
42 break ;
43
44 case 1:
45 for( ii=0 ; ii < npix ; ii++ ){
46 nar[ii] = CJMULT( far[ii] , gar[ii] ) ;
47 }
48 break ;
49
50 default:
51 fprintf( stderr , "mri_multiply_complex illegal mode %d\n" , mode ) ;
52 MRI_FATAL_ERROR ;
53 }
54 return newImg ;
55 }
56
57 /****************************************************************************/
58
mri_complex_phase(MRI_IMAGE * im)59 MRI_IMAGE *mri_complex_phase( MRI_IMAGE *im )
60 {
61 register int ii , npix ;
62 MRI_IMAGE *newImg ;
63 float *nar ; complex *iar ;
64
65 if( im->kind != MRI_complex ){
66 fprintf( stderr , "mri_complex_phase illegal image type!\n" ) ;
67 MRI_FATAL_ERROR ;
68 }
69
70 npix = im->nvox ;
71 newImg = mri_new_conforming( im , MRI_float ) ;
72 MRI_COPY_AUX( newImg , im ) ;
73 iar = MRI_COMPLEX_PTR(im) ; nar = MRI_FLOAT_PTR(newImg) ;
74
75 for( ii=0 ; ii < npix ; ii++ )
76 nar[ii] = atan2( iar[ii].i , iar[ii].r ) ;
77
78 return newImg ;
79 }
80
81 /*--------------------------------------------------------------------------*/
82
complex_abs(complex z)83 float complex_abs( complex z ) /* 24 Aug 2009 */
84 {
85 float x , y , val ;
86
87 x = fabsf(z.r) ; floatfix(x) ;
88 y = fabsf(z.i) ; floatfix(y) ;
89 if( x > y && x > 0.0f )
90 val = x*sqrtf(1.0f+(y*y)/(x*x)) ;
91 else if( y > x && y > 0.0f )
92 val = y*sqrtf(1.0f+(x*x)/(y*y)) ;
93 else
94 val = x*1.414214f ; /* case x==y */
95
96 floatfix(val) ; return val ;
97 }
98
99 /***************************************************************************/
100
mri_complex_abs(MRI_IMAGE * im)101 MRI_IMAGE *mri_complex_abs( MRI_IMAGE *im )
102 {
103 register int ii , npix ;
104 MRI_IMAGE *newImg ;
105 float *nar ; complex *iar ;
106
107 if( im->kind != MRI_complex ){
108 fprintf( stderr , "mri_complex_abs illegal type!\n" ) ;
109 MRI_FATAL_ERROR ;
110 }
111
112 npix = im->nvox ;
113 newImg = mri_new_conforming( im , MRI_float ) ;
114 MRI_COPY_AUX( newImg , im ) ;
115 iar = MRI_COMPLEX_PTR(im) ; nar = MRI_FLOAT_PTR(newImg) ;
116
117 for( ii=0 ; ii < npix ; ii++ ) nar[ii] = complex_abs(iar[ii]) ;
118
119 return newImg ;
120 }
121
122 /***************************************************************************/
123
mri_complex_real(MRI_IMAGE * im)124 MRI_IMAGE *mri_complex_real( MRI_IMAGE *im )
125 {
126 register int ii , npix ;
127 MRI_IMAGE *newImg ;
128 float *nar ; complex *iar ;
129
130 if( im->kind != MRI_complex ){
131 fprintf( stderr , "mri_complex_abs illegal type!\n" ) ;
132 MRI_FATAL_ERROR ;
133 }
134
135 npix = im->nvox ;
136 newImg = mri_new_conforming( im , MRI_float ) ;
137 MRI_COPY_AUX( newImg , im ) ;
138 iar = MRI_COMPLEX_PTR(im) ; nar = MRI_FLOAT_PTR(newImg) ;
139
140 for( ii=0 ; ii < npix ; ii++ ) nar[ii] = iar[ii].r ;
141
142 return newImg ;
143 }
144
145 /***************************************************************************/
146
mri_complex_imag(MRI_IMAGE * im)147 MRI_IMAGE *mri_complex_imag( MRI_IMAGE *im )
148 {
149 register int ii , npix ;
150 MRI_IMAGE *newImg ;
151 float *nar ; complex *iar ;
152
153 if( im->kind != MRI_complex ){
154 fprintf( stderr , "mri_complex_abs illegal type!\n" ) ;
155 MRI_FATAL_ERROR ;
156 }
157
158 npix = im->nvox ;
159 newImg = mri_new_conforming( im , MRI_float ) ;
160 MRI_COPY_AUX( newImg , im ) ;
161 iar = MRI_COMPLEX_PTR(im) ; nar = MRI_FLOAT_PTR(newImg) ;
162
163 for( ii=0 ; ii < npix ; ii++ ) nar[ii] = iar[ii].i ;
164
165 return newImg ;
166 }
167