aboutsummaryrefslogtreecommitdiff
path: root/noao/twodspec/multispec/exstrip.x
blob: a114b5a816b92c782f447464d339b56a114b02b5 (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
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
include	<imhdr.h>
include	<math/interp.h>
include "ms.h"

# EX_STRIP  -- Simple strip extraction of spectra.
# EX_STRIP1 -- Extract integrated spectra.
# EX_STRIP2 -- Extract two dimensional strip spectra.


# EX_STRIP -- Simple strip extraction of spectra.
#
# This procedure is called either by t_extract to extract spectra (either
# integrated or strip) or by t_newimage to extract a new image.
# Since there is no modeling only data spectra or image lines are extracted.
# It outputs the extracted spectra or image lines to the output image file.

procedure ex_strip (ms, im_in, im_out, spectra, lines, lower, upper,
    ex_spectra, ex_model, ex_integral)

pointer	ms				# MULTISPEC pointer
pointer	im_in				# Input image descriptor
pointer	im_out				# Output image descriptor
int	spectra[ARB]			# Spectra range list
int	lines[ARB]			# Line range list
real	lower				# Lower limit of strips
real	upper				# Upper limit of strips
bool	ex_spectra			# Extract spectra or image line
bool	ex_model			# Extract model or data
bool	ex_integral			# Extract integrated spectra or strip

int	line_in, line_out
pointer	data_in, data_out

int	get_next_number()
pointer	imgl2r(), impl2r()

begin
	if (ex_model)
	    call error (MS_ERROR, "Can't extract model")

	if (ex_spectra) {
	    # Extract spectra using ex_strip1 for integrated spectra and
	    # ex_strip2 for strip spectra.
	    if (ex_integral)
	        call ex_strip1 (ms, im_in, im_out, spectra, lines, lower,
		    upper)
	    else
	        call ex_strip2 (ms, im_in, im_out, spectra, lines, lower,
		    upper)
	} else {
	    # Create a new multi-spectra image by copying the selected
	    # input image lines to the output image.
	    line_in = 0
	    line_out = 0
	    while (get_next_number (lines, line_in) != EOF) {
		line_out = line_out + 1
		data_in = imgl2r (im_in, line_in)
		data_out = impl2r (im_out, line_out)
		call amovr (Memr[data_in], Memr[data_out], IM_LEN(im_out, 1))
	    }
	}
end

# EX_STRIP1 -- Extract integrated spectra.
#
# For each spectrum in the spectra range list and for each line in
# the line range list the pixels between lower and upper (relative
# to the spectrum center) are summed.
# The spectra positions are obtained from the MULTISPEC database.

procedure ex_strip1 (ms, im_in, im_out, spectra, lines, lower, upper)

pointer	ms				# MULTISPEC pointer
pointer	im_in				# Input image descriptor
pointer	im_out				# Output image descriptor
int	spectra[ARB]			# Spectra range list
int	lines[ARB]			# Line range list
real	lower				# Lower limit of strips
real	upper				# Upper limit of strips

int	line_in, line_out, spectrum_in, spectrum_out
real	x_center, x_start, x_end
pointer	buf_in, buf_out

real	sum_pixels(), cveval()
int	get_next_number()
pointer	imgl2r(), impl3r()

begin
	# Get fit functions for spectra positions. 
	call msgfits (ms, X0_FIT)

	# Loop through the input lines and write integrated spectra out.
	line_in = 0
	line_out = 0
	while (get_next_number (lines, line_in) != EOF) {
	    line_out = line_out + 1

	    # Get the input data line.
	    buf_in = imgl2r (im_in, line_in)

	    # Loop the the spectra, calculate the integrated luminosity and
	    # write it to the output image.
	    spectrum_in = 0
	    spectrum_out = 0
	    while (get_next_number (spectra, spectrum_in) != EOF) {
		spectrum_out = spectrum_out + 1

		buf_out = impl3r (im_out, line_out, spectrum_out)

		# Determine the spectrum limits from spectrum center position.
		x_center = cveval (CV(ms, X0_FIT, spectrum_in), real (line_in))
		x_start = max (1., x_center + lower)
		x_end = min (real (IM_LEN(im_in, 1)), x_center + upper)
		Memr[buf_out] =
		    sum_pixels (Memr[buf_in], x_start, x_end)
	    }
	}
end

# EX_STRIP2 -- Extract two dimensional strip spectra.
#
# Each line in the range list is fit by an image interpolator and then for
# each spectrum in spectra range list the interpolator values between lower
# and upper (relative to the spectrum center) are written to a three
# dimensional image.  There is one band for each spectrum.  The spectra
# positions are obtained from the MULTISPEC database.
# The procedure requests the interpolator type using CLIO.

procedure ex_strip2 (ms, im_in, im_out, spectra, lines, lower, upper)

pointer	ms				# MULTISPEC pointer
pointer	im_in				# Input image descriptor
pointer	im_out				# Output image descriptor
int	spectra[ARB]			# Spectra range list
int	lines[ARB]			# Line range list
real	lower				# Lower limit of strip
real	upper				# Upper limit of strip

int	interpolator			# Array interpolar type

int	i, len_in, len_out, line_in, line_out, spectrum_in, spectrum_out
real	x, x_start
pointer	buf_in, buf_out
pointer	sp, coeff

int	get_next_number(), clginterp()
real	asival(), cveval()
pointer	imgl2r(), impl3r()
errchk	salloc, imgl2r, impl3r
errchk	asiset, asifit, asival, clginterp

begin
	# Get the image interpolator type.
	interpolator = clginterp ("interpolator")

	len_in = IM_LEN (im_in, 1)
	len_out = nint (upper - lower + 1)

	# Set up the interpolator coefficient array.
	call smark (sp)
	call salloc (coeff, 2 * len_in + SZ_ASI, TY_REAL)
	call asiset (Memr[coeff], interpolator)

	# Get the spectra position functions from the database.
	call msgfits (ms, X0_FIT)

	# Loop through the input lines, do the image interpolation and write
	# the strip spectra to the output.
	line_in = 0
	line_out = 0
	while (get_next_number (lines, line_in) != EOF) {
	    line_out = line_out + 1

	    # Get the input data and fit an interpolation function.
	    buf_in = imgl2r (im_in, line_in)
	    call asifit (Memr[buf_in], len_in, Memr[coeff])

	    # Loop through the spectra writing the strip spectra.
	    spectrum_in = 0
	    spectrum_out = 0
	    while (get_next_number (spectra, spectrum_in) != EOF) {
		spectrum_out = spectrum_out + 1
		buf_out = impl3r (im_out, line_out, spectrum_out)

		# Determine the starting position for the strips and
		# evaluate the interpolation function at each point in
		# the strip.
		x_start = cveval (CV(ms, X0_FIT, spectrum_in), real (line_in)) +
		    lower
		do i = 1, len_out {
		    x = x_start + i - 1
		    if ((x < 1) || (x > len_in))
			Memr[buf_out + i - 1] = 0.
		    else
		        Memr[buf_out + i - 1] = asival (x, Memr[coeff])
		}
	    }
	}

	# Free interpolator memory.
	call sfree (sp)
end