CR2RE Pipeline Reference Manual 1.6.11
cr2res_detlin.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#ifdef HAVE_CONFIG_H
21#include <config.h>
22#endif
23
24/*-----------------------------------------------------------------------------
25 Includes
26 -----------------------------------------------------------------------------*/
27
28#include <math.h>
29#include <cpl.h>
30#include "cr2res_detlin.h"
31#include "cr2res_pfits.h"
32#include "cr2res_qc.h"
33
34/*-----------------------------------------------------------------------------
35 Defines
36 -----------------------------------------------------------------------------*/
37
38#define pow2(x) (x)*(x)
39#define pow3(x) (x)*(x)*(x)
40
41/*-----------------------------------------------------------------------------
42 Functions prototypes
43 -----------------------------------------------------------------------------*/
44
45static int cr2res_detlin_frames_dit_compare(
46 const cpl_frame * in1,
47 const cpl_frame * in2) ;
48static void cr2res_matrix_fill_normal_vandermonde(cpl_matrix * self,
49 cpl_matrix * mx,
50 const cpl_vector * xhat,
51 cpl_boolean is_eqdist,
52 cpl_size mindeg,
53 const cpl_vector * values);
54
55/*----------------------------------------------------------------------------*/
59/*----------------------------------------------------------------------------*/
60
62
63/*----------------------------------------------------------------------------*/
70/*----------------------------------------------------------------------------*/
71int cr2res_detlin_correct(
72 hdrl_image * in,
73 const hdrl_imagelist * detlin)
74{
75 const cpl_image * ima ;
76 const cpl_image * erra ;
77 const cpl_image * imb ;
78 const cpl_image * errb ;
79 const cpl_image * imc ;
80 const cpl_image * errc ;
81 const double * pima ;
82 const double * perra ;
83 const double * pimb ;
84 const double * perrb ;
85 const double * pimc ;
86 const double * perrc ;
87 cpl_image * cur_ima ;
88 double * pdata ;
89 double * perr ;
90 int nx, ny ;
91 double correction_factor;
92 int i ;
93
94 /* Test entries */
95 if (!in || !detlin) return -1 ;
96
97 /* Initialise */
98 pdata = cpl_image_get_data_double(hdrl_image_get_image(in)) ;
99 perr = cpl_image_get_data_double(hdrl_image_get_error(in)) ;
100
101
102 /* Load the 3 coeffs images */
103 ima = hdrl_image_get_image_const(hdrl_imagelist_get_const(detlin, 0)) ;
104 erra = hdrl_image_get_error_const(hdrl_imagelist_get_const(detlin, 0)) ;
105 imb = hdrl_image_get_image_const(hdrl_imagelist_get_const(detlin, 1)) ;
106 errb = hdrl_image_get_error_const(hdrl_imagelist_get_const(detlin, 1)) ;
107 imc = hdrl_image_get_image_const(hdrl_imagelist_get_const(detlin, 2)) ;
108 errc = hdrl_image_get_error_const(hdrl_imagelist_get_const(detlin, 2)) ;
109
110 if (!ima || !imb || !imc ) {
111 cpl_msg_error(cpl_func, "Cannot access the detlin images") ;
112 return -1 ;
113 }
114 pima = cpl_image_get_data_double_const(ima) ;
115 pimb = cpl_image_get_data_double_const(imb) ;
116 pimc = cpl_image_get_data_double_const(imc) ;
117 perra = cpl_image_get_data_double_const(erra);
118 perrb = cpl_image_get_data_double_const(errb);
119 perrc = cpl_image_get_data_double_const(errc);
120
121 /* Test sizes */
122 cur_ima = hdrl_image_get_image(in) ;
123 nx = cpl_image_get_size_x(cur_ima) ;
124 ny = cpl_image_get_size_y(cur_ima) ;
125 if ((cpl_image_get_size_x(ima) != nx) ||
126 (cpl_image_get_size_x(imb) != nx) ||
127 (cpl_image_get_size_x(imc) != nx) ||
128 (cpl_image_get_size_y(ima) != ny) ||
129 (cpl_image_get_size_y(imb) != ny) ||
130 (cpl_image_get_size_y(imc) != ny)) {
131 cpl_msg_error(cpl_func, "Incompatible sizes") ;
132 return -1 ;
133 }
134
135 /* Loop on pixels */
136 for (i=0 ; i<nx*ny ; i++) {
137
138 // Only correct non-linear regime
139 if (pdata[i] < CR2RES_DETLIN_THRESHOLD) continue ;
140
141 // for each pixel p' = (a + b * p + c * p * p) * p
142 perr[i] = pow2(perra[i] * pdata[i]) + pow2(perrb[i] * pow2(pdata[i]))
143 + pow2(perrc[i] * pow3(pdata[i]))
144 + pow2(perr[i] * (pima[i] + 2. * pimb[i] * pdata[i]
145 + 3. * pimc[i] * pow2(pdata[i])));
146 perr[i] = sqrt(perr[i]);
147
148 correction_factor = pima[i] + \
149 ((pimb[i] + pimc[i] * pdata[i]) * pdata[i]);
150 pdata[i] = pdata[i] * correction_factor;
151 }
152 /* return */
153 return 0 ;
154}
155
156/*----------------------------------------------------------------------------*/
174/*----------------------------------------------------------------------------*/
175int cr2res_detlin_compute(
176 const cpl_vector * dits,
177 const cpl_vector * adus,
178 cpl_size max_degree,
179 cpl_polynomial ** fitted,
180 cpl_vector ** error)
181{
182 cpl_matrix * samppos ;
183 cpl_polynomial * fitted_local ;
184 cpl_vector * error_local ;
185 cpl_vector * adusPsec ;
186 cpl_vector * y_tofit, *tmp;
187 cpl_vector * adus_loc, *dits_loc;
188 double aduPsec;
189 cpl_size i=0, first_satur=-1 ;
190 int count_linear=0;
191
192 /* Test entries */
193 if (fitted == NULL || dits == NULL || adus == NULL) return -1 ;
194 if (cpl_vector_get_size(dits) != cpl_vector_get_size(adus))
195 return -1 ;
196
197 /* Determine true ADU/s by assuming it is linear up to threshold */
198 if (cpl_vector_get_min(adus) > CR2RES_DETLIN_THRESHOLD) return -1;
199 if (cpl_vector_get_max(adus) < CR2RES_DETLIN_THRESHOLD) return -1;
200 for (i = 0; i < cpl_vector_get_size(adus); i++) {
201 if (cpl_vector_get(adus,i) < CR2RES_DETLIN_THRESHOLD ) count_linear++;
202 if (cpl_vector_get(adus,i)>CR2RES_DETLIN_MAXFIT && first_satur==-1)
203 first_satur=i;
204 }
205
206 // in case no adus above threshold
207 if (first_satur==-1) first_satur= cpl_vector_get_size(dits);
208 // in case, we have more linear points than the saturation
209 if (count_linear > first_satur-1) count_linear = first_satur-1;
210
211 adus_loc = cpl_vector_extract(adus, 0, first_satur-1, 1);
212 dits_loc = cpl_vector_extract(dits, 0, first_satur-1, 1);
213
214
215 adusPsec = cpl_vector_duplicate(adus_loc);
216 cpl_vector_divide(adusPsec, dits_loc);
217 tmp = cpl_vector_extract(adusPsec,0,count_linear-1,1);
218 aduPsec = cpl_vector_get_median(tmp);
219 cpl_vector_delete(tmp);
220
221 samppos = cpl_matrix_wrap(1,
222 cpl_vector_get_size(adus_loc),
223 cpl_vector_get_data(adus_loc)) ;
224
225 y_tofit = cpl_vector_new(cpl_vector_get_size(dits_loc));
226 for(i = 0; i < cpl_vector_get_size(dits_loc); i++) {
227 double y;
228 // We fit the ratio of true ADU/s over the measured ones.
229 y = aduPsec / cpl_vector_get(adusPsec,i);
230 cpl_vector_set(y_tofit, i, y);
231 }
232 cpl_vector_delete(adusPsec);
233
234 /* Fit */
235 fitted_local = cpl_polynomial_new(1);
236 if (
237 cpl_polynomial_fit(fitted_local, samppos, NULL, y_tofit, NULL,
238 CPL_FALSE, NULL, &max_degree) != CPL_ERROR_NONE) {
239
240 /* Failed Fit - Fill the coefficients */
241 cpl_vector_delete(adus_loc);
242 cpl_vector_delete(dits_loc);
243 cpl_matrix_unwrap(samppos) ;
244 cpl_vector_delete(y_tofit);
245 cpl_polynomial_delete(fitted_local) ;
246 cpl_error_reset() ;
247 return -1 ;
248 }
249
250 /* Compute the error */
251 error_local = cpl_vector_new(max_degree+1) ;
252 cpl_size nc = max_degree + 1;
253 cpl_size ndata = cpl_vector_get_size(y_tofit);
254 if (nc >= ndata){
255 // No uncertainty as fit should be perfectly aligned with data points
256 for (i = 0; i < max_degree + 1; i++) {
257 cpl_vector_set(error_local, i, 0);
258 }
259 } else {
260 // lhs = vandermode(x, order)
261 // hankel = dot(lhs.T, lhs)
262 // cov = inv(hankel)
263 // cov *= sum(abs(resids)) / (len(x) - order)
264 // error_local = sqrt(diag(cov))
265
266 cpl_matrix * hankel = cpl_matrix_new(nc, nc);
267 cpl_matrix * mx = cpl_matrix_new(nc, 1); // just a temporary matrix
268
269 // this actually returns the hankel matrix, not vandermode
270 // also directly copied from the cpl source code (cpl_polynomial.c)
271 cr2res_matrix_fill_normal_vandermonde(hankel, mx, adus_loc, CPL_FALSE,
272 0, y_tofit);
273 cpl_matrix * inverse = cpl_matrix_invert_create(hankel);
274 cpl_vector * resids = cpl_vector_new(ndata);
275
276 cpl_vector_fill_polynomial_fit_residual(resids, y_tofit, NULL,
277 fitted_local, samppos, NULL);
278 for (i = 0; i < ndata; i++) {
279 cpl_vector_set(resids, i, fabs(cpl_vector_get(resids, i)));
280 }
281
282 cpl_matrix_multiply_scalar(inverse,
283 cpl_vector_get_sum(resids) / (double)(ndata - nc));
284
285 for (i = 0; i < max_degree + 1; i++) {
286 cpl_vector_set(error_local, i, sqrt(cpl_matrix_get(inverse, i, i)));
287 }
288 cpl_matrix_delete(hankel);
289 cpl_matrix_delete(mx);
290 cpl_matrix_delete(inverse);
291 cpl_vector_delete(resids);
292 }
293 cpl_vector_delete(y_tofit);
294 cpl_matrix_unwrap(samppos) ;
295 cpl_vector_delete(dits_loc);
296 cpl_vector_delete(adus_loc);
297
298 /* Check Result - Polynomial coefficients are NaN sometimes */
299 for (i=0 ; i<=max_degree ; i++) {
300 double cur_coeff;
301 cur_coeff = cpl_polynomial_get_coeff(fitted_local, &i) ;
302 if (isnan(cur_coeff)) {
303 cpl_polynomial_delete(fitted_local) ;
304 cpl_vector_delete(error_local) ;
305 *fitted = NULL ;
306 *error = NULL ;
307 return -1 ;
308 }
309 }
310
311 /* Catch an error in CPL */
312 if (cpl_error_get_code()) {
313 cpl_error_reset() ;
314 cpl_polynomial_delete(fitted_local) ;
315 cpl_vector_delete(error_local) ;
316 *fitted = NULL ;
317 *error = NULL ;
318 return -1 ;
319 }
320
321 /* Return */
322 *fitted = fitted_local ;
323 if (error != NULL) *error = error_local ;
324 else cpl_vector_delete(error_local) ;
325 return 0 ;
326}
327
328/*----------------------------------------------------------------------------*/
334/*----------------------------------------------------------------------------*/
335cpl_frameset * cr2res_detlin_sort_frames(
336 const cpl_frameset * in)
337{
338 cpl_frameset * sorted ;
339
340 /* Check Inputs */
341 if (in == NULL) return NULL ;
342
343 sorted = cpl_frameset_duplicate(in) ;
344 if (cpl_frameset_sort(sorted,
345 cr2res_detlin_frames_dit_compare) != CPL_ERROR_NONE) {
346 cpl_frameset_delete(sorted) ;
347 return NULL ;
348 }
349 return sorted ;
350}
351
353
354/*----------------------------------------------------------------------------*/
361/*----------------------------------------------------------------------------*/
362static int cr2res_detlin_frames_dit_compare(
363 const cpl_frame * in1,
364 const cpl_frame * in2)
365{
366 cpl_propertylist * plist1 ;
367 cpl_propertylist * plist2 ;
368 double dit1, dit2 ;
369
370 /* Test entries */
371 if (in1==NULL || in2==NULL) return 0 ;
372
373 /* Get property lists */
374 plist1=cpl_propertylist_load(cpl_frame_get_filename(in1),0) ;
375 plist2=cpl_propertylist_load(cpl_frame_get_filename(in2),0) ;
376
377 /* Get DITs */
378 dit1 = cr2res_pfits_get_dit(plist1) ;
379 dit2 = cr2res_pfits_get_dit(plist2) ;
380 if (plist1 != NULL) cpl_propertylist_delete(plist1) ;
381 if (plist2 != NULL) cpl_propertylist_delete(plist2) ;
382 if (cpl_error_get_code()) return 0 ;
383
384 if (dit1<dit2) return -1 ;
385 else if (dit1>dit2) return 1 ;
386 return 0 ;
387}
388
389/*----------------------------------------------------------------------------*/
403/*----------------------------------------------------------------------------*/
404static void cr2res_matrix_fill_normal_vandermonde(cpl_matrix * self,
405 cpl_matrix * mx,
406 const cpl_vector * xhat,
407 cpl_boolean is_eqdist,
408 cpl_size mindeg,
409 const cpl_vector * values)
410{
411
412
413 const double * dval = cpl_vector_get_data_const(values);
414 const double * xval = cpl_vector_get_data_const(xhat);
415 cpl_vector * phat = cpl_vector_duplicate(xhat); /* Powers of xhat */
416 cpl_vector * qhat = NULL; /* mindeg Power of xhat */
417 double * dhat = cpl_vector_get_data(phat);
418 double * ehat = NULL;
419 const cpl_size nc = cpl_matrix_get_ncol(self);
420 const cpl_size np = cpl_vector_get_size(xhat);
421 cpl_size i,j;
422
423 /* Fill Hankel matrix from top-left to main skew diagonal
424 - on and above (non-skew) main diagonal */
425 /* Also compute transpose(V) * b */
426 /* Peel off 1st iteration */
427 if (mindeg > 0) {
428 double hsum = 0.0;
429 cpl_size k;
430
431 qhat = mindeg == 1 ? cpl_vector_duplicate(xhat) : cpl_vector_new(np);
432 ehat = cpl_vector_get_data(qhat);
433
434 /* Raise xhat to the power of mindeg */
435 for (k=0; k < np; k++) {
436 const double x = xval[k];
437
438 if (mindeg > 1) ehat[k] = pow(x, (int)mindeg);
439 dhat[k] *= ehat[k];
440
441 hsum += ehat[k] * ehat[k];
442 }
443 cpl_matrix_set(self, 0, 0, hsum);
444 } else {
445 cpl_matrix_set(self, 0, 0, (double)np);
446 }
447 /* qhat is xhat to the power of mindeg, iff mindeg > 0 */
448 /* dhat is xhat to the power of 1+mindeg, iff mindeg > 0 */
449 for (j=1; j < 2; j++) {
450 double vsum0 = 0.0;
451 double hsum = 0.0;
452 double vsum = 0.0;
453 cpl_size k;
454
455 for (k=0; k < np; k++) {
456 const double y = dval[k];
457
458 hsum += mindeg > 0 ? ehat[k] * dhat[k] : dhat[k];
459 vsum += y * dhat[k];
460 vsum0 += mindeg > 0 ? ehat[k] * y : y;
461 }
462 cpl_matrix_set(mx, 0, 0, vsum0);
463 cpl_matrix_set(mx, j, 0, vsum);
464 if (is_eqdist) continue;
465 k = j;
466 for (i=0; i <= k; i++, k--) {
467 cpl_matrix_set(self, i, k, hsum);
468 }
469 }
470 for (; j < nc; j++) {
471 double hsum = 0.0;
472 double vsum = 0.0;
473 cpl_size k;
474
475 for (k=0; k < np; k++) {
476 const double x = xval[k];
477 const double y = dval[k];
478
479 dhat[k] *= x;
480 hsum += mindeg > 0 ? ehat[k] * dhat[k] : dhat[k];
481 vsum += y * dhat[k];
482 }
483 cpl_matrix_set(mx, j, 0, vsum);
484 if (is_eqdist && (j&1)) continue;
485 k = j;
486 for (i=0; i <= k; i++, k--) {
487 cpl_matrix_set(self, i, k, hsum);
488 }
489 }
490 /* Fill remaining Hankel matrix - on and above (non-skew) main diagonal */
491
492 if (mindeg > 0) {
493 cpl_vector_multiply(phat, qhat);
494 cpl_vector_delete(qhat);
495 }
496
497 for (i = 1; i < nc; i++) {
498 cpl_size k;
499 double hsum = 0.0;
500
501 if (is_eqdist && ((i+nc)&1)==0) {
502 cpl_vector_multiply(phat, xhat);
503 continue;
504 }
505
506 for (k=0; k < np; k++) {
507 const double x = xval[k];
508
509 dhat[k] *= x;
510 hsum += dhat[k];
511 }
512 k = i;
513 for (j = nc-1; k <= j; k++, j--) {
514 cpl_matrix_set(self, k, j, hsum);
515 }
516 }
517
518 cpl_vector_delete(phat);
519
520}
521
cr2res_detlin_sort_frames
cpl_frameset * cr2res_detlin_sort_frames(const cpl_frameset *in)
Sort the frames by increasing DIT.
Definition cr2res_detlin.c:335
cr2res_detlin_compute
int cr2res_detlin_compute(const cpl_vector *dits, const cpl_vector *adus, cpl_size max_degree, cpl_polynomial **fitted, cpl_vector **error)
Fits the response of a given pixel to the illumination increase.
Definition cr2res_detlin.c:175
cr2res_detlin_correct
int cr2res_detlin_correct(hdrl_image *in, const hdrl_imagelist *detlin)
Apply the detector linearity correction.
Definition cr2res_detlin.c:71
cr2res_pfits_get_dit
double cr2res_pfits_get_dit(const cpl_propertylist *plist)
find out the DIT value
Definition cr2res_pfits.c:199
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_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_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_imagelist_get_const
const hdrl_image * hdrl_imagelist_get_const(const hdrl_imagelist *himlist, cpl_size inum)
Get an image from a list of images.
Definition hdrl_imagelist_io.c:229
Generated on Tue Jun 16 2026 00:11:34 for CR2RE Pipeline Reference Manual by doxygen 1.13.2