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 "mrilib.h"
8 #include "thd.h"
9
10 static int oblique_report_index = 0;
11 static int oblique_report_repeat = 20;
12 static int oblique_report_repeat2 = 100;
13 static int first_oblique = 1;
14 static int oblique_update = 0;
15 static int OBL_report=1;
set_obliquity_report(int v)16 void set_obliquity_report(int v) { OBL_report=v; }
17
18
19 /*====================================================================
20 3D coordinate conversion routines;
21 tags for coordinate systems:
22 fdind = FD_brick voxel indices (ints)
23 fdfind = FD_brick voxel indices (floats - added 30 Aug 2001)
24 3dind = THD_3dim_dataset voxel indices (ints)
25 3dfind = THD_3dim_dataset floating point voxel indices
26 (used for subvoxel resolution)
27 3dmm = THD_3dim_dataset millimetric coordinates (floats)
28 dicomm = DICOM 3.0 millimetric coordinates (floats)
29
30 The 3dmm coordinate measurements are oriented the same way
31 as the dicomm coordinates (which means that the ??del values
32 can be negative if the voxel axes are reflected), but the 3dmm
33 coordinates are in general a permutation of the DICOM coordinates.
34
35 These routines all take as input and produce as output
36 THD_?vec3 (i=int,f=float) structs.
37 ======================================================================*/
38
THD_3dfind_to_3dmm(THD_3dim_dataset * dset,THD_fvec3 iv)39 THD_fvec3 THD_3dfind_to_3dmm( THD_3dim_dataset *dset , THD_fvec3 iv )
40 {
41 THD_dataxes *daxes ;
42 THD_fvec3 fv ;
43
44 daxes = CURRENT_DAXES(dset) ;
45
46 fv.xyz[0] = daxes->xxorg + iv.xyz[0] * daxes->xxdel ;
47 fv.xyz[1] = daxes->yyorg + iv.xyz[1] * daxes->yydel ;
48 fv.xyz[2] = daxes->zzorg + iv.xyz[2] * daxes->zzdel ;
49 return fv ;
50 }
51
52 /*------------------------------------------------------------------*/
53
THD_3dind_to_3dmm(THD_3dim_dataset * dset,THD_ivec3 iv)54 THD_fvec3 THD_3dind_to_3dmm( THD_3dim_dataset *dset , THD_ivec3 iv )
55 {
56 THD_dataxes *daxes ;
57 THD_fvec3 fv ;
58
59 daxes = CURRENT_DAXES(dset) ;
60
61 fv.xyz[0] = daxes->xxorg + iv.ijk[0] * daxes->xxdel ;
62 fv.xyz[1] = daxes->yyorg + iv.ijk[1] * daxes->yydel ;
63 fv.xyz[2] = daxes->zzorg + iv.ijk[2] * daxes->zzdel ;
64 return fv ;
65 }
66
67 /*--------------------------------------------------------------------*/
68 /* this version is without using wod dataxes 7 Mar 2005 [rickr] */
69
THD_3dind_to_3dmm_no_wod(THD_3dim_dataset * dset,THD_ivec3 iv)70 THD_fvec3 THD_3dind_to_3dmm_no_wod( THD_3dim_dataset *dset , THD_ivec3 iv )
71 {
72 THD_dataxes *daxes ;
73 THD_fvec3 fv ;
74
75 daxes = dset->daxes ;
76
77 fv.xyz[0] = daxes->xxorg + iv.ijk[0] * daxes->xxdel ;
78 fv.xyz[1] = daxes->yyorg + iv.ijk[1] * daxes->yydel ;
79 fv.xyz[2] = daxes->zzorg + iv.ijk[2] * daxes->zzdel ;
80 return fv ;
81 }
82
83 /*--------------------------------------------------------------------*/
84 /* 04 Apr 2012 - RWCox */
85
THD_3dind_to_dicomm_no_wod(THD_3dim_dataset * dset,THD_ivec3 iv)86 THD_fvec3 THD_3dind_to_dicomm_no_wod( THD_3dim_dataset *dset , THD_ivec3 iv )
87 {
88 THD_fvec3 fv , gv ;
89 fv = THD_3dind_to_3dmm_no_wod( dset , iv ) ;
90 gv = THD_3dmm_to_dicomm ( dset , fv ) ;
91 return gv ;
92 }
93
94 /*--------------------------------------------------------------------*/
95
THD_3dmm_to_3dfind(THD_3dim_dataset * dset,THD_fvec3 fv)96 THD_fvec3 THD_3dmm_to_3dfind( THD_3dim_dataset *dset , THD_fvec3 fv )
97 {
98 THD_dataxes *daxes ;
99 THD_fvec3 iv ;
100
101 daxes = CURRENT_DAXES(dset) ;
102
103 iv.xyz[0] = (fv.xyz[0] - daxes->xxorg) / daxes->xxdel ;
104 iv.xyz[1] = (fv.xyz[1] - daxes->yyorg) / daxes->yydel ;
105 iv.xyz[2] = (fv.xyz[2] - daxes->zzorg) / daxes->zzdel ;
106
107 if( iv.xyz[0] < 0 ) iv.xyz[0] = 0 ;
108 else if( iv.xyz[0] > daxes->nxx-1 ) iv.xyz[0] = daxes->nxx-1 ;
109
110 if( iv.xyz[1] < 0 ) iv.xyz[1] = 0 ;
111 else if( iv.xyz[1] > daxes->nyy-1 ) iv.xyz[1] = daxes->nyy-1 ;
112
113 if( iv.xyz[2] < 0 ) iv.xyz[2] = 0 ;
114 else if( iv.xyz[2] > daxes->nzz-1 ) iv.xyz[2] = daxes->nzz-1 ;
115
116 return iv ;
117 }
118
119 /*--------------------------------------------------------------------*/
120
THD_3dmm_to_3dind_warn(THD_3dim_dataset * dset,THD_fvec3 fv,int * out)121 THD_ivec3 THD_3dmm_to_3dind_warn( THD_3dim_dataset* dset ,
122 THD_fvec3 fv, int *out )
123 {
124 THD_dataxes *daxes ;
125 THD_ivec3 iv ;
126
127 *out = 0;
128 daxes = CURRENT_DAXES(dset) ;
129
130 iv.ijk[0] = (fv.xyz[0] - daxes->xxorg) / daxes->xxdel + 0.49f ;
131 iv.ijk[1] = (fv.xyz[1] - daxes->yyorg) / daxes->yydel + 0.49f ;
132 iv.ijk[2] = (fv.xyz[2] - daxes->zzorg) / daxes->zzdel + 0.49f ;
133
134 if( iv.ijk[0] < 0 ) { iv.ijk[0] = 0 ; *out = 1; }
135 else if( iv.ijk[0] > daxes->nxx-1 ) { iv.ijk[0] = daxes->nxx-1 ; *out = 1; }
136
137 if( iv.ijk[1] < 0 ) { iv.ijk[1] = 0 ; *out = 1; }
138 else if( iv.ijk[1] > daxes->nyy-1 ) { iv.ijk[1] = daxes->nyy-1 ; *out = 1; }
139
140 if( iv.ijk[2] < 0 ) { iv.ijk[2] = 0 ; *out = 1; }
141 else if( iv.ijk[2] > daxes->nzz-1 ) { iv.ijk[2] = daxes->nzz-1 ; *out = 1; }
142
143 return iv ;
144 }
145
THD_3dmm_to_3dind(THD_3dim_dataset * dset,THD_fvec3 fv)146 THD_ivec3 THD_3dmm_to_3dind( THD_3dim_dataset *dset , THD_fvec3 fv )
147 {
148 THD_dataxes *daxes ;
149 THD_ivec3 iv ;
150
151 daxes = CURRENT_DAXES(dset) ;
152
153 iv.ijk[0] = (fv.xyz[0] - daxes->xxorg) / daxes->xxdel + 0.49f ;
154 iv.ijk[1] = (fv.xyz[1] - daxes->yyorg) / daxes->yydel + 0.49f ;
155 iv.ijk[2] = (fv.xyz[2] - daxes->zzorg) / daxes->zzdel + 0.49f ;
156
157 if( iv.ijk[0] < 0 ) iv.ijk[0] = 0 ;
158 else if( iv.ijk[0] > daxes->nxx-1 ) iv.ijk[0] = daxes->nxx-1 ;
159
160 if( iv.ijk[1] < 0 ) iv.ijk[1] = 0 ;
161 else if( iv.ijk[1] > daxes->nyy-1 ) iv.ijk[1] = daxes->nyy-1 ;
162
163 if( iv.ijk[2] < 0 ) iv.ijk[2] = 0 ;
164 else if( iv.ijk[2] > daxes->nzz-1 ) iv.ijk[2] = daxes->nzz-1 ;
165
166 #if defined(USE_TRACING) && 0
167 INFO_message("THD_3dmm_to_3dind: fv=%f %f %f iv=%d %d %d from %s",
168 fv.xyz[0],fv.xyz[1],fv.xyz[2] , iv.ijk[0],iv.ijk[1],iv.ijk[2] , DBROUT ) ;
169 #endif
170
171 return iv ;
172 }
173
174 /*--------------------------------------------------------------------*/
175
176 /* this version is without using wod dataxes 28 Sep 2004 [rickr] */
177
THD_3dmm_to_3dind_no_wod(THD_3dim_dataset * dset,THD_fvec3 fv)178 THD_ivec3 THD_3dmm_to_3dind_no_wod( THD_3dim_dataset *dset , THD_fvec3 fv )
179 {
180 THD_dataxes *daxes ;
181 THD_ivec3 iv ;
182
183 daxes = dset->daxes ;
184
185 iv.ijk[0] = (fv.xyz[0] - daxes->xxorg) / daxes->xxdel + 0.49f ;
186 iv.ijk[1] = (fv.xyz[1] - daxes->yyorg) / daxes->yydel + 0.49f ;
187 iv.ijk[2] = (fv.xyz[2] - daxes->zzorg) / daxes->zzdel + 0.49f ;
188
189 if( iv.ijk[0] < 0 ) iv.ijk[0] = 0 ;
190 else if( iv.ijk[0] > daxes->nxx-1 ) iv.ijk[0] = daxes->nxx-1 ;
191
192 if( iv.ijk[1] < 0 ) iv.ijk[1] = 0 ;
193 else if( iv.ijk[1] > daxes->nyy-1 ) iv.ijk[1] = daxes->nyy-1 ;
194
195 if( iv.ijk[2] < 0 ) iv.ijk[2] = 0 ;
196 else if( iv.ijk[2] > daxes->nzz-1 ) iv.ijk[2] = daxes->nzz-1 ;
197
198 #if defined(USE_TRACING) && 0
199 INFO_message("THD_3dmm_to_3dind_no_wod: fv=%f %f %f iv=%d %d %d from %s",
200 fv.xyz[0],fv.xyz[1],fv.xyz[2] , iv.ijk[0],iv.ijk[1],iv.ijk[2] , DBROUT ) ;
201 #endif
202
203 return iv ;
204 }
205
206 /*---------------------------------------------------------------------
207 convert from input image oriented x,y,z to Dicom x,y,z
208 (x axis = R->L , y axis = A->P , z axis = I->S)
209
210 N.B.: image distances are oriented the same as Dicom,
211 just in a permuted order.
212 -----------------------------------------------------------------------*/
213
THD_3dmm_to_dicomm(THD_3dim_dataset * dset,THD_fvec3 imv)214 THD_fvec3 THD_3dmm_to_dicomm( THD_3dim_dataset *dset , THD_fvec3 imv )
215 {
216 THD_fvec3 dicv ;
217 float xim,yim,zim , xdic=0,ydic=0,zdic=0 ;
218
219 xim = imv.xyz[0] ; yim = imv.xyz[1] ; zim = imv.xyz[2] ;
220
221 switch( dset->daxes->xxorient ){
222 case ORI_R2L_TYPE:
223 case ORI_L2R_TYPE: xdic = xim ; break ;
224 case ORI_P2A_TYPE:
225 case ORI_A2P_TYPE: ydic = xim ; break ;
226 case ORI_I2S_TYPE:
227 case ORI_S2I_TYPE: zdic = xim ; break ;
228
229 default: THD_FATAL_ERROR("illegal xxorient code") ;
230 }
231
232 switch( dset->daxes->yyorient ){
233 case ORI_R2L_TYPE:
234 case ORI_L2R_TYPE: xdic = yim ; break ;
235 case ORI_P2A_TYPE:
236 case ORI_A2P_TYPE: ydic = yim ; break ;
237 case ORI_I2S_TYPE:
238 case ORI_S2I_TYPE: zdic = yim ; break ;
239
240 default: THD_FATAL_ERROR("illegal yyorient code") ;
241 }
242
243 switch( dset->daxes->zzorient ){
244 case ORI_R2L_TYPE:
245 case ORI_L2R_TYPE: xdic = zim ; break ;
246 case ORI_P2A_TYPE:
247 case ORI_A2P_TYPE: ydic = zim ; break ;
248 case ORI_I2S_TYPE:
249 case ORI_S2I_TYPE: zdic = zim ; break ;
250
251 default: THD_FATAL_ERROR("illegal zzorient code") ;
252 }
253
254 dicv.xyz[0] = xdic ; dicv.xyz[1] = ydic ; dicv.xyz[2] = zdic ;
255 return dicv ;
256 }
257
258 /*---------------------------------------------------------------------
259 convert to input image oriented x,y,z from Dicom x,y,z
260 -----------------------------------------------------------------------*/
261
THD_dicomm_to_3dmm(THD_3dim_dataset * dset,THD_fvec3 dicv)262 THD_fvec3 THD_dicomm_to_3dmm( THD_3dim_dataset *dset , THD_fvec3 dicv )
263 {
264 THD_fvec3 imv ;
265 float xim,yim,zim , xdic,ydic,zdic ;
266
267 xdic = dicv.xyz[0] ; ydic = dicv.xyz[1] ; zdic = dicv.xyz[2] ;
268
269 switch( dset->daxes->xxorient ){
270 case ORI_R2L_TYPE:
271 case ORI_L2R_TYPE: xim = xdic ; break ;
272 case ORI_P2A_TYPE:
273 case ORI_A2P_TYPE: xim = ydic ; break ;
274 case ORI_I2S_TYPE:
275 case ORI_S2I_TYPE: xim = zdic ; break ;
276
277 default: THD_FATAL_ERROR("illegal xxorient code") ;
278 }
279
280 switch( dset->daxes->yyorient ){
281 case ORI_R2L_TYPE:
282 case ORI_L2R_TYPE: yim = xdic ; break ;
283 case ORI_P2A_TYPE:
284 case ORI_A2P_TYPE: yim = ydic ; break ;
285 case ORI_I2S_TYPE:
286 case ORI_S2I_TYPE: yim = zdic ; break ;
287
288 default: THD_FATAL_ERROR("illegal yyorient code") ;
289 }
290
291 switch( dset->daxes->zzorient ){
292 case ORI_R2L_TYPE:
293 case ORI_L2R_TYPE: zim = xdic ; break ;
294 case ORI_P2A_TYPE:
295 case ORI_A2P_TYPE: zim = ydic ; break ;
296 case ORI_I2S_TYPE:
297 case ORI_S2I_TYPE: zim = zdic ; break ;
298
299 default: THD_FATAL_ERROR("illegal zzorient code") ;
300 }
301
302 imv.xyz[0] = xim ; imv.xyz[1] = yim ; imv.xyz[2] = zim ;
303 return imv ;
304 }
305
306 /*---------------------------------------------------------------------*/
307
THD_fdind_to_3dind(FD_brick * br,THD_ivec3 ib)308 THD_ivec3 THD_fdind_to_3dind( FD_brick *br , THD_ivec3 ib )
309 {
310 THD_ivec3 id ;
311 int qq , ax ;
312
313 for( qq=0 ; qq < 3 ; qq++ ){
314 ax = abs( br->a123.ijk[qq] ) - 1 ; /* 0,1,2, for x,y,z */
315
316 if( br->a123.ijk[qq] > 0 ) id.ijk[ax] = ib.ijk[qq] ;
317 else id.ijk[ax] = br->sxyz.ijk[ax] - ib.ijk[qq];
318 }
319
320 return id ;
321 }
322
THD_3dind_to_fdind(FD_brick * br,THD_ivec3 id)323 THD_ivec3 THD_3dind_to_fdind( FD_brick *br , THD_ivec3 id )
324 {
325 THD_ivec3 ib ;
326 int qq , ax ;
327
328 for( qq=0 ; qq < 3 ; qq++ ){
329 ax = abs( br->a123.ijk[qq] ) - 1 ;
330
331 if( br->a123.ijk[qq] > 0 ) ib.ijk[qq] = id.ijk[ax] ;
332 else ib.ijk[qq] = br->sxyz.ijk[ax] - id.ijk[ax];
333 }
334
335 return ib ;
336 }
337
338 /*---------------------------------------------------------------------*/
339
THD_fdfind_to_3dfind(FD_brick * br,THD_fvec3 ib)340 THD_fvec3 THD_fdfind_to_3dfind( FD_brick *br , THD_fvec3 ib ) /* 30 Aug 2001 */
341 {
342 THD_fvec3 id ;
343 int qq , ax ;
344
345 for( qq=0 ; qq < 3 ; qq++ ){
346 ax = abs( br->a123.ijk[qq] ) - 1 ; /* 0,1,2, for x,y,z */
347
348 if( br->a123.ijk[qq] > 0 ) id.xyz[ax] = ib.xyz[qq] ;
349 else id.xyz[ax] = br->sxyz.ijk[ax] - ib.xyz[qq];
350 }
351
352 return id ;
353 }
354
THD_3dfind_to_fdfind(FD_brick * br,THD_fvec3 id)355 THD_fvec3 THD_3dfind_to_fdfind( FD_brick *br , THD_fvec3 id ) /* 30 Aug 2001 */
356 {
357 THD_fvec3 ib ;
358 int qq , ax ;
359
360 for( qq=0 ; qq < 3 ; qq++ ){
361 ax = abs( br->a123.ijk[qq] ) - 1 ;
362
363 if( br->a123.ijk[qq] > 0 ) ib.xyz[qq] = id.xyz[ax] ;
364 else ib.xyz[qq] = br->sxyz.ijk[ax] - id.xyz[ax];
365 }
366
367 return ib ;
368 }
369
370 /*********************************************************************
371 Convenience is bliss. Coordinate inter-conversion routines.
372 **********************************************************************/
373
AFNI_ijk_to_xyz(THD_3dim_dataset * dset,int ii,int jj,int kk,float * xx,float * yy,float * zz)374 void AFNI_ijk_to_xyz( THD_3dim_dataset * dset ,
375 int ii , int jj , int kk ,
376 float * xx , float * yy , float * zz )
377 {
378 THD_fvec3 fv ;
379
380 if( ! ISVALID_DSET(dset) ) return ;
381
382 fv = THD_3dind_to_3dmm( dset , TEMP_IVEC3(ii,jj,kk) ) ;
383 *xx = fv.xyz[0] ;
384 *yy = fv.xyz[1] ;
385 *zz = fv.xyz[2] ;
386 return ;
387 }
388
AFNI_xyz_to_ijk(THD_3dim_dataset * dset,float xx,float yy,float zz,int * ii,int * jj,int * kk)389 void AFNI_xyz_to_ijk( THD_3dim_dataset * dset ,
390 float xx , float yy , float zz ,
391 int * ii , int * jj , int * kk )
392 {
393 THD_ivec3 iv ;
394
395 if( ! ISVALID_DSET(dset) ) return ;
396
397 iv = THD_3dmm_to_3dind ( dset , TEMP_FVEC3(xx,yy,zz) ) ;
398 *ii = iv.ijk[0] ;
399 *jj = iv.ijk[1] ;
400 *kk = iv.ijk[2] ;
401 return ;
402 }
403
AFNI_xyz_to_dicomm(THD_3dim_dataset * dset,float xx,float yy,float zz,float * xd,float * yd,float * zd)404 void AFNI_xyz_to_dicomm( THD_3dim_dataset * dset ,
405 float xx , float yy , float zz ,
406 float * xd , float * yd , float * zd )
407 {
408 THD_fvec3 fv ;
409
410 if( ! ISVALID_DSET(dset) ) return ;
411
412 fv = THD_3dmm_to_dicomm( dset , TEMP_FVEC3(xx,yy,zz) ) ;
413 *xd = fv.xyz[0] ;
414 *yd = fv.xyz[1] ;
415 *zd = fv.xyz[2] ;
416 return ;
417 }
418
AFNI_dicomm_to_xyz(THD_3dim_dataset * dset,float xd,float yd,float zd,float * xx,float * yy,float * zz)419 void AFNI_dicomm_to_xyz( THD_3dim_dataset * dset ,
420 float xd , float yd , float zd ,
421 float * xx , float * yy , float * zz )
422 {
423 THD_fvec3 fv ;
424
425 if( ! ISVALID_DSET(dset) ) return ;
426
427 fv = THD_3dmm_to_dicomm( dset , TEMP_FVEC3(xd,yd,zd) ) ;
428 *xx = fv.xyz[0] ;
429 *yy = fv.xyz[1] ;
430 *zz = fv.xyz[2] ;
431 return ;
432 }
433
434 /*-------------------------------------------------------------------*/
435
THD_coorder_fill(char * in_orcode,THD_coorder * cord)436 void THD_coorder_fill( char *in_orcode , THD_coorder *cord )
437 {
438 char acod , orcode[4] ;
439 int xx,yy,zz , ss1,ss2,ss3 , ii,ll ;
440
441 if( cord == NULL ) return ;
442
443 /* default values */
444
445 cord->xxsign = cord->yysign = cord->zzsign = 1 ;
446 cord->first = 0 ;
447 cord->second = 1 ;
448 cord->third = 2 ;
449 cord->xxor = ORI_R2L_TYPE ;
450 cord->yyor = ORI_A2P_TYPE ;
451 cord->zzor = ORI_I2S_TYPE ;
452 strcpy(cord->orcode,"RAI") ;
453
454 /* check string for OKness */
455
456 if( in_orcode == NULL ) return ;
457 strncpy(orcode,in_orcode,3) ; orcode[3] = '\0' ;
458 ll = strlen(orcode) ; if( ll != 3 ) return ;
459 for( ii=0 ; ii < 3 ; ii++ ) orcode[ii] = toupper(orcode[ii]) ;
460 if( strncmp(orcode,"FLI",3) == 0 ) strcpy(orcode,"LPI") ;
461
462 /* extract direction codes */
463
464 acod = orcode[0] ; xx = ORCODE(acod) ;
465 acod = orcode[1] ; yy = ORCODE(acod) ;
466 acod = orcode[2] ; zz = ORCODE(acod) ;
467
468 /* check direction codes for OKness */
469
470 if( xx<0 || yy<0 || zz<0 || ! OR3OK(xx,yy,zz) ) return ;
471
472 /* all is OK. get signs of orientations */
473
474 ss1 = (ORIENT_sign[xx] == '-') ? -1 : 1 ;
475 ss2 = (ORIENT_sign[yy] == '-') ? -1 : 1 ;
476 ss3 = (ORIENT_sign[zz] == '-') ? -1 : 1 ;
477
478 /* whose on first? */
479
480 cord->first = xx / 2 ;
481 cord->second = yy / 2 ;
482 cord->third = zz / 2 ;
483
484 cord->xxsign = (cord->first ==0) ? ss1
485 :(cord->second==0) ? ss2 : ss3 ;
486
487 cord->yysign = (cord->first ==1) ? ss1
488 :(cord->second==1) ? ss2 : ss3 ;
489
490 cord->zzsign = (cord->first ==2) ? ss1
491 :(cord->second==2) ? ss2 : ss3 ;
492
493 cord->xxor = xx ;
494 cord->yyor = yy ;
495 cord->zzor = zz ;
496
497 strcpy(cord->orcode,orcode) ;
498 return ;
499 }
500
501 /*---------------------------------------------------------------------
502 convert to output order x,y,z from Dicom x,y,z
503 -----------------------------------------------------------------------*/
504
THD_dicom_to_coorder(THD_coorder * cord,float * xx,float * yy,float * zz)505 void THD_dicom_to_coorder( THD_coorder *cord ,
506 float *xx , float *yy , float *zz )
507 {
508 float xval , yval , zval ;
509
510 if( cord == NULL ) return ;
511
512 /* changes signs first */
513
514 xval = cord->xxsign * (*xx) ;
515 yval = cord->yysign * (*yy) ;
516 zval = cord->zzsign * (*zz) ;
517
518 /* scramble order */
519
520 *xx = (cord->first == 0) ? xval
521 :(cord->first == 1) ? yval : zval ;
522
523 *yy = (cord->second == 0) ? xval
524 :(cord->second == 1) ? yval : zval ;
525
526 *zz = (cord->third == 0) ? xval
527 :(cord->third == 1) ? yval : zval ;
528
529 return ;
530 }
531
532 /*-------------------------------------------------------------------
533 convert to Dicom x,y,z from output order x,y,z
534 ---------------------------------------------------------------------*/
535
THD_coorder_to_dicom(THD_coorder * cord,float * xx,float * yy,float * zz)536 void THD_coorder_to_dicom( THD_coorder *cord ,
537 float *xx , float *yy , float *zz )
538 {
539 float xval , yval , zval ;
540
541 if( cord == NULL ) return ;
542
543 /* unscramble order */
544
545 xval = (cord->first == 0) ? (*xx)
546 :(cord->second == 0) ? (*yy) : (*zz) ;
547
548 yval = (cord->first == 1) ? (*xx)
549 :(cord->second == 1) ? (*yy) : (*zz) ;
550
551 zval = (cord->first == 2) ? (*xx)
552 :(cord->second == 2) ? (*yy) : (*zz) ;
553
554 /* change signs */
555
556 *xx = cord->xxsign * xval ;
557 *yy = cord->yysign * yval ;
558 *zz = cord->zzsign * zval ;
559
560 return ;
561 }
562
563 /*---------------------------------------------------------------------
564 Return rotation and shift param. to go from i, j, k to Cardinal
565 Dicom x,y,z (non-oblique)
566 -----------------------------------------------------------------------*/
THD_dicom_card_xform(THD_3dim_dataset * dset,THD_dmat33 * tmat,THD_dfvec3 * dics)567 void THD_dicom_card_xform (THD_3dim_dataset *dset ,
568 THD_dmat33 *tmat, THD_dfvec3 *dics )
569 {
570
571 THD_dfvec3 dicr;
572
573 /* rotation business */
574 switch( dset->daxes->xxorient ){
575 case ORI_R2L_TYPE:
576 case ORI_L2R_TYPE:
577 tmat->mat[0][0] = dset->daxes->xxdel ; tmat->mat[0][1] = tmat->mat[0][2] = 0.0;
578 dics->xyz[0] = dset->daxes->xxorg;
579 break ;
580 case ORI_P2A_TYPE:
581 case ORI_A2P_TYPE:
582 tmat->mat[1][0] = dset->daxes->xxdel ; tmat->mat[1][1] = tmat->mat[1][2] = 0.0;
583 dics->xyz[1] = dset->daxes->xxorg;
584 break ;
585 case ORI_I2S_TYPE:
586 case ORI_S2I_TYPE:
587 tmat->mat[2][0] = dset->daxes->xxdel ; tmat->mat[2][1] = tmat->mat[2][2] = 0.0;
588 dics->xyz[2] = dset->daxes->xxorg;
589 break ;
590 default: THD_FATAL_ERROR("illegal xxorient code") ;
591 }
592
593 switch( dset->daxes->yyorient ){
594 case ORI_R2L_TYPE:
595 case ORI_L2R_TYPE:
596 tmat->mat[0][1] = dset->daxes->yydel ; tmat->mat[0][0] = tmat->mat[0][2] = 0.0;
597 dics->xyz[0] = dset->daxes->yyorg;
598 break ;
599 case ORI_P2A_TYPE:
600 case ORI_A2P_TYPE:
601 tmat->mat[1][1] = dset->daxes->yydel ; tmat->mat[1][0] = tmat->mat[1][2] = 0.0;
602 dics->xyz[1] = dset->daxes->yyorg;
603 break ;
604 case ORI_I2S_TYPE:
605 case ORI_S2I_TYPE:
606 tmat->mat[2][1] = dset->daxes->yydel ; tmat->mat[2][0] = tmat->mat[2][2] = 0.0;
607 dics->xyz[2] = dset->daxes->yyorg;
608 break ;
609
610 default: THD_FATAL_ERROR("illegal yyorient code") ;
611 }
612
613 switch( dset->daxes->zzorient ){
614 case ORI_R2L_TYPE:
615 case ORI_L2R_TYPE:
616 tmat->mat[0][2] = dset->daxes->zzdel ; tmat->mat[0][0] = tmat->mat[0][1] = 0.0;
617 dics->xyz[0] = dset->daxes->zzorg;
618 break ;
619 case ORI_P2A_TYPE:
620 case ORI_A2P_TYPE:
621 tmat->mat[1][2] = dset->daxes->zzdel ; tmat->mat[1][0] = tmat->mat[1][1] = 0.0;
622 dics->xyz[1] = dset->daxes->zzorg;
623 break ;
624 case ORI_I2S_TYPE:
625 case ORI_S2I_TYPE:
626 tmat->mat[2][2] = dset->daxes->zzdel ; tmat->mat[2][0] = tmat->mat[2][1] = 0.0;
627 dics->xyz[2] = dset->daxes->zzorg;
628 break ;
629 default: THD_FATAL_ERROR("illegal zzorient code") ;
630 }
631
632 return ;
633 }
634
635 /* -------------------------------------------------------------------------
636 * THD_dicom_real_to_card - convert ijk_to_dicom_real coords to ijk_to_dicom
637 * version (both input and output coords are dicom)
638 *
639 * Convert input coords->xyz from ijk_to_dicom_real to ijk_to_dicom.
640 * This function essentially just modifies the coords vector.
641 * If rnd is set, round to the nearest voxel (i.e., round ijk indices).
642 *
643 *
644 * method: inverse warp dicom_real_xyz to indices
645 * possibly round
646 * warp to dicom_xyz
647 *
648 * Note: it is okay if the input or output coordinates are outside of the
649 * dataset's field of view.
650 *
651 * return 0 on success 23 Mar 2020 [rickr]
652 * -------------------------------------------------------------------------
653 */
THD_dicom_real_to_card(THD_3dim_dataset * dset,THD_fvec3 * coords,int rnd)654 int THD_dicom_real_to_card(THD_3dim_dataset *dset, THD_fvec3 * coords, int rnd)
655 {
656 THD_ivec3 iv={{0,0,0}};
657 float fi, fj, fk;
658 mat44 r2i; /* real xyz to ijk (inverse of ijk_to_dicom_real */
659
660 ENTRY("THD_dicom_real_to_card");
661
662 /* check for bad inputs */
663 if ( !dset || !dset->daxes ) {
664 WARNING_message("THD_dicom_real_to_card: null input");
665 RETURN(1);
666 }
667 if ( !ISVALID_MAT44(dset->daxes->ijk_to_dicom) ) {
668 WARNING_message("THD_dicom_real_to_card: invalid ijk_to_dicom");
669 RETURN(1);
670 }
671 if ( !ISVALID_MAT44(dset->daxes->ijk_to_dicom_real) ) {
672 WARNING_message("THD_dicom_real_to_card: invalid ijk_to_dicom_real");
673 RETURN(1);
674 }
675
676 /* get dicom_to_ijk mat */
677 r2i = nifti_mat44_inverse(dset->daxes->ijk_to_dicom_real);
678
679 /* apply r2i to input coords */
680 MAT44_VEC(r2i, coords->xyz[0], coords->xyz[1], coords->xyz[2], fi, fj, fk);
681
682 /* possibly round to nearest voxel center (okay if outside dset grid) */
683 if( rnd ) { fi = roundf(fi); fj = roundf(fj); fk = roundf(fk); }
684
685 /* apply ijk_to_dicom to index coords */
686 MAT44_VEC(dset->daxes->ijk_to_dicom, fi, fj, fk,
687 coords->xyz[0], coords->xyz[1], coords->xyz[2]);
688
689 RETURN(0);
690 }
691
692 /*---------------------------------------------------------------------
693 Return rotation and shift param. to go from i, j, k to Real
694 Dicom x,y,z
695 -----------------------------------------------------------------------*/
THD_dicom_real_xform(THD_3dim_dataset * dset,THD_dmat33 * tmat,THD_dfvec3 * dics)696 void THD_dicom_real_xform(THD_3dim_dataset *dset ,
697 THD_dmat33 *tmat, THD_dfvec3 *dics )
698 {
699 if ( !dset || !dset->daxes ||
700 !ISVALID_MAT44(dset->daxes->ijk_to_dicom_real)) {
701 THD_FATAL_ERROR("null input or no valid ijk_to_dicom_real") ;
702 }
703
704 UNLOAD_MAT44(dset->daxes->ijk_to_dicom_real,
705 tmat->mat[0][0], tmat->mat[0][1], tmat->mat[0][2], dics->xyz[0], \
706 tmat->mat[1][0], tmat->mat[1][1], tmat->mat[1][2], dics->xyz[1], \
707 tmat->mat[2][0], tmat->mat[2][1], tmat->mat[2][2], dics->xyz[2] );
708 return;
709 }
710 #define MAXNUM(a,b) ( (a) > (b) ? (a):(b))
711 #define MAX3(a,b,c) ( (MAXNUM(a,b)) > (MAXNUM(a,c)) ? (MAXNUM(a,b)):(MAXNUM(a,c)))
712 #define MINNUM(a,b) ( (a) < (b) ? (a):(b))
713 #define MIN3(a,b,c) ( (MINNUM(a,b)) < (MINNUM(a,c)) ? (MINNUM(a,b)):(MINNUM(a,c)))
714
715 /* compute angle of greatest obliquity given transformation matrix */
716
THD_compute_oblique_angle(mat44 ijk_to_dicom44,int verbose)717 float THD_compute_oblique_angle(mat44 ijk_to_dicom44, int verbose)
718 {
719 float dxtmp, dytmp, dztmp ;
720 float xmax, ymax, zmax ;
721 float fig_merit, ang_merit ;
722
723
724 dxtmp = sqrt ( ijk_to_dicom44.m[0][0] * ijk_to_dicom44.m[0][0] +
725 ijk_to_dicom44.m[1][0] * ijk_to_dicom44.m[1][0] +
726 ijk_to_dicom44.m[2][0] * ijk_to_dicom44.m[2][0] ) ;
727
728 xmax = MAX3(fabs(ijk_to_dicom44.m[0][0]),fabs(ijk_to_dicom44.m[1][0]),fabs(ijk_to_dicom44.m[2][0])) / dxtmp ;
729
730 dytmp = sqrt ( ijk_to_dicom44.m[0][1] * ijk_to_dicom44.m[0][1] +
731 ijk_to_dicom44.m[1][1] * ijk_to_dicom44.m[1][1] +
732 ijk_to_dicom44.m[2][1] * ijk_to_dicom44.m[2][1] ) ;
733
734 ymax = MAX3(fabs(ijk_to_dicom44.m[0][1]),
735 fabs(ijk_to_dicom44.m[1][1]),
736 fabs(ijk_to_dicom44.m[2][1])) / dytmp ;
737
738 dztmp = sqrt ( ijk_to_dicom44.m[0][2] * ijk_to_dicom44.m[0][2] +
739 ijk_to_dicom44.m[1][2] * ijk_to_dicom44.m[1][2] +
740 ijk_to_dicom44.m[2][2] * ijk_to_dicom44.m[2][2] ) ;
741
742 zmax = MAX3(fabs(ijk_to_dicom44.m[0][2]),
743 fabs(ijk_to_dicom44.m[1][2]),
744 fabs(ijk_to_dicom44.m[2][2])) / dztmp ;
745
746 fig_merit = MIN3(xmax,ymax,zmax) ;
747 ang_merit = acos (fig_merit) * 180.0 / 3.141592653 ;
748
749 if (fabs(ang_merit) > .01) {
750 if ( verbose ) INFO_message("%f degrees from plumb.\n",ang_merit ) ;
751 }
752 else
753 ang_merit = 0.0;
754 return(ang_merit);
755 }
756
757 /*
758 [PT: Nov 3, 2020] more funcs for dealing with obliquity matrix
759 (dset->daxes->ijk_to_dicom_real)
760
761 [PT: Nov 12, 2020] adopting the following language/terminology for
762 systematically approaching these mats, summarizing conversation with
763 RCR:
764
765 aform_real : mat44 that is the AFNI equivalent of NIFTI sform matrix
766 ('a'form for AFNI), what we typically call
767 ijk_to_dicom_real in the AFNI extension/header; likely
768 only difference is that aform_* are RAI, while sform
769 matrix is LPI (sigh). aform_real, like sform, is not
770 restricted to being orthogonal -- can represent any
771 linear affine transform.
772
773 aform_orth : mat44 that is the orthogonalized version of aform_real,
774 which should correspond to the NIFTI quaternion, since
775 the NIFTI quaternion is inherently orthogonalized when
776 created from the NIFTI sform matrix; we even use the
777 NIFTI2 mat44->quaternion->mat44 functions to calc
778 aform_orth, so it really should agree, EXCEPT that,
779 like all good things, it is RAI. While orthogonal,
780 this matrix might still contain obliquity information
781 (i.e., rotation information: it need not be diagonal).
782
783 In most cases in AFNI-generated programs, aform_orth
784 should be the same as aform_real (as of the writing of
785 this note). But we allow for reading in more general
786 matrix info.
787
788 aform_card : mat44 that is the unrotated (diagonalized) form of
789 aform_orth, what we typically call ijk_to_dicom in the
790 AFNI extension/header. This is the matrix mostly used
791 to display data in the GUI, and when we ignore
792 obliquity information in transformations, etc. Is also
793 RAI. We aim now to create this in a way that the (i*,
794 j*, k*) location where (x, y, z) = (0, 0, 0) in
795 aform_orth is preserved. That is, we preserve the
796 coordinate origin (hoping it is good!).
797 */
798
799 // get the integer form of orientation in RLPAIS ordering. So, a dset
800 // with "RAI" -> {0, 3, 4}.
THD_orient_to_int_rlpais(char ochar[4],int oint[3])801 void THD_orient_to_int_rlpais( char ochar[4], int oint[3] )
802 {
803 int i, j;
804
805 ENTRY("THD_orient_to_int_rlpais");
806
807 // translate new orient to int form
808 for( i=0 ; i<3 ; i++)
809 for( j=0 ; j<6 ; j++ ) {
810 if( strncmp(&ochar[i], &ORIENT_first[j], 1) == 0 ) {
811 oint[i] = j;
812 break;
813 }
814 if( j == 6 )
815 ERROR_message("Perm calc, couldn't find orient char match\n"
816 "\tCheck inp orient string: %s", ochar);
817 }
818
819 EXRETURN;
820 }
821
822 // get the char form of orientation from int array in RLPAIS ordering.
823 // So, a dset with {0, 3, 4} -> "RAI".
THD_int_to_orient_rlpais(int oint[3],char ochar[4])824 void THD_int_to_orient_rlpais( int oint[3], char ochar[4])
825 {
826 int i, j;
827
828 ENTRY("THD_int_to_orient_rlpais");
829
830 for( i=0 ; i<3 ; i++)
831 ochar[i] = ORIENT_first[oint[i]];
832 ochar[3] = '\0';
833
834 EXRETURN;
835 }
836
837 /* Check if orient is valid: must contain exactly one member of each
838 of the following pairs: LR, AP, IS.
839 Return 0 if valid, 1 otherwise.
840 */
is_valid_orient_char(char ochar[3])841 int is_valid_orient_char(char ochar[3])
842 {
843 int ii;
844 int oint[3] = {0, 0, 0};
845
846 ENTRY("is_valid_orient_char");
847
848 THD_orient_to_int_rlpais(ochar, oint); // get int form of orient
849 ii = is_valid_orient_int(oint); // ... and check this
850
851 return ii;
852 }
853
854 /* Check if orient is valid: must contain exactly one member of each
855 of the following pairs: 01, 23, 45.
856 Return 0 if valid, 1 otherwise.
857 */
is_valid_orient_int(int oint[3])858 int is_valid_orient_int(int oint[3])
859 {
860 int i, j;
861 int score[3] = {0, 0, 0};
862 int score_sum = 0;
863 char ochar[4];
864
865 ENTRY("is_valid_orient_int");
866
867 for( i=0 ; i<3 ; i++ )
868 score[ ORIENT_xyzint[ oint[i]]-1 ] = 1;
869
870 for( i=0 ; i<3 ; i++ )
871 score_sum+= score[i];
872
873 THD_int_to_orient_rlpais( oint, ochar );
874
875 if( score_sum == 3 ) return 1; // valid
876 else return 0; // invalid
877 }
878
879 // Calc permutation mat33 needed when changing dset orientation from
880 // orient A (e.g., "SLP") to orient B (e.g., "RAI").
THD_char_reorient_perm_mat33(char * ocharA,char * ocharB)881 mat33 THD_char_reorient_perm_mat33(char *ocharA, char *ocharB)
882 {
883 int i, j;
884 int ointA[3] = {-1, -1, -1}; // int form of orient
885 int ointB[3] = {-1, -1, -1};
886 mat33 P33;
887
888 ENTRY("THD_char_reorient_perm_mat33");
889
890 INFO_message("CHECKING %s -> %s", ocharA, ocharB);
891 if( !is_valid_orient_char(ocharA) || !is_valid_orient_char(ocharB) )
892 ERROR_exit("Invalid orientation for permuting: %s -> %s",
893 ocharA, ocharB);
894
895 // translate new orient to int form
896 THD_orient_to_int_rlpais(ocharA, ointA);
897 THD_orient_to_int_rlpais(ocharB, ointB);
898
899 P33 = THD_int_reorient_perm_mat33(ointA, ointB);
900
901 return P33;
902 }
903
904 // Calc permutation mat33 needed when changing dset orientation from
905 // orient A (e.g., "531") to orient B (e.g., "024").
THD_int_reorient_perm_mat33(int * ointA,int * ointB)906 mat33 THD_int_reorient_perm_mat33(int *ointA, int *ointB)
907 {
908 int i, j;
909 mat33 PA, PB, PAINV, PBINV; // intermediate
910 mat33 P33;
911 char ocharA[4], ocharB[4];
912
913 char ochar_rai[4] = "RAI\0";
914 int oint_rai[3] = {-1, -1, -1}; // int form of orient
915
916 ENTRY("THD_int_reorient_perm_mat33");
917
918 THD_orient_to_int_rlpais(ochar_rai, oint_rai);
919
920 LOAD_ZERO_MAT33(PA); // init mat
921 LOAD_ZERO_MAT33(PB); // init mat
922 LOAD_ZERO_MAT33(P33); // init mat
923
924 if( !is_valid_orient_int(ointA) ) {
925 THD_int_to_orient_rlpais( ointA, ocharA );
926 ERROR_exit("Dset has invalid orientation for permuting: %s", ocharA);
927 }
928
929 if( !is_valid_orient_int(ointB) ) {
930 THD_int_to_orient_rlpais( ointB, ocharB );
931 ERROR_exit("Specified orientation is invalid for permuting: %s", ocharB);
932 }
933
934 // make permutation matrix; will be applied as first arg in
935 // multiplications, such as: MAT33_MUL(P33, dset_mat33)
936 for( i=0 ; i<3 ; i++ )
937 for( j=0 ; j<3 ; j++ ){
938 if( oint_rai[i] == ointA[j] )
939 PA.m[i][j] = 1;
940 else if ( ORIENT_OPPOSITE(oint_rai[i]) == ointA[j] )
941 PA.m[i][j] = -1;
942 }
943 for( i=0 ; i<3 ; i++ )
944 for( j=0 ; j<3 ; j++ ){
945 if( oint_rai[i] == ointB[j] )
946 PB.m[i][j] = 1;
947 else if ( ORIENT_OPPOSITE(oint_rai[i]) == ointB[j] )
948 PB.m[i][j] = -1;
949 }
950
951 PAINV = MAT33_INV(PA);
952 P33 = MAT33_MUL(PB, PAINV);
953
954 return P33;
955 }
956
957
958 /* OLD, apparently broken way of thinking about this
959
960
961 // Calc permutation mat33 needed when changing dset orientation from
962 // orient A (e.g., "531") to orient B (e.g., "024").
963 mat33 THD_int_reorient_perm_mat33(int *ointA, int *ointB)
964 {
965 int i, j;
966 mat33 P33;
967 char ocharA[4], ocharB[4];
968
969 ENTRY("THD_int_reorient_perm_mat33");
970
971 LOAD_ZERO_MAT33(P33); // init mat
972
973 if( !is_valid_orient_int(ointA) ) {
974 THD_int_to_orient_rlpais( ointA, ocharA );
975 ERROR_exit("Dset has invalid orientation for permuting: %s", ocharA);
976 }
977
978 if( !is_valid_orient_int(ointB) ) {
979 THD_int_to_orient_rlpais( ointB, ocharB );
980 ERROR_exit("Specified orientation is invalid for permuting: %s", ocharB);
981 }
982
983 // make permutation matrix; will be applied as first arg in
984 // multiplications, such as: MAT33_MUL(P33, dset_mat33)
985 for( i=0 ; i<3 ; i++ )
986 for( j=0 ; j<3 ; j++ ){
987 if( ointA[i] == ointB[j] )
988 P33.m[i][j] = 1;
989 else if ( ORIENT_OPPOSITE(ointA[i]) == ointB[j] )
990 P33.m[i][j] = -1;
991 }
992
993 return P33;
994 }
995 */
996
997
998 // Calc permutation mat33 needed for ijk_to_dicom_real when changing
999 // orientation from what a current dset has to some new_ori. 'P33' is
1000 // essentially calc'ed here for input dset and new_ori.
THD_dset_reorient_perm_mat33(THD_3dim_dataset * dsetA,char * ocharB)1001 mat33 THD_dset_reorient_perm_mat33( THD_3dim_dataset *dsetA, char *ocharB)
1002 {
1003 int i, j;
1004 int ointA[3] = {-1, -1, -1}; // current dset orient (int form)
1005 int ointB[3] = {-1, -1, -1}; // new orient (int form)
1006 mat33 P33;
1007
1008 ENTRY("THD_dset_reorient_perm_mat33");
1009
1010 LOAD_ZERO_MAT33(P33); // init mat
1011
1012 if( !ISVALID_DSET(dsetA) ) return(P33);
1013
1014 THD_fill_orient_int_3_rlpais( dsetA->daxes, ointA ); // dset orient as ints
1015 THD_orient_to_int_rlpais(ocharB, ointB); // new orient as ints
1016
1017 P33 = THD_int_reorient_perm_mat33(ointA, ointB);
1018
1019 return P33;
1020 }
1021
1022 // apply permutation+new orientation to dset, calc+return new
1023 // ijk_to_dicom_real mat
THD_refit_orient_ijk_to_dicom_real(THD_3dim_dataset * dsetA,char * ocharB)1024 mat44 THD_refit_orient_ijk_to_dicom_real( THD_3dim_dataset *dsetA,
1025 char *ocharB )
1026 {
1027 int i, j;
1028 float origA[3], origB[3];
1029 mat33 P33;
1030 mat33 dsetA_mat33, dsetA_mat33_P;
1031 mat44 dsetA_mat44_P;
1032
1033 ENTRY("THD_refit_orient_ijk_to_dicom_real");
1034
1035 ZERO_MAT44(dsetA_mat44_P); // init mat
1036
1037 // checks
1038 if( !ISVALID_DSET(dsetA) || oblique_report_repeat==0 ) return(dsetA_mat44_P);
1039 THD_check_oblique_field(dsetA);
1040
1041 MAT44_TO_MAT33(dsetA->daxes->ijk_to_dicom_real, dsetA_mat33); // get current
1042
1043 P33 = THD_dset_reorient_perm_mat33( dsetA, ocharB); // perm from old to new
1044
1045 // INITIAL
1046 //INFO_message("OLD mats");
1047 //DUMP_MAT44("IJK_TO_DICOM_REAL", dsetA->daxes->ijk_to_dicom_real);
1048 //DUMP_MAT44("IJK_TO_DICOM", dsetA->daxes->ijk_to_dicom);
1049
1050 for( i=0 ; i<3 ; i++ )
1051 origA[i] = dsetA->daxes->ijk_to_dicom_real.m[i][3];
1052 MAT33_VEC(P33, origA[0], origA[1], origA[2],
1053 origB[0], origB[1], origB[2]);
1054
1055 //INFO_message("BEFORE AND AFTER");
1056 //INFO_message("%f %f %f", origA[0], origA[1], origA[2]);
1057 //INFO_message("%f %f %f", origB[0], origB[1], origB[2]);
1058
1059 // apply permutation to mat33, and then make M44
1060 dsetA_mat33_P = MAT33_MUL(P33, dsetA_mat33);
1061 MAT33_TO_MAT44(dsetA_mat33_P, dsetA_mat44_P);
1062 LOAD_MAT44_VEC(dsetA_mat44_P,
1063 origB[0], origB[1], origB[2]);
1064
1065 // FINAL
1066 //INFO_message("NEW mat");
1067 //DUMP_MAT44("", dsetA_mat44_P);
1068
1069 return dsetA_mat44_P;
1070 }
1071
1072 /*
1073 [PT: Nov 16, 2020] Use the fact that nifti_mat44_to_quatern produces
1074 a quaternion representation of the input matrix that is
1075 orthogonalized, even if the input mat isn't. Thus, when converting
1076 quatern -> mat44, the output mat should be orthogonalized.
1077
1078 Use this to calculate aform_orth (Mout) from aform_real (Min).
1079 */
nifti_orthogonalize_mat44(mat44 Min)1080 mat44 nifti_orthogonalize_mat44( mat44 Min )
1081 {
1082 mat44 Mout;
1083 float qb, qc, qd;
1084 float qx, qy, qz;
1085 float dx, dy, dz, qfac;
1086
1087 nifti_mat44_to_quatern( Min,
1088 &qb, &qc, &qd,
1089 &qx, &qy, &qz,
1090 &dx, &dy, &dz, &qfac );
1091
1092 Mout = nifti_quatern_to_mat44( qb, qc, qd,
1093 qx, qy, qz,
1094 dx, dy, dz, qfac );
1095
1096 return Mout;
1097 }
1098
1099 /*
1100 [PT: Nov 16, 2020] Test orthogonality by orthogonalizing the mat44,
1101 and then comparing sum of elementwise diffs. Could be done
1102 differently, but typically this should be fine.
1103 */
is_mat44_orthogonal(mat44 A)1104 int is_mat44_orthogonal(mat44 A)
1105 {
1106 mat44 B;
1107
1108 B = nifti_orthogonalize_mat44( A );
1109
1110 return MAT44_FLEQ(A, B);
1111 }
1112
1113
1114 /*
1115 Calc the cardinal 3x3 matrix of a dset in RAI ordering based on
1116 voxelsize
1117
1118 mat33 THD_dset_card_mat33( THD_3dim_dataset *dset )
1119 {
1120 int i,j;
1121 mat33 M;
1122 THD_dmat33 tmat;
1123 THD_dfvec3 dics;
1124
1125 THD_dicom_card_xform( dset,
1126 &tmat, &dics );
1127
1128
1129
1130 return M;
1131 }
1132
1133 */
1134
1135
1136
1137
THD_report_obliquity(THD_3dim_dataset * dset)1138 void THD_report_obliquity(THD_3dim_dataset *dset)
1139 {
1140 double angle;
1141
1142 ENTRY("THD_report_obliquity");
1143 if(AFNI_yesenv("AFNI_NO_OBLIQUE_WARNING") || !OBL_report) EXRETURN;
1144
1145 if( !ISVALID_DSET(dset) || oblique_report_repeat==0 ) EXRETURN;
1146
1147 THD_check_oblique_field(dset); /* make sure oblique field is available*/
1148 angle = THD_compute_oblique_angle(dset->daxes->ijk_to_dicom_real, 0);
1149 if(angle == 0.0) EXRETURN;
1150
1151 if(oblique_report_index<oblique_report_repeat) {
1152 if(first_oblique) {
1153 WARNING_message(
1154 " If you are performing spatial transformations on an oblique dset,\n"
1155 " such as %s,\n"
1156 " or viewing/combining it with volumes of differing obliquity,\n"
1157 " you should consider running: \n"
1158 " 3dWarp -deoblique \n"
1159 " on this and other oblique datasets in the same session.\n"
1160 " See 3dWarp -help for details.\n", DSET_BRIKNAME(dset));
1161 first_oblique = 0;
1162 }
1163
1164 INFO_message("Oblique dataset:%s is %f degrees from plumb.\n",
1165 DSET_BRIKNAME(dset), angle ) ;
1166
1167 }
1168
1169
1170 oblique_report_index++;
1171
1172 if(oblique_report_repeat2==-1) { /* report obliquity n times, stop */
1173 if(oblique_report_index>oblique_report_repeat)
1174 oblique_report_index = oblique_report_repeat;
1175 EXRETURN;
1176 }
1177
1178 /* reset counter if needed*/
1179 if(oblique_report_index>=(oblique_report_repeat+oblique_report_repeat2))
1180 oblique_report_index = 0;
1181
1182 EXRETURN;
1183 }
1184
1185 /* set the number of times to report obliquity and
1186 the number of oblique datasets to skip.
1187 If the first number is 0, don't report at all.
1188 If the second number is 0 (and the first number is not), always report.
1189 If the second number is -1, stop reporting after reporting the first n1
1190 */
THD_set_oblique_report(int n1,int n2)1191 void THD_set_oblique_report(int n1, int n2)
1192 {
1193 oblique_report_repeat = n1;
1194 oblique_report_repeat2 = n2;
1195 }
1196
1197
THD_get_oblique_report()1198 int THD_get_oblique_report()
1199 {
1200 return(oblique_report_repeat);
1201 }
1202
THD_reset_oblique_report_index()1203 void THD_reset_oblique_report_index()
1204 {
1205 oblique_report_index = 0;
1206 }
1207
1208 /* make daxes structure set to cardinal orientation versus oblique */
1209 /* usually done if not already set (in older datasets) */
THD_make_cardinal(THD_3dim_dataset * dset)1210 void THD_make_cardinal(THD_3dim_dataset *dset)
1211 {
1212 THD_dmat33 tmat ;
1213 THD_dfvec3 tvec ;
1214 mat44 Tc, Tr;
1215
1216 THD_dicom_card_xform(dset, &tmat, &tvec);
1217 LOAD_MAT44(Tc,
1218 tmat.mat[0][0], tmat.mat[0][1], tmat.mat[0][2], tvec.xyz[0],
1219 tmat.mat[1][0], tmat.mat[1][1], tmat.mat[1][2], tvec.xyz[1],
1220 tmat.mat[2][0], tmat.mat[2][1], tmat.mat[2][2], tvec.xyz[2]);
1221 dset->daxes->ijk_to_dicom_real = Tc;
1222 /* DUMP_MAT44("thd_coords: set ijk_to_dicom_real",Tc) ; */
1223 }
1224
1225 /* check the obliquity transformation field to see if it's valid */
1226 /* if not, reset it to the cardinal field */
THD_check_oblique_field(THD_3dim_dataset * dset)1227 void THD_check_oblique_field(THD_3dim_dataset *dset)
1228 {
1229 if( !ISVALID_MAT44(dset->daxes->ijk_to_dicom_real) )
1230 THD_make_cardinal(dset);
1231 }
1232
1233 /* allow for updating of obliquity - IJK_TO_DICOM_REAL attribute and
1234 structure */
THD_updating_obliquity(int update)1235 void THD_updating_obliquity(int update)
1236 {
1237 oblique_update = update;
1238 }
1239
1240 /* accessor function to get current status - can we update
1241 ijk_to_dicom_real*/
THD_update_obliquity_status()1242 int THD_update_obliquity_status()
1243 {
1244 return(oblique_update);
1245 }
1246