From octave-maintainers-request at bevo dot che dot wisc dot edu Wed Jan 28 16:03:21 2004 Subject: contour plots and doing graphics in octave From: Donald J Bindner To: octave-maintainers at bevo dot che dot wisc dot edu Date: Wed, 28 Jan 2004 15:48:03 -0600 --2fHTh5uZTiUOsy+g Content-Type: text/plain; charset=us-ascii Content-Disposition: inline I've looked at the code for doing contour plots in octave (both version 2.0.x and 2.1.x) and both seem to be somewhat problematic because they place what is effectively a 2-D object into a 3-D context. One deficiency of this is that you can't really draw on a contour plot. For example in a Calculus class it might desirable to draw gradient vectors against the contours or show the path a particle would take if it followed the gradient up to a maximum. I have a modest alternative that generates good 2-D contour plots from within Octave rather than leaving that jobs to gnuplot. It is not as good as the current model in some cases, because I used linear pieces, and the current code specifies bsplines, but it is better because it isn't 3-D. What I have done is the most obvious simple thing. I interpolate the function on a mesh of points (like returned my meshgrid) and connect up the places that a contour would cross. I dump these line segments to files (one for each contour) and then ask gnuplot to plot the file for each contour in a different color. It works, so I know this can work pretty well and can be done quickly. But, at the moment, I'm calling gnuplot myself rather than tying in with the graphics primitives already in Octave. That means it is 2-D but I still can't draw on it. My question in a nutshell is this. From a dynamically linked function, how do I create graphics? There is a send_to_plot_stream() in pt-plot.cc but its symbol is not exported, except through the graw() function. I see that there are structures, tree_plot_command, etc., but I don't understand what to do with them. And I can't seem to trace from plot.m -> __plt__.m --> __plt2__.m --> to the point where it interfaces with the C++ code. I attach my working-but-not-finished prototype so you can see what context I come from. Example usage: octave:1> r = linspace(0,1,50); t = linspace(0,2*pi,50); octave:2> [rr,tt] = meshgrid(r,t); xx = rr .* cos(tt); yy = rr .* sin(tt); zz = sqrt(1-rr.^2); octave:3> contour(xx,yy,zz) -- Don Bindner --2fHTh5uZTiUOsy+g Content-Type: text/x-c++src; charset=us-ascii Content-Disposition: attachment; filename="__cntr__.cc" #include #include #include #include #include #include static std::string Vgnuplot_command_plot; static std::string Vgnuplot_command_replot; static oprocstream *plot_stream = 0; static int set_string_var (std::string& var, const char *nm) { int retval = 0; std::string s = builtin_string_variable (nm); if (s.empty ()) { gripe_invalid_value_specified (nm); retval = -1; } else var = s; return retval; } static int gnuplot_command_plot (void) { return set_string_var (Vgnuplot_command_plot, "gnuplot_command_plot"); } static int gnuplot_command_replot (void) { return set_string_var (Vgnuplot_command_replot, "gnuplot_command_replot"); } static double mmin( const Matrix &m ) { double d; int i,j,nr,nc; nr = m.rows(); nc = m.columns(); if( nr==0 || nc==0 ) return 0.0; d = m(0,0); for(i = 0; i < nr; i++ ) { for( j = 0; j < nc; j++ ) { d = (m(i,j)d) ? m(i,j) : d; } } return d; } DEFUN_DLD( __cntr__, args, , "__cntr__ (x,y,z,levels): plot contour lines of a function\n\n" "Plot contour lines by interpolating the values of a function\n" "on a grid. If x,y,z are matrices with the same dimensions\n" "the the corresponding elements represent points on the function.\n" "levels contains a vector of heights at which to draw level curves.\n" ) { octave_value_list retval; int nargin = args.length(); int nr, nc, nl, i, j, k; int flag=0; double x11, x21, x12, x22, y11, y21, y12, y22, z11, z21, z12, z22; double c, x, y; double minx, miny, maxx, maxy; std::string fname; FILE *f; char cmd[256]; if( nargin != 4 ) { print_usage( "cnt" ); return retval; } // Initialize these strings gnuplot_command_plot(); gnuplot_command_replot(); Matrix xx = args(0).matrix_value(); Matrix yy = args(1).matrix_value(); Matrix zz = args(2).matrix_value(); ColumnVector level = ColumnVector( args(3).vector_value()); nl = level.length(); nr = xx.rows(); nc = xx.columns(); minx = mmin(xx); maxx = mmax(xx); miny = mmin(yy); maxy = mmax(yy); if( plot_stream && *plot_stream ) { delete plot_stream; plot_stream = NULL; } plot_stream = new oprocstream( "gnuplot" ); if( !plot_stream || ! *plot_stream ) { error( "Could not open pipe." ); return retval; } for( k = 0; k < nl; k++ ) { c = level(k); fname = file_ops::tempnam ("", "oct-"); f = fopen( fname.c_str(), "w" ); mark_for_deletion( fname ); for( i = 0; i < nr-1; i++ ) { for( j = 0; j < nc-1; j++ ) { x11 = xx(i,j ); x21 = xx(i,j+1); x12 = xx(i+1,j); x22 = xx(i+1,j+1); y11 = yy(i,j ); y21 = yy(i,j+1); y12 = yy(i+1,j); y22 = yy(i+1,j+1); z11 = zz(i,j ); z21 = zz(i,j+1); z12 = zz(i+1,j); z22 = zz(i+1,j+1); // find bottom intersection with contour if(( z11>=c && z21=c )) { x = (c - z11)*(x21-x11)/(z21-z11) + x11; y = (c - z11)*(y21-y11)/(z21-z11) + y11; fprintf( f, "%lf %lf\n", x, y ); flag = 1; } // find top intersection with contour line if(( z12>=c && z22=c )) { x = (c - z12)*(x22-x12)/(z22-z12) + x12; y = (c - z12)*(y22-y12)/(z22-z12) + y12; fprintf( f, "%f %f\n", x, y ); flag = 1; } // find left intersection with contour line if(( z11>=c && z12=c )) { x = (c - z11)*(x12-x11)/(z12-z11) + x11; y = (c - z11)*(y12-y11)/(z12-z11) + y11; fprintf( f, "%f %f\n", x, y ); flag = 1; } // find right intersection with contour line if(( z21>=c && z22=c )) { x = (c - z21)*(x22-x21)/(z22-z21) + x21; y = (c - z21)*(y22-y21)/(z22-z21) + y21; fprintf( f, "%f %f\n", x, y ); flag = 1; } if( flag == 1 ) { fprintf( f, "\n" ); flag = 0; } } } fclose( f ); if( k == 0 ) { snprintf( cmd, sizeof(cmd), "%s [%lf:%lf] [%lf:%lf] \"%s\" title \"%g\" with lines 1\n", Vgnuplot_command_plot.c_str(), minx, maxx, miny, maxy, fname.c_str(), c ); *plot_stream << cmd; } else { snprintf( cmd, sizeof(cmd), "%s \"%s\" title \"%g\" with lines %d\n", Vgnuplot_command_replot.c_str(), fname.c_str(), c, k+1 ); *plot_stream << cmd; } } plot_stream->flush(); return retval; } --2fHTh5uZTiUOsy+g Content-Type: text/plain; charset=us-ascii Content-Disposition: attachment; filename="contour.m" function contour( x, y, z, levels ) # contour (x,y,z,levels): plot contour lines of a function # # Plot contour lines by interpolating the values of a function # on a grid. If x,y,z are matrices with the same dimensions # the the corresponding elements represent points on the function. # levels contains a vector of heights at which to draw level curves, # or a scalar containing the number of equally spaced curves to # draw. # # If x is a vector and y is a vector, then z must be a matrix # with the same number of rows as x and the same number of # columns as y. The corresponding points are x(i),y(j),z(i,y) # are points of the function. if( nargin < 3 ) usage( "contour" ); return; endif if( nargin == 3 ) levels = 5; endif if( is_scalar( levels ) ) minz = min(min(z)); maxz = max(max(z)); diff = (maxz-minz)/(2*levels); levels = linspace(minz+diff,maxz-diff,levels ); endif if( !is_vector( levels )) error( "levels must be a vector or a scalar" ); endif if( is_vector(x) && is_vector(y)) if( length(x) != rows(z) || length(y) != columns(z) ) usage( "contour" ); return; endif [x,y] = meshgrid(x,y); endif if( rows(x) != rows(z) || rows(y) != rows(z) || columns(x) != columns(z) || columns(y) != columns(z) ) usage( "contour" ); return; endif __cntr__(x,y,z,levels); endfunction --2fHTh5uZTiUOsy+g--