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 #include "sinfo_distortion.h"
00041 #include "sinfo_functions.h"
00042 #include "sinfo_msg.h"
00043 #include "sinfo_error.h"
00044 #include "sinfo_utils_wrappers.h"
00045
00046
00047
00048
00049
00050 #define ARC_NBSAMPLES 20
00051 #define ARC_THRESHFACT (1.0/3.0)
00052 #define ARC_MINGOODPIX 100
00053 #define ARC_MINARCLENFACT 1.19
00054 #define ARC_MINNBARCS 32
00055 #define ARC_RANGE_FACT 3.0
00056 #define ARC_WINDOWSIZE 10
00057
00058 #define TRESH_MEDIAN_MIN 0.0
00059 #define TRESH_SIGMA_MAX 200.0
00060
00061
00065
00066
00067
00068
00069
00070 static cpl_apertures *
00071 sinfo_distortion_detect_arcs_new(cpl_image* ,cpl_image **,
00072 int,int,double,int,int,int,int,double,int);
00073
00074 static
00075 cpl_apertures * sinfo_distortion_detect_arcs(cpl_image *,
00076 cpl_image **, int, int, int, int, int, int) ;
00077 static int
00078 sinfo_distortion_fill_badzones(cpl_image *, int, int, int, int, double) ;
00079 static int
00080 sinfo_distortion_threshold1d(cpl_image *, double, cpl_image *, double) ;
00081 static int
00082 sinfo_distortion_purge_arcs(cpl_image *, cpl_apertures **,
00083 cpl_image **, int, int, double) ;
00084 static cpl_bivector **
00085 sinfo_distortion_get_arc_positions(cpl_image *,
00086 cpl_image *,
00087 cpl_apertures *, int, double **) ;
00088 static double sinfo_distortion_fine_pos(cpl_image *, cpl_image *, int, int) ;
00089 static int sinfo_distortion_sub_hor_lowpass(cpl_image *, int) ;
00090 static cpl_image * sinfo_distortion_remove_ramp(const cpl_image *) ;
00091 static cpl_image *
00092 sinfo_distortion_smooth(cpl_image* inp,const int r,const int d);
00093
00094
00095
00096
00097
00098
00099
00100
00103
00113 static cpl_image *
00114 sinfo_distortion_smooth(cpl_image* inp,const int r,const int d)
00115 {
00116
00117 int sx=0;
00118 int sy=0;
00119 int i=0;
00120 int j=0;
00121 int z=0;
00122
00123 float sum;
00124 cpl_image* out=NULL;
00125 float* pi=NULL;
00126 float* po=NULL;
00127 int min=0;
00128
00129 cknull(inp,"Null input image!");
00130 check_nomsg(sx=cpl_image_get_size_x(inp));
00131 check_nomsg(sy=cpl_image_get_size_y(inp));
00132 check_nomsg(out=cpl_image_duplicate(inp));
00133 check_nomsg(pi=cpl_image_get_data_float(inp));
00134 check_nomsg(po=cpl_image_get_data_float(out));
00135 min = r/2;
00136 switch (d) {
00137 case 0:
00138 for(j=0;j<sy;j++) {
00139 for(i=min;i<sx-min;i++) {
00140 sum=0;
00141 for(z=i-min;z<i+min+1;z++) {
00142 sum+=pi[z+j*sx];
00143 }
00144 po[i+j*sx]=sum/r;
00145 }
00146 }
00147 break;
00148
00149 case 1:
00150 for(i=0;i<sx;i++) {
00151 for(j=min;j<sy-min;j++) {
00152 sum=0;
00153 for(z=j-min;z<j+min+1;z++) {
00154 sum+=pi[i+z*sx];
00155 }
00156 po[i+j*sx]=sum;
00157 }
00158 }
00159 break;
00160
00161 default:
00162 sinfo_msg_error("case not supported");
00163 goto cleanup;
00164
00165 }
00166 check_nomsg(cpl_image_delete(inp));
00167 return out;
00168 cleanup:
00169 return NULL;
00170 }
00171
00172
00185 cpl_image *
00186 sinfo_distortion_image_restore(const cpl_image* inp,
00187 const int r,
00188 const int d,
00189 const double kappa,
00190 const int ks_method,
00191 const int n)
00192 {
00193
00194 int sx=0;
00195 int sy=0;
00196 int i=0;
00197 int j=0;
00198 int z=0;
00199 int k=0;
00200
00201
00202 cpl_image* out=NULL;
00203 float* pi=NULL;
00204 float* po=NULL;
00205 int min=0;
00206 cpl_vector* vec=NULL;
00207 double* pv=NULL;
00208 double mean=0;
00209 double median=0;
00210
00211 cknull(inp,"Null input image!");
00212 check_nomsg(sx=cpl_image_get_size_x(inp));
00213 check_nomsg(sy=cpl_image_get_size_y(inp));
00214 check_nomsg(out=cpl_image_duplicate(inp));
00215 check_nomsg(pi=cpl_image_get_data_float(inp));
00216 check_nomsg(po=cpl_image_get_data_float(out));
00217 min = r/2;
00218 check_nomsg(vec=cpl_vector_new(r));
00219 check_nomsg(pv=cpl_vector_get_data(vec));
00220 switch (d) {
00221 case 0:
00222 for(j=0;j<sy;j++) {
00223 for(i=min;i<sx-min;i++) {
00224 k=0;
00225 for(z=i-min;z<i+min+1;z++) {
00226 pv[k]=(double)pi[z+j*sx];
00227 k++;
00228 }
00229 cknull_nomsg(vec=sinfo_vector_clip(vec,kappa,n,ks_method));
00230 check_nomsg(mean=cpl_vector_get_mean(vec));
00231 check_nomsg(median=cpl_vector_get_mean(vec));
00232 po[i+j*sx]+=(mean-median);
00233 }
00234 }
00235 break;
00236
00237 case 1:
00238 for(i=0;i<sx;i++) {
00239 for(j=min;j<sy-min;j++) {
00240 k=0;
00241 for(z=j-min;z<j+min+1;z++) {
00242 pv[k]=(double)pi[i+z*sx];
00243 k++;
00244 }
00245 cknull_nomsg(vec=sinfo_vector_clip(vec,kappa,n,ks_method));
00246 check_nomsg(mean=cpl_vector_get_mean(vec));
00247 check_nomsg(median=cpl_vector_get_mean(vec));
00248 po[i+j*sx]+=(mean-median);
00249 }
00250 }
00251 break;
00252
00253 default:
00254 sinfo_msg_error("case not supported");
00255 goto cleanup;
00256
00257 }
00258 check_nomsg(cpl_image_delete((cpl_image*)inp));
00259 return out;
00260 cleanup:
00261 return NULL;
00262 }
00263
00264
00288
00289 cpl_polynomial * sinfo_distortion_estimate_new(
00290 const cpl_image * org,
00291 int xmin,
00292 int ymin,
00293 int xmax,
00294 int ymax,
00295 int auto_ramp_sub,
00296 int arc_sat,
00297 int max_arc_width,
00298 double kappa,
00299 double arcs_min_arclen_factor,
00300 int arcs_window_size,
00301 int smooth_rad,
00302 int degree,
00303 double offset,
00304 cpl_apertures ** arcs)
00305 {
00306 cpl_image * local_im ;
00307 cpl_image * label_image ;
00308 double rightmost, leftmost ;
00309 cpl_bivector ** arcs_pos ;
00310 double * parc_posx ;
00311 double * parc_posy ;
00312 double * lines_pos ;
00313 cpl_bivector * grid ;
00314 double * pgridx ;
00315 double * pgridy ;
00316 cpl_vector * values_to_fit ;
00317 double * pvalues_to_fit ;
00318 int min_arc_range ;
00319 int n_calib ;
00320 int n_arcs ;
00321 cpl_polynomial * poly2d ;
00322 int nx ;
00323 int i, j ;
00324
00325
00326 cpl_vector * lines_pos_tmp ;
00327 cpl_bivector * grid_tmp ;
00328 cpl_vector* grid_tot=0;
00329 double* pgrid_tmp_x=NULL;
00330 double* pgrid_tmp_y=NULL;
00331 double* pgrid_tot=NULL;
00332 double* plines_pos_tmp=NULL;
00333 int n_lines=0;
00334 int k=0;
00335
00336
00337
00338 if (org == NULL) return NULL ;
00339 if (kappa < 0.0) return NULL ;
00340
00341
00342 n_calib = ARC_NBSAMPLES ;
00343 nx = cpl_image_get_size_x(org) ;
00344
00345 if (auto_ramp_sub) {
00346 local_im = sinfo_distortion_remove_ramp(org) ;
00347 } else {
00348
00349 local_im = cpl_image_duplicate(org) ;
00350 }
00351 if (local_im == NULL) {
00352 cpl_msg_error(cpl_func, "Cannot clean the image") ;
00353 return NULL ;
00354 }
00355 if(smooth_rad > 1) {
00356 local_im=sinfo_distortion_smooth(local_im,smooth_rad,1);
00357
00358
00359
00360
00361 }
00362
00363 cpl_msg_info(cpl_func, "Detect arcs") ;
00364 if ((*arcs = sinfo_distortion_detect_arcs_new(local_im,
00365 &label_image,
00366 arc_sat, max_arc_width, kappa,
00367 xmin, ymin, xmax, ymax,
00368 arcs_min_arclen_factor,arcs_window_size)) == NULL) {
00369 cpl_image_delete(local_im) ;
00370 cpl_msg_error(cpl_func, "Cannot detect the arcs") ;
00371 return NULL ;
00372 }
00373 n_arcs = cpl_apertures_get_size(*arcs) ;
00374 cpl_msg_info(cpl_func, "%d detected arcs", n_arcs) ;
00375
00376
00377 rightmost = leftmost = cpl_apertures_get_max_x(*arcs, 1) ;
00378 for (i=1 ; i<n_arcs ; i++) {
00379 if (cpl_apertures_get_max_x(*arcs, i+1) < leftmost)
00380 leftmost = cpl_apertures_get_max_x(*arcs, i+1) ;
00381 if (cpl_apertures_get_max_x(*arcs, i+1) > rightmost)
00382 rightmost = cpl_apertures_get_max_x(*arcs, i+1) ;
00383 }
00384 min_arc_range = (int)(nx / ARC_RANGE_FACT) ;
00385 if ((int)(rightmost-leftmost) < min_arc_range) {
00386 cpl_msg_error(cpl_func, "too narrow range (%g-%g)<%d",
00387 rightmost, leftmost, min_arc_range) ;
00388 cpl_apertures_delete(*arcs) ;
00389 cpl_image_delete(local_im) ;
00390 cpl_image_delete(label_image) ;
00391 return NULL ;
00392 }
00393
00394
00395 cpl_msg_info(cpl_func, "Create deformation grid") ;
00396 lines_pos = cpl_malloc(n_arcs * sizeof(double)) ;
00397 if ((arcs_pos = sinfo_distortion_get_arc_positions(local_im,
00398 label_image, *arcs, n_calib, &lines_pos))==NULL){
00399 cpl_msg_error(cpl_func, "cannot get arcs positions") ;
00400 cpl_apertures_delete(*arcs) ;
00401 cpl_image_delete(local_im) ;
00402 cpl_free(lines_pos) ;
00403 cpl_image_delete(label_image) ;
00404 return NULL ;
00405 }
00406 cpl_image_delete(label_image) ;
00407 cpl_image_delete(local_im) ;
00408
00409
00410 lines_pos_tmp=cpl_vector_new(n_arcs);
00411 plines_pos_tmp=cpl_vector_get_data(lines_pos_tmp);
00412
00413
00414 sinfo_msg("Fit the 2d polynomial") ;
00415 grid = cpl_bivector_new(n_arcs * n_calib) ;
00416 pgridx = cpl_bivector_get_x_data(grid) ;
00417 pgridy = cpl_bivector_get_y_data(grid) ;
00418 values_to_fit = cpl_vector_new(n_arcs * n_calib) ;
00419 pvalues_to_fit = cpl_vector_get_data(values_to_fit) ;
00420
00421 for (i=0 ; i<n_arcs ; i++) {
00422 parc_posx = cpl_bivector_get_x_data(arcs_pos[i]) ;
00423 parc_posy = cpl_bivector_get_y_data(arcs_pos[i]) ;
00424 for (j=0 ; j<n_calib ; j++) {
00425 plines_pos_tmp[i]=lines_pos[i] ;
00426 pgridx[j+i*n_calib] = lines_pos[i] ;
00427 pgridy[j+i*n_calib] = parc_posy[j] ;
00428 pvalues_to_fit[j+i*n_calib] = parc_posx[j] ;
00429 }
00430 }
00431
00432 n_lines= n_arcs/32.0;
00433 if(n_lines < 1) {
00434 n_lines=1;
00435 }
00436 cpl_vector_sort(lines_pos_tmp,1);
00437 plines_pos_tmp=cpl_vector_get_data(lines_pos_tmp);
00438 grid_tmp=cpl_bivector_duplicate(grid);
00439 grid_tot=cpl_vector_new(n_calib);
00440 pgrid_tmp_x = cpl_bivector_get_x_data(grid_tmp) ;
00441 pgrid_tmp_y = cpl_bivector_get_y_data(grid_tmp) ;
00442 pgrid_tot = cpl_vector_get_data(grid_tot);
00443 for(j=0;j<n_calib;j++) {
00444 pgrid_tot[j]=0;
00445 for(i=n_lines ;i<n_arcs;i=i+n_lines)
00446 {
00447 for(k=0;k<n_lines;k++) {
00448 pgrid_tot[j] += (plines_pos_tmp[i+k]-
00449 plines_pos_tmp[k]);
00450
00451
00452
00453
00454 }
00455 }
00456
00457
00458
00459 }
00460
00461 for(j=0;j<n_calib;j++) {
00462 for (i=0 ; i<n_arcs ; i++) {
00463 pgridx[j+i*n_calib]=pgridx[j+i*n_calib]*
00464 ((nx/32.0)*n_lines*(31*32/2))/pgrid_tot[j]-offset;
00465
00466
00467
00468
00469 pgrid_tmp_x[j+i*n_calib]=pgrid_tmp_x[j+i*n_calib]*
00470 ((nx/32.0)*n_lines*(31*32/2))/pgrid_tot[j]-
00471 offset;
00472
00473 }
00474 }
00475 cpl_vector_delete(lines_pos_tmp);
00476 cpl_bivector_delete(grid_tmp);
00477 cpl_vector_delete(grid_tot);
00478
00479
00480
00481 for (i=0 ; i<n_arcs ; i++) cpl_bivector_delete(arcs_pos[i]) ;
00482 cpl_free(arcs_pos) ;
00483 cpl_free(lines_pos) ;
00484
00485
00486 if ((poly2d = cpl_polynomial_fit_2d_create(grid, values_to_fit,
00487 degree, NULL))==NULL) {
00488 cpl_msg_error(cpl_func, "cannot apply the 2d fit") ;
00489 cpl_bivector_delete(grid) ;
00490 cpl_vector_delete(values_to_fit) ;
00491 cpl_apertures_delete(*arcs) ;
00492 return NULL ;
00493 }
00494
00495
00496 cpl_bivector_delete(grid) ;
00497 cpl_vector_delete(values_to_fit) ;
00498 return poly2d ;
00499 }
00500
00501
00502
00503
00519
00520 static cpl_apertures * sinfo_distortion_detect_arcs_new(
00521 cpl_image * im,
00522 cpl_image ** label_im,
00523 int arc_sat,
00524 int max_arc_width,
00525 double kappa,
00526 int xmin,
00527 int ymin,
00528 int xmax,
00529 int ymax,
00530 double arcs_min_arclen_factor,
00531 int arcs_window_size)
00532 {
00533 cpl_image * filt_im ;
00534 cpl_matrix * filter ;
00535 cpl_image * collapsed ;
00536 cpl_mask * bin_im ;
00537 double threshold, fillval, median_val, sigma ;
00538 int min_arclen = 0 ;
00539 cpl_apertures * det ;
00540 int nobj ;
00541 int ngoodpix ;
00542 int ny ;
00543
00544 ny = cpl_image_get_size_y(im) ;
00545
00546 *label_im = NULL ;
00547
00548
00549 median_val = cpl_image_get_median_dev(im, &sigma) ;
00550 fillval = median_val-sigma/2.0 ;
00551 if (sinfo_distortion_fill_badzones(im, xmin, ymin, xmax, ymax,
00552 fillval) == -1) {
00553 cpl_msg_error(cpl_func, "cannot fill bad zones") ;
00554 return NULL ;
00555 }
00556
00557 filter = cpl_matrix_new(3, 1) ;
00558 cpl_matrix_fill(filter, 1.0) ;
00559
00560 filt_im = cpl_image_duplicate(im) ;
00561 cpl_matrix_delete(filter) ;
00562
00563
00564
00565 if (sinfo_distortion_sub_hor_lowpass(filt_im, arcs_window_size) == -1) {
00566 cpl_image_delete(filt_im) ;
00567 return NULL ;
00568 }
00569
00570
00571
00572
00573 median_val = cpl_image_get_median_dev(filt_im, &sigma) ;
00574
00575
00576 if (median_val < TRESH_MEDIAN_MIN) median_val = TRESH_MEDIAN_MIN ;
00577 if (sigma > TRESH_SIGMA_MAX) sigma = TRESH_SIGMA_MAX ;
00578
00579
00580 threshold = median_val + sigma * kappa ;
00581
00582
00583 collapsed = cpl_image_collapse_median_create(filt_im, 0, 0, 0) ;
00584
00585
00586 if (sinfo_distortion_threshold1d(filt_im, median_val,
00587 collapsed, 0.0)==-1) {
00588 cpl_msg_error(cpl_func, "cannot threshold the filtered image") ;
00589 cpl_image_delete(filt_im) ;
00590 cpl_image_delete(collapsed) ;
00591 return NULL ;
00592 }
00593 cpl_image_delete(collapsed) ;
00594
00595
00596 bin_im = cpl_mask_threshold_image_create(filt_im, threshold,
00597 CPL_PIXEL_MAXVAL);
00598 cpl_image_delete(filt_im) ;
00599 if (bin_im == NULL) {
00600 cpl_msg_error(cpl_func, "cannot binarise the image") ;
00601 return NULL ;
00602 }
00603
00604
00605 ngoodpix = cpl_mask_count(bin_im) ;
00606 if (ngoodpix < ARC_MINGOODPIX) {
00607 cpl_msg_error(cpl_func, "Too few (%d) white pixels", ngoodpix) ;
00608 cpl_mask_delete(bin_im) ;
00609 return NULL ;
00610 }
00611
00612
00613 filter = cpl_matrix_new(3, 3) ;
00614 cpl_matrix_fill(filter, 1.0) ;
00615 cpl_mask_closing(bin_im, filter) ;
00616 cpl_matrix_delete(filter) ;
00617
00618
00619 *label_im = cpl_image_labelise_mask_create(bin_im, &nobj) ;
00620 cpl_mask_delete(bin_im) ;
00621
00622
00623
00624
00625 if ((det = cpl_apertures_new_from_image(im, *label_im)) == NULL) {
00626 cpl_msg_error(cpl_func, "Cannot compute arcs stats") ;
00627 cpl_image_delete(*label_im) ;
00628 *label_im = NULL ;
00629 return NULL ;
00630 }
00631
00632 min_arclen = (int)(ny /arcs_min_arclen_factor) ;
00633
00634
00635
00636
00637 if (sinfo_distortion_purge_arcs(im, &det, label_im, min_arclen,
00638 max_arc_width, arc_sat) == -1) {
00639 cpl_msg_error(cpl_func, "Cannot purge the arcs") ;
00640 cpl_image_delete(*label_im) ;
00641 *label_im = NULL ;
00642 cpl_apertures_delete(det) ;
00643 return NULL ;
00644 }
00645
00646 if (cpl_apertures_get_size(det) < ARC_MINNBARCS) {
00647 cpl_msg_error(cpl_func, "Not enough valid arcs (%d < %d)",
00648 cpl_apertures_get_size(det), ARC_MINNBARCS) ;
00649 cpl_image_delete(*label_im) ;
00650 *label_im = NULL ;
00651 cpl_apertures_delete(det) ;
00652 return NULL ;
00653 }
00654
00655
00656 return det ;
00657 }
00658
00659
00660
00661
00685
00686 cpl_polynomial * sinfo_distortion_estimate(
00687 const cpl_image * org,
00688 int xmin,
00689 int ymin,
00690 int xmax,
00691 int ymax,
00692 int auto_ramp_sub,
00693 int arc_sat,
00694 int max_arc_width,
00695 int degree,
00696 double offset,
00697 cpl_apertures ** arcs)
00698 {
00699 const char * fctid = "sinfo_distortion_estimate" ;
00700 cpl_image * local_im ;
00701 cpl_image * label_image ;
00702 double rightmost, leftmost ;
00703 cpl_bivector ** arcs_pos ;
00704 double * parc_posx ;
00705 double * parc_posy ;
00706 double * lines_pos ;
00707 cpl_bivector * grid ;
00708 double * pgridx ;
00709 double * pgridy ;
00710 cpl_vector * values_to_fit ;
00711 double * pvalues_to_fit ;
00712 int min_arc_range ;
00713 int n_calib ;
00714 int n_arcs ;
00715 cpl_polynomial * poly2d ;
00716 int nx ;
00717 int i, j ;
00718
00719
00720 cpl_vector * lines_pos_tmp ;
00721 cpl_bivector * grid_tmp ;
00722 int n_lines=0;
00723 int k=0;
00724 cpl_vector* grid_tot=0;
00725 double* pgrid_tmp_x=NULL;
00726 double* pgrid_tmp_y=NULL;
00727 double* pgrid_tot=NULL;
00728 double* plines_pos_tmp=NULL;
00729
00730
00731 if (org == NULL) return NULL ;
00732
00733
00734 n_calib = ARC_NBSAMPLES ;
00735 nx = cpl_image_get_size_x(org) ;
00736
00737 if (auto_ramp_sub) {
00738 local_im = sinfo_distortion_remove_ramp(org) ;
00739 } else {
00740
00741 local_im = cpl_image_duplicate(org) ;
00742 }
00743 if (local_im == NULL) {
00744 cpl_msg_error(fctid, "Cannot clean the image") ;
00745 return NULL ;
00746 }
00747
00748
00749 cpl_msg_info(fctid, "Detect arcs") ;
00750 if ((*arcs = sinfo_distortion_detect_arcs(local_im,
00751 &label_image,
00752 arc_sat, max_arc_width,
00753 xmin, ymin, xmax, ymax)) == NULL) {
00754 cpl_image_delete(local_im) ;
00755 cpl_msg_error(fctid, "Cannot detect the arcs") ;
00756 return NULL ;
00757 }
00758 n_arcs = cpl_apertures_get_size(*arcs) ;
00759 cpl_msg_info(fctid, "%d detected arcs", n_arcs) ;
00760
00761
00762 rightmost = leftmost = cpl_apertures_get_max_x(*arcs, 1) ;
00763 for (i=1 ; i<n_arcs ; i++) {
00764 if (cpl_apertures_get_max_x(*arcs, i+1) < leftmost)
00765 leftmost = cpl_apertures_get_max_x(*arcs, i+1) ;
00766 if (cpl_apertures_get_max_x(*arcs, i+1) > rightmost)
00767 rightmost = cpl_apertures_get_max_x(*arcs, i+1) ;
00768 }
00769 min_arc_range = (int)(nx / ARC_RANGE_FACT) ;
00770 if ((int)(rightmost-leftmost) < min_arc_range) {
00771 cpl_msg_error(fctid, "too narrow range (%g-%g)<%d",
00772 rightmost, leftmost, min_arc_range) ;
00773 cpl_apertures_delete(*arcs) ;
00774 cpl_image_delete(local_im) ;
00775 cpl_image_delete(label_image) ;
00776 return NULL ;
00777 }
00778
00779
00780 cpl_msg_info(fctid, "Create deformation grid") ;
00781 lines_pos = cpl_malloc(n_arcs * sizeof(double)) ;
00782 if ((arcs_pos = sinfo_distortion_get_arc_positions(local_im,
00783 label_image, *arcs, n_calib, &lines_pos))==NULL){
00784 cpl_msg_error(fctid, "cannot get arcs positions") ;
00785 cpl_apertures_delete(*arcs) ;
00786 cpl_image_delete(local_im) ;
00787 cpl_free(lines_pos) ;
00788 cpl_image_delete(label_image) ;
00789 return NULL ;
00790 }
00791 cpl_image_delete(label_image) ;
00792 cpl_image_delete(local_im) ;
00793
00794
00795 lines_pos_tmp=cpl_vector_new(n_arcs);
00796 plines_pos_tmp=cpl_vector_get_data(lines_pos_tmp);
00797
00798 cpl_msg_info(fctid, "Fit the 2d polynomial") ;
00799 grid = cpl_bivector_new(n_arcs * n_calib) ;
00800 pgridx = cpl_bivector_get_x_data(grid) ;
00801 pgridy = cpl_bivector_get_y_data(grid) ;
00802 values_to_fit = cpl_vector_new(n_arcs * n_calib) ;
00803 pvalues_to_fit = cpl_vector_get_data(values_to_fit) ;
00804 for (i=0 ; i<n_arcs ; i++) {
00805 parc_posx = cpl_bivector_get_x_data(arcs_pos[i]) ;
00806 parc_posy = cpl_bivector_get_y_data(arcs_pos[i]) ;
00807 for (j=0 ; j<n_calib ; j++) {
00808 plines_pos_tmp[i]=lines_pos[i] ;
00809 pgridx[j+i*n_calib] = lines_pos[i] ;
00810 pgridy[j+i*n_calib] = parc_posy[j] ;
00811 pvalues_to_fit[j+i*n_calib] = parc_posx[j];
00812
00813
00814
00815
00816
00817 }
00818 }
00819
00820
00821
00822 n_lines= n_arcs/32.0;
00823 if(n_lines < 1) {
00824 n_lines=1;
00825 }
00826 cpl_vector_sort(lines_pos_tmp,1);
00827 plines_pos_tmp=cpl_vector_get_data(lines_pos_tmp);
00828
00829 grid_tmp=cpl_bivector_duplicate(grid);
00830 grid_tot=cpl_vector_new(n_calib);
00831 pgrid_tmp_x = cpl_bivector_get_x_data(grid_tmp) ;
00832 pgrid_tmp_y = cpl_bivector_get_y_data(grid_tmp) ;
00833 pgrid_tot = cpl_vector_get_data(grid_tot);
00834 for(j=0;j<n_calib;j++) {
00835 pgrid_tot[j]=0;
00836 for(i=n_lines ;i<n_arcs;i=i+n_lines)
00837 {
00838 for(k=0;k<n_lines;k++) {
00839 pgrid_tot[j] += (plines_pos_tmp[i+k]-
00840 plines_pos_tmp[k]);
00841
00842
00843
00844
00845 }
00846 }
00847
00848
00849
00850 }
00851
00852 for(j=0;j<n_calib;j++) {
00853 for (i=0 ; i<n_arcs ; i++) {
00854 pgridx[j+i*n_calib]=pgridx[j+i*n_calib]*
00855 ((nx/32.0)*n_lines*(31*32/2))/pgrid_tot[j]-offset;
00856
00857
00858
00859
00860 pgrid_tmp_x[j+i*n_calib]=pgrid_tmp_x[j+i*n_calib]*
00861 ((nx/32.0)*n_lines*(31*32/2))/pgrid_tot[j]-
00862 offset;
00863
00864 }
00865 }
00866
00867
00868
00869 for (i=0 ; i<n_arcs ; i++) cpl_bivector_delete(arcs_pos[i]) ;
00870 cpl_free(arcs_pos) ;
00871 cpl_free(lines_pos) ;
00872
00873
00874 if ((poly2d = cpl_polynomial_fit_2d_create(grid, values_to_fit,
00875 degree, NULL))==NULL) {
00876 cpl_msg_error(fctid, "cannot apply the 2d fit") ;
00877 cpl_bivector_delete(grid) ;
00878 cpl_vector_delete(values_to_fit) ;
00879 cpl_apertures_delete(*arcs) ;
00880 return NULL ;
00881 }
00882
00883
00884 cpl_bivector_delete(grid) ;
00885 cpl_vector_delete(values_to_fit) ;
00886 return poly2d ;
00887 }
00888
00891
00906
00907 static cpl_apertures * sinfo_distortion_detect_arcs(
00908 cpl_image * im,
00909 cpl_image ** label_im,
00910 int arc_sat,
00911 int max_arc_width,
00912 int xmin,
00913 int ymin,
00914 int xmax,
00915 int ymax)
00916 {
00917 const char * fctid = "sinfo_distortion_detect_arcs" ;
00918 cpl_image * filt_im ;
00919 cpl_matrix * filter ;
00920 cpl_image * collapsed ;
00921 cpl_mask * bin_im ;
00922 double threshold, fillval, median_val, sigma ;
00923 int min_arclen = 0 ;
00924 cpl_apertures * det ;
00925 int nobj ;
00926 int ngoodpix ;
00927 int ny ;
00928
00929 ny = cpl_image_get_size_y(im) ;
00930
00931
00932 *label_im = NULL ;
00933
00934
00935 median_val = cpl_image_get_median_dev(im, &sigma) ;
00936 fillval = median_val-sigma/2.0 ;
00937 if (sinfo_distortion_fill_badzones(im, xmin, ymin, xmax, ymax,
00938 fillval) == -1) {
00939 cpl_msg_error(fctid, "cannot fill bad zones") ;
00940 return NULL ;
00941 }
00942
00943
00944 filter = cpl_matrix_new(3, 1) ;
00945 cpl_matrix_fill(filter, 1.0) ;
00946
00947 filt_im = cpl_image_duplicate(im) ;
00948 cpl_matrix_delete(filter) ;
00949
00950
00951 if (sinfo_distortion_sub_hor_lowpass(filt_im, ARC_WINDOWSIZE) == -1) {
00952 cpl_image_delete(filt_im) ;
00953 return NULL ;
00954 }
00955
00956
00957 median_val = cpl_image_get_median_dev(filt_im, &sigma) ;
00958
00959
00960 if (median_val < TRESH_MEDIAN_MIN) median_val = TRESH_MEDIAN_MIN ;
00961 if (sigma > TRESH_SIGMA_MAX) sigma = TRESH_SIGMA_MAX ;
00962
00963
00964 threshold = median_val + sigma * ARC_THRESHFACT ;
00965
00966
00967 collapsed = cpl_image_collapse_median_create(filt_im, 0, 0, 0) ;
00968
00969
00970 if (sinfo_distortion_threshold1d(filt_im, median_val,
00971 collapsed, 0.0)==-1) {
00972 cpl_msg_error(fctid, "cannot threshold the filtered image") ;
00973 cpl_image_delete(filt_im) ;
00974 cpl_image_delete(collapsed) ;
00975 return NULL ;
00976 }
00977 cpl_image_delete(collapsed) ;
00978
00979
00980 bin_im = cpl_mask_threshold_image_create(filt_im, threshold,
00981 CPL_PIXEL_MAXVAL);
00982 cpl_image_delete(filt_im) ;
00983 if (bin_im == NULL) {
00984 cpl_msg_error(fctid, "cannot binarise the image") ;
00985 return NULL ;
00986 }
00987
00988
00989 ngoodpix = cpl_mask_count(bin_im) ;
00990 if (ngoodpix < ARC_MINGOODPIX) {
00991 cpl_msg_error(fctid, "Too few (%d) white pixels", ngoodpix) ;
00992 cpl_mask_delete(bin_im) ;
00993 return NULL ;
00994 }
00995
00996
00997 filter = cpl_matrix_new(3, 3) ;
00998 cpl_matrix_fill(filter, 1.0) ;
00999 cpl_mask_closing(bin_im, filter) ;
01000 cpl_matrix_delete(filter) ;
01001
01002
01003 *label_im = cpl_image_labelise_mask_create(bin_im, &nobj) ;
01004 cpl_mask_delete(bin_im) ;
01005
01006
01007 if ((det = cpl_apertures_new_from_image(im, *label_im)) == NULL) {
01008 cpl_msg_error(fctid, "Cannot compute arcs stats") ;
01009 cpl_image_delete(*label_im) ;
01010 *label_im = NULL ;
01011 return NULL ;
01012 }
01013
01014
01015 min_arclen = (int)(ny / ARC_MINARCLENFACT) ;
01016
01017
01018 if (sinfo_distortion_purge_arcs(im, &det, label_im, min_arclen,
01019 max_arc_width, arc_sat) == -1) {
01020 cpl_msg_error(fctid, "Cannot purge the arcs") ;
01021 cpl_image_delete(*label_im) ;
01022 *label_im = NULL ;
01023 cpl_apertures_delete(det) ;
01024 return NULL ;
01025 }
01026 if (cpl_apertures_get_size(det) < ARC_MINNBARCS) {
01027 cpl_msg_error(fctid, "Not enough valid arcs (%d < %d)",
01028 cpl_apertures_get_size(det), ARC_MINNBARCS) ;
01029 cpl_image_delete(*label_im) ;
01030 *label_im = NULL ;
01031 cpl_apertures_delete(det) ;
01032 return NULL ;
01033 }
01034
01035
01036 return det ;
01037 }
01038
01039 static int sinfo_distortion_fill_badzones(
01040 cpl_image * im,
01041 int xmin,
01042 int ymin,
01043 int xmax,
01044 int ymax,
01045 double fillval)
01046 {
01047 float * pfi ;
01048 int nx, ny ;
01049 int i, j ;
01050
01051
01052 if (im == NULL) return -1 ;
01053 if (cpl_image_get_type(im) != CPL_TYPE_FLOAT) return -1 ;
01054
01055
01056 pfi = cpl_image_get_data_float(im) ;
01057 nx = cpl_image_get_size_x(im) ;
01058 ny = cpl_image_get_size_y(im) ;
01059
01060
01061 for (i=0 ; i<nx ; i++) {
01062 for (j=0 ; j<ny ; j++) {
01063 if ((i<xmin-1) || (i>xmax-1) || (j<ymin-1) || (j>ymax-1)) {
01064 pfi[i+j*nx] = (float)fillval ;
01065 }
01066 }
01067 }
01068 return 0 ;
01069 }
01070
01071 static int sinfo_distortion_threshold1d(
01072 cpl_image * im,
01073 double threshold,
01074 cpl_image * im1d,
01075 double newval)
01076 {
01077 float * pim ;
01078 float * pim1d ;
01079 int nx, ny ;
01080 int i, j ;
01081
01082
01083 if (im == NULL) return -1 ;
01084 if (im1d == NULL) return -1 ;
01085 if (cpl_image_get_type(im) != CPL_TYPE_FLOAT) return -1 ;
01086 if (cpl_image_get_type(im1d) != CPL_TYPE_FLOAT) return -1 ;
01087
01088
01089 pim = cpl_image_get_data_float(im) ;
01090 pim1d = cpl_image_get_data_float(im1d) ;
01091 nx = cpl_image_get_size_x(im) ;
01092 ny = cpl_image_get_size_y(im) ;
01093
01094
01095 for (i=0 ; i<nx ; i++)
01096 if (pim1d[i] < threshold) {
01097 for (j=0 ; j<ny ; j++) pim[i+j*nx] = (float)newval ;
01098 }
01099
01100
01101 return 0 ;
01102 }
01103
01104 static int sinfo_distortion_sub_hor_lowpass(
01105 cpl_image * im,
01106 int filt_size)
01107 {
01108 cpl_vector * linehi ;
01109 cpl_vector * linelo ;
01110 cpl_vector * avglinehi ;
01111 cpl_vector * avglinelo ;
01112 double * pavglinehi ;
01113 float * pim ;
01114 int lopos, hipos, nx, ny ;
01115 int i, j ;
01116
01117
01118 if (im == NULL) return -1 ;
01119 if (filt_size <= 0) return -1 ;
01120
01121
01122 nx = cpl_image_get_size_x(im) ;
01123 ny = cpl_image_get_size_y(im) ;
01124 lopos = (int)(ny/4) ;
01125 hipos = (int)(3*ny/4) ;
01126
01127
01128 if ((linehi = cpl_vector_new_from_image_row(im, hipos)) == NULL) {
01129 return -1 ;
01130 }
01131 if ((linelo = cpl_vector_new_from_image_row(im, lopos)) == NULL) {
01132 cpl_vector_delete(linehi) ;
01133 return -1 ;
01134 }
01135
01136
01137 if ((avglinehi = cpl_vector_filter_median_create(linehi,
01138 filt_size)) == NULL) {
01139 cpl_vector_delete(linehi) ;
01140 cpl_vector_delete(linelo) ;
01141 return -1 ;
01142 }
01143 cpl_vector_delete(linehi) ;
01144
01145 if ((avglinelo = cpl_vector_filter_median_create(linelo,
01146 filt_size)) == NULL) {
01147 cpl_vector_delete(linelo) ;
01148 cpl_vector_delete(avglinehi) ;
01149 return -1 ;
01150 }
01151 cpl_vector_delete(linelo) ;
01152
01153
01154 cpl_vector_add(avglinehi, avglinelo) ;
01155 cpl_vector_delete(avglinelo) ;
01156 cpl_vector_divide_scalar(avglinehi, 2.0) ;
01157
01158
01159 pavglinehi = cpl_vector_get_data(avglinehi) ;
01160 pim = cpl_image_get_data_float(im) ;
01161 for (i=0 ; i<nx ; i++) {
01162 for (j=0 ; j<ny ; j++) {
01163 pim[i+j*nx] -= pavglinehi[i] ;
01164 }
01165 }
01166 cpl_vector_delete(avglinehi) ;
01167
01168 return 0 ;
01169 }
01170
01171
01172
01173
01174
01175
01176
01177
01178 static int sinfo_distortion_purge_arcs(
01179 cpl_image * im,
01180 cpl_apertures ** arcs,
01181 cpl_image ** lab_im,
01182 int min_arclen,
01183 int max_arcwidth,
01184 double arc_sat)
01185 {
01186 const char * fctid = "sinfo_distortion_purge_arcs" ;
01187 int nb_arcs ;
01188 int * selection ;
01189 int arclen, arcwidth, edge ;
01190 double mean ;
01191 int * plabim ;
01192 cpl_mask * bin_im ;
01193 int nx, ny ;
01194 int i, j ;
01195
01196
01197 if (arcs == NULL) return -1 ;
01198 if (*arcs == NULL) return -1 ;
01199 if (*lab_im == NULL) return -1 ;
01200
01201
01202 nb_arcs = cpl_apertures_get_size(*arcs) ;
01203 nx = cpl_image_get_size_x(*lab_im) ;
01204 ny = cpl_image_get_size_y(*lab_im) ;
01205
01206
01207 selection = cpl_malloc(nb_arcs * sizeof(int)) ;
01208
01209
01210 for (i=0 ; i<nb_arcs ; i++) {
01211 arclen = cpl_apertures_get_top(*arcs, i+1) -
01212 cpl_apertures_get_bottom(*arcs, i+1) + 1 ;
01213 arcwidth = cpl_apertures_get_right(*arcs, i+1) -
01214 cpl_apertures_get_left(*arcs, i+1) + 1 ;
01215 edge = cpl_apertures_get_left_y(*arcs, i+1) ;
01216 mean = cpl_apertures_get_mean(*arcs, i+1) ;
01217
01218
01219
01220 if (
01221 (arclen>min_arclen) &&
01222 (arcwidth<max_arcwidth) &&
01223 (edge>0) &&
01224 (mean < arc_sat)) {
01225
01226
01227
01228
01229
01230 selection[i] = 1 ;
01231 } else {
01232
01233
01234
01235
01236
01237 selection[i] = 0 ;
01238 }
01239 }
01240
01241
01242 for (i=0 ; i<nb_arcs ; i++) {
01243 if (selection[i] == 0) {
01244 plabim = cpl_image_get_data_int(*lab_im) ;
01245 for (j=0 ; j<nx*ny ; j++) {
01246 if (plabim[j] == i+1) plabim[j] = 0 ;
01247 }
01248 }
01249 }
01250 cpl_free(selection) ;
01251
01252
01253 bin_im = cpl_mask_threshold_image_create(*lab_im, 0.5, CPL_PIXEL_MAXVAL) ;
01254 cpl_image_delete(*lab_im) ;
01255 *lab_im = cpl_image_labelise_mask_create(bin_im, NULL) ;
01256 cpl_mask_delete(bin_im) ;
01257
01258
01259 cpl_apertures_delete(*arcs) ;
01260 *arcs = cpl_apertures_new_from_image(im, *lab_im) ;
01261
01262
01263 if (cpl_apertures_get_size(*arcs) <= 0) {
01264 cpl_msg_error(fctid, "No valid arc found") ;
01265 return -1 ;
01266 }
01267
01268 return 0 ;
01269 }
01270
01271 static cpl_bivector **
01272 sinfo_distortion_get_arc_positions(
01273 cpl_image * in,
01274 cpl_image * label_im,
01275 cpl_apertures * det,
01276 int nb_samples,
01277 double ** lines_pos)
01278 {
01279 const char * fctid = "sinfo_distortion_get_arc_positions" ;
01280 int n_arcs ;
01281 cpl_image * filt_img ;
01282 cpl_matrix * kernel ;
01283 cpl_bivector ** pos ;
01284 double * biv_x ;
01285 double * biv_y ;
01286 double x_finepos ;
01287 int * plabel_im ;
01288 int * arcs_samples_y ;
01289 int * computed ;
01290 double arclen ;
01291 int use_this_arc ;
01292 int obj ;
01293 int nx, ny ;
01294 int i, j, k ;
01295
01296
01297
01298
01299 n_arcs = cpl_apertures_get_size(det) ;
01300 nx = cpl_image_get_size_x(label_im) ;
01301 ny = cpl_image_get_size_y(label_im) ;
01302
01303
01304 pos = cpl_calloc(n_arcs, sizeof(cpl_bivector*)) ;
01305 for (i=0 ; i<n_arcs ; i++) pos[i] = cpl_bivector_new(nb_samples) ;
01306
01307
01308 kernel = cpl_matrix_new(3, 3) ;
01309 cpl_matrix_fill(kernel, 1.0) ;
01310 filt_img = cpl_image_filter_median(in, kernel) ;
01311 cpl_matrix_delete(kernel) ;
01312
01313
01314 arcs_samples_y = cpl_malloc(n_arcs * nb_samples * sizeof(int)) ;
01315 computed = cpl_calloc(n_arcs*nb_samples, sizeof(int)) ;
01316
01317
01318 for (j=0 ; j<n_arcs ; j++) {
01319 arclen = cpl_apertures_get_top(det,j+1) -
01320 cpl_apertures_get_bottom(det,j+1) + 1 ;
01321 for (i=0 ; i<nb_samples ; i++) {
01322 arcs_samples_y[i+j*nb_samples] =
01323 (int)(cpl_apertures_get_bottom(det, j+1) +
01324 (arclen * (i + 0.5)) / (double)nb_samples) ;
01325 }
01326 }
01327
01328
01329 plabel_im = cpl_image_get_data_int(label_im) ;
01330 for (i=0 ; i<nx ; i++) {
01331 for (j=0 ; j<ny ; j++) {
01332
01333
01334 obj = plabel_im[i + j * nx] ;
01335
01336 if (obj==0) continue ;
01337
01338 else obj-- ;
01339
01340 use_this_arc = 0 ;
01341 for (k=0 ; k<nb_samples ; k++) {
01342 if (arcs_samples_y[k+obj*nb_samples] == j) {
01343 use_this_arc = 1 ;
01344 break ;
01345 }
01346 }
01347 if ((use_this_arc) && (computed[k+obj*nb_samples] == 0)) {
01348
01349 if ((x_finepos = sinfo_distortion_fine_pos(filt_img,
01350 label_im, i, j)) < 0.0) {
01351 cpl_msg_error(fctid, "cannot find fine arc position") ;
01352 cpl_image_delete(filt_img) ;
01353 cpl_free(arcs_samples_y);
01354 cpl_free(computed) ;
01355 for (i=0 ; i<n_arcs ; i++) cpl_bivector_delete(pos[i]);
01356 cpl_free(pos) ;
01357 return NULL ;
01358 } else {
01359 biv_x = cpl_bivector_get_x_data(pos[obj]) ;
01360 biv_y = cpl_bivector_get_y_data(pos[obj]) ;
01361 biv_x[k] = x_finepos ;
01362 biv_y[k] = j ;
01363 (*lines_pos)[obj] = cpl_apertures_get_centroid_x(det,obj+1);
01364 computed[k+obj*nb_samples] = 1 ;
01365 }
01366 }
01367 }
01368 }
01369
01370
01371 cpl_image_delete(filt_img) ;
01372 cpl_free(arcs_samples_y) ;
01373 cpl_free(computed) ;
01374 return pos ;
01375 }
01376
01377 static double
01378 sinfo_distortion_fine_pos(
01379 cpl_image * im,
01380 cpl_image * label_im,
01381 int x,
01382 int y)
01383 {
01384 float * pim ;
01385 int * plabel_im ;
01386 int objnum ;
01387 int curr_obj ;
01388 int start_pos ;
01389 double grav_c ;
01390 double sum ;
01391 double max ;
01392 double val ;
01393 int maxpos ;
01394 int im_extrem ;
01395 double arc_pos ;
01396 int nx ;
01397
01398
01399 nx = cpl_image_get_size_x(im) ;
01400 grav_c = 0.0 ;
01401 sum = 0.0 ;
01402 start_pos = x ;
01403 maxpos = start_pos ;
01404 pim = cpl_image_get_data_float(im) ;
01405 max = (double)pim[start_pos + y * nx] ;
01406 plabel_im = cpl_image_get_data_int(label_im) ;
01407 objnum = plabel_im[start_pos + y * nx] ;
01408 im_extrem = nx ;
01409
01410
01411 do {
01412 val = (double)pim[start_pos + y * nx] ;
01413 if (start_pos == 0) grav_c = 0.0 ;
01414 else grav_c += start_pos * val ;
01415 sum += val ;
01416 if (val > max) {
01417 max = val ;
01418 maxpos = start_pos ;
01419 }
01420
01421
01422 start_pos++ ;
01423
01424 curr_obj = plabel_im[start_pos + y * nx] ;
01425 } while (curr_obj == objnum) ;
01426
01427
01428 if ((fabs(grav_c) < 1.0e-40) || (fabs(sum) < 1.0e-40)) {
01429 arc_pos = maxpos ;
01430 } else {
01431 arc_pos = grav_c / sum ;
01432 if (fabs(arc_pos) >= start_pos) arc_pos = maxpos ;
01433 }
01434
01435
01436 return arc_pos ;
01437 }
01438
01439
01445
01446 #define IS_NB_TESTPOINTS 8
01447 #define IS_MIN_SLOPE 0.01
01448 #define IS_MAX_SLOPE_DIF 0.075
01449 #define IS_MAX_FIT_EDGE_DIF 0.05
01450 #define IS_MIN_RAMP 10.0
01451 #define IS_MAX_MNERR 13.0
01452 #define IS_MAX_MNERR_DIF 8.0
01453 #define IS_MAX_INTER_DIF 20.0
01454 #define IS_SKIPZONE 2.5
01455 #define SQR(x) ((x)*(x))
01456 static cpl_image * sinfo_distortion_remove_ramp(const cpl_image * in)
01457 {
01458 const char * fctid = "sinfo_distortion_remove_ramp" ;
01459 int ramp_present ;
01460 int nx, ny ;
01461 int y, yhi, ylo;
01462 cpl_vector * tmp_vector ;
01463 cpl_bivector * testpointlo ;
01464 double * testpointlo_x ;
01465 double * testpointlo_y ;
01466 cpl_bivector * testpointhi ;
01467 double * testpointhi_x ;
01468 double * testpointhi_y ;
01469 int spacing;
01470 double rampdif, fitslope;
01471 double * pol_coefhi,
01472 * pol_coeflo ;
01473 cpl_vector * median ;
01474 double * median_data ;
01475 double medianerrlo, medianerrhi;
01476 double slope ;
01477 cpl_image * out ;
01478 float * pout ;
01479 float val ;
01480 int i, j ;
01481
01482
01483 nx = cpl_image_get_size_x(in) ;
01484 ny = cpl_image_get_size_y(in) ;
01485
01486
01487 if (in==NULL) return NULL ;
01488
01489 if (ny<IS_SKIPZONE*IS_NB_TESTPOINTS){
01490 cpl_msg_error(fctid, "image has %d lines, min=%d ",
01491 ny, (int)(IS_SKIPZONE*IS_NB_TESTPOINTS*2));
01492 return NULL ;
01493 }
01494
01495 slope=0.0 ;
01496 spacing= ny / (IS_SKIPZONE*IS_NB_TESTPOINTS) ;
01497 yhi = (int)(ny/2) ;
01498 ylo = yhi - 1 ;
01499
01500 testpointhi = cpl_bivector_new(IS_NB_TESTPOINTS) ;
01501 testpointhi_x = cpl_bivector_get_x_data(testpointhi) ;
01502 testpointhi_y = cpl_bivector_get_y_data(testpointhi) ;
01503 testpointlo = cpl_bivector_new(IS_NB_TESTPOINTS) ;
01504 testpointlo_x = cpl_bivector_get_x_data(testpointlo) ;
01505 testpointlo_y = cpl_bivector_get_y_data(testpointlo) ;
01506 for (i=0 ; i<IS_NB_TESTPOINTS ; i++) {
01507 y = yhi + i * spacing;
01508 tmp_vector = cpl_vector_new_from_image_row(in, y+1) ;
01509 testpointhi_x[i] = y - ny / 2;
01510 testpointhi_y[i] = cpl_vector_get_median(tmp_vector) ;
01511 cpl_vector_delete(tmp_vector) ;
01512 y = ylo - i * spacing;
01513 tmp_vector = cpl_vector_new_from_image_row(in, y+1) ;
01514 testpointlo_x[IS_NB_TESTPOINTS-i-1] = y ;
01515 testpointlo_y[IS_NB_TESTPOINTS-i-1]=cpl_vector_get_median(tmp_vector) ;
01516 cpl_vector_delete(tmp_vector) ;
01517 }
01518
01519
01520 pol_coefhi = irplib_flat_fit_slope_robust(testpointhi_x,
01521 testpointhi_y, IS_NB_TESTPOINTS) ;
01522 pol_coeflo = irplib_flat_fit_slope_robust(testpointlo_x,
01523 testpointlo_y, IS_NB_TESTPOINTS) ;
01524
01525
01526 median = cpl_vector_new(IS_NB_TESTPOINTS) ;
01527 median_data = cpl_vector_get_data(median) ;
01528 for (i=0 ; i<IS_NB_TESTPOINTS ; i++) {
01529 median_data[i]=SQR(testpointhi_y[i]
01530 - pol_coefhi[0] - pol_coefhi[1] * testpointhi_x[i]);
01531 }
01532 medianerrhi = cpl_vector_get_median(median) ;
01533 for (i=0; i<IS_NB_TESTPOINTS; i++) {
01534 median_data[i]=SQR(testpointlo_y[i]
01535 - pol_coeflo[0] - pol_coeflo[1] * testpointlo_x[i]);
01536 }
01537 medianerrlo = cpl_vector_get_median(median) ;
01538 cpl_vector_delete(median) ;
01539 rampdif = testpointlo_y[IS_NB_TESTPOINTS-1] - testpointhi_y[0];
01540 slope = rampdif / (ny/2.0) ;
01541 fitslope = (pol_coefhi[1] + pol_coeflo[1]) / 2.0 ;
01542
01543 cpl_bivector_delete(testpointlo);
01544 cpl_bivector_delete(testpointhi);
01545
01546
01547 if (fabs(rampdif)<IS_MIN_RAMP ||
01548 fabs(pol_coefhi[1]) < IS_MIN_SLOPE ||
01549 fabs(pol_coeflo[1]) < IS_MIN_SLOPE ||
01550 pol_coefhi[1]/pol_coeflo[1]<0.5 ||
01551 pol_coefhi[1]/pol_coeflo[1]>2.0 ||
01552 fabs(pol_coefhi[1]-pol_coeflo[1])>IS_MAX_SLOPE_DIF ||
01553 fabs(pol_coefhi[0]-pol_coeflo[0]) > IS_MAX_INTER_DIF ||
01554 medianerrlo> IS_MAX_MNERR ||
01555 medianerrhi> IS_MAX_MNERR ||
01556 fabs(medianerrlo-medianerrhi) >IS_MAX_MNERR_DIF ||
01557 fabs(slope-fitslope) > IS_MAX_FIT_EDGE_DIF ||
01558 slope/fitslope<0.5 ||
01559 slope/fitslope>2.0) ramp_present = 0 ;
01560 else ramp_present = 1 ;
01561
01562 cpl_free(pol_coeflo) ;
01563 cpl_free(pol_coefhi) ;
01564
01565
01566 out = cpl_image_duplicate(in) ;
01567 pout = cpl_image_get_data_float(out) ;
01568 if (ramp_present == 1) {
01569 for (j=0 ; j<ny/2 ; j++) {
01570 val = slope * (j-ny/2) ;
01571 for (i=0 ; i<nx ; i++)
01572 pout[i+j*nx] -= val ;
01573 }
01574 for (j=ny/2 ; j<ny ; j++) {
01575 val = slope * (j-ny) ;
01576 for (i=0 ; i<nx ; i++)
01577 pout[i+j*nx] -= val ;
01578 }
01579
01580 }
01581
01582 return out ;
01583 }
01584