1 /* minccmp.c */
2 /* */
3 /* Copyright Andrew Janke - a.janke@gmail.com */
4 /* Permission to use, copy, modify, and distribute this software and its */
5 /* documentation for any purpose and without fee is hereby granted, */
6 /* provided that the above copyright notice appear in all copies. The */
7 /* author makes no representations about the */
8 /* suitability of this software for any purpose. It is provided "as is" */
9 /* without express or implied warranty. */
10 /* */
11 /* calculates measures of similarity/difference between 2 or more volumes */
12 /* */
13 /* Measures used (sum(x) denotes the sum of x over a volume): */
14 /* RMSE - Root Mean Squared Error */
15 /* = sqrt( 1/n * sum((a-b)^2)) */
16 /* xcorr - Cross Correlation */
17 /* = sum((a*b)^2) / (sqrt(sum(a^2)) * sqrt(sum(b^2)) */
18 /* zscore - z-score differences */
19 /* = sum( |((a - mean(a)) / stdev(a)) - */
20 /* ((b - mean(b)) / stdev(b))| ) / nvox */
21 /* */
22 /* Tue Jun 17 11:31:10 EST 2003 - initial version inspired by voldiff and */
23 /* peter's compare_volumes */
24
25 #include <stdio.h>
26 #include <stdlib.h>
27 #include <unistd.h>
28 #include <math.h>
29 #include <float.h>
30 #include <voxel_loop.h>
31 #include <ParseArgv.h>
32
33 #ifndef FALSE
34 # define FALSE 0
35 #endif
36 #ifndef TRUE
37 # define TRUE 1
38 #endif
39
40 #define SQR2(x) ((x) * (x))
41
42 typedef struct {
43 double nvox;
44 double sum; /* sum of valid voxels */
45 double ssum; /* squared sum of valid voxels */
46 double min;
47 double max;
48
49 double mean;
50 double var;
51 double sd;
52
53 double sum_prd0; /* sum of product of file[x] with file[0] */
54 double ssum_add0; /* squared sum of addition of file[x] with file[0] */
55 double ssum_dif0; /* squared sum of difference of file[x] with file[0] */
56 double ssum_prd0; /* squared sum of product of file[x] with file[0] */
57
58 double sum_zdif0; /* sum of zscore differences of file[x] and file[0] */
59
60 /* result stores */
61 double rmse;
62 double xcorr;
63 double zscore;
64 double vratio;
65 } Vol_Data;
66
67 typedef struct {
68 int n_datafiles;
69
70 int mask;
71 int mask_idx;
72
73 /* individual volume data */
74 Vol_Data *vd;
75 } Loop_Data;
76
77 /* Function prototypes */
78 void pass_0(void *caller_data, long num_voxels,
79 int input_num_buffers, int input_vector_length,
80 double *input_data[],
81 int output_num_buffers, int output_vector_length,
82 double *output_data[], Loop_Info * loop_info);
83 void pass_1(void *caller_data, long num_voxels,
84 int input_num_buffers, int input_vector_length,
85 double *input_data[],
86 int output_num_buffers, int output_vector_length,
87 double *output_data[], Loop_Info * loop_info);
88 void print_result(char *title, double result);
89 void dump_stats(Loop_Data * ld);
90 void do_int_calcs(Loop_Data * ld);
91 void do_final_calcs(Loop_Data * ld);
92
93 /* Argument variables and table */
94 static int verbose = FALSE;
95 static int debug = FALSE;
96 static int quiet = FALSE;
97 static int clobber = FALSE;
98 static int max_buffer_size_in_kb = 4 * 1024;
99 static int check_dim_info = TRUE;
100 static char *mask_fname = NULL;
101 static double valid_range[2] = { -DBL_MAX, DBL_MAX };
102
103 static int do_all = FALSE;
104 static int do_ssq = FALSE;
105 static int do_rmse = FALSE;
106 static int do_xcorr = FALSE;
107 static int do_zscore = FALSE;
108 static int do_vratio = FALSE;
109
110 ArgvInfo argTable[] = {
111 {"-verbose", ARGV_CONSTANT, (char *)TRUE, (char *)&verbose,
112 "be verbose"},
113 {"-debug", ARGV_CONSTANT, (char *)TRUE, (char *)&debug,
114 "dump all stats info"},
115 {"-quiet", ARGV_CONSTANT, (char *)TRUE, (char *)&quiet,
116 "print requested values only"},
117 {"-clobber", ARGV_CONSTANT, (char *)TRUE, (char *)&clobber,
118 "clobber existing files"},
119 {"-max_buffer_size_in_kb", ARGV_INT, (char *)1, (char *)&max_buffer_size_in_kb,
120 "maximum size of internal buffers."},
121 {"-check_dimensions", ARGV_CONSTANT, (char *) TRUE, (char *) &check_dim_info,
122 "Check that files have matching dimensions (default)."},
123 {"-nocheck_dimensions", ARGV_CONSTANT, (char *) FALSE, (char *) &check_dim_info,
124 "Do not check that files have matching dimensions."},
125
126 {NULL, ARGV_HELP, (char *)NULL, (char *)NULL,
127 "\nVoxel selection options (applies to first volume ONLY):"},
128 {"-floor", ARGV_FLOAT, (char *)1, (char *)&valid_range[0],
129 "Ignore voxels below this value. (incl)"},
130 {"-ceil", ARGV_FLOAT, (char *)1, (char *)&valid_range[1],
131 "Ignore voxels above this value. (incl)"},
132 {"-range", ARGV_FLOAT, (char *)2, (char *)&valid_range,
133 "Ignore voxels outside the range. (incl)"},
134
135 {"-mask", ARGV_STRING, (char *)1, (char *)&mask_fname,
136 "Use <mask.mnc> for calculations."},
137
138 {NULL, ARGV_HELP, (char *)NULL, (char *)NULL,
139 "\nImage Statistics (printed in this order)"},
140 {"-all", ARGV_CONSTANT, (char *)TRUE, (char *)&do_all,
141 "all statistics (default)."},
142 {"-ssq", ARGV_CONSTANT, (char *)TRUE, (char *)&do_ssq,
143 "sum of squared difference (2 volumes)"},
144 {"-rmse", ARGV_CONSTANT, (char *)TRUE, (char *)&do_rmse,
145 "root mean squared error (2 volumes)"},
146 {"-xcorr", ARGV_CONSTANT, (char *)TRUE, (char *)&do_xcorr,
147 "cross correlation (2 volumes)"},
148 {"-zscore", ARGV_CONSTANT, (char *)TRUE, (char *)&do_zscore,
149 "z-score (2 volumes)"},
150 // {"-vr", ARGV_CONSTANT, (char *)TRUE, (char *)&do_vratio,
151 // "variance ratio (2 volumes)"},
152
153 // {NULL, ARGV_HELP, (char *)NULL, (char *)NULL,
154 // "\nBinary Image Only Statistics"},
155 // {"-kappa", ARGV_CONSTANT, (char *)TRUE, (char *)&do_kappa,
156 // "all statistics (default)."},
157
158 {NULL, ARGV_HELP, NULL, NULL, ""},
159 {NULL, ARGV_END, NULL, NULL, NULL}
160 };
161
main(int argc,char * argv[])162 int main(int argc, char *argv[]){
163 char **infiles;
164 int n_infiles;
165
166 Loop_Options *loop_opt;
167 Loop_Data ld;
168
169 int i;
170
171 /* Get arguments */
172 if(ParseArgv(&argc, argv, argTable, 0) || (argc < 3)){
173 fprintf(stderr, "\nUsage: %s [options] <in1.mnc> <in2.mnc> [<inn.mnc>]\n", argv[0]);
174 fprintf(stderr, " %s -help\n\n", argv[0]);
175 return (EXIT_FAILURE);
176 }
177 n_infiles = argc - 1;
178 infiles = &argv[1];
179
180 /* check arguments */
181 if(!do_rmse && !do_xcorr && !do_zscore && !do_vratio){
182 do_all = TRUE;
183 }
184 if(do_all){
185 do_ssq = do_rmse = do_xcorr = do_zscore = do_vratio = TRUE;
186 }
187
188 /* check for infiles */
189 if(verbose){
190 fprintf(stderr, "\n+++ infiles +++\n");
191 }
192 for(i = 0; i < n_infiles; i++){
193 if(verbose){
194 fprintf(stderr, " | [%02d]: %s\n", i, infiles[i]);
195 }
196 if(access(infiles[i], F_OK) != 0){
197 fprintf(stderr, "%s: Couldn't find %s\n", argv[0], infiles[i]);
198 exit(EXIT_FAILURE);
199 }
200 }
201
202 /* set up Loop_Data struct and mask file */
203 ld.n_datafiles = n_infiles;
204 if(mask_fname != NULL){
205 if(verbose){
206 fprintf(stderr, " | mask: %s\n", mask_fname);
207 }
208 if(access(mask_fname, F_OK) != 0){
209 fprintf(stderr, "%s: Couldn't find mask file: %s\n", argv[0], mask_fname);
210 exit(EXIT_FAILURE);
211 }
212
213 ld.mask = TRUE;
214 ld.mask_idx = n_infiles;
215
216 infiles[n_infiles] = mask_fname;
217 n_infiles++;
218 }
219 else {
220 ld.mask = FALSE;
221 ld.mask_idx = 0;
222 }
223
224 /* allocate space and initialise volume stats data */
225 ld.vd = (Vol_Data *) malloc(sizeof(Vol_Data) * ld.n_datafiles);
226 for(i = 0; i < ld.n_datafiles; i++){
227 ld.vd[i].nvox = 0;
228 ld.vd[i].sum = 0;
229 ld.vd[i].ssum = 0;
230 ld.vd[i].min = DBL_MAX;
231 ld.vd[i].max = -DBL_MAX;
232
233 ld.vd[i].sum_prd0 = 0;
234 ld.vd[i].ssum_add0 = 0;
235 ld.vd[i].ssum_dif0 = 0;
236 ld.vd[i].ssum_prd0 = 0;
237
238 ld.vd[i].mean = 0;
239 ld.vd[i].var = 0;
240 ld.vd[i].sd = 0;
241
242 ld.vd[i].rmse = 0.0;
243 ld.vd[i].xcorr = 0.0;
244 ld.vd[i].zscore = 0.0;
245 ld.vd[i].vratio = 0.0;
246 }
247
248 /* set up and do voxel_loop(s) */
249 loop_opt = create_loop_options();
250 set_loop_verbose(loop_opt, verbose);
251 set_loop_buffer_size(loop_opt, (long)1024 * max_buffer_size_in_kb);
252 set_loop_check_dim_info(loop_opt, check_dim_info);
253
254 /* first pass */
255 voxel_loop(n_infiles, infiles, 0, NULL, NULL, loop_opt, pass_0, (void *)&ld);
256
257 /* intermediate calculations */
258 do_int_calcs(&ld);
259
260 /* run the second pass if we have to */
261 if(do_zscore){
262 voxel_loop(n_infiles, infiles, 0, NULL, NULL, loop_opt, pass_1, (void *)&ld);
263 }
264
265 /* final calculations */
266 do_final_calcs(&ld);
267
268 free_loop_options(loop_opt);
269
270 if(debug){
271 dump_stats(&ld);
272 }
273
274 /* calculate and print result(s) */
275 if(do_all && !quiet){
276 fprintf(stdout, "file[0]: %s\n", infiles[0]);
277 fprintf(stdout, "file[1]: %s\n", infiles[1]);
278 fprintf(stdout, "mask file: %s\n", mask_fname);
279 }
280 if(do_ssq){
281 print_result("ssq: ", ld.vd[1].ssum_dif0);
282 }
283 if(do_rmse){
284 print_result("rmse: ", ld.vd[1].rmse);
285 }
286 if(do_xcorr){
287 print_result("xcorr: ", ld.vd[1].xcorr);
288 }
289 if(do_zscore){
290 print_result("zscore: ", ld.vd[1].zscore);
291 }
292 if(!quiet){
293 fprintf(stdout, "\n");
294 }
295
296 return EXIT_SUCCESS;
297 }
298
299 /* voxel loop function for first pass through data */
pass_0(void * caller_data,long num_voxels,int input_num_buffers,int input_vector_length,double * input_data[],int output_num_buffers,int output_vector_length,double * output_data[],Loop_Info * loop_info)300 void pass_0(void *caller_data, long num_voxels,
301 int input_num_buffers, int input_vector_length,
302 double *input_data[],
303 int output_num_buffers, int output_vector_length,
304 double *output_data[], Loop_Info * loop_info){
305 long ivox;
306 double valuei, value0;
307 int i;
308
309 /* get pointer to loop data */
310 Loop_Data *ld = (Loop_Data *)caller_data;
311
312 /* shut the compiler up - yes I _know_ I don't use these */
313 (void)output_num_buffers;
314 (void)output_vector_length;
315 (void)output_data;
316 (void)loop_info;
317
318 /* sanity check */
319 if((input_num_buffers < 2) || (output_num_buffers != 0)){
320 fprintf(stderr, "Bad arguments to pass_0\n");
321 exit(EXIT_FAILURE);
322 }
323
324 /* for each voxel */
325 for(ivox = num_voxels * input_vector_length; ivox--;){
326
327 /* skip voxels out of the mask region */
328 if(ld->mask && !(int)input_data[ld->mask_idx][ivox]){
329 continue;
330 }
331
332 value0 = input_data[0][ivox];
333 if(value0 >= valid_range[0] && value0 <= valid_range[1]){
334
335 /* for each volume */
336 for(i = 0; i < ld->n_datafiles; i++){
337
338 valuei = input_data[i][ivox];
339
340 /* various voxel sums */
341 ld->vd[i].nvox++;
342 ld->vd[i].sum += valuei;
343 ld->vd[i].ssum += SQR2(valuei);
344
345 if(i != 0){
346 ld->vd[i].sum_prd0 += valuei * value0;
347
348 ld->vd[i].ssum_add0 += SQR2(valuei + value0);
349 ld->vd[i].ssum_dif0 += SQR2(valuei - value0);
350 ld->vd[i].ssum_prd0 += SQR2(valuei * value0);
351 }
352
353 /* min and max */
354 if(valuei < ld->vd[i].min){
355 ld->vd[i].min = valuei;
356 }
357 else if(valuei > ld->vd[i].max){
358 ld->vd[i].max = valuei;
359 }
360 }
361 }
362 }
363
364 return;
365 }
366
367 /* intermediate calculations */
do_int_calcs(Loop_Data * ld)368 void do_int_calcs(Loop_Data * ld){
369 int i;
370 double denom;
371
372 for(i = 0; i < ld->n_datafiles; i++){
373 /* mean */
374 ld->vd[i].mean = ld->vd[i].sum / ld->vd[i].nvox;
375
376 /* variance */
377 ld->vd[i].var = ((ld->vd[i].nvox * ld->vd[i].ssum) - SQR2(ld->vd[i].sum)) /
378 (ld->vd[i].nvox * (ld->vd[i].nvox - 1));
379
380 /* sd */
381 ld->vd[i].sd = sqrt(ld->vd[i].var);
382
383 /* RMSE */
384 ld->vd[i].rmse = sqrt((1.0 / ld->vd[0].nvox) * ld->vd[i].ssum_dif0);
385
386 /* xcorr */
387 denom = sqrt(ld->vd[0].ssum * ld->vd[i].ssum);
388 ld->vd[i].xcorr = (denom == 0.0) ? 0.0 : ld->vd[i].sum_prd0 / denom;
389 }
390 }
391
392 /* voxel loop function for second pass through data */
pass_1(void * caller_data,long num_voxels,int input_num_buffers,int input_vector_length,double * input_data[],int output_num_buffers,int output_vector_length,double * output_data[],Loop_Info * loop_info)393 void pass_1(void *caller_data, long num_voxels,
394 int input_num_buffers, int input_vector_length,
395 double *input_data[],
396 int output_num_buffers, int output_vector_length,
397 double *output_data[], Loop_Info * loop_info){
398 long ivox;
399 double valuei, value0;
400 int i;
401
402 /* get pointer to loop data */
403 Loop_Data *ld = (Loop_Data *)caller_data;
404
405 /* shut the compiler up - yes I _know_ I don't use these */
406 (void)output_num_buffers;
407 (void)output_vector_length;
408 (void)output_data;
409 (void)loop_info;
410
411 /* sanity check */
412 if((input_num_buffers < 2) || (output_num_buffers != 0)){
413 fprintf(stderr, "Bad arguments to pass_1\n");
414 exit(EXIT_FAILURE);
415 }
416
417 /* for each voxel */
418 for(ivox = num_voxels * input_vector_length; ivox--;){
419
420 /* skip voxels out of the mask region */
421 if(ld->mask && !(int)input_data[ld->mask_idx][ivox]){
422 continue;
423 }
424
425 value0 = input_data[0][ivox];
426 if(value0 >= valid_range[0] && value0 <= valid_range[1]){
427
428 /* for each volume */
429 for(i = 0; i < ld->n_datafiles; i++){
430
431 valuei = input_data[i][ivox];
432
433 if(i != 0){
434
435 /* zscore total */
436 ld->vd[i].sum_zdif0 +=
437 fabs(((value0 - ld->vd[0].mean) / ld->vd[0].sd) -
438 ((valuei - ld->vd[i].mean) / ld->vd[i].sd));
439
440 }
441 }
442 }
443 }
444
445 return;
446 }
447
448 /* final calculations */
do_final_calcs(Loop_Data * ld)449 void do_final_calcs(Loop_Data * ld){
450 int i;
451
452 for(i = 0; i < ld->n_datafiles; i++){
453 /* zscore */
454 ld->vd[i].zscore = ld->vd[i].sum_zdif0 / ld->vd[i].nvox;
455 }
456 }
457
458 /* dirty little function to print out results */
print_result(char * title,double result)459 void print_result(char *title, double result){
460 if(!quiet){
461 fprintf(stdout, "%s", title);
462 }
463 fprintf(stdout, "%.10g\n", result);
464 }
465
466 /* debug function to dump stats structure */
dump_stats(Loop_Data * ld)467 void dump_stats(Loop_Data * ld){
468 int i;
469
470 fprintf(stdout, " + Main Loop data structure\n");
471 fprintf(stdout, " | n_datafiles %d\n", ld->n_datafiles);
472 fprintf(stdout, " | mask %d\n", ld->mask);
473 fprintf(stdout, " | mask_idx %d\n", ld->mask_idx);
474
475 fprintf(stdout, " +++ volume data stats\n");
476 for(i = 0; i < ld->n_datafiles; i++){
477 fprintf(stdout, " |---------------------------------\n");
478 fprintf(stdout, " | [%02d] nvox %10g\n", i, ld->vd[i].nvox);
479 fprintf(stdout, " | [%02d] sum %.10g\n", i, ld->vd[i].sum);
480 fprintf(stdout, " | [%02d] ssum %.10g\n", i, ld->vd[i].ssum);
481 fprintf(stdout, " | [%02d] min %.10g\n", i, ld->vd[i].min);
482 fprintf(stdout, " | [%02d] max %.10g\n", i, ld->vd[i].max);
483 fprintf(stdout, " | [%02d] mean %.10g\n", i, ld->vd[i].mean);
484 fprintf(stdout, " | [%02d] var %.10g\n", i, ld->vd[i].var);
485 fprintf(stdout, " | [%02d] sd %.10g\n", i, ld->vd[i].sd);
486 fprintf(stdout, " | [%02d] sum_prd0 %.10g\n", i, ld->vd[i].sum_prd0);
487 fprintf(stdout, " | [%02d] ssum_add0 %.10g\n", i, ld->vd[i].ssum_add0);
488 fprintf(stdout, " | [%02d] ssum_dif0 %.10g\n", i, ld->vd[i].ssum_dif0);
489 fprintf(stdout, " | [%02d] ssum_prd0 %.10g\n", i, ld->vd[i].ssum_prd0);
490 fprintf(stdout, " | [%02d] sum_zdif0 %.10g\n", i, ld->vd[i].sum_zdif0);
491
492 fprintf(stdout, " | [%02d] rmse %.10g\n", i, ld->vd[i].rmse);
493 fprintf(stdout, " | [%02d] xcorr %.10g\n", i, ld->vd[i].xcorr);
494 fprintf(stdout, " | [%02d] zscore %.10g\n", i, ld->vd[i].zscore);
495 fprintf(stdout, " | [%02d] vratio %.10g\n", i, ld->vd[i].vratio);
496 }
497 }
498