PRO cafe_plot_model2, env, group, position=position, $ color=color, distinct=distinct, $ linestyle=linestyle, $ psym=psym, $ range=range, quiet=quiet, $ undefined=undefined, $ noselect=noselect, $ selectcolor=selectcolor, $ undefcolor=undefcolor, $ deltacolor=deltacolor, $ deltalinestyle=deltalinestyle, $ deltapsym=deltapsym, $ norange=norange, $ polyshade=polyshade, $ shade=shade, $ surf=surf, $ contour=contour, $ isolines=isolines, $ gridspacing=gridspacing, $ gridlimits=gridlimits, $ trisurf=trisurf, $ tv=tv,zbuff=zbuff, $ shading=shading, $ skipmargin=skipmargin, $ zvalue=zvalue, $ modelrange=modelrange, $ xrange=xrange, $ yrange=yrange, $ zrange=zrange, $ _EXTRA=ex, $ help=help,shorthelp=shorthelp ;+ ; NAME: ; plot_model2 ; ; PURPOSE: ; plot 2-d model over existing data ; ; CATEGORY: ; cafe ; ; SUBCATEGORY: ; plot type ; ; INPUTS: ; ; group - (optional) Define the data group to plot. ; Default is the current group. Must be in ; range [0..29]. ; ; PLOT OUTPUT: ; ; model2 - plots the data of group given . X/Y ranges are ; kept. Undefined points are not plotted by default. ; No subgroup distinction will be performed (the ; data points are plotted together). ; ; SETPLOT KEYWORDS: ; Apart from IDL shade_surf/surface keywords following special ; are defined: ; shade - Plot a shaded surface instead surface mesh. ; USING THIS OPTION WILL CLEAR ALL FORMER ; DRAWINGS! ; ; polyshade - Plot a shaded surface. Instead of shade the data ; will be plot as polygons with their ; positions. Usefull for very irregulary gridded ; data. Uses polyshade internally. ; USING THIS OPTION WILL CLEAR ALL FORMER ; DRAWINGS! ; ; surf - Plot a surface mesh. (default) ; ; trisurf - Plot a surface mesh consisting of ; triangles. The advantage in respect of surf is ; that (as in polyshade) the actual data point position is ; taken into account. The disadvantage is that no ; hidden line removal is performed, and plotting ; takes a long time for larger data sets. ; Triangles with orientations normal pointing not ; to the observer are not displayed (simplest ; form of hidden line removal) ; ; contour - Plot a contour plot. ; ; isolines - Plot on shade/surface isolevel contour lines ; ; tv - Plot pixel in a plane. Internally the procedure ; polyfill is used. This keyword works with zbuff ; only. Use the shading keyword to define the ; color distribution. ; ; gridspacing- Two component vector defining the grid spacing ; in x and y direction (refer TRIGRID, keyword ; GS). ; ; gridlimits - Four component vector defining the grid ; boundary as [x0,y0,x1,y1] (refer TRIGRID, ; keyword limits). ; ; shading - Expression string which defines how to color ; encode shaded plots. Applicable for ; shade/polyshade modes only. The expression ; string may be any arithmetic ; expression. Special keywords are 'x', 'y' and ; 'z' which will be used. ; The expression will be applied to all (x,y) ; pairs used in given grid. Default is the ; built-in light source shading (refer shades ; keyword for shade_surf procedure). ; The result of the expression should be within ; 0..1 (assuming that x,y,z are in this range ; also). ; Examples: ; > plot,model2[shading=z] ; -> color encode z values ; > plot,model2[shading=y*z] ; -> color encode diagonal to y/z (less usefull) ; > plot,model2[shading=200] ; -> monochrome red shading ; ; undefined- Plot model for undefined data points also. Useful for ; closer inspection. (No color distinction to ; defined data points) ; ; selected - Plot model for selected datapoints only. ; ; norange - Do not compute range (auxiliary data). Should ; not be used if single plot type. ; ; modelrange - Perform range determination for model size also ; (if range is not required anyway because still ; not defined). ; ; deltacolor - Color difference for distinct data sets. ; ; deltalinestyle - Increase of line styles for distinct data ; sets. Default: 0 (no change). ; ; deltapsym - Increase of point symbols for distinct data ; sets. May be positive or negative ; integer number. Default: 0 (no change). ; ; If psym exceeds -7..7 it will be set at 0 (no ; symbol). Therefore it is recomended to ; plot data sets with lines with a negative ; deltapsym (to decrease and restart from 0). ; ; selectcolor - Color for selected data points. ; ; ; SIDE EFFECTS: ; None. ; ; EXAMPLE: ; > data, bla.dat,dat2 ; > plot, model2 ; -> 2-dim data is displayed as shaded. ; ; HISTORY: ; $Id: cafe_plot_model2.pro,v 1.18 2004/07/30 17:17:26 goehler Exp $ ;- ; ; $Log: cafe_plot_model2.pro,v $ ; Revision 1.18 2004/07/30 17:17:26 goehler ; fix: allow negative deltacolor also without exceeding 0..255 range ; ; Revision 1.17 2004/06/21 17:23:40 goehler ; fix: modelrange was not visible from outside and thus without effect ; ; Revision 1.16 2004/06/08 08:28:37 goehler ; fix: range determination only when "modelrange" is set. ; ; Revision 1.15 2004/01/08 10:48:39 goehler ; - update concerning zvalue ; - model grid now without trigrid (faster..) ; ; Revision 1.14 2004/01/07 16:31:55 goehler ; - incorporate benefits from data2 ; - actually create model grid. this is still an intermediate version because ; still using trigrid function. ; ; Revision 1.13 2003/11/06 10:25:11 goehler ; update: allow psym/linestyle to vary for distinct data groups. ; documentation is updated. ; ; Revision 1.12 2003/11/05 09:03:55 goehler ; replace x/y/zrange with range expressions ; ; Revision 1.11 2003/11/04 14:27:37 goehler ; added help with shorthelp for all subtopics ; added qdp x binsize info (for plot also) ; ; Revision 1.10 2003/11/04 08:42:28 goehler ; implementation of trisurf hidden triangle removal ; ; Revision 1.9 2003/11/03 17:59:54 goehler ; added shading keyword ; ; Revision 1.8 2003/10/31 15:18:45 goehler ; added polyshade ; ; Revision 1.7 2003/05/08 09:20:13 goehler ; renamed color setting keywords/undefined (avoid name clash) ; ; Revision 1.6 2003/03/17 14:11:33 goehler ; review/documentation updated. ; ; Revision 1.5 2003/03/03 11:18:24 goehler ; major change: environment struct has become a pointer -> support of wplot/command line ; in common. ; Branch to be able to maintain the former line also. ; ; Revision 1.4 2003/02/18 13:56:22 goehler ; added isolines option for z-varying contour levels ; improved model2 ; ; Revision 1.3 2003/02/12 12:45:22 goehler ; surf->shade; added contour keyword ; ; Revision 1.2 2003/02/11 17:25:59 goehler ; improved axis display ; ; Revision 1.1 2003/02/11 15:03:35 goehler ; added gauss model, and method to plot models. ; ; ; ;; command name of this source (needed for automatic help) name="plot_model2" ;; ------------------------------------------------------------ ;; HELP ;; ------------------------------------------------------------ ;; if help given -> print the specification above (from this file) IF keyword_set(help) THEN BEGIN cafe_help,env, name return ENDIF ;; ------------------------------------------------------------ ;; SHORT HELP ;; ------------------------------------------------------------ IF keyword_set(shorthelp) THEN BEGIN cafereport,env, "model2 - plot 2-d model" return ENDIF ;; ------------------------------------------------------------ ;; SETUP ;; ------------------------------------------------------------ ;; set contour to enable iso lines IF keyword_set(isolines) THEN contour=1 ;; define default plotting: surface: IF NOT keyword_set(shade) AND $ NOT keyword_set(surf) AND $ NOT keyword_set(trisurf) AND $ NOT keyword_set(polyshade) AND $ NOT (keyword_set(tv) $ AND keyword_set(zbuff)) AND $ NOT keyword_set(contour) THEN surf=1 ;; ------------------------------------------------------------ ;; COLORS ;; ------------------------------------------------------------ ;; set default color IF n_elements(color) EQ 0 THEN color = 255 ;; define color delta IF n_elements(deltacolor) EQ 0 THEN deltacolor = 23 ;; define linestyle delta IF n_elements(deltalinestyle) EQ 0 THEN deltalinestyle = 0 ;; define psym delta IF n_elements(deltapsym) EQ 0 THEN deltapsym = 0 ;; color used for selected datapoints: IF n_elements(selectcolor) EQ 0 THEN selectcolor = 100 ;; ------------------------------------------------------------ ;; DEFINE RANGE ;; ------------------------------------------------------------ IF (n_elements(range) NE 0) AND keyword_set(modelrange) THEN BEGIN ;; for all subgroups: FOR i = 0, n_elements((*env).groups[group].data)-1 DO BEGIN ;; skip empty subgroups: IF NOT PTR_VALID((*env).groups[group].data[i].x) THEN CONTINUE ;; look only for defined datapoints: def_index = where(*(*env).groups[group].data[i].def) ;; skip undefined datasets IF def_index[0] EQ -1 THEN continue ;; replace range by superior value: IF def_index[0] NE -1 THEN BEGIN range[0] = range[0] < min((*(*env).groups[group].data[i].x)[def_index,0]) range[1] = range[1] > max((*(*env).groups[group].data[i].x)[def_index,0]) range[2] = range[2] < min((*(*env).groups[group].data[i].x)[def_index,1]) range[3] = range[3] > max((*(*env).groups[group].data[i].x)[def_index,1]) x = (*(*env).groups[group].data[i].x)[def_index,*] z = cafemodel(env,x,group) range[4] = range[4] < min(z) range[5] = range[5] > max(z) ENDIF ENDFOR ;; 2.) apply x/y/zrange if given IF n_elements(xrange) EQ 2 THEN range[0:1] = xrange IF n_elements(yrange) EQ 2 THEN range[2:3] = yrange IF n_elements(zrange) EQ 2 THEN range[4:5] = zrange ENDIF ;; do not plot if quiet: ;; (needed for range determination) IF keyword_set(quiet) THEN RETURN ;; ------------------------------------------------------------ ;; COLLECT DATA ;; ------------------------------------------------------------ ;; data to collect: x=0.D0 y=0.D0 z=0.D0 ;; plot all subgroups: FOR i = 0, n_elements((*env).groups[group].data)-1 DO BEGIN ;; skip empty subgroups: IF NOT PTR_VALID((*env).groups[group].data[i].x) THEN CONTINUE ;; plot only defined/selected datapoints: IF NOT keyword_set(selected) THEN $ def_index = where(*(*env).groups[group].data[i].def) $ ELSE $ def_index = where(*(*env).groups[group].data[i].selected) ;; check which data are not defined: undef_index = where(*(*env).groups[group].data[i].def EQ 0) x = [x,(*(*env).groups[group].data[i].x)[def_index,0]] y = [y,(*(*env).groups[group].data[i].x)[def_index,1]] z = [z,cafemodel(env,(*(*env).groups[group].data[i].x)[def_index,*],group)] ;; add undefined measure points ;; if requested IF (undef_index[0] NE -1) AND $ keyword_set(undef) THEN BEGIN x = [x,(*(*env).groups[group].data[i].x)[undef_index,0]] y = [y,(*(*env).groups[group].data[i].x)[undef_index,1]] z = [z,cafemodel(env,(*(*env).groups[group].data[i].x)[undef_index,*],group)] ENDIF ENDFOR ;; remove first dummy element: IF n_elements(x) GT 1 THEN BEGIN x =x[1:n_elements(x)-1] y =y[1:n_elements(y)-1] z =z[1:n_elements(z)-1] ENDIF ;; ------------------------------------------------------------ ;; CREATE GRID ;; ------------------------------------------------------------ ;; define default grid limits (idl complains if not doing so): IF n_elements(gridlimits) EQ 0 THEN $ gridlimits=[min(x),min(y),max(x),max(y)] ;; define default grid spacing (idl complains if not doing so): IF n_elements(gridspacing) EQ 0 THEN $ gridspacing = [(gridlimits[2]-gridlimits[0])/50, $ (gridlimits[3]-gridlimits[1])/50] ;; create triangulated list needed for trigrid/polyshade triangulate,x,y,tr,b ;; create grid axis: nxgrid = ceil((gridlimits[2]-gridlimits[0])/gridspacing[0])+1 ; number of x grid points nygrid = ceil((gridlimits[3]-gridlimits[1])/gridspacing[1])+1 ; number of y grid points xgrid = dindgen(nxgrid)*gridspacing[0]+gridlimits[0] ygrid = dindgen(nygrid)*gridspacing[1]+gridlimits[1] ;; create gridded data in dat: xgridpoints=reform(xgrid#make_array(nygrid,value=1.D0),nxgrid*nygrid) ygridpoints=reform(ygrid##make_array(nxgrid,value=1.D0),nxgrid*nygrid) gridpoints = [[xgridpoints],[ygridpoints]] dat = cafemodel(env,gridpoints,group) dat = reform(dat,nxgrid,nygrid,/overwrite) ; reform to make 2-d array ;; remove margin if desired: IF n_elements(skipmargin) NE 0 THEN BEGIN datsize=size(dat,/dimensions) dat=dat[skipmargin:datsize[0]-skipmargin,$ skipmargin:datsize[1]-skipmargin] xgrid=xgrid[skipmargin:datsize[0]-skipmargin] ygrid=ygrid[skipmargin:datsize[0]-skipmargin] ENDIF xrange=range[0:1] yrange=range[2:3] zrange=range[4:5] ;; make shure styles present: IF n_elements(xstyle) EQ 0 THEN xstyle=0 IF n_elements(ystyle) EQ 0 THEN ystyle=0 IF n_elements(zstyle) EQ 0 THEN zstyle=0 ;; ------------------------------------------------------------ ;; EVALUATE SHADING EXPRESSION ;; ------------------------------------------------------------ ;; create shades array: IF n_elements(shading) NE 0 THEN BEGIN ;; convert to string: aux_shading = string(shading) ;; replace z -> data in [0..1] aux_shading=strepex(aux_shading,"z", $ "((dat-zrange[0])/(zrange[1]-zrange[0]))",/all) ;; replace x -> xgrid in [0..1] aux_shading=strepex(aux_shading,"x", $ "((XGRID#MAKE_ARRAY(N_ELEMENTS(YGRID),VALUE=1.)-xrange[0])/(xrange[1]-xrange[0]))",/all) ;; replace y -> ygrid in [0..1] aux_shading=strepex(aux_shading,"y", $ "((YGRID##MAKE_ARRAY(N_ELEMENTS(XGRID),VALUE=1.)-yrange[0])/(yrange[1]-yrange[0]))",/all) ;; predefine shade array matching in dimensions: shades = make_array(dimension=size(dat,/dimension)) ;; execute full expression: IF NOT execute("shades[*,*] ="+aux_shading) THEN $ cafereport,env, !ERR_STRING ; plotting failed ;; cut to [0..255]: shades = ((shades < 1.) > 0.)*255. ENDIF ;; ------------------------------------------------------------ ;; PLOT POLYGON SHADED WITH IMPROVED OBJECT DISPLAY ;; ------------------------------------------------------------ IF keyword_set(polyshade) THEN BEGIN ;; add number of vertices for triangles :-) tr1 = make_array(4,(size(tr))[2],value=3) tr1[1:3,*] = tr ;; create shades array for shaded polygons: IF n_elements(shading) NE 0 THEN BEGIN ;; convert to string: aux_shading = string(shading) ;; replace x,y,z -> scaled versions of x,y,z aux_shading=strepex(aux_shading,"(x|y|z)", $ "((&0-&0range[0])/(&0range[1]-&0range[0]))",/all) ;; predefine shade array matching in dimensions: polyshades = make_array(dimension=size(z,/dimension)) ;; execute full expression: IF NOT execute("polyshades[*] ="+aux_shading) THEN $ cafereport,env, !ERR_STRING ; plotting failed ;; cut to [0..255]: polyshades = ((polyshades < 1.) > 0.)*255. ENDIF ;; create image (distinguish shades/no shades): IF n_elements(shading) NE 0 THEN $ img = polyshade(x,y,z,tr1,/t3d,shades=polyshades) $ ELSE img = polyshade(x,y,z,tr1,/t3d) ;; display image if no z-buffer: IF !D.name NE 'Z' THEN tvscl,img ;; add axis surface,dat,xgrid,ygrid, /noerase, /nodata, $ position=position,/save, $ xrange=xrange,yrange=yrange,zrange=zrange, $ xstyle=xstyle,ystyle=ystyle,zstyle=zstyle, $ zvalue=zvalue,_EXTRA=ex ENDIF ;; ------------------------------------------------------------ ;; PLOT SHADED ;; ------------------------------------------------------------ IF keyword_set(shade) THEN BEGIN ;; shaded drawing ;; distinction of shading due idl 5.4 bug: IF n_elements(shading) NE 0 THEN $ shade_surf,dat,xgrid,ygrid, /noerase,position=position,/save,$ xrange=xrange,yrange=yrange,zrange=zrange, $ shades=shades,color=color, $ zvalue=zvalue, _EXTRA=ex $ ELSE $ shade_surf,dat,xgrid,ygrid, /noerase,position=position,/save,$ xrange=xrange,yrange=yrange,zrange=zrange, $ zvalue=zvalue,color=color,_EXTRA=ex ENDIF ;; ------------------------------------------------------------ ;; PLOT SURFACE MESH ;; ------------------------------------------------------------ IF keyword_set(surf) THEN BEGIN ;; surface drawing: ;; disable axis: xstyle=xstyle OR 4 ystyle=ystyle OR 4 zstyle=zstyle OR 4 surface,dat,xgrid,ygrid, /noerase, position=position,/save,$ xrange=xrange,yrange=yrange,zrange=zrange, $ xstyle=xstyle,ystyle=ystyle,zstyle=zstyle,color=color, $ zvalue=zvalue,_EXTRA=ex ENDIF ;; ------------------------------------------------------------ ;; PLOT TRIANGLE MESH ;; ------------------------------------------------------------ IF keyword_set(trisurf) THEN BEGIN ;; surface drawing: ;; disable axis: xstyle=xstyle OR 4 ystyle=ystyle OR 4 zstyle=zstyle OR 4 FOR i = 0, n_elements(tr[0,*])-1 DO BEGIN ;;----------------------------------------------- ;; remove triangles with orientation making them invisible: ;; convert triangle position to homogenuous vectors ;; transformed into viewer system: p0 = [x[tr[0,i]],y[tr[0,i]],z[tr[0,i]],1] # !P.T p1 = [x[tr[1,i]],y[tr[1,i]],z[tr[1,i]],1] # !P.T p2 = [x[tr[2,i]],y[tr[2,i]],z[tr[2,i]],1] # !P.T ;; get triangle vertice vectors (differences) d1 = p1[0:2]/p1[3] - p0[0:2]/p0[3] d2 = p2[0:2]/p2[3] - p0[0:2]/p0[3] ;; positive triangle orientation -> show triangle IF determ([[0,0,1],[d1],[d2]]) gt 0 THEN BEGIN plots,x[tr[0:1,i]],y[tr[0:1,i]],z[tr[0:1,i]],/t3d,color=color,_EXTRA=ex plots,x[tr[1:2,i]],y[tr[1:2,i]],z[tr[1:2,i]],/t3d,color=color,_EXTRA=ex plots,x[tr[[0,2],i]],y[tr[[0,2],i]],z[tr[[0,2],i]],color=color,/t3d,_EXTRA=ex ENDIF ENDFOR ENDIF ;; ------------------------------------------------------------ ;; PLOT TV (PIXEL) ;; ------------------------------------------------------------ IF keyword_set(tv) AND keyword_set(zbuff) THEN BEGIN ;; define image range: xr = [xgrid[0],xgrid[n_elements(xgrid)-1]] yr = [ygrid[0],ygrid[n_elements(ygrid)-1]] xpoly=[xr[0],xr[1],xr[1],xr[0]] ypoly=[yr[0],yr[0],yr[1],yr[1]] IF n_elements(zvalue) NE 0 THEN BEGIN zpoly=make_array(4,value=zvalue*(zrange[1]-zrange[0])+zrange[0]) ENDIF ELSE BEGIN zpoly=make_array(4,value=zrange[1]) ENDELSE ;; use shading if existing: IF n_elements(shades) NE 0 THEN tvdat = shades $ ELSE tvdat=(dat-min(dat))/(max(dat)-min(dat))*255. ;; define image size: tvsize=size(tvdat,/dimensions)-1 ;; fill the image: polyfill, xpoly,ypoly,zpoly,/t3d,/data, $ _EXTRA=ex,pattern=tvdat, $ image_coord=[[0,0],[tvsize[0],0], $ [tvsize[0],tvsize[1]],[0,[tvsize[1]]]] ENDIF ;; ------------------------------------------------------------ ;; PLOT DATA ;; ------------------------------------------------------------ ;; plot all subgroups: FOR i = 0, n_elements((*env).groups[group].data)-1 DO BEGIN ;; skip empty subgroups: IF NOT PTR_VALID((*env).groups[group].data[i].x) THEN CONTINUE ;; plot only defined datapoints properly: def_index = where(*(*env).groups[group].data[i].def) undef_index = where(*(*env).groups[group].data[i].def EQ 0) select_index = where(*(*env).groups[group].data[i].selected $ AND *(*env).groups[group].data[i].def) ;; add undefined measure points, but do not connect with lines ;; if requested IF (undef_index[0] NE -1) AND $ keyword_set(undef) THEN $ plots, (*(*env).groups[group].data[i].x)[undef_index,0], $ (*(*env).groups[group].data[i].x)[undef_index,1], $ (*(*env).groups[group].data[i].y)[undef_index], $ /t3d,psym=4, color = undefcolor, _EXTRA=ex ;; show selected measure points, but do not connect with lines ;; if not rejected IF (select_index[0] NE -1) AND $ (NOT keyword_set(noselect)) THEN BEGIN plots, (*(*env).groups[group].data[i].x)[select_index,0], $ (*(*env).groups[group].data[i].x)[select_index,1], $ (*(*env).groups[group].data[i].y)[select_index], $ /t3d,psym=4, color = selectcolor, _EXTRA=ex ENDIF ;; change style when each data set is plotted different: IF keyword_set(distinct) THEN BEGIN ;; change color if some given: color = color - deltacolor ;; wrap around if too much colors: IF color LT 0 THEN color = color + 255 IF color GT 255 THEN color = color MOD 256 ;; change linestyle if different: linestyle = linestyle + deltalinestyle ;; wrap around linestyles: IF linestyle LT 0 THEN linestyle = 5 IF linestyle GT 5 THEN linestyle = 0 ;; change psym if different: psym = psym + deltapsym ;; wrap around psym, respect psym sign, allow psym > 7 if ;; not changing: IF deltapsym NE 0 AND abs(psym) GT 7 THEN psym = 0 ENDIF ENDFOR ;; ------------------------------------------------------------ ;; PLOT CONTOUR ;; ------------------------------------------------------------ IF keyword_set(contour) THEN BEGIN ;; set z value if not plotting iso lines: IF NOT keyword_set(isolines) AND n_elements(zvalue) EQ 0 THEN zvalue =1.0 contour,dat,xgrid,ygrid,/t3d, zvalue=zvalue, $ /noerase,position=position,color=color, $ xrange=xrange,yrange=yrange,zrange=zrange, $ _EXTRA=ex ENDIF END