00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031
00032
00033
00034
00035
00036 #include <cpl.h>
00037 #include <math.h>
00038
00039 #include "irplib_flat.h"
00040
00041 #include "sinfoni_distortion.h"
00042
00043
00044
00045
00046
00047 #define ARC_NBSAMPLES 20
00048 #define ARC_THRESHFACT (1.0/3.0)
00049 #define ARC_MINGOODPIX 100
00050 #define ARC_MINARCLENFACT 2.0
00051 #define ARC_MAXARCWIDTH 33
00052 #define ARC_MINNBARCS 4
00053 #define ARC_RANGE_FACT 3.0
00054
00055 #define TRESH_MEDIAN_MIN 0.0
00056 #define TRESH_SIGMA_MAX 200.0
00057
00058
00062
00063
00064
00065
00066
00067
00068 static int sinfoni_distortion_fill_badzones(cpl_image *, int, int, int, int,
00069 double) ;
00070 static int sinfoni_distortion_threshold1d(cpl_image *, double, cpl_image *,
00071 double) ;
00072 static int sinfoni_distortion_purge_arcs(cpl_image *, cpl_apertures **,
00073 cpl_image **, int, int, double) ;
00074 static cpl_bivector ** sinfoni_distortion_get_arc_positions(cpl_image *,
00075 cpl_image *, cpl_apertures *, int, double **) ;
00076 static double sinfoni_distortion_fine_pos(cpl_image *, cpl_image *, int, int) ;
00077 static cpl_image * sinfoni_distortion_remove_ramp(const cpl_image *) ;
00078
00079
00080
00081
00082
00085
00106
00107 cpl_polynomial * sinfoni_distortion_estimate(
00108 const cpl_image * org,
00109 int xmin,
00110 int ymin,
00111 int xmax,
00112 int ymax,
00113 int auto_ramp_sub,
00114 int arc_sat,
00115 cpl_apertures ** arcs)
00116 {
00117 const char * fctid = "sinfoni_distortion_estimate" ;
00118 cpl_image * local_im ;
00119 cpl_image * label_image ;
00120 double rightmost, leftmost ;
00121 cpl_bivector ** arcs_pos ;
00122 double * parc_posx ;
00123 double * parc_posy ;
00124 double * lines_pos ;
00125 cpl_bivector * grid ;
00126 double * pgridx ;
00127 double * pgridy ;
00128 cpl_vector * values_to_fit ;
00129 double * pvalues_to_fit ;
00130 int min_arc_range ;
00131 int n_calib ;
00132 int n_arcs ;
00133 cpl_polynomial * poly2d ;
00134 int nx ;
00135 int i, j ;
00136
00137
00138 if (org == NULL) return NULL ;
00139
00140
00141 n_calib = ARC_NBSAMPLES ;
00142 nx = cpl_image_get_size_x(org) ;
00143
00144 if (auto_ramp_sub) {
00145 local_im = sinfoni_distortion_remove_ramp(org) ;
00146 } else {
00147
00148 local_im = cpl_image_duplicate(org) ;
00149 }
00150 if (local_im == NULL) {
00151 cpl_msg_error(fctid, "Cannot clean the image") ;
00152 return NULL ;
00153 }
00154
00155
00156 cpl_msg_info(fctid, "Detect arcs") ;
00157 if ((*arcs = sinfoni_distortion_detect_arcs(local_im,
00158 &label_image,
00159 arc_sat,
00160 xmin, ymin, xmax, ymax)) == NULL) {
00161 cpl_image_delete(local_im) ;
00162 return NULL ;
00163 }
00164 n_arcs = cpl_apertures_get_size(*arcs) ;
00165
00166
00167 rightmost = leftmost = cpl_apertures_get_max_x(*arcs, 1) ;
00168 for (i=1 ; i<n_arcs ; i++) {
00169 if (cpl_apertures_get_max_x(*arcs, i+1) < leftmost)
00170 leftmost = cpl_apertures_get_max_x(*arcs, i+1) ;
00171 if (cpl_apertures_get_max_x(*arcs, i+1) > rightmost)
00172 rightmost = cpl_apertures_get_max_x(*arcs, i+1) ;
00173 }
00174 min_arc_range = (int)(nx / ARC_RANGE_FACT) ;
00175 if ((int)(rightmost-leftmost) < min_arc_range) {
00176 cpl_msg_error(fctid, "too narrow range (%g-%g)<%d",
00177 rightmost, leftmost, min_arc_range) ;
00178 cpl_apertures_delete(*arcs) ;
00179 cpl_image_delete(local_im) ;
00180 cpl_image_delete(label_image) ;
00181 return NULL ;
00182 }
00183
00184
00185 lines_pos = cpl_malloc(n_arcs * sizeof(double)) ;
00186 if ((arcs_pos = sinfoni_distortion_get_arc_positions(local_im,
00187 label_image, *arcs, n_calib, &lines_pos))==NULL){
00188 cpl_msg_error(fctid, "cannot get arcs positions") ;
00189 cpl_apertures_delete(*arcs) ;
00190 cpl_image_delete(local_im) ;
00191 cpl_free(lines_pos) ;
00192 cpl_image_delete(label_image) ;
00193 return NULL ;
00194 }
00195 cpl_image_delete(label_image) ;
00196 cpl_image_delete(local_im) ;
00197
00198
00199 grid = cpl_bivector_new(n_arcs * n_calib) ;
00200 pgridx = cpl_bivector_get_x_data(grid) ;
00201 pgridy = cpl_bivector_get_y_data(grid) ;
00202 values_to_fit = cpl_vector_new(n_arcs * n_calib) ;
00203 pvalues_to_fit = cpl_vector_get_data(values_to_fit) ;
00204 for (i=0 ; i<n_arcs ; i++) {
00205 parc_posx = cpl_bivector_get_x_data(arcs_pos[i]) ;
00206 parc_posy = cpl_bivector_get_y_data(arcs_pos[i]) ;
00207 for (j=0 ; j<n_calib ; j++) {
00208 pgridx[j+i*n_calib] = lines_pos[i] ;
00209 pgridy[j+i*n_calib] = parc_posy[j] ;
00210 pvalues_to_fit[j+i*n_calib] = parc_posx[j] ;
00211 }
00212 }
00213 for (i=0 ; i<n_arcs ; i++) cpl_bivector_delete(arcs_pos[i]) ;
00214 cpl_free(arcs_pos) ;
00215 cpl_free(lines_pos) ;
00216
00217
00218 if ((poly2d = cpl_polynomial_fit_2d_create(grid, values_to_fit, 2,
00219 NULL))==NULL) {
00220 cpl_msg_error(fctid, "cannot apply the 2d fit") ;
00221 cpl_bivector_delete(grid) ;
00222 cpl_vector_delete(values_to_fit) ;
00223 cpl_apertures_delete(*arcs) ;
00224 return NULL ;
00225 }
00226
00227
00228 cpl_bivector_delete(grid) ;
00229 cpl_vector_delete(values_to_fit) ;
00230 return poly2d ;
00231 }
00232
00233
00247
00248 cpl_apertures * sinfoni_distortion_detect_arcs(
00249 cpl_image * im,
00250 cpl_image ** label_im,
00251 int arc_sat,
00252 int xmin,
00253 int ymin,
00254 int xmax,
00255 int ymax)
00256 {
00257 const char * fctid = "sinfoni_distortion_detect_arcs" ;
00258 cpl_image * filt_im ;
00259 cpl_matrix * filter ;
00260 cpl_image * collapsed ;
00261 cpl_mask * bin_im ;
00262 double threshold, fillval, median_val, sigma ;
00263 int min_arclen = 0 ;
00264 cpl_apertures * det ;
00265 int nobj ;
00266 int ngoodpix ;
00267 int ny ;
00268
00269 ny = cpl_image_get_size_y(im) ;
00270
00271
00272 *label_im = NULL ;
00273
00274
00275 median_val = cpl_image_get_median_dev(im, &sigma) ;
00276 fillval = median_val-sigma/2.0 ;
00277 if (sinfoni_distortion_fill_badzones(im, xmin, ymin, xmax, ymax,
00278 fillval) == -1) {
00279 cpl_msg_error(fctid, "cannot fill bad zones") ;
00280 return NULL ;
00281 }
00282
00283
00284 filter = cpl_matrix_new(3, 3) ;
00285 cpl_matrix_fill(filter, 1.0) ;
00286 filt_im = cpl_image_filter_median(im, filter) ;
00287 cpl_matrix_delete(filter) ;
00288
00289
00290 median_val = cpl_image_get_median_dev(filt_im, &sigma) ;
00291
00292
00293 if (median_val < TRESH_MEDIAN_MIN) median_val = TRESH_MEDIAN_MIN ;
00294 if (sigma > TRESH_SIGMA_MAX) sigma = TRESH_SIGMA_MAX ;
00295
00296
00297 threshold = median_val + sigma * ARC_THRESHFACT ;
00298
00299
00300 collapsed = cpl_image_collapse_median_create(filt_im, 0, 0, 0) ;
00301
00302
00303 if (sinfoni_distortion_threshold1d(filt_im, median_val, collapsed, 0.0)==-1) {
00304 cpl_image_delete(filt_im) ;
00305 cpl_image_delete(collapsed) ;
00306 return NULL ;
00307 }
00308 cpl_image_delete(collapsed) ;
00309
00310
00311 bin_im = cpl_mask_threshold_image_create(filt_im, threshold,
00312 CPL_PIXEL_MAXVAL);
00313 cpl_image_delete(filt_im) ;
00314 if (bin_im == NULL) return NULL ;
00315
00316
00317 ngoodpix = cpl_mask_count(bin_im) ;
00318 if (ngoodpix < ARC_MINGOODPIX) {
00319 cpl_msg_error(fctid, "Too few (%d) white pixels", ngoodpix) ;
00320 cpl_mask_delete(bin_im) ;
00321 return NULL ;
00322 }
00323
00324
00325 filter = cpl_matrix_new(3, 3) ;
00326 cpl_matrix_fill(filter, 1.0) ;
00327 cpl_mask_closing(bin_im, filter) ;
00328 cpl_matrix_delete(filter) ;
00329
00330
00331 *label_im = cpl_image_labelise_mask_create(bin_im, &nobj) ;
00332 cpl_mask_delete(bin_im) ;
00333
00334
00335 if ((det = cpl_apertures_new_from_image(im, *label_im)) == NULL) {
00336 cpl_msg_error(fctid, "Cannot compute arcs stats") ;
00337 cpl_image_delete(*label_im) ;
00338 *label_im = NULL ;
00339 return NULL ;
00340 }
00341
00342
00343 min_arclen = (int)(ny / ARC_MINARCLENFACT) ;
00344
00345
00346 if (sinfoni_distortion_purge_arcs(im, &det, label_im, min_arclen,
00347 ARC_MAXARCWIDTH, arc_sat) == -1) {
00348 cpl_image_delete(*label_im) ;
00349 *label_im = NULL ;
00350 cpl_apertures_delete(det) ;
00351 return NULL ;
00352 }
00353 if (cpl_apertures_get_size(det) < ARC_MINNBARCS) {
00354 cpl_image_delete(*label_im) ;
00355 *label_im = NULL ;
00356 cpl_apertures_delete(det) ;
00357 return NULL ;
00358 }
00359
00360
00361 return det ;
00362 }
00363
00366 static int sinfoni_distortion_fill_badzones(
00367 cpl_image * im,
00368 int xmin,
00369 int ymin,
00370 int xmax,
00371 int ymax,
00372 double fillval)
00373 {
00374 float * pfi ;
00375 int nx, ny ;
00376 int i, j ;
00377
00378
00379 if (im == NULL) return -1 ;
00380 if (cpl_image_get_type(im) != CPL_TYPE_FLOAT) return -1 ;
00381
00382
00383 pfi = cpl_image_get_data_float(im) ;
00384 nx = cpl_image_get_size_x(im) ;
00385 ny = cpl_image_get_size_y(im) ;
00386
00387
00388 for (i=0 ; i<nx ; i++) {
00389 for (j=0 ; j<ny ; j++) {
00390 if ((i<xmin-1) || (i>xmax-1) || (j<ymin-1) || (j>ymax-1)) {
00391 pfi[i+j*nx] = (float)fillval ;
00392 }
00393 }
00394 }
00395 return 0 ;
00396 }
00397
00398 static int sinfoni_distortion_threshold1d(
00399 cpl_image * im,
00400 double threshold,
00401 cpl_image * im1d,
00402 double newval)
00403 {
00404 float * pim ;
00405 float * pim1d ;
00406 int nx, ny ;
00407 int i, j ;
00408
00409
00410 if (im == NULL) return -1 ;
00411 if (im1d == NULL) return -1 ;
00412 if (cpl_image_get_type(im) != CPL_TYPE_FLOAT) return -1 ;
00413 if (cpl_image_get_type(im1d) != CPL_TYPE_FLOAT) return -1 ;
00414
00415
00416 pim = cpl_image_get_data_float(im) ;
00417 pim1d = cpl_image_get_data_float(im1d) ;
00418 nx = cpl_image_get_size_x(im) ;
00419 ny = cpl_image_get_size_y(im) ;
00420
00421
00422 for (i=0 ; i<nx ; i++)
00423 if (pim1d[i] < threshold) {
00424 for (j=0 ; j<ny ; j++) pim[i+j*nx] = (float)newval ;
00425 }
00426
00427
00428 return 0 ;
00429 }
00430
00431 static int sinfoni_distortion_purge_arcs(
00432 cpl_image * im,
00433 cpl_apertures ** arcs,
00434 cpl_image ** lab_im,
00435 int min_arclen,
00436 int max_arcwidth,
00437 double arc_sat)
00438 {
00439 const char * fctid = "sinfoni_distortion_purge_arcs" ;
00440 int nb_arcs ;
00441 int * selection ;
00442 int arclen, arcwidth, edge ;
00443 double mean ;
00444 int * plabim ;
00445 cpl_mask * bin_im ;
00446 int nx, ny ;
00447 int i, j ;
00448
00449
00450 if (arcs == NULL) return -1 ;
00451 if (*arcs == NULL) return -1 ;
00452 if (*lab_im == NULL) return -1 ;
00453
00454
00455 nb_arcs = cpl_apertures_get_size(*arcs) ;
00456 nx = cpl_image_get_size_x(*lab_im) ;
00457 ny = cpl_image_get_size_y(*lab_im) ;
00458
00459
00460 selection = cpl_malloc(nb_arcs * sizeof(int)) ;
00461
00462
00463 for (i=0 ; i<nb_arcs ; i++) {
00464 arclen = cpl_apertures_get_top(*arcs, i+1) -
00465 cpl_apertures_get_bottom(*arcs, i+1) + 1 ;
00466 arcwidth = cpl_apertures_get_right(*arcs, i+1) -
00467 cpl_apertures_get_left(*arcs, i+1) + 1 ;
00468 edge = cpl_apertures_get_left_y(*arcs, i+1) ;
00469 mean = cpl_apertures_get_mean(*arcs, i+1) ;
00470
00471
00472 if ((arclen>min_arclen) && (arcwidth<max_arcwidth)
00473 && (edge>0) && (mean < arc_sat)) {
00474 selection[i] = 1 ;
00475 } else {
00476 selection[i] = 0 ;
00477 }
00478 }
00479
00480
00481 for (i=0 ; i<nb_arcs ; i++) {
00482 if (selection[i] == 0) {
00483 plabim = cpl_image_get_data_int(*lab_im) ;
00484 for (j=0 ; j<nx*ny ; j++) {
00485 if (plabim[j] == i+1) plabim[j] = 0 ;
00486 }
00487 }
00488 }
00489 cpl_free(selection) ;
00490
00491
00492 bin_im = cpl_mask_threshold_image_create(*lab_im, 0.5, CPL_PIXEL_MAXVAL) ;
00493 cpl_image_delete(*lab_im) ;
00494 *lab_im = cpl_image_labelise_mask_create(bin_im, NULL) ;
00495 cpl_mask_delete(bin_im) ;
00496
00497
00498 cpl_apertures_delete(*arcs) ;
00499 *arcs = cpl_apertures_new_from_image(im, *lab_im) ;
00500
00501
00502 if (cpl_apertures_get_size(*arcs) <= 0) {
00503 cpl_msg_error(fctid, "No valid arc found") ;
00504 return -1 ;
00505 }
00506
00507 return 0 ;
00508 }
00509
00510 static cpl_bivector ** sinfoni_distortion_get_arc_positions(
00511 cpl_image * in,
00512 cpl_image * label_im,
00513 cpl_apertures * det,
00514 int nb_samples,
00515 double ** lines_pos)
00516 {
00517 const char * fctid = "sinfoni_distortion_get_arc_positions" ;
00518 int n_arcs ;
00519 cpl_image * filt_img ;
00520 cpl_matrix * kernel ;
00521 cpl_bivector ** pos ;
00522 double * biv_x ;
00523 double * biv_y ;
00524 double x_finepos ;
00525 int * plabel_im ;
00526 int * arcs_samples_y ;
00527 int * computed ;
00528 double arclen ;
00529 int use_this_arc ;
00530 int obj ;
00531 int nx, ny ;
00532 int i, j, k ;
00533
00534
00535
00536
00537 n_arcs = cpl_apertures_get_size(det) ;
00538 nx = cpl_image_get_size_x(label_im) ;
00539 ny = cpl_image_get_size_y(label_im) ;
00540
00541
00542 pos = cpl_calloc(n_arcs, sizeof(cpl_bivector*)) ;
00543 for (i=0 ; i<n_arcs ; i++) pos[i] = cpl_bivector_new(nb_samples) ;
00544
00545
00546 kernel = cpl_matrix_new(3, 3) ;
00547 cpl_matrix_fill(kernel, 1.0) ;
00548 filt_img = cpl_image_filter_median(in, kernel) ;
00549 cpl_matrix_delete(kernel) ;
00550
00551
00552 arcs_samples_y = cpl_malloc(n_arcs * nb_samples * sizeof(int)) ;
00553 computed = cpl_calloc(n_arcs*nb_samples, sizeof(int)) ;
00554
00555
00556 for (j=0 ; j<n_arcs ; j++) {
00557 arclen = cpl_apertures_get_top(det,j+1) -
00558 cpl_apertures_get_bottom(det,j+1) + 1 ;
00559 for (i=0 ; i<nb_samples ; i++) {
00560 arcs_samples_y[i+j*nb_samples] =
00561 (int)(cpl_apertures_get_bottom(det, j+1) +
00562 (arclen * (i + 0.5)) / (double)nb_samples) ;
00563 }
00564 }
00565
00566
00567 plabel_im = cpl_image_get_data_int(label_im) ;
00568 for (i=0 ; i<nx ; i++) {
00569 for (j=0 ; j<ny ; j++) {
00570
00571
00572 obj = plabel_im[i + j * nx] ;
00573
00574 if (obj==0) continue ;
00575
00576 else obj-- ;
00577
00578 use_this_arc = 0 ;
00579 for (k=0 ; k<nb_samples ; k++) {
00580 if (arcs_samples_y[k+obj*nb_samples] == j) {
00581 use_this_arc = 1 ;
00582 break ;
00583 }
00584 }
00585 if ((use_this_arc) && (computed[k+obj*nb_samples] == 0)) {
00586
00587 if ((x_finepos = sinfoni_distortion_fine_pos(filt_img,
00588 label_im, i, j)) < 0.0) {
00589 cpl_msg_error(fctid, "cannot find fine arc position") ;
00590 cpl_image_delete(filt_img) ;
00591 cpl_free(arcs_samples_y);
00592 cpl_free(computed) ;
00593 for (i=0 ; i<n_arcs ; i++) cpl_bivector_delete(pos[i]);
00594 cpl_free(pos) ;
00595 return NULL ;
00596 } else {
00597 biv_x = cpl_bivector_get_x_data(pos[obj]) ;
00598 biv_y = cpl_bivector_get_y_data(pos[obj]) ;
00599 biv_x[k] = x_finepos ;
00600 biv_y[k] = j ;
00601 (*lines_pos)[obj] = cpl_apertures_get_centroid_x(det,obj+1);
00602 computed[k+obj*nb_samples] = 1 ;
00603 }
00604 }
00605 }
00606 }
00607
00608
00609 cpl_image_delete(filt_img) ;
00610 cpl_free(arcs_samples_y) ;
00611 cpl_free(computed) ;
00612 return pos ;
00613 }
00614
00615 static double sinfoni_distortion_fine_pos(
00616 cpl_image * im,
00617 cpl_image * label_im,
00618 int x,
00619 int y)
00620 {
00621 float * pim ;
00622 int * plabel_im ;
00623 int objnum ;
00624 int curr_obj ;
00625 int start_pos ;
00626 double grav_c ;
00627 double sum ;
00628 double max ;
00629 double val ;
00630 int maxpos ;
00631 int im_extrem ;
00632 double arc_pos ;
00633 int nx ;
00634
00635
00636 nx = cpl_image_get_size_x(im) ;
00637 grav_c = 0.0 ;
00638 sum = 0.0 ;
00639 start_pos = x ;
00640 maxpos = start_pos ;
00641 pim = cpl_image_get_data_float(im) ;
00642 max = (double)pim[start_pos + y * nx] ;
00643 plabel_im = cpl_image_get_data_int(label_im) ;
00644 objnum = plabel_im[start_pos + y * nx] ;
00645 im_extrem = nx ;
00646
00647
00648 do {
00649 val = (double)pim[start_pos + y * nx] ;
00650 if (start_pos == 0) grav_c = 0.0 ;
00651 else grav_c += start_pos * val ;
00652 sum += val ;
00653 if (val > max) {
00654 max = val ;
00655 maxpos = start_pos ;
00656 }
00657
00658
00659 start_pos++ ;
00660
00661 curr_obj = plabel_im[start_pos + y * nx] ;
00662 } while (curr_obj == objnum) ;
00663
00664
00665 if ((fabs(grav_c) < 1.0e-40) || (fabs(sum) < 1.0e-40)) {
00666 arc_pos = maxpos ;
00667 } else {
00668 arc_pos = grav_c / sum ;
00669 if (fabs(arc_pos) >= start_pos) arc_pos = maxpos ;
00670 }
00671
00672
00673 return arc_pos ;
00674 }
00675
00676
00682
00683 #define IS_NB_TESTPOINTS 8
00684 #define IS_MIN_SLOPE 0.01
00685 #define IS_MAX_SLOPE_DIF 0.075
00686 #define IS_MAX_FIT_EDGE_DIF 0.05
00687 #define IS_MIN_RAMP 10.0
00688 #define IS_MAX_MNERR 13.0
00689 #define IS_MAX_MNERR_DIF 8.0
00690 #define IS_MAX_INTER_DIF 20.0
00691 #define IS_SKIPZONE 2.5
00692 #define SQR(x) ((x)*(x))
00693 static cpl_image * sinfoni_distortion_remove_ramp(const cpl_image * in)
00694 {
00695 const char * fctid = "sinfoni_distortion_remove_ramp" ;
00696 int ramp_present ;
00697 int nx, ny ;
00698 int y, yhi, ylo;
00699 cpl_vector * tmp_vector ;
00700 cpl_bivector * testpointlo ;
00701 double * testpointlo_x ;
00702 double * testpointlo_y ;
00703 cpl_bivector * testpointhi ;
00704 double * testpointhi_x ;
00705 double * testpointhi_y ;
00706 int spacing;
00707 double rampdif, fitslope;
00708 double * pol_coefhi,
00709 * pol_coeflo ;
00710 cpl_vector * median ;
00711 double * median_data ;
00712 double medianerrlo, medianerrhi;
00713 double slope ;
00714 cpl_image * out ;
00715 float * pout ;
00716 float val ;
00717 int i, j ;
00718
00719
00720 nx = cpl_image_get_size_x(in) ;
00721 ny = cpl_image_get_size_y(in) ;
00722
00723
00724 if (in==NULL) return NULL ;
00725
00726 if (ny<IS_SKIPZONE*IS_NB_TESTPOINTS){
00727 cpl_msg_error(fctid, "image has %d lines, min=%d ",
00728 ny, (int)(IS_SKIPZONE*IS_NB_TESTPOINTS*2));
00729 return NULL ;
00730 }
00731
00732 slope=0.0 ;
00733 spacing= ny / (IS_SKIPZONE*IS_NB_TESTPOINTS) ;
00734 yhi = (int)(ny/2) ;
00735 ylo = yhi - 1 ;
00736
00737 testpointhi = cpl_bivector_new(IS_NB_TESTPOINTS) ;
00738 testpointhi_x = cpl_bivector_get_x_data(testpointhi) ;
00739 testpointhi_y = cpl_bivector_get_y_data(testpointhi) ;
00740 testpointlo = cpl_bivector_new(IS_NB_TESTPOINTS) ;
00741 testpointlo_x = cpl_bivector_get_x_data(testpointlo) ;
00742 testpointlo_y = cpl_bivector_get_y_data(testpointlo) ;
00743 for (i=0 ; i<IS_NB_TESTPOINTS ; i++) {
00744 y = yhi + i * spacing;
00745 tmp_vector = cpl_vector_new_from_image_row(in, y+1) ;
00746 testpointhi_x[i] = y - ny / 2;
00747 testpointhi_y[i] = cpl_vector_get_median(tmp_vector) ;
00748 cpl_vector_delete(tmp_vector) ;
00749 y = ylo - i * spacing;
00750 tmp_vector = cpl_vector_new_from_image_row(in, y+1) ;
00751 testpointlo_x[IS_NB_TESTPOINTS-i-1] = y ;
00752 testpointlo_y[IS_NB_TESTPOINTS-i-1]=cpl_vector_get_median(tmp_vector) ;
00753 cpl_vector_delete(tmp_vector) ;
00754 }
00755
00756
00757 pol_coefhi = irplib_flat_fit_slope_robust(testpointhi_x,
00758 testpointhi_y, IS_NB_TESTPOINTS) ;
00759 pol_coeflo = irplib_flat_fit_slope_robust(testpointlo_x,
00760 testpointlo_y, IS_NB_TESTPOINTS) ;
00761
00762
00763 median = cpl_vector_new(IS_NB_TESTPOINTS) ;
00764 median_data = cpl_vector_get_data(median) ;
00765 for (i=0 ; i<IS_NB_TESTPOINTS ; i++) {
00766 median_data[i]=SQR(testpointhi_y[i]
00767 - pol_coefhi[0] - pol_coefhi[1] * testpointhi_x[i]);
00768 }
00769 medianerrhi = cpl_vector_get_median(median) ;
00770 for (i=0; i<IS_NB_TESTPOINTS; i++) {
00771 median_data[i]=SQR(testpointlo_y[i]
00772 - pol_coeflo[0] - pol_coeflo[1] * testpointlo_x[i]);
00773 }
00774 medianerrlo = cpl_vector_get_median(median) ;
00775 cpl_vector_delete(median) ;
00776 rampdif = testpointlo_y[IS_NB_TESTPOINTS-1] - testpointhi_y[0];
00777 slope = rampdif / (ny/2.0) ;
00778 fitslope = (pol_coefhi[1] + pol_coeflo[1]) / 2.0 ;
00779
00780 cpl_bivector_delete(testpointlo);
00781 cpl_bivector_delete(testpointhi);
00782
00783
00784 if (fabs(rampdif)<IS_MIN_RAMP ||
00785 fabs(pol_coefhi[1]) < IS_MIN_SLOPE ||
00786 fabs(pol_coeflo[1]) < IS_MIN_SLOPE ||
00787 pol_coefhi[1]/pol_coeflo[1]<0.5 ||
00788 pol_coefhi[1]/pol_coeflo[1]>2.0 ||
00789 fabs(pol_coefhi[1]-pol_coeflo[1])>IS_MAX_SLOPE_DIF ||
00790 fabs(pol_coefhi[0]-pol_coeflo[0]) > IS_MAX_INTER_DIF ||
00791 medianerrlo> IS_MAX_MNERR ||
00792 medianerrhi> IS_MAX_MNERR ||
00793 fabs(medianerrlo-medianerrhi) >IS_MAX_MNERR_DIF ||
00794 fabs(slope-fitslope) > IS_MAX_FIT_EDGE_DIF ||
00795 slope/fitslope<0.5 ||
00796 slope/fitslope>2.0) ramp_present = 0 ;
00797 else ramp_present = 1 ;
00798
00799 cpl_free(pol_coeflo) ;
00800 cpl_free(pol_coefhi) ;
00801
00802
00803 out = cpl_image_duplicate(in) ;
00804 pout = cpl_image_get_data_float(out) ;
00805 if (ramp_present == 1) {
00806 for (j=0 ; j<ny/2 ; j++) {
00807 val = slope * (j-ny/2) ;
00808 for (i=0 ; i<nx ; i++)
00809 pout[i+j*nx] -= val ;
00810 }
00811 for (j=ny/2 ; j<ny ; j++) {
00812 val = slope * (j-ny) ;
00813 for (i=0 ; i<nx ; i++)
00814 pout[i+j*nx] -= val ;
00815 }
00816
00817 }
00818
00819 return out ;
00820 }
00821