1 #include "mrilib.h"
2
3 /*-------------------------------------------------------------------*/
4
5 #undef OV
6 #undef EP
7 #undef AN
8
9 #define OV(i,j,k) ov[(i)+(j)*ovx+(k)*ovxy]
10 #define EP(i,j,k) epiar[(i)+(j)*nxepi+(k)*nxyepi]
11 #define AN(i,j,k) antar[(i)+(j)*nxant+(k)*nxyant]
12
THD_autonudge(THD_3dim_dataset * dsepi,int ivepi,THD_3dim_dataset * dsant,int ivant,float step,int xstep,int ystep,int zstep,int code)13 THD_fvec3 THD_autonudge( THD_3dim_dataset *dsepi, int ivepi,
14 THD_3dim_dataset *dsant, int ivant,
15 float step, int xstep, int ystep, int zstep, int code )
16 {
17 THD_fvec3 fv1,fv2 , fvorg_old,fvorg_new , dxorg ;
18 THD_ivec3 iv1,iv2 ;
19 float *tar ;
20 byte *epiar , *antar ;
21 int ii,jj,kk , nxepi,nyepi,nzepi , nxyepi,nxyzepi ;
22 int nxant,nyant,nzant , nxyant,nxyzant ;
23 MRI_IMAGE *tim ;
24 float xorgepi , yorgepi , zorgepi , xx1,xx2,yy1,yy2,zz1,zz2 ;
25 float epiclip , xorgant,yorgant,zorgant , f1,f2,g1,g2,h1,h2 , f,g,h ;
26 float tx,ty,tz , dxepi,dyepi,dzepi , dxant,dyant,dzant , z1,z2,y1,y2,x1,x2 ;
27 int *ov,*ovp , ovx=2*xstep+1 , ovy=2*ystep+1 , ovz=2*zstep+1 , ovxy=ovx*ovy ;
28 int xant,yant,zant , pp,qq,rr , i,j,k , ip,jp,kp , ovtop , kstep ;
29 float dxyz_ratio , vsum_thresh , vsum , sx,sy,sz ;
30 int verb = ((code & 1) != 0) ;
31
32 /*-- start the action! --*/
33
34 ENTRY("THD_autonudge") ;
35
36 /*-- sanity checks --*/
37
38 if( !ISVALID_DSET(dsepi) ||
39 !ISVALID_DSET(dsant) ||
40 ivepi < 0 || ivepi >= DSET_NVALS(dsepi) ||
41 ivant < 0 || ivant >= DSET_NVALS(dsant) ||
42 step <= 0.0 || ovx < 1 || ovy < 1 || ovz < 1 || ovx*ovy*ovz < 3 ){
43
44 fprintf(stderr,"THD_autonudge: bad inputs!\n") ; EXIT(1) ;
45 }
46
47 /*-- load chosen sub-brick of epi into local float array --*/
48
49 if( DSET_ARRAY(dsepi,ivepi) == NULL ){
50 DSET_load(dsepi) ;
51 if( !DSET_LOADED(dsepi) ){
52 fprintf(stderr,"THD_autonudge: can't load %s\n",DSET_HEADNAME(dsepi));
53 EXIT(1) ;
54 }
55 }
56
57 nxepi = DSET_NX(dsepi) ;
58 nyepi = DSET_NY(dsepi) ; nxyepi = nxepi * nyepi ;
59 nzepi = DSET_NZ(dsepi) ; nxyzepi = nxyepi * nzepi ;
60
61 tar = (float *) malloc( sizeof(float) * nxyzepi ) ;
62 if( tar == NULL ){
63 fprintf(stderr,"THD_autonudge: malloc failure for epiar\n"); EXIT(1);
64 }
65
66 EDIT_coerce_scale_type( nxyzepi ,
67 DSET_BRICK_FACTOR(dsepi,ivepi) ,
68 DSET_BRICK_TYPE(dsepi,ivepi) ,
69 DSET_ARRAY(dsepi,ivepi) , MRI_float , tar ) ;
70 DSET_unload(dsepi) ;
71
72 /*-- clip epi array values --*/
73
74 tim = mri_new_vol_empty( nxepi , nyepi , nzepi , MRI_float ) ;
75 mri_fix_data_pointer( tar , tim ) ;
76 epiclip = THD_cliplevel( tim , 0.5 ) ; /* get clip value */
77 mri_clear_data_pointer(tim) ; mri_free(tim) ;
78
79 if( epiclip <= 0.0 ){
80 fprintf(stderr,"THD_autonudge: can't compute epiclip\n"); EXIT(1);
81 }
82
83 if( verb )
84 fprintf(stderr,"THD_autonudge: epi clip level=%g\n",epiclip) ;
85
86 epiar = (byte *) malloc(sizeof(byte)*nxyzepi) ;
87 if( epiar == 0 ){
88 fprintf(stderr,"THD_autonudge: malloc failed for epiar\n"); EXIT(1);
89 }
90
91 for( ii=0 ; ii < nxyzepi ; ii++ ) /* mask of supra-clip voxels */
92 epiar[ii] = (tar[ii] > epiclip) ;
93
94 free(tar) ;
95
96 /*-- load chosen sub-brick of ant into local float array --*/
97
98 if( DSET_ARRAY(dsant,ivant) == NULL ){
99 DSET_load(dsant) ;
100 if( !DSET_LOADED(dsant) ){
101 fprintf(stderr,"THD_autonudge: can't load %s\n",DSET_HEADNAME(dsant));
102 EXIT(1) ;
103 }
104 }
105
106 nxant = DSET_NX(dsant) ;
107 nyant = DSET_NY(dsant) ; nxyant = nxant * nyant ;
108 nzant = DSET_NZ(dsant) ; nxyzant = nxyant * nzant ;
109
110 tar = (float *) malloc( sizeof(float) * nxyzant ) ;
111 if( tar == NULL ){
112 fprintf(stderr,"THD_autonudge: malloc failure for antar\n"); EXIT(1);
113 }
114
115 EDIT_coerce_scale_type( nxyzant ,
116 DSET_BRICK_FACTOR(dsant,ivant) ,
117 DSET_BRICK_TYPE(dsant,ivant) ,
118 DSET_ARRAY(dsant,ivant) , MRI_float , tar ) ;
119 DSET_unload(dsant) ;
120
121 antar = (byte *) malloc(sizeof(byte)*nxyzant) ;
122 if( antar == NULL ){
123 fprintf(stderr,"THD_autonudge: malloc failure for antar\n"); EXIT(1);
124 }
125
126 for( ii=0 ; ii < nxyzant ; ii++ ) /* make mask */
127 antar[ii] = (tar[ii] > 0.0) ;
128
129 free(tar) ;
130
131 /*-- find axis in ant that corresponds to x-axis in epi --*/
132
133 LOAD_FVEC3(fv1,0,0,0) ;
134 fv1 = THD_3dfind_to_3dmm( dsepi, fv1 ) ; /* coords in dsepi */
135 fv1 = THD_3dmm_to_dicomm( dsepi, fv1 ) ; /* DICOM in dsepi */
136 fv1 = THD_dicomm_to_3dmm( dsant, fv1 ) ; /* coords in dsant */
137 iv1 = THD_3dmm_to_3dind ( dsant, fv1 ) ; /* index in dsant */
138
139 LOAD_FVEC3(fv2,nxepi-1,0,0) ;
140 fv2 = THD_3dfind_to_3dmm( dsepi , fv2 ) ; /* coords in dsepi */
141 fv2 = THD_3dmm_to_dicomm( dsepi , fv2 ) ; /* DICOM in dsepi */
142 fv2 = THD_dicomm_to_3dmm( dsant, fv2 ) ; /* coords in dsant */
143 iv2 = THD_3dmm_to_3dind ( dsant, fv2 ) ; /* index in dsant */
144
145 if( iv1.ijk[0] != iv2.ijk[0] ) xant = 0 ; /* epi x-axis */
146 else if( iv1.ijk[1] != iv2.ijk[1] ) xant = 1 ; /* in ant */
147 else if( iv1.ijk[2] != iv2.ijk[2] ) xant = 2 ;
148 else {
149 fprintf(stderr,"THD_autonudge: incoherent x slicing?!\n");
150 DUMP_IVEC3("iv1",iv1); DUMP_IVEC3("iv2",iv2); EXIT(1);
151 }
152
153 /*-- find axis in ant that corresponds to y-axis in epi --*/
154
155 LOAD_FVEC3(fv2,0,nyepi-1,0) ;
156 fv2 = THD_3dfind_to_3dmm( dsepi, fv2 ) ; /* coords in dsepi */
157 fv2 = THD_3dmm_to_dicomm( dsepi, fv2 ) ; /* DICOM in dsepi */
158 fv2 = THD_dicomm_to_3dmm( dsant, fv2 ) ; /* coords in dsant */
159 iv2 = THD_3dmm_to_3dind ( dsant, fv2 ) ; /* index in dsant */
160
161 if( iv1.ijk[0] != iv2.ijk[0] ) yant = 0 ; /* epi y-axis */
162 else if( iv1.ijk[1] != iv2.ijk[1] ) yant = 1 ; /* in ant */
163 else if( iv1.ijk[2] != iv2.ijk[2] ) yant = 2 ;
164 else {
165 fprintf(stderr,"THD_autonudge: incoherent y slicing?!\n");
166 DUMP_IVEC3("iv1",iv1); DUMP_IVEC3("iv2",iv2); EXIT(1);
167 }
168
169 /*-- find axis in ant that corresponds to z-axis in epi --*/
170
171 LOAD_FVEC3(fv2,0,0,nzepi-1) ;
172 fv2 = THD_3dfind_to_3dmm( dsepi, fv2 ) ; /* coords in dsepi */
173 fv2 = THD_3dmm_to_dicomm( dsepi, fv2 ) ; /* DICOM in dsepi */
174 fv2 = THD_dicomm_to_3dmm( dsant, fv2 ) ; /* coords in dsant */
175 iv2 = THD_3dmm_to_3dind ( dsant, fv2 ) ; /* index in dsant */
176
177 if( iv1.ijk[0] != iv2.ijk[0] ) zant = 0 ; /* epi z-axis */
178 else if( iv1.ijk[1] != iv2.ijk[1] ) zant = 1 ; /* in ant */
179 else if( iv1.ijk[2] != iv2.ijk[2] ) zant = 2 ;
180 else {
181 fprintf(stderr,"THD_autonudge: incoherent z slicing?!\n");
182 DUMP_IVEC3("iv1",iv1); DUMP_IVEC3("iv2",iv2); EXIT(1);
183 }
184
185 #if 0
186 if( verb )
187 fprintf(stderr," xant=%d yant=%d zant=%d\n",
188 xant,yant,zant) ;
189 #endif
190
191 if( ((1<<xant) | (1<<yant) | (1<<zant)) != 7 ){
192 fprintf(stderr,"THD_autonudge: incoherent xyz slicing!\n"
193 " xant=%d yant=%d zant=%d\n",
194 xant,yant,zant) ;
195 EXIT(1) ;
196 }
197
198 /*-- allocate space for array of overlap counts --*/
199
200 ov = (int *) calloc( sizeof(int) , ovx*ovy*ovz ) ;
201 if( ov == NULL ){
202 fprintf(stderr,"THD_autonudge: can't malloc space for overlap counts!\n");
203 EXIT(1) ;
204 }
205
206 /*-- for each origin shift in the xant,yant directions,
207 shift origin of ant dataset, then compare to epi dataset --*/
208
209 xorgepi = DSET_XORG(dsepi) ; dxepi = DSET_DX(dsepi) ;
210 yorgepi = DSET_YORG(dsepi) ; dyepi = DSET_DY(dsepi) ;
211 zorgepi = DSET_ZORG(dsepi) ; dzepi = DSET_DZ(dsepi) ;
212
213 xorgant = DSET_XORG(dsant) ; dxant = DSET_DX(dsant) ;
214 yorgant = DSET_YORG(dsant) ; dyant = DSET_DY(dsant) ;
215 zorgant = DSET_ZORG(dsant) ; dzant = DSET_DZ(dsant) ;
216
217 dxyz_ratio = fabs( (dxepi*dyepi*dzepi)/(dxant*dyant*dzant) ) ;
218 vsum_thresh = 0.5*dxyz_ratio ;
219
220 LOAD_FVEC3( fvorg_old,xorgant,yorgant,zorgant) ;
221 LOAD_FVEC3( dxorg ,dxant ,dyant ,dzant ) ;
222
223 ovtop = 0 ;
224 if( nzepi < 11 ){
225 kstep = 2 ;
226 } else {
227 kstep = (int)(0.2*nzepi+0.5) ;
228 }
229
230 for( pp=0 ; pp < ovx ; pp++ ){ /* loop over shifts in 3 directions */
231 for( qq=0 ; qq < ovy ; qq++ ){
232 for( rr=0 ; rr < ovz ; rr++ ){
233
234 if( verb ) fprintf(stderr," starting shift %2d %2d %2d ",
235 pp-xstep,qq-ystep,rr-zstep) ;
236
237 fvorg_new = fvorg_old ; /* old ant origin */
238
239 fvorg_new.xyz[xant] += (pp-xstep)*step*dxorg.xyz[xant]; /* shift epi-x */
240 fvorg_new.xyz[yant] += (qq-ystep)*step*dxorg.xyz[yant]; /* shift epi-y */
241 fvorg_new.xyz[zant] += (rr-zstep)*step*dxorg.xyz[zant]; /* shift epi-z */
242
243 xorgant = fvorg_new.xyz[0] ; /* load new ant origin */
244 yorgant = fvorg_new.xyz[1] ;
245 zorgant = fvorg_new.xyz[2] ;
246
247 ovp = ov + (pp + qq*ovx + rr*ovxy) ; /* place to store result */
248
249 /*-- foreach voxel in epi dataset,
250 find how much of it is filled by nonzero ant voxels --*/
251
252 for( kk=0 ; kk < nzepi ; kk++ ){
253 z1 = zorgepi + dzepi*(kk-0.5) ; z2 = zorgepi + dzepi*(kk+0.49999) ;
254
255 if( verb && kk%kstep == 0 ) fprintf(stderr,".") ;
256
257 for( jj=0 ; jj < nyepi ; jj++ ){
258 y1 = yorgepi + dyepi*(jj-0.5) ; y2 = yorgepi + dyepi*(jj+0.49999) ;
259
260 for( ii=0 ; ii < nxepi ; ii++ ){
261 if( EP(ii,jj,kk) == 0 ) continue ; /* skip voxel */
262
263 x1 = xorgepi + dxepi*(ii-0.5) ; x2 = xorgepi + dxepi*(ii+0.49999) ;
264
265 /* epi voxel covers coords [x1..x2] X [y1..y2] X [z1..z2] */
266
267 /* transform these to ant dataset coords */
268
269 LOAD_FVEC3(fv1,x1,y1,z1) ; /* coords in epi */
270 fv1 = THD_3dmm_to_dicomm( dsepi, fv1 ) ; /* DICOM coords */
271 fv1 = THD_dicomm_to_3dmm( dsant, fv1 ) ; /* coords in ant */
272 UNLOAD_FVEC3(fv1,xx1,yy1,zz1) ;
273
274 LOAD_FVEC3(fv2,x2,y2,z2) ; /* coords in epi */
275 fv2 = THD_3dmm_to_dicomm( dsepi, fv2 ) ; /* DICOM coords */
276 fv2 = THD_dicomm_to_3dmm( dsant, fv2 ) ; /* coords in ant */
277 UNLOAD_FVEC3(fv2,xx2,yy2,zz2) ;
278
279 /* epi voxel spans ant coords [xx1..xx2] X [yy1..yy2] X [zz1..zz2] */
280
281 /* compute indices into ant dataset voxels */
282
283 f1 = (xx1-xorgant)/dxant + 0.49999 ; f2 = (xx2-xorgant)/dxant + 0.49999 ;
284 if( f1 > f2 ){ tx = f1 ; f1 = f2 ; f2 = tx ; }
285 if( f1 >= nxant || f2 <= 0.0 ) continue ;
286 if( f1 < 0.0 ) f1 = 0.0 ; if( f2 >= nxant ) f2 = nxant - 0.001 ;
287
288 g1 = (yy1-yorgant)/dyant + 0.49999 ; g2 = (yy2-yorgant)/dyant + 0.49999 ;
289 if( g1 > g2 ){ ty = g1 ; g1 = g2 ; g2 = ty ; }
290 if( g1 >= nyant || g2 <= 0.0 ) continue ;
291 if( g1 < 0.0 ) g1 = 0.0 ; if( g2 >= nyant ) g2 = nyant - 0.001 ;
292
293 h1 = (zz1-zorgant)/dzant + 0.49999 ; h2 = (zz2-zorgant)/dzant + 0.49999 ;
294 if( h1 > h2 ){ tz = h1 ; h1 = h2 ; h2 = tz ; }
295 if( h1 >= nzant || h2 <= 0.0 ) continue ;
296 if( h1 < 0.0 ) h1 = 0.0 ; if( h2 >= nzant ) h2 = nzant - 0.001 ;
297
298 /* epi voxel covers voxels [f1..f2] X [g1..g2] X [h1..h2] in ant */
299
300 /* loop over these, and count how much is nonzero */
301
302 vsum = 0.0 ;
303 for( f=f1 ; f < f2 ; f = ip ){
304 i = (int) f ; ip = i+1 ; tx = MIN(ip,f2) ; sx = tx - f ;
305 for( g=g1 ; g < g2 ; g = jp ){
306 j = (int) g ; jp = j+1 ; ty = MIN(jp,g2) ; sy = ty - g ;
307 for( h=h1 ; h < h2 ; h = kp ){
308 k = (int) h ; kp = k+1 ; tz = MIN(kp,h2) ; sz = tz - h ;
309 if( AN(i,j,k) ) vsum += sx * sy * sz ;
310 }}} /* end of loop over ant voxels */
311
312 #if 0
313 fprintf(stderr," epi=%d %d %d ant=%6.2f..%6.2f %6.2f..%6.2f %6.2f..%6.2f vsum=%6.2f\n",
314 ii,jj,kk , f1,f2 , g1,g2 , h1,h2 , vsum) ;
315 #endif
316
317 /* add to results for this shift */
318
319 if( vsum > vsum_thresh ) (*ovp)++ ;
320
321 }}} /* end of loop over epi voxels */
322
323 if( verb ) fprintf(stderr," overlap=%d",*ovp) ;
324 if( *ovp > ovtop ){
325 ovtop = *ovp ; if( verb ) fprintf(stderr," *") ;
326 }
327 if( verb ) fprintf(stderr,"\n") ;
328
329 }}} /* end of loop over shifts */
330
331 /*-- find best shift in list --*/
332
333 ii = jj = kk = ip = 0 ;
334
335 for( pp=0 ; pp < ovx ; pp++ ){ /* loop over shifts in 3 directions */
336 for( qq=0 ; qq < ovy ; qq++ ){
337 for( rr=0 ; rr < ovz ; rr++ ){
338 if( OV(pp,qq,rr) > ip ){
339 ii = pp ; jj = qq ; kk = rr ; ip = OV(pp,qq,rr) ;
340 }
341 }}}
342
343 fvorg_new.xyz[xant] = (ii-xstep)*step*dxorg.xyz[xant]; /* shift epi-x */
344 fvorg_new.xyz[yant] = (jj-ystep)*step*dxorg.xyz[yant]; /* shift epi-y */
345 fvorg_new.xyz[zant] = (kk-zstep)*step*dxorg.xyz[zant]; /* shift epi-z */
346
347 if( verb ){
348 fprintf(stderr," best shift: %d %d %d overlap=%d\n",
349 ii-xstep,jj-ystep,kk-zstep,ip) ;
350 DUMP_FVEC3(" best shift",fvorg_new) ;
351 }
352
353 free(ov) ; free(antar) ; free(epiar) ;
354
355 RETURN( fvorg_new ) ;
356 }
357