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 "afni.h"
8
9 #ifndef ALLOW_PLUGINS
10 # error "Plugins not properly set up -- see machdep.h"
11 #endif
12
13 /***********************************************************************
14 Plugin to do what 3dvolreg.c does. Most of the code is
15 adapted from plug_imreg.c, 3dvolreg.c, and plug_realtime.c.
16 ************************************************************************/
17
18 /*--------------------- string to 'help' the user --------------------*/
19
20 static char helpstring[] =
21 "Purpose: 3D (volumetric) registration of 3D+time dataset\n"
22 "[See also program '3dvolreg']\n"
23 "\n"
24 "Input items to this plugin are:\n"
25 " Datasets: Input = 3D+time dataset to process\n"
26 " Output = Prefix for new dataset\n"
27 "\n"
28 " Parameters: Base Brick = Time index for base image\n"
29 " Resampling = Resampling method\n"
30 "\n"
31 " Basis: Dataset = Dataset from which to take base brick\n"
32 " [if not selected, uses input dataset]\n"
33 "\n"
34 " Outputs: dfile = If not blank, name of file in which to\n"
35 " save estimated movement parameters plus\n"
36 " before-and-after residuals;\n"
37 " see '3dvolreg -help' for details\n"
38 " 1Dfile = If not blank, name of file in which to save\n"
39 " ONLY the estimated movement parameters;\n"
40 " this file will also be added to the AFNI\n"
41 " timeseries list, so that it can be used\n"
42 " as an 'ort' with FIM.\n"
43 " Graph = If 'Yes', will plot the estimated movement\n"
44 " parameters when the registration is done.\n"
45 "\n"
46 "WARNING: This plugin can be slow, particularly when registering\n"
47 " large datasets using Fourier resampling. The algorithm\n"
48 " used will only work to correct for small movements.\n"
49 "\n"
50 "Author -- RW Cox, November 1998"
51 ;
52
53 /*----------------- prototypes for internal routines -----------------*/
54
55 char * VOLREG_main( PLUGIN_interface * ) ; /* the entry point */
56
57 /*---------------------------- global data ---------------------------*/
58
59 static PLUGIN_interface * global_plint = NULL ;
60
61 #define NRESAM 4
62 static char * REG_resam_strings[NRESAM] = {
63 "Cubic" , "Quintic" , "Heptic" , "Fourier"
64 } ;
65
66 static int REG_resam_ints[NRESAM] = {
67 MRI_CUBIC , MRI_QUINTIC , MRI_HEPTIC , MRI_FOURIER
68 } ;
69
70 #define NGRAPH 2
71 static char * GRAPH_strings[NGRAPH] = { "No" , "Yes" } ;
72
73 static int VL_nbase = 3 ;
74 static int VL_resam = 3 ;
75 static int VL_clipit = 1 ;
76 static int VL_intern = 1 ;
77 static int VL_graph = 0 ;
78 static MRI_IMAGE * VL_imbase = NULL ;
79
80 static char VL_dfile[256] = "\0" ;
81 static char VL_1Dfile[256] = "\0" ; /* 14 Apr 2000 */
82 static char VL_prefix[THD_MAX_PREFIX] = "\0" ;
83
84 static THD_3dim_dataset * VL_dset = NULL ;
85
86 static int VL_maxite = 9 ;
87 static float VL_dxy = 0.05 ; /* voxels */
88 static float VL_dph = 0.07 ; /* degrees */
89 static float VL_del = 0.70 ; /* voxels */
90
91 /***********************************************************************
92 Set up the interface to the user:
93 1) Create a new interface using "PLUTO_new_interface";
94
95 2) For each line of inputs, create the line with "PLUTO_add_option"
96 (this line of inputs can be optional or mandatory);
97
98 3) For each item on the line, create the item with
99 "PLUTO_add_dataset" for a dataset chooser,
100 "PLUTO_add_string" for a string chooser,
101 "PLUTO_add_number" for a number chooser.
102 ************************************************************************/
103
104
105 DEFINE_PLUGIN_PROTOTYPE
106
PLUGIN_init(int ncall)107 PLUGIN_interface * PLUGIN_init( int ncall )
108 {
109 PLUGIN_interface * plint ; /* will be the output of this routine */
110
111 if( ncall > 0 ) return NULL ; /* only one interface */
112 CHECK_IF_ALLOWED("3DREGISTRATION","3D Registration") ; /* 30 Sep 2016 */
113
114 /*---------------- set titles and call point ----------------*/
115
116 plint = PLUTO_new_interface( "3D Registration" ,
117 "3D Registration of a 3D+time Dataset" ,
118 helpstring ,
119 PLUGIN_CALL_VIA_MENU , VOLREG_main ) ;
120
121 PLUTO_add_hint( plint , "3D Registration of a 3D+time Dataset" ) ;
122
123 global_plint = plint ; /* make global copy */
124
125 PLUTO_set_sequence( plint , "A:newdset:reg" ) ;
126
127 /*--------- 1st line ---------*/
128
129 PLUTO_add_option( plint ,
130 "Datasets" , /* label at left of input line */
131 "DAtasets" , /* tag to return to plugin */
132 TRUE /* is this mandatory? */
133 ) ;
134
135 PLUTO_add_dataset( plint ,
136 "Input" , /* label next to button */
137 ANAT_ALL_MASK , /* take any anat datasets */
138 FUNC_FIM_MASK , /* only allow fim funcs */
139 DIMEN_ALL_MASK | BRICK_ALLREAL_MASK
140 ) ;
141
142 PLUTO_add_string( plint ,
143 "Output" , /* label next to textfield */
144 0,NULL , /* no fixed strings to choose among */
145 19 /* 19 spaces for typing in value */
146 ) ;
147
148 /*---------- 2nd line --------*/
149
150 PLUTO_add_option( plint ,
151 "Parameters" , /* label at left of input line */
152 "Parameters" , /* tag to return to plugin */
153 TRUE /* is this mandatory? */
154 ) ;
155
156 PLUTO_add_number( plint ,
157 "Base Brick" , /* label next to chooser */
158 0 , /* smallest possible value */
159 98 , /* largest possible value */
160 0 , /* decimal shift (none in this case) */
161 VL_nbase , /* default value */
162 FALSE /* allow user to edit value? */
163 ) ;
164
165 PLUTO_add_string( plint, "Resampling", NRESAM, REG_resam_strings, VL_resam ) ;
166
167 /*---------- 3rd line --------*/
168
169 PLUTO_add_option( plint , "Basis" , "Basis" , FALSE ) ;
170
171 PLUTO_add_dataset( plint ,
172 "Dataset" , /* label next to button */
173 ANAT_ALL_MASK , /* take any anat datasets */
174 FUNC_FIM_MASK , /* only allow fim funcs */
175 DIMEN_ALL_MASK | BRICK_ALLREAL_MASK
176 ) ;
177
178 /*---------- 4th line ---------*/
179
180 PLUTO_add_option( plint , "Outputs" , "Outputs" , FALSE ) ;
181
182 PLUTO_add_string( plint ,
183 "dfile" , /* label next to textfield */
184 0,NULL , /* no fixed strings to choose among */
185 19 /* 19 spaces for typing in value */
186 ) ;
187
188 /* 14 Apr 2000: add 1Dfile */
189
190 PLUTO_add_string( plint ,
191 "1Dfile" , /* label next to textfield */
192 0,NULL , /* no fixed strings to choose among */
193 19 /* 19 spaces for typing in value */
194 ) ;
195
196 PLUTO_add_string( plint , "Graph" , NGRAPH , GRAPH_strings , VL_graph ) ;
197
198 /*--------- done with interface setup ---------*/
199
200 return plint ;
201 }
202
203 /***************************************************************************
204 Main routine for this plugin (will be called from AFNI).
205 If the return string is not NULL, some error transpired, and
206 AFNI will popup the return string in a message box.
207 ****************************************************************************/
208
209 #undef FREEUP
210 #define FREEUP(x) do{if((x) != NULL){free((x)); (x)=NULL;}}while(0)
211
VOLREG_main(PLUGIN_interface * plint)212 char * VOLREG_main( PLUGIN_interface * plint )
213 {
214 MCW_idcode * idc ;
215 char * str ;
216 MRI_3dalign_basis * albase ;
217 THD_3dim_dataset * new_dset ;
218 MRI_IMAGE * qim , * tim , * fim ;
219 float *dx, *dy, *dz, *roll, *yaw, *pitch, *rmsnew, *rmsold, *imb, *tar ;
220 float ddx,ddy,ddz , sum ;
221 float dxtop=0.0,dytop=0.0,dztop=0.0 , rolltop=0.0,yawtop=0.0,pitchtop=0.0 ;
222 float dxbot=0.0,dybot=0.0,dzbot=0.0 , rollbot=0.0,yawbot=0.0,pitchbot=0.0 ;
223 float dxbar,dybar,dzbar , rollbar,yawbar,pitchbar ;
224 int kim,ii , imcount , iha , ax1,ax2,ax3 , hax1,hax2,hax3 ;
225
226 double cputim = COX_cpu_time();
227
228 /*--------------------------------------------------------------------*/
229 /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
230
231 /*--------- go to first input line ---------*/
232
233 PLUTO_next_option(plint) ;
234
235 idc = PLUTO_get_idcode(plint) ; /* get dataset item */
236 VL_dset = PLUTO_find_dset(idc) ; /* get ptr to dataset */
237 if( VL_dset == NULL )
238 return "*************************\n"
239 "Cannot find Input Dataset\n"
240 "*************************" ;
241
242 ii = DSET_NVALS_PER_TIME(VL_dset) ;
243 if( ii > 1 )
244 return "************************************\n"
245 "Dataset has > 1 value per time point\n"
246 "************************************" ;
247
248 str = PLUTO_get_string(plint) ; /* get the output prefix */
249 if( ! PLUTO_prefix_ok(str) ) /* check if it is OK */
250 return "************************\n"
251 "Output Prefix is illegal\n"
252 "************************" ;
253 strcpy( VL_prefix , str ) ;
254
255 /*--------- go to next input line ---------*/
256
257 PLUTO_next_option(plint) ;
258
259 VL_nbase = PLUTO_get_number(plint) ;
260
261 str = PLUTO_get_string(plint) ;
262 VL_resam = PLUTO_string_index( str , NRESAM , REG_resam_strings ) ;
263 VL_resam = REG_resam_ints[VL_resam] ;
264
265 /*--------- see if the other (non-mandatory) option line are present --------*/
266
267 VL_imbase = NULL ;
268 VL_dfile[0] = '\0' ;
269 VL_1Dfile[0]= '\0' ; /* 14 Apr 2000 */
270 VL_graph = 0 ;
271
272 while( 1 ){
273 str = PLUTO_get_optiontag( plint ) ; if( str == NULL ) break ;
274
275 if( strcmp(str,"Basis") == 0 ){ /* Alternate basis dataset */
276 THD_3dim_dataset * bset ;
277
278 idc = PLUTO_get_idcode(plint) ; /* get dataset item */
279 bset = PLUTO_find_dset(idc) ; /* get ptr to dataset */
280 if( bset == NULL )
281 return "*************************\n"
282 "Cannot find Basis Dataset\n"
283 "*************************" ;
284
285 if( VL_nbase >= DSET_NVALS(bset) || VL_nbase < 0 )
286 return "*****************************************\n"
287 "Base index out of range for Basis dataset\n"
288 "*****************************************" ;
289
290 if( DSET_NX(bset) != DSET_NX(VL_dset) ||
291 DSET_NY(bset) != DSET_NY(VL_dset) || DSET_NZ(bset) != DSET_NZ(VL_dset) )
292 return "***************************\n"
293 "Brick size mismatch between\n"
294 " Input and Basis datasets! \n"
295 "***************************" ;
296
297 DSET_load(bset) ;
298 if( ! DSET_LOADED(bset) ) return "************************\n"
299 "Can't load Basis dataset\n"
300 "************************" ;
301
302 VL_imbase = mri_to_float( DSET_BRICK(bset,VL_nbase) ) ; /* copy this */
303 VL_intern = 0 ;
304 DSET_unload(bset) ;
305
306 } else if( strcmp(str,"Outputs") == 0 ){ /* Optional outputs */
307
308 str = PLUTO_get_string(plint) ;
309 if( str != NULL && str[0] != '\0' ){
310 if( THD_filename_ok(str) ) strcpy( VL_dfile , str ) ;
311 else return "*************************\n"
312 "dfile name not acceptable\n"
313 "*************************" ;
314 }
315
316 /* 14 Apr 2000: get 1Dfile */
317
318 str = PLUTO_get_string(plint) ;
319 if( str != NULL && str[0] != '\0' ){
320 if( THD_filename_ok(str) ) strcpy( VL_1Dfile , str ) ;
321 else return "**************************\n"
322 "1Dfile name not acceptable\n"
323 "**************************" ;
324 }
325
326 str = PLUTO_get_string(plint) ;
327 VL_graph = PLUTO_string_index( str , NGRAPH , GRAPH_strings ) ;
328 }
329
330 }
331
332 /*-- must manufacture base image --*/
333
334 if( VL_imbase == NULL ){
335
336 if( DSET_NVALS(VL_dset) < 2 )
337 return "******************************************\n"
338 "Can't register a 1 brick dataset to itself\n"
339 "******************************************" ;
340
341 if( VL_nbase >= DSET_NVALS(VL_dset) || VL_nbase < 0 )
342 return "*****************************************\n"
343 "Base index out of range for Input dataset\n"
344 "*****************************************" ;
345
346 DSET_load(VL_dset) ;
347 if( ! DSET_LOADED(VL_dset) ) return "************************\n"
348 "Can't load Input dataset\n"
349 "************************" ;
350
351 VL_imbase = mri_to_float( DSET_BRICK(VL_dset,VL_nbase) ) ; /* copy this */
352 VL_intern = 1 ;
353 }
354
355 VL_imbase->dx = fabs( DSET_DX(VL_dset) ) ; /* must set the voxel dimensions */
356 VL_imbase->dy = fabs( DSET_DY(VL_dset) ) ;
357 VL_imbase->dz = fabs( DSET_DZ(VL_dset) ) ;
358 imb = MRI_FLOAT_PTR( VL_imbase ) ; /* need this to compute rms */
359
360 /*------------- set up to compute new dataset -----------*/
361
362 DSET_load( VL_dset ) ;
363 if( ! DSET_LOADED(VL_dset) ){
364 mri_free( VL_imbase ) ;
365 return "************************\n"
366 "Can't load Input dataset\n"
367 "************************" ;
368 }
369
370 imcount = DSET_NVALS( VL_dset ) ;
371 dx = (float *) malloc( sizeof(float) * imcount ) ;
372 dy = (float *) malloc( sizeof(float) * imcount ) ;
373 dz = (float *) malloc( sizeof(float) * imcount ) ;
374 roll = (float *) malloc( sizeof(float) * imcount ) ;
375 pitch = (float *) malloc( sizeof(float) * imcount ) ;
376 yaw = (float *) malloc( sizeof(float) * imcount ) ;
377 rmsnew = (float *) malloc( sizeof(float) * imcount ) ;
378 rmsold = (float *) malloc( sizeof(float) * imcount ) ;
379
380 iha = THD_handedness( VL_dset ) ; /* LH or RH? */
381 ax1 = THD_axcode( VL_dset , 'I' ) ; hax1 = ax1 * iha ; /* roll */
382 ax2 = THD_axcode( VL_dset , 'R' ) ; hax2 = ax2 * iha ; /* pitch */
383 ax3 = THD_axcode( VL_dset , 'A' ) ; hax3 = ax3 * iha ; /* yaw */
384
385 mri_3dalign_params( VL_maxite , VL_dxy , VL_dph , VL_del ,
386 abs(ax1)-1 , abs(ax2)-1 , abs(ax3)-1 , -1 ) ;
387
388 mri_3dalign_method( VL_resam , 0 , 0 , VL_clipit ) ;
389
390 new_dset = EDIT_empty_copy( VL_dset ) ;
391 EDIT_dset_items( new_dset , ADN_prefix , VL_prefix , ADN_none ) ;
392
393 { char * his = PLUTO_commandstring(plint) ;
394 tross_Copy_History( VL_dset , new_dset ) ;
395 tross_Append_History( new_dset , his ) ; free(his) ;
396 }
397
398 albase = mri_3dalign_setup( VL_imbase , NULL ) ;
399
400 /*-- loop over sub-bricks and register them --*/
401
402 PLUTO_popup_meter(plint) ;
403
404 dxbar = dybar = dzbar = rollbar = yawbar = pitchbar = 0.0 ;
405
406 for( kim=0 ; kim < imcount ; kim++ ){
407
408 qim = DSET_BRICK( VL_dset , kim ) ; /* the sub-brick in question */
409 fim = mri_to_float( qim ) ; /* make a float copy */
410 fim->dx = fabs( DSET_DX(VL_dset) ) ; /* must set voxel dimensions */
411 fim->dy = fabs( DSET_DY(VL_dset) ) ;
412 fim->dz = fabs( DSET_DZ(VL_dset) ) ;
413
414 /*-- the actual registration --*/
415
416 if( !VL_intern || kim != VL_nbase ){
417 tim = mri_3dalign_one( albase , fim ,
418 roll+kim , pitch+kim , yaw+kim ,
419 &ddx , &ddy , &ddz ) ;
420 } else {
421 tim = mri_to_float( VL_imbase ) ;
422 roll[kim] = pitch[kim] = yaw[kim] = ddx = ddy = ddz = 0.0 ;
423 }
424
425 /*-- need to massage output parameters --*/
426
427 roll[kim] *= (180.0/PI) ; if( hax1 < 0 ) roll[kim] = -roll[kim] ;
428 pitch[kim] *= (180.0/PI) ; if( hax2 < 0 ) pitch[kim] = -pitch[kim] ;
429 yaw[kim] *= (180.0/PI) ; if( hax3 < 0 ) yaw[kim] = -yaw[kim] ;
430
431 switch( new_dset->daxes->xxorient ){
432 case ORI_R2L_TYPE: dy[kim] = ddx ; break ;
433 case ORI_L2R_TYPE: dy[kim] = -ddx ; break ;
434 case ORI_P2A_TYPE: dz[kim] = -ddx ; break ;
435 case ORI_A2P_TYPE: dz[kim] = ddx ; break ;
436 case ORI_I2S_TYPE: dx[kim] = ddx ; break ;
437 case ORI_S2I_TYPE: dx[kim] = -ddx ; break ;
438 }
439
440 switch( new_dset->daxes->yyorient ){
441 case ORI_R2L_TYPE: dy[kim] = ddy ; break ;
442 case ORI_L2R_TYPE: dy[kim] = -ddy ; break ;
443 case ORI_P2A_TYPE: dz[kim] = -ddy ; break ;
444 case ORI_A2P_TYPE: dz[kim] = ddy ; break ;
445 case ORI_I2S_TYPE: dx[kim] = ddy ; break ;
446 case ORI_S2I_TYPE: dx[kim] = -ddy ; break ;
447 }
448
449 switch( new_dset->daxes->zzorient ){
450 case ORI_R2L_TYPE: dy[kim] = ddz ; break ;
451 case ORI_L2R_TYPE: dy[kim] = -ddz ; break ;
452 case ORI_P2A_TYPE: dz[kim] = -ddz ; break ;
453 case ORI_A2P_TYPE: dz[kim] = ddz ; break ;
454 case ORI_I2S_TYPE: dx[kim] = ddz ; break ;
455 case ORI_S2I_TYPE: dx[kim] = -ddz ; break ;
456 }
457
458 /*-- collect statistics --*/
459
460 if( kim == 0 ){
461 dxtop = dxbot = dx[kim] ;
462 dytop = dybot = dy[kim] ;
463 dztop = dzbot = dz[kim] ;
464 rolltop = rollbot = roll[kim] ;
465 pitchtop = pitchbot = pitch[kim] ;
466 yawtop = yawbot = yaw[kim] ;
467 } else {
468 dxtop = MAX(dxtop , dx[kim] ); dxbot = MIN(dxbot , dx[kim] );
469 dytop = MAX(dytop , dy[kim] ); dybot = MIN(dybot , dy[kim] );
470 dztop = MAX(dztop , dz[kim] ); dzbot = MIN(dzbot , dz[kim] );
471 rolltop = MAX(rolltop , roll[kim] ); rollbot = MIN(rollbot , roll[kim] );
472 pitchtop = MAX(pitchtop, pitch[kim]); pitchbot = MIN(pitchbot, pitch[kim]);
473 yawtop = MAX(yawtop , yaw[kim] ); yawbot = MIN(yawbot , yaw[kim] );
474 }
475
476 dxbar += dx[kim] ; dybar += dy[kim] ; dzbar += dz[kim] ;
477 rollbar += roll[kim] ; pitchbar += pitch[kim] ; yawbar += yaw[kim] ;
478
479 sum = 0.0 ;
480 tar = MRI_FLOAT_PTR(tim) ;
481 for( ii=0 ; ii < tim->nvox ; ii++ ) sum += SQR( imb[ii] - tar[ii] ) ;
482 rmsnew[kim] = sqrt( sum / tim->nvox ) ;
483
484 sum = 0.0 ;
485 tar = MRI_FLOAT_PTR(fim) ;
486 for( ii=0 ; ii < fim->nvox ; ii++ ) sum += SQR( imb[ii] - tar[ii] ) ;
487 rmsold[kim] = sqrt( sum / fim->nvox ) ;
488
489 mri_free(fim) ; /* only needed this to compute rmsold */
490
491 /*-- Attach the registered brick to output dataset,
492 converting it to the correct type, if necessary
493 (the new registered brick in "tim" is stored as floats). --*/
494
495 switch( qim->kind ){
496
497 case MRI_float:
498 EDIT_substitute_brick( new_dset , kim , MRI_float , MRI_FLOAT_PTR(tim) ) ;
499 mri_fix_data_pointer( NULL , tim ) ; mri_free( tim ) ;
500 break ;
501
502 case MRI_short:
503 fim = mri_to_short(1.0,tim) ; mri_free( tim ) ;
504 EDIT_substitute_brick( new_dset , kim , MRI_short , MRI_SHORT_PTR(fim) ) ;
505 mri_fix_data_pointer( NULL , fim ) ; mri_free( fim ) ;
506 break ;
507
508 case MRI_byte:
509 fim = mri_to_byte(tim) ; mri_free( tim ) ;
510 EDIT_substitute_brick( new_dset , kim , MRI_byte , MRI_BYTE_PTR(fim) ) ;
511 mri_fix_data_pointer( NULL , fim ) ; mri_free( fim ) ;
512 break ;
513
514 /*-- should not ever get here, but who knows? --*/
515
516 default:
517 mri_free( VL_imbase ) ;
518 FREEUP(dx) ; FREEUP(dy) ; FREEUP(dz) ;
519 FREEUP(roll) ; FREEUP(yaw) ; FREEUP(pitch) ;
520 FREEUP(rmsnew) ; FREEUP(rmsold) ;
521
522 fprintf(stderr,"\n*** Can't register bricks of type %s\n",
523 MRI_TYPE_name[qim->kind] ) ;
524 return("*** Can't register bricks of this type");
525 /* EXIT(1) ;*/
526 }
527
528 DSET_unload_one( VL_dset , kim ) ; /* don't need this anymore */
529
530 PLUTO_set_meter(plint, (100*(kim+1))/imcount ) ;
531
532 } /* end of loop over sub-bricks */
533
534 /*-- done with registration --*/
535
536 mri_3dalign_cleanup( albase ) ;
537 mri_free( VL_imbase ) ;
538 DSET_unload( VL_dset ) ;
539
540 /*-- show some statistics of the results --*/
541
542 dxbar /= imcount ; dybar /= imcount ; dzbar /= imcount ;
543 rollbar /= imcount ; pitchbar /= imcount ; yawbar /= imcount ;
544
545 { char buf[512] ; int nbuf ;
546
547 cputim = COX_cpu_time() - cputim ;
548 sprintf(buf," \nCPU time for realignment=%.3g s\n\n" , cputim) ;
549 nbuf = strlen(buf) ;
550
551 sprintf(buf+nbuf,"Min : roll=%+.3f pitch=%+.3f yaw=%+.3f"
552 " dS=%+.3f dL=%+.3f dP=%+.3f\n\n" ,
553 rollbot , pitchbot , yawbot , dxbot , dybot , dzbot ) ;
554 nbuf = strlen(buf) ;
555
556 sprintf(buf+nbuf,"Mean: roll=%+.3f pitch=%+.3f yaw=%+.3f"
557 " dS=%+.3f dL=%+.3f dP=%+.3f\n\n" ,
558 rollbar , pitchbar , yawbar , dxbar , dybar , dzbar ) ;
559 nbuf = strlen(buf) ;
560
561 sprintf(buf+nbuf,"Max : roll=%+.3f pitch=%+.3f yaw=%+.3f"
562 " dS=%+.3f dL=%+.3f dP=%+.3f\n" ,
563 rolltop , pitchtop , yawtop , dxtop , dytop , dztop ) ;
564
565 PLUTO_popup_message( plint , buf ) ;
566 }
567
568 /*-- save to a file --*/
569
570 if( VL_dfile[0] != '\0' ){
571 FILE * fp ;
572
573 if( THD_is_file(VL_dfile) )
574 PLUTO_popup_transient( plint , "** Warning:\n"
575 "** Overwriting dfile" ) ;
576
577 fp = fopen( VL_dfile , "w" ) ;
578 for( kim=0 ; kim < imcount ; kim++ )
579 fprintf(fp , "%4d %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %11.4g %11.4g\n" ,
580 kim , roll[kim], pitch[kim], yaw[kim],
581 dx[kim], dy[kim], dz[kim], rmsold[kim] , rmsnew[kim] ) ;
582 fclose(fp) ;
583 }
584
585 if( VL_1Dfile[0] != '\0' ){ /* 14 Apr 2000 */
586 FILE * fp ;
587 char fn[256] ;
588 MRI_IMAGE * tsim ;
589 float * tsar ;
590
591 strcpy(fn,VL_1Dfile) ;
592 if( strstr(VL_1Dfile,"1D") == NULL ) strcat(fn,".1D") ;
593
594 if( THD_is_file(fn) )
595 PLUTO_popup_transient( plint , "** Warning:\n"
596 "** Overwriting 1Dfile" ) ;
597
598 tsim = mri_new( imcount , 6 , MRI_float ) ;
599 tsar = MRI_FLOAT_PTR(tsim) ;
600 fp = fopen( fn , "w" ) ;
601 for( kim=0 ; kim < imcount ; kim++ ){
602 fprintf(fp , "%7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n" ,
603 roll[kim], pitch[kim], yaw[kim],
604 dx[kim] , dy[kim] , dz[kim] ) ;
605 tsar[0*imcount+kim] = roll[kim] ;
606 tsar[1*imcount+kim] = pitch[kim] ;
607 tsar[2*imcount+kim] = yaw[kim] ;
608 tsar[3*imcount+kim] = dx[kim] ;
609 tsar[4*imcount+kim] = dy[kim] ;
610 tsar[5*imcount+kim] = dz[kim] ;
611 }
612 fclose(fp) ;
613 PLUTO_register_timeseries( fn , tsim ) ;
614 mri_free(tsim) ;
615 }
616
617 /*-- graph --*/
618
619 if( VL_graph && imcount > 1 ){
620 float * yar[7] , dt ;
621 int nn = imcount ;
622 static char * nar[6] = {
623 "\\Delta I-S [mm]" , "\\Delta R-L [mm]" , "\\Delta A-P [mm]" ,
624 "Roll [\\degree]" , "Pitch [\\degree]" , "Yaw [\\degree]" } ;
625
626 yar[0] = (float *) malloc( sizeof(float) * nn ) ;
627 dt = DSET_TIMESTEP(VL_dset) ; if( dt <= 0.0 ) dt = 1.0 ;
628 for( ii=0 ; ii < nn ; ii++ ) yar[0][ii] = ii * dt ;
629
630 yar[1] = dx ; yar[2] = dy ; yar[3] = dz ;
631 yar[4] = roll ; yar[5] = pitch ; yar[6] = yaw ;
632
633 plot_ts_lab( GLOBAL_library.dc->display ,
634 nn , yar[0] , -6 , yar+1 ,
635 "time" , NULL , DSET_FILECODE(VL_dset) , nar , NULL ) ;
636
637 free(yar[0]) ;
638 }
639
640 /* done with these statistics */
641
642 FREEUP(dx) ; FREEUP(dy) ; FREEUP(dz) ;
643 FREEUP(roll) ; FREEUP(yaw) ; FREEUP(pitch) ;
644 FREEUP(rmsnew) ; FREEUP(rmsold) ;
645
646 /*-- save new dataset to disk --*/
647
648 PLUTO_add_dset( plint , new_dset , DSET_ACTION_MAKE_CURRENT ) ;
649 return NULL ;
650 }
651