CR2RE Pipeline Reference Manual 1.7.0
irplib_spectrum.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_wlxcorr.h"
29#include "irplib_spectrum.h"
30
31#include <math.h>
32#include <float.h>
33#include <cpl.h>
34
35/*-----------------------------------------------------------------------------
36 Define
37 -----------------------------------------------------------------------------*/
38
39#define SPECTRUM_HW 16
40#define MIN_THRESH_FACT 0.9
41#define MAX_THRESH_FACT 1.1
42#define SPEC_SHADOW_FACT 30.0 /* Negative spectrum intensity*/
43#define SPEC_MAXWIDTH 48
44
45/*-----------------------------------------------------------------------------
46 Functions prototypes
47 -----------------------------------------------------------------------------*/
48
49static int select_valid_spectra(cpl_image *,
50 cpl_apertures *,
51 int,
52 spec_shadows,
53 int,
54 int *,
55 int **);
56static int
57valid_spectrum(cpl_image *, cpl_apertures *, int, spec_shadows, int, int);
58
59/*----------------------------------------------------------------------------*/
63/*----------------------------------------------------------------------------*/
64
66
67/*----------------------------------------------------------------------------*/
82/*----------------------------------------------------------------------------*/
83int
84irplib_spectrum_find_brightest(const cpl_image *in,
85 int offset,
86 spec_shadows shadows,
87 double min_bright,
88 int orient,
89 double *pos)
90{
91 cpl_image *loc_ima;
92 cpl_image *filt_image;
93 cpl_image *collapsed;
94 float *pcollapsed;
95 cpl_vector *line;
96 double *pline;
97 cpl_vector *line_filt;
98 double threshold;
99 double median, stdev, max, mean;
100 cpl_mask *mask;
101 cpl_image *labels;
102 cpl_size nlabels;
103 cpl_apertures *aperts;
104 int n_valid_specs;
105 int *valid_specs;
106 double brightness;
107 int i;
108
109 /* Test entries */
110 if (in == NULL)
111 return -1;
112 if (orient != 0 && orient != 1)
113 return -1;
114
115 /* Flip the image if necessary */
116 if (orient == 1) {
117 loc_ima = cpl_image_duplicate(in);
118 cpl_image_flip(loc_ima, 1);
119 }
120 else {
121 loc_ima = cpl_image_duplicate(in);
122 }
123
124 /* Median vertical filtering 3x3 */
125 mask = cpl_mask_new(3, 3);
126 cpl_mask_not(mask);
127 filt_image = cpl_image_new(cpl_image_get_size_x(loc_ima),
128 cpl_image_get_size_y(loc_ima),
129 cpl_image_get_type(loc_ima));
130 if (cpl_image_filter_mask(filt_image, loc_ima, mask, CPL_FILTER_MEDIAN,
131 CPL_BORDER_FILTER) != CPL_ERROR_NONE) {
132 cpl_msg_error(__func__, "Cannot filter the image");
133 cpl_mask_delete(mask);
134 cpl_image_delete(filt_image);
135 return -1;
136 }
137 cpl_mask_delete(mask);
138 cpl_image_delete(loc_ima);
139
140 /* Collapse the image */
141 if ((collapsed = cpl_image_collapse_median_create(filt_image, 1, 0, 0)) ==
142 NULL) {
143 cpl_msg_error(cpl_func,
144 "collapsing image: aborting spectrum detection");
145 cpl_image_delete(filt_image);
146 return -1;
147 }
148 cpl_image_delete(filt_image);
149
150 /* Subtract low frequency signal */
151 line = cpl_vector_new_from_image_column(collapsed, 1);
152 cpl_image_delete(collapsed);
153 line_filt = cpl_vector_filter_median_create(line, SPECTRUM_HW);
154 cpl_vector_subtract(line, line_filt);
155 cpl_vector_delete(line_filt);
156
157 /* Get relevant stats for thresholding */
158 median = cpl_vector_get_median_const(line);
159 stdev = cpl_vector_get_stdev(line);
160 max = cpl_vector_get_max(line);
161 mean = cpl_vector_get_mean(line);
162
163 /* Set the threshold */
164 threshold = median + stdev;
165 if (threshold > MIN_THRESH_FACT * max)
166 threshold = MIN_THRESH_FACT * max;
167 if (threshold < MAX_THRESH_FACT * mean)
168 threshold = MAX_THRESH_FACT * mean;
169
170 /* Recreate the image */
171 collapsed = cpl_image_new(1, cpl_vector_get_size(line), CPL_TYPE_FLOAT);
172 pcollapsed = cpl_image_get_data_float(collapsed);
173 pline = cpl_vector_get_data(line);
174 for (i = 0; i < cpl_vector_get_size(line); i++)
175 pcollapsed[i] = (float)pline[i];
176 cpl_vector_delete(line);
177
178 /* Binarise the image */
179 if ((mask = cpl_mask_threshold_image_create(collapsed, threshold,
180 DBL_MAX)) == NULL) {
181 cpl_msg_error(cpl_func, "cannot binarise");
182 cpl_image_delete(collapsed);
183 return -1;
184 }
185 if (cpl_mask_count(mask) < 1) {
186 cpl_msg_error(cpl_func, "not enough signal to detect spectra");
187 cpl_image_delete(collapsed);
188 cpl_mask_delete(mask);
189 return -1;
190 }
191 /* Labelise the different detected apertures */
192 if ((labels = cpl_image_labelise_mask_create(mask, &nlabels)) == NULL) {
193 cpl_msg_error(cpl_func, "cannot labelise");
194 cpl_image_delete(collapsed);
195 cpl_mask_delete(mask);
196 return -1;
197 }
198 cpl_mask_delete(mask);
199
200 /* Create the detected apertures list */
201 if ((aperts = cpl_apertures_new_from_image(collapsed, labels)) == NULL) {
202 cpl_msg_error(cpl_func, "cannot compute apertures");
203 cpl_image_delete(collapsed);
204 cpl_image_delete(labels);
205 return -1;
206 }
207 cpl_image_delete(labels);
208
209 /* Select only relevant specs, create corresponding LUT's */
210 if (select_valid_spectra(collapsed, aperts, offset, shadows, SPEC_MAXWIDTH,
211 &n_valid_specs, &valid_specs) == -1) {
212 cpl_msg_debug(
213 cpl_func,
214 "Could not select valid spectra from the %" CPL_SIZE_FORMAT
215 " apertures in %" CPL_SIZE_FORMAT "-col 1D-image, offset=%d"
216 ", min_bright=%d",
217 cpl_apertures_get_size(aperts), cpl_image_get_size_y(collapsed),
218 offset, SPEC_MAXWIDTH);
219 if (cpl_msg_get_level() <= CPL_MSG_DEBUG)
220 cpl_apertures_dump(aperts, stderr);
221 cpl_image_delete(collapsed);
222 cpl_apertures_delete(aperts);
223 return -1;
224 }
225 cpl_image_delete(collapsed);
226 if (n_valid_specs < 1) {
227 cpl_msg_error(cpl_func, "no valid spectrum detected");
228 cpl_free(valid_specs);
229 cpl_apertures_delete(aperts);
230 return -1;
231 }
232
233 /* Look for the brightest, among the detected spectra */
234 *pos = cpl_apertures_get_centroid_y(aperts, valid_specs[0] + 1);
235 brightness = cpl_apertures_get_flux(aperts, valid_specs[0] + 1);
236 for (i = 0; i < n_valid_specs; i++) {
237 if (cpl_apertures_get_flux(aperts, valid_specs[i] + 1) > brightness) {
238 *pos = cpl_apertures_get_centroid_y(aperts, valid_specs[i] + 1);
239 brightness = cpl_apertures_get_flux(aperts, valid_specs[i] + 1);
240 }
241 }
242 cpl_apertures_delete(aperts);
243 cpl_free(valid_specs);
244
245 /* Minimum brightness required */
246 if (brightness < min_bright) {
247 cpl_msg_error(cpl_func, "brightness %f too low <%f", brightness,
248 min_bright);
249 return -1;
250 }
251
252 /* Return */
253 return 0;
254}
255
256/*----------------------------------------------------------------------------*/
268/*----------------------------------------------------------------------------*/
269cpl_vector *
270irplib_spectrum_detect_peaks(const cpl_vector *in,
271 int fwhm,
272 double sigma,
273 int display,
274 cpl_vector **fwhms_out,
275 cpl_vector **areas_out)
276{
277 cpl_vector *filtered;
278 cpl_vector *spec_clean;
279 cpl_vector *spec_convolved;
280 double *pspec_convolved;
281 int filt_size;
282 cpl_vector *big_detected;
283 cpl_vector *big_fwhms;
284 cpl_vector *big_area;
285 double *pbig_detected;
286 double *pbig_fwhms;
287 double *pbig_area;
288 cpl_vector *detected;
289 cpl_vector *fwhms;
290 cpl_vector *area;
291 double max, med, stdev;
292 double x0, sig, norm, offset;
293 int nb_det, nb_samples, hwidth, start, stop;
294 int i, j;
295
296 /* Test entries */
297 if (in == NULL)
298 return NULL;
299
300 /* Initialise */
301 nb_samples = cpl_vector_get_size(in);
302 filt_size = 5;
303 hwidth = 5;
304
305 /* Subtract the low frequency part */
306 cpl_msg_debug(__func__, "Low Frequency signal removal");
307 if ((filtered = cpl_vector_filter_median_create(in, filt_size)) == NULL) {
308 cpl_msg_error(__func__, "Cannot filter the spectrum");
309 return NULL;
310 }
311 spec_clean = cpl_vector_duplicate(in);
312 cpl_vector_subtract(spec_clean, filtered);
313 cpl_vector_delete(filtered);
314
315 /* Display if requested */
316 if (display) {
317 cpl_plot_vector("set grid;set xlabel 'Position (pixels)';set ylabel "
318 "'Intensity (ADU)';",
319 "t 'Filtered extracted spectrum' w lines", "",
320 spec_clean);
321 }
322
323 /* Convolve */
324 spec_convolved = cpl_vector_duplicate(spec_clean);
325 if (fwhm > 0) {
326 cpl_vector *conv_kernel;
327 cpl_msg_debug(__func__, "Spectrum convolution");
328 /* Create convolution kernel */
329 if ((conv_kernel = irplib_wlxcorr_convolve_create_kernel(fwhm, fwhm)) ==
330 NULL) {
331 cpl_msg_error(cpl_func, "Cannot create convolution kernel");
332 cpl_vector_delete(spec_clean);
333 cpl_vector_delete(spec_convolved);
334 return NULL;
335 }
336
337 /* Smooth the instrument resolution */
338 if (irplib_wlxcorr_convolve(spec_convolved, conv_kernel)) {
339 cpl_msg_error(cpl_func, "Cannot smoothe the signal");
340 cpl_vector_delete(spec_clean);
341 cpl_vector_delete(spec_convolved);
342 cpl_vector_delete(conv_kernel);
343 return NULL;
344 }
345 cpl_vector_delete(conv_kernel);
346
347 /* Display if requested */
348 if (display) {
349 cpl_plot_vector("set grid;set xlabel 'Position (pixels)';set "
350 "ylabel 'Intensity (ADU)';",
351 "t 'Convolved extracted spectrum' w lines", "",
352 spec_convolved);
353 }
354 }
355
356 /* Apply the detection */
357 big_detected = cpl_vector_duplicate(spec_convolved);
358 big_fwhms = cpl_vector_duplicate(spec_convolved);
359 big_area = cpl_vector_duplicate(spec_convolved);
360 pbig_detected = cpl_vector_get_data(big_detected);
361 pbig_fwhms = cpl_vector_get_data(big_fwhms);
362 pbig_area = cpl_vector_get_data(big_area);
363
364 pspec_convolved = cpl_vector_get_data(spec_convolved);
365
366 /* To avoid detection on the side */
367 pspec_convolved[0] = pspec_convolved[nb_samples - 1] = 0.0;
368
369 /* Compute stats */
370 max = cpl_vector_get_max(spec_convolved);
371 stdev = cpl_vector_get_stdev(spec_convolved);
372 med = cpl_vector_get_median_const(spec_convolved);
373
374 /* Loop on the detected lines */
375 nb_det = 0;
376 while (max > med + stdev * sigma) {
377 cpl_vector *extract;
378 cpl_vector *extract_x;
379 double cur_val;
380
381 /* Compute the position */
382 i = 0;
383 while (pspec_convolved[i] < max)
384 i++;
385 if (i <= 0 || i >= nb_samples - 1)
386 break;
387
388 /* Extract the line */
389 if (i - hwidth >= 0)
390 start = i - hwidth;
391 else
392 start = 0;
393 if (i + hwidth <= nb_samples - 1)
394 stop = i + hwidth;
395 else
396 stop = nb_samples - 1;
397 extract = cpl_vector_extract(spec_clean, start, stop, 1);
398 extract_x = cpl_vector_duplicate(extract);
399 for (j = 0; j < cpl_vector_get_size(extract_x); j++) {
400 cpl_vector_set(extract_x, j, (double)j + 1);
401 }
402 /* Fit the gaussian */
403 if (cpl_vector_fit_gaussian(extract_x, NULL, extract, NULL, CPL_FIT_ALL,
404 &x0, &sig, &norm, &offset, NULL, NULL,
405 NULL) != CPL_ERROR_NONE) {
406 cpl_msg_debug(__func__, "Cannot fit a gaussian at [%d, %d]", start,
407 stop);
408 cpl_error_reset();
409 }
410 else {
411 pbig_detected[nb_det] = x0 + start;
412 pbig_area[nb_det] = norm;
413 pbig_fwhms[nb_det] = 2 * sig * sqrt(2 * log(2));
414 cpl_msg_debug(__func__, "Line nb %d at position %g", nb_det + 1,
415 pbig_detected[nb_det]);
416 nb_det++;
417 }
418 cpl_vector_delete(extract);
419 cpl_vector_delete(extract_x);
420
421 /* Cancel out the line on the left */
422 j = i - 1;
423 cur_val = pspec_convolved[i];
424 while (j >= 0 && pspec_convolved[j] < cur_val) {
425 cur_val = pspec_convolved[j];
426 pspec_convolved[j] = 0.0;
427 j--;
428 }
429 /* Cancel out the line on the right */
430 j = i + 1;
431 cur_val = pspec_convolved[i];
432 while (j <= nb_samples - 1 && pspec_convolved[j] < cur_val) {
433 cur_val = pspec_convolved[j];
434 pspec_convolved[j] = 0.0;
435 j++;
436 }
437 /* Cancel out the line on center */
438 pspec_convolved[i] = 0.0;
439
440 /* Recompute the stats */
441 max = cpl_vector_get_max(spec_convolved);
442 stdev = cpl_vector_get_stdev(spec_convolved);
443 med = cpl_vector_get_median_const(spec_convolved);
444 }
445 cpl_vector_delete(spec_convolved);
446 cpl_vector_delete(spec_clean);
447
448 /* Create the output vector */
449 if (nb_det == 0) {
450 detected = NULL;
451 area = NULL;
452 fwhms = NULL;
453 }
454 else {
455 double *pdetected;
456 double *pfwhms;
457 double *parea;
458 detected = cpl_vector_new(nb_det);
459 area = cpl_vector_new(nb_det);
460 fwhms = cpl_vector_new(nb_det);
461 pdetected = cpl_vector_get_data(detected);
462 parea = cpl_vector_get_data(area);
463 pfwhms = cpl_vector_get_data(fwhms);
464 for (i = 0; i < nb_det; i++) {
465 pdetected[i] = pbig_detected[i];
466 parea[i] = pbig_area[i];
467 pfwhms[i] = pbig_fwhms[i];
468 }
469 }
470 cpl_vector_delete(big_detected);
471 cpl_vector_delete(big_area);
472 cpl_vector_delete(big_fwhms);
473
474 /* Return */
475 if (fwhms_out == NULL)
476 cpl_vector_delete(fwhms);
477 else
478 *fwhms_out = fwhms;
479 if (areas_out == NULL)
480 cpl_vector_delete(area);
481 else
482 *areas_out = area;
483 return detected;
484}
485
487
488/*----------------------------------------------------------------------------*/
500/*----------------------------------------------------------------------------*/
501static int
502select_valid_spectra(cpl_image *in,
503 cpl_apertures *aperts,
504 int offset,
505 spec_shadows shadows,
506 int max_spec_width,
507 int *n_valid_specs,
508 int **valid_specs)
509{
510 int nb_aperts;
511 int i, j;
512
513 /* Initialise */
514 *valid_specs = NULL;
515 nb_aperts = cpl_apertures_get_size(aperts);
516 *n_valid_specs = 0;
517
518 /* Test entries */
519 if (nb_aperts < 1)
520 return -1;
521
522 /* Count nb of valid specs */
523 j = 0;
524 for (i = 0; i < nb_aperts; i++)
525 if (valid_spectrum(in, aperts, offset, shadows, max_spec_width, i + 1))
526 (*n_valid_specs)++;
527
528 /* Associate to each spectrum, its object number */
529 if (*n_valid_specs) {
530 *valid_specs = cpl_calloc(*n_valid_specs, sizeof(int));
531 j = 0;
532 for (i = 0; i < nb_aperts; i++)
533 if (valid_spectrum(in, aperts, offset, shadows, max_spec_width,
534 i + 1)) {
535 (*valid_specs)[j] = i;
536 j++;
537 }
538 }
539 else
540 return -1;
541
542 return 0;
543}
544
545/*---------------------------------------------------------------------------*/
556/*----------------------------------------------------------------------------*/
557static int
558valid_spectrum(cpl_image *in,
559 cpl_apertures *aperts,
560 int offset,
561 spec_shadows shadows,
562 int max_spec_width,
563 int objnum)
564{
565 int objwidth;
566 double valover, valunder, valcenter;
567
568 /* Find objwidth */
569 objwidth = cpl_apertures_get_top(aperts, objnum) -
570 cpl_apertures_get_bottom(aperts, objnum) + 1;
571 if (objwidth > max_spec_width) {
572 cpl_msg_error(cpl_func, "object is too wide");
573 return 0;
574 }
575
576 /* Object is too small */
577 if (cpl_apertures_get_npix(aperts, objnum) < 2)
578 return 0;
579
580 /* no shadow required */
581 if (shadows == NO_SHADOW)
582 return 1;
583
584 /* Get the median of the object (valcenter) */
585 valcenter = cpl_apertures_get_median(aperts, objnum);
586
587 /* Get the black shadows medians (valunder and valover) */
588 if (cpl_apertures_get_bottom(aperts, objnum) - offset < 1)
589 valunder = 0.0;
590 else
591 valunder = cpl_image_get_median_window(
592 in, 1, cpl_apertures_get_bottom(aperts, objnum) - offset, 1,
593 cpl_apertures_get_top(aperts, objnum) - offset);
594
595 if (cpl_apertures_get_top(aperts, objnum) + offset > 1024)
596 valover = 0.0;
597 else
598 valover = cpl_image_get_median_window(
599 in, 1, cpl_apertures_get_bottom(aperts, objnum) + offset, 1,
600 cpl_apertures_get_top(aperts, objnum) + offset);
601
602 switch (shadows) {
603 case TWO_SHADOWS:
604 if ((valunder < -fabs(valcenter / SPEC_SHADOW_FACT)) &&
605 (valover < -fabs(valcenter / SPEC_SHADOW_FACT)) &&
606 (valunder / valover > 0.5) && (valunder / valover < 2.0))
607 return 1;
608 break;
609
610 case ONE_SHADOW:
611 if ((valunder < -fabs(valcenter / SPEC_SHADOW_FACT)) ||
612 (valover < -fabs(valcenter / SPEC_SHADOW_FACT)))
613 return 1;
614 break;
615
616 case NO_SHADOW:
617 return 1;
618
619 default:
620 cpl_msg_error(cpl_func, "unknown spec_detect_mode");
621 break;
622 }
623
624 cpl_msg_debug(cpl_func, "No spectrum(%d): under=%g, center=%g, over=%g",
625 shadows, valunder, valcenter, valover);
626
627 return 0;
628}
irplib_spectrum_find_brightest
int irplib_spectrum_find_brightest(const cpl_image *in, int offset, spec_shadows shadows, double min_bright, int orient, double *pos)
Finds the brightest spectrum in an image.
Definition irplib_spectrum.c:84
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:270
Generated on Tue Oct 13 2026 12:27:26 for CR2RE Pipeline Reference Manual by doxygen 1.13.2