CR2RE Pipeline Reference Manual 1.7.0
cr2res_extract.c
1/*
2 * This file is part of the CR2RES Pipeline
3 * Copyright (C) 2002,2003 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/* Removing this pragma to give significant speedup to cal_flat
21 * Don't think its actually been needed since the original development of
22 * curved extraction. It also isn't recognized by clang and some
23 * versions of gcc.
24 * pragma GCC optimize("O0")
25*/
26
27#ifdef HAVE_CONFIG_H
28#include <config.h>
29#endif
30
31/*-----------------------------------------------------------------------------
32 Includes
33 -----------------------------------------------------------------------------*/
34#include <limits.h>
35#include <math.h>
36#include <cpl.h>
37
38#include "irplib_utils.h"
39
40#include "cr2res_dfs.h"
41#include "cr2res_trace.h"
42#include "cr2res_extract.h"
43#include "cr2res_wave.h"
44#include "cr2res_bpm.h"
45
46/*-----------------------------------------------------------------------------
47 Defines
48 -----------------------------------------------------------------------------*/
49
50typedef unsigned char byte;
51#define min(a,b) (((a)<(b))?(a):(b))
52#define max(a,b) (((a)>(b))?(a):(b))
53#define signum(a) (((a)>0)?1:((a)<0)?-1:0)
54/* The zeta entries of one detector pixel are contiguous in memory, as
55 required by the pixel-centric SLE fill loops */
56#define zeta_index(x, y, z) \
57 ((z) + (y) * (3 * (osample + 1)) + (x) * (3 * (osample + 1)) * nrows)
58#define mzeta_index(x, y) ((y) + (x) * nrows)
59/* Band matrices are stored row-major (band entries for one row are
60 contiguous): this matches the access pattern of both the SLE fill
61 loops and the row-major bandsol */
62#define laij_index(x, y) ((x) * (4 * osample + 1) + (y))
63#define paij_index(x, y) ((x) * nx + (y))
64/* 6 polynomial coefficients per column (degrees up to 5 supported) */
65#define curve_index(x, y) ((x) * 6 + (y))
66
67typedef struct {
68 int x ;
69 int iy ; /* Contributing subpixel x,iy */
70 double w; /* Contribution weight <= 1/osample */
71} zeta_ref;
72
73/* Key ranges of one pixel's zeta list, maintained by the zeta build.
74 The SLE fill loops merge the list by subpixel index iy (sL system) or
75 source column x (sP system) into a small dense window [min, max]
76 instead of searching a list of unique keys. */
77typedef struct {
78 int min_iy, max_iy ;
79 int min_x, max_x ;
80} zeta_rng;
81
82/*-----------------------------------------------------------------------------
83 Functions prototypes
84 -----------------------------------------------------------------------------*/
85
86static int cr2res_extract_slit_func_curved(
87 double error_factor,
88 int ncols,
89 int nrows,
90 int osample,
91 double pclip,
92 double * im,
93 double * pix_unc,
94 int * mask,
95 double * ycen,
96 int * ycen_offset,
97 int y_lower_lim,
98 const double * slitcurve,
99 int delta_x,
100 double * sL,
101 double * sP,
102 double * model,
103 double * unc,
104 double lambda_sP,
105 double lambda_sL,
106 double sP_stop,
107 int maxiter,
108 double kappa,
109 const double * slit_func_in,
110 double * sP_old,
111 double * l_Aij,
112 double * p_Aij,
113 double * l_bj,
114 double * p_bj,
115 double * zw,
116 int * zk,
117 zeta_ref * zeta,
118 int * m_zeta,
119 zeta_rng * z_rng) ;
120
121static int cr2res_extract_zeta_tensors(
122 int ncols,
123 int nrows,
124 const double * ycen,
125 const int * ycen_offset,
126 int y_lower_lim,
127 int osample,
128 const double * slitcurve,
129 zeta_ref * zeta,
130 int * m_zeta,
131 zeta_rng * z_rng) ;
132
133static int cr2res_extract_bandsol_rowmajor(
134 double * a,
135 double * r,
136 int n,
137 int nd) ;
138
139static int cr2res_extract_slitdec_adjust_swath(
140 cpl_vector * ycen,
141 int height,
142 int leny,
143 int sw,
144 int lenx,
145 int dx,
146 cpl_vector ** bins_begin,
147 cpl_vector ** bins_end) ;
148
149/*----------------------------------------------------------------------------*/
153/*----------------------------------------------------------------------------*/
154
156
157/*----------------------------------------------------------------------------*/
185/*----------------------------------------------------------------------------*/
186int cr2res_extract_traces(
187 const hdrl_image * img,
188 const cpl_table * traces,
189 const cpl_table * slit_func_in,
190 const cpl_table * blaze_table_in,
191 float blaze_norm,
192 int reduce_order,
193 int reduce_trace,
194 cr2res_extr_method extr_method,
195 int extr_height,
196 int swath_width,
197 int oversample,
198 double pclip,
199 double smooth_slit,
200 double smooth_spec,
201 int niter,
202 double kappa,
203 double error_factor,
204 int display,
205 int disp_order_idx,
206 int disp_trace,
207 cpl_table ** extracted,
208 cpl_table ** slit_func,
209 hdrl_image ** model_master)
210{
211 cpl_bivector ** spectrum ;
212 cpl_vector ** slit_func_vec ;
213 hdrl_image ** model_one ;
214 cpl_table * slit_func_loc ;
215 cpl_table * extract_loc ;
216 hdrl_image * model_loc ;
217 int nb_traces, i;
218
219 /* Check Entries */
220 if (img == NULL || traces == NULL) return -1 ;
221
222 /* Initialise */
223 nb_traces = cpl_table_get_nrow(traces) ;
224
225 /* Allocate Data containers */
226 spectrum = cpl_malloc(nb_traces * sizeof(cpl_bivector *)) ;
227 slit_func_vec = cpl_malloc(nb_traces * sizeof(cpl_vector *)) ;
228 model_one = cpl_malloc(nb_traces * sizeof(hdrl_image *)) ;
229 model_loc = hdrl_image_duplicate(img) ;
230 hdrl_image_mul_scalar(model_loc, (hdrl_value){0.0, 0.0}) ;
231
232 /* Loop over the traces and extract them. The traces are independent
233 of each other: with OpenMP they are extracted in parallel, each
234 trace entirely within one thread (so per-trace results do not
235 depend on the threading), and the per-trace models are merged
236 sequentially in trace order below. */
237#ifdef _OPENMP
238#pragma omp parallel for schedule(dynamic, 1)
239#endif
240 for (i=0 ; i<nb_traces ; i++) {
241 /* Initialise */
242 cpl_vector * slit_func_in_vec ;
243 hdrl_image * model_loc_one ;
244 int trace_id;
245 int order;
246 slit_func_vec[i] = NULL ;
247 spectrum[i] = NULL ;
248 model_loc_one = NULL ;
249 model_one[i] = NULL ;
250
251 /* Get Order and trace id */
252 order = cpl_table_get(traces, CR2RES_COL_ORDER, i, NULL) ;
253 trace_id = cpl_table_get(traces, CR2RES_COL_TRACENB, i, NULL) ;
254
255 /* Check if this order needs to be skipped */
256 if (reduce_order > -1 && order != reduce_order) continue ;
257
258 /* Check if this trace needs to be skipped */
259 if (reduce_trace > -1 && trace_id != reduce_trace) continue ;
260
261 cpl_msg_info(__func__, "Process Order %d/Trace %d",order,trace_id) ;
262 cpl_msg_indent_more() ;
263
264 /* Get the input slit_func if available */
265 if (slit_func_in != NULL) {
266 /* Load the proper slit function vector */
267 cr2res_extract_SLIT_FUNC_get_vector(slit_func_in, order,
268 trace_id, &slit_func_in_vec) ;
269 } else {
270 slit_func_in_vec = NULL ;
271 }
272
273 /* Call the Extraction */
274 if (extr_method == CR2RES_EXTR_SUM) {
275 if (cr2res_extract_sum_vert(img, traces, order,
276 trace_id, extr_height, &(slit_func_vec[i]),
277 &(spectrum[i]), &model_loc_one) != 0) {
278 cpl_msg_error(__func__, "Cannot (sum-)extract the trace") ;
279 if (slit_func_in_vec != NULL)
280 cpl_vector_delete(slit_func_in_vec) ;
281 slit_func_vec[i] = NULL ;
282 spectrum[i] = NULL ;
283 model_loc_one = NULL ;
284 cpl_error_reset() ;
285 cpl_msg_indent_less() ;
286 continue ;
287 }
288 } else if (extr_method == CR2RES_EXTR_MEDIAN) {
289 if (cr2res_extract_median(img, traces, order,
290 trace_id, extr_height, &(slit_func_vec[i]),
291 &(spectrum[i]), &model_loc_one) != 0) {
292 cpl_msg_error(__func__, "Cannot (median-)extract the trace") ;
293 if (slit_func_in_vec != NULL)
294 cpl_vector_delete(slit_func_in_vec) ;
295 slit_func_vec[i] = NULL ;
296 spectrum[i] = NULL ;
297 model_loc_one = NULL ;
298 cpl_error_reset() ;
299 cpl_msg_indent_less() ;
300 continue ;
301 }
302 } else if (extr_method == CR2RES_EXTR_TILTSUM) {
303 if (cr2res_extract_sum_tilt(img, traces, order,
304 trace_id, extr_height, &(slit_func_vec[i]),
305 &(spectrum[i]), &model_loc_one) != 0) {
306 cpl_msg_error(__func__, "Cannot (tiltsum-)extract the trace") ;
307 if (slit_func_in_vec != NULL)
308 cpl_vector_delete(slit_func_in_vec) ;
309 slit_func_vec[i] = NULL ;
310 spectrum[i] = NULL ;
311 model_loc_one = NULL ;
312 cpl_error_reset() ;
313 cpl_msg_indent_less() ;
314 continue ;
315 }
316 } else if (extr_method == CR2RES_EXTR_OPT_CURV) {
317 if (cr2res_extract_slitdec_curved(img, traces, slit_func_in_vec,
318 order, trace_id, extr_height, swath_width,
319 oversample, pclip, smooth_slit, smooth_spec,
320 niter, kappa, error_factor,
321 &(slit_func_vec[i]),
322 &(spectrum[i]), &model_loc_one) != 0) {
323 cpl_msg_error(__func__,
324 "Cannot extract order %d, trace %d", order, trace_id) ;
325 if (slit_func_in_vec != NULL)
326 cpl_vector_delete(slit_func_in_vec) ;
327 slit_func_vec[i] = NULL ;
328 spectrum[i] = NULL ;
329 model_loc_one = NULL ;
330 cpl_error_reset() ;
331 cpl_msg_indent_less() ;
332 continue ;
333 }
334 }
335 if (slit_func_in_vec != NULL) cpl_vector_delete(slit_func_in_vec) ;
336
337 /* Correct the blaze if passed */
338 if (blaze_table_in != NULL) {
339 cpl_bivector * blaze_biv ;
340 cpl_bivector * blaze_err_biv ;
341 if (cr2res_extract_EXTRACT1D_get_spectrum(blaze_table_in, order,
342 trace_id, &blaze_biv, &blaze_err_biv)) {
343 cpl_msg_warning(__func__,
344 "Cannot Get the Blaze for order/trace:%d/%d - skip",
345 order, trace_id) ;
346 } else {
347 double * pblaze ;
348 double * pblaze_err ;
349 double first_nonzero_value, first_nonzero_error ;
350 int j ;
351 /* The Blaze needs to be 'cleaned from 0s before division */
352 pblaze =
353 cpl_vector_get_data(cpl_bivector_get_y(blaze_biv)) ;
354 pblaze_err =
355 cpl_vector_get_data(cpl_bivector_get_y(blaze_err_biv)) ;
356 first_nonzero_value = first_nonzero_error = 0.0 ;
357 for (j=0 ; j<cpl_bivector_get_size(blaze_biv) ; j++) {
358 if (fabs(pblaze[j])>1e-3) {
359 first_nonzero_value = pblaze[j] ;
360 first_nonzero_error = pblaze_err[j] ;
361 break ;
362 }
363 }
364 if (fabs(first_nonzero_value)<1e-3) {
365 cpl_msg_warning(__func__, "Blaze filled with zeros - skip");
366 } else {
367 double * pspec ;
368 double * pspec_err ;
369 double norm_factor, val, err ;
370 for (j=0 ; j<cpl_bivector_get_size(blaze_biv) ; j++) {
371 if (fabs(pblaze[j])<1e-3) {
372 pblaze[j] = first_nonzero_value ;
373 pblaze_err[j] = first_nonzero_error ;
374 }
375 }
376
377 /* Normalize the Blaze, prefer the normalization factor given in
378 * QC FLAT BLAZE NORM if present otherwise revert to trace-wise
379 * 95th percentile normalization
380 */
381
382 if (blaze_norm > 0){
383 norm_factor = blaze_norm;
384 } else {
385 cpl_vector * tmp_vec ;
386 cpl_size kth ;
387 tmp_vec=cpl_vector_duplicate(cpl_bivector_get_y(blaze_biv));
388 kth = (cpl_size)(cpl_bivector_get_size(blaze_biv)*0.95) ;
389 irplib_vector_get_kth(tmp_vec, kth) ;
390 norm_factor = cpl_vector_get(tmp_vec, kth) ;
391 cpl_vector_delete(tmp_vec) ;
392 }
393
394 /* Divide by the Blaze */
395 pspec = cpl_bivector_get_x_data(spectrum[i]) ;
396 pspec_err = cpl_bivector_get_y_data(spectrum[i]);
397 for (j=0 ; j<cpl_bivector_get_size(blaze_biv) ; j++) {
398 if (fabs(pspec[j]) > 1e-3) {
399 /* Apply division by normalized blaze */
400 val = pspec[j] / (pblaze[j]/norm_factor) ;
401
402 /* Error */
403 /* err(a/b)=abs(a/b)sqrt((err_a/a)^2+(err_b/b)^2) */
404 err = fabs(val) * sqrt(((pspec_err[j]*pspec_err[j])/
405 (pspec[j]*pspec[j]))+
406 ((pblaze_err[j]*pblaze_err[j])/
407 (pblaze[j]*pblaze[j])));
408 } else {
409 val = err = 0.0 ;
410 }
411
412 /* Set Results */
413 pspec[j] = val ;
414 pspec_err[j] = err ;
415 }
416 if (cpl_error_get_code()) {
417 cpl_error_reset();
418 cpl_msg_warning(__func__,
419 "Cannot Correct Blaze for order/trace:%d/%d - skip",
420 order, trace_id) ;
421 }
422 }
423 cpl_bivector_delete(blaze_biv) ;
424 cpl_bivector_delete(blaze_err_biv) ;
425 }
426 }
427
428 /* Keep the model for the in-order merge below */
429 model_one[i] = model_loc_one ;
430
431 cpl_msg_indent_less() ;
432 }
433
434 /* Update the model global image: merge the per-trace models in trace
435 order, and plot if requested. Equivalent to setting every pixel
436 where the trace model is non-zero and not rejected (cpl_image_set
437 semantics: the written pixel is accepted). Direct buffer access:
438 the per-pixel hdrl calls dominated the runtime of this function. */
439 for (i=0 ; i<nb_traces ; i++) {
440 if (model_one[i] != NULL) {
441 const cpl_image * src_data =
442 hdrl_image_get_image_const(model_one[i]) ;
443 const cpl_image * src_err =
444 hdrl_image_get_error_const(model_one[i]) ;
445 const double * sd = cpl_image_get_data_double_const(src_data) ;
446 const double * se = cpl_image_get_data_double_const(src_err) ;
447 const cpl_mask * src_bpm = cpl_image_get_bpm_const(src_data) ;
448 const cpl_binary * sb =
449 src_bpm ? cpl_mask_get_data_const(src_bpm) : NULL ;
450 cpl_image * dst_data = hdrl_image_get_image(model_loc) ;
451 cpl_image * dst_err = hdrl_image_get_error(model_loc) ;
452 double * dd = cpl_image_get_data_double(dst_data) ;
453 double * de = cpl_image_get_data_double(dst_err) ;
454 /* clearing flags is only needed where a map exists */
455 cpl_binary * dbd = cpl_image_get_bpm_const(dst_data) ?
456 cpl_mask_get_data(cpl_image_get_bpm(dst_data)) : NULL ;
457 cpl_binary * dbe = cpl_image_get_bpm_const(dst_err) ?
458 cpl_mask_get_data(cpl_image_get_bpm(dst_err)) : NULL ;
459 cpl_size npix = hdrl_image_get_size_x(model_loc)
460 * hdrl_image_get_size_y(model_loc) ;
461 cpl_size k ;
462 for (k = 0 ; k < npix ; k++) {
463 if (sd[k] != 0 && (sb == NULL || !sb[k])) {
464 dd[k] = sd[k] ;
465 de[k] = se[k] ;
466 if (dbd != NULL) dbd[k] = CPL_BINARY_0 ;
467 if (dbe != NULL) dbe[k] = CPL_BINARY_0 ;
468 }
469 }
470 hdrl_image_delete(model_one[i]) ;
471 model_one[i] = NULL ;
472 }
473
474 /* Plot the Spectrum */
475 if (display && spectrum[i] != NULL &&
476 disp_order_idx ==
477 cpl_table_get(traces, CR2RES_COL_ORDER, i, NULL) &&
478 disp_trace ==
479 cpl_table_get(traces, CR2RES_COL_TRACENB, i, NULL)) {
480 cpl_plot_vector(
481 "set grid;set xlabel 'pixels';set ylabel 'Flux (ADU)';",
482 "t 'Extracted Spectrum' w lines", "",
483 cpl_bivector_get_x_const(spectrum[i])) ;
484 }
485 }
486 cpl_free(model_one) ;
487
488 /* Create the slit_func_tab for the current detector */
489 if ((slit_func_loc = cr2res_extract_SLITFUNC_create(slit_func_vec,
490 traces)) == NULL) {
491 cpl_msg_error(__func__, "Cannot compute the slit function") ;
492 for (i=0 ; i<nb_traces ; i++) {
493 if (slit_func_vec[i] != NULL) cpl_vector_delete(slit_func_vec[i]) ;
494 if (spectrum[i] != NULL) cpl_bivector_delete(spectrum[i]) ;
495 }
496 cpl_free(spectrum) ;
497 cpl_free(slit_func_vec) ;
498 hdrl_image_delete(model_loc) ;
499 return -1;
500 }
501
502 /* Create the extracted_tab for the current detector */
503 if ((extract_loc = cr2res_extract_EXTRACT1D_create(spectrum, traces))
504 == NULL) {
505 for (i=0 ; i<nb_traces ; i++) {
506 if (slit_func_vec[i] != NULL) cpl_vector_delete(slit_func_vec[i]) ;
507 if (spectrum[i] != NULL) cpl_bivector_delete(spectrum[i]) ;
508 }
509 cpl_free(spectrum) ;
510 cpl_free(slit_func_vec) ;
511 hdrl_image_delete(model_loc) ;
512 cpl_table_delete(slit_func_loc);
513 return -1;
514 }
515
516 /* Deallocate Vectors */
517 for (i=0 ; i<nb_traces ; i++) {
518 if (slit_func_vec[i] != NULL) cpl_vector_delete(slit_func_vec[i]) ;
519 if (spectrum[i] != NULL) cpl_bivector_delete(spectrum[i]) ;
520 }
521 cpl_free(spectrum) ;
522 cpl_free(slit_func_vec) ;
523
524 /* Return */
525 *extracted = extract_loc ;
526 *slit_func = slit_func_loc ;
527 *model_master = model_loc;
528 return 0 ;
529}
530
531/*----------------------------------------------------------------------------*/
554/*----------------------------------------------------------------------------*/
555int cr2res_extract_sum_vert(
556 const hdrl_image * hdrl_in,
557 const cpl_table * trace_tab,
558 int order,
559 int trace_id,
560 int height,
561 cpl_vector ** slit_func,
562 cpl_bivector ** spec,
563 hdrl_image ** model)
564{
565 int * ycen_int;
566 cpl_vector * ycen ;
567 cpl_image * img_tmp;
568 cpl_image * img_1d;
569 const cpl_image * img_in;
570 const cpl_image * err_in;
571 cpl_vector * spc;
572 cpl_vector * slitfu;
573 cpl_vector * sigma;
574 cpl_size lenx, leny;
575 int i, j, y;
576 double trace_cen, trace_height ;
577 cpl_type imtyp;
578
579 /* Check Entries */
580 if (hdrl_in == NULL || trace_tab == NULL) return -1 ;
581
582 img_in = hdrl_image_get_image_const(hdrl_in);
583 err_in = hdrl_image_get_error_const(hdrl_in);
584
585 /* use the same type as input for temp images below */
586 imtyp = cpl_image_get_type(img_in);
587 lenx = cpl_image_get_size_x(img_in);
588 leny = cpl_image_get_size_y(img_in);
589
590 /* Compute height if not given */
591 if (height <= 0) {
592 height = cr2res_trace_get_height(trace_tab, order, trace_id);
593 if (height <= 0) {
594 cpl_msg_error(__func__, "Cannot compute height");
595 return -1;
596 }
597 }
598 /* Get ycen */
599 if ((ycen = cr2res_trace_get_ycen(trace_tab, order,
600 trace_id, lenx)) == NULL) {
601 cpl_msg_error(__func__, "Cannot get ycen");
602 return -1 ;
603 }
604 trace_cen = cpl_vector_get(ycen, cpl_vector_get_size(ycen)/2) ;
605 trace_height = (double)cr2res_trace_get_height(trace_tab, order, trace_id) ;
606 cpl_msg_info(__func__, "Y position of the trace: %g -> %g",
607 trace_cen-(trace_height/2), trace_cen+(trace_height/2)) ;
608
609 img_tmp = cr2res_image_cut_rectify(img_in, ycen, height);
610 if (img_tmp == NULL) {
611 cpl_msg_error(__func__, "Cannot rectify order");
612 cpl_vector_delete(ycen);
613 return -1;
614 }
615
616 if (cpl_msg_get_level() == CPL_MSG_DEBUG) {
617 cpl_image_save(img_tmp, "debug_rectorder.fits", imtyp,
618 NULL, CPL_IO_CREATE);
619 }
620 img_1d = cpl_image_collapse_create(img_tmp, 0); // sum of rows
621 spc = cpl_vector_new_from_image_row(img_1d, 1);
622 cpl_image_delete(img_1d);
623
624 img_1d = cpl_image_collapse_create(img_tmp, 1); // sum of columns
625 slitfu = cpl_vector_new_from_image_column(img_1d, 1);
626 cpl_vector_divide_scalar(slitfu, cpl_vector_get_sum(slitfu));
627 cpl_image_delete(img_1d);
628 cpl_image_delete(img_tmp);
629
630 img_tmp = cr2res_image_cut_rectify(err_in, ycen, height);
631 if (img_tmp == NULL) {
632 cpl_msg_error(__func__, "Cannot rectify error");
633 cpl_vector_delete(ycen);
634 return -1;
635 }
636 cpl_image_multiply(img_tmp, img_tmp);
637 img_1d = cpl_image_collapse_create(img_tmp, 0);
638 sigma = cpl_vector_new_from_image_row(img_1d, 1);
639 cpl_vector_sqrt(sigma);
640 cpl_image_delete(img_tmp);
641 cpl_image_delete(img_1d);
642
643 // reconstruct the 2d image with the "model"
644 img_tmp = cpl_image_new(lenx, leny, imtyp);
645 ycen_int = cr2res_vector_get_int(ycen);
646 for (i=1;i<=lenx;i++){
647 for (j=1;j<=height;j++){
648 y = ycen_int[i-1]-(height/2)+j;
649 if ((y <=0) || (y > leny)){ continue; }
650 cpl_image_set(img_tmp, i, y,
651 cpl_vector_get(spc,i-1)*cpl_vector_get(slitfu,j-1) );
652 }
653 }
654 cpl_vector_delete(ycen);
655 cpl_free(ycen_int);
656
657 if (cpl_msg_get_level() == CPL_MSG_DEBUG) {
658 cpl_image_save(img_tmp, "debug_model.fits", imtyp,
659 NULL, CPL_IO_CREATE);
660 }
661
662
663 *slit_func = slitfu;
664 *spec = cpl_bivector_wrap_vectors(spc, sigma);
665 *model = hdrl_image_create(img_tmp, NULL);
666 cpl_image_delete(img_tmp);
667
668 if (cpl_error_get_code() != CPL_ERROR_NONE){
669 cpl_msg_error(__func__, "Error in the vertical sum extraction %s",
670 cpl_error_get_where());
671 cpl_error_reset();
672 return -1;
673 } else {
674 return 0;
675 }
676}
677
678/*----------------------------------------------------------------------------*/
702/*----------------------------------------------------------------------------*/
703int cr2res_extract_median(
704 const hdrl_image * hdrl_in,
705 const cpl_table * trace_tab,
706 int order,
707 int trace_id,
708 int height,
709 cpl_vector ** slit_func,
710 cpl_bivector ** spec,
711 hdrl_image ** model)
712{
713 int * ycen_int;
714 cpl_vector * ycen ;
715 cpl_image * img_tmp;
716 cpl_image * img_1d;
717 const cpl_image * img_in;
718 const cpl_image * err_in;
719 cpl_vector * spc;
720 cpl_vector * slitfu;
721 cpl_vector * sigma;
722 cpl_size lenx, leny;
723 int i, j;
724 double trace_cen, trace_height ;
725 cpl_type imtyp;
726
727 /* Check Entries */
728 if (hdrl_in == NULL || trace_tab == NULL) return -1 ;
729
730 img_in = hdrl_image_get_image_const(hdrl_in);
731 err_in = hdrl_image_get_error_const(hdrl_in);
732
733 /* use the same type as input for temp images below */
734 imtyp = cpl_image_get_type(img_in);
735 lenx = cpl_image_get_size_x(img_in);
736 leny = cpl_image_get_size_y(img_in);
737
738 /* Compute height if not given */
739 if (height <= 0) {
740 height = cr2res_trace_get_height(trace_tab, order, trace_id);
741 if (height <= 0) {
742 cpl_msg_error(__func__, "Cannot compute height");
743 return -1;
744 }
745 }
746 /* Get ycen */
747 if ((ycen = cr2res_trace_get_ycen(trace_tab, order,
748 trace_id, lenx)) == NULL) {
749 cpl_msg_error(__func__, "Cannot get ycen");
750 return -1 ;
751 }
752 trace_cen = cpl_vector_get(ycen, cpl_vector_get_size(ycen)/2) ;
753 trace_height = (double)cr2res_trace_get_height(trace_tab, order, trace_id) ;
754 cpl_msg_info(__func__, "Y position of the trace: %g -> %g",
755 trace_cen-(trace_height/2), trace_cen+(trace_height/2)) ;
756
757 img_tmp = cr2res_image_cut_rectify(img_in, ycen, height);
758 if (img_tmp == NULL) {
759 cpl_msg_error(__func__, "Cannot rectify order");
760 cpl_vector_delete(ycen);
761 return -1;
762 }
763
764 if (cpl_msg_get_level() == CPL_MSG_DEBUG) {
765 cpl_image_save(img_tmp, "debug_rectorder.fits", imtyp,
766 NULL, CPL_IO_CREATE);
767 }
768 img_1d = cpl_image_collapse_median_create(img_tmp, 0, 0, 0); // rows
769 spc = cpl_vector_new_from_image_row(img_1d, 1);
770
771 // Multiply so that scaling matches a vertical sum
772 cpl_vector_multiply_scalar(spc,(double)height);
773 cpl_image_delete(img_1d);
774
775 img_1d = cpl_image_collapse_median_create(img_tmp, 1, 0, 0); // cols
776 slitfu = cpl_vector_new_from_image_column(img_1d, 1);
777 cpl_vector_divide_scalar(slitfu, cpl_vector_get_sum(slitfu));
778 cpl_image_delete(img_1d);
779 cpl_image_delete(img_tmp);
780
781 img_tmp = cr2res_image_cut_rectify(err_in, ycen, height);
782 if (img_tmp == NULL)
783 {
784 cpl_msg_error(__func__, "Cannot rectify error");
785 cpl_vector_delete(ycen);
786 return -1;
787 }
788 cpl_image_multiply(img_tmp, img_tmp);
789 img_1d = cpl_image_collapse_create(img_tmp, 0);
790 sigma = cpl_vector_new_from_image_row(img_1d, 1);
791 cpl_vector_sqrt(sigma);
792 cpl_image_delete(img_tmp);
793 cpl_image_delete(img_1d);
794
795 // reconstruct the 2d image with the "model"
796 img_tmp = cpl_image_new(lenx, leny, imtyp);
797 ycen_int = cr2res_vector_get_int(ycen);
798 for (i=1;i<=lenx;i++){
799 for (j=1;j<=height;j++){
800 cpl_image_set(img_tmp, i, ycen_int[i-1]-(height/2)+j,
801 cpl_vector_get(spc,i-1)*cpl_vector_get(slitfu,j-1) );
802 }
803 }
804 cpl_vector_delete(ycen);
805 cpl_free(ycen_int);
806
807 if (cpl_msg_get_level() == CPL_MSG_DEBUG) {
808 cpl_image_save(img_tmp, "debug_model.fits", imtyp,
809 NULL, CPL_IO_CREATE);
810 }
811
812
813 *slit_func = slitfu;
814 *spec = cpl_bivector_wrap_vectors(spc, sigma);
815 *model = hdrl_image_create(img_tmp, NULL);
816 cpl_image_delete(img_tmp);
817
818 return 0;
819}
820
821/*----------------------------------------------------------------------------*/
845/*----------------------------------------------------------------------------*/
846int cr2res_extract_sum_tilt(
847 const hdrl_image * hdrl_in,
848 const cpl_table * trace_tab,
849 int order,
850 int trace_id,
851 int height,
852 cpl_vector ** slit_func,
853 cpl_bivector ** spec,
854 hdrl_image ** model)
855{
856 int * ycen_int;
857 cpl_vector * ycen ;
858 cpl_image * img_tmp;
859 cpl_image * img_1d;
860 const cpl_image * img_in;
861 const cpl_image * err_in;
862 cpl_vector * spc;
863 cpl_vector * slitfu;
864 cpl_vector * sigma;
865 cpl_size lenx, leny;
866 int i, j;
867 double trace_cen, trace_height ;
868 cpl_type imtyp;
869
870 int badpix;
871 double a, b, c, value;
872 cpl_polynomial * slitcurve_A, * slitcurve_B, *slitcurve_C;
873 cpl_bivector * xi, *xt;
874
875 /* Check Entries */
876 if (hdrl_in == NULL || trace_tab == NULL) return -1 ;
877
878 img_in = hdrl_image_get_image_const(hdrl_in);
879 err_in = hdrl_image_get_error_const(hdrl_in);
880
881 /* use the same type as input for temp images below */
882 imtyp = cpl_image_get_type(img_in);
883 lenx = cpl_image_get_size_x(img_in);
884 leny = cpl_image_get_size_y(img_in);
885
886 /* Compute height if not given */
887 if (height <= 0) {
888 height = cr2res_trace_get_height(trace_tab, order, trace_id);
889 if (height <= 0) {
890 cpl_msg_error(__func__, "Cannot compute height");
891 return -1;
892 }
893 }
894 /* Get ycen */
895 if ((ycen = cr2res_trace_get_ycen(trace_tab, order,
896 trace_id, lenx)) == NULL) {
897 cpl_msg_error(__func__, "Cannot get ycen");
898 return -1 ;
899 }
900 trace_cen = cpl_vector_get(ycen, cpl_vector_get_size(ycen)/2) ;
901 trace_height = (double)cr2res_trace_get_height(trace_tab, order, trace_id) ;
902 cpl_msg_info(__func__, "Y position of the trace: %g -> %g",
903 trace_cen-(trace_height/2), trace_cen+(trace_height/2)) ;
904
905 img_tmp = cr2res_image_cut_rectify(img_in, ycen, height);
906 if (img_tmp == NULL) {
907 cpl_msg_error(__func__, "Cannot rectify order");
908 cpl_vector_delete(ycen);
909 return -1;
910 }
911
912 if (cpl_msg_get_level() == CPL_MSG_DEBUG) {
913 cpl_image_save(img_tmp, "debug_rectorder.fits", imtyp,
914 NULL, CPL_IO_CREATE);
915 }
916
917 // Fix curvature if existing
918 // def correct_for_curvature(img_order, tilt, shear, xwd):
919 // img_order = np.ma.filled(img_order, 0)
920 // xt = np.arange(img_order.shape[1])
921 // for y, yt in zip(range(xwd[0] + xwd[1]), range(-xwd[0], xwd[1])):
922 // xi = xt + yt * tilt + yt ** 2 * shear
923 // img_order[y] = np.interp(xi, xt, img_order[y])
924 // return img_order
925
926 // Read the slitcurvature from the trace table
927 slitcurve_A = cr2res_get_trace_wave_poly(trace_tab, CR2RES_COL_SLIT_CURV_A,
928 order, trace_id);
929 slitcurve_B = cr2res_get_trace_wave_poly(trace_tab, CR2RES_COL_SLIT_CURV_B,
930 order, trace_id);
931 slitcurve_C = cr2res_get_trace_wave_poly(trace_tab, CR2RES_COL_SLIT_CURV_C,
932 order, trace_id);
933
934 if ((slitcurve_A == NULL) || (slitcurve_B == NULL) || (slitcurve_C == NULL))
935 {
936 cpl_msg_error(__func__,
937 "No (or incomplete) slitcurve data found in trace table");
938 cpl_vector_delete(ycen);
939 cpl_polynomial_delete(slitcurve_A);
940 cpl_polynomial_delete(slitcurve_B);
941 cpl_polynomial_delete(slitcurve_C);
942 return -1;
943 }
944
945
946 // xi is the regular coordinates
947 // xt is the shifted coordinates
948 xi = cpl_bivector_new(lenx);
949 xt = cpl_bivector_new(lenx);
950 for (i = 0; i < lenx; i++){
951 cpl_vector_set(cpl_bivector_get_x(xi), i, i);
952 cpl_vector_set(cpl_bivector_get_x(xt), i, i);
953 cpl_vector_set(cpl_bivector_get_y(xi), i, 0);
954 cpl_vector_set(cpl_bivector_get_y(xt), i, 0);
955 }
956
957 for (i = 0; i < height; i++){
958 int yt = i - height / 2 + 0.5;
959 int yc = cpl_vector_get(ycen, i);
960
961 for (j = 1; j < lenx - 1; j++){
962 //a = cpl_polynomial_eval_1d(slitcurve_A, j, NULL); this is ignored apparently?
963 b = cpl_polynomial_eval_1d(slitcurve_B, j, NULL);
964 c = cpl_polynomial_eval_1d(slitcurve_C, j, NULL);
965
966 // shift polynomial to local frame
967 // a = a - j + yc * b + yc * yc * c;
968 a = 0;
969 b += 2 * yc * c;
970
971 value = j - a - yt * b - yt * yt * c;
972 value = max(min(value, lenx-1), 0);
973 cpl_vector_set(cpl_bivector_get_x(xt), j, value);
974 value = cpl_image_get(img_tmp, j + 1, i + 1, &badpix);
975 if (badpix) value = NAN;
976 cpl_vector_set(cpl_bivector_get_y(xt), j, value);
977 }
978
979 cpl_bivector_interpolate_linear(xi, xt);
980
981 for (j = 0; j < lenx; j++){
982 value = cpl_vector_get(cpl_bivector_get_y(xi), j);
983 cpl_image_set(img_tmp, j+1, i+1, value);
984 if (isnan(value)){
985 cpl_image_set(img_tmp, j+1, i+1, 0);
986 cpl_image_reject(img_tmp, j+1, i+1);
987 }
988 }
989 }
990
991 cpl_bivector_delete(xi);
992 cpl_bivector_delete(xt);
993 cpl_polynomial_delete(slitcurve_A);
994 cpl_polynomial_delete(slitcurve_B);
995 cpl_polynomial_delete(slitcurve_C);
996
997 if (cpl_msg_get_level() == CPL_MSG_DEBUG) {
998 cpl_image_save(img_tmp, "debug_img_shifted.fits", imtyp,
999 NULL, CPL_IO_CREATE);
1000 }
1001
1002 // Fill bad pixels with linear approximation
1003 // xi is the input coordinates
1004 // xt is the target coordinates
1005 xi = cpl_bivector_new(height);
1006 xt = cpl_bivector_new(height);
1007
1008 for (j = 0; j < lenx; j++){
1009 // We need to set the first and last point in case they are bad pixels and need to be interpolated
1010 cpl_vector_set(cpl_bivector_get_x(xi), 0, 0);
1011 cpl_vector_set(cpl_bivector_get_y(xi), 0, 0);
1012 cpl_vector_set(cpl_bivector_get_x(xi), height-1, height-1);
1013 cpl_vector_set(cpl_bivector_get_y(xi), height-1, 0);
1014 for (i = 0; i < height; i++){
1015 cpl_vector_set(cpl_bivector_get_x(xt), i, i);
1016 cpl_vector_set(cpl_bivector_get_y(xt), i, 0);
1017
1018 value = cpl_image_get(img_tmp, j+1, i+1, &badpix);
1019 if (!badpix){
1020 cpl_vector_set(cpl_bivector_get_x(xi), i, i);
1021 cpl_vector_set(cpl_bivector_get_y(xi), i, value);
1022 }
1023 }
1024 cpl_bivector_interpolate_linear(xt, xi);
1025 for (i = 0; i < height; i++){
1026 value = cpl_vector_get(cpl_bivector_get_y(xt), i);
1027 cpl_image_set(img_tmp, j+1, i+1, value);
1028 }
1029 }
1030
1031 cpl_bivector_delete(xi);
1032 cpl_bivector_delete(xt);
1033
1034 if (cpl_msg_get_level() == CPL_MSG_DEBUG) {
1035 cpl_image_save(img_tmp, "debug_img_shifted_corrected.fits", imtyp,
1036 NULL, CPL_IO_CREATE);
1037 }
1038
1039 img_1d = cpl_image_collapse_create(img_tmp, 0); // sum of rows
1040 spc = cpl_vector_new_from_image_row(img_1d, 1);
1041 cpl_image_delete(img_1d);
1042
1043 img_1d = cpl_image_collapse_create(img_tmp, 1); // sum of columns
1044 slitfu = cpl_vector_new_from_image_column(img_1d, 1);
1045 cpl_vector_divide_scalar(slitfu, cpl_vector_get_sum(slitfu));
1046 cpl_image_delete(img_1d);
1047 cpl_image_delete(img_tmp);
1048
1049 img_tmp = cr2res_image_cut_rectify(err_in, ycen, height);
1050 if (img_tmp == NULL)
1051 {
1052 cpl_msg_error(__func__, "Cannot rectify error");
1053 cpl_vector_delete(ycen);
1054 return -1;
1055 }
1056 cpl_image_multiply(img_tmp, img_tmp);
1057 img_1d = cpl_image_collapse_create(img_tmp, 0);
1058 sigma = cpl_vector_new_from_image_row(img_1d, 1);
1059 cpl_vector_sqrt(sigma);
1060 cpl_image_delete(img_tmp);
1061 cpl_image_delete(img_1d);
1062
1063 // reconstruct the 2d image with the "model"
1064 img_tmp = cpl_image_new(lenx, leny, imtyp);
1065 ycen_int = cr2res_vector_get_int(ycen);
1066 for (i=1;i<=lenx;i++){
1067 for (j=1;j<=height;j++){
1068 cpl_image_set(img_tmp, i, ycen_int[i-1]-(height/2)+j,
1069 cpl_vector_get(spc,i-1)*cpl_vector_get(slitfu,j-1) );
1070 }
1071 }
1072 cpl_vector_delete(ycen);
1073 cpl_free(ycen_int);
1074
1075 if (cpl_msg_get_level() == CPL_MSG_DEBUG) {
1076 cpl_image_save(img_tmp, "debug_model.fits", imtyp,
1077 NULL, CPL_IO_CREATE);
1078 }
1079
1080
1081 *slit_func = slitfu;
1082 *spec = cpl_bivector_wrap_vectors(spc, sigma);
1083 *model = hdrl_image_create(img_tmp, NULL);
1084 cpl_image_delete(img_tmp);
1085
1086 return 0;
1087}
1088
1089/*----------------------------------------------------------------------------*/
1096/*----------------------------------------------------------------------------*/
1097cpl_table * cr2res_extract_EXTRACT1D_create(
1098 cpl_bivector ** spectrum,
1099 const cpl_table * trace_table)
1100{
1101 cpl_table * out ;
1102 char * col_name ;
1103 const double * pspec ;
1104 const double * perr ;
1105 cpl_vector * wave_vec ;
1106 const double * pwl ;
1107 int nrows, all_null, i, order, trace_id, nb_traces ;
1108
1109 /* Check entries */
1110 if (spectrum == NULL || trace_table == NULL) return NULL ;
1111
1112 /* Initialise */
1113 nb_traces = cpl_table_get_nrow(trace_table) ;
1114
1115 /* Check if all vectors are not null */
1116 all_null = 1 ;
1117 for (i=0 ; i<nb_traces ; i++)
1118 if (spectrum[i] != NULL) {
1119 nrows = cpl_bivector_get_size(spectrum[i]) ;
1120 all_null = 0 ;
1121 }
1122 if (all_null == 1) return NULL ;
1123
1124 /* Check the sizes */
1125 for (i=0 ; i<nb_traces ; i++)
1126 if (spectrum[i] != NULL && cpl_bivector_get_size(spectrum[i]) != nrows)
1127 return NULL ;
1128
1129 /* Create the table */
1130 out = cpl_table_new(nrows);
1131 for (i=0 ; i<nb_traces ; i++) {
1132 order = cpl_table_get(trace_table, CR2RES_COL_ORDER, i, NULL) ;
1133 trace_id = cpl_table_get(trace_table, CR2RES_COL_TRACENB, i, NULL) ;
1134 /* Create SPEC column */
1135 col_name = cr2res_dfs_SPEC_colname(order, trace_id) ;
1136 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
1137 cpl_free(col_name) ;
1138 /* Create SPEC_ERR column */
1139 col_name = cr2res_dfs_SPEC_ERR_colname(order, trace_id) ;
1140 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
1141 cpl_free(col_name) ;
1142 /* Create WAVELENGTH column */
1143 col_name = cr2res_dfs_WAVELENGTH_colname(order, trace_id) ;
1144 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
1145 cpl_free(col_name) ;
1146 }
1147
1148 /* Fill the table */
1149 for (i=0 ; i<nb_traces ; i++) {
1150 order = cpl_table_get(trace_table, CR2RES_COL_ORDER, i, NULL) ;
1151 trace_id = cpl_table_get(trace_table, CR2RES_COL_TRACENB, i, NULL) ;
1152 if (spectrum[i] != NULL) {
1153 pspec = cpl_bivector_get_x_data_const(spectrum[i]) ;
1154 perr = cpl_bivector_get_y_data_const(spectrum[i]);
1155 /* Fill SPEC column */
1156 col_name = cr2res_dfs_SPEC_colname(order, trace_id) ;
1157 cpl_table_copy_data_double(out, col_name, pspec) ;
1158 cpl_free(col_name) ;
1159 /* Fill SPEC_ERR column */
1160 col_name = cr2res_dfs_SPEC_ERR_colname(order, trace_id) ;
1161 cpl_table_copy_data_double(out, col_name, perr) ;
1162 cpl_free(col_name) ;
1163
1164 /* Compute Wavelength column */
1165 wave_vec = cr2res_trace_get_wl(trace_table, order, trace_id,
1166 CR2RES_DETECTOR_SIZE);
1167 if (wave_vec == NULL) {
1168 wave_vec = cpl_vector_new(CR2RES_DETECTOR_SIZE) ;
1169 cpl_vector_fill(wave_vec, 0.0) ;
1170 }
1171 pwl = cpl_vector_get_data(wave_vec) ;
1172
1173 /* Fill WAVELENGTH column */
1174 col_name = cr2res_dfs_WAVELENGTH_colname(order, trace_id) ;
1175 cpl_table_copy_data_double(out, col_name, pwl) ;
1176 cpl_free(col_name) ;
1177 cpl_vector_delete(wave_vec) ;
1178 } else {
1179 /* Fill SPEC column */
1180 col_name = cr2res_dfs_SPEC_colname(order, trace_id) ;
1181 cpl_table_fill_column_window_double(out, col_name, 0, CR2RES_DETECTOR_SIZE, NAN);
1182 cpl_table_set_column_invalid(out, col_name, 0, CR2RES_DETECTOR_SIZE);
1183 cpl_free(col_name) ;
1184 /* Fill SPEC_ERR column */
1185 col_name = cr2res_dfs_SPEC_ERR_colname(order, trace_id) ;
1186 cpl_table_fill_column_window_double(out, col_name, 0, CR2RES_DETECTOR_SIZE, NAN);
1187 cpl_table_set_column_invalid(out, col_name, 0, CR2RES_DETECTOR_SIZE);
1188 cpl_free(col_name) ;
1189 /* Fill WAVELENGTH column */
1190 col_name = cr2res_dfs_WAVELENGTH_colname(order, trace_id) ;
1191 cpl_table_fill_column_window_double(out, col_name, 0, CR2RES_DETECTOR_SIZE, NAN);
1192 cpl_table_set_column_invalid(out, col_name, 0, CR2RES_DETECTOR_SIZE);
1193 cpl_free(col_name) ;
1194 }
1195 }
1196 return out ;
1197}
1198
1199/*----------------------------------------------------------------------------*/
1206/*----------------------------------------------------------------------------*/
1207cpl_table * cr2res_extract_SLITFUNC_create(
1208 cpl_vector ** slit_func,
1209 const cpl_table * trace_table)
1210{
1211 cpl_table * out ;
1212 char * col_name ;
1213 const double * pslit ;
1214 int nrows, nrows_max, all_null, i, order, trace_id,
1215 nb_traces ;
1216
1217 /* Check entries */
1218 if (slit_func == NULL || trace_table == NULL) return NULL ;
1219
1220 /* Initialise */
1221 nrows_max = -1 ;
1222 nb_traces = cpl_table_get_nrow(trace_table) ;
1223
1224 /* Check that all vectors are not null */
1225 all_null = 1 ;
1226 for (i=0 ; i<nb_traces ; i++)
1227 if (slit_func[i] != NULL) {
1228 nrows = cpl_vector_get_size(slit_func[i]) ;
1229 if (nrows > nrows_max) nrows_max = nrows ;
1230 all_null = 0 ;
1231 }
1232 if (all_null == 1) return NULL ;
1233
1234 /* Resize all vectors to nrows_max, otherwise cpl_table_copy_data_double
1235 accesses values outside the vector*/
1236 for (i = 0; i < nb_traces; i++)
1237 if (slit_func[i] != NULL)
1238 {
1239 nrows = cpl_vector_get_size(slit_func[i]);
1240 cpl_vector_set_size(slit_func[i], nrows_max);
1241 for (int j = nrows; j < nrows_max; j++)
1242 {
1243 cpl_vector_set(slit_func[i], j, 0.0);
1244 }
1245 }
1246
1247 /* Create the table */
1248 out = cpl_table_new(nrows_max);
1249 for (i=0 ; i<nb_traces ; i++) {
1250 order = cpl_table_get(trace_table, CR2RES_COL_ORDER, i, NULL) ;
1251 trace_id = cpl_table_get(trace_table, CR2RES_COL_TRACENB, i, NULL) ;
1252 col_name = cr2res_dfs_SLIT_FUNC_colname(order, trace_id) ;
1253 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
1254 cpl_free(col_name) ;
1255 }
1256
1257 /* Fill the table */
1258 for (i=0 ; i<nb_traces ; i++) {
1259 if (slit_func[i] != NULL) {
1260 order = cpl_table_get(trace_table, CR2RES_COL_ORDER, i, NULL) ;
1261 trace_id = cpl_table_get(trace_table, CR2RES_COL_TRACENB, i, NULL) ;
1262 pslit = cpl_vector_get_data_const(slit_func[i]) ;
1263 col_name = cr2res_dfs_SLIT_FUNC_colname(order, trace_id) ;
1264 cpl_table_copy_data_double(out, col_name, pslit) ;
1265 cpl_free(col_name) ;
1266 } else {
1267 order = cpl_table_get(trace_table, CR2RES_COL_ORDER, i, NULL) ;
1268 trace_id = cpl_table_get(trace_table, CR2RES_COL_TRACENB, i, NULL) ;
1269 col_name = cr2res_dfs_SLIT_FUNC_colname(order, trace_id) ;
1270 cpl_table_fill_column_window_double(out, col_name, 0, nrows_max, NAN);
1271 cpl_table_set_column_invalid(out, col_name, 0, nrows_max);
1272 cpl_free(col_name) ;
1273 }
1274 }
1275 return out ;
1276}
1277
1278/*----------------------------------------------------------------------------*/
1288/*----------------------------------------------------------------------------*/
1289int cr2res_extract_EXTRACT1D_get_spectrum(
1290 const cpl_table * tab,
1291 int order,
1292 int trace_nb,
1293 cpl_bivector ** spec,
1294 cpl_bivector ** spec_err)
1295{
1296 char * spec_name ;
1297 char * wave_name ;
1298 char * spec_err_name ;
1299 const double * pspec ;
1300 const double * pwave ;
1301 const double * pspec_err ;
1302 double * pxspec ;
1303 double * pyspec ;
1304 double * pxspec_err ;
1305 double * pyspec_err ;
1306 int i, tab_size ;
1307
1308 /* Check entries */
1309 if (tab == NULL || spec == NULL || spec_err == NULL) return -1 ;
1310
1311 /* Get the Spectrum */
1312 spec_name = cr2res_dfs_SPEC_colname(order, trace_nb) ;
1313 if ((pspec = cpl_table_get_data_double_const(tab, spec_name)) == NULL) {
1314 cpl_msg_error(__func__, "Cannot find the spectrum") ;
1315 cpl_free(spec_name) ;
1316 return -1 ;
1317 }
1318 cpl_free(spec_name) ;
1319
1320 /* Get the Wavelength */
1321 wave_name = cr2res_dfs_WAVELENGTH_colname(order, trace_nb) ;
1322 if ((pwave = cpl_table_get_data_double_const(tab, wave_name)) == NULL) {
1323 cpl_msg_error(__func__, "Cannot find the wavelength") ;
1324 cpl_free(wave_name) ;
1325 return -1 ;
1326 }
1327 cpl_free(wave_name) ;
1328
1329 /* Get the Spectrum Error */
1330 spec_err_name = cr2res_dfs_SPEC_ERR_colname(order, trace_nb) ;
1331 if ((pspec_err = cpl_table_get_data_double_const(tab, spec_err_name))
1332 == NULL) {
1333 cpl_msg_error(__func__, "Cannot find the spectrum error") ;
1334 cpl_free(spec_err_name) ;
1335 return -1 ;
1336 }
1337 cpl_free(spec_err_name) ;
1338
1339 /* Create the output */
1340 tab_size = cpl_table_get_nrow(tab) ;
1341 *spec = cpl_bivector_new(tab_size) ;
1342 *spec_err = cpl_bivector_new(tab_size) ;
1343 pxspec = cpl_bivector_get_x_data(*spec) ;
1344 pyspec = cpl_bivector_get_y_data(*spec) ;
1345 pxspec_err = cpl_bivector_get_x_data(*spec_err) ;
1346 pyspec_err = cpl_bivector_get_y_data(*spec_err) ;
1347 for (i=0 ; i<tab_size ; i++) {
1348 pxspec[i] = pwave[i] ;
1349 pyspec[i] = pspec[i] ;
1350 pxspec_err[i] = pwave[i] ;
1351 pyspec_err[i] = pspec_err[i] ;
1352 }
1353 return 0 ;
1354}
1355
1356/*----------------------------------------------------------------------------*/
1365/*----------------------------------------------------------------------------*/
1366int cr2res_extract_SLIT_FUNC_get_vector(
1367 const cpl_table * tab,
1368 int order,
1369 int trace_nb,
1370 cpl_vector ** vec)
1371{
1372 char * vec_name ;
1373 const double * ptab ;
1374 double * pvec ;
1375 int i, tab_size ;
1376
1377 /* Check entries */
1378 if (tab == NULL || vec == NULL) return -1 ;
1379
1380 /* Get the Vector */
1381 vec_name = cr2res_dfs_SLIT_FUNC_colname(order, trace_nb) ;
1382 if ((ptab = cpl_table_get_data_double_const(tab, vec_name)) == NULL) {
1383 cpl_msg_error(__func__, "Cannot find the slit_func") ;
1384 cpl_free(vec_name) ;
1385 cpl_error_reset() ;
1386 *vec = NULL ;
1387 return -1 ;
1388 }
1389 cpl_free(vec_name) ;
1390
1391 /* Create the output */
1392 tab_size = cpl_table_get_nrow(tab) ;
1393 *vec = cpl_vector_new(tab_size) ;
1394 pvec = cpl_vector_get_data(*vec) ;
1395 for (i=0 ; i<tab_size ; i++) pvec[i] = ptab[i] ;
1396
1397 return 0 ;
1398}
1399
1400/*----------------------------------------------------------------------------*/
1437/*----------------------------------------------------------------------------*/
1438int cr2res_extract_slitdec_curved(
1439 const hdrl_image * img_hdrl,
1440 const cpl_table * trace_tab,
1441 const cpl_vector * slit_func_vec_in,
1442 int order,
1443 int trace_id,
1444 int height,
1445 int swath,
1446 int oversample,
1447 double pclip,
1448 double smooth_slit,
1449 double smooth_spec,
1450 int niter,
1451 double kappa,
1452 double error_factor,
1453 cpl_vector ** slit_func,
1454 cpl_bivector ** spec,
1455 hdrl_image ** model)
1456{
1457 double * ycen_rest;
1458 double * ycen_sw;
1459 int * ycen_offset_sw;
1460 double * slitfu_sw_data;
1461 double * model_sw;
1462 const double * slit_func_in;
1463 int * mask_sw;
1464 const cpl_image * img_in;
1465 const cpl_image * err_in;
1466 cpl_image * img_sw;
1467 cpl_image * err_sw;
1468 cpl_image * img_rect;
1469 cpl_image * err_rect;
1470 cpl_image * model_rect;
1471 cpl_vector * ycen ;
1472 cpl_image * img_out;
1473 cpl_vector * slitfu_sw;
1474 cpl_vector * unc_sw;
1475 cpl_vector * spc;
1476 cpl_vector * slitfu;
1477 cpl_vector * weights_sw;
1478 cpl_vector * tmp_vec;
1479 cpl_vector * bins_begin;
1480 cpl_vector * bins_end;
1481 cpl_vector * unc_decomposition;
1482 cpl_size lenx, leny;
1483 cpl_type imtyp;
1484 cpl_polynomial *slitcurve_A, *slitcurve_B, *slitcurve_C;
1485 double * slitcurve_sw;
1486 hdrl_image * model_out;
1487 cpl_bivector * spectrum_loc;
1488 double * sP_old;
1489 double * l_Aij;
1490 double * p_Aij;
1491 double * l_bj;
1492 double * p_bj;
1493 double * zw;
1494 int * zk;
1495 zeta_ref * zeta;
1496 int * m_zeta;
1497 zeta_rng * z_rng;
1498 char * path;
1499 double pixval, errval;
1500 double trace_cen, trace_height;
1501 int i, j, k, nswaths, col, x, y, ny_os,
1502 badpix, delta_x;
1503 int ny, nx;
1504
1505
1506 /* Check Entries */
1507 if (img_hdrl == NULL || trace_tab == NULL) return -1 ;
1508
1509 if (smooth_slit == 0.0) {
1510 cpl_msg_error(__func__, "Slit-smoothing cannot be 0.0");
1511 return -1;
1512 } else if (smooth_slit < 0.1) {
1513 cpl_msg_warning(__func__, "Slit-smoothing unreasonably small");
1514 } else if (smooth_slit > 100.0) {
1515 cpl_msg_warning(__func__, "Slit-smoothing unreasonably big");
1516 }
1517 if (oversample < 3){
1518 cpl_msg_error(__func__, "Oversampling too small");
1519 return -1;
1520 } else if (oversample > 15) {
1521 cpl_msg_warning(__func__, "Large oversampling, runtime will be long");
1522 }
1523 if (niter < 5){
1524 cpl_msg_warning(__func__,
1525 "Allowing at least 5 iterations is recommended.");
1526 }
1527 if (kappa < 4){
1528 cpl_msg_warning(__func__,
1529 "Rejecting outliers < 4 sigma risks making good data.");
1530 }
1531 img_in = hdrl_image_get_image_const(img_hdrl);
1532 err_in = hdrl_image_get_error_const(img_hdrl);
1533
1534 /* Initialise */
1535 imtyp = cpl_image_get_type(img_in);
1536 lenx = cpl_image_get_size_x(img_in);
1537 leny = cpl_image_get_size_y(img_in);
1538
1539 /* Compute height if not given */
1540 if (height <= 0) {
1541 height = cr2res_trace_get_height(trace_tab, order, trace_id);
1542 if (height <= 0) {
1543 cpl_msg_error(__func__, "Cannot compute height");
1544 return -1;
1545 }
1546 }
1547 if (height > leny){
1548 height = leny;
1549 cpl_msg_warning(__func__,
1550 "Given height larger than image, clipping height");
1551 }
1552
1553 /* Get ycen */
1554 if ((ycen = cr2res_trace_get_ycen(trace_tab, order,
1555 trace_id, lenx)) == NULL) {
1556 cpl_msg_error(__func__, "Cannot get ycen");
1557 return -1 ;
1558 }
1559 trace_cen = cpl_vector_get(ycen, cpl_vector_get_size(ycen)/2) ;
1560 trace_height = (double)cr2res_trace_get_height(trace_tab, order, trace_id) ;
1561 cpl_msg_info(__func__, "Y position of the trace: %g -> %g",
1562 trace_cen-(trace_height/2), trace_cen+(trace_height/2)) ;
1563 if (trace_cen-(height/2) < 0.0 ||
1564 trace_cen+(height/2) > CR2RES_DETECTOR_SIZE) {
1565 cpl_msg_error(__func__, "Extraction outside detector edges impossible");
1566 cpl_vector_delete(ycen);
1567 return -1;
1568 }
1569
1570 // Get cut-out rectified order
1571 img_rect = cr2res_image_cut_rectify(img_in, ycen, height);
1572 if (img_rect == NULL){
1573 cpl_msg_error(__func__, "Cannot rectify order");
1574 cpl_vector_delete(ycen);
1575 return -1;
1576 }
1577 if (cpl_msg_get_level() == CPL_MSG_DEBUG) {
1578 cpl_image_save(img_rect, "debug_rectorder_curved.fits", imtyp,
1579 NULL, CPL_IO_CREATE);
1580 }
1581 err_rect = cr2res_image_cut_rectify(err_in, ycen, height);
1582 ycen_rest = cr2res_vector_get_rest(ycen);
1583
1584 /* Retrieve the polynomials that describe the slit tilt and curvature*/
1585 slitcurve_A = cr2res_get_trace_wave_poly(trace_tab, CR2RES_COL_SLIT_CURV_A,
1586 order, trace_id);
1587 slitcurve_B = cr2res_get_trace_wave_poly(trace_tab, CR2RES_COL_SLIT_CURV_B,
1588 order, trace_id);
1589 slitcurve_C = cr2res_get_trace_wave_poly(trace_tab, CR2RES_COL_SLIT_CURV_C,
1590 order, trace_id);
1591 if ((slitcurve_A == NULL) || (slitcurve_B == NULL) || (slitcurve_C == NULL))
1592 {
1593 cpl_msg_error(__func__,
1594 "No (or incomplete) slitcurve data found in trace table");
1595 cpl_vector_delete(ycen);
1596 cpl_free(ycen_rest) ;
1597 cpl_image_delete(err_rect) ;
1598 cpl_image_delete(img_rect) ;
1599 cpl_polynomial_delete(slitcurve_A);
1600 cpl_polynomial_delete(slitcurve_B);
1601 cpl_polynomial_delete(slitcurve_C);
1602 return -1;
1603 }
1604
1605 /* Maximum horizontal shift in detector pixels due to slit image curv. */
1606 delta_x=0;
1607 for (i=1; i<=lenx; i+=swath/2){
1608 double delta_tmp, a, b, c, yc;
1609 /* Do a coarse sweep through the order and evaluate the slitcurve */
1610 /* polynomials at +- height/2, update the value. */
1611 /* Note: The index i is subtracted from a because the polys have */
1612 /* their origin at the edge of the full frame */
1613 //a = cpl_polynomial_eval_1d(slitcurve_A, i, NULL); this is ignored apparently?
1614 b = cpl_polynomial_eval_1d(slitcurve_B, i, NULL);
1615 c = cpl_polynomial_eval_1d(slitcurve_C, i, NULL);
1616 yc = cpl_vector_get(ycen, i-1);
1617
1618 // Shift polynomial to local frame
1619 // We fix a to 0, see comment later on, when we create the
1620 // polynomials for the extraction
1621 a = 0;
1622 b += 2 * yc * c;
1623
1624 delta_tmp = max( fabs(a + (c*height/2. + b)*height/2.),
1625 fabs(a + (c*height/-2. + b)*height/-2.));
1626 if (delta_tmp > delta_x) delta_x = (int)ceil(delta_tmp);
1627 }
1628 delta_x += 1;
1629 cpl_msg_debug(__func__, "Max delta_x from slit curv: %d pix.", delta_x);
1630
1631 if (delta_x >= swath / 4){
1632 cpl_msg_error(__func__,
1633 "Curvature is larger than the swath, try again with a larger swath size");
1634 cpl_vector_delete(ycen);
1635 cpl_free(ycen_rest) ;
1636 cpl_image_delete(err_rect) ;
1637 cpl_image_delete(img_rect) ;
1638 cpl_polynomial_delete(slitcurve_A);
1639 cpl_polynomial_delete(slitcurve_B);
1640 cpl_polynomial_delete(slitcurve_C);
1641 return -1;
1642 }
1643
1644 /* Number of rows after oversampling */
1645 ny_os = oversample*(height+1) +1;
1646 if ((swath = cr2res_extract_slitdec_adjust_swath(ycen, height, leny, swath,
1647 lenx, delta_x, &bins_begin, &bins_end)) == -1){
1648 cpl_msg_error(__func__, "Cannot calculate swath size");
1649 cpl_vector_delete(ycen);
1650 cpl_free(ycen_rest) ;
1651 cpl_image_delete(err_rect) ;
1652 cpl_image_delete(img_rect) ;
1653 cpl_polynomial_delete(slitcurve_A);
1654 cpl_polynomial_delete(slitcurve_B);
1655 cpl_polynomial_delete(slitcurve_C);
1656 return -1;
1657 }
1658 nswaths = cpl_vector_get_size(bins_begin);
1659
1660 /* Use existing slitfunction if given */
1661 slit_func_in = NULL;
1662 if (slit_func_vec_in != NULL) {
1663 cpl_size size;
1664 size = cpl_vector_get_size(slit_func_vec_in);
1665 if (size == ny_os){
1666 slit_func_in = cpl_vector_get_data_const(slit_func_vec_in);
1667 } else {
1668 cpl_msg_warning(__func__, "Ignoring the given slit_func since it is"
1669 " of the wrong size, expected %i but got %lli points.",
1670 ny_os, size);
1671 }
1672 }
1673
1674 /* Allocate */
1675 mask_sw = cpl_malloc(height * swath*sizeof(int));
1676 model_sw = cpl_malloc(height * swath*sizeof(double));
1677 unc_sw = cpl_vector_new(swath);
1678 img_sw = cpl_image_new(swath, height, CPL_TYPE_DOUBLE);
1679 err_sw = cpl_image_new(swath, height, CPL_TYPE_DOUBLE);
1680 ycen_sw = cpl_malloc(swath*sizeof(double));
1681 ycen_offset_sw = cpl_malloc(swath * sizeof(int));
1682
1683 slitcurve_sw = cpl_malloc(swath * 6 * sizeof(double));
1684
1685 // Local versions of return data
1686 slitfu = cpl_vector_new(ny_os);
1687 spectrum_loc = cpl_bivector_new(lenx);
1688 spc = cpl_bivector_get_x(spectrum_loc);
1689 unc_decomposition = cpl_bivector_get_y(spectrum_loc);
1690 for (j=0; j<lenx ; j++){
1691 cpl_vector_set(spc, j, 0.);
1692 cpl_vector_set(unc_decomposition, j, 0.);
1693 }
1694 model_out = hdrl_image_new(lenx, leny);
1695 img_out = hdrl_image_get_image(model_out);
1696 model_rect = cpl_image_new(lenx, height, CPL_TYPE_DOUBLE);
1697
1698 // Work vectors
1699 slitfu_sw = cpl_vector_new(ny_os);
1700 for (j=0; j < ny_os; j++) cpl_vector_set(slitfu_sw, j, 0);
1701 slitfu_sw_data = cpl_vector_get_data(slitfu_sw);
1702 weights_sw = cpl_vector_new(swath);
1703 for (i = 0; i < swath; i++) cpl_vector_set(weights_sw, i, 0);
1704
1705 /* Pre-calculate the weights for overlapping swaths*/
1706 for (i=delta_x; i < swath/2; i++) {
1707 j = i - delta_x + 1;
1708 cpl_vector_set(weights_sw, i, j);
1709 cpl_vector_set(weights_sw, swath - i - 1, j);
1710 }
1711 // normalize such that max(w)=1
1712 cpl_vector_divide_scalar(weights_sw, swath/2 - delta_x + 1);
1713
1714 // assert cpl_vector_get_sum(weights_sw) == swath / 2 - delta_x
1715 // Assign memory for extract_curved algorithm
1716 // Since the arrays always have the same size, we can reuse allocated memory
1717 ny = oversample * (height + 1) + 1;
1718 nx = 4 * delta_x + 1;
1719 if(nx < 3) nx = 3;
1720
1721 sP_old = cpl_malloc(swath * sizeof(double));
1722 l_Aij = cpl_malloc(ny * (4*oversample+1) * sizeof(double));
1723 p_Aij = cpl_malloc(swath * nx * sizeof(double));
1724 l_bj = cpl_malloc(ny * sizeof(double));
1725 p_bj = cpl_malloc(swath * sizeof(double));
1726
1727 /* Scratch buffers for per-pixel merged zeta weights: large enough for
1728 both the slit-function window (2*oversample+1) and the spectrum
1729 window (2*delta_x+1 <= nx) */
1730 i = 3 * (oversample + 1);
1731 if (i < nx) i = nx;
1732 zw = cpl_malloc(i * sizeof(double));
1733 zk = cpl_malloc(i * sizeof(int));
1734
1735 /* Convolution tensor telling the coordinates of subpixels {x, iy}
1736 contributing to detector pixel {x, y}. [ncols][nrows][3*(osample+1)]
1737 */
1738 zeta = cpl_malloc(swath * height * 3 * (oversample + 1)
1739 * sizeof(zeta_ref));
1740
1741 /* The actual number of contributing elements in zeta [ncols][nrows]
1742 and the key ranges of each pixel's zeta list */
1743 m_zeta = cpl_malloc(swath * height * sizeof(int));
1744 z_rng = cpl_malloc(swath * height * sizeof(zeta_rng));
1745
1746 for (i = 0; i < nswaths; i++) {
1747 double *img_sw_data;
1748 double *err_sw_data;
1749 double *spec_sw_data;
1750 double *unc_sw_data;
1751
1752 cpl_image *img_tmp;
1753 cpl_vector *spec_sw;
1754 cpl_vector *spec_tmp;
1755
1756 int sw_start, sw_end, y_lower_limit;
1757
1758 double img_sum;
1759
1760 sw_start = cpl_vector_get(bins_begin, i);
1761 sw_end = cpl_vector_get(bins_end, i);
1762
1763 /* Prepare swath cut-outs and auxiliary data */
1764 for(col=1; col<=swath; col++){ // col is x-index in swath
1765 x = sw_start + col; // coords in large image
1766
1767 /* prepare signal, error and mask */
1768 for(y=1;y<=height;y++){
1769 errval = cpl_image_get(err_rect, x, y, &badpix);
1770 if (isnan(errval) || badpix){
1771 // default to errval of 1 instead of 0
1772 // this avoids division by 0
1773 errval = 1;
1774 }
1775 pixval = cpl_image_get(img_rect, x, y, &badpix);
1776 if (isnan(pixval) || badpix){
1777 // We set bad pixels to neg. infinity, to make sure they are
1778 // rejected in the extraction
1779 // The algorithm does not like NANs!
1780 badpix = 1;
1781 pixval = -DBL_MAX;
1782 errval = 1;
1783 }
1784 cpl_image_set(img_sw, col, y, pixval);
1785 cpl_image_set(err_sw, col, y, errval);
1786 if (badpix){
1787 // Reject the pixel here, so it is not used for the initial
1788 // guess of the spectrum
1789 cpl_image_reject(img_sw, col, y);
1790 }
1791
1792 // raw index for mask, start with 0!
1793 j = (y-1)*swath + (col-1) ;
1794 // The mask value is inverted for the extraction
1795 // 1 for good pixel and 0 for bad pixel
1796 mask_sw[j] = !badpix;
1797 }
1798
1799 /* set slit curvature coefficients, shifted to the local frame */
1800 // The slit curvature has been determined in the global reference
1801 // frame, with the a coefficient set to 0 in the local frame.
1802 // Shifting the polynomial by ycen moves it into the local frame
1803 // again and should result in a = 0:
1804 // a - x + yc * b + yc * yc * c
1805 // However this only works, as long as ycen
1806 // is the same ycen that was used for the slitcurvature. If e.g. we
1807 // switch traces, then ycen will change and a will be unequal 0.
1808 // in fact a will be the offset due to the curvature between the
1809 // old ycen and the new. This will then cause an offset in the
1810 // pixels used for the extraction, so that all traces will have the
1811 // same spectrum, with no relative offsets.
1812 // Which would be great, if we didn't have an offset in the
1813 // wavelength calibration of the different traces.
1814 // Therefore we force a to be 0 in the local frame regardless of
1815 // ycen. For the extraction we only need the b and c coefficient
1816 // anyways.
1817 // Note that this means, we use the curvature a few pixels offset.
1818 // Usually this is no problem, since it only varies slowly over the
1819 // order.
1820 {
1821 double b, c, yc;
1822 b = cpl_polynomial_eval_1d(slitcurve_B, x, NULL);
1823 c = cpl_polynomial_eval_1d(slitcurve_C, x, NULL);
1824 yc = cpl_vector_get(ycen, x-1);
1825 slitcurve_sw[curve_index(col-1, 0)] = 0;
1826 slitcurve_sw[curve_index(col-1, 1)] = b + 2 * yc * c;
1827 slitcurve_sw[curve_index(col-1, 2)] = c;
1828 slitcurve_sw[curve_index(col-1, 3)] = 0;
1829 slitcurve_sw[curve_index(col-1, 4)] = 0;
1830 slitcurve_sw[curve_index(col-1, 5)] = 0;
1831 }
1832 }
1833
1834 for (j=0; j< height * swath; j++) model_sw[j] = 0;
1835 img_sw_data = cpl_image_get_data_double(img_sw);
1836 err_sw_data = cpl_image_get_data_double(err_sw);
1837 unc_sw_data = cpl_vector_get_data(unc_sw);
1838 // First guess for the spectrum
1839 // img_tmp = cpl_image_collapse_median_create(img_sw, 0, 0, 0);
1840 img_tmp = cpl_image_collapse_median_create(img_sw, 0, 0, 0);
1841 spec_tmp = cpl_vector_new_from_image_row(img_tmp, 1);
1842 cpl_vector_multiply_scalar(spec_tmp,
1843 (double)cpl_image_get_size_y(img_sw));
1844 spec_sw = cpl_vector_filter_median_create(spec_tmp, 1);
1845 cpl_vector_delete(spec_tmp);
1846 cpl_image_delete(img_tmp);
1847 spec_sw_data = cpl_vector_get_data(spec_sw);
1848
1849 for (j=sw_start;j<sw_end;j++){
1850 ycen_sw[j-sw_start] = ycen_rest[j];
1851 ycen_offset_sw[j-sw_start] = (int) cpl_vector_get(ycen, j);
1852 }
1853 y_lower_limit = height / 2;
1854
1855 img_tmp = cpl_image_wrap_int(swath, height, mask_sw);
1856 img_sum = cpl_image_get_flux(img_tmp);
1857 if (img_sum < 0.5 * swath*height){
1858 cpl_msg_error(__func__,
1859 "Only %.0f %% of pixels not masked, cannot extract",
1860 100*img_sum/(swath*height));
1861 cpl_image_unwrap(img_tmp);
1862 cpl_vector_delete(spec_sw);
1863 break;
1864 }
1865 if (cpl_msg_get_level() == CPL_MSG_DEBUG)
1866 {
1867 cpl_image_save(img_tmp, "debug_mask_before_sw.fits", CPL_TYPE_INT, NULL, CPL_IO_CREATE);
1868 cpl_vector_save(spec_sw, "debug_spc_initial_guess.fits",
1869 CPL_TYPE_DOUBLE, NULL, CPL_IO_CREATE);
1870 }
1871 cpl_image_unwrap(img_tmp);
1872
1873 /* Finally ready to call the slit-decomp */
1874 cr2res_extract_slit_func_curved(error_factor, swath, height, oversample, pclip,
1875 img_sw_data, err_sw_data, mask_sw, ycen_sw, ycen_offset_sw,
1876 y_lower_limit, slitcurve_sw, delta_x, slitfu_sw_data,
1877 spec_sw_data, model_sw, unc_sw_data, smooth_spec, smooth_slit,
1878 5.e-5, niter, kappa, slit_func_in, sP_old, l_Aij, p_Aij, l_bj,
1879 p_bj, zw, zk, zeta, m_zeta, z_rng);
1880
1881 // add up slit-functions, divide by nswaths below to get average
1882 if (i==0) cpl_vector_copy(slitfu,slitfu_sw);
1883 else cpl_vector_add(slitfu,slitfu_sw);
1884
1885 if (cpl_msg_get_level() == CPL_MSG_DEBUG) {
1886 path = cpl_sprintf("debug_spc_%i.fits", i);
1887 cpl_vector_save(spec_sw, path , CPL_TYPE_DOUBLE, NULL,
1888 CPL_IO_CREATE);
1889 cpl_free(path);
1890
1891 path = cpl_sprintf("debug_mask_%i.fits", i);
1892 img_tmp = cpl_image_wrap_int(swath, height, mask_sw);
1893 cpl_image_save(img_tmp, path, CPL_TYPE_INT, NULL, CPL_IO_CREATE);
1894 cpl_free(path);
1895 cpl_image_unwrap(img_tmp);
1896
1897 tmp_vec = cpl_vector_wrap(swath, ycen_sw);
1898 path = cpl_sprintf("debug_ycen_%i.fits", i);
1899 cpl_vector_save(tmp_vec, path, CPL_TYPE_DOUBLE, NULL,
1900 CPL_IO_CREATE);
1901 cpl_vector_unwrap(tmp_vec);
1902 cpl_free(path);
1903
1904 cpl_vector_save(weights_sw, "debug_weights.fits", CPL_TYPE_DOUBLE,
1905 NULL, CPL_IO_CREATE);
1906 path = cpl_sprintf("debug_slitfu_%i.fits", i);
1907 cpl_vector_save(slitfu_sw, path, CPL_TYPE_DOUBLE,
1908 NULL, CPL_IO_CREATE);
1909 cpl_free(path);
1910
1911 path = cpl_sprintf("debug_model_%i.fits", i);
1912 img_tmp = cpl_image_wrap_double(swath, height, model_sw);
1913 cpl_image_save(img_tmp, path, CPL_TYPE_DOUBLE,
1914 NULL, CPL_IO_CREATE);
1915 cpl_image_unwrap(img_tmp);
1916 cpl_free(path);
1917
1918 path = cpl_sprintf("debug_img_sw_%i.fits", i);
1919 cpl_image_save(img_sw, path, CPL_TYPE_DOUBLE, NULL,
1920 CPL_IO_CREATE);
1921 cpl_free(path);
1922 }
1923
1924 // The last bins are shifted, overwriting the first k values
1925 // this is the same amount the bin was shifted to the front before
1926 // (when creating the bins)
1927 // The duplicate values in the vector will not matter as they are
1928 // not used below
1929 if ((i == nswaths - 1) && (i != 0)){
1930 k = cpl_vector_get(bins_end, i-1) -
1931 cpl_vector_get(bins_begin, i) - swath / 2 - delta_x;
1932
1933 for (j = 0; j < swath - k; j++){
1934 cpl_vector_set(spec_sw, j, cpl_vector_get(spec_sw, j + k));
1935 cpl_vector_set(unc_sw, j, cpl_vector_get(unc_sw, j + k));
1936 for (y = 0; y < height; y++)
1937 model_sw[y * swath + j] = model_sw[y * swath + j + k];
1938 }
1939 sw_start = cpl_vector_get(bins_begin, i-1) + swath / 2 - delta_x;
1940 cpl_vector_set(bins_begin, i, sw_start);
1941 // for the following k's
1942 cpl_vector_set(bins_end, i, lenx);
1943 }
1944
1945 if (nswaths==1) ; // no weighting if only one swath
1946 else if (i==0){ // first and last half swath are not weighted
1947 for (j = 0; j < delta_x; j++)
1948 {
1949 cpl_vector_set(spec_sw, j, 0);
1950 cpl_vector_set(unc_sw, j, 0.);
1951 for (y = 0; y < height; y++) model_sw[y * swath + j] = 0;
1952 }
1953 for (j = swath/2; j < swath; j++) {
1954 cpl_vector_set(spec_sw, j,
1955 cpl_vector_get(spec_sw,j) * cpl_vector_get(weights_sw,j));
1956 cpl_vector_set(unc_sw, j,
1957 cpl_vector_get(unc_sw, j) * cpl_vector_get(weights_sw, j));
1958 for (y = 0; y < height; y++) {
1959 model_sw[y * swath + j] *= cpl_vector_get(weights_sw, j);
1960 }
1961 }
1962 } else if (i == nswaths - 1) {
1963 for (j = sw_end-sw_start-1; j >= sw_end-sw_start-delta_x-1; j--)
1964 {
1965 cpl_vector_set(spec_sw, j, 0);
1966 cpl_vector_set(unc_sw, j, 0);
1967 for (y = 0; y < height; y++) model_sw[y * swath + j] = 0;
1968 }
1969 for (j = 0; j < swath / 2; j++) {
1970 cpl_vector_set(spec_sw, j,
1971 cpl_vector_get(spec_sw,j) * cpl_vector_get(weights_sw,j));
1972 cpl_vector_set(unc_sw, j,
1973 cpl_vector_get(unc_sw,j) * cpl_vector_get(weights_sw,j));
1974 for (y = 0; y < height; y++) {
1975 model_sw[y * swath + j] *= cpl_vector_get(weights_sw,j);
1976 }
1977 }
1978 } else {
1979 /* Multiply by weights and add to output array */
1980 cpl_vector_multiply(spec_sw, weights_sw);
1981 cpl_vector_multiply(unc_sw, weights_sw);
1982 for (y = 0; y < height; y++) {
1983 for (j = 0; j < swath; j++){
1984 model_sw[y * swath + j] *= cpl_vector_get(weights_sw,j);
1985 }
1986 }
1987 }
1988
1989 if (cpl_msg_get_level() == CPL_MSG_DEBUG) {
1990 img_tmp = cpl_image_wrap_double(swath, height, model_sw);
1991 cpl_image_save(img_tmp, "debug_model_after_sw.fits", CPL_TYPE_DOUBLE,
1992 NULL, CPL_IO_CREATE);
1993 cpl_image_unwrap(img_tmp);
1994 }
1995
1996 // Save swath to output vector
1997 for (j=sw_start;j<sw_end;j++) {
1998 cpl_vector_set(spc, j,
1999 cpl_vector_get(spec_sw, j-sw_start) + cpl_vector_get(spc, j));
2000 // just add weighted errors (instead of squared sum)
2001 // as they are not independent
2002 cpl_vector_set(unc_decomposition, j,
2003 cpl_vector_get(unc_sw, j - sw_start)
2004 + cpl_vector_get(unc_decomposition, j));
2005
2006 for(y = 0; y < height; y++){
2007 cpl_image_set(model_rect, j+1, y+1,
2008 cpl_image_get(model_rect, j+1, y+1, &badpix)
2009 + model_sw[y * swath + j - sw_start]);
2010 if (badpix) cpl_image_reject(model_rect, j+1, y+1);
2011 }
2012 }
2013
2014 if (cpl_msg_get_level() == CPL_MSG_DEBUG) {
2015 cpl_image_save(model_rect, "debug_model_after_merge.fits",
2016 CPL_TYPE_DOUBLE, NULL, CPL_IO_CREATE);
2017 }
2018
2019 cpl_vector_delete(spec_sw);
2020 } // End loop over swaths
2021
2022 // divide by nswaths to make the slitfu into the average over all swaths.
2023 cpl_vector_divide_scalar(slitfu, nswaths);
2024
2025 // Deallocate loop memory
2026 cpl_free(sP_old);
2027 cpl_free(l_Aij);
2028 cpl_free(p_Aij);
2029 cpl_free(l_bj);
2030 cpl_free(p_bj);
2031 cpl_free(zw);
2032 cpl_free(zk);
2033
2034 cpl_free(zeta);
2035 cpl_free(m_zeta);
2036 cpl_free(z_rng);
2037
2038 cpl_image_delete(img_rect);
2039 cpl_image_delete(err_rect);
2040 cpl_image_delete(img_sw);
2041 cpl_image_delete(err_sw);
2042
2043 cpl_free(mask_sw);
2044 cpl_free(model_sw);
2045 cpl_vector_delete(unc_sw);
2046 cpl_free(ycen_rest);
2047 cpl_free(ycen_sw);
2048 cpl_free(ycen_offset_sw);
2049
2050 cpl_vector_delete(bins_begin);
2051 cpl_vector_delete(bins_end);
2052 cpl_vector_delete(slitfu_sw);
2053 cpl_vector_delete(weights_sw);
2054
2055 cpl_polynomial_delete(slitcurve_A);
2056 cpl_polynomial_delete(slitcurve_B);
2057 cpl_polynomial_delete(slitcurve_C);
2058 cpl_free(slitcurve_sw);
2059
2060 // insert model_rect into large frame
2061 if (cr2res_image_insert_rect(model_rect, ycen, img_out) == -1) {
2062 // Cancel
2063 cpl_msg_error(__func__, "failed to reinsert model swath into model image");
2064 cpl_image_delete(model_rect);
2065 hdrl_image_delete(model_out);
2066 cpl_vector_delete(ycen);
2067 cpl_bivector_delete(spectrum_loc);
2068 cpl_vector_delete(slitfu);
2069 return -1;
2070 }
2071
2072 if (cpl_msg_get_level() == CPL_MSG_DEBUG) {
2073 cpl_image_save(model_rect, "debug_model_rect.fits", CPL_TYPE_DOUBLE,
2074 NULL, CPL_IO_CREATE);
2075 cpl_image_save(img_out, "debug_model_all.fits", CPL_TYPE_DOUBLE,
2076 NULL, CPL_IO_CREATE);
2077 cpl_vector_save(spc, "debug_spc_all.fits", CPL_TYPE_DOUBLE,
2078 NULL, CPL_IO_CREATE);
2079 }
2080
2081 cpl_image_delete(model_rect);
2082 cpl_vector_delete(ycen);
2083
2084 if (cpl_error_get_code() != CPL_ERROR_NONE){
2085 cpl_msg_error(__func__,
2086 "Something went wrong in the extraction. Error Code: %i, loc: %s",
2087 cpl_error_get_code(), cpl_error_get_where());
2088 cpl_error_reset();
2089 cpl_vector_delete(slitfu);
2090 cpl_bivector_delete(spectrum_loc);
2091 hdrl_image_delete(model_out);
2092 return -1;
2093 }
2094
2095 *slit_func = slitfu;
2096 *spec = spectrum_loc;
2097 *model = model_out;
2098 return 0;
2099}
2100
2101
2102/*-------------------------------------------------------------------------*/
2103/*-------------------- EXTRACT 2d ------------------------------*/
2104/*-------------------------------------------------------------------------*/
2105
2106/*-------------------------------------------------------------------------*/
2118/*--------------------------------------------------------------------------*/
2119int cr2res_extract2d_traces(
2120 const hdrl_image * img,
2121 const cpl_table * traces,
2122 int reduce_order,
2123 int reduce_trace,
2124 cpl_table ** extracted)
2125{
2126 cpl_bivector ** spectrum ;
2127 cpl_bivector ** position ;
2128 cpl_vector ** wavelength ;
2129 cpl_vector ** slit_fraction ;
2130 cpl_table * extract_loc ;
2131 cpl_image * wavemap;
2132 cpl_image * slitmap;
2133 int nb_traces, i, npoints ;
2134
2135 /* Check Entries */
2136 if (img == NULL || traces == NULL) return -1 ;
2137
2138 /* Initialise */
2139 nb_traces = cpl_table_get_nrow(traces) ;
2140 npoints = CR2RES_DETECTOR_SIZE * CR2RES_DETECTOR_SIZE / nb_traces;
2141
2142 /* Allocate Data containers */
2143 spectrum = cpl_malloc(nb_traces * sizeof(cpl_bivector *)) ;
2144 position = cpl_malloc(nb_traces * sizeof(cpl_bivector *)) ;
2145 wavelength = cpl_malloc(nb_traces * sizeof(cpl_vector *)) ;
2146 slit_fraction = cpl_malloc(nb_traces * sizeof(cpl_vector *)) ;
2147
2148 // Calculate wavelength and slitfunction map once
2149 if (cr2res_slit_pos_image(traces, &slitmap, &wavemap) != 0)
2150 {
2151 cpl_msg_error(__func__,
2152 "Could not create wavelength / slit_fraction image");
2153 cpl_free(spectrum);
2154 cpl_free(position);
2155 cpl_free(wavelength);
2156 cpl_free(slit_fraction);
2157 return -1;
2158 }
2159
2160 /* Loop over the traces and extract them */
2161 for (i=0 ; i<nb_traces ; i++) {
2162 /* Initialise */
2163 spectrum[i] = NULL ;
2164 position[i] = NULL ;
2165 wavelength[i] = NULL ;
2166 slit_fraction[i] = NULL ;
2167
2168 int order, trace_id;
2169
2170 /* Get Order and trace id */
2171 order = cpl_table_get(traces, CR2RES_COL_ORDER, i, NULL) ;
2172 trace_id = cpl_table_get(traces, CR2RES_COL_TRACENB, i, NULL) ;
2173
2174 /* Check if this order needs to be skipped */
2175 if (reduce_order > -1 && order != reduce_order) continue ;
2176
2177 /* Check if this trace needs to be skipped */
2178 if (reduce_trace > -1 && trace_id != reduce_trace) continue ;
2179
2180 cpl_msg_info(__func__, "Process Order %d/Trace %d",order,trace_id) ;
2181 cpl_msg_indent_more() ;
2182
2183 /* Call the Extraction */
2184 if (cr2res_extract2d_trace(img, traces, order, trace_id,
2185 npoints, wavemap, slitmap,
2186 &(spectrum[i]), &(position[i]), &(wavelength[i]),
2187 &(slit_fraction[i])) != 0) {
2188 cpl_msg_error(__func__, "Cannot extract2d the trace") ;
2189 spectrum[i] = NULL ;
2190 position[i] = NULL ;
2191 wavelength[i] = NULL ;
2192 slit_fraction[i] = NULL ;
2193 cpl_error_reset() ;
2194 cpl_msg_indent_less() ;
2195 continue ;
2196 }
2197 cpl_msg_indent_less() ;
2198 }
2199
2200 /* Create the extracted_tab for the current detector */
2201 extract_loc = cr2res_extract_EXTRACT2D_create(spectrum, position,
2202 wavelength, slit_fraction, traces) ;
2203
2204 /* Deallocate Vectors */
2205 for (i=0 ; i<nb_traces ; i++) {
2206 if (spectrum[i] != NULL) cpl_bivector_delete(spectrum[i]) ;
2207 if (position[i] != NULL) cpl_bivector_delete(position[i]) ;
2208 if (wavelength[i] != NULL) cpl_vector_delete(wavelength[i]) ;
2209 if (slit_fraction[i] != NULL) cpl_vector_delete(slit_fraction[i]) ;
2210 }
2211 cpl_free(spectrum) ;
2212 cpl_free(position) ;
2213 cpl_free(wavelength) ;
2214 cpl_free(slit_fraction) ;
2215 cpl_image_delete(wavemap);
2216 cpl_image_delete(slitmap);
2217
2218 /* Return */
2219 *extracted = extract_loc ;
2220 return 0 ;
2221}
2222
2223/*----------------------------------------------------------------------------*/
2242/*----------------------------------------------------------------------------*/
2243int cr2res_extract2d_trace(
2244 const hdrl_image * in,
2245 const cpl_table * trace_tab,
2246 int order,
2247 int trace_id,
2248 int npoints,
2249 const cpl_image * wavemap,
2250 const cpl_image * slitmap,
2251 cpl_bivector ** spectrum,
2252 cpl_bivector ** position,
2253 cpl_vector ** wavelength,
2254 cpl_vector ** slit_fraction)
2255{
2256 cpl_bivector * spectrum_local ;
2257 cpl_vector * spectrum_flux ;
2258 cpl_vector * spectrum_error ;
2259 cpl_bivector * position_local ;
2260 cpl_vector * position_x;
2261 cpl_vector * position_y;
2262 cpl_vector * wavelength_local ;
2263 cpl_vector * slit_fraction_local ;
2264 const cpl_array * lower_array;
2265 const cpl_array * upper_array;
2266 cpl_polynomial * lower_poly;
2267 cpl_polynomial * upper_poly;
2268 cpl_vector * lower;
2269 cpl_vector * upper;
2270 cpl_vector * x;
2271
2272 int k, bad_pix;
2273 cpl_size i, j, row;
2274 double flux, err;
2275
2276 /* Check Entries */
2277 if (in==NULL || trace_tab==NULL || spectrum==NULL || position==NULL
2278 || wavelength==NULL || slit_fraction==NULL) return -1 ;
2279
2280 if ((k = cr2res_get_trace_table_index(trace_tab, order, trace_id)) == -1)
2281 {
2282 cpl_msg_error(__func__, "Order and/or Trace not found in trace table");
2283 return -1;
2284 }
2285
2286 // Step 0: Initialise output arrays
2287 wavelength_local = cpl_vector_new(npoints);
2288 slit_fraction_local = cpl_vector_new(npoints);
2289 position_x = cpl_vector_new(npoints);
2290 position_y = cpl_vector_new(npoints);
2291 spectrum_flux = cpl_vector_new(npoints);
2292 spectrum_error = cpl_vector_new(npoints);
2293
2294 // Step 1: Figure out pixels in the current trace
2295 // i.e. everything between upper and lower in trace_wave
2296
2297 // The x coordinate of the detector
2298 x = cpl_vector_new(CR2RES_DETECTOR_SIZE);
2299 for (i = 0; i < CR2RES_DETECTOR_SIZE; i++) cpl_vector_set(x, i, i + 1);
2300
2301 lower_array = cpl_table_get_array(trace_tab, CR2RES_COL_LOWER, k);
2302 upper_array = cpl_table_get_array(trace_tab, CR2RES_COL_UPPER, k);
2303 lower_poly = cr2res_convert_array_to_poly(lower_array);
2304 upper_poly = cr2res_convert_array_to_poly(upper_array);
2305 lower = cr2res_polynomial_eval_vector(lower_poly, x);
2306 upper = cr2res_polynomial_eval_vector(upper_poly, x);
2307
2308 // Step 2: Iterate over pixels in the given trace
2309 // and fill the vectors
2310 row = -1;
2311 for (i = 0; i < CR2RES_DETECTOR_SIZE; i++)
2312 {
2313 for (j = cpl_vector_get(lower, i); j < cpl_vector_get(upper, i); j++)
2314 {
2315 /* Protect the case where the trace goes out of the det */
2316 if (j<1 || j>CR2RES_DETECTOR_SIZE) continue ;
2317 row++;
2318 cpl_vector_set(position_x, row, i +1);
2319 cpl_vector_set(position_y, row, j);
2320 flux = cpl_image_get(hdrl_image_get_image_const(in), i + 1, j,
2321 &bad_pix);
2322 err = cpl_image_get(hdrl_image_get_error_const(in), i + 1, j,
2323 &bad_pix);
2324 cpl_vector_set(spectrum_flux, row, flux);
2325 cpl_vector_set(spectrum_error, row, err);
2326
2327 // Set wavelength
2328 cpl_vector_set(wavelength_local, row, cpl_image_get(
2329 wavemap, i + 1, j, &bad_pix));
2330 // Set Slitfraction
2331 cpl_vector_set(slit_fraction_local, row, cpl_image_get(
2332 slitmap, i + 1, j, &bad_pix));
2333 }
2334 }
2335
2336 // Step 3: resize output
2337 for (i = row; i < npoints; i++){
2338 cpl_vector_set(position_x, i, NAN);
2339 cpl_vector_set(position_y, i, NAN);
2340 cpl_vector_set(spectrum_flux, i, NAN);
2341 cpl_vector_set(spectrum_error, i, NAN);
2342 cpl_vector_set(wavelength_local, i, NAN);
2343 cpl_vector_set(slit_fraction_local, i, NAN);
2344 }
2345
2346 position_local = cpl_bivector_wrap_vectors(position_x, position_y);
2347 spectrum_local = cpl_bivector_wrap_vectors(spectrum_flux, spectrum_error);
2348
2349 cpl_polynomial_delete(upper_poly);
2350 cpl_polynomial_delete(lower_poly);
2351 cpl_vector_delete(upper);
2352 cpl_vector_delete(lower);
2353 cpl_vector_delete(x);
2354
2355 *spectrum = spectrum_local ;
2356 *position = position_local ;
2357 *wavelength = wavelength_local ;
2358 *slit_fraction = slit_fraction_local ;
2359 return 0;
2360}
2361
2362/*----------------------------------------------------------------------------*/
2372/*----------------------------------------------------------------------------*/
2373cpl_table * cr2res_extract_EXTRACT2D_create(
2374 cpl_bivector ** spectrum,
2375 cpl_bivector ** position,
2376 cpl_vector ** wavelength,
2377 cpl_vector ** slit_fraction,
2378 const cpl_table * trace_table)
2379{
2380 cpl_table * out ;
2381 char * col_name ;
2382 const double * pspec ;
2383 const double * perr ;
2384 const double * pxposition ;
2385 const double * pyposition ;
2386 const double * pwave ;
2387 const double * pslit_frac ;
2388 int nrows, all_null, i, order, trace_id, nb_traces ;
2389
2390 /* Check entries */
2391 if (spectrum==NULL || trace_table==NULL || position==NULL ||
2392 wavelength==NULL || slit_fraction==NULL)
2393 return NULL ;
2394
2395 /* Initialise */
2396 nb_traces = cpl_table_get_nrow(trace_table) ;
2397
2398 /* Check if all vectors are not null */
2399 all_null = 1 ;
2400 for (i=0 ; i<nb_traces ; i++)
2401 if (spectrum[i] != NULL) {
2402 nrows = cpl_bivector_get_size(spectrum[i]) ;
2403 all_null = 0 ;
2404 }
2405 if (all_null == 1) return NULL ;
2406
2407 /* Check the sizes */
2408 for (i=0 ; i<nb_traces ; i++)
2409 if (spectrum[i] != NULL && cpl_bivector_get_size(spectrum[i]) != nrows)
2410 return NULL ;
2411
2412 /* Create the table */
2413 out = cpl_table_new(nrows);
2414 for (i=0 ; i<nb_traces ; i++) {
2415 order = cpl_table_get(trace_table, CR2RES_COL_ORDER, i, NULL) ;
2416 trace_id = cpl_table_get(trace_table, CR2RES_COL_TRACENB, i, NULL) ;
2417 /* Create SPEC column */
2418 col_name = cr2res_dfs_SPEC_colname(order, trace_id) ;
2419 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
2420 cpl_free(col_name) ;
2421 /* Create SPEC_ERR column */
2422 col_name = cr2res_dfs_SPEC_ERR_colname(order, trace_id) ;
2423 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
2424 cpl_free(col_name) ;
2425 /* Create WAVELENGTH column */
2426 col_name = cr2res_dfs_WAVELENGTH_colname(order, trace_id) ;
2427 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
2428 cpl_free(col_name) ;
2429 /* Create POSITIONX column */
2430 col_name = cr2res_dfs_POSITIONX_colname(order, trace_id) ;
2431 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
2432 cpl_free(col_name) ;
2433 /* Create POSITIONY column */
2434 col_name = cr2res_dfs_POSITIONY_colname(order, trace_id) ;
2435 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
2436 cpl_free(col_name) ;
2437 /* Create SLIT_FRACTION column */
2438 col_name = cr2res_dfs_SLIT_FRACTION_colname(order, trace_id) ;
2439 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
2440 cpl_free(col_name) ;
2441 }
2442
2443 /* Fill the table */
2444 for (i=0 ; i<nb_traces ; i++) {
2445 if (spectrum[i]!=NULL && position[i]!=NULL &&
2446 wavelength[i]!=NULL && slit_fraction[i]!=NULL) {
2447 order = cpl_table_get(trace_table, CR2RES_COL_ORDER, i, NULL) ;
2448 trace_id = cpl_table_get(trace_table, CR2RES_COL_TRACENB, i, NULL) ;
2449 pspec = cpl_bivector_get_x_data_const(spectrum[i]) ;
2450 perr = cpl_bivector_get_y_data_const(spectrum[i]);
2451 pxposition = cpl_bivector_get_x_data_const(position[i]) ;
2452 pyposition = cpl_bivector_get_y_data_const(position[i]) ;
2453 pwave = cpl_vector_get_data_const(wavelength[i]) ;
2454 pslit_frac = cpl_vector_get_data_const(slit_fraction[i]) ;
2455 /* Fill SPEC column */
2456 col_name = cr2res_dfs_SPEC_colname(order, trace_id) ;
2457 cpl_table_copy_data_double(out, col_name, pspec) ;
2458 cpl_free(col_name) ;
2459 /* Fill SPEC_ERR column */
2460 col_name = cr2res_dfs_SPEC_ERR_colname(order, trace_id) ;
2461 cpl_table_copy_data_double(out, col_name, perr) ;
2462 cpl_free(col_name) ;
2463 /* Fill WAVELENGTH column */
2464 col_name = cr2res_dfs_WAVELENGTH_colname(order, trace_id) ;
2465 cpl_table_copy_data_double(out, col_name, pwave) ;
2466 cpl_free(col_name) ;
2467 /* Fill POSITIONX column */
2468 col_name = cr2res_dfs_POSITIONX_colname(order, trace_id) ;
2469 cpl_table_copy_data_double(out, col_name, pxposition) ;
2470 cpl_free(col_name) ;
2471 /* Fill POSITIONY column */
2472 col_name = cr2res_dfs_POSITIONY_colname(order, trace_id) ;
2473 cpl_table_copy_data_double(out, col_name, pyposition) ;
2474 cpl_free(col_name) ;
2475 /* Fill SLIT_FRACTION column */
2476 col_name = cr2res_dfs_SLIT_FRACTION_colname(order, trace_id) ;
2477 cpl_table_copy_data_double(out, col_name, pslit_frac) ;
2478 cpl_free(col_name) ;
2479 }
2480 }
2481 return out ;
2482}
2483
2485
2486
2487/*----------------------------------------------------------------------------*/
2491/*----------------------------------------------------------------------------*/
2492static inline void cr2res_extract_zeta_add(
2493 zeta_ref * zeta,
2494 int * m_zeta,
2495 zeta_rng * z_rng,
2496 int ncols,
2497 int nrows,
2498 int osample,
2499 int x,
2500 int iy,
2501 int xx,
2502 int yy,
2503 double w)
2504{
2505 if (xx >= 0 && xx < ncols && yy >= 0 && yy < nrows && w > 0)
2506 {
2507 const int m = m_zeta[mzeta_index(xx, yy)];
2508 zeta_rng * zr = &z_rng[mzeta_index(xx, yy)];
2509 zeta[zeta_index(xx, yy, m)].x = x;
2510 zeta[zeta_index(xx, yy, m)].iy = iy;
2511 zeta[zeta_index(xx, yy, m)].w = w;
2512 m_zeta[mzeta_index(xx, yy)]++;
2513 if (iy < zr->min_iy) zr->min_iy = iy;
2514 if (iy > zr->max_iy) zr->max_iy = iy;
2515 if (x < zr->min_x) zr->min_x = x;
2516 if (x > zr->max_x) zr->max_x = x;
2517 }
2518}
2519
2520/*----------------------------------------------------------------------------*/
2552/*----------------------------------------------------------------------------*/
2553static int cr2res_extract_zeta_tensors(
2554 int ncols,
2555 int nrows,
2556 const double * ycen,
2557 const int * ycen_offset,
2558 int y_lower_lim,
2559 int osample,
2560 const double * slitcurve,
2561 zeta_ref * zeta,
2562 int * m_zeta,
2563 zeta_rng * z_rng)
2564{
2565 int x, xx, y, yy, ix1, ix2, iy, iy1, iy2;
2566 double step, delta, dy, w, d1, d2;
2567
2568 step = 1.e0 / osample;
2569
2570 /* Clean zeta counts. The zeta entries themselves need no initialization:
2571 only the first m_zeta[x, y] entries of each list are ever read.
2572 Same for the key ranges: only read where m_zeta[x, y] > 0. */
2573 for (x = 0; x < ncols; x++)
2574 for (y = 0; y < nrows; y++) {
2575 zeta_rng * zr = &z_rng[mzeta_index(x, y)];
2576 m_zeta[mzeta_index(x, y)] = 0;
2577 zr->min_iy = INT_MAX;
2578 zr->max_iy = INT_MIN;
2579 zr->min_x = INT_MAX;
2580 zr->max_x = INT_MIN;
2581 }
2582
2583 /*
2584 Construct the zeta tensor. It contains pixel references and contribution
2585 values coming from subpixels to a given detector pixel.
2586 Note that zeta is used in the equations for sL, sP and for the model but it
2587 does not involve the data, only the geometry. Thus it can be pre-computed once.
2588 */
2589 for (x = 0; x < ncols; x++)
2590 {
2591 /*
2592 I promised to reconsider the initial offset. Here it is. For the original layout
2593 (no column shifts and discontinuities in ycen) there is pixel y that contains the
2594 central line yc. There are two options here (by construction of ycen that can be 0
2595 but cannot be 1): (1) yc is inside pixel y and (2) yc falls at the boundary between
2596 pixels y and y-1. yc cannot be at the boundary of pixels y+1 and y because we would
2597 select y+1 to be pixel y in that case.
2598
2599 Next we need to define starting and ending indices iy for sL subpixels that contribute
2600 to pixel y. I call them iy1 and iy2. For both cases we assume osample+1 subpixels covering
2601 pixel y (weird). So for case 1 iy1 will be (y-1)*osample and iy2 == y*osample. Special
2602 treatment of the boundary subpixels will compensate for introducing extra subpixel in
2603 case 1. In case 2 things are more logical: iy1=(yc-y)*osample+(y-1)*osample;
2604 iy2=(y+1-yc)*osample)+(y-1)*osample. ycen is yc-y making things simpler. Note also that
2605 the same pattern repeats for all rows: we only need to initialize iy1 and iy2 and keep
2606 incrementing them by osample.
2607 */
2608
2609 iy2 = osample - floor(ycen[x] * osample);
2610 iy1 = iy2 - osample;
2611
2612 /*
2613 Handling partial subpixels cut by detector pixel rows is again tricky. Here we have three
2614 cases (mostly because of the decision to assume that we always have osample+1 subpixels
2615 per one detector pixel). Here d1 is the fraction of the subpixel iy1 inside detector pixel y.
2616 d2 is then the fraction of subpixel iy2 inside detector pixel y. By definition d1+d2==step.
2617 Case 1: ycen falls on the top boundary of each detector pixel (ycen == 1). Here we conclude
2618 that the first subpixel is fully contained inside pixel y and d1 is set to step.
2619 Case 2: ycen falls on the bottom boundary of each detector pixel (ycen == 0). Here we conclude
2620 that the first subpixel is totally outside of pixel y and d1 is set to 0.
2621 Case 3: ycen falls inside of each pixel (0>ycen>1). In this case d1 is set to the fraction of
2622 the first step contained inside of each pixel.
2623 And BTW, this also means that central line coincides with the upper boundary of subpixel iy2
2624 when the y loop reaches pixel y_lower_lim. In other words:
2625
2626 dy=(iy-(y_lower_lim+ycen[x])*osample)*step-0.5*step
2627 */
2628
2629 d1 = fmod(ycen[x], step);
2630 if (d1 == 0)
2631 d1 = step;
2632 d2 = step - d1;
2633
2634 /* Define initial distance from ycen */
2635 /* It is given by the center of the first */
2636 /* subpixel falling into pixel y_lower_lim */
2637 dy = ycen[x] - floor((y_lower_lim + ycen[x]) / step) * step - step;
2638
2639 /*
2640 Now we go detector pixels x and y incrementing subpixels looking for their contributions
2641 to the current and adjacent pixels. Note that the curvature/tilt of the projected slit
2642 image could be so large that subpixel iy may not contribute to column x at all. On the
2643 other hand, subpixels around ycen by definition must contribute to pixel x,y.
2644
2645 Each subpixel is assumed to be exactly 1 detector pixel wide; a horizontal shift delta
2646 divides its weight w between columns ix1=int(delta) and ix2=ix1+signum(delta) as
2647 (1-|delta-ix1|)*w and |delta-ix1|*w. The yy offset is required because the iy subpixel
2648 contributes to the yy row in the xx column of detector pixels where yy and y are in the
2649 same row. In the packed array this is not necessarily true. Instead, what we know is:
2650 y+ycen_offset[x] == yy+ycen_offset[xx]
2651 */
2652
2653 for (y = 0; y < nrows; y++)
2654 {
2655 iy1 += osample; // Bottom subpixel falling in row y
2656 iy2 += osample; // Top subpixel falling in row y
2657 dy -= step;
2658 for (iy = iy1; iy <= iy2; iy++)
2659 {
2660 double t;
2661 if (iy == iy1)
2662 w = d1;
2663 else if (iy == iy2)
2664 w = d2;
2665 else
2666 w = step;
2667 dy += step;
2668 t = dy - ycen[x];
2669 delta = t * (slitcurve[curve_index(x, 1)] +
2670 t * (slitcurve[curve_index(x, 2)] +
2671 t * (slitcurve[curve_index(x, 3)] +
2672 t * (slitcurve[curve_index(x, 4)] +
2673 t * slitcurve[curve_index(x, 5)]))));
2674 ix1 = delta;
2675 ix2 = ix1 + signum(delta);
2676
2677 if (ix1 < ix2) /* Subpixel iy shifts to the right from column x */
2678 {
2679 if (x + ix1 >= 0 && x + ix2 < ncols)
2680 {
2681 xx = x + ix1;
2682 yy = y + ycen_offset[x] - ycen_offset[xx];
2683 cr2res_extract_zeta_add(zeta, m_zeta, z_rng, ncols, nrows,
2684 osample, x, iy, xx, yy,
2685 w - fabs(delta - ix1) * w);
2686 xx = x + ix2;
2687 yy = y + ycen_offset[x] - ycen_offset[xx];
2688 cr2res_extract_zeta_add(zeta, m_zeta, z_rng, ncols, nrows,
2689 osample, x, iy, xx, yy,
2690 fabs(delta - ix1) * w);
2691 }
2692 }
2693 else if (ix1 > ix2) /* Subpixel iy shifts to the left from column x */
2694 {
2695 if (x + ix2 >= 0 && x + ix1 < ncols)
2696 {
2697 xx = x + ix2;
2698 yy = y + ycen_offset[x] - ycen_offset[xx];
2699 cr2res_extract_zeta_add(zeta, m_zeta, z_rng, ncols, nrows,
2700 osample, x, iy, xx, yy,
2701 fabs(delta - ix1) * w);
2702 xx = x + ix1;
2703 yy = y + ycen_offset[x] - ycen_offset[xx];
2704 cr2res_extract_zeta_add(zeta, m_zeta, z_rng, ncols, nrows,
2705 osample, x, iy, xx, yy,
2706 w - fabs(delta - ix1) * w);
2707 }
2708 }
2709 else /* Subpixel iy stays inside column x */
2710 {
2711 xx = x + ix1;
2712 yy = y + ycen_offset[x] - ycen_offset[xx];
2713 cr2res_extract_zeta_add(zeta, m_zeta, z_rng, ncols, nrows,
2714 osample, x, iy, xx, yy, w);
2715 }
2716 }
2717 }
2718 }
2719 return 0;
2720}
2721
2722/*----------------------------------------------------------------------------*/
2760/*----------------------------------------------------------------------------*/
2761static int cr2res_extract_slit_func_curved(
2762 double error_factor,
2763 int ncols,
2764 int nrows,
2765 int osample,
2766 double pclip,
2767 double * im,
2768 double * pix_unc,
2769 int * mask,
2770 double * ycen,
2771 int * ycen_offset,
2772 int y_lower_lim,
2773 const double * slitcurve,
2774 int delta_x,
2775 double * sL,
2776 double * sP,
2777 double * model,
2778 double * unc,
2779 double lambda_sP,
2780 double lambda_sL,
2781 double sP_stop,
2782 int maxiter,
2783 double kappa,
2784 const double * slit_func_in,
2785 double * sP_old,
2786 double * l_Aij,
2787 double * p_Aij,
2788 double * l_bj,
2789 double * p_bj,
2790 double * zw,
2791 int * zk,
2792 zeta_ref * zeta,
2793 int * m_zeta,
2794 zeta_rng * z_rng)
2795{
2796 int x, xx, y, yy, iy, n, m, nk, mz, ny, nx;
2797 double norm, lambda, diag_tot, ww, dev, cost, tmp, sum;
2798 double sP_change, sP_med;
2799 cpl_vector * tmp_vec;
2800 int nclip, iter, isum;
2801
2802 /* The size of the sL array. */
2803 /* Extra osample is because ycen can be between 0 and 1. */
2804 ny = osample * (nrows + 1) + 1;
2805 nx = 4 * delta_x + 1;
2806 if (nx < 3)
2807 nx = 3;
2808
2809 cr2res_extract_zeta_tensors(ncols, nrows, ycen, ycen_offset,
2810 y_lower_lim, osample, slitcurve, zeta, m_zeta,
2811 z_rng);
2812
2813 // If a slit func is given, use that instead of recalculating it
2814 if (slit_func_in != NULL) {
2815 // Normalize the input just in case
2816 norm = 0.e0;
2817 for (iy = 0; iy < ny; iy++) {
2818 sL[iy] = slit_func_in[iy];
2819 norm += sL[iy];
2820 }
2821 norm /= osample;
2822 for (iy = 0; iy < ny; iy++)
2823 sL[iy] /= norm;
2824 }
2825
2826 /* Pre-clip the most extreme pixels from the mask */
2827 nclip = (int)(ncols * nrows * pclip / 100);
2828 if (nclip > 0)
2829 {
2830 if (nclip > ncols * nrows / 2)
2831 nclip = (ncols * nrows / 2) - 1;
2832 cpl_vector *im_vec = cpl_vector_wrap(ncols * nrows, im);
2833 cpl_vector *pos_vec = cpl_vector_new(ncols * nrows);
2834 for (x = 0; x < ncols; x++)
2835 {
2836 for (y = 0; y < nrows; y++)
2837 {
2838 cpl_vector_set(pos_vec, y * ncols + x, y * ncols + x);
2839 }
2840 }
2841 cpl_bivector *im_in_bvec = cpl_bivector_wrap_vectors(pos_vec, im_vec);
2842 cpl_bivector *im_sort_bvec = cpl_bivector_new(ncols * nrows);
2843 cpl_bivector_sort(im_sort_bvec, im_in_bvec, CPL_SORT_ASCENDING, CPL_SORT_BY_Y);
2844 cpl_bivector_unwrap_vectors(im_in_bvec);
2845 cpl_vector_unwrap(im_vec);
2846 cpl_vector_delete(pos_vec);
2847
2848 cpl_vector *pos_sort_vec = cpl_bivector_get_x(im_sort_bvec);
2849
2850 int ii = 0, mm = 0;
2851 while (mm < nclip)
2852 {
2853 int pos = (int)cpl_vector_get(pos_sort_vec, ii);
2854 if (mask[pos] != 0)
2855 {
2856 mask[pos] = 0;
2857 mm++;
2858 }
2859 ii++;
2860 }
2861
2862 ii = ncols * nrows - 1;
2863 mm = 0;
2864 while (mm < nclip)
2865 {
2866 int pos = (int)cpl_vector_get(pos_sort_vec, ii);
2867 if (mask[pos] != 0)
2868 {
2869 mask[pos] = 0;
2870 mm++;
2871 }
2872 ii--;
2873 }
2874 cpl_bivector_delete(im_sort_bvec);
2875 pos_sort_vec = NULL;
2876 }
2877
2878 /* Loop through sL , sP reconstruction until convergence is reached */
2879 iter = 0;
2880 do {
2881 /* Save the spectrum from the previous iteration */
2882 for (x = 0; x < ncols; x++)
2883 sP_old[x] = sP[x];
2884
2885 if (slit_func_in == NULL) {
2886 /* Compute slit function sL */
2887
2888 /* Prepare the RHS and the matrix */
2889 for (iy = 0; iy < ny; iy++)
2890 l_bj[iy] = 0.e0; /* Clean RHS */
2891 for (iy = 0; iy < ny * (4 * osample + 1); iy++)
2892 l_Aij[iy] = 0.e0;
2893
2894 /* Fill in SLE arrays for slit function */
2895 for (xx = 0; xx < ncols; xx++) {
2896 for (yy = 0; yy < nrows; yy++) {
2897 const zeta_ref *zrow;
2898 const zeta_rng *zr;
2899 double imv;
2900 int k0, rng;
2901 mz = m_zeta[mzeta_index(xx, yy)];
2902 if (mz <= 0 || !mask[yy * ncols + xx])
2903 continue;
2904 zrow = &zeta[zeta_index(xx, yy, 0)];
2905 imv = im[yy * ncols + xx];
2906 /* Merge entries sharing the same subpixel index iy: only
2907 the summed weight enters both the matrix and the RHS.
2908 The iy of one pixel span at most 2*osample+1 indices
2909 (the band width assumed by the matrix), so merge into
2910 a dense window zw[iy - k0] instead of searching a
2911 list of unique keys. */
2912 zr = &z_rng[mzeta_index(xx, yy)];
2913 k0 = zr->min_iy;
2914 rng = zr->max_iy - k0;
2915 if (rng <= 2 * osample) {
2916 for (n = 0; n <= rng; n++)
2917 zw[n] = 0.e0;
2918 for (m = 0; m < mz; m++)
2919 zw[zrow[m].iy - k0] += sP[zrow[m].x] * zrow[m].w;
2920 /* The matrix is symmetric: accumulate each unordered
2921 pair once into the upper bands (mirrored below
2922 after the fill). Walking the window row-wise makes
2923 the inner loop contiguous in both operands. Window
2924 entries between actual keys are zero and add
2925 exactly nothing. */
2926 for (m = 0; m <= rng; m++) {
2927 const double um = zw[m];
2928 const double *restrict uv = zw + m;
2929 double *restrict arow =
2930 &l_Aij[laij_index(k0 + m, 2 * osample)];
2931 const int dmax = rng - m;
2932 /* element-wise independent fmas: vectorization
2933 does not reorder per-element arithmetic */
2934#ifdef _OPENMP
2935#pragma omp simd
2936#endif
2937 for (n = 0; n <= dmax; n++)
2938 arow[n] += um * uv[n];
2939 l_bj[k0 + m] += imv * um;
2940 }
2941 continue;
2942 }
2943 /* Over-wide list (extreme geometry): merge by searching
2944 unique keys, as before */
2945 nk = 0;
2946 for (m = 0; m < mz; m++) {
2947 const int key = zrow[m].iy;
2948 const double v = sP[zrow[m].x] * zrow[m].w;
2949 for (n = 0; n < nk; n++) {
2950 if (zk[n] == key) {
2951 zw[n] += v;
2952 break;
2953 }
2954 }
2955 if (n == nk) {
2956 zk[nk] = key;
2957 zw[nk++] = v;
2958 }
2959 }
2960 for (m = 0; m < nk; m++) {
2961 const double um = zw[m];
2962 iy = zk[m];
2963 l_Aij[laij_index(iy, 2 * osample)] += um * um;
2964 for (n = m + 1; n < nk; n++) {
2965 const int iyn = zk[n];
2966 const int lo = min(iy, iyn);
2967 const int d = iyn > iy ? iyn - iy : iy - iyn;
2968 l_Aij[laij_index(lo, d + 2 * osample)] +=
2969 zw[n] * um;
2970 }
2971 l_bj[iy] += imv * um;
2972 }
2973 }
2974 }
2975
2976 /* Mirror the upper bands into the lower bands:
2977 A[r+d, 2o-d] = A[r, 2o+d] */
2978 for (m = 1; m <= 2 * osample; m++)
2979 for (iy = 0; iy < ny - m; iy++)
2980 l_Aij[laij_index(iy + m, 2 * osample - m)] =
2981 l_Aij[laij_index(iy, 2 * osample + m)];
2982
2983 diag_tot = 0.e0;
2984 for (iy = 0; iy < ny; iy++)
2985 diag_tot += l_Aij[laij_index(iy, 2 * osample)];
2986
2987 /* Scale regularization parameters */
2988 lambda = lambda_sL * diag_tot / ny;
2989
2990 /* Add regularization parts for the SLE matrix */
2991 l_Aij[laij_index(0, 2 * osample)] += lambda; /* Main diag */
2992 l_Aij[laij_index(0, 2 * osample + 1)] -= lambda; /* Upper diag */
2993 for (iy = 1; iy < ny - 1; iy++) {
2994 l_Aij[laij_index(iy, 2 * osample - 1)] -= lambda;
2995 l_Aij[laij_index(iy, 2 * osample)] += lambda * 2.e0;
2996 l_Aij[laij_index(iy, 2 * osample + 1)] -= lambda;
2997 }
2998 l_Aij[laij_index(ny - 1, 2 * osample - 1)] -= lambda;
2999 l_Aij[laij_index(ny - 1, 2 * osample)] += lambda;
3000
3001 /* Regularize diagonal to prevent singular matrix from fully
3002 masked rows */
3003 {
3004 double max_diag = 0.0;
3005 for (iy = 0; iy < ny; iy++)
3006 if (l_Aij[laij_index(iy, 2 * osample)] > max_diag)
3007 max_diag = l_Aij[laij_index(iy, 2 * osample)];
3008 if (max_diag > 0.0) {
3009 const double min_diag = max_diag * 1.0e-10;
3010 for (iy = 0; iy < ny; iy++)
3011 if (l_Aij[laij_index(iy, 2 * osample)] < min_diag)
3012 l_Aij[laij_index(iy, 2 * osample)] = min_diag;
3013 }
3014 }
3015
3016 /* Solve the system of equations */
3017 cr2res_extract_bandsol_rowmajor(l_Aij, l_bj, ny, 4 * osample + 1);
3018
3019 /* Normalize the slit function; sum of |sL| (not plain sum) keeps
3020 the flux scale of the historic cr2res convention when sL has
3021 negative parts (e.g. background residuals in nodding A-B) */
3022 norm = 0.e0;
3023 for (iy = 0; iy < ny; iy++) {
3024 sL[iy] = l_bj[iy];
3025 norm += fabs(sL[iy]);
3026 }
3027 norm /= osample;
3028 for (iy = 0; iy < ny; iy++)
3029 sL[iy] /= norm;
3030 }
3031
3032 /* Compute spectrum sP */
3033 for (x = 0; x < ncols; x++)
3034 p_bj[x] = 0.e0;
3035 for (x = 0; x < ncols * nx; x++)
3036 p_Aij[x] = 0.e0;
3037
3038 /* Pixel-centric fill, see the slit function SLE above */
3039 for (xx = 0; xx < ncols; xx++) {
3040 for (yy = 0; yy < nrows; yy++) {
3041 const zeta_ref *zrow;
3042 const zeta_rng *zr;
3043 double imv;
3044 int k0, rng;
3045 mz = m_zeta[mzeta_index(xx, yy)];
3046 if (mz <= 0 || !mask[yy * ncols + xx])
3047 continue;
3048 zrow = &zeta[zeta_index(xx, yy, 0)];
3049 imv = im[yy * ncols + xx];
3050 /* Merge entries sharing the same source column x; with small
3051 curvature this collapses the list to just a few entries.
3052 Sources span at most 2*delta_x+1 columns (the band width
3053 of the matrix), so merge into a dense window zw[x - k0]. */
3054 zr = &z_rng[mzeta_index(xx, yy)];
3055 k0 = zr->min_x;
3056 rng = zr->max_x - k0;
3057 if (rng <= 2 * delta_x) {
3058 for (n = 0; n <= rng; n++)
3059 zw[n] = 0.e0;
3060 for (m = 0; m < mz; m++)
3061 zw[zrow[m].x - k0] += sL[zrow[m].iy] * zrow[m].w;
3062 /* Symmetric matrix: upper bands only, mirrored after
3063 the fill */
3064 for (m = 0; m <= rng; m++) {
3065 const double um = zw[m];
3066 const double *restrict uv = zw + m;
3067 double *restrict arow =
3068 &p_Aij[paij_index(k0 + m, 2 * delta_x)];
3069 const int dmax = rng - m;
3070#ifdef _OPENMP
3071#pragma omp simd
3072#endif
3073 for (n = 0; n <= dmax; n++)
3074 arow[n] += um * uv[n];
3075 p_bj[k0 + m] += imv * um;
3076 }
3077 continue;
3078 }
3079 /* Over-wide list: merge by searching unique keys */
3080 nk = 0;
3081 for (m = 0; m < mz; m++) {
3082 const int key = zrow[m].x;
3083 const double v = sL[zrow[m].iy] * zrow[m].w;
3084 for (n = 0; n < nk; n++) {
3085 if (zk[n] == key) {
3086 zw[n] += v;
3087 break;
3088 }
3089 }
3090 if (n == nk) {
3091 zk[nk] = key;
3092 zw[nk++] = v;
3093 }
3094 }
3095 for (m = 0; m < nk; m++) {
3096 const double um = zw[m];
3097 x = zk[m];
3098 p_Aij[paij_index(x, 2 * delta_x)] += um * um;
3099 for (n = m + 1; n < nk; n++) {
3100 const int xn = zk[n];
3101 const int lo = min(x, xn);
3102 const int d = xn > x ? xn - x : x - xn;
3103 p_Aij[paij_index(lo, d + 2 * delta_x)] += zw[n] * um;
3104 }
3105 p_bj[x] += imv * um;
3106 }
3107 }
3108 }
3109
3110 /* Mirror the upper bands into the lower bands */
3111 for (m = 1; m <= 2 * delta_x; m++)
3112 for (x = 0; x < ncols - m; x++)
3113 p_Aij[paij_index(x + m, 2 * delta_x - m)] =
3114 p_Aij[paij_index(x, 2 * delta_x + m)];
3115
3116 if (lambda_sP > 0.e0) {
3117 lambda = lambda_sP; /* Scale regularization parameter */
3118 p_Aij[paij_index(0, 2 * delta_x)] += lambda; /* Main diag */
3119 p_Aij[paij_index(0, 2 * delta_x + 1)] -= lambda; /* Upper diag */
3120 for (x = 1; x < ncols - 1; x++) {
3121 p_Aij[paij_index(x, 2 * delta_x - 1)] -= lambda;
3122 p_Aij[paij_index(x, 2 * delta_x)] += lambda * 2.e0;
3123 p_Aij[paij_index(x, 2 * delta_x + 1)] -= lambda;
3124 }
3125 p_Aij[paij_index(ncols - 1, 2 * delta_x - 1)] -= lambda;
3126 p_Aij[paij_index(ncols - 1, 2 * delta_x)] += lambda;
3127 }
3128
3129 /* Regularize diagonal to prevent singular matrix from fully masked
3130 columns. When a column has no valid data (all pixels masked), the
3131 corresponding row of the matrix is zero, causing division by zero
3132 in bandsol. The resulting spectrum value for masked columns will
3133 be ~0 (from p_bj[x]/diag). */
3134 {
3135 double max_diag = 0.0;
3136 for (x = 0; x < ncols; x++)
3137 if (p_Aij[paij_index(x, 2 * delta_x)] > max_diag)
3138 max_diag = p_Aij[paij_index(x, 2 * delta_x)];
3139 if (max_diag > 0.0) {
3140 const double min_diag = max_diag * 1.0e-10;
3141 for (x = 0; x < ncols; x++)
3142 if (p_Aij[paij_index(x, 2 * delta_x)] < min_diag)
3143 p_Aij[paij_index(x, 2 * delta_x)] = min_diag;
3144 }
3145 }
3146
3147 /* Solve the system of equations */
3148 cr2res_extract_bandsol_rowmajor(p_Aij, p_bj, ncols, nx);
3149
3150 for (x = 0; x < ncols; x++)
3151 sP[x] = p_bj[x]; /* New Spectrum vector */
3152
3153 /* Compute median value of the spectrum for normalisation purpose */
3154 tmp_vec = cpl_vector_wrap(ncols, sP);
3155 sP_med = fabs(cpl_vector_get_median_const(tmp_vec));
3156 cpl_vector_unwrap(tmp_vec);
3157
3158 /* Compute the change in the spectrum */
3159 sP_change = 0.e0;
3160 for (x = 0; x < ncols; x++) {
3161 if (fabs(sP[x] - sP_old[x]) > sP_change)
3162 sP_change = fabs(sP[x] - sP_old[x]);
3163 }
3164
3165 /* Compute the model. x is the outer loop so that the zeta tensor,
3166 by far the largest array, is read sequentially */
3167 for (x = 0; x < ncols; x++) {
3168 for (y = 0; y < nrows; y++) {
3169 const zeta_ref *zrow = &zeta[zeta_index(x, y, 0)];
3170 double acc = 0.;
3171 mz = m_zeta[mzeta_index(x, y)];
3172 for (m = 0; m < mz; m++) {
3173 xx = zrow[m].x;
3174 iy = zrow[m].iy;
3175 ww = zrow[m].w;
3176 acc += sP[xx] * sL[iy] * ww;
3177 }
3178 model[y * ncols + x] = acc;
3179 }
3180 }
3181
3182 /* Compare model and data: reduced chi-square as convergence
3183 criterion and RMS of residuals for outlier rejection */
3184 cost = 0.e0;
3185 sum = 0.e0;
3186 isum = 0;
3187 for (y = 0; y < nrows; y++) {
3188 for (x = delta_x; x < ncols - delta_x; x++) {
3189 if (mask[y * ncols + x]) {
3190 tmp = model[y * ncols + x] - im[y * ncols + x];
3191 sum += tmp * tmp;
3192 tmp /= max(pix_unc[y * ncols + x], 1);
3193 cost += tmp * tmp;
3194 isum++;
3195 }
3196 }
3197 }
3198 cost /= (isum - (ncols + ny));
3199 dev = sqrt(sum / isum);
3200
3201 /* Adjust the mask marking outliers */
3202 if (kappa > 0) {
3203 for (y = 0; y < nrows; y++) {
3204 for (x = delta_x; x < ncols - delta_x; x++) {
3205 if (fabs(model[y * ncols + x] - im[y * ncols + x]) <
3206 kappa * dev)
3207 mask[y * ncols + x] = 1;
3208 else
3209 mask[y * ncols + x] = 0;
3210 }
3211 }
3212 }
3213
3214 cpl_msg_debug(
3215 __func__,
3216 "Iter: %i, Sigma: %g, Cost: %g, sP_change: %g, sP_lim: %g",
3217 iter, dev, cost, sP_change, sP_stop * sP_med);
3218
3219 iter++;
3220
3221 /* Check for convergence: largest change in the spectrum between
3222 two iterations, relative to its median (historic criterion,
3223 empirically robust on a large variety of data). On NaNs the
3224 comparison is false and the loop exits, as in the old code. */
3225 } while (iter == 1 ||
3226 (iter <= maxiter && sP_change > sP_stop * sP_med));
3227
3228 if (iter > maxiter && sP_change > sP_stop * sP_med)
3229 cpl_msg_warning(
3230 __func__,
3231 "Maximum number of %d iterations reached without converging.",
3232 maxiter);
3233
3234 /* Flip sign if converged in negative direction */
3235 sum = 0.0;
3236 for (y = 0; y < ny; y++)
3237 sum += sL[y];
3238 if (sum < 0.0) {
3239 for (y = 0; y < ny; y++)
3240 sL[y] *= -1.0;
3241 for (x = 0; x < ncols; x++)
3242 sP[x] *= -1.0;
3243 }
3244
3245 if (error_factor == -1)
3246 {
3247 // Uncertainty calculation, following Horne 1986.
3248 for (x = 0; x < ncols; x++)
3249 {
3250 double num_sum;
3251 double den_sum;
3252 double model_sum = 0.0;
3253
3254 unc[x] = 0.0;
3255 num_sum = 0.0;
3256 den_sum = 0.0;
3257 for (y = 0; y < nrows; y++)
3258 {
3259 model_sum += model[y * ncols + x];
3260 }
3261 for (y = 0; y < nrows; y++)
3262 {
3263 double model_norm = model[y * ncols + x] / model_sum;
3264 num_sum += model_norm * mask[y * ncols + x];
3265 den_sum += (model_norm * model_norm) * mask[y * ncols + x] / (pix_unc[y * ncols + x] * pix_unc[y * ncols + x]);
3266 }
3267 if (den_sum != 0)
3268 {
3269 unc[x] = sqrt(fabs(num_sum / den_sum));
3270 }
3271 else
3272 {
3273 unc[x] = NAN;
3274 }
3275 }
3276 }
3277 else
3278 {
3279 // Uncertainty calculation only using total object flux, needs later correction.
3280 for (x = 0; x < ncols; x++)
3281 {
3282 double msum;
3283
3284 unc[x] = 0.0;
3285 msum = 0.0;
3286 sum = 0.0;
3287 for (y = 0; y < nrows; y++)
3288 {
3289 if (mask[y * ncols + x])
3290 {
3291 msum += (im[y * ncols + x] * model[y * ncols + x]) *
3292 mask[y * ncols + x];
3293 sum += (model[y * ncols + x] * model[y * ncols + x]) *
3294 mask[y * ncols + x];
3295 }
3296 }
3297 if (msum != 0)
3298 {
3299 // This can give NaNs if m/sum is less than zero, i.e. low/no flux
3300 // due to ignoring background flux.
3301 unc[x] = sqrt(fabs(sP[x]) * fabs(sum) / fabs(msum) / error_factor);
3302 }
3303 else
3304 {
3305 // Fix bad value to NaN as Phase3 doesn't allow Inf.
3306 unc[x] = NAN;
3307 }
3308 }
3309 }
3310 return 0;
3311}
3312
3313/*----------------------------------------------------------------------------*/
3328/*----------------------------------------------------------------------------*/
3329static int cr2res_extract_bandsol_rowmajor(
3330 double * a,
3331 double * r,
3332 int n,
3333 int nd)
3334{
3335 double aa;
3336 int i, j, k;
3337
3338 /* Forward sweep */
3339 for (i = 0; i < n - 1; i++)
3340 {
3341 aa = a[i * nd + nd / 2];
3342 r[i] /= aa;
3343 for (j = 0; j < nd; j++)
3344 a[i * nd + j] /= aa;
3345 for (j = 1; j < min(nd / 2 + 1, n - i); j++)
3346 {
3347 aa = a[(i + j) * nd + nd / 2 - j];
3348 r[i + j] -= r[i] * aa;
3349 for (k = 0; k < nd - j; k++)
3350 a[(i + j) * nd + k] -= a[i * nd + k + j] * aa;
3351 }
3352 }
3353
3354 /* Backward sweep */
3355 r[n - 1] /= a[(n - 1) * nd + nd / 2];
3356 for (i = n - 1; i > 0; i--)
3357 {
3358 for (j = 1; j <= min(nd / 2, i); j++)
3359 r[i - j] -= r[i] * a[(i - j) * nd + nd / 2 + j];
3360 r[i - 1] /= a[(i - 1) * nd + nd / 2];
3361 }
3362
3363 r[0] /= a[nd / 2];
3364 return 0;
3365}
3366
3367/*----------------------------------------------------------------------------*/
3393/*----------------------------------------------------------------------------*/
3394int cr2res_extract_slitdec_bandsol(
3395 double * a,
3396 double * r,
3397 int n,
3398 int nd,
3399 double lambda)
3400{
3401 double aa;
3402 int i, j, k;
3403
3404 //if(fmod(nd,2)==0) return -1;
3405
3406 /* Forward sweep */
3407 for(i=0; i<n-1; i++)
3408 {
3409 aa=a[i+n*(nd/2)];
3410 if(aa==0.e0) aa = lambda; //return -3;
3411 r[i]/=aa;
3412 for(j=0; j<nd; j++) a[i+j*n]/=aa;
3413 for(j=1; j<min(nd/2+1,n-i); j++)
3414 {
3415 aa=a[i+j+n*(nd/2-j)];
3416 r[i+j]-=r[i]*aa;
3417 for(k=0; k<n*(nd-j); k+=n) a[i+j+k]-=a[i+k+n*j]*aa;
3418 }
3419 }
3420
3421 /* Backward sweep */
3422 aa = a[n-1+n*(nd/2)];
3423 if (aa == 0) aa = lambda; //return -4;
3424 r[n-1]/=aa;
3425 for(i=n-1; i>0; i--)
3426 {
3427 for(j=1; j<=min(nd/2,i); j++){
3428 r[i-j]-=r[i]*a[i-j+n*(nd/2+j)];
3429 }
3430 aa = a[i-1+n*(nd/2)];
3431 if(aa==0.e0) aa = lambda; //return -5;
3432
3433 r[i-1]/=aa;
3434 }
3435
3436 aa = a[n*(nd/2)];
3437 if(aa==0.e0) aa = lambda; //return -6;
3438 r[0]/=aa;
3439 return 0;
3440}
3441
3442/*----------------------------------------------------------------------------*/
3453/*----------------------------------------------------------------------------*/
3454static int cr2res_extract_slitdec_adjust_swath(
3455 cpl_vector * ycen,
3456 int height,
3457 int leny,
3458 int sw,
3459 int lenx,
3460 int dx,
3461 cpl_vector ** bins_begin,
3462 cpl_vector ** bins_end)
3463{
3464 if (sw <= 0 || lenx <= 0) return -1;
3465 if (bins_begin == NULL || bins_end == NULL) return -1;
3466 if (ycen == NULL) return -1;
3467
3468 // special case only one swath
3469 if (sw==CR2RES_DETECTOR_SIZE){
3470 *bins_begin = cpl_vector_new(1);
3471 *bins_end = cpl_vector_new(1);
3472 cpl_vector_set(*bins_begin, 0, 0);
3473 cpl_vector_set(*bins_end, 0, CR2RES_DETECTOR_SIZE);
3474 return CR2RES_DETECTOR_SIZE;
3475 }
3476
3477 int nbin, nx, i = 0;
3478 double step = 0;
3479 int start = 0, end = lenx;
3480 // Setting them as int, makes this comparable to using ycen_int
3481 int ymin, ymax;
3482 int * ycen_int = cr2res_vector_get_int(ycen);
3483
3484 for (i=0; i < lenx; i++){
3485 ymin = ycen_int[i] - (height/2);
3486 ymax = ycen_int[i] + (height/2) + height%2 ;
3487 if (!((ymax <= 1) || (ymin > leny))){
3488 start = i;
3489 break;
3490 }
3491 start = lenx;
3492 }
3493 for (i = lenx - 1; i >= 0; i--){
3494 ymin = ycen_int[i] - (height/2);
3495 ymax = ycen_int[i] + (height/2) + height%2 ;
3496 if (!((ymax <= 1) || (ymin > leny))){
3497 end = i + 1;
3498 break;
3499 }
3500 end = 0;
3501 }
3502 cpl_free(ycen_int);
3503 nx = end - start;
3504 if (nx <= 0){
3505 // No valid points in this order
3506 return -1;
3507 }
3508
3509 if (sw > nx - 2 * dx) {
3510 sw = nx - 2 * dx;
3511 if (sw % 2 == 1) sw -= 1;
3512 } else if (sw % 2 == 1) sw += 1;
3513
3514 // Calculate number of bins
3515 nbin = 2 * ((nx - 2 * dx) / sw);
3516 if ((nx - 2 * dx) % sw > sw / 2) nbin++;
3517 if (nbin < 1) nbin = 1;
3518
3519 // Step / 2, to get half width swaths
3520 step = sw / 2;
3521 *bins_begin = cpl_vector_new(nbin);
3522 *bins_end = cpl_vector_new(nbin);
3523
3524 // boundaries of bins
3525 for(i = 0; i < nbin; i++)
3526 {
3527 int bin;
3528 bin = start + min(i * step, nx - sw - 2 * dx);
3529 cpl_vector_set(*bins_begin, i, bin);
3530 cpl_vector_set(*bins_end, i, bin + sw + 2 * dx);
3531 cpl_msg_debug(__func__, "Swath %d goes from %d to %d.", i, bin,
3532 bin + sw + 2 * dx);
3533 }
3534 return sw + 2 * dx;
3535}
cr2res_dfs_SLIT_FUNC_colname
char * cr2res_dfs_SLIT_FUNC_colname(int order_idx, int trace)
Get the SLIT_FUNC table column name for a given order/trace.
Definition cr2res_dfs.c:432
cr2res_dfs_SPEC_ERR_colname
char * cr2res_dfs_SPEC_ERR_colname(int order_idx, int trace)
Get the ERR column name for a given order/trace.
Definition cr2res_dfs.c:414
cr2res_dfs_WAVELENGTH_colname
char * cr2res_dfs_WAVELENGTH_colname(int order_idx, int trace)
Get the WAVELENGTH column name for a given order/trace.
Definition cr2res_dfs.c:396
cr2res_dfs_POSITIONY_colname
char * cr2res_dfs_POSITIONY_colname(int order_idx, int trace)
Get the POSITIONY table column name for a given order/trace.
Definition cr2res_dfs.c:468
cr2res_dfs_SLIT_FRACTION_colname
char * cr2res_dfs_SLIT_FRACTION_colname(int order_idx, int trace)
Get the SLIT_FRACTION table column name for a given order/trace.
Definition cr2res_dfs.c:486
cr2res_dfs_SPEC_colname
char * cr2res_dfs_SPEC_colname(int order_idx, int trace)
Get the SPEC column name for a given order/trace.
Definition cr2res_dfs.c:378
cr2res_dfs_POSITIONX_colname
char * cr2res_dfs_POSITIONX_colname(int order_idx, int trace)
Get the POSITIONX table column name for a given order/trace.
Definition cr2res_dfs.c:450
cr2res_extract_EXTRACT1D_create
cpl_table * cr2res_extract_EXTRACT1D_create(cpl_bivector **spectrum, const cpl_table *trace_table)
Create the extract 1D table to be saved.
Definition cr2res_extract.c:1097
cr2res_extract_SLIT_FUNC_get_vector
int cr2res_extract_SLIT_FUNC_get_vector(const cpl_table *tab, int order, int trace_nb, cpl_vector **vec)
Get a vector from the SLIT_FUNC table.
Definition cr2res_extract.c:1366
cr2res_extract_median
int cr2res_extract_median(const hdrl_image *hdrl_in, const cpl_table *trace_tab, int order, int trace_id, int height, cpl_vector **slit_func, cpl_bivector **spec, hdrl_image **model)
Simple extraction function with the median.
Definition cr2res_extract.c:703
cr2res_extract_sum_tilt
int cr2res_extract_sum_tilt(const hdrl_image *hdrl_in, const cpl_table *trace_tab, int order, int trace_id, int height, cpl_vector **slit_func, cpl_bivector **spec, hdrl_image **model)
Simple extraction function with curvature correction.
Definition cr2res_extract.c:846
cr2res_extract_EXTRACT2D_create
cpl_table * cr2res_extract_EXTRACT2D_create(cpl_bivector **spectrum, cpl_bivector **position, cpl_vector **wavelength, cpl_vector **slit_fraction, const cpl_table *trace_table)
Create the extract 2D table to be saved.
Definition cr2res_extract.c:2373
cr2res_extract_EXTRACT1D_get_spectrum
int cr2res_extract_EXTRACT1D_get_spectrum(const cpl_table *tab, int order, int trace_nb, cpl_bivector **spec, cpl_bivector **spec_err)
Get a Spectrum and its error from the EXTRACT_1D table.
Definition cr2res_extract.c:1289
cr2res_extract_sum_vert
int cr2res_extract_sum_vert(const hdrl_image *hdrl_in, const cpl_table *trace_tab, int order, int trace_id, int height, cpl_vector **slit_func, cpl_bivector **spec, hdrl_image **model)
Simple extraction function.
Definition cr2res_extract.c:555
cr2res_extract2d_trace
int cr2res_extract2d_trace(const hdrl_image *in, const cpl_table *trace_tab, int order, int trace_id, int npoints, const cpl_image *wavemap, const cpl_image *slitmap, cpl_bivector **spectrum, cpl_bivector **position, cpl_vector **wavelength, cpl_vector **slit_fraction)
Extraction2d function.
Definition cr2res_extract.c:2243
cr2res_extract_SLITFUNC_create
cpl_table * cr2res_extract_SLITFUNC_create(cpl_vector **slit_func, const cpl_table *trace_table)
Create the slit functions table to be saved.
Definition cr2res_extract.c:1207
cr2res_extract_slitdec_curved
int cr2res_extract_slitdec_curved(const hdrl_image *img_hdrl, const cpl_table *trace_tab, const cpl_vector *slit_func_vec_in, int order, int trace_id, int height, int swath, int oversample, double pclip, double smooth_slit, double smooth_spec, int niter, double kappa, double error_factor, cpl_vector **slit_func, cpl_bivector **spec, hdrl_image **model)
Extract optimally (slit-decomposition) with polynomial slit.
Definition cr2res_extract.c:1438
cr2res_extract_traces
int cr2res_extract_traces(const hdrl_image *img, const cpl_table *traces, const cpl_table *slit_func_in, const cpl_table *blaze_table_in, float blaze_norm, int reduce_order, int reduce_trace, cr2res_extr_method extr_method, int extr_height, int swath_width, int oversample, double pclip, double smooth_slit, double smooth_spec, int niter, double kappa, double error_factor, int display, int disp_order_idx, int disp_trace, cpl_table **extracted, cpl_table **slit_func, hdrl_image **model_master)
Extracts all the passed traces at once.
Definition cr2res_extract.c:186
cr2res_extract2d_traces
int cr2res_extract2d_traces(const hdrl_image *img, const cpl_table *traces, int reduce_order, int reduce_trace, cpl_table **extracted)
Extract2d all the passed traces at once.
Definition cr2res_extract.c:2119
cr2res_trace_get_ycen
cpl_vector * cr2res_trace_get_ycen(const cpl_table *trace, int order_idx, int trace_nb, int size)
Retrieves the middle (All) polynomial from trace table and evaluates.
Definition cr2res_trace.c:769
cr2res_trace_get_height
int cr2res_trace_get_height(const cpl_table *trace, cpl_size order_idx, cpl_size trace_nb)
Computes the average height (pix) of an order, from trace polys.
Definition cr2res_trace.c:805
cr2res_get_trace_table_index
cpl_size cr2res_get_trace_table_index(const cpl_table *trace_wave, int order_idx, int trace_nb)
Get the index in a TRACE_WAVE table.
Definition cr2res_trace.c:652
cr2res_get_trace_wave_poly
cpl_polynomial * cr2res_get_trace_wave_poly(const cpl_table *trace_wave, const char *poly_column, int order_idx, int trace_nb)
Get a polynomial from a TRACE_WAVE table.
Definition cr2res_trace.c:685
cr2res_trace_get_wl
cpl_vector * cr2res_trace_get_wl(const cpl_table *trace_wave, int order_idx, int trace_nb, int size)
Get the Wavelength vector from a TRACE_WAVE table.
Definition cr2res_trace.c:730
cr2res_image_cut_rectify
cpl_image * cr2res_image_cut_rectify(const cpl_image *img_in, const cpl_vector *ycen, int height)
Cut a bent order into a rectangle, shifting columns.
Definition cr2res_utils.c:907
cr2res_vector_get_int
int * cr2res_vector_get_int(const cpl_vector *ycen)
Definition cr2res_utils.c:230
cr2res_convert_array_to_poly
cpl_polynomial * cr2res_convert_array_to_poly(const cpl_array *arr)
Convert an array to polynomial.
Definition cr2res_utils.c:1305
cr2res_vector_get_rest
double * cr2res_vector_get_rest(const cpl_vector *ycen)
Definition cr2res_utils.c:253
cr2res_polynomial_eval_vector
cpl_vector * cr2res_polynomial_eval_vector(const cpl_polynomial *poly, const cpl_vector *vec)
Evaluate a polynomial on a vector.
Definition cr2res_utils.c:1059
cr2res_slit_pos_image
int cr2res_slit_pos_image(const cpl_table *trace_wave, cpl_image **slitpos, cpl_image **wavelength)
get a image of the slitposition (and wavelength) along the slit
Definition cr2res_utils.c:841
cr2res_image_insert_rect
int cr2res_image_insert_rect(const cpl_image *rect_in, const cpl_vector *ycen, cpl_image *img_out)
Re-insert a rectangular cut-out of an order into the full frame.
Definition cr2res_utils.c:984
hdrl_image_mul_scalar
cpl_error_code hdrl_image_mul_scalar(hdrl_image *self, hdrl_value value)
Elementwise multiplication of an image with a scalar.
Definition hdrl_image_math.c:245
hdrl_image_duplicate
hdrl_image * hdrl_image_duplicate(const hdrl_image *himg)
copy hdrl_image
Definition hdrl_image.c:391
hdrl_image_get_error
cpl_image * hdrl_image_get_error(hdrl_image *himg)
get error as cpl image
Definition hdrl_image.c:131
hdrl_image_get_size_y
cpl_size hdrl_image_get_size_y(const hdrl_image *self)
return size of Y dimension of image
Definition hdrl_image.c:540
hdrl_image_get_size_x
cpl_size hdrl_image_get_size_x(const hdrl_image *self)
return size of X dimension of image
Definition hdrl_image.c:525
hdrl_image_get_error_const
const cpl_image * hdrl_image_get_error_const(const hdrl_image *himg)
get error as cpl image
Definition hdrl_image.c:144
hdrl_image_create
hdrl_image * hdrl_image_create(const cpl_image *image, const cpl_image *error)
create a new hdrl_image from to existing images by copying them
Definition hdrl_image.c:295
hdrl_image_get_image
cpl_image * hdrl_image_get_image(hdrl_image *himg)
get data as cpl image
Definition hdrl_image.c:105
hdrl_image_get_image_const
const cpl_image * hdrl_image_get_image_const(const hdrl_image *himg)
get data as cpl image
Definition hdrl_image.c:118
hdrl_image_new
hdrl_image * hdrl_image_new(cpl_size nx, cpl_size ny)
create new zero filled hdrl image
Definition hdrl_image.c:311
hdrl_image_delete
void hdrl_image_delete(hdrl_image *himg)
delete hdrl_image
Definition hdrl_image.c:379
Generated on Fri Jun 19 2026 14:54:49 for CR2RE Pipeline Reference Manual by doxygen 1.13.2