1 #include "mrilib.h"
2
main(int argc,char * argv[])3 int main( int argc , char * argv[] )
4 {
5 THD_3dim_dataset *dset1=NULL,*dset2=NULL , *mset=NULL , *oset ;
6 MRI_IMAGE *dbr1=NULL,*dbr2=NULL;
7 float fac1=0.0, fac2=0.0;
8 char *prefix = "cmplx" ;
9 byte *mask=NULL ; float *rr,*gg ; complex *cmx ;
10 int iarg=1 , dofim=0 , ii,i2,nvox , mode=0 ;
11 int ivol, nvol=0 ; /* 16 Oct 2019 rickr */
12 int ninput=0 ; /* have 1 or 2 datasets as input */
13
14 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
15 printf("Usage #1: 3dTwotoComplex [options] dataset\n"
16 "Usage #2: 3dTwotoComplex [options] dataset1 dataset2\n"
17 "\n"
18 "Converts 2 sub-bricks of input to a complex-valued dataset.\n"
19 "* If you have 1 input dataset, then sub-bricks [0..1] are\n"
20 " used to form the 2 components of the output.\n"
21 "* If you have 2 input datasets, then the [0] sub-brick of\n"
22 " each is used to form the components.\n"
23 "* Complex datasets have two 32-bit float components per voxel.\n"
24 "\n"
25 "Options:\n"
26 " -prefix ppp = Write output into dataset with prefix 'ppp'.\n"
27 " [default='cmplx']\n"
28 " -RI = The 2 inputs are real and imaginary parts.\n"
29 " [this is the default]\n"
30 " -MP = The 2 inputs are magnitude and phase.\n"
31 " [phase is in radians, please!]\n"
32 " -mask mset = Only output nonzero values where the mask\n"
33 " dataset 'mset' is nonzero.\n"
34 "Notes:\n"
35 "* Input datasets must be byte-, short-, or float-valued.\n"
36 "* You might calculate the component datasets using 3dcalc.\n"
37 "* At present, there is limited support for complex datasets.\n"
38 " About the only thing you can do is display them in 2D\n"
39 " slice windows in AFNI.\n"
40 "\n"
41 "-- RWCox - March 2006\n"
42 ) ;
43 PRINT_COMPILE_DATE ; exit(0) ;
44 }
45
46 mainENTRY("3dTwotoComplex main"); machdep();
47 AFNI_logger("3dTwotoComplex",argc,argv);
48
49 /*-- options --*/
50
51 #define GOOD_TYPE(tt) ((tt)==MRI_short || (tt)==MRI_byte || (tt)==MRI_float)
52
53 while( iarg < argc && argv[iarg][0] == '-' ){
54
55 if( strcmp(argv[iarg],"-prefix") == 0 ){
56 prefix = argv[++iarg] ;
57 if( !THD_filename_ok(prefix) )
58 ERROR_exit("-prefix %s is illegal!\n",prefix) ;
59 iarg++ ; continue ;
60 }
61
62 if( strcmp(argv[iarg],"-mask") == 0 ){
63 if( mset != NULL )
64 ERROR_exit("can't have 2 -mask options!\n");
65 mset = THD_open_dataset( argv[++iarg] ) ;
66 if( !ISVALID_DSET(mset) )
67 ERROR_exit("can't open -mset %s\n",argv[iarg]);
68 if( !GOOD_TYPE( DSET_BRICK_TYPE(mset,0) ) )
69 ERROR_exit("-mset %s has invalid data type\n",argv[iarg]);
70 DSET_load(mset) ; CHECK_LOAD_ERROR(mset) ;
71 iarg++ ; continue ;
72 }
73
74 if( strcmp(argv[iarg],"-RI") == 0 ){
75 mode = 0 ; iarg++ ; continue ;
76 }
77 if( strcmp(argv[iarg],"-MP") == 0 ){
78 mode = 1 ; iarg++ ; continue ;
79 }
80
81 ERROR_exit("ILLEGAL option: %s\n",argv[iarg]) ;
82 }
83
84 if( iarg >= argc )
85 ERROR_exit("No datasets on command line!?\n");
86
87 /*-- read input datasets --*/
88
89 #define OPENIT(ds,aa) \
90 do{ ds = THD_open_dataset(aa) ; \
91 if( !ISVALID_DSET(ds) ) ERROR_exit("Can't open dataset %s\n",aa); \
92 DSET_load(ds) ; CHECK_LOAD_ERROR(ds) ; \
93 } while(0)
94
95 if( iarg == argc-1 ){ /* one dataset */
96
97 OPENIT(dset1,argv[iarg]) ;
98
99 nvol = DSET_NVALS(dset1) ;
100 if( (nvol % 2) || nvol < 2 )
101 ERROR_exit("Dataset %s needs an even number (>0) of sub-bricks!\n",
102 argv[iarg]);
103 nvol /= 2; /* convert to the number of output volumes */
104
105 if( !GOOD_TYPE( DSET_BRICK_TYPE(dset1,0) ) ||
106 !GOOD_TYPE( DSET_BRICK_TYPE(dset1,1) ) )
107 ERROR_exit("ILLEGAL dataset type in %s\n",argv[iarg]);
108
109 /* use dset2 in any case, to generalize in loop */
110 dset2 = dset1 ;
111
112 /* only 1 input */
113 ninput = 1 ;
114
115 } else if( iarg+1 < argc ){ /* two datasets */
116
117 OPENIT(dset1,argv[iarg]) ;
118 if( !GOOD_TYPE( DSET_BRICK_TYPE(dset1,0) ) )
119 ERROR_exit("ILLEGAL dataset type in %s\n",argv[iarg]) ;
120
121 OPENIT(dset2,argv[iarg+1]) ;
122 if( !GOOD_TYPE( DSET_BRICK_TYPE(dset2,0) ) )
123 ERROR_exit("ILLEGAL dataset type in %s\n",argv[iarg+1]);
124
125 nvol = DSET_NVALS(dset1) ;
126 if( nvol != DSET_NVALS(dset2) )
127 ERROR_exit("Datasets must have the same number of volumes");
128
129 if( ! EQUIV_GRIDS(dset1, dset2) )
130 ERROR_exit("Dataset grids do not match\n") ;
131
132 if( iarg+1 < argc-1 )
133 WARNING_message("extra arguments on command line: %s ...\n",
134 argv[iarg+2] ) ;
135
136 /* have 2 inputs */
137 ninput = 2 ;
138 } else { /* this is bad */
139 ERROR_exit("Need either 1 or 2 datasets on command line!?\n");
140 }
141
142 /* make a mask, if any */
143
144 if( mset != NULL ){
145 if( ! EQUIV_GRIDS_NXYZ(dset1, mset) )
146 ERROR_exit("Mask and dataset nx/y/z don't match!\n");
147 mask = THD_makemask( mset , 0 , 1.0,0.0 ) ;
148 if( mask == NULL )
149 ERROR_exit("Can't make mask!\n");
150 DSET_unload(mset) ;
151 }
152
153 /* convert input sub-bricks to floats */
154
155 nvox = DSET_NVOX(dset1) ;
156
157 /* allocate computational memory */
158 rr = (float *) calloc(sizeof(float),nvox) ;
159 gg = (float *) calloc(sizeof(float),nvox) ;
160
161 /* make output dataset */
162
163 oset = EDIT_empty_copy( dset1 ) ;
164 EDIT_dset_items( oset ,
165 ADN_prefix , prefix ,
166 ADN_datum_all, MRI_complex ,
167 ADN_ntt , nvol ,
168 ADN_type , (dofim) ? HEAD_FUNC_TYPE : HEAD_ANAT_TYPE,
169 ADN_func_type, (dofim) ? FUNC_FIM_TYPE : ANAT_SPGR_TYPE,
170 ADN_none ) ;
171
172 for( ivol=0; ivol < nvol; ivol++ ) {
173
174 /* choose input indices, and note that if ninput==1, dset2==dset1 */
175 if( ninput == 2 ) {
176 ii = ivol;
177 i2 = ivol;
178 } else {
179 ii = 2*ivol;
180 i2 = ii + 1;
181 }
182
183 /* set each data block and scale factor */
184 dbr1 = DSET_BRICK(dset1,ii) ; fac1 = DSET_BRICK_FACTOR(dset1,ii) ;
185 dbr2 = DSET_BRICK(dset2,i2) ; fac2 = DSET_BRICK_FACTOR(dset2,i2) ;
186
187 /* explicitly clear memory, though it shouldn't really matter */
188 memset(rr, '\0', nvox*sizeof(float));
189 memset(gg, '\0', nvox*sizeof(float));
190
191 /* convert to float */
192 if( fac1 == 0.0 ) fac1 = 1.0 ;
193 EDIT_coerce_scale_type( nvox , fac1 ,
194 DSET_BRICK_TYPE(dset1,ii), mri_data_pointer(dbr1),
195 MRI_float , rr ) ;
196
197 if( fac2 == 0.0 ) fac2 = 1.0 ;
198 EDIT_coerce_scale_type( nvox , fac2 ,
199 DSET_BRICK_TYPE(dset2,i2), mri_data_pointer(dbr2),
200 MRI_float , gg ) ;
201
202 /* allocate space for the new volume, to be attached to new dset */
203 cmx = (complex *) calloc(sizeof(complex),nvox) ;
204
205 #define T2C(a,b) \
206 (mode==0) ? CMPLX((a),(b)) : CMPLX((a)*cos(b),(a)*sin(b))
207
208 /* merge inputs to output */
209 if( mask != NULL ){
210 for( ii=0 ; ii < nvox ; ii++ )
211 if( mask[ii] ) cmx[ii] = T2C(rr[ii],gg[ii]) ;
212 } else {
213 for( ii=0 ; ii < nvox ; ii++ )
214 cmx[ii] = T2C(rr[ii],gg[ii]) ;
215 }
216
217 EDIT_substitute_brick( oset , ivol , MRI_complex , cmx ) ;
218 }
219
220 /* fly, be free! */
221 DSET_unload(dset1) ;
222 if( ninput == 2 ) DSET_unload(dset2) ;
223 free(mask) ; free(rr) ; free(gg) ;
224
225 /* make history */
226
227 tross_Copy_History( oset , dset1 ) ;
228 tross_Make_History( "3dTwotoComplex", argc,argv, oset ) ;
229
230 DSET_write( oset ) ;
231 WROTE_DSET( oset ) ;
232 DSET_unload( oset ) ;
233
234 exit(0) ;
235 }
236