;+ ; NAME: ; SCATTER3D ; ; PURPOSE: ; Show a scatterplot of two arrays with color coding to ; indicate density of points. ; ; AUTHOR: ; ; Thomas Berger ; Lockheed Martin Solar and Astrophysics Laboratory ; berger@lmsal.com ; ; CATEGORY: ; Data analysis, graphics. ; ; CALLING SEQUENCE: ; SCATTER3D,array1,array2 ; ; INPUTS: ; ARRAY1: data array of any dimension > 0 ; ARRAY2: data array of the same dimension as ARRAY1 ; ; OUTPUTS: ; BINS: The binsize actually used in HIST2D (determined from ; data ranges and NBIN parameter. ; ; R: The HIST2D "image" used to display the scatterplot ; ; CENTER: The center of the principle moment ellipse. Only if ; keyword ELLIPSE is set. ; ; SEMI_AIXS: FLTARR[2] containing the semi-axes of the principle ; moment ellipse in data units. Only if keyword ELLIPSE ; is set. ; ; ROTANG: Rotation angle of the principle moment ellipse in ; degrees counterclockwise from x-axis. Only if ; keyword ELLIPSE is set. ; ; KEYWORD PARAMETERS: ; ; NBIN: The number of bins for HIST2D - effectively the ; number of "pixels" in the scatterplot image. ; Default is 100. ; ; ELLIPSE: If set, routine calculates and plots the principle moment ; ellipse fit to the scatterplot "image". ; ; XTITLE: Passed directly to PLOT.PRO for axis labelling. ; YTITLE ; ; PLOT45: If set, plots the 45-degree, slope=1, line for comparison ; to the scatterplot data. ; ; COMMON BLOCKS: ; None. ; ; SIDE EFFECTS: ; 720x540 pixel direct graphics window is opened in upper right of display. ; ; RESTRICTIONS: ; Requires HIST2D.PRO (built-in IDL function) and ; COLORBAR.PRO (Fanning Software IDL procedure, also found in SolarSoft). ; ; EXAMPLE: ; To make a scatterplot of two images IMAGE1 and IMAGE2 with ; the priniciple moment ellipse and a 45-degree line ; ; SCATTER3D, image1, image2, bins, r, center, semi, rotang,$ ; /ELLIPSE,/PLOT45,XTITLE='Image 1', YTITLE='Image 2' ; ; MODIFICATION HISTORY: ; Written by: Thomas Berger, LMSAL, 10/2001. ; 2-November-2001 - S.L.Freeland; ; change 'ellipse' to scatter3d_ellipse ; 'diagonal' to scatter3d_diagonal ; to avoid current (and future) naming conflicts ;- ; PRO scatter3d, dat1, dat2, bins, r, center, semi_axis, rotang, $ NBIN=nbin, ELLIPSE=ellipse, $ XTITLE=xtitle, YTITLE=ytitle, $ PLOT45=plot45 ON_ERROR,2 s1 = SIZE(dat1) s2 = SIZE(dat1) if s2(0) ne s1(0) OR s2(1) ne s1(1) OR s2(2) ne s1(2) then begin MESSAGE,'Dimensions of datasets unequal - returning.' RETURN end if KEYWORD_SET(xtitle) then xt = xtitle else xt = 'DATA1' if KEYWORD_SET(ytitle) then yt = ytitle else yt = 'DATA2' mx1 = FLOAT(MAX(dat1)) mn1 = FLOAT(MIN(dat1)) if KEYWORD_SET(nbin) then ng = nbin else ng=100. rng1 = mx1-mn1 b1 = ROUND(ng*rng1)/FLOAT(ng)/ng mx2 = FLOAT(MAX(dat2)) mn2 = FLOAT(MIN(dat2)) rng2 = mx2-mn2 b2 = ROUND(ng*rng2)/FLOAT(ng)/ng mins = [mn1,mn2] r = HIST_2D(dat1,dat2,BIN1=b1,BIN2=b2,MIN1=mn1,MAX1=mx1,MIN2=mn2,MAX2=mx2) szr = SIZE(r) ;actual binsize of x and y pixels in r: bins = [rng1/szr[1],rng2/szr[2]] ;Make the FIXED ASPECT RATIO plot WINX,720,540 LOADCT,2 ar = FLOAT(680)/640 x0 = 0.15 y0 = 0.10 x1 = 0.76 y1 = y0 + (x1-x0)*ar !p.position=[x0,y0,x1,y1] y1 = ((x1-x0)*!d.x_vsize)/((y1-y0)*!d.y_vsize)*y1 !p.position=[x0,y0,x1,y1] !x.margin=[0,4] !y.margin=[4,4] PLOT,[mn1,mx1],[mn2,mx2],/nodata,$ xrange=xrang,yrange=yrang,/xsty,/ysty px = !X.WINDOW*!D.X_VSIZE py = !Y.WINDOW*!D.Y_VSIZE sx = px[1]-px[0]+1 sy = py[1]-py[0]+1 TVSCL,CONGRID(ALOG(r+1),sx,sy),px[0],py[0] PLOT,[mn1,mx1],[mn2,mx2],/noerase,$ /nodata,$ xrange=xrang,yrange=yrang,/xsty,/ysty,$ title='Log(1+N) Histogram Scatterplot',xtitle=xt,ytitle=yt PLOTS,[mn1,mx1],[0,0] PLOTS,[0,0],[mn2,mx2] ;Plot 45-degree line for reference: if KEYWORD_SET(plot45) then PLOTS,[mn1,mx1],[mn1,mx1] ;Add colorbar. colorbar,position=[.90,.1,.94,y1],$ range=[Min(r),Max(r)],/vertical,$ format='(I6)',title='N Values',charsize=0.7 ;Ellipse fit if KEYWORD_SET(ellipse) then begin center = 0 semi_axis = 0 rotang = 0 scatter3d_ellipse,r,center,semi_axis,rotang center = center*bins-ABS(mins) PLOTS,center[0],center[1],$ psym=4,symsize=1,thick=3,color=0 a = semi_axis[0] b = semi_axis[1] xp1 = FINDGEN(4*a)/2 - a xp1 = [xp1,-xp1[0]] yp = b*SQRT(1.0 - (xp1/a)^2.) yn = -yp cosp = COS(!DTOR*rotang) sinp = SIN(!DTOR*rotang) rotmat = [[sinp,cosp],[sinp,-cosp]] xypr = rotmat#TRANSPOSE([[xp1],[yp]]) xynr = rotmat#TRANSPOSE([[xp1],[yn]]) ;Draw the ellipse in the plot box only. Should use vectors but got lazy... xpr = REFORM(xypr[0,*])*bins[0]+center[0] ypr = REFORM(xypr[1,*])*bins[1]+center[1] goodx = WHERE(xpr ge mn1 and xpr le mx1) xpr = xpr[goodx] ypr = ypr[goodx] goody = WHERE(ypr ge mn2 and ypr le mx2) xpr = xpr[goody] ypr = ypr[goody] PLOTS,xpr,ypr, col=0, thick=2, lines=0 goodx = WHERE(xynr[0,*] gt mn1 and xynr[0,*] lt mx1) xnr = REFORM(xynr[0,*])*bins[0]+center[0] ynr = REFORM(xynr[1,*])*bins[1]+center[1] goodx = WHERE(xnr ge mn1 and xnr le mx1) xnr = xnr[goodx] ynr = ynr[goodx] goody = WHERE(ynr ge mn2 and ynr le mx2) xnr = xnr[goody] ynr = ynr[goody] PLOTS,xnr,ynr, col=0, thick=2, lines=0 semi_axis = semi_axis*bins end RETURN END PRO scatter3d_ellipse,eabt,center,semi_axes,rotation ; ; eabt - 2-dimensional object array ; center - center of computed ellipse (pixels) (centroid of eabt) ; semi_axes - semi-axes of computed ellipse (pixels) ; rotation - rotation of computed ellipse wrt input array (degrees) ; ;Written by Bill Rosenberg 3/95. sz=size(eabt) if sz(0) ne 2 then begin print,' Dimension of input array, ',sz(0),' not equal to 2' retall endif xd=sz(1) & yd=sz(2) xr=findgen(xd,yd) mod xd yr=findgen(yd,xd) mod yd yr = transpose(yr) ;print,xr,format='(i3)' ;print,yr,format='(i3)' s=total(eabt) xbar=total(xr*eabt)/s ;print,xbar ybar=total(yr*eabt)/s ;print,ybar sxx=total(eabt*(xr-xbar)^2)/s ;print,sxx syy=total(eabt*(yr-ybar)^2)/s ;print,syy sxy=total(eabt*(xr-xbar)*(yr-ybar))/s ;print,sxy ;Covariance matrix: kk=[[sxx,sxy],[sxy,syy]] scatter3d_diagonal,kk,dd,t1 ;print,2.*sqrt(dd > 0) ;print,t1 rotation=t1 semi_axes=2.*sqrt([dd(0,0),dd(1,1)] > 0) center=[xbar,ybar] ;print,' Ellipse center: ',center ;print,' Ellipse semi-axes: ',semi_axes ;print,' Ellipse rotation: ',rotation,' degrees' RETURN END PRO scatter3d_diagonal,kk,dd,t1 ;Matrix diagonalization routine. Used by ELLIPSE.PRO ; Rosenberg 3/95. ;print,kk ; 50.0000 22.0000 ; 22.0000 10.0000 a=kk(0,0) & b=kk(0,1) & d=kk(1,1) ;print,a*d-b^2 ; 16.0000 l1=((a+d)+sqrt((a-d)^2+4.*b^2))/2. ;print,l1 ; 59.7321 l2=((a+d)-sqrt((a-d)^2+4.*b^2))/2. ;print,l2 ; 0.267862 t1=rtand(-b/(a-l1)) ;print,t1 ; 66.1368 t2=rtand(-b/(a-l2)) ;print,t2 ; -23.8632 x1=[sind(t1),cosd(t1)] x2=[sind(t2),cosd(t2)] pp=[[x1],[x2]] ;print,pp ; 0.914514 0.404554 ; -0.404554 0.914514 ppi=transpose(pp) dd=ppi#kk#pp ;print,dd ; 59.7321 2.942979e-06 ; 3.814697e-06 0.267863 RETURN END