sinfo_image_ops.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 * M.P.E. - SPIFFI project
00021 *
00022 *
00023 *
00024 * who       when      what
00025 * --------  --------  ----------------------------------------------
00026 * rabuter 2004-12-03 support one dimensional image in sinfo_shiftImage
00027 * schreib  23/05/00  created
00028 */
00029 
00030 /************************************************************************
00031 *   NAME
00032 *        sinfo_image_ops.c -
00033 *        image arithmetic routines 
00034 *
00035 *   SYNOPSIS
00036 *   #include "sinfo_image_ops.h"
00037 *
00038 *   1) Vector * sinfo_new_mean_of_columns( cpl_image *im )
00039 *   2) Vector * sinfo_new_clean_mean_of_columns( cpl_image *im, 
00040 *                                   double lo_reject, 
00041 *                                   double hi_reject)
00042 *   3) cpl_image * sinfo_new_div_image_by_row( cpl_image *im, Vector *row )
00043 *   4) cpl_image * sinfo_new_mult_row_to_image( cpl_image *im, Vector *row )
00044 *   5) cpl_image * sinfo_new_col_tilt ( cpl_image * image, float sigmaFactor )
00045 *   6) cpl_image * sinfo_new_median_image( cpl_image * im, float fmedian )
00046 *   7) cpl_image * sinfo_new_compare_images( cpl_image * im1, 
00047                                              cpl_image * im2, 
00048                                              cpl_image * origim )
00049 *   8) cpl_image * sinfo_new_thresh_image ( cpl_image * im, 
00050                                             float lo_cut, float hi_cut )
00051 *   9) pixel_map * sinfo_new_promote_image_to_pixelmap ( cpl_image * im )
00052 *  10) cpl_image * sinfo_new_promote_image_to_mask ( cpl_image * im, 
00053                                                      int * n_badpixels )
00054 *  11) cpl_image * sinfo_new_mult_image_by_mask ( cpl_image * im, 
00055                                                   cpl_image * mask )
00056 *  12) cpl_image * sinfo_new_interpol_image ( cpl_image * im,
00057 *                                 cpl_image * mask,
00058 *                                 int        max_radius,
00059 *                                 int        n_pixels )
00060 *  13) cpl_image * sinfo_interpol_source_image ( cpl_image * im,
00061 *                                       cpl_image * mask,
00062 *                                       int        max_rad,
00063 *                                       float   ** slit_edges )
00064 *  14) cpl_image * sinfo_new_stack_row_to_image ( Vector * row, int ly )
00065 *  15) Stats * sinfo_new_image_stats_on_rectangle ( cpl_image * im, 
00066 *                                      float      loReject,
00067 *                                      float      hiReject,
00068 *                                      int        llx, 
00069 *                                      int        lly, 
00070 *                                      int        urx, 
00071 *                                      int        ury )
00072 *  16) cpl_image * sinfo_new_normalize_to_central_pixel ( cpl_image * image )
00073 *  17) cpl_image *
00074 *      sinfo_new_shift_image(
00075 *            cpl_image    *    image_in,
00076 *            double           shift_x,
00077 *            double           shift_y,
00078 *            double       *   interp_kernel)
00079 *  18) cpl_image * sinfo_new_combine_masks ( cpl_image * firstMask, 
00080                                              cpl_image * secondMask ) 
00081 *  19) cpl_image * sinfo_new_slice_cube (cpl_imagelist * cube, int x, int y ) 
00082 *  20) cpl_image * sinfo_new_div_images_robust ( cpl_image * im1, 
00083                                                  cpl_image * im2 ) 
00084 *
00085 *
00086 *   DESCRIPTION
00087 *   1) takes the average of each image column
00088 *   2) takes the average of each image column by sorting the
00089 *      column values and rejecting the given percentage of 
00090 *      the highest and lowest values  [0...1]
00091 *   3) divides each image column by a row value
00092 *   4) multiplies each image column with a row value
00093 *   5) first calculates statistics for each column of an image.
00094 *      sinfo_median value and standard deviation of columns are de-
00095 *      termined, blank values are excluded. Fits a straight 
00096 *      line through the pixel values of each column and subtracts
00097 *      the fit in order to remove the tilt of each column.
00098 *      Only those pixels are used for the fit that are within 
00099 *      a defined factor of sigma noise limit. The noise is 
00100 *      calculated from pixels between the 10percentil and 
00101 *      90percentil points. 
00102 *      if the straight line could not be determined, the sinfo_median
00103 *      of the column is subtracted from the column
00104 *   6) sinfo_median filter, calculates the sinfo_median for an image
00105 *      by using the 8 closest pixels to each pixel.
00106 *      The values in the output image are determined according
00107 *      to the values of the input parameter. 
00108 *      If fmedian = 0: always replace by sinfo_median
00109 *      if fmedian < 0: replace by sinfo_median if |pixel - sinfo_median| >
00110 *                      -fmedian
00111 *      if fmedian > 0: replace by sinfo_median (fmedian as a factor of 
00112 *                      the square root of the sinfo_median itself)          
00113 *                      if |pixel - median| >= fmedian * sqrt ( median ) 
00114 *                      This can be used to consider photon noise.
00115 *                      This considers a dependence of the differences on the 
00116 *                      pixel values themselves.
00117 *   7) if a pixel value of one image (im1) equals 
00118 *      the pixel value of the other image keep the 
00119 *      pixel value of the original image otherwise replace 
00120 *      it with ZEROs
00121 *   8) simple search for static bad pixels for a flat field
00122 *      or sinfo_dark frame, values below and above the threshold
00123 *      values are set to ZERO.
00124 *   9) changes an image with ZERO indicated bad pixels to
00125 *      a bad pixel map.
00126 *  10) changes an image with ZERO indicated bad pixels to
00127 *      a bad pixel mask image, that means the returned
00128 *      image has values 1 at positions of good pixels and 
00129 *      ZEROs at positions of bad pixels.
00130 *  11) changes an image to an image that has ZERO indicated 
00131 *      static bad pixels 
00132 *  12) interpolates all bad pixels indicated by the bad pixel mask. 
00133 *      Therefore, the mean of at least 2 valid values of 
00134 *      the nearest 8 neighbors is taken. If too much 
00135 *      neighbors are also bad pixels 
00136 *      the neighbor radius is increased to a maximum of 
00137 *      max_radius until n_pixels valid pixels are found.
00138 *      The valid neighbors are searched by going through
00139 *      the columns and rows around the central square that
00140 *      was already searched.
00141 *      The bad pixel is interpolated by the mean of these 
00142 *      valid pixels (less than 9) or by the sinfo_median of them
00143 *      (more than 8).
00144 *  13) interpolates all bad pixels indicated by the bad pixel mask. 
00145 *      Therefore, the mean of the nearest 4 neighbors is taken, 
00146 *      two in spectral direction and 2 in spatial direction.
00147 *      The routine cares about the image and slitlet edges.
00148 *      If there are no good pixel found within the nearest neighbors,
00149 *      the next 4 nearest neighbors in spatial and spectral direction
00150 *      are searched for valid pixels until a limit of max_rad.   
00151 *      A maximum of 4 valid pixels are used for interpolation by their mean.
00152 *  14) stack a given image row to build a whole image
00153 *  15) computes the mean and standard deviation of 
00154 *      a given rectangle on an image by leaving the extreme
00155 *      intensity values.
00156 *  16) normalizes a raw flatfield image by dividing by the median of the 
00157        central spectral pixels to produce a master flatfield
00158 *  17) This function is a conversion to CPL of the ECLIPSE function 
00159        shift_image()
00160 *      but slightly changed. If a blank (ZERO) pixel appears the blank pixel
00161 *      is shifted but preserved as blank.
00162 *      If a blank (ZERO) pixel appears within the
00163 *      interpolation kernel the blank pixel is set to 0. 
00164 * 
00165 *      This function shifts an image by a non-integer offset, using
00166 *      interpolation. You can either generate an interpolation kernel once and
00167 *      pass it to this function, or let it generate a default kernel. In the
00168 *      former case, use sinfo_generate_interpolation_kernel() to generate an
00169 *      appropriate kernel. In the latter case, pass NULL as last argument. A
00170 *      default interpolation kernel is then generated then discarded 
00171        before this function returns.
00172 *
00173 *      The returned image is a newly allocated object, it must be deallocated
00174 *      using cpl_image_delete().
00175 *  18) combines two bad pixel mask to one using an or relation
00176 *  19) slices a data cube in x or y direction
00177 *  20) divides two images by considering blanks and
00178 *      calculating first 1/im2 by
00179 *      cutting the very high values and setting to 1,
00180 *      then multiplying im1 * 1/im2.  
00181 *
00182 *   FILES
00183 *
00184 *   ENVIRONMENT
00185 *
00186 *   RETURN VALUES
00187 *
00188 *   CAUTIONS
00189 *
00190 *   EXAMPLES
00191 *
00192 *   SEE ALSO
00193 *
00194 *   BUGS
00195 *
00196 *------------------------------------------------------------------------
00197 */
00198 
00199 #ifdef HAVE_CONFIG_H
00200 #  include <config.h>
00201 #endif
00202 #include "sinfo_vltPort.h"
00203 
00204 /*
00205  * System Headers
00206  */
00207 
00208 /*
00209  * Local Headers
00210  */
00211 
00212 #include "sinfo_image_ops.h"
00213 #include "sinfo_resampling.h"
00214 #include "sinfo_local_types.h"
00223 /*----------------------------------------------------------------------------
00224  *                          Function codes
00225  *--------------------------------------------------------------------------*/
00226 
00235 double sinfo_new_my_median_image(cpl_image* im)
00236 {
00237   double m=0;
00238   register int i=0;
00239   int n=0;
00240   pixelvalue* pv=0;
00241   int ilx=0;
00242   int ily=0;
00243   float* pidata=NULL;
00244 
00245 
00246   if(im==NULL) sinfo_msg_error("Null Image");
00247   ilx=cpl_image_get_size_x(im);
00248   ily=cpl_image_get_size_y(im);
00249   pidata=cpl_image_get_data_float(im);
00250   
00251    for ( i = 0 ; i < (int) ilx*ily ; i++ )
00252     {
00253       if ( isnan(pidata[i]) )
00254         {
00255 
00256     } else {
00257       n++;
00258     }
00259     }
00260    pv = cpl_calloc(n,sizeof(pixelvalue));
00261    n=0;
00262    for ( i = 0 ; i < (int) ilx*ily ; i++ )
00263     {
00264       if ( isnan(pidata[i]) )
00265         {
00266 
00267     } else {
00268       pv[n]=pidata[i];
00269           n++;
00270     }
00271     }
00272    if(pv == NULL || n == 0) {
00273      m=0;
00274    } else {
00275      m=sinfo_new_median(pv,n);
00276    }
00277    cpl_free(pv);
00278    if(isnan(m)){
00279      m=0;
00280    }
00281   return m;
00282 }
00283 
00292 Vector * sinfo_new_mean_of_columns( cpl_image *im )
00293 {
00294     Vector * row=NULL ;
00295     int i=0;
00296     int j=0;
00297     int ilx=0;
00298     int ily=0;
00299     float* pidata=NULL;
00300 
00301     if ( im == NULL )
00302     {
00303         sinfo_msg_error ("null image") ;
00304         return NullVector ;
00305     }
00306     ilx=cpl_image_get_size_x(im);
00307     ily=cpl_image_get_size_y(im);
00308     pidata=cpl_image_get_data_float(im);
00309     
00310     /* allocate memory for a row with the length of the image x-axis */
00311     if ( NULL == (row = sinfo_new_vector (ilx)) )
00312     {
00313         sinfo_msg_error ("not able to allocate a sinfo_vector" ) ;
00314         return NullVector ;
00315     }
00316 
00317     for ( i = 0 ; i < ilx ; i++ )
00318     {
00319         for ( j = 0 ; j < ily ; j++ )
00320         {
00321             if (!isnan(pidata[i+j*ilx]))
00322             {
00323                 row->data[i] += pidata[i + j*(ilx)] ;
00324             }
00325         }
00326 
00327         row->data[i] /= ily ;
00328     }
00329     return row ;
00330 }
00342 cpl_image * sinfo_new_clean_mean_of_columns( cpl_image *im, 
00343                              float lo_reject, 
00344                              float hi_reject)
00345 {
00346     cpl_image     * row=NULL ;
00347     pixelvalue*   buffer=NULL ;
00348     int          i=0;
00349     int          j=0;
00350     int          k=0;
00351     int          nv=0;
00352     int          lo_n=0;
00353     int          hi_n=0;
00354     int ilx=0;
00355     int ily=0;
00356     float* pidata=NULL;
00357     float* podata=NULL;
00358 
00359     if ( im == NULL )
00360     {
00361         sinfo_msg_error ("null image") ;
00362         return NULL ;
00363     }
00364     ilx=cpl_image_get_size_x(im);
00365     ily=cpl_image_get_size_y(im);
00366     pidata=cpl_image_get_data_float(im);
00367 
00368     if ((lo_reject + hi_reject) > 0.9)
00369     {
00370         sinfo_msg_error("illegal rejection thresholds: [%f] and [%f]", 
00371                         lo_reject, hi_reject) ;
00372         sinfo_msg_error("threshold sum should not be over "
00373                         "0.90 aborting average") ;
00374         return NULL ;
00375     }
00376 
00377     lo_n = (int) (ily * lo_reject + 0.5) ;
00378     hi_n = (int) (ily * hi_reject + 0.5) ;
00379     if (lo_n + hi_n >= ily)
00380     {
00381         sinfo_msg_error ("everything would be rejected") ;
00382         return NULL ;
00383     }
00384 
00385     /* allocate memory for a row with the length of the image x-axis */
00386     if ( NULL == (row = cpl_image_new (ilx, 1,CPL_TYPE_FLOAT)) ) 
00387     {
00388         sinfo_msg_error ("cannot allocate new image") ;
00389         return NULL ;
00390     }
00391     podata=cpl_image_get_data_float(row);
00392 
00393     buffer=(pixelvalue*) cpl_calloc(ily,sizeof(pixelvalue)) ;
00394 
00395     for ( i = 0 ; i < ilx ; i++ )
00396     {
00397         for ( j = 0 ; j < ily ; j++ )
00398         {
00399             buffer[j] = pidata[i + j*(ilx)] ;
00400         }
00401         sinfo_pixel_qsort (buffer, ily) ;
00402         
00403         nv = 0 ;
00404         for (k = lo_n ; k < ily - hi_n ; k ++)
00405         {
00406             if ( !isnan(buffer[k]) )
00407             {
00408                 podata[i] += buffer[k] ;
00409                 nv ++ ;
00410             }
00411         }
00412         podata[i] /= nv ;
00413  
00414     }
00415     cpl_free(buffer);
00416     return row ;
00417 }
00418 
00419 
00429 cpl_image * sinfo_new_div_image_by_row( cpl_image *im, Vector *row )
00430 {
00431     cpl_image *image=NULL ;
00432     int         i=0;
00433     int         j=0;
00434     int ilx=0;
00435     int ily=0;
00436     float* pidata=NULL;
00437     float* podata=NULL;
00438 
00439     if ( im == NULL || row == NULL )
00440     {
00441         sinfo_msg_error ("null image or null row") ;
00442         return NULL ;
00443     }
00444     ilx=cpl_image_get_size_x(im);
00445     ily=cpl_image_get_size_y(im);
00446     pidata=cpl_image_get_data_float(im);
00447 
00448     if ( ilx != row -> n_elements )
00449     {
00450         sinfo_msg_error("image and row size not compatible") ;
00451         return NULL ;
00452     }
00453     
00454     if ( NULL == (image = cpl_image_duplicate (im)) ) 
00455     {
00456         sinfo_msg_error ("cannot copy image") ;
00457         return NULL ;
00458     }
00459     podata=cpl_image_get_data_float(image);
00460     
00461     for (i = 0 ; i < ilx ; i++ )
00462     {
00463         for (j = 0 ; j < ily ; j++)
00464         {
00465             if ( !isnan(pidata[i + j*(ilx)]) )
00466             {
00467                 podata[i + j*(ilx)] = pidata[i + j*(ilx)] / row -> data[i] ;
00468             }
00469         }
00470     }
00471     return image ;
00472 }
00473 
00474 
00484 cpl_image * sinfo_new_mult_row_to_image( cpl_image *im, Vector *row )
00485 {
00486     cpl_image *image=NULL;
00487     int         i=0;
00488     int         j=0;
00489     int ilx=0;
00490     int ily=0;
00491     float* pidata=NULL;
00492     float* podata=NULL;
00493 
00494 
00495 
00496 
00497     if ( im == NULL || row == NULL )
00498     {
00499         sinfo_msg_error ("null image or null row") ;
00500         return NULL ;
00501     }
00502     ilx=cpl_image_get_size_x(im);
00503     ily=cpl_image_get_size_y(im);
00504     pidata=cpl_image_get_data_float(im);
00505 
00506     if ( ilx != row -> n_elements )
00507     {
00508         sinfo_msg_error("image and row size not compatible") ;
00509         return NULL ;
00510     }
00511 
00512     if ( NULL == (image = cpl_image_duplicate (im)) ) 
00513     {
00514         sinfo_msg_error ("cannot copy image") ;
00515         return NULL ;
00516     }
00517     podata=cpl_image_get_data_float(image);
00518     
00519     for (i = 0 ; i < ilx ; i++ )
00520     {
00521         for (j = 0 ; j < ily ; j++)
00522         {
00523             if ( !isnan(pidata[i + j*(ilx)]) )
00524             {
00525                 podata[i + j*(ilx)] = pidata[i + j*(ilx)] * row -> data[i] ;
00526             }
00527         }
00528     }
00529     return image ;
00530 }
00531 
00532 
00533 
00534 
00535 
00536 
00560 cpl_image * sinfo_new_col_tilt ( cpl_image * image, float sigmaFactor )
00561 {
00562     cpl_image   * im=NULL;
00563     float      * column=NULL ;
00564     double       sum=0;
00565     double  sum2=0;
00566     double  mean=0;
00567     float   sinfo_median=0;
00568     float   noise=0 ;
00569     float      * sig=NULL;
00570     float    * dat=NULL; 
00571     float        a=0;
00572     float        b=0;
00573     float        siga=0;
00574     float        sigb=0;
00575     float        chi2=0;
00576     float        q=0;
00577     int          i=0;
00578     int          j=0;
00579     int          colnum=0;
00580     int         npix=0;
00581     int         mwt=0 ;
00582     int lx=0;
00583     int ly=0;
00584     float* p_in_data=NULL;
00585     float* p_ou_data=NULL;
00586 
00587 
00588     if ( image == NULL )
00589     {
00590         sinfo_msg_error ("no image given" ) ;     
00591         return NULL ;
00592     }
00593     
00594     if ( sigmaFactor <= 0. )
00595     {
00596         sinfo_msg_error ("no or negative sigma factor") ;
00597         return NULL ;
00598     }
00599     lx = cpl_image_get_size_x(image);
00600     ly = cpl_image_get_size_y(image);
00601 
00602     
00603     /* allocate memory */
00604     if ( NULL == (im = cpl_image_new (lx,ly,CPL_TYPE_FLOAT )) ) 
00605     {
00606         sinfo_msg_error ("cannot allocate new image" ) ;
00607         return NULL ;
00608     }
00609 
00610     /* go through the columns */
00611     p_in_data = cpl_image_get_data_float(image);
00612     p_ou_data = cpl_image_get_data_float(im);
00613     for ( i = 0 ; i < lx ; i ++ )
00614     {
00615         /* initialize the buffer variables for each column */
00616         colnum = 0 ;
00617         column = (float *) cpl_calloc ( ly , sizeof (float *) ) ;
00618         sig    = (float *) cpl_calloc ( ly , sizeof (float *) ) ;
00619         dat    = (float *) cpl_calloc ( ly , sizeof (float *) ) ;
00620 
00621         /*select only non-ZERO values of one column*/
00622         for ( j = 0 ; j < ly ; j++ )    
00623         {
00624             if ( !isnan(p_in_data[i + j*lx]) )
00625             {
00626                 column[j] = p_in_data[i + j*lx] ;
00627                 colnum ++ ;
00628             }
00629         }
00630         if ( colnum < 10 ) 
00631         {
00632             /*sinfo_msg_warning ("sinfo_new_col_tilt:",
00633           "column %d has almost only blank pixels and is set to blank", i+1) ;*/
00634             for ( j = 0 ; j < ly ; j++ )    
00635             {
00636                 p_ou_data[i + j*lx] = ZERO;
00637             }
00638             /*
00639             cpl_free (column) ;
00640             cpl_free (sig);
00641             cpl_free (dat) ;
00642             continue ; 
00643             */
00644         }
00645 
00646         /*-------------------------------------------------------------------
00647          * sort the data, clip off the extremes, determine the noise
00648          * and get the range for the valid data. It is assumed here
00649          * that most pixels are o.k. 
00650          */
00651 
00652         sinfo_pixel_qsort (column, colnum) ;
00653         
00654         sum   = 0. ;
00655         sum2  = 0. ;
00656         npix  = 0  ;
00657 
00658         for ( j = 0.1*colnum + 1 ; j <= 0.9*colnum ; j++ )
00659         {
00660             sum  += column[j] ;
00661             sum2 += column[j] * column[j] ;
00662             npix ++ ;
00663         }
00664 
00665         if (npix <= 1)
00666         {
00667             noise = sigmaFactor * 1000.;
00668         }
00669         else
00670         {
00671             mean   = sum/(float)npix ;
00672             noise  = sqrt( (sum2 - sum*mean)/(double)(npix -1) ) ;
00673             noise *= sigmaFactor ;
00674         }
00675 
00676         /* -------------------------------------------------------------
00677          * determine sinfo_median if colnum is odd, sinfo_median will be the 
00678            colnum/2 th value, otherwise
00679          * sinfo_median is the mean of colnum/2-1 th and colnum/2 th value.
00680          */       
00681 
00682         if ( colnum % 2 == 1 )
00683         {
00684             sinfo_median = column[colnum/2] ;
00685         }
00686         else
00687         {
00688             sinfo_median = (column[colnum/2 - 1] + column[colnum/2])/2. ;
00689         }
00690 
00691         /* now select the pixels for the tilt calculation */
00692         
00693         colnum = 0 ;
00694         for ( j = 0; j < ly ; j ++ )
00695         {
00696             if ( !isnan(p_in_data[i+j*lx]) &&
00697                  fabs ( (p_in_data[i+j*lx]) - sinfo_median) <= noise )
00698             {
00699                 column[colnum] = p_in_data[i+j*lx] ;
00700                 dat[colnum] = (float) j ;
00701                 sig[colnum] = 1. ;
00702                 colnum ++ ;
00703             }
00704         }
00705 
00706         if ( colnum == 0 )
00707         {
00708             /*for ( j = 0; j < ly; j++ )
00709             {
00710                 p_ou_data[i+j*lx] -= sinfo_median ;
00711             }
00712             cpl_free (column) ;
00713             cpl_free (sig)    ;
00714             cpl_free (dat)    ;
00715             continue ;*/
00716         a=0./0.;
00717         b=0./0.;
00718         }
00719     else
00720     {
00721         mwt = 0 ;
00722         sinfo_my_fit ( dat, column, colnum, sig, mwt, &a, 
00723                        &b, &siga, &sigb, &chi2, &q ) ;
00724     }
00725         if ( fabs(b) >= SLOPE || fabs(a) >= SATURATION  || 
00726              isnan(b) || isnan(a))
00727         {
00728             sinfo_msg_warning ("linear fit: slope is greater than limit: %f"
00729                                " saturation level is reached: %f in column"
00730                                " number %d ", b, a , i+1) ; 
00731         }
00732         
00733         /* subtract fit or sinfo_median from data */
00734         for ( j = 0; j < ly; j++ )
00735         {
00736             if ( !isnan(p_in_data[i+j*lx]) &&
00737                  fabs(b) < SLOPE && fabs(a) < SATURATION )
00738             {
00739                 p_ou_data[i+j*lx] = p_in_data[i+j*lx] - (a + b * (float)j) ;
00740             }
00741             else if ( isnan(p_in_data[i+j*lx]) )
00742             {
00743                 p_ou_data[i+j*lx] = ZERO ;
00744             }
00745             else if ( (fabs(b) >= SLOPE || 
00746                        fabs(a) >= SATURATION || isnan(a) || isnan(b)) && 
00747                       !isnan(p_in_data[i+j*lx]) )
00748             {
00749                 p_ou_data[i+j*lx] -= sinfo_median ;
00750             }
00751             else
00752             {
00753                 sinfo_msg_error (" case is not possible! %f %f", b,a) ;
00754                 /*cpl_free (column) ;
00755                 cpl_free (sig)    ;
00756                 cpl_free (dat)    ;
00757                 cpl_image_delete(im) ;
00758                 return NULL ;*/
00759             }
00760         }
00761         cpl_free (column) ;
00762         cpl_free (sig)    ;
00763         cpl_free (dat)    ;
00764     }
00765 
00766     return im     ;
00767 }
00768 
00769 
00770 
00771 
00772 
00796 cpl_image * sinfo_new_median_image( cpl_image * im, float fmedian )
00797 {
00798     cpl_image *   image=NULL       ;
00799     pixelvalue * value=NULL       ;
00800     pixelvalue   sinfo_median=0      ;
00801     int        * position=NULL    ;
00802     int          nposition=0   ;
00803     int          n=0;
00804     int          i=0;
00805     int          j=0;
00806     int lx=0;
00807     int ly=0;
00808     float* p_in_data=NULL;
00809     float* p_ou_data=NULL;
00810 
00811     if ( im == NULL )
00812     {
00813         sinfo_msg_error ("no image input") ;
00814         return NULL ;
00815     }
00816     
00817     image = cpl_image_duplicate ( im ) ;
00818     lx=cpl_image_get_size_x(im);
00819     ly=cpl_image_get_size_y(im);
00820     p_in_data=cpl_image_get_data_float(im);
00821     p_ou_data=cpl_image_get_data_float(image);
00822 
00823     /*---------------------------------------------------------------------- 
00824      * go through all pixels 
00825      */
00826 
00827     for ( i = 0 ; i < (int) lx*ly ; i++ )
00828     {
00829         /* blank pixels are not replaced */ 
00830         if ( isnan(p_in_data[i]) )
00831         {
00832             continue ;
00833         }
00834 
00835         /* initialize the buffer variables for the 8 nearest neighbors */
00836         value = (pixelvalue * )cpl_calloc ( 8, sizeof ( pixelvalue * ) ) ;
00837         position = ( int * ) cpl_calloc ( 8, sizeof ( int * ) ) ;
00838 
00839         /*-------------------------------------------------------------------- 
00840          * determine the pixel position of the 8 nearest neighbors  
00841          */
00842 
00843         position[0] = i + lx - 1 ; /* upper left  */
00844         position[1] = i + lx     ; /* upper       */
00845         position[2] = i + lx + 1 ; /* upper right */
00846         position[3] = i + 1      ; /* right       */
00847         position[4] = i - lx + 1 ; /* lower right */
00848         position[5] = i - lx     ; /* lower       */
00849         position[6] = i - lx - 1 ; /* lower left  */
00850         position[7] = i - 1      ; /* left        */
00851        
00852         /*------------------------------------------------------------------- 
00853          * determine the positions of the image margins, top positions are 
00854            changed to low positions and vice versa. Right positions are 
00855            changed to left positions and vice versa.
00856          */
00857 
00858         if ( i >= 0 && i < lx )    /* bottom line */
00859         {
00860             position[4] += 2 * lx ;  
00861             position[5] += 2 * lx ;  
00862             position[6] += 2 * lx ;  
00863         }
00864         else if ( i >= ((int) lx*ly - lx ) && i < (int) lx*ly ) /* top line */
00865         {
00866             position[0] -= 2 * lx ;  
00867             position[1] -= 2 * lx ;  
00868             position[2] -= 2 * lx ;  
00869         }
00870         else if ( i % lx == 0 )    /* left side */
00871         {
00872             position[0] += 2 ;
00873             position[6] += 2 ;
00874             position[7] += 2 ;
00875         }
00876         else if ( i % lx == lx - 1 )    /* right side */
00877         {
00878             position[2] -= 2 ;
00879             position[3] -= 2 ;
00880             position[4] -= 2 ;
00881         }
00882 
00883         /* --------------------------------------------------------------------
00884          * read the pixel values of the neighboring pixels,
00885          * blanks are not considered 
00886          */
00887 
00888         nposition = 8 ;
00889         n = 0 ;
00890         for ( j = 0 ; j < nposition ; j ++ )
00891         {
00892             if ( !isnan(p_in_data[position[j]]) )
00893             {
00894                 value[n] = p_in_data[position[j]] ;    
00895                 n ++ ;
00896             }
00897         }
00898         nposition = n ;
00899 
00900         if ( nposition <= 1 )  /* almost all neighbors are blank */
00901         {
00902             p_ou_data[i] = ZERO ;
00903             cpl_free(value) ;
00904             cpl_free(position) ;
00905             continue ;
00906         }
00907         
00908         /* sort the values and determine the sinfo_median */
00909 
00910         sinfo_pixel_qsort ( value, nposition ) ;
00911         if ( nposition % 2 == 1 )
00912         {
00913             sinfo_median = value [ nposition/2 ] ;
00914         }
00915         else
00916         {
00917             sinfo_median = ( value [nposition/2 - 1] + 
00918                              value [nposition/2] ) / 2. ;
00919         }
00920 
00921         /* -----------------------------------------------------------------
00922          * replace the pixel value by the sinfo_median on conditions:
00923          * fmedian = 0: always replace with sinfo_median.
00924          * fmedian < 0: interpret as absolute condition:
00925          *              if |pixel - sinfo_median| > -fmedian
00926          *              replace with sinfo_median.
00927          * fmedian > 0: replace by sinfo_median (fmedian as a factor of 
00928          *              the square root of the sinfo_median itself)            
00929          *              if |pixel - sinfo_median| >= fmedian * 
00930                                                      sqrt ( sinfo_median ) 
00931          *              considers a dependence on the pixel value.
00932          *              This can be used to consider photon noise.
00933          */
00934 
00935         if ( fmedian == 0 )
00936         {
00937             p_ou_data[i] = sinfo_median ;
00938         }
00939         else if ( fmedian < 0 &&
00940                   fabs ( sinfo_median - p_in_data[i] ) >= -fmedian )
00941         {
00942             p_ou_data[i] = sinfo_median ;
00943         }
00944         else if ( fmedian > 0 &&
00945                   fabs ( sinfo_median - p_in_data[i] ) >= fmedian * 
00946                                                       sqrt(fabs(sinfo_median)) )
00947         {
00948             p_ou_data[i] = sinfo_median ;
00949         }   
00950         else
00951         {
00952             cpl_free (value) ; 
00953             cpl_free (position) ; 
00954             continue ;
00955         }
00956 
00957         cpl_free (value) ;
00958         cpl_free (position) ;
00959     }       
00960     return image ;
00961 }
00962 
00963 
00964 
00965 
00976 cpl_image * 
00977 sinfo_new_compare_images(cpl_image * im1,cpl_image * im2,cpl_image * origim ) 
00978 {
00979     cpl_image * image=NULL ;
00980     int            i=0 ;
00981     int lx1=0;
00982     int ly1=0;
00983     int lx2=0;
00984     int ly2=0;
00985     float* p_in1_data=NULL;
00986     float* p_in2_data=NULL;
00987     float* p_ou_data=NULL;
00988     float* p_org_data=NULL;
00989 
00990  
00991     if ( im1 == NULL || im2 == NULL || origim == NULL )
00992     {
00993         sinfo_msg_error ("Null images as input" ) ;
00994         return NULL ;
00995     }
00996     lx1=cpl_image_get_size_x(im1);
00997     ly1=cpl_image_get_size_y(im1);
00998 
00999     lx2=cpl_image_get_size_x(im2);
01000     ly2=cpl_image_get_size_y(im2);
01001 
01002     p_in1_data=cpl_image_get_data_float(im1);
01003     p_in2_data=cpl_image_get_data_float(im2);
01004     p_org_data=cpl_image_get_data_float(origim);
01005     if ( lx1 != lx2 || ly1 != ly2 )
01006     {
01007         sinfo_msg_error ("incompatible image sizes" ) ;
01008         return NULL ;
01009     }
01010     
01011     /* allocate memory */
01012     if ( NULL == (image = cpl_image_new ( lx1, ly1, CPL_TYPE_FLOAT )) ) 
01013     {
01014         sinfo_msg_error ("cannot allocate new image" ) ;
01015         return NULL ;
01016     }
01017     p_ou_data=cpl_image_get_data_float(image);
01018     for ( i = 0 ; i < (int) lx1*ly1 ; i ++ )
01019     {
01020         if ( isnan(p_in1_data[i]) && isnan(p_in2_data[i]) )
01021         {
01022             p_ou_data[i] = ZERO ;
01023         }
01024         else
01025         {
01026             if ( p_in1_data[i] == p_in2_data[i] )
01027             {
01028                 p_ou_data[i] = p_org_data[i] ;
01029             }
01030             else
01031             {
01032                 p_ou_data[i] = ZERO ;
01033             }
01034         }
01035     }
01036     return image ;
01037 }
01038 
01039     
01040 
01052 cpl_image * 
01053 sinfo_new_promote_image_to_mask (cpl_image * im, int * n_badpixels )
01054 {
01055     cpl_image * reImage=NULL ;
01056     int        i=0 ;
01057     int lx=0;
01058     int ly=0;
01059     float* p_in_data=NULL;
01060     float* p_ou_data=NULL;
01061 
01062     if ( NULL == im )
01063     {
01064         sinfo_msg_error("no input image given!") ;
01065         return NULL ;
01066     }
01067     lx=cpl_image_get_size_x(im);
01068     ly=cpl_image_get_size_y(im);
01069     p_in_data=cpl_image_get_data_float(im);
01070 
01071     /* allocate memory for the returned image */
01072     if ( NULL == (reImage = cpl_image_new (lx,ly,CPL_TYPE_FLOAT )) )
01073     {
01074         sinfo_msg_error ("cannot allocate new image!") ;
01075         return NULL ;
01076     } 
01077     p_ou_data=cpl_image_get_data_float(reImage);
01078     
01079     *n_badpixels = 0 ;
01080     for ( i = 0 ; i < (int) lx*ly ; i ++ )
01081     {
01082         if ( isnan(p_in_data[i]) )
01083         {
01084             p_ou_data[i] = 0. ;
01085             (*n_badpixels)++ ;
01086         }
01087         else
01088         {
01089             p_ou_data[i] = 1. ;
01090         }
01091     }
01092     return reImage ;
01093 }
01094 
01095 
01106 cpl_image * sinfo_new_mult_image_by_mask (cpl_image * im,cpl_image * mask )
01107 {
01108     cpl_image * reImage=NULL ;
01109     int        i=0 ;
01110     int ix=0;
01111     int iy=0;    
01112     int mx=0;
01113     int my=0;
01114 
01115 
01116     float* pmdata=NULL;
01117     float* podata=NULL;
01118 
01119     if ( NULL == im )
01120     {
01121         sinfo_msg_error("no input image given!") ;
01122         return NULL ;
01123     }
01124     if ( NULL == mask )
01125     {
01126         sinfo_msg_error("no mask image given!") ;
01127         return NULL ;
01128     }
01129     ix=cpl_image_get_size_x(im);
01130     iy=cpl_image_get_size_y(im);
01131     mx=cpl_image_get_size_x(mask);
01132     my=cpl_image_get_size_y(mask);
01133 
01134     if ( ix != mx || iy != my)
01135     {
01136         sinfo_msg_error("image sizes are not correspondent!") ;
01137         return NULL ;
01138     }
01139 
01140     reImage = cpl_image_duplicate( im ) ;
01141     podata=cpl_image_get_data_float(reImage);
01142     pmdata=cpl_image_get_data_float(mask);
01143 
01144     for ( i = 0 ; i < (int) ix*iy ; i ++ )
01145     {
01146         if ( pmdata[i] == 0. )
01147         {
01148             podata[i] = ZERO ; 
01149         }
01150     }
01151 
01152     return reImage ;
01153 }
01154 
01155 
01156 
01166 cpl_image * 
01167 sinfo_new_thresh_image (cpl_image * im, float lo_cut, float hi_cut )
01168 {
01169     cpl_image * image=NULL ;
01170     float* p_inp_data=NULL;
01171     float* p_out_data=NULL;
01172     int lx=0;
01173     int ly=0;
01174 
01175     int            i=0 ;
01176     
01177     if (im == NULL)
01178     {
01179         sinfo_msg_error ("null image given") ;
01180         return NULL ;
01181     }
01182     lx=cpl_image_get_size_x(im);
01183     ly=cpl_image_get_size_y(im);
01184 
01185     image = cpl_image_duplicate(im) ;
01186     p_inp_data=cpl_image_get_data(im);
01187     p_out_data=cpl_image_get_data(image);
01188     for ( i = 0 ; i < (int) lx*ly ; i ++ )
01189     {
01190         if ( p_inp_data[i] > (pixelvalue) hi_cut || 
01191              p_inp_data[i] < (pixelvalue) lo_cut )
01192         {
01193              p_out_data[i] = ZERO ;
01194         }
01195     }
01196     return image ;
01197 }
01198 
01199 
01200 
01201 
01226 cpl_image * sinfo_new_interpol_image ( cpl_image * im,
01227                            cpl_image * mask,
01228                            int        max_radius,
01229                            int        n_pixels )
01230 {
01231     cpl_image * returnImage=NULL ;
01232     float* neighbors=NULL ;
01233     float sum=0;
01234     float mean=0;
01235     int i=0;
01236     int j=0;
01237     int k=0;
01238     int row=0;
01239     int col=0;
01240     int n_valid=0;
01241     int agreed=0;
01242 
01243     int ilx=0;
01244     int ily=0;
01245     int mlx=0;
01246     int mly=0;
01247     float* pidata=NULL;
01248     float* podata=NULL;
01249     float* pmdata=NULL;
01250 
01251     if ( NULL == im )
01252     {
01253         sinfo_msg_error("sorry, no input image given!") ;
01254         return NULL ;
01255     }
01256     ilx=cpl_image_get_size_x(im);
01257     ily=cpl_image_get_size_y(im);
01258     pidata=cpl_image_get_data_float(im);
01259 
01260     if ( NULL == mask )
01261     {
01262         sinfo_msg_error("sorry, no mask image given!") ;
01263         return NULL ;
01264     }
01265 
01266     mlx=cpl_image_get_size_x(mask);
01267     mly=cpl_image_get_size_y(mask);
01268     pmdata=cpl_image_get_data_float(mask);
01269 
01270     if ( mlx != ilx || mly != mly )
01271     {
01272         sinfo_msg_error("images not compatible !") ;
01273         return NULL ;
01274     }
01275 
01276     if ( max_radius <= 0 )
01277     {
01278         sinfo_msg_error("wrong number of pixels for maximal "
01279                         "search radius given!") ;
01280         return NULL ;
01281     }
01282 
01283     if ( n_pixels <= 2 )
01284     {
01285         sinfo_msg_error("wrong number of pixels used for interpolation given!") ;
01286         return NULL ;
01287     }
01288 
01289     returnImage = cpl_image_duplicate ( im ) ;
01290     podata=cpl_image_get_data_float(returnImage);
01291 
01292     /* go through the columns and rows of the input and mask image */
01293     neighbors=cpl_calloc(4*max_radius*max_radius,sizeof(float)) ;
01294 
01295     for ( col = 0 ; col < ilx ; col++ )
01296     {
01297         for ( row = 0 ; row < ily ; row++ )
01298         {
01299             /* look for the ZEROS that means the detected bad pixels */
01300             if ( isnan(pmdata[col+row*ilx]) || pmdata[col+row*ilx] == 0. )
01301             {
01302                 /* now the neighbors must be considered */
01303                 n_valid = 0 ;
01304                 agreed  = 0 ;
01305                 for ( j = 1 ; j <= max_radius ; j++ )
01306                 {
01307 
01308                     /* go through the left column */
01309                     for ( k = -j ; k < j ; k++ )
01310                     {
01311                         if ( col-j >= 0 && row+k < ily && row+k >= 0 )
01312                         {
01313                             if ( !isnan(pmdata[col-j+(row+k)*mlx]) || 
01314                                  pmdata[col-j+(row+k)*mlx] != 0 )
01315                             {
01316                                 neighbors[n_valid]=pidata[col-j+(row+k)*ilx] ;
01317                                 n_valid++ ;
01318                             }
01319                         }
01320                     }
01321 
01322                     /* go through the upper row */
01323                     for ( k = -j ; k < j ; k++ )
01324                     {
01325                         if ( col+k < ilx && col+k >= 0 && row+j < ily )
01326                         {
01327                             if ( !isnan(pmdata[col+k+(row+j)*mlx]) || 
01328                                  pmdata[col+k+(row+j)*mlx] != 0. )
01329                             {
01330                                 neighbors[n_valid]=pidata[col+k+(row+j)*ilx] ;
01331                                 n_valid++ ;
01332                             }
01333                         }
01334                     }
01335 
01336                     /* go through the right column */
01337                     for ( k = -j ; k < j ; k++ )
01338                     {
01339                         if ( col+j < ilx  && row-k >= 0 && row-k < ily )
01340                         {
01341                             if ( !isnan(pmdata[col+j+(row-k)*mlx]) || 
01342                                  pmdata[col+j+(row-k)*mlx] != 0. )
01343                             {
01344                                 neighbors[n_valid]=pidata[col+j+(row-k)*ilx] ;
01345                                 n_valid++ ;
01346                             }
01347                         }
01348                     }
01349 
01350                     /* go through the lower row */
01351                     for ( k = -j ; k < j ; k++ )
01352                     {
01353                         if ( col-k >= 0 && col-k < ilx && row-j < ily )
01354                         {
01355                             if ( !isnan(pmdata[col-k+(row-j)*mlx]) ||
01356                                  pmdata[col-k+(row-j)*mlx] != 0. )
01357                             {
01358                                 neighbors[n_valid]=pidata[col-k+(row-j)*ilx] ;
01359                                 n_valid++ ;
01360                             }
01361                         }
01362                     }
01363 
01364                     /* control if the breaking criteria is fullfilled */
01365                     if ( n_valid >= n_pixels )
01366                     {
01367                         agreed = 1 ;
01368                         break ;
01369                     }
01370                     /* do a break if more than 2 nearest neighbors are found */
01371                     if ( j == 1 && n_valid >= 2 )
01372                     {
01373                         agreed = 1 ;
01374                         break ;
01375                     }
01376                 }
01377                 if ( n_valid < n_pixels && agreed == 0 )
01378                 {
01379                     sinfo_msg_error("not enough valid neighbors found to "
01380                                     "interpolate bad pixel in col: "
01381                                     "%d, row: %d", col, row ) ;
01382                     return NULL ;
01383                 }
01384                 else
01385                 {
01386                     /* ------------------------------------------------------ 
01387                      * take the mean of the valid neighboring pixels if less 
01388                      * than 9 valid pixels are available else take the 
01389                        sinfo_median.
01390                      */
01391                     if ( n_valid <= 8 )
01392                     {
01393                         sum = 0. ;
01394                          
01395                         for ( i = 0 ; i < n_valid ; i++ )
01396                         {
01397                             sum += neighbors[i] ;
01398                         }
01399                         mean = sum / n_valid ;
01400 
01401                         podata[col+row*ilx] = mean ;
01402                     }
01403                     else
01404                     {
01405                        podata[col+row*ilx]=sinfo_new_median(neighbors,n_valid);
01406                     }
01407                 }
01408             }
01409         }
01410     }
01411     cpl_free(neighbors);
01412     return returnImage ;
01413 }
01414 
01415 
01438 cpl_image * sinfo_interpol_source_image ( cpl_image * im,
01439                                  cpl_image * mask,
01440                                  int        max_rad,
01441                                  float   ** slit_edges )
01442 {
01443     cpl_image * returnImage=NULL ;
01444     float validpixel[6] ;
01445     float sum=0 ;
01446     int n=0;
01447     int row=0;
01448     int col=0;
01449     int i=0;
01450     int k=0;
01451     int slitlet=0;
01452     int n_slitlets=0;
01453     int ilx=0;
01454     int ily=0;
01455     int mlx=0;
01456     int mly=0;
01457   
01458     float* pidata=NULL;
01459     float* podata=NULL;
01460     float* pmdata=NULL;
01461 
01462 
01463     if ( NULL == im )
01464     {
01465         sinfo_msg_error("sorry, no input image given!") ;
01466         return NULL ;
01467     }
01468     ilx=cpl_image_get_size_x(im);
01469     ily=cpl_image_get_size_y(im);
01470     pidata=cpl_image_get_data_float(im);
01471 
01472     if ( NULL == mask )
01473     {
01474         sinfo_msg_error("sorry, no bad pixel mask image given!") ;
01475         return NULL ;
01476     }
01477     mlx=cpl_image_get_size_x(mask);
01478     mly=cpl_image_get_size_y(mask);
01479     pmdata=cpl_image_get_data_float(mask);
01480 
01481     if ( mlx != ilx || mly != ily )
01482     {
01483         sinfo_msg_error("images not compatible in size!") ;
01484         return NULL ;
01485     }
01486 
01487     if ( max_rad < 1 )
01488     {
01489         sinfo_msg_error("sorry, wrong maximum distance given!") ;
01490         return NULL ;
01491     }
01492 
01493     if ( slit_edges == NULL )
01494     {
01495         sinfo_msg_error("sorry, array slit_edges is empty!") ;
01496         return NULL ;
01497     }
01498  
01499     /* determine the number of slitlets */
01500     n_slitlets = N_SLITLETS ;
01501  
01502     /* copy the original image in the image that will be returned */
01503     returnImage = cpl_image_duplicate( im ) ;
01504     podata=cpl_image_get_data_float(returnImage);
01505     
01506     /* go through the rows and columns of the image and search for 
01507       the bad pixels */
01508     for ( row = 0 ; row < ily ; row++ )
01509     {
01510         for ( col = 0 ; col < ilx ; col++ )
01511         {
01512             n = 0 ;
01513             if ( isnan(pmdata[col + row*mlx]) || 
01514                  pmdata[col + row*mlx] == 0. ||
01515                  isnan(pidata[col + row*mlx]) )
01516             {
01517                 /* look for the slitlet where the bad pixel is found */
01518                 slitlet = -1000 ;
01519                 for ( k = 0 ; k < n_slitlets ; k++ )
01520                 {
01521                     if ( sinfo_new_nint(slit_edges[k][0]) <= col && 
01522                          sinfo_new_nint(slit_edges[k][1]) >= col )
01523                     {
01524                         slitlet = k ;
01525                     }
01526 /* The following else statement is wrong, because in the 
01527      end slitlet will always be -1000
01528             else
01529                     {
01530                         slitlet = -1000 ;
01531                     } 
01532 */
01533                 }      
01534                 for ( i = 0 ; i < 6 ; i++ )
01535                 {
01536                     validpixel[i] = 0. ;
01537                 }
01538                 /* look for the valid nearest neighbors 
01539                    and collect them but only a maximum of 4 */
01540                 for ( i = 1 ; i <= max_rad ; i++ )
01541                 {
01542                     if ( row + i < ily)
01543                     {
01544                         if ( !isnan(pmdata[col + (row+i) * mlx])
01545                              && pmdata[col + (row+i) * mlx] != 0. &&
01546                                 !isnan(pidata[col + (row+i) * ilx]) )
01547                         {
01548                             validpixel[n] = pidata[col + (row+i) * ilx] ;
01549                             n++ ;
01550                         }
01551                     }
01552                     if ( row - i >= 0 )
01553                     {
01554                         if ( !isnan(pmdata[col + (row-i) * mlx]) 
01555                              && pmdata[col + (row-i) * mlx] != 0. &&
01556                                 !isnan(pidata[col + (row-i) * ilx]) )
01557                         {
01558                             validpixel[n] = pidata[col + (row-i) * ilx] ;
01559                             n++ ;
01560                         }
01561                     }
01562 
01563                     /* be aware of the slitlet edges in the 
01564                        spatial direction */
01565                     if ( col + i < ilx )  
01566                     {
01567                       if (  slitlet != -1000 )
01568                       {
01569                          if (col+i <= sinfo_new_nint(slit_edges[slitlet][1]) &&
01570                              !isnan(pmdata[col + i + row * mlx]) && 
01571                              pmdata[col + i + row * mlx] != 0. &&
01572                              !isnan(pidata[col + i + row * ilx]) )
01573                          {
01574                              validpixel[n] = pidata[col + i + row * ilx] ;
01575                              n++ ;
01576                          }
01577                       }
01578                     }
01579                     if ( col - i >= 0 )  
01580                     {
01581                       if ( slitlet != -1000 )
01582                       {
01583                          if (col-i >= sinfo_new_nint(slit_edges[slitlet][0]) &&
01584                              !isnan(pmdata[col - i + row * mlx]) &&
01585                              pmdata[col - i + row * mlx] != 0. && 
01586                              !isnan(pidata[col - i + row * ilx]) )  
01587                          {
01588                              validpixel[n] = pidata[col - i + row * ilx] ;
01589                              n++ ;
01590                          }
01591                       }
01592                     }
01593 
01594                     if ( i == 1 && n > 1 )
01595                     {
01596                         break ;
01597                     }
01598                     if ( n > 2 )
01599                     {
01600                         break ;
01601                     }
01602                 }
01603 
01604                 if ( n == 0 )
01605                 {
01606                     podata[col + row*ilx] = ZERO ; 
01607             /*sinfo_msg_warning("sinfo_interpolSourceImage:",
01608                                         "bad pixel in column: %d and row: %d"
01609                                         " could not be interpolated!",col,row);
01610                              */
01611                 }
01612                 else
01613                 {
01614                     /* now compute the mean and replace 
01615                        the bad pixel value by the mean */ 
01616                     sum = 0. ;
01617                     for ( i = 0 ; i < n ; i++ )
01618                     {
01619                         sum += validpixel[i] ;
01620                     }
01621                     podata[col + row*ilx] = sum/n ; 
01622                 }
01623             }
01624         }
01625     }
01626 
01627     return returnImage ;
01628 }
01629 
01639 cpl_image * sinfo_new_stack_row_to_image ( Vector * row, int ly )
01640 {
01641     cpl_image * image=NULL;
01642     int        col=0;
01643     int        ro=0;
01644     float* podata=NULL;
01645 
01646     if ( row == NullVector )
01647     {
01648         sinfo_msg_error ("Null sinfo_vector as input" ) ;
01649         return NULL ;
01650     }
01651     if ( ly <= 1 )
01652     {
01653         sinfo_msg_error ("wrong image length given" ) ;
01654         return NULL ;
01655     }
01656 
01657     /* allocate memory */
01658     if (NULL == (image = cpl_image_new(row->n_elements ,ly,CPL_TYPE_FLOAT )) ) 
01659     {
01660         sinfo_msg_error ("cannot allocate new image" ) ;
01661         return NULL ;
01662     }
01663     podata=cpl_image_get_data_float(image);
01664 
01665     for ( col = 0 ; col < row -> n_elements ; col++ )
01666     {
01667         for ( ro = 0 ; ro < ly ; ro++ )
01668         {
01669             podata[col + ro*ly] = row -> data[col] ;
01670         }
01671     }
01672     return image ;
01673 }
01674 
01690 Stats * sinfo_new_image_stats_on_rectangle ( cpl_image * im, 
01691                                 float      loReject,
01692                                 float      hiReject,
01693                                 int        llx, 
01694                                 int        lly, 
01695                                 int        urx, 
01696                                 int        ury )
01697 {
01698     Stats * retstats=NULL;
01699     int i=0 ;
01700     int row=0;
01701     int col=0;
01702     int n=0;
01703     int npix=0;
01704     int lo_n=0;
01705     int hi_n=0;
01706     double pix_sum=0;
01707     double sqr_sum=0;
01708     float * pix_array=NULL;
01709     int im_lx=0;
01710     int im_ly=0;
01711     float* pim=NULL;
01712 
01713     if ( NULL == im )
01714     {
01715         sinfo_msg_error("sorry, no input image given!") ;
01716         return NULL ;
01717     }
01718     if ( loReject+hiReject >= 100. )
01719     {
01720         sinfo_msg_error("sorry, too much pixels rejected!") ;
01721         return NULL ;
01722     }
01723     if ( loReject < 0. || loReject >= 100. || 
01724          hiReject < 0. || hiReject >= 100. )
01725     {
01726         sinfo_msg_error("sorry, negative reject values!") ;
01727         return NULL ;
01728     }
01729 
01730     im_lx=cpl_image_get_size_x(im);
01731     im_ly=cpl_image_get_size_x(im);
01732 
01733     if ( llx < 0 || lly < 0 || urx < 0 || ury < 0 ||
01734          llx >= im_lx || lly >= im_ly || urx >= im_lx || 
01735          ury >= im_ly || ury <= lly || urx <= llx )
01736     {
01737         sinfo_msg_error("sorry, wrong pixel coordinates of rectangle!") ;
01738         return NULL ;
01739     }
01740 
01741      /* allocate memory */
01742     retstats = (Stats*) cpl_calloc(1, sizeof(Stats)) ;
01743     npix = (urx - llx + 1) * (ury - lly + 1) ;
01744     pix_array = (float*) cpl_calloc ( npix, sizeof(float) ) ;
01745 
01746     /*------------------------------------------------------------------------- 
01747      * go through the rectangle and copy the pixel values into an array.
01748      */
01749     n = 0 ;
01750     pim = cpl_image_get_data_float(im);
01751     for ( row = lly ; row <= ury ; row++ )
01752     {
01753         for ( col = llx ; col <= urx ; col++ )
01754         {
01755             if ( !isnan(pim[col + row*im_lx]) )
01756             {
01757                 pix_array[n] = pim[col + row*im_lx] ;
01758                 n++ ;
01759             }
01760     }
01761     }
01762     
01763     npix = n;
01764     /*if (n != npix)
01765     {
01766         sinfo_msg_error("the computed number of pixel equals "
01767                         "not the counted number, impossible!") ;
01768         cpl_free(retstats) ;
01769         cpl_free(pix_array) ;
01770         return NULL ;
01771     }*/
01772 
01773     /* determining the clean mean is already done in the recipes */
01774     if ( FLT_MAX == (retstats->cleanmean = sinfo_new_clean_mean(pix_array, 
01775                                            npix, loReject, hiReject)) )
01776     {    
01777         sinfo_msg_error("sinfo_new_clean_mean() did not work!") ;
01778         cpl_free(retstats) ;
01779         cpl_free(pix_array) ;
01780         return NULL ;
01781     } 
01782 
01783     /* now the clean standard deviation must be calculated */
01784     /* initialize sums */
01785     lo_n = (int) (loReject / 100. * (float)npix) ;
01786     hi_n = (int) (hiReject / 100. * (float)npix) ;
01787     pix_sum = 0. ;
01788     sqr_sum = 0. ;
01789     n = 0 ;
01790     for ( i = lo_n ; i <= npix - hi_n ; i++ )
01791     {
01792         pix_sum += (double)pix_array[i] ;
01793         sqr_sum += ((double)pix_array[i] * (double)pix_array[i]) ;
01794         n++ ;
01795     }
01796 
01797     if ( n == 0 )
01798     {    
01799         sinfo_msg_error("number of clean pixels is zero!") ;
01800         cpl_free(retstats) ;
01801         cpl_free(pix_array) ;
01802         return NULL ;
01803     } 
01804     retstats -> npix = n ;
01805     pix_sum /= (double) n ;
01806     sqr_sum /= (double) n ;
01807     retstats -> cleanstdev = (float)sqrt(sqr_sum - pix_sum * pix_sum) ;
01808     cpl_free (pix_array) ;
01809     return retstats ;
01810 }
01811 
01812 
01813 
01822 cpl_image * sinfo_new_normalize_to_central_pixel ( cpl_image * image )
01823 {
01824     int col=0;
01825     int row=0;
01826     int i=0;
01827     int n=0;
01828     float* array=NULL ;
01829     float divisor=0;
01830     cpl_image * retImage=NULL;
01831     int ilx=0;
01832     int ily=0;
01833     float* pidata=NULL;
01834     float* podata=NULL;
01835 
01836     if ( image == NULL )
01837     {
01838         sinfo_msg_error("no image given!") ;
01839         return NULL ;
01840     }
01841     ilx=cpl_image_get_size_x(image);
01842     ily=cpl_image_get_size_y(image);
01843     pidata=cpl_image_get_data_float(image);
01844 
01845     retImage = cpl_image_duplicate(image) ;
01846     podata=cpl_image_get_data_float(retImage);
01847 
01848     n = 0 ;
01849     /* go through the central two image rows and store 
01850        the values in an array */
01851     array=cpl_calloc(2*ilx,sizeof(float)) ;
01852 
01853     for ( row = ily/2 ; row < ily/2+1 ; row++ )
01854     {
01855         for ( col = 0 ; col < ilx ; col++ )
01856         {
01857             if ( !isnan(pidata[col+ilx*row]) )
01858             {
01859                 array[n] = pidata[col+ilx*row] ;
01860                 n++ ;
01861             }
01862         }
01863     }
01864     /* compute the sinfo_median of the central 2 spectral 
01865        values of all spatial pixels*/
01866     if ( isnan(divisor = sinfo_new_median(array, n) ) )
01867     {
01868         sinfo_msg_error("no sinfo_median possible!") ;
01869         return NULL ;
01870     }
01871     if ( 0 == divisor )
01872     {
01873         sinfo_msg_error("cannot divide by 0") ;
01874         return NULL ;
01875     }
01876 
01877     for ( i = 0 ; i < (int) ilx*ily ; i++ )
01878     {
01879         if ( isnan(pidata[i]) )
01880         {
01881             podata[i] = ZERO ;
01882         }
01883         else
01884         {
01885             podata[i] = pidata[i]/divisor ;
01886         }
01887     }
01888     cpl_free(array);
01889     return retImage ;
01890 }
01891 
01892 /*-------------------------------------------------------------------------*/
01921 /*--------------------------------------------------------------------------*/
01922 
01923 cpl_image *
01924 sinfo_new_mpe_shift_image(
01925     cpl_image    *    image_in,
01926     double           shift_x,
01927     double           shift_y,
01928     double       *   interp_kernel)
01929 {
01930     cpl_image    *       shifted=NULL ;
01931     pixelvalue  *       first_pass=NULL ;
01932     pixelvalue  *       second_pass=NULL ;
01933     int             samples = KERNEL_SAMPLES ;
01934     int          i=0, j=0 ;
01935     double           fx=0, fy=0 ;
01936     double           rx=0, ry=0 ;
01937     int             px=0, py=0 ;
01938     int             tabx=0, taby=0 ;
01939     double           value=0 ;
01940     size_t          pos ;
01941     register pixelvalue     *   pix ;
01942     register pixelvalue     *   pixint ;
01943     int             mid=0;
01944     double          norm=0 ;
01945     double       *      ker=NULL ;
01946     int                         freeKernel = 1 ;
01947 
01948     int ilx=0;
01949     int ily=0;
01950     float* pidata=NULL;
01951     float* psdata=NULL;
01952 
01953 
01954     /* error handling: test entries */
01955     if (image_in==NULL) return NULL ;
01956 
01957     /* Shifting by a zero offset returns a copy of the input image */
01958     if ((fabs(shift_x)<1e-2) && (fabs(shift_y)<1e-2))
01959       return cpl_image_duplicate(image_in) ;
01960     ilx=cpl_image_get_size_x(image_in);
01961     ily=cpl_image_get_size_y(image_in);
01962     pidata=cpl_image_get_data_float(image_in);
01963 
01964 
01965         /* See if a kernel needs to be generated */
01966     if (interp_kernel == NULL) {
01967         ker = sinfo_generate_interpolation_kernel("default") ;
01968         if (ker == NULL) {
01969             sinfo_msg_error("kernel generation failure:aborting resampling") ;
01970             return NULL ;
01971         }
01972     } else {
01973         ker = interp_kernel ;
01974         freeKernel = 0 ;
01975     }
01976 
01977     mid = (int)samples/(int)2 ;
01978     first_pass = cpl_calloc(ilx, ily*sizeof(pixelvalue)) ;
01979     shifted = cpl_image_new(ilx, ily,CPL_TYPE_FLOAT) ;
01980     psdata=cpl_image_get_data_float(shifted);
01981     
01982     second_pass = psdata ;
01983 
01984     pix = pidata ;
01985     if ( ilx != 1 )
01986     { 
01987     for (j=0 ; j<ily ; j++) 
01988         {
01989          for (i=0 ; i<ilx ; i++) {
01990            fx = (double)i-shift_x ;
01991            px = (int)fx ;
01992            rx = fx - (double)px ;
01993            pos = px + j * ilx ;
01994 
01995            if ((px>1) && (px<(ilx-2))) 
01996            {
01997            tabx = (int)(fabs((double)mid * rx)) ;
01998            /* exclude blank (ZERO) pixels from interpolation */
01999            if (isnan(pix[pos]))
02000                {
02001                value = ZERO ;
02002                }
02003            else
02004                {
02005                if (isnan(pix[pos-1]))
02006                {
02007                pix[pos-1] = 0. ;
02008                }
02009                if (isnan(pix[pos+1]))
02010                {
02011                pix[pos+1] = 0. ;
02012                }
02013                if (isnan(pix[pos+2]))
02014                {
02015                pix[pos+2] = 0. ;
02016                }
02017                
02018                /*
02019             * Sum up over 4 closest pixel values,
02020             * weighted by interpolation kernel values
02021             */
02022                value =     (double)pix[pos-1] * ker[mid+tabx] +
02023                (double)pix[pos] * ker[tabx] +
02024                (double)pix[pos+1] * ker[mid-tabx] +
02025                (double)pix[pos+2] * ker[samples-tabx-1] ;
02026                /*
02027             * Also sum up interpolation kernel coefficients
02028             * for further normalization
02029             */
02030                norm =      (double)ker[mid+tabx] +
02031                (double)ker[tabx] +
02032                (double)ker[mid-tabx] +
02033                (double)ker[samples-tabx-1] ;
02034                if (fabs(norm) > 1e-4) {
02035                value /= norm ;
02036                }
02037                }
02038            } else {
02039            value = ZERO ;
02040            }  
02041            /*
02042         * There may be a problem of rounding here if pixelvalue
02043         * has not enough bits to sustain the accuracy.
02044         */
02045            if ( isnan(value) )
02046            {
02047            first_pass[i+j*ilx] = ZERO ;
02048            }
02049            else
02050            {
02051            first_pass[i+j*ilx] = (pixelvalue)value ;
02052            }
02053          }
02054         }
02055     } 
02056     else
02057     {
02058     memcpy(first_pass,pix,ily*sizeof(pixelvalue));
02059     }
02060 
02061     pixint = first_pass ;
02062     for (i=0 ; i<ilx ; i++) {
02063         for (j=0 ; j<ily ; j++) {
02064             fy = (double)j - shift_y ;
02065             py = (int)fy ;
02066             ry = fy - (double)py ;
02067             pos = i + py * ilx ;
02068 
02069             taby = (int)(fabs((double)mid * ry)) ;
02070 
02071             if ((py>(int)1) && (py<(ily-2))) {
02072                 /* exclude blank (ZERO) pixels from interpolation */
02073                 if (isnan(pixint[pos]) && ilx != 1 )
02074                 {
02075                     value = ZERO ;
02076                 }
02077                 else
02078                 {
02079                     if (isnan(pixint[pos-ilx]))
02080                     {
02081                         pixint[pos-ilx] = 0. ;
02082                     }
02083                     if (isnan(pixint[pos+ilx]))
02084                     {
02085                         pixint[pos+ilx] = 0. ;
02086                     }
02087                     if (isnan(pixint[pos+2*ilx]))
02088                     {
02089                         pixint[pos+2*ilx] = 0. ;
02090                     }
02091                     /*
02092                      * Sum up over 4 closest pixel values,
02093                      * weighted by interpolation kernel values
02094                      */
02095                     value = (double)pixint[pos-ilx] * ker[mid+taby] +
02096                             (double)pixint[pos] * ker[taby] +
02097                             (double)pixint[pos+ilx] * ker[mid-taby] +
02098                             (double)pixint[pos+2*ilx]*ker[samples-taby-1];
02099                     /*
02100                      * Also sum up interpolation kernel coefficients
02101                      * for further normalization
02102                      */
02103                     norm =      (double)ker[mid+taby] +
02104                                 (double)ker[taby] +
02105                                 (double)ker[mid-taby] +
02106                                 (double)ker[samples-taby-1] ;
02107 
02108                     if (fabs(norm) > 1e-4) {
02109                         value /= norm ;
02110                     }
02111                 }
02112             } else {
02113                 value = ZERO ;
02114             }  
02115             if (isnan(value))
02116             {
02117                 second_pass[i+j*ilx] = ZERO ;
02118             }
02119             else
02120             {
02121                 second_pass[i+j*ilx] = (pixelvalue)value ;
02122             }
02123         }
02124     }
02125 
02126     cpl_free(first_pass) ;
02127     if (freeKernel)
02128         cpl_free(ker) ;
02129     return shifted ;
02130 }
02131 
02132 
02133 /*-------------------------------------------------------------------------*/
02164 /*--------------------------------------------------------------------------*/
02165 
02166 void
02167 sinfo_new_shift_image_in_cube(
02168     cpl_image     *   image_in,
02169     double           shift_x,
02170     double           shift_y,
02171     double       *   interp_kernel,
02172     cpl_image     *   shifted,
02173     pixelvalue   *   first_pass)
02174 {  
02175     pixelvalue  *       second_pass=NULL ;
02176     int             samples = KERNEL_SAMPLES ;
02177     int          i=0, j=0 ;
02178     double           fx=0, fy=0 ;
02179     double           rx=0, ry=0 ;
02180     int             px=0, py=0 ;
02181     int             tabx=0, taby=0 ;
02182     double           value=0 ;
02183     size_t          pos ;
02184     register pixelvalue     *   pix ;
02185     register pixelvalue     *   pixint ;
02186     int             mid=0;
02187     double          norm=0 ;
02188     double       *      ker=NULL ;
02189     int                         freeKernel = 1 ;
02190 
02191     int ilx=0;
02192     int ily=0;
02193     int slx=0;
02194     int sly=0;
02195     float* pidata=NULL;
02196     float* psdata=NULL;
02197 
02198     /* error handling: test entries */
02199         if (image_in==NULL) shifted = NULL ;
02200         pidata=cpl_image_get_data_float(image_in);
02201         ilx=cpl_image_get_size_x(image_in);
02202         ily=cpl_image_get_size_y(image_in);
02203 
02204         shifted=cpl_image_new(ilx,ily,CPL_TYPE_FLOAT);
02205         slx=ilx;
02206         sly=ily;
02207 
02208         psdata=cpl_image_get_data_float(shifted);
02209 
02210         /* Shifting by a zero offset returns a copy of the input image */
02211         if ((fabs(shift_x)<1e-2) && (fabs(shift_y)<1e-2))
02212                 memcpy(psdata,pidata, (size_t) slx*sly * sizeof(pixelvalue)) ;
02213 
02214         /* See if a kernel needs to be generated */
02215     if (interp_kernel == NULL) {
02216         ker = sinfo_generate_interpolation_kernel("default") ;
02217         if (ker == NULL) {
02218             sinfo_msg_error("kernel generation failure:aborting resampling") ;
02219             shifted = NULL ;
02220         }
02221     } else {
02222         ker = interp_kernel ;
02223         freeKernel = 0 ;
02224     }
02225 
02226     mid = (int)samples/(int)2 ;
02227     second_pass = psdata ;
02228 
02229     pix = pidata ;
02230     for (j=0 ; j<ily ; j++) {
02231         for (i=1 ; i<ilx-2 ; i++) {
02232             fx = (double)i-shift_x ;
02233             px = (int)fx ;
02234             rx = fx - (double)px ;
02235 
02236             pos = px + j * ilx ;
02237 
02238             if ((px>1) && (px<(ilx-2))) {
02239                 tabx = (int)(fabs((double)mid * rx)) ;
02240                 /* exclude blank (ZERO) pixels from interpolation */
02241                 if (isnan(pix[pos]))
02242                 {
02243                     value = ZERO ;
02244                 }
02245                 else
02246                 {
02247                     if (isnan(pix[pos-1]))
02248                     {
02249                         pix[pos-1] = 0. ;
02250                     }
02251                     if (isnan(pix[pos+1]))
02252                     {
02253                         pix[pos+1] = 0. ;
02254                     }
02255                     if (isnan(pix[pos+2]))
02256                     {
02257                         pix[pos+2] = 0. ;
02258                     }
02259                    
02260                     /*
02261                      * Sum up over 4 closest pixel values,
02262                      * weighted by interpolation kernel values
02263                      */
02264                     value =     (double)pix[pos-1] * ker[mid+tabx] +
02265                                 (double)pix[pos] * ker[tabx] +
02266                                 (double)pix[pos+1] * ker[mid-tabx] +
02267                                 (double)pix[pos+2] * ker[samples-tabx-1] ;
02268                     /*
02269                      * Also sum up interpolation kernel coefficients
02270                      * for further normalization
02271                      */
02272                     norm =      (double)ker[mid+tabx] +
02273                                 (double)ker[tabx] +
02274                                 (double)ker[mid-tabx] +
02275                                 (double)ker[samples-tabx-1] ;
02276                     if (fabs(norm) > 1e-4) {
02277                         value /= norm ;
02278                     }
02279                 }
02280             } else {
02281                 value = 0.0 ;
02282             }  
02283             /*
02284              * There may be a problem of rounding here if pixelvalue
02285              * has not enough bits to sustain the accuracy.
02286              */
02287             if ( isnan(value) )
02288             {
02289                 first_pass[i+j*ilx] = ZERO ;
02290             }
02291             else
02292             {
02293                 first_pass[i+j*ilx] = (pixelvalue)value ;
02294             }
02295         }
02296     }
02297     pixint = first_pass ;
02298     for (i=0 ; i< ilx ; i++) {
02299         for (j=1 ; j< ily-2 ; j++) {
02300             fy = (double)j - shift_y ;
02301             py = (int)fy ;
02302             ry = fy - (double)py ;
02303             pos = i + py * ilx ;
02304 
02305             taby = (int)(fabs((double)mid * ry)) ;
02306 
02307             if ((py>(int)1) && (py<(ily-2))) {
02308                 /* exclude blank (ZERO) pixels from interpolation */
02309                 if (isnan(pixint[pos]))
02310                 {
02311                     value = ZERO ;
02312                 }
02313                 else
02314                 {
02315                     if (isnan(pixint[pos-ilx]))
02316                     {
02317                         pixint[pos-ilx] = 0. ;
02318                     }
02319                     if (isnan(pixint[pos+ilx]))
02320                     {
02321                         pixint[pos+ilx] = 0. ;
02322                     }
02323                     if (isnan(pixint[pos+2*ilx]))
02324                     {
02325                         pixint[pos+2*ilx] = 0. ;
02326                     }
02327                     /*
02328                      * Sum up over 4 closest pixel values,
02329                      * weighted by interpolation kernel values
02330                      */
02331                     value = (double)pixint[pos-ilx] * ker[mid+taby] +
02332                             (double)pixint[pos] * ker[taby] +
02333                             (double)pixint[pos+ilx] * ker[mid-taby] +
02334                             (double)pixint[pos+2*ilx]*ker[samples-taby-1];
02335                     /*
02336                      * Also sum up interpolation kernel coefficients
02337                      * for further normalization
02338                      */
02339                     norm =      (double)ker[mid+taby] +
02340                                 (double)ker[taby] +
02341                                 (double)ker[mid-taby] +
02342                                 (double)ker[samples-taby-1] ;
02343 
02344                     if (fabs(norm) > 1e-4) {
02345                         value /= norm ;
02346                     }
02347                 }
02348             } else {
02349           /* value = 0.0 ; AMo: This affect slitlet #1 */
02350             }  
02351             if (isnan(value))
02352             {
02353                 second_pass[i+j*ilx] = ZERO ;
02354             }
02355             else
02356             {
02357                 second_pass[i+j*ilx] = (pixelvalue)value ;
02358             }
02359         }
02360     }
02361 
02362     if (freeKernel)
02363         cpl_free(ker) ;
02364 }
02365 
02366 /* function to delete the image statistics within python */
02367 void sinfo_new_del_Stats( Stats * st)
02368 {
02369     cpl_free (st) ;
02370 }
02371 
02378 cpl_image * 
02379 sinfo_new_combine_masks ( cpl_image * firstMask, cpl_image * secondMask ) 
02380 {
02381     cpl_image * retMask=NULL ;
02382     int n=0 ;
02383     int olx=0;
02384     int oly=0;
02385     float* podata=NULL;
02386     float* pm1data=NULL;
02387     float* pm2data=NULL;
02388 
02389     if ( firstMask == NULL || secondMask == NULL )
02390     {
02391         sinfo_msg_error ("no input mask image given!") ;
02392         return NULL ;
02393     }
02394     retMask = cpl_image_duplicate (firstMask) ;
02395     podata = cpl_image_get_data_float(retMask);
02396     pm1data = cpl_image_get_data_float(firstMask);
02397     pm2data = cpl_image_get_data_float(secondMask);
02398     olx=cpl_image_get_size_x(retMask);
02399     oly=cpl_image_get_size_y(retMask);
02400 
02401     for ( n = 0 ; n < (int) olx*oly ; n++ )
02402     {
02403        if ( podata[n] == 0. || pm2data[n] == 0. )
02404        {
02405            podata[n] = 0. ;
02406        }
02407        else
02408        {
02409            podata[n] = 1. ;
02410        }
02411     }
02412     return retMask ;
02413 }
02414 
02423 cpl_image * sinfo_new_slice_cube (cpl_imagelist * cube, int x, int y ) 
02424 {
02425     cpl_image * retImage=NULL ;
02426     int col=0, row=0, z=0 ;
02427     int inp=0;
02428     int ilx=0;
02429     int ily=0;
02430     cpl_image* img=NULL;
02431     float* podata=NULL;
02432     float* pidata=NULL;
02433 
02434     if ( cube == NULL )
02435     {
02436         sinfo_msg_error("no cube given!") ;
02437         return NULL ;
02438     }
02439     if ( x > 31 || y > 31 )
02440     {
02441         sinfo_msg_warning("wrong x or y values!") ;
02442     }
02443 
02444     img=cpl_imagelist_get(cube,0);
02445     ilx=cpl_image_get_size_x(img);
02446     ily=cpl_image_get_size_y(img);
02447     inp=cpl_imagelist_get_size(cube);
02448     if ( x < 0 )
02449     {
02450         /* allocate memory */
02451         if ( NULL == (retImage = cpl_image_new(ilx, inp, CPL_TYPE_FLOAT)) ) 
02452         {
02453             sinfo_msg_error("could not allocate memory!") ;
02454             return NULL ;
02455         }
02456         podata=cpl_image_get_data_float(retImage);
02457         for ( z = 0 ; z < inp ; z++ )
02458         {
02459             
02460             pidata=cpl_image_get_data_float(cpl_imagelist_get(cube,z));
02461             for ( col = 0 ; col < ilx ; col++ )  
02462             {
02463                 podata[col+z*ilx] = pidata[col+y*ilx] ;
02464             }
02465         }
02466     }
02467     else if ( y < 0 )
02468     {
02469         /* allocate memory */
02470         if ( NULL == (retImage = cpl_image_new(ily, inp,CPL_TYPE_FLOAT)) ) 
02471         {
02472             sinfo_msg_error("could not allocate memory!") ;
02473             return NULL ;
02474         }
02475         podata=cpl_image_get_data_float(retImage);
02476 
02477         for ( z = 0 ; z < inp ; z++ )
02478         {
02479             pidata=cpl_image_get_data_float(cpl_imagelist_get(cube,z));
02480             for ( row = 0 ; row < ily ; row++ )  
02481             {
02482                 podata[row+z*ily] = pidata[x+row*ily] ;
02483             }
02484         }
02485     }
02486     else 
02487     {
02488         sinfo_msg_error("wrong input!") ;
02489         return NULL ;
02490     }
02491     return retImage ;
02492 }
02493 
02505 cpl_image * sinfo_new_div_images_robust ( cpl_image * im1, cpl_image * im2 ) 
02506 {
02507     cpl_image * retIm=NULL ;
02508     float help=0 ;
02509     int i=0 ;
02510     int lx1=0;
02511     int ly1=0;
02512     int lx2=0;
02513     int ly2=0;
02514 
02515     float* p1data=NULL;
02516     float* p2data=NULL;
02517     float* podata=NULL;
02518   
02519     if ( im1 == NULL || im2 == NULL ) 
02520     {
02521         sinfo_msg_error("no input images given!") ;
02522         return NULL ;
02523     }
02524     lx1=cpl_image_get_size_x(im1);
02525     ly1=cpl_image_get_size_y(im1);
02526     lx2=cpl_image_get_size_x(im2);
02527     ly2=cpl_image_get_size_y(im2);
02528     p1data=cpl_image_get_data_float(im1);
02529     p2data=cpl_image_get_data_float(im2);
02530 
02531     if ( lx1 != lx2 || ly1 != ly2 ) 
02532     {
02533         sinfo_msg_error("images not compatible!") ;
02534         return NULL ;
02535     }
02536     if ( NULL == (retIm = cpl_image_new(lx1, ly1, CPL_TYPE_FLOAT)) ) 
02537     {
02538         sinfo_msg_error("could not allocate memory!") ;
02539         return NULL ;
02540     }
02541     podata=cpl_image_get_data_float(retIm);
02542 
02543     for ( i = 0 ; i < (int) lx2*ly2 ; i++ )
02544     {
02545         if ( !isnan(p2data[i]) )
02546         {
02547             help = 1./p2data[i] ;
02548             if (fabs( help )> THRESH )
02549             {
02550                help = 1. ;
02551             }
02552         }
02553         else
02554         {
02555             help = ZERO ;
02556         }
02557         if ( isnan(help) || isnan(p1data[i]) )
02558         {
02559             podata[i] = ZERO ;
02560         }
02561         else
02562         {
02563             podata[i] = p1data[i] * help ;
02564         }
02565     }
02566     return retIm ;
02567 }
02568 
02569 cpl_image * sinfo_new_null_edges ( cpl_image * image) 
02570 {
02571     cpl_image * new=NULL ;
02572     int i=0,j=0 ;
02573     int ilx=0;
02574     int ily=0;
02575     int olx=0;
02576     int oly=0;
02577 
02578     float* pidata=NULL;
02579     float* podata=NULL;
02580 
02581     if ( image == NULL )
02582     {
02583         sinfo_msg_error ("no input image given!\n") ;
02584         return NULL ;
02585     }
02586 
02587 
02588     new = cpl_image_duplicate (image) ;
02589     ilx=cpl_image_get_size_x(image);
02590     ily=cpl_image_get_size_y(image);
02591     olx=cpl_image_get_size_x(new);
02592     oly=cpl_image_get_size_y(new);
02593     pidata=cpl_image_get_data_float(image);
02594     podata=cpl_image_get_data_float(new);
02595     
02596     for ( i = 0 ; i < olx ; i++ )
02597     {
02598         for ( j = 0 ; j < 4 ; j++)
02599     {
02600         podata[i+j*olx]=0;
02601         podata[i+(oly-j-1)*olx]=0;
02602     }
02603     }
02604     for ( i = 0 ; i < oly ; i++ )
02605     {
02606         for ( j = 0 ; j < 4 ; j++)
02607     {
02608         podata[j+i*olx]=0;
02609         podata[(olx-j-1)+i*olx]=0;
02610     }
02611     }
02612     return new ;
02613 }
02614 
02615 
02616 void sinfo_new_used_cor_map( cpl_image *im, cpl_image *map)
02617 {
02618     int i=0,j=0,index=0;
02619     float temp_array[2048];
02620     int lx=cpl_image_get_size_x(im);
02621     int ly=cpl_image_get_size_y(im);
02622     float* pidata=cpl_image_get_data_float(im);
02623     float* pmdata=cpl_image_get_data_float(map);
02624 
02625     for( j=0; j<ly; j++)
02626     {
02627     for( i=0;i<lx;i++)
02628         {
02629           index = (int)pmdata[i+j*lx];
02630           temp_array[i] = pidata[index+j*lx];
02631         }
02632     for( i=0;i<lx;i++)
02633         {
02634           pidata[i+j*lx]= temp_array[i];
02635         }
02636     }
02637 }
02638 
02639 
02640 
02641 
02642 /*-------------------------------------------------------------------------*/
02665 /*--------------------------------------------------------------------------*/
02666 
02667 cpl_image *
02668 sinfo_new_shift_image(
02669     cpl_image   *   image_in,
02670     double           shift_x,
02671     double           shift_y,
02672     double       *  interp_kernel)
02673 {
02674     cpl_image    *  shifted=NULL ;
02675     float  *    first_pass=NULL ;
02676     float  *    second_pass=NULL ;
02677     int             samples = KERNEL_SAMPLES ;
02678     int          i=0, j=0 ;
02679     double           fx=0, fy=0 ;
02680     double           rx=0, ry=0 ;
02681     int             px=0, py=0 ;
02682     int             tabx=0, taby=0 ;
02683     double           value=0 ;
02684     size_t          pos ;
02685     register float     *    pix=NULL ;
02686     register float     *    pixint=NULL ;
02687     int             mid=0;
02688     double          norm=0 ;
02689     double       *  ker=NULL ;
02690     int             freeKernel = 1 ;
02691     int ilx=0;
02692     int ily=0;
02693 
02694     /* error handling: test entries */
02695     if (image_in==NULL) return NULL ;
02696 
02697     /* Shifting by a zero offset returns a copy of the input image */
02698     if ((fabs(shift_x)<1e-2) && (fabs(shift_y)<1e-2))
02699         return cpl_image_duplicate(image_in) ;
02700 
02701     /* See if a kernel needs to be generated */
02702     if (interp_kernel == NULL) {
02703         ker = sinfo_generate_interpolation_kernel("default") ;
02704         if (ker == NULL) {
02705             sinfo_msg_error("kernel generation failure: aborting resampling") ;
02706             return NULL ;
02707         }
02708     } else {
02709         ker = interp_kernel ;
02710         freeKernel = 0 ;
02711     }
02712 
02713     ilx=cpl_image_get_size_x(image_in);
02714     ily=cpl_image_get_size_y(image_in);
02715 
02716     mid = (int)samples/(int)2 ;
02717     first_pass = cpl_calloc(ilx, ily*sizeof(float)) ;
02718     shifted = cpl_image_new(ilx, ily,CPL_TYPE_FLOAT) ;
02719     second_pass = cpl_image_get_data_float(shifted);
02720 
02721     pix = cpl_image_get_data_float(image_in);
02722     for (j=0 ; j<ily ; j++) {
02723         for (i=1 ; i<ilx-2 ; i++) {
02724             fx = (double)i-shift_x ;
02725             px = (int)fx ;
02726             rx = fx - (double)px ;
02727 
02728             pos = px + j * ilx ;
02729 
02730             if ((px>1) && (px<(ilx-3))) {
02731                 tabx = (int)(fabs((double)mid * rx)) ;
02732                 /* 
02733                  * Sum up over 4 closest pixel values, 
02734                  * weighted by interpolation kernel values
02735                  */
02736                 value =     (double)pix[pos-1] * ker[mid+tabx] +
02737                             (double)pix[pos] * ker[tabx] + 
02738                             (double)pix[pos+1] * ker[mid-tabx] +
02739                             (double)pix[pos+2] * ker[samples-tabx-1] ;
02740                 /*
02741                  * Also sum up interpolation kernel coefficients
02742                  * for further normalization
02743                  */
02744                 norm =      (double)ker[mid+tabx] + 
02745                             (double)ker[tabx] + 
02746                             (double)ker[mid-tabx] + 
02747                             (double)ker[samples-tabx-1] ;
02748                 if (fabs(norm) > 1e-4) {
02749                     value /= norm ;
02750                 }
02751             } else {
02752                 value = 0.0 ;
02753             }   
02754             /* 
02755              * There may be a problem of rounding here if pixelvalue
02756              * has not enough bits to sustain the accuracy.
02757              */
02758             first_pass[i+j*ilx] = (float)value ;
02759         }
02760     }
02761     pixint = first_pass ;
02762     for (i=0 ; i<ilx ; i++) {
02763         for (j=1 ; j<ily-3 ; j++) {
02764             fy = (double)j - shift_y ;
02765             py = (int)fy ;
02766             ry = fy - (double)py ;
02767             pos = i + py * ilx ;
02768 
02769             taby = (int)(fabs((double)mid * ry)) ;
02770 
02771             if ((py>(int)1) && (py<(ily-2))) {
02772                 /* 
02773                  * Sum up over 4 closest pixel values, 
02774                  * weighted by interpolation kernel values
02775                  */
02776                 value = (double)pixint[pos-ilx] * ker[mid+taby] + 
02777                         (double)pixint[pos] * ker[taby] + 
02778                         (double)pixint[pos+ilx] * ker[mid-taby] + 
02779                         (double)pixint[pos+2*ilx]*ker[samples-taby-1];
02780                 /*
02781                  * Also sum up interpolation kernel coefficients
02782                  * for further normalization
02783                  */
02784                 norm =      (double)ker[mid+taby] + 
02785                             (double)ker[taby] + 
02786                             (double)ker[mid-taby] +
02787                             (double)ker[samples-taby-1] ;
02788 
02789                 if (fabs(norm) > 1e-4) {
02790                     value /= norm ;
02791                 }
02792             } else {
02793                 value = 0.0 ;
02794             }   
02795             second_pass[i+j*ilx] = (float)value ;
02796         }
02797     }
02798 
02799     cpl_free(first_pass) ;
02800     if (freeKernel)
02801         cpl_free(ker) ;
02802     return shifted ;
02803 }
02804 
02805 /*--------------------------------------------------------------------------*/

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