slang-histogram 0.3.2a-3.1 / usr / share / doc / slang-histogram / examples / evt2img

#!/usr/bin/env isis-script
% -*- mode: slang; mode:fold -*-
%add_to_isis_load_path ("/nfs/cxc/a1/src/slgtk/devel/share/slsh/local-packages");

require ("histogram");
require ("rgbwidget");
require("gtkplot");

implements ("rgbselect"); %{{{ 

#ifnexists rgbselect

typedef struct { r0, r1, g0, g1, b0, b1 } RGBLimits;
static variable Scale_Funcs = [ &gtk_plot_set_xscale, &gtk_plot_set_yscale];
static variable Range_Funcs = [ &_gtk_plot_set_xrange, &_gtk_plot_set_yrange];
static variable Region = struct {color, lslider, rslider, canvas, widget, sd};
static variable Slider = struct {lneighbor,rneighbor,widget,parent,inverted};
static variable SelectorData = struct { 
   	x, y,
	x_non_zeroes,
	angstroms,
	npoints,
	scale,				% x/y log/lin
	rgblims,
	wavelength_plotd,
	energy_plotd,
	current_plotd,
	plots_area,
	apply_cb,
	apply_cb_args,
	sliders,
	Red,
	Green,
	Blue,
	status };

#ifnexists hist_bsearch
static define hist_bsearch (x, xs)
{
   variable n0 = 0;
   variable n = length (xs);
   variable n1 = n;
   
   while (n1 > n0 + 1)
     {
	variable n2 = (n1 + n0)/2;
	if (xs[n2] >= x)
	  {
	     n1 = n2;
	     continue;
	  }
	n0 = n2;
     }
   return n0;
}
#endif

static define index_of(X, value) { hist_bsearch(value,X); }

static define keV_to_A (e)
{
   return 12.3984185734/e;
}

static define data_to_canvas(pd,dx,dy)
{
   % Map values in (input) data space to positions in canvas (drawing) space
   variable px, py, cx, cy;
   gtk_plot_get_pixel(pd.plot, dx, dy, &px, &py);
   gtk_plot_canvas_get_position(pd.canvas, int(px), int(py), &cx, &cy);
   return (cx,cy);
}

static define region_draw(r)
{
   if (r.widget != NULL) {
	() = gtk_plot_canvas_remove_child(r.canvas,r.widget);
	r.widget = NULL;
   }

   variable span = [ int ( gtk_range_get_value (r.lslider.widget)) :
			int( gtk_range_get_value(r.rslider.widget)) ];
   if (length(span) == 1) return;

   variable pd = r.sd.current_plotd;
   variable yaxis = pd.dsets[0].axes[1], x = pd.dsets[0].axes[0].data;
   variable rx = [ x[ span[0] ],  x[span], x[ span[-1] ] ];
   variable ry = [ yaxis.min,  yaxis.data[span], yaxis.min ];

   % Would like to use array_map, but can't
   _for(0,length(rx)-1,1) {
	variable i = ();
	(rx[i], ry[i]) = data_to_canvas(pd,rx[i],ry[i]);
   }

   r.canvas = pd.canvas;
   r.widget = _gtk_plot_canvas_put_polygon(r.canvas, rx, ry,
	 		GTK_PLOT_LINE_SOLID, 0, gdk_black,  r.color, TRUE);
}

static define update_regions(sd)
{
   region_draw(sd.Red);
   region_draw(sd.Green);
   region_draw(sd.Blue);
   _gtk_plot_redraw(sd.current_plotd);
}

static define update_slider(slider, sd)
{
   variable value = int(gtk_range_get_value(slider.widget));

   variable neighbor = sd.sliders[slider.lneighbor];
   if (neighbor != NULL) {
	if (value <= int(gtk_range_get_value(neighbor.widget)))
		gtk_range_set_value(neighbor.widget,value);
   }

   neighbor = sd.sliders[slider.rneighbor];
   if (neighbor != NULL) {
	if (value >= int(gtk_range_get_value(neighbor.widget)))
	   gtk_range_set_value(neighbor.widget,value);
   }

   update_regions(sd);
}

static define format_slider_value(slider, value, sd)
{
   variable X = sd.current_plotd.dsets[0].axes[0].data;
   return sprintf("%3.2f", X [ int(value) ]);
}

static define set_all_slider_values(sd,r0,r1,g0,g1,b0,b1)
{
   gtk_range_set_value(sd.sliders[1].widget, r0);
   gtk_range_set_value(sd.sliders[2].widget, r1);
   gtk_range_set_value(sd.sliders[3].widget, g0);
   gtk_range_set_value(sd.sliders[4].widget, g1);
   gtk_range_set_value(sd.sliders[5].widget, b0);
   gtk_range_set_value(sd.sliders[6].widget, b1);
   update_regions(sd);
}

static define make_slider(sd, container, value, lneighbor, rneighbor)
{
   variable s  = @Slider;
   s.lneighbor = lneighbor;
   s.rneighbor = rneighbor;
   s.parent    = container;
   s.inverted  = FALSE;

   variable adj = gtk_adjustment_new(value, 0, sd.npoints - 1, 1, 1, 0.0);
   s.widget = gtk_hscale_new(adj);
   () = g_signal_connect (s.widget, "format_value", &format_slider_value, sd);
   gtk_range_set_update_policy (s.widget, GTK_UPDATE_DISCONTINUOUS);
   gtk_container_add(container,s.widget);
   return s;
}

static define make_region(selector_data, color,left_slider,right_slider)
{ 
   variable r	= @Region;
   r.lslider	= selector_data.sliders[left_slider];
   r.rslider	= selector_data.sliders[right_slider];
   r.color	= color;
   r.sd         = selector_data;

   gtk_widget_modify_bg(r.lslider.widget,GTK_STATE_NORMAL,color);
   gtk_widget_modify_bg(r.rslider.widget,GTK_STATE_NORMAL,color);

   () = g_signal_connect_swapped(r.lslider.widget,"value-changed",
   				&update_slider, r.lslider, selector_data);
   () = g_signal_connect_swapped(r.rslider.widget,"value-changed",
   				&update_slider, r.rslider, selector_data);
   return r;
}

static define set_rgb_values(rgblims, sd)
{
   rgblims.r0 = sd.x [ int( gtk_range_get_value(sd.sliders[1].widget) ) ];
   rgblims.r1 = sd.x [ int( gtk_range_get_value(sd.sliders[2].widget) ) ];
   rgblims.g0 = sd.x [ int( gtk_range_get_value(sd.sliders[3].widget) ) ];
   rgblims.g1 = sd.x [ int( gtk_range_get_value(sd.sliders[4].widget) ) ];
   rgblims.b0 = sd.x [ int( gtk_range_get_value(sd.sliders[5].widget) ) ];
   rgblims.b1 = sd.x [ int( gtk_range_get_value(sd.sliders[6].widget) ) ];
}

static define apply(win, sd, dismiss)
{
   if (dismiss) {
	sd.status = 0;
	set_rgb_values(sd.rgblims,sd);
	gtk_widget_destroy(win);
	gtk_main_quit();
	return;
   }

   if (sd.apply_cb != NULL) {
	variable r = @RGBLimits;
	set_rgb_values(r,sd);
	(@sd.apply_cb)(r, __push_args(sd.apply_cb_args) );
   }
}

static define equalize(sd)
{
   % Compute total area of spectrum covered by color regions,
   % then assign each color region an equal slice. 

   variable rmin = int(gtk_range_get_value(sd.sliders[1].widget));
   variable bmax = int(gtk_range_get_value(sd.sliders[6].widget));
   if ( bmax - rmin <= 3) {
	rmin = 0;
	bmax = sd.npoints - 1;
   }

   variable counts = sd.y[ [ rmin : bmax ] ];
   variable total_color = cumsum(counts), slice = total_color [-1] / 3.0;

   variable rmax = rmin + hist_bsearch(slice, total_color);
   variable gmax = rmin + hist_bsearch(2*slice, total_color);

   set_all_slider_values(sd,rmin,rmax,rmax,gmax,gmax,bmax);
}

static define revert(sd)
{
   variable r = sd.rgblims, x = sd.x;
   set_all_slider_values(sd, index_of(x,r.r0), index_of(x,r.r1),
	index_of(x,r.g0), index_of(x,r.g1), index_of(x,r.b0), index_of(x,r.b1));
}

static define set_scale(which_axis,what_scale,sd)
{
   variable func, plotd = sd.current_plotd;
   variable axis = plotd.dsets[0].axes[which_axis];

   if (what_scale == GTK_PLOT_SCALE_LOG10) {
	if (axis.min <= 0.0) {
	   () = printf("Log scaling disabled: some values are non-positive\n");
	   return FALSE;
	}
      gtk_plot_axis_set_ticks(plotd.plot, which_axis , 1, 8);
   }
   else
     {
	func = Range_Funcs[which_axis];
	@func(plotd, axis.min, axis.max);       % this will reset ticks, too
     }

   sd.scale[which_axis] = what_scale;
   func = Scale_Funcs[which_axis];		% this requires 2 LOC due
   (@func)(plotd.plot, what_scale);		% to a bug in S-Lang
   update_regions(sd);

   return TRUE;
}

static define set_scale_cb(radiob, which_axis, what_scale, default_radiob, sd)
{
   if (0 == gtk_toggle_button_get_active (radiob)) return;

   !if (set_scale(which_axis,what_scale,sd))
	gtk_toggle_button_set_active (default_radiob,TRUE);
}

static define plot_expose(event,sd) { update_regions(sd); return FALSE; }

static define customize_plotd(plotd,xlabel,sd)
{
   gtk_plot_canvas_unset_flags(plotd.canvas,GTK_PLOT_CANVAS_CAN_SELECT);
   gtk_plot_data_set_connector(plotd.dset.widget, GTK_PLOT_CONNECT_STRAIGHT);
   gtk_plot_axis_set_title(plotd.plot,GTK_PLOT_AXIS_LEFT,"Counts");
   gtk_plot_axis_set_title(plotd.plot,GTK_PLOT_AXIS_BOTTOM,xlabel);
   () = g_signal_connect_swapped(plotd.plot,"expose_event",&plot_expose,sd);
}

static define set_units(radiob, default_radiob, sd)
{
   if (0 == gtk_toggle_button_get_active (radiob)) return;

   variable newpd, label = gtk_button_get_label (radiob);
   switch(label)
     { case "Energy" : newpd = sd.energy_plotd; }

     { case "Wavelength" :

	if (sd.x_non_zeroes == NULL)
	   sd.x_non_zeroes = where(sd.x != 0.0);

	if (length(sd.x_non_zeroes) == 0) {
	   () = printf("Cannot plot Wavelength, Energy contains only zeroes\n");
	   gtk_toggle_button_set_active (default_radiob,TRUE);
	   return;
	}

	if (sd.wavelength_plotd == NULL) {
	   sd.angstroms = keV_to_A ( sd.x [ sd.x_non_zeroes ] );
	   variable pd = _gtk_plot( sd.angstroms, sd.y [sd.x_non_zeroes] );
	   customize_plotd(pd,label,sd);
	   gtk_box_pack_start(sd.plots_area, pd.canvas, TRUE, TRUE, 1);
	   sd.wavelength_plotd = pd;
	}
	newpd = sd.wavelength_plotd;
   }

   if (newpd == sd.current_plotd) return;

   gtk_widget_hide(sd.current_plotd.canvas);
   sd.current_plotd = newpd;
   ()=set_scale(0,sd.scale[0],sd);
   ()=set_scale(1,sd.scale[1],sd);
   gtk_widget_show_all(newpd.canvas);

   _for (1,6,1) {
	variable i = (), slider = sd.sliders[i];
	gtk_box_reorder_child (slider.parent, slider.widget,
					(i + slider.inverted) mod 2);
	slider.inverted = not(slider.inverted);
	gtk_range_set_inverted(slider.widget,slider.inverted);
	gtk_widget_queue_draw(slider.widget);
   }
}

static define make_radiob()
{
   variable box, widget, label, callback, callback_args;
   callback_args = __pop_args (_NARGS - 4);
   (box, widget, label, callback) = ();
   
   variable rb = gtk_radio_button_new_with_label_from_widget (widget, label);
   () = g_signal_connect (rb, "clicked", callback, __push_args(callback_args));
   gtk_box_pack_start (box, rb, 1, 1, 0);

   return rb;
}

static define make_aligned_button()
{
   variable vbox, text, cb, cbargs = __pop_args(_NARGS - 3);
   (vbox, text, cb) = ();

   variable align = gtk_alignment_new(0.5, 0.5, 0, 0);
   gtk_box_pack_start(vbox,align,FALSE,FALSE,5);
   variable button = gtk_button_new_with_label(text);
   gtk_container_add(align,button);
   if (cb != NULL)
      () = g_signal_connect_swapped(button,"clicked", cb, __push_args(cbargs) );
   else
      gtk_widget_set_sensitive(button,FALSE);
}

static define make_axis_scaler(vbox, which_axis, sd)
{
   variable frame = gtk_frame_new( char(88 + which_axis) + " Axis Scale");
   gtk_box_pack_start(vbox,frame,FALSE,FALSE,5);

   variable rvbox = gtk_vbox_new(FALSE,5);
   gtk_container_add(frame,rvbox);

   variable linrb = make_radiob(rvbox, NULL, "Linear", &set_scale_cb,
   			which_axis, GTK_PLOT_SCALE_LINEAR,  NULL, sd);

   () = make_radiob(rvbox, linrb,"Log10", &set_scale_cb, which_axis,
			GTK_PLOT_SCALE_LOG10,  linrb, sd);
}

public define rgbselect()
{
   variable sd = @SelectorData;

   switch(_NARGS) 
        { case 0 or case 1 or case 2:
		usage("Integer_Type = rgbselect(x, y, rgb_struct "+
			"[,apply_callback_func [, apply_callback_args] );");
	}
	{ case 3: (sd.x, sd.y, sd.rgblims) = (); }
	{
	   	sd.apply_cb_args = __pop_args(_NARGS - 4);
		(sd.x, sd.y, sd.rgblims, sd.apply_cb) = ();
	}

   if ( andelse {sd.apply_cb != NULL} {not(_is_callable(sd.apply_cb))} )
	error("Specified Apply callback is not a callable function");

   sd.x = typecast(sd.x,Double_Type);
   sd.y = typecast(sd.y,Double_Type);
   sd.status = -1;

   variable plotd = _gtk_plot(sd.x, sd.y);	% plot energy vs counts
   customize_plotd(plotd,"Energy",sd);

   sd.npoints = length(sd.x);
   sd.energy_plotd  = plotd;
   sd.current_plotd = plotd;

   variable win = gtk_window_new (GTK_WINDOW_TOPLEVEL);
   () = g_signal_connect(win,"destroy",&gtk_main_quit);
   gtk_window_set_title(win,"rgbselect 0.91");

   variable main_vbox = gtk_vbox_new(FALSE,5);
   gtk_container_add(win,main_vbox);

   sd.plots_area = gtk_hbox_new(FALSE,5);
   gtk_box_pack_start(main_vbox,sd.plots_area,TRUE,TRUE,0);
   gtk_box_pack_start(sd.plots_area,plotd.canvas,TRUE,TRUE,1);

   variable button_vbox = gtk_vbox_new(FALSE,20);
   gtk_box_pack_end(sd.plots_area,button_vbox,FALSE,FALSE,5);

   make_axis_scaler(button_vbox, GTK_PLOT_AXIS_X, sd);
   make_axis_scaler(button_vbox, GTK_PLOT_AXIS_Y, sd);

   variable frame = gtk_frame_new("Units");
   gtk_box_pack_start(button_vbox,frame,FALSE,FALSE,0);
   variable vbox = gtk_vbox_new(FALSE,10);
   gtk_container_add(frame,vbox);

   variable button = make_radiob(vbox,NULL,"Energy", &set_units, NULL, sd);
   button = make_radiob(vbox,button,"Wavelength", &set_units, button, sd);

   frame = gtk_frame_new("RGB Selectors");
   gtk_box_pack_start(main_vbox,frame,FALSE,FALSE,0);

   variable slider_vbox = gtk_vbox_new(FALSE,10);
   gtk_container_add(frame,slider_vbox);

   variable hbox = gtk_hbox_new(FALSE,10);
   gtk_box_pack_start(slider_vbox,hbox,FALSE,FALSE,0);

   sd.scale = [ GTK_PLOT_SCALE_LINEAR, GTK_PLOT_SCALE_LINEAR];

   sd.sliders = Struct_Type[8];
   sd.sliders[0] = NULL;
   sd.sliders[1] = make_slider(sd,hbox,index_of(sd.x,sd.rgblims.r0),0,2);
   sd.sliders[2] = make_slider(sd,hbox,index_of(sd.x,sd.rgblims.r1),1,3);

   hbox = gtk_hbox_new(FALSE,10);
   gtk_container_add(slider_vbox,hbox);
   sd.sliders[3] = make_slider(sd,hbox,index_of(sd.x,sd.rgblims.g0),2,4);
   sd.sliders[4] = make_slider(sd,hbox,index_of(sd.x,sd.rgblims.g1),3,5);

   hbox = gtk_hbox_new(FALSE,10);
   gtk_container_add(slider_vbox,hbox);
   sd.sliders[5] = make_slider(sd,hbox,index_of(sd.x,sd.rgblims.b0),4,6);
   sd.sliders[6] = make_slider(sd,hbox,index_of(sd.x,sd.rgblims.b1),5,7);
   sd.sliders[7] = NULL;

   hbox = gtk_hbox_new(FALSE,5);
   gtk_box_pack_end(main_vbox,hbox,FALSE,FALSE,5);
   button = gtk_button_new_from_stock(GTK_STOCK_APPLY);
   () = g_signal_connect_swapped(button,"clicked", &apply, win, sd, 0);
   gtk_box_pack_start(hbox,button,FALSE,FALSE,5);

   button = gtk_button_new_from_stock(GTK_STOCK_OK);
   () = g_signal_connect_swapped(button,"clicked", &apply, win, sd, 1);
   gtk_box_pack_start(hbox,button,FALSE,FALSE,5);

   button = gtk_button_new_from_stock(GTK_STOCK_CANCEL);
   () = g_signal_connect_swapped(button,"clicked",&gtk_widget_destroy, win);
   gtk_box_pack_end(hbox,button,FALSE,FALSE,5);

   sd.Red = make_region(sd, gdk_red, 1, 2);
   sd.Green = make_region(sd, gdk_green, 3, 4);
   sd.Blue = make_region(sd, gdk_blue, 5, 6);

   frame = gtk_frame_new("Color Ratios");
   gtk_box_pack_start(button_vbox,frame,FALSE,FALSE,0);
   vbox = gtk_vbox_new(FALSE,5);
   gtk_container_add(frame,vbox);

   make_aligned_button(vbox,"Equalize",&equalize,sd);
   make_aligned_button(vbox," Revert ",&revert,sd);

   gtk_widget_show_all(win);
   gtk_main();

   return sd.status;
}

#ifexists provide
provide("rgbselect");
#endif

#endif

%}}}

