1 #include "mrilib.h"
2
3 /*---------------------------------------------------------------------------*/
4 /* Functions to deal with the mat44 elements in an AFNI dataset dataxes
5 structure. Note that the mat44 (4x4 matrix) type is defined in
6 nifti1_io.h, and that the last row of the matrix is ALWAYS [0 0 0 1].
7 It is to be applied to column vectors of the form [ x y z 1 ]', to
8 define an affine transformation of (x,y,z).
9
10 The library nifti1_io.c includes nifti1_mat44_inverse() and some
11 other utility functions.
12
13 -- RWCox - 07 Dec 2005
14 -----------------------------------------------------------------------------*/
15
16 /*---------------------------------------------------------------------------*/
17 /*! Multiply 2 mat44 matrices (a utility missing from nifti1_io.c).
18 -----------------------------------------------------------------------------*/
19
THD_mat44_mul(mat44 A,mat44 B)20 mat44 THD_mat44_mul( mat44 A , mat44 B )
21 {
22 mat44 C ; int i,j ;
23 for( i=0 ; i < 3 ; i++ )
24 for( j=0 ; j < 4 ; j++ )
25 C.m[i][j] = A.m[i][0] * B.m[0][j] + A.m[i][1] * B.m[1][j]
26 + A.m[i][2] * B.m[2][j] + A.m[i][3] * B.m[3][j] ;
27
28 C.m[3][0] = C.m[3][1] = C.m[3][2] = 0.0f ; C.m[3][3] = 1.0f ;
29 return C ;
30 }
31
32 /*---------------------------------------------------------------------------*/
33 /*! Compute the mat44 elements of a dataset dataxes struct, given the
34 other elements that define the coordinate system. Also see
35 function THD_daxes_from_mat44(), which does the reverse.
36
37 On input to this function, the XXdel, XXorg, nXX, and to_dicomm elements
38 of dax must have been set already, for XX=xx,yy,zz. Return value is
39 -1 if an error occurs, 0 if all is cool.
40 -----------------------------------------------------------------------------*/
41
THD_daxes_to_mat44(THD_dataxes * dax)42 int THD_daxes_to_mat44( THD_dataxes *dax )
43 {
44 mat44 ijk_to_dxyz , dxyz_to_dicom ;
45
46 if( dax == NULL ) return -1 ;
47
48 /* ijk_to_dxyz: transforms (i,j,k) to dataset (x,y,z) coords */
49
50 LOAD_MAT44( ijk_to_dxyz ,
51 dax->xxdel , 0.0f , 0.0f , dax->xxorg ,
52 0.0f , dax->yydel , 0.0f , dax->yyorg ,
53 0.0f , 0.0f , dax->zzdel , dax->zzorg ) ;
54
55 /* set to_dicomm (takes dataset order xyz to DICOM order xyz) */
56
57 THD_set_daxes_to_dicomm(dax) ;
58
59 /* dxyz_to_dicom: transforms dataset (x,y,z) coords to DICOM coords */
60
61 LOAD_MAT44( dxyz_to_dicom ,
62 dax->to_dicomm.mat[0][0] , dax->to_dicomm.mat[0][1] ,
63 dax->to_dicomm.mat[0][2] , 0.0f ,
64 dax->to_dicomm.mat[1][0] , dax->to_dicomm.mat[1][1] ,
65 dax->to_dicomm.mat[1][2] , 0.0f ,
66 dax->to_dicomm.mat[2][0] , dax->to_dicomm.mat[2][1] ,
67 dax->to_dicomm.mat[2][2] , 0.0f ) ;
68
69 /* dax->ijk_to_dicom: transforms (i,j,k) to DICOM (x,y,z) */
70
71 dax->ijk_to_dicom = THD_mat44_mul( dxyz_to_dicom , ijk_to_dxyz ) ;
72
73 /* and the inverse transformation: DICOM (x,y,z) to indexes (i,j,k) */
74
75 dax->dicom_to_ijk = nifti_mat44_inverse( dax->ijk_to_dicom ) ;
76
77 THD_set_dicom_box( dax ) ;
78
79 #if 0
80 { THD_dataxes new_daxes ;
81 new_daxes.type = DATAXES_TYPE;
82 THD_edit_dataxes(1.0, dax, &new_daxes);
83 printf("edited dataxes:\n"
84 " nxx=%d xxorg=%11.4g xxdel=%11.4g\n"
85 " nyy=%d yyorg=%11.4g yydel=%11.4g\n"
86 " nzz=%d zzorg=%11.4g zzdel=%11.4g\n" ,
87 new_daxes.nxx , new_daxes.xxorg , new_daxes.xxdel ,
88 new_daxes.nyy , new_daxes.yyorg , new_daxes.yydel ,
89 new_daxes.nzz , new_daxes.zzorg , new_daxes.zzdel ) ;
90 DUMP_MAT44("ijk_to_dicom",new_daxes.ijk_to_dicom) ;
91 }
92 #endif
93
94 return 0 ;
95 }
96
97 /*------------------------------------------------------------------------*/
98 /*! Set the min and max DICOM coords that can occur:
99 scan all 8 possible corners of the box the matrix defines.
100 --------------------------------------------------------------------------*/
101
THD_set_dicom_box(THD_dataxes * dax)102 void THD_set_dicom_box( THD_dataxes *dax )
103 {
104 float x,y,z , nx1,ny1,nz1 ;
105 float xbot,ybot,zbot , xtop,ytop,ztop ;
106
107 if( dax == NULL || !ISVALID_MAT44(dax->ijk_to_dicom) ) return ;
108
109 nx1 = dax->nxx - 1.0f; ny1 = dax->nyy - 1.0f; nz1 = dax->nzz - 1.0f;
110
111 MAT44_VEC(dax->ijk_to_dicom , 0,0,0 , x,y,z ) ; /* (0,0,0) corner */
112 xbot = xtop = x ; ybot = ytop = y ; zbot = ztop = z ;
113
114 /* this macro checks the (a,b,c) corner and updates the bot/top values */
115
116 #undef BT
117 #define BT(a,b,c) \
118 do{ MAT44_VEC(dax->ijk_to_dicom , a,b,c , x,y,z ) ; \
119 xbot = MIN(xbot,x); xtop = MAX(xtop,x) ; \
120 ybot = MIN(ybot,y); ytop = MAX(ytop,y) ; \
121 zbot = MIN(zbot,z); ztop = MAX(ztop,z) ; } while(0)
122
123 BT(nx1, 0 , 0 ) ; BT( 0 ,ny1, 0 ) ; BT(nx1,ny1, 0 ) ;
124 BT( 0 , 0 ,nz1) ; BT(nx1, 0 ,nz1) ; BT( 0 ,ny1,nz1) ; BT(nx1,ny1,nz1) ;
125
126 dax->dicom_xxmin = xbot ; dax->dicom_xxmax = xtop ;
127 dax->dicom_yymin = ybot ; dax->dicom_yymax = ytop ;
128 dax->dicom_zzmin = zbot ; dax->dicom_zzmax = ztop ;
129
130 #undef BT
131 return ;
132 }
133
134 /*---------------------------------------------------------------------------*/
135 /*! Given the ijk_to_dicom index to DICOM transformation in the header, AND
136 the nxx,nyy,nzz values, load the legacy dataxes information:
137 - xxorient = Orientation code
138 - yyorient = Orientation code
139 - zzorient = Orientation code
140 - xxorg = Center of (0,0,0) voxel
141 - yyorg = Center of (0,0,0) voxel
142 - zzorg = Center of (0,0,0) voxel
143 - xxdel = Spacings between voxel centers (mm) - may be negative
144 - yydel = Spacings between voxel centers (mm) - may be negative
145 - zzdel = Spacings between voxel centers (mm) - may be negative
146 - xxmin = Min value of x (etc.)
147 - xxmax = Min value of x (etc.)
148 - to_dicomm = Orthogonal matrix transforming from dataset coordinates
149 to DICOM coordinates
150
151 Return value is -1 if there is an error, 0 if not.
152 -----------------------------------------------------------------------------*/
153
THD_daxes_from_mat44(THD_dataxes * dax)154 int THD_daxes_from_mat44( THD_dataxes *dax )
155 {
156 int icod , jcod , kcod ;
157 mat44 nmat ;
158 float xx,yy,zz , ss , aa,bb,cc ;
159
160 /* table to convert NIfTI-1 orientation codes (1..6)
161 into AFNI orientation codes (0..5, and in a different order) */
162
163 static int orient_nifti2afni[7] =
164 { -1 , ORI_L2R_TYPE, ORI_R2L_TYPE, ORI_P2A_TYPE,
165 ORI_A2P_TYPE, ORI_I2S_TYPE, ORI_S2I_TYPE } ;
166
167 if( dax == NULL ) return -1 ;
168 if( dax->nxx < 1 || dax->nyy < 1 || dax->nzz < 1 ) return -1 ;
169
170 /* use the NIfTI-1 library function to determine best orientation;
171 but, must remember that NIfTI-1 x and y are reversed from AFNI's,
172 so we must negate the x and y rows of the matrix before func call */
173
174 nmat = dax->ijk_to_dicom ; XYINVERT_MAT44(nmat) ;
175
176 nifti_mat44_to_orientation( nmat , &icod, &jcod, &kcod ) ;
177
178 if( icod == 0 || jcod == 0 || kcod == 0 ) return -1 ;
179
180 dax->xxorient = orient_nifti2afni[icod] ;
181 dax->yyorient = orient_nifti2afni[jcod] ;
182 dax->zzorient = orient_nifti2afni[kcod] ;
183
184 /* grid offset int i-direction is projection of last
185 column of ijk_to_dicom matrix (the shifts) along the
186 direction of the first column of the matrix (the i-column) */
187
188 aa = dax->ijk_to_dicom.m[0][3] ; /* extract last column */
189 bb = dax->ijk_to_dicom.m[1][3] ;
190 cc = dax->ijk_to_dicom.m[2][3] ;
191
192 xx = dax->ijk_to_dicom.m[0][0] ; /* extract 1st column */
193 yy = dax->ijk_to_dicom.m[1][0] ;
194 zz = dax->ijk_to_dicom.m[2][0] ;
195 ss = sqrt(xx*xx+yy*yy+zz*zz) ; if( ss == 0.0f ) ss = 1.0f ;
196 dax->xxorg = (xx*aa+yy*bb+zz*cc) / ss ;
197 if( ORIENT_sign[dax->xxorient] == '-' ) dax->xxorg = -dax->xxorg ;
198
199 xx = dax->ijk_to_dicom.m[0][1] ; /* extract 2nd column */
200 yy = dax->ijk_to_dicom.m[1][1] ;
201 zz = dax->ijk_to_dicom.m[2][1] ;
202 ss = sqrt(xx*xx+yy*yy+zz*zz) ; if( ss == 0.0f ) ss = 1.0f ;
203 dax->yyorg = (xx*aa+yy*bb+zz*cc) / ss ;
204 if( ORIENT_sign[dax->yyorient] == '-' ) dax->yyorg = -dax->yyorg ;
205
206 xx = dax->ijk_to_dicom.m[0][2] ; /* extract 3rd column */
207 yy = dax->ijk_to_dicom.m[1][2] ;
208 zz = dax->ijk_to_dicom.m[2][2] ;
209 ss = sqrt(xx*xx+yy*yy+zz*zz) ; if( ss == 0.0f ) ss = 1.0f ;
210 dax->zzorg = (xx*aa+yy*bb+zz*cc) / ss ;
211 if( ORIENT_sign[dax->zzorient] == '-' ) dax->zzorg = -dax->zzorg ;
212
213 /* grid spacing along i-direction is length of 1st column of matrix */
214
215 dax->xxdel = sqrt( SQR(dax->ijk_to_dicom.m[0][0])
216 +SQR(dax->ijk_to_dicom.m[1][0])
217 +SQR(dax->ijk_to_dicom.m[2][0]) ) ;
218 if( ORIENT_sign[dax->xxorient] == '-' ) dax->xxdel = -dax->xxdel ;
219
220 /* mutatis mutandis for j- and k-directions */
221
222 dax->yydel = sqrt( SQR(dax->ijk_to_dicom.m[0][1])
223 +SQR(dax->ijk_to_dicom.m[1][1])
224 +SQR(dax->ijk_to_dicom.m[2][1]) ) ;
225 if( ORIENT_sign[dax->yyorient] == '-' ) dax->yydel = -dax->yydel ;
226
227 dax->zzdel = sqrt( SQR(dax->ijk_to_dicom.m[0][2])
228 +SQR(dax->ijk_to_dicom.m[1][2])
229 +SQR(dax->ijk_to_dicom.m[2][2]) ) ;
230 if( ORIENT_sign[dax->zzorient] == '-' ) dax->zzdel = -dax->zzdel ;
231
232 /* to_dicomm orthogonal matrix:
233 we make an orthogonal matrix out of the columns of ijk_to_dicom */
234
235 nmat = nifti_make_orthog_mat44(
236 dax->ijk_to_dicom.m[0][0], dax->ijk_to_dicom.m[1][0], dax->ijk_to_dicom.m[2][0],
237 dax->ijk_to_dicom.m[0][1], dax->ijk_to_dicom.m[1][1], dax->ijk_to_dicom.m[2][1],
238 dax->ijk_to_dicom.m[0][2], dax->ijk_to_dicom.m[1][2], dax->ijk_to_dicom.m[2][2] );
239
240 dax->to_dicomm.mat[0][0] = nmat.m[0][0] ;
241 dax->to_dicomm.mat[0][1] = nmat.m[1][0] ;
242 dax->to_dicomm.mat[0][2] = nmat.m[2][0] ;
243 dax->to_dicomm.mat[1][0] = nmat.m[0][1] ;
244 dax->to_dicomm.mat[1][1] = nmat.m[1][1] ;
245 dax->to_dicomm.mat[1][2] = nmat.m[2][1] ;
246 dax->to_dicomm.mat[2][0] = nmat.m[0][2] ;
247 dax->to_dicomm.mat[2][1] = nmat.m[1][2] ;
248 dax->to_dicomm.mat[2][2] = nmat.m[2][2] ;
249
250 /** min and max values of (x,y,z) **/
251
252 dax->xxmin = dax->xxorg ;
253 dax->xxmax = dax->xxorg + (dax->nxx-1) * dax->xxdel ;
254 if( dax->xxmin > dax->xxmax ){
255 float temp = dax->xxmin ;
256 dax->xxmin = dax->xxmax ; dax->xxmax = temp ;
257 }
258
259 dax->yymin = dax->yyorg ;
260 dax->yymax = dax->yyorg + (dax->nyy-1) * dax->yydel ;
261 if( dax->yymin > dax->yymax ){
262 float temp = dax->yymin ;
263 dax->yymin = dax->yymax ; dax->yymax = temp ;
264 }
265
266 dax->zzmin = dax->zzorg ;
267 dax->zzmax = dax->zzorg + (dax->nzz-1) * dax->zzdel ;
268 if( dax->zzmin > dax->zzmax ){
269 float temp = dax->zzmin ;
270 dax->zzmin = dax->zzmax ; dax->zzmax = temp ;
271 }
272
273 #ifdef EXTEND_BBOX
274 dax->xxmin -= 0.5 * fabsf(dax->xxdel) ; /* pushes edges back by 1/2 */
275 dax->xxmax += 0.5 * fabsf(dax->xxdel) ; /* voxel dimensions (the box */
276 dax->yymin -= 0.5 * fabsf(dax->yydel) ; /* defined above is based on */
277 dax->yymax += 0.5 * fabsf(dax->yydel) ; /* voxel centers, not edges) */
278 dax->zzmin -= 0.5 * fabsf(dax->zzdel) ;
279 dax->zzmax += 0.5 * fabsf(dax->zzdel) ;
280 #endif
281
282 return 0 ;
283 }
284
285 /*-----------------------------------------------------------------------*/
286 /*! Inputs:
287 - matrix old_mat: transforms (i,j,k) indexes to (x,y,z) coords
288 - ni,nj,nk: number of grid points along each direction
289 - ri,rj,rk: new grid spacing along each direction
290 Outputs:
291 - ninew,njnew,nknew: number of grid points along each new direction
292 - return value: new matrix; the [3][3] element will be 0 on error
293 -------------------------------------------------------------------------*/
294
THD_resample_mat44(mat44 old_mat,int ni,int nj,int nk,float ri,float rj,float rk,int * ninew,int * njnew,int * nknew)295 mat44 THD_resample_mat44( mat44 old_mat ,
296 int ni , int nj , int nk ,
297 float ri , float rj , float rk ,
298 int *ninew , int *njnew , int *nknew )
299 {
300 mat44 new_mat ;
301 float di,dj,dk , fi,fj,fk , aa,bb,cc , pp,qq,rr ;
302
303 ZZME(new_mat) ;
304 new_mat.m[3][3] = 0.0f ; /* mark output as bad; will fix later */
305
306 if( !ISVALID_MAT44(old_mat) ||
307 ninew == NULL || njnew == NULL || nknew == NULL ||
308 ni <= 0 || nj <= 0 || nk <= 0 ) return new_mat ;
309
310 /* spacing along input (i,j,k) directions */
311
312 di = MAT33_CLEN(old_mat,0) ; if( di == 0.0f ) di = 1.0f ;
313 dj = MAT33_CLEN(old_mat,1) ; if( dj == 0.0f ) dj = di ;
314 dk = MAT33_CLEN(old_mat,2) ; if( dk == 0.0f ) dk = di ;
315
316 if( ri <= 0.0f ) ri = 1.0f ; /* don't allow negative spacing */
317 if( rj <= 0.0f ) rj = ri ; /* along output (i,j,k) directions */
318 if( rk <= 0.0f ) rk = ri ;
319
320 fi = ri/di ; fj = rj/dj ; fk = rk/dk ;
321
322 /* copy input matrix to output,
323 then scale the upperleft 3x3 to allow for new spacings */
324
325 new_mat = old_mat ;
326
327 new_mat.m[0][0] *= fi ; /* first column */
328 new_mat.m[1][0] *= fi ;
329 new_mat.m[2][0] *= fi ;
330 new_mat.m[0][1] *= fj ; /* second column */
331 new_mat.m[1][1] *= fj ;
332 new_mat.m[2][1] *= fj ;
333 new_mat.m[0][2] *= fk ; /* third column */
334 new_mat.m[1][2] *= fk ;
335 new_mat.m[2][2] *= fk ;
336
337 /* number of points along each axis in new grid */
338
339 *ninew = (int)rint(ni/fi) ; /* so that ninew*ri == ni*di */
340 *njnew = (int)rint(nj/fj) ;
341 *nknew = (int)rint(nk/fk) ;
342
343 /* now compute offset (fourth column) in new_mat so
344 that the input and output grids have the same (x,y,z) centers;
345 center of input is at (i,j,k) = 0.5*(ni-1,nj-1,nk-1), et cetera */
346
347 di = 0.5*( ni -1) ; dj = 0.5*( nj -1) ; dk = 0.5*( nk -1) ;
348 fi = 0.5*(*ninew-1) ; fj = 0.5*(*njnew-1) ; fk = 0.5*(*nknew-1) ;
349
350 MAT33_VEC( old_mat , di,dj,dk , aa,bb,cc ) ;
351 MAT33_VEC( new_mat , fi,fj,fk , pp,qq,rr ) ;
352
353 /* adjust so that MAT44_VEC(new_mat,fi,fj,fk,...) is
354 same result as MAT44_VEC(old_mat,di,dj,dk,...) */
355
356 new_mat.m[0][3] += (aa-pp) ;
357 new_mat.m[1][3] += (bb-qq) ;
358 new_mat.m[2][3] += (cc-rr) ;
359
360 return new_mat ;
361 }
362