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 /************************************************************************
11 July 1997: moved warp manipulations routines from afni_warp.c to here.
12 *************************************************************************/
13
14 /*------------------------------------------------------------------------
15 Inputs: DICOM coordinate warp from old_dset to new_dset
16 Output: voxel index coordinates warp that can actually be used
17 --------------------------------------------------------------------------*/
18
AFNI_make_voxwarp(THD_warp * inwarp,THD_3dim_dataset * old_dset,THD_3dim_dataset * new_dset)19 THD_warp * AFNI_make_voxwarp( THD_warp * inwarp ,
20 THD_3dim_dataset * old_dset ,
21 THD_3dim_dataset * new_dset )
22 {
23 THD_warp * newwarp ;
24 THD_linear_mapping * map ;
25 THD_dataxes * new_daxes ;
26
27 newwarp = myRwcNew( THD_warp ) ;
28 newwarp->type = inwarp->type ;
29 new_daxes = CURRENT_DAXES(new_dset) ;
30
31 switch( inwarp->type ){
32
33 default:{
34 fprintf(stderr,"\a\n*** ILLEGAL warp code!!! %d\n",inwarp->type) ;
35 sleep(1) ; EXIT(1) ;
36 }
37 break ;
38
39 case WARP_AFFINE_TYPE:{
40
41 map = AFNI_make_voxmap( &(inwarp->rig_bod.warp),
42 old_dset->daxes , new_daxes ) ;
43
44 /* load the (inclusive) voxel index ranges into the affine map */
45
46 LOAD_FVEC3( map->bot, 0,0,0 ) ;
47 LOAD_FVEC3( map->top, new_daxes->nxx-1,
48 new_daxes->nyy-1, new_daxes->nzz-1 ) ;
49
50
51 newwarp->rig_bod.warp = *map ;
52
53 myRwcFree( map ) ;
54 }
55 break ;
56
57 case WARP_TALAIRACH_12_TYPE:{
58 int iw ;
59 for( iw=0 ; iw < 12 ; iw++ ){
60 map = AFNI_make_voxmap( &(inwarp->tal_12.warp[iw]) ,
61 old_dset->daxes , new_daxes ) ;
62
63 map->bot = THD_dicomm_to_3dmm(new_dset,inwarp->tal_12.warp[iw].bot);
64 map->top = THD_dicomm_to_3dmm(new_dset,inwarp->tal_12.warp[iw].top);
65
66 map->bot = THD_3dmm_to_3dfind( new_dset , map->bot ) ;
67 map->top = THD_3dmm_to_3dfind( new_dset , map->top ) ;
68
69 newwarp->tal_12.warp[iw] = *map ;
70
71 myRwcFree( map ) ;
72 }
73 }
74 break ;
75 }
76
77 return newwarp ;
78 }
79
80 /*----------------------------------------------------------------*/
81
82 #ifdef SOLARIS
do_nothing(int iii)83 void do_nothing(int iii){return ;} /* for compiler optimization error on Sun */
84 #else
85 #define do_nothing(iii)
86 #endif
87
88 /*-----------------------------------------------------------------
89 This routine takes as input a linear mapping from Dicom to Dicom
90 coordinates, and computes a linear mapping from voxel-index
91 to voxel-index coordinates as output.
92 Pointers to the old and new dataset axes structures are also input.
93
94 x3dmm_new = [new_dicomm_to_3dmm]
95 ( [dd_trans] [old_3dmm_to_dicomm] x3dmm_old - dd_base )
96
97 is the conversion between true (3dmm) coordinates. In turn, we also
98 have
99
100 x3dfind_new = [new_scale] ( x3dmm_new - new_origin )
101 x3dmm_old = [old_scale] x3dfind_old + old_origin
102
103 as the transformations between voxel-index (3dfind) coordinates
104 and true (3dmm) coordinates.
105 --------------------------------------------------------------------*/
106
AFNI_make_voxmap(THD_linear_mapping * inmap,THD_dataxes * old_daxes,THD_dataxes * new_daxes)107 THD_linear_mapping * AFNI_make_voxmap( THD_linear_mapping * inmap ,
108 THD_dataxes * old_daxes ,
109 THD_dataxes * new_daxes )
110 {
111 THD_mat33 old_scale , old_3dmm_to_dicomm ,
112 dd_trans , new_scale , new_dicomm_to_3dmm ,
113 mt ; /* temp matrix */
114
115 THD_fvec3 dd_base , new_origin , old_origin , vt0,vt1,vt2 ;
116
117 THD_linear_mapping * newmap ;
118
119 /*--- set up the elements of the transformation ---*/
120
121 dd_trans = inmap->mfor ;
122
123 LOAD_DIAG_MAT( old_scale , old_daxes->xxdel ,
124 old_daxes->yydel ,
125 old_daxes->zzdel ) ;
126
127 LOAD_DIAG_MAT( new_scale , 1.0/new_daxes->xxdel , /* notice */
128 1.0/new_daxes->yydel , /* inversion */
129 1.0/new_daxes->zzdel ) ;
130
131 old_3dmm_to_dicomm = old_daxes->to_dicomm ;
132 new_dicomm_to_3dmm = TRANSPOSE_MAT(new_daxes->to_dicomm) ; /* inversion */
133
134 /* vector elements */
135
136 dd_base = inmap->bvec ; /* in new dicomm */
137
138 LOAD_FVEC3( new_origin , new_daxes->xxorg , /* in old 3dmm */
139 new_daxes->yyorg ,
140 new_daxes->zzorg ) ;
141
142 LOAD_FVEC3( old_origin , old_daxes->xxorg , /* in new 3dmm */
143 old_daxes->yyorg ,
144 old_daxes->zzorg ) ;
145
146 /* multiply the matrices together */
147
148 mt = MAT_MUL( old_3dmm_to_dicomm , old_scale ) ; do_nothing(0) ;
149 mt = MAT_MUL( dd_trans , mt ) ; do_nothing(0) ;
150 mt = MAT_MUL( new_dicomm_to_3dmm , mt ) ; do_nothing(0) ;
151 mt = MAT_MUL( new_scale , mt ) ; do_nothing(0) ;
152
153 /* compute the new bvec */
154
155 vt0 = MATVEC( old_3dmm_to_dicomm , old_origin ) ;
156 vt0 = MATVEC( dd_trans , vt0 ) ;
157 vt0 = MATVEC( new_dicomm_to_3dmm , vt0 ) ;
158 vt0 = MATVEC( new_scale , vt0 ) ;
159
160 vt1 = MATVEC( new_dicomm_to_3dmm , dd_base ) ;
161 vt1 = MATVEC( new_scale , vt1 ) ;
162
163 vt2 = MATVEC( new_scale , new_origin ) ; /* want vt1 + vt2 - vt0 */
164
165 vt2 = ADD_FVEC3( vt1 , vt2 ) ;
166 vt2 = SUB_FVEC3( vt2 , vt0 ) ;
167
168 /* make the output map */
169
170 newmap = myRwcNew( THD_linear_mapping ) ;
171
172 newmap->type = MAPPING_LINEAR_TYPE ;
173 newmap->mfor = mt ;
174 newmap->mbac = MAT_INV(mt) ;
175 newmap->bvec = vt2 ;
176
177 newmap->svec = MATVEC(newmap->mbac,newmap->bvec) ;
178 NEGATE_FVEC3(newmap->svec) ;
179
180 return newmap ;
181 }
182
183 /*--------------------------------------------------------------------
184 this routine takes as input two warps, A and B, and produces the
185 warp A*B; B is the warp that occurs before A.
186 Restriction: A and B cannot both be Talairach_12 warps!
187 The output is placed into A (*warp_in).
188
189 Note that a NULL warp pointer for warp_prior will have the same
190 effect as an identity warp, since the routine will return immediately.
191 ----------------------------------------------------------------------*/
192
AFNI_concatenate_warp(THD_warp * warp_in,THD_warp * warp_prior)193 void AFNI_concatenate_warp( THD_warp * warp_in , THD_warp * warp_prior )
194 {
195 THD_linear_mapping * prior_map , * new_map ;
196
197 if( warp_in == NULL || warp_prior == NULL ) return ;
198
199 switch( warp_in->type + 100*warp_prior->type ){
200
201 default:
202 warp_in->type = -1 ; /* set error flag! */
203 return ;
204
205 /*-- 2 affine warps ==> a new affine warp --*/
206
207 case WARP_AFFINE_TYPE + 100*WARP_AFFINE_TYPE:{
208 prior_map = &(warp_prior->rig_bod.warp) ;
209 new_map = AFNI_concatenate_lmap(
210 &(warp_in->rig_bod.warp) , prior_map ) ;
211
212 warp_in->rig_bod.warp = *new_map ; /* write over input warp */
213 myRwcFree( new_map ) ;
214 }
215 break ;
216
217 /*--- Talairach preceeded by affine ==> new Talairach --*/
218
219 case WARP_TALAIRACH_12_TYPE + 100*WARP_AFFINE_TYPE:{
220 int iw ;
221 prior_map = &(warp_prior->rig_bod.warp) ;
222 for( iw=0 ; iw < 12 ; iw++ ){
223
224 new_map = AFNI_concatenate_lmap(
225 &(warp_in->tal_12.warp[iw]) , prior_map ) ;
226
227 warp_in->tal_12.warp[iw] = *new_map ; /* write over input warp */
228 myRwcFree( new_map ) ;
229 }
230 }
231 break ;
232
233 /*-- affine preceeded by Talairach ==> new Talairach
234 [this case is not currently used, since there are no warps
235 AFTER a Talairach warp, but it may be useful in the future]
236 -- RWCox, November 1994 A.D. --*/
237
238 case WARP_AFFINE_TYPE + 100*WARP_TALAIRACH_12_TYPE:{
239 int iw ;
240 THD_talairach_12_warp * new_warp = myRwcNew( THD_talairach_12_warp ) ;
241
242 new_warp->type = WARP_TALAIRACH_12_TYPE ;
243 for( iw=0 ; iw < 12 ; iw++ ){
244 prior_map = &(warp_prior->tal_12.warp[iw]) ;
245 new_map = AFNI_concatenate_lmap(
246 &(warp_in->rig_bod.warp) , prior_map ) ;
247 new_warp->warp[iw] = *new_map ;
248 myRwcFree( new_map ) ;
249 }
250
251 warp_in->tal_12 = *new_warp ; /* write over input warp */
252 myRwcFree( new_warp ) ;
253 }
254 break ;
255
256 } /* end of switch on warp types */
257
258 return ;
259 }
260
261 /*-----------------------------------------------------------------------
262 Concatenate 2 linear maps (including allowing for shift of origin)
263 -------------------------------------------------------------------------*/
264
AFNI_concatenate_lmap(THD_linear_mapping * map_2,THD_linear_mapping * map_1)265 THD_linear_mapping * AFNI_concatenate_lmap( THD_linear_mapping * map_2 ,
266 THD_linear_mapping * map_1 )
267 {
268 THD_linear_mapping * map_out ;
269 THD_fvec3 tvec ;
270
271 /* make a new linear mapping */
272
273 map_out = myRwcNew(THD_linear_mapping) ;
274 map_out->type = MAPPING_LINEAR_TYPE ;
275
276 /* matrix */
277
278 map_out->mfor = MAT_MUL( map_2->mfor , map_1->mfor ) ;
279 map_out->mbac = MAT_INV( map_out->mfor ) ;
280
281 /* vector */
282
283 tvec = MATVEC( map_2->mfor , map_1->bvec ) ;
284 map_out->bvec = ADD_FVEC3( tvec , map_2->bvec ) ;
285 map_out->svec = MATVEC( map_out->mbac , map_out->bvec ) ;
286 NEGATE_FVEC3(map_out->svec) ;
287
288 map_out->bot = map_2->bot ;
289 map_out->top = map_2->top ;
290
291 return map_out ;
292 }
293
294 /*--------------------------------------------------------------------------*/
295 /*! Make an affine warp from 12 input numbers:
296 - [ a11 a12 a13 ] [ s1 ]
297 - x_map = [ a21 a22 a23 ] x_in + [ s2 ]
298 - [ a31 a32 a33 ] [ s3 ]
299
300 27 Aug 2002 - RWCox.
301 ----------------------------------------------------------------------------*/
302
AFNI_make_affwarp_12(float a11,float a12,float a13,float s1,float a21,float a22,float a23,float s2,float a31,float a32,float a33,float s3)303 THD_warp * AFNI_make_affwarp_12( float a11, float a12, float a13, float s1 ,
304 float a21, float a22, float a23, float s2 ,
305 float a31, float a32, float a33, float s3 )
306 {
307 THD_warp *warp ;
308 THD_linear_mapping map ;
309 float dd , nn ;
310
311 warp = myRwcNew( THD_warp ) ;
312 warp->type = WARP_AFFINE_TYPE ;
313
314 ZZME(map) ;
315 map.type = MAPPING_LINEAR_TYPE ;
316
317 LOAD_MAT(map.mfor,a11,a12,a13,a21,a22,a23,a31,a32,a33) ;
318 dd = MAT_DET(map.mfor) ; nn = MAT_FNORM(map.mfor) ;
319 if( fabs(dd) < 1.e-5*nn*nn*nn ) return NULL ; /* bad input */
320 LOAD_FVEC3(map.bvec,-s1,-s2,-s3) ;
321 LOAD_INVERSE_LMAP(map) ;
322
323 warp->rig_bod.warp = map ;
324
325 return warp ;
326 }
327
328 /*-------------------------------------------------------------------------*/
329
AFNI_make_affwarp_mat(THD_mat33 mmm)330 THD_warp * AFNI_make_affwarp_mat( THD_mat33 mmm )
331 {
332 return AFNI_make_affwarp_12( mmm.mat[0][0], mmm.mat[0][1], mmm.mat[0][2], 0.0 ,
333 mmm.mat[1][0], mmm.mat[1][1], mmm.mat[1][2], 0.0 ,
334 mmm.mat[2][0], mmm.mat[2][1], mmm.mat[2][2], 0.0 ) ;
335 }
336
337 /*-------------------------------------------------------------------------*/
338
AFNI_make_affwarp_matvec(THD_mat33 mmm,THD_fvec3 vvv)339 THD_warp * AFNI_make_affwarp_matvec( THD_mat33 mmm , THD_fvec3 vvv )
340 {
341 return AFNI_make_affwarp_12( mmm.mat[0][0], mmm.mat[0][1], mmm.mat[0][2], vvv.xyz[0] ,
342 mmm.mat[1][0], mmm.mat[1][1], mmm.mat[1][2], vvv.xyz[1] ,
343 mmm.mat[2][0], mmm.mat[2][1], mmm.mat[2][2], vvv.xyz[2] ) ;
344 }
345