Initial commit

1 files changed, 209 insertions, 0 deletions

diff --git a/pkg/plot/crtpict/plothgms.x b/pkg/plot/crtpict/plothgms.x
new file mode 100644
index 00000000..803c709f
--- /dev/null
+++ b/pkg/plot/crtpict/plothgms.x
@@ -0,0 +1,209 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<gset.h>
+include	<imhdr.h>
+include	<mach.h>
+include "wdes.h"
+include	"crtpict.h"
+
+# CRT_PLOT_HISTOGRAMS -- Calculate and plot three histograms describing the
+# intensity to greyscale mapping.
+
+procedure crt_plot_histograms (gp, cl, im, wdes, xs, xe, ys, ye)
+
+pointer	gp
+pointer	cl		# Pointer to cl structure.
+pointer im
+pointer	wdes
+real	xs, xe, ys, ye, z1, z2
+
+pointer	w1, inten_hgram, greys_hgram, sp, text, syv, greys, hgram
+pointer real_greys, xval, yval
+int	nsig_bits, i, zrange, mask
+real	plot_ys, plot_ye, plot_width, plot_spacing, x, y, delta_inten, inten
+real	plot1_xs, plot1_xe, plot2_xs, plot2_xe, plot3_xs, plot3_xe 
+real	dz1, dz2, gio_char_x
+real	major_length, minor_length, ux1, ux2, uy1, uy2, y_pos, label_y
+real	wx1, wx2, wy1, wy2
+
+bool	ggetb()
+real 	ggetr()
+int	ggeti(), and()
+errchk	ggetb, ggeti, ggetr, crt_calc_hgrams, ggwind, gswind, gsview, gploto
+errchk	gsetr, gseti, glabax, amapr, gpline, gvline, achtir, gtext
+
+begin
+	call smark (sp)
+	call salloc (inten_hgram, NBINS, TY_INT)
+	call salloc (greys_hgram, NBINS, TY_INT)
+	call salloc (text, SZ_LINE, TY_CHAR)
+	call salloc (syv, NBINS, TY_SHORT)
+	call salloc (greys, NBINS, TY_INT)
+	call salloc (hgram, NBINS, TY_REAL)
+	call salloc (real_greys, NBINS, TY_REAL)
+	call salloc (xval, NBINS, TY_REAL)
+	call salloc (yval, NBINS, TY_REAL)
+
+	# First, get pointer to WCS 1 and some device parameters.
+	w1 = W_WC(wdes, 1)
+	z1 = W_ZS(w1)
+	z2 = W_ZE(w1)
+
+	# If z1 and z2 not in graphcap, set to some reasonable numbers for
+	# plots to be generated.
+
+	if (ggetb (gp, "z1") && ggetb (gp, "z2")) {
+	    dz1 = ggetr (gp, "z1")
+	    dz2 = ggetr (gp, "z2")
+	} else {
+	    dz1 = 0.
+	    dz2 = 255.
+	}
+
+	# To allow room for annotation, the y limits of each plot are
+	# drawn in by 5%.  The y limits are the same for each plot.
+
+	plot_ys = ys + (0.05 * (ye - ys))
+	plot_ye = ye - (0.05 * (ye - ys))
+	label_y = plot_ys - (plot_ye - plot_ys) * 0.20
+
+	# Now calculate the x limits.  Each plot occupys a fourth of the
+	# available space in x.  The distance between plot centers is a
+	# third of the available space.
+
+	plot_width = (xe - xs) / 4.0
+	plot_spacing = (xe - xs) / 3.0
+
+	plot1_xs = xs + (plot_spacing / 2.0) - (plot_width / 2.0)
+	plot1_xe = plot1_xs + plot_width
+	plot2_xs = plot1_xs + plot_spacing
+	plot2_xe = plot2_xs + plot_width
+	plot3_xs = plot2_xs + plot_spacing
+	plot3_xe = plot3_xs + plot_width
+
+	# Calculate the histograms for both the untransformed (intensity) and
+	# transformed (greyscale) image in a single procedure. A separate
+	# path is taken for linear or user transformations:
+
+	if (W_ZT(w1) == W_USER)
+	    call crt_user_hgram (im, gp, z1, z2, Mems[LUT(cl)], 
+		Memi[inten_hgram], Memi[greys_hgram])
+	else
+	    call crt_linear_hgram (im, gp, z1, z2, W_ZT(w1), Memi[inten_hgram], 
+	        Memi[greys_hgram])
+
+	# Each histogram plot is a separate mapping in WCS 3
+	call gseti (gp, G_WCS, 3)
+
+	# The first histogram shows the number of pixels at a given
+	# intensity versus intensity for the original image.
+
+	call gsview (gp, plot1_xs, plot1_xe, plot_ys, plot_ye)
+	gio_char_x = ((plot1_xe - plot1_xs) / 50.) / ggetr (gp, "cw")
+	major_length = gio_char_x * (ggetr (gp, "cw"))
+	minor_length = gio_char_x * (0.5 * (ggetr (gp, "cw")))
+
+	call gseti (gp, G_YTRAN, GW_LOG)
+	call gseti (gp, G_ROUND, YES)
+	call gsetr (gp, G_MAJORLENGTH, major_length)
+	call gsetr (gp, G_MINORLENGTH, minor_length)
+	call gseti (gp, G_LABELAXIS, NO)
+	call gsetr (gp, G_TICKLABELSIZE, 0.25)
+	call achtir (Memi[inten_hgram], Memr[hgram], NBINS)
+	call gploto (gp, Memr[hgram], NBINS, IM_MIN(im), IM_MAX(im), "")
+
+	# Now to label the plot:
+	call ggwind (gp, ux1, ux2, uy1, uy2)
+	y_pos = uy1 - ((uy2 - uy1) * 0.20)  #y_pos below yaxis by 20% of height
+	call gseti (gp, G_YTRAN, GW_LINEAR)
+	call gtext (gp, (ux1+ux2)/2.0, y_pos, "LOG10(N(DN)) VS DN", 
+	    "v=t;h=c;s=.25")
+
+	# The third plot shows the number of pixels at a given greyscale
+	# level versus greyscale level for the range of intensities
+	# transformed. 
+
+	call gsview (gp, plot3_xs, plot3_xe, plot_ys, plot_ye)
+	call achtir (Memi[greys_hgram], Memr[hgram], NBINS)
+
+	if (z2 > z1)
+	    call gploto (gp, Memr[hgram], NBINS, dz1, dz2, "")
+	else
+	    call gploto (gp, Memr[hgram], NBINS, dz2, dz1, "")
+
+	# Now to label the plot:
+	call ggwind (gp, ux1, ux2, uy1, uy2)
+	call gseti (gp, G_YTRAN, GW_LINEAR)
+	y_pos = uy1 - ((uy2 - uy1) * 0.20) # y_pos below yaxis by 20% of height
+	call gtext (gp, (ux1+ux2)/2.0, y_pos, "TRANSFORMED HISTOGRAM", 
+	    "v=t;h=c;s=.25")
+
+	# The second plot shows how the dynamic range of the transformed
+	# image maps to the dynamic range of the output device.
+
+	call gsview (gp, plot2_xs, plot2_xe, plot_ys, plot_ye)
+	call gswind (gp, IM_MIN(im), IM_MAX(im), real (dz1), real (dz2))
+	call gseti (gp, G_YTRAN, GW_LINEAR)
+	call glabax (gp, "", "", "")
+
+	if (W_ZT(w1) != W_UNITARY) {
+	    do i = 1, NBINS 
+		Memr[xval+i-1] = IM_MIN(im) + (i-1) * (IM_MAX(im) - 
+		    IM_MIN(im))/ (NBINS-1)
+
+	    if (W_ZT(w1) == W_USER) {
+	        call sprintf (Memc[text], SZ_LINE, 
+		    "USER DEFINED FUNCTION: FROM %g TO %g")
+		call pargr (z1)
+		call pargr (z2)
+		call amapr (Memr[xval], Memr[yval], NBINS, z1, z2, STARTPT, 
+		    ENDPT)
+		call achtrs (Memr[yval], Mems[syv], NBINS)
+		call aluts (Mems[syv], Mems[syv], NBINS, Mems[LUT[cl]])
+		call achtsr (Mems[syv], Memr[yval], NBINS)
+	    } else {
+	        call sprintf (Memc[text], SZ_LINE, 
+		    "TRANSFER FUNCTION: LINEAR FROM %g TO %g")
+		call pargr (z1)
+		call pargr (z2)
+	        call amapr (Memr[xval], Memr[yval], NBINS, z1, z2, real (dz1), 
+		    real (dz2))
+	    }
+
+	    call gpline (gp, Memr[xval], Memr[yval], NBINS)
+	    call ggwind (gp, wx1, wx2, wy1, wy2)
+	    x = (wx2 + wx1) / 2.0
+	    y = wy1 - (wy2 - wy1) * 0.20
+	    call gtext (gp, x, y, Memc[text], "h=c;v=t;s=0.25")
+	    call printf ("%s\n")
+		call pargstr (Memc[text])
+
+	} else {
+	    # Calculate number of bits depth in output device
+	    zrange = ggeti (gp, "zr")
+	    for (nsig_bits = 0; ; nsig_bits = nsig_bits + 1) {
+	        zrange = zrange / 2
+	        if (zrange == 0)
+	            break
+	    }
+
+	    # Truncate intensities to dynamic range of output device.
+	    delta_inten = (IM_MAX(im) - IM_MIN(im)) / (NBINS - 1)
+	    mask = 2**(nsig_bits) - 1
+	    do i = 1, NBINS {
+		inten =  IM_MIN(im) + ((i-1) * delta_inten)
+		Memi[greys+i-1] = and (int (inten), mask)
+	    }
+
+	    call achtir (Memi[greys], Memr[real_greys], NBINS)
+	    call gvline (gp, Memr[real_greys], NBINS, IM_MIN(im), IM_MAX(im))
+	    call ggwind (gp, wx1, wx2, wy1, wy2)
+	    x = (wx2 + wx1) / 2.0
+	    y = wy1 - (wy2 - wy1) * 0.20
+	    call gtext (gp, x, y, "TRANSFER FUNCTION: UNITARY","h=c;v=t;s=0.25")
+	    call printf ("Unitary Transfer Function; Lowest %d bits output.\n")
+		call pargi (nsig_bits)
+	}
+
+	call sfree (sp)
+end