aboutsummaryrefslogtreecommitdiff
path: root/pkg/plot/crtpict/crtulut.x
blob: 43609c55a610655fcdd2aa79e2798a4159c31ee4 (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.

include	<error.h>
include	<ctype.h>
include	"crtpict.h"

# CRT_ULUT -- Generates a look up table from data supplied by user.  The
# data is read from a two column text file of intensity, greyscale values.
# The input data are sorted, then mapped to the x range [0-4096].  A 
# piecewise linear look up table of 4096 values is then constructed from 
# the (x,y) pairs given.  A pointer to the look up table, as well as the z1 
# and z2 intensity endpoints, is returned.

procedure crt_ulut (fname, z1, z2, lut)

char	fname[SZ_FNAME]		# Name of file with intensity, greyscale values
real	z1			# Intensity mapped to minimum gs value
real	z2			# Intensity mapped to maximum gs value
pointer	lut			# Look up table - pointer is returned

pointer	sp, x, y
int	nvalues, i, j, x1, x2, y1
real	delta_gs, delta_xv, slope
errchk	crt_rlut, crt_sort, malloc

begin
	call smark (sp)
	call salloc (x, SZ_BUF, TY_REAL)
	call salloc (y, SZ_BUF, TY_REAL)

	# Read intensities and greyscales from the user's input file.  The
	# intensity range is then mapped into a standard range and the 
	# values sorted.

	call crt_rlut (fname, Memr[x], Memr[y], nvalues)
	call alimr (Memr[x], nvalues, z1, z2)
	call amapr (Memr[x], Memr[x], nvalues, z1, z2, STARTPT, ENDPT)
	call crt_sort (Memr[x], Memr[y], nvalues)

	# Fill lut in straight line segments - piecewise linear
	call malloc (lut, SZ_BUF, TY_SHORT)
	do i = 1, nvalues-1 {
	    delta_gs = Memr[y+i] - Memr[y+i-1]
	    delta_xv = Memr[x+i] - Memr[x+i-1]
	    slope = delta_gs / delta_xv
	    x1 = int (Memr[x+i-1]) 
	    x2 = int (Memr[x+i])
	    y1 = int (Memr[y+i-1])
	    do j = x1, x2-1
		Mems[lut+j-1] = y1 + slope * (j-x1)
	}

	call sfree (sp)
end


# CRT_RLUT -- Read text file of x, y, values.

procedure crt_rlut (utab, x, y, nvalues)

char	utab[SZ_FNAME]		# Name of list file
real	x[SZ_BUF]		# Array of x values, filled on return
real	y[SZ_BUF]		# Array of y values, filled on return
int	nvalues			# Number of values in x, y vectors - returned

int	n, fd
pointer	sp, lbuf, ip
real	xval, yval
int	getline(), open()
errchk	open, sscan, getline, malloc

begin
	call smark (sp)
	call salloc (lbuf, SZ_LINE, TY_CHAR)

	iferr (fd = open (utab, READ_ONLY, TEXT_FILE))
	    call error (0, "Error opening user table")

	n = 0

	while (getline (fd, Memc[lbuf]) != EOF) {
	    # Skip comment lines and blank lines.
	    if (Memc[lbuf] == '#')
		next
	    for (ip=lbuf;  IS_WHITE(Memc[ip]);  ip=ip+1)
		;
	    if (Memc[ip] == '\n' || Memc[ip] == EOS)
		next

	    # Decode the points to be plotted.
	    call sscan (Memc[ip])
		call gargr (xval)
		call gargr (yval)

	    n = n + 1
	    if (n > SZ_BUF) 
	        call error (0,
		    "Intensity transformation table cannot exceed 4096 values")

	    x[n] = xval
	    y[n] = yval
	}

	nvalues = n
	call close (fd)
	call sfree (sp)
end


# CRT_SORT -- Bubble sort of paired arrays.  

procedure crt_sort (xvals, yvals, nvals)

real	xvals[nvals]		# Array of x values
real	yvals[nvals]		# Array of y values
int	nvals			# Number of values in each array

int	i, j
real	temp
define	swap	{temp=$1;$1=$2;$2=temp}

begin
	for (i = nvals; i > 1; i = i - 1)
	    for (j = 1; j < i; j = j + 1) 
		if (xvals[j] > xvals[j+1]) {
		    # Out of order; exchange y values
		    swap (xvals[j], xvals[j+1])
		    swap (yvals[j], yvals[j+1])
		}
end