diff options
Diffstat (limited to 'noao/twodspec/multispec/exstrip.x')
| -rw-r--r-- | noao/twodspec/multispec/exstrip.x | 203 |
1 files changed, 203 insertions, 0 deletions
diff --git a/noao/twodspec/multispec/exstrip.x b/noao/twodspec/multispec/exstrip.x new file mode 100644 index 00000000..a114b5a8 --- /dev/null +++ b/noao/twodspec/multispec/exstrip.x @@ -0,0 +1,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 |