1 /*****************************************************************************
2    Major portions of this software are copyrighted by the Medical College
3    of Wisconsin, 1994-2000, and are released under the Gnu General Public
4    License, Version 2.  See the file README.Copyright for details.
5 ******************************************************************************/
6 
7 #include "afni.h"
8 
9 #ifndef ALLOW_PLUGINS
10 #  error "Plugins not properly set up -- see machdep.h"
11 #endif
12 
13 /***********************************************************************
14   Plugin to plot a scatterplot of 2 bricks.
15   [Adapted from plug_histog.c - RWCox - 13 Jan 2000]
16   11 Aug 2001: modified to print correlation coefficient
17 ************************************************************************/
18 
19 char * SCAT_main( PLUGIN_interface * ) ;
20 
21 static char helpstring[] =
22    " \n"
23    " Purpose:  Scatterplot data from 2 bricks.\n"
24    "\n"
25    " Source x: Dataset   = dataset for x-axis values\n"
26    "           Sub-brick = which one to use\n"
27    "           Bottom    = minimum value to include\n"
28    "           Top       = maximum value to include\n"
29    "\n"
30    " Source y: Dataset   = dataset for y-axis values\n"
31    "           Sub-brick = which one to use\n"
32    "           Bottom    = minimum value to include\n"
33    "           Top       = maximum value to include\n"
34    "\n"
35    " Mask:     Dataset   = masking dataset\n"
36    "           Sub-brick = which one to use\n"
37    "           Bottom    = min value from mask dataset to use\n"
38    "           Top       = max value from mask dataset to use\n"
39    "\n"
40    "In the above definitions:\n"
41    "  if Bottom >  Top, then all voxels will be used\n"
42    "  if Bottom <= Top, then only voxels in this range will be used\n"
43    "\n"
44    "Note: A maximum of 4,000,000 (four million) points can be plotted.\n\n"
45    " Author -- RW Cox - 13 January 2000\n"
46 ;
47 
48 /***********************************************************************
49    Set up the interface to the user
50 ************************************************************************/
51 
52 
53 DEFINE_PLUGIN_PROTOTYPE
54 
PLUGIN_init(int ncall)55 PLUGIN_interface * PLUGIN_init( int ncall )
56 {
57    PLUGIN_interface * plint ;
58 
59    if( ncall > 0 ) return NULL ;  /* only one interface */
60 
61    /*-- set titles and call point --*/
62 
63    plint = PLUTO_new_interface( "ScatterPlot" ,
64                                 "ScatterPlot of 2 Bricks" ,
65                                 helpstring ,
66                                 PLUGIN_CALL_VIA_MENU , SCAT_main  ) ;
67 
68    PLUTO_add_hint( plint , "Of 2 Bricks" ) ;
69 
70    PLUTO_set_sequence( plint , "A:afniinfo:dset" ) ;
71 
72    PLUTO_set_runlabels( plint , "Plot+Keep" , "Plot+Close" ) ;  /* 04 Nov 2003 */
73 
74    /*-- first line of input --*/
75 
76    PLUTO_add_option( plint , "Source x" , "Source x" , TRUE ) ;
77    PLUTO_add_dataset(plint , "Dataset" ,
78                                     ANAT_ALL_MASK , FUNC_ALL_MASK ,
79                                     DIMEN_ALL_MASK | BRICK_ALLREAL_MASK ) ;
80    PLUTO_add_number( plint , "Sub-brick" , 0,9999,0 , 0,1 ) ;
81    PLUTO_add_number( plint , "Bottom" , -99999,99999, 0,  1,1 ) ;
82    PLUTO_add_number( plint , "Top"    , -99999,99999, 0, -1,1 ) ;
83 
84    /*-- second line of input --*/
85 
86    PLUTO_add_option( plint , "Source y" , "Source y" , TRUE ) ;
87    PLUTO_add_dataset(plint , "Dataset" ,
88                                     ANAT_ALL_MASK , FUNC_ALL_MASK ,
89                                     DIMEN_ALL_MASK | BRICK_ALLREAL_MASK ) ;
90    PLUTO_add_number( plint , "Sub-brick" , 0,9999,0 , 0,1 ) ;
91    PLUTO_add_number( plint , "Bottom" , -99999,99999, 0,  1,1 ) ;
92    PLUTO_add_number( plint , "Top"    , -99999,99999, 0, -1,1 ) ;
93 
94    /*-- third line of input --*/
95 
96    PLUTO_add_option( plint , "Mask" , "Mask" , FALSE ) ;
97    PLUTO_add_dataset( plint , "Dataset" ,
98                                     ANAT_ALL_MASK , FUNC_ALL_MASK ,
99                                     DIMEN_ALL_MASK | BRICK_ALLREAL_MASK ) ;
100    PLUTO_add_number( plint , "Sub-brick" , 0,9999,0 , 0,1 ) ;
101    PLUTO_add_number( plint , "Bottom" , -99999,99999, 0,  1,1 ) ;
102    PLUTO_add_number( plint , "Top"    , -99999,99999, 0, -1,1 ) ;
103 
104    /*-- fourth line of input [30 Oct 2006] --*/
105 
106    PLUTO_add_option( plint , "Aboot" , "Aboot" , FALSE ) ;
107    PLUTO_add_number( plint , "Radius" , 2,100,0,10,1 ) ;
108 
109    /*-- fifth and sixth lines of input [08 Aug 2007] --*/
110 
111    PLUTO_add_option( plint , "Percent" , "Percent" , FALSE ) ;
112    PLUTO_add_number( plint , "Percent" , 1,100,0,10,1 ) ;
113 
114    PLUTO_add_option( plint , "Label" , "Label" , FALSE ) ;
115    PLUTO_add_string( plint , "Label" , 0,NULL , 16 ) ;
116 
117    return plint ;
118 }
119 
120 /***************************************************************************
121   Main routine for this plugin (will be called from AFNI).
122 ****************************************************************************/
123 
SCAT_main(PLUGIN_interface * plint)124 char * SCAT_main( PLUGIN_interface *plint )
125 {
126    MCW_idcode *idc ;
127    THD_3dim_dataset *xdset, *ydset , *mask_dset=NULL ;
128    int ivx,ivy , mcount , nvox , ii,jj , nbin=-1 , did_corr=0 , zzc=0 ;
129    float mask_bot=666.0 , mask_top=-666.0 ;
130    float xbot,xtop, ybot,ytop, pcor=0.0f, a=0.0f,b=0.0f, p025=0.0f,p975=0.0f, mi=0.0f ;
131    char *tag , *str ;
132    char xlab[THD_MAX_NAME],ylab[THD_MAX_NAME],tlab[THD_MAX_NAME+64] ;
133    char ab[16] , bb[16] , *pab,*pbb ;
134    byte *mmm ;
135    float *xar , *yar ;
136    float scor=0.0f,s025=0.0f,s975=0.0f , sa=0.0f,sb=0.0f ;  /* 23 Dec 2014 */
137    float alin[2] , blin[2] ; float_triple clin[2] ;         /* 24 Dec 2014 */
138    char *eplot=NULL ;
139 
140    char *cname=NULL ;  /* 06 Aug 1998 */
141    int miv=0 ;
142 
143    float hrad=0.0f ;   /* 30 Oct 2006 */
144    float perc=0.0f ;   /* 08 Aug 2007 */
145    char *label=NULL ;  /* 09 Aug 2007 */
146 
147    /*--------------------------------------------------------------------*/
148    /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
149 
150    if( plint == NULL )
151       return "**********************\n"
152              "SCAT_main:  NULL input\n"
153              "**********************"  ;
154 
155    /*-- read x dataset line --*/
156 
157    PLUTO_next_option(plint) ;
158    idc   = PLUTO_get_idcode(plint) ;
159    xdset = PLUTO_find_dset(idc) ;
160    if( xdset == NULL )
161       return "*************************\n"
162              "SCAT_main:  bad x dataset\n"
163              "*************************"  ;
164 
165    ivx = (int) PLUTO_get_number(plint) ;
166    if( ivx >= DSET_NVALS(xdset) || ivx < 0 )
167       return "***************************\n"
168              "SCAT_main:  bad x sub-brick\n"
169              "***************************" ;
170 
171    DSET_load(xdset) ;
172    if( DSET_ARRAY(xdset,ivx) == NULL )
173       return "********************************\n"
174              "SCAT_main:  can't load x dataset\n"
175              "********************************"  ;
176    nvox = DSET_NVOX(xdset) ;
177 
178    xbot = PLUTO_get_number(plint) ;
179    xtop = PLUTO_get_number(plint) ;
180 
181    /*-- read y dataset line --*/
182 
183    PLUTO_next_option(plint) ;
184    idc   = PLUTO_get_idcode(plint) ;
185    ydset = PLUTO_find_dset(idc) ;
186    if( ydset == NULL )
187       return "*************************\n"
188              "SCAT_main:  bad y dataset\n"
189              "*************************"  ;
190 
191    ivy = (int) PLUTO_get_number(plint) ;
192    if( ivy >= DSET_NVALS(ydset) || ivy < 0 )
193       return "***************************\n"
194              "SCAT_main:  bad y sub-brick\n"
195              "***************************" ;
196 
197    if( DSET_NVOX(ydset) != nvox )
198       return "************************************************\n"
199              "SCAT_main:  x and y datasets don't match in size\n"
200              "************************************************" ;
201 
202    DSET_load(ydset) ;
203    if( DSET_ARRAY(ydset,ivy) == NULL )
204       return "********************************\n"
205              "SCAT_main:  can't load y dataset\n"
206              "********************************"  ;
207 
208    ybot = PLUTO_get_number(plint) ;
209    ytop = PLUTO_get_number(plint) ;
210 
211    /*-- read optional line(s) --*/
212 
213    while( (tag=PLUTO_get_optiontag(plint)) != NULL ){
214 
215       /*-- Mask dataset --*/
216 
217       if( strcmp(tag,"Mask") == 0 ){
218 
219          idc       = PLUTO_get_idcode(plint) ;
220          mask_dset = PLUTO_find_dset(idc) ;
221 
222          if( mask_dset == NULL )
223             return "****************************\n"
224                    "SCAT_main:  bad mask dataset\n"
225                    "****************************"  ;
226 
227          if( DSET_NVOX(mask_dset) != nvox )
228            return "**********************************************************\n"
229                   "SCAT_main: mask input dataset doesn't match source dataset\n"
230                   "**********************************************************" ;
231 
232          miv = (int) PLUTO_get_number(plint) ;  /* 06 Aug 1998 */
233          if( miv >= DSET_NVALS(mask_dset) || miv < 0 )
234             return "**************************************************\n"
235                    "SCAT_main: mask dataset sub-brick index is illegal\n"
236                    "**************************************************"  ;
237 
238          DSET_load(mask_dset) ;
239          if( DSET_ARRAY(mask_dset,miv) == NULL )
240             return "***********************************\n"
241                    "SCAT_main:  can't load mask dataset\n"
242                    "***********************************"  ;
243 
244          mask_bot = PLUTO_get_number(plint) ;
245          mask_top = PLUTO_get_number(plint) ;
246          continue ;
247       }
248 
249       /*-- 30 Oct 2006: Aboot --*/
250 
251       if( strcmp(tag,"Aboot") == 0 ){
252          hrad = PLUTO_get_number(plint) ;
253          continue ;
254       }
255 
256       /*-- 08 Aug 2007: Percent --*/
257 
258       if( strcmp(tag,"Percent") == 0 ){
259          perc = PLUTO_get_number(plint) ;
260          continue ;
261       }
262 
263       /*-- 09 Aug 2007: Label --*/
264 
265       if( strcmp(tag,"Label") == 0 ){
266         char *ll = PLUTO_get_string(plint) ;
267         if( ll != NULL ) label = strdup(ll) ;
268         continue ;
269       }
270    }
271 
272    /*------------------------------------------------------*/
273    /*---------- At this point, the inputs are OK ----------*/
274 
275    /*-- build the byte mask array --*/
276 
277    if( mask_dset == NULL ){
278       mmm = (byte *) malloc( sizeof(byte) * nvox ) ;
279       if( mmm == NULL )
280          return " \n*** Can't malloc workspace for mask! ***\n" ;
281       memset( mmm , 1, nvox ) ; mcount = nvox ;
282    } else {
283 
284       mmm = THD_makemask( mask_dset , miv , mask_bot , mask_top ) ;
285       if( mmm == NULL )
286          return " \n*** Can't make mask for some reason! ***\n" ;
287       mcount = THD_countmask( nvox , mmm ) ;
288 
289 #if 0
290       if( !EQUIV_DSETS(mask_dset,xdset) &&
291           !EQUIV_DSETS(mask_dset,ydset)   ) DSET_unload(mask_dset) ;
292 #endif
293 
294       if( mcount < 3 ){
295         free(mmm) ;
296         return " \n*** Less than 3 voxels survive the mask! ***\n" ;
297       }
298    }
299 
300    /*-- 30 Oct 2006: modify mask via Aboot Radius, if present --*/
301 
302    if( hrad > 0.0f ){
303       MCW_cluster *cl ;
304       short *di,*dj,*dk ;
305       int nd , xx,yy,zz , dd , nx,ny,nz,nxy, nx1,ny1,nz1 , ip,jp,kp ;
306 
307       cl = MCW_build_mask( fabs(DSET_DX(xdset)) ,
308                            fabs(DSET_DY(xdset)) ,
309                            fabs(DSET_DZ(xdset)) , hrad ) ;
310 
311       if( cl == NULL || cl->num_pt < 6 ){
312          KILL_CLUSTER(cl);
313          PLUTO_popup_transient(plint, " \n"
314                                       " Aboot Radius too small\n"
315                                       " for this dataset!\n"     ) ;
316       } else {
317          ADDTO_CLUSTER(cl,0,0,0,0) ;
318          di = cl->i ; dj = cl->j ; dk = cl->k ; nd = cl->num_pt ;
319          nx = DSET_NX(xdset) ; nx1 = nx-1 ;
320          ny = DSET_NY(xdset) ; ny1 = ny-1 ; nxy  = nx*ny ;
321          nz = DSET_NZ(xdset) ; nz1 = nz-1 ;
322          xx = plint->im3d->vinfo->i1 ;
323          yy = plint->im3d->vinfo->j2 ;
324          zz = plint->im3d->vinfo->k3 ;
325          for( dd=0 ; dd < nd ; dd++ ){
326            ip = xx+di[dd] ; if( ip < 0 || ip > nx1 ) continue ;
327            jp = yy+dj[dd] ; if( jp < 0 || jp > ny1 ) continue ;
328            kp = zz+dk[dd] ; if( kp < 0 || kp > nz1 ) continue ;
329            mmm[ip+jp*nx+kp*nxy]++ ;
330          }
331          KILL_CLUSTER(cl) ;
332          for( dd=0 ; dd < nvox ; dd++ ) if( mmm[dd] == 1 ) mmm[dd] = 0 ;
333          mcount = THD_countmask( nvox , mmm ) ;
334 
335          if( mcount < 3 ){
336            free(mmm) ;
337            return " \n*** Less than 3 voxels survive the mask+radius! ***\n" ;
338          }
339          if( perc == 0.0f || perc == 100.0f ){
340            char buf[256] ;
341            sprintf(buf," \n"
342                        " %d voxels in the mask+radius\n"
343                        " out of %d dataset voxels\n ",mcount,nvox) ;
344            PLUTO_popup_transient(plint,buf) ;
345          }
346       }
347    }
348 
349    /* 08 Aug 2007: modify mask by randomly selecting a percentage of it */
350 
351    if( perc > 0.0f && perc < 100.0f ){
352      double pp = 0.01 * perc ; char buf[256] ;
353      for( ii=0 ; ii < nvox ; ii++ ){
354        if( mmm[ii] && drand48() > pp ) mmm[ii] = 0 ;
355      }
356      mcount = THD_countmask( nvox , mmm ) ;
357      if( mcount < 3 ){
358        free(mmm) ;
359        return " \n*** Less than 3 voxels survive the percenting! ***\n" ;
360      }
361      sprintf(buf," \n"
362                  " %d voxels in mask after percenting\n"
363                  " out of %d dataset voxels (randomly)\n ",mcount,nvox) ;
364      PLUTO_popup_transient(plint,buf) ;
365    }
366 
367    /*-- allocate the space for the data to be plotted --*/
368 
369    xar = (float *) malloc(sizeof(float)*mcount) ;
370    yar = (float *) malloc(sizeof(float)*mcount) ;
371 
372    /*-- load values into x array --*/
373 
374    switch( DSET_BRICK_TYPE(xdset,ivx) ){
375       case MRI_short:{
376          short *bar = (short *) DSET_ARRAY(xdset,ivx) ;
377          float mfac = DSET_BRICK_FACTOR(xdset,ivx) ;
378          if( mfac == 0.0 ) mfac = 1.0 ;
379          for( ii=jj=0 ; ii < nvox ; ii++ )
380             if( mmm[ii] ) xar[jj++] = mfac*bar[ii] ;
381       }
382       break ;
383 
384       case MRI_byte:{
385          byte *bar = (byte *) DSET_ARRAY(xdset,ivx) ;
386          float mfac = DSET_BRICK_FACTOR(xdset,ivx) ;
387          if( mfac == 0.0 ) mfac = 1.0 ;
388          for( ii=jj=0 ; ii < nvox ; ii++ )
389             if( mmm[ii] ) xar[jj++] = mfac*bar[ii] ;
390       }
391       break ;
392 
393       case MRI_float:{
394          float *bar = (float *) DSET_ARRAY(xdset,ivx) ;
395          float mfac = DSET_BRICK_FACTOR(xdset,ivx) ;
396          if( mfac == 0.0 ) mfac = 1.0 ;
397          for( ii=jj=0 ; ii < nvox ; ii++ )
398             if( mmm[ii] ) xar[jj++] = mfac*bar[ii] ;
399       }
400       break ;
401 
402       default: break ;
403    }
404 
405    /*-- load values into y array --*/
406 
407    switch( DSET_BRICK_TYPE(ydset,ivy) ){
408       case MRI_short:{
409          short *bar = (short *) DSET_ARRAY(ydset,ivy) ;
410          float mfac = DSET_BRICK_FACTOR(ydset,ivy) ;
411          if( mfac == 0.0 ) mfac = 1.0 ;
412          for( ii=jj=0 ; ii < nvox ; ii++ )
413             if( mmm[ii] ) yar[jj++] = mfac*bar[ii] ;
414       }
415       break ;
416 
417       case MRI_byte:{
418          byte *bar = (byte *) DSET_ARRAY(ydset,ivy) ;
419          float mfac = DSET_BRICK_FACTOR(ydset,ivy) ;
420          if( mfac == 0.0 ) mfac = 1.0 ;
421          for( ii=jj=0 ; ii < nvox ; ii++ )
422             if( mmm[ii] ) yar[jj++] = mfac*bar[ii] ;
423       }
424       break ;
425 
426       case MRI_float:{
427          float *bar = (float *) DSET_ARRAY(ydset,ivy) ;
428          float mfac = DSET_BRICK_FACTOR(ydset,ivy) ;
429          if( mfac == 0.0 ) mfac = 1.0 ;
430          for( ii=jj=0 ; ii < nvox ; ii++ )
431             if( mmm[ii] ) yar[jj++] = mfac*bar[ii] ;
432       }
433       break ;
434 
435       default: break ;
436    }
437 
438    /* remove those voxels that aren't in the data ranges for both datasets */
439 
440    if( xbot < xtop || ybot < ytop ){
441       int nm ; float *tar ;
442 
443       /* make the mask of those that survive the x range */
444 
445       if( xbot < xtop ){
446          for( jj=0 ; jj < mcount ; jj++ )
447             mmm[jj] = ( xar[jj] >= xbot && xar[jj] <= xtop ) ;
448       } else {
449          memset( mmm , 1 , mcount ) ;
450       }
451 
452       /* and hit that mask with the y range */
453 
454       if( ybot < ytop ){
455          for( jj=0 ; jj < mcount ; jj++ )
456             if( mmm[jj] )
457                mmm[jj] = ( yar[jj] >= ybot && yar[jj] <= ytop ) ;
458       }
459 
460       nm = THD_countmask( mcount , mmm ) ;  /* how many are left */
461 
462       if( nm < 1 ){
463          free(mmm) ; free(xar) ; free(yar) ;
464          return " \n*** No values survive to be plotted! ***\n" ;
465       }
466 
467       /* copy survivors into new lists and free old lists */
468 
469       if( nm < mcount ){  /* we actually lost some points */
470          tar = (float *) malloc(sizeof(float)*nm) ;
471          for( ii=jj=0 ; ii < mcount ; ii++ )
472             if( mmm[ii] ) tar[jj++] = xar[ii] ;
473          free(xar) ; xar = tar ;
474 
475          tar = (float *) malloc(sizeof(float)*nm) ;
476          for( ii=jj=0 ; ii < mcount ; ii++ )
477             if( mmm[ii] ) tar[jj++] = yar[ii] ;
478          free(yar) ; yar = tar ;
479 
480          mcount = nm ;
481       }
482    }
483 
484    free(mmm) ;  /* don't need this no more */
485 
486 #define MAX_SPLOT 500000
487 
488    if( mcount > MAX_SPLOT ){
489       static char msg[128] ;
490       free(xar) ; free(yar) ;
491       sprintf(msg," \n*** Attempt to scatterplot %d points!\n"
492                      "*** Maximum allowed is     %d.\n"
493                      "*** Use the Percent function :-(\n" , mcount , MAX_SPLOT ) ;
494       return msg ;
495    }
496 
497    /* check for multiple (0,0) points [24 Dec 2014] */
498 
499    for( zzc=ii=0 ; ii < mcount ; ii++ ){
500      if( xar[ii] == 0.0f && yar[ii] == 0.0f ) zzc++ ;
501    }
502    if( zzc > 2+mcount/64 )
503      WARNING_message("Scatterplot has %d (0,0) pairs out of %d data points",zzc,mcount) ;
504 
505    /* compute the labels for the plot */
506 
507    if( xbot >= xtop ){
508       sprintf( xlab , "\\noesc %s[%d]" , DSET_FILECODE(xdset),ivx ) ;
509    } else {
510       AV_fval_to_char(xbot,ab) ; AV_fval_to_char(xtop,bb) ;
511       pab = ab ; if( *pab == ' ' ) pab++ ;
512       pbb = bb ; if( *pbb == ' ' ) pbb++ ;
513       sprintf( xlab , "\\noesc %s[%d]<%s..%s>" , DSET_FILECODE(xdset),ivx,pab,pbb ) ;
514    }
515 
516    if( ybot >= ytop ){
517       sprintf( ylab , "\\noesc %s[%d]" , DSET_FILECODE(ydset),ivy ) ;
518    } else {
519       AV_fval_to_char(ybot,ab) ; AV_fval_to_char(ytop,bb) ;
520       pab = ab ; if( *pab == ' ' ) pab++ ;
521       pbb = bb ; if( *pbb == ' ' ) pbb++ ;
522       sprintf( ylab , "\\noesc %s[%d]<%s..%s>" , DSET_FILECODE(ydset),ivy,pab,pbb ) ;
523    }
524 
525    /*- 11 Aug 2001: compute correlation coefficient -*/
526 
527    eplot = getenv("AFNI_SCATPLOT_LINES") ;
528    if( eplot == NULL ) eplot = "BOTH" ;
529 
530    if( mcount >= 5 ){           /* 02 Mar 2011: the new way [IPad-2 day!] */
531      float_triple aaa,bbb,rrr ; float_pair sab ;
532      if( mcount < 6666 ){
533        THD_pearson_corr_boot( mcount,xar,yar , &rrr,&aaa,&bbb ) ;
534        pcor = rrr.a ; p025 = rrr.b ; p975 = rrr.c ; a = aaa.a ; b = bbb.a ;
535        INFO_message("Scatterplot: R=%.3f  a=%.3f  b=%.3f  [y=ax+b fit]",pcor,a,b) ;
536        THD_spearman_corr_boot( mcount,xar,yar , &rrr ) ;
537        scor = rrr.a ; s025 = rrr.b ; s975 = rrr.c ;  /* 23 Dec 2014 */
538        did_corr = 2 ;
539      } else {
540        rrr = THD_pearson_indexed( mcount,NULL , xar,yar ) ;
541        pcor = rrr.c ; a = rrr.a ; b = rrr.b ; p025 = 6.66f ; p975 = -6.66f ;
542        mi = THD_mutual_info( mcount , xar , yar ) ;
543        INFO_message("Scatterplot: R=%.3f  a=%.3f  b=%.3f  [y=ax+b fit]  MI=%.3f",pcor,a,b,mi) ;
544        scor = THD_spearman_corr_nd(mcount,xar,yar) ;  /* 23 Dec 2014 */
545        did_corr = 1 ;
546      }
547      if( strcasestr(eplot,"NOL1") == NULL ){        /* turn on or off some line plots */
548        sab = THD_l1_fit_to_line( mcount,xar,yar ) ;
549        sa  = sab.a ; sb = sab.b ;
550      }
551      if( strcasestr(eplot,"NOL2") != NULL ){
552        a = b = 0.0f ;
553      }
554    }
555 
556    if( did_corr == 2 ) strcpy(tlab,"\\small ") ;
557    else                tlab[0] = '\0' ;
558 
559    if( label != NULL && *label != '\0' ){
560      strcpy(tlab+strlen(tlab),label) ;
561    } else {
562      sprintf(tlab+strlen(tlab),"[%d vox]",mcount) ;
563    }
564    if( did_corr ){
565      sprintf(tlab+strlen(tlab),"\\esc\\red  R\\approx %.2f",pcor) ;
566      if( p025 < pcor && p975 > pcor )
567        sprintf(tlab+strlen(tlab),"\\in[%.2f..%.2f]_{95%%}",p025,p975) ;
568      sprintf(tlab+strlen(tlab)," \\blue\\rho\\approx %.2f",scor) ;
569      if( s025 < scor && s975 > scor )
570        sprintf(tlab+strlen(tlab),"\\in[%.2f..%.2f]_{95%%}",s025,s975) ;
571      if( mi != 0.0f )
572        sprintf(tlab+strlen(tlab)," \\green MI\\approx %.2f",mi) ;
573      sprintf(tlab+strlen(tlab),"\\black") ;
574 
575      if( strlen(xlab) < 60 ){
576        if( a != 0.0f || b != 0.0f ){
577          sprintf(xlab+strlen(xlab),
578                  "\\esc\\red  \\{y\\approx %.3g*x%s%.3g\\}\\black" ,
579                  a , (b<0.0f)?"-":"+" , fabs(b) ) ;
580        } else if( sa != 0.0f || sb != 0.0f ){
581          sprintf(xlab+strlen(xlab),
582                  "\\esc\\blue  \\{y\\approx %.3g*x%s%.3g\\}\\black" ,
583                  sa , (sb<0.0f)?"-":"+" , fabs(b) ) ;
584        }
585      }
586    }
587 
588    /*-- actually plot data (cf. afni_plugin.c) --*/
589 
590    alin[0] = a ; blin[0] = b ; clin[0].a = 0.7f; clin[0].b = 0.0f; clin[0].c = 0.0f;
591    alin[1] = sa; blin[1] = sb; clin[1].a = 0.0f; clin[1].b = 0.0f; clin[1].c = 0.7f;
592 
593    PLUTO_scatterplot_NEW( mcount , xar,yar , xlab,ylab,tlab , 2,alin,blin,clin ) ;
594 
595    /*-- go home to papa --*/
596 
597    if( label != NULL ) free(label) ;
598    free(xar) ; free(yar) ; return NULL ;
599 }
600