/*----------------------------------------------------------------------
 * afni_vol2surf.c   - afni's interface to the vol2surf library
 *
 * AFNI_vol2surf_func_overlay() - similar to AFNI_vnlist_func_overlay()
 *----------------------------------------------------------------------
 */

/*----------------------------------------------------------------------
 * R. Reynolds    September, 2004
 *
 * history:
 *
 * 08 Oct 2004 [rickr]:
 *   - AFNI_vol2surf_func_overlay() has new params, surfA,surfB,use_defaults
 *   - pass use_defaults to afni_vol2surf()
 * 25 Oct 2004 [rickr]:
 *   - accept Rdata and Rthr pointers, for optionally returning the data
 *     and global threshold
 *   - make threshold masks absolute
 * 09 Aug 2006 [rickr]:
 *  - set the surface volume dataset in the global v2s_plugin_opts struct
 *  - also store the index and threshold value of the threshold sub-brick
 *----------------------------------------------------------------------
 */

#include "afni.h"
#include "vol2surf.h"

static int map_v2s_results(v2s_results *res, Three_D_View *im3d,
                           SUMA_irgba **map, int debug, int dnode );
/*-----------------------------------------------------------------------*/
/*! Create a nodal color overlay from a voxel map.
    - Surface index(es) come from global v2s_plugin_opts struct
    - Return value is number of nodes overlaid:
      -  0 return ==> no overlay was computed
      - -1 return ==> some error (e.g., no surface nodes on this dataset)
      = positive  ==> can use *map
    - *map is set to a newly malloc()-ed array (if return > 0)
    - if (Rdata) return data from results structure
    - if (Rthr) return overlay threshold from afni

    Sample usage:
    - SUMA_irgba * map;
    - int          nmap;
    - nmap = AFNI_vnlist_func_overlay( im3d, &map );
    -      if( nmap <  0 ){ ** error ** }
    - else if( nmap == 0 ){ ** nothing to show ** }
    - else                { ** show map[0..nmap-1] ** }

    * based on AFNI_vnlist_func_overlay()
-------------------------------------------------------------------------*/
int AFNI_vol2surf_func_overlay(Three_D_View *im3d, SUMA_irgba **map,
         int surfA, int surfB, int use_defaults, float ** Rdata, float * Rthr )
{
    THD_3dim_dataset * oset;            /* overlay dataset */
    THD_slist_find     find;
    THD_session      * ss;
    SUMA_surface     * sA, * sB;
    MRI_IMAGE        * im_oim;
    v2s_plugin_opts  * go;
    v2s_results      * results;
    byte             * cmask;
    int                nout, oind, debug;

ENTRY("AFNI_vol2surf_func_overlay") ;

    /* check inputs */
    if ( map == NULL || !IM3D_VALID(im3d) ) RETURN(-1);

    go = &gv2s_plug_opts;
    if ( ! use_defaults && ! go->ready ) RETURN(-1);

    debug = go->sopt.debug;         /* because I'm lazy */

    ss = im3d->ss_now;              /* session must have needed surface(s) */
    if( ss                 == NULL      ||
        ss->su_num         <= 0         ||
        surfA              <  0         ||
        ss->su_num         <= surfA     ||
        ss->su_num         <= surfB     ||
        ss->su_surf[surfA] == NULL      ||
        (surfB >= 0 && ss->su_surf[surfB] == NULL) )
    {
        if ( debug > 1 )
        {
            if ( !ss ) fprintf(stderr,"** v2s: NULL session\n");
            else
                fprintf(stderr,"** v2s: bad session data:\n"
                        "   su_num, surfA, surfB = %d, %d, %d\n",
                        ss->su_num, surfA, surfB);
        }
        RETURN(-1);
    }

    /* init return values */
    if ( Rdata ) *Rdata = NULL;
    if ( Rthr )  *Rthr  = 0.0;

    /* set surface pointers */
    sA = ss->su_surf[surfA];
    sB = ( surfB >= 0 ) ? ss->su_surf[surfB] : NULL;

    /* store the surface volume pointer            9 Aug 2006 [rickr] */
    find = PLUTO_dset_finder(sA->idcode_dset);
    go->sv_dset = find.dset;

    if ( debug )
    {
        fprintf(stderr,"++ v2s overlay: sa,sb = %d,%d\n", surfA, surfB);
        if ( debug > 1 )
            fprintf(stderr,"  surfA is %s, surfB is %s\n",
                sA->label[0] ? sA->label : "<no label>",
                sB ? (sB->label[0] ? sB->label : "<no label>") : "<not used>");
    }

    /*-------------------- overlay image --------------------*/
    oset = im3d->fim_now;
    oind = im3d->vinfo->fim_index;  /* overlay sub-brick index */

    if( oset == NULL )
    {
        if ( debug > 1 ) fprintf(stderr,"** v2s: no overlay dset\n");
        RETURN(-1);
    }
    im_oim = DSET_BRICK(oset,oind);

    if( ! im_oim || ! AFNI_GOOD_FUNC_DTYPE(im_oim->kind) )
    {
        if ( debug > 1 ) fprintf(stderr,"** v2s: no overlay image\n");
        RETURN(-1);
    }
    DSET_load(oset);                    /* to be sure */
    if( !DSET_LOADED(oset) ) RETURN(-1);

    /*-------------------- mask from threshold --------------------*/
    go->gpt_index = -1; go->gpt_thresh = 0.0;   /* init threshold options */
    cmask = NULL;

    /* if a new overlay has been created from clustering, get it as a mask
       (for now - since it's an MRI_IMAGE, and we can't pass to v2s)
                                                       11 Jan 2008 [rickr] */

    if( debug ) {
        fprintf(stderr,"-- CLUST: ");
        if( ! im3d->vinfo )
            fprintf(stderr,"no vinfo to cluster on\n");
        else if( ! oset->dblk->vedim )
            fprintf(stderr,"no dblk->vedim to cluster on\n");
        else {
            fprintf(stderr,"fim_index %d, thr_index %d, VEDIT_IVAL %d\n",
                   im3d->vinfo->fim_index,
                   im3d->vinfo->thr_index,
                   DSET_VEDIT_IVAL(oset));
            fprintf(stderr,"   (kind %d, brick type %d)\n",
                   oset->dblk->vedim->kind, DSET_BRICK_TYPE(oset,oind));
        }
    }

    if( im3d->vinfo && oset->dblk->vedim                      &&
            (im3d->vinfo->fim_index == DSET_VEDIT_IVAL(oset)) &&
            (oset->dblk->vedim->kind == DSET_BRICK_TYPE(oset,oind)) ) {
        cmask = mri_to_bytemask(oset->dblk->vedim, 1.0, 0.0);
        if( debug > 1 ) fprintf(stderr,"++ applying mask from edited image\n");
    }
    else if( im3d->vinfo->func_threshold > 0.0 &&
             im3d->vinfo->thr_onoff )               /* then want a threshold */
    {
        MRI_IMAGE * im_thr;
        float       thresh;
        int         tind = im3d->vinfo->thr_index;

        im_thr = DSET_BRICK(oset,tind);

        /* note real threshold */
        thresh = get_3Dview_func_thresh(im3d,1) ;

        /* maybe we want to return this */
        if ( Rthr ) *Rthr = thresh;

        if( im_thr && !AFNI_GOOD_FUNC_DTYPE(im_thr->kind) )
            im_thr = NULL;

        /* create the mask, if this fails, just continue... */
        if( im_thr != NULL )
        {
            int nset;
            if ( debug > 1 )
                fprintf(stderr,"++ mask from index %d and thresh %f\n",
                        tind,thresh);
            nset = thd_mask_from_brick(oset, tind, thresh, &cmask, 1);
            if ( debug > 1 )
            {
                if ( ! cmask )
                    fprintf(stderr,"-- no mask created\n");
                else
                    fprintf(stderr,"++ mask has %d set voxels\n", nset);
            }
        }
        else if ( debug > 1 )
            fprintf(stderr,"-- no threshold mask\n");

        /* note mask options for v2s command output  9 Aug 2006 [rickr] */
        if( cmask ) { go->gpt_index = tind; go->gpt_thresh = thresh; }
    } else {    /* then make the mask from non-zero overlay values */
        cmask = THD_makemask(oset, oind, 1.0, 0.0);
        /* we are applying a cmask, so make it clear in any printed command
                                                         8 Dec 2016 [rickr] */
        go->gpt_index = oind; go->gpt_thresh = 0;
    }

    /*-------------------- vol2surf computation --------------------*/
    results = afni_vol2surf(oset, oind, sA, sB, cmask, use_defaults);

    if ( cmask ) free(cmask);   /* we're done with the mask */
    if ( ! results )
    {
        if ( debug ) fprintf(stderr,"-- vol2surf failure\n");
        RETURN(-1);             /* failure? */
    }

    /*-------------------- lose old vnlist --------------------*/
    /* vol2surf has no method for counting voxels (right now)  */
    if( sA->vn )
    {
        if ( debug > 1 ) fprintf(stderr,"-- destroying vnlist\n");
        SUMA_destroy_vnlist( sA->vn ) ;
        sA->vn = NULL;
    }

    /*-------------------- set overlay colors --------------------*/
    nout = map_v2s_results(results, im3d, map, debug, go->sopt.dnode);

    /* before free()ing results, check whether we want to return the values */
    if ( Rdata )
    {
        *Rdata = results->vals[0];
        results->vals[0] = NULL;   /* do not let free_v2s_results() free it */
    }

    free_v2s_results(results);

    if ( debug > 1 ) fprintf(stderr,"++ map_v2s_results: nout = %d\n", nout);

    /* check failure */
    if ( ! *map ) RETURN(-1);

    RETURN(nout);
}

