1 #include "mrilib.h"
2
3 /*-------------------------------------------------------------------
4 08 Feb 2001: read a matrix+vector from a file - RWCox
5 14 Feb 2001: modified to add "-rotate ..." mode
6 10 Aug 2005: add MATRIX(...)
7 ---------------------------------------------------------------------*/
8
THD_read_dvecmat(char * fname,int invert)9 THD_dvecmat THD_read_dvecmat( char *fname , int invert )
10 {
11 THD_dvecmat dvm ;
12 THD_dmat33 tmat ;
13 THD_dfvec3 tvec ;
14 int nn , loaded=0 ;
15
16 ENTRY("THD_read_dvecmat") ;
17
18 LOAD_DIAG_DMAT(tmat,0.0,0.0,0.0) ; /* initialize to all 0s */
19 LOAD_DFVEC3(tvec,0.0,0.0,0.0) ;
20
21 if( fname == NULL || fname[0] == '\0' ){
22 dvm.vv = tvec ; dvm.mm = tmat ; RETURN(dvm) ;
23 }
24
25 /*-- filename is "dataset::attribute" --*/
26
27 if( strstr(fname,"::") != NULL ){ /*== read from dataset header ==*/
28 char *dname = strdup(fname) ;
29 char *cc = strstr(dname,"::") ;
30 char *ss ; int iss=0 ;
31 THD_3dim_dataset *dset ;
32 ATR_float *atr, *atrc ;
33 int incl;
34 THD_dfvec3 tvec_co, tvec_cb;
35
36 *cc = '\0' ; cc += 2 ; /* dname = dataset name ; cc = attribute name */
37
38 dset = THD_open_one_dataset( dname ) ;
39 if( !ISVALID_DSET(dset) ){
40 ERROR_message("THD_read_dvecmat: can't open dataset %s\n",dname) ;
41 dvm.vv = tvec ; dvm.mm = tmat ; free(dname) ; RETURN(dvm) ;
42 }
43
44 if (strstr(cc,"IJK_TO_CARD_DICOM")) { /* return the matrix that turns i,j,k to mm RAI ZSS May 07*/
45 THD_dicom_card_xform(dset, &tmat, &tvec);
46 } if (strstr(cc,"IJK_TO_DICOM_REAL")) {
47 THD_dicom_real_xform(dset, &tmat, &tvec);
48 } else {
49 ss = strstr(cc,"[") ;
50 if( ss != NULL ){
51 *ss = '\0'; ss++; iss = (int)strtol(ss,NULL,10); if( iss < 0 ) iss = 0;
52 }
53 atr = THD_find_float_atr( dset->dblk , cc ) ;
54 if( atr == NULL ){
55 ERROR_message("THD_read_dvecmat: can't find attribute %s in dataset %s\n",
56 cc,dname) ;
57 dvm.vv = tvec ; dvm.mm = tmat ; free(dname) ; RETURN(dvm) ;
58 }
59
60 if( strcmp(cc,"WARP_DATA") == 0 ){ /* 10 Aug 2005 */
61
62 LOAD_DMAT(tmat,atr->fl[0],atr->fl[1],atr->fl[2],
63 atr->fl[3],atr->fl[4],atr->fl[5],
64 atr->fl[6],atr->fl[7],atr->fl[8] ) ;
65 LOAD_DFVEC3(tvec,-atr->fl[18],-atr->fl[19],-atr->fl[20]) ;
66 loaded = 1 ;
67 }
68
69 if( !loaded ){
70 switch( atr->nfl ){
71
72 default:
73 ERROR_message(
74 "THD_read_dvecmat: can't read matrix+vector from %s::%s\n",
75 dname,cc) ;
76 dvm.vv = tvec ; dvm.mm = tmat ; free(dname) ; RETURN(dvm) ;
77 break ; /* unreachable */
78
79 case 12: LOAD_DMAT(tmat,atr->fl[0],atr->fl[1],atr->fl[2],
80 atr->fl[4],atr->fl[5],atr->fl[6],
81 atr->fl[8],atr->fl[9],atr->fl[10] ) ;
82 LOAD_DFVEC3(tvec,atr->fl[3],atr->fl[7],atr->fl[11]) ;
83 break ;
84
85 case 9: LOAD_DMAT(tmat,atr->fl[0],atr->fl[1],atr->fl[2],
86 atr->fl[3],atr->fl[4],atr->fl[5],
87 atr->fl[6],atr->fl[7],atr->fl[8] ) ;
88 LOAD_DFVEC3(tvec,0,0,0) ;
89 break ;
90 }
91 }
92
93 {
94 /* ZSS: Dec 27, the year of our Lord when the War on Christmas was raging */
95 /* Need to include VOLREG_CENTER_OLD, VOLREG_CENTER_BASE,
96 ROTATE_CENTER_OLD, ROTATE_CENTER_BASE,
97 if applicable. */
98 /* Do we have a ROTATE or VOLREG attribute ?*/
99 incl = 0;
100 if (strstr(cc,"VOLREG_MATVEC")) {
101 atrc = THD_find_float_atr( dset->dblk , "VOLREG_CENTER_OLD");
102 if ( atrc ) {
103 LOAD_DFVEC3(tvec_co,atrc->fl[0],atrc->fl[1],atrc->fl[2]) ;
104 incl = 1;
105 } else {
106 LOAD_DFVEC3(tvec_co,0,0,0) ;
107 }
108
109 atrc = THD_find_float_atr( dset->dblk , "VOLREG_CENTER_BASE");
110 if ( atrc ) {
111 LOAD_DFVEC3(tvec_cb,atrc->fl[0],atrc->fl[1],atrc->fl[2]) ;
112 incl = 1;
113 } else {
114 LOAD_DFVEC3(tvec_cb,0,0,0) ;
115 }
116
117 if (incl == 1)
118 INFO_message(
119 "THD_read_dvecmat:\n"
120 " Including VOLREG_CENTER_BASE and VOLREG_CENTER_OLD\n"
121 " attributes in final transform\n");
122 else INFO_message(
123 "THD_read_dvecmat:\n"
124 " No VOLREG_CENTER_BASE or VOLREG_CENTER_OLD\n"
125 " attributes found with VOLREG_MATVEC\n");
126 } else if (strstr(cc,"ROTATE_MATVEC")) {
127 atrc = THD_find_float_atr( dset->dblk , "ROTATE_CENTER_OLD");
128 if ( atrc ) {
129 LOAD_DFVEC3(tvec_co,atrc->fl[0],atrc->fl[1],atrc->fl[2]) ;
130 incl = 1;
131 } else {
132 LOAD_DFVEC3(tvec_co,0,0,0) ;
133 }
134
135 atrc = THD_find_float_atr( dset->dblk , "ROTATE_CENTER_BASE");
136 if ( atrc ) {
137 LOAD_DFVEC3(tvec_cb,atrc->fl[0],atrc->fl[1],atrc->fl[2]) ;
138 incl = 1;
139 } else {
140 LOAD_DFVEC3(tvec_cb,0,0,0) ;
141 }
142
143 if (incl == 1) INFO_message("THD_read_dvecmat:\n"
144 " Including ROTATE_CENTER_BASE and ROTATE_CENTER_OLD\n"
145 " attributes in final transform\n");
146 else INFO_message("THD_read_dvecmat:\n"
147 " No ROTATE_CENTER_BASE or ROTATE_CENTER_OLD\n"
148 " attributes found with ROTATE_MATVEC\n");
149 }
150 if (incl == 1) {
151 tvec_co = DMATVEC(tmat, tvec_co);
152 NEGATE_DFVEC3(tvec_co);
153 tvec.xyz[0] += tvec_cb.xyz[0] + tvec_co.xyz[0];
154 tvec.xyz[1] += tvec_cb.xyz[1] + tvec_co.xyz[1];
155 tvec.xyz[2] += tvec_cb.xyz[2] + tvec_co.xyz[2];
156 }
157 }
158 }
159 free(dname) ; DSET_delete(dset) ;
160
161 /*-- 14 Feb 2001: filename is "-rotate a b c -[ab]shift x y z" string --*/
162
163 } else if( strstr(fname,"-rotate") != NULL ){ /*== compute directly ==*/
164
165 dvm = THD_rotcom_to_matvec( NULL , fname ) ; /* thd_rotangles.c */
166 tvec = dvm.vv ; tmat = dvm.mm ;
167
168 /*-- MATRIX(...) --*/
169
170 } else if( strncmp(fname,"MATRIX(",7) == 0 ){
171 double matar[12] ; int nn ;
172 nn = sscanf(fname,"MATRIX(%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf",
173 matar+0 , matar+1 , matar+2 , matar+3 ,
174 matar+4 , matar+5 , matar+6 , matar+7 ,
175 matar+8 , matar+9 , matar+10, matar+11 ) ;
176 if( nn < 12 ){
177 ERROR_message("badly formatted: %s",fname) ;
178 } else {
179 LOAD_DMAT(tmat,matar[0],matar[1],matar[2],
180 matar[4],matar[5],matar[6],
181 matar[8],matar[9],matar[10] ) ;
182 LOAD_DFVEC3(tvec,matar[3],matar[7],matar[11]) ;
183 }
184
185 /*-- just a normal filename --*/
186
187 } else { /*== read numbers from file ==*/
188
189 MRI_IMAGE *dim; double *dar;
190 dim = mri_read_double_ascii( fname ) ;
191 if( dim == NULL || dim->nvox < 9 ){
192 ERROR_message("THD_read_dvecmat: can't read matrix+vector from '%s'\n",
193 fname) ;
194 dvm.vv = tvec ; dvm.mm = tmat ; RETURN(dvm) ;
195 }
196 dar = MRI_DOUBLE_PTR(dim) ;
197 if( dim->nvox < 12 ){
198 LOAD_DMAT(tmat,dar[0],dar[1],dar[2], /* 9 == matrix only */
199 dar[3],dar[4],dar[5],
200 dar[6],dar[7],dar[8] ) ;
201 } else {
202 LOAD_DMAT(tmat,dar[0],dar[1],dar[2], /* 12 ==> matrix+vector */
203 dar[4],dar[5],dar[6],
204 dar[8],dar[9],dar[10] ) ;
205 LOAD_DFVEC3(tvec,dar[3],dar[7],dar[11]) ;
206 }
207 }
208
209 /*-- invert the transformation we just read? --*/
210 /*-- [y] = [R][x] + [v] is transformation, so --*/
211 /*-- [x] = inv[R] - inv[R][v] is inverse transformation --*/
212
213 if( invert ){
214 THD_dmat33 imat ; THD_dfvec3 ivec ;
215 imat = DMAT_INV(tmat) ; /* matrix inverse */
216 ivec = DMATVEC(imat,tvec) ; /* multiply inverse into vector */
217 tmat = imat ;
218 tvec = ivec ; NEGATE_DFVEC3(tvec) ; /* negate vector */
219 }
220
221 /*-- store results and get outta here, dude! --*/
222
223 dvm.vv = tvec ; dvm.mm = tmat ; RETURN(dvm) ;
224 }
225
226 /*-------------------------------------------------------------------------------*/
227
invert_dvecmat(THD_dvecmat avm)228 THD_dvecmat invert_dvecmat( THD_dvecmat avm ) /* 24 Jul 2007 */
229 {
230 THD_dvecmat bvm ;
231 THD_dmat33 imat ; THD_dfvec3 ivec ;
232 imat = DMAT_INV(avm.mm) ; /* matrix inverse */
233 ivec = DMATVEC(imat,avm.vv) ; /* multiply inverse into vector */
234 NEGATE_DFVEC3(ivec) ;
235 bvm.mm = imat ; bvm.vv = ivec ; return bvm ;
236 }
237
238 /*-------------------------------------------------------------------------------*/
239
sqrt_dvecmat(THD_dvecmat avm)240 THD_dvecmat sqrt_dvecmat( THD_dvecmat avm ) /* 30 Jul 2007 */
241 {
242 mat44 A , B ; THD_dvecmat bvm ;
243 VECMAT_TO_MAT44( avm , A ) ;
244 B = MAT44_SQRT( A ) ;
245 MAT44_TO_VECMAT( B , bvm ) ;
246 return bvm ;
247 }
248