1 #include "mrilib.h"
2
3 /*---------------------------------------------------------------------
4 01 Feb 2001: determine if 2 datasets don't match in some way for
5 voxel-by-voxel comparison.
6 return = 0 if they are good in all ways
7 return = binary OR of MISMATCH_ codes (cf. 3ddata.h) if
8 something in their dataxes don't match
9 A generalization of EQUIV_GRIDS
10 -----------------------------------------------------------------------*/
11
THD_dataset_mismatch(THD_3dim_dataset * ds1,THD_3dim_dataset * ds2)12 int THD_dataset_mismatch( THD_3dim_dataset *ds1 , THD_3dim_dataset *ds2 )
13 {
14 THD_dataxes * dax1 , * dax2 ;
15 THD_fvec3 fv1 , fv2 , dv ;
16 int code ;
17 float cd,c1,c2 ;
18 double angle;
19
20 ENTRY("THD_dataset_mismatch") ;
21
22 if( !ISVALID_DSET(ds1) || !ISVALID_DSET(ds2) ) RETURN(-1) ;
23
24 dax1 = ds1->daxes ;
25 dax2 = ds2->daxes ;
26 code = 0 ; /* will be return value */
27
28 /* check if the number of voxels in each direction is the same */
29
30 if( dax1->nxx != dax2->nxx ||
31 dax1->nyy != dax2->nyy ||
32 dax1->nzz != dax2->nzz ) code |= MISMATCH_DIMEN ;
33
34 /* check if the grid spacings are the same */
35
36 if( fabs(dax1->xxdel-dax2->xxdel) > 0.01*fabs(dax1->xxdel) ||
37 fabs(dax1->yydel-dax2->yydel) > 0.01*fabs(dax1->yydel) ||
38 fabs(dax1->zzdel-dax2->zzdel) > 0.01*fabs(dax1->zzdel) ) code |= MISMATCH_DELTA ;
39
40 /* check if the orientations are the same */
41
42 if( dax1->xxorient != dax2->xxorient ||
43 dax1->yyorient != dax2->yyorient ||
44 dax1->zzorient != dax2->zzorient ) code |= MISMATCH_ORIENT ;
45
46 /* check if they have the same centers */
47
48 fv1 = THD_dataset_center( ds1 ) ;
49 fv2 = THD_dataset_center( ds2 ) ;
50 dv = SUB_FVEC3(fv1,fv2) ; cd = SIZE_FVEC3(dv) ;
51
52 LOAD_FVEC3(fv1,dax1->xxdel,dax1->yydel,dax1->zzdel) ; c1 = SIZE_FVEC3(fv1) ;
53 LOAD_FVEC3(fv2,dax2->xxdel,dax2->yydel,dax2->zzdel) ; c2 = SIZE_FVEC3(fv2) ;
54
55 /* change default limit from 0.1 * vox diag sum to 0.001 * ave_vox_diag */
56 /* done to appease the mighty ZSS [9 Dec 2021 rickr] */
57 if( cd > 0.0005*(c1+c2) ) code |= MISMATCH_CENTER ;
58
59 /* check if the obliquity is the same */
60 /* (allow for truncation differnces) 22 May 2015 [rickr] */
61 angle = dset_obliquity_angle_diff(ds1, ds2, OBLIQ_ANGLE_THRESH);
62 if (angle > 0.0) code |= MISMATCH_OBLIQ ;
63
64 RETURN(code) ;
65 }
66
67 /*---------------------------------------------------------------------
68 23 Feb 2012: Return the absolute value of the difference between
69 two volumes, divided by the number of voxels
70 and the number of sub-bricks. Voxels that are zero
71 in both sets are not counted.
72 Comparisons are done after conversion of data to double
73 return = -1.0 ERROR
74 = 0.0 Exactly the same
75 -----------------------------------------------------------------------*/
THD_diff_vol_vals(THD_3dim_dataset * d1,THD_3dim_dataset * d2,int scl)76 double THD_diff_vol_vals(THD_3dim_dataset *d1, THD_3dim_dataset *d2, int scl) {
77 double dd=0.0, denom=0.0;
78 int i=0, k=0;
79 double *a1=NULL, *a2=NULL;
80 MRI_IMAGE *b1 = NULL , *b2 = NULL;
81
82 ENTRY("THD_diff_vol_vals");
83
84 if (!d1 && !d2) RETURN(dd);
85 if (!d1 || !d2) RETURN(-1.0);
86
87 if (!EQUIV_GRIDS(d1,d2)) RETURN(-1.0);
88 if (DSET_NVALS(d1) != DSET_NVALS(d2)) RETURN(-1.0);
89
90 DSET_mallocize(d1) ; DSET_load(d1) ;
91 DSET_mallocize(d2) ; DSET_load(d2) ;
92 dd = 0.0; denom = 0;
93 for (i=0; i<DSET_NVALS(d1); ++i) {
94 b1 = THD_extract_double_brick(i, d1);
95 b2 = THD_extract_double_brick(i, d2);
96 a1 = MRI_DOUBLE_PTR(b1);
97 a2 = MRI_DOUBLE_PTR(b2);
98 for (k=0; k<DSET_NVOX(d1); ++k) {
99 dd += ABS(a1[k]-a2[k]);
100 if (a1[k]!=0.0 || a2[k]!=0.0) ++denom;
101 }
102 mri_clear_data_pointer(b1); mri_free(b1) ;
103 mri_clear_data_pointer(b2); mri_free(b2) ;
104 }
105 if (scl && denom>0.0) dd /= denom;
106
107 RETURN(dd);
108 }
109