CR2RE Pipeline Reference Manual 1.6.11
irplib_ppm.c
1/*
2 * This file is part of the irplib package
3 * Copyright (C) 2002,2003,2014 European Southern Observatory
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02111-1307 USA
18 */
19
20#ifdef HAVE_CONFIG_H
21#include <config.h>
22#endif
23
24/*-----------------------------------------------------------------------------
25 Includes
26 -----------------------------------------------------------------------------*/
27
28#include "irplib_ppm.h"
29#include "irplib_wlxcorr.h"
30#include "irplib_spectrum.h"
31
32#include <math.h>
33#include <cpl.h>
34
35/*-----------------------------------------------------------------------------
36 Private functions
37 -----------------------------------------------------------------------------*/
38#ifdef IRPLIB_PPM_USE_METHOD2
39static cpl_vector * irplib_ppm_convolve_line(const cpl_vector *, double,double);
40static cpl_vector * irplib_ppm_detect_lines(const cpl_vector *, double) ;
41#endif
42
43/*----------------------------------------------------------------------------*/
47/*----------------------------------------------------------------------------*/
48
50
51/*----------------------------------------------------------------------------*/
68/*----------------------------------------------------------------------------*/
69cpl_polynomial * irplib_ppm_engine(
70 const cpl_vector * spectrum,
71 const cpl_bivector * lines_catalog,
72 const cpl_polynomial * poly_init,
73 double slitw,
74 double fwhm,
75 double thresh,
76 int degree,
77 int doplot,
78 cpl_table ** tab_infos)
79{
80#ifdef IRPLIB_PPM_USE_METHOD2
81 cpl_vector * spec_conv ;
82#endif
83 int spec_sz ;
84 cpl_vector * det_lines ;
85 cpl_vector * cat_lines ;
86 double * pcat_lines ;
87 double wmin, wmax ;
88 double disp_min, disp_max, disp ;
89 int nlines_cat, nlines ;
90 const double * plines_catalog_x ;
91 const double * plines_catalog_y ;
92 cpl_bivector * matched ;
93 cpl_matrix * matchedx;
94 int match_sz;
95 cpl_polynomial * fitted ;
96 cpl_table * spc_table ;
97 const cpl_vector* vectors_plot[3];
98 cpl_vector * plot_y ;
99 int start_ind, stop_ind ;
100 double fill_val ;
101 cpl_size deg_loc ;
102 int i ;
103 cpl_error_code error;
104
105 /* Check entries */
106 if (spectrum == NULL) return NULL ;
107 if (lines_catalog == NULL) return NULL ;
108 if (poly_init == NULL) return NULL ;
109
110 /* Initialise */
111 spec_sz = cpl_vector_get_size(spectrum) ;
112 deg_loc = (cpl_size)degree ;
113
114#ifdef IRPLIB_PPM_USE_METHOD2
115 /* METHOD 2 */
116 /* Correlate the spectrum with the line profile */
117 if ((spec_conv = irplib_ppm_convolve_line(spectrum, slitw, fwhm)) == NULL) {
118 cpl_msg_error(cpl_func, "Cannot convolve the signal") ;
119 return NULL ;
120 }
121
122 /* Apply the lines detection */
123 if ((det_lines = irplib_ppm_detect_lines(spec_conv, 0.9)) == NULL) {
124 cpl_msg_error(cpl_func, "Cannot detect lines") ;
125 cpl_vector_delete(spec_conv) ;
126 return NULL ;
127 }
128 cpl_vector_delete(spec_conv) ;
129#else
130 /* METHOD 1 */
131 if ((det_lines = irplib_spectrum_detect_peaks(spectrum, fwhm,
132 thresh, 0, NULL, NULL)) == NULL) {
133 cpl_msg_error(cpl_func, "Cannot convolve the signal") ;
134 return NULL ;
135 }
136#endif
137 cpl_msg_info(cpl_func, "Detected %"CPL_SIZE_FORMAT" lines",
138 cpl_vector_get_size(det_lines));
139
140 /* Get the catalog lines */
141 wmin = cpl_polynomial_eval_1d(poly_init, 1.0, NULL) ;
142 wmax = cpl_polynomial_eval_1d(poly_init, spec_sz, NULL) ;
143 plines_catalog_x = cpl_bivector_get_x_data_const(lines_catalog) ;
144 plines_catalog_y = cpl_bivector_get_y_data_const(lines_catalog) ;
145 nlines = cpl_bivector_get_size(lines_catalog) ;
146 nlines_cat = 0 ;
147 start_ind = stop_ind = -1 ;
148 for (i=0 ; i<nlines ; i++) {
149 if (plines_catalog_x[i] > wmin && plines_catalog_x[i] < wmax &&
150 plines_catalog_y[i] > 0.0) {
151 nlines_cat++ ;
152 if (start_ind<0) start_ind = i ;
153 stop_ind = i ;
154 }
155 }
156 if (nlines_cat == 0) {
157 cpl_msg_error(cpl_func, "No lines in catalog") ;
158 cpl_vector_delete(det_lines) ;
159 return NULL ;
160 }
161 cat_lines = cpl_vector_new(nlines_cat) ;
162 pcat_lines = cpl_vector_get_data(cat_lines) ;
163 nlines_cat = 0 ;
164 for (i=0 ; i<nlines ; i++) {
165 if (plines_catalog_x[i] > wmin && plines_catalog_x[i] < wmax &&
166 plines_catalog_y[i] > 0.0) {
167 pcat_lines[nlines_cat] = plines_catalog_x[i] ;
168 nlines_cat++ ;
169 }
170 }
171
172 /* Plot inputs */
173 if (doplot) {
174 double * pdet_lines ;
175
176 /* Catalog */
177 irplib_wlxcorr_catalog_plot(lines_catalog, wmin, wmax) ;
178
179 /* Spectrum with detected lines */
180 fill_val = cpl_vector_get_max(spectrum) ;
181 plot_y = cpl_vector_new(spec_sz);
182 cpl_vector_fill(plot_y, 0.0) ;
183 pdet_lines = cpl_vector_get_data(det_lines) ;
184 for (i=0 ; i<cpl_vector_get_size(det_lines) ; i++) {
185 cpl_vector_set(plot_y, (int)pdet_lines[i], fill_val) ;
186 }
187 vectors_plot[0] = NULL ;
188 vectors_plot[1] = spectrum ;
189 vectors_plot[2] = plot_y ;
190
191 cpl_plot_vectors("set grid;set xlabel 'Position (Pixel)';set ylabel "
192 "'Intensity (ADU/sec)';",
193 "t 'Spectrum with detected lines' w lines", "",
194 vectors_plot, 3);
195 cpl_vector_delete(plot_y) ;
196 }
197
198 /* Apply the point pattern matching */
199 disp = (wmax-wmin) / spec_sz ;
200 disp_min = disp - (disp/10) ;
201 disp_max = disp + (disp/10) ;
202 matched = cpl_ppm_match_positions(det_lines, cat_lines, disp_min,
203 disp_max, 0.05, NULL, NULL);
204 cpl_vector_delete(det_lines) ;
205 cpl_vector_delete(cat_lines) ;
206
207 if (matched == NULL) {
208 cpl_msg_error(cpl_func, "Cannot apply the point pattern matching") ;
209 return NULL ;
210 }
211
212 match_sz = cpl_bivector_get_size(matched);
213
214 cpl_msg_info(cpl_func, "Matched %d lines", match_sz) ;
215
216 if (match_sz <= deg_loc) {
217 cpl_msg_error(cpl_func, "Not enough match for the fit") ;
218 cpl_bivector_delete(matched) ;
219 return NULL ;
220 }
221
222 /* Plot if requested */
223 if (doplot) {
224 const double * pmatched ;
225 cpl_bivector * biplot ;
226 cpl_vector * plot_cat_x ;
227 cpl_vector * plot_cat_y ;
228 /* Spectrum with matched lines */
229 fill_val = cpl_vector_get_max(spectrum) ;
230 plot_y = cpl_vector_new(spec_sz);
231 cpl_vector_fill(plot_y, 0.0) ;
232 pmatched = cpl_bivector_get_x_data_const(matched) ;
233 for (i=0 ; i < match_sz; i++) {
234 cpl_vector_set(plot_y, (int)pmatched[i], fill_val) ;
235 }
236 vectors_plot[0] = NULL ;
237 vectors_plot[1] = spectrum ;
238 vectors_plot[2] = plot_y ;
239
240 cpl_plot_vectors("set grid;set xlabel 'Position (Pixel)';set ylabel "
241 "'Intensity (ADU/sec)';",
242 "t 'Spectrum with matched lines' w lines", "",
243 vectors_plot, 3);
244 cpl_vector_delete(plot_y) ;
245
246 /* Catalog with matched lines */
247 plot_cat_x=cpl_vector_extract(cpl_bivector_get_x_const(lines_catalog),
248 start_ind, stop_ind, 1) ;
249 plot_cat_y=cpl_vector_extract(cpl_bivector_get_y_const(lines_catalog),
250 start_ind, stop_ind, 1) ;
251 biplot = cpl_bivector_wrap_vectors(plot_cat_x, plot_cat_y) ;
252 cpl_plot_bivector("set grid;set xlabel 'Wavelength';set ylabel "
253 "'Emission';", "t 'Catalog' w impulses", "",
254 biplot);
255 cpl_bivector_unwrap_vectors(biplot) ;
256
257 plot_y = cpl_vector_duplicate(plot_cat_y) ;
258 cpl_vector_fill(plot_y, 0.0) ;
259 pmatched = cpl_bivector_get_y_data_const(matched) ;
260 fill_val=cpl_vector_get_mean(plot_cat_y) ;
261 for (i=0 ; i < match_sz; i++) {
262 int wl_ind = 0 ;
263 while (pmatched[i] > cpl_vector_get(plot_cat_x, wl_ind)
264 && wl_ind < spec_sz) wl_ind++ ;
265 if (wl_ind < spec_sz) cpl_vector_set(plot_y, wl_ind, fill_val) ;
266 }
267 biplot = cpl_bivector_wrap_vectors(plot_cat_x, plot_y) ;
268 cpl_plot_bivector("set grid;set xlabel 'Wavelength';set ylabel "
269 "'Emission';", "t 'Catalog (matched lines)' w "
270 "impulses", "", biplot) ;
271 cpl_bivector_unwrap_vectors(biplot) ;
272 cpl_vector_delete(plot_cat_x) ;
273 cpl_vector_delete(plot_cat_y) ;
274 cpl_vector_delete(plot_y) ;
275 }
276
277 /* Apply the fit */
278 matchedx = cpl_matrix_wrap(1, match_sz, cpl_bivector_get_x_data(matched));
279 fitted = cpl_polynomial_new(1);
280 error = cpl_polynomial_fit(fitted, matchedx, NULL,
281 cpl_bivector_get_y_const(matched), NULL,
282 CPL_FALSE, NULL, &deg_loc);
283 cpl_bivector_delete(matched);
284 (void)cpl_matrix_unwrap(matchedx);
285 if (error) {
286 cpl_msg_error(cpl_func, "Cannot fit the polynomial") ;
287 cpl_polynomial_delete(fitted);
288 return NULL ;
289 }
290
291 /* Create the infos table */
292 if ((spc_table = irplib_wlxcorr_gen_spc_table(spectrum,
293 lines_catalog, slitw, fwhm, poly_init, fitted)) == NULL) {
294 cpl_msg_error(cpl_func, "Cannot generate the infos table") ;
295 cpl_polynomial_delete(fitted) ;
296 return NULL ;
297 }
298 if (tab_infos != NULL) *tab_infos = spc_table ;
299 else cpl_table_delete(spc_table) ;
300 return fitted ;
301}
302
304
305#ifdef IRPLIB_PPM_USE_METHOD2
306/*----------------------------------------------------------------------------*/
317/*----------------------------------------------------------------------------*/
318static cpl_vector * irplib_ppm_convolve_line(
319 const cpl_vector * spectrum,
320 double slitw,
321 double fwhm)
322{
323 cpl_vector * conv_kernel ;
324 cpl_vector * line_profile ;
325 cpl_vector * xcorrs ;
326 cpl_vector * xc_single ;
327 int hs, line_sz, sp_sz ;
328 int i ;
329
330 /* Test entries */
331 if (spectrum == NULL) return NULL ;
332
333 /* Create the convolution kernel */
334 if ((conv_kernel = irplib_wlxcorr_convolve_create_kernel(slitw,
335 fwhm)) == NULL) {
336 cpl_msg_error(cpl_func, "Cannot create kernel") ;
337 return NULL ;
338 }
339 hs = cpl_vector_get_size(conv_kernel) ;
340 line_sz = 2 * hs + 1 ;
341
342 /* Create the line profile */
343 line_profile = cpl_vector_new(line_sz) ;
344 cpl_vector_fill(line_profile, 0.0) ;
345 cpl_vector_set(line_profile, hs, 1.0) ;
346 if (irplib_wlxcorr_convolve(line_profile, conv_kernel) != 0) {
347 cpl_msg_error(cpl_func, "Cannot create line profile") ;
348 cpl_vector_delete(line_profile) ;
349 cpl_vector_delete(conv_kernel) ;
350 return NULL ;
351 }
352 cpl_vector_delete(conv_kernel) ;
353
354 /* Create the correlations values vector */
355 sp_sz = cpl_vector_get_size(spectrum) ;
356 xcorrs = cpl_vector_new(sp_sz) ;
357 cpl_vector_fill(xcorrs, 0.0) ;
358 xc_single = cpl_vector_new(1) ;
359
360 /* Loop on the pixels of the spectrum */
361 for (i=hs ; i<sp_sz-hs ; i++) {
362 cpl_vector * spec_ext ;
363 /* Extract the current spectrum part */
364 if ((spec_ext = cpl_vector_extract(spectrum, i-hs, i+hs, 1)) == NULL) {
365 cpl_msg_error(cpl_func, "Cannot extract spectrum") ;
366 cpl_vector_delete(xc_single) ;
367 cpl_vector_delete(line_profile) ;
368 return NULL ;
369 }
370 if (cpl_vector_correlate(xc_single, spec_ext, line_profile) < 0) {
371 cpl_msg_error(cpl_func, "Cannot correlate") ;
372 cpl_vector_delete(xc_single) ;
373 cpl_vector_delete(line_profile) ;
374 cpl_vector_delete(spec_ext) ;
375 return NULL ;
376 }
377 cpl_vector_set(xcorrs, i, cpl_vector_get(xc_single, 0)) ;
378 cpl_vector_delete(spec_ext) ;
379 }
380 cpl_vector_delete(xc_single) ;
381 cpl_vector_delete(line_profile) ;
382
383 return xcorrs ;
384}
385
386/*----------------------------------------------------------------------------*/
395/*----------------------------------------------------------------------------*/
396static cpl_vector * irplib_ppm_detect_lines(
397 const cpl_vector * spec,
398 double threshold)
399{
400 cpl_vector * spec_loc ;
401 double * pspec_loc ;
402 cpl_vector * lines ;
403 double * plines ;
404 int spec_loc_sz, nlines ;
405 double max ;
406 int i ;
407
408 /* Test inputs */
409 if (spec == NULL) return NULL ;
410
411 /* Local spectrum */
412 spec_loc = cpl_vector_duplicate(spec) ;
413 pspec_loc = cpl_vector_get_data(spec_loc) ;
414 spec_loc_sz = cpl_vector_get_size(spec_loc) ;
415
416 /* Threshold the local spectrum */
417 for (i=0 ; i<spec_loc_sz ; i++)
418 if (pspec_loc[i] < threshold) pspec_loc[i] = 0.0 ;
419
420 /* Allocate lines container */
421 lines = cpl_vector_new(spec_loc_sz) ;
422 plines = cpl_vector_get_data(lines) ;
423 nlines = 0 ;
424
425 /* Loop as long as there are lines */
426 while ((max = cpl_vector_get_max(spec_loc)) > threshold) {
427 /* Find the max position */
428 int max_ind = 0 ;
429 while (max_ind < spec_loc_sz && pspec_loc[max_ind] < max) max_ind++ ;
430 if (max_ind == spec_loc_sz) {
431 cpl_msg_error(cpl_func, "Cannot find maximum") ;
432 cpl_vector_delete(spec_loc) ;
433 cpl_vector_delete(lines) ;
434 return NULL ;
435 }
436 if (max_ind == 0 || max_ind == spec_loc_sz-1) {
437 pspec_loc[max_ind] = 0 ;
438 continue ;
439 }
440
441 /* Get the precise position from the neighbours values */
442 plines[nlines] = pspec_loc[max_ind] * max_ind +
443 pspec_loc[max_ind-1] * (max_ind-1) +
444 pspec_loc[max_ind+1] * (max_ind+1) ;
445 plines[nlines] /= pspec_loc[max_ind] + pspec_loc[max_ind+1] +
446 pspec_loc[max_ind-1] ;
447 plines[nlines] ++ ;
448 nlines ++ ;
449
450 /* Clean the line */
451 i = max_ind ;
452 while (i>=0 && pspec_loc[i] > threshold) {
453 pspec_loc[i] = 0.0 ;
454 i-- ;
455 }
456 i = max_ind+1 ;
457 while (i<spec_loc_sz && pspec_loc[i] > threshold) {
458 pspec_loc[i] = 0.0 ;
459 i++ ;
460 }
461 }
462 cpl_vector_delete(spec_loc) ;
463
464 /* Check if there are lines */
465 if (nlines == 0) {
466 cpl_msg_error(cpl_func, "Cannot detect any line") ;
467 cpl_vector_delete(lines) ;
468 return NULL ;
469 }
470
471 /* Resize the vector */
472 cpl_vector_set_size(lines, nlines) ;
473
474 /* Sort the lines */
475 cpl_vector_sort(lines, 1) ;
476
477 return lines ;
478}
479
480#endif
irplib_ppm_engine
cpl_polynomial * irplib_ppm_engine(const cpl_vector *spectrum, const cpl_bivector *lines_catalog, const cpl_polynomial *poly_init, double slitw, double fwhm, double thresh, int degree, int doplot, cpl_table **tab_infos)
The Wavelength Calibration using PPM.
Definition irplib_ppm.c:69
irplib_spectrum_detect_peaks
cpl_vector * irplib_spectrum_detect_peaks(const cpl_vector *in, int fwhm, double sigma, int display, cpl_vector **fwhms_out, cpl_vector **areas_out)
Detect the brightest features in a spectrum.
Definition irplib_spectrum.c:258
Generated on Tue Jun 16 2026 00:11:35 for CR2RE Pipeline Reference Manual by doxygen 1.13.2