/*! given vol2surf results, fill SUMA_irgba struct */
static int map_v2s_results(v2s_results *res, Three_D_View *im3d,
                           SUMA_irgba **map, int debug, int dnode )
{
    SUMA_irgba * mptr;
    MCW_pbar   * pbar;
    rgbyte     * cmap;
    float      * result_vals;
    float        bbot, btop, fac, fval;
    float        pane_scale = 1.0;
    byte         r, g, b;
    int          nindex, node, npanes, ival;
    int          map_index = 0;  /* counter for results */

ENTRY("map_v2s_results");

    /*-------------------- create output node list --------------------*/
    /* from malloc    12 Feb 2009 [lesstif patrol] */
    mptr = (SUMA_irgba *) calloc(res->nused, sizeof(SUMA_irgba));
    if ( ! mptr )
    {
        fprintf(stderr,"** malloc failure: %d irgba's\n", res->nused);
        RETURN(-1);
    }
    *map = mptr;  /* and store the address */

    pbar       = im3d->vwid->func->inten_pbar;
    npanes     = pbar->num_panes;

    /* with AFNI_PBAR_FULLRANGE, pane_scale = 1.0  [20 May 2019 rickr] */
    /*
       pane_scale = im3d->vinfo->fim_range;
       if ( pane_scale == 0.0 ) pane_scale = im3d->vinfo->fim_autorange;
       if ( pane_scale == 0.0 ) pane_scale = 1.0;
    */
    if( PBAR_FULLRANGE ) {
       pane_scale = 1.0;
       if( debug > 2 )
          fprintf(stderr,"+d mvr: have PBAR_FULLRANGE, resetting pane_scale\n");
    } else {
       pane_scale = FIM_RANGE(im3d);
       if( debug > 2 )
          fprintf(stderr,"+d mvr: no PBAR_FULLRANGE, pane_scale = FIM_RANGE\n");
    }

    if ( debug > 1 ) {
        fprintf(stderr,"+d mvr: npanes = %d, pane_scale = %f\n",
                pbar->bigmode ? NPANE_BIG : npanes, pane_scale);
        fprintf(stderr,"   fim_range = %f, fim_autorange = %f\n",
                im3d->vinfo->fim_range, im3d->vinfo->fim_autorange);
    }

    result_vals = res->vals[0]; /* for typing and potential speed */

    if( pbar->bigmode )         /* colorscale */
    {
        int zbot;

        bbot = pane_scale*pbar->bigbot;
        btop = pane_scale*pbar->bigtop;
        fac  = NPANE_BIG / (btop-bbot);
        cmap = pbar->bigcolor;
        zbot = (bbot == 0.0);

        if ( debug > 1 ) {
            fprintf(stderr,"+d bigmode: bbot,btop,fac, zbot = %f,%f,%f, %d\n",
                    bbot, btop, fac, zbot);
            fprintf(stderr,"   NPANE_BIG = %d, BIGGEST = %d\n",
                            NPANE_BIG, NPANE_BIGGEST);
        }

        for ( nindex = 0; nindex < res->nused; nindex++ )
        {
            if ( res->nodes ) node = res->nodes[nindex];
            else              node = nindex;

            /* note value, and check to ignore */
            fval = result_vals[nindex];

            if ( debug > 1 && node == dnode )
                fprintf(stderr, "+d dnode %d, fval %f\n", dnode, fval);

            if ( zbot && fval < 0 ) continue;

            /* note the color panel index, and bound it in [0,NPANE_BIG-1] */
            if( fval >= btop ) ival = 0;        /* guard against overflow */
            else if ( fval <= bbot ) ival = NPANE_BIG - 1;
            /* trunc: matching afni_func.c    [20 May 2019 rickr] */
            else ival = (int)(fac * (btop - fval));

            if ( ival < 0 ) ival = 0;
            if ( ival >= NPANE_BIG ) ival = NPANE_BIG - 1;

            /* get color map, and check if non-zero */
            r = cmap[ival].r;  g = cmap[ival].g;  b = cmap[ival].b;

            if ( debug > 1 && node == dnode )
                fprintf(stderr, "+d pane, r,g,b : %d, %d,%d,%d\n",ival,r,g,b);

            if ( r == 0 && g == 0 && b == 0 ) continue;

            /* we are set, fill the SUMA_irgba struct  */
            mptr->id = res->nodes[nindex];
            mptr->r = r;  mptr->g = g;  mptr->b = b;  mptr->a = 255;

            /* increment pointer and counter (okay, so counter is unneeded) */
            mptr++;  map_index++;
        }
    }
    else        /* indexed colors */
    {
        float othr[NPANE_MAX+1];        /* threshold */
        short ovc[NPANE_MAX+1]; /* color */
        byte  ovc_r[NPANE_MAX+1], ovc_g[NPANE_MAX+1], ovc_b[NPANE_MAX+1];

        /* set the overlay color indices */
        for( ival=0 ; ival < npanes ; ival++ )
            ovc[ival] = pbar->ov_index[ival];   /* from top of pbar down */
        /* npanes->ival to remove compile warning    [12 Jan 2022 rickr] */
        ovc[ival] = im3d->vinfo->use_posfunc ? 0 : ovc[ival-1];

        /* get the actual RGB colors of each pane on the pbar */
        for( ival=0 ; ival <= npanes ; ival++ ) /* include npanes 27 Aug 2007 */
        {
            ovc_r[ival] = DCOV_REDBYTE  (im3d->dc,ovc[ival]);
            ovc_g[ival] = DCOV_GREENBYTE(im3d->dc,ovc[ival]);
            ovc_b[ival] = DCOV_BLUEBYTE (im3d->dc,ovc[ival]);
        }

        /* compute the thresholds */
        for( ival=0 ; ival < npanes ; ival++ )
            othr[ival] = pane_scale * pbar->pval[ival+1];
        /* npanes->ival to remove compile warning    [12 Jan 2022 rickr] */
        othr[ival] = othr[ival-1] ;

        if ( debug > 2 )
        {
            for( ival=0 ; ival <= npanes ; ival++ )
                fprintf(stderr,"+d pane #%2d, ovc = %d, othr = %f\n",
                        ival, ovc[ival], othr[ival]);
        }

        for ( nindex = 0; nindex < res->nused; nindex++ )
        {
            if ( res->nodes ) node = res->nodes[nindex];
            else              node = nindex;

            /* note value, and check to ignore */
            fval = result_vals[nindex];

            if ( debug > 1 && node == dnode )
                fprintf(stderr, "+d dnode %d, fval %f\n", dnode, fval);

            if ( fval == 0.0 ) continue;

            for ( ival = 0; ival < npanes && fval < othr[ival]; ival++ )
                ;
            if ( ovc[ival] == 0 ) continue;  /* no color in this pane */
            r = ovc_r[ival];  g = ovc_g[ival];  b = ovc_b[ival]; 

            /* we are set, fill the SUMA_irgba struct  */
            mptr->id = res->nodes[nindex];
            mptr->r = r;  mptr->g = g;  mptr->b = b;  mptr->a = 255;

            if ( debug > 1 && node == dnode )
                fprintf(stderr, "+d pane, r,g,b = %d, %d,%d,%d\n",ival,r,g,b);

            if ( r == 0 && g == 0 && b == 0 ) continue;

            /* increment pointer and counter (okay, so counter is unneeded) */
            mptr++;  map_index++;
        }
   }

    if ( debug > 0 )
        fprintf(stderr,"+d mvr v2s map size = %d\n", map_index);

   /* forfeit unused memory */
   *map = (SUMA_irgba *)realloc(*map, sizeof(SUMA_irgba)*map_index);
   if ( ! *map )
   {
        fprintf(stderr,"** mvr: failed realloc of map, %d irgba's\n",map_index);
        map_index = -1;
   }

   RETURN(map_index) ;  /* number of entries in map */
}


