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 "parser.h"
8 #include "mrilib.h"
9
10 /*-------------------------- global data --------------------------*/
11
12 static int CALC_datum = -1 ;
13 static int CALC_nvox = -1 ;
14 static PARSER_code * CALC_code = NULL ;
15
16 static int CALC_nx=-1 , CALC_ny=-1 ; /* 27 Jul 2000 */
17
18 static MRI_IMAGE * CALC_im[26] ;
19 static int CALC_type[26] ;
20 static byte * CALC_byte[26] ;
21 static short * CALC_short[26] ;
22 static float * CALC_float[26] ;
23
24 static char CALC_output_name[256] = "imcalc.out" ;
25
26 /*--------------------------- prototypes ---------------------------*/
27 void CALC_read_opts( int , char ** ) ;
28 void CALC_Syntax(void) ;
29
30 /*--------------------------------------------------------------------
31 read the arguments, load the global variables
32 ----------------------------------------------------------------------*/
33
CALC_read_opts(int argc,char * argv[])34 void CALC_read_opts( int argc , char * argv[] )
35 {
36 int nopt = 1 ;
37 int ids ;
38
39 for( ids=0 ; ids < 26 ; ids++ ){
40 CALC_type[ids] = -1 ;
41 }
42
43 while( nopt < argc && argv[nopt][0] == '-' ){
44
45 /**** -nxy nx ny ****/
46
47 if( strncmp(argv[nopt],"-nxy",4) == 0 ){
48 if( ++nopt >= argc-1 ){
49 fprintf(stderr,"need 2 arguments after -nxy!\n"); exit(1);
50 }
51 CALC_nx = strtol(argv[nopt++],NULL,10) ;
52 CALC_ny = strtol(argv[nopt++],NULL,10) ;
53 continue ; /* go to next arg */
54 }
55
56 /**** -datum type ****/
57
58 if( strncmp(argv[nopt],"-datum",6) == 0 ){
59 if( ++nopt >= argc ){
60 fprintf(stderr,"need an argument after -datum!\n") ; exit(1) ;
61 }
62 if( strcmp(argv[nopt],"short") == 0 ){
63 CALC_datum = MRI_short ;
64 } else if( strcmp(argv[nopt],"float") == 0 ){
65 CALC_datum = MRI_float ;
66 } else if( strcmp(argv[nopt],"byte") == 0 ){
67 CALC_datum = MRI_byte ;
68 } else {
69 fprintf(stderr,"-datum of type '%s' is not supported in 3dmerge!\n",
70 argv[nopt] ) ;
71 exit(1) ;
72 }
73 nopt++ ; continue ; /* go to next arg */
74 }
75
76 /**** -output name ****/
77
78 if( strncmp(argv[nopt],"-output",6) == 0 ){
79 nopt++ ;
80 if( nopt >= argc ){
81 fprintf(stderr,"need argument after -output!\n") ; exit(1) ;
82 }
83 strncpy( CALC_output_name , argv[nopt++] , 255 ) ; CALC_output_name[255] = '\0' ;
84 continue ;
85 }
86
87 /**** -expr expression ****/
88
89 if( strncmp(argv[nopt],"-expr",4) == 0 ){
90 if( CALC_code != NULL ){
91 fprintf(stderr,"cannot have 2 -expr options!\n") ; exit(1) ;
92 }
93 nopt++ ;
94 if( nopt >= argc ){
95 fprintf(stderr,"need argument after -expr!\n") ; exit(1) ;
96 }
97 CALC_code = PARSER_generate_code( argv[nopt++] ) ;
98 if( CALC_code == NULL ){
99 fprintf(stderr,"illegal expression!\n") ; exit(1) ;
100 }
101 continue ;
102 }
103
104 /**** -[a-z] filename ****/
105
106 ids = strlen(argv[nopt]) ;
107
108 if( (argv[nopt][1] >= 'a' && argv[nopt][1] <= 'z') && (ids == 2) ){
109
110 int ival , nxyz ;
111 MRI_IMAGE * im ;
112
113 ival = argv[nopt][1] - 'a' ;
114 if( CALC_im[ival] != NULL ){
115 fprintf(stderr,"can't open duplicate %s images!\n",argv[nopt]) ;
116 exit(1) ;
117 }
118
119 nopt++ ;
120 if( nopt >= argc ){
121 fprintf(stderr,"need argument after %s!\n",argv[nopt-1]) ; exit(1) ;
122 }
123
124 im = mri_read_just_one( argv[nopt++] ) ;
125 if( im == NULL ){
126 fprintf(stderr,"can't open image %s\n",argv[nopt-1]) ; exit(1) ;
127 }
128
129 nxyz = im->nvox ;
130 if( CALC_nvox < 0 ){
131 CALC_nvox = nxyz ;
132 } else if( nxyz != CALC_nvox ){
133 fprintf(stderr,"image %s differs in size from others\n",argv[nopt-1]);
134 exit(1) ;
135 }
136
137 CALC_type[ival] = im->kind ;
138 CALC_im[ival] = im ;
139
140 switch( CALC_type[ival] ){
141 default: fprintf(stderr,"illegal datum in image %s\n",argv[nopt-1]);
142 exit(1) ;
143
144 case MRI_short: CALC_short[ival] = MRI_SHORT_PTR(im) ; break ;
145 case MRI_float: CALC_float[ival] = MRI_FLOAT_PTR(im) ; break ;
146 case MRI_byte : CALC_byte [ival] = MRI_BYTE_PTR(im) ; break ;
147 }
148
149 if( CALC_datum < 0 ) CALC_datum = CALC_type[ival] ;
150 continue ;
151 }
152
153 fprintf(stderr,"Unknown option: %s\n",argv[nopt]) ; exit(1) ;
154
155 } /* end of loop over options */
156
157 /*** cleanup ***/
158
159 if( nopt < argc ){
160 fprintf(stderr,"Extra command line arguments puzzle me! %s ...\n",argv[nopt]) ;
161 exit(1) ;
162 }
163
164 for( ids=0 ; ids < 26 ; ids++ ) if( CALC_im[ids] != NULL ) break ;
165
166 if( ids == 26 && (CALC_nx < 2 || CALC_ny < 2) ){
167 fprintf(stderr,"No input images given!\n") ; exit(1) ;
168 }
169
170 if( CALC_code == NULL ){
171 fprintf(stderr,"No expression given!\n") ; exit(1) ;
172 }
173
174 return ;
175 }
176
177 /*------------------------------------------------------------------*/
178
CALC_Syntax(void)179 void CALC_Syntax(void)
180 {
181 printf(
182 "Do arithmetic on 2D images, pixel-by-pixel.\n"
183 "Usage: imcalc options\n"
184 "where the options are:\n"
185 ) ;
186
187 printf(
188 " -datum type = Coerce the output data to be stored as the given type,\n"
189 " which may be byte, short, or float.\n"
190 " [default = datum of first input image on command line]\n"
191 " -a dname = Read image 'dname' and call the voxel values 'a'\n"
192 " in the expression. 'a' may be any letter from 'a' to 'z'.\n"
193 " ** If some letter name is used in the expression, but not\n"
194 " present in one of the image options here, then that\n"
195 " variable is set to 0.\n"
196 " -expr \"expression\"\n"
197 " Apply the expression within quotes to the input images,\n"
198 " one voxel at a time, to produce the output image.\n"
199 " (\"sqrt(a*b)\" to compute the geometric mean, for example)\n"
200 " -output name = Use 'name' for the output image filename.\n"
201 " [default='imcalc.out']\n"
202 "\n"
203 "See the output of '3dcalc -help' for details on what kinds of expressions\n"
204 "are possible. Note that complex-valued images cannot be processed (byte,\n"
205 "short, and float are OK).\n"
206 ) ;
207 exit(0) ;
208 }
209
210 /*------------------------------------------------------------------*/
211
main(int argc,char * argv[])212 int main( int argc , char * argv[] )
213 {
214 double atoz[26] ;
215 int ii , ids ;
216 MRI_IMAGE * new_im=NULL ;
217 float * fnew ;
218
219 /*** read input options ***/
220
221 if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ) CALC_Syntax() ;
222
223 machdep() ;
224
225 CALC_read_opts( argc , argv ) ;
226
227 /*** make output image (always float datum) ***/
228
229 if( THD_filesize(CALC_output_name) > 0 ){
230 fprintf(stderr,
231 "*** Output file %s already exists -- cannot continue!\n",
232 CALC_output_name ) ;
233 exit(1) ;
234 }
235
236 for( ids=0 ; ids < 26 ; ids++ ) if( CALC_im[ids] != NULL ) break ;
237 if( ids < 26 ){
238 new_im = mri_new_conforming( CALC_im[ids] , MRI_float ) ;
239 } else if( CALC_nx > 1 && CALC_ny > 1 ){
240 new_im = mri_new( CALC_nx , CALC_ny , MRI_float ) ;
241 CALC_nvox = new_im->nvox ;
242 }
243 fnew = MRI_FLOAT_PTR(new_im) ;
244
245 for( ids=0 ; ids < 26 ; ids++ ) atoz[ids] = 0.0 ;
246
247 /*** loop over voxels and compute result ***/
248
249 for( ii=0 ; ii < CALC_nvox ; ii++ ){
250
251 for( ids=0 ; ids < 26 ; ids++ ){
252 switch( CALC_type[ids] ){
253 case MRI_short: atoz[ids] = CALC_short[ids][ii] ; break;
254 case MRI_float: atoz[ids] = CALC_float[ids][ii] ; break;
255 case MRI_byte : atoz[ids] = CALC_byte [ids][ii] ; break;
256 }
257 }
258
259 fnew[ii] = PARSER_evaluate_one( CALC_code , atoz ) ;
260 }
261
262 /*** toss input images ***/
263
264 for( ids=0 ; ids < 26 ; ids++ )
265 if( CALC_im[ids] != NULL ) mri_free( CALC_im[ids] ) ;
266
267 /*** scale to output image, if needed ***/
268
269 switch( CALC_datum ){
270
271 case MRI_short:{
272 float top ;
273 MRI_IMAGE * shim ;
274
275 top = mri_maxabs( new_im ) ;
276 if( top < 32767.0 ) shim = mri_to_short( 1.0 , new_im ) ;
277 else shim = mri_to_short_scl( 0.0 , 10000.0 , new_im ) ;
278
279 mri_free(new_im) ; new_im = shim ;
280 }
281 break ;
282
283 case MRI_byte:{
284 float top ;
285 MRI_IMAGE * bim ;
286
287 top = mri_maxabs( new_im ) ;
288 if( top < 255.0 ) bim = mri_to_byte_scl( 1.0 , 0.0 , new_im ) ;
289 else bim = mri_to_byte_scl( 0.0 , 255.0 , new_im ) ;
290
291 mri_free(new_im) ; new_im = bim ;
292 }
293 break ;
294 }
295
296 /*** done ***/
297
298 mri_write( CALC_output_name , new_im ) ;
299 exit(0) ;
300 }
301