1 #include "mrilib.h"
2
main(int argc,char * argv[])3 int main( int argc , char *argv[] )
4 {
5 THD_3dim_dataset *dset=NULL ;
6 int nx,ny,nz,nxy,nxyz ;
7 int ii,jj,kk, mm, ip,jp,kp, im,jm,km, nmsk=0 , qq , verb=0 , iarg ;
8 byte *msk=NULL , *dsk=NULL ;
9 THD_ivec3 iv ;
10 THD_fvec3 *dispvec ; int num_dispvec=0 ;
11 MRI_IMAGE *matim=NULL; float *matim_far=NULL; int matim_nx=0,matim_ny=0;
12 mat44 qmat; float xx,yy,zz , dd,dmax,dbar ;
13
14 mainENTRY("3dmaxdisp") ; machdep() ;
15 PRINT_VERSION("3dmaxdisp"); AUTHOR("Zhark the Displacer");
16
17 /*-----------------------------------------------------------------------*/
18 if( argc < 2 || strcasecmp(argv[1],"-help") == 0 ){
19 printf("\n"
20 "Program 3dmaxdisp!\n"
21 "\n"
22 " * Reads in a 3D dataset and a DICOM-based affine matrix\n"
23 " * Outputs the average and maximum displacement that the matrix produces\n"
24 " when applied to the edge voxels of the 3D dataset's automask.\n"
25 " * The motivation for this program was to check if two\n"
26 " affine transformation matrices are 'close' -- but of course,\n"
27 " you can use this program for anything else you like.\n"
28 " * Suppose you have two affine transformation matrices that\n"
29 " transform a dataset XX.nii to MNI space, stored in files\n"
30 " AA.aff12.1D and BB.aff12.1D\n"
31 " and they aren't identical but they are close. How close?\n"
32 " * If these matrices are from 3dAllineate (-1Dmatrix_save),\n"
33 " then each matrix transforms DICOM-order coordinates\n"
34 " in XX.nii to MNI-space.\n"
35 " * So Inverse(AA) transforms MNI-space to XX-space\n"
36 " * So Inverse(AA)*BB transforms MNI-space to MNI-space,\n"
37 " going back to XX-space via matrix Inverse(AA) and then forward\n"
38 " to MNI-space via BB.\n"
39 " * This program (3dmaxdisp) can compute the average and maximum\n"
40 " displacement of Inverse(AA)*BB over the edges of the MNI template,\n"
41 " which will give you the answer to 'How close?' are the matrices.\n"
42 " If these displacements are on the order of a voxel size\n"
43 " (e.g., 1 mm), then the two matrices are for practical purposes\n"
44 " 'close enough' (in Zhark's opinion).\n"
45 " * How to do this?\n"
46 " cat_matvec AA.aff12.1D -I BB.aff12.1D > AinvB.aff12.1D\n"
47 " 3dmaxdisp -dset ~/abin/MNI152_2009_template_SSW.nii.gz'[0]' -matrix AinvB.aff12.1D\n"
48 " * Results are sent to stdout, two numbers per row (average and maximum),\n"
49 " one row of output for each matrix row given. Usually you will want to\n"
50 " capture stdout to a file with '>' or '| tee', depending on your further plans.\n"
51 "\n"
52 "-------\n"
53 "OPTIONS:\n"
54 "-------\n"
55 " -inset ddd }= The input dataset is 'ddd', which is used only to form\n"
56 " *OR* }= the mask over which the displacements will be computed.\n"
57 " -dset ddd }=\n"
58 "\n"
59 " -matrix mmm = File 'mmm' has 12 numbers per row, which are assembled\n"
60 " into the 3x4 affine transformation matrix to be applied\n"
61 " to the coordinates of the voxels in the dataset mask.\n"
62 " * As a special case, you can use the word 'IDENTITY'\n"
63 " for the matrix filename, which should result in\n"
64 " a max displacement of 0 mm.\n"
65 " * If there is more than 1 row in 'mmm', then each row\n"
66 " is treated as a separate matrix, and the max displacement\n"
67 " will be computed separately for each matrix.\n"
68 "\n"
69 " -verb = Print a few progress reports (to stderr).\n"
70 "\n"
71 "------\n"
72 "Author: Zhark the Displacer (AKA Bob the Inverted) -- June 2021\n"
73 "------\n"
74 "\n"
75 ) ;
76 exit(0) ;
77 }
78 /*-----------------------------------------------------------------------*/
79
80 /*---------- scan input options ----------*/
81
82 iarg = 1 ;
83 while( iarg < argc && argv[iarg][0] == '-' ){
84
85 if( strcmp(argv[iarg],"-dset") == 0 || strcmp(argv[iarg],"-inset") == 0 ){
86 if( dset != NULL ) ERROR_exit("Can't have multiple %s options :-(",argv[iarg]) ;
87 if( ++iarg >= argc ) ERROR_exit("no argument after '%s' :-(",argv[iarg-1]) ;
88 dset = THD_open_dataset( argv[iarg] ) ;
89 if( dset == NULL ) ERROR_exit("can't open input dataset '%s' :-(",argv[iarg]);
90 iarg++ ; continue ;
91 }
92
93 /*-----*/
94
95 if( strcmp(argv[iarg],"-matrix") == 0 ){
96 char *fname ; MRI_IMAGE *qim ;
97 if( matim != NULL ) ERROR_exit("Can't use '%s' more than once :(",argv[iarg]) ;
98 if( ++iarg >= argc ) ERROR_exit("no argument after '%s' :-(",argv[iarg-1]) ;
99 fname = argv[iarg] ;
100 if( strcasecmp(fname,"IDENTITY")==0 ) /* special case */
101 fname = "1D: 1 0 0 0 0 1 0 0 0 0 1 0" ;
102 qim = mri_read_1D(fname) ;
103 if( qim == NULL ) ERROR_exit("Can't read -matrix '%s' :-(",fname) ;
104 matim = mri_transpose(qim); mri_free(qim);
105 matim_far = MRI_FLOAT_PTR(matim) ;
106 matim_nx = matim->nx ; /* # of values per row */
107 matim_ny = matim->ny ; /* number of rows */
108 if( matim_nx < 12 && matim->nvox == 12 ){ /* special case of a 3x4 array */
109 matim_nx = 12 ; matim_ny = 1 ;
110 if(verb)INFO_message("-matrix: converting input 3x4 array to 1 row of 12 numbers") ;
111 }
112 if( matim_nx < 12 )
113 ERROR_exit("%d = Less than 12 numbers per row in -matrix '%s' :-(" ,matim_nx,fname) ;
114 else if( matim_nx > 12 )
115 WARNING_message("%d = More than 12 numbers per row in -matrix '%s'",matim_ny,fname) ;
116 iarg++ ; continue ;
117 }
118
119 /*-----*/
120
121 if( strcmp(argv[iarg],"-verb") == 0 ){
122 verb++ ; iarg++ ; continue ;
123 }
124
125 /*-----*/
126
127 ERROR_message("Unknown and Illegal option '%s' :-( :-( :-(",argv[iarg]) ;
128 suggest_best_prog_option(argv[0], argv[iarg]);
129 exit(1);
130
131 }
132
133 /*----- inputs without formal options? -----*/
134
135 if( matim == NULL )
136 ERROR_exit("No -matrix option ???") ;
137
138 if( dset == NULL && iarg >= argc )
139 ERROR_exit("No input dataset ???") ;
140
141 if( dset == NULL ){
142 dset = THD_open_dataset( argv[iarg] ) ;
143 if( dset == NULL ) ERROR_exit("can't open input dataset '%s' :-(",argv[iarg]);
144 iarg++ ;
145 }
146
147 nx = DSET_NX(dset) ;
148 ny = DSET_NY(dset) ;
149 nz = DSET_NZ(dset) ; nxy = nx*ny ; nxyz = nxy*nz ;
150
151 /*---------- create the xyz coords for automask edge voxels ----------*/
152
153 #undef DSK
154 #define DSK(i,j,k) dsk[(i)+(j)*nx+(k)*nxy]
155 if( verb ) INFO_message("creating automask") ;
156 THD_automask_set_clipfrac(0.333f) ;
157 DSET_load(dset) ;
158 dsk = THD_automask( dset ) ; /* volume automask */
159 if( dsk == NULL )
160 ERROR_exit("Cannot create automask for input dataset :(") ;
161
162 if( verb ) ININFO_message(" %s voxels in automask",
163 commaized_integer_string(mask_count(nxyz,dsk)) ) ;
164
165 msk = (byte *)calloc(1,nxyz) ;
166 for( nmsk=mm=0 ; mm < nxyz ; mm++ ){ /* create edges of the automask */
167 if( dsk[mm] == 0 ) continue ; /* not in the mask */
168 ii = mm % nx ; kk = mm / nxy ; jj = (mm%nxy) / nx ; /* voxel indexes */
169 ip=ii+1; im=ii-1; if(ip>=nx || im<0){ msk[mm]=1; nmsk++; continue; }
170 jp=jj+1; jm=jj-1; if(jp>=ny || jm<0){ msk[mm]=1; nmsk++; continue; }
171 kp=kk+1; km=kk-1; if(kp>=nz || km<0){ msk[mm]=1; nmsk++; continue; }
172 if( DSK(ip,jj,kk) && DSK(im,jj,kk) && /* if all 6 nbhrs */
173 DSK(ii,jp,kk) && DSK(ii,jm,kk) && /* are in automask */
174 DSK(ii,jj,kp) && DSK(ii,jj,km) ) continue ; /* skip this voxel */
175 msk[mm] = 1 ; nmsk++ ;
176 }
177 free(dsk) ;
178 if( nmsk == 0 )
179 ERROR_exit("Cannot find edge voxels for input dataset automask :(") ;
180 dispvec = (THD_fvec3 *)malloc(sizeof(THD_fvec3)*nmsk) ;
181 for( qq=mm=0 ; mm < nxyz ; mm++ ){
182 if( msk[mm] == 0 ) continue ;
183 ii = mm % nx ; kk = mm / nxy ; jj = (mm%nxy) / nx ;
184 iv.ijk[0] = ii ; iv.ijk[1] = jj ; iv.ijk[2] = kk ; /* convert 3D */
185 dispvec[qq++] = THD_3dind_to_dicomm_no_wod( dset , iv ) ; /* index to xyz */
186 }
187 free(msk) ;
188 num_dispvec = qq ;
189 if( verb ) ININFO_message(" %s edge voxels in automask",commaized_integer_string(num_dispvec)) ;
190
191 /*--------- Loop over matrix rows and compute displacments ----------*/
192
193 for( qq=0 ; qq < matim_ny ; qq++ ){
194 LOAD_MAT44_AR( qmat , matim_far+(12*qq) ) ; /* get matrix */
195 qmat.m[0][0] -= 1.0f ; /* subtract identity */
196 qmat.m[1][1] -= 1.0f ; /* matrix */
197 qmat.m[2][2] -= 1.0f ;
198
199 for( dbar=dmax=0.0f,mm=0 ; mm < num_dispvec ; mm++ ){ /* loop over voxels */
200 MAT44_VEC( qmat,
201 dispvec[mm].xyz[0], dispvec[mm].xyz[1], dispvec[mm].xyz[2],
202 xx , yy , zz ) ;
203 dd = sqrtf(xx*xx + yy*yy + zz*zz) ; /* length of this displacement */
204 if( dd > dmax ) dmax = dd ;
205 dbar += dd ;
206 }
207 dbar /= num_dispvec ;
208 printf(" %.6g %.6g\n",dbar,dmax) ;
209 }
210
211 exit(0) ;
212 }
213