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
|
include <imhdr.h>
include <fset.h>
include "ms.h"
# T_FIND_PEAKS -- Find the spectra peaks in a MULTISPEC image and record
# their positions in the database.
#
# An average of naverage lines from the MULTISPEC image is searched
# for peaks satisfying constraints on the minimum and maximum number,
# columns, peak values, and separation between peaks. The positions
# of the peaks satisfying these constraints is entered in the database.
# It is an error if fewer than the minimum number of peaks is found
# or if the number of peaks differs from a previously determined number.
# The peak finding is done by the function FIND_PEAKS which is numerical
# and may be used outside the MULTISPEC package.
procedure t_find_peaks ()
# CL parameters:
char image[SZ_FNAME] # Image to be searched
int lines[3, MAX_RANGES] # Image lines in which to find spectra
int min_npeaks # Minimum number of spectra to be found
int max_npeaks # Maximum number of spectra to be accepted
int separation # Minimum pixel separation between spectra
int edge # Minimum distance to edge of image
real threshold # Minimum peak value
real contrast # Max contrast between strongest and weakest
int columns[3, MAX_RANGES] # Spectra positions limited to these columns
int naverage # Number of image lines to average
bool debug # Print debugging information
char comment[SZ_LINE]
int i, j, k, line, sample, nsamples, npoints, nspectra
pointer ms, im
pointer sp, data, x, samples
int find_peaks(), get_sample_lines()
int clgeti(), clgranges()
real clgetr()
bool clgetb(), is_in_range()
pointer msmap(), immap()
begin
# Get task parameters and access files.
call clgstr ("image", image, SZ_FNAME)
ms = msmap (image, READ_WRITE, 0)
im = immap (image, READ_ONLY, 0)
i = clgranges ("lines", 1, IM_LEN(im, 2), lines, MAX_RANGES)
min_npeaks = clgeti ("min_npeaks")
max_npeaks = clgeti ("max_npeaks")
separation = clgeti ("separation")
edge = clgeti ("edge")
threshold = clgetr ("threshold")
contrast = clgetr ("contrast")
i = clgranges ("columns", 1, IM_LEN(im, 1), columns, MAX_RANGES)
naverage = clgeti ("naverage")
debug = clgetb ("debug")
call fseti (STDOUT, F_FLUSHNL, YES)
# Allocate working memory.
npoints = IM_LEN(im, 1)
call smark (sp)
call salloc (samples, MS_NSAMPLES(ms), TY_INT)
call salloc (data, npoints, TY_REAL)
call salloc (x, npoints, TY_REAL)
# Get the sample lines.
nsamples = get_sample_lines (ms, lines, Memi[samples])
# Loop through each sample line.
do i = 1, nsamples {
sample = Memi[samples + i - 1]
line = LINE(ms, sample)
# Get the image data with averaging.
call msgimage (im, line, naverage, Memr[data])
# Mark columns which are to be ignored with INDEFR.
do j = 1, npoints
if (!is_in_range (columns, j))
Memr[data + j - 1] = INDEFR
# Find the peaks.
nspectra = find_peaks (Memr[data], Memr[x], npoints,
contrast, separation, edge, max_npeaks, threshold, debug)
if (debug) {
call printf (" Number of spectra found in line %d = %d.\n")
call pargi (line)
call pargi (nspectra)
}
if (nspectra < min_npeaks)
call error (MS_ERROR, "Too few spectra found")
# Enter the spectra found in the database. If the number of
# spectra has not been previously set in the database then
# enter the number of spectra and make entries in the
# database. Otherwise check that the number of spectra found
# agrees with that already in the database.
if (MS_NSPECTRA(ms) == 0) {
if (nspectra == 0)
next
MS_NSPECTRA(ms) = nspectra
call dbenter (MS_DB(ms), NAME(ms, I0), nspectra * SZ_REAL,
MS_NSAMPLES(ms))
call dbenter (MS_DB(ms), NAME(ms, X0), nspectra * SZ_REAL,
MS_NSAMPLES(ms))
} else if (MS_NSPECTRA(ms) != nspectra)
call error (MS_ERROR, "Attempt to change the number of spectra")
call msgparam (ms, X0, sample)
call amovr (Memr[x], PARAMETER(ms, X0, 1), nspectra)
call mspparam (ms, X0, sample)
# The peak scale is taken and the pixel value at the peak.
call msgparam (ms, I0, sample)
do j = 1, nspectra {
k = PARAMETER(ms, X0, j)
PARAMETER(ms, I0, j) = Memr[data + k - 1]
}
call mspparam (ms, I0, sample)
# Enter a comment in the database.
call sprintf (comment, SZ_LINE,
"Spectra located in sample line %d.")
call pargi (sample)
call history (ms, comment)
}
# Update the database and close the database and image.
call msphdr (ms)
call msunmap (ms)
call imunmap (im)
call sfree (sp)
end
|
