sinfo_new_find_distortions.c

00001 /*
00002  * This file is part of the ESO SINFONI Pipeline
00003  * Copyright (C) 2004,2005 European Southern Observatory
00004  *
00005  * This program is free software; you can redistribute it and/or modify
00006  * it under the terms of the GNU General Public License as published by
00007  * the Free Software Foundation; either version 2 of the License, or
00008  * (at your option) any later version.
00009  *
00010  * This program is distributed in the hope that it will be useful,
00011  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013  * GNU General Public License for more details.
00014  *
00015  * You should have received a copy of the GNU General Public License
00016  * along with this program; if not, write to the Free Software
00017  * Foundation, 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
00018  */
00019 /*----------------------------------------------------------------------------
00020  
00021    File name    :       sinfo_find_distortions.c
00022    Author       :   A. Modigliani
00023    Created on   :   Aug 12, 2004
00024    Description  : 
00025 
00026  ---------------------------------------------------------------------------*/
00027 #ifdef HAVE_CONFIG_H
00028 #  include <config.h>
00029 #endif
00030 
00031 /*----------------------------------------------------------------------------
00032                                 Includes
00033  ---------------------------------------------------------------------------*/
00034 #include <irplib_distortion.h>
00035 #include "sinfo_new_find_distortions.h"
00036 #include "sinfo_pro_save.h"
00037 #include "sinfo_pro_types.h"
00038 #include "sinfo_finddist_ini_by_cpl.h"
00039 #include "sinfo_wave_calibration.h"
00040 #include "sinfo_cube_construct.h"
00041 #include "sinfo_absolute.h"
00042 #include "sinfo_distortion.h"
00043 #include "sinfo_utilities.h"
00044 #include "sinfo_recipes.h"
00045 #include "sinfo_utils_wrappers.h"
00046 #include "sinfo_error.h"
00047 #include "sinfo_globals.h"
00048 
00049 /*----------------------------------------------------------------------------
00050                                 Defines
00051  ---------------------------------------------------------------------------*/
00052 #define SINFO_ARC_SATURATION                 100000
00053 #define SINFO_ARC_MAX_WIDTH                  64
00054 
00055 
00064 /*----------------------------------------------------------------------------
00065                              Function Definitions
00066  ---------------------------------------------------------------------------*/
00067 
00068 /*----------------------------------------------------------------------------
00069    Function     :       sinfo_find_distortions()
00070    In           :       
00071    Out          :        
00072    Job          :
00073 
00074  ---------------------------------------------------------------------------*/
00075 static double
00076 new_compute_shift(double x,
00077           double y,
00078           double c00,
00079           double c10,
00080           double c01,
00081           double c11,
00082           double c20,
00083           double c02,
00084           double c21,
00085           double c12,
00086           double c30,
00087           double c03);
00088 
00089 int 
00090 sinfo_new_find_distortions(const char* plugin_id, 
00091                            cpl_parameterlist* config, 
00092                            cpl_frameset* sof)
00093 {
00094 
00095   finddist_config * cfg=NULL ;
00096   cpl_image * imonind=NULL ;
00097   cpl_image * impoly=NULL ;
00098   cpl_image * im=NULL ;
00099   cpl_image *  map=NULL ;
00100   cpl_image *  mask=NULL ;
00101 
00102   int*            degx=NULL;
00103   int*            degy=NULL;
00104   int* n_found_lines=NULL;
00105   int* sum_pointer=NULL;
00106   int** row_clean=NULL;
00107 
00108 
00109   int x_l=SIZEX/4*1;
00110   int x_u=SIZEX/4*3;
00111   int y_l=SIZEY/4*1;
00112   int y_u=SIZEY/4*3;
00113   int x_c=SIZEX/4*2;
00114   int y_c=SIZEY/4*2;
00115   const int pdx=4;
00116   const int pdy=4;
00117   int check = 0;
00118   int i = 0;
00119   int lx=0;
00120   int ly=0;
00121 
00122 
00123   int sum=0;
00124   int n_lines=0;
00125   int fit=0;
00126   int n=0;
00127   int j=0;
00128   int coef_pow[2];
00129   int arcs_window_size=0;
00130 
00131   float value=0;
00132   float newval=0;
00133   float total_dist=0;
00134   float shift=0;
00135 
00136   float* distances=NULL;
00137   float** wavelength_clean=NULL;
00138   float** slit_pos=NULL;
00139   float** acoefs=NULL;
00140   float* wave=NULL;
00141   float* intens=NULL;
00142   float* first =NULL;
00143 
00144   double*         coef=NULL;
00145 
00146   double xshift=0.;
00147   double arcs_thres_factor=0;
00148 
00149   double arcs_min_arclen_factor=0;
00150   double pcf[pdx][pdy];
00151 
00152 
00153   char tbl_name[FILE_NAME_SZ];
00154   char key_name[FILE_NAME_SZ];
00155 
00156   cpl_table* poly_tbl=NULL;
00157   cpl_table* tbl_spos=NULL;
00158   cpl_table* tbl_line_list=NULL;
00159 
00160   FitParams** par=NULL;
00161 
00162   cpl_frameset* stk=NULL;
00163   cpl_table* qclog_tbl=NULL;
00164 
00165   cpl_polynomial* distor_poly=NULL;
00166   cpl_apertures       *   arcs=NULL ;
00167   cpl_parameter* p=NULL;
00168   int smooth_rad=0;
00169 
00170   /* 
00171     parse the file names and parameters to the finddist_config 
00172     data structure cfg 
00173    */
00174   check_nomsg(stk=cpl_frameset_new());
00175   cknull(cfg = sinfo_parse_cpl_input_finddist(config,sof,&stk),
00176      "could not parse CPL input!");
00177 
00178   if(sinfo_is_fits_file (cfg->inFrame) != 1) {
00179     sinfo_msg_error("Input file %s is not FITS",cfg->inFrame);
00180     goto cleanup;
00181   }
00182 
00183   /* load the emission line frame--- */
00184   check(im = cpl_image_load(cfg->inFrame,CPL_TYPE_FLOAT,0,0),
00185     "could not load arc image");
00186 
00187   /* load the fake on - fake off frame--- */
00188   check(imonind = cpl_image_load(cfg->nsFrame,CPL_TYPE_FLOAT,0,0),
00189     "could not load on-off fake image");
00190 
00191   check(impoly  = cpl_image_load(cfg->nsFrame,CPL_TYPE_FLOAT,0,0),
00192     "could not load on-off fake image");
00193 
00194   /* load the fake on - fake off frame--- */
00195   check(mask = cpl_image_load(cfg->mask,CPL_TYPE_FLOAT,0,0),
00196     "could not load mask image");
00197 
00198   check(cpl_image_multiply (im,mask),
00199     "Failing to correct arc ima by bp mask");
00200 
00201   check(cpl_image_multiply (imonind,mask),
00202     "Failing to correct on-off fake ima by bp mask");
00203 
00204   check(cpl_image_multiply(impoly,mask),
00205     "Failing to correct on-off fake ima by bp mask");
00206 
00207   sinfo_free_image(&mask);
00208 
00209   check_nomsg(lx = cpl_image_get_size_x(im));
00210   check_nomsg(ly = cpl_image_get_size_y(im));
00211   /* TO BE CHENGED THE FOLLOWING INPUT */
00212   /* open the line list and read the number of lines */
00213   check(tbl_line_list = cpl_table_load(cfg->lineList,1,0),
00214     "problems loading table %s",cfg->lineList);
00215 
00216   check_nomsg(n_lines = cpl_table_get_nrow(tbl_line_list));
00217   check_nomsg(wave=cpl_table_get_data_float(tbl_line_list,"wave"));
00218   check_nomsg(intens = cpl_table_get_data_float(tbl_line_list,"int"));
00219 
00220   /*
00221        #---------------------------------------------------------------------
00222        #---------------------------FINDLINES---------------------------------
00223        #---------------------------------------------------------------------
00224   */
00225 
00226   sinfo_msg("Find Lines");
00227   /* 
00228        do the wavelength calibration, that means: 
00229        find the dispersion relation and parameterize its coefficients
00230   */
00231 
00232   acoefs  = sinfo_new_2Dfloatarray(cfg->nrDispCoefficients, lx);
00233   /* allocate memory */
00234   n_found_lines    = sinfo_new_intarray(lx); 
00235   row_clean        = sinfo_new_2Dintarray(lx, n_lines);
00236   wavelength_clean = sinfo_new_2Dfloatarray(lx, n_lines);
00237   sum_pointer      = sinfo_new_intarray(1);
00238  
00239   /* find the emission lines in each image column */
00240   sinfo_new_intarray_set_value(sum_pointer, 0, 0);
00241 
00242   ck0(check = sinfo_new_find_lines(im, 
00243                    wave, 
00244                    intens,
00245                    n_lines,
00246                    row_clean,
00247                    wavelength_clean,
00248                    cfg->guessBeginWavelength,
00249                    cfg->guessDispersion1,
00250                    cfg->guessDispersion2,
00251                    cfg->mindiff,
00252                    cfg->halfWidth,
00253                    n_found_lines,
00254                    cfg->sigma,
00255                    sum_pointer ),
00256       "FindLines failed!");
00257 
00258   sinfo_free_table(&tbl_line_list);
00259  
00260   /*
00261    #-------------------------------------------------------------------------
00262    #---------------------------WAVECALIB-------------------------------------
00263    #-------------------------------------------------------------------------
00264     */
00265   sinfo_msg("Do wave calib");
00266   sum = sinfo_new_intarray_get_value(sum_pointer,0);
00267 
00268   /* #allocate memory for the fit parameters */
00269   cknull(par = sinfo_new_fit_params( sum ),
00270      "sinfo_new_fit_params failed!");
00271 
00272   /* 
00273       fit each line, make a polynomial fit and fit the resulting fit 
00274        coefficients across the columns of the slitlet
00275   */
00276   cknull(map = sinfo_new_wave_cal(im, 
00277                           par,
00278                       acoefs,
00279                       cfg->nslitlets,
00280                       row_clean,
00281                       wavelength_clean,
00282                       n_found_lines,
00283                       cfg->guessDispersion1,
00284                       cfg->halfWidth,
00285                       cfg->minAmplitude,
00286                       cfg->maxResidual,
00287                       cfg->fwhm,
00288                       cfg->nrDispCoefficients, 
00289                       cfg->nrCoefCoefficients,
00290                       cfg->sigmaFactor,
00291                       cfg->pixeldist,
00292                       cfg->pixel_tolerance),
00293      "sinfo_waveCal failed!");
00294   /* here mem leak */
00295 
00296   sinfo_msg("Check line positions");
00297   shift = sinfo_new_check_line_positions (im, acoefs, 
00298                                           cfg->nrDispCoefficients, 
00299                                           cfg->guessDispersion1,par);
00300 
00301   if (FLAG == shift) {
00302     sinfo_msg_error ( "checkForLinePositions failed!");
00303     goto cleanup;
00304   }
00305   sinfo_free_image(&map);
00306 
00307 
00308  /* free memory */
00309   sinfo_new_destroy_2Dfloatarray ( &wavelength_clean, lx );
00310   sinfo_new_destroy_2Dintarray (&row_clean, lx);
00311   sinfo_new_destroy_intarray(&n_found_lines );
00312   sinfo_new_destroy_intarray(&sum_pointer );
00313   sinfo_new_destroy_2Dfloatarray ( &acoefs, cfg->nrDispCoefficients );
00314   /*
00315   sinfo_new_destroy_array (& wave );
00316   sinfo_new_destroy_array (& intens );
00317   */
00318   /*
00319   #----------------------------------------------------------------------------
00320   #-------------------------SLITFITS-------------------------------------------
00321   #----------------------------------------------------------------------------
00322   #--fit the slitlet sinfo_edge positions if desired--
00323   */
00324 
00325   sinfo_msg("Find slit pos");
00326   /*allocate memory for the slitlet position array */
00327   slit_pos = sinfo_new_2Dfloatarray(32,2);
00328 
00329   fit = sinfo_new_fit_slits_boltz( im, par, slit_pos, cfg->boxLength, 
00330                       cfg->yBox, cfg->diffTol );
00331 
00332 
00333   if (fit < 0) {
00334     sinfo_msg_error( "sinfo_fitSlitsBoltz failed" );
00335     goto cleanup;
00336   }
00337   sinfo_new_destroy_fit_params(&par); 
00338 
00339   check_nomsg(tbl_spos=cpl_table_new(32));
00340   check_nomsg(cpl_table_new_column(tbl_spos,"pos1", CPL_TYPE_DOUBLE));
00341   check_nomsg(cpl_table_new_column(tbl_spos,"pos2", CPL_TYPE_DOUBLE));
00342   check_nomsg(cpl_table_set_column_format(tbl_spos,"pos1", "15.9f"));
00343   check_nomsg(cpl_table_set_column_format(tbl_spos,"pos2", "15.9f"));
00344 
00345   for (i =0; i< 32; i++) {
00346       check_nomsg(cpl_table_set_double(tbl_spos,"pos1",i,
00347                                 sinfo_new_array2D_get_value(slit_pos,i,0)));
00348       check_nomsg(cpl_table_set_double(tbl_spos,"pos2",i,
00349                                 sinfo_new_array2D_get_value(slit_pos,i,1)));
00350      
00351   }
00352   strcpy(tbl_name,"out_slitlets_pos_predist.fits");
00353 
00354   ck0(sinfo_pro_save_tbl(tbl_spos,stk,sof,tbl_name,
00355              PRO_SLITLETS_POS_PREDIST,NULL,plugin_id,config),
00356       "cannot save tbl %s", tbl_name);
00357 
00358   sinfo_free_table(&tbl_spos);
00359 
00360   /*
00361   #---------------------------------------------------------
00362   # do the north - south - test
00363   #---------------------------------------------------------
00364   */
00365   sinfo_msg("Do north south test");
00366 
00367   sinfo_msg("cfg->nslits =%d\n", cfg->nslits);
00368   sinfo_msg("cfg->nshalfWidth =%d\n", cfg->nshalfWidth);
00369   sinfo_msg("cfg->nsfwhm=%f\n",cfg->nsfwhm ); 
00370   sinfo_msg("cfg->minDiff=%f\n", cfg->minDiff);
00371   sinfo_msg("cfg->estimated_dist=%f\n", cfg->estimated_dist);
00372   sinfo_msg("cfg->devtol=%f\n", cfg->devtol);
00373   sinfo_msg("cfg->loPos=%d\n", cfg->loPos);
00374   sinfo_msg("cfg->hiPos=%d\n",cfg->hiPos);
00375 
00376   cknull(distances = sinfo_north_south_test(imonind, 
00377                       cfg->nslits,
00378                       cfg->nshalfWidth,
00379                       cfg->nsfwhm ,
00380                       cfg->minDiff,
00381                       cfg->estimated_dist,
00382                       cfg->devtol,
00383                       cfg->loPos,
00384                       cfg->hiPos),
00385      "north_south_test failed");
00386 
00387   sinfo_free_image(&imonind);
00388 
00389   /*
00390   #---------------------------------------------------------
00391   # get firstcol
00392   #---------------------------------------------------------
00393   */
00394   sinfo_msg("get first col");
00395   first = sinfo_new_floatarray(61);
00396   total_dist=0;
00397   n=0;
00398 
00399   for (i=0; i< 31; i++) {
00400     total_dist=total_dist+sinfo_new_array_get_value(distances,i);
00401 
00402     for (j=0; j< 2; j++) {
00403       if (j == 0) {
00404         if (i != 30) {
00405       newval = sinfo_new_array2D_get_value(slit_pos,i+1,j) - total_dist;
00406       sinfo_new_array_set_value(first, newval, n);
00407       n = n+1;
00408         }
00409       }
00410       if (j == 1) {
00411         newval = sinfo_new_array2D_get_value(slit_pos,i,j) - total_dist;
00412         sinfo_new_array_set_value(first, newval, n);
00413         n = n+1;
00414       }
00415     }
00416   }
00417   value = sinfo_new_f_median(first,61);
00418   sinfo_msg("Firstcol =%f", value);
00419   sinfo_new_destroy_array(&first);
00420   sinfo_new_destroy_array(&distances);
00421   sinfo_new_destroy_2Dfloatarray ( &slit_pos, 32 );
00422  
00423   /*
00424   #---------------------------------------------------------
00425   # Determine distortions
00426   #---------------------------------------------------------
00427   */
00428   sinfo_msg("Determine distortions");
00429   degx=cpl_calloc(8,sizeof(int));
00430   degy=cpl_calloc(8,sizeof(int));
00431   coef=cpl_calloc(8,sizeof(double));
00432 
00433   check_nomsg(p=cpl_parameterlist_find(config,
00434                 "sinfoni.distortion.arcs_thresh_factor"));
00435   check_nomsg(arcs_thres_factor=cpl_parameter_get_double(p));
00436 
00437   check_nomsg(p=cpl_parameterlist_find(config,
00438                 "sinfoni.distortion.arcs_min_arclen_factor"));
00439   check_nomsg(arcs_min_arclen_factor=cpl_parameter_get_double(p));
00440 
00441   check_nomsg(p=cpl_parameterlist_find(config,
00442                 "sinfoni.distortion.arcs_window_size"));
00443   check_nomsg(arcs_window_size=cpl_parameter_get_int(p));
00444 
00445   check_nomsg(p=cpl_parameterlist_find(config,
00446                 "sinfoni.distortion.smooth_rad"));
00447   check_nomsg(smooth_rad=cpl_parameter_get_int(p));
00448 
00449 
00450   cknull(distor_poly=sinfo_distortion_estimate_new(impoly, 
00451                            0, 
00452                            0,
00453                                                    cpl_image_get_size_x(impoly), 
00454                                                    cpl_image_get_size_y(impoly), 
00455                                                    FALSE, 
00456                                                    SINFO_ARC_SATURATION,
00457                            SINFO_ARC_MAX_WIDTH, 
00458                            arcs_thres_factor, 
00459                                                    arcs_min_arclen_factor,
00460                                                    arcs_window_size,
00461                                                    smooth_rad,
00462                                                    3,
00463                                                    (double)value,
00464                                                    &arcs),
00465      "cannot estimage distortion") ;
00466   /*AMo: additional mem leaks */
00467 
00468   sinfo_free_apertures(&arcs);
00469   coef_pow[0]=0;
00470   coef_pow[1]=0;
00471   check_nomsg(pcf[0][0]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
00472   sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",0,0,pcf[0][0]);
00473   /*
00474   pcf[0][0]+=value;
00475   */
00476   sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",0,0,pcf[0][0]);
00477   check_nomsg(cpl_polynomial_set_coeff(distor_poly,coef_pow,pcf[0][0]));
00478 
00479   /* up to here ok */
00480   /* To check mem leaks */
00481   sinfo_free_image(&impoly);
00482 
00483   sinfo_msg("Polynomial fit results: deg=%d",
00484         cpl_polynomial_get_degree(distor_poly));
00485 
00486  
00487   coef_pow[0]=0;
00488   coef_pow[1]=0;
00489   check_nomsg(pcf[0][0]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
00490   sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",0,0,pcf[0][0]);
00491 
00492   coef_pow[0]=1;
00493   coef_pow[1]=0;
00494   check_nomsg(pcf[1][0]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
00495   sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",1,0,pcf[1][0]);
00496 
00497   coef_pow[0]=0;
00498   coef_pow[1]=1;
00499   check_nomsg(pcf[0][1]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
00500   sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",0,1,pcf[0][1]);
00501 
00502   coef_pow[0]=1;
00503   coef_pow[1]=1;
00504   check_nomsg(pcf[1][1]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
00505   sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",1,1,pcf[1][1]);
00506 
00507   coef_pow[0]=2;
00508   coef_pow[1]=0;
00509   check_nomsg(pcf[2][0]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
00510   sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",2,0,pcf[2][0]);
00511 
00512   coef_pow[0]=0;
00513   coef_pow[1]=2;
00514   check_nomsg(pcf[0][2]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
00515   sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",0,2,pcf[0][2]);
00516 
00517   coef_pow[0]=2;
00518   coef_pow[1]=1;
00519   check_nomsg(pcf[2][1]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
00520   sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",2,1,pcf[2][1]);
00521 
00522   coef_pow[0]=1;
00523   coef_pow[1]=2;
00524   check_nomsg(pcf[1][2]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
00525   sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",1,2,pcf[1][2]);
00526 
00527   coef_pow[0]=3;
00528   coef_pow[1]=0;
00529   check_nomsg(pcf[3][0]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
00530   sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",3,0,pcf[3][0]);
00531 
00532   coef_pow[0]=0;
00533   coef_pow[1]=3;
00534   check_nomsg(pcf[0][3]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
00535   sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",0,3,pcf[0][3]);
00536   
00537   sinfo_msg("Save results");
00538 
00539   check_nomsg(poly_tbl=cpl_table_new(10));
00540   check_nomsg(cpl_table_new_column(poly_tbl,"degx", CPL_TYPE_INT));
00541   check_nomsg(cpl_table_new_column(poly_tbl,"degy", CPL_TYPE_INT));
00542   check_nomsg(cpl_table_new_column(poly_tbl,"coeff", CPL_TYPE_DOUBLE));
00543   check_nomsg(cpl_table_set_int(poly_tbl,"degx",0,0));
00544   check_nomsg(cpl_table_set_int(poly_tbl,"degx",1,1));
00545   check_nomsg(cpl_table_set_int(poly_tbl,"degx",2,0));
00546   check_nomsg(cpl_table_set_int(poly_tbl,"degx",3,1));
00547   check_nomsg(cpl_table_set_int(poly_tbl,"degx",4,2));
00548   check_nomsg(cpl_table_set_int(poly_tbl,"degx",5,0));
00549   check_nomsg(cpl_table_set_int(poly_tbl,"degx",6,2));
00550   check_nomsg(cpl_table_set_int(poly_tbl,"degx",7,1));
00551   check_nomsg(cpl_table_set_int(poly_tbl,"degx",8,3));
00552   check_nomsg(cpl_table_set_int(poly_tbl,"degx",9,0));
00553 
00554   check_nomsg(cpl_table_set_int(poly_tbl,"degy",0,0));
00555   check_nomsg(cpl_table_set_int(poly_tbl,"degy",1,0));
00556   check_nomsg(cpl_table_set_int(poly_tbl,"degy",2,1));
00557   check_nomsg(cpl_table_set_int(poly_tbl,"degy",3,1));
00558   check_nomsg(cpl_table_set_int(poly_tbl,"degy",4,0));
00559   check_nomsg(cpl_table_set_int(poly_tbl,"degy",5,2));
00560   check_nomsg(cpl_table_set_int(poly_tbl,"degy",6,1));
00561   check_nomsg(cpl_table_set_int(poly_tbl,"degy",7,2));
00562   check_nomsg(cpl_table_set_int(poly_tbl,"degy",8,0));
00563   check_nomsg(cpl_table_set_int(poly_tbl,"degy",9,3));
00564 
00565   check_nomsg(cpl_table_set_double(poly_tbl,"coeff",0,pcf[0][0]));
00566   check_nomsg(cpl_table_set_double(poly_tbl,"coeff",1,pcf[1][0]));
00567   check_nomsg(cpl_table_set_double(poly_tbl,"coeff",2,pcf[0][1]));
00568   check_nomsg(cpl_table_set_double(poly_tbl,"coeff",3,pcf[1][1]));
00569   check_nomsg(cpl_table_set_double(poly_tbl,"coeff",4,pcf[2][0]));
00570   check_nomsg(cpl_table_set_double(poly_tbl,"coeff",5,pcf[0][2]));
00571   check_nomsg(cpl_table_set_double(poly_tbl,"coeff",6,pcf[2][1]));
00572   check_nomsg(cpl_table_set_double(poly_tbl,"coeff",7,pcf[1][2]));
00573   check_nomsg(cpl_table_set_double(poly_tbl,"coeff",8,pcf[3][0]));
00574   check_nomsg(cpl_table_set_double(poly_tbl,"coeff",9,pcf[0][3]));
00575 
00576   sinfo_msg("Polynomial distortion coefficients \n"); 
00577   sinfo_msg("c= %g   %g  %g  %g  %g  %g  %g  %g %g %g\n",
00578          pcf[0][0],pcf[1][0],pcf[0][1],pcf[1][1],
00579              pcf[2][0],pcf[0][2],pcf[2][1],pcf[1][2],pcf[3][0],pcf[0][3]);
00580 
00581   /* QC LOG */
00582   check_nomsg(qclog_tbl = sinfo_qclog_init());
00583   snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%d","QC COEFF",0,0);
00584   ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,pcf[0][0],
00585              "Polynomial distortion coefficient","%g"));
00586 
00587   snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%d","QC COEFF",1,0);
00588   ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,pcf[1][0],
00589              "Polynomial distortion coefficient","%g"));
00590 
00591   snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%d","QC COEFF",0,1);
00592   ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,pcf[0][1],
00593              "Polynomial distortion coefficient","%g"));
00594 
00595   snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%d","QC COEFF",1,1);
00596   ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,pcf[1][1],
00597              "Polynomial distortion coefficient","%g"));
00598 
00599   snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%d","QC COEFF",2,0);
00600   ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,pcf[2][0],
00601              "Polynomial distortion coefficient","%g"));
00602 
00603   snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%d","QC COEFF",0,2);
00604   ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,pcf[0][2],
00605              "Polynomial distortion coefficient","%g"));
00606 
00607   snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%d","QC COEFF",2,1);
00608   ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,pcf[2][1],
00609              "Polynomial distortion coefficient","%g"));
00610 
00611   snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%d","QC COEFF",1,2);
00612   ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,pcf[1][2],
00613              "Polynomial distortion coefficient","%g"));
00614 
00615   snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%d","QC COEFF",3,0);
00616   ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,pcf[3][0],
00617              "Polynomial distortion coefficient","%g"));
00618 
00619   snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%d","QC COEFF",0,3);
00620   ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,pcf[0][3],
00621              "Polynomial distortion coefficient","%g"));
00622 
00623   snprintf(key_name,MAX_NAME_SIZE-1,"%s","QC OFFSET");
00624   ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,value,
00625              "Polynomial distortion coefficient","%g"));
00626  
00627   xshift=new_compute_shift(x_c,y_c,pcf[0][0],pcf[1][0],pcf[0][1],
00628                                    pcf[1][1],pcf[2][0],pcf[0][2],
00629                                    pcf[2][1],pcf[1][2],pcf[3][0],pcf[0][3]);
00630  
00631   ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC XSHIFT CC",xshift,
00632                                    "X shift in x_c,y_c","%g"));
00633 
00634 
00635   xshift=new_compute_shift(x_l,y_l,pcf[0][0],pcf[1][0],pcf[0][1],
00636                                    pcf[1][1],pcf[2][0],pcf[0][2],
00637                                    pcf[2][1],pcf[1][2],pcf[3][0],pcf[0][3]);
00638 
00639   ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC XSHIFT LL",xshift,
00640                                    "X shift in x_l,y_l","%g"));
00641 
00642   xshift=new_compute_shift(x_l,y_u,pcf[0][0],pcf[1][0],pcf[0][1],
00643                                    pcf[1][1],pcf[2][0],pcf[0][2],
00644                                    pcf[2][1],pcf[1][2],pcf[3][0],pcf[0][3]);
00645   
00646   ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC XSHIFT UL",xshift,
00647                                    "X shift in x_l,y_u","%g"));
00648 
00649   xshift=new_compute_shift(x_u,y_u,pcf[0][0],pcf[1][0],pcf[0][1],
00650                                    pcf[1][1],pcf[2][0],pcf[0][2],
00651                                    pcf[2][1],pcf[1][2],pcf[3][0],pcf[0][3]);
00652   
00653   ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC XSHIFT UR",xshift,
00654                                    "X shift in x_u,y_u","%g"));
00655 
00656 
00657   xshift=new_compute_shift(x_u,y_l,pcf[0][0],pcf[1][0],pcf[0][1],
00658                                    pcf[1][1],pcf[2][0],pcf[0][2],
00659                                    pcf[2][1],pcf[1][2],pcf[3][0],pcf[0][3]);
00660 
00661   ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC XSHIFT LR",xshift,
00662                                    "X shift in x_u,y_l","%g"));
00663 
00664 
00665   ck0(sinfo_pro_save_tbl(poly_tbl,stk,sof,cfg->outName,
00666              PRO_DISTORTION,qclog_tbl,plugin_id,config),
00667       "cannot dump tbl %s", cfg->outName);
00668 
00669   sinfo_free_table(&poly_tbl);
00670   sinfo_free_table(&qclog_tbl);
00671   sinfo_free_polynomial(&distor_poly);
00672   sinfo_free_int(&degx);
00673   sinfo_free_int(&degy);
00674   sinfo_free_double(&coef);
00675   sinfo_free_image(&im);
00676   sinfo_free_frameset(&stk);
00677   sinfo_finddist_free (&cfg);
00678  
00679   return 0;
00680 
00681  cleanup:
00682   sinfo_free_table(&poly_tbl);
00683   sinfo_free_table(&qclog_tbl);
00684   sinfo_free_polynomial(&distor_poly);
00685   sinfo_free_int(&degx);
00686   sinfo_free_int(&degy);
00687   sinfo_free_double(&coef);
00688   sinfo_free_apertures(&arcs);
00689 
00690 
00691 
00692   /*if(wave != NULL) sinfo_new_destroy_array (& wave );*/
00693   /*if(intens != NULL) sinfo_new_destroy_array (& intens );*/ 
00694   if(first != NULL) sinfo_new_destroy_array(&first);
00695   sinfo_free_table(&tbl_spos);
00696   if(slit_pos != NULL) sinfo_new_destroy_2Dfloatarray (&slit_pos,32);
00697   if(distances != NULL) sinfo_new_destroy_array(&distances); 
00698   if(par!=NULL) sinfo_new_destroy_fit_params(&par); 
00699   if(n_found_lines != NULL) sinfo_new_destroy_intarray(&n_found_lines );
00700   if(row_clean != NULL) sinfo_new_destroy_2Dintarray(&row_clean, lx);
00701   if(wavelength_clean != NULL) {
00702       sinfo_new_destroy_2Dfloatarray(&wavelength_clean,lx);
00703   }
00704   if(sum_pointer != NULL) sinfo_new_destroy_intarray(&sum_pointer);
00705   if(acoefs != NULL) {
00706       sinfo_new_destroy_2Dfloatarray(&acoefs, cfg->nrDispCoefficients );
00707   }
00708   sinfo_free_table(&tbl_line_list);
00709   sinfo_free_image(&mask);
00710   sinfo_free_image(&impoly);
00711   sinfo_free_image(&imonind);
00712   sinfo_free_image(&map);
00713   sinfo_free_image(&im);
00714   sinfo_finddist_free (&cfg);
00715   sinfo_free_frameset(&stk);
00716   return -1;
00717 
00718 }
00719 
00720 
00721 static double
00722 new_compute_shift(double x,
00723               double y,
00724               double c00,
00725               double c10,
00726               double c01,
00727               double c11,
00728               double c20,
00729               double c02,
00730               double c21,
00731               double c12,
00732               double c30,
00733               double c03)
00734 {
00735 
00736   double x_shift=0;
00737   double shift=0;
00738     sinfo_msg("c= %g   %g  %g  %g  %g  %g  %g  %g %g %g\n",
00739          c00,c10,c01,c11,
00740              c20,c02,c21,c12,c30,c03);
00741 
00742   shift=c00+c10*x+c01*y+
00743         c11*x*y+c20*x*x+c02*y*y+
00744         c21*x*x*y+c12*x*y*y+c30*x*x*x+c03*y*y*y;
00745   x_shift=x-shift;
00746   return x_shift;
00747 
00748 
00749 
00750 }

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