implements ("viewrgb"); %{{{
%require ("histogram");
%require ("gtkdevel");

%() = evalfile ("rgbwidget.sl");
%() = evalfile ("./rgbselect.sl");

_debug_info = 1;
_traceback=1;


static define create_rgb (spectrum_info)
{
   variable xgrid = spectrum_info.xgrid;
   variable ygrid = spectrum_info.ygrid;
   variable energy = spectrum_info.energy;
   variable x = spectrum_info.x;
   variable y = spectrum_info.y;
   variable rgb_limits = spectrum_info.rgb_limits;

   variable i_low = where ((energy >= rgb_limits.r0) and (energy < rgb_limits.r1));
   variable i_med = where ((energy >= rgb_limits.g0) and (energy < rgb_limits.g1));
   variable i_hi = where ((energy >= rgb_limits.b0) and (energy < rgb_limits.b1));

   variable r = hist2d (y[i_low], x[i_low], ygrid, xgrid);
   variable g = hist2d (y[i_med], x[i_med], ygrid, xgrid);
   variable b = hist2d (y[i_hi], x[i_hi], ygrid, xgrid);

   return (r, g, b);
}

static define save_to_png (rgb_widget, file)
{
   gdk_pixbuf_save (rgb_widget.zoomed_pixbuf, file, "png");
   %gdk_pixbuf_save (rgb_widget.pixbuf, file+".unzoomed", "png");
}

static define save_to_jpg (rgb_widget, file)
{
   gdk_pixbuf_save (rgb_widget.zoomed_pixbuf, file, "jpeg");
}

% Usage: make_radio_button (box, group, label, callback, ...)
static define make_radio_button ()
{
   variable box, group, label, callback, callback_args;
   callback_args = __pop_args (_NARGS-4);
   (box, group, label, callback) = ();
   
   variable button = gtk_radio_button_new_with_label (group, label);
   () = g_signal_connect (button, "clicked",
			  callback, __push_args(callback_args));
   gtk_box_pack_start (box, button, 1, 1, 0);
   group = gtk_radio_button_get_group (button);

   return (button, group);
}

static define set_save_filename_ext (w, ext)
{
   variable file = gtk_file_selection_get_filename (w);
   
   if (file == NULL)
     return;
   
   variable st = stat_file (file);
   if (st != NULL)
     {
	if (stat_is ("dir", st.st_mode))
	  return;
     }

   variable dir = path_dirname (file);
   file = path_sans_extname (path_basename (file));

   if (file == "")
     return;

   file += ext;
   file = path_concat (dir, file);
   gtk_file_selection_set_filename (w, file);
}

static define save_as_png_callback (w, save_info)
{
   set_save_filename_ext (save_info.window, ".png");
   save_info.is_png = 1;
}

static define save_as_jpg_callback (w, save_info)
{
   set_save_filename_ext (save_info.window, ".jpg");
   save_info.is_png = 0;
}

static define file_selection_ok (save_info)
{
   variable file = gtk_file_selection_get_filename (save_info.window);
   variable ext = path_extname (file);
   
   if (ext == "")
     {
	if (save_info.is_png)
	  ext = ".png";
	else
	  ext = ".jpg";
	
	file += ext;
     }

   if (ext == ".png")
     save_info.is_png = 1;
   else if ((ext == ".jpg") or (ext == ".jpeg"))
     save_info.is_png = 0;

   save_info.filename = file;
   gtk_widget_destroy (save_info.window);
   save_info.window = NULL;
}

static define save_destroy_callback (w)
{
   gtk_main_quit ();
}

static define save_callback (window, rgb_widget)
{
   variable save_info = struct
     { 
	save_ok, is_png, filename, window
     };
   
   window = gtk_file_selection_new ("Save Image Dialog");
   save_info.window = window;

   gtk_file_selection_hide_fileop_buttons (window);
   gtk_window_set_position (window, GTK_WIN_POS_MOUSE);
   gtk_window_set_modal (window,1);

   () = g_signal_connect (window,"destroy", &save_destroy_callback);
   () = g_signal_connect_swapped (gtk_file_selection_get_ok_button(window),
				  "clicked", &file_selection_ok, save_info);

   variable frame, button, group, box;
   
   button = gtk_file_selection_get_cancel_button (window);
   () = g_signal_connect_swapped (button,"clicked",
				  &gtk_widget_destroy, window);

   frame = gtk_frame_new ("File Type");

   box = gtk_vbox_new (1, 0);
   gtk_container_add (frame, box);

   group = NULL;
   (button, group) = make_radio_button (box, NULL, "png",
   					&save_as_png_callback, save_info);

   (button, group) = make_radio_button (box, group, "jpeg",
   					&save_as_jpg_callback, save_info);
   gtk_toggle_button_set_active (button, TRUE);
   save_info.is_png = 0;

   box = gtk_file_selection_get_action_area (window);
   gtk_box_pack_start (box, frame, 1, 1, 10);

   gtk_widget_show_all (window);
   gtk_main ();

   variable file = save_info.filename;
   if (file == NULL)
     return;

   if (save_info.is_png)
     save_to_png (rgb_widget, file);
   else
     save_to_jpg (rgb_widget, file);
}

static define adjust_rgb_callback (rgb_limits, spectrum_info, rgb_widget)
{
   spectrum_info.rgb_limits = @rgb_limits;
   variable r, g, b;
   (r, g, b) = create_rgb (spectrum_info);
   
   rgb_display_widget_set_rgb (rgb_widget, r, g, b);
}

static define adjust_rgb (win, rgb_widget, spectrum_info)
{

   if (-1 == rgbselect (spectrum_info.engrid, spectrum_info.counts, 
			spectrum_info.rgb_limits, &adjust_rgb_callback, spectrum_info, rgb_widget))
     return;

   variable r, g, b;
   (r, g, b) = create_rgb (spectrum_info);
   
   rgb_display_widget_set_rgb (rgb_widget, r, g, b);
}

static define quit_callback (win)
{
   gtk_main_quit ();
}

static define wcs_callback_function (x, y, xgrid, ygrid)
{
   x *= 1.0;
   y *= 1.0;
   variable nx = length(xgrid);
   variable ny = length(ygrid);
   if ((nx <= 1) or (ny <= 1))
     return "";

   variable px = (xgrid[-1] - xgrid[0])/(nx-1.0);
   variable py = (ygrid[-1] - ygrid[0])/(ny-1.0);
   x = xgrid[0] + px*x + 0.5;
   y = ygrid[0] + py*y + 0.5;
   return sprintf ("(%g,%g)",x,y);
}

public define view_rgb (r, g, b, spectrum_info)
{
   variable rgb_widget = rgb_display_widget_new (r, g, b);

   variable save_button = gtk_button_new_with_label ("Save");
   () = g_signal_connect (save_button, "clicked", &save_callback, rgb_widget);
   gtk_box_pack_end (rgb_widget.button_hbox, save_button, 0, 0, 5);

   if (spectrum_info != NULL)
     {
	variable adjust_rgb_button = gtk_button_new_with_label ("Adjust RGB Levels");
	() = g_signal_connect (adjust_rgb_button, "clicked", &adjust_rgb, rgb_widget, spectrum_info);
	gtk_box_pack_start (rgb_widget.button_hbox, adjust_rgb_button, 0, 0, 5);
	gtk_box_reorder_child (rgb_widget.button_hbox, adjust_rgb_button, 0);
     }

   % Use convention where nth pixel runs from n-0.5 to n+0.5
   rgb_display_widget_set_wcs_callback (rgb_widget, &wcs_callback_function,
					spectrum_info.xgrid, spectrum_info.ygrid);
                                        %0.5+spectrum_info.xgrid[0],
                                        %0.5+spectrum_info.ygrid[0]);

   gtk_widget_show_all(rgb_widget.window);
   gtk_main ();   
}


public define create_and_display_rgb_image (spectrum_info)
{
   variable r, g, b;
   (r,g,b) = create_rgb (spectrum_info);
   view_rgb (r, g, b, spectrum_info);
}


%}}}

