sinfo_distortion.c

00001 /* $Id: sinfo_distortion.c,v 1.16 2007/01/10 14:26:05 amodigli Exp $
00002  *
00003  * This file is part of the irplib package
00004  * Copyright (C) 2002,2003 European Southern Observatory
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00019  */
00020 
00021 /*
00022  * $Author: amodigli $
00023  * $Date: 2007/01/10 14:26:05 $
00024  * $Revision: 1.16 $
00025  * $Name:  $
00026  */
00027 
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031 
00032 /*----------------------------------------------------------------------------
00033                                    Includes
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                                    Define
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   /* 1.1-2 */
00054 #define ARC_MINNBARCS       32     /* 4-32 */
00055 #define ARC_RANGE_FACT      3.0
00056 #define ARC_WINDOWSIZE      10  /* 32 */
00057 
00058 #define TRESH_MEDIAN_MIN    0.0
00059 #define TRESH_SIGMA_MAX     200.0
00060 
00061 /*---------------------------------------------------------------------------*/
00065 /*---------------------------------------------------------------------------*/
00066 
00067 /*----------------------------------------------------------------------------
00068                                 Functions prototypes
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                                 Functions code
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     /* AMO added to use offset */
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     /* Check entries */
00338     if (org == NULL) return NULL ;
00339     if (kappa < 0.0) return NULL ;
00340 
00341     /* Initialise */
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         /* Local copy of input image */
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       //cpl_image_save(local_im,"out_local_im.fits",CPL_BPP_DEFAULT,NULL,CPL_IO_DEFAULT);
00358       //local_im=sinfo_distortion_image_restore(local_im,smooth_rad,1,2,0,2);
00359       //cpl_image_save(local_im,"out_local_im_post.fits",CPL_BPP_DEFAULT,NULL,CPL_IO_DEFAULT);
00360     
00361     }
00362     /* Detect the arcs in the input image */
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     /* Check that the arcs are not concentrated in the same zone */
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     /* Create a 2-D deformation grid with detected arcs */
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     /* Prepare the fitting */
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     /* AMO new to use offset */
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             sinfo_msg("diff=%g",(plines_pos_tmp[i+k]-
00452                  plines_pos_tmp[k]));
00453         */
00454       }
00455     }
00456       /*
00457       sinfo_msg("j=%d pgrid_tot=%g",j,pgrid_tot[j]);
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          sinfo_msg_error("AMo after corr grid[%d,%d]=%g",
00467                           i,k,pgridx[k+i*n_calib]);
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     /* end AMO: to use the offset */    
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     /* Apply the fitting */
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     /* Free and return */
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     /* Default values for output parameters */
00546     *label_im = NULL ;
00547 
00548     /* Clear zones to be ignored (to avoid false detections) */
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     /* Median vertical filter */
00557     filter = cpl_matrix_new(3, 1) ;
00558     cpl_matrix_fill(filter, 1.0) ;
00559     /* filt_im = cpl_image_filter_median(im, filter) ; */
00560     filt_im = cpl_image_duplicate(im) ;
00561     cpl_matrix_delete(filter) ;
00562 
00563     /* Subtract a low-pass */
00564     /* AMO: suppressed as may remove arcs */
00565     if (sinfo_distortion_sub_hor_lowpass(filt_im, arcs_window_size) == -1) {
00566         cpl_image_delete(filt_im) ;
00567         return NULL ;
00568     }
00569     //cpl_image_save(filt_im,"out_filt_im_lp.fits",CPL_BPP_DEFAULT,
00570     //               NULL,CPL_IO_DEFAULT);
00571 
00572     /* Get relevant stats for thresholding */
00573     median_val = cpl_image_get_median_dev(filt_im, &sigma) ;
00574 
00575     /* Correct median_val and sigma if necessary */
00576     if (median_val < TRESH_MEDIAN_MIN) median_val = TRESH_MEDIAN_MIN ;
00577     if (sigma > TRESH_SIGMA_MAX) sigma = TRESH_SIGMA_MAX ;
00578 
00579     /* Set the threshold */
00580     threshold = median_val + sigma * kappa ;
00581 
00582     /* Collapse the image */
00583     collapsed = cpl_image_collapse_median_create(filt_im, 0, 0, 0) ;
00584 
00585     /* Threshold to keep only the arcs - use of the collapsed image */
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     /* Binarize the image */
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     /* Test if there are enough good pixels */
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     /* Apply a morphological closing to clean the isolated pixels */
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     /* Labelize pixel map to a label image */
00619     *label_im = cpl_image_labelise_mask_create(bin_im, &nobj) ;
00620     cpl_mask_delete(bin_im) ;
00621     //cpl_image_save(*label_im,"out_label_im.fits",CPL_BPP_DEFAULT,
00622     //               NULL,CPL_IO_DEFAULT);
00623 
00624     /* Compute statistics on objects */
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     /* Set min_arclen */
00632     min_arclen = (int)(ny /arcs_min_arclen_factor) ;
00633     //cpl_image_save(im,"out_im.fits",CPL_BPP_DEFAULT,NULL,CPL_IO_DEFAULT);
00634 
00635     /* Purge non-relevant arcs */
00636     /* cpl_apertures_dump(det,stdout); */
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     /* cpl_apertures_dump(det,stdout); */ 
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     /* Return  */
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     /* AMO added to use offset */
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     /* Check entries */
00731     if (org == NULL) return NULL ;
00732 
00733     /* Initialise */
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         /* Local copy of input image */
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     /* Detect the arcs in the input image */
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     /* Check that the arcs are not concentrated in the same zone */
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     /* Create a 2-D deformation grid with detected arcs */
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     /* Prepare the fitting */
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       sinfo_msg("pgridx=%g pgridy=%g pvalues=%g",
00815           pgridx[j+i*n_calib],pgridy[j+i*n_calib],pvalues_to_fit[j+i*n_calib]);
00816 */     
00817         }
00818     }
00819 
00820 
00821     /* AMO new to use offset */
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             sinfo_msg("diff=%g",(plines_pos_tmp[i+k]-
00843                  plines_pos_tmp[k]));
00844         */
00845       }
00846     }
00847       /*
00848       sinfo_msg("j=%d pgrid_tot=%g",j,pgrid_tot[j]);
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          sinfo_msg_error("AMo after corr grid[%d,%d]=%g",
00858                           i,k,pgridx[k+i*n_calib]);
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     /* end AMO: to use the offset */    
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     /* Apply the fitting */
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     /* Free and return */
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     /* Default values for output parameters */
00932     *label_im = NULL ;
00933 
00934     /* Clear zones to be ignored (to avoid false detections) */
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     /* Median vertical filter */
00944     filter = cpl_matrix_new(3, 1) ;
00945     cpl_matrix_fill(filter, 1.0) ;
00946     /* filt_im = cpl_image_filter_median(im, filter) ; */
00947     filt_im = cpl_image_duplicate(im) ;
00948     cpl_matrix_delete(filter) ;
00949 
00950     /* Subtract a low-pass */
00951     if (sinfo_distortion_sub_hor_lowpass(filt_im, ARC_WINDOWSIZE) == -1) {
00952         cpl_image_delete(filt_im) ;
00953         return NULL ;
00954     }
00955     
00956     /* Get relevant stats for thresholding */
00957     median_val = cpl_image_get_median_dev(filt_im, &sigma) ;
00958 
00959     /* Correct median_val and sigma if necessary */
00960     if (median_val < TRESH_MEDIAN_MIN) median_val = TRESH_MEDIAN_MIN ;
00961     if (sigma > TRESH_SIGMA_MAX) sigma = TRESH_SIGMA_MAX ;
00962 
00963     /* Set the threshold */
00964     threshold = median_val + sigma * ARC_THRESHFACT ;
00965 
00966     /* Collapse the image */
00967     collapsed = cpl_image_collapse_median_create(filt_im, 0, 0, 0) ;
00968 
00969     /* Threshold to keep only the arcs - use of the collapsed image */
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     /* Binarize the image */
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     /* Test if there are enough good pixels */
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     /* Apply a morphological closing to clean the isolated pixels */
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     /* Labelize pixel map to a label image */
01003     *label_im = cpl_image_labelise_mask_create(bin_im, &nobj) ;
01004     cpl_mask_delete(bin_im) ;
01005 
01006     /* Compute statistics on objects */
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     /* Set min_arclen */
01015     min_arclen = (int)(ny / ARC_MINARCLENFACT) ;
01016 
01017     /* Purge non-relevant arcs */
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     /* Return  */
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     /* Check entries */
01052     if (im == NULL) return -1 ;
01053     if (cpl_image_get_type(im) != CPL_TYPE_FLOAT) return -1 ;
01054 
01055     /* Get the data */
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     /* Fill the zone */
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     /* Check entries */
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     /* Get access to the im / im1d data */
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     /* Apply the thresholding */
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     /* Return */
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     /* Test entries */
01118     if (im == NULL) return -1 ;
01119     if (filt_size <= 0) return -1 ;
01120     
01121     /* Initialise */
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     /* Get the vectors out of the image */
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     /* Filter the vectors */
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     /* Average the filtered vectors to get the low freq signal */
01154     cpl_vector_add(avglinehi, avglinelo) ;
01155     cpl_vector_delete(avglinelo) ;
01156     cpl_vector_divide_scalar(avglinehi, 2.0) ;
01157 
01158     /* Subtract the low frequency signal */
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     /* Check entries */
01197     if (arcs == NULL) return -1 ;
01198     if (*arcs == NULL) return -1 ;
01199     if (*lab_im == NULL) return -1 ;
01200 
01201     /* Get number of arcs */
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     /* Allocate selection array */
01207     selection = cpl_malloc(nb_arcs * sizeof(int)) ;
01208     /* Loop on the different arcs candidates */
01209     /* sinfo_msg("min_arclen=%d max_arcwidth=%d",min_arclen,max_arcwidth); */
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         /* Test if the current object is a valid arc */
01219       
01220         if (
01221             (arclen>min_arclen) && 
01222         (arcwidth<max_arcwidth) && 
01223             (edge>0) &&
01224             (mean < arc_sat)) {
01225       /*
01226         sinfo_msg_warning("Take Pos=%5.4d len=%d width=%d edge=%d mean=%f ",
01227     (cpl_apertures_get_right(*arcs, i+1)+cpl_apertures_get_left(*arcs, i+1))/2,
01228      arclen,arcwidth,edge,mean);
01229       */
01230             selection[i] = 1 ;
01231         } else {
01232       /*
01233     sinfo_msg_warning("Rej Pos=%5.4d len=%d width=%d edge=%d mean=%f i=%d",
01234          (cpl_apertures_get_right(*arcs, i+1)+
01235           cpl_apertures_get_left(*arcs, i+1))/2,arclen,arcwidth,edge,mean,i);
01236       */
01237             selection[i] = 0 ;
01238         }
01239     }
01240 
01241     /* Update the labelised image by erasing non valid arcs */
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     /* Reset the labels to have consecutive ones */
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     /* Purge the bad arcs */
01259     cpl_apertures_delete(*arcs) ;
01260     *arcs = cpl_apertures_new_from_image(im, *lab_im) ;
01261 
01262     /* Check if there are some valid arcs */
01263     if (cpl_apertures_get_size(*arcs) <= 0) {
01264         cpl_msg_error(fctid, "No valid arc found") ;
01265         return -1 ;
01266     }
01267     /* Return  */
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     /* Check entries */
01297 
01298     /* Initialise */
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     /* Allocate positions (pos. of n_arcs*nb_samples pts on the arcs) */
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     /* Median filter on input image */
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     /* Measured Arcs coordinates along curvature */
01314     arcs_samples_y = cpl_malloc(n_arcs * nb_samples * sizeof(int)) ;
01315     computed = cpl_calloc(n_arcs*nb_samples, sizeof(int)) ;
01316 
01317     /* Find out the Y coordinates along the arcs  */
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     /* Find out the X coord. at nb_samples Y positions on all arcs */
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             /* use_this_arc is set to 1 if we are on the arc at a y */
01333             /* coordinate where the x coord should be found */
01334             obj = plabel_im[i + j * nx] ;
01335             /* Handle background */
01336             if (obj==0) continue ;
01337             /* Decrease by one to index the array from 0 */
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                 /* Find x coordinate of obj at the Y coord. */
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     /* Free and return */
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     /* Initialize */
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     /* While we stay in the same object... */
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         /* Next point */
01422         start_pos++ ;
01423 
01424         curr_obj = plabel_im[start_pos + y * nx] ;
01425     } while (curr_obj == objnum) ;
01426 
01427     /* Returned position is the gravity center or the max in bad cases */
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     /* Return */
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     /* Initialise */
01483     nx = cpl_image_get_size_x(in) ;
01484     ny = cpl_image_get_size_y(in) ;
01485                     
01486     /* Check entries */
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     /* Fill the vectors */
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     /* Apply the fit */
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     /* Compute the errors */
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     /* Decide if there is a ramp or not  */
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     /* Correct the ramp if it is there */
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 

Generated on Wed Jan 17 08:33:42 2007 for SINFONI Pipeline Reference Manual by  doxygen 1.4.4