1 #include "mrilib.h"
2
3 /*----------------------------------------------------------------------------*/
4 #undef TRIPROD
5 #define TRIPROD(ax,ay,az,bx,by,bz,cx,cy,cz) ( (ax)*((by)*(cz)-(bz)*(cy)) \
6 +(bx)*((cy)*(az)-(cz)*(ay)) \
7 +(cx)*((ay)*(bz)-(az)*(by)) )
8
9 #undef DA
10 #define DA(p,q) (p.a-q.a)
11 #undef DB
12 #define DB(p,q) (p.b-q.b)
13 #undef DC
14 #define DC(p,q) (p.c-q.c)
15
hexahedron_volume(float_triple x0,float_triple x1,float_triple x2,float_triple x3,float_triple x4,float_triple x5,float_triple x6,float_triple x7)16 static float hexahedron_volume( float_triple x0 , float_triple x1 , float_triple x2 ,
17 float_triple x3 , float_triple x4 , float_triple x5 ,
18 float_triple x6 , float_triple x7 )
19 {
20 float xa,ya,za , xb,yb,zb , xc,yc,zc , vol ;
21
22 xa = DA(x7,x1)+DA(x6,x0); ya = DB(x7,x1)+DB(x6,x0); za = DC(x7,x1)+DC(x6,x0);
23 xb = DA(x7,x2) ; yb = DB(x7,x2) ; zb = DC(x7,x2) ;
24 xc = DA(x3,x0) ; yc = DB(x3,x0) ; zc = DC(x3,x0) ;
25 vol = TRIPROD(xa,ya,za,xb,yb,zb,xc,yc,zc) ;
26 xa = DA(x6,x0) ; ya = DB(x6,x0) ; za = DC(x6,x0) ;
27 xb = DA(x7,x2)+DA(x5,x0); yb = DB(x7,x2)+DB(x5,x0); zb = DC(x7,x2)+DC(x5,x0);
28 xc = DA(x7,x4) ; yc = DB(x7,x4) ; zc = DC(x7,x4) ;
29 vol += TRIPROD(xa,ya,za,xb,yb,zb,xc,yc,zc) ;
30 xa = DA(x7,x1) ; ya = DB(x7,x1) ; za = DC(x7,x1) ;
31 xb = DA(x5,x0) ; yb = DB(x5,x0) ; zb = DC(x5,x0) ;
32 xc = DA(x7,x4)+DA(x3,x0); yc = DB(x7,x4)+DB(x3,x0); zc = DC(x7,x4)+DC(x3,x0);
33 vol += TRIPROD(xa,ya,za,xb,yb,zb,xc,yc,zc) ;
34 return vol ;
35 }
36
37 #undef MINBLUR
38 #define MINBLUR 1.234567f
39
40 /*----------------------------------------------------------------------------*/
41 /* fim = image to warp
42 wimar = contains 3 images with index warps
43 mode = interpolation method
44 */
45
GW_warpim(MRI_IMAGE * fim,MRI_IMARR * wimar,int mode)46 MRI_IMAGE * GW_warpim( MRI_IMAGE *fim, MRI_IMARR *wimar , int mode )
47 {
48 float *iar,*jar,*kar,*gar ; int npt ; MRI_IMAGE *gim ;
49
50 /* ININFO_message("Enter GW_warpim") ; */
51
52 iar = MRI_FLOAT_PTR( IMARR_SUBIM(wimar,0) ) ;
53 jar = MRI_FLOAT_PTR( IMARR_SUBIM(wimar,1) ) ;
54 kar = MRI_FLOAT_PTR( IMARR_SUBIM(wimar,2) ) ;
55 npt = fim->nvox ;
56 gim = mri_new_vol( fim->nx , fim->ny , fim->nz , MRI_float ) ;
57 gar = MRI_FLOAT_PTR(gim) ;
58
59 switch( mode ){
60 case MRI_NN: GA_interp_NN (fim,npt,iar,jar,kar,gar) ; break ;
61
62 default:
63 case MRI_LINEAR: GA_interp_linear (fim,npt,iar,jar,kar,gar) ; break ;
64
65 case MRI_CUBIC: GA_interp_cubic (fim,npt,iar,jar,kar,gar) ; break ;
66
67 case MRI_QUINTIC: GA_interp_quintic(fim,npt,iar,jar,kar,gar) ; break ;
68 }
69
70 /* ININFO_message(" exit GW_warpim") ; */
71 return gim ;
72 }
73
74 /*----------------------------------------------------------------------------*/
75
GW_blurim(MRI_IMAGE * qim,float blur)76 void GW_blurim( MRI_IMAGE *qim , float blur )
77 {
78 blur = MAX(blur,MINBLUR) ; blur = FWHM_TO_SIGMA(blur) ;
79 FIR_blur_volume( qim->nx,qim->ny,qim->nz , 1.0f,1.0f,1.0f ,
80 MRI_FLOAT_PTR(qim) , blur ) ;
81 return ;
82 }
83
84 /*----------------------------------------------------------------------------*/
85
GW_differentiate(MRI_IMAGE * iim,float blur)86 MRI_IMARR * GW_differentiate( MRI_IMAGE *iim , float blur )
87 {
88 MRI_IMARR *outar ;
89 MRI_IMAGE *qim , *qxim , *qyim , *qzim ;
90 float *qar , *qxar , *qyar , *qzar ;
91 int nx,ny,nz,nxy,nxyz , ii,jj,kk , nx1,ny1,nz1 , km,kp,jm,jp,im,ip ;
92
93 /* INFO_message("Enter GW_differentiate") ; */
94
95 nx = iim->nx ; nx1 = nx-1 ;
96 ny = iim->ny ; ny1 = ny-1 ;
97 nz = iim->nz ; nz1 = nz-1 ; nxy = nx*ny ; nxyz = nx*ny*nz ;
98
99 qim = mri_copy(iim) ; GW_blurim(qim,blur) ; qar = MRI_FLOAT_PTR(qim) ;
100
101 qxim = mri_new_vol(nx,ny,nz,MRI_float) ; qxar = MRI_FLOAT_PTR(qxim) ;
102 qyim = mri_new_vol(nx,ny,nz,MRI_float) ; qyar = MRI_FLOAT_PTR(qyim) ;
103 qzim = mri_new_vol(nx,ny,nz,MRI_float) ; qzar = MRI_FLOAT_PTR(qzim) ;
104
105 #undef IJK
106 #define IJK(i,j,k) ((i)+(j)*nx+(k)*nxy)
107 for( kk=0 ; kk < nz ; kk++ ){
108 km = kk-1 ; if( km < 0 ) km = 0 ;
109 kp = kk+1 ; if( kp > nz1 ) kp = nz1 ;
110 for( jj=0 ; jj < ny ; jj++ ){
111 jm = jj-1 ; if( jm < 0 ) jm = 0 ;
112 jp = jj+1 ; if( jp > ny1 ) jp = ny1 ;
113 for( ii=0 ; ii < nx ; ii++ ){
114 im = ii-1 ; if( im < 0 ) im = 0 ;
115 ip = ii+1 ; if( ip > nx1 ) ip = nx1 ;
116 qxar[IJK(ii,jj,kk)] = 0.5f * ( qar[IJK(ip,jj,kk)] - qar[IJK(im,jj,kk)] ) ;
117 qyar[IJK(ii,jj,kk)] = 0.5f * ( qar[IJK(ii,jp,kk)] - qar[IJK(ii,jm,kk)] ) ;
118 qzar[IJK(ii,jj,kk)] = 0.5f * ( qar[IJK(ii,jj,kp)] - qar[IJK(ii,jj,km)] ) ;
119 }}}
120 #undef IJK
121
122 INIT_IMARR(outar) ;
123 ADDTO_IMARR(outar,qim) ; ADDTO_IMARR(outar,qxim) ;
124 ADDTO_IMARR(outar,qyim) ; ADDTO_IMARR(outar,qzim) ; return outar ;
125 }
126
127 /*----------------------------------------------------------------------------*/
128 /* jim = base image
129 iimar = 4 images: [0] = source image, [1..3] = derivatives of source
130 wimar = 3 images with current index warp
131 dblur = blurring to apply to deltim
132 output = 3 images with additional index warp
133 */
134
GW_deltim(MRI_IMAGE * jim,MRI_IMARR * iimar,MRI_IMARR * wimar,float blur)135 MRI_IMARR * GW_deltim( MRI_IMAGE *jim , MRI_IMARR *iimar ,
136 MRI_IMARR *wimar, float blur )
137 {
138 MRI_IMARR *outar;
139 MRI_IMAGE *iim , *ixim , *iyim , *izim , *tim ;
140 float *iwar, *ixar , *iyar , *izar ;
141 float *iar , *jar ;
142 int nx,ny,nz,nxy,nxyz , ii,jj,kk , nx1,ny1,nz1 , qoff ;
143 float val,vmax,gg , sij,sjj ;
144
145 /* INFO_message("Enter GW_deltim") ; */
146
147 nx = jim->nx ; ny = jim->ny ; nz = jim->nz ; nxy = nx*ny ; nxyz = nxy*nz ;
148 nx1 = nx-1 ; ny1 = ny-1 ; nz1 = nz-1 ;
149
150 iim = GW_warpim( IMARR_SUBIM(iimar,0) , wimar , MRI_LINEAR ) ;
151 ixim = GW_warpim( IMARR_SUBIM(iimar,1) , wimar , MRI_LINEAR ) ;
152 iyim = GW_warpim( IMARR_SUBIM(iimar,2) , wimar , MRI_LINEAR ) ;
153 izim = GW_warpim( IMARR_SUBIM(iimar,3) , wimar , MRI_LINEAR ) ;
154
155 iar = MRI_FLOAT_PTR(iim) ;
156 ixar = MRI_FLOAT_PTR(ixim) ;
157 iyar = MRI_FLOAT_PTR(iyim) ;
158 izar = MRI_FLOAT_PTR(izim) ;
159
160 jar = MRI_FLOAT_PTR(jim) ;
161
162 /* compute scale factor to convert jar values to iar values */
163
164 sij = sjj = 0.0f ;
165 for( ii=0 ; ii < nxyz ; ii++ ){
166 sij += jar[ii]*iar[ii] ; sjj += jar[ii]*jar[ii] ;
167 }
168 if( sjj == 0.0f ) sjj = 1.0f ;
169 sij = sij / sjj ;
170
171 for( ii=0 ; ii < nxyz ; ii++ ){
172 gg = ixar[ii]*ixar[ii] + iyar[ii]*iyar[ii] + izar[ii]*izar[ii] ;
173 if( gg > 0.0f ){
174 val = sij*jar[ii] - iar[ii] ;
175 val /= ( gg + val*val ) ;
176 ixar[ii] *= val ; iyar[ii] *= val ; izar[ii] *= val ;
177 }
178 }
179 mri_free(iim) ;
180
181 tim=GW_warpim(ixim,wimar,MRI_LINEAR); mri_free(ixim); ixim=tim; ixar=MRI_FLOAT_PTR(ixim);
182 tim=GW_warpim(iyim,wimar,MRI_LINEAR); mri_free(iyim); iyim=tim; iyar=MRI_FLOAT_PTR(iyim);
183 tim=GW_warpim(izim,wimar,MRI_LINEAR); mri_free(izim); izim=tim; izar=MRI_FLOAT_PTR(izim);
184
185 GW_blurim(ixim,blur) ; GW_blurim(iyim,blur) ; GW_blurim(izim,blur) ;
186
187 vmax = 0.0f ;
188 for( ii=0 ; ii < nxyz ; ii++ ){
189 val = ixar[ii]*ixar[ii] + iyar[ii]*iyar[ii] + izar[ii]*izar[ii] ;
190 if( val > vmax ) vmax = val ;
191 }
192 if( vmax == 0.0f ) vmax = 1.0f ; /* should not happen */
193 val = 1.0f / sqrtf(vmax) ;
194 for( ii=0 ; ii < nxyz ; ii++ ){
195 ixar[ii] *= val ; iyar[ii] *= val ; izar[ii] *= val ;
196 }
197
198 INIT_IMARR(outar) ;
199 ADDTO_IMARR(outar,ixim) ; ADDTO_IMARR(outar,iyim) ; ADDTO_IMARR(outar,izim) ;
200 return outar ;
201 }
202
203 /*----------------------------------------------------------------------------*/
204 /* iim = source image
205 jim = base image
206 wimar = base index warp [3]
207 deltar = additional index warp [3]
208 ds = amount of deltar to apply
209 */
210
GW_diffval(MRI_IMAGE * iim,MRI_IMAGE * jim,MRI_IMARR * wimar,MRI_IMARR * deltar,float ds)211 float GW_diffval( MRI_IMAGE *iim , MRI_IMAGE *jim ,
212 MRI_IMARR *wimar , MRI_IMARR *deltar , float ds )
213 {
214 int ii , nxyz ; float vv,val , sij,sjj ;
215 MRI_IMAGE *wim ; float *war , *jar ;
216 MRI_IMARR *qimar ;
217
218 nxyz = iim->nx * iim->ny * iim->nz ;
219 jar = MRI_FLOAT_PTR(jim) ;
220
221 /* INFO_message("Enter GW_diffval") ; */
222
223 if( ds == 0.0f ){
224 qimar = wimar ;
225 } else {
226 MRI_IMAGE *qxim,*qyim,*qzim ;
227 float *qxar,*qyar,*qzar , *wxar,*wyar,*wzar , *dxar,*dyar,*dzar ;
228 INIT_IMARR(qimar) ;
229 qxim = mri_new_vol(iim->nx,iim->ny,iim->nz,MRI_float); qxar = MRI_FLOAT_PTR(qxim);
230 qyim = mri_new_vol(iim->nx,iim->ny,iim->nz,MRI_float); qyar = MRI_FLOAT_PTR(qyim);
231 qzim = mri_new_vol(iim->nx,iim->ny,iim->nz,MRI_float); qzar = MRI_FLOAT_PTR(qzim);
232 wxar = MRI_FLOAT_PTR(IMARR_SUBIM(wimar ,0)) ;
233 wyar = MRI_FLOAT_PTR(IMARR_SUBIM(wimar ,1)) ;
234 wzar = MRI_FLOAT_PTR(IMARR_SUBIM(wimar ,2)) ;
235 dxar = MRI_FLOAT_PTR(IMARR_SUBIM(deltar,0)) ;
236 dyar = MRI_FLOAT_PTR(IMARR_SUBIM(deltar,1)) ;
237 dzar = MRI_FLOAT_PTR(IMARR_SUBIM(deltar,2)) ;
238 for( ii=0 ; ii < nxyz ; ii++ ){
239 qxar[ii] = wxar[ii] + ds*dxar[ii] ;
240 qyar[ii] = wyar[ii] + ds*dyar[ii] ;
241 qzar[ii] = wzar[ii] + ds*dzar[ii] ;
242 }
243 ADDTO_IMARR(qimar,qxim) ; ADDTO_IMARR(qimar,qyim) ; ADDTO_IMARR(qimar,qzim) ;
244 }
245
246 wim = GW_warpim( iim , qimar , MRI_LINEAR ) ; war = MRI_FLOAT_PTR(wim) ;
247 if( qimar != wimar ) DESTROY_IMARR(qimar) ;
248
249 sij = sjj = 0.0f ;
250 for( ii=0 ; ii < nxyz ; ii++ ){
251 sij += jar[ii]*war[ii] ; sjj += jar[ii]*jar[ii] ;
252 }
253 /* ININFO_message(" compute sij=%g,sjj=%g",sij,sjj) ; */
254 if( sjj == 0.0f ) sjj = 1.0f ;
255 sij = sij / sjj ;
256
257 vv = 0.0f ; jar = MRI_FLOAT_PTR(jim) ; war = MRI_FLOAT_PTR(wim) ;
258 for( ii=0 ; ii < nxyz ; ii++ ){
259 val = sij*jar[ii] - war[ii] ; vv += val*val ;
260 }
261 mri_free(wim) ;
262
263 /* ININFO_message("GW_diffval(%f) = %g",ds,vv) ; */
264 return vv ;
265 }
266
267 /*----------------------------------------------------------------------------*/
268
GW_best_step(MRI_IMAGE * iim,MRI_IMAGE * jim,MRI_IMARR * wimar,MRI_IMARR * deltar)269 float GW_best_step( MRI_IMAGE *iim , MRI_IMAGE *jim ,
270 MRI_IMARR *wimar , MRI_IMARR *deltar )
271 {
272 float vbot , vtop , vmid , db=0.0f, dt=0.16f , dm , vbm,vmt ;
273 int ii ;
274
275 /* INFO_message("Enter GW_best_step") ; */
276
277 vbot = GW_diffval( iim , jim , wimar , deltar , db ) ;
278 vtop = GW_diffval( iim , jim , wimar , deltar , dt ) ;
279
280 #if 0
281 for( ii=0 ; ii < 2 ; ii++ ){
282 dm = 0.5f*(db+dt) ;
283 vmid = GW_diffval( iim , jim , wimar , deltar , dm ) ;
284 vbm = MIN(vbot,vmid) ;
285 vmt = MIN(vmid,vtop) ;
286 if( vbm <= vmt ){ dt = dm ; vtop = vmid ; }
287 else { db = dm ; vbot = vmid ; }
288 }
289 #endif
290
291 return ((vtop < vbot) ? dt : db) ;
292 }
293
294 /*----------------------------------------------------------------------------*/
295 /* Scale jim (in place) to match iim */
296
GW_scaleim(MRI_IMAGE * iim,MRI_IMAGE * jim)297 float GW_scaleim( MRI_IMAGE *iim , MRI_IMAGE *jim )
298 {
299 float iclip , jclip , aa ; int ii , nvox ;
300 float *jar ;
301
302 iclip = THD_cliplevel( iim , 0.4f ) ;
303 jclip = THD_cliplevel( iim , 0.4f ) ;
304 if( jclip <= 0.0f || iclip <= 0.0f ) return 0.0f ;
305
306 aa = iclip / jclip ;
307 nvox = jim->nvox ;
308 jar = MRI_FLOAT_PTR(jim) ;
309 for( ii=0 ; ii < nvox ; ii++ ) jar[ii] *= aa ;
310 return aa ;
311 }
312
313 /*----------------------------------------------------------------------------*/
314
315 static float blur = 0.0f ;
316 static MRI_IMARR *iimar = NULL ;
317 static MRI_IMARR *wimar = NULL ;
318 static MRI_IMAGE *jim = NULL ;
319
320 static THD_3dim_dataset *qset = NULL ;
321
GW_setup(MRI_IMAGE * inim,MRI_IMAGE * bsim,float iblur,float jblur)322 void GW_setup( MRI_IMAGE *inim , MRI_IMAGE *bsim , float iblur , float jblur )
323 {
324 MRI_IMAGE *wim ; float *war ;
325 int nx,ny,nz,nxy,nxyz , ii,jj,kk,qq ;
326
327 nx = inim->nx ;
328 ny = inim->ny ;
329 nz = inim->nz ; nxy = nx*ny ; nxyz = nx*ny*nz ;
330
331 if( iimar != NULL ) DESTROY_IMARR(iimar) ;
332 if( wimar != NULL ) DESTROY_IMARR(wimar) ;
333 if( jim != NULL ) mri_free(jim) ;
334
335 blur = MAX(iblur,MINBLUR) ;
336 iimar = GW_differentiate( inim , blur ) ;
337 jim = mri_copy(bsim) ;
338 blur = MAX(jblur,MINBLUR) ;
339 GW_blurim(jim,blur) ;
340 (void)GW_scaleim( IMARR_SUBIM(iimar,0) , jim ) ;
341
342 { MRI_IMAGE *tim = mri_copy(jim) ;
343 EDIT_dset_items( qset , ADN_prefix,"zzbase" , ADN_none ) ;
344 EDIT_substitute_brick( qset , 0 , MRI_float , MRI_FLOAT_PTR(tim) ) ;
345 DSET_write(qset) ; WROTE_DSET(qset) ;
346 }
347
348 INIT_IMARR(wimar) ;
349 wim = mri_new_vol(jim->nx,jim->ny,jim->nz,MRI_float); ADDTO_IMARR(wimar,wim);
350 war = MRI_FLOAT_PTR(wim) ;
351 for( qq=kk=0 ; kk < nz ; kk++ )
352 for( jj=0 ; jj < ny ; jj++ )
353 for( ii=0 ; ii < nx ; ii++,qq++ ) war[qq] = ii ;
354
355 wim = mri_new_vol(jim->nx,jim->ny,jim->nz,MRI_float); ADDTO_IMARR(wimar,wim);
356 war = MRI_FLOAT_PTR(wim) ;
357 for( qq=kk=0 ; kk < nz ; kk++ )
358 for( jj=0 ; jj < ny ; jj++ )
359 for( ii=0 ; ii < nx ; ii++,qq++ ) war[qq] = jj ;
360
361 wim = mri_new_vol(jim->nx,jim->ny,jim->nz,MRI_float); ADDTO_IMARR(wimar,wim);
362 war = MRI_FLOAT_PTR(wim) ;
363 for( qq=kk=0 ; kk < nz ; kk++ )
364 for( jj=0 ; jj < ny ; jj++ )
365 for( ii=0 ; ii < nx ; ii++,qq++ ) war[qq] = kk ;
366
367 return ;
368 }
369
370 /*----------------------------------------------------------------------------*/
371
GW_iterate(void)372 float GW_iterate(void)
373 {
374 MRI_IMARR *deltar ; float ds , qblur ;
375
376 /* INFO_message("Enter GW_iterate") ; */
377
378 qblur = 1.444f*blur ; qblur = MAX(qblur,14.444f) ;
379 deltar = GW_deltim( jim , iimar , wimar , qblur ) ;
380 ds = GW_best_step( IMARR_SUBIM(iimar,0) , jim , wimar , deltar ) ;
381
382 ININFO_message("best_step = %f",ds) ;
383
384 if( ds != 0.0f ){
385 float *wxar = MRI_FLOAT_PTR(IMARR_SUBIM(wimar ,0)) ;
386 float *wyar = MRI_FLOAT_PTR(IMARR_SUBIM(wimar ,1)) ;
387 float *wzar = MRI_FLOAT_PTR(IMARR_SUBIM(wimar ,2)) ;
388 float *dxar = MRI_FLOAT_PTR(IMARR_SUBIM(deltar,0)) ;
389 float *dyar = MRI_FLOAT_PTR(IMARR_SUBIM(deltar,1)) ;
390 float *dzar = MRI_FLOAT_PTR(IMARR_SUBIM(deltar,2)) ;
391 int ii , nxyz = jim->nvox ;
392 for( ii=0 ; ii < nxyz ; ii++ ){
393 wxar[ii] += ds*dxar[ii] ;
394 wyar[ii] += ds*dyar[ii] ;
395 wzar[ii] += ds*dzar[ii] ;
396 }
397
398 /* debug output */
399
400 { MRI_IMAGE *tim = GW_warpim(IMARR_SUBIM(iimar,0),wimar,MRI_LINEAR) ;
401 char prefix[64] ;
402 static int nite=0 ;
403 nite++ ; sprintf(prefix,"zziter%04d",nite) ;
404 EDIT_dset_items( qset , ADN_prefix,prefix , ADN_none ) ;
405 EDIT_substitute_brick( qset , 0 , MRI_float , MRI_FLOAT_PTR(tim) ) ;
406 DSET_write(qset) ; WROTE_DSET(qset) ; mri_clear_data_pointer(tim) ; mri_free(tim) ;
407 }
408 }
409
410 DESTROY_IMARR(deltar) ; return ds ;
411 }
412
413 /*----------------------------------------------------------------------------*/
414
main(int argc,char * argv[])415 int main( int argc , char *argv[] )
416 {
417 THD_3dim_dataset *iset , *bset , *oset ;
418 MRI_IMAGE *iim , *bim ; int ii,iarg=1 ; float ds , bblur=0.0f,sblur=4.9f ;
419
420 if( argc < 3 ){ printf("gwarp [-blur b s] base source\n") ; exit(0) ; }
421
422 INFO_message("gwarp is NOT a production program: it is merely for testing!") ;
423
424 if( strcmp(argv[iarg],"-blur") == 0 ){
425 bblur = strtod(argv[++iarg],NULL) ;
426 sblur = strtod(argv[++iarg],NULL) ; iarg++ ;
427 }
428
429 bset = THD_open_dataset(argv[iarg++]) ; if( !ISVALID_DSET(bset) ) ERROR_exit("bset") ;
430 DSET_load(bset) ; bim = mri_to_float(DSET_BRICK(bset,0)) ; DSET_unload(bset) ;
431
432 iset = THD_open_dataset(argv[iarg++]) ; if( !ISVALID_DSET(iset) ) ERROR_exit("iset") ;
433 DSET_load(iset) ; iim = mri_to_float(DSET_BRICK(iset,0)) ; DSET_unload(iset) ;
434
435 qset = EDIT_empty_copy(iset) ;
436 EDIT_dset_items( qset , ADN_nvals,1 , ADN_none ) ;
437
438 GW_setup(iim,bim,sblur,bblur) ; mri_free(bim) ;
439
440 for( ii=0 ; ii < 499 ; ii++ ){
441 ds = GW_iterate() ; if( ds == 0.0f ) break ;
442 }
443
444 bim = GW_warpim(iim,wimar,MRI_LINEAR) ;
445
446 oset = EDIT_empty_copy(iset) ;
447 EDIT_dset_items( oset , ADN_prefix,"gwarp" , ADN_nvals,1 , ADN_none ) ;
448 EDIT_substitute_brick( oset , 0 , MRI_float , MRI_FLOAT_PTR(bim) ) ;
449 DSET_write(oset) ; WROTE_DSET(oset) ; exit(0) ;
450 }
451