00001 /*----------------------------------------------------------------------
00002  * afni_vol2surf.c   - afni's interface to the vol2surf library
00003  *
00004  * AFNI_vol2surf_func_overlay() - similar to AFNI_vnlist_func_overlay()
00005  *----------------------------------------------------------------------
00006  */
00007 
00008 /*----------------------------------------------------------------------
00009  * R. Reynolds    September, 2004
00010  *
00011  * history:
00012  *
00013  * 08 Oct 2004 [rickr]:
00014  *   - AFNI_vol2surf_func_overlay() has new params, surfA,surfB,use_defaults
00015  *   - pass use_defaults to afni_vol2surf()
00016  *
00017  * 25 Oct 2004 [rickr]:
00018  *   - accept Rdata and Rthr pointers, for optionally returning the data
00019  *     and global threshold
00020  *   - make threshold masks absolute
00021  *----------------------------------------------------------------------
00022  */
00023 
00024 #include "afni.h"
00025 #include "vol2surf.h"
00026 
00027 static int map_v2s_results(v2s_results *res, Three_D_View *im3d,
00028                            SUMA_irgba **map, int debug, int dnode );
00029 /*-----------------------------------------------------------------------*/
00030 /*! Create a nodal color overlay from a voxel map.
00031     - Surface index(es) come from global v2s_plugin_opts struct
00032     - Return value is number of nodes overlaid:
00033       -  0 return ==> no overlay was computed
00034       - -1 return ==> some error (e.g., no surface nodes on this dataset)
00035       = positive  ==> can use *map
00036     - *map is set to a newly malloc()-ed array (if return > 0)
00037     - if (Rdata) return data from results structure
00038     - if (Rthr) return overlay threshold from afni
00039 
00040     Sample usage:
00041     - SUMA_irgba * map;
00042     - int          nmap;
00043     - nmap = AFNI_vnlist_func_overlay( im3d, &map );
00044     -      if( nmap <  0 ){ ** error ** }
00045     - else if( nmap == 0 ){ ** nothing to show ** }
00046     - else                { ** show map[0..nmap-1] ** }
00047 
00048     * based on AFNI_vnlist_func_overlay()
00049 -------------------------------------------------------------------------*/
00050 int AFNI_vol2surf_func_overlay(Three_D_View *im3d, SUMA_irgba **map,
00051          int surfA, int surfB, int use_defaults, float ** Rdata, float * Rthr )
00052 {
00053     THD_3dim_dataset * oset;            /* overlay dataset */
00054     THD_session      * ss;
00055     SUMA_surface     * sA, * sB;
00056     MRI_IMAGE        * im_oim;
00057     v2s_plugin_opts  * go;
00058     v2s_results      * results;
00059     byte             * cmask;
00060     int                nout, oind, debug;
00061 
00062 ENTRY("AFNI_vol2surf_func_overlay") ;
00063 
00064     /* check inputs */
00065     if ( map == NULL || !IM3D_VALID(im3d) ) RETURN(-1);
00066 
00067     go = &gv2s_plug_opts;
00068     if ( ! use_defaults && ! go->ready ) RETURN(-1);
00069 
00070     debug = go->sopt.debug;         /* because I'm lazy */
00071 
00072     ss = im3d->ss_now;              /* session must have needed surface(s) */
00073     if( ss                 == NULL      ||
00074         ss->su_num         <= 0         ||
00075         surfA              <  0         ||
00076         ss->su_num         <= surfA     ||
00077         ss->su_num         <= surfB     ||
00078         ss->su_surf[surfA] == NULL      ||
00079         (surfB >= 0 && ss->su_surf[surfB] == NULL) )
00080     {
00081         if ( debug > 1 )
00082         {
00083             if ( !ss ) fprintf(stderr,"** v2s: NULL session\n");
00084             else
00085                 fprintf(stderr,"** v2s: bad session data:\n"
00086                         "   su_num, surfA, surfB = %d, %d, %d\n",
00087                         ss->su_num, surfA, surfB);
00088         }
00089         RETURN(-1);
00090     }
00091 
00092     /* init return values */
00093     if ( Rdata ) *Rdata = NULL;
00094     if ( Rthr )  *Rthr  = 0.0;
00095 
00096     /* set surface pointers */
00097     sA = ss->su_surf[surfA];
00098     sB = ( surfB >= 0 ) ? ss->su_surf[surfB] : NULL;
00099 
00100     if ( debug )
00101     {
00102         fprintf(stderr,"++ v2s overlay: sa,sb = %d,%d\n", surfA, surfB);
00103         if ( debug > 1 )
00104             fprintf(stderr,"  surfA is %s, surfB is %s\n",
00105                 sA->label[0] ? sA->label : "<no label>",
00106                 sB ? (sB->label[0] ? sB->label : "<no label>") : "<not used>");
00107     }
00108 
00109     /*-------------------- overlay image --------------------*/
00110     oset = im3d->fim_now;
00111     oind = im3d->vinfo->fim_index;  /* overlay sub-brick index */
00112 
00113     if( oset == NULL )
00114     {
00115         if ( debug > 1 ) fprintf(stderr,"** v2s: no overlay dset\n");
00116         RETURN(-1);
00117     }
00118     im_oim = DSET_BRICK(oset,oind);
00119 
00120     if( ! im_oim || ! AFNI_GOOD_FUNC_DTYPE(im_oim->kind) )
00121     {
00122         if ( debug > 1 ) fprintf(stderr,"** v2s: no overlay image\n");
00123         RETURN(-1);
00124     }
00125     DSET_load(oset);                    /* to be sure */
00126     if( !DSET_LOADED(oset) ) RETURN(-1);
00127 
00128     /*-------------------- mask from threshold --------------------*/
00129     cmask = NULL;
00130     if( im3d->vinfo->func_threshold > 0.0 )     /* then want a threshold */
00131     {
00132         MRI_IMAGE * im_thr;
00133         float       thresh;
00134         int         tind = im3d->vinfo->thr_index;
00135 
00136         im_thr = DSET_BRICK(oset,tind);
00137 
00138         /* note real threshold */
00139         thresh = im3d->vinfo->func_threshold * im3d->vinfo->func_thresh_top;
00140 
00141         /* maybe we want to return this */
00142         if ( Rthr ) *Rthr = thresh;
00143 
00144         if( im_thr && !AFNI_GOOD_FUNC_DTYPE(im_thr->kind) )
00145             im_thr = NULL;
00146 
00147         /* create the mask, if this fails, just continue... */
00148         if( im_thr != NULL )
00149         {
00150             int nset;
00151             if ( debug > 1 )
00152                 fprintf(stderr,"++ mask from index %d and thresh %f\n",
00153                         tind,thresh);
00154             nset = thd_mask_from_brick(oset, tind, thresh, &cmask, 1);
00155             if ( debug > 1 )
00156             {
00157                 if ( ! cmask )
00158                     fprintf(stderr,"-- no mask created\n");
00159                 else
00160                     fprintf(stderr,"++ mask has %d set voxels\n", nset);
00161             }
00162         }
00163         else if ( debug > 1 )
00164             fprintf(stderr,"-- no threshold mask\n");
00165     }
00166 
00167     /*-------------------- vol2surf computation --------------------*/
00168     results = afni_vol2surf(oset, oind, sA, sB, cmask, use_defaults);
00169 
00170     if ( cmask ) free(cmask);   /* we're done with the mask */
00171     if ( ! results )
00172     {
00173         if ( debug ) fprintf(stderr,"-- vol2surf failure\n");
00174         RETURN(-1);             /* failure? */
00175     }
00176 
00177     /*-------------------- lose old vnlist --------------------*/
00178     /* vol2surf has no method for counting voxels (right now)  */
00179     if( sA->vn )
00180     {
00181         if ( debug > 1 ) fprintf(stderr,"-- destroying vnlist\n");
00182         SUMA_destroy_vnlist( sA->vn ) ;
00183         sA->vn = NULL;
00184     }
00185 
00186     /*-------------------- set overlay colors --------------------*/
00187     nout = map_v2s_results(results, im3d, map, debug, go->sopt.dnode);
00188 
00189     /* before free()ing results, check whether we want to return the values */
00190     if ( Rdata )
00191     {
00192         *Rdata = results->vals[0];
00193         results->vals[0] = NULL;   /* do not let free_v2s_results() free it */
00194     }
00195 
00196     free_v2s_results(results);
00197 
00198     if ( debug > 1 ) fprintf(stderr,"++ map_v2s_results: nout = %d\n", nout);
00199 
00200     /* check failure */
00201     if ( ! *map ) RETURN(-1);
00202 
00203     RETURN(nout);
00204 }
00205 
00206 /*! given vol2surf results, fill SUMA_irgba struct */
00207 static int map_v2s_results(v2s_results *res, Three_D_View *im3d,
00208                            SUMA_irgba **map, int debug, int dnode )
00209 {
00210     SUMA_irgba * mptr;
00211     MCW_pbar   * pbar;
00212     rgbyte     * cmap;
00213     float      * result_vals;
00214     float        bbot, btop, fac, fval;
00215     float        pane_scale = 1.0;
00216     byte         r, g, b;
00217     int          nindex, node, npanes, ival;
00218     int          map_index = 0;  /* counter for results */
00219 
00220 ENTRY("map_v2s_results");
00221 
00222     /*-------------------- create output node list --------------------*/
00223     mptr = (SUMA_irgba *) malloc(sizeof(SUMA_irgba) * res->nused);
00224     if ( ! mptr )
00225     {
00226         fprintf(stderr,"** malloc failure: %d irgba's\n", res->nused);
00227         RETURN(-1);
00228     }
00229     *map = mptr;  /* and store the address */
00230 
00231     pbar       = im3d->vwid->func->inten_pbar;
00232     npanes     = pbar->num_panes;
00233     pane_scale = im3d->vinfo->fim_range;
00234     if ( pane_scale == 0.0 ) pane_scale = im3d->vinfo->fim_autorange;
00235     if ( pane_scale == 0.0 ) pane_scale = 1.0;
00236 
00237     if ( debug > 1 )
00238         fprintf(stderr,"+d mvr: npanes = %d, pane_scale = %f\n",
00239                 pbar->bigmode ? NPANE_BIG : npanes, pane_scale);
00240 
00241     result_vals = res->vals[0]; /* for typing and potential speed */
00242 
00243     if( pbar->bigmode )         /* colorscale */
00244     {
00245         int zbot;
00246 
00247         bbot = pane_scale*pbar->bigbot;
00248         btop = pane_scale*pbar->bigtop;
00249         fac  = NPANE_BIG / (btop-bbot);
00250         cmap = pbar->bigcolor;
00251         zbot = (bbot == 0.0);
00252 
00253         if ( debug > 1 )
00254             fprintf(stderr,"+d bigmode: bbot,btop,fac, zbot = %f,%f,%f, %d\n",
00255                     bbot, btop, fac, zbot);
00256 
00257         for ( nindex = 0; nindex < res->nused; nindex++ )
00258         {
00259             if ( res->nodes ) node = res->nodes[nindex];
00260             else              node = nindex;
00261 
00262             /* note value, and check to ignore */
00263             fval = result_vals[nindex];
00264 
00265             if ( debug > 1 && node == dnode )
00266                 fprintf(stderr, "+d dnode %d, fval %f\n", dnode, fval);
00267 
00268             if ( zbot && fval < 0 ) continue;
00269 
00270             /* note the color panel index, and bound it in [0,NPANE_BIG-1] */
00271             ival = (int)(fac * (btop - fval) + 0.49);
00272             if ( ival < 0 ) ival = 0;
00273             if ( ival >= NPANE_BIG ) ival = NPANE_BIG - 1;
00274 
00275             /* get color map, and check if non-zero */
00276             r = cmap[ival].r;  g = cmap[ival].g;  b = cmap[ival].b;
00277 
00278             if ( debug > 1 && node == dnode )
00279                 fprintf(stderr, "+d pane, r,g,b = %d, %d,%d,%d\n",ival,r,g,b);
00280 
00281             if ( r == 0 && g == 0 && b == 0 ) continue;
00282 
00283             /* we are set, fill the SUMA_irgba struct  */
00284             mptr->id = res->nodes[nindex];
00285             mptr->r = r;  mptr->g = g;  mptr->b = b;  mptr->a = 255;
00286 
00287             /* increment pointer and counter (okay, so counter is unneeded) */
00288             mptr++;  map_index++;
00289         }
00290     }
00291     else        /* indexed colors */
00292     {
00293         float othr[NPANE_MAX];  /* threshold */
00294         short ovc[NPANE_MAX+1]; /* color */
00295         byte  ovc_r[NPANE_MAX+1], ovc_g[NPANE_MAX+1], ovc_b[NPANE_MAX+1];
00296 
00297         /* set the overlay color indices */
00298         for( ival=0 ; ival < npanes ; ival++ )
00299             ovc[ival] = pbar->ov_index[ival];   /* from top of pbar down */
00300         ovc[npanes] = im3d->vinfo->use_posfunc ? 0 : ovc[npanes-1];
00301 
00302         /* get the actual RGB colors of each pane on the pbar */
00303         for( ival=0 ; ival < npanes ; ival++ )
00304         {
00305             ovc_r[ival] = DCOV_REDBYTE  (im3d->dc,ovc[ival]);
00306             ovc_g[ival] = DCOV_GREENBYTE(im3d->dc,ovc[ival]);
00307             ovc_b[ival] = DCOV_BLUEBYTE (im3d->dc,ovc[ival]);
00308         }
00309 
00310         /* compute the thresholds */
00311         for( ival=0 ; ival < npanes ; ival++ )
00312             othr[ival] = pane_scale * pbar->pval[ival+1];
00313 
00314         if ( debug > 2 )
00315         {
00316             for( ival=0 ; ival <= npanes ; ival++ )
00317                 fprintf(stderr,"+d pane #%2d, ovc = %d, othr = %f\n",
00318                         ival, ovc[ival], othr[ival]);
00319         }
00320 
00321         for ( nindex = 0; nindex < res->nused; nindex++ )
00322         {
00323             if ( res->nodes ) node = res->nodes[nindex];
00324             else              node = nindex;
00325 
00326             /* note value, and check to ignore */
00327             fval = result_vals[nindex];
00328 
00329             if ( debug > 1 && node == dnode )
00330                 fprintf(stderr, "+d dnode %d, fval %f\n", dnode, fval);
00331 
00332             if ( fval == 0.0 ) continue;
00333 
00334             for ( ival = 0; ival < npanes && fval < othr[ival]; ival++ )
00335                 ;
00336             if ( ovc[ival] == 0 ) continue;  /* no color in this pane */
00337             r = ovc_r[ival];  g = ovc_g[ival];  b = ovc_b[ival]; 
00338 
00339             /* we are set, fill the SUMA_irgba struct  */
00340             mptr->id = res->nodes[nindex];
00341             mptr->r = r;  mptr->g = g;  mptr->b = b;  mptr->a = 255;
00342 
00343             if ( debug > 1 && node == dnode )
00344                 fprintf(stderr, "+d pane, r,g,b = %d, %d,%d,%d\n",ival,r,g,b);
00345 
00346             if ( r == 0 && g == 0 && b == 0 ) continue;
00347 
00348             /* increment pointer and counter (okay, so counter is unneeded) */
00349             mptr++;  map_index++;
00350         }
00351    }
00352 
00353     if ( debug > 0 )
00354         fprintf(stderr,"+d mvr v2s map size = %d\n", map_index);
00355 
00356    /* forfeit unused memory */
00357    *map = (SUMA_irgba *)realloc(*map, sizeof(SUMA_irgba)*map_index);
00358    if ( ! *map )
00359    {
00360         fprintf(stderr,"** mvr: failed realloc of map, %d irgba's\n",map_index);
00361         map_index = -1;
00362    }
00363 
00364    RETURN(map_index) ;  /* number of entries in map */
00365 }
00366