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 
00009 #ifndef ALLOW_PLUGINS
00010 #  error "Plugins not properly set up -- see machdep.h"
00011 #endif
00012 
00013 /***********************************************************************
00014   Plugin to plot a scatterplot of 2 bricks.
00015   [Adapted from plug_histog.c - RWCox - 13 Jan 2000]
00016   11 Aug 2001: modified to print correlation coefficient
00017 ************************************************************************/
00018 
00019 char * SCAT_main( PLUGIN_interface * ) ;
00020 
00021 static char helpstring[] =
00022    " \n"
00023    " Purpose:  Scatterplot data from 2 bricks.\n"
00024    "\n"
00025    " Source x: Dataset   = dataset for x-axis values\n"
00026    "           Sub-brick = which one to use\n"
00027    "           Bottom    = minimum value to include\n"
00028    "           Top       = maximum value to include\n"
00029    "\n"
00030    " Source y: Dataset   = dataset for y-axis values\n"
00031    "           Sub-brick = which one to use\n"
00032    "           Bottom    = minimum value to include\n"
00033    "           Top       = maximum value to include\n"
00034    "\n"
00035    " Mask:     Dataset   = masking dataset\n"
00036    "           Sub-brick = which one to use\n"
00037    "           Bottom    = min value from mask dataset to use\n"
00038    "           Top       = max value from mask dataset to use\n"
00039    "\n"
00040    "In the above definitions:\n"
00041    "  if Bottom >  Top, then all voxels will be used\n"
00042    "  if Bottom <= Top, then only voxels in this range will be used\n"
00043    "\n"
00044    "Note: A maximum of 4,000,000 (four million) points can be plotted.\n\n"
00045    " Author -- RW Cox - 13 January 2000\n"
00046 ;
00047 
00048 /***********************************************************************
00049    Set up the interface to the user
00050 ************************************************************************/
00051 
00052 
00053 DEFINE_PLUGIN_PROTOTYPE
00054 
00055 PLUGIN_interface * PLUGIN_init( int ncall )
00056 {
00057    PLUGIN_interface * plint ;
00058 
00059    if( ncall > 0 ) return NULL ;  /* only one interface */
00060 
00061    /*-- set titles and call point --*/
00062 
00063    plint = PLUTO_new_interface( "ScatterPlot" ,
00064                                 "ScatterPlot of 2 Bricks" ,
00065                                 helpstring ,
00066                                 PLUGIN_CALL_VIA_MENU , SCAT_main  ) ;
00067 
00068    PLUTO_add_hint( plint , "Of 2 Bricks" ) ;
00069 
00070    PLUTO_set_sequence( plint , "A:afniinfo:dset" ) ;
00071 
00072    PLUTO_set_runlabels( plint , "Plot+Keep" , "Plot+Close" ) ;  /* 04 Nov 2003 */
00073 
00074    /*-- first line of input --*/
00075 
00076    PLUTO_add_option( plint , "Source x" , "Source x" , TRUE ) ;
00077    PLUTO_add_dataset(plint , "Dataset" ,
00078                                     ANAT_ALL_MASK , FUNC_ALL_MASK ,
00079                                     DIMEN_ALL_MASK | BRICK_ALLREAL_MASK ) ;
00080    PLUTO_add_number( plint , "Sub-brick" , 0,9999,0 , 0,1 ) ;
00081    PLUTO_add_number( plint , "Bottom" , -99999,99999, 0,  1,1 ) ;
00082    PLUTO_add_number( plint , "Top"    , -99999,99999, 0, -1,1 ) ;
00083 
00084    /*-- second line of input --*/
00085 
00086    PLUTO_add_option( plint , "Source y" , "Source y" , TRUE ) ;
00087    PLUTO_add_dataset(plint , "Dataset" ,
00088                                     ANAT_ALL_MASK , FUNC_ALL_MASK ,
00089                                     DIMEN_ALL_MASK | BRICK_ALLREAL_MASK ) ;
00090    PLUTO_add_number( plint , "Sub-brick" , 0,9999,0 , 0,1 ) ;
00091    PLUTO_add_number( plint , "Bottom" , -99999,99999, 0,  1,1 ) ;
00092    PLUTO_add_number( plint , "Top"    , -99999,99999, 0, -1,1 ) ;
00093 
00094    /*-- third line of input --*/
00095 
00096    PLUTO_add_option( plint , "Mask" , "Mask" , FALSE ) ;
00097    PLUTO_add_dataset( plint , "Dataset" ,
00098                                     ANAT_ALL_MASK , FUNC_ALL_MASK ,
00099                                     DIMEN_ALL_MASK | BRICK_ALLREAL_MASK ) ;
00100    PLUTO_add_number( plint , "Sub-brick" , 0,9999,0 , 0,1 ) ;
00101    PLUTO_add_number( plint , "Bottom" , -99999,99999, 0,  1,1 ) ;
00102    PLUTO_add_number( plint , "Top"    , -99999,99999, 0, -1,1 ) ;
00103 
00104    return plint ;
00105 }
00106 
00107 /***************************************************************************
00108   Main routine for this plugin (will be called from AFNI).
00109 ****************************************************************************/
00110 
00111 char * SCAT_main( PLUGIN_interface *plint )
00112 {
00113    MCW_idcode *idc ;
00114    THD_3dim_dataset *xdset, *ydset , *mask_dset=NULL ;
00115    int ivx,ivy , mcount , nvox , ii,jj , nbin=-1 ;
00116    float mask_bot=666.0 , mask_top=-666.0 ;
00117    float xbot,xtop , ybot,ytop , pcor=0 , a=0,b=0 ;
00118    char *tag , *str ;
00119    char xlab[THD_MAX_NAME],ylab[THD_MAX_NAME],tlab[THD_MAX_NAME] ;
00120    char ab[16] , bb[16] , *pab,*pbb ;
00121    byte *mmm ;
00122    float *xar , *yar ;
00123 
00124    char *cname=NULL ;  /* 06 Aug 1998 */
00125    int miv=0 ;
00126 
00127    /*--------------------------------------------------------------------*/
00128    /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
00129 
00130    if( plint == NULL )
00131       return "**********************\n"
00132              "SCAT_main:  NULL input\n"
00133              "**********************"  ;
00134 
00135    /*-- read x dataset line --*/
00136 
00137    PLUTO_next_option(plint) ;
00138    idc   = PLUTO_get_idcode(plint) ;
00139    xdset = PLUTO_find_dset(idc) ;
00140    if( xdset == NULL )
00141       return "*************************\n"
00142              "SCAT_main:  bad x dataset\n"
00143              "*************************"  ;
00144 
00145    ivx = (int) PLUTO_get_number(plint) ;
00146    if( ivx >= DSET_NVALS(xdset) || ivx < 0 )
00147       return "***************************\n"
00148              "SCAT_main:  bad x sub-brick\n"
00149              "***************************" ;
00150 
00151    DSET_load(xdset) ;
00152    if( DSET_ARRAY(xdset,ivx) == NULL )
00153       return "********************************\n"
00154              "SCAT_main:  can't load x dataset\n"
00155              "********************************"  ;
00156    nvox = DSET_NVOX(xdset) ;
00157 
00158    xbot = PLUTO_get_number(plint) ;
00159    xtop = PLUTO_get_number(plint) ;
00160 
00161    /*-- read y dataset line --*/
00162 
00163    PLUTO_next_option(plint) ;
00164    idc   = PLUTO_get_idcode(plint) ;
00165    ydset = PLUTO_find_dset(idc) ;
00166    if( ydset == NULL )
00167       return "*************************\n"
00168              "SCAT_main:  bad y dataset\n"
00169              "*************************"  ;
00170 
00171    ivy = (int) PLUTO_get_number(plint) ;
00172    if( ivy >= DSET_NVALS(ydset) || ivy < 0 )
00173       return "***************************\n"
00174              "SCAT_main:  bad y sub-brick\n"
00175              "***************************" ;
00176 
00177    if( DSET_NVOX(ydset) != nvox )
00178       return "************************************************\n"
00179              "SCAT_main:  x and y datasets don't match in size\n"
00180              "************************************************" ;
00181 
00182    DSET_load(ydset) ;
00183    if( DSET_ARRAY(ydset,ivy) == NULL )
00184       return "********************************\n"
00185              "SCAT_main:  can't load y dataset\n"
00186              "********************************"  ;
00187 
00188    ybot = PLUTO_get_number(plint) ;
00189    ytop = PLUTO_get_number(plint) ;
00190 
00191    /*-- read optional line(s) --*/
00192 
00193    while( (tag=PLUTO_get_optiontag(plint)) != NULL ){
00194 
00195       /*-- Mask dataset --*/
00196 
00197       if( strcmp(tag,"Mask") == 0 ){
00198 
00199          idc       = PLUTO_get_idcode(plint) ;
00200          mask_dset = PLUTO_find_dset(idc) ;
00201 
00202          if( mask_dset == NULL )
00203             return "****************************\n"
00204                    "SCAT_main:  bad mask dataset\n"
00205                    "****************************"  ;
00206 
00207          if( DSET_NVOX(mask_dset) != nvox )
00208            return "**********************************************************\n"
00209                   "SCAT_main: mask input dataset doesn't match source dataset\n"
00210                   "**********************************************************" ;
00211 
00212          miv = (int) PLUTO_get_number(plint) ;  /* 06 Aug 1998 */
00213          if( miv >= DSET_NVALS(mask_dset) || miv < 0 )
00214             return "**************************************************\n"
00215                    "SCAT_main: mask dataset sub-brick index is illegal\n"
00216                    "**************************************************"  ;
00217 
00218          DSET_load(mask_dset) ;
00219          if( DSET_ARRAY(mask_dset,miv) == NULL )
00220             return "***********************************\n"
00221                    "SCAT_main:  can't load mask dataset\n"
00222                    "***********************************"  ;
00223 
00224          mask_bot = PLUTO_get_number(plint) ;
00225          mask_top = PLUTO_get_number(plint) ;
00226          continue ;
00227       }
00228    }
00229 
00230    /*------------------------------------------------------*/
00231    /*---------- At this point, the inputs are OK ----------*/
00232 
00233    /*-- build the byte mask array --*/
00234 
00235    if( mask_dset == NULL ){
00236       mmm = (byte *) malloc( sizeof(byte) * nvox ) ;
00237       if( mmm == NULL )
00238          return " \n*** Can't malloc workspace for mask! ***\n" ;
00239       memset( mmm , 1, nvox ) ; mcount = nvox ;
00240    } else {
00241 
00242       mmm = THD_makemask( mask_dset , miv , mask_bot , mask_top ) ;
00243       if( mmm == NULL )
00244          return " \n*** Can't make mask for some reason! ***\n" ;
00245       mcount = THD_countmask( nvox , mmm ) ;
00246 
00247       if( !EQUIV_DSETS(mask_dset,xdset) &&
00248           !EQUIV_DSETS(mask_dset,ydset)   ) DSET_unload(mask_dset) ;
00249 
00250       if( mcount < 3 ){
00251          free(mmm) ;
00252          return " \n*** Less than 3 voxels survive the mask! ***\n" ;
00253       }
00254    }
00255 
00256    /*-- allocate the space for the data to be plotted --*/
00257 
00258    xar = (float *) malloc(sizeof(float)*mcount) ;
00259    yar = (float *) malloc(sizeof(float)*mcount) ;
00260 
00261    /*-- load values into x array --*/
00262 
00263    switch( DSET_BRICK_TYPE(xdset,ivx) ){
00264       case MRI_short:{
00265          short *bar = (short *) DSET_ARRAY(xdset,ivx) ;
00266          float mfac = DSET_BRICK_FACTOR(xdset,ivx) ;
00267          if( mfac == 0.0 ) mfac = 1.0 ;
00268          for( ii=jj=0 ; ii < nvox ; ii++ )
00269             if( mmm[ii] ) xar[jj++] = mfac*bar[ii] ;
00270       }
00271       break ;
00272 
00273       case MRI_byte:{
00274          byte *bar = (byte *) DSET_ARRAY(xdset,ivx) ;
00275          float mfac = DSET_BRICK_FACTOR(xdset,ivx) ;
00276          if( mfac == 0.0 ) mfac = 1.0 ;
00277          for( ii=jj=0 ; ii < nvox ; ii++ )
00278             if( mmm[ii] ) xar[jj++] = mfac*bar[ii] ;
00279       }
00280       break ;
00281 
00282       case MRI_float:{
00283          float *bar = (float *) DSET_ARRAY(xdset,ivx) ;
00284          float mfac = DSET_BRICK_FACTOR(xdset,ivx) ;
00285          if( mfac == 0.0 ) mfac = 1.0 ;
00286          for( ii=jj=0 ; ii < nvox ; ii++ )
00287             if( mmm[ii] ) xar[jj++] = mfac*bar[ii] ;
00288       }
00289       break ;
00290    }
00291 
00292    /*-- load values into y array --*/
00293 
00294    switch( DSET_BRICK_TYPE(ydset,ivy) ){
00295       case MRI_short:{
00296          short *bar = (short *) DSET_ARRAY(ydset,ivy) ;
00297          float mfac = DSET_BRICK_FACTOR(ydset,ivy) ;
00298          if( mfac == 0.0 ) mfac = 1.0 ;
00299          for( ii=jj=0 ; ii < nvox ; ii++ )
00300             if( mmm[ii] ) yar[jj++] = mfac*bar[ii] ;
00301       }
00302       break ;
00303 
00304       case MRI_byte:{
00305          byte *bar = (byte *) DSET_ARRAY(ydset,ivy) ;
00306          float mfac = DSET_BRICK_FACTOR(ydset,ivy) ;
00307          if( mfac == 0.0 ) mfac = 1.0 ;
00308          for( ii=jj=0 ; ii < nvox ; ii++ )
00309             if( mmm[ii] ) yar[jj++] = mfac*bar[ii] ;
00310       }
00311       break ;
00312 
00313       case MRI_float:{
00314          float *bar = (float *) DSET_ARRAY(ydset,ivy) ;
00315          float mfac = DSET_BRICK_FACTOR(ydset,ivy) ;
00316          if( mfac == 0.0 ) mfac = 1.0 ;
00317          for( ii=jj=0 ; ii < nvox ; ii++ )
00318             if( mmm[ii] ) yar[jj++] = mfac*bar[ii] ;
00319       }
00320       break ;
00321    }
00322 
00323    /* remove those voxels that aren't in the data ranges for both datasets */
00324 
00325    if( xbot < xtop || ybot < ytop ){
00326       int nm ; float *tar ;
00327 
00328       /* make the mask of those that survive the x range */
00329 
00330       if( xbot < xtop ){
00331          for( jj=0 ; jj < mcount ; jj++ )
00332             mmm[jj] = ( xar[jj] >= xbot && xar[jj] <= xtop ) ;
00333       } else {
00334          memset( mmm , 1 , mcount ) ;
00335       }
00336 
00337       /* and hit that mask with the y range */
00338 
00339       if( ybot < ytop ){
00340          for( jj=0 ; jj < mcount ; jj++ )
00341             if( mmm[jj] )
00342                mmm[jj] = ( yar[jj] >= ybot && yar[jj] <= ytop ) ;
00343       }
00344 
00345       nm = THD_countmask( mcount , mmm ) ;  /* how many are left */
00346 
00347       if( nm < 1 ){
00348          free(mmm) ; free(xar) ; free(yar) ;
00349          return " \n*** No values survive to be plotted! ***\n" ;
00350       }
00351 
00352       /* copy survivors into new lists and free old lists */
00353 
00354       if( nm < mcount ){  /* we actually lost some points */
00355          tar = (float *) malloc(sizeof(float)*nm) ;
00356          for( ii=jj=0 ; ii < mcount ; ii++ )
00357             if( mmm[ii] ) tar[jj++] = xar[ii] ;
00358          free(xar) ; xar = tar ;
00359 
00360          tar = (float *) malloc(sizeof(float)*nm) ;
00361          for( ii=jj=0 ; ii < mcount ; ii++ )
00362             if( mmm[ii] ) tar[jj++] = yar[ii] ;
00363          free(yar) ; yar = tar ;
00364 
00365          mcount = nm ;
00366       }
00367    }
00368 
00369    free(mmm) ;  /* don't need this no more */
00370 
00371    if( mcount > 4000000 ){
00372       static char msg[128] ;
00373       free(xar) ; free(yar) ;
00374       sprintf(msg," \n*** Attempt to scatterplot %d points!\n"
00375                      "*** Maximum allowed is     4000000.\n" , mcount ) ;
00376       return msg ;
00377    }
00378 
00379    /* compute the labels for the plot */
00380 
00381    if( xbot >= xtop ){
00382       sprintf( xlab , "\\noesc %s[%d]" , DSET_FILECODE(xdset),ivx ) ;
00383    } else {
00384       AV_fval_to_char(xbot,ab) ; AV_fval_to_char(xtop,bb) ;
00385       pab = ab ; if( *pab == ' ' ) pab++ ;
00386       pbb = bb ; if( *pbb == ' ' ) pbb++ ;
00387       sprintf( xlab , "\\noesc %s[%d]<%s..%s>" , DSET_FILECODE(xdset),ivx,pab,pbb ) ;
00388    }
00389 
00390    if( ybot >= ytop ){
00391       sprintf( ylab , "%s[%d]" , DSET_FILECODE(ydset),ivy ) ;
00392    } else {
00393       AV_fval_to_char(ybot,ab) ; AV_fval_to_char(ytop,bb) ;
00394       pab = ab ; if( *pab == ' ' ) pab++ ;
00395       pbb = bb ; if( *pbb == ' ' ) pbb++ ;
00396       sprintf( ylab , "\\noesc %s[%d]<%s..%s>" , DSET_FILECODE(ydset),ivy,pab,pbb ) ;
00397    }
00398 
00399    /*- 11 Aug 2001: compute correlation coefficient -*/
00400 
00401    if( mcount > 1 ){
00402      float xbar=0,ybar=0 , xq=0,yq=0,xyq=0 ;
00403      for( ii=0 ; ii < mcount ; ii++ ){ xbar += xar[ii]; ybar += yar[ii]; }
00404      xbar /= mcount ; ybar /= mcount ;
00405      for( ii=0 ; ii < mcount ; ii++ ){
00406         xq  += (xar[ii]-xbar)*(xar[ii]-xbar) ;
00407         yq  += (yar[ii]-ybar)*(yar[ii]-ybar) ;
00408         xyq += (xar[ii]-xbar)*(yar[ii]-ybar) ;
00409      }
00410      if( xq > 0.0 && yq > 0.0 ){
00411        pcor = xyq/sqrt(xq*yq); a = xyq/xq; b = (xq*ybar-xbar*xyq)/xq; }
00412    }
00413 
00414    if( mask_dset == NULL ){
00415       sprintf(tlab,"Scatter Plot: %d Voxels",mcount) ;
00416    } else {
00417       sprintf(tlab,"\\noesc Scatter Plot: %d Voxels (%s)",
00418               mcount , DSET_FILECODE(mask_dset)   ) ;
00419    }
00420    if( pcor != 0.0 ){
00421      char abuf[16] ; AV_fval_to_char(a,abuf) ;
00422      sprintf(tlab+strlen(tlab)," R=%.3f a=%s",pcor,abuf) ;
00423    }
00424 
00425    /*-- actually plot data --*/
00426 
00427    PLUTO_scatterplot( mcount , xar,yar , xlab,ylab,tlab , a,b ) ;
00428 
00429    /*-- go home to papa --*/
00430 
00431    free(xar) ; free(yar) ; return NULL ;
00432 }