#include "afni.h"

#ifndef ALLOW_PLUGINS
#  error "Plugins not properly set up -- see machdep.h"
#endif

/***********************************************************************
  Plugin to plot averages over ROIs.
************************************************************************/

static int                tsgraph_num=0 ;
static MEM_topshell_data *tsgraph_mtd=NULL ;

static void ROIPLOT_tsgraph_draw   (void) ;
static void ROIPLOT_tsgraph_mtdkill( MEM_topshell_data * mp ) ;

static char * ROIPLOT_main( PLUGIN_interface * ) ;

static char helpstring[] =
   " Purpose: Plot averages over ROIs.\n"
   "\n"
   " INPUTS\n"
   " ------\n"
   " Source: Defines where the data to be averaged over the ROIs comes from.\n"
   "   Dataset = Where to get time series values from\n"
   "   Start   = 1st sub-brick index to extract\n"
   "   Stop    = Last sub-brick index to extract\n"
   "\n"
   " Two methods are available to define the ROIs from which dataset values\n"
   " will be extracted.  Only one of these can be used at a time:\n"
   "\n"
   " Cluster: Clusters of contiguous voxels are extracted (as in 3dclust),\n"
   " -------- and these are used for the ROIs.\n"
   "   From   = Dataset to define the clusters in\n"
   "   Index  = Sub-brick index to define the clusters in\n"
   "   Thresh = Minimum value to accept in the sub-brick\n"
   "   rmm   }= Clustering radius and\n"
   "   vmul  }= minimum cluster volume\n"
   "\n"
   " Values: ROIs are defined by common values stored in a dataset, as\n"
   " ------- might be generated by the 'Draw Dataset' plugin.\n"
   "   From   = Dataset to define the ROIs\n"
   "   Index  = Sub-brick index to define the ROIs\n"
   "   Min   }= Range of values to\n"
   "   Max   }= use for the ROIS\n"
   " *** N.B.: Values selection is NOT YET IMPLEMENTED ***\n"
   "\n"
   " Author -- RW Cox - April 2002"
;

/***********************************************************************
   Set up the interface to the user
************************************************************************/


DEFINE_PLUGIN_PROTOTYPE

PLUGIN_interface * PLUGIN_init( int ncall )
{
   PLUGIN_interface * plint ;

   if( ncall > 0 ) return NULL ;  /* only one interface */

   /*-- set titles and call point --*/

   plint = PLUTO_new_interface( "ROI Plot" , "Plot Average Timeseries over ROI" ,
                                 helpstring ,
                                 PLUGIN_CALL_VIA_MENU , ROIPLOT_main  ) ;

   PLUTO_short_choose(plint) ;
   PLUTO_short_number(plint) ;

   PLUTO_add_hint( plint , "Plot Average Timeseries over ROI" ) ;

   PLUTO_set_sequence( plint , "A:afniinfo:dset" ) ;

   PLUTO_set_runlabels( plint , "Plot+Keep" , "Plot+Close" ) ;  /* 04 Nov 2003 */

   /*-- first line of input --*/

   PLUTO_add_option( plint , "Source" , "Source" , TRUE ) ;
   PLUTO_add_dataset(plint , "Timeseries" ,
                                    ANAT_ALL_MASK , FUNC_ALL_MASK ,
                                    DIMEN_4D_MASK | BRICK_ALLREAL_MASK ) ;

   PLUTO_add_number( plint , "Start" , 0,99999,0 , 0    ,1 ) ;
   PLUTO_add_number( plint , "Stop"  , 0,99999,0 , 99999,1 ) ;

   /*-- second line of input --*/

   PLUTO_add_option( plint  , "Cluster" , "Cluster" , FALSE ) ;
   PLUTO_add_dataset( plint , "From" ,
                                    ANAT_ALL_MASK , FUNC_ALL_MASK ,
                                    DIMEN_ALL_MASK | BRICK_ALLREAL_MASK ) ;
   PLUTO_add_number( plint , "Index" , 0,99999,0 ,  0,1 ) ;
   PLUTO_add_number( plint , "Thresh", 0,9999 ,1 ,  1,1 ) ;
   PLUTO_add_number( plint , "rmm"   , 0,99   ,1 , 11,1 ) ;
   PLUTO_add_number( plint , "vmul"  , 0,9999,-1 , 20,1 ) ;

   /*-- third line of input --*/

   PLUTO_add_option( plint  , "Values" , "Values" , FALSE ) ;
   PLUTO_add_dataset( plint , "From" ,
                                    ANAT_ALL_MASK , FUNC_ALL_MASK ,
                                    DIMEN_ALL_MASK | BRICK_ALLREAL_MASK ) ;
   PLUTO_add_number( plint , "Index" , 0,99999,0 ,  0,1 ) ;
   PLUTO_add_number( plint , "Min"   , 0,99999,0 , 0    ,1 ) ;
   PLUTO_add_number( plint , "Max"   , 0,99999,0 , 99999,1 ) ;

   return plint ;
}

/***************************************************************************
  Main routine for this plugin (will be called from AFNI).
****************************************************************************/

