00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #include "afni.h"
00008 #include "parser.h"
00009 
00010 #if 0
00011 # define VMCHECK do{ if(verbose == 2) MCHECK; } while(0)
00012 #else
00013 # define VMCHECK /* nada */
00014 #endif
00015 
00016 #define TCP_CONTROL "tcp:*:7954"      /* control channel specification */
00017 #define INFO_SIZE  (32*1024)          /* change this ==> change SHM_CHILD below */
00018 #define SHM_CHILD  "shm:afnibahn:32K" /* for data from the child */
00019 
00020 #define SHORT_DELAY      1            /* msec */
00021 #define LONG_DELAY      10
00022 
00023 #ifndef ALLOW_PLUGINS
00024 #  error "Plugins not properly set up -- see machdep.h"
00025 #endif
00026 
00027 #define ALLOW_REGISTRATION
00028 
00029 #ifdef ALLOW_REGISTRATION
00030 #  include "coxplot.h"
00031 #endif
00032 
00033 /***********************************************************************
00034   Plugin to accept data from an external process and assemble
00035   it into a dataset.  Initial implementation for the Bruker 3T/60
00036   system at the Medical College of Wisconsin.
00037 ************************************************************************/
00038 
00039 /** 24 Jun 2002: modified to allow nzz=1 for UCSD trolls                     **/
00040 /** 27 Jun 2003: added BYTEORDER command for automatic byte swapping [rickr] **/
00041 /** 30 Jun 2003: allow MRI_complex data type when using BYTEORDER    [rickr] **/
00042 /** 30 Oct 2003: if possible, compute function on registered data    [rickr] **/
00043 /** 29 Jan 2004: allow 100 chars in root_prefix via PREFIX (from 31) [rickr]
00044                * x-axis of 3-D motion graphs changed from time to reps
00045                * plot_ts_... functions now use reg_rep for x-axis values
00046                * reg_graph_xr is no longer scaled by TR
00047                * added (float *)reg_rep, for graphing with x == rep num
00048                * added RT_set_grapher_pinnums(), to call more than once
00049                * added GRAPH_XRANGE and GRAPH_YRANGE command strings for
00050                      control over the scales of the motion graph
00051                * if GRAPH_XRANGE and GRAPH_YRANGE commands are both passed,
00052                      do not display the final (scaled) motion graph          **/
00053 /** 13 Feb 2004: added RT_MAX_PREFIX for incoming PREFIX command     [rickr]
00054                * if GRAPH_?RANGE is given, disable 'pushing'
00055                      (see plot_ts_xypush())
00056                * added GRAPH_EXPR command, to compute and display a single
00057                      motion curve, instead of the normal six
00058                      (see p_code, etc., reg_eval, and RT_parser_init())      **/
00059 /** 31 Mar 2004: added ability to send registration parameters       [rickr]
00060                * If the AFNI_REALTIME_MP_HOST_PORT environment variable is
00061                      set (as HOST:PORT, e.g. localhost:53214), then the six
00062                      registration correction parameters will be sent to that
00063                      host/port via a tcp socket.  This is done only in the
00064                      case of graphing the 3D registration parameters.
00065                * added RT_input variables to manage the new socket
00066                * added RT_mp_comm_...() functions
00067                * modified yar[] logic in RT_registration_3D_realtime() to
00068                      pass the registration parameters before adjusting the
00069                      base pointers for plot_ts_addto()                       **/
00070 /** 02 Apr 2004: move RT_mp_comm_close() from last plot check  [tross/rickr] **/
00071 /** 10 May 2005: added TPATTERN command to set timing pattern        [rickr] **/
00072 
00073 
00074 /**************************************************************************/
00075 /*********************** struct for reading data **************************/
00076 
00077 #include <sys/types.h>
00078 #include <unistd.h>
00079 #include <sys/wait.h>
00080 #include <signal.h>
00081 #include <ctype.h>
00082 
00083 #define DTYPE_2DZ   77
00084 #define DTYPE_2DZT  78
00085 #define DTYPE_3D    79
00086 #define DTYPE_3DT   80
00087 
00088 #define ZORDER_ALT  33
00089 #define ZORDER_SEQ  34
00090 #define ZORDER_EXP  39
00091 #define NZMAX       1024
00092 
00093 #define RT_MAX_EXPR     1024    /* max size for parser expression  */
00094 #define RT_MAX_PREFIX    100    /* max size for output file prefix */
00095 
00096 #define RT_MP_DEF_PORT 53214    /* default port for sending motion params */
00097 
00098 #define DEFAULT_XYFOV    240.0
00099 #define DEFAULT_XYMATRIX  64
00100 #define DEFAULT_ZDELTA     5.0
00101 #define DEFAULT_ZNUM       2
00102 #define DEFAULT_TR         1.0
00103 
00104 #define RT_NBUF          INFO_SIZE
00105 #define NNAME            128
00106 
00107 #define ORCODE(aa) \
00108   ( (aa)=='R' ? ORI_R2L_TYPE : (aa)=='L' ? ORI_L2R_TYPE : \
00109     (aa)=='P' ? ORI_P2A_TYPE : (aa)=='A' ? ORI_A2P_TYPE : \
00110     (aa)=='I' ? ORI_I2S_TYPE : (aa)=='S' ? ORI_S2I_TYPE : ILLEGAL_TYPE )
00111 
00112 #define MAX_CHAN 32    /* 30 Jul 2002 */
00113 
00114 static int        num_open_controllers ;                       /* 02 Aug 2002 */
00115 static int            open_controller_index[MAX_CONTROLLERS] ;
00116 static Three_D_View * open_controller      [MAX_CONTROLLERS] ;
00117 
00118 typedef struct {
00119 
00120    int info_ok , no_data ;        /* status flags */
00121 
00122    int         image_mode ;       /* 28 Apr 2000 */
00123    void      * image_handle ;
00124    MRI_IMAGE * image_space ;
00125 
00126    char     name_data[NNAME] ;    /* where to get image data */
00127    IOCHAN * ioc_data ;            /* IO channel for image data */
00128 
00129    char     name_info[NNAME] ;    /* where to get control data */
00130    IOCHAN * ioc_info ;            /* IO channel for control data */
00131    pid_t    child_info ;          /* pid of child that will give me control information */
00132    double   child_start_time ;    /* 10 Dec 1998: elapsed time when child process was forked */
00133 
00134    int   nxx  ,nyy  ,nzz  ,       /* matrix */
00135          orcxx,orcyy,orczz ;      /* orientations */
00136 
00137    float xxfov , yyfov , zzfov ;  /* field of view */
00138    float dxx   , dyy   , dzz ,    /* voxel sizes */
00139          xxorg,yyorg,zzorg ;      /* grid origins */
00140    int   xxdcode,yydcode,zzdcode; /* direction code for x,y,z offsets */
00141 
00142    float xxoff , yyoff , zzoff ;  /* offsets (instead of xxorg, etc.) [18 Dec 2002] */
00143 
00144    float zgap ;                   /* extra gap between slices [18 Dec 2002] */
00145                                   /* (only used if zzfov is used instead of dzz) */
00146 
00147    int   xcen ,ycen ,zcen  ;      /* centering of axes? */
00148 
00149    float tr ;                     /* TR */
00150    int   nr ;                     /* expected # volumes */
00151 
00152    int   dtype ;                  /* dataset type: a DTYPE code */
00153    int   zorder ;                 /* 2D slice ordering: a ZORDER code */
00154    int   zorder_lock ;            /* 22 Feb 1999: lock zorder value */
00155    int   tpattern ;               /* 2D slice timing     10 May 2005 [rickr] */
00156    int   nzseq , zseq[NZMAX] ;    /* slice input order */
00157    int   datum ;                  /* a MRI_type code from mrilib.h */
00158    int   swap_on_read ;           /* flag: swap bytes?   26 Jun 2003 [rickr] */
00159 
00160    int   nbuf ;                   /* current buffer size */
00161    char  buf[RT_NBUF] ;           /* buffer for reading command strings */
00162 
00163    char root_prefix[THD_MAX_PREFIX] ;  /* name of dataset (sort of) */
00164 
00165    /* 01 Aug 2002: these items are now indexed by MAX_CHAN
00166                    to allow for multiple channel acquisitions */
00167 
00168    int num_chan ;                               /* number of dataset channels */
00169    int cur_chan ;                               /* current channel index */
00170 
00171    int afni_status[MAX_CHAN] ;                  /* does AFNI know about this dataset yet? */
00172    THD_3dim_dataset * dset[MAX_CHAN] ;          /* AFNI dataset under construction */
00173    Three_D_View * im3d[MAX_CHAN] ;              /* AFNI controller which gets it */
00174    THD_session * sess[MAX_CHAN] ;               /* AFNI session which gets it */
00175    int           sess_num[MAX_CHAN] ;           /* AFNI session index */
00176    char * sbr[MAX_CHAN] ;                       /* sub-brick under construction */
00177    char * im[MAX_CHAN] ;                        /* place for next slice to be read */
00178    int nvol[MAX_CHAN] ;                         /* # volumes read so far */
00179    int nsl[MAX_CHAN] ;                          /* # slices read so far */
00180 
00181    int sbr_size ;                 /* # bytes per sub-brick */
00182    int imsize ;                   /* # bytes per image */
00183    MRI_IMARR * bufar ;            /* buffer for input images */
00184 
00185    THD_3dim_dataset * func_dset ; /* functional dataset, if any */
00186    int func_status ;              /* does AFNI know about this dataset yet? */
00187    int func_code ;                /* how to compute function */
00188    int func_condit ;              /* condition that function computation is in now */
00189    int_func * func_func ;         /* function to compute the function */
00190 
00191 #ifdef ALLOW_REGISTRATION
00192    /*-- 07 Apr and 01 Aug 1998: real-time image registration stuff --*/
00193 
00194    THD_3dim_dataset * reg_dset ;      /* registered dataset, if any */
00195    MRI_2dalign_basis ** reg_2dbasis ; /* stuff for each slice */
00196    int reg_base_index ;               /* where to start? */
00197    int reg_mode ;                     /* how to register? */
00198    int reg_status ;                   /* does AFNI know about this dataset yet? */
00199    int reg_nvol ;                     /* number of volumes registered so far */
00200    int reg_graph ;                    /* 17 Aug 1998: to graph, or not to graph */
00201    int reg_nest ;                     /* number of estimated parameters */
00202    float * reg_tim , * reg_dx  ,
00203          * reg_dy  , * reg_phi  ;     /* estimated motion parameters */
00204 
00205    /*--  Oct 1998: more stuff for 3D registration --*/
00206 
00207    float * reg_dz , * reg_theta , * reg_psi , * reg_rep ;
00208    MRI_3dalign_basis * reg_3dbasis ;
00209    int iha , ax1,hax1 , ax2,hax2 , ax3,hax3 ;
00210    MEM_topshell_data * mp ;
00211 
00212    int   reg_resam , reg_final_resam ;
00213    int   reg_graph_xnew, reg_graph_ynew ;
00214    float reg_graph_xr , reg_graph_yr ;
00215 
00216    /*-- Feb 2004 [rickr]: parser fields for GRAPH_EXPR command --*/
00217    PARSER_code * p_code ;                       /* parser expression code    */
00218    char          p_expr   [RT_MAX_EXPR+1] ;     /* user's parser expression  */
00219    double        p_atoz   [26] ;                /* source values to evaluate */
00220    int           p_has_sym[26] ;                /* symbol indices            */
00221    int           p_max_sym ;                    /* max index+1 of p_has_sym  */
00222    float       * reg_eval ;                     /* EXPR evaluation results   */
00223 
00224    /*-- Mar 2004 [rickr]: tcp comm fields for motion params (for Tom Ross) --*/
00225    int           mp_tcp_use ;     /* are we using tcp comm for motion params */
00226    int           mp_tcp_sd ;      /* socket descriptor                       */
00227    int           mp_port ;        /* destination port for motion params      */
00228    char          mp_host[128] ;   /* destination host for motion params      */
00229    int           mp_nmsg ;        /* count the number of sent messages       */
00230    int           mp_npsets ;      /* count the number of sent data lists     */
00231 #endif
00232 
00233    double elapsed , cpu ;         /* times */
00234    double last_elapsed ;
00235    int    last_nvol ;
00236 
00237    int    num_note ;              /* 01 Oct 2002 */
00238    char **note ;
00239 
00240    int    marked_for_death ;      /* 10 Dec 2002 */
00241 
00242 } RT_input ;
00243 
00244 #define SHOW_TIMES                                              \
00245   fprintf(stderr,"RT: cpu time = %.2f  elapsed time = %.2f\n" , \
00246           PLUTO_cpu_time()-rtin->cpu , PLUTO_elapsed_time()-rtin->elapsed )
00247 
00248 /**************************************************************************/
00249 
00250 static char helpstring[] =
00251    " Purpose: Controlling realtime dataset acquisitions.\n"
00252    "\n"
00253    " USAGE:\n"
00254    " Set the controls to the state you want BEFORE realtime image input\n"
00255    " begins, then press one of the 'Set' buttons to send the control\n"
00256    " information to AFNI.\n"
00257    "\n"
00258    " INPUTS:\n"
00259    " Images Only = If 'No', then the input images will be used\n"
00260    "                 to assemble a dataset.  If 'Yes', then the\n"
00261    "                 input images will be displayed but not made\n"
00262    "                 into a dataset or otherwise saved.\n"
00263    "                 (They will only be stored in the image viewer\n"
00264    "                  that is opened - you can use 'Save:bkg' to\n"
00265    "                  write them to disk, if desired, before\n"
00266    "                  closing the 'Realtime Images' viewer window.)\n"
00267    "               If this input is 'Yes', then none of the other\n"
00268    "                 inputs below have much meaning!\n"
00269    "\n"
00270    " Root     = Root prefix to be used for new datasets.\n"
00271    "              The actual prefix will be of the form Root#001\n"
00272    "\n"
00273    " Update   = How often AFNI's display will be updated\n"
00274    "             = How many 3D bricks are gathered before\n"
00275    "               the images and graphs are redrawn:\n"
00276    "                0 --> don't update until all data is acquired\n"
00277    "                1 --> redraw images and graphs every time\n"
00278    "                2 --> redraw every other brick, etc.\n"
00279    "\n"
00280    " Function = Controls what kind of function analysis is done\n"
00281    "              during realtime input of time dependent datasets.\n"
00282    "\n"
00283    " Verbose  = If set to 'Yes', the plugin will print out progress\n"
00284    "              reports as information flows into it.\n"
00285    "              If set to 'Very', prints out LOTS of information.\n"
00286 #ifdef ALLOW_REGISTRATION
00287    "\n"
00288    " Registration = If activated, image registration will take place\n"
00289    "                  either in realtime or at the end of image acquisition.\n"
00290    "                * The un-registered dataset will be named something like\n"
00291    "                  Root#001 and the aligned dataset Root#001%reg2D.\n"
00292    "                * You can choose either 2D (slice-wise) registration,\n"
00293    "                  or 3D (volume-wise) registration.\n"
00294    "                * The 3D case also allows a quicker mode where the volumes\n"
00295    "                  aren't actually registered, but estimates of the motion\n"
00296    "                  parameters are computed for graphing purposes.  This is\n"
00297    "                  done in realtime.\n"
00298    "                * If 3D realtime or estimate modes are chosen, the motion\n"
00299    "                  parameters can be graphed in realtime - see 'Graph' below.\n"
00300    " Base Image   = The value sets the time index in the new dataset to which\n"
00301    "                  the alignment will take place.\n"
00302    " Resampling   = Determines the interpolation method used:\n"
00303    "                  Cubic     = fastest, least accurate interpolation\n"
00304    "                  Quintic   = intermediate in speed and accuracy\n"
00305    "                  Heptic    = another intermediate interpolation method\n"
00306    "                  Fourier   = slowest, most accurate\n"
00307    "                  Hept+Four = use Heptic for estimation, Fourier for\n"
00308    "                              the final alignment\n"
00309    "              N.B.: Linear interpolation is always used for '3D: estimate'.\n"
00310    "                    Quintic/Heptic interpolation is only available for 3D.\n"
00311    "\n"
00312    " Graph        = If set to 'Yes', will display a graph of the estimated\n"
00313    "                  motion parameters at the end of the run.\n"
00314    "                * If set to 'Realtime', will also show a cruder graph that\n"
00315    "                  alters in realtime (if realtime registration is active).\n"
00316    " NR [x-axis]  = If a realtime graph is generated, this entry specifies\n"
00317    "                  how many repetitions (TR's) to expect.\n"
00318    " YR [y-axis]  = If a realtime graph is generated, this entry specifies\n"
00319    "                  the vertical range for each motion parameter.\n"
00320 #endif
00321    "\n"
00322    "MULTICHANNEL ACQUISITION [Aug 2002]:\n"
00323    " Multiple image channels can be acquired (e.g., from multi-coil and/or\n"
00324    " multi-echo sequences).  Each image channel goes into a separate dataset.\n"
00325    "  * These datasets will have names like Root#007_03 (for the 3rd channel\n"
00326    "    in the 7th acquisition).\n"
00327    "  * Functional activation cannot be computed with multichannel acquisition.\n"
00328    "  * Registration cannot be computed with multichannel acquisition.\n"
00329    "\n"
00330    "HOW TO SEND DATA:\n"
00331    " See file README.realtime for details; see program rtfeedme.c for an example.\n"
00332    "\n"
00333 ;
00334 
00335 /** variables encoding the state of options from the plugin interface **/
00336 
00337 static char root[THD_MAX_PREFIX] = "rt." ;  /* default prefix */
00338 static int  update               = 1  ;     /* default update frequency */
00339 
00340 #define NFUNC  2
00341 static char * FUNC_strings[NFUNC] = { "None" , "FIM" } ;
00342 static int func_code = 0 ;
00343 
00344 static int image_mode = 0 ;  /* 28 Apr 2000 */
00345 
00346 #define FUNC_NONE  0  /* symbolic values for the func_code */
00347 #define FUNC_RTFIM 1
00348 
00349 #define INIT_MODE     -1  /* modes for RT_fim_recurse (AKA func_func) */
00350 #define UPDATE_MODE   -2
00351 #define FINAL_MODE    -3
00352 
00353 int RT_fim_recurse( RT_input * , int ) ;
00354 
00355 int_func * FUNC_funcs[NFUNC] = { NULL , RT_fim_recurse } ;
00356 
00357 #define NYESNO 2
00358 #define NVERB  3
00359 static char * VERB_strings[NVERB] = { "No" , "Yes" , "Very" } ;
00360 static int verbose = 1 ;
00361 
00362 #ifdef ALLOW_REGISTRATION
00363   static int regtime  = 3 ;  /* index of base time for registration */
00364   static int regmode  = 0 ;  /* index into REG_strings */
00365   static int reggraph = 0 ;  /* graphing mode */
00366 
00367   static int reg_resam = 1 ;    /* index into REG_resam_strings */
00368   static int   reg_nr  = 100 ;  /* units of TR */
00369   static float reg_yr  = 1.0 ;  /* mm and degrees */
00370 
00371 #define NGRAPH 3
00372 static char * GRAPH_strings[NGRAPH] = { "No" , "Yes" , "Realtime" } ;
00373 
00374 #define NRESAM 5
00375   static char * REG_resam_strings[NRESAM] = {
00376      "Cubic" , "Quintic" , "Heptic" , "Fourier" , "Hept+Four" } ;
00377   static int REG_resam_ints[NRESAM] = {
00378      MRI_CUBIC , MRI_QUINTIC , MRI_HEPTIC , MRI_FOURIER , -666 } ;
00379 
00380 # define NREG 6
00381   static char * REG_strings[NREG] = {
00382     "None" , "2D: realtime" , "2D: at end"
00383            , "3D: realtime" , "3D: at end" , "3D: estimate" } ;
00384 
00385   static char * REG_strings_ENV[NREG] = {  /* ' ' -> '_'  17 May 2005 [rickr */
00386     "None" , "2D:_realtime" , "2D:_at_end"
00387            , "3D:_realtime" , "3D:_at_end" , "3D:_estimate" } ;
00388 
00389 # define REGMODE_NONE      0
00390 # define REGMODE_2D_RTIME  1
00391 # define REGMODE_2D_ATEND  2
00392 # define REGMODE_3D_RTIME  3
00393 # define REGMODE_3D_ATEND  4
00394 # define REGMODE_3D_ESTIM  5
00395 
00396 # define REG_IS_2D(mm) ( (mm) == REGMODE_2D_RTIME || (mm) == REGMODE_2D_ATEND )
00397 
00398 # define REG_IS_3D(mm) \
00399   ( (mm)==REGMODE_3D_RTIME || (mm)==REGMODE_3D_ATEND || (mm)==REGMODE_3D_ESTIM )
00400 
00401 # define REG_MAKE_DSET(mm)                                   \
00402   ( (mm) == REGMODE_2D_RTIME || (mm) == REGMODE_2D_ATEND ||  \
00403     (mm) == REGMODE_3D_RTIME || (mm) == REGMODE_3D_ATEND   )
00404 #endif
00405 
00406 /************ global data for reading data *****************/
00407 
00408 static IOCHAN * ioc_control = NULL ;     /* control channel */
00409 static RT_input * rtinp = NULL ;         /* only read 1 stream at a time */
00410 static PLUGIN_interface * plint = NULL ; /* AFNI plugin structure */
00411 
00412 /************ prototypes ***********/
00413 
00414 DEFINE_PLUGIN_PROTOTYPE
00415 
00416 char * RT_main( PLUGIN_interface * ) ;
00417 Boolean RT_worker( XtPointer ) ;
00418 RT_input * new_RT_input( IOCHAN * ) ;
00419 int RT_check_listen(void) ;
00420 int RT_acquire_info( char * ) ;
00421 int RT_process_info( int , char * , RT_input * ) ;
00422 void RT_start_dataset( RT_input * ) ;
00423 void RT_read_image( RT_input * , char * ) ;
00424 int RT_process_data( RT_input * ) ;
00425 void RT_process_image( RT_input * ) ;
00426 void RT_finish_dataset( RT_input * ) ;
00427 void RT_tell_afni( RT_input * , int ) ;
00428 void RT_start_child( RT_input * ) ;
00429 void RT_check_info( RT_input * , int ) ;
00430 void RT_process_xevents( RT_input * ) ;  /* 13 Oct 2000 */
00431 
00432 void RT_tell_afni_one( RT_input * , int , int ) ;  /* 01 Aug 2002 */
00433 
00434 #ifdef ALLOW_REGISTRATION
00435   void RT_registration_2D_atend( RT_input * rtin ) ;
00436   void RT_registration_2D_setup( RT_input * rtin ) ;
00437   void RT_registration_2D_close( RT_input * rtin ) ;
00438   void RT_registration_2D_onevol( RT_input * rtin , int tt ) ;
00439   void RT_registration_2D_realtime( RT_input * rtin ) ;
00440 
00441   void RT_registration_3D_atend( RT_input * rtin ) ;
00442   void RT_registration_3D_setup( RT_input * rtin ) ;
00443   void RT_registration_3D_close( RT_input * rtin ) ;
00444   void RT_registration_3D_onevol( RT_input * rtin , int tt ) ;
00445   void RT_registration_3D_realtime( RT_input * rtin ) ;
00446 
00447   int  RT_mp_comm_close       ( RT_input * rtin );
00448   int  RT_mp_comm_init        ( RT_input * rtin );
00449   int  RT_mp_comm_init_vars   ( RT_input * rtin );
00450   int  RT_mp_comm_send_data   ( RT_input * rtin, float * mp[6], int nt );
00451   int  RT_parser_init         ( RT_input * rtin );
00452   void RT_set_grapher_pinnums ( int pinnum );
00453 #endif
00454 
00455 #define TELL_NORMAL  0
00456 #define TELL_WRITE   1
00457 #define TELL_FINAL   2
00458 
00459 #define WRITE_INTERVAL 37.954  /* seconds */
00460 
00461 #define REAP_CHILDREN
00462 #ifdef REAP_CHILDREN
00463 #  define GRIM_REAPER waitpid(-1,NULL,WNOHANG)
00464 #else
00465 #  define GRIM_REAPER /* nada */
00466 #endif
00467 
00468 #define USE_RT_STARTUP
00469 #ifdef  USE_RT_STARTUP
00470 /********************** Register a work process **********************/
00471 
00472 static void RT_startup( XtPointer junk )
00473 {
00474    PLUTO_register_workproc( RT_worker , NULL ) ;
00475    return ;
00476 }
00477 #endif
00478 
00479 /***********************************************************************
00480    Set up the interface to the user
00481    09 Oct 2000: allow all options to be initialized via the environment
00482 ************************************************************************/
00483 
00484 PLUGIN_interface * PLUGIN_init( int ncall )
00485 {
00486    char * ept ; /* 09 Oct 2000 */
00487 
00488    if( ncall > 0 )         return NULL ;  /* only one interface */
00489    if( ! ALLOW_real_time ) return NULL ;  /* do nothing if not allowed */
00490 
00491    /*-- set titles and call point --*/
00492 
00493    plint = PLUTO_new_interface( "RT Options" , "Set Real-Time Acquisition Options" ,
00494                                 helpstring , PLUGIN_CALL_VIA_MENU , RT_main  ) ;
00495 
00496    PLUTO_add_hint( plint , "Set Real-Time Acquisition Options" ) ;
00497 
00498    PLUTO_set_sequence( plint , "A:AArealtime" ) ;
00499    PLUTO_set_butcolor( plint , "hot" ) ;
00500 
00501    PLUTO_set_runlabels( plint , "Set+Keep" , "Set+Close" ) ;  /* 04 Nov 2003 */
00502 
00503    /*-- 28 Apr 2000: Images Only mode --*/
00504 
00505    ept = getenv("AFNI_REALTIME_Images_Only") ;  /* 09 Oct 2000 */
00506    if( ept != NULL ){
00507       int ii = PLUTO_string_index( ept , NYESNO , VERB_strings ) ;
00508       if( ii >= 0 && ii < NYESNO ) image_mode = ii ;
00509    }
00510 
00511    PLUTO_add_option( plint , "" , "Mode" , FALSE ) ;
00512    PLUTO_add_string( plint , "Images Only" , NYESNO,VERB_strings,image_mode ) ;
00513 
00514    /*-- next line of input: Prefix for output dataset --*/
00515 
00516    ept = getenv("AFNI_REALTIME_Root") ;        /* 09 Oct 2000 */
00517    if( !THD_filename_pure(ept) ) ept = NULL ;
00518    if( ept != NULL ) MCW_strncpy(root,ept,THD_MAX_PREFIX) ;
00519 
00520    PLUTO_add_option( plint , "" , "Root" , FALSE ) ;
00521    PLUTO_add_string( plint , "Root" , 0, (ept!=NULL) ? &ept : NULL , 19 ) ;
00522 
00523    /*-- next line of input: Update frequency --*/
00524 
00525    ept = getenv("AFNI_REALTIME_Update") ;      /* 09 Oct 2000 */
00526    if( ept != NULL ){
00527       int ii = (int) rint(strtod(ept,NULL)) ;
00528       if( ii >= 0 && ii <= 19 ) update = ii ;
00529    }
00530 
00531    PLUTO_add_option( plint , "" , "Update" , FALSE ) ;
00532    PLUTO_add_number( plint , "Update" , 0,19,0 , update , FALSE ) ;
00533 
00534    /*-- next line of input: Function computation --*/
00535 
00536    ept = getenv("AFNI_REALTIME_Function") ;   /* 09 Oct 2000 */
00537    if( ept != NULL ){
00538       int ii = PLUTO_string_index( ept , NFUNC , FUNC_strings ) ;
00539       if( ii >= 0 && ii < NFUNC ) func_code = ii ;
00540    }
00541 
00542    PLUTO_add_option( plint , "" , "Function" , FALSE ) ;
00543    PLUTO_add_string( plint , "Function" , NFUNC , FUNC_strings , func_code ) ;
00544 
00545    /*-- next line of input: Verbosity flag --*/
00546 
00547    ept = getenv("AFNI_REALTIME_Verbose") ;   /* 09 Oct 2000 */
00548    if( ept != NULL ){
00549       int ii = PLUTO_string_index( ept , NVERB , VERB_strings ) ;
00550       if( ii >= 0 && ii < NVERB ) verbose = ii ;
00551    }
00552 
00553    PLUTO_add_option( plint , "" , "Verbose" , FALSE ) ;
00554    PLUTO_add_string( plint , "Verbose" , NVERB , VERB_strings , verbose ) ;
00555 
00556 #ifdef ALLOW_REGISTRATION
00557    /*-- next line of input: registration mode --*/
00558 
00559    ept = getenv("AFNI_REALTIME_Registration") ;  /* 09 Oct 2000 */
00560    if( ept != NULL ){
00561       int ii = PLUTO_string_index( ept , NREG , REG_strings ) ;
00562       /* on failure, try alternate forms               17 May 2005 [rickr] */
00563       if ( ii < 0 ) ii = PLUTO_string_index( ept , NREG , REG_strings_ENV ) ;
00564       if( ii >= 0 && ii < NREG ) regmode = ii ;
00565    }
00566 
00567    ept = getenv("AFNI_REALTIME_Base_Image") ;     /* 09 Oct 2000 */
00568    if( ept != NULL ){
00569       int ii = (int) rint(strtod(ept,NULL)) ;
00570       if( ii >= 0 && ii <= 59 ) regtime = ii ;
00571    }
00572 
00573    ept = getenv("AFNI_REALTIME_Resampling") ;    /* 09 Oct 2000 */
00574    if( ept != NULL ){
00575       int ii = PLUTO_string_index( ept , NRESAM , REG_resam_strings ) ;
00576       if( ii >= 0 && ii < NRESAM ) reg_resam = ii ;
00577    }
00578 
00579    PLUTO_add_option( plint , "" , "Registration" , FALSE ) ;
00580    PLUTO_add_string( plint , "Registration" , NREG , REG_strings , regmode ) ;
00581    PLUTO_add_number( plint , "Base Image" , 0,59,0 , regtime , FALSE ) ;
00582    PLUTO_add_string( plint , "Resampling" , NRESAM , REG_resam_strings , reg_resam ) ;
00583 
00584    /*-- next line of input: registration graphing --*/
00585 
00586    ept = getenv("AFNI_REALTIME_Graph")  ;  /* 09 Oct 2000 */
00587    if( ept != NULL ){
00588       int ii = PLUTO_string_index( ept , NGRAPH , GRAPH_strings ) ;
00589       if( ii >= 0 && ii < NGRAPH ) reggraph = ii ;
00590    }
00591 
00592    ept = getenv("AFNI_REALTIME_NR") ;     /* 09 Oct 2000 */
00593    if( ept != NULL ){
00594       int ii = (int) rint(strtod(ept,NULL)) ;
00595       if( ii >= 5 && ii <= 9999 ) reg_nr = ii ;
00596    }
00597 
00598    ept = getenv("AFNI_REALTIME_YR") ;     /* 09 Oct 2000 */
00599    if( ept != NULL ){
00600       float ff = strtod(ept,NULL) ;
00601       if( ff > 0.0 ) reg_yr = ff ;
00602    }
00603 
00604    PLUTO_add_option( plint , "" , "Graphing" , FALSE ) ;
00605    PLUTO_add_string( plint , "Graph" , NGRAPH , GRAPH_strings , reggraph ) ;
00606    PLUTO_add_number( plint , "NR [x-axis]" , 5,9999,0 , reg_nr , TRUE ) ;
00607    PLUTO_add_number( plint , "YR [y-axis]" , 1,100,1 , (int)(reg_yr*10.0) , TRUE ) ;
00608 #endif
00609 
00610    /***** Register a work process *****/
00611 
00612 #ifndef USE_RT_STARTUP
00613    PLUTO_register_workproc( RT_worker , NULL ) ;
00614 #else
00615    PLUTO_register_timeout( 1954 , RT_startup , NULL ) ;
00616 #endif
00617 
00618    /***** 12 Oct 2000: set the graphing geometry, if any *****/
00619 
00620    ept = getenv("AFNI_REALTIME_volreg_graphgeom") ;
00621    if( ept != NULL ){
00622       char *str = malloc(strlen(ept)+20) ;
00623       sprintf(str,"AFNI_tsplotgeom=%s",ept) ;
00624       putenv(str) ;
00625    }
00626 
00627    /***** go home to mama (i.e., AFNI) *****/
00628 
00629    ALLOW_real_time = 1 ;  /* flag to AFNI that realtime work is started */
00630    return plint ;
00631 }
00632 
00633 /***************************************************************************
00634   Main routine for this plugin (will be called from AFNI).
00635 ****************************************************************************/
00636 
00637 char * RT_main( PLUGIN_interface * plint )
00638 {
00639    char * tag , * str , * new_prefix ;
00640    static char buf[256] ;
00641 
00642    if( plint == NULL )
00643       return "*********************\n"
00644              "RT_main:  NULL input\n"
00645              "*********************"  ;
00646 
00647    /** loop over input from AFNI **/
00648 
00649 #if 0                         /* 09 Oct 2000: turn this "feature" off */
00650 #ifdef ALLOW_REGISTRATION
00651    regmode = REGMODE_NONE ;   /* no registration if not ordered explicitly */
00652 #endif
00653 #endif
00654 
00655    while( (tag=PLUTO_get_optiontag(plint)) != NULL ){
00656 
00657       /* 28 Apr 2000: added "Images Only" mode */
00658 
00659       if( strcmp(tag,"Mode") == 0 ){
00660          str        = PLUTO_get_string(plint) ;
00661          image_mode = PLUTO_string_index( str , NYESNO , VERB_strings ) ;
00662          continue ;
00663       }
00664 
00665       if( strcmp(tag,"Root") == 0 ){
00666          new_prefix = PLUTO_get_string(plint) ;
00667          if( ! THD_filename_pure(new_prefix) )
00668             return "**************************\n"
00669                    "RT_main:  bad root prefix\n"
00670                    "**************************"  ;
00671          strcpy(root,new_prefix) ;
00672          continue ;
00673       }
00674 
00675       if( strcmp(tag,"Update") == 0 ){
00676          update = PLUTO_get_number(plint) ;
00677          continue ;
00678       }
00679 
00680       if( strcmp(tag,"Function") == 0 ){
00681          str       = PLUTO_get_string(plint) ;
00682          func_code = PLUTO_string_index( str , NFUNC , FUNC_strings ) ;
00683          continue ;
00684       }
00685 
00686       if( strcmp(tag,"Verbose") == 0 ){
00687          str     = PLUTO_get_string(plint) ;
00688          verbose = PLUTO_string_index( str , NVERB , VERB_strings ) ;
00689          continue ;
00690       }
00691 
00692 #ifdef ALLOW_REGISTRATION
00693       if( strcmp(tag,"Registration") == 0 ){
00694          str       = PLUTO_get_string(plint) ;
00695          regmode   = PLUTO_string_index( str , NREG , REG_strings ) ;
00696 
00697          regtime   = PLUTO_get_number(plint) ;
00698 
00699          str       = PLUTO_get_string(plint) ;
00700          reg_resam = PLUTO_string_index( str , NRESAM , REG_resam_strings ) ;
00701          continue ;
00702       }
00703 
00704       if( strcmp(tag,"Graphing") == 0 ){
00705          str      = PLUTO_get_string(plint) ;
00706          reggraph = PLUTO_string_index( str , NGRAPH , GRAPH_strings ) ;
00707 
00708          reg_nr   = PLUTO_get_number(plint) ;
00709          reg_yr   = PLUTO_get_number(plint) ;
00710 
00711          /* 12 Oct 2000: set pin_num on all graphs now open */
00712 
00713          if( reg_nr >= MIN_PIN && reg_nr <= MAX_PIN && IM3D_OPEN(plint->im3d) )
00714              RT_set_grapher_pinnums(reg_nr);
00715 
00716          continue ;
00717       }
00718 #endif
00719 
00720       /** How the hell did this happen? **/
00721 
00722       sprintf(buf,"*****************\n"
00723                   "Illegal optiontag: %s\n"
00724                   "*****************"      , tag ) ;
00725       return buf ;
00726 
00727    }  /* end of loop over active input options */
00728 
00729    /*-- 28 Apr 2000: if in Image Only mode, turn some stuff off --*/
00730 
00731    if( image_mode ){
00732       func_code = 0 ;
00733       regmode   = 0 ;
00734       reggraph  = 0 ;
00735    }
00736 
00737    PLUTO_turnoff_options( plint ) ;  /* 21 Feb 2001 */
00738 
00739    return NULL ;  /* nothing bad happened */
00740 }
00741 
00742 /***************************************************************************
00743    Listen for an incoming real-time data control connection.
00744    Return 1 if the connection is established,
00745    return 0 if no connection is ready, return -1 if an error occurs.
00746 ****************************************************************************/
00747 
00748 int RT_check_listen(void)
00749 {
00750    int jj ;
00751    static int newcon = 1 ;
00752 
00753    /** see if we need to start listening **/
00754 
00755    if( ioc_control == NULL ){
00756       if( verbose )
00757          fprintf(stderr,"RT: starting to listen for control stream.\n") ;
00758       ioc_control = iochan_init( TCP_CONTROL , "accept" ) ;
00759       newcon      = 1 ;
00760       if( ioc_control == NULL ){
00761          fprintf(stderr,"RT: can't listen for control stream\a\n") ;
00762          return -1 ;
00763       }
00764    }
00765 
00766    /** Check if the network channel is active **/
00767 
00768    jj = iochan_goodcheck(ioc_control,SHORT_DELAY) ;
00769 
00770    if( jj == 1 ){  /* external connection! */
00771 
00772       if( newcon ){
00773          fprintf(stderr,"RT:---------------------------------------\n") ;
00774          fprintf(stderr,"RT: connected to control stream %s\n",ioc_control->name) ;
00775          newcon = 0 ;
00776       } else {
00777 /**
00778          if( verbose == 2 )
00779             fprintf(stderr,"RT: still waiting for control data.\n") ;
00780 **/
00781       }
00782 
00783       if( ! TRUST_host(ioc_control->name) ){
00784          fprintf(stderr,"RT: untrusted host connection - closing!\a\n") ;
00785          IOCHAN_CLOSENOW(ioc_control) ;
00786          return 0 ;
00787       }
00788 
00789       jj = iochan_readcheck(ioc_control,0) ;  /* is something ready to read? */
00790 
00791       if( jj > 0 && verbose == 2 ) fprintf(stderr,"RT: control data is present!\n") ;
00792 
00793       return jj ;
00794 
00795    } else if( jj == -1 ){  /* something bad! */
00796 
00797       fprintf(stderr,"RT: failure while listening for control stream!\a\n") ;
00798       IOCHAN_CLOSENOW(ioc_control) ;
00799       return 0 ;
00800    }
00801 
00802    return 0 ;  /* no external connection yet */
00803 }
00804 
00805 /***************************************************************************/
00806 /**             macro to close down the realtime input stream             **/
00807 
00808 #define CLEANUP(n) cleanup_rtinp(n)
00809 
00810 /**------------ 01 Aug 1998: replace the macro with a function -----------**/
00811 
00812 #undef  FREEUP
00813 #define FREEUP(x) do{ if( (x) != NULL ){free((x)); (x)=NULL;} } while(0)
00814 
00815 /* 10 Dec 2002: add keep_ioc_data flag,
00816                 so as not to close the data channel;
00817                 note that the pointer to it will be lost when rtinp is freed */
00818 
00819 void cleanup_rtinp( int keep_ioc_data )
00820 {
00821    int cc ;
00822 
00823    if( !keep_ioc_data )
00824      IOCHAN_CLOSENOW(rtinp->ioc_data) ;    /* close open I/O channels */
00825 
00826    IOCHAN_CLOSENOW(rtinp->ioc_info) ;
00827 
00828    if( rtinp->child_info > 0 )             /* destroy child process */
00829       kill( rtinp->child_info , SIGTERM ) ;
00830 
00831    DESTROY_IMARR(rtinp->bufar) ;           /* destroy any buffered images */
00832 
00833    for( cc=0 ; cc < MAX_CHAN ; cc++ ){     /* 01 Aug 2002: loop over channels */
00834      if( rtinp->sbr[cc] != NULL )
00835        free( rtinp->sbr[cc] ) ;            /* destroy buffered data */
00836    }
00837 
00838 #ifdef ALLOW_REGISTRATION
00839    if( rtinp->reg_2dbasis != NULL ){       /* destroy registration setup */
00840       int kk ;
00841       for( kk=0 ; kk < rtinp->nzz ; kk++ )
00842          mri_2dalign_cleanup( rtinp->reg_2dbasis[kk] ) ;
00843       free( rtinp->reg_2dbasis ) ;
00844    }
00845 
00846    if( rtinp->reg_3dbasis != NULL ){
00847       mri_3dalign_cleanup( rtinp->reg_3dbasis ) ;
00848    }
00849 
00850    FREEUP( rtinp->reg_tim   ) ; FREEUP( rtinp->reg_dx    ) ;
00851    FREEUP( rtinp->reg_dy    ) ; FREEUP( rtinp->reg_dz    ) ;
00852    FREEUP( rtinp->reg_phi   ) ; FREEUP( rtinp->reg_psi   ) ;
00853    FREEUP( rtinp->reg_theta ) ; FREEUP( rtinp->reg_rep   ) ;
00854    FREEUP( rtinp->reg_eval  ) ;
00855 #endif
00856 
00857    if( rtinp->image_handle != NULL )
00858       PLUTO_imseq_rekill( rtinp->image_handle,NULL,NULL ) ;  /* 28 Apr 2000 */
00859 
00860    if( rtinp->image_space != NULL ){
00861       mri_clear_data_pointer(rtinp->image_space) ; mri_free(rtinp->image_space) ;
00862    }
00863 
00864    /* 01 Oct 2002: free stored notes */
00865 
00866    if( rtinp->num_note > 0 && rtinp->note != NULL ){
00867      int kk ;
00868      for( kk=0 ; kk < rtinp->num_note ; kk++ ) FREEUP( rtinp->note[kk] ) ;
00869      FREEUP(rtinp->note) ;
00870    }
00871 
00872    free(rtinp) ; rtinp = NULL ;            /* destroy data structure */
00873    ioc_control = NULL ;                    /* ready to listen again */
00874    GRIM_REAPER ;                           /* reap dead child, if any */
00875 }
00876 
00877 /*********************************************************************
00878   Real time routine called every so often by Xt.
00879   This handles connections from other processes that want to
00880   send AFNI data.
00881 
00882   This is a work process -- it returns False if it wants to be
00883   called again, and return True if it is done.  For real-time
00884   AFNI-izing, this should always return False.
00885 **********************************************************************/
00886 
00887 Boolean RT_worker( XtPointer elvis )
00888 {
00889    int jj ;
00890    static int first=1 ;
00891 
00892 #if 0
00893    if( first ){
00894       if( verbose ) fprintf(stderr,"RT: first call to RT_worker()\n") ;
00895       first = 0 ;
00896    }
00897 #endif
00898 
00899    /**--------------------------------------------------------------------**/
00900    /** if we are waiting for a connection, check to see if it is good yet **/
00901 
00902    if( rtinp == NULL ){         /* not doing real-time input at this time */
00903       static int nerr = 0 ;
00904       jj = RT_check_listen() ;             /* see if someone wants to talk */
00905       if( jj < 0 ){                        /* quit if there's an error */
00906          nerr++ ; iochan_sleep(1) ;
00907          return (nerr > 9) ? True : False; /* 12 Oct 2000: if a lot of errors */
00908       }
00909       nerr = 0 ;                           /* reset error count */
00910       if( jj == 0 ) return False ;         /* try later if no connection now */
00911       rtinp = new_RT_input(NULL) ;         /* try to make a new input struct */
00912       IOCHAN_CLOSENOW( ioc_control ) ;     /* not needed any more */
00913       if( rtinp == NULL ) return False ;   /* try later (?) */
00914    }
00915 
00916    /**---------------------------------------------------------------**/
00917    /**---- at this point, the rtinp struct is good,              ----**/
00918    /**---- meaning that we are in the business of acquiring data ----**/
00919    /**---------------------------------------------------------------**/
00920 
00921    /**-----------------------------------------------------------------------**/
00922    /** if input data stream has gone bad, close down the real-time operation **/
00923 
00924    if( iochan_goodcheck(rtinp->ioc_data,0) != 1 ){
00925 
00926       if( rtinp->sbr[0] != NULL ){     /* if we started acquisition */
00927          if( verbose == 2 )
00928             fprintf(stderr,"RT: data stream closed down.\n") ;
00929          VMCHECK ;
00930          RT_finish_dataset( rtinp ) ;  /* then we can finish it */
00931       } else {
00932          fprintf(stderr,"RT: data stream closed before dataset was fully defined!\a\n") ;
00933       }
00934 
00935       CLEANUP(0) ; return False ;
00936    }
00937 
00938    /**-----------------------------------------------------------------**/
00939    /**---- See if any data is waiting to be read from the child   ----**/
00940    /**---- process that is supposed to supply control information ----**/
00941 
00942 #define CHILD_MAX_WAIT 66.6  /* seconds */
00943 
00944    if( rtinp->child_info > 0 ){
00945 
00946       jj = iochan_readcheck( rtinp->ioc_info , SHORT_DELAY ) ;  /** data ready? **/
00947 
00948       if( jj == 0 ){       /** 10 Dec 1998: child not sending data? **/
00949          double et = PLUTO_elapsed_time() - rtinp->child_start_time ;
00950          double mw = CHILD_MAX_WAIT ;
00951          char *eee = getenv("AFNI_REALTIME_CHILDWAIT") ;
00952          if( eee != NULL ){
00953             double val=strtod(eee,NULL); if(val >= 1.0) mw = val;
00954          }
00955 
00956          if( et > mw ){  /* don't wait any more, give up */
00957             fprintf(stderr,
00958                     "RT: no data from child after %f seconds!  Giving up.\a\n",et) ;
00959             CLEANUP(0) ; return False ;
00960          }
00961 
00962       } else if( jj < 0 ){   /** something bad happened? **/
00963 
00964          fprintf(stderr,"RT: child info stream closed prematurely!\a\n") ;
00965          CLEANUP(0) ; return False ;
00966 
00967       } else if( jj > 0 ){   /** something good happened? **/
00968 
00969          char * info = (char *) malloc( sizeof(char) * INFO_SIZE ) ;
00970          int   ninfo , ii ;
00971 
00972          /* read all data possible from the info channel */
00973 
00974          if( verbose )
00975             fprintf(stderr,"RT: receiving data from child process\n") ;
00976          VMCHECK ;
00977 
00978          ninfo = 0 ;
00979          while(1){
00980             ii = iochan_recv( rtinp->ioc_info , info+ninfo , INFO_SIZE-ninfo ) ;
00981             if( ii < 1 ) break ;
00982             ninfo += ii ;
00983             if( ninfo >= INFO_SIZE ){
00984                fprintf(stderr,"RT: child process sent too much info!\a\n") ;
00985                break ;
00986             }
00987             ii = iochan_readcheck( rtinp->ioc_info , SHORT_DELAY ) ;
00988             if( ii < 1 ) break ;
00989          }
00990          IOCHAN_CLOSENOW( rtinp->ioc_info ) ; rtinp->child_info = 0 ;
00991 
00992          if( ninfo <= 0 ){
00993             fprintf(stderr,"RT: child info stream returned no data!\a\n") ;
00994             CLEANUP(0) ; return False ;
00995          }
00996 
00997          if( verbose == 2 )
00998             fprintf(stderr,"RT: child info stream returned %d bytes.\n",ninfo) ;
00999          VMCHECK ;
01000 
01001          /* process the info and store it in the real-time struct */
01002 
01003          jj = RT_process_info( ninfo , info , rtinp ) ; free(info) ;
01004 
01005          /* check if data was processed OK */
01006 
01007          if( jj <= 0 ){
01008             fprintf(stderr,"RT: child info was badly formatted!\a\n") ;
01009             CLEANUP(0) ; return False ;
01010          }
01011 
01012          /* if we already received the first image data,
01013             then all the setup info should be OK.  if not, quit */
01014 
01015          if( ! rtinp->no_data && ! rtinp->info_ok ){
01016             fprintf(stderr,"RT: child info was incomplete or erroneous!\a\n") ;
01017             RT_check_info( rtinp , 1 ) ;
01018             PLUTO_beep() ;
01019             PLUTO_popup_transient( plint , " \n"
01020                                            "      Heads down!\n"
01021                                            "Realtime header was bad!\n" ) ;
01022             CLEANUP(0) ; return FALSE ;
01023          }
01024 
01025          /* if all the setup info is OK, we can create the dataset now */
01026 
01027          if( rtinp->sbr[0] == NULL && rtinp->info_ok ){
01028             if( verbose == 2 )
01029                fprintf(stderr,"RT: info complete --> creating dataset.\n") ;
01030             VMCHECK ;
01031             RT_start_dataset( rtinp ) ;
01032          }
01033       }
01034 
01035       /** if nothing happened (good or bad),
01036           will try to read it again on the next time thru **/
01037    }
01038 
01039    /**---------------------------------------------**/
01040    /** See if any image data is waiting to be read **/
01041 
01042    jj = iochan_readcheck( rtinp->ioc_data , SHORT_DELAY ) ;
01043 
01044    if( jj == 0 ){     /** no data available now **/
01045 
01046 #ifdef WRITE_INTERVAL
01047       double delt ;
01048       double mw = WRITE_INTERVAL ;
01049       char *eee = getenv("AFNI_REALTIME_WRITEWAIT") ;
01050       if( eee != NULL ){
01051          double val=strtod(eee,NULL); if(val >= 1.0) mw = val;
01052       }
01053 
01054       /** if there is no image data for a long time,
01055           and there is something new to write out,
01056           then write the datasets out again.         **/
01057 
01058       if( !rtinp->no_data                    &&
01059           rtinp->nvol[0] > rtinp->last_nvol  &&
01060           !rtinp->image_mode                 &&  /* 28 Apr 2000 */
01061           (delt=PLUTO_elapsed_time()-rtinp->last_elapsed) > mw  ){
01062 
01063          int cmode , cc ;
01064 
01065          fprintf(stderr,"RT: no image data for %g seconds --> saving to disk.\n",delt) ;
01066 
01067          PLUTO_popup_transient( plint , " \n"
01068                                         " Pause in input data stream:\n"
01069                                         " Saving current dataset(s) to disk.\n" ) ;
01070 
01071          RT_tell_afni(rtinp,TELL_NORMAL) ;
01072 
01073          cmode = THD_get_write_compression() ;      /* 20 Mar 1998 */
01074          THD_set_write_compression(COMPRESS_NONE) ;
01075          SHOW_AFNI_PAUSE ;
01076 
01077          for( cc=0 ; cc < rtinp->num_chan ; cc++ )
01078            THD_write_3dim_dataset( NULL,NULL , rtinp->dset[cc] , True ) ;
01079 
01080          if( rtinp->func_dset != NULL )
01081             THD_write_3dim_dataset( NULL,NULL , rtinp->func_dset , True ) ;
01082 
01083          THD_set_write_compression(cmode) ;  sync() ; /* 08 Mar 2000: sync disk */
01084          SHOW_AFNI_READY ;
01085 
01086          rtinp->last_nvol    = rtinp->nvol[0] ;
01087          rtinp->last_elapsed = PLUTO_elapsed_time() ;
01088       }
01089 #endif
01090 
01091       return False ;
01092    }
01093 
01094    if( jj < 0 ){               /** something bad happened to data channel **/
01095       if( verbose == 2 )
01096          fprintf(stderr,"RT: data stream closed down.\n") ;
01097       VMCHECK ;
01098       if( rtinp->sbr[0] != NULL ) RT_finish_dataset( rtinp ) ;
01099       CLEANUP(0) ; return False ;
01100    }
01101 
01102    /************************************************/
01103    /** If here, data is ready on the data channel **/
01104 
01105    /**-----------------------------------------------------**/
01106    /** If this is the first time reading the data channel, **/
01107    /** read and process any header info present there.     **/
01108    /** (Will read the image data a little farther below.)  **/
01109 
01110    if( rtinp->no_data ){
01111       int ii , nb ;
01112 
01113       /** if so ordered, start a child process to get header info **/
01114 
01115       if( strlen(rtinp->name_info) > 0 ) RT_start_child( rtinp ) ;
01116 
01117       /** read initial stuff from the data channel **/
01118 
01119       fprintf(stderr,"RT: receiving image metadata") ; fflush(stderr) ;
01120 
01121       nb = iochan_recv( rtinp->ioc_data , rtinp->buf , RT_NBUF ) ;
01122 
01123       if( nb <= 0 ){
01124          fprintf(stderr,"\nRT: recv on data stream fails on first try!\a\n") ;
01125          CLEANUP(0) ; return False ;
01126       }
01127 
01128       /* read any more that comes available very quickly [06 Aug 2002] */
01129 
01130       while( nb < RT_NBUF-1 ){
01131         ii = iochan_readcheck( rtinp->ioc_data , SHORT_DELAY ) ;
01132         if( ii <= 0 ) break ;
01133         ii = iochan_recv( rtinp->ioc_data , rtinp->buf+nb , RT_NBUF-nb ) ;
01134         if( ii <= 0 ) break ;
01135         nb += ii ;
01136       }
01137       rtinp->nbuf = nb ;
01138 
01139       fprintf(stderr,"=%d bytes\n",rtinp->nbuf) ;
01140 
01141       PLUTO_beep() ;
01142       PLUTO_popup_transient( plint , " \n"
01143                                      "***************************\n"
01144                                      "*       Heads Up!         *\n"
01145                                      "* Incoming realtime data! *\n"
01146                                      "***************************\n" ) ;
01147 
01148       /** process header info in the data channel;
01149           note that there must be SOME header info this first time **/
01150 
01151       jj = RT_process_info( rtinp->nbuf , rtinp->buf , rtinp ) ;
01152 
01153       if( jj <= 0 ){                  /* header info not OK */
01154 
01155          fprintf(stderr,"RT: initial image metadata is badly formatted!\a\n") ;
01156          CLEANUP(0) ; return False ;
01157 
01158       } else if( jj < rtinp->nbuf ){  /* data bytes left after header info */
01159 
01160          memmove( rtinp->buf , rtinp->buf + jj , rtinp->nbuf - jj ) ;
01161          rtinp->nbuf = rtinp->nbuf - jj ;
01162 
01163       } else {                        /* no data left after header info */
01164          rtinp->nbuf = 0 ;
01165       }
01166 
01167       if( verbose == 2 )
01168          fprintf(stderr,
01169                  "RT: processed %d bytes of header info from data stream\n",jj) ;
01170       VMCHECK ;
01171 
01172       rtinp->no_data = 0 ;  /* can't say we never received data */
01173 
01174       /* if we have already received the child process info,
01175          or if there is no child process info, then all should be OK now */
01176 
01177       if( rtinp->child_info == 0 && ! rtinp->info_ok ){
01178          fprintf(stderr,"RT: image header info was incomplete!\a\n") ;
01179          RT_check_info( rtinp , 1 ) ;
01180          PLUTO_beep() ;
01181          PLUTO_popup_transient( plint , " \n"
01182                                         "      Heads down!\n"
01183                                         "Realtime header was bad!\n" ) ;
01184          CLEANUP(0) ; return FALSE ;
01185       }
01186    }
01187 
01188    /**-------------------------------**/
01189    /** At last -- process image data **/
01190 
01191    jj = RT_process_data( rtinp ) ;
01192 
01193    if( jj < 0 ){
01194       fprintf(stderr,"RT: data stream aborted during image read!\n") ;
01195       if( rtinp->sbr[0] != NULL ) RT_finish_dataset( rtinp ) ;
01196       CLEANUP(0) ; return False ;
01197    }
01198 
01199    /* 10 Dec 2002: if data stream marked itself for death,
01200                    then close down the dataset,
01201                    but don't close down the data stream itself;
01202                    instead, create a new RT_input context for
01203                    reading new metadata+image data from the
01204                    same data stream                          */
01205 
01206    if( rtinp->marked_for_death ){
01207       RT_input *new_rtinp ;
01208       fprintf(stderr,"RT: data stream says to close dataset.\n") ;
01209       if( rtinp->sbr[0] != NULL ) RT_finish_dataset( rtinp ) ;
01210       fprintf(stderr,"RT: starting to read from existing data stream.\n") ;
01211       new_rtinp = new_RT_input( rtinp->ioc_data ) ; /* new RT_input context */
01212       CLEANUP(1) ;                                 /* will delete old rtinp */
01213       rtinp = new_rtinp ; return False ;
01214    }
01215 
01216    /*** continue normally! ***/
01217 
01218    rtinp->last_elapsed = PLUTO_elapsed_time() ;  /* record this time */
01219    return False ;
01220 }
01221 
01222 /*---------------------------------------------------------------------
01223    13 Oct 2000: process some X events, to make AFNI less sluggish
01224 -----------------------------------------------------------------------*/
01225 
01226 #define MAX_NEV 6
01227 #if MAX_NEV > 0
01228 void RT_process_xevents( RT_input * rtin )
01229 {
01230    Display * dis = THE_DISPLAY ;
01231    XEvent ev ; int nev=0 ;
01232 
01233    XSync( dis , False ) ;
01234    while( nev++ < MAX_NEV &&
01235           XCheckMaskEvent(dis ,
01236                           ButtonMotionMask   |PointerMotionMask|
01237                           ButtonPressMask    |ButtonReleaseMask|
01238                           KeyPressMask       |KeyReleaseMask   |
01239                           StructureNotifyMask|ExposureMask      ,&ev) ){
01240 
01241           XtDispatchEvent( &ev ) ;  /* do the actual work for this event */
01242    }
01243    XmUpdateDisplay(THE_TOPSHELL) ;
01244    if( verbose == 2 && nev > 1 )
01245       fprintf(stderr,"RT: processed %d events\n",nev-1);
01246    return ;
01247 }
01248 #else
01249 void RT_process_xevents( RT_input * rtin ){}  /* doesn't do much */
01250 #endif
01251 
01252 /*--------------------------------------------------------------------------
01253    Initialize a new RT_input structure; returns pointer to new structure.
01254    Will read from the control channel, which will tell it what data
01255    channel to open, and if it needs to fork a child process to gather
01256    header info.
01257 
01258    10 Dec 2002: new input ioc_data is a pointer to an IOCHAN that
01259                 should be used for data acquisition:
01260                   - if this is not NULL, then the control channel is NOT used
01261                   - instead, the RT_input struct is setup to read
01262                      directly from ioc_data (which should already be good)
01263                   - in this case, you can't use an info_command; all the
01264                      metadata must come from the ioc_data channel
01265 ----------------------------------------------------------------------------*/
01266 
01267 RT_input * new_RT_input( IOCHAN *ioc_data )
01268 {
01269    RT_input *rtin ;
01270    int ii , cc ;
01271    Three_D_View *im3d ;
01272 
01273    if( ioc_data == NULL ){ /*** THE OLD WAY: get info from control channel ***/
01274 
01275      int ncon ;
01276      char *con , *ptr ;
01277 
01278      /** wait until data can be read, or something terrible happens **/
01279 
01280      if( iochan_readcheck(ioc_control,-1) <= 0 ){
01281         fprintf(stderr,"RT: control stream fails readcheck!\a\n") ;
01282         return NULL ;
01283      }
01284 
01285      /** make new structure **/
01286 
01287      rtin = (RT_input *) calloc( 1 , sizeof(RT_input) ) ;
01288      con  = (char *)     malloc( INFO_SIZE ) ;
01289 
01290      if( rtin == NULL || con == NULL ){
01291        fprintf(stderr,"RT: malloc fails in new_RT_input!\a\n") ; EXIT(1) ;
01292      }
01293 
01294      /** read all data possible from control channel **/
01295 
01296      ncon = 0 ;  /* read all data possible from the control channel */
01297      while(1){
01298         ii = iochan_recv( ioc_control , con+ncon , INFO_SIZE-ncon ) ;
01299         if( ii < 1 ) break ;
01300         ncon += ii ;
01301         if( ncon >= INFO_SIZE ){
01302            fprintf(stderr,"RT: control stream buffer overflow!\a\n") ;
01303            break ;
01304         }
01305         iochan_sleep( SHORT_DELAY ) ;
01306      }
01307 
01308      if( ncon < 1 ){
01309         fprintf(stderr,"RT: control stream sends no data!\a\n") ;
01310         free(rtin) ; free(con) ; return NULL ;
01311      }
01312 
01313      /** The data we now have in 'con' should be of the form
01314             image_channel_spec\n
01315             info_command\n\0
01316          where 'image_channel_spec' is an IOCHAN specifier saying
01317            how AFNI should read the image data
01318          and 'info_command' is the name of a command to run that
01319            will write the dataset control info to stdout (where
01320            a child of AFNI will get it using the 'popen' routine).
01321            If no 'info_command' is needed (all the dataset info will
01322            come thru on the image channel), then the info_command
01323            should just be skipped.
01324          Each string should be less than NNAME characters long!
01325      **/
01326 
01327      ncon = MIN( ncon , INFO_SIZE-2 ) ;
01328      con[ncon+1] = '\0' ;               /* make sure it is NUL terminated */
01329 
01330      /** copy first string -- this is image_channel_spec **/
01331 
01332      ii = 0 ;
01333      for( ptr=con ; ii < NNAME && *ptr != '\0' && *ptr != '\n' ; ptr++ ){
01334         rtin->name_data[ii++] = *ptr ;
01335      }
01336      if( ii >= NNAME ){
01337         fprintf(stderr,"RT: control image_channel_spec buffer overflow!\a\n") ;
01338         ii = NNAME - 1 ;
01339      }
01340      rtin->name_data[ii] = '\0' ;
01341 
01342      /** open data channel **/
01343 
01344      rtin->ioc_data = iochan_init( rtin->name_data , "accept" ) ;
01345      if( rtin->ioc_data == NULL ){
01346         fprintf(stderr,"RT: failure to open data IOCHAN %s\a\n",rtin->name_data) ;
01347         free(rtin) ; free(con) ; return NULL ;
01348      }
01349 
01350      if( verbose == 2 )
01351         fprintf(stderr,"RT: opened data stream %s\n",rtin->name_data) ;
01352      VMCHECK ;
01353 
01354      /** if more follows, that is a command for a child process
01355          (which will not be started until the first image data arrives) **/
01356 
01357      ii = 0 ;
01358      if( *ptr != '\0' ){
01359         for( ptr++ ; ii < NNAME && *ptr != '\0' && *ptr != '\n' ; ptr++ ){
01360            rtin->name_info[ii++] = *ptr ;
01361         }
01362         if( ii >= NNAME ){
01363            fprintf(stderr,"RT: control info_command buffer overflow!\a\n") ;
01364            ii = NNAME - 1 ;
01365         }
01366      }
01367      rtin->name_info[ii] = '\0' ;
01368 
01369      if( verbose == 2 ){
01370         if( strlen(rtin->name_info) > 0 )
01371            fprintf(stderr,"RT: info command for child will be '%s'\n",rtin->name_info) ;
01372         else
01373            fprintf(stderr,"RT: no info command given.\n") ;
01374      }
01375 
01376      free(con) ;
01377      VMCHECK ;
01378 
01379      /** the command (if any) will be run later **/
01380 
01381    } else {  /*** THE NEW WAY: directly connect to ioc_data channel [10 DEC 2002] ***/
01382 
01383      /** make new structure **/
01384 
01385      rtin = (RT_input *) calloc( 1 , sizeof(RT_input) ) ;
01386      if( rtin == NULL ){
01387        fprintf(stderr,"RT: malloc fails in new_RT_input!\a\n") ; EXIT(1) ;
01388      }
01389 
01390      MCW_strncpy( rtin->name_data , ioc_data->name , NNAME ) ;  /* not really needed */
01391      rtin->ioc_data     = ioc_data ;
01392      rtin->name_info[0] = '\0' ;                                /* no child */
01393 
01394    } /*** end of THE NEW WAY of starting up ***/
01395 
01396    rtin->ioc_info   = NULL ;
01397    rtin->child_info = 0 ;
01398 
01399    /** wait until the data channel is good **/
01400 
01401    if( verbose )
01402       fprintf(stderr,"RT: waiting for data stream to become good.\n") ;
01403    VMCHECK ;
01404 
01405    while(1){
01406       ii = iochan_goodcheck(rtin->ioc_data,1000) ;             /* wait up to 1000 msec */
01407            if( ii >  0 )                 break ;               /* is good! */
01408       else if( ii == 0 && verbose == 2 ) fprintf(stderr,".") ; /* not good yet */
01409       else {                                                   /* is bad! */
01410          fprintf(stderr,"RT: data stream fails to become good!\a\n") ;
01411          IOCHAN_CLOSENOW(rtin->ioc_data) ; free(rtin) ;
01412          return NULL ;
01413       }
01414    }
01415 
01416    if( verbose == 2 )
01417       fprintf(stderr,"RT: data stream is now bodaciously good.\n") ;
01418    VMCHECK ;
01419 
01420    /** initialize internal constants in the struct **/
01421 
01422    rtin->info_ok = 0 ; rtin->no_data = 1 ;
01423 
01424    rtin->image_mode   = image_mode ;  /* 28 Apr 2000 */
01425    rtin->image_handle = NULL ;
01426    rtin->image_space  = NULL ;
01427 
01428 #if 0                               /*** THE VERY OLD WAY to set things up ***/
01429    rtin->nxx = DEFAULT_XYMATRIX ;
01430    rtin->nyy = DEFAULT_XYMATRIX ;
01431    rtin->nzz = DEFAULT_ZNUM ;
01432 
01433    rtin->orcxx = ORI_S2I_TYPE ;
01434    rtin->orcyy = ORI_A2P_TYPE ;
01435    rtin->orczz = ORI_L2R_TYPE ;
01436 
01437    rtin->xxfov = DEFAULT_XYFOV ;
01438    rtin->yyfov = DEFAULT_XYFOV ;
01439    rtin->zzfov = -1.0 ;
01440    rtin->dxx   = DEFAULT_XYFOV / DEFAULT_XYMATRIX ;
01441    rtin->dyy   = DEFAULT_XYFOV / DEFAULT_XYMATRIX ;
01442    rtin->dzz   = DEFAULT_ZDELTA ;
01443 
01444    rtin->zgap  = 0.0 ;                  /* 18 Dec 2002 */
01445    rtin->xxoff = rtin->yyoff = rtin->zzoff = 0.0 ;
01446 
01447    rtin->xxorg = 0.5 * (rtin->nxx - 1) * rtin->dxx ; rtin->xcen = 1 ;
01448    rtin->yyorg = 0.5 * (rtin->nyy - 1) * rtin->dyy ; rtin->ycen = 1 ;
01449    rtin->zzorg = 0.5 * (rtin->nzz - 1) * rtin->dzz ; rtin->zcen = 1 ;
01450    rtin->xxdcode = ILLEGAL_TYPE ;
01451    rtin->yydcode = ILLEGAL_TYPE ;
01452    rtin->zzdcode = ILLEGAL_TYPE ;
01453 
01454    rtin->tr     = DEFAULT_TR ;
01455    rtin->nr     = 1 ;
01456    rtin->dtype  = DTYPE_2DZT ;
01457    rtin->datum  = MRI_short  ;
01458 
01459    rtin->zorder      = ZORDER_ALT ;
01460    rtin->zorder_lock = 0 ;              /* 22 Feb 1999 */
01461 #else
01462    rtin->nxx = -1 ;
01463    rtin->nyy = -1 ;
01464    rtin->nzz = -1 ;
01465 
01466    rtin->orcxx = -1 ;
01467    rtin->orcyy = -1 ;
01468    rtin->orczz = -1 ;
01469 
01470    rtin->xxfov = 0.0 ;
01471    rtin->yyfov = 0.0 ;
01472    rtin->zzfov = 0.0 ;
01473    rtin->dxx   = 0.0 ;
01474    rtin->dyy   = 0.0 ;
01475    rtin->dzz   = 0.0 ;
01476 
01477    rtin->zgap  = 0.0 ;                  /* 18 Dec 2002 */
01478    rtin->xxoff = rtin->yyoff = rtin->zzoff = 0.0 ;
01479 
01480    rtin->xxorg = 0.0 ; rtin->xcen = 1 ; rtin->xxdcode = ILLEGAL_TYPE ;
01481    rtin->yyorg = 0.0 ; rtin->ycen = 1 ; rtin->yydcode = ILLEGAL_TYPE ;
01482    rtin->zzorg = 0.0 ; rtin->zcen = 1 ; rtin->zzdcode = ILLEGAL_TYPE ;
01483 
01484    rtin->tr     = DEFAULT_TR ;
01485    rtin->nr     = 1 ;
01486    rtin->dtype  = DTYPE_2DZT ;
01487    rtin->datum  = MRI_short  ;
01488    rtin->nzseq  = 0 ;
01489 
01490    rtin->zorder      = ZORDER_ALT ;
01491    rtin->zorder_lock = 0 ;              /* 22 Feb 1999 */
01492    rtin->tpattern    = ZORDER_ALT ;     /* 10 May 2005 [rickr] */
01493 #endif
01494 
01495    rtin->swap_on_read = 0 ;  /* wait for BYTEORDER     26 Jun 2003 [rickr] */
01496 
01497    rtin->nbuf   = 0    ;     /* no input buffer data yet */
01498    rtin->imsize = 0 ;        /* don't know how many bytes/image yet */
01499    rtin->bufar  = NULL ;     /* no input image buffer yet */
01500 
01501    im3d = plint->im3d ;                           /* default controller */
01502    if( !IM3D_OPEN(im3d) )                         /* may not be open */
01503       im3d = AFNI_find_open_controller() ;
01504 
01505    /** 01 Aug 2002: initialize dataset dependent items */
01506 
01507    for( cc=0 ; cc < MAX_CHAN ; cc++ ){
01508       rtin->afni_status[cc] = 0 ;      /* AFNI is ignorant of the dataset */
01509       rtin->dset[cc]        = NULL ;   /* no dataset yet */
01510       rtin->sbr[cc]         = NULL ;   /* no brick space yet */
01511       rtin->im[cc]          = NULL ;
01512       rtin->nvol[cc]        = 0 ;
01513       rtin->nsl[cc]         = 0 ;
01514       rtin->im3d[cc]        = im3d ;
01515       rtin->sess[cc]        = GLOBAL_library.sslist->ssar[im3d->vinfo->sess_num] ;
01516       rtin->sess_num[cc]    = im3d->vinfo->sess_num ;
01517    }
01518 
01519    rtin->num_chan = 1 ;   /* default number of channels */
01520    rtin->cur_chan = 0 ;   /* current channel index */
01521 
01522    strcpy( rtin->root_prefix , root ) ;
01523 
01524    /** setup realtime function evaluation **/
01525 
01526    rtin->func_dset   = NULL ;  /* no function yet (if ever) */
01527    rtin->func_status = 0 ;     /* AFNI is ignorant */
01528    rtin->func_condit = 0 ;     /* no condition yet */
01529 
01530    /* function evaluation is not done on time-independent datasets */
01531 
01532    rtin->func_code   = (rtin->dtype == DTYPE_2DZ || rtin->dtype == DTYPE_3D)
01533                        ? FUNC_NONE : func_code ;
01534 
01535    rtin->func_func   = FUNC_funcs[rtin->func_code] ; /* where to evaluate */
01536 
01537 #ifdef ALLOW_REGISTRATION
01538    rtin->reg_base_index = regtime ;  /* save these now, in case the evil */
01539    rtin->reg_mode       = regmode ;  /* user changes them on the fly!    */
01540    rtin->reg_dset       = NULL ;
01541    rtin->reg_2dbasis    = NULL ;
01542    rtin->reg_status     = 0 ;        /* AFNI knows nothing. NOTHING.     */
01543    rtin->reg_nvol       = 0 ;        /* number volumes registered so far */
01544 
01545    rtin->reg_nest  = 0 ;
01546    rtin->reg_tim   = (float *) malloc( sizeof(float) ) ;
01547    rtin->reg_dx    = (float *) malloc( sizeof(float) ) ;
01548    rtin->reg_dy    = (float *) malloc( sizeof(float) ) ;
01549    rtin->reg_phi   = (float *) malloc( sizeof(float) ) ;
01550 
01551    rtin->reg_dz       = (float *) malloc( sizeof(float) ) ;
01552    rtin->reg_theta    = (float *) malloc( sizeof(float) ) ;
01553    rtin->reg_psi      = (float *) malloc( sizeof(float) ) ;
01554    rtin->reg_rep      = (float *) malloc( sizeof(float) ) ;
01555    rtin->reg_eval     = (float *) malloc( sizeof(float) ) ;
01556    rtin->reg_3dbasis  = NULL ;
01557    rtin->mp           = NULL ;    /* no plot yet */
01558 
01559    rtin->reg_graph_xnew = 0 ;     /* did the user update reg_graph_xr */
01560    rtin->reg_graph_ynew = 0 ;     /* did the user update reg_graph_yr */
01561    rtin->reg_graph_xr = reg_nr ;  /* will scale by TR when we know it */
01562    rtin->reg_graph_yr = reg_yr ;
01563 
01564    rtin->p_code = NULL ;          /* init parser code    12 Feb 2004 [rickr] */
01565 
01566    rtin->mp_tcp_use = 0 ;         /* tcp motion param    30 Mar 2004 [rickr] */
01567    rtin->mp_tcp_sd  = 0 ;
01568    rtin->mp_port    = RT_MP_DEF_PORT ;
01569    rtin->mp_nmsg    = 0 ;
01570    rtin->mp_npsets  = 0 ;
01571    strcpy(rtin->mp_host, "localhost") ;
01572 
01573    rtin->reg_resam = REG_resam_ints[reg_resam] ;
01574    if( rtin->reg_resam < 0 ){                    /* 20 Nov 1998: */
01575       rtin->reg_resam       = MRI_HEPTIC ;       /* special case */
01576       rtin->reg_final_resam = MRI_FOURIER ;
01577    } else {
01578       rtin->reg_final_resam = -1 ;
01579    }
01580 
01581    rtin->reg_graph = reggraph ;
01582    if( regmode==REGMODE_3D_ESTIM && reggraph==0 ) rtin->reg_graph = 1;
01583 #endif
01584 
01585    /** record the times for later reportage **/
01586 
01587    rtin->elapsed = PLUTO_elapsed_time() ; rtin->last_elapsed = rtin->elapsed ;
01588    rtin->cpu     = PLUTO_cpu_time() ;     rtin->last_nvol = 0 ;
01589 
01590    rtin->num_note = 0 ;      /* 01 Oct 2002: notes list */
01591    rtin->note     = NULL ;
01592 
01593    rtin->marked_for_death = 0 ;  /* 10 Dec 2002 */
01594 
01595    return rtin ;
01596 }
01597 
01598 /*-------------------------------------------------------------------------
01599    Create a child process to run a command that will send header info
01600    back to the parent AFNI process.
01601 ---------------------------------------------------------------------------*/
01602 
01603 void RT_start_child( RT_input * rtin )
01604 {
01605    pid_t child_pid ;
01606 
01607    if( rtin == NULL || strlen(rtin->name_info) == 0 ) return ;  /* sanity */
01608 
01609    child_pid = fork() ;             /* AKA bifurcation */
01610 
01611    if( child_pid == (pid_t)(-1) ){  /* real bad news */
01612      fprintf(stderr,"RT: can't fork child process!\a\n") ; EXIT(1) ;
01613    }
01614 
01615    if( child_pid > 0 ){              /** I'm the parent **/
01616 
01617       if( verbose == 2 )
01618          fprintf(stderr,"RT: forked a child process to execute '%s'\n",rtin->name_info) ;
01619       VMCHECK ;
01620 
01621       /** open a channel to communicate with the child **/
01622 
01623       rtin->child_info = child_pid ;
01624       rtin->ioc_info   = iochan_init( SHM_CHILD , "accept" ) ;
01625       if( rtinp->ioc_info == NULL ){
01626         kill( child_pid , SIGTERM ) ;
01627         fprintf(stderr,"RT: can't create read stream from child!\a\n") ;
01628         EXIT(1) ;
01629       }
01630 
01631       rtin->child_start_time = PLUTO_elapsed_time() ;  /* 10 Dec 1998 */
01632 
01633       /** next time thru RT_worker will try to read data from this channel **/
01634 
01635    } else {                          /** I'm the child **/
01636 
01637       RT_acquire_info( rtin->name_info ) ;  /* get info, send back to parent */
01638       _exit(0) ;
01639    }
01640 
01641    return ;
01642 }
01643 
01644 /*-------------------------------------------------------------------------
01645 ---------------------------------------------------------------------------
01646    This routine is only called from the child process:
01647      Run a command, get its output, send it back to the parent, exit.
01648      The entire existence of the child is confined to this routine.
01649 ---------------------------------------------------------------------------
01650 ---------------------------------------------------------------------------*/
01651 
01652 int RT_acquire_info( char * command )
01653 {
01654    FILE * fp ;
01655    char * info = (char *) malloc( sizeof(char) * INFO_SIZE ) ; int ninfo = 0 ;
01656    IOCHAN * ioc ;
01657    int jj ;
01658 
01659    /** create channel back to parent **/
01660 
01661    ioc = iochan_init( SHM_CHILD , "create" ) ;
01662    if( ioc == NULL ){
01663       fprintf(stderr,"RT: child fails to open stream back to parent!\a\n") ;
01664       _exit(1) ;
01665    }
01666 
01667    /** run command, read its output into this pipe **/
01668 
01669    fp = popen( command , "r" ) ;
01670    if( fp == NULL ){
01671       fprintf(stderr,"RT: child fails to open pipe to command=%s\a\n",command) ;
01672       IOCHAN_CLOSENOW(ioc) ; _exit(1) ;
01673    }
01674 
01675    /** read pipe until nothing more **/
01676 
01677    while( fgets(info+ninfo,INFO_SIZE-ninfo,fp) != NULL ){
01678       ninfo = strlen(info) ;
01679    }
01680    pclose(fp) ;
01681 
01682    /** send output back to parent **/
01683 
01684    jj = iochan_writecheck(ioc,-1) ;  /* wait until ready */
01685    if( jj < 0 ){
01686       fprintf(stderr,"RT: child can't write IOCHAN to parent!\a\n") ;
01687       IOCHAN_CLOSENOW(ioc) ; _exit(1) ;
01688    }
01689 
01690    iochan_sendall( ioc , info , ninfo+1 ) ;        /* include the NUL character */
01691    iochan_sleep(LONG_DELAY) ;                      /* wait a bit */
01692    while( ! iochan_clearcheck(ioc,LONG_DELAY) )    /* loop until cleared */
01693       iochan_sleep(LONG_DELAY) ;
01694 
01695    iochan_sleep(LONG_DELAY) ;                      /* once more, just for luck */
01696 
01697                                                    /**************************/
01698    free(info); IOCHAN_CLOSENOW(ioc); _exit(0);     /** END OF CHILD PROCESS **/
01699 }                                                  /**************************/
01700 
01701 /****************************************************************************
01702    Check the realtime input structure for OK-ness, and set its internal flag.
01703    If if it bad, and prt isn't 0, print a message about the error(s).
01704 *****************************************************************************/
01705 
01706 #define EPR(s) fprintf(stderr,"RT: HEADER DATA ERROR - %s\a\n",(s))
01707 
01708 #define OR3OK(x,y,z) ( ((x)&6) + ((y)&6) + ((z)&6) == 6 )
01709 
01710 void RT_check_info( RT_input * rtin , int prt )
01711 {
01712    if( rtin == NULL ) return ;
01713 
01714    /*-- 28 Apr 2000: if in Image Only mode, do fewer checks --*/
01715 
01716    if( rtin->image_mode ){
01717 
01718       rtin->info_ok = ( rtin->nxx > 1 )                         &&
01719                       ( rtin->nyy > 1 )                         &&
01720                       ( AFNI_GOOD_DTYPE(rtin->datum) ) ;
01721 
01722       if( rtin->info_ok || !prt ) return ;  /* if good, or if no print */
01723 
01724       if( !(rtin->nxx > 1)                ) EPR("Image x-dimen not > 1") ;
01725       if( !(rtin->nyy > 1)                ) EPR("Image y-dimen not > 1") ;
01726       if( !(AFNI_GOOD_DTYPE(rtin->datum)) ) EPR("Bad datum") ;
01727       return ;
01728    }
01729 
01730    /*-- below here: must construct dataset, so do all necessary checks --*/
01731 
01732    rtin->info_ok = ( rtin->dtype > 0 )                            &&
01733                    ( THD_filename_pure(rtin->root_prefix) )       &&
01734                    ( strlen(rtin->root_prefix) < THD_MAX_PREFIX ) &&
01735                    ( rtin->tr > 0 )                               &&
01736                    ( rtin->dzz > 0 || rtin->zzfov > 0 )           &&
01737                    ( rtin->xxfov > 0 )                            &&
01738                    ( rtin->yyfov > 0 )                            &&
01739                    ( rtin->nxx > 1 )                              &&
01740                    ( rtin->nyy > 1 )                              &&
01741                    ( rtin->nzz >= 1 )                             &&
01742                    ( AFNI_GOOD_DTYPE(rtin->datum) )               &&
01743                    ( rtin->zorder > 0 )                           &&
01744                    ( rtin->tpattern > 0 )                         &&
01745                    ( rtin->orcxx >= 0 )                           &&
01746                    ( rtin->orcyy >= 0 )                           &&
01747                    ( rtin->orczz >= 0 )                           &&
01748                    ( OR3OK(rtin->orcxx,rtin->orcyy,rtin->orczz) )    ;
01749 
01750    if( rtin->info_ok || !prt ) return ;  /* if good, or if no print */
01751 
01752    /* print error messages */
01753 
01754    if( !(rtin->dtype > 0)                            ) EPR("Bad acquisition type") ;
01755    if( !(THD_filename_pure(rtin->root_prefix))       ) EPR("Bad prefix") ;
01756    if( !(strlen(rtin->root_prefix) < THD_MAX_PREFIX) ) EPR("Overlong prefix") ;
01757    if( !(rtin->tr > 0)                               ) EPR("TR is not positive") ;
01758    if( !(rtin->dzz > 0 || rtin->zzfov > 0)           ) EPR("Slice thickness not positive") ;
01759    if( !(rtin->xxfov > 0)                            ) EPR("x-FOV not positive") ;
01760    if( !(rtin->yyfov > 0)                            ) EPR("y-FOV not positive") ;
01761    if( !(rtin->nxx > 1)                              ) EPR("Image x-dimen not > 1") ;
01762    if( !(rtin->nyy > 1)                              ) EPR("Image y-dimen not > 1") ;
01763    if( !(rtin->nzz >= 1)                             ) EPR("Slice count (z-dimen) not >= 1") ;
01764    if( !(AFNI_GOOD_DTYPE(rtin->datum))               ) EPR("Bad datum") ;
01765    if( !(rtin->zorder > 0)                           ) EPR("Slice ordering illegal") ;
01766    if( !(rtin->tpattern > 0)                         ) EPR("Timing pattern illegal") ;
01767    if( !(rtin->orcxx >= 0)                           ) EPR("x-orientation illegal") ;
01768    if( !(rtin->orcyy >= 0)                           ) EPR("y-orientation illegal") ;
01769    if( !(rtin->orczz >= 0)                           ) EPR("z-orientation illegal") ;
01770    if( !(OR3OK(rtin->orcxx,rtin->orcyy,rtin->orczz)) ) EPR("Inconsistent xyz-orientations") ;
01771 
01772    return ;
01773 }
01774 
01775 #ifdef ALLOW_REGISTRATION
01776 /*---------------------------------------------------------------------------
01777    Close the socket connection.                        30 Mar 2004 [rickr]
01778 
01779    return   0 : on success
01780           < 0 : on error
01781 -----------------------------------------------------------------------------*/
01782 int RT_mp_comm_close( RT_input * rtin )
01783 {
01784     char magic_bye[] = { 0xde, 0xad, 0xde, 0xad, 0 };
01785 
01786     if ( rtin->mp_tcp_use != 1 || rtin->mp_tcp_sd <= 0 )
01787         return 0;
01788 
01789     if ( (tcp_writecheck(rtin->mp_tcp_sd, 1)   == -1) ||
01790          (send(rtin->mp_tcp_sd, magic_bye, 4, 0) == -1 ) )
01791         fprintf(stderr,"** closing: our MP socket has gone bad?\n");
01792 
01793     fprintf(stderr,"RT: MP: closing motion param socket, "
01794                    "sent %d param sets over %d messages\n",
01795                    rtin->mp_npsets, rtin->mp_nmsg);
01796 
01797     /* in any case, close the socket */
01798     close(rtin->mp_tcp_sd);
01799     rtin->mp_tcp_sd  = 0;
01800     rtin->mp_tcp_use = 0;
01801     rtin->mp_npsets  = 0;
01802     rtin->mp_nmsg    = 0;
01803 
01804     return 0;
01805 }
01806 
01807 
01808 /*---------------------------------------------------------------------------
01809    Send the current motion params.                        30 Mar 2004 [rickr]
01810 
01811    return   0 : on success
01812           < 0 : on error
01813 -----------------------------------------------------------------------------*/
01814 int RT_mp_comm_send_data( RT_input * rtin, float * mp[6], int nt )
01815 {
01816     float data[600];            /* max transfer is nt == 100 */
01817     int   rv, nvals, remain;
01818     int   c, c2;
01819 
01820     if ( rtin->mp_tcp_use != 1 || nt <= 0 )
01821         return 0;
01822 
01823     if ( rtin->mp_tcp_sd <= 0 )
01824         return -1;
01825 
01826     /* hmmmm, Bob has a good function to test the socket... */
01827     if ( (rv = tcp_writecheck(rtin->mp_tcp_sd, 1)) == -1 )
01828     {
01829         fprintf(stderr,"** our MP socket has gone bad?\n");
01830         close(rtin->mp_tcp_sd);
01831         rtin->mp_tcp_sd  = 0;
01832         rtin->mp_tcp_use = 0;   /* allow a later re-try... */
01833         return -1;
01834     }
01835 
01836     remain = nt;
01837     while ( remain > 0 )
01838     {
01839         nvals = MIN(remain, 100);
01840 
01841         /* copy floats to 'data' */
01842         for ( c = 0; c < nvals; c++ )
01843             for ( c2 = 0; c2 < 6; c2++ )
01844                 data[6*c+c2] = mp[c2][c];
01845 
01846         if ( send(rtin->mp_tcp_sd, data, 6*nvals*sizeof(float), 0) == -1 )
01847         {
01848             fprintf(stderr,"** failed to send %d floats, closing socket...\n",
01849                     6*nvals);
01850             close(rtin->mp_tcp_sd);
01851             rtin->mp_tcp_sd  = 0;
01852             rtin->mp_tcp_use = 0;  /* allow a later re-try... */
01853             return -1;
01854         }
01855 
01856         /* keep track of num messages and num param sets */
01857         rtin->mp_nmsg++;
01858         rtin->mp_npsets += nvals;
01859 
01860         remain -= nvals;
01861     }
01862 
01863     return 0;
01864 }
01865 
01866 
01867 /*---------------------------------------------------------------------------
01868    Initialize the motion parameter communications.        30 Mar 2004 [rickr]
01869 
01870    return   0 : on success
01871           < 0 : on error
01872 -----------------------------------------------------------------------------*/
01873 int RT_mp_comm_init( RT_input * rtin )
01874 {
01875     struct sockaddr_in   sin;
01876     struct hostent     * hostp;
01877     char                 magic_hi[] = { 0xab, 0xcd, 0xef, 0xab };
01878     int                  sd;
01879 
01880     if ( rtin->mp_tcp_sd != 0 )
01881         fprintf(stderr,"** warning, did we not close the MP socket?\n");
01882 
01883     if ( (hostp = gethostbyname(rtin->mp_host)) == NULL )
01884     {
01885         fprintf(stderr,"** cannot lookup host '%s'\n", rtin->mp_host);
01886         rtin->mp_tcp_use = -1;
01887         return -1;
01888     }
01889 
01890     /* fill the sockaddr_in struct */
01891     memset(&sin, 0, sizeof(sin));
01892     sin.sin_family      = AF_INET;
01893     sin.sin_addr.s_addr = ((struct in_addr *)(hostp->h_addr))->s_addr;
01894     sin.sin_port        = htons(rtin->mp_port);
01895 
01896     /* get a socket */
01897     if ( (sd = socket(AF_INET, SOCK_STREAM, 0)) == -1 )
01898     {
01899         perror("pe: socket");
01900         rtin->mp_tcp_use = -1;   /* let us not try, try again */
01901         return -1;
01902     }
01903 
01904     if ( connect(sd, (struct sockaddr *)&sin, sizeof(sin)) == -1 )
01905     {
01906         perror("pe: connect");
01907         rtin->mp_tcp_use = -1;
01908         return -1;
01909     }
01910 
01911     /* send the hello message */
01912     if ( send(sd, magic_hi, 4*sizeof(char), 0) == -1 )
01913     {
01914         perror("pe: send hello");
01915         rtin->mp_tcp_use = -1;
01916         return -1;
01917     }
01918 
01919     fprintf(stderr,"RT: MP: opened motion param socket to %s:%d\n",
01920             rtin->mp_host, rtin->mp_port);
01921 
01922     /* everything worked out, we're good to van Gogh */
01923 
01924     rtin->mp_tcp_sd = sd;
01925 
01926     return 0;
01927 }
01928 
01929 
01930 /*---------------------------------------------------------------------------
01931    Initialize the motion parameter communication variables. 30 Mar 2004 [rickr]
01932 
01933    We are expecting AFNI_REALTIME_MP_HOST_PORT to read "hostname:port".
01934 
01935    return   0 : on success
01936           < 0 : on error
01937 -----------------------------------------------------------------------------*/
01938 int RT_mp_comm_init_vars( RT_input * rtin )
01939 {
01940     char * ept, * cp;
01941     int    len;
01942 
01943     if ( rtin->mp_tcp_use < 0 )     /* we've failed out, do not try again */
01944         return 0;
01945 
01946     if ( rtin->mp_tcp_sd != 0 )
01947         fprintf(stderr,"** warning, did we not close the MP socket?\n");
01948     rtin->mp_tcp_sd   = 0;
01949 
01950     /* for now, we will only init this if the HOST:PORT env var exists */
01951     ept = getenv("AFNI_REALTIME_MP_HOST_PORT") ;  /* 09 Oct 2000 */
01952     if( ept == NULL )
01953         return 0;
01954 
01955     cp = strchr(ept, ':');      /* find ':' seperator */
01956 
01957     if ( cp == NULL || !isdigit(*(cp+1)) )
01958     {
01959         fprintf(stderr,"** env var AFNI_REALTIME_MP_HOST_PORT must be in the "
01960                        "form hostname:port_num\n   (var is '%s')\n", ept);
01961         return -1;
01962     }
01963 
01964     len = cp - ept;     /* length of hostname */
01965     if ( len > 127 )
01966     {
01967         fprintf(stderr,"** motion param hostname restricted to 127 bytes,\n"
01968                        "   found %d from host in %s\n", len, ept);
01969         return -1;
01970     }
01971 
01972     fprintf(stderr,"RT: MP: found motion param env var '%s'\n", ept);
01973 
01974     rtin->mp_port = atoi(cp+1);
01975     strncpy(rtin->mp_host, ept, len);
01976     rtin->mp_host[len] = '\0';
01977     rtin->mp_tcp_use = 1;
01978 
01979     return 0;
01980 }
01981 
01982 
01983 /*---------------------------------------------------------------------------
01984    Initialize the parser fields from the user expression.  12 Feb 2004 [rickr]
01985 
01986    return   0 : on success
01987           < 0 : on error
01988 -----------------------------------------------------------------------------*/
01989 int RT_parser_init( RT_input * rtin )
01990 {
01991     PARSER_set_printout(1);
01992     rtin->p_code = PARSER_generate_code( rtin->p_expr );
01993 
01994     if ( ! rtin->p_code )
01995     {
01996         fprintf(stderr,"** cannot parse expression '%s'\n", rtin->p_expr);
01997         return -1;
01998     }
01999 
02000     /* p_max_sym will be 0 if nothing is marked, 26 if z is, etc. */
02001     PARSER_mark_symbols( rtin->p_code, rtin->p_has_sym );
02002     for ( rtin->p_max_sym = 26; rtin->p_max_sym > 0; rtin->p_max_sym-- )
02003         if ( rtin->p_has_sym[rtin->p_max_sym - 1] )
02004             break;
02005 
02006     if ( rtin->p_max_sym > 6 )
02007     {
02008         fprintf(stderr,"** parser expression may only contain symbols a-f\n");
02009         return -2;
02010     }
02011 
02012     return 0;
02013 }
02014 #endif
02015 
02016 /*---------------------------------------------------------------------------
02017    Process command strings in the buffer, up to the first '\0' character.
02018    Returns the number of characters processed, which will be one more than
02019    the index of the first '\0'.  If an error occurs, -1 is returned.
02020 -----------------------------------------------------------------------------*/
02021 
02022 #define BADNEWS     fprintf(stderr,"RT: illegal header info=%s\a\n",buf)
02023 #define STARTER(st) (strncmp(buf,st,strlen(st)) == 0)
02024 #define NBUF        1024
02025 
02026 int RT_process_info( int ninfo , char * info , RT_input * rtin )
02027 {
02028    int ii , jj , nstart,nend , nuse , nbuf ;
02029    char buf[NBUF] ;
02030 
02031    if( rtin == NULL || info == NULL || ninfo == 0 ) return -1 ;
02032 
02033    for( nend=0 ; nend < ninfo && info[nend] != '\0' ; nend++ ) ; /* nada */
02034    if( nend == ninfo ){
02035       fprintf(stderr,"RT: info string not NUL-terminated!\a\n") ;
02036       return -1 ;
02037    }
02038 
02039    /******************************************/
02040    /** Scan ahead until next command string **/
02041 
02042    nstart = 0 ;
02043    while( nstart < nend ){
02044 
02045       /** skip whitespace **/
02046 
02047       for( ; nstart < nend && isspace(info[nstart]) ; nstart++ ) ; /* nada */
02048       if( nstart >= nend ) break ;
02049 
02050       /** copy characters into buf until the next '\n' or '\0' **/
02051 
02052       for( nbuf=0,jj=nstart ; nbuf < NBUF && info[jj] != '\n' && info[jj] != '\0' ; ){
02053          buf[nbuf++] = info[jj++] ;
02054       }
02055       if( nbuf == NBUF ){
02056          fprintf(stderr,"RT: line buffer overflow in control information!\a\n") ;
02057          nbuf-- ;
02058       }
02059       buf[nbuf] = '\0' ; nstart = jj ;
02060 
02061       if( verbose == 2 )
02062          fprintf(stderr,"RT: info line buffer=%s\n",buf) ;
02063       VMCHECK ;
02064 
02065       /************************************/
02066       /*** Scan for legal input strings ***/
02067 
02068       if( STARTER("ACQUISITION_TYPE") ){
02069          char typ[32] ;
02070 
02071          sscanf( buf , "ACQUISITION_TYPE %31s" , typ ) ;
02072 
02073               if( strcmp(typ,"2D+z")  == 0 ) rtin->dtype = DTYPE_2DZ ;
02074          else if( strcmp(typ,"2D+zt") == 0 ) rtin->dtype = DTYPE_2DZT ;
02075          else if( strcmp(typ,"3D")    == 0 ) rtin->dtype = DTYPE_3D ;
02076          else if( strcmp(typ,"3D+t")  == 0 ) rtin->dtype = DTYPE_3DT ;
02077          else
02078               BADNEWS ;
02079 
02080 #ifdef ALLOW_REGISTRATION
02081 
02082       /* Some troublesome physicist, let's just call him "Tom" (possibly
02083        * named by his mother, Mrs. Ross), wants control over the ranges
02084        * in the motion correction graph window.
02085        *                                                29 Jan 2004 [rickr] */
02086 
02087       } else if( STARTER("GRAPH_XRANGE") ){
02088          float fval = 0.0 ;
02089          sscanf( buf , "GRAPH_XRANGE %f" , &fval ) ;
02090          if( fval >= MIN_PIN && fval <= MAX_PIN ) {
02091              rtin->reg_graph_xnew = 1;
02092              rtin->reg_graph_xr   = fval;
02093 
02094              if( rtin->reg_graph && REG_IS_3D(rtin->reg_mode) &&
02095                  IM3D_OPEN(plint->im3d) )
02096              {
02097                  plot_ts_xypush(1-rtin->reg_graph_xnew, 1-rtin->reg_graph_ynew);
02098                  RT_set_grapher_pinnums((int)(fval+0.5));
02099              }
02100 
02101          } else
02102               BADNEWS ;
02103 
02104       } else if( STARTER("GRAPH_YRANGE") ){
02105          float fval = 0.0 ;
02106          sscanf( buf , "GRAPH_YRANGE %f" , &fval ) ;
02107          if( fval > 0.0 ) {
02108             rtin->reg_graph_ynew = 1;
02109             rtin->reg_graph_yr   = fval ;
02110 
02111             /* if the user sets scales, don't 'push'    11 Feb 2004 [rickr] */
02112             if( rtin->reg_graph && REG_IS_3D(rtin->reg_mode) )
02113                 plot_ts_xypush(1-rtin->reg_graph_xnew, 1-rtin->reg_graph_ynew);
02114          } else
02115             BADNEWS ;
02116 
02117       /* Allow the user to specify an expression, to evalue the six motion
02118        * parameters into one.
02119        *                                                12 Feb 2004 [rickr] */
02120       } else if( STARTER("GRAPH_EXPR") ){
02121          sscanf( buf , "GRAPH_EXPR %1024s" , rtin->p_expr ) ;
02122          rtin->p_expr[RT_MAX_EXPR] = '\0';
02123          if ( RT_parser_init(rtin) != 0 )
02124             BADNEWS ;
02125 #endif
02126 
02127       } else if( STARTER("NAME") ){
02128          char npr[THD_MAX_PREFIX] = "\0" ;
02129          /* RT_MAX_PREFIX is used here, currently 100 */
02130          sscanf( buf , "NAME %100s" , npr ) ; /* 31->100  29 Jan 2004 [rickr] */
02131          if( THD_filename_pure(npr) ) strcpy( rtin->root_prefix , npr ) ;
02132          else
02133               BADNEWS ;
02134 
02135       } else if( STARTER("PREFIX") ){           /* 01 Aug 2002 */
02136          char npr[THD_MAX_PREFIX] = "\0" ;
02137          /* RT_MAX_PREFIX is used here, currently 100 */
02138          sscanf( buf , "PREFIX %100s" , npr ) ;
02139          if( THD_filename_pure(npr) ) strcpy( rtin->root_prefix , npr ) ;
02140          else
02141               BADNEWS ;
02142 
02143       } else if( STARTER("NOTE") ) {            /* 01 Oct 2002: notes list */
02144          int nn = rtin->num_note ;
02145          if( nbuf > 6 ){
02146            int ii ;
02147            rtin->note = realloc( rtin->note , sizeof(char *)*(nn+1) ); /* extend array */
02148            rtin->note[nn] = strdup(buf+5) ;                            /* skip "NOTE "  */
02149            for( ii=0 ; rtin->note[nn][ii] != '\0' ; ii++ )             /* '\n' insertion */
02150              if( rtin->note[nn][ii] == '\a' || rtin->note[nn][ii] == '\f' )
02151                rtin->note[nn][ii] = '\n';
02152            rtin->num_note ++ ;
02153          } else
02154            BADNEWS ;
02155 
02156       } else if( STARTER("NUMVOL") ){
02157          int val = 0 ;
02158          sscanf( buf , "NUMVOL %d" , &val ) ;
02159          if( val > 0 ) rtin->nr = val ;
02160          else
02161               BADNEWS ;
02162 
02163       } else if( STARTER("TR") ){          /* units are seconds */
02164          float val = 0.0 ;
02165          sscanf( buf , "TR %f" , &val ) ;
02166          if( val > 0.0 ) rtin->tr = val ;
02167          else
02168               BADNEWS ;
02169 
02170       } else if( STARTER("ZDELTA") ){
02171          float val = 0.0 ;
02172          sscanf( buf , "ZDELTA %f" , &val ) ;
02173          if( val > 0.0 ) rtin->dzz = val ;
02174          else
02175               BADNEWS ;
02176          if( verbose == 2 )
02177             fprintf(stderr,"RT: dzz = %g\n",rtin->dzz) ;
02178          VMCHECK ;
02179 
02180       } else if( STARTER("ZGAP") ){               /* 18 Dec 2002 */
02181          float val = 0.0 ;
02182          sscanf( buf , "ZGAP %f" , &val ) ;
02183          if( val >= 0.0 ) rtin->zgap = val ;
02184          else
02185               BADNEWS ;
02186          if( verbose == 2 )
02187             fprintf(stderr,"RT: zgap = %g\n",rtin->zgap) ;
02188          VMCHECK ;
02189 
02190       } else if( STARTER("XYZOFF") ){             /* 18 Dec 2002 */
02191          float xval = 0.0 , yval = 0.0 , zval = 0.0 ;
02192          sscanf( buf , "XYZOFF %f %f %f" , &xval , &yval , &zval ) ;
02193          rtin->xxoff = xval ; rtin->xcen = 1 ;
02194          rtin->yyoff = yval ; rtin->ycen = 1 ;
02195          rtin->zzoff = zval ; rtin->zcen = 1 ;
02196          if( verbose == 2 )
02197             fprintf(stderr,"RT: offset = %g %g %g\n",rtin->xxoff,rtin->yyoff,rtin->zzoff) ;
02198          VMCHECK ;
02199 
02200       } else if( STARTER("ZFIRST") ){   /* set z origin */
02201          float val = 0.0 ;
02202          char dcode = ' ' ;
02203          sscanf( buf , "ZFIRST %f%c" , &val,&dcode ) ;
02204          rtin->zzorg   = val ;
02205          rtin->zcen    = 0 ;
02206          rtin->zzdcode = ORCODE(dcode) ;
02207          if( verbose == 2 )
02208             fprintf(stderr,"RT: zzorg = %g%c\n" ,
02209                     rtin->zzorg ,
02210                     (rtin->zzdcode < 0) ? ' ' : ORIENT_first[rtin->zzdcode] ) ;
02211          VMCHECK ;
02212 
02213       } else if( STARTER("XYZFIRST") ){  /* 10 Dec 2002: set all origins */
02214          float xf=0.0,yf=0.0,zf=0.0 ;
02215          char  xc=' ',yc=' ',zc=' ' ;
02216          sscanf( buf , "XYZFIRST %f%c%f%c%f%c" , &xf,&xc,&yf,&yc,&zf,&zc ) ;
02217          rtin->xxorg = xf ; rtin->xcen = 0 ; rtin->xxdcode = ORCODE(xc) ;
02218          rtin->yyorg = yf ; rtin->ycen = 0 ; rtin->yydcode = ORCODE(yc) ;
02219          rtin->zzorg = zf ; rtin->zcen = 0 ; rtin->zzdcode = ORCODE(zc) ;
02220          if( verbose == 2 )
02221             fprintf(stderr,"RT: xxorg=%g%c yyorg=%g%c zzorg=%g%c\n" ,
02222                     rtin->xxorg ,
02223                     (rtin->xxdcode < 0) ? ' ' : ORIENT_first[rtin->xxdcode] ,
02224                     rtin->yyorg ,
02225                     (rtin->yydcode < 0) ? ' ' : ORIENT_first[rtin->yydcode] ,
02226                     rtin->zzorg ,
02227                     (rtin->zzdcode < 0) ? ' ' : ORIENT_first[rtin->zzdcode]
02228                    ) ;
02229          VMCHECK ;
02230 
02231       } else if( STARTER("XYFOV") ){
02232          float xval = 0.0 , yval = 0.0 , zval = 0.0 ;
02233          sscanf( buf , "XYFOV %f %f %f" , &xval , &yval , &zval ) ;
02234          if( xval > 0.0 ){
02235             rtin->xxfov = xval ;
02236             rtin->yyfov = (yval > 0.0) ? yval : xval ;
02237             if( zval > 0.0 ) rtin->zzfov = zval ;
02238          } else
02239                 BADNEWS ;
02240          if( verbose == 2 )
02241             fprintf(stderr,"RT: fov = %g %g %g\n",rtin->xxfov,rtin->yyfov,rtin->zzfov) ;
02242          VMCHECK ;
02243 
02244       } else if( STARTER("XYMATRIX") ){
02245          int xval = 0 , yval = 0 , zval = 0 ;
02246          sscanf( buf , "XYMATRIX %d %d %d" , &xval , &yval , &zval ) ;
02247          if( xval > 1 ){
02248             rtin->nxx = xval ;
02249             rtin->nyy = (yval > 1) ? yval : xval ;
02250             if( zval > 0 ){
02251                rtin->nzz = zval ;
02252                if( rtin->nzz < 1 ) fprintf(stderr,"RT: # slices = %d!\a\n",zval) ;
02253             }
02254          } else
02255                 BADNEWS ;
02256          if( verbose == 2 )
02257             fprintf(stderr,"RT: matrix = %d %d %d\n",rtin->nxx,rtin->nyy,rtin->nzz) ;
02258          VMCHECK ;
02259 
02260       } else if( STARTER("ZNUM") ){
02261          int zval = 0 ;
02262          sscanf( buf , "ZNUM %d" , &zval ) ;
02263          if( zval > 0 ){
02264              rtin->nzz = zval ;
02265              if( rtin->nzz < 1 ) fprintf(stderr,"RT: # slices = %d!\a\n",zval) ;
02266          }
02267          else
02268               BADNEWS ;
02269          if( verbose == 2 && rtin->nzz >= 1 )
02270             fprintf(stderr,"RT: # slices = %d\n",rtin->nzz) ;
02271          VMCHECK ;
02272 
02273       } else if( STARTER("DATUM") ){
02274          int ii ;
02275          char tstr[32] = "\0" ;
02276          sscanf( buf , "DATUM %31s" , tstr ) ;
02277          for( ii=0 ; ii <= LAST_MRI_TYPE ; ii++ )
02278             if( strcmp(tstr,MRI_TYPE_name[ii]) == 0 ) break ;
02279 
02280          if( AFNI_GOOD_DTYPE(ii) ) rtin->datum = ii ;
02281          else
02282               BADNEWS ;
02283          if( verbose == 2 )
02284             fprintf(stderr,"RT: datum code = %d\n",rtin->datum) ;
02285          VMCHECK ;
02286 
02287       } else if( STARTER("BYTEORDER") ){    /* 27 Jun 2003:           [rickr] */
02288          int bo = 0 ;
02289          char tstr[10] = "\0" ;
02290          sscanf( buf, "BYTEORDER %9s", tstr ) ;
02291 
02292          /* first, note the incoming endian */
02293          if      ( strncmp(tstr,"LSB_FIRST",9) == 0 ) bo = LSB_FIRST ;
02294          else if ( strncmp(tstr,"MSB_FIRST",9) == 0 ) bo = MSB_FIRST ;
02295          else
02296              BADNEWS ;
02297 
02298          /* if different from the local endian, we will swap bytes */
02299          if ( bo != 0 ) {
02300             int local_bo, one = 1;
02301             local_bo = (*(char *)&one == 1) ? LSB_FIRST : MSB_FIRST ;
02302 
02303             /* if we are informed, and the orders differ, we will swap */
02304             if ( bo != local_bo )
02305                 rtin->swap_on_read = 1 ;
02306          }
02307 
02308          if( verbose > 1 )
02309             fprintf(stderr,"RT: BYTEORDER string = '%s', swap_on_read = %d\n",
02310                     BYTE_ORDER_STRING(bo), rtin->swap_on_read) ;
02311       } else if( STARTER("LOCK_ZORDER") ){  /* 22 Feb 1999:                   */
02312          rtin->zorder_lock = 1 ;            /* allow program to 'lock' zorder,*/
02313                                             /* so that later changes are      */
02314       } else if( STARTER("ZORDER") ){       /* ineffective.                   */
02315          if( ! rtin->zorder_lock ){
02316            char str[32] = "\0" ; int nord=0 , nb = 0 , nq ;
02317            sscanf( buf , "ZORDER %31s%n" , str , &nb ) ;
02318                 if( strcmp(str,"alt") == 0 ) rtin->zorder = ZORDER_ALT ;
02319            else if( strcmp(str,"seq") == 0 ) rtin->zorder = ZORDER_SEQ ;
02320            else if( strcmp(str,"explicit") == 0 ){
02321               rtin->zorder = ZORDER_EXP ;
02322               do{                                  /* get all numbers */
02323                  rtin->zseq[nord] = -1 ; nq = -1 ;
02324                  sscanf(buf+nb , "%d%n" , &(rtin->zseq[nord]) , &nq) ;
02325                  if( nq < 1 ) break ;              /* failed to get */
02326                  nb += nq ; nord++ ;               /* otherwise, increment */
02327               } while( nb < nbuf ) ;               /* until end of buffer */
02328               rtin->nzseq = nord ;                 /* number found */
02329               if( nord < 1 || nord > NZMAX ) BADNEWS ;
02330            }
02331            else
02332                 BADNEWS ;
02333          }
02334 
02335       } else if( STARTER("TPATTERN") ){               /* get timing pattern  */
02336          if( ! rtin->zorder_lock ){                   /* 10 May 2005 [rickr] */
02337            char str[32] = "\0" ;
02338            sscanf( buf , "TPATTERN %31s" , str ) ;
02339                 if( strcmp(str,"alt+z") == 0 ) rtin->tpattern = ZORDER_ALT ;
02340            else if( strcmp(str,"seq+z") == 0 ) rtin->tpattern = ZORDER_SEQ ;
02341            else
02342                 BADNEWS ;
02343          }
02344 
02345       } else if( STARTER("XYZAXES") ){
02346          int ii , orx=-1,ory=-1,orz=-1 ;
02347          char xstr[32] = "\0" , ystr[32] = "\0" , zstr[32] = "\0" ;
02348 
02349          sscanf( buf , "XYZAXES %31s %31s %31s" , xstr,ystr,zstr ) ;
02350 
02351          for( ii=0 ; ii < 6 ; ii++ ){
02352             if( strcmp(xstr,ORIENT_shortstr[ii]) == 0 ||
02353                 strcmp(xstr,ORIENT_typestr[ii])  == 0 ||
02354                 strcmp(xstr,ORIENT_tinystr[ii])  == 0   ) orx = ii ;
02355 
02356             if( strcmp(ystr,ORIENT_shortstr[ii]) == 0 ||
02357                 strcmp(ystr,ORIENT_typestr[ii])  == 0 ||
02358                 strcmp(ystr,ORIENT_tinystr[ii])  == 0   ) ory = ii ;
02359 
02360             if( strcmp(zstr,ORIENT_shortstr[ii]) == 0 ||
02361                 strcmp(zstr,ORIENT_typestr[ii])  == 0 ||
02362                 strcmp(zstr,ORIENT_tinystr[ii])  == 0   ) orz = ii ;
02363          }
02364 
02365          if( orx >= 0 && ory >= 0 && orz >= 0 && OR3OK(orx,ory,orz) ){
02366             rtin->orcxx = orx ;
02367             rtin->orcyy = ory ;
02368             rtin->orczz = orz ;
02369          } else
02370                 BADNEWS ;
02371 
02372       } else if( STARTER("NUM_CHAN") ){     /* 01 Aug 2002 */
02373          int nn=0 ;
02374          sscanf( buf , "NUM_CHAN %d",&nn) ;
02375          if( nn >= 1 && nn <= MAX_CHAN )
02376            rtin->num_chan = nn ;
02377          else
02378                 BADNEWS ;
02379 
02380       } else if( STARTER("DRIVE_AFNI") ){   /* 30 Jul 2002 */
02381          char cmd[1024]="\0" ;
02382          int ii ;
02383          if( strlen(buf) < 11 ){
02384             fprintf(stderr,"RT: DRIVE_AFNI lacks command\n") ;
02385          } else {  /* the command is everything after "DRIVE_AFNI " */
02386             MCW_strncpy(cmd,buf+11,1024) ;
02387             if( verbose == 2 )
02388                fprintf(stderr,"RT: command DRIVE_AFNI %s\n",cmd) ;
02389             ii = AFNI_driver( cmd ) ;  /* just do it */
02390             if( ii < 0 )
02391                fprintf(stderr,"RT: command DRIVE_AFNI %s **FAILS**\n",cmd) ;
02392          }
02393 
02394       } else {                              /* this is bad news */
02395          BADNEWS ;
02396       }
02397    }  /* end of loop over command buffers */
02398 
02399    /** now, determine if enough information exists to create a dataset **/
02400 
02401    if( rtin->image_mode ){
02402       rtin->nzz      = 1 ;  /* 28 Apr 2000 */
02403       rtin->num_chan = 1 ;  /* 01 Aug 2002 */
02404    }
02405 
02406    /** 01 Aug 2002: turn some things off in multi-channel mode **/
02407 
02408    if( rtin->num_chan > 1 ){
02409 
02410      if( rtin->reg_mode > 0 && verbose )
02411        fprintf(stderr,"RT: %d channel acquisition => no registration!\n",rtin->num_chan) ;
02412 
02413      if( rtin->func_code > 0 && verbose )
02414        fprintf(stderr,"RT: %d channel acquisition => no function!\n"    ,rtin->num_chan) ;
02415 
02416      rtin->reg_mode  = REGMODE_NONE ;  /* no registration */
02417      rtin->func_code = FUNC_NONE ;     /* no function */
02418      rtin->func_func = NULL ;
02419      rtin->reg_graph = 0 ;
02420    }
02421 
02422    if( rtin->nzz == 1 ){                   /* 24 Jun 2002: 1 slice only? */
02423      rtin->zorder = ZORDER_SEQ ;
02424      rtin->tpattern = ZORDER_SEQ ;
02425      if( REG_IS_3D(rtin->reg_mode) ){
02426        rtin->reg_mode = REGMODE_NONE ;
02427        fprintf(stderr,"RT: can't do 3D registration on 2D dataset!\n") ;
02428      }
02429    }
02430 
02431    RT_check_info( rtin , 0 ) ;
02432 
02433    /** if possible, now compute the number of bytes in each input image **/
02434 
02435    if( AFNI_GOOD_DTYPE(rtin->datum) && rtin->nxx > 0 && rtin->nyy > 0 ){
02436       int n1 = mri_datum_size( rtin->datum ) ;
02437 
02438       if( rtin->dtype == DTYPE_2DZT || rtin->dtype == DTYPE_2DZ )
02439          rtin->imsize = rtin->nxx * rtin->nyy * n1 ;
02440       else if( rtin->nzz > 0 )
02441          rtin->imsize = rtin->nxx * rtin->nyy * rtin->nzz * n1 ;
02442    }
02443 
02444    /* validate data type when swapping */
02445    if ( rtin->swap_on_read == 1 ) {
02446 
02447       if( (rtin->datum != MRI_short) &&         /* if the type is not okay, */
02448           (rtin->datum != MRI_int)   &&         /* then turn off swapping   */
02449           (rtin->datum != MRI_float) &&
02450           (rtin->datum != MRI_complex) )
02451       {
02452          if( rtin->datum != MRI_byte )          /* don't complain about bytes */
02453              fprintf(stderr,"RT: BYTEORDER applies only to short, int, float "
02454                             "or complex\n");
02455          rtin->swap_on_read = 0;
02456       }
02457    }
02458 
02459    /** return the number of characters processed **/
02460 
02461    return (nend+1) ;
02462 }
02463 
02464 /*--------------------------------------------------------------------
02465    Create the dataset and prepare to receive image data.  Some of this
02466    code is taken shamelessly from to3d.c (I wrote it, I can steal it).
02467    This routine should only be called when rtin->info_ok is true.
02468 ----------------------------------------------------------------------*/
02469 
02470 void RT_start_dataset( RT_input * rtin )
02471 {
02472    THD_ivec3 nxyz , orixyz ;
02473    THD_fvec3 dxyz , orgxyz ;
02474    int nvox , npix , n1 , ii , cc ;
02475    char npr[THD_MAX_PREFIX] , ccpr[THD_MAX_PREFIX] ;
02476 
02477    /*********************************************/
02478    /** 28 Apr 2000: Image Only mode is simpler **/
02479 
02480    if( rtin->image_mode ){
02481       nvox = rtin->nxx * rtin->nyy ;
02482       n1   = mri_datum_size( rtin->datum ) ;
02483 
02484       rtin->sbr_size = nvox * n1 ;                /* size of image space */
02485       rtin->sbr[0]   = malloc( rtin->sbr_size ) ; /* image space */
02486       rtin->imsize   = rtin->sbr_size ;
02487       rtin->im[0]    = rtin->sbr[0] ;             /* image space again */
02488       rtin->nzz      = 1 ;
02489 
02490       rtin->image_space = mri_new_vol_empty( rtin->nxx, rtin->nyy, 1, rtin->datum ) ;
02491       mri_fix_data_pointer( rtin->sbr[0] , rtin->image_space ) ;
02492 
02493       return ;
02494    }
02495 
02496    /*******************************************/
02497    /** Must create a dataset and other stuff **/
02498 
02499    /*******************************************************************/
02500    /** 02 Aug 2002: in multi-channel mode, need multiple controllers **/
02501 
02502    if( rtin->num_chan > 1 ){
02503      int ic , tc , nc=AFNI_count_controllers() , sn ;
02504 
02505      Three_D_View *im3d = plint->im3d ;             /* default controller */
02506      if( !IM3D_OPEN(im3d) )                         /* may not be open */
02507        im3d = AFNI_find_open_controller() ;
02508 
02509      sn = im3d->vinfo->sess_num ;                   /* session for all channels */
02510 
02511      /* open new controllers if needed */
02512 
02513      if( nc < rtin->num_chan && nc < MAX_CONTROLLERS ){
02514        nc = MIN( rtin->num_chan , MAX_CONTROLLERS ) - nc ;  /* number to add */
02515 
02516        fprintf(stderr,
02517                "RT: %d channel data requires opening %d more AFNI controllers",
02518                rtin->num_chan , nc ) ;
02519 
02520        for( ic=0 ; ic < nc ; ic++ )                         /* add them */
02521          AFNI_clone_controller_CB(NULL,NULL,NULL) ;
02522      }
02523 
02524      /* make list of open controllers */
02525 
02526      num_open_controllers = AFNI_count_controllers() ;
02527      for( ic=nc=0 ; ic < MAX_CONTROLLERS ; ic++ ){
02528        if( IM3D_OPEN(GLOBAL_library.controllers[ic]) ){
02529          open_controller[nc]       = GLOBAL_library.controllers[ic] ;
02530          open_controller_index[nc] = ic ;
02531          nc++ ;
02532        }
02533      }
02534      nc = num_open_controllers ;  /* nugatory or moot */
02535 
02536      /* now assign channels to controllers */
02537 
02538      tc = MIN( nc , rtin->num_chan ) ;
02539      for( cc=0 ; cc < rtin->num_chan ; cc++ ){
02540        if( cc < tc ) rtin->im3d[cc] = open_controller[cc] ; /* next available */
02541        else          rtin->im3d[cc] = NULL ;                /* none available */
02542 
02543        /* all channels go to the 1st session directory */
02544 
02545        rtin->sess_num[cc] = sn ;
02546        rtin->sess[cc]     = GLOBAL_library.sslist->ssar[sn] ;
02547      }
02548 
02549    /* single channel mode: make sure calling controller is still open */
02550 
02551    } else {
02552      Three_D_View *im3d = plint->im3d ;             /* default controller */
02553      if( !IM3D_OPEN(im3d) )                         /* may not be open */
02554        im3d = AFNI_find_open_controller() ;
02555 
02556      rtin->im3d[0]     = im3d ;
02557      rtin->sess_num[0] = rtin->im3d[0]->vinfo->sess_num ;
02558      rtin->sess[0]     = GLOBAL_library.sslist->ssar[rtin->sess_num[0]] ;
02559    }
02560 
02561    /***************************/
02562    /** Let there be dataset! **/  /* 01 Aug 2002: or datasets */
02563 
02564    for( cc=0 ; cc < rtin->num_chan ; cc++ ){
02565      rtin->dset[cc] = EDIT_empty_copy(NULL) ;
02566      tross_Append_History( rtin->dset[cc] , "plug_realtime: creation" ) ;
02567 
02568      if( rtin->num_note > 0 && rtin->note != NULL ){  /* 01 Oct 2002 */
02569        for( ii=0 ; ii < rtin->num_note ; ii++ )
02570          tross_Add_Note( rtin->dset[cc] , rtin->note[ii] ) ;
02571      }
02572    }
02573 
02574    /********************************/
02575    /** make a good dataset prefix **/
02576 
02577    for( ii=1 ; ; ii++ ){
02578       sprintf( npr , "%.*s#%03d" , RT_MAX_PREFIX, rtin->root_prefix , ii ) ;
02579 
02580       if( rtin->num_chan == 1 ){                  /* the old way */
02581 
02582          if( PLUTO_prefix_ok(npr) ) break ;
02583 
02584       } else {                                    /* 02 Aug 2002: the multichannel way */
02585 
02586         for( cc=0 ; cc < rtin->num_chan ; cc++ ){ /* check each channel prefix */
02587           sprintf(ccpr, "%s_%02d", npr,cc+1 ) ;
02588           if( !PLUTO_prefix_ok(ccpr) ) break ;    /* not good? quit cc loop */
02589         }
02590         if( cc == rtin->num_chan ) break ;        /* all were good? we are done */
02591       }
02592    }
02593 
02594    /**********************************************************/
02595    /** correct the spatial information for axes orientation **/
02596 
02597    rtin->dxx = rtin->xxfov / rtin->nxx ;
02598    rtin->dyy = rtin->yyfov / rtin->nyy ;
02599 
02600    /* 18 Dec 2002: add zgap here, per UCSD request */
02601 
02602    if( rtin->zzfov > 0 ) rtin->dzz = rtin->zzfov / rtin->nzz + rtin->zgap ;
02603 
02604    /* 18 Dec 2002: add xxoff (etc.) here, per UCSD request */
02605 
02606    if( rtin->xcen ) rtin->xxorg = 0.5 * (rtin->nxx - 1) * rtin->dxx + rtin->xxoff ;
02607    if( rtin->ycen ) rtin->yyorg = 0.5 * (rtin->nyy - 1) * rtin->dyy + rtin->yyoff ;
02608    if( rtin->zcen ) rtin->zzorg = 0.5 * (rtin->nzz - 1) * rtin->dzz + rtin->zzoff ;
02609 
02610    /* if axis direction is a 'minus', voxel increments are negative */
02611 
02612    if( ORIENT_sign[rtin->orcxx] == '-' ) rtin->dxx = - rtin->dxx ;
02613    if( ORIENT_sign[rtin->orcyy] == '-' ) rtin->dyy = - rtin->dyy ;
02614    if( ORIENT_sign[rtin->orczz] == '-' ) rtin->dzz = - rtin->dzz ;
02615 
02616    /* 10 Dec 2002:
02617       We now allow direction codes input for x and y as well,
02618       so must duplicate the complicated zzorg logic for setting
02619       the xxorg and yyorg signs.  This duplicated logic follows
02620       the zzdcode stuff, below.                                 */
02621 
02622 #if 0
02623    /* if axis direction is a 'plus',
02624       then the origin is in a 'minus' direction,
02625       so the origin offset must have its sign changed */  /*** THE OLD WAY ***/
02626 
02627    if( ORIENT_sign[rtin->orcxx] == '+' ) rtin->xxorg = - rtin->xxorg ;
02628    if( ORIENT_sign[rtin->orcyy] == '+' ) rtin->yyorg = - rtin->yyorg ;
02629 #endif
02630 
02631    /* z-origin is complex, because the remote program might
02632       have given a direction code, or it might not have.
02633       We must set the sign of the z-origin based on the
02634       sign of the axis direction if no specific
02635       direction code is given, otherwise it should be
02636       based on the specific direction code (rtin->zzdcode). */
02637 
02638    if( rtin->zzdcode < 0 ){  /* no specific code */
02639 
02640       if( ORIENT_sign[rtin->orczz] == '+' ) rtin->zzorg = - rtin->zzorg ;
02641 
02642    } else {  /* have code */
02643 
02644       /* check if code matches the axis direction
02645          (or is the opposite, which is also OK).  */
02646 
02647       if( rtin->orczz != rtin->zzdcode                 &&
02648           rtin->orczz != ORIENT_OPPOSITE(rtin->zzdcode)  ){  /* this is bad */
02649 
02650          fprintf(stderr,"RT: ZFIRST direction code = %c but Z axis = %s!\a\n",
02651                  ORIENT_first[rtin->zzdcode] , ORIENT_shortstr[rtin->orczz] ) ;
02652 
02653          if( ORIENT_sign[rtin->orczz] == '+' ) rtin->zzorg = - rtin->zzorg ;
02654 
02655       } else {  /* things are OK */
02656 
02657          if( ORIENT_sign[rtin->zzdcode] == '+' ) rtin->zzorg = - rtin->zzorg ;
02658       }
02659    }
02660 
02661    /* 10 Dec 2002: duplicate logic for xxorg */
02662 
02663    if( rtin->xxdcode < 0 ){
02664       if( ORIENT_sign[rtin->orcxx] == '+' ) rtin->xxorg = - rtin->xxorg ;
02665    } else {
02666       if( rtin->orcxx != rtin->xxdcode                 &&
02667           rtin->orcxx != ORIENT_OPPOSITE(rtin->xxdcode)  ){
02668          fprintf(stderr,"RT: XFIRST direction code = %c but X axis = %s!\a\n",
02669                  ORIENT_first[rtin->xxdcode] , ORIENT_shortstr[rtin->orcxx] ) ;
02670          if( ORIENT_sign[rtin->orcxx] == '+' ) rtin->xxorg = - rtin->xxorg ;
02671       } else {  /* things are OK */
02672          if( ORIENT_sign[rtin->xxdcode] == '+' ) rtin->xxorg = - rtin->xxorg ;
02673       }
02674    }
02675 
02676    /* 10 Dec 2002: duplicate logic for yyorg */
02677 
02678    if( rtin->yydcode < 0 ){
02679       if( ORIENT_sign[rtin->orcyy] == '+' ) rtin->yyorg = - rtin->yyorg ;
02680    } else {
02681       if( rtin->orcyy != rtin->yydcode                 &&
02682           rtin->orcyy != ORIENT_OPPOSITE(rtin->yydcode)  ){
02683          fprintf(stderr,"RT: YFIRST direction code = %c but Y axis = %s!\a\n",
02684                  ORIENT_first[rtin->yydcode] , ORIENT_shortstr[rtin->orcyy] ) ;
02685          if( ORIENT_sign[rtin->orcyy] == '+' ) rtin->yyorg = - rtin->yyorg ;
02686       } else {  /* things are OK */
02687          if( ORIENT_sign[rtin->yydcode] == '+' ) rtin->yyorg = - rtin->yyorg ;
02688       }
02689    }
02690 
02691    /************************************************/
02692    /** add the spatial information to the dataset **/  /* 01 Aug 2002: datasets */
02693 
02694    nxyz.ijk[0]   = rtin->nxx   ; dxyz.xyz[0]   = rtin->dxx ;
02695    nxyz.ijk[1]   = rtin->nyy   ; dxyz.xyz[1]   = rtin->dyy ;
02696    nxyz.ijk[2]   = rtin->nzz   ; dxyz.xyz[2]   = rtin->dzz ;
02697 
02698    orixyz.ijk[0] = rtin->orcxx ; orgxyz.xyz[0] = rtin->xxorg ;
02699    orixyz.ijk[1] = rtin->orcyy ; orgxyz.xyz[1] = rtin->yyorg ;
02700    orixyz.ijk[2] = rtin->orczz ; orgxyz.xyz[2] = rtin->zzorg ;
02701 
02702    for( cc=0 ; cc < rtin->num_chan ; cc++ ){  /* 01 Aug 2002: loop over datasets */
02703 
02704      if( rtin->num_chan == 1 )
02705        strcpy( ccpr , npr ) ;                      /* 1 channel prefix */
02706      else
02707        sprintf( ccpr , "%s_%02d" , npr , cc+1 ) ;  /* new prefix for multi-channel */
02708 
02709      EDIT_dset_items( rtin->dset[cc] ,
02710                          ADN_prefix      , ccpr ,
02711                          ADN_datum_all   , rtin->datum ,
02712                          ADN_nxyz        , nxyz ,
02713                          ADN_xyzdel      , dxyz ,
02714                          ADN_xyzorg      , orgxyz ,
02715                          ADN_xyzorient   , orixyz ,
02716                          ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
02717                          ADN_nvals       , 1 ,
02718                          ADN_type        , HEAD_ANAT_TYPE ,
02719                          ADN_view_type   , VIEW_ORIGINAL_TYPE ,
02720                          ADN_func_type   , ANAT_EPI_TYPE ,
02721                       ADN_none ) ;
02722 
02723       /******************************************/
02724       /** add time axis information, if needed **/
02725 
02726       if( rtin->dtype == DTYPE_3DT ||                      /* all slices   */
02727          (rtin->dtype == DTYPE_2DZT && rtin->nzz == 1) ){  /* simultaneous */
02728 
02729          EDIT_dset_items( rtin->dset[cc] ,
02730                              ADN_ntt      , 1 ,
02731                              ADN_ttorg    , 0.0 ,
02732                              ADN_ttdel    , rtin->tr ,
02733                              ADN_ttdur    , 0.0 ,
02734                              ADN_tunits   , UNITS_SEC_TYPE ,
02735                           ADN_none ) ;
02736 
02737       } else if( rtin->dtype == DTYPE_2DZT ){  /* slices at different times */
02738 
02739          float * tpattern  = (float *) malloc( sizeof(float) * rtin->nzz ) ;
02740          float   tframe    = rtin->tr / rtin->nzz ;
02741          float   tsl ;
02742          int     ii ;
02743 
02744          if( rtin->zorder == ZORDER_ALT || rtin->tpattern == ZORDER_ALT ){
02745             /* alternating +z direction */
02746             tsl = 0.0 ;
02747             for( ii=0 ; ii < rtin->nzz ; ii+=2 ){
02748                tpattern[ii] = tsl ; tsl += tframe ;
02749             }
02750             for( ii=1 ; ii < rtin->nzz ; ii+=2 ){
02751                tpattern[ii] = tsl ; tsl += tframe ;
02752             }
02753          }
02754          /* make sequential a terminating case     10 May 2005 [rickr] */
02755          /* else if( rtin->zorder == ZORDER_SEQ ){                     */
02756          /* ... more tpatterns may be added above here                 */
02757          else {
02758             tsl = 0.0 ;
02759             for( ii=0 ; ii < rtin->nzz ; ii++ ){
02760                tpattern[ii] = tsl ; tsl += tframe ;
02761             }
02762          }
02763 
02764          EDIT_dset_items( rtin->dset[cc] ,
02765                              ADN_ntt      , 1 ,
02766                              ADN_ttorg    , 0.0 ,
02767                              ADN_ttdel    , rtin->tr ,
02768                              ADN_ttdur    , 0.0 ,
02769                              ADN_tunits   , UNITS_SEC_TYPE ,
02770                              ADN_nsl      , rtin->nzz ,
02771                              ADN_zorg_sl  , rtin->zzorg ,
02772                              ADN_dz_sl    , rtin->dzz ,
02773                              ADN_toff_sl  , tpattern ,
02774                           ADN_none ) ;
02775 
02776          free( tpattern ) ;
02777       }
02778 
02779       rtin->afni_status[cc] = 0 ;  /* uninformed at this time */
02780       DSET_lock(rtin->dset[cc]) ;  /* 20 Mar 1998 */
02781    }
02782 
02783 #ifdef ALLOW_REGISTRATION
02784    /*---- Make a dataset for registration, if need be ----*/
02785 
02786    if( REG_MAKE_DSET(rtin->reg_mode) &&
02787        ((rtin->dtype==DTYPE_2DZT) || (rtin->dtype==DTYPE_3DT)) ){
02788 
02789       rtin->reg_dset = EDIT_empty_copy( rtin->dset[0] ) ;
02790       tross_Append_History( rtin->reg_dset , "plug_realtime: registration" ) ;
02791 
02792       strcpy(ccpr,npr) ;
02793            if( REG_IS_2D(rtin->reg_mode) ) strcat(ccpr,"%reg2D") ;
02794       else if( REG_IS_3D(rtin->reg_mode) ) strcat(ccpr,"%reg3D") ;
02795       else                                 strcat(ccpr,"%reg"  ) ;
02796 
02797       EDIT_dset_items( rtin->reg_dset , ADN_prefix , ccpr , ADN_none ) ;
02798       DSET_lock(rtin->reg_dset) ;
02799 
02800       if( rtin->num_note > 0 && rtin->note != NULL ){  /* 01 Oct 2002 */
02801         for( ii=0 ; ii < rtin->num_note ; ii++ )
02802           tross_Add_Note( rtin->reg_dset , rtin->note[ii] ) ;
02803       }
02804    }
02805 
02806    rtin->reg_status    = 0 ;
02807    rtin->reg_nvol      = 0 ;
02808 
02809    /* No longer scale to time units as x-axis is now in terms of reps.
02810     *                                      *** 29 Jan 2004 [rickr] ***
02811     *  rtin->reg_graph_xr *= rtin->tr ;    *** scale to time units ***
02812     */
02813 
02814 #endif
02815 
02816    /***********************************************/
02817    /** now prepare space for incoming image data **/
02818 
02819    nvox = rtin->nxx * rtin->nyy * rtin->nzz ;
02820    n1   = mri_datum_size( rtin->datum ) ;
02821 
02822    rtin->sbr_size = nvox * n1 ;                /* size of sub-brick */
02823    for( cc=0 ; cc < rtin->num_chan ; cc++ ){
02824      rtin->sbr[cc] = malloc( rtin->sbr_size ) ; /* space to hold sub-brick */
02825      if( rtin->sbr[cc] == NULL ){
02826        fprintf(stderr,
02827                "RT: can't malloc data space for real-time channel %02d!\a\n",
02828                cc+1) ;
02829        EXIT(1) ;
02830      }
02831    }
02832 
02833    if( rtin->dtype == DTYPE_2DZT || rtin->dtype == DTYPE_2DZ )
02834       rtin->imsize = rtin->nxx * rtin->nyy * n1 ;
02835    else
02836       rtin->imsize = nvox * n1 ;
02837 
02838    for( cc=0 ; cc < rtin->num_chan ; cc++ ){
02839      rtin->im[cc]   = rtin->sbr[cc] ;  /* place to put first image in sub-brick */
02840      rtin->nsl[cc]  = 0 ;              /* number of slices gathered so far */
02841      rtin->nvol[cc] = 0 ;              /* number of volumes gathered so far */
02842    }
02843 
02844    /** if there is already data stored in the temporary buffer
02845        (acquired before this routine was called), then we want
02846        to place it into the dataset now.                       **/
02847 
02848    if( rtin->bufar != NULL ){
02849       int ii , upsave ;
02850       MRI_IMAGE * ibb ;
02851       char * bbb ;
02852 
02853       if( verbose == 2 )
02854          fprintf(stderr,"RT: putting %d buffered images into dataset\n" ,
02855                         IMARR_COUNT(rtin->bufar) ) ;
02856       VMCHECK ;
02857 
02858       upsave = update ;  /* don't do updates to AFNI during this step */
02859       update = 0 ;
02860 
02861       for( ii=0 ; ii < IMARR_COUNT(rtin->bufar) ; ii++ ){
02862         ibb = IMARR_SUBIMAGE(rtin->bufar,ii) ;                    /* next buffered image */
02863         bbb = (char *) MRI_BYTE_PTR( ibb ) ;                      /* pointer to its data */
02864         memcpy( rtin->im[rtin->cur_chan] , bbb , rtin->imsize ) ; /* copy into sub-brick */
02865         mri_free( ibb ) ;                                         /* free buffered image */
02866         RT_process_image( rtin ) ;                                /* send into dataset   */
02867       }
02868       FREE_IMARR( rtin->bufar ) ;   /* throw this away like a used Kleenex */
02869 
02870       update = upsave ;
02871 
02872       if( verbose == 2 )
02873          fprintf(stderr,"RT: buffered images all placed into dataset\n") ;
02874       VMCHECK ;
02875    }
02876 
02877    /** dataset(s) now created and ready to boogie! **/
02878    /** popup a message to keep the user occupied.  **/
02879 
02880    { char str[1024] ;
02881      char acq[128] , *sli ;
02882 
02883      switch( rtin->dtype ){
02884        case DTYPE_2DZ:  strcpy(acq,"2D+z (1 volume, by slice")        ; break ;
02885        case DTYPE_2DZT: strcpy(acq,"2D+zt (multivolume, by slice")    ; break ;
02886        case DTYPE_3D:   strcpy(acq,"3D (1 volume, all at once)")      ; break ;
02887        case DTYPE_3DT:  strcpy(acq,"3D+t (multivolume, by 3D array)") ; break ;
02888        default:         strcpy(acq,"Bizarro world")                   ; break ;
02889      }
02890 
02891      if( rtin->dtype == DTYPE_2DZ || rtin->dtype == DTYPE_2DZT ){
02892        switch( rtin->zorder ){
02893          case ZORDER_ALT: strcat(acq," - interleaved order)") ; break ;
02894          case ZORDER_SEQ: strcat(acq," - sequential order)")  ; break ;
02895          case ZORDER_EXP: strcat(acq," - explicit order)")    ; break ;
02896        }
02897      }
02898 
02899      switch( rtin->orczz ){
02900        case ORI_I2S_TYPE:
02901        case ORI_S2I_TYPE: sli = "(Axial)"    ; break ;
02902 
02903        case ORI_R2L_TYPE:
02904        case ORI_L2R_TYPE: sli = "(Sagittal)" ; break ;
02905 
02906        case ORI_P2A_TYPE:
02907        case ORI_A2P_TYPE: sli = "(Coronal)"  ; break ;
02908 
02909        default:           sli = "\0"         ; break ;  /* say what? */
02910      }
02911 
02912      sprintf(str," \n"
02913                  " ** Realtime Header Information **\n"
02914                  "\n"
02915                  " Dataset prefix  : %s\n"
02916                  " Brick Dimensions: %d x %d x %d\n"
02917                  " Voxel Grid Size : %.4f x %.4f x %.4f (mm)\n"
02918                  " Grid Orientation: %s x %s x %s %s\n"
02919                  " Grid Offset     : %.1f x %.1f x %.1f (mm)\n"
02920                  " Datum Type      : %s\n"
02921                  " Number Channels : %d\n"
02922                  " Acquisition Type: %s\n" ,
02923              npr ,
02924              rtin->nxx , rtin->nyy , rtin->nzz ,
02925              fabs(rtin->dxx) , fabs(rtin->dyy) , fabs(rtin->dzz) ,
02926              ORIENT_shortstr[rtin->orcxx],ORIENT_shortstr[rtin->orcyy],
02927                ORIENT_shortstr[rtin->orczz] , sli ,
02928              rtin->xxorg , rtin->yyorg , rtin->zzorg ,
02929              MRI_TYPE_name[rtin->datum] ,
02930              rtin->num_chan ,
02931              acq
02932           ) ;
02933 
02934      PLUTO_popup_transient(plint,str);
02935    }
02936 
02937    return ;
02938 }
02939 
02940 /*--------------------------------------------------------------------
02941   Read one image (or volume) into the space pointed to by im.
02942 ----------------------------------------------------------------------*/
02943 
02944 #define COMMAND_MARKER        "Et Earello Endorenna utulien!!"
02945 #define COMMAND_MARKER_LENGTH 30
02946 
02947 void RT_read_image( RT_input * rtin , char * im )
02948 {
02949    int need , have , nbuffed ;
02950 
02951    /** sanity checks **/
02952 
02953    if( rtin == NULL || im == NULL ){
02954      fprintf(stderr,"RT: illegal inputs to RT_read_image!\a\n") ;
02955      EXIT(1) ;
02956    }
02957 
02958    if( rtin->imsize <= 0 ){
02959      fprintf(stderr,"RT: image data present, but don't know its size!\a\n") ;
02960      EXIT(1) ;
02961    }
02962 
02963    /** see if any data in buffer already **/
02964 
02965    have = rtin->nbuf ;
02966 
02967    /** if we have data already, use it first **/
02968 
02969    if( have > 0 ){
02970 
02971       nbuffed = MIN( have , rtin->imsize ) ;  /* how much to read from buffer */
02972       memcpy( im , rtin->buf , nbuffed ) ;    /* copy it into image */
02973 
02974       if( nbuffed < have ){                   /* didn't use up all of buffer */
02975          memmove( rtin->buf , rtin->buf + nbuffed , rtin->nbuf - nbuffed ) ;
02976          rtin->nbuf = rtin->nbuf - nbuffed ;
02977       } else {                                /* used up all of buffer */
02978          rtin->nbuf = 0 ;
02979       }
02980    } else {
02981       nbuffed = 0 ;
02982    }
02983 
02984    /** at this point, nbuffed is the number of bytes read from the buffer.
02985        if we need to read more bytes from the I/O channel, then do so now. **/
02986 
02987    need = rtin->imsize - nbuffed ;  /* number of bytes left to fill image */
02988 
02989    /* read this many bytes for sure (waiting if need be) */
02990 
02991    if( need > 0 )
02992       iochan_recvall( rtin->ioc_data , im + nbuffed , need ) ;
02993 
02994    /* 10 Dec 2002:
02995       Check if the command string is present at the start of im.
02996       If it is, then mark rtin for death.                       */
02997 
02998    if( memcmp(im,COMMAND_MARKER,COMMAND_MARKER_LENGTH) == 0 )
02999      rtin->marked_for_death = 1 ;
03000    else {
03001       /* we have a complete image, check for byte swapping */
03002       if ( rtin->swap_on_read != 0 ) {
03003          if( rtin->datum == MRI_short )
03004             mri_swap2( rtin->imsize / 2, (short *)im );
03005          else
03006             mri_swap4( rtin->imsize / 4, (int *)im );
03007       }
03008    }
03009 
03010 
03011    return ;
03012 }
03013 
03014 /*--------------------------------------------------------------------
03015    Read images and put them into place.  Note that if there is any
03016    data left in the rtin buffer, it will go into the image first.
03017 ----------------------------------------------------------------------*/
03018 
03019 int RT_process_data( RT_input * rtin )
03020 {
03021    int vdone ;
03022 
03023    /** can we create a dataset yet? **/
03024 
03025    if( rtin->sbr[0] == NULL && rtin->info_ok ){
03026       if( verbose == 2 )
03027          fprintf(stderr,"RT: info complete --> creating dataset.\n") ;
03028       VMCHECK ;
03029       RT_start_dataset( rtin ) ;
03030    }
03031 
03032    /** read images as long as there is data to suck up **/
03033 
03034    while( rtin->nbuf > 0 || iochan_readcheck(rtin->ioc_data,0) > 0 ){
03035 
03036       if( rtin->im[0] != NULL ){  /** process data into dataset directly **/
03037 
03038          RT_read_image( rtin , rtin->im[rtin->cur_chan] ) ; /* read into dataset buffer */
03039 
03040          if( rtin->marked_for_death ) return 0 ;            /* 10 Dec 2002 */
03041 
03042          RT_process_image( rtin ) ;                         /* process it for the dataset */
03043 
03044       } else {                 /** read data into temporary buffer space **/
03045                                /** (will have to process it later, dude) **/
03046          MRI_IMAGE * newim ;
03047          char * newbuf ;
03048 
03049          if( rtin->imsize <= 0 ){
03050            fprintf(stderr,"RT: image data present, but don't know its size!\a\n") ;
03051            EXIT(1) ;
03052          }
03053 
03054          if( rtin->bufar == NULL )    /* initialize buffer for input images */
03055            INIT_IMARR(rtin->bufar) ;
03056 
03057          if( verbose == 2 && rtin->bufar->num % 10 == 0 ){
03058            fprintf(stderr,"RT: reading image into buffer[%d]\n",rtin->bufar->num) ;
03059            VMCHECK ;
03060          }
03061 
03062          newim  = mri_new( rtin->imsize , 1 , MRI_byte ) ; /* make space for next image */
03063          newbuf = (char *) MRI_BYTE_PTR(newim) ;           /* pointer to image data */
03064          ADDTO_IMARR( rtin->bufar , newim ) ;              /* add to saved image list */
03065          RT_read_image( rtin , newbuf ) ;                  /* read image data */
03066          if( rtin->marked_for_death ) return 0 ;           /* 10 Dec 2002 */
03067       }
03068 
03069       RT_process_xevents( rtinp ) ;
03070 
03071    }  /* end of loop over reading images as long as we can */
03072 
03073    /** return an error code if the input data channel has gone bad **/
03074 
03075 /**
03076    if( iochan_goodcheck(rtin->ioc_data,0) <= 0 ) return -1 ;
03077 **/
03078 
03079    return 1 ;
03080 }
03081 
03082 /*--------------------------------------------------------------------
03083    Function to be called if the user kills the Image Only viewer
03084 ----------------------------------------------------------------------*/
03085 
03086 static void RT_image_kfun(void * kdata)                /* 28 Apr 2000 */
03087 {
03088    RT_input * rtin = (RT_input *) kdata ;
03089    if( rtin != NULL ) rtin->image_handle = NULL ;
03090    return ;
03091 }
03092 
03093 /*-------------------------------------------------------------------
03094    Given image data stored in rtin->im, process it into the dataset
03095    that is under construction.  This routine should only be called
03096    after RT_start_dataset has been executed.
03097 ---------------------------------------------------------------------*/
03098 
03099 #define MIN_TO_GRAPH 2
03100 
03101 void RT_process_image( RT_input * rtin )
03102 {
03103    int vdone , cc = rtin->cur_chan ;
03104 
03105    /** 28 Apr 2000: Image Only mode stuff **/
03106 
03107    if( rtin->image_mode ){
03108       if( rtin->image_handle != NULL ){
03109          PLUTO_imseq_addto( rtin->image_handle , rtin->image_space ) ;
03110       } else {
03111          rtin->image_handle = PLUTO_imseq_popim( rtin->image_space, RT_image_kfun,rtin ) ;
03112          PLUTO_imseq_retitle( rtin->image_handle , "Realtime Images" ) ;
03113       }
03114       return ;
03115    }
03116 
03117    /** check if new image completes the current volume **/
03118 
03119    if( rtin->dtype == DTYPE_2DZT || rtin->dtype == DTYPE_2DZ ){
03120 
03121       if( verbose == 2 )
03122          fprintf(stderr,"RT: read image into dataset brick %d slice %d.\n",
03123                  rtin->nvol[cc],rtin->nsl[cc]) ;
03124 
03125       rtin->nsl[cc] ++ ;                     /* 1 more slice */
03126       vdone = (rtin->nsl[cc] == rtin->nzz) ; /* have all slices? */
03127 
03128    } else if( rtin->dtype == DTYPE_3DT || rtin->dtype == DTYPE_3D ){
03129 
03130 #if 0
03131       if( verbose == 2 )
03132         fprintf(stderr,"RT: read image into dataset brick %d\n",rtin->nvol[cc]) ;
03133 #endif
03134 
03135       vdone = 1 ;                        /* 3D gets all data at once */
03136    }
03137 
03138    /** if volume is done, add it to the dataset **/
03139 
03140    if( vdone ){
03141 
03142       rtin->nvol[cc] ++ ;        /* 1 more volume is acquired! */
03143 
03144       if( verbose == 2 )
03145          fprintf(stderr,"RT: now have %d complete sub-bricks in channel %02d.\n",
03146                  rtin->nvol[cc],cc+1) ;
03147       VMCHECK ;
03148 
03149       /* first time: put this volume in as "substitute" for empty 1st brick
03150          later:      add new volume at end of chain                         */
03151 
03152       if( rtin->nvol[cc] == 1 )
03153          EDIT_substitute_brick( rtin->dset[cc] , 0 , rtin->datum , rtin->sbr[cc] ) ;
03154       else
03155          EDIT_add_brick( rtin->dset[cc] , rtin->datum , 0.0 , rtin->sbr[cc] ) ;
03156 
03157       VMCHECK ;
03158       if( verbose == 2 )
03159          fprintf(stderr,"RT: added brick to dataset in channel %02d\n",cc+1) ;
03160       VMCHECK ;
03161 
03162       /* must also change the number of times recorded
03163          [EDIT_add_brick does 'nvals' correctly, but not 'ntt'] */
03164 
03165       if( rtin->dtype == DTYPE_3DT || rtin->dtype == DTYPE_2DZT ){
03166          EDIT_dset_items( rtin->dset[cc] , ADN_ntt , rtin->nvol[cc] , ADN_none ) ;
03167          if( verbose == 2 )
03168             fprintf(stderr,"RT: altered ntt in dataset header in channel %02d\n",cc+1) ;
03169          VMCHECK ;
03170       } else if( rtin->nvol[cc] > 1 ){
03171          fprintf(stderr,"RT: have %d bricks for time-independent dataset!\a\n",
03172                  rtin->nvol[cc]) ;
03173          VMCHECK ;
03174       }
03175 
03176 #ifdef ALLOW_REGISTRATION
03177       /* do registration before function computation   - 30 Oct 2003 [rickr] */
03178       switch( rtin->reg_mode ){
03179            case REGMODE_2D_RTIME: RT_registration_2D_realtime( rtin ) ;
03180            break ;
03181 
03182            case REGMODE_3D_RTIME:
03183            case REGMODE_3D_ESTIM: RT_registration_3D_realtime( rtin ) ;
03184            break ;
03185       }
03186 #endif /* ALLOW_REGISTRATION */
03187 
03188       /** compute function, maybe? **/
03189 
03190       if( rtin->func_code > 0 ){
03191          int jj ;
03192 
03193          /** if needed, initialize the function computations **/
03194 
03195          if( rtin->func_condit == 0 ){
03196 #if 0
03197             jj = rtin->func_func( rtin , INIT_MODE ) ;
03198 #else
03199             AFNI_CALL_VALU_2ARG( rtin->func_func , int,jj ,
03200                                  RT_input *,rtin , int,INIT_MODE ) ;
03201 #endif
03202             if( jj < 0 ){ rtin->func_code = 0 ; rtin->func_func = NULL ; }
03203             rtin->func_condit = 1 ;  /* initialized */
03204          }
03205 
03206          /** do the function computations for this volume **/
03207 
03208          if( rtin->func_code > 0 )
03209 #if 0
03210             jj = rtin->func_func( rtin , rtin->nvol[cc] - 1 ) ;
03211 #else
03212             AFNI_CALL_VALU_2ARG( rtin->func_func , int,jj ,
03213                                  RT_input *,rtin , int,rtin->nvol[cc]-1 ) ;
03214 #endif
03215       }
03216 
03217       /** make space for next sub-brick to arrive **/
03218 
03219       if( verbose == 2 )
03220          fprintf(stderr,"RT: malloc-ing %d bytes for next volume in channel %02d\n",
03221                  rtin->sbr_size,cc+1) ;
03222       VMCHECK ;
03223 
03224       rtin->sbr[cc] = malloc( rtin->sbr_size ) ;
03225       if( rtin->sbr[cc] == NULL ){
03226         fprintf(stderr,"RT: can't malloc real-time brick %d for channel %02d\a\n",
03227                 rtin->nvol[cc]+1,cc+1) ;
03228         EXIT(1) ;
03229       }
03230       if( verbose == 2 )
03231          fprintf(stderr,"RT: malloc succeeded\n") ;
03232       VMCHECK ;
03233 
03234       rtin->im[cc]  = rtin->sbr[cc] ;  /* location of slice #0 within sub-brick */
03235       rtin->nsl[cc] = 0 ;              /* number of slices gathered so far */
03236 
03237       /** tell AFNI about this dataset, maybe **/
03238 
03239       if( update > 0 ){  /* if we want updates every so often */
03240          int doit ;
03241 
03242          if( verbose == 2 )
03243             fprintf(stderr,"RT: checking for update status\n") ;
03244          VMCHECK ;
03245 
03246          doit = ( (rtin->dtype==DTYPE_3DT || rtin->dtype==DTYPE_2DZT) &&
03247                   (cc+1 == rtin->num_chan)                            && /* 01 Aug 2002 */
03248                   (rtin->nvol[cc] == MIN_TO_GRAPH ||
03249                    (rtin->nvol[cc] > MIN_TO_GRAPH && rtin->nvol[cc] % update == 0)) ) ;
03250 
03251          if( doit ){
03252             if( verbose == 2 )
03253                fprintf(stderr,"RT: about to tell AFNI about dataset.\n") ;
03254             VMCHECK ;
03255             RT_tell_afni(rtin,TELL_NORMAL) ;
03256          }
03257       }
03258 
03259    } else {  /** need to add more slices before volume is done **/
03260 
03261       int noff ;  /* slice position in output brick */
03262 
03263       if( rtin->zorder == ZORDER_SEQ ){  /* sequential:       */
03264                                          /* slice position is */
03265          noff = rtin->nsl[cc] ;          /* just slice number */
03266 
03267       } else if ( rtin->zorder == ZORDER_ALT ){  /* alternating: */
03268 
03269          int nhalf = (rtin->nzz + 1)/2 ;         /* number in first 1/2 */
03270 
03271          if( rtin->nsl[cc] < nhalf )
03272             noff = 2 * rtin->nsl[cc] ;               /* first half of slices */
03273          else
03274             noff = 1 + 2 * (rtin->nsl[cc] - nhalf) ; /* second half of slices */
03275       }
03276 
03277       rtin->im[cc] = rtin->sbr[cc] + (noff * rtin->imsize) ; /* where next slice goes */
03278 
03279    }
03280 
03281    /* 01 Aug 2002: update which channel gets the next image */
03282 
03283    if( rtin->num_chan > 1 )
03284      rtin->cur_chan = (cc+1) % rtin->num_chan ;
03285 
03286    return ;
03287 }
03288 
03289 /*-------------------------------------------------------------------
03290    Tell AFNI about all datasets we are building.
03291 ---------------------------------------------------------------------*/
03292 
03293 void RT_tell_afni( RT_input *rtin , int mode )
03294 {
03295    int cc ;
03296 
03297    if( rtin == NULL ) return ;
03298 
03299    /*** tell separately for each one ***/
03300 
03301    for( cc=0 ; cc < rtin->num_chan ; cc++ )
03302      RT_tell_afni_one( rtin , mode , cc ) ;
03303 
03304    /*** at the end of acquisition, show messages ***/
03305 
03306    if( mode == TELL_FINAL && ISVALID_DSET(rtin->dset[0]) ){
03307      char qbuf[256*(MAX_CHAN+1)] , zbuf[256] ;
03308      sprintf( qbuf ,
03309                 " \n"
03310                 " Acquisition Terminated\n\n"
03311                 " Brick Dimensions: %d x %d x %d  Datum: %s\n\n" ,
03312                rtin->nxx,rtin->nyy,rtin->nzz , MRI_TYPE_name[rtin->datum] );
03313 
03314      for( cc=0 ; cc < rtin->num_chan ; cc++ ){
03315        if( ISVALID_DSET(rtin->dset[cc]) )
03316          sprintf( zbuf ,
03317                   " Channel %02d: dataset %s has %d sub-bricks\n",
03318                   cc+1 , DSET_FILECODE(rtin->dset[cc]) , rtin->nvol[cc] ) ;
03319        else
03320          sprintf( zbuf ,
03321                   " Channel %d: INVALID DATASET?!\n",cc) ;
03322        strcat( qbuf , zbuf ) ;
03323      }
03324      strcat(qbuf,"\n") ;
03325 
03326      PLUTO_beep(); PLUTO_popup_transient(plint,qbuf);
03327 
03328      if( verbose == 2 ) SHOW_TIMES ;
03329 
03330      sync() ;  /* 08 Mar 2000: sync disk */
03331    }
03332 
03333    return ;
03334 }
03335 
03336 /*-------------------------------------------------------------------
03337    Tell AFNI about one dataset we are building.
03338 ---------------------------------------------------------------------*/
03339 
03340 void RT_tell_afni_one( RT_input *rtin , int mode , int cc )
03341 {
03342    Three_D_View * im3d ;
03343    THD_session * sess ;
03344    THD_3dim_dataset * old_anat ;
03345    int ii , id , afni_init=0 ;
03346    char clll , *cstr ;
03347    MCW_arrowval * tav ;
03348 
03349    /** sanity check **/
03350 
03351    if( rtin == NULL || !ISVALID_DSET(rtin->dset[cc]) ) return ;
03352 
03353    im3d = rtin->im3d[cc] ;                      /* AFNI controller */
03354    sess = rtin->sess[cc] ;                      /* AFNI session */
03355 
03356    if( im3d != NULL ){   /* 05 Aug 2002: allow for NULL im3d controller */
03357 
03358      tav  = im3d->vwid->imag->time_index_av ;   /* time index widget */
03359      ii   = AFNI_controller_index( im3d ) ;
03360      cstr = AFNI_controller_label( im3d ) ; clll = cstr[1] ;
03361 
03362      old_anat = im3d->anat_now ;  /* the default */
03363    }
03364 
03365    /**--- deal with the input dataset ---**/
03366 
03367    if( rtin->afni_status[cc] == 0 ){  /** first time for this dataset **/
03368 
03369       EDIT_dset_items( rtin->dset[cc],
03370                          ADN_directory_name,sess->sessname,
03371                        ADN_none ) ;
03372 
03373       if( im3d != NULL ){
03374         if( verbose )
03375            fprintf(stderr , "RT: sending dataset %s with %d bricks\n"
03376                             "    to AFNI[%c], session %s\n" ,
03377                    DSET_FILECODE(rtin->dset[cc]) , rtin->nvol[cc] ,
03378                    clll , sess->sessname ) ;
03379 
03380       }
03381       THD_load_statistics( rtin->dset[cc] ) ;
03382 
03383       /** put it into the current session in the current controller **/
03384 
03385       if( GLOBAL_library.have_dummy_dataset ) UNDUMMYIZE ;
03386 
03387       id = sess->num_dsset ;
03388 
03389       if( id >= THD_MAX_SESSION_SIZE ){
03390         fprintf(stderr,"RT: max number of anat datasets exceeded!\a\n") ;
03391         EXIT(1) ;
03392       }
03393       sess->dsset[id][VIEW_ORIGINAL_TYPE] = rtin->dset[cc] ;
03394       sess->num_dsset = id+1 ;
03395       POPDOWN_strlist_chooser ;
03396 
03397       if( ISFUNC(rtin->dset[cc]) )
03398         AFNI_force_adoption( sess , False ) ;
03399 
03400       /** tell AFNI controller to jump to this dataset and session **/
03401 
03402       if( im3d != NULL ){
03403         if( ISANAT(rtin->dset[cc]) )
03404            im3d->vinfo->anat_num = sess->num_dsset - 1 ;
03405         else
03406            im3d->vinfo->func_num = sess->num_dsset - 1 ;
03407 
03408         im3d->vinfo->sess_num = rtin->sess_num[cc] ;
03409 
03410         im3d->vinfo->time_index = 0 ;
03411         AV_assign_ival( tav , 0 ) ;
03412       }
03413 
03414       afni_init             = 1 ; /* below: will initialize AFNI to see this dataset */
03415       rtin->afni_status[cc] = 1 ; /* mark dataset to show that AFNI knows about it now */
03416 
03417    } else {  /** 2nd or later call for this dataset **/
03418 
03419       THD_update_statistics( rtin->dset[cc] ) ;
03420 
03421       if( im3d != NULL ){
03422         if( mode != TELL_FINAL && verbose )
03423           fprintf(stderr,"RT: update with %d bricks in channel %02d to AFNI[%c]\n",
03424                   rtin->nvol[cc],cc+1,clll) ;
03425 
03426         tav->fmax = tav->imax = im3d->vinfo->time_index = rtin->nvol[cc] - 1  ;
03427         AV_assign_ival( tav , tav->imax ) ;
03428       }
03429 
03430    }
03431 
03432    /**--- Deal with the computed function, if any ---**/
03433 
03434    if( rtin->func_dset != NULL ){
03435       int jj ;
03436 
03437 #if 0
03438       jj = rtin->func_func( rtin , UPDATE_MODE ) ;  /* put data into dataset */
03439 #else
03440       AFNI_CALL_VALU_2ARG( rtin->func_func , int,jj ,
03441                            RT_input *,rtin , int,UPDATE_MODE ) ;
03442 #endif
03443 
03444       if( rtin->func_condit > 1 ){             /* data actually inside */
03445 
03446          THD_load_statistics( rtin->func_dset ) ; /* statistickizification */
03447 
03448          if( rtin->func_status == 0 ){  /** first time for this dataset **/
03449 
03450             if( im3d != NULL && verbose )
03451                fprintf(stderr , "RT: sending dataset %s with %d bricks\n"
03452                                 "    to AFNI controller [%c] session %s\n" ,
03453                        DSET_FILECODE(rtin->func_dset) , DSET_NVALS(rtin->func_dset) ,
03454                        clll , sess->sessname ) ;
03455 
03456             EDIT_dset_items( rtin->func_dset, ADN_directory_name,sess->sessname, ADN_none ) ;
03457             id = sess->num_dsset ;
03458             if( id >= THD_MAX_SESSION_SIZE ){
03459               fprintf(stderr,"RT: max number of datasets exceeded!\a\n") ;
03460               EXIT(1) ;
03461             }
03462             sess->dsset[id][VIEW_ORIGINAL_TYPE] = rtin->func_dset ; sess->num_dsset++ ;
03463             AFNI_force_adoption( sess , False ) ;
03464             POPDOWN_strlist_chooser ;
03465 
03466             if( im3d != NULL ){
03467               im3d->vinfo->func_num = sess->num_dsset - 1 ;
03468               AFNI_SETUP_FUNC_ON(im3d) ;
03469             }
03470 
03471             rtin->func_status = 1 ;   /* AFNI knows now */
03472             afni_init = 1 ;           /* below: tell AFNI to look at this function */
03473 
03474          } else {  /** 2nd or later call for this dataset **/
03475 
03476             /**--- actually, there's nothing to do here ---**/
03477 
03478          }
03479       }  /* end of if functional dataset actually has data */
03480    }  /* end of if functional dataset exists */
03481 
03482 #ifdef ALLOW_REGISTRATION
03483    /**--- Deal with the registered dataset, if any ---**/
03484 
03485    if( rtin->reg_dset != NULL && rtin->reg_nvol > 0 ){
03486 
03487       if( rtin->reg_status == 0 ){  /** first time for this dataset **/
03488 
03489          THD_load_statistics( rtin->reg_dset ) ;
03490 
03491          if( im3d != NULL ){
03492            if( verbose )
03493                  fprintf(stderr , "RT: sending dataset %s with %d bricks\n"
03494                                   "    to AFNI controller [%c] session %s\n" ,
03495                          DSET_FILECODE(rtin->reg_dset) , DSET_NVALS(rtin->reg_dset) ,
03496                          clll , sess->sessname ) ;
03497          }
03498 
03499          EDIT_dset_items( rtin->reg_dset, ADN_directory_name,sess->sessname, ADN_none ) ;
03500 
03501          id = sess->num_dsset ;
03502          if( id >= THD_MAX_SESSION_SIZE ){
03503            fprintf(stderr,"RT: max number of datasets exceeded!\a\n") ;
03504            EXIT(1) ;
03505          }
03506          sess->dsset[id][VIEW_ORIGINAL_TYPE] = rtin->reg_dset ; sess->num_dsset = id+1 ;
03507          POPDOWN_strlist_chooser ;
03508 
03509          rtin->reg_status = 1 ;   /* AFNI knows about this dataset now */
03510 
03511       } else {  /** 2nd or later call for this dataset **/
03512 
03513          THD_update_statistics( rtin->reg_dset ) ;
03514       }
03515    }
03516 #endif
03517 
03518    /**--- actually talk to AFNI now ---**/
03519 
03520    if( afni_init ){  /* tell AFNI to view new dataset(s) */
03521 
03522      if( im3d != NULL ){
03523        AFNI_SETUP_VIEW(im3d,VIEW_ORIGINAL_TYPE) ;
03524        if( EQUIV_DSETS(rtin->dset[cc],old_anat) ) THD_update_statistics( rtin->dset[cc] ) ;
03525        else                                       THD_load_statistics  ( rtin->dset[cc] ) ;
03526 
03527        AFNI_initialize_view( old_anat , im3d ) ;  /* Geronimo! */
03528      }
03529 
03530    } else {          /* just tell AFNI to refresh the images/graphs */
03531 
03532       Three_D_View * qq3d ;
03533       int review ;
03534 
03535       /* check all controllers to see if they are looking at the same datasets */
03536 
03537       for( ii=0 ; ii < MAX_CONTROLLERS ; ii++ ){
03538          qq3d = GLOBAL_library.controllers[ii] ;
03539          if( !IM3D_OPEN(qq3d) ) break ;  /* skip this one */
03540 
03541          review = (qq3d == im3d)                             ||   /* same viewer?  */
03542                   EQUIV_DSETS(rtin->dset[cc],qq3d->anat_now) ||   /* or same anat? */
03543                   ( rtin->func_dset != NULL   &&                  /* or same func? */
03544                     qq3d->vinfo->func_visible &&
03545                     EQUIV_DSETS(rtin->func_dset,qq3d->fim_now) ) ;
03546 
03547 #ifdef ALLOW_REGISTRATION
03548          review = review || ( rtin->reg_dset != NULL &&       /* or same reg?  */
03549                               rtin->reg_nvol > 0     &&
03550                               EQUIV_DSETS(rtin->reg_dset,qq3d->anat_now) ) ;
03551 #endif
03552 
03553          if( review ){
03554             AFNI_modify_viewing( qq3d , False ) ;  /* Crazy Horse! */
03555          }
03556       }
03557    }
03558 
03559    if( im3d != NULL ) XmUpdateDisplay( THE_TOPSHELL ) ;
03560 
03561    /**--- if this is the final call, do some cleanup stuff ---**/
03562 
03563    if( mode == TELL_FINAL ){
03564 
03565       int cmode ;
03566 
03567       if( verbose == 2 )
03568          fprintf(stderr,"RT: finalizing dataset to AFNI (including disk output).\n") ;
03569 
03570 #if 0
03571       if( im3d != NULL )
03572         fprintf(stderr , "RT: sending dataset %s with %d bricks\n"
03573                          "    to AFNI controller [%c] session %s\n" ,
03574                 DSET_FILECODE(rtin->dset[cc]) , rtin->nvol[cc] ,
03575                 clll , sess->sessname ) ;
03576 #endif
03577 
03578 #if 0
03579       THD_load_statistics( rtin->dset[cc] ) ;
03580 #endif
03581 
03582       /* 20 Mar 1998: write in uncompressed mode for speed */
03583 
03584       cmode = THD_get_write_compression() ;
03585       THD_set_write_compression(COMPRESS_NONE) ;
03586       SHOW_AFNI_PAUSE ;
03587 
03588       THD_write_3dim_dataset( NULL,NULL , rtin->dset[cc] , True ) ;
03589       DSET_unlock( rtin->dset[cc] ) ;  /* 20 Mar 1998 */
03590 
03591       if( rtin->func_dset != NULL ){
03592          int jj ;
03593 #if 0
03594          jj = rtin->func_func( rtin , FINAL_MODE ) ;
03595 #else
03596          AFNI_CALL_VALU_2ARG( rtin->func_func , int,jj ,
03597                               RT_input *,rtin , int,FINAL_MODE ) ;
03598 #endif
03599          THD_write_3dim_dataset( NULL,NULL , rtin->func_dset , True ) ;
03600          DSET_unlock( rtin->func_dset ) ;  /* 20 Mar 1998 */
03601          THD_force_malloc_type( rtin->func_dset->dblk , DATABLOCK_MEM_ANY ) ;
03602       }
03603 
03604 #ifdef ALLOW_REGISTRATION
03605       if( rtin->reg_dset != NULL && rtin->reg_nvol > 0 ){
03606          THD_write_3dim_dataset( NULL,NULL , rtin->reg_dset , True ) ;
03607          DSET_unlock( rtin->reg_dset ) ;
03608          THD_force_malloc_type( rtin->reg_dset->dblk , DATABLOCK_MEM_ANY ) ;
03609       }
03610 #endif
03611 
03612       THD_set_write_compression(cmode) ;  /* restore compression mode */
03613 
03614       AFNI_force_adoption( sess , GLOBAL_argopt.warp_4D ) ;
03615       AFNI_make_descendants( GLOBAL_library.sslist ) ;
03616       THD_force_malloc_type( rtin->dset[cc]->dblk , DATABLOCK_MEM_ANY ) ;
03617 
03618       AFNI_purge_unused_dsets() ;
03619       SHOW_AFNI_READY ;
03620    }
03621 
03622    return ;
03623 }
03624 
03625 /*--------------------------------------------------------------------------
03626   This routine is called when no more data will be added to the
03627   incoming dataset(s).
03628 ----------------------------------------------------------------------------*/
03629 
03630 void RT_finish_dataset( RT_input * rtin )
03631 {
03632    int ii , cc , nbad=0 ;
03633 
03634    if( rtin->image_mode ){
03635       if( verbose == 2 ) SHOW_TIMES ;
03636       return ;
03637    }
03638 
03639    for( cc=0 ; cc < rtin->num_chan ; cc++ ){
03640 
03641       if( ! ISVALID_3DIM_DATASET(rtin->dset[cc]) ){
03642         fprintf(stderr,"RT: attempt to finish channel %02d with incomplete dataset!\a\n",cc+1) ;
03643         nbad++ ; continue ;
03644       }
03645 
03646       if( rtin->nvol[cc] < 1 ){
03647         fprintf(stderr,"RT: attempt to finish channel %02d with 0 completed bricks!\a\n",cc+1) ;
03648         DSET_delete( rtin->dset[cc] ) ; rtin->dset[cc] = NULL ;
03649         if( rtin->func_dset != NULL ){
03650            DSET_delete( rtin->func_dset ) ; rtin->func_dset = NULL ;
03651         }
03652 #ifdef ALLOW_REGISTRATION
03653         if( rtin->reg_dset != NULL ){
03654            DSET_delete( rtin->reg_dset ) ; rtin->reg_dset = NULL ;
03655         }
03656 #endif
03657         nbad++ ;
03658       }
03659 
03660       if( rtin->nsl[cc] > 0 )
03661          fprintf(stderr,"RT: finish channel %02d with %d slices unused!\a\n",
03662                  cc+1,rtin->nsl[cc]);
03663 
03664       fprintf(stderr,"RT: finish channel %02d with %d bricks completed.\n",
03665               cc+1,rtin->nvol[cc]) ;
03666    }
03667 
03668    if( verbose ) SHOW_TIMES ;
03669    if( nbad    ) return ;
03670 
03671 #ifdef ALLOW_REGISTRATION
03672    /*--- Do the "at end" registration, if ordered and if possible ---*/
03673 
03674    switch( rtin->reg_mode ){
03675       case REGMODE_2D_ATEND: RT_registration_2D_atend( rtin ) ; break ;
03676       case REGMODE_3D_ATEND: RT_registration_3D_atend( rtin ) ; break ;
03677    }
03678 
03679    /*--- 17 Aug 1998: Graph the motion parameters, if desired ---*/
03680 
03681    if( rtin->reg_graph && rtin->reg_nest > 1 && REG_IS_2D(rtin->reg_mode) ){
03682       float * yar[4] ;  /* not Tasha */
03683       int * iar , ii , nn = rtin->reg_nest ;
03684       static char * nar[3] = { "\\Delta x [mm]" , "\\Delta y [mm]" , "\\phi   [\\degree]" } ;
03685 
03686       if( verbose == 2 )
03687          fprintf(stderr,"RT: graphing estimated 2D motion parameters\n") ;
03688 
03689       /* sort the arrays by time */
03690 
03691       iar    = (int *)   malloc( sizeof(int)   * nn ) ;  /* index */
03692       yar[0] = (float *) malloc( sizeof(float) * nn ) ;  /* time  */
03693       yar[1] = (float *) malloc( sizeof(float) * nn ) ;  /* dx    */
03694       yar[2] = (float *) malloc( sizeof(float) * nn ) ;  /* dy    */
03695       yar[3] = (float *) malloc( sizeof(float) * nn ) ;  /* phi   */
03696 
03697       for( ii=0 ; ii < nn ; ii++ ){
03698          iar[ii] = ii ; yar[0][ii] = rtin->reg_tim[ii] ;
03699       }
03700 
03701       qsort_floatint( nn , yar[0] , iar ) ;     /* sort by time, carry index */
03702 
03703       for( ii=0 ; ii < nn ; ii++ ){             /* use index to reorder params */
03704          yar[1][ii] = rtin->reg_dx [ iar[ii] ] ;
03705          yar[2][ii] = rtin->reg_dy [ iar[ii] ] ;
03706          yar[3][ii] = rtin->reg_phi[ iar[ii] ] ;
03707       }
03708 
03709       plot_ts_lab( THE_DISPLAY ,
03710                    nn , yar[0] , -3 , yar+1 ,
03711                    "time" , NULL , DSET_FILECODE(rtin->dset[0]) , nar , NULL ) ;
03712 
03713       free(iar) ; free(yar[0]) ; free(yar[1]) ; free(yar[2]) ; free(yar[3]) ;
03714    }
03715 
03716    /*--- Oct 1998: do 3D graphing ---*/
03717    /*--- Jan 2004: if both reg_graph_xnew and reg_graph_ynew, don't plot */
03718 
03719    if( rtin->reg_graph && rtin->reg_nest > 1 && REG_IS_3D(rtin->reg_mode) &&
03720        ( rtin->reg_graph_xnew == 0 || rtin->reg_graph_ynew == 0 ) ){
03721       float * yar[7] ;
03722       int     ycount = -6 ;
03723       int     nn = rtin->reg_nest ;
03724       static char * nar[6] = {
03725          "\\Delta I-S [mm]" , "\\Delta R-L [mm]" , "\\Delta A-P [mm]" ,
03726          "Roll [\\degree]" , "Pitch [\\degree]" , "Yaw [\\degree]"  } ;
03727 
03728       char *ttl = malloc( strlen(DSET_FILECODE(rtin->dset[0])) + 32 ) ;
03729       strcpy(ttl,"\\noesc ") ;
03730       strcat(ttl,DSET_FILECODE(rtin->dset[0])) ;
03731       if( rtin->reg_mode == REGMODE_3D_ESTIM ) strcat(ttl," [Estimate]") ;
03732 
03733       if( verbose == 2 )
03734          fprintf(stderr,"RT: graphing estimated 3D motion parameters\n") ;
03735 
03736       /* arrays are already sorted by time */
03737 
03738       yar[0] = rtin->reg_rep   ;  /* repetition numbers */
03739       yar[1] = rtin->reg_dx    ;  /* dx    */
03740       yar[2] = rtin->reg_dy    ;  /* dy    */
03741       yar[3] = rtin->reg_dz    ;  /* dz    */
03742       yar[4] = rtin->reg_phi   ;  /* roll  */
03743       yar[5] = rtin->reg_psi   ;  /* pitch */
03744       yar[6] = rtin->reg_theta ;  /* yaw   */
03745 
03746       if ( rtin->p_code )
03747       {
03748          ycount = 1;
03749          yar[1] = rtin->reg_eval;
03750       }
03751 
03752       plot_ts_lab( THE_DISPLAY ,
03753                    nn , yar[0] , ycount , yar+1 ,
03754                    "reps" , NULL , ttl , nar , NULL ) ;
03755 
03756       free(ttl) ;
03757    }
03758 
03759    /* close any open tcp connection */
03760    if ( rtin->mp_tcp_use )
03761       RT_mp_comm_close( rtin );
03762 
03763    /* if we have a parser expression, free it */
03764    if ( rtin->p_code )
03765    {
03766       free( rtin->p_code ) ;
03767       rtin->p_code = NULL ;
03768    }
03769 #endif
03770 
03771    /** tell afni about it one last time **/
03772 
03773    RT_tell_afni(rtin,TELL_FINAL) ;
03774 
03775    return ;
03776 }
03777 
03778 #ifdef ALLOW_REGISTRATION
03779 /*****************************************************************************
03780   The collection of functions for handling slice registration in this plugin.
03781 ******************************************************************************/
03782 
03783 /*---------------------------------------------------------------------------
03784   Set the sizes of any open graph windows.
03785 
03786   Moved Bob's code to its own function.                  2004 Jan 28  [rickr]
03787 -----------------------------------------------------------------------------*/
03788 
03789 void RT_set_grapher_pinnums( int pinnum )
03790 {
03791     /* 12 Oct 2000: set pin_num on all graphs now open */
03792 
03793     if( pinnum < MIN_PIN || pinnum > MAX_PIN || !IM3D_OPEN(plint->im3d) )
03794         return;
03795 
03796     drive_MCW_grapher( plint->im3d->g123, graDR_setpinnum, (XtPointer) pinnum );
03797     drive_MCW_grapher( plint->im3d->g231, graDR_setpinnum, (XtPointer) pinnum );
03798     drive_MCW_grapher( plint->im3d->g312, graDR_setpinnum, (XtPointer) pinnum );
03799 }
03800 
03801 /*---------------------------------------------------------------------------
03802   Do pieces of 2D registration during realtime.
03803 -----------------------------------------------------------------------------*/
03804 
03805 void RT_registration_2D_realtime( RT_input * rtin )
03806 {
03807    int tt , ntt ;
03808 
03809    if( rtin->reg_dset == NULL ) return ;
03810 
03811    /*-- check to see if we need to setup first --*/
03812 
03813    if( rtin->reg_2dbasis == NULL ){  /* need to setup */
03814 
03815       /* check if enough data to setup */
03816 
03817       if( rtin->reg_base_index >= rtin->nvol[0] ) return ;  /* can't setup */
03818 
03819       /* setup the registration process */
03820 
03821       if( verbose )
03822          fprintf(stderr,"RT: setting up 2D registration 'realtime'\n") ;
03823 
03824       SHOW_AFNI_PAUSE ;
03825       RT_registration_2D_setup( rtin ) ;  /* may take a little while */
03826 
03827       if( rtin->reg_2dbasis == NULL ){
03828          fprintf(stderr,"RT: can't setup %s registration!\a\n",
03829                  REG_strings[REGMODE_2D_RTIME] ) ;
03830          DSET_delete( rtin->reg_dset ) ; rtin->reg_dset = NULL ;
03831          rtin->reg_mode = REGMODE_NONE ;
03832          SHOW_AFNI_READY ; return ;
03833       }
03834    }
03835 
03836    /*-- register all sub-bricks that aren't done yet --*/
03837 
03838    ntt = DSET_NUM_TIMES( rtin->dset[0] ) ;
03839    for( tt=rtin->reg_nvol ; tt < ntt ; tt++ )
03840       RT_registration_2D_onevol( rtin , tt ) ;
03841 
03842    /*-- my work here is done --*/
03843 
03844    XmUpdateDisplay( THE_TOPSHELL ) ;
03845    SHOW_AFNI_READY ; return ;
03846 }
03847 
03848 /*---------------------------------------------------------------------------
03849   Do 2D registration of all slices at once, when the realtime dataset
03850   is all present and accounted for.
03851 -----------------------------------------------------------------------------*/
03852 
03853 void RT_registration_2D_atend( RT_input * rtin )
03854 {
03855    int tt , ntt ;
03856 
03857    /* check if have enough data to register as ordered */
03858 
03859    if( rtin->reg_base_index >= rtin->nvol[0] ){
03860       fprintf(stderr,"RT: can't do %s registration: not enough 3D volumes!\a\n",
03861               REG_strings[REGMODE_2D_ATEND] ) ;
03862       DSET_delete( rtin->reg_dset ) ; rtin->reg_dset = NULL ;
03863       rtin->reg_mode = REGMODE_NONE ;
03864       return ;
03865    }
03866 
03867    /* set up the registration process */
03868 
03869    if( verbose )
03870       fprintf(stderr,"RT: starting 2D registration 'at end'\n") ;
03871 
03872    SHOW_AFNI_PAUSE ;
03873    RT_registration_2D_setup( rtin ) ;
03874 
03875    if( rtin->reg_2dbasis == NULL ){
03876       fprintf(stderr,"RT: can't setup %s registration!\a\n",
03877               REG_strings[REGMODE_2D_ATEND] ) ;
03878       DSET_delete( rtin->reg_dset ) ; rtin->reg_dset = NULL ;
03879       rtin->reg_mode = REGMODE_NONE ;
03880       SHOW_AFNI_READY ; return ;
03881    }
03882 
03883    /* register each volume into the new dataset */
03884 
03885    ntt = DSET_NUM_TIMES( rtin->dset[0] ) ;
03886    for( tt=0 ; tt < ntt ; tt++ ){
03887       XmUpdateDisplay( THE_TOPSHELL ) ;
03888       RT_registration_2D_onevol( rtin , tt ) ;
03889       if( verbose == 1 ) fprintf(stderr,"%d",tt%10) ;
03890    }
03891    if( verbose == 1 ) fprintf(stderr,"\n") ;
03892 
03893    /* un-setup the registration process */
03894 
03895    RT_registration_2D_close( rtin ) ;
03896    if( verbose ) SHOW_TIMES ;
03897    SHOW_AFNI_READY ; return ;
03898 }
03899 
03900 /*-----------------------------------------------------------------------
03901    Setup the 2D registration data for each slice
03902 -------------------------------------------------------------------------*/
03903 
03904 void RT_registration_2D_setup( RT_input * rtin )
03905 {
03906    int ibase = rtin->reg_base_index ;
03907    int kk , nx,ny,nz , kind , nbar ;
03908    MRI_IMAGE * im ;
03909    char * bar ;
03910 
03911    nx   = DSET_NX( rtin->dset[0] ) ;
03912    ny   = DSET_NY( rtin->dset[0] ) ;
03913    nz   = DSET_NZ( rtin->dset[0] ) ;
03914    kind = DSET_BRICK_TYPE( rtin->dset[0] , ibase ) ;
03915 
03916    rtin->reg_nvol  = 0 ;
03917 
03918    rtin->reg_2dbasis = (MRI_2dalign_basis **)
03919                          malloc( sizeof(MRI_2dalign_basis *) * nz ) ;
03920 
03921    im   = mri_new_vol_empty( nx,ny,1 , kind ) ;        /* fake image for slices */
03922    bar  = DSET_BRICK_ARRAY( rtin->dset[0] , ibase ) ;  /* ptr to base volume    */
03923    nbar = im->nvox * im->pixel_size ;               /* offset for each slice */
03924 
03925    for( kk=0 ; kk < nz ; kk++ ){                         /* loop over slices */
03926       mri_fix_data_pointer( bar + kk*nbar , im ) ;             /* base image */
03927       rtin->reg_2dbasis[kk] = mri_2dalign_setup( im , NULL ) ;
03928    }
03929 
03930    kk = rtin->reg_resam ;
03931    if( kk == MRI_QUINTIC || kk == MRI_HEPTIC )
03932       kk = MRI_BICUBIC ; /* quintic & heptic not available in 2D */
03933 
03934    mri_2dalign_method( MRI_BILINEAR, MRI_BICUBIC, kk ) ;
03935 
03936    mri_fix_data_pointer( NULL , im ) ; mri_free( im ) ;   /* get rid of this */
03937    return ;
03938 }
03939 
03940 /*---------------------------------------------------------------------------
03941     Undo the routine just above!
03942 -----------------------------------------------------------------------------*/
03943 
03944 void RT_registration_2D_close( RT_input * rtin )
03945 {
03946    int kk , nz ;
03947 
03948    nz = DSET_NZ( rtin->dset[0] ) ;
03949    for( kk=0 ; kk < nz ; kk++ )
03950       mri_2dalign_cleanup( rtin->reg_2dbasis[kk] ) ;
03951 
03952    free( rtin->reg_2dbasis ) ; rtin->reg_2dbasis = NULL ;
03953    return ;
03954 }
03955 
03956 /*-----------------------------------------------------------------------
03957    Carry out the 2D registration for each slice in the tt-th sub-brick
03958 -------------------------------------------------------------------------*/
03959 
03960 #ifndef D2RFAC
03961 #  define D2RFAC (PI/180.0)
03962 #  define R2DFAC (180.0/PI)
03963 #endif
03964 
03965 void RT_registration_2D_onevol( RT_input * rtin , int tt )
03966 {
03967    int kk , nx,ny,nz , kind , nbar ;
03968    MRI_IMAGE * im , * rim , * qim ;
03969    char * bar , * rar , * qar ;
03970    float dx,dy,phi ;        /* 17 Aug 1998 */
03971    int   nest , nxy ;
03972 
03973    /*-- sanity check --*/
03974 
03975    if( rtin->dset[0] == NULL || rtin->reg_dset == NULL ) return ;
03976 
03977    nx   = DSET_NX( rtin->dset[0] ) ;
03978    ny   = DSET_NY( rtin->dset[0] ) ; nxy = nx * ny ;
03979    nz   = DSET_NZ( rtin->dset[0] ) ;
03980    kind = DSET_BRICK_TYPE( rtin->dset[0] , 0 ) ;
03981 
03982    im   = mri_new_vol_empty( nx,ny,1 , kind ) ;     /* fake image for slices */
03983    bar  = DSET_BRICK_ARRAY( rtin->dset[0] , tt ) ;  /* ptr to input volume   */
03984    nbar = im->nvox * im->pixel_size ;               /* offset for each slice */
03985 
03986    /* make space for new sub-brick in reg_dset */
03987 
03988    if( verbose == 2 )
03989       fprintf(stderr,"RT: 2D registering sub-brick %d",tt) ;
03990 
03991    rar = (char *) malloc( sizeof(char) * nx*ny*nz * im->pixel_size ) ;
03992 
03993    if( rar == NULL ){
03994       fprintf(stderr,"RT: can't malloc space for registered dataset!\a\n") ;
03995       DSET_delete( rtin->reg_dset ) ; rtin->reg_dset = NULL ;
03996       rtin->reg_mode = REGMODE_NONE ;
03997       return ;
03998    }
03999 
04000    /*-- loop over slices --*/
04001 
04002    for( kk=0 ; kk < nz ; kk++ ){
04003 
04004       if( verbose == 2 ) fprintf(stderr,".") ;
04005 
04006       mri_fix_data_pointer( bar + kk*nbar , im ) ;  /* image to register */
04007 
04008       /* registration! */
04009 
04010       rim = mri_2dalign_one( rtin->reg_2dbasis[kk] , im , &dx , &dy , &phi ) ;
04011 
04012       /* 17 Aug 1998: save estimated motion parameters */
04013 
04014       nest = rtin->reg_nest ;
04015       rtin->reg_tim = (float *) realloc( (void *) rtin->reg_tim ,
04016                                          sizeof(float) * (nest+1) ) ;
04017       rtin->reg_dx  = (float *) realloc( (void *) rtin->reg_dx  ,
04018                                          sizeof(float) * (nest+1) ) ;
04019       rtin->reg_dy  = (float *) realloc( (void *) rtin->reg_dy  ,
04020                                          sizeof(float) * (nest+1) ) ;
04021       rtin->reg_phi = (float *) realloc( (void *) rtin->reg_phi ,
04022                                          sizeof(float) * (nest+1) ) ;
04023 
04024       rtin->reg_tim[nest] = THD_timeof_vox( tt , kk*nxy , rtin->dset[0] ) ;
04025       rtin->reg_dx [nest] = dx * DSET_DX(rtin->dset[0]) ;
04026       rtin->reg_dy [nest] = dy * DSET_DY(rtin->dset[0]) ;
04027       rtin->reg_phi[nest] = phi * R2DFAC             ; rtin->reg_nest ++ ;
04028 
04029       /* convert output image to desired type;
04030          set qar to point to data in the converted registered image */
04031 
04032       switch( kind ){
04033          case MRI_float:                        /* rim is already floats */
04034             qar = (char *) MRI_FLOAT_PTR(rim) ;
04035          break ;
04036 
04037          case MRI_short:
04038             qim = mri_to_short(1.0,rim) ; mri_free(rim) ; rim = qim ;
04039             qar = (char *) MRI_SHORT_PTR(rim) ;
04040          break ;
04041 
04042          case MRI_byte:
04043             qim = mri_to_byte(rim) ; mri_free(rim) ; rim = qim ;
04044             qar = (char *) MRI_BYTE_PTR(rim) ;
04045          break ;
04046 
04047          default:
04048             fprintf(stderr,"RT: can't do registration on %s images!\a\n",
04049                     MRI_TYPE_name[kind] ) ;
04050             DSET_delete( rtin->reg_dset ) ; rtin->reg_dset = NULL ;
04051             rtin->reg_mode = REGMODE_NONE ;
04052             free(rar) ;
04053             mri_free(rim) ; mri_fix_data_pointer(NULL,im) ; mri_free(im) ;
04054          return ;
04055       }
04056 
04057       /* copy data from registered image into output sub-brick */
04058 
04059       memcpy( rar + kk*nbar , qar , nbar ) ;
04060 
04061       mri_free(rim) ; /* don't need this no more no how */
04062    }
04063 
04064    mri_fix_data_pointer(NULL,im) ; mri_free(im) ;   /* get rid of this */
04065 
04066    /*-- attach rar to reg_dset --*/
04067 
04068    if( tt == 0 )
04069       EDIT_substitute_brick( rtin->reg_dset , 0 , rtin->datum , rar ) ;
04070    else
04071       EDIT_add_brick( rtin->reg_dset , rtin->datum , 0.0 , rar ) ;
04072 
04073    rtin->reg_nvol = tt+1 ;
04074 
04075    EDIT_dset_items( rtin->reg_dset , ADN_ntt , rtin->reg_nvol ,  ADN_none ) ;
04076 
04077    if( verbose == 2 ) fprintf(stderr,"\n") ;
04078    return ;
04079 }
04080 /*---------------------------------------------------------------------------*/
04081 #endif /* ALLOW_REGISTRATION */
04082 
04083 
04084 #ifdef ALLOW_REGISTRATION
04085 /*****************************************************************************
04086   The collection of functions for handling volume registration in this plugin.
04087 ******************************************************************************/
04088 
04089 /*---------------------------------------------------------------------------
04090    Called when the user kills the realtime graph of motion parameters
04091 -----------------------------------------------------------------------------*/
04092 
04093 void MTD_killfunc( MEM_topshell_data * mp )
04094 {
04095    /* if mp is for active input, then set the active one to nada */
04096 
04097    if( mp == NULL ) return ;
04098    if( rtinp != NULL && mp == rtinp->mp ){
04099       if( verbose ) fprintf(stderr,"RT: user killed active realtime graph\n") ;
04100       rtinp->mp = NULL ;
04101    } else {
04102       if( verbose ) fprintf(stderr,"RT: user killed inactive realtime graph\n") ;
04103    }
04104 
04105    if( mp->userdata != NULL ){ free(mp->userdata) ; mp->userdata = NULL ; }
04106    return ;
04107 }
04108 
04109 /*---------------------------------------------------------------------------
04110   Do pieces of 3D registration during realtime.
04111 -----------------------------------------------------------------------------*/
04112 
04113 void RT_registration_3D_realtime( RT_input * rtin )
04114 {
04115    int tt , ntt , ttbot ;
04116 
04117    /*-- check to see if we need to setup first --*/
04118 
04119    if( rtin->reg_3dbasis == NULL ){  /* need to setup */
04120 
04121       /* check if enough data to setup */
04122 
04123       if( rtin->reg_base_index >= rtin->nvol[0] ) return ;  /* can't setup */
04124 
04125       /* setup the registration process */
04126 
04127       if( verbose ){
04128          if( rtin->reg_mode == REGMODE_3D_ESTIM )
04129             fprintf(stderr,"RT: setting up 3D registration 'estimate'\n") ;
04130          else
04131             fprintf(stderr,"RT: setting up 3D registration 'realtime'\n") ;
04132       }
04133 
04134       SHOW_AFNI_PAUSE ;
04135       RT_registration_3D_setup( rtin ) ;  /* may take a while */
04136 
04137       if( rtin->reg_3dbasis == NULL ){
04138          fprintf(stderr,"RT: can't setup %s registration!\a\n",
04139                  REG_strings[rtin->reg_mode] ) ;
04140          if( rtin->reg_dset != NULL ){
04141             DSET_delete( rtin->reg_dset ) ; rtin->reg_dset = NULL ;
04142          }
04143          rtin->reg_mode = REGMODE_NONE ;
04144          SHOW_AFNI_READY ; return ;
04145       }
04146 
04147       /* realtime graphing? */
04148 
04149       if( rtin->reg_graph == 2 ){
04150          int    ycount = -6 ;  /* default number of graphs */
04151          static char * nar[6] = {
04152             "\\Delta I-S [mm]" , "\\Delta R-L [mm]" , "\\Delta A-P [mm]" ,
04153             "Roll [\\degree]" , "Pitch [\\degree]" , "Yaw [\\degree]"  } ;
04154 
04155          char *ttl = malloc( strlen(DSET_FILECODE(rtin->dset[0])) + 32 ) ;
04156          strcpy(ttl,"\\noesc ") ;
04157          strcat(ttl,DSET_FILECODE(rtin->dset[0])) ;
04158          if( rtin->reg_mode == REGMODE_3D_ESTIM ) strcat(ttl," [Estimate]") ;
04159 
04160          /* if p_code, only plot the reg_eval */
04161          if ( rtin->p_code )
04162             ycount = 1;
04163 
04164          rtin->mp = plot_ts_init( GLOBAL_library.dc->display ,
04165                                   0.0,rtin->reg_graph_xr-1 ,
04166                                   ycount,-rtin->reg_graph_yr,rtin->reg_graph_yr,
04167                                   "reps", NULL, ttl, nar , NULL ) ;
04168 
04169          if( rtin->mp != NULL ) rtin->mp->killfunc = MTD_killfunc ;
04170 
04171          free(ttl) ;
04172 
04173          /* set up comm for motion params    30 Mar 2004 [rickr] */
04174          RT_mp_comm_init_vars( rtin ) ;
04175          if ( rtin->mp_tcp_use ) RT_mp_comm_init( rtin ) ;
04176       }
04177    }
04178 
04179    /*-- register all sub-bricks that aren't done yet --*/
04180 
04181    ntt   = DSET_NUM_TIMES( rtin->dset[0] ) ;
04182    ttbot = rtin->reg_nvol ;
04183    for( tt=ttbot ; tt < ntt ; tt++ )
04184       RT_registration_3D_onevol( rtin , tt ) ;
04185 
04186    if( rtin->mp != NULL && ntt > ttbot ){
04187       float        * yar[7] ;
04188       int            ycount = -6 ;
04189 
04190       yar[0] = rtin->reg_rep   + ttbot ;
04191       yar[1] = rtin->reg_dx    + ttbot ;
04192       yar[2] = rtin->reg_dy    + ttbot ;
04193       yar[3] = rtin->reg_dz    + ttbot ;
04194       yar[4] = rtin->reg_phi   + ttbot ;
04195       yar[5] = rtin->reg_psi   + ttbot ;
04196       yar[6] = rtin->reg_theta + ttbot ;
04197 
04198       /* send the data off over tcp connection          30 Mar 2004 [rickr] */
04199       if ( rtin->mp_tcp_use )
04200           RT_mp_comm_send_data( rtin, yar+1, ntt-ttbot );
04201 
04202       if( ttbot > 0 )  /* modify yar after send_data, when ttbot > 0 */
04203       {
04204           int c;
04205           for ( c = 0; c < 7; c++ )  /* apply the old ttbot-- */
04206               yar[c]--;
04207           ttbot-- ;
04208       }
04209 
04210       /* if p_code, only plot the reg_eval */
04211       if ( rtin->p_code )
04212       {
04213          ycount = 1;
04214          yar[1] = rtin->reg_eval + ttbot;
04215       }
04216 
04217       plot_ts_addto( rtin->mp , ntt-ttbot , yar[0] , ycount , yar+1 ) ;
04218    }
04219 
04220    /*-- my work here is done --*/
04221 
04222    XmUpdateDisplay( THE_TOPSHELL ) ;
04223    SHOW_AFNI_READY ; return ;
04224 }
04225 
04226 /*---------------------------------------------------------------------------
04227   Do 3D registration of all slices at once, when the realtime dataset
04228   is all present and accounted for.
04229 -----------------------------------------------------------------------------*/
04230 
04231 void RT_registration_3D_atend( RT_input * rtin )
04232 {
04233    int tt , ntt ;
04234 
04235    /* check if have enough data to register as ordered */
04236 
04237    if( rtin->reg_base_index >= rtin->nvol[0] ){
04238       fprintf(stderr,"RT: can't do %s registration: not enough 3D volumes!\a\n",
04239               REG_strings[rtin->reg_mode] ) ;
04240       DSET_delete( rtin->reg_dset ) ; rtin->reg_dset = NULL ;
04241       rtin->reg_mode = REGMODE_NONE ;
04242       return ;
04243    }
04244 
04245    /* set up the registration process */
04246 
04247    if( verbose )
04248       fprintf(stderr,"RT: starting 3D registration 'at end'\n") ;
04249 
04250    SHOW_AFNI_PAUSE ;
04251    RT_registration_3D_setup( rtin ) ;
04252 
04253    if( rtin->reg_3dbasis == NULL ){
04254       fprintf(stderr,"RT: can't setup %s registration!\a\n",
04255               REG_strings[rtin->reg_mode] ) ;
04256       DSET_delete( rtin->reg_dset ) ; rtin->reg_dset = NULL ;
04257       rtin->reg_mode = REGMODE_NONE ;
04258       SHOW_AFNI_READY ; return ;
04259    }
04260 
04261    /* register each volume into the new dataset */
04262 
04263    ntt = DSET_NUM_TIMES( rtin->dset[0] ) ;
04264    if( verbose == 1 ) fprintf(stderr,"RT: ") ;
04265    for( tt=0 ; tt < ntt ; tt++ ){
04266       XmUpdateDisplay( THE_TOPSHELL ) ;
04267       RT_registration_3D_onevol( rtin , tt ) ;
04268       if( verbose == 1 ) fprintf(stderr,"%d",tt%10) ;
04269    }
04270    if( verbose == 1 ) fprintf(stderr,"\n") ;
04271 
04272    /* un-setup the registration process */
04273 
04274    RT_registration_3D_close( rtin ) ;
04275    if( verbose ) SHOW_TIMES ;
04276    SHOW_AFNI_READY ; return ;
04277 }
04278 
04279 /*-----------------------------------------------------------------------
04280    Setup the 3D registration data
04281 -------------------------------------------------------------------------*/
04282 
04283 #define VL_MIT  9     /* iterations */
04284 #define VL_DXY  0.05  /* voxels */
04285 #define VL_DPH  0.07  /* degrees */
04286 #define VL_DEL  0.70  /* voxels */
04287 
04288 void RT_registration_3D_setup( RT_input * rtin )
04289 {
04290    int ibase = rtin->reg_base_index ;
04291    int kk , nx,ny,nz , kind , nbar ;
04292    MRI_IMAGE * im ;
04293    char * bar , * ept ;
04294 
04295    /*-- extract info about coordinate axes of dataset --*/
04296 
04297    rtin->iha  = THD_handedness( rtin->dset[0] )   ;  /* LH or RH? */
04298 
04299    rtin->ax1  = THD_axcode( rtin->dset[0] , 'I' ) ;
04300    rtin->hax1 = rtin->ax1 * rtin->iha      ;      /* roll */
04301 
04302    rtin->ax2  = THD_axcode( rtin->dset[0] , 'R' ) ;
04303    rtin->hax2 = rtin->ax2 * rtin->iha      ;      /* pitch */
04304 
04305    rtin->ax3  = THD_axcode( rtin->dset[0] , 'A' ) ;
04306    rtin->hax3 = rtin->ax3 * rtin->iha      ;      /* yaw */
04307 
04308    im     = DSET_BRICK(rtin->dset[0],ibase) ;
04309    im->dx = fabs( DSET_DX(rtin->dset[0]) ) ;  /* must set voxel dimensions */
04310    im->dy = fabs( DSET_DY(rtin->dset[0]) ) ;
04311    im->dz = fabs( DSET_DZ(rtin->dset[0]) ) ;
04312 
04313    switch( rtin->reg_mode ){
04314       default: rtin->reg_3dbasis = NULL ;   /* should not occur */
04315       return ;
04316 
04317       case REGMODE_3D_RTIME:                /* actual registration */
04318       case REGMODE_3D_ATEND:
04319          ept = getenv("AFNI_REALTIME_volreg_maxite") ;
04320          kk  = VL_MIT ;
04321          if( ept != NULL ){
04322             kk = strtol(ept,NULL,10) ; if( kk <= 0 ) kk = VL_MIT ;
04323          }
04324          mri_3dalign_params( kk , VL_DXY , VL_DPH , VL_DEL ,
04325                              abs(rtin->ax1)-1 , abs(rtin->ax2)-1 , abs(rtin->ax3)-1 , -1 ) ;
04326 
04327                                                            /* noreg , clipit */
04328          mri_3dalign_method( rtin->reg_resam , verbose==2 ,   0     , 1        ) ;
04329 
04330          mri_3dalign_final_regmode( rtin->reg_final_resam ) ;
04331 
04332          rtin->reg_3dbasis = mri_3dalign_setup( im ,  NULL ) ;
04333       break ;
04334 
04335       case REGMODE_3D_ESTIM:               /* just estimate motion */
04336          ept = getenv("AFNI_REALTIME_volreg_maxite_est") ;
04337          kk  = 1 ;
04338          if( ept != NULL ){
04339             kk = strtol(ept,NULL,10) ; if( kk <= 0 ) kk = 1 ;
04340          }
04341 
04342          mri_3dalign_params( kk , VL_DXY , VL_DPH , 2.0*VL_DEL ,
04343                              abs(rtin->ax1)-1 , abs(rtin->ax2)-1 , abs(rtin->ax3)-1 , -1 ) ;
04344 
04345          kk = MRI_CUBIC ;                     /* noreg , clipit */
04346          mri_3dalign_method( kk , verbose==2 ,   1     , 0        ) ;
04347 
04348          rtin->reg_3dbasis = mri_3dalign_setup( im ,  NULL ) ;
04349 #if 0
04350          kk = MRI_LINEAR ;
04351          mri_3dalign_method( kk , verbose==2 ,   1     , 0        ) ;
04352 #endif
04353       break ;
04354    }
04355 
04356    rtin->reg_nvol  = 0 ;
04357    return ;
04358 }
04359 
04360 /*---------------------------------------------------------------------------
04361     Undo the routine just above!
04362 -----------------------------------------------------------------------------*/
04363 
04364 void RT_registration_3D_close( RT_input * rtin )
04365 {
04366    mri_3dalign_cleanup( rtin->reg_3dbasis ) ;
04367    rtin->reg_3dbasis = NULL ;
04368    return ;
04369 }
04370 
04371 /*-----------------------------------------------------------------------
04372    Carry out the 3D registration for the tt-th sub-brick
04373 -------------------------------------------------------------------------*/
04374 
04375 #ifndef D2RFAC
04376 #  define D2RFAC (PI/180.0)
04377 #  define R2DFAC (180.0/PI)
04378 #endif
04379 
04380 void RT_registration_3D_onevol( RT_input * rtin , int tt )
04381 {
04382    MRI_IMAGE * rim , * qim ;
04383    char * qar ;
04384    float dx,dy,dz , roll,pitch,yaw , ddx,ddy,ddz ;
04385    int   nest ;
04386 
04387    /*-- sanity check --*/
04388 
04389    if( rtin->dset[0] == NULL ) return ;
04390 
04391    /*-- actual registration --*/
04392 
04393    if( verbose == 2 )
04394       fprintf(stderr,"RT: 3D registering sub-brick %d\n",tt) ;
04395 
04396    qim     = DSET_BRICK(rtin->dset[0],tt) ;
04397    qim->dx = fabs( DSET_DX(rtin->dset[0]) ) ;  /* must set voxel dimensions */
04398    qim->dy = fabs( DSET_DY(rtin->dset[0]) ) ;
04399    qim->dz = fabs( DSET_DZ(rtin->dset[0]) ) ;
04400 
04401    rim = mri_3dalign_one( rtin->reg_3dbasis , qim ,
04402                           &roll , &pitch , &yaw , &dx , &dy , &dz ) ;
04403 
04404    /*-- massage and store movement parameters --*/
04405 
04406    roll  *= R2DFAC ; if( rtin->hax1 < 0 ) roll  = -roll  ;
04407    pitch *= R2DFAC ; if( rtin->hax2 < 0 ) pitch = -pitch ;
04408    yaw   *= R2DFAC ; if( rtin->hax3 < 0 ) yaw   = -yaw   ;
04409 
04410    switch( rtin->dset[0]->daxes->xxorient ){
04411       case ORI_R2L_TYPE: ddy =  dx; break ; case ORI_L2R_TYPE: ddy = -dx; break ;
04412       case ORI_P2A_TYPE: ddz = -dx; break ; case ORI_A2P_TYPE: ddz =  dx; break ;
04413       case ORI_I2S_TYPE: ddx =  dx; break ; case ORI_S2I_TYPE: ddx = -dx; break ;
04414    }
04415 
04416    switch( rtin->dset[0]->daxes->yyorient ){
04417       case ORI_R2L_TYPE: ddy =  dy; break ; case ORI_L2R_TYPE: ddy = -dy; break ;
04418       case ORI_P2A_TYPE: ddz = -dy; break ; case ORI_A2P_TYPE: ddz =  dy; break ;
04419       case ORI_I2S_TYPE: ddx =  dy; break ; case ORI_S2I_TYPE: ddx = -dy; break ;
04420    }
04421 
04422    switch( rtin->dset[0]->daxes->zzorient ){
04423       case ORI_R2L_TYPE: ddy =  dz; break ; case ORI_L2R_TYPE: ddy = -dz; break ;
04424       case ORI_P2A_TYPE: ddz = -dz; break ; case ORI_A2P_TYPE: ddz =  dz; break ;
04425       case ORI_I2S_TYPE: ddx =  dz; break ; case ORI_S2I_TYPE: ddx = -dz; break ;
04426    }
04427 
04428    nest = rtin->reg_nest ;
04429    rtin->reg_tim = (float *) realloc( (void *) rtin->reg_tim ,
04430                                       sizeof(float) * (nest+1) ) ;
04431    rtin->reg_dx  = (float *) realloc( (void *) rtin->reg_dx  ,
04432                                       sizeof(float) * (nest+1) ) ;
04433    rtin->reg_dy  = (float *) realloc( (void *) rtin->reg_dy  ,
04434                                       sizeof(float) * (nest+1) ) ;
04435    rtin->reg_dz  = (float *) realloc( (void *) rtin->reg_dz  ,
04436                                       sizeof(float) * (nest+1) ) ;
04437    rtin->reg_phi = (float *) realloc( (void *) rtin->reg_phi ,
04438                                       sizeof(float) * (nest+1) ) ;
04439    rtin->reg_psi = (float *) realloc( (void *) rtin->reg_psi ,
04440                                       sizeof(float) * (nest+1) ) ;
04441    rtin->reg_theta = (float *) realloc( (void *) rtin->reg_theta ,
04442                                       sizeof(float) * (nest+1) ) ;
04443    rtin->reg_rep   = (float *) realloc( (void *) rtin->reg_rep ,
04444                                       sizeof(float) * (nest+1) ) ;
04445    rtin->reg_eval   = (float *) realloc( (void *) rtin->reg_eval ,
04446                                       sizeof(float) * (nest+1) ) ;
04447 
04448    rtin->reg_tim[nest]   = THD_timeof_vox( tt , 0 , rtin->dset[0] ) ;
04449    rtin->reg_dx [nest]   = ddx   ;
04450    rtin->reg_dy [nest]   = ddy   ;
04451    rtin->reg_dz [nest]   = ddz   ;
04452    rtin->reg_phi[nest]   = roll  ;
04453    rtin->reg_psi[nest]   = pitch ;
04454    rtin->reg_theta[nest] = yaw   ;
04455    rtin->reg_rep[nest]   = tt    ;
04456 
04457    /* evaluate parser expression, since all input values are assigned here */
04458    if ( rtin->p_code )
04459    {
04460        int c ;
04461 
04462        /* clear extras, just to be safe */
04463        for ( c = 6; c < 26; c++ ) rtin->p_atoz[c] = 0;
04464 
04465        rtin->p_atoz[0] = ddx ; rtin->p_atoz[1] = ddy  ; rtin->p_atoz[2] = ddz;
04466        rtin->p_atoz[3] = roll; rtin->p_atoz[4] = pitch; rtin->p_atoz[5] = yaw;
04467 
04468        /* actually set the value via parser evaluation */
04469        rtin->reg_eval[nest] = PARSER_evaluate_one(rtin->p_code, rtin->p_atoz);
04470    }
04471 
04472    rtin->reg_nest ++ ; rtin->reg_nvol = tt+1 ;
04473 
04474    /*-- convert output image to desired type;
04475         set qar to point to data in the converted registered image --*/
04476 
04477    if( rim != NULL && rtin->reg_dset != NULL ){
04478       switch( rtin->datum ){
04479          case MRI_float:                        /* rim is already floats */
04480             qar = (char *) MRI_FLOAT_PTR(rim) ;
04481          break ;
04482 
04483          case MRI_short:
04484             qim = mri_to_short(1.0,rim) ; mri_free(rim) ; rim = qim ;
04485             qar = (char *) MRI_SHORT_PTR(rim) ;
04486          break ;
04487 
04488          case MRI_byte:
04489             qim = mri_to_byte(rim) ; mri_free(rim) ; rim = qim ;
04490             qar = (char *) MRI_BYTE_PTR(rim) ;
04491          break ;
04492 
04493          /* this case should not occur: */
04494 
04495          default:
04496             fprintf(stderr,"RT: can't do registration on %s images!\a\n",
04497                     MRI_TYPE_name[rtin->datum] ) ;
04498             DSET_delete( rtin->reg_dset ) ; rtin->reg_dset = NULL ;
04499             rtin->reg_mode = REGMODE_NONE ;
04500             mri_free(rim) ;
04501          return ;
04502       }
04503 
04504       /* attach to reg_dset */
04505 
04506       if( tt == 0 )
04507          EDIT_substitute_brick( rtin->reg_dset , 0 , rtin->datum , qar ) ;
04508       else
04509          EDIT_add_brick( rtin->reg_dset , rtin->datum , 0.0 , qar ) ;
04510 
04511       EDIT_dset_items( rtin->reg_dset , ADN_ntt , rtin->reg_nvol ,  ADN_none ) ;
04512 
04513       mri_fix_data_pointer(NULL,rim) ; mri_free(rim) ;   /* get rid of this */
04514    }
04515 
04516    return ;
04517 }
04518 /*---------------------------------------------------------------------------*/
04519 #endif /* ALLOW_REGISTRATION */
04520 
04521 #ifndef FIM_THR
04522 #define FIM_THR  0.0999
04523 #endif
04524 
04525 #define MIN_UPDT 5
04526 
04527 /*--------------------------------------------------------------------------
04528   Recursively update the functional image computation.
04529     mode = (if >= 0) index of sub-brick to process this time
04530            (if <  0) one of the MODE codes above:
04531              INIT_MODE:   create a new dataset, and initialize workspaces
04532              UPDATE_MODE: put values into dataset bricks
04533              FINAL_MODE:  clean up workspaces
04534 
04535   The return value is 1 if all is OK, -1 if something bad happens.
04536 
04537   Adapted 27 Apr 1997 from AFNI_fimmer_compute in afni_fimmer.c -- only
04538   the computational stuff is left; the display stuff is relegated to
04539   another routine.
04540 -----------------------------------------------------------------------------*/
04541 
04542 #define IFree(x) do{ if( (x)!=NULL ){ free(x) ; (x)=NULL ; } } while(0)
04543 
04544 int RT_fim_recurse( RT_input *rtin , int mode )
04545 {
04546    static Three_D_View * im3d ;
04547    static THD_3dim_dataset * dset_time ;
04548    static MRI_IMAGE * ref_ts ;
04549    static MRI_IMAGE * ort_ts ;
04550 
04551    static THD_3dim_dataset * new_dset ;
04552    static int nvox , ngood_ref , it1 , dtyp , nxyz , nupdt ;
04553    static float * vval=NULL , * aval=NULL , * rbest=NULL , * abest=NULL ;
04554    static int   * indx=NULL ;
04555    static PCOR_references ** pc_ref=NULL ;
04556    static PCOR_voxel_corr ** pc_vc=NULL ;
04557    static int nx_ref , ny_ref ;
04558 
04559    static int fim_nref , nx_ort , ny_ort , internal_ort ;
04560    static float * ortar ;
04561    static float * ref_vec = NULL ;
04562    static int    nref_vec = -666 ;
04563 
04564    static int polort = 1 ;  /* 02 Jun 1999 */
04565 
04566    static char new_prefix[THD_MAX_PREFIX] ;
04567    static char old_prefix[THD_MAX_PREFIX] ;
04568 
04569    static int first_pass = 1;     /* flag for first FIM computation   30 Oct 2003 [rickr] */
04570    int        start_index;        /* in case of loop FIM computation  30 Oct 2003 [rickr] */
04571 
04572    int ifim, it, iv, ivec, nnow, itbot , ip ;
04573    float fthr , topval , stataux[MAX_STAT_AUX] ;
04574    float * tsar ;
04575    short * bar ;
04576    void  * ptr ;
04577 
04578    char strbuf[128] ;  /* 03 Jun 1999 */
04579 
04580    /******************/
04581    /** sanity check **/
04582 
04583    if( rtin == NULL ) return -1 ;
04584 
04585    /****************************/
04586    /***** initialize stuff *****/
04587 
04588    if( mode == INIT_MODE ){
04589 
04590       dset_time = rtin->dset[0] ;       /* assign dset_time first   30 Oct 2003 [rickr] */
04591 
04592 #ifdef ALLOW_REGISTRATION
04593       /* now check for use of registered dataset                    30 Oct 2003 [rickr] */
04594       if( (rtin->reg_mode == REGMODE_2D_RTIME) || (rtin->reg_mode == REGMODE_3D_RTIME) )
04595          dset_time = rtin->reg_dset;
04596 #endif /* ALLOW_REGISTRATION */
04597 
04598       if( dset_time == NULL ) return -1 ;
04599 
04600       im3d      = rtin->im3d[0] ;         if( im3d      == NULL ) return -1 ;
04601       ref_ts    = im3d->fimdata->fimref ; if( ref_ts    == NULL ) return -1 ;
04602       ort_ts    = im3d->fimdata->fimort ; /* NULL is legal here */
04603       nupdt     = 0 ;                     /* number of updates done yet */
04604 
04605       polort    = im3d->fimdata->polort ; /* 02 Jun 1999 */
04606 
04607       if( ort_ts != NULL ){               /* load some dimensions */
04608          nx_ort = ort_ts->nx ;
04609          ny_ort = ort_ts->ny ;
04610          ortar  = MRI_FLOAT_PTR(ort_ts) ;
04611          internal_ort = 0 ;
04612       } else {
04613          nx_ort = ny_ort = 0 ;
04614          ortar  = NULL ;
04615          internal_ort = 1 ;
04616       }
04617       fim_nref = (internal_ort) ? (polort+2) : (ny_ort+polort+2) ;
04618 
04619       if( nref_vec < fim_nref ){  /* allocate space for ref work vector */
04620           ref_vec = (float *) XtRealloc( (char *)ref_vec, sizeof(float)*fim_nref ) ;
04621          nref_vec = fim_nref ;
04622       }
04623 
04624       itbot     = im3d->fimdata->init_ignore ;  /* certainly don't use these points! */
04625       nx_ref    = ref_ts->nx ;                  /* load dimensions */
04626       ny_ref    = ref_ts->ny ;
04627       ngood_ref =  0 ;
04628       it1       = -1 ;
04629       for( ivec=0 ; ivec < ny_ref ; ivec++ ){           /* scan ref vectors */
04630          tsar = MRI_FLOAT_PTR(ref_ts) + (ivec*nx_ref) ; /* ptr to ivec-th vector */
04631          ifim = 0 ;                                     /* count # good entries */
04632          for( it=itbot ; it < nx_ref ; it++ ){
04633             if( tsar[it] < WAY_BIG ){ ifim++ ; if( it1 < 0 ) it1 = it ; }
04634          }
04635 
04636          if( ifim < MIN_UPDT ){
04637             fprintf(stderr,"RTfim: ideal timeseries has too few good entries!\a\n") ;
04638             return -1 ;
04639          }
04640 
04641          ngood_ref = MAX( ifim , ngood_ref ) ;
04642       }
04643 
04644       /** at this point, ngood_ref = max number of good reference points,
04645           and                  it1 = index of first point used in first reference **/
04646 
04647       dtyp = DSET_BRICK_TYPE(dset_time,0) ;  /* type of each voxel */
04648       if( ! AFNI_GOOD_FUNC_DTYPE(dtyp) ){
04649          fprintf(stderr,"RTfim: input dataset has unsupported datum %s!\a\n",
04650                  MRI_TYPE_name[dtyp] ) ;
04651          return -1 ;
04652       }
04653 
04654       if( nx_ort > 0 && nx_ort < nx_ref ){  /* check for compatibility */
04655          fprintf(stderr,
04656                  "RTfim: WARNING -- ort length=%d  ideal length=%d\n",
04657                  nx_ort , nx_ref ) ;
04658       }
04659 
04660       /*--- Create a new prefix ---*/
04661 
04662       MCW_strncpy( old_prefix , DSET_PREFIX(dset_time) , THD_MAX_PREFIX-3 ) ;
04663 
04664       for( ifim=1 ; ifim < 99 ; ifim++ ){
04665          sprintf( new_prefix , "%s@%d" , old_prefix , ifim ) ;
04666          if( PLUTO_prefix_ok(new_prefix) ) break ;
04667       }
04668       if( ifim == 99 ){
04669          fprintf(stderr,"RTfim: can't create new prefix!\a\n") ;
04670          return -1 ;
04671       }
04672 
04673       nxyz =  dset_time->dblk->diskptr->dimsizes[0]    /* total voxel count */
04674             * dset_time->dblk->diskptr->dimsizes[1]
04675             * dset_time->dblk->diskptr->dimsizes[2] ;
04676 
04677       return 1 ;
04678    }
04679 
04680    /**************************/
04681    /***** finalize stuff *****/
04682 
04683    if( mode == FINAL_MODE ){
04684 
04685       /*--- End of recursive updates; now free temporary workspaces ---*/
04686 
04687       if( pc_ref != NULL ){
04688          for( ivec=0 ; ivec < ny_ref ; ivec++ ){
04689             free_PCOR_references(pc_ref[ivec]) ;
04690             free_PCOR_voxel_corr(pc_vc[ivec]) ;
04691          }
04692       }
04693       IFree(vval) ; IFree(indx)  ; IFree(pc_ref) ; IFree(pc_vc)  ;
04694       IFree(aval) ; IFree(rbest) ; IFree(abest)  ;
04695 
04696       first_pass = 1;                /* reset for future datasets    30 Oct 2003 [rickr] */
04697 
04698       return 1 ;
04699    }
04700 
04701    /*********************************************/
04702    /***** update the dataset under creation *****/
04703 
04704    if( mode == UPDATE_MODE ){
04705 
04706       if( nupdt < MIN_UPDT ) return 1 ;  /* can't do anything yet */
04707 
04708       rtin->func_condit = 2 ;  /* updated at least once */
04709 
04710       /*--- set the statistical parameters ---*/
04711 
04712       stataux[0] = nupdt ;               /* number of points used */
04713       stataux[1] = (ny_ref==1) ? 1 : 2 ; /* number of references  */
04714       stataux[2] = fim_nref - 1 ;        /* number of orts        */
04715       for( iv=3 ; iv < MAX_STAT_AUX ; iv++ ) stataux[iv] = 0.0 ;
04716 
04717       (void) EDIT_dset_items( new_dset, ADN_stat_aux, stataux, ADN_none ) ;
04718 
04719       /*** Compute brick arrays for new dataset ***/
04720 
04721       if( ny_ref == 1 ){
04722 
04723       /*** Just 1 ref vector --> load values directly into dataset ***/
04724 
04725          /*--- get alpha (coef) into vval,
04726                find max value, scale into brick array ---*/
04727 
04728          PCOR_get_coef( pc_ref[0] , pc_vc[0] , vval ) ;
04729 
04730          topval = 0.0 ;
04731          for( iv=0 ; iv < nvox ; iv++ )
04732             if( fabs(vval[iv]) > topval ) topval = fabs(vval[iv]) ;
04733 
04734          bar = DSET_ARRAY( new_dset , FUNC_ival_fim[FUNC_COR_TYPE] ) ;
04735 
04736 #ifdef DONT_USE_MEMCPY
04737          for( iv=0 ; iv < nxyz ; iv++ ) bar[iv] = 0 ;
04738 #else
04739          memset( bar , 0 , sizeof(short)*nxyz ) ;
04740 #endif
04741 
04742          if( topval > 0.0 ){
04743             topval = MRI_TYPE_maxval[MRI_short] / topval ;
04744             for( iv=0 ; iv < nvox ; iv++ )
04745                bar[indx[iv]] = (short)(topval * vval[iv] + 0.499) ;
04746 
04747             stataux[0] = 1.0/topval ;
04748          } else {
04749             stataux[0] = 0.0 ;
04750          }
04751 
04752          /*--- get correlation coefficient (pcor) into vval,
04753                scale into brick array (with fixed scaling factor) ---*/
04754 
04755          PCOR_get_pcor( pc_ref[0] , pc_vc[0] , vval ) ;
04756 
04757          bar = DSET_ARRAY( new_dset , FUNC_ival_thr[FUNC_COR_TYPE] ) ;
04758 
04759 #ifdef DONT_USE_MEMCPY
04760          for( iv=0 ; iv < nxyz ; iv++ ) bar[iv] = 0 ;
04761 #else
04762          memset( bar , 0 , sizeof(short)*nxyz ) ;
04763 #endif
04764 
04765          for( iv=0 ; iv < nvox ; iv++ )
04766             bar[indx[iv]] = (short)(FUNC_COR_SCALE_SHORT * vval[iv] + 0.499) ;
04767 
04768          stataux[1] = 1.0 / FUNC_COR_SCALE_SHORT ;
04769 
04770       } else {
04771 
04772       /*** Multiple references --> find best correlation at each voxel ***/
04773 
04774          /*--- get first ref results into abest and rbest (best so far) ---*/
04775 
04776          PCOR_get_coef( pc_ref[0] , pc_vc[0] , abest ) ;
04777          PCOR_get_pcor( pc_ref[0] , pc_vc[0] , rbest ) ;
04778 
04779          /*--- for each succeeding ref vector,
04780                get results into aval and vval,
04781                if |vval| > |rbest|, then use that result instead ---*/
04782 
04783          for( ivec=1 ; ivec < ny_ref ; ivec++ ){
04784 
04785             PCOR_get_coef( pc_ref[ivec] , pc_vc[ivec] , aval ) ;
04786             PCOR_get_pcor( pc_ref[ivec] , pc_vc[ivec] , vval ) ;
04787 
04788             for( iv=0 ; iv < nvox ; iv++ ){
04789                if( fabs(vval[iv]) > fabs(rbest[iv]) ){
04790                   rbest[iv] = vval[iv] ;
04791                   abest[iv] = aval[iv] ;
04792                }
04793             }
04794          }
04795 
04796          /*--- at this point, abest and rbest are the best
04797                results, so scale them into the dataset bricks ---*/
04798 
04799          topval = 0.0 ;
04800          for( iv=0 ; iv < nvox ; iv++ )
04801             if( fabs(abest[iv]) > topval ) topval = fabs(abest[iv]) ;
04802 
04803          bar = DSET_ARRAY( new_dset , FUNC_ival_fim[FUNC_COR_TYPE] ) ;
04804 
04805 #ifdef DONT_USE_MEMCPY
04806          for( iv=0 ; iv < nxyz ; iv++ ) bar[iv] = 0 ;
04807 #else
04808          memset( bar , 0 , sizeof(short)*nxyz ) ;
04809 #endif
04810 
04811          if( topval > 0.0 ){
04812             topval = MRI_TYPE_maxval[MRI_short] / topval ;
04813             for( iv=0 ; iv < nvox ; iv++ )
04814                bar[indx[iv]] = (short)(topval * abest[iv] + 0.499) ;
04815 
04816             stataux[0] = 1.0/topval ;
04817          } else {
04818             stataux[0] = 0.0 ;
04819          }
04820 
04821          bar = DSET_ARRAY( new_dset , FUNC_ival_thr[FUNC_COR_TYPE] ) ;
04822 
04823 #ifdef DONT_USE_MEMCPY
04824          for( iv=0 ; iv < nxyz ; iv++ ) bar[iv] = 0 ;
04825 #else
04826          memset( bar , 0 , sizeof(short)*nxyz ) ;
04827 #endif
04828 
04829          for( iv=0 ; iv < nvox ; iv++ )
04830             bar[indx[iv]] = (short)(FUNC_COR_SCALE_SHORT * rbest[iv] + 0.499) ;
04831 
04832          stataux[1] = 1.0 / FUNC_COR_SCALE_SHORT ;
04833       }
04834 
04835       /* set scaling factors for each sub-brick */
04836 
04837       (void) EDIT_dset_items( new_dset, ADN_brick_fac, stataux, ADN_none ) ;
04838       return 1 ;
04839    }
04840 
04841    /*********************************************************/
04842    /***** at this point, must process time point # mode *****/
04843 
04844    if( mode < it1 ) return 1 ;  /* nothing to do before first time */
04845 
04846    /*---- at first time, find values above threshold to fim:
04847             find the mean of the first array,
04848             compute the threshold (fthr) from it,
04849             make indx[i] be the 3D index of the i-th voxel above threshold ----*/
04850 
04851    /* be sure dset_time has enough sub-bricks for this time point   30 Oct 2003 [rickr] */
04852    if ( (DSET_NVALS(dset_time) <= mode) || (DSET_ARRAY(dset_time,mode) == NULL) )
04853       return 1;
04854 
04855 #if 0                           /* change to first_pass test */
04856    if( mode == it1 )
04857 #endif
04858 
04859    /* check for first fim computation with this dataset            30 Oct 2003 [rickr] */
04860    if( first_pass == 1 ){
04861 
04862       switch( dtyp ){  /* process each datum type separately */
04863 
04864          case MRI_short:{
04865             short * dar = (short *) DSET_ARRAY(dset_time,it1) ;
04866             for( iv=0,fthr=0.0 ; iv < nxyz ; iv++ ) fthr += abs(dar[iv]) ;
04867             fthr = FIM_THR * fthr / nxyz ;
04868             for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
04869                if( abs(dar[iv]) > fthr ) nvox++ ;
04870             indx = (int *) malloc( sizeof(int) * nvox ) ;
04871             if( indx == NULL ){
04872                fprintf(stderr,"RTfim: indx malloc failure!\a\n") ;
04873                EXIT(1) ;
04874             }
04875             for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
04876                if( abs(dar[iv]) > fthr ) indx[nvox++] = iv ;
04877          }
04878          break ;
04879 
04880          case MRI_float:{
04881             float * dar = (float *) DSET_ARRAY(dset_time,it1) ;
04882             for( iv=0,fthr=0.0 ; iv < nxyz ; iv++ ) fthr += fabs(dar[iv]) ;
04883             fthr = FIM_THR * fthr / nxyz ;
04884             for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
04885                if( fabs(dar[iv]) > fthr ) nvox++ ;
04886             indx = (int *) malloc( sizeof(int) * nvox ) ;
04887             if( indx == NULL ){
04888               fprintf(stderr,"RTfim: indx malloc failure!\a\n") ;
04889               EXIT(1) ;
04890             }
04891             for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
04892                if( fabs(dar[iv]) > fthr ) indx[nvox++] = iv ;
04893          }
04894          break ;
04895 
04896          case MRI_byte:{
04897             byte * dar = (byte *) DSET_ARRAY(dset_time,it1) ;
04898             for( iv=0,fthr=0.0 ; iv < nxyz ; iv++ ) fthr += dar[iv] ;
04899             fthr = FIM_THR * fthr / nxyz ;
04900             for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
04901                if( dar[iv] > fthr ) nvox++ ;
04902             indx = (int *) malloc( sizeof(int) * nvox ) ;
04903             if( indx == NULL ){
04904               fprintf(stderr,"RTfim: indx malloc failure!\a\n") ;
04905               EXIT(1) ;
04906             }
04907             for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
04908                if( dar[iv] > fthr ) indx[nvox++] = iv ;
04909          }
04910          break ;
04911       }
04912 
04913       /** allocate space for voxel values that are above the threshold **/
04914       /** (nvox = # voxels above threshold;  nxyz = total # voxels)    **/
04915 
04916       if( verbose == 2 )
04917          fprintf(stderr,"RTfim: %d/%d voxels being FIMmed.\n",nvox,nxyz) ;
04918       VMCHECK ;
04919 
04920       vval = (float *) malloc( sizeof(float) * nvox) ;
04921       if( vval == NULL ){
04922         fprintf(stderr,"RTfim: vval malloc failure!\a\n") ;
04923         EXIT(1) ;
04924       }
04925 
04926       /** allocate extra space for comparing results from multiple ref vectors **/
04927 
04928       if( ny_ref > 1 ){
04929          aval  = (float *) malloc( sizeof(float) * nvox) ;
04930          rbest = (float *) malloc( sizeof(float) * nvox) ;
04931          abest = (float *) malloc( sizeof(float) * nvox) ;
04932          if( aval==NULL || rbest==NULL || abest==NULL ){
04933            fprintf(stderr,"RTfim: abest malloc failure!\a\n") ;
04934            EXIT(1) ;
04935          }
04936       } else {
04937          aval = rbest = abest = NULL ;
04938       }
04939 
04940       /*--- initialize recursive updates ---*/
04941 
04942       pc_ref = (PCOR_references **) malloc( sizeof(PCOR_references *) * ny_ref ) ;
04943       pc_vc  = (PCOR_voxel_corr **) malloc( sizeof(PCOR_voxel_corr *) * ny_ref ) ;
04944 
04945       if( pc_ref == NULL || pc_vc == NULL ){
04946         fprintf(stderr,"RTfim: can't malloc recursion space!\a\n") ;
04947         EXIT(1) ;
04948       }
04949 
04950       for( ivec=0 ; ivec < ny_ref ; ivec++ ){
04951          pc_ref[ivec] = new_PCOR_references( fim_nref ) ;
04952          pc_vc[ivec]  = new_PCOR_voxel_corr( nvox , fim_nref ) ;
04953          if( pc_ref[ivec] == NULL || pc_vc[ivec] == NULL ){
04954            fprintf(stderr,"RTfim: can't malloc refs and corr!\a\n") ;
04955            EXIT(1) ;
04956          }
04957       }
04958 
04959       /*--- Make a new dataset to hold the output ---*/
04960 
04961       rtin->func_dset = new_dset = EDIT_empty_copy( dset_time ) ;
04962       if( new_dset == NULL ){
04963         fprintf(stderr,"RTfim: can't create empty dataset!\a\n") ;
04964         EXIT(1) ;
04965       }
04966       tross_Append_History( new_dset , "plug_realtime: FIM" ) ;
04967 
04968       it = EDIT_dset_items( new_dset ,
04969                                ADN_prefix      , new_prefix ,
04970                                ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
04971                                ADN_type        , ISHEAD(dset_time)
04972                                                  ? HEAD_FUNC_TYPE : GEN_FUNC_TYPE ,
04973                                ADN_func_type   , FUNC_COR_TYPE ,
04974                                ADN_nvals       , FUNC_nvals[FUNC_COR_TYPE] ,
04975                                ADN_datum_all   , MRI_short ,
04976                                ADN_ntt         , 0 ,
04977                             ADN_none ) ;
04978 
04979       if( it > 0 )
04980          fprintf(stderr,"RTfim: WARNING -- %d errors in dataset creation!\a\n",it) ;
04981 
04982       /** attach bricks to this dataset **/
04983 
04984       for( iv=0 ; iv < new_dset->dblk->nvals ; iv++ ){
04985          ptr = malloc( DSET_BRICK_BYTES(new_dset,iv) ) ;
04986          if( ptr == NULL ){
04987            fprintf(stderr,"RTfim: can't malloc dataset bricks!\a\n") ;
04988            EXIT(1) ;
04989          }
04990          mri_fix_data_pointer( ptr ,  DSET_BRICK(new_dset,iv) ) ;
04991 
04992          sprintf(strbuf,"#%d",iv) ;              /* 03 Jun 1999 */
04993          EDIT_BRICK_LABEL(new_dset,iv,strbuf) ;  /* add a label to each brick */
04994       }
04995 
04996       DSET_lock( rtin->func_dset ) ; /* 20 Mar 1998 */
04997 
04998    }  /*---- end of initializing at first important time point ----*/
04999 
05000    /*---- do recursive update for this time point ----*/
05001 
05002    /* process as loop, in case this is the first time (and we need to backtrack) */
05003    /*                                                        30 Oct 2003 [rickr] */
05004    if ( first_pass == 1 )
05005    {
05006       first_pass  = 0;
05007       start_index = it1;      /* start from timepoint of first FIM computation */
05008    }
05009    else
05010       start_index = mode;     /* if not first pass, the "loop" will be this one iteration */
05011 
05012    for ( it = start_index; it <= mode; it++ )
05013    {
05014       nnow = 0 ;     /* number of refs used this time */
05015 
05016       if( it >= nx_ref ) return 1 ;  /* can't update without references! */
05017 
05018       for( ivec=0 ; ivec < ny_ref ; ivec++ ){           /* for each ideal time series */
05019          tsar = MRI_FLOAT_PTR(ref_ts) + (ivec*nx_ref) ;
05020          if( tsar[it] >= WAY_BIG ) continue ;           /* skip this time series */
05021 
05022          ref_vec[0] = 1.0 ;         /* we always supply ort for constant */
05023          for( ip=1 ; ip <= polort ; ip++ )              /* 02 Jun 1999:    */
05024             ref_vec[ip] = ref_vec[ip-1] * ((float)it) ; /* and polynomials */
05025 
05026          if( internal_ort ){
05027             ref_vec[ip] = tsar[it] ;                    /* ideal */
05028          } else {
05029             for( iv=0 ; iv < ny_ort ; iv++ )            /* any other orts? */
05030                ref_vec[iv+ip] = (it < nx_ort) ? ortar[it+iv*nx_ort]
05031                                               : 0.0 ;
05032 
05033             ref_vec[ny_ort+ip] = tsar[it] ;             /* ideal */
05034          }
05035 
05036          update_PCOR_references( ref_vec , pc_ref[ivec] ) ;  /* Choleski refs */
05037 
05038          /* load data from dataset into local float array vval */
05039 
05040          switch( dtyp ){
05041             case MRI_short:{
05042                short * dar = (short *) DSET_ARRAY(dset_time,it) ;
05043                for( iv=0 ; iv < nvox ; iv++ ) vval[iv] = (float) dar[indx[iv]] ;
05044             }
05045             break ;
05046 
05047             case MRI_float:{
05048                float * dar = (float *) DSET_ARRAY(dset_time,it) ;
05049                for( iv=0 ; iv < nvox ; iv++ ) vval[iv] = (float) dar[indx[iv]] ;
05050             }
05051             break ;
05052 
05053             case MRI_byte:{
05054                byte * dar = (byte *) DSET_ARRAY(dset_time,it) ;
05055                for( iv=0 ; iv < nvox ; iv++ ) vval[iv] = (float) dar[indx[iv]] ;
05056             }
05057             break ;
05058          }
05059 
05060          PCOR_update_float( vval , pc_ref[ivec] , pc_vc[ivec] ) ;  /* Choleski data */
05061          nnow++ ;
05062 
05063       }
05064       if( nnow > 0 ) nupdt++ ;  /* at least one ideal vector was used */
05065    }
05066 
05067    return 1 ;
05068 }