00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008 #include "coxplot.h"
00009 #include "display.h"
00010
00011
00012
00013
00014
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
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
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
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
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
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
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
00142
00143 xobot = 0.15 ; xotop = 1.27 ;
00144 yobot = 0.1 ; yotop = 0.95 ;
00145
00146 if( STGOOD(tlab) ){ yotop -= 0.02 ; yobot -= 0.01 ; }
00147
00148
00149
00150 if( STGOOD(xlab) )
00151 plotpak_pwritf( 0.5*(xobot+xotop) , yobot-0.06 , xlab , 16 , 0 , 0 ) ;
00152
00153
00154
00155 if( STGOOD(ylab) )
00156 plotpak_pwritf( xobot-0.12 , 0.5*(yobot+yotop) , ylab , 16 , 90 , 0 ) ;
00157
00158
00159
00160 if( STGOOD(tlab) )
00161 plotpak_pwritf( xobot+0.01 , yotop+0.01 , tlab , 18 , 0 , -2 ) ;
00162
00163
00164
00165 plotpak_set( xobot,xotop , yobot,yotop , xbot,xtop , ybot,ytop , 1 ) ;
00166 plotpak_periml( nnax,mmax , nnay,mmay ) ;
00167
00168
00169
00170 #define NELL 64
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
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
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
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") ;
00294
00295
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 ){
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 ;
00397
00398 if( use > 1 && nx > use ) nx = use ;
00399
00400
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 }