/*-----------------------------------------------------------------------*/
/*! Return the vol2surf time series data at a single node.
 *
 *      sess    : session to get surfaces from (e.g. im3d->ss_now)
 *      dset    : dataset to get values from (e.g. im3d->anat_now)
 *      surf    : surface index (will use LDP of this surface)
 *                              (e.g. in {0..from im3d->ss_now->su_num-1})
 *      node    : node index in surface
 *
 *  Note: surf/node might come from AFNI_get_xhair_node(), say.
 *
 *  returns: a pointer to the float data or NULL on failure.
 *
 *  Array will be of length DSET_NVALS(dset), and should be free'd by
 *  the calling function.
 * ----------------------------------------------------------------------*/
float * AFNI_v2s_node_timeseries(THD_session * sess, THD_3dim_dataset * dset,
                                 int surfA, int surfB, int node, int use_v2s)
{
   v2s_results     * results = NULL;
   v2s_opts_t        sopt;
   SUMA_surface    * sA, * sB;
   float           * ts = NULL;  /* time series to return */
   int               ind, len, verb = gv2s_plug_opts.sopt.debug;

   ENTRY("AFNI_v2s_node_timeseries");

   if( !sess || !dset || surfA < 0 ){
      fprintf(stderr,"** v2s_node_ts - bad inputs: %p,%p,%d\n",sess,dset,surfA);
      RETURN(NULL);
   }

   sA = (surfA >= 0) ? sess->su_surf[surfA] : NULL;
   sB = (surfB >= 0) ? sess->su_surf[surfB] : NULL;

   /* check and apply the node index */
   if( node < 0 || node >= sA->num_ixyz ) {
      fprintf(stderr,"** v2s_NTS: node %d outside [0,%d)\n",node,sA->num_ixyz);
      RETURN(NULL);
   } else if( verb > 1 )
      fprintf(stderr,"-- v2s_NTS: getting time series at node %d\n", node);

   /* fill options struct (from defaults or plugin opts) */
   if( use_v2s ) {
      /* is using plugin, nuke any gp_index, filenames or structs */
      sopt = gv2s_plug_opts.sopt;

      sopt.gp_index = -1;  /* get all sub-bricks */
      sopt.outfile_1D = sopt.outfile_niml = sopt.segc_file = NULL;
      memset(&sopt.cmd, 0, sizeof(v2s_cmd_t));
      memset(&sopt.oob, 0, sizeof(v2s_oob_t));
      memset(&sopt.oom, 0, sizeof(v2s_oob_t));
   } else
      v2s_fill_sopt_default(&sopt, sB ? 2 : 1);

   sopt.skip_cols = V2S_SKIP_ALL ^ V2S_SKIP_VALS;  /* get all data */
   sopt.first_node = sopt.last_node = node;        /* get only 1 node */
   if( verb > 2 ) sopt.dnode = node;               /* maybe make verbose */

   /* get the results */
   results = opt_vol2surf(dset, &sopt, sA, sB, NULL);
   if( !results ) RETURN(NULL);

   /* extract single time series */
   len = DSET_NVALS(dset);
   if( len != results->nlab ){
      fprintf(stderr,"** v2s_NTS: nvals != nlab (%d, %d)\n",len,results->nlab);
      free_v2s_results(results); RETURN(NULL);
   }

   ts = (float *)malloc(len*sizeof(float));
   if( !ts ) {
      fprintf(stderr,"** v2s_NTS: failed to alloc %d floats\n", len);
      free_v2s_results(results); RETURN(NULL);
   }

   /* actually fill the array, woohoo! */
   if(verb > 3) fprintf(stderr,"-- v2s_NTS: copying %d vals for return\n",len);
   for( ind = 0; ind < len; ind++ )
      ts[ind] = results->vals[ind][0];

   if(verb > 3) fprintf(stderr,"-- v2s_NTS: freeing results... ");
   free_v2s_results(results);

   if(verb > 3) fprintf(stderr,"done\n");
   RETURN(ts);
}

