sinfoni_distortion.c

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

Generated on Wed Oct 26 13:08:55 2005 for SINFONI Pipeline Reference Manual by doxygen1.2.13.1 written by Dimitri van Heesch, © 1997-2001