static char * ROIPLOT_main( PLUGIN_interface * plint )
{
   MCW_idcode *idc ;
   THD_3dim_dataset *input_dset , *mask_dset ;
   int iv_start,iv_stop , nvox , iv_mask ;
   float thr,rmm,vmul    , dx,dy,dz ;
   int val_min , val_max , nx,ny,nz , ii,jj , nt,nc ;
   char *tag ;
   MCW_cluster_array *clustar ;
   MRI_IMAGE *flim ;
   float     *flar , **yar , *xar , xcm,ycm,zcm ;
   char      **nam_yyy ;
   MEM_plotdata * mp ;

#define MAX_NC 7

   /*--------------------------------------------------------------------*/
   /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/

   if( plint == NULL )
      return "*************************\n"
             "ROIPLOT_main:  NULL input\n"
             "*************************"  ;

   /*-- read 1st line --*/

   PLUTO_next_option(plint) ;
   idc        = PLUTO_get_idcode(plint) ;
   input_dset = PLUTO_find_dset(idc) ;
   if( input_dset == NULL )
      return "*********************************\n"
             "ROIPLOT_main:  bad Source dataset\n"
             "*********************************"  ;

   iv_start = (int) PLUTO_get_number(plint) ;
   iv_stop  = (int) PLUTO_get_number(plint) ;
   if( iv_stop >= DSET_NVALS(input_dset) ) iv_stop = DSET_NVALS(input_dset)-1 ;
   if( iv_start >= iv_stop )
      return "***********************************\n"
             "ROIPLOT_main:  bad Start/Stop range\n"
             "***********************************" ;
   DSET_load(input_dset) ;
   if( DSET_ARRAY(input_dset,iv_stop) == NULL )
      return "****************************************\n"
             "ROIPLOT_main:  can't load Source dataset\n"
             "****************************************"  ;

   nx   = input_dset->daxes->nxx ; dx = fabs(input_dset->daxes->xxdel) ;
   ny   = input_dset->daxes->nyy ; dy = fabs(input_dset->daxes->yydel) ;
   nz   = input_dset->daxes->nzz ; dz = fabs(input_dset->daxes->zzdel) ;
   nvox = nx*ny*nz ;

   /*-- read next lines --*/

   tag = PLUTO_get_optiontag(plint) ;
   if( tag == NULL )
      return "***********************************"
             "ROIPLOT_main: no ROI selection made\n"
             "***********************************"   ;

   /*--------------------------------------*/

   if( strcmp(tag,"Cluster") == 0 ){
     float *vfim ;
     MCW_cluster_array *clbig ;
     MCW_cluster *cl ;
     int ptmin ;

     idc       = PLUTO_get_idcode(plint) ;
     mask_dset = PLUTO_find_dset(idc) ;
     if( mask_dset == NULL )
       return "**********************************\n"
              "ROIPLOT_main:  bad Cluster dataset\n"
              "**********************************"  ;
     if( DSET_NVOX(mask_dset) != nvox )
       return "******************************************************\n"
              "ROIPLOT_main:  Source and Cluster datasets don't match\n"
              "******************************************************"  ;
     iv_mask = (int) PLUTO_get_number(plint) ;
     if( iv_mask >= DSET_NVALS(mask_dset) )
       return "********************************\n"
              "ROIPLOT_main:  bad Cluster index\n"
              "*********************************"  ;
     DSET_load(mask_dset) ;
     if( DSET_ARRAY(mask_dset,iv_mask) == NULL )
        return "*****************************************\n"
               "ROIPLOT_main:  can't load Cluster dataset\n"
               "*****************************************\n" ;

     thr  = PLUTO_get_number(plint) ;
     rmm  = PLUTO_get_number(plint) ;
     vmul = PLUTO_get_number(plint) ;

     vfim = (float *) malloc(sizeof(float)*nx*ny*nz) ;
     EDIT_coerce_scale_type( nvox ,
                             DSET_BRICK_FACTOR(mask_dset,iv_mask) ,
                             DSET_BRICK_TYPE(mask_dset,iv_mask) ,
                             DSET_ARRAY(mask_dset,iv_mask) ,
                             MRI_float , vfim ) ;
     if( thr > 0.0 ){
       for( ii=0 ; ii < nvox ; ii++ )
         if( fabs(vfim[ii]) < thr ) vfim[ii] = 0.0 ;
     }

     clustar = MCW_find_clusters( nx,ny,nz , dx,dy,dz ,
                                  MRI_float , vfim , rmm ) ;
     free(vfim) ;

     if( clustar == NULL )
        return "********************************\n"
               "ROIPLOT_main:  no Clusters found\n"
               "********************************\n" ;

     /* edit for volume */

     ptmin = (int) (vmul / (dx*dy*dz) + 0.99) ;
     INIT_CLARR(clbig) ;
     for( ii=0 ; ii < clustar->num_clu ; ii++ ){
       cl = clustar->clar[ii] ;
       if( cl != NULL && cl->num_pt >= ptmin ){ /* big enough */
          ADDTO_CLARR(clbig,cl) ;               /* copy pointer */
          clustar->clar[ii] = NULL ;            /* null out original */
       }
     }
     DESTROY_CLARR(clustar) ;
     clustar = clbig ;
     if( clustar == NULL || clustar->num_clu == 0 ){
       printf("** NO CLUSTERS FOUND ***\n") ;
       if( clustar != NULL ) DESTROY_CLARR(clustar) ;
       return "********************************\n"
              "ROIPLOT_main:  no Clusters found\n"
              "********************************\n" ;
     }
     SORT_CLARR(clustar) ;

   /*--------------------------------------*/

   } else if( strcmp(tag,"Values") == 0 ){
     idc       = PLUTO_get_idcode(plint) ;
     mask_dset = PLUTO_find_dset(idc) ;
     if( mask_dset == NULL )
       return "*********************************\n"
              "ROIPLOT_main:  bad Values dataset\n"
              "*********************************"  ;
     if( DSET_NVOX(mask_dset) != nvox )
       return "*****************************************************\n"
              "ROIPLOT_main:  Source and Values datasets don't match\n"
              "*****************************************************"  ;
     iv_mask = (int) PLUTO_get_number(plint) ;
     if( iv_mask >= DSET_NVALS(mask_dset) )
       return "*******************************\n"
              "ROIPLOT_main:  bad Values index\n"
              "********************************"  ;
     DSET_load(mask_dset) ;
     if( DSET_ARRAY(mask_dset,iv_mask) == NULL )
        return "****************************************\n"
               "ROIPLOT_main:  can't load Values dataset\n"
               "****************************************\n" ;

     val_min = (int) PLUTO_get_number(plint) ;
     val_max = (int) PLUTO_get_number(plint) ;
     if( val_max < val_min ) val_max = val_min ;

     /* must build cluster here */

     return "****************************************\n"
            "ROIPLOT_main:  Values NOT IMPLEMENT YET!\n"
            "****************************************\n" ;

   } else {
     return "*******************************\n"
            "ROIPLOT_main: no ROI selection?\n"
            "********************************"  ;
   }

   /* check for too many options */

   tag = PLUTO_get_optiontag(plint) ;
   if( tag != NULL )
      PLUTO_popup_transient(plint,"\n"
                                  "+++++++ WARNING ++++++++\n"
                                  " Using Cluster option;\n"
                                  " Ignoring Values option\n"  ) ;

   /* extract average time series */

   flim = THD_average_timeseries( clustar , input_dset ) ;

   if( flim == NULL ){
     DESTROY_CLARR(clustar) ;
     return "***********************************************\n"
            "ROIPLOT_main:  Can't extract average timeseries\n"
            "***********************************************\n" ;
   }

   nc   = flim->ny ; if( nc > MAX_NC ) nc = MAX_NC ;
   nt   = flim->nx ;
   flar = MRI_FLOAT_PTR(flim) ;
   yar  = (float **) malloc(sizeof(float *)*nc) ;
   for( ii=0 ; ii < nc ; ii++ ) yar[ii] = flar + (ii*nt+iv_start) ;
   xar  = (float *) malloc(sizeof(float)*nt) ;
   for( ii=0 ; ii < nt ; ii++ ) xar[ii] = ii ;
   nam_yyy = (char **) malloc(sizeof(char *)*nc) ;
   for( ii=0 ; ii < nc ; ii++ ){
      xcm = ycm = zcm = 0.0 ;
      for( jj=0 ; jj < clustar->clar[ii]->num_pt ; jj++ ){
        xcm += clustar->clar[ii]->i[ii] ;
        ycm += clustar->clar[ii]->j[ii] ;
        zcm += clustar->clar[ii]->k[ii] ;
      }
      xcm /= clustar->clar[ii]->num_pt ;
      ycm /= clustar->clar[ii]->num_pt ;
      zcm /= clustar->clar[ii]->num_pt ;
      nam_yyy[ii] = malloc(256) ;
      sprintf(nam_yyy[ii],"%d voxels\n"
                          "_{ijk: %d %d %d}" ,
              clustar->clar[ii]->num_pt ,
              (int)(xcm+0.499) , (int)(ycm+0.499) , (int)(zcm+0.499) ) ;
   }
   DESTROY_CLARR(clustar);

   mp = plot_ts_mem( iv_stop-iv_start+1 , xar+iv_start ,
                     nc,TSP_SEPARATE_YBOX,yar ,
                     NULL,NULL,NULL,nam_yyy ) ;

   for( ii=0 ; ii < nc ; ii++ ) free(nam_yyy[ii]) ;
   free(nam_yyy);
   free(xar); free(yar); mri_free(flim);

   if( mp == NULL )
     return "******************************************\n"
            "ROIPLOT_main:  Can't create in-memory plot\n"
            "******************************************\n" ;

   memplot_to_topshell( plint->im3d->dc->display , mp , NULL ) ;

   return NULL ;
}

/*------------------------------------------------------------------*/

/*------------------------------------------------------------------*/