Skip to content

AFNI/NIfTI Server

Sections
Personal tools
You are here: Home » AFNI » Documentation

Doxygen Source Code Documentation


Main Page   Alphabetical List   Data Structures   File List   Data Fields   Globals   Search  

1dsigplot.c

Go to the documentation of this file.
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 "mrilib.h"
00008 #include "coxplot.h"
00009 #include "display.h"
00010 
00011 /*----------------------------------------------------------------------
00012   Return p10 as a power of 10 such that
00013     p10 <= fabs(x) < 10*p10
00014   unless x == 0, in which case return 0.
00015 ------------------------------------------------------------------------*/
00016 
00017 static float p10( float x )
00018 {
00019    double y ;
00020    if( x == 0.0 ) return 0.0 ;
00021    if( x <  0.0 ) x = -x ;
00022    y = floor(log10(x)+0.000001) ; y = pow( 10.0 , y ) ;
00023    return (float) y ;
00024 }
00025 
00026 #define STGOOD(s) ( (s) != NULL && (s)[0] != '\0' )
00027 
00028 /*-----------------------------------------------------------------
00029    Plot a scatterplot with ellipses:
00030      npt  = # of points in x[] and y[]
00031      x    = x-axis values array
00032      y    = y-axis values array
00033      xsig = x standard deviation array
00034      ysig = y standard deviation array
00035      corr = correlation coefficient array
00036      pell = CDF for bivariate normal to be inside ellipse
00037      xlab } labels for x-axis,
00038      ylab }            y-axis
00039      tlab }        and top of graph (NULL => skip this label)
00040 -------------------------------------------------------------------*/
00041 
00042 MEM_plotdata * PLOT_scatterellipse( int npt ,
00043                                     float xlo,float xhi , float ylo,float yhi ,
00044                                     float * x , float * y ,
00045                                     float * xsig , float * ysig , float * corr ,
00046                                     float pell ,
00047                                     char * xlab , char * ylab , char * tlab )
00048 {
00049    int ii,jj , np , nnax,mmax , nnay,mmay ;
00050    float xbot,xtop , ybot,ytop , pbot,ptop ,
00051          xobot,xotop,yobot,yotop , xa,xb,ya,yb , dx,dy,dt ;
00052    char str[32] ;
00053    MEM_plotdata * mp ;
00054    float xu,yu , xd,yd ;
00055 
00056    if( npt < 2 || x == NULL || y == NULL ) return NULL ;
00057    if( xsig == NULL || ysig == NULL || corr == NULL ) return NULL ;
00058 
00059    if( pell <= 0.0 || pell >= 1.0 ) pell = 0.5 ;
00060 
00061    pell = -2.0 * log(1.0-pell) ;
00062 
00063    /* find range of data */
00064 
00065    xbot = xtop = x[0] ; ybot = ytop = y[0] ;
00066    for( ii=0 ; ii < npt ; ii++ ){
00067       xu = x[ii] + pell*xsig[ii] ;
00068       xd = x[ii] - pell*xsig[ii] ;
00069       if( xd < xbot ) xbot = xd ;
00070       if( xu > xtop ) xtop = xu ;
00071 
00072       yu = y[ii] + pell*ysig[ii] ;
00073       yd = y[ii] - pell*ysig[ii] ;
00074       if( yd < ybot ) ybot = yd ;
00075       if( yu > ytop ) ytop = yu ;
00076    }
00077    if( xbot >= xtop || ybot >= ytop ){
00078       fprintf(stderr,"*** Data has no range in PLOT_scatterellipse!\n\a");
00079       return NULL ;
00080    }
00081 
00082    if( xlo < xhi ){ xbot = xlo ; xtop = xhi ; }
00083    if( ylo < yhi ){ ybot = ylo ; ytop = yhi ; }
00084 
00085    /*-- push range of x outwards --*/
00086 
00087 #if 0
00088    pbot = p10(xbot) ; ptop = p10(xtop) ; if( ptop < pbot ) ptop = pbot ;
00089    if( ptop != 0.0 ){
00090       np = (xtop-xbot) / ptop + 0.5 ;
00091       switch( np ){
00092          case 1:  ptop *= 0.1  ; break ;
00093          case 2:  ptop *= 0.2  ; break ;
00094          case 3:  ptop *= 0.25 ; break ;
00095          case 4:
00096          case 5:  ptop *= 0.5  ; break ;
00097       }
00098       xbot = floor( xbot/ptop ) * ptop ;
00099       xtop =  ceil( xtop/ptop ) * ptop ;
00100       nnax = floor( (xtop-xbot) / ptop + 0.5 ) ;
00101       mmax = (nnax < 3) ? 10
00102                         : (nnax < 6) ? 5 : 2 ;
00103    } else {
00104       nnax = 1 ; mmax = 10 ;
00105    }
00106 #else
00107       nnax = 1 ; mmax = 10 ;
00108 #endif
00109 
00110    /*-- push range of y outwards --*/
00111 
00112 #if 0
00113    pbot = p10(ybot) ; ptop = p10(ytop) ; if( ptop < pbot ) ptop = pbot ;
00114    if( ptop != 0.0 ){
00115       np = (ytop-ybot) / ptop + 0.5 ;
00116       switch( np ){
00117          case 1:  ptop *= 0.1  ; break ;
00118          case 2:  ptop *= 0.2  ; break ;
00119          case 3:  ptop *= 0.25 ; break ;
00120          case 4:
00121          case 5:  ptop *= 0.5  ; break ;
00122       }
00123       ybot = floor( ybot/ptop ) * ptop ;
00124       ytop =  ceil( ytop/ptop ) * ptop ;
00125       nnay = floor( (ytop-ybot) / ptop + 0.5 ) ;
00126       mmay = (nnay < 3) ? 10
00127                         : (nnay < 6) ? 5 : 2 ;
00128    } else {
00129       nnay = 1 ; mmay = 10 ;
00130    }
00131 #else
00132       nnay = 1 ; mmay = 10 ;
00133 #endif
00134 
00135    /*-- setup to plot --*/
00136 
00137    create_memplot_surely( "SellPlot" , 1.3 ) ;
00138    set_color_memplot( 0.0 , 0.0 , 0.0 ) ;
00139    set_thick_memplot( 0.0 ) ;
00140 
00141    /*-- plot labels, if any --*/
00142 
00143    xobot = 0.15 ; xotop = 1.27 ;  /* set objective size of plot */
00144    yobot = 0.1  ; yotop = 0.95 ;
00145 
00146    if( STGOOD(tlab) ){ yotop -= 0.02 ; yobot -= 0.01 ; }
00147 
00148    /* x-axis label? */
00149 
00150    if( STGOOD(xlab) )
00151       plotpak_pwritf( 0.5*(xobot+xotop) , yobot-0.06 , xlab , 16 , 0 , 0 ) ;
00152 
00153    /* y-axis label? */
00154 
00155    if( STGOOD(ylab) )
00156       plotpak_pwritf( xobot-0.12 , 0.5*(yobot+yotop) , ylab , 16 , 90 , 0 ) ;
00157 
00158    /* label at top? */
00159 
00160    if( STGOOD(tlab) )
00161       plotpak_pwritf( xobot+0.01 , yotop+0.01 , tlab , 18 , 0 , -2 ) ;
00162 
00163    /* plot axes */
00164 
00165    plotpak_set( xobot,xotop , yobot,yotop , xbot,xtop , ybot,ytop , 1 ) ;
00166    plotpak_periml( nnax,mmax , nnay,mmay ) ;
00167 
00168    /* plot data */
00169 
00170 #define NELL 64                /* should be divisible by 4 */
00171 #define DTH  (2.0*PI/NELL)
00172 
00173    for( ii=0 ; ii < npt ; ii++ ){
00174       dx = pell * xsig[ii] ;
00175       dy = pell * ysig[ii] ;
00176       dt = asin(corr[ii]) ;
00177       xb = x[ii] + dx ;
00178       yb = y[ii] + dy*sin(dt) ;
00179       for( jj=1 ; jj <= NELL ; jj++ ){
00180          xa = xb ; ya = yb ;
00181          xb = x[ii] + dx*cos(jj*DTH) ;
00182          yb = y[ii] + dy*sin(jj*DTH+dt) ;
00183          plotpak_line( xa,ya , xb,yb ) ;
00184       }
00185    }
00186 
00187    set_color_memplot( 1.0 , 0.0 , 0.0 ) ;
00188    for( ii=0 ; ii < npt-1 ; ii++ )
00189      plotpak_line( x[ii],y[ii] , x[ii+1],y[ii+1] ) ;
00190 
00191 #define DSQ 0.005
00192 
00193    dx = DSQ*(xtop-xbot) ;
00194    dy = DSQ*(ytop-ybot) * (xotop-xobot)/(yotop-yobot) ;
00195    set_color_memplot( 0.0 , 0.0 , 1.0 ) ;
00196    for( ii=0 ; ii < npt ; ii++ ){
00197       xa = x[ii] - dx ; xb = x[ii] + dx ;
00198       ya = y[ii] - dy ; yb = y[ii] + dy ;
00199       plotpak_line( xa,ya , xa,yb ) ;
00200       plotpak_line( xa,yb , xb,yb ) ;
00201       plotpak_line( xb,yb , xb,ya ) ;
00202       plotpak_line( xb,ya , xa,ya ) ;
00203    }
00204 
00205    set_color_memplot( 0.0 , 0.0 , 0.0 ) ;
00206 
00207    mp = get_active_memplot() ;
00208    return mp ;
00209 }
00210 
00211 /*---- quickie program to look at some graphs - RWCox - Feb 1999 ----*/
00212 
00213 #define DEFAULT_NCOLOVR 20
00214 
00215 static char * INIT_colovr[DEFAULT_NCOLOVR] = {
00216    "#ffff00" , "#ffcc00"   , "#ff9900"  , "#ff6900" , "#ff4400" , "#ff0000" ,
00217    "#0000ff" , "#0044ff"   , "#0069ff"  , "#0099ff" , "#00ccff" , "#00ffff" ,
00218    "green"   , "limegreen" , "violet"   , "hotpink" ,
00219    "white"   , "#dddddd"   , "#bbbbbb"  , "black"
00220 } ;
00221 
00222 static char * INIT_labovr[DEFAULT_NCOLOVR] = {
00223    "yellow" , "yell-oran" , "oran-yell" , "orange"   , "oran-red" , "red"   ,
00224    "dk-blue", "blue"      , "lt-blue1"  , "lt-blue2" , "blue-cyan", "cyan"  ,
00225    "green"  , "limegreen" , "violet"    , "hotpink"  ,
00226    "white"  , "gry-dd"    , "gry-bb"    , "black"
00227 } ;
00228 
00229 static int nx ;
00230 static float *xx, *yy, *xsig, *ysig, *xycor , pell=0.5 ;
00231 
00232 static float xlo=1.0,xhi=-1.0 , ylo=1.0,yhi=-1.0 ;
00233 
00234 static MCW_DC * dc ;
00235 static char * title = NULL , * xlabel = NULL , * ylabel = NULL ;
00236 
00237 void startup_timeout_CB( XtPointer client_data , XtIntervalId * id ) ;
00238 
00239 int main( int argc , char * argv[] )
00240 {
00241    int iarg , ii , ny , ignore=0 , use=0 , install=0 ;
00242    char * tsfile , * cpt ;
00243    char dname[THD_MAX_NAME] , subv[THD_MAX_NAME] ;
00244    MRI_IMAGE * flim ;
00245    float * far ;
00246    XtAppContext app ;
00247    Widget shell ;
00248    int cxx=0 , cyy=1 , cxsig=2 , cysig=3 , crho=4 ;
00249 
00250    /*-- help? --*/
00251 
00252    if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00253      printf("Usage: 1dsigplot [options] infile\n"
00254             "Scatterplots a 5 column *.1D file to the screen.\n"
00255             "The columns of the input file must be\n"
00256             "   x  y  sigma_x sigma_y rho_xy\n"
00257             "\n"
00258             "Options:\n"
00259             " -install   = Install a new X11 colormap.\n"
00260             " -ignore nn = Skip first 'nn' rows in the input file\n"
00261             "                [default = 0]\n"
00262             " -use mm    = Plot 'mm' points [default = all of them]\n"
00263             " -xlabel aa = Put string 'aa' below the x-axis\n"
00264             "                [default = no axis label]\n"
00265             " -ylabel aa = Put string 'aa' to the left of the y-axis\n"
00266             "                [default = no axis label]\n"
00267             " -pell p    = Set CDF probability contour level for ellipses.\n"
00268             " -col abcde = 'abcde' is a permutation of '01234', and\n"
00269             "                specifies the column order for the data,\n"
00270             "                where a=column index for x\n"
00271             "                      b=column index for y\n"
00272             "                      c=column index for sigma_x\n"
00273             "                      d=column index for sigma_y\n"
00274             "                      e=column index for rho_xy\n"
00275             "                [default = 01234]\n"
00276             "\n"
00277             " -xrange x1 x2 = Range of x-axis\n"
00278             " -yrange y1 y2 = Range of y-axis\n"
00279            ) ;
00280       exit(0) ;
00281    }
00282 
00283    machdep() ;
00284 
00285    /* open X11 */
00286 
00287    shell = XtVaAppInitialize(
00288               &app , "AFNI" , NULL , 0 , &argc , argv , NULL , NULL ) ;
00289    if( shell == NULL ){
00290       fprintf(stderr,"** Cannot initialize X11!\n") ; exit(1) ;
00291    }
00292 
00293    cpt = my_getenv("TMPDIR") ;  /* just for fun */
00294 
00295    /*-- scan arguments that X11 didn't eat --*/
00296 
00297    iarg = 1 ;
00298    while( iarg < argc && argv[iarg][0] == '-' ){
00299 
00300      if( strcmp(argv[iarg],"-xrange") == 0 ){
00301         xlo = strtod(argv[++iarg],NULL) ;
00302         xhi = strtod(argv[++iarg],NULL) ;
00303         iarg++ ; continue ;
00304      }
00305 
00306      if( strcmp(argv[iarg],"-yrange") == 0 ){
00307         ylo = strtod(argv[++iarg],NULL) ;
00308         yhi = strtod(argv[++iarg],NULL) ;
00309         iarg++ ; continue ;
00310      }
00311 
00312      if( strcmp(argv[iarg],"-pell") == 0 ){
00313         pell = strtod(argv[++iarg],NULL) ;
00314         iarg++ ; continue ;
00315      }
00316 
00317      if( strcmp(argv[iarg],"-col") == 0 ){
00318         int ierr=0 ;
00319         iarg++ ;
00320         cxx   = argv[iarg][0] - '0' ; if( cxx   < 0 || cxx   > 4 ) ierr++ ;
00321         cyy   = argv[iarg][1] - '0' ; if( cyy   < 0 || cyy   > 4 ) ierr++ ;
00322         cxsig = argv[iarg][2] - '0' ; if( cxsig < 0 || cxsig > 4 ) ierr++ ;
00323         cysig = argv[iarg][3] - '0' ; if( cysig < 0 || cysig > 4 ) ierr++ ;
00324         crho  = argv[iarg][4] - '0' ; if( crho  < 0 || crho  > 4 ) ierr++ ;
00325         if( ierr || cxx+cyy+cxsig+cysig+crho != 10 ){
00326            fprintf(stderr,"*** Illegal argument after -ord!\n");exit(1);
00327         }
00328         iarg++ ; continue ;
00329      }
00330 
00331      if( strcmp(argv[iarg],"-install") == 0 ){
00332         install++ ; iarg++ ; continue ;
00333      }
00334 
00335      if( strcmp(argv[iarg],"-") == 0 ){  /* skip */
00336         iarg++ ; continue ;
00337      }
00338 
00339      if( strcmp(argv[iarg],"-title") == 0 ){
00340         title = argv[++iarg] ;
00341         iarg++ ; continue ;
00342      }
00343 
00344      if( strcmp(argv[iarg],"-xlabel") == 0 ){
00345         xlabel = argv[++iarg] ;
00346         iarg++ ; continue ;
00347      }
00348 
00349      if( strcmp(argv[iarg],"-ylabel") == 0 ){
00350         ylabel = argv[++iarg] ;
00351         iarg++ ; continue ;
00352      }
00353 
00354      if( strcmp(argv[iarg],"-ignore") == 0 ){
00355         ignore = strtod( argv[++iarg] , NULL ) ;
00356         if( ignore < 0 ){fprintf(stderr,"** Illegal -ignore value!\n");exit(1);}
00357         iarg++ ; continue ;
00358      }
00359 
00360      if( strcmp(argv[iarg],"-use") == 0 ){
00361         use = strtod( argv[++iarg] , NULL ) ;
00362         if( use < 2 ){fprintf(stderr,"** Illegal -use value!\n");exit(1);}
00363         iarg++ ; continue ;
00364      }
00365 
00366      fprintf(stderr,"** Unknown option: %s\n",argv[iarg]) ; exit(1) ;
00367    }
00368 
00369    if( iarg >= argc ){
00370       fprintf(stderr,"** No tsfile on command line!\n") ; exit(1) ;
00371    }
00372 
00373    dc = MCW_new_DC( shell , 16 ,
00374                     DEFAULT_NCOLOVR , INIT_colovr , INIT_labovr ,
00375                     1.0 , install ) ;
00376 
00377    flim = mri_read_1D( argv[iarg] ) ;
00378    if( flim == NULL ){
00379       fprintf(stderr,"** Can't read input file %s\n",argv[iarg]) ;
00380       exit(1);
00381    }
00382 
00383    if( flim->ny != 5 ){
00384       fprintf(stderr,"** Input file doesn't have exactly 5 columns!\n") ;
00385       exit(1) ;
00386    }
00387    far  = MRI_FLOAT_PTR(flim) ;
00388    nx   = flim->nx ;
00389 
00390    xx    = far + cxx  *nx ;
00391    yy    = far + cyy  *nx ;
00392    xsig  = far + cxsig*nx ;
00393    ysig  = far + cysig*nx ;
00394    xycor = far + crho *nx ;
00395 
00396    nx = nx - ignore ;  /* cut off the ignored points */
00397 
00398    if( use > 1 && nx > use ) nx = use ;
00399 
00400    /* start X11 */
00401 
00402    (void) XtAppAddTimeOut( app , 123 , startup_timeout_CB , NULL ) ;
00403 
00404    XtAppMainLoop(app) ;
00405    exit(0) ;
00406 }
00407 
00408 /*-----------------------------------------------------------------*/
00409 void killfunc(void * fred){ exit(0) ; }
00410 /*-----------------------------------------------------------------*/
00411 
00412 void startup_timeout_CB( XtPointer client_data , XtIntervalId * id )
00413 {
00414    MEM_plotdata * mp ;
00415 
00416    mp = PLOT_scatterellipse( nx , xlo,xhi,ylo,yhi ,
00417                              xx,yy,xsig,ysig,xycor ,
00418                              pell , xlabel,ylabel,title ) ;
00419 
00420    if( mp != NULL )
00421       (void) memplot_to_topshell( dc->display , mp , killfunc ) ;
00422 
00423    return ;
00424 }
 

Powered by Plone

This site conforms to the following standards: