CR2RE Pipeline Reference Manual 1.7.0
irplib_distortion.c
1/* $Id: irplib_distortion.c,v 1.52 2013-01-29 08:43:33 jtaylor Exp $
2 *
3 * This file is part of the irplib package
4 * Copyright (C) 2002,2003 European Southern Observatory
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02111-1307 USA
19 */
20
21/*
22 * $Author: jtaylor $
23 * $Date: 2013-01-29 08:43:33 $
24 * $Revision: 1.52 $
25 * $Name: not supported by cvs2svn $
26 */
27
28#ifdef HAVE_CONFIG_H
29#include <config.h>
30#endif
31
32/*-----------------------------------------------------------------------------
33 Includes
34 -----------------------------------------------------------------------------*/
35
36#include "irplib_distortion.h"
37
38#include "irplib_flat.h"
39#include "irplib_utils.h"
40#include "irplib_polynomial.h"
41
42#include <math.h>
43#include <float.h>
44
45/*-----------------------------------------------------------------------------
46 Define
47 -----------------------------------------------------------------------------*/
48
49#define IRPLIB_MAX(A, B) ((A) > (B) ? (A) : (B))
50#define IRPLIB_MIN(A, B) ((A) < (B) ? (A) : (B))
51
52#define ARC_MINGOODPIX 100
53#define ARC_MINARCLENFACT 2.0
54#define ARC_MINNBARCS 4
55#define ARC_RANGE_FACT 3.0
56#define ARC_WINDOWSIZE 32
57
58#define TRESH_MEDIAN_MIN 0.0
59#define TRESH_SIGMA_MAX 200.0
60
61/*----------------------------------------------------------------------------*/
65/*----------------------------------------------------------------------------*/
66
67/*-----------------------------------------------------------------------------
68 Functions prototypes
69 -----------------------------------------------------------------------------*/
70
71static cpl_apertures *irplib_distortion_detect_arcs(cpl_image *,
72 cpl_image **,
73 int,
74 int,
75 double,
76 int,
77 int,
78 int,
79 int);
80static cpl_error_code
81irplib_distortion_fill_border(cpl_image *, int, int, int, int, double);
82static int
83irplib_distortion_threshold1d(cpl_image *, double, cpl_image *, double);
84static cpl_error_code irplib_distortion_purge_arcs(cpl_apertures **,
85 cpl_image *,
86 const cpl_image *,
87 int,
88 int,
89 double);
90static cpl_error_code
91irplib_distortion_fill_arc_positions(cpl_bivector *,
92 cpl_vector *,
93 const cpl_image *,
94 const cpl_image *,
95 const cpl_apertures *);
96
97static double irplib_distortion_get_row_centroid(const cpl_image *,
98 const cpl_image *,
99 int,
100 int);
101
102static int irplib_distortion_sub_hor_lowpass(cpl_image *, int);
103static cpl_image *irplib_distortion_remove_ramp(const cpl_image *);
104
105static cpl_error_code irplib_image_filter_background_line(cpl_image *,
106 const cpl_image *,
107 int,
108 cpl_boolean);
109
110static cpl_error_code irplib_polynomial_fit_2d(cpl_polynomial *,
111 const cpl_bivector *,
112 const cpl_vector *,
113 int,
114 double,
115 double *);
116
117static cpl_matrix *irplib_matrix_product_normal_create(const cpl_matrix *);
118
119/*-----------------------------------------------------------------------------
120 Functions code
121 -----------------------------------------------------------------------------*/
122
124
125/*----------------------------------------------------------------------------*/
154/*----------------------------------------------------------------------------*/
155cpl_polynomial *
156irplib_distortion_estimate(const cpl_image *org,
157 int xmin,
158 int ymin,
159 int xmax,
160 int ymax,
161 int auto_ramp_sub,
162 int arc_sat,
163 int max_arc_width,
164 double kappa,
165 int degree,
166 cpl_apertures **arcs)
167{
168 cpl_image *local_im;
169 cpl_image *filtered;
170 cpl_image *label_image;
171 double rightmost, leftmost;
172 cpl_bivector *grid;
173 cpl_vector *values_to_fit;
174 int n_arcs;
175 cpl_polynomial *poly2d;
176 double mse = 0.0;
177 const int nx = cpl_image_get_size_x(org);
178 const int ny = cpl_image_get_size_y(org);
179 const int min_arc_range = (int)(nx / ARC_RANGE_FACT);
180 int i;
181
182 /* Check entries */
183 cpl_ensure(org != NULL, CPL_ERROR_NULL_INPUT, NULL);
184 cpl_ensure(kappa >= 0.0, CPL_ERROR_ILLEGAL_INPUT, NULL);
185 cpl_ensure(max_arc_width > 0, CPL_ERROR_ILLEGAL_INPUT, NULL);
186
187 /* The background may vary strongly along the vertical line. */
188 /* Detect and rm background with a 1+2*max_arc_width x 1 median filter */
189
190 filtered = cpl_image_new(nx, ny, cpl_image_get_type(org));
191
192 irplib_image_filter_background_line(filtered, org, max_arc_width, CPL_TRUE);
193
194 if (auto_ramp_sub) {
195 local_im = irplib_distortion_remove_ramp(filtered);
196 cpl_image_delete(filtered);
197 }
198 else {
199 local_im = filtered;
200 }
201
202 cpl_error_ensure(local_im != NULL, cpl_error_get_code(), return (NULL),
203 "Cannot clean the image");
204
205 /* Detect the arcs in the input image */
206 *arcs = irplib_distortion_detect_arcs(local_im, &label_image, arc_sat,
207 max_arc_width, kappa, xmin, ymin,
208 xmax, ymax);
209 if (*arcs == NULL) {
210 cpl_image_delete(local_im);
211 cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
212 "Cannot detect the arcs");
213 return NULL;
214 }
215 n_arcs = cpl_apertures_get_size(*arcs);
216 cpl_msg_info(cpl_func, "%d detected arcs", n_arcs);
217
218 /* Check that the arcs are not concentrated in the same zone */
219 rightmost = leftmost = cpl_apertures_get_pos_x(*arcs, 1);
220 for (i = 1; i < n_arcs; i++) {
221 if (cpl_apertures_get_pos_x(*arcs, i + 1) < leftmost)
222 leftmost = cpl_apertures_get_pos_x(*arcs, i + 1);
223 if (cpl_apertures_get_pos_x(*arcs, i + 1) > rightmost)
224 rightmost = cpl_apertures_get_pos_x(*arcs, i + 1);
225 }
226 if ((int)(rightmost - leftmost) < min_arc_range) {
227#if defined CPL_HAVE_VA_ARGS && CPL_HAVE_VA_ARGS != 0
228 cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
229 "too narrow range (%g-%g)<%d", rightmost,
230 leftmost, min_arc_range);
231#else
232 cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
233 "too narrow range");
234#endif
235 cpl_apertures_delete(*arcs);
236 cpl_image_delete(local_im);
237 cpl_image_delete(label_image);
238 *arcs = NULL;
239 return NULL;
240 }
241
242 /* Create a 2-D deformation grid with detected arcs */
243 cpl_msg_info(cpl_func, "Create deformation grid");
244 grid = cpl_bivector_new(n_arcs * ny);
245 values_to_fit = cpl_vector_new(n_arcs * ny);
246
247 if (irplib_distortion_fill_arc_positions(grid, values_to_fit, local_im,
248 label_image, *arcs)) {
249 cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
250 "cannot get arcs positions");
251 cpl_apertures_delete(*arcs);
252 cpl_image_delete(local_im);
253 cpl_image_delete(label_image);
254 *arcs = NULL;
255 return NULL;
256 }
257 cpl_image_delete(label_image);
258 cpl_image_delete(local_im);
259
260 /* Apply the fitting */
261 poly2d = cpl_polynomial_new(2);
262 if (irplib_polynomial_fit_2d(poly2d, grid, values_to_fit, degree,
263 0.5 * (ny + 1), &mse)) {
264 cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
265 "cannot apply the 2d fit");
266 cpl_bivector_delete(grid);
267 cpl_vector_delete(values_to_fit);
268 cpl_apertures_delete(*arcs);
269 *arcs = NULL;
270 return NULL;
271 }
272
273 cpl_msg_info(cpl_func,
274 "Fitted a %d. degree 2D-polynomial to %" CPL_SIZE_FORMAT
275 " points "
276 "with mean-square error: %g",
277 degree, cpl_vector_get_size(values_to_fit), mse);
278
279 /* Free and return */
280 cpl_bivector_delete(grid);
281 cpl_vector_delete(values_to_fit);
282 return poly2d;
283}
284
286
287/*----------------------------------------------------------------------------*/
303/*----------------------------------------------------------------------------*/
304static cpl_apertures *
305irplib_distortion_detect_arcs(cpl_image *im,
306 cpl_image **label_im,
307 int arc_sat,
308 int max_arc_width,
309 double kappa,
310 int xmin,
311 int ymin,
312 int xmax,
313 int ymax)
314{
315 const int ny = cpl_image_get_size_y(im);
316 /* Set min_arclen */
317 const int min_arclen = (int)(ny / ARC_MINARCLENFACT);
318 cpl_image *filt_im;
319 cpl_mask *filter;
320 cpl_image *collapsed;
321 cpl_mask *bin_im;
322 double threshold, fillval, median_val, sigma;
323 cpl_apertures *det;
324 cpl_size nobj;
325 int ngoodpix;
326
327 /* Default values for output parameters */
328 *label_im = NULL;
329
330 /* Clear zones to be ignored (to avoid false detections) */
331 median_val = cpl_image_get_median_dev(im, &sigma);
332 fillval = median_val - sigma / 2.0;
333 if (irplib_distortion_fill_border(im, xmin, ymin, xmax, ymax, fillval)) {
334 cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
335 "cannot fill bad zones");
336 return NULL;
337 }
338
339 /* Subtract a low-pass */
340 filt_im = cpl_image_duplicate(im);
341 if (irplib_distortion_sub_hor_lowpass(filt_im, ARC_WINDOWSIZE) == -1) {
342 cpl_image_delete(filt_im);
343 return NULL;
344 }
345
346 /* Get relevant stats for thresholding */
347 median_val = cpl_image_get_median_dev(filt_im, &sigma);
348
349 /* Correct median_val and sigma if necessary */
350 if (median_val < TRESH_MEDIAN_MIN)
351 median_val = TRESH_MEDIAN_MIN;
352 if (sigma > TRESH_SIGMA_MAX)
353 sigma = TRESH_SIGMA_MAX;
354
355 /* Set the threshold */
356 threshold = median_val + sigma * kappa;
357
358 /* Collapse the image */
359 collapsed = cpl_image_collapse_median_create(filt_im, 0, 0, 0);
360
361 /* Threshold to keep only the arcs - use of the collapsed image */
362 if (irplib_distortion_threshold1d(filt_im, median_val, collapsed, 0.0) ==
363 -1) {
364 cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
365 "cannot threshold the filtered image");
366 cpl_image_delete(filt_im);
367 cpl_image_delete(collapsed);
368 return NULL;
369 }
370 cpl_image_delete(collapsed);
371
372 /* Binarize the image */
373 bin_im = cpl_mask_threshold_image_create(filt_im, threshold, DBL_MAX);
374 cpl_image_delete(filt_im);
375 if (bin_im == NULL) {
376 cpl_error_set_message(cpl_func, CPL_ERROR_DATA_NOT_FOUND,
377 "cannot binarise the image");
378 return NULL;
379 }
380
381 /* Test if there are enough good pixels */
382 ngoodpix = cpl_mask_count(bin_im);
383 if (ngoodpix < ARC_MINGOODPIX) {
384#if defined CPL_HAVE_VA_ARGS && CPL_HAVE_VA_ARGS != 0
385 cpl_error_set_message(cpl_func, CPL_ERROR_DATA_NOT_FOUND,
386 "Too few (%d) white pixels", ngoodpix);
387#else
388 cpl_error_set_message(cpl_func, CPL_ERROR_DATA_NOT_FOUND,
389 "Too few white pixels");
390#endif
391 cpl_mask_delete(bin_im);
392 return NULL;
393 }
394
395 /* Apply a morphological opening to clean the isolated pixels */
396 filter = cpl_mask_new(3, 3);
397 cpl_mask_not(filter);
398 cpl_mask_filter(bin_im, bin_im, filter, CPL_FILTER_OPENING,
399 CPL_BORDER_ZERO);
400 cpl_mask_delete(filter);
401
402 /* Labelize pixel map to a label image */
403 *label_im = cpl_image_labelise_mask_create(bin_im, &nobj);
404 cpl_mask_delete(bin_im);
405
406 /* Compute statistics on objects */
407 if ((det = cpl_apertures_new_from_image(im, *label_im)) == NULL) {
408 cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
409 "Cannot compute arcs stats");
410 cpl_image_delete(*label_im);
411 *label_im = NULL;
412 return NULL;
413 }
414
415 /* Purge non-relevant arcs */
416 if (irplib_distortion_purge_arcs(&det, *label_im, im, min_arclen,
417 max_arc_width, arc_sat)) {
418 cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
419 "Cannot purge the arcs");
420 cpl_image_delete(*label_im);
421 *label_im = NULL;
422 cpl_apertures_delete(det);
423 return NULL;
424 }
425 if (cpl_apertures_get_size(det) < ARC_MINNBARCS) {
426#if defined CPL_HAVE_VA_ARGS && CPL_HAVE_VA_ARGS != 0
427 cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
428 "Not enough valid arcs (%" CPL_SIZE_FORMAT
429 " < %d)",
430 cpl_apertures_get_size(det), ARC_MINNBARCS);
431#else
432 cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
433 "Not enough valid arcs, min=" CPL_STRINGIFY(ARC_MINNBARCS));
434#endif
435 cpl_image_delete(*label_im);
436 *label_im = NULL;
437 cpl_apertures_delete(det);
438 return NULL;
439 }
440
441 /* Return */
442 return det;
443}
444
445/*----------------------------------------------------------------------------*/
455/*----------------------------------------------------------------------------*/
456static cpl_error_code
457irplib_distortion_fill_border(cpl_image *self,
458 int xmin,
459 int ymin,
460 int xmax,
461 int ymax,
462 double fillval)
463{
464 const int nx = cpl_image_get_size_x(self);
465 const int ny = cpl_image_get_size_y(self);
466 float *pfi = cpl_image_get_data_float(self);
467 const float fvalue = (float)fillval;
468 int i, j;
469
470
471 cpl_ensure_code(pfi != NULL, cpl_error_get_code());
472
473 /* Ensure validity of pixel buffer access */
474 xmin = IRPLIB_MIN(xmin, nx + 1);
475 ymax = IRPLIB_MIN(ymax, ny);
476
477 /* - and avoid double access */
478 xmax = IRPLIB_MAX(xmax, xmin - 1);
479 ymin = IRPLIB_MIN(ymin, ymax + 1);
480
481 /* Fill the zone */
482
483 for (j = 0; j < ymin - 1; j++) {
484 for (i = 0; i < nx; i++) {
485 pfi[i + j * nx] = fvalue;
486 }
487 }
488 /* assert( j == IRPLIB_MAX(0, ymin-1) ); */
489
490 for (; j < ymax; j++) {
491 for (i = 0; i < xmin - 1; i++) {
492 pfi[i + j * nx] = fvalue;
493 }
494 for (i = xmax; i < nx; i++) {
495 pfi[i + j * nx] = fvalue;
496 }
497 }
498 /* assert( j == IRPLIB_MAX(0, ymax) ); */
499
500 for (; j < ny; j++) {
501 for (i = 0; i < nx; i++) {
502 pfi[i + j * nx] = fvalue;
503 }
504 }
505
506 return CPL_ERROR_NONE;
507}
508
509static int
510irplib_distortion_threshold1d(cpl_image *im,
511 double threshold,
512 cpl_image *im1d,
513 double newval)
514{
515 float *pim;
516 float *pim1d;
517 int nx, ny;
518 int i, j;
519
520 /* Check entries */
521 if (im == NULL)
522 return -1;
523 if (im1d == NULL)
524 return -1;
525 if (cpl_image_get_type(im) != CPL_TYPE_FLOAT)
526 return -1;
527 if (cpl_image_get_type(im1d) != CPL_TYPE_FLOAT)
528 return -1;
529
530 /* Get access to the im / im1d data */
531 pim = cpl_image_get_data_float(im);
532 pim1d = cpl_image_get_data_float(im1d);
533 nx = cpl_image_get_size_x(im);
534 ny = cpl_image_get_size_y(im);
535
536 /* Apply the thresholding */
537 for (i = 0; i < nx; i++)
538 if (pim1d[i] < threshold) {
539 for (j = 0; j < ny; j++)
540 pim[i + j * nx] = (float)newval;
541 }
542
543 /* Return */
544 return 0;
545}
546
547static int
548irplib_distortion_sub_hor_lowpass(cpl_image *im, int filt_size)
549{
550 cpl_vector *linehi;
551 cpl_vector *linelo;
552 cpl_vector *avglinehi;
553 cpl_vector *avglinelo;
554 double *pavglinehi;
555 float *pim;
556 int lopos, hipos, nx, ny;
557 int i, j;
558
559 /* Test entries */
560 if (im == NULL)
561 return -1;
562 if (filt_size <= 0)
563 return -1;
564
565 /* Initialise */
566 nx = cpl_image_get_size_x(im);
567 ny = cpl_image_get_size_y(im);
568 lopos = (int)(ny / 4);
569 hipos = (int)(3 * ny / 4);
570
571 /* Get the vectors out of the image */
572 if ((linehi = cpl_vector_new_from_image_row(im, hipos)) == NULL) {
573 return -1;
574 }
575 if ((linelo = cpl_vector_new_from_image_row(im, lopos)) == NULL) {
576 cpl_vector_delete(linehi);
577 return -1;
578 }
579
580 /* Filter the vectors */
581 if ((avglinehi = cpl_vector_filter_median_create(linehi, filt_size)) ==
582 NULL) {
583 cpl_vector_delete(linehi);
584 cpl_vector_delete(linelo);
585 return -1;
586 }
587 cpl_vector_delete(linehi);
588
589 if ((avglinelo = cpl_vector_filter_median_create(linelo, filt_size)) ==
590 NULL) {
591 cpl_vector_delete(linelo);
592 cpl_vector_delete(avglinehi);
593 return -1;
594 }
595 cpl_vector_delete(linelo);
596
597 /* Average the filtered vectors to get the low freq signal */
598 cpl_vector_add(avglinehi, avglinelo);
599 cpl_vector_delete(avglinelo);
600 cpl_vector_divide_scalar(avglinehi, 2.0);
601
602 /* Subtract the low frequency signal */
603 pavglinehi = cpl_vector_get_data(avglinehi);
604 pim = cpl_image_get_data_float(im);
605 for (i = 0; i < nx; i++) {
606 for (j = 0; j < ny; j++) {
607 pim[i + j * nx] -= pavglinehi[i];
608 }
609 }
610 cpl_vector_delete(avglinehi);
611
612 return 0;
613}
614
615/*----------------------------------------------------------------------------*/
626/*----------------------------------------------------------------------------*/
627static cpl_error_code
628irplib_distortion_purge_arcs(cpl_apertures **self,
629 cpl_image *lab_im,
630 const cpl_image *arc_im,
631 int min_arclen,
632 int max_arcwidth,
633 double arc_sat)
634{
635 const double ycenter = 0.5 * (1 + cpl_image_get_size_y(arc_im));
636 int narcs;
637 int nkeep = 0;
638 int ifirst = 1;
639 int *relabel;
640 int i;
641
642 /* Check entries */
643 cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT);
644
645 /* Get number of arcs */
646 narcs = cpl_apertures_get_size(*self);
647
648 cpl_ensure_code(narcs > 0, CPL_ERROR_DATA_NOT_FOUND);
649 cpl_ensure_code(cpl_image_get_type(lab_im) == CPL_TYPE_INT,
650 CPL_ERROR_ILLEGAL_INPUT);
651
652 /* Allocate relabel array with default relabelling to zero */
653 relabel = cpl_calloc(narcs, sizeof(int));
654
655 /* Loop on the different arcs candidates */
656 for (i = 0; i < narcs; i++) {
657 /* Test if the current object is a valid arc */
658 const int arclen = 1 + cpl_apertures_get_top(*self, i + 1) -
659 cpl_apertures_get_bottom(*self, i + 1);
660
661 if (cpl_apertures_get_top(*self, i + 1) < ycenter)
662 continue;
663 if (cpl_apertures_get_bottom(*self, i + 1) > ycenter)
664 continue;
665
666 if (arclen > min_arclen) {
667 const int arcwidth = 1 + cpl_apertures_get_right(*self, i + 1) -
668 cpl_apertures_get_left(*self, i + 1);
669 if (arcwidth < max_arcwidth) {
670 const int edge = cpl_apertures_get_left_y(*self, i + 1);
671 if (edge > 0) {
672 const double mean = cpl_apertures_get_mean(*self, i + 1);
673 if (mean < arc_sat) {
674 relabel[i] = ++nkeep;
675 /* Relabeling, if any, starts with ifirst */
676 if (nkeep == i + 1)
677 ifirst = nkeep;
678 }
679 }
680 }
681 }
682 }
683
684 if (nkeep < narcs) {
685 /* Update the labelised image by erasing non valid arcs */
686 int *plabim = cpl_image_get_data_int(lab_im);
687 const int npix =
688 cpl_image_get_size_x(lab_im) * cpl_image_get_size_y(lab_im);
689
690 if (nkeep == 0) {
691 cpl_free(relabel);
692#if defined CPL_HAVE_VA_ARGS && CPL_HAVE_VA_ARGS != 0
693 return cpl_error_set_message(cpl_func, CPL_ERROR_DATA_NOT_FOUND,
694 "All %d arc(s) are invalid", narcs);
695#else
696 return cpl_error_set_message(cpl_func, CPL_ERROR_DATA_NOT_FOUND,
697 "All arcs are invalid");
698#endif
699 }
700
701 for (i = 0; i < npix; i++) {
702 const int label = plabim[i];
703
704 if (label < 0 || label > narcs)
705 break;
706 if (label >= ifirst)
707 plabim[i] = relabel[label - 1];
708 }
709
710 if (i < npix) {
711 /* lab_im is not a valid label image */
712 cpl_free(relabel);
713 return cpl_error_set(cpl_func, plabim[i] < 0
714 ? CPL_ERROR_ILLEGAL_INPUT
715 : CPL_ERROR_INCOMPATIBLE_INPUT);
716 }
717
718 /* Purge the bad arcs */
719 cpl_apertures_delete(*self);
720 *self = cpl_apertures_new_from_image(arc_im, lab_im);
721 }
722
723 cpl_free(relabel);
724
725 cpl_msg_info(cpl_func, "Purged %d of %d arcs (1st purged=%d)",
726 narcs - nkeep, narcs, ifirst);
727
728 /* arc_im may be invalid */
729 cpl_ensure_code(*self != NULL, cpl_error_get_code());
730
731 return CPL_ERROR_NONE;
732}
733
734
735/*----------------------------------------------------------------------------*/
749/*----------------------------------------------------------------------------*/
750static cpl_error_code
751irplib_distortion_fill_arc_positions(cpl_bivector *grid,
752 cpl_vector *fitvalues,
753 const cpl_image *in,
754 const cpl_image *label_im,
755 const cpl_apertures *det)
756{
757 const int narcs = cpl_apertures_get_size(det);
758 int nfitvals = cpl_vector_get_size(fitvalues);
759 const int nx = cpl_image_get_size_x(label_im);
760 const int ny = cpl_image_get_size_y(label_im);
761 cpl_image *filt_img;
762 cpl_mask *kernel;
763 cpl_vector *gridx = cpl_bivector_get_x(grid);
764 cpl_vector *gridy = cpl_bivector_get_y(grid);
765 cpl_polynomial *dist1d;
766 cpl_matrix *dist1dx = NULL;
767 cpl_vector *dist1dy = NULL;
768 double *dgridx;
769 double *dgridy;
770 double *dfitv;
771 int ndone = 0;
772 int i, obj;
773
774 cpl_ensure_code(nfitvals > 0, CPL_ERROR_DATA_NOT_FOUND);
775 cpl_ensure_code(narcs > 0, CPL_ERROR_DATA_NOT_FOUND);
776 cpl_ensure_code(cpl_image_get_type(label_im) == CPL_TYPE_INT,
777 CPL_ERROR_TYPE_MISMATCH);
778
779 /* Ensure space for output */
780 if (nfitvals < narcs * ny) {
781 nfitvals = narcs * ny;
782 cpl_vector_set_size(fitvalues, nfitvals);
783 }
784 if (cpl_vector_get_size(gridx) < nfitvals ||
785 cpl_vector_get_size(gridy) < nfitvals) {
786 cpl_vector_set_size(gridx, nfitvals);
787 cpl_vector_set_size(gridy, nfitvals);
788 }
789
790 /* Get data after resizing */
791 dgridx = cpl_vector_get_data(gridx);
792 dgridy = cpl_vector_get_data(gridy);
793 dfitv = cpl_vector_get_data(fitvalues);
794
795 /* Median filter on input image */
796 kernel = cpl_mask_new(3, 3);
797 cpl_mask_not(kernel);
798 filt_img = cpl_image_new(nx, ny, cpl_image_get_type(in));
799 cpl_image_filter_mask(filt_img, in, kernel, CPL_FILTER_MEDIAN,
800 CPL_BORDER_FILTER);
801 cpl_mask_delete(kernel);
802
803 dist1d = cpl_polynomial_new(1);
804
805 for (obj = 0; obj < narcs; obj++) {
806 /* Find the reference X-coordinate for the arc */
807 const int *plabel_im = cpl_image_get_data_int_const(label_im);
808 const int ndist1d = cpl_apertures_get_top(det, obj + 1) -
809 cpl_apertures_get_bottom(det, obj + 1) + 1;
810 cpl_boolean sampsym = CPL_TRUE;
811 int j, prevj = 0;
812 int k = 0;
813
814 (void)cpl_matrix_unwrap(dist1dx);
815 (void)cpl_vector_unwrap(dist1dy);
816 dist1dx = cpl_matrix_wrap(1, ndist1d, dgridy + ndone);
817 dist1dy = cpl_vector_wrap(ndist1d, dfitv + ndone);
818
819 /* Find out the X coord. at all Y positions on the arc */
820
821 for (j = cpl_apertures_get_bottom(det, obj + 1) - 1;
822 j < cpl_apertures_get_top(det, obj + 1); j++) {
823 for (i = 0; i < nx; i++) {
824 if (plabel_im[i + j * nx] == obj + 1)
825 break;
826 }
827 if (i < nx) {
828 /* Found 1st pixel of aperture obj+1 in row j+1 */
829 cpl_errorstate prestate = cpl_errorstate_get();
830
831 const double x_finepos =
832 irplib_distortion_get_row_centroid(filt_img, label_im, i,
833 j);
834 if (!cpl_errorstate_is_equal(prestate)) {
835 irplib_error_recover(prestate,
836 "Could not find X-position "
837 "for line %d at y=%d (x=%d)",
838 obj + 1, j + 1, i + 1);
839 }
840 else if (x_finepos >= 0.0) {
841 cpl_matrix_set(dist1dx, 0, k, 1.0 + j);
842 cpl_vector_set(dist1dy, k, 1.0 + x_finepos);
843 if (k > 0 && j != 1 + prevj)
844 sampsym = CPL_FALSE;
845 prevj = j;
846 k++;
847 }
848 }
849 }
850 if (k > 0) {
851 double ref_xpos, grad;
852 cpl_error_code error;
853 const cpl_boolean did_drop = k != ndist1d;
854 const cpl_size mindeg = 0;
855 const cpl_size maxdeg = 2;
856
857 if (did_drop) {
858 /* Set correct size */
859 dist1dx = cpl_matrix_wrap(1, k, cpl_matrix_unwrap(dist1dx));
860 dist1dy = cpl_vector_wrap(k, cpl_vector_unwrap(dist1dy));
861 }
862
863 error = cpl_polynomial_fit(dist1d, dist1dx, &sampsym, dist1dy, NULL,
864 CPL_FALSE, &mindeg, &maxdeg);
865 if (error) {
866 cpl_msg_error(cpl_func, "1D-fit failed");
867 break;
868 }
869
870 ref_xpos = cpl_polynomial_eval_1d(dist1d, 0.5 * (ny + 1), &grad);
871
872 for (j = cpl_apertures_get_bottom(det, obj + 1) - 1;
873 j < cpl_apertures_get_top(det, obj + 1); j++) {
874 const double xpos =
875 cpl_polynomial_eval_1d(dist1d, j + 1.0, NULL);
876
877 dfitv[ndone] = xpos;
878 dgridx[ndone] = ref_xpos;
879 /* Wrapping dist1dx does _not_ take care of dgridy,
880 in case of "Could not find X-position " */
881 if (did_drop)
882 dgridy[ndone] = 1.0 + j;
883 ndone++;
884 }
885 cpl_msg_info(cpl_func, "Line %d has center gradient %g", obj + 1,
886 grad);
887 }
888 }
889
890 cpl_image_delete(filt_img);
891 cpl_polynomial_delete(dist1d);
892 (void)cpl_matrix_unwrap(dist1dx);
893 (void)cpl_vector_unwrap(dist1dy);
894
895 cpl_msg_info(cpl_func,
896 "Found %d fitting points ("
897 "expected up to %d points)",
898 ndone, nfitvals);
899
900 cpl_ensure_code(obj == narcs, cpl_error_get_code());
901
902 cpl_ensure_code(ndone > 0, CPL_ERROR_DATA_NOT_FOUND);
903
904 cpl_vector_set_size(fitvalues, ndone);
905 cpl_vector_set_size(gridx, ndone);
906 cpl_vector_set_size(gridy, ndone);
907
908 return CPL_ERROR_NONE;
909}
910
911/*----------------------------------------------------------------------------*/
921/*----------------------------------------------------------------------------*/
922static double
923irplib_distortion_get_row_centroid(const cpl_image *im,
924 const cpl_image *label_im,
925 int x,
926 int y)
927{
928 const int nx = cpl_image_get_size_x(im);
929 const int ny = cpl_image_get_size_y(im);
930 const int ynx = y * nx;
931 const float *pim = cpl_image_get_data_float_const(im);
932 const int *plabel_im = cpl_image_get_data_int_const(label_im);
933 int firstpos = -1;
934 int lastpos = -1;
935 int maxpos = x;
936 int objnum;
937 double wsum = 0.0;
938 double sum = 0.0;
939 double max = 0.0;
940
941 cpl_ensure(pim != NULL, cpl_error_get_code(), -1.0);
942 cpl_ensure(plabel_im != NULL, cpl_error_get_code(), -2.0);
943 cpl_ensure(x >= 0, CPL_ERROR_ILLEGAL_INPUT, -3.0);
944 cpl_ensure(y >= 0, CPL_ERROR_ILLEGAL_INPUT, -4.0);
945 cpl_ensure(x < nx, CPL_ERROR_ILLEGAL_INPUT, -5.0);
946 cpl_ensure(y < ny, CPL_ERROR_ILLEGAL_INPUT, -6.0);
947
948 max = (double)pim[x + ynx];
949 objnum = plabel_im[x + ynx];
950
951 /* While we stay in the same object... */
952 do {
953 const double val = (double)pim[x + ynx];
954
955 if (val > 0.0) { /* FIXME: Handle this exception better */
956 wsum += x * val;
957 sum += val;
958
959 if (firstpos < 0)
960 firstpos = x;
961 lastpos = x;
962
963 if (val > max) {
964 max = val;
965 maxpos = x;
966 }
967 }
968
969
970 /* Next point */
971 x++;
972
973 } while (x < nx && objnum == plabel_im[x + ynx]);
974
975 cpl_ensure(sum > 0.0, CPL_ERROR_DATA_NOT_FOUND, -7.0);
976
977 /*
978 assert( 0 <= maxpos && maxpos < nx );
979 assert( objnum == plabel_im[maxpos + ynx] );
980 assert( wsum >= 0.0 );
981 */
982
983 return (wsum < sum * firstpos || wsum > sum * lastpos) ? maxpos
984 : wsum / sum;
985}
986
987/*----------------------------------------------------------------------------*/
993/*----------------------------------------------------------------------------*/
994#define IS_NB_TESTPOINTS 8
995#define IS_MIN_SLOPE 0.01
996#define IS_MAX_SLOPE_DIF 0.075
997#define IS_MAX_FIT_EDGE_DIF 0.05
998#define IS_MIN_RAMP 10.0
999#define IS_MAX_MNERR 13.0
1000#define IS_MAX_MNERR_DIF 8.0
1001#define IS_MAX_INTER_DIF 20.0
1002#define IS_SKIPZONE 2.5
1003#define SQR(x) ((x) * (x))
1004static cpl_image *
1005irplib_distortion_remove_ramp(const cpl_image *in)
1006{
1007 int ramp_present;
1008 const int nx = cpl_image_get_size_x(in);
1009 const int ny = cpl_image_get_size_y(in);
1010 const int yhi = (int)(ny / 2);
1011 const int ylo = yhi - 1;
1012 cpl_bivector *testpointlo;
1013 double *testpointlo_x;
1014 double *testpointlo_y;
1015 cpl_bivector *testpointhi;
1016 double *testpointhi_x;
1017 double *testpointhi_y;
1018 const int spacing = ny / (IS_SKIPZONE * IS_NB_TESTPOINTS);
1019 double rampdif, fitslope;
1020 double *pol_coefhi, *pol_coeflo;
1021 cpl_vector *median;
1022 double *median_data;
1023 double medianerrlo, medianerrhi;
1024 double slope;
1025 cpl_image *out;
1026 float *pout;
1027 int i;
1028
1029 cpl_ensure(cpl_image_get_type(in) == CPL_TYPE_FLOAT,
1030 CPL_ERROR_UNSUPPORTED_MODE, NULL);
1031
1032 if (ny < IS_SKIPZONE * IS_NB_TESTPOINTS) {
1033#if defined CPL_HAVE_VA_ARGS && CPL_HAVE_VA_ARGS != 0
1034 cpl_error_set_message(
1035 cpl_func, CPL_ERROR_ILLEGAL_INPUT,
1036 "image has %d lines, min=" CPL_STRINGIFY(IS_SKIPZONE) "*" CPL_STRINGIFY(IS_NB_TESTPOINTS), ny);
1037#else
1038 cpl_error_set_message(
1039 cpl_func, CPL_ERROR_ILLEGAL_INPUT,
1040 "image has too few lines, min=" CPL_STRINGIFY(IS_SKIPZONE) "*" CPL_STRINGIFY(IS_NB_TESTPOINTS));
1041#endif
1042 return NULL;
1043 }
1044
1045 /* Fill the vectors */
1046 testpointhi = cpl_bivector_new(IS_NB_TESTPOINTS);
1047 testpointhi_x = cpl_bivector_get_x_data(testpointhi);
1048 testpointhi_y = cpl_bivector_get_y_data(testpointhi);
1049 testpointlo = cpl_bivector_new(IS_NB_TESTPOINTS);
1050 testpointlo_x = cpl_bivector_get_x_data(testpointlo);
1051 testpointlo_y = cpl_bivector_get_y_data(testpointlo);
1052 for (i = 0; i < IS_NB_TESTPOINTS; i++) {
1053 int y;
1054 cpl_vector *tmp_vector;
1055 y = yhi + i * spacing;
1056 tmp_vector = cpl_vector_new_from_image_row(in, y + 1);
1057 testpointhi_x[i] = y - ny / 2;
1058 testpointhi_y[i] = cpl_vector_get_median_const(tmp_vector);
1059 cpl_vector_delete(tmp_vector);
1060 y = ylo - i * spacing;
1061 tmp_vector = cpl_vector_new_from_image_row(in, y + 1);
1062 testpointlo_x[IS_NB_TESTPOINTS - i - 1] = y;
1063 testpointlo_y[IS_NB_TESTPOINTS - i - 1] =
1064 cpl_vector_get_median_const(tmp_vector);
1065 cpl_vector_delete(tmp_vector);
1066 }
1067
1068 /* Apply the fit */
1069 pol_coefhi = irplib_flat_fit_slope_robust(testpointhi_x, testpointhi_y,
1070 IS_NB_TESTPOINTS);
1071 pol_coeflo = irplib_flat_fit_slope_robust(testpointlo_x, testpointlo_y,
1072 IS_NB_TESTPOINTS);
1073
1074 /* Compute the errors */
1075 median = cpl_vector_new(IS_NB_TESTPOINTS);
1076 median_data = cpl_vector_get_data(median);
1077 for (i = 0; i < IS_NB_TESTPOINTS; i++) {
1078 median_data[i] = SQR(testpointhi_y[i] - pol_coefhi[0] -
1079 pol_coefhi[1] * testpointhi_x[i]);
1080 }
1081 medianerrhi = cpl_vector_get_median(median);
1082 for (i = 0; i < IS_NB_TESTPOINTS; i++) {
1083 median_data[i] = SQR(testpointlo_y[i] - pol_coeflo[0] -
1084 pol_coeflo[1] * testpointlo_x[i]);
1085 }
1086 medianerrlo = cpl_vector_get_median(median);
1087 cpl_vector_delete(median);
1088 rampdif = testpointlo_y[IS_NB_TESTPOINTS - 1] - testpointhi_y[0];
1089 slope = rampdif / (ny / 2.0);
1090 fitslope = (pol_coefhi[1] + pol_coeflo[1]) / 2.0;
1091
1092 cpl_bivector_delete(testpointlo);
1093 cpl_bivector_delete(testpointhi);
1094
1095 /* Decide if there is a ramp or not */
1096 if (fabs(rampdif) < IS_MIN_RAMP || fabs(pol_coefhi[1]) < IS_MIN_SLOPE ||
1097 fabs(pol_coeflo[1]) < IS_MIN_SLOPE ||
1098 pol_coefhi[1] / pol_coeflo[1] < 0.5 ||
1099 pol_coefhi[1] / pol_coeflo[1] > 2.0 ||
1100 fabs(pol_coefhi[1] - pol_coeflo[1]) > IS_MAX_SLOPE_DIF ||
1101 fabs(pol_coefhi[0] - pol_coeflo[0]) > IS_MAX_INTER_DIF ||
1102 medianerrlo > IS_MAX_MNERR || medianerrhi > IS_MAX_MNERR ||
1103 fabs(medianerrlo - medianerrhi) > IS_MAX_MNERR_DIF ||
1104 fabs(slope - fitslope) > IS_MAX_FIT_EDGE_DIF ||
1105 slope / fitslope < 0.5 || slope / fitslope > 2.0)
1106 ramp_present = 0;
1107 else
1108 ramp_present = 1;
1109
1110 cpl_free(pol_coeflo);
1111 cpl_free(pol_coefhi);
1112
1113 /* Correct the ramp if it is there */
1114 out = cpl_image_duplicate(in);
1115 pout = cpl_image_get_data_float(out);
1116 if (ramp_present == 1) {
1117 float val;
1118 int j;
1119 for (j = 0; j < ny / 2; j++) {
1120 val = slope * (j - ny / 2);
1121 for (i = 0; i < nx; i++)
1122 pout[i + j * nx] -= val;
1123 }
1124 for (j = ny / 2; j < ny; j++) {
1125 val = slope * (j - ny);
1126 for (i = 0; i < nx; i++)
1127 pout[i + j * nx] -= val;
1128 }
1129 }
1130
1131 return out;
1132}
1133
1134/*----------------------------------------------------------------------------*/
1148/*----------------------------------------------------------------------------*/
1149static cpl_error_code
1150irplib_image_filter_background_line(cpl_image *self,
1151 const cpl_image *other,
1152 int hsize,
1153 cpl_boolean vertical)
1154{
1155 const int nx = cpl_image_get_size_x(self);
1156 const int ny = cpl_image_get_size_y(self);
1157 const int msize = 1 + 2 * hsize;
1158 cpl_mask *mask;
1159 cpl_image *background;
1160 cpl_error_code error = CPL_ERROR_NONE;
1161
1162 cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT);
1163 cpl_ensure_code(hsize >= 0, CPL_ERROR_ILLEGAL_INPUT);
1164
1165 if (other == NULL)
1166 other = self;
1167
1168 mask = vertical ? cpl_mask_new(msize, 1) : cpl_mask_new(1, msize);
1169
1170 error |= cpl_mask_not(mask);
1171
1172 background = cpl_image_new(nx, ny, cpl_image_get_type(other));
1173
1174 error |= cpl_image_filter_mask(background, other, mask, CPL_FILTER_MEDIAN,
1175 CPL_BORDER_FILTER);
1176 cpl_mask_delete(mask);
1177
1178 if (self != other) {
1179 error |= cpl_image_copy(self, other, 1, 1);
1180 }
1181
1182 error |= cpl_image_subtract(self, background);
1183 cpl_image_delete(background);
1184
1185 return error ? cpl_error_set_where(cpl_func) : CPL_ERROR_NONE;
1186}
1187
1188
1213
1214static cpl_matrix *
1215irplib_matrix_product_normal_create(const cpl_matrix *self)
1216{
1217 double sum;
1218 cpl_matrix *product;
1219 const double *ai = cpl_matrix_get_data_const(self);
1220 double *bwrite;
1221 const int m = cpl_matrix_get_nrow(self);
1222 const int n = cpl_matrix_get_ncol(self);
1223 int i, j, k;
1224
1225
1226 cpl_ensure(self != NULL, CPL_ERROR_NULL_INPUT, NULL);
1227
1228#if 0
1229 /* Initialize all values to zero.
1230 This is done to avoid access of uninitilized memory, in case
1231 someone passes the matrix to for example cpl_matrix_dump(). */
1232 product = cpl_matrix_new(m, m);
1233 bwrite = cpl_matrix_get_data(product);
1234#else
1235 bwrite = (double *)cpl_malloc(m * m * sizeof(double));
1236 product = cpl_matrix_wrap(m, m, bwrite);
1237#endif
1238
1239 /* The result at (i,j) is the dot-product of i'th and j'th row */
1240 for (i = 0; i < m; i++, bwrite += m, ai += n) {
1241 const double *aj;
1242 aj = ai; /* aj points to first entry in j'th row */
1243 for (j = i; j < m; j++, aj += n) {
1244 sum = 0.0;
1245 for (k = 0; k < n; k++) {
1246 sum += ai[k] * aj[k];
1247 }
1248 bwrite[j] = sum;
1249 }
1250 }
1251
1252 return product;
1253}
1254
1255/*----------------------------------------------------------------------------*/
1269/*----------------------------------------------------------------------------*/
1270static cpl_error_code
1271irplib_polynomial_fit_2d(cpl_polynomial *self,
1272 const cpl_bivector *xy_pos,
1273 const cpl_vector *values,
1274 int degree,
1275 double fixy,
1276 double *mse)
1277{
1278 const int np = cpl_bivector_get_size(xy_pos);
1279 /* Number of unknowns to determine in one dimension */
1280 const int nc1 = 1 + degree;
1281 /* Number of unknowns to determine */
1282 /* P_{i,0} = 0, except P_{1,0} = 1 */
1283 const int nc = nc1 * (1 + nc1) / 2 - nc1;
1284 cpl_matrix *mv; /* The transpose of the Vandermonde matrix */
1285 cpl_matrix *mh; /* Block-Hankel matrix, V'*V */
1286 cpl_matrix *mb;
1287 cpl_matrix *mx;
1288#ifdef IRPLIB_DISTORTION_ASSERT
1289 /*const double * coeffs1d;*/
1290#endif
1291 double *dmv;
1292 cpl_vector *xhat;
1293 cpl_vector *yhat;
1294 cpl_vector *zhat;
1295 cpl_size powers[2];
1296 int degx, degy;
1297 int i, j;
1298 cpl_error_code error;
1299
1300
1301 cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT);
1302 cpl_ensure_code(cpl_polynomial_get_dimension(self) == 2,
1303 CPL_ERROR_INVALID_TYPE);
1304 cpl_ensure_code(np > 0, cpl_error_get_code());
1305 cpl_ensure_code(values != NULL, CPL_ERROR_NULL_INPUT);
1306
1307 cpl_ensure_code(cpl_vector_get_size(values) == np,
1308 CPL_ERROR_INCOMPATIBLE_INPUT);
1309
1310 cpl_ensure_code(degree > 0, CPL_ERROR_ILLEGAL_INPUT);
1311 cpl_ensure_code(np >= nc, CPL_ERROR_DATA_NOT_FOUND);
1312
1313 /* transform zero-point to fixy */
1314 yhat = cpl_vector_duplicate(cpl_bivector_get_y_const(xy_pos));
1315 cpl_vector_subtract_scalar(yhat, fixy);
1316
1317 /* - and ensure P(y) = y on center line */
1318 xhat = cpl_vector_duplicate(cpl_bivector_get_x_const(xy_pos));
1319 zhat = cpl_vector_duplicate(values);
1320 cpl_vector_subtract(zhat, xhat);
1321
1322 /* Initialize matrices */
1323 /* mv contains the polynomial terms in the order described */
1324 /* above in each row, for each input point. */
1325 dmv = (double *)cpl_malloc(nc * np * sizeof(double));
1326 mv = cpl_matrix_wrap(nc, np, dmv);
1327
1328 /* Has redundant FLOPs, appears to improve accuracy */
1329 for (i = 0; i < np; i++) {
1330 const double x = cpl_vector_get(xhat, i);
1331 const double y = cpl_vector_get(yhat, i);
1332 double yvalue = y;
1333 j = 0;
1334 for (degy = 1; degy <= degree; degy++) {
1335 double xvalue = 1;
1336 for (degx = 0; degx <= degree - degy; degx++, j++) {
1337 dmv[np * j + i] = xvalue * yvalue;
1338 xvalue *= x;
1339 }
1340 yvalue *= y;
1341 }
1342 /* cx_assert( j == nc ); */
1343 }
1344 cpl_vector_delete(xhat);
1345 cpl_vector_delete(yhat);
1346
1347 /* mb contains the values, it is not modified */
1348 mb = cpl_matrix_wrap(np, 1, cpl_vector_get_data(zhat));
1349
1350 /* Form the right hand side of the normal equations */
1351 mx = cpl_matrix_product_create(mv, mb);
1352
1353 cpl_matrix_unwrap(mb);
1354 cpl_vector_delete(zhat);
1355
1356 /* Form the matrix of the normal equations */
1357 mh = irplib_matrix_product_normal_create(mv);
1358 cpl_matrix_delete(mv);
1359
1360 /* Solve XA=B by a least-square solution (aka pseudo-inverse). */
1361 error = cpl_matrix_decomp_chol(mh) || cpl_matrix_solve_chol(mh, mx);
1362
1363 cpl_matrix_delete(mh);
1364
1365 if (error) {
1366 cpl_matrix_delete(mx);
1367 cpl_ensure_code(0, error);
1368 }
1369
1370 /* Store coefficients for output */
1371
1372#ifdef IRPLIB_DISTORTION_ASSERT
1373 /*coeffs1d = cpl_matrix_get_data(mx);*/
1374#endif
1375
1376 j = 0;
1377 for (degy = 1; degy <= degree; degy++) {
1378 powers[1] = degy;
1379 for (degx = 0; degx <= degree - degy; degx++, j++) {
1380 powers[0] = degx;
1381 /* cx_assert( coeffs1d[j] == cpl_matrix_get(mx, j, 0) ); */
1382 cpl_polynomial_set_coeff(self, powers, cpl_matrix_get(mx, j, 0));
1383 }
1384 }
1385 /* cx_assert( j == nc ); */
1386
1387 cpl_matrix_delete(mx);
1388
1389 /* P_{1,0} = 1 */
1390 powers[0] = 1;
1391 powers[1] = 0;
1392 cpl_polynomial_set_coeff(self, powers, 1.0);
1393
1394 /* Transform the polynomial back in Y */
1395 cpl_polynomial_shift_1d(self, 1, -fixy);
1396
1397 /* If requested, compute mean squared error */
1398 if (mse != NULL) {
1399 const cpl_vector *x_pos = cpl_bivector_get_x_const(xy_pos);
1400 const cpl_vector *y_pos = cpl_bivector_get_y_const(xy_pos);
1401 cpl_vector *x_val = cpl_vector_new(2);
1402
1403 *mse = 0;
1404 for (i = 0; i < np; i++) {
1405 double residue;
1406 cpl_vector_set(x_val, 0, cpl_vector_get(x_pos, i));
1407 cpl_vector_set(x_val, 1, cpl_vector_get(y_pos, i));
1408 /* Subtract from the true value, square, accumulate */
1409 residue =
1410 cpl_vector_get(values, i) - cpl_polynomial_eval(self, x_val);
1411 *mse += residue * residue;
1412 }
1413 cpl_vector_delete(x_val);
1414 /* Average the error term */
1415 *mse /= np;
1416 }
1417
1418 return CPL_ERROR_NONE;
1419}
irplib_flat_fit_slope_robust
double * irplib_flat_fit_slope_robust(double *x, double *y, int np)
Fit a slope to a list of points (robust fit).
Definition irplib_flat.c:194
Generated on Tue Oct 13 2026 12:27:26 for CR2RE Pipeline Reference Manual by doxygen 1.13.2