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 "mrilib.h"
8 #include "coxplot.h"
9 #include "display.h"
10
11 /*----------------------------------------------------------------------
12 Return p10 as a power of 10 such that
13 p10 <= fabs(x) < 10*p10
14 unless x == 0, in which case return 0.
15 ------------------------------------------------------------------------*/
16
p10(float x)17 static float p10( float x )
18 {
19 double y ;
20 if( x == 0.0 ) return 0.0 ;
21 if( x < 0.0 ) x = -x ;
22 y = floor(log10(x)+0.000001) ; y = pow( 10.0 , y ) ;
23 return (float) y ;
24 }
25
26 #define STGOOD(s) ( (s) != NULL && (s)[0] != '\0' )
27
28 /*-----------------------------------------------------------------
29 Plot a scatterplot with ellipses:
30 npt = # of points in x[] and y[]
31 x = x-axis values array
32 y = y-axis values array
33 xsig = x standard deviation array
34 ysig = y standard deviation array
35 corr = correlation coefficient array
36 pell = CDF for bivariate normal to be inside ellipse
37 xlab } labels for x-axis,
38 ylab } y-axis
39 tlab } and top of graph (NULL => skip this label)
40 -------------------------------------------------------------------*/
41
PLOT_scatterellipse(int npt,float xlo,float xhi,float ylo,float yhi,float * x,float * y,float * xsig,float * ysig,float * corr,float pell,char * xlab,char * ylab,char * tlab)42 MEM_plotdata * PLOT_scatterellipse( int npt ,
43 float xlo,float xhi , float ylo,float yhi ,
44 float * x , float * y ,
45 float * xsig , float * ysig , float * corr ,
46 float pell ,
47 char * xlab , char * ylab , char * tlab )
48 {
49 int ii,jj , np , nnax,mmax , nnay,mmay ;
50 float xbot,xtop , ybot,ytop , pbot,ptop ,
51 xobot,xotop,yobot,yotop , xa,xb,ya,yb , dx,dy,dt ;
52 char str[32] ;
53 MEM_plotdata * mp ;
54 float xu,yu , xd,yd ;
55
56 if( npt < 2 || x == NULL || y == NULL ) return NULL ;
57 if( xsig == NULL || ysig == NULL || corr == NULL ) return NULL ;
58
59 if( pell <= 0.0 || pell >= 1.0 ) pell = 0.5 ;
60
61 pell = -2.0 * log(1.0-pell) ;
62
63 /* find range of data */
64
65 xbot = xtop = x[0] ; ybot = ytop = y[0] ;
66 for( ii=0 ; ii < npt ; ii++ ){
67 xu = x[ii] + pell*xsig[ii] ;
68 xd = x[ii] - pell*xsig[ii] ;
69 if( xd < xbot ) xbot = xd ;
70 if( xu > xtop ) xtop = xu ;
71
72 yu = y[ii] + pell*ysig[ii] ;
73 yd = y[ii] - pell*ysig[ii] ;
74 if( yd < ybot ) ybot = yd ;
75 if( yu > ytop ) ytop = yu ;
76 }
77 if( xbot >= xtop || ybot >= ytop ){
78 fprintf(stderr,"*** Data has no range in PLOT_scatterellipse!\n\a");
79 return NULL ;
80 }
81
82 if( xlo < xhi ){ xbot = xlo ; xtop = xhi ; }
83 if( ylo < yhi ){ ybot = ylo ; ytop = yhi ; }
84
85 /*-- push range of x outwards --*/
86
87 #if 0
88 pbot = p10(xbot) ; ptop = p10(xtop) ; if( ptop < pbot ) ptop = pbot ;
89 if( ptop != 0.0 ){
90 np = (xtop-xbot) / ptop + 0.5 ;
91 switch( np ){
92 case 1: ptop *= 0.1 ; break ;
93 case 2: ptop *= 0.2 ; break ;
94 case 3: ptop *= 0.25 ; break ;
95 case 4:
96 case 5: ptop *= 0.5 ; break ;
97 }
98 xbot = floor( xbot/ptop ) * ptop ;
99 xtop = ceil( xtop/ptop ) * ptop ;
100 nnax = floor( (xtop-xbot) / ptop + 0.5 ) ;
101 mmax = (nnax < 3) ? 10
102 : (nnax < 6) ? 5 : 2 ;
103 } else {
104 nnax = 1 ; mmax = 10 ;
105 }
106 #else
107 nnax = 1 ; mmax = 10 ;
108 #endif
109
110 /*-- push range of y outwards --*/
111
112 #if 0
113 pbot = p10(ybot) ; ptop = p10(ytop) ; if( ptop < pbot ) ptop = pbot ;
114 if( ptop != 0.0 ){
115 np = (ytop-ybot) / ptop + 0.5 ;
116 switch( np ){
117 case 1: ptop *= 0.1 ; break ;
118 case 2: ptop *= 0.2 ; break ;
119 case 3: ptop *= 0.25 ; break ;
120 case 4:
121 case 5: ptop *= 0.5 ; break ;
122 }
123 ybot = floor( ybot/ptop ) * ptop ;
124 ytop = ceil( ytop/ptop ) * ptop ;
125 nnay = floor( (ytop-ybot) / ptop + 0.5 ) ;
126 mmay = (nnay < 3) ? 10
127 : (nnay < 6) ? 5 : 2 ;
128 } else {
129 nnay = 1 ; mmay = 10 ;
130 }
131 #else
132 nnay = 1 ; mmay = 10 ;
133 #endif
134
135 /*-- setup to plot --*/
136
137 create_memplot_surely( "SellPlot" , 1.3 ) ;
138 set_color_memplot( 0.0 , 0.0 , 0.0 ) ;
139 set_thick_memplot( 0.0 ) ;
140
141 /*-- plot labels, if any --*/
142
143 xobot = 0.15 ; xotop = 1.27 ; /* set objective size of plot */
144 yobot = 0.1 ; yotop = 0.95 ;
145
146 if( STGOOD(tlab) ){ yotop -= 0.02 ; yobot -= 0.01 ; }
147
148 /* x-axis label? */
149
150 if( STGOOD(xlab) )
151 plotpak_pwritf( 0.5*(xobot+xotop) , yobot-0.06 , xlab , 16 , 0 , 0 ) ;
152
153 /* y-axis label? */
154
155 if( STGOOD(ylab) )
156 plotpak_pwritf( xobot-0.12 , 0.5*(yobot+yotop) , ylab , 16 , 90 , 0 ) ;
157
158 /* label at top? */
159
160 if( STGOOD(tlab) )
161 plotpak_pwritf( xobot+0.01 , yotop+0.01 , tlab , 18 , 0 , -2 ) ;
162
163 /* plot axes */
164
165 plotpak_set( xobot,xotop , yobot,yotop , xbot,xtop , ybot,ytop , 1 ) ;
166 plotpak_periml( nnax,mmax , nnay,mmay ) ;
167
168 /* plot data */
169
170 #define NELL 64 /* should be divisible by 4 */
171 #define DTH (2.0*PI/NELL)
172
173 for( ii=0 ; ii < npt ; ii++ ){
174 dx = pell * xsig[ii] ;
175 dy = pell * ysig[ii] ;
176 dt = asin(corr[ii]) ;
177 xb = x[ii] + dx ;
178 yb = y[ii] + dy*sin(dt) ;
179 for( jj=1 ; jj <= NELL ; jj++ ){
180 xa = xb ; ya = yb ;
181 xb = x[ii] + dx*cos(jj*DTH) ;
182 yb = y[ii] + dy*sin(jj*DTH+dt) ;
183 plotpak_line( xa,ya , xb,yb ) ;
184 }
185 }
186
187 set_color_memplot( 1.0 , 0.0 , 0.0 ) ;
188 for( ii=0 ; ii < npt-1 ; ii++ )
189 plotpak_line( x[ii],y[ii] , x[ii+1],y[ii+1] ) ;
190
191 #define DSQ 0.005
192
193 dx = DSQ*(xtop-xbot) ;
194 dy = DSQ*(ytop-ybot) * (xotop-xobot)/(yotop-yobot) ;
195 set_color_memplot( 0.0 , 0.0 , 1.0 ) ;
196 for( ii=0 ; ii < npt ; ii++ ){
197 xa = x[ii] - dx ; xb = x[ii] + dx ;
198 ya = y[ii] - dy ; yb = y[ii] + dy ;
199 plotpak_line( xa,ya , xa,yb ) ;
200 plotpak_line( xa,yb , xb,yb ) ;
201 plotpak_line( xb,yb , xb,ya ) ;
202 plotpak_line( xb,ya , xa,ya ) ;
203 }
204
205 set_color_memplot( 0.0 , 0.0 , 0.0 ) ;
206
207 mp = get_active_memplot() ;
208 return mp ;
209 }
210
211 /*---- quickie program to look at some graphs - RWCox - Feb 1999 ----*/
212
213 #define DEFAULT_NCOLOVR 20
214
215 static char * INIT_colovr[DEFAULT_NCOLOVR] = {
216 "#ffff00" , "#ffcc00" , "#ff9900" , "#ff6900" , "#ff4400" , "#ff0000" ,
217 "#0000ff" , "#0044ff" , "#0069ff" , "#0099ff" , "#00ccff" , "#00ffff" ,
218 "green" , "limegreen" , "violet" , "hotpink" ,
219 "white" , "#dddddd" , "#bbbbbb" , "black"
220 } ;
221
222 static char * INIT_labovr[DEFAULT_NCOLOVR] = {
223 "yellow" , "yell-oran" , "oran-yell" , "orange" , "oran-red" , "red" ,
224 "dk-blue", "blue" , "lt-blue1" , "lt-blue2" , "blue-cyan", "cyan" ,
225 "green" , "limegreen" , "violet" , "hotpink" ,
226 "white" , "gry-dd" , "gry-bb" , "black"
227 } ;
228
229 static int nx ;
230 static float *xx, *yy, *xsig, *ysig, *xycor , pell=0.5 ;
231
232 static float xlo=1.0,xhi=-1.0 , ylo=1.0,yhi=-1.0 ;
233
234 static MCW_DC * dc ;
235 static char * title = NULL , * xlabel = NULL , * ylabel = NULL ;
236
237 void startup_timeout_CB( XtPointer client_data , XtIntervalId * id ) ;
238
main(int argc,char * argv[])239 int main( int argc , char * argv[] )
240 {
241 int iarg , ii , ny , ignore=0 , use=0 , install=0 ;
242 char * tsfile , * cpt ;
243 char dname[THD_MAX_NAME] , subv[THD_MAX_NAME] ;
244 MRI_IMAGE * flim ;
245 float * far ;
246 XtAppContext app ;
247 Widget shell ;
248 int cxx=0 , cyy=1 , cxsig=2 , cysig=3 , crho=4 ;
249
250 /*-- help? --*/
251
252 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
253 printf("Usage: 1dsigplot [options] infile\n"
254 "Scatterplots a 5 column *.1D file to the screen.\n"
255 "The columns of the input file must be\n"
256 " x y sigma_x sigma_y rho_xy\n"
257 "\n"
258 "Options:\n"
259 " -install = Install a new X11 colormap.\n"
260 " -ignore nn = Skip first 'nn' rows in the input file\n"
261 " [default = 0]\n"
262 " -use mm = Plot 'mm' points [default = all of them]\n"
263 " -xlabel aa = Put string 'aa' below the x-axis\n"
264 " [default = no axis label]\n"
265 " -ylabel aa = Put string 'aa' to the left of the y-axis\n"
266 " [default = no axis label]\n"
267 " -pell p = Set CDF probability contour level for ellipses.\n"
268 " -col abcde = 'abcde' is a permutation of '01234', and\n"
269 " specifies the column order for the data,\n"
270 " where a=column index for x\n"
271 " b=column index for y\n"
272 " c=column index for sigma_x\n"
273 " d=column index for sigma_y\n"
274 " e=column index for rho_xy\n"
275 " [default = 01234]\n"
276 "\n"
277 " -xrange x1 x2 = Range of x-axis\n"
278 " -yrange y1 y2 = Range of y-axis\n"
279 ) ;
280 PRINT_COMPILE_DATE ; exit(0) ;
281 }
282
283 machdep() ;
284
285 /* open X11 */
286
287 shell = XtVaAppInitialize(
288 &app , "AFNI" , NULL , 0 , &argc , argv , NULL , NULL ) ;
289 if( shell == NULL ){
290 fprintf(stderr,"** Cannot initialize X11!\n") ; exit(1) ;
291 }
292
293 cpt = my_getenv("TMPDIR") ; /* just for fun */
294
295 /*-- scan arguments that X11 didn't eat --*/
296
297 iarg = 1 ;
298 while( iarg < argc && argv[iarg][0] == '-' ){
299
300 if( strcmp(argv[iarg],"-xrange") == 0 ){
301 xlo = strtod(argv[++iarg],NULL) ;
302 xhi = strtod(argv[++iarg],NULL) ;
303 iarg++ ; continue ;
304 }
305
306 if( strcmp(argv[iarg],"-yrange") == 0 ){
307 ylo = strtod(argv[++iarg],NULL) ;
308 yhi = strtod(argv[++iarg],NULL) ;
309 iarg++ ; continue ;
310 }
311
312 if( strcmp(argv[iarg],"-pell") == 0 ){
313 pell = strtod(argv[++iarg],NULL) ;
314 iarg++ ; continue ;
315 }
316
317 if( strcmp(argv[iarg],"-col") == 0 ){
318 int ierr=0 ;
319 iarg++ ;
320 cxx = argv[iarg][0] - '0' ; if( cxx < 0 || cxx > 4 ) ierr++ ;
321 cyy = argv[iarg][1] - '0' ; if( cyy < 0 || cyy > 4 ) ierr++ ;
322 cxsig = argv[iarg][2] - '0' ; if( cxsig < 0 || cxsig > 4 ) ierr++ ;
323 cysig = argv[iarg][3] - '0' ; if( cysig < 0 || cysig > 4 ) ierr++ ;
324 crho = argv[iarg][4] - '0' ; if( crho < 0 || crho > 4 ) ierr++ ;
325 if( ierr || cxx+cyy+cxsig+cysig+crho != 10 ){
326 fprintf(stderr,"*** Illegal argument after -ord!\n");exit(1);
327 }
328 iarg++ ; continue ;
329 }
330
331 if( strcmp(argv[iarg],"-install") == 0 ){
332 install++ ; iarg++ ; continue ;
333 }
334
335 if( strcmp(argv[iarg],"-") == 0 ){ /* skip */
336 iarg++ ; continue ;
337 }
338
339 if( strcmp(argv[iarg],"-title") == 0 ){
340 title = argv[++iarg] ;
341 iarg++ ; continue ;
342 }
343
344 if( strcmp(argv[iarg],"-xlabel") == 0 ){
345 xlabel = argv[++iarg] ;
346 iarg++ ; continue ;
347 }
348
349 if( strcmp(argv[iarg],"-ylabel") == 0 ){
350 ylabel = argv[++iarg] ;
351 iarg++ ; continue ;
352 }
353
354 if( strcmp(argv[iarg],"-ignore") == 0 ){
355 ignore = strtod( argv[++iarg] , NULL ) ;
356 if( ignore < 0 ){fprintf(stderr,"** Illegal -ignore value!\n");exit(1);}
357 iarg++ ; continue ;
358 }
359
360 if( strcmp(argv[iarg],"-use") == 0 ){
361 use = strtod( argv[++iarg] , NULL ) ;
362 if( use < 2 ){fprintf(stderr,"** Illegal -use value!\n");exit(1);}
363 iarg++ ; continue ;
364 }
365
366 fprintf(stderr,"** Unknown option: %s\n",argv[iarg]) ; exit(1) ;
367 }
368
369 if( iarg >= argc ){
370 fprintf(stderr,"** No tsfile on command line!\n") ; exit(1) ;
371 }
372
373 dc = MCW_new_DC( shell , 16 ,
374 DEFAULT_NCOLOVR , INIT_colovr , INIT_labovr ,
375 1.0 , install ) ;
376
377 flim = mri_read_1D( argv[iarg] ) ;
378 if( flim == NULL ){
379 fprintf(stderr,"** Can't read input file %s\n",argv[iarg]) ;
380 exit(1);
381 }
382
383 if( flim->ny != 5 ){
384 fprintf(stderr,"** Input file doesn't have exactly 5 columns!\n") ;
385 exit(1) ;
386 }
387 far = MRI_FLOAT_PTR(flim) ;
388 nx = flim->nx ;
389
390 xx = far + cxx *nx ;
391 yy = far + cyy *nx ;
392 xsig = far + cxsig*nx ;
393 ysig = far + cysig*nx ;
394 xycor = far + crho *nx ;
395
396 nx = nx - ignore ; /* cut off the ignored points */
397
398 if( use > 1 && nx > use ) nx = use ;
399
400 /* start X11 */
401
402 (void) XtAppAddTimeOut( app , 123 , startup_timeout_CB , NULL ) ;
403
404 XtAppMainLoop(app) ;
405 exit(0) ;
406 }
407
408 /*-----------------------------------------------------------------*/
killfunc(void * fred)409 void killfunc(void * fred){ exit(0) ; }
410 /*-----------------------------------------------------------------*/
411
startup_timeout_CB(XtPointer client_data,XtIntervalId * id)412 void startup_timeout_CB( XtPointer client_data , XtIntervalId * id )
413 {
414 MEM_plotdata * mp ;
415
416 mp = PLOT_scatterellipse( nx , xlo,xhi,ylo,yhi ,
417 xx,yy,xsig,ysig,xycor ,
418 pell , xlabel,ylabel,title ) ;
419
420 if( mp != NULL )
421 (void) memplot_to_topshell( dc->display , mp , killfunc ) ;
422
423 return ;
424 }
425