CR2RE Pipeline Reference Manual 1.7.0
irplib_wavecal.c
1/*
2 * This file is part of the IRPLIB Pipeline
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_wavecal_impl.h"
29
30/* Needed for irplib_errorstate_dump_debug() */
31#include "irplib_utils.h"
32
33#include <string.h>
34#include <math.h>
35
36#ifdef HAVE_GSL
37#include <gsl/gsl_multimin.h>
38#endif
39
40/*-----------------------------------------------------------------------------
41 Private types
42 -----------------------------------------------------------------------------*/
43
44#ifdef HAVE_GSL
45
46typedef struct
47{
48 const cpl_vector *observed;
49 cpl_polynomial *disp1d;
50 cpl_vector *spectrum;
51 irplib_base_spectrum_model *param;
52 cpl_error_code (*filler)(cpl_vector *,
53 const cpl_polynomial *,
54 irplib_base_spectrum_model *);
55 cpl_vector *vxc;
56 double xc;
57 int maxxc;
58 double mxc;
59 cpl_polynomial *mdisp;
60 int ishift;
61
62} irplib_multimin;
63
64#endif /* HAVE_GSL */
65
66/*-----------------------------------------------------------------------------
67 Defines
68 -----------------------------------------------------------------------------*/
69
70#ifndef inline
71#define inline /* inline */
72#endif
73
74#define IRPLIB_MAX(A, B) ((A) > (B) ? (A) : (B))
75#define IRPLIB_MIN(A, B) ((A) < (B) ? (A) : (B))
76
77/*-----------------------------------------------------------------------------
78 Private functions
79 -----------------------------------------------------------------------------*/
80
81#ifdef HAVE_GSL
82static double irplib_gsl_correlation(const gsl_vector *, void *);
83#endif
84
85static cpl_error_code irplib_polynomial_find_1d_from_correlation_(
86 cpl_polynomial *,
87 int,
88 const cpl_vector *,
89 irplib_base_spectrum_model *,
90 cpl_error_code (*)(cpl_vector *,
91 const cpl_polynomial *,
92 irplib_base_spectrum_model *),
93 double,
94 double,
95 int,
96 int,
97 double *,
98 cpl_boolean *);
99
100
101/*----------------------------------------------------------------------------*/
105/*----------------------------------------------------------------------------*/
106
108
109
110/*----------------------------------------------------------------------------*/
118/*----------------------------------------------------------------------------*/
119int
120irplib_bivector_count_positive(const cpl_bivector *self,
121 double x_min,
122 double x_max)
123{
124 const int nself = cpl_bivector_get_size(self);
125 const double *px = cpl_bivector_get_x_data_const(self);
126 const double *py = cpl_bivector_get_y_data_const(self);
127 int npos = 0;
128 int i = 0;
129
130 cpl_ensure(self != NULL, CPL_ERROR_NULL_INPUT, -1);
131 cpl_ensure(x_min <= x_max, CPL_ERROR_ILLEGAL_INPUT, -2);
132
133 /* FIXME: Use cpl_vector_find() */
134 while (i < nself && px[i] < x_min)
135 i++;
136 while (i < nself && px[i] < x_max)
137 if (py[i++] > 0)
138 npos++;
139
140 return npos;
141}
142
143/*----------------------------------------------------------------------------*/
153/*----------------------------------------------------------------------------*/
154cpl_error_code
155irplib_polynomial_fit_2d_dispersion(cpl_polynomial *self,
156 const cpl_image *imgwave,
157 int fitdeg,
158 double *presid)
159{
160 const int nx = cpl_image_get_size_x(imgwave);
161 const int ny = cpl_image_get_size_y(imgwave);
162 const int nbad = cpl_image_count_rejected(imgwave);
163 const int nsamp = nx * ny - nbad;
164 cpl_matrix *xy_pos;
165 double *xdata;
166 double *ydata;
167 cpl_vector *wlen;
168 double *dwlen;
169 const cpl_size nfitdeg = (cpl_size)fitdeg;
170 int i, j;
171 int k = 0;
172
173 cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT);
174 cpl_ensure_code(imgwave != NULL, CPL_ERROR_NULL_INPUT);
175 cpl_ensure_code(presid != NULL, CPL_ERROR_NULL_INPUT);
176 cpl_ensure_code(fitdeg > 0, CPL_ERROR_ILLEGAL_INPUT);
177
178 cpl_ensure_code(cpl_polynomial_get_dimension(self) == 2,
179 CPL_ERROR_ILLEGAL_INPUT);
180
181 xy_pos = cpl_matrix_new(2, nsamp);
182 xdata = cpl_matrix_get_data(xy_pos);
183 ydata = xdata + nsamp;
184
185 dwlen = (double *)cpl_malloc(nsamp * sizeof(double));
186 wlen = cpl_vector_wrap(nsamp, dwlen);
187
188 for (i = 1; i <= nx; i++) {
189 for (j = 1; j <= ny; j++) {
190 int is_bad;
191 const double value = cpl_image_get(imgwave, i, j, &is_bad);
192 if (!is_bad) {
193 xdata[k] = i;
194 ydata[k] = j;
195 dwlen[k] = value;
196 k++;
197 }
198 }
199 }
200
201 cpl_msg_info(cpl_func,
202 "Fitting 2D polynomial to %d X %d image, ignoring "
203 "%d poorly calibrated pixels",
204 nx, ny, nbad);
205
206 if (cpl_polynomial_fit(self, xy_pos, NULL, wlen, NULL, CPL_FALSE, NULL,
207 &nfitdeg) == CPL_ERROR_NONE &&
208 presid != NULL) {
209 cpl_vector_fill_polynomial_fit_residual(wlen, wlen, NULL, self, xy_pos,
210 NULL);
211 *presid = cpl_vector_product(wlen, wlen) / nsamp;
212 }
213 cpl_matrix_delete(xy_pos);
214 cpl_vector_delete(wlen);
215
216 cpl_ensure_code(k == nsamp, CPL_ERROR_UNSPECIFIED);
217
218 return CPL_ERROR_NONE;
219}
220
221
222/*----------------------------------------------------------------------------*/
240/*----------------------------------------------------------------------------*/
241cpl_error_code
242irplib_polynomial_find_1d_from_correlation(
243 cpl_polynomial *self,
244 int maxdeg,
245 const cpl_vector *obs,
246 irplib_base_spectrum_model *model,
247 cpl_error_code (*filler)(cpl_vector *,
248 const cpl_polynomial *,
249 irplib_base_spectrum_model *),
250 double pixtol,
251 double pixstep,
252 int hsize,
253 int maxite,
254 double *pxc)
255{
256 cpl_boolean restart = CPL_FALSE;
257 const cpl_error_code error =
258 irplib_polynomial_find_1d_from_correlation_(self, maxdeg, obs, model,
259 filler, pixtol, pixstep,
260 hsize, maxite, pxc,
261 &restart);
262
263 return error ? cpl_error_set_where(cpl_func)
264 : (restart ? cpl_error_set(cpl_func, CPL_ERROR_CONTINUE)
265 : CPL_ERROR_NONE);
266}
267
268/*----------------------------------------------------------------------------*/
289/*----------------------------------------------------------------------------*/
290static cpl_error_code
291irplib_polynomial_find_1d_from_correlation_(
292 cpl_polynomial *self,
293 int maxdeg,
294 const cpl_vector *obs,
295 irplib_base_spectrum_model *model,
296 cpl_error_code (*filler)(cpl_vector *,
297 const cpl_polynomial *,
298 irplib_base_spectrum_model *),
299 double pixtol,
300 double pixstep,
301 int hsize,
302 int maxite,
303 double *pxc,
304 cpl_boolean *prestart)
305{
306#ifdef HAVE_GSL
307 const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
308 gsl_multimin_fminimizer *minimizer;
309 gsl_multimin_function my_func;
310 irplib_multimin data;
311 gsl_vector *dispgsl;
312 gsl_vector *stepsize;
313 gsl_vector *dispprev;
314 int status = GSL_CONTINUE;
315 const int nobs = cpl_vector_get_size(obs);
316 const cpl_size nfit = maxdeg + 1;
317 cpl_errorstate prestate = cpl_errorstate_get();
318 /* Convert pixel step to wavelength step on detector center */
319 const double wlstep =
320 cpl_polynomial_eval_1d_diff(self, 0.5 * (nobs + pixstep),
321 0.5 * (nobs - pixstep), NULL);
322 double wlstepi = wlstep;
323 int iter;
324 cpl_size i;
325
326#endif
327
328 cpl_ensure_code(prestart != NULL, CPL_ERROR_NULL_INPUT);
329 *prestart = CPL_FALSE;
330 cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT);
331 cpl_ensure_code(obs != NULL, CPL_ERROR_NULL_INPUT);
332 cpl_ensure_code(model != NULL, CPL_ERROR_NULL_INPUT);
333 cpl_ensure_code(filler != NULL, CPL_ERROR_NULL_INPUT);
334 cpl_ensure_code(pxc != NULL, CPL_ERROR_NULL_INPUT);
335
336 cpl_ensure_code(cpl_polynomial_get_dimension(self) == 1,
337 CPL_ERROR_ILLEGAL_INPUT);
338
339 cpl_ensure_code(cpl_polynomial_get_degree(self) > 0,
340 CPL_ERROR_ILLEGAL_INPUT);
341
342 cpl_ensure_code(maxdeg >= 0, CPL_ERROR_ILLEGAL_INPUT);
343 cpl_ensure_code(pixtol > 0.0, CPL_ERROR_ILLEGAL_INPUT);
344 cpl_ensure_code(pixstep > 0.0, CPL_ERROR_ILLEGAL_INPUT);
345 cpl_ensure_code(hsize >= 0, CPL_ERROR_ILLEGAL_INPUT);
346 cpl_ensure_code(maxite >= 0, CPL_ERROR_ILLEGAL_INPUT);
347
348#ifndef HAVE_GSL
349 return cpl_error_set_message(cpl_func, CPL_ERROR_UNSUPPORTED_MODE,
350 "GSL is not available");
351#else
352
353 minimizer = gsl_multimin_fminimizer_alloc(T, (size_t)nfit);
354
355 cpl_ensure_code(minimizer != NULL, CPL_ERROR_ILLEGAL_OUTPUT);
356
357 dispgsl = gsl_vector_alloc((size_t)nfit);
358 stepsize = gsl_vector_alloc((size_t)nfit);
359 dispprev = gsl_vector_alloc((size_t)nfit);
360
361 for (i = 0; i < nfit; i++) {
362 const double value = cpl_polynomial_get_coeff(self, &i);
363 gsl_vector_set(dispgsl, (size_t)i, value);
364 gsl_vector_set(stepsize, (size_t)i, wlstepi);
365 wlstepi /= (double)nobs;
366 }
367
368 my_func.n = nfit;
369 my_func.f = &irplib_gsl_correlation;
370 my_func.params = (void *)(&data);
371
372 data.observed = obs;
373 data.disp1d = self;
374 data.spectrum = cpl_vector_new(nobs + 2 * hsize);
375 data.vxc = cpl_vector_new(1 + 2 * hsize);
376 data.xc = 0;
377 data.param = model;
378 data.filler = filler;
379 data.maxxc = 0; /* Output */
380 data.ishift = 0; /* Output */
381 data.mxc = -1.0; /* Output */
382 data.mdisp = NULL; /* Output */
383
384 gsl_multimin_fminimizer_set(minimizer, &my_func, dispgsl, stepsize);
385
386 for (iter = 0; status == GSL_CONTINUE && iter < maxite; iter++) {
387 double size;
388 const double fprev = minimizer->fval;
389
390 gsl_vector_memcpy(dispprev, minimizer->x);
391 status = gsl_multimin_fminimizer_iterate(minimizer);
392
393 if (status || !cpl_errorstate_is_equal(prestate))
394 break;
395
396 size = gsl_multimin_fminimizer_size(minimizer);
397 status = gsl_multimin_test_size(size, pixtol);
398
399 if (status == GSL_SUCCESS) {
400 cpl_msg_debug(cpl_func, "converged to minimum at");
401
402 if (nfit == 0) {
403 cpl_msg_debug(cpl_func, "%5d %g df() = %g size = %g", iter,
404 gsl_vector_get(minimizer->x, 0) -
405 gsl_vector_get(dispprev, 0),
406 minimizer->fval - fprev, size);
407 }
408 else if (nfit == 1) {
409 cpl_msg_debug(cpl_func, "%5d %g %g df() = %g size = %g", iter,
410 gsl_vector_get(minimizer->x, 0) -
411 gsl_vector_get(dispprev, 0),
412 gsl_vector_get(minimizer->x, 1) -
413 gsl_vector_get(dispprev, 1),
414 minimizer->fval - fprev, size);
415 }
416 else {
417 cpl_msg_debug(cpl_func, "%5d %g %g %g df() = %g size = %g",
418 iter,
419 gsl_vector_get(minimizer->x, 0) -
420 gsl_vector_get(dispprev, 0),
421 gsl_vector_get(minimizer->x, 1) -
422 gsl_vector_get(dispprev, 1),
423 gsl_vector_get(minimizer->x, 2) -
424 gsl_vector_get(dispprev, 2),
425 minimizer->fval - fprev, size);
426 }
427 }
428 }
429
430 if (status == GSL_SUCCESS && cpl_errorstate_is_equal(prestate)) {
431 if (data.mxc > -minimizer->fval) {
432 *pxc = data.mxc;
433 cpl_msg_warning(cpl_func, "Local maximum: %g(%d) > %g", data.mxc,
434 data.ishift, -minimizer->fval);
435 cpl_polynomial_shift_1d(data.mdisp, 0, (double)data.ishift);
436 cpl_polynomial_copy(self, data.mdisp);
437 *prestart = CPL_TRUE;
438 }
439 else {
440 *pxc = -minimizer->fval;
441 for (i = 0; i < nfit; i++) {
442 const double value = gsl_vector_get(minimizer->x, i);
443 cpl_polynomial_set_coeff(self, &i, value);
444 }
445 }
446 }
447
448 cpl_vector_delete(data.spectrum);
449 cpl_vector_delete(data.vxc);
450 cpl_polynomial_delete(data.mdisp);
451 gsl_multimin_fminimizer_free(minimizer);
452 gsl_vector_free(dispgsl);
453 gsl_vector_free(dispprev);
454 gsl_vector_free(stepsize);
455
456 cpl_ensure_code(status != GSL_CONTINUE, CPL_ERROR_CONTINUE);
457 cpl_ensure_code(status == GSL_SUCCESS, CPL_ERROR_DATA_NOT_FOUND);
458 cpl_ensure_code(cpl_errorstate_is_equal(prestate), cpl_error_get_code());
459
460 return CPL_ERROR_NONE;
461#endif
462}
463
464
465/*----------------------------------------------------------------------------*/
493/*----------------------------------------------------------------------------*/
494cpl_error_code
495irplib_vector_fill_line_spectrum(cpl_vector *self,
496 const cpl_polynomial *disp,
497 irplib_base_spectrum_model *lsslamp)
498{
499 irplib_line_spectrum_model *arclamp = (irplib_line_spectrum_model *)lsslamp;
500 cpl_error_code error;
501
502 cpl_ensure_code(arclamp != NULL, CPL_ERROR_NULL_INPUT);
503
504 arclamp->cost++;
505
506 error =
507 irplib_vector_fill_line_spectrum_model(self, arclamp->linepix,
508 arclamp->erftmp, disp,
509 arclamp->lines, arclamp->wslit,
510 arclamp->wfwhm, arclamp->xtrunc,
511 0, CPL_FALSE, CPL_FALSE,
512 &(arclamp->ulines));
513 cpl_ensure_code(!error, error);
514
515 arclamp->xcost++;
516
517 return CPL_ERROR_NONE;
518}
519
520/*----------------------------------------------------------------------------*/
533/*----------------------------------------------------------------------------*/
534cpl_error_code
535irplib_vector_fill_logline_spectrum(cpl_vector *self,
536 const cpl_polynomial *disp,
537 irplib_base_spectrum_model *lsslamp)
538{
539 irplib_line_spectrum_model *arclamp = (irplib_line_spectrum_model *)lsslamp;
540 cpl_error_code error;
541
542 cpl_ensure_code(arclamp != NULL, CPL_ERROR_NULL_INPUT);
543
544 arclamp->cost++;
545
546 error =
547 irplib_vector_fill_line_spectrum_model(self, arclamp->linepix,
548 arclamp->erftmp, disp,
549 arclamp->lines, arclamp->wslit,
550 arclamp->wfwhm, arclamp->xtrunc,
551 0, CPL_FALSE, CPL_TRUE,
552 &(arclamp->ulines));
553 cpl_ensure_code(!error, error);
554
555 arclamp->xcost++;
556
557 return CPL_ERROR_NONE;
558}
559
560
561/*----------------------------------------------------------------------------*/
574/*----------------------------------------------------------------------------*/
575cpl_error_code
576irplib_vector_fill_line_spectrum_fast(cpl_vector *self,
577 const cpl_polynomial *disp,
578 irplib_base_spectrum_model *lsslamp)
579{
580 irplib_line_spectrum_model *arclamp = (irplib_line_spectrum_model *)lsslamp;
581 cpl_error_code error;
582
583 cpl_ensure_code(arclamp != NULL, CPL_ERROR_NULL_INPUT);
584
585 arclamp->cost++;
586
587 error =
588 irplib_vector_fill_line_spectrum_model(self, arclamp->linepix,
589 arclamp->erftmp, disp,
590 arclamp->lines, arclamp->wslit,
591 arclamp->wfwhm, arclamp->xtrunc,
592 0, CPL_TRUE, CPL_FALSE,
593 &(arclamp->ulines));
594 cpl_ensure_code(!error, error);
595
596 arclamp->xcost++;
597
598 return CPL_ERROR_NONE;
599}
600
601/*----------------------------------------------------------------------------*/
614/*----------------------------------------------------------------------------*/
615cpl_error_code
616irplib_vector_fill_logline_spectrum_fast(cpl_vector *self,
617 const cpl_polynomial *disp,
618 irplib_base_spectrum_model *lsslamp)
619{
620 irplib_line_spectrum_model *arclamp = (irplib_line_spectrum_model *)lsslamp;
621 cpl_error_code error;
622
623 cpl_ensure_code(arclamp != NULL, CPL_ERROR_NULL_INPUT);
624
625 arclamp->cost++;
626
627 error =
628 irplib_vector_fill_line_spectrum_model(self, arclamp->linepix,
629 arclamp->erftmp, disp,
630 arclamp->lines, arclamp->wslit,
631 arclamp->wfwhm, arclamp->xtrunc,
632 0, CPL_TRUE, CPL_TRUE,
633 &(arclamp->ulines));
634 cpl_ensure_code(!error, error);
635
636 arclamp->xcost++;
637
638 return CPL_ERROR_NONE;
639}
640
641/*----------------------------------------------------------------------------*/
652/*----------------------------------------------------------------------------*/
653cpl_error_code
654irplib_plot_spectrum_and_model(
655 const cpl_vector *self,
656 const cpl_polynomial *disp1d,
657 irplib_base_spectrum_model *model,
658 cpl_error_code (*filler)(cpl_vector *,
659 const cpl_polynomial *,
660 irplib_base_spectrum_model *))
661{
662 cpl_errorstate prestate = cpl_errorstate_get();
663 cpl_vector *wl;
664 cpl_vector *spectrum;
665 cpl_vector *vxc;
666 const int len = cpl_vector_get_size(self);
667 double maxval, xc;
668 int ixc;
669 int error = 0;
670
671 cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT);
672 cpl_ensure_code(disp1d != NULL, CPL_ERROR_NULL_INPUT);
673 cpl_ensure_code(model != NULL, CPL_ERROR_NULL_INPUT);
674 cpl_ensure_code(filler != NULL, CPL_ERROR_NULL_INPUT);
675
676 cpl_ensure_code(cpl_polynomial_get_dimension(disp1d) == 1,
677 CPL_ERROR_ILLEGAL_INPUT);
678
679 cpl_ensure_code(cpl_polynomial_get_degree(disp1d) > 0,
680 CPL_ERROR_ILLEGAL_INPUT);
681
682 wl = cpl_vector_new(len);
683 spectrum = cpl_vector_new(len);
684 vxc = cpl_vector_new(1);
685
686 error |= (int)cpl_vector_fill_polynomial(wl, disp1d, 1.0, 1.0);
687 error |= filler(spectrum, disp1d, model);
688
689 ixc = cpl_vector_correlate(vxc, self, spectrum);
690 xc = cpl_vector_get(vxc, ixc);
691
692 maxval = cpl_vector_get_max(spectrum);
693 if (maxval != 0.0)
694 error |= cpl_vector_multiply_scalar(spectrum,
695 cpl_vector_get_max(self) / maxval);
696 if (!error) {
697 const cpl_vector *spair[] = {wl, self, spectrum};
698 char *pre =
699 cpl_sprintf("set grid;set xlabel 'Wavelength (%g -> %g)'; "
700 "set ylabel 'Intensity';",
701 cpl_vector_get(wl, 0), cpl_vector_get(wl, len - 1));
702 char *title = cpl_sprintf("t 'Observed and modelled spectra (%d pixel "
703 "XC=%g) ' w linespoints",
704 len, xc);
705
706 (void)cpl_plot_vectors(pre, title, "", spair, 3);
707 cpl_free(pre);
708 cpl_free(title);
709 }
710
711 cpl_vector_delete(wl);
712 cpl_vector_delete(spectrum);
713 cpl_vector_delete(vxc);
714
715 cpl_errorstate_set(prestate);
716
717 return CPL_ERROR_NONE;
718}
719
720/*----------------------------------------------------------------------------*/
740/*----------------------------------------------------------------------------*/
741cpl_error_code
742irplib_bivector_find_shift_from_correlation(
743 cpl_bivector *self,
744 const cpl_polynomial *disp,
745 const cpl_vector *obs,
746 irplib_base_spectrum_model *model,
747 cpl_error_code (*filler)(cpl_vector *,
748 const cpl_polynomial *,
749 irplib_base_spectrum_model *),
750 int hsize,
751 cpl_boolean doplot,
752 double *pxc)
753{
754 const int nobs = cpl_vector_get_size(obs);
755 const int nmodel = 2 * hsize + nobs;
756 cpl_polynomial *shdisp;
757 cpl_vector *xself = cpl_bivector_get_x(self);
758 cpl_vector *yself = cpl_bivector_get_y(self);
759 cpl_vector *mspec1d;
760 cpl_vector *xcorr;
761 cpl_error_code error = CPL_ERROR_NONE;
762 double xcprev, xcnext;
763 int ixc, imax = 0;
764 int i;
765
766 cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT);
767 cpl_ensure_code(disp != NULL, CPL_ERROR_NULL_INPUT);
768 cpl_ensure_code(obs != NULL, CPL_ERROR_NULL_INPUT);
769 cpl_ensure_code(model != NULL, CPL_ERROR_NULL_INPUT);
770 cpl_ensure_code(filler != NULL, CPL_ERROR_NULL_INPUT);
771 cpl_ensure_code(hsize > 0, CPL_ERROR_ILLEGAL_INPUT);
772
773 shdisp = cpl_polynomial_duplicate(disp);
774
775 /* Shift reference by -hsize so filler can be used without offset */
776 if (cpl_polynomial_shift_1d(shdisp, 0, -hsize)) {
777 cpl_polynomial_delete(shdisp);
778 return cpl_error_set_where(cpl_func);
779 }
780
781 mspec1d = cpl_vector_new(nmodel);
782
783 if (filler(mspec1d, shdisp, model)) {
784 cpl_vector_delete(mspec1d);
785 return cpl_error_set_where(cpl_func);
786 }
787
788 /* Should not be able to fail now */
789 xcorr = cpl_vector_new(1 + 2 * hsize);
790 ixc = cpl_vector_correlate(xcorr, mspec1d, obs);
791
792#ifdef IRPLIB_SPC_DUMP
793 /* Need irplib_wavecal.c rev. 1.12 through 1.15 */
794 irplib_polynomial_dump_corr_step(shdisp, xcorr, "Shift");
795#endif
796
797 cpl_vector_delete(mspec1d);
798 cpl_polynomial_delete(shdisp);
799
800 /* Find local maxima. */
801 /* FIXME(?): Also include stationary points */
802 i = 0;
803 xcprev = cpl_vector_get(xcorr, i);
804 xcnext = cpl_vector_get(xcorr, i + 1);
805
806 if (xcprev >= xcnext) {
807 /* 1st data point is an extreme */
808 /* FIXME: This could also be an error, recoverable by caller by
809 increasing hsize */
810 imax++;
811
812 cpl_vector_set(xself, 0, i - hsize);
813 cpl_vector_set(yself, 0, xcprev);
814 }
815
816 for (i = 1; i < 2 * hsize; i++) {
817 const double xc = xcnext;
818 xcnext = cpl_vector_get(xcorr, i + 1);
819 if (xc >= xcprev && xc >= xcnext) {
820 /* Found (local) maximum at shift i - hsize */
821 int j;
822
823 imax++;
824
825 if (cpl_bivector_get_size(self) < imax) {
826 cpl_vector_set_size(xself, imax);
827 cpl_vector_set_size(yself, imax);
828 }
829
830 for (j = imax - 1; j > 0; j--) {
831 if (xc <= cpl_vector_get(yself, j - 1))
832 break;
833 cpl_vector_set(xself, j, cpl_vector_get(xself, j - 1));
834 cpl_vector_set(yself, j, cpl_vector_get(yself, j - 1));
835 }
836 cpl_vector_set(xself, j, i - hsize);
837 cpl_vector_set(yself, j, xc);
838 }
839 xcprev = xc;
840 }
841
842 /* assert( i == 2 * hsize ); */
843
844 if (xcnext >= xcprev) {
845 /* Last data point is an extreme */
846 /* FIXME: This could also be an error, recoverable by caller by
847 increasing hsize */
848 int j;
849
850 imax++;
851
852 if (cpl_bivector_get_size(self) < imax) {
853 cpl_vector_set_size(xself, imax);
854 cpl_vector_set_size(yself, imax);
855 }
856
857 for (j = imax - 1; j > 0; j--) {
858 if (xcnext <= cpl_vector_get(yself, j - 1))
859 break;
860 cpl_vector_set(xself, j, cpl_vector_get(xself, j - 1));
861 cpl_vector_set(yself, j, cpl_vector_get(yself, j - 1));
862 }
863 cpl_vector_set(xself, j, i - hsize);
864 cpl_vector_set(yself, j, xcnext);
865 }
866
867 if (doplot) {
868 /* Vector of -hsize, 1-hsize, 2-hsize, ..., 0, ..., hsize */
869 cpl_vector *xvals = cpl_vector_new(1 + 2 * hsize);
870 cpl_bivector *bcorr = cpl_bivector_wrap_vectors(xvals, xcorr);
871 double x = (double)-hsize;
872 char *title =
873 cpl_sprintf("t 'Cross-correlation of shifted %d-pixel "
874 "spectrum (XCmax=%g at %d)' w linespoints",
875 nobs, cpl_vector_get(xcorr, ixc), ixc - hsize);
876
877 for (i = 0; i < 1 + 2 * hsize; i++, x += 1.0) {
878 cpl_vector_set(xvals, i, x);
879 }
880
881 cpl_plot_bivector("set grid;set xlabel 'Offset [pixel]';", title, "",
882 bcorr);
883 cpl_bivector_unwrap_vectors(bcorr);
884 cpl_vector_delete(xvals);
885 cpl_free(title);
886 }
887
888 if (pxc != NULL)
889 *pxc = cpl_vector_get(xcorr, hsize);
890
891 cpl_vector_delete(xcorr);
892
893 if (imax < 1) {
894 error = CPL_ERROR_DATA_NOT_FOUND;
895 }
896 else if (cpl_bivector_get_size(self) > imax) {
897 cpl_vector_set_size(xself, imax);
898 cpl_vector_set_size(yself, imax);
899 }
900
901 /* Propagate error, if any */
902 return cpl_error_set(cpl_func, error);
903}
904
905/*----------------------------------------------------------------------------*/
918/*----------------------------------------------------------------------------*/
919cpl_error_code
920irplib_polynomial_shift_1d_from_correlation(
921 cpl_polynomial *self,
922 const cpl_vector *obs,
923 irplib_base_spectrum_model *model,
924 cpl_error_code (*filler)(cpl_vector *,
925 const cpl_polynomial *,
926 irplib_base_spectrum_model *),
927 int hsize,
928 cpl_boolean doplot,
929 double *pxc)
930{
931 const int nobs = cpl_vector_get_size(obs);
932 const int nmodel = 2 * hsize + nobs;
933 cpl_vector *mspec1d;
934 cpl_vector *xcorr;
935 cpl_error_code error;
936 int ixc, xxc;
937 double xc;
938
939 cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT);
940 cpl_ensure_code(obs != NULL, CPL_ERROR_NULL_INPUT);
941 cpl_ensure_code(model != NULL, CPL_ERROR_NULL_INPUT);
942 cpl_ensure_code(filler != NULL, CPL_ERROR_NULL_INPUT);
943 cpl_ensure_code(hsize > 0, CPL_ERROR_ILLEGAL_INPUT);
944
945 /* Shift reference by -hsize so filler can be used without offset */
946 cpl_ensure_code(!cpl_polynomial_shift_1d(self, 0, -hsize),
947 cpl_error_get_code());
948
949 mspec1d = cpl_vector_new(nmodel);
950
951 if (filler(mspec1d, self, model)) {
952 cpl_vector_delete(mspec1d);
953 return cpl_error_set_where(cpl_func);
954 }
955
956 /* Should not be able to fail now */
957 xcorr = cpl_vector_new(1 + 2 * hsize);
958 ixc = cpl_vector_correlate(xcorr, mspec1d, obs);
959
960#ifdef IRPLIB_SPC_DUMP
961 /* Need irplib_wavecal.c rev. 1.12 through 1.15 */
962 irplib_polynomial_dump_corr_step(self, xcorr, "Shift");
963#endif
964
965 cpl_vector_delete(mspec1d);
966
967 error = cpl_polynomial_shift_1d(self, 0, (double)ixc);
968
969 xc = cpl_vector_get(xcorr, ixc);
970
971 xxc = ixc - hsize; /* The effect of the two shifts */
972
973 cpl_msg_info(cpl_func, "Shifting %d pixels (%g < %g)", xxc,
974 cpl_vector_get(xcorr, hsize), xc);
975
976 if (doplot) {
977 cpl_vector *xvals = cpl_vector_new(1 + 2 * hsize);
978 cpl_bivector *bcorr = cpl_bivector_wrap_vectors(xvals, xcorr);
979 int i;
980 double x = (double)-hsize;
981 char *title = cpl_sprintf("t 'Cross-correlation of shifted %d-pixel "
982 "spectrum (XCmax=%g at %d)' w linespoints",
983 nobs, cpl_vector_get(xcorr, ixc), xxc);
984
985 for (i = 0; i < 1 + 2 * hsize; i++, x += 1.0) {
986 cpl_vector_set(xvals, i, x);
987 }
988
989 cpl_plot_bivector("set grid;set xlabel 'Offset [pixel]';", title, "",
990 bcorr);
991 cpl_bivector_unwrap_vectors(bcorr);
992 cpl_vector_delete(xvals);
993 cpl_free(title);
994 }
995
996 cpl_vector_delete(xcorr);
997
998 cpl_ensure_code(!error, error);
999
1000 if (pxc != NULL)
1001 *pxc = xc;
1002
1003 return CPL_ERROR_NONE;
1004}
1005
1006
1007/*----------------------------------------------------------------------------*/
1028/*----------------------------------------------------------------------------*/
1029cpl_error_code
1030irplib_vector_fill_line_spectrum_model(cpl_vector *self,
1031 cpl_vector *linepix,
1032 cpl_vector *erftmp,
1033 const cpl_polynomial *disp,
1034 const cpl_bivector *lines,
1035 double wslit,
1036 double wfwhm,
1037 double xtrunc,
1038 int hsize,
1039 cpl_boolean dofast,
1040 cpl_boolean dolog,
1041 cpl_size *pulines)
1042{
1043 cpl_errorstate prestate;
1044 const double sigma = wfwhm * CPL_MATH_SIG_FWHM;
1045 const cpl_vector *xlines = cpl_bivector_get_x_const(lines);
1046 const double *dxlines = cpl_vector_get_data_const(xlines);
1047 const double *dylines = cpl_bivector_get_y_data_const(lines);
1048 double *plinepix = linepix ? cpl_vector_get_data(linepix) : NULL;
1049 const int nlines = cpl_vector_get_size(xlines);
1050 const int nself = cpl_vector_get_size(self);
1051 double *dself = cpl_vector_get_data(self);
1052 cpl_polynomial *dispi;
1053 double *profile = NULL;
1054 const cpl_size i0 = 0;
1055 const double p0 = cpl_polynomial_get_coeff(disp, &i0);
1056 double wl;
1057 double xpos = (double)(1 - hsize) - xtrunc;
1058 const double xmax = (double)(nself - hsize) + xtrunc;
1059 double xderiv, xextreme;
1060 cpl_error_code error = CPL_ERROR_NONE;
1061 int iline;
1062 cpl_size ulines = 0;
1063
1064 cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT);
1065 cpl_ensure_code(disp != NULL, CPL_ERROR_NULL_INPUT);
1066 cpl_ensure_code(lines != NULL, CPL_ERROR_NULL_INPUT);
1067
1068 cpl_ensure_code(wslit > 0.0, CPL_ERROR_ILLEGAL_INPUT);
1069 cpl_ensure_code(wfwhm > 0.0, CPL_ERROR_ILLEGAL_INPUT);
1070 cpl_ensure_code(hsize >= 0, CPL_ERROR_ILLEGAL_INPUT);
1071 cpl_ensure_code(xtrunc > 0.0, CPL_ERROR_ILLEGAL_INPUT);
1072 cpl_ensure_code(nself > 2 * hsize, CPL_ERROR_ILLEGAL_INPUT);
1073
1074 cpl_ensure_code(cpl_polynomial_get_dimension(disp) == 1,
1075 CPL_ERROR_ILLEGAL_INPUT);
1076 cpl_ensure_code(cpl_polynomial_get_degree(disp) > 0,
1077 CPL_ERROR_ILLEGAL_INPUT);
1078
1079 /* The smallest wavelength contributing to the spectrum. */
1080 wl = cpl_polynomial_eval_1d(disp, xpos, &xderiv);
1081
1082 if (wl <= 0.0)
1083 return cpl_error_set_message_macro(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
1084 __FILE__, __LINE__,
1085 "Non-positive wavelength at x=%g: "
1086 "P(x)=%g, P'(x)=%g",
1087 xpos, wl, xderiv);
1088
1089 if (xderiv <= 0.0)
1090 return cpl_error_set_message_macro(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
1091 __FILE__, __LINE__,
1092 "Non-increasing dispersion at "
1093 "x=%g: P'(x)=%g, P(x)=%g",
1094 xpos, xderiv, wl);
1095
1096 /* Find the 1st line */
1097 iline = cpl_vector_find(xlines, wl);
1098
1099 /* The first line must be at least at wl */
1100 if (dxlines[iline] < wl)
1101 iline++;
1102
1103 if (iline >= nlines)
1104 return cpl_error_set_message_macro(cpl_func, CPL_ERROR_DATA_NOT_FOUND,
1105 __FILE__, __LINE__,
1106 "The %d-line catalogue has only "
1107 "lines below P(%g)=%g > %g",
1108 nlines, xpos, wl,
1109 dxlines[nlines - 1]);
1110
1111 memset(dself, 0, nself * sizeof(double));
1112
1113 dispi = cpl_polynomial_duplicate(disp);
1114
1115 /* Verify monotony of dispersion */
1116 cpl_polynomial_derivative(dispi, 0);
1117
1118 prestate = cpl_errorstate_get();
1119
1120 if (cpl_polynomial_solve_1d(dispi, 0.5 * (nlines + 1), &xextreme, 1)) {
1121 cpl_errorstate_set(prestate);
1122 }
1123 else if (xpos < xextreme && xextreme < xmax) {
1124 cpl_polynomial_delete(dispi);
1125 return cpl_error_set_message_macro(
1126 cpl_func, CPL_ERROR_ILLEGAL_INPUT, __FILE__, __LINE__,
1127 "Non-monotone "
1128 "dispersion at x=%g: P'(x)=0, "
1129 "P(x)=%g",
1130 xextreme, cpl_polynomial_eval_1d(disp, xextreme, NULL));
1131 }
1132
1133 if (dofast) {
1134 const int npix = 1 + (int)xtrunc;
1135
1136 if (erftmp != NULL && cpl_vector_get_size(erftmp) == npix &&
1137 cpl_vector_get(erftmp, 0) > 0.0) {
1138 profile = cpl_vector_get_data(erftmp);
1139 }
1140 else {
1141 const double yval = 0.5 / wslit;
1142 const double x0p = 0.5 * wslit + 0.5;
1143 const double x0n = -0.5 * wslit + 0.5;
1144 double x1diff = irplib_erf_antideriv(x0p, sigma) -
1145 irplib_erf_antideriv(x0n, sigma);
1146 int ipix;
1147
1148 if (erftmp == NULL) {
1149 profile = (double *)cpl_malloc(sizeof(double) * (size_t)npix);
1150 }
1151 else {
1152 cpl_vector_set_size(erftmp, npix);
1153 profile = cpl_vector_get_data(erftmp);
1154 }
1155
1156 profile[0] = 2.0 * yval * x1diff;
1157
1158 for (ipix = 1; ipix < npix; ipix++) {
1159 const double x1 = (double)ipix;
1160 const double x1p = x1 + 0.5 * wslit + 0.5;
1161 const double x1n = x1 - 0.5 * wslit + 0.5;
1162 const double x0diff = x1diff;
1163
1164 x1diff = irplib_erf_antideriv(x1p, sigma) -
1165 irplib_erf_antideriv(x1n, sigma);
1166
1167 profile[ipix] = yval * (x1diff - x0diff);
1168 }
1169 }
1170 }
1171
1172 cpl_polynomial_copy(dispi, disp);
1173
1174 /* FIXME: A custom version of cpl_polynomial_solve_1d() which returns
1175 P'(xpos) can be used for the 1st NR-iteration. */
1176 /* Further, the sign of P'(xpos) could be checked for all lines. */
1177 /* Perform 1st NR-iteration in solving for P(xpos) = dxlines[iline] */
1178 xpos -= (wl - dxlines[iline]) / xderiv;
1179
1180 /* Iterate through the lines */
1181 for (; !error && iline < nlines; iline++) {
1182 /* Lines may have a non-physical intensity (e.g. zero) to indicate some
1183 property of the line, e.g. unknown intensity due to blending */
1184 if (dylines[iline] <= 0.0)
1185 continue;
1186
1187 /* Use 1st guess, if available (Use 0.0 to flag unavailable) */
1188 if (plinepix != NULL && plinepix[iline] > 0.0)
1189 xpos = plinepix[iline];
1190
1191 if (xpos > xmax)
1192 xpos = xmax; /* FIXME: Better to limit xpos ? */
1193
1194 /* Find the (sub-) pixel position of the line */
1195 error = cpl_polynomial_set_coeff(dispi, &i0, p0 - dxlines[iline]) ||
1196 cpl_polynomial_solve_1d(dispi, xpos, &xpos, 1);
1197
1198 if (xpos > xmax) {
1199 if (error) {
1200 error = 0;
1201 cpl_msg_debug(cpl_func,
1202 "Stopping spectrum fill at line %d/%d "
1203 "at xpos=%g > xmax=%g",
1204 iline, nlines, xpos, xmax);
1205 cpl_errorstate_dump(prestate, CPL_FALSE,
1206 irplib_errorstate_dump_debug);
1207 cpl_errorstate_set(prestate);
1208 }
1209 break;
1210 }
1211 else if (error) {
1212 if (linepix != NULL && ulines)
1213 (void)cpl_vector_fill(linepix, 0.0);
1214 (void)cpl_error_set_message_macro(cpl_func, cpl_error_get_code(),
1215 __FILE__, __LINE__,
1216 "Could not find pixel-position "
1217 "of line %d/%d at wavelength=%g."
1218 " xpos=%g, xmax=%g",
1219 iline, nlines, dxlines[iline],
1220 xpos, xmax);
1221 break;
1222 }
1223 else if (dofast) {
1224 const double frac = fabs(xpos - floor(xpos));
1225#ifdef IRPLIB_WAVECAL_FAST_FAST
1226 const double frac0 = 1.0 - frac; /* Weight opposite of distance */
1227#else
1228 /* Center intensity correctly */
1229 const double ep1pw =
1230 irplib_erf_antideriv(frac + 0.5 * wslit, sigma);
1231 const double en1pw =
1232 irplib_erf_antideriv(frac + 0.5 * wslit - 1.0, sigma);
1233 const double ep1nw =
1234 irplib_erf_antideriv(frac - 0.5 * wslit, sigma);
1235 const double en1nw =
1236 irplib_erf_antideriv(frac - 0.5 * wslit - 1.0, sigma);
1237 const double frac0 =
1238 (en1nw - en1pw) / (ep1pw - en1pw - ep1nw + en1nw);
1239
1240#endif
1241 const double frac1 = 1.0 - frac0;
1242 const double yval0 = frac0 * dylines[iline];
1243 const double yval1 = frac1 * dylines[iline];
1244 const int npix = 1 + (int)xtrunc;
1245 int ipix;
1246 int i0n = hsize - 1 + floor(xpos);
1247 int i0p = i0n;
1248 int i1n = i0n + 1;
1249 int i1p = i1n;
1250 cpl_boolean didline = CPL_FALSE;
1251
1252
1253 /* Update 1st guess for next time, if available */
1254 if (plinepix != NULL)
1255 plinepix[iline] = xpos;
1256
1257 if (frac0 < 0.0) {
1258 (void)cpl_error_set_message_macro(cpl_func,
1259 CPL_ERROR_UNSPECIFIED,
1260 __FILE__, __LINE__,
1261 "Illegal split at x=%g: %g + "
1262 "%g = 1",
1263 xpos, frac0, frac1);
1264#ifdef IRPLIB_WAVEVAL_DEBUG
1265 }
1266 else {
1267 cpl_msg_warning(cpl_func, "profile split at x=%g: %g + %g = 1",
1268 xpos, frac0, frac1);
1269#endif
1270 }
1271
1272 for (ipix = 0; ipix < npix; ipix++, i0n--, i0p++, i1n--, i1p++) {
1273 if (i0n >= 0 && i0n < nself) {
1274 dself[i0n] += yval0 * profile[ipix];
1275 didline = CPL_TRUE;
1276 }
1277 if (i1n >= 0 && i1n < nself && ipix + 1 < npix) {
1278 dself[i1n] += yval1 * profile[ipix + 1];
1279 didline = CPL_TRUE;
1280 }
1281
1282 if (ipix == 0)
1283 continue;
1284
1285 if (i0p >= 0 && i0p < nself) {
1286 dself[i0p] += yval0 * profile[ipix];
1287 didline = CPL_TRUE;
1288 }
1289 if (i1p >= 0 && i1p < nself && ipix + 1 < npix) {
1290 dself[i1p] += yval1 * profile[ipix + 1];
1291 didline = CPL_TRUE;
1292 }
1293 }
1294
1295 if (didline)
1296 ulines++;
1297 }
1298 else {
1299 const double yval = 0.5 * dylines[iline] / wslit;
1300 const int ifirst =
1301 IRPLIB_MAX((int)(xpos - xtrunc + 0.5), 1 - hsize);
1302 const int ilast = IRPLIB_MIN((int)(xpos + xtrunc), nself - hsize);
1303 int ipix;
1304 const double x0 = (double)ifirst - xpos;
1305 const double x0p = x0 + 0.5 * wslit - 0.5;
1306 const double x0n = x0 - 0.5 * wslit - 0.5;
1307 double x1diff = irplib_erf_antideriv(x0p, sigma) -
1308 irplib_erf_antideriv(x0n, sigma);
1309
1310 /* Update 1st guess for next time, if available */
1311 if (plinepix != NULL)
1312 plinepix[iline] = xpos;
1313
1314 if (ilast >= ifirst)
1315 ulines++;
1316
1317 for (ipix = ifirst; ipix <= ilast; ipix++) {
1318 const double x1 = (double)ipix - xpos;
1319 const double x1p = x1 + 0.5 * wslit + 0.5;
1320 const double x1n = x1 - 0.5 * wslit + 0.5;
1321 const double x0diff = x1diff;
1322
1323 x1diff = irplib_erf_antideriv(x1p, sigma) -
1324 irplib_erf_antideriv(x1n, sigma);
1325
1326 dself[ipix + hsize - 1] += yval * (x1diff - x0diff);
1327 }
1328 }
1329 }
1330
1331 cpl_polynomial_delete(dispi);
1332 if (erftmp == NULL)
1333 cpl_free(profile);
1334
1335 cpl_ensure_code(!error, cpl_error_get_code());
1336
1337 if (dolog) {
1338 int i;
1339 for (i = 0; i < nself; i++) {
1340 dself[i] = dself[i] > 0.0 ? log(1.0 + dself[i]) : 0.0;
1341 }
1342 }
1343
1344 if (!ulines)
1345 return cpl_error_set_message_macro(
1346 cpl_func, CPL_ERROR_DATA_NOT_FOUND, __FILE__, __LINE__,
1347 "The %d-line "
1348 "catalogue has no lines in the range "
1349 "%g -> P(%g)=%g",
1350 nlines, wl, xmax, cpl_polynomial_eval_1d(disp, xmax, NULL));
1351
1352 if (pulines != NULL)
1353 *pulines = ulines;
1354
1355 return CPL_ERROR_NONE;
1356}
1357
1358/*----------------------------------------------------------------------------*/
1367/*----------------------------------------------------------------------------*/
1368inline double
1369irplib_erf_antideriv(double x, double sigma)
1370{
1371 return x * erf(x / (sigma * CPL_MATH_SQRT2)) +
1372 2.0 * sigma / CPL_MATH_SQRT2PI * exp(-0.5 * x * x / (sigma * sigma));
1373}
1374
1375
1376#ifdef HAVE_GSL
1377
1378/*----------------------------------------------------------------------------*/
1385/*----------------------------------------------------------------------------*/
1386static double
1387irplib_gsl_correlation(const gsl_vector *self, void *data)
1388{
1389 irplib_multimin *mindata = (irplib_multimin *)data;
1390 cpl_errorstate prestate = cpl_errorstate_get();
1391 int nobs, nmodel, ndiff;
1392 cpl_size i;
1393
1394 cpl_ensure(self != NULL, CPL_ERROR_NULL_INPUT, GSL_NAN);
1395 cpl_ensure(data != NULL, CPL_ERROR_NULL_INPUT, GSL_NAN);
1396
1397 cpl_ensure(mindata->filler != NULL, CPL_ERROR_NULL_INPUT, GSL_NAN);
1398 cpl_ensure(mindata->observed != NULL, CPL_ERROR_NULL_INPUT, GSL_NAN);
1399 cpl_ensure(mindata->spectrum != NULL, CPL_ERROR_NULL_INPUT, GSL_NAN);
1400
1401 nobs = cpl_vector_get_size(mindata->observed);
1402 nmodel = cpl_vector_get_size(mindata->spectrum);
1403 ndiff = nmodel - nobs;
1404
1405 cpl_ensure((ndiff & 1) == 0, CPL_ERROR_ILLEGAL_INPUT, GSL_NAN);
1406
1407 cpl_ensure(cpl_vector_get_size(mindata->vxc) == 1 + ndiff,
1408 CPL_ERROR_ILLEGAL_INPUT, GSL_NAN);
1409
1410 ndiff /= 2;
1411
1412 for (i = 0; i < (cpl_size)self->size; i++) {
1413 const double value = gsl_vector_get(self, (size_t)i);
1414 cpl_polynomial_set_coeff(mindata->disp1d, &i, value);
1415 }
1416
1417 /* Shift reference by -ndiff so filler can be used without offset.
1418 The subsequent polynomial shift is reduced by -ndiff. */
1419 cpl_ensure_code(!cpl_polynomial_shift_1d(mindata->disp1d, 0, -ndiff),
1420 cpl_error_get_code());
1421
1422 if (mindata->filler(mindata->spectrum, mindata->disp1d, mindata->param) ||
1423 !cpl_errorstate_is_equal(prestate)) {
1424 /* The fill failed. Ensure the discarding of this candidate by
1425 setting the cross-correlation to its minimum possible value. */
1426
1427 (void)cpl_vector_fill(mindata->vxc, -1.0);
1428
1429 mindata->maxxc = ndiff;
1430
1431 if (!cpl_errorstate_is_equal(prestate)) {
1432 cpl_msg_debug(cpl_func, "Spectrum fill failed:");
1433 cpl_errorstate_dump(prestate, CPL_FALSE,
1434 irplib_errorstate_dump_debug);
1435 cpl_errorstate_set(prestate);
1436 }
1437 }
1438 else {
1439 mindata->maxxc = cpl_vector_correlate(mindata->vxc, mindata->spectrum,
1440 mindata->observed);
1441 }
1442
1443#ifdef IRPLIB_SPC_DUMP
1444 /* Need irplib_wavecal.c rev. 1.12 through 1.15 */
1445 irplib_polynomial_dump_corr_step(mindata->disp1d, mindata->vxc, "Optimize");
1446#endif
1447
1448 mindata->xc = cpl_vector_get(mindata->vxc, ndiff);
1449
1450 if (mindata->maxxc != ndiff &&
1451 cpl_vector_get(mindata->vxc, mindata->maxxc) > mindata->mxc) {
1452 const irplib_base_spectrum_model *arclamp =
1453 (const irplib_base_spectrum_model *)mindata->param;
1454
1455 if (mindata->mdisp == NULL) {
1456 mindata->mdisp = cpl_polynomial_duplicate(mindata->disp1d);
1457 }
1458 else {
1459 cpl_polynomial_copy(mindata->mdisp, mindata->disp1d);
1460 }
1461 mindata->mxc = cpl_vector_get(mindata->vxc, mindata->maxxc);
1462 mindata->ishift = mindata->maxxc; /* Offset -ndiff pre-shifted above */
1463 cpl_msg_debug(cpl_func,
1464 "Local maximum: %g(%d) > %g(%d) (cost=%u:%u. "
1465 "lines=%u)",
1466 mindata->mxc, mindata->maxxc, mindata->xc, ndiff,
1467 (unsigned)arclamp->cost, (unsigned)arclamp->xcost,
1468 (unsigned)arclamp->ulines);
1469 }
1470
1471 return -mindata->xc;
1472}
1473
1474#endif
1475
1476/*----------------------------------------------------------------------------*/
1499/*----------------------------------------------------------------------------*/
1500cpl_error_code
1501irplib_polynomial_find_1d_from_correlation_all(
1502 cpl_polynomial *self,
1503 int maxdeg,
1504 const cpl_vector *obs,
1505 int nmaxima,
1506 int linelim,
1507 irplib_base_spectrum_model *model,
1508 cpl_error_code (*filler)(cpl_vector *,
1509 const cpl_polynomial *,
1510 irplib_base_spectrum_model *),
1511 double pixtol,
1512 double pixstep,
1513 int hsize,
1514 int maxite,
1515 int maxfail,
1516 int maxcont,
1517 cpl_boolean doplot,
1518 double *pxc)
1519{
1520#ifdef HAVE_GSL
1521
1522 cpl_errorstate prestate = cpl_errorstate_get();
1523 cpl_polynomial *start;
1524 cpl_polynomial *cand;
1525 cpl_polynomial *backup;
1526 cpl_error_code error = CPL_ERROR_NONE;
1527 double xc;
1528 cpl_bivector *xtshift = cpl_bivector_new(nmaxima ? nmaxima : 1);
1529 const cpl_vector *xtshiftx = cpl_bivector_get_x_const(xtshift);
1530 const cpl_vector *xtshifty = cpl_bivector_get_y_const(xtshift);
1531 int nshift;
1532 int imaximum = -1;
1533 int imaxima;
1534
1535#endif
1536
1537 cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT);
1538 cpl_ensure_code(obs != NULL, CPL_ERROR_NULL_INPUT);
1539 cpl_ensure_code(model != NULL, CPL_ERROR_NULL_INPUT);
1540 cpl_ensure_code(filler != NULL, CPL_ERROR_NULL_INPUT);
1541 cpl_ensure_code(pxc != NULL, CPL_ERROR_NULL_INPUT);
1542
1543 cpl_ensure_code(cpl_polynomial_get_dimension(self) == 1,
1544 CPL_ERROR_ILLEGAL_INPUT);
1545
1546 cpl_ensure_code(cpl_polynomial_get_degree(self) > 0,
1547 CPL_ERROR_ILLEGAL_INPUT);
1548
1549 cpl_ensure_code(maxdeg >= 0, CPL_ERROR_ILLEGAL_INPUT);
1550 cpl_ensure_code(pixtol > 0.0, CPL_ERROR_ILLEGAL_INPUT);
1551 cpl_ensure_code(pixstep > 0.0, CPL_ERROR_ILLEGAL_INPUT);
1552 cpl_ensure_code(hsize >= 0, CPL_ERROR_ILLEGAL_INPUT);
1553 cpl_ensure_code(maxite >= 0, CPL_ERROR_ILLEGAL_INPUT);
1554 cpl_ensure_code(nmaxima >= 0, CPL_ERROR_ILLEGAL_INPUT);
1555 cpl_ensure_code(maxfail > 0, CPL_ERROR_ILLEGAL_INPUT);
1556 cpl_ensure_code(maxcont > 0, CPL_ERROR_ILLEGAL_INPUT);
1557 cpl_ensure_code(linelim >= 0, CPL_ERROR_ILLEGAL_INPUT);
1558
1559#ifndef HAVE_GSL
1560 /* Avoid unused variable warning */
1561 cpl_ensure_code(doplot == CPL_TRUE || doplot == CPL_FALSE,
1562 CPL_ERROR_ILLEGAL_INPUT);
1563 return cpl_error_set_message(cpl_func, CPL_ERROR_UNSUPPORTED_MODE,
1564 "GSL is not available");
1565#else
1566
1567 if (irplib_bivector_find_shift_from_correlation(xtshift, self, obs, model,
1568 filler, hsize, doplot,
1569 &xc)) {
1570 cpl_bivector_delete(xtshift);
1571 return cpl_error_set_where(cpl_func);
1572 }
1573
1574 if (model->ulines > (cpl_size)linelim) {
1575 /* The initial, optimal (integer) shift */
1576 const double xxc = cpl_vector_get(xtshiftx, 0);
1577 const double xc0 = cpl_vector_get(xtshifty, 0);
1578
1579 cpl_msg_warning(cpl_func,
1580 "Doing only shift=%g pixels with lines=%u > "
1581 "%d and XC=%g",
1582 xxc, (unsigned)model->ulines, linelim, xc0);
1583
1584 cpl_polynomial_shift_1d(self, 0, xxc);
1585
1586 *pxc = xc0;
1587
1588 cpl_bivector_delete(xtshift);
1589
1590 return CPL_ERROR_NONE;
1591 }
1592
1593 start = cpl_polynomial_duplicate(self);
1594 cand = cpl_polynomial_new(1);
1595 backup = cpl_polynomial_new(1);
1596
1597 /* Number of (local) maxima to use as starting point of the optimization */
1598 nshift = cpl_bivector_get_size(xtshift);
1599 if (nmaxima == 0 || nmaxima > nshift)
1600 nmaxima = nshift;
1601
1602 cpl_msg_info(cpl_func,
1603 "Optimizing %d/%d local shift-maxima "
1604 "(no-shift xc=%g. linelim=%d)",
1605 nmaxima, nshift, xc, linelim);
1606 if (cpl_msg_get_level() <= CPL_MSG_DEBUG)
1607 cpl_bivector_dump(xtshift, stdout);
1608
1609 for (imaxima = 0; imaxima < nmaxima; imaxima++) {
1610 /* The initial, optimal (integer) shift */
1611 const double xxc = cpl_vector_get(xtshiftx, imaxima);
1612 double xtpixstep = pixstep;
1613 double xtpixtol = pixtol;
1614 double xtxc;
1615 cpl_boolean ok = CPL_FALSE;
1616 int nfail;
1617
1618
1619 cpl_polynomial_copy(cand, start);
1620 cpl_polynomial_shift_1d(cand, 0, xxc);
1621 cpl_polynomial_copy(backup, cand);
1622
1623 /* Increase tolerance until convergence */
1624 for (nfail = 0; nfail < maxfail;
1625 nfail++, xtpixtol *= 2.0, xtpixstep *= 2.0) {
1626 int restart = maxcont;
1627 cpl_boolean redo;
1628
1629 do {
1630 if (error) {
1631 cpl_errorstate_dump(prestate, CPL_FALSE,
1632 irplib_errorstate_dump_debug);
1633 cpl_errorstate_set(prestate);
1634 }
1635 error = irplib_polynomial_find_1d_from_correlation_(
1636 cand, maxdeg, obs, model, filler, xtpixtol, xtpixstep, 2,
1637 maxite, &xtxc, &redo);
1638 if (redo && !error)
1639 error = CPL_ERROR_CONTINUE;
1640 } while (((!error && redo) || error == CPL_ERROR_CONTINUE) &&
1641 --restart);
1642
1643 if (!error && !redo) {
1644 cpl_msg_debug(cpl_func,
1645 "XC(imax=%d/%d:xtpixtol=%g): %g "
1646 "(cost=%u:%u)",
1647 1 + imaxima, nmaxima, xtpixtol, xtxc,
1648 (unsigned)model->cost, (unsigned)model->xcost);
1649 break;
1650 }
1651 cpl_msg_warning(cpl_func,
1652 "Increasing xtpixtol from %g (%g, imax="
1653 "%d/%d)",
1654 xtpixtol, xtpixstep, 1 + imaxima, nmaxima);
1655 if (model->ulines > (cpl_size)linelim) {
1656 cpl_msg_warning(cpl_func,
1657 "Stopping search-refinement via "
1658 "catalogue with %u lines > %d",
1659 (unsigned)model->ulines, linelim);
1660 break;
1661 }
1662 cpl_polynomial_copy(cand, start);
1663 }
1664
1665 /* Decrease tolerance until divergence, keep previous */
1666 for (; !error && xtpixtol > 0.0; xtpixtol *= 0.25, xtpixstep *= 0.5) {
1667 int restart = maxcont;
1668 cpl_boolean redo;
1669
1670 cpl_polynomial_copy(backup, cand);
1671 do {
1672 if (error) {
1673 cpl_errorstate_dump(prestate, CPL_FALSE,
1674 irplib_errorstate_dump_debug);
1675 cpl_errorstate_set(prestate);
1676 }
1677 error = irplib_polynomial_find_1d_from_correlation_(
1678 cand, maxdeg, obs, model, filler, xtpixtol, xtpixstep, 2,
1679 maxite, &xtxc, &redo);
1680 if (redo && !error)
1681 error = CPL_ERROR_CONTINUE;
1682 } while (((!error && redo) || error == CPL_ERROR_CONTINUE) &&
1683 --restart);
1684 if (error)
1685 break;
1686 ok = CPL_TRUE;
1687 if (redo)
1688 break;
1689 cpl_msg_debug(cpl_func,
1690 "XC(imax=%d/%d:xtpixtol=%g): %g (cost=%u:%u"
1691 ". ulines=%u)",
1692 1 + imaxima, nmaxima, xtpixtol, xtxc,
1693 (unsigned)model->cost, (unsigned)model->xcost,
1694 (unsigned)model->ulines);
1695 if (model->ulines > (cpl_size)linelim) {
1696 cpl_msg_info(cpl_func,
1697 "Stopping search-refinement via "
1698 "catalogue with %u lines > %u",
1699 (unsigned)model->ulines, linelim);
1700 break;
1701 }
1702 }
1703
1704 if (error) {
1705 error = 0;
1706 cpl_errorstate_dump(prestate, CPL_FALSE,
1707 irplib_errorstate_dump_debug);
1708 cpl_errorstate_set(prestate);
1709 cpl_polynomial_copy(cand, backup);
1710 }
1711 if (ok && xtxc > xc) {
1712 imaximum = imaxima;
1713 cpl_polynomial_copy(self, cand);
1714 xc = xtxc;
1715
1716 cpl_msg_info(cpl_func,
1717 "XC(imax=%d/%d): %g -> %g (initial-shift=%g. "
1718 "cost=%u:%u. lines=%u)",
1719 1 + imaxima, nmaxima,
1720 cpl_vector_get(xtshifty, imaxima), xtxc,
1721 cpl_vector_get(xtshiftx, imaxima),
1722 (unsigned)model->cost, (unsigned)model->xcost,
1723 (unsigned)model->ulines);
1724 }
1725 else {
1726 cpl_msg_info(cpl_func,
1727 "xc(imax=%d/%d): %g -> %g (initial-shift=%g. "
1728 "cost=%u:%u. lines=%u)",
1729 1 + imaxima, nmaxima,
1730 cpl_vector_get(xtshifty, imaxima), xtxc,
1731 cpl_vector_get(xtshiftx, imaxima),
1732 (unsigned)model->cost, (unsigned)model->xcost,
1733 (unsigned)model->ulines);
1734 }
1735 }
1736
1737 cpl_polynomial_delete(start);
1738 cpl_polynomial_delete(backup);
1739 cpl_polynomial_delete(cand);
1740
1741 if (imaximum < 0) {
1742 /* The initial, optimal (integer) shift */
1743 const double xxc = cpl_vector_get(xtshiftx, 0);
1744 const double xc0 = cpl_vector_get(xtshifty, 0);
1745
1746 error = cpl_error_set_message(cpl_func, CPL_ERROR_DATA_NOT_FOUND,
1747 "Could not improve XC=%g over %d "
1748 "local shift-maxima, best at shift %g",
1749 xc0, nmaxima, xxc);
1750 }
1751 else {
1752 cpl_msg_info(cpl_func,
1753 "Maximal XC=%g (up from %g, with initial pixel-"
1754 "shift of %g) at %d/%d local shift-maximi",
1755 xc, cpl_vector_get(xtshifty, imaximum),
1756 cpl_vector_get(xtshiftx, imaximum), 1 + imaximum, nmaxima);
1757
1758 if (doplot) {
1759 irplib_plot_spectrum_and_model(obs, self, model, filler);
1760 }
1761
1762 *pxc = xc;
1763 }
1764
1765 cpl_bivector_delete(xtshift);
1766
1767 return error;
1768
1769#endif
1770}
1771
irplib_errorstate_dump_debug
void irplib_errorstate_dump_debug(unsigned self, unsigned first, unsigned last)
Dump a single CPL error at the CPL debug level.
Definition irplib_utils.c:156
irplib_bivector_count_positive
int irplib_bivector_count_positive(const cpl_bivector *self, double x_min, double x_max)
Count the positive Y-entries in a given X-range.
Definition irplib_wavecal.c:120
irplib_vector_fill_line_spectrum
cpl_error_code irplib_vector_fill_line_spectrum(cpl_vector *self, const cpl_polynomial *disp, irplib_base_spectrum_model *lsslamp)
Generate a 1D spectrum from a model and a dispersion relation.
Definition irplib_wavecal.c:495
irplib_bivector_find_shift_from_correlation
cpl_error_code irplib_bivector_find_shift_from_correlation(cpl_bivector *self, const cpl_polynomial *disp, const cpl_vector *obs, irplib_base_spectrum_model *model, cpl_error_code(*filler)(cpl_vector *, const cpl_polynomial *, irplib_base_spectrum_model *), int hsize, cpl_boolean doplot, double *pxc)
Find shift(s) that maximizes (locally) the cross-correlation.
Definition irplib_wavecal.c:742
irplib_erf_antideriv
double irplib_erf_antideriv(double x, double sigma)
The antiderivative of erx(x/sigma/sqrt(2)) with respect to x.
Definition irplib_wavecal.c:1369
irplib_polynomial_fit_2d_dispersion
cpl_error_code irplib_polynomial_fit_2d_dispersion(cpl_polynomial *self, const cpl_image *imgwave, int fitdeg, double *presid)
Fit a 2D-dispersion from an image of wavelengths.
Definition irplib_wavecal.c:155
irplib_vector_fill_line_spectrum_fast
cpl_error_code irplib_vector_fill_line_spectrum_fast(cpl_vector *self, const cpl_polynomial *disp, irplib_base_spectrum_model *lsslamp)
Generate a 1D spectrum from a model and a dispersion relation.
Definition irplib_wavecal.c:576
irplib_vector_fill_logline_spectrum_fast
cpl_error_code irplib_vector_fill_logline_spectrum_fast(cpl_vector *self, const cpl_polynomial *disp, irplib_base_spectrum_model *lsslamp)
Generate a 1D spectrum from a model and a dispersion relation.
Definition irplib_wavecal.c:616
irplib_plot_spectrum_and_model
cpl_error_code irplib_plot_spectrum_and_model(const cpl_vector *self, const cpl_polynomial *disp1d, irplib_base_spectrum_model *model, cpl_error_code(*filler)(cpl_vector *, const cpl_polynomial *, irplib_base_spectrum_model *))
Plot a 1D spectrum and one from a model.
Definition irplib_wavecal.c:654
irplib_polynomial_find_1d_from_correlation_all
cpl_error_code irplib_polynomial_find_1d_from_correlation_all(cpl_polynomial *self, int maxdeg, const cpl_vector *obs, int nmaxima, int linelim, irplib_base_spectrum_model *model, cpl_error_code(*filler)(cpl_vector *, const cpl_polynomial *, irplib_base_spectrum_model *), double pixtol, double pixstep, int hsize, int maxite, int maxfail, int maxcont, cpl_boolean doplot, double *pxc)
Modify self by maximizing the cross-correlation across all maxima.
Definition irplib_wavecal.c:1501
irplib_vector_fill_line_spectrum_model
cpl_error_code irplib_vector_fill_line_spectrum_model(cpl_vector *self, cpl_vector *linepix, cpl_vector *erftmp, const cpl_polynomial *disp, const cpl_bivector *lines, double wslit, double wfwhm, double xtrunc, int hsize, cpl_boolean dofast, cpl_boolean dolog, cpl_size *pulines)
Generate a 1D spectrum from (arc) lines and a dispersion relation.
Definition irplib_wavecal.c:1030
irplib_polynomial_find_1d_from_correlation
cpl_error_code irplib_polynomial_find_1d_from_correlation(cpl_polynomial *self, int maxdeg, const cpl_vector *obs, irplib_base_spectrum_model *model, cpl_error_code(*filler)(cpl_vector *, const cpl_polynomial *, irplib_base_spectrum_model *), double pixtol, double pixstep, int hsize, int maxite, double *pxc)
Modify self by maximizing the cross-correlation.
Definition irplib_wavecal.c:242
irplib_vector_fill_logline_spectrum
cpl_error_code irplib_vector_fill_logline_spectrum(cpl_vector *self, const cpl_polynomial *disp, irplib_base_spectrum_model *lsslamp)
Generate a 1D spectrum from a model and a dispersion relation.
Definition irplib_wavecal.c:535
irplib_polynomial_shift_1d_from_correlation
cpl_error_code irplib_polynomial_shift_1d_from_correlation(cpl_polynomial *self, const cpl_vector *obs, irplib_base_spectrum_model *model, cpl_error_code(*filler)(cpl_vector *, const cpl_polynomial *, irplib_base_spectrum_model *), int hsize, cpl_boolean doplot, double *pxc)
Shift self by the amount that maximizes the cross-correlation.
Definition irplib_wavecal.c:920
Generated on Tue Oct 13 2026 12:27:26 for CR2RE Pipeline Reference Manual by doxygen 1.13.2