implements ("evt2img"); %{{{
static define read_events (file, xgrid, ygrid)
{
   ERROR_BLOCK 
     {
	_clear_error ();
	exit (1);
     }
   variable t = _traceback; _traceback = 0;
   % filtering is done this way to minimize memory use
   variable en, x, y, i, j;
   x = fits_read_col (file, "X");
   
   if (xgrid != NULL)
     {
	i = where ((x >= xgrid[0]) and (x < xgrid[-1]));
	x = x[i];
     }
   else i = [0:length(x)-1];

   y = fits_read_col (file, "Y")[i];   
   if (ygrid != NULL)
     {
	j = where ((y >= ygrid[0]) and (y < ygrid[-1]));
	y = y[j];
	x = x[j];
	i = i[j];
     }

   if (fits_binary_table_column_exists (file, "ENERGY"))
     {
	en = fits_read_col (file, "ENERGY")[i];
	en *= 0.001;
     }
   else
     {
	() = fprintf (stderr, "*** WARNING: ENERGY column does not exist.  Trying PI and PHA\n");
	if (fits_binary_table_column_exists (file, "PI"))
	  en = fits_read_col (file, "PI")[i];
	else
	  en = fits_read_col (file, "PHA")[i];
	
	() = fprintf (stderr, "*** WARNING: Arbitrarily assuming energies span from 1 to 10 keV\n");
	en = (10.0*en)/max(en);
     }

   _traceback = t;
   return x, y, en;
}

static define set_rgblimits (r0, r1, g0, g1, b0, b1)
{
   variable rgblimits = struct
     {
	r0, r1, g0, g1, b0, b1
     };
   rgblimits.r0 = r0;
   rgblimits.r1 = r1;
   rgblimits.g0 = g0;
   rgblimits.g1 = g1;
   rgblimits.b0 = b0;
   rgblimits.b1 = b1;
   return rgblimits;
}

static define main (file, xgrid, ygrid, rgb_limits)
{
   variable x, y, energy;
   (x,y,energy) = read_events (file, xgrid, ygrid);

   if ((xgrid == NULL) or (ygrid == NULL))
     {
	variable min_y = min(y);
	variable min_x = min(x);
	variable max_y = max(y);
	variable max_x = max(x);
	variable dy = max_y - min_y;
	variable dx = max_x - min_x;
	if (dy > dx) dx = dy;
	ygrid = min_y+[0:511]*(dx/512.0);
	xgrid = min_x+[0:511]*(dx/512.0);
     }
	
   variable engrid = [min(energy):max(energy):0.01];
   variable counts = hist1d (energy, engrid);

   if (rgb_limits == NULL)
     {
	variable s = cumsum (counts);
	variable r1 = engrid[hist_bsearch (s[-1]/3.0, s)];
	variable g1 = engrid[hist_bsearch ((2.0/3.0)*s[-1], s)];
	rgb_limits = set_rgblimits (NULL, r1, r1, g1, g1, NULL);
     }
   
   if (rgb_limits.r0 == NULL)
     rgb_limits.r0 = engrid[0];
   if (rgb_limits.b1 == NULL)
     rgb_limits.b1 = engrid[-1];
   
   vmessage ("Using: Red=%g-%g, Green=%g-%g, Blue=%g-%g",
	     rgb_limits.r0, rgb_limits.r1,
	     rgb_limits.g0, rgb_limits.g1,
	     rgb_limits.b0, rgb_limits.b1);

   variable spectrum_info = struct 
     {
	x, y, xgrid, ygrid,
	  energy, engrid, counts,
	  rgb_limits, wcs
     };
   spectrum_info.x = x;
   spectrum_info.y = y;
   spectrum_info.xgrid = xgrid;
   spectrum_info.ygrid = ygrid;
   spectrum_info.energy = energy;
   spectrum_info.engrid = engrid;
   spectrum_info.counts = counts;
   spectrum_info.rgb_limits = rgb_limits;

   create_and_display_rgb_image (spectrum_info);
}

%   XGrid = 4096-256 + [0:511];
%   YGrid = 4096-256 + [0:511];
%main ("./o120_g02346Es_bp_dstrk-s4_evt1.fits");

static define usage ()
{
   () = fprintf (stderr, "Usage: %s evt-file [ X/YGrid [energy-bands] ]\n", __argv[0]);
   () = fprintf (stderr, "The X/Y grids are specified using one of the forms\n");
   () = fprintf (stderr, "   xlo:xhi,ylo:yhi\n");
   () = fprintf (stderr, "   xlo:xhi:xinc,ylo:yhi:yinc\n");
   () = fprintf (stderr, "The default X/Y grid is 3840:4351,3840:4351\n");
   () = fprintf (stderr, "The RGB energy bounds are given either by one of the forms\n");
   () = fprintf (stderr, "   Rhi,Blo\n");
   () = fprintf (stderr, "   Rhi,Glo,Ghi,Blo\n");
   () = fprintf (stderr, "   Rlo,Rhi,Glo,Ghi,Blo,Bhi\n");
   () = fprintf (stderr, "in units of keV\n");


   exit (1);
}

static define parse_grid (grid)
{
   variable lo, hi, inc, n;
   
   inc = 1;
   n = sscanf (grid, "%g:%g:%g", &lo, &hi, &inc);
   if (n >= 2)
     {
	variable g = [lo:hi:inc];
	if (length (g))
	  return g;
     }
   () = fprintf (stderr, "**** Invalid XYGrid specification ****\n");
   exit (1);
}
   
static define parse_xygrid (xygrid)
{
   xygrid = strtok (xygrid, ",");
   if (length (xygrid) == 1)
     return parse_grid (xygrid[0]), NULL;

   return parse_grid (xygrid[0]), parse_grid (xygrid[1]);
}

static define parse_rgblimits (rgbspec)
{	
   variable p1 = -1, p2 = -1, p3 = -1, p4 = -1, p5 = -1, p6 = -1;

   switch (sscanf (rgbspec, "%g,%g,%g,%g,%g,%g", &p1,&p2,&p3,&p4,&p5,&p6))
     {
      case 6: % ok
     }
     {
      case 4:
	(p1,p2,p3,p4,p5,p6) = (NULL,p1,p2,p3,p4,NULL);
     }
     {
      case 2:
	(p1,p2,p3,p4,p5,p6) = (NULL,p1,p1,p2,p2,NULL);
     }
     {
	() = fprintf (stderr, "*** Invalid RGB grid specification ***\n");
	usage ();
     }
   return set_rgblimits (p1,p2,p3,p4,p5,p6);
}

static define parse_argv ()
{
   variable p1,p2,p3,p4,p5,p6;
   variable i;
   variable evtfile = NULL;
   variable xgrid = NULL, ygrid = NULL;
   variable rgblimits = NULL;

   i = 1;
   while (i < __argc)
     {
	variable arg = __argv[i];
	i++;

	if (evtfile == NULL)
	  {
	     evtfile = arg;
	     continue;
	  }
	if (xgrid == NULL)
	  {
	     (xgrid, ygrid) = parse_xygrid (arg);
	     continue;
	  }
	if (ygrid == NULL)
	  {
	     ygrid = parse_grid (arg);
	     continue;
	  }
	if (rgblimits == NULL)
	  {
	     rgblimits = parse_rgblimits (arg);
	     continue;
	  }
	usage ();
     }
   if (evtfile == NULL)
     usage ();

#iffalse
   if (xgrid == NULL)
     xgrid = 4096-256 + [0:511];
   
   if (ygrid == NULL)
     ygrid = 4096-256 + [0:511];
#endif

   return evtfile, xgrid, ygrid, rgblimits;
}

main (parse_argv ());
exit (0);


%}}}