image_ops.c

00001 /*******************************************************************************
00002 * M.P.E. - SPIFFI project
00003 *
00004 *
00005 *
00006 * who       when      what
00007 * --------  --------  ----------------------------------------------
00008 * rabuter 2004-12-03 support one dimensional image in shiftImage
00009 * schreib  23/05/00  created
00010 */
00011 
00012 /************************************************************************
00013 *   NAME
00014 *        image_ops.c -
00015 *        image arithmetic routines 
00016 *
00017 *   SYNOPSIS
00018 *   #include "image_ops.h"
00019 *
00020 *   1) Vector * meanOfColumns( OneImage *im )
00021 *   2) Vector * cleanMeanOfColumns( OneImage *im, 
00022 *                                   double lo_reject, 
00023 *                                   double hi_reject)
00024 *   3) OneImage * divImageByRow( OneImage *im, Vector *row )
00025 *   4) OneImage * multRowToImage( OneImage *im, Vector *row )
00026 *   5) OneImage * colTilt ( OneImage * image, float sigmaFactor )
00027 *   6) OneImage * medianImage( OneImage * im, float fmedian )
00028 *   7) OneImage * compareImages( OneImage * im1, OneImage * im2, OneImage * origim )
00029 *   8) OneImage * threshImage ( OneImage * im, float lo_cut, float hi_cut )
00030 *   9) pixel_map * promoteImageToPixelmap ( OneImage * im )
00031 *  10) OneImage * promoteImageToMask ( OneImage * im, int * n_badpixels )
00032 *  11) OneImage * multImageByMask ( OneImage * im, OneImage * mask )
00033 *  12) OneImage * interpolImage ( OneImage * im,
00034 *                                 OneImage * mask,
00035 *                                 int        max_radius,
00036 *                                 int        n_pixels )
00037 *  13) OneImage * interpolSourceImage ( OneImage * im,
00038 *                                       OneImage * mask,
00039 *                                       int        max_rad,
00040 *                                       float   ** slit_edges )
00041 *  14) OneImage * stackRowToImage ( Vector * row, int ly )
00042 *  15) Stats * imageStatsOnRectangle ( OneImage * im, 
00043 *                                      float      loReject,
00044 *                                      float      hiReject,
00045 *                                      int        llx, 
00046 *                                      int        lly, 
00047 *                                      int        urx, 
00048 *                                      int        ury )
00049 *  16) OneImage * normalize_to_central_pixel ( OneImage * image )
00050 *  17) OneImage *
00051 *      shiftImage(
00052 *            OneImage    *    image_in,
00053 *            double           shift_x,
00054 *            double           shift_y,
00055 *            double       *   interp_kernel)
00056 *  18) OneImage * combineMasks ( OneImage * firstMask, OneImage * secondMask ) 
00057 *  19) OneImage * sliceCube ( OneCube * cube, int x, int y ) 
00058 *  20) OneImage * divImagesRobust ( OneImage * im1, OneImage * im2 ) 
00059 *
00060 *
00061 *   DESCRIPTION
00062 *   1) takes the average of each image column
00063 *   2) takes the average of each image column by sorting the
00064 *      column values and rejecting the given percentage of 
00065 *      the highest and lowest values  [0...1]
00066 *   3) divides each image column by a row value
00067 *   4) multiplies each image column with a row value
00068 *   5) first calculates statistics for each column of an image.
00069 *      median value and standard deviation of columns are de-
00070 *      termined, blank values are excluded. Fits a straight 
00071 *      line through the pixel values of each column and subtracts
00072 *      the fit in order to remove the tilt of each column.
00073 *      Only those pixels are used for the fit that are within 
00074 *      a defined factor of sigma noise limit. The noise is 
00075 *      calculated from pixels between the 10percentil and 
00076 *      90percentil points. 
00077 *      if the straight line could not be determined, the median
00078 *      of the column is subtracted from the column
00079 *   6) median filter, calculates the median for an image
00080 *      by using the 8 closest pixels to each pixel.
00081 *      The values in the output image are determined according
00082 *      to the values of the input parameter. 
00083 *      If fmedian = 0: always replace by median
00084 *      if fmedian < 0: replace by median if |pixel - median| >
00085 *                      -fmedian
00086 *      if fmedian > 0: replace by median (fmedian as a factor of 
00087 *                      the square root of the median itself)                                  
00088 *                      if |pixel - median| >= fmedian * sqrt ( median ) 
00089 *                      This can be used to consider photon noise.
00090 *                      This considers a dependence of the differences on the 
00091 *                      pixel values themselves.
00092 *   7) if a pixel value of one image (im1) equals 
00093 *      the pixel value of the other image keep the 
00094 *      pixel value of the original image otherwise replace 
00095 *      it with ZEROs
00096 *   8) simple search for static bad pixels for a flat field
00097 *      or dark frame, values below and above the threshold
00098 *      values are set to ZERO.
00099 *   9) changes an image with ZERO indicated bad pixels to
00100 *      a bad pixel map.
00101 *  10) changes an image with ZERO indicated bad pixels to
00102 *      a bad pixel mask image, that means the returned
00103 *      image has values 1 at positions of good pixels and 
00104 *      ZEROs at positions of bad pixels.
00105 *  11) changes an image to an image that has ZERO indicated 
00106 *      static bad pixels 
00107 *  12) interpolates all bad pixels indicated by the bad pixel mask. 
00108 *      Therefore, the mean of at least 2 valid values of 
00109 *      the nearest 8 neighbors is taken. If too much 
00110 *      neighbors are also bad pixels 
00111 *      the neighbor radius is increased to a maximum of 
00112 *      max_radius until n_pixels valid pixels are found.
00113 *      The valid neighbors are searched by going through
00114 *      the columns and rows around the central square that
00115 *      was already searched.
00116 *      The bad pixel is interpolated by the mean of these 
00117 *      valid pixels (less than 9) or by the median of them
00118 *      (more than 8).
00119 *  13) interpolates all bad pixels indicated by the bad pixel mask. 
00120 *      Therefore, the mean of the nearest 4 neighbors is taken, 
00121 *      two in spectral direction and 2 in spatial direction.
00122 *      The routine cares about the image and slitlet edges.
00123 *      If there are no good pixel found within the nearest neighbors,
00124 *      the next 4 nearest neighbors in spatial and spectral direction
00125 *      are searched for valid pixels until a limit of max_rad.   
00126 *      A maximum of 4 valid pixels are used for interpolation by their mean.
00127 *  14) stack a given image row to build a whole image
00128 *  15) computes the mean and standard deviation of 
00129 *      a given rectangle on an image by leaving the extreme
00130 *      intensity values.
00131 *  16) normalizes a raw flatfield image by dividing by the median of the central
00132 *      spectral pixels to produce a master flatfield
00133 *  17) This function is a copy of the ECLIPSE function shift_image()
00134 *      but slightly changed. If a blank (ZERO) pixel appears the blank pixel
00135 *      is shifted but preserved as blank.
00136 *      If a blank (ZERO) pixel appears within the
00137 *      interpolation kernel the blank pixel is set to 0. 
00138 * 
00139 *      This function shifts an image by a non-integer offset, using
00140 *      interpolation. You can either generate an interpolation kernel once and
00141 *      pass it to this function, or let it generate a default kernel. In the
00142 *      former case, use generate_interpolation_kernel() to generate an
00143 *      appropriate kernel. In the latter case, pass NULL as last argument. A
00144 *      default interpolation kernel is then generated then discarded before this
00145 *      function returns.
00146 *
00147 *      The returned image is a newly allocated object, it must be deallocated
00148 *      using destroy_image().
00149 *  18) combines two bad pixel mask to one using an or relation
00150 *  19) slices a data cube in x or y direction
00151 *  20) divides two images by considering blanks and
00152 *      calculating first 1/im2 by
00153 *      cutting the very high values and setting to 1,
00154 *      then multiplying im1 * 1/im2.  
00155 *
00156 *   FILES
00157 *
00158 *   ENVIRONMENT
00159 *
00160 *   RETURN VALUES
00161 *
00162 *   CAUTIONS
00163 *
00164 *   EXAMPLES
00165 *
00166 *   SEE ALSO
00167 *
00168 *   BUGS
00169 *
00170 *------------------------------------------------------------------------
00171 */
00172 
00173 #include "vltPort.h"
00174 
00175 /*
00176  * System Headers
00177  */
00178 
00179 /*
00180  * Local Headers
00181  */
00182 
00183 #include "image_ops.h"
00184 
00185 /*----------------------------------------------------------------------------
00186  *                          Function codes
00187  *--------------------------------------------------------------------------*/
00188 
00189 
00190 
00191 
00192 double median_image(OneImage* im)
00193 {
00194   double m=0;
00195   register int i=0;
00196   int n=0;
00197   pixelvalue* pv=0;
00198   
00199    for ( i = 0 ; i < (int) im -> nbpix ; i++ )
00200     {
00201       if ( qfits_isnan(im -> data[i]) )
00202         {
00203 
00204     } else {
00205       n++;
00206     }
00207     }
00208    pv = cpl_calloc(n,sizeof(pixelvalue));
00209    n=0;
00210    for ( i = 0 ; i < (int) im -> nbpix ; i++ )
00211     {
00212       if ( qfits_isnan(im -> data[i]) )
00213         {
00214 
00215     } else {
00216       pv[n]=im->data[i];
00217           n++;
00218     }
00219     }
00220    if(pv == NULL || n == 0) {
00221      m=0;
00222    } else {
00223      m=median(pv,n);
00224    }
00225    cpl_free(pv);
00226    if(qfits_isnan(m)){
00227      m=0;
00228    }
00229   return m;
00230 }
00231 /*----------------------------------------------------------------------------
00232    Function     :       meanOfColumns()
00233    In           :       image
00234    Out          :       resulting row array
00235    Job          :       takes the average of each image column
00236  ---------------------------------------------------------------------------*/
00237 
00238 Vector * meanOfColumns( OneImage *im )
00239 {
00240     Vector * row ;
00241     int         i, j ;
00242 
00243     if ( im == NullImage )
00244     {
00245         cpl_msg_error ("meanOfColumns:","null image") ;
00246         return NullVector ;
00247     }
00248     
00249     /* allocate memory for a row with the length of the image x-axis */
00250     if ( NULL == (row = newVector (im -> lx)) )
00251     {
00252         cpl_msg_error ("meanOfColumns:","not able to allocate a vector" ) ;
00253         return NullVector ;
00254     }
00255 
00256     for ( i = 0 ; i < im -> lx ; i++ )
00257     {
00258         for ( j = 0 ; j < im -> ly ; j++ )
00259         {
00260             if (!qfits_isnan(im->data[i+j*im->lx]))
00261             {
00262                 row -> data[i] += im -> data[i + j*(im->lx)] ;
00263             }
00264         }
00265 
00266         row -> data[i] /= im -> ly ;
00267     }
00268     return row ;
00269 }
00270 
00271 
00272 /*----------------------------------------------------------------------------
00273    Function     :       cleanMeanOfColumns()
00274    In           :       image , percentage of lowest and highest values to 
00275                         reject
00276    Out          :       resulting row image
00277    Job          :       takes the average of each image column by sorting the
00278                         column values and rejecting the given percentage of 
00279                         the highest and lowest values  [0...1]
00280  ---------------------------------------------------------------------------*/
00281 
00282 OneImage * cleanMeanOfColumns( OneImage *im, 
00283                              float lo_reject, 
00284                              float hi_reject)
00285 {
00286     OneImage     * row ;
00287     pixelvalue   buffer[im->ly] ;
00288     int          i, j, k, nv ;
00289     int          lo_n, hi_n ;
00290 
00291     if ( im == NULL )
00292     {
00293         cpl_msg_error ("cleanMeanOfColumns:","null image") ;
00294         return NullImage ;
00295     }
00296 
00297     if ((lo_reject + hi_reject) > 0.9)
00298     {
00299         cpl_msg_error("cleanMeanOfColumns:","illegal rejection thresholds: [%f] and [%f]", lo_reject, hi_reject) ;
00300         cpl_msg_error("cleanMeanOfColumns:","threshold sum should not be over 0.90  aborting average") ;
00301         return NullImage ;
00302     }
00303 
00304     lo_n = (int) (im -> ly * lo_reject + 0.5) ;
00305     hi_n = (int) (im -> ly * hi_reject + 0.5) ;
00306     if (lo_n + hi_n >= im -> ly)
00307     {
00308         cpl_msg_error ("cleanMeanOfColumns:","everything would be rejected") ;
00309         return NullImage ;
00310     }
00311 
00312     /* allocate memory for a row with the length of the image x-axis */
00313     if ( NULL == (row = new_image (im -> lx, 1)) ) 
00314     {
00315         cpl_msg_error ("cleanMeanOfColumns:","cannot allocate new image") ;
00316         return NullImage ;
00317     }
00318 
00319     for ( i = 0 ; i < im -> lx ; i++ )
00320     {
00321         for ( j = 0 ; j < im -> ly ; j++ )
00322         {
00323             buffer[j] = im -> data[i + j*(im->lx)] ;
00324         }
00325         pixel_qsort (buffer, im -> ly) ;
00326         
00327         nv = 0 ;
00328         for (k = lo_n ; k < im->ly - hi_n ; k ++)
00329         {
00330             if ( !qfits_isnan(buffer[k]) )
00331             {
00332                 row -> data[i] += buffer[k] ;
00333                 nv ++ ;
00334             }
00335         }
00336         row -> data[i] /= nv ;
00337  
00338     }
00339     return row ;
00340 }
00341 
00342 
00343 /*----------------------------------------------------------------------------
00344    Function     :       divImageByRow()
00345    In           :       image, row array
00346    Out          :       resulting image
00347    Job          :       divides each image column by a row value
00348  ---------------------------------------------------------------------------*/
00349 
00350 OneImage * divImageByRow( OneImage *im, Vector *row )
00351 {
00352     OneImage *image ;
00353     int         i,j ;
00354 
00355     if ( im == NULL || row == NULL )
00356     {
00357         cpl_msg_error ("divImageByRow:","null image or null row") ;
00358         return NullImage ;
00359     }
00360 
00361     if ( im -> lx != row -> n_elements )
00362     {
00363         cpl_msg_error("divImageByRow:","image and row size not compatible") ;
00364         return NullImage ;
00365     }
00366     
00367     if ( NullImage == (image = copy_image (im)) ) 
00368     {
00369         cpl_msg_error ("divImageByRow:","cannot copy image") ;
00370         return NullImage ;
00371     }
00372     
00373     for (i = 0 ; i < im -> lx ; i++ )
00374     {
00375         for (j = 0 ; j < im -> ly ; j++)
00376         {
00377             if ( !qfits_isnan(im -> data[i + j*(im->lx)]) )
00378             {
00379                 image -> data[i + j*(im->lx)] = im -> data[i + j*(im->lx)] / row -> data[i] ;
00380             }
00381         }
00382     }
00383     return image ;
00384 }
00385 
00386 /*----------------------------------------------------------------------------
00387    Function     :       multRowToImage()
00388    In           :       image, row array
00389    Out          :       resulting image
00390    Job          :       multiplies each image column with a row value
00391  ---------------------------------------------------------------------------*/
00392 
00393 OneImage * multRowToImage( OneImage *im, Vector *row )
00394 {
00395     OneImage *image ;
00396     int         i,j ;
00397 
00398     if ( im == NULL || row == NULL )
00399     {
00400         cpl_msg_error ("multRowToImage:","null image or null row") ;
00401         return NullImage ;
00402     }
00403 
00404     if ( im -> lx != row -> n_elements )
00405     {
00406         cpl_msg_error("multRowToImage:","image and row size not compatible") ;
00407         return NullImage ;
00408     }
00409 
00410     if ( NULL == (image = copy_image (im)) ) 
00411     {
00412         cpl_msg_error ("multRowToImage:","cannot copy image") ;
00413         return NullImage ;
00414     }
00415     
00416     for (i = 0 ; i < im -> lx ; i++ )
00417     {
00418         for (j = 0 ; j < im -> ly ; j++)
00419         {
00420             if ( !qfits_isnan(im -> data[i + j*(im->lx)]) )
00421             {
00422                 image -> data[i + j*(im->lx)] = im -> data[i + j*(im->lx)] * row -> data[i] ;
00423             }
00424         }
00425     }
00426     return image ;
00427 }
00428 
00429 
00430 /*----------------------------------------------------------------------------
00431    Function     :       colTilt()
00432    In           :       image , factor of sigma noise limit to determine 
00433                         pixels that are used for the fit.
00434    Out          :       image
00435    Job          :       first calculates statistics for each column of an image.
00436                         median value and standard deviation of columns are de-
00437                         termined, blank values are excluded. Fits a straight 
00438                         line through the pixel values of each column and subtracts
00439                         the fit in order to remove the tilt of each column.
00440                         Only those pixels are used for the fit that are within 
00441                         a defined factor of sigma noise limit. The noise is 
00442                         calculated from pixels between the 10percentil and 
00443                         90percentil points. 
00444                         if the straight line could not be determined, the median
00445                         of the column is subtracted from the column
00446    Notice       :       works only for raw or averaged raw images
00447  ---------------------------------------------------------------------------*/
00448 
00449 OneImage * colTilt ( OneImage * image, float sigmaFactor )
00450 {
00451     OneImage   * im ;
00452     float      * column ;
00453     double       sum, sum2, mean ;
00454     float        median, noise ;
00455     float      * sig, * dat ; 
00456     float        a, b, siga, sigb, chi2, q ;
00457     int          i, j, colnum, npix, mwt ;
00458 
00459 
00460 
00461     if ( image == NullImage )
00462     {
00463         cpl_msg_error ( "colTilt:","no image given" ) ;     
00464         return NullImage ;
00465     }
00466     
00467     if ( sigmaFactor <= 0. )
00468     {
00469         cpl_msg_error ( "colTilt:","no or negative sigma factor") ;
00470         return NullImage ;
00471     }
00472     
00473     /* allocate memory */
00474     if ( NULL == (im = new_image ( image -> lx, image -> ly )) ) 
00475     {
00476         cpl_msg_error ( "colTilt:","cannot allocate new image" ) ;
00477         return NullImage ;
00478     }
00479 
00480     /* go through the columns */
00481     for ( i = 0 ; i < image -> lx ; i ++ )
00482     {
00483         /* initialize the buffer variables for each column */
00484         colnum = 0 ;
00485         column = (float *) cpl_calloc ( image -> ly , sizeof (float *) ) ;
00486         sig    = (float *) cpl_calloc ( image -> ly , sizeof (float *) ) ;
00487         dat    = (float *) cpl_calloc ( image -> ly , sizeof (float *) ) ;
00488 
00489         /*select only non-ZERO values of one column*/
00490         for ( j = 0 ; j < image -> ly ; j++ )    
00491         {
00492             if ( !qfits_isnan(image -> data[i + j*image -> lx]) )
00493             {
00494                 column[j] = image -> data[i + j*image -> lx] ;
00495                 colnum ++ ;
00496             }
00497         }
00498         if ( colnum < 10 ) 
00499         {
00500             /*cpl_msg_warning ("colTilt:","column %d has almost only blank pixels and is set to blank", i+1) ;*/
00501             for ( j = 0 ; j < image -> ly ; j++ )    
00502             {
00503                 im -> data[i + j*image -> lx] = ZERO;
00504             }
00505             /*
00506             cpl_free (column) ;
00507             cpl_free (sig);
00508             cpl_free (dat) ;
00509             continue ; 
00510             */
00511         }
00512 
00513         /*-------------------------------------------------------------------
00514          * sort the data, clip off the extremes, determine the noise
00515          * and get the range for the valid data. It is assumed here
00516          * that most pixels are o.k. 
00517          */
00518 
00519         pixel_qsort (column, colnum) ;
00520         
00521         sum   = 0. ;
00522         sum2  = 0. ;
00523         npix  = 0  ;
00524 
00525         for ( j = 0.1*colnum + 1 ; j <= 0.9*colnum ; j++ )
00526         {
00527             sum  += column[j] ;
00528             sum2 += column[j] * column[j] ;
00529             npix ++ ;
00530         }
00531 
00532         if (npix <= 1)
00533         {
00534             noise = sigmaFactor * 1000.;
00535         }
00536         else
00537         {
00538             mean   = sum/(float)npix ;
00539             noise  = sqrt( (sum2 - sum*mean)/(double)(npix -1) ) ;
00540             noise *= sigmaFactor ;
00541         }
00542 
00543         /* -------------------------------------------------------------
00544          * determine median
00545          * if colnum is odd, median will be the colnum/2 th value, otherwise
00546          * median is the mean of colnum/2-1 th and colnum/2 th value.
00547          */       
00548 
00549         if ( colnum % 2 == 1 )
00550         {
00551             median = column[colnum/2] ;
00552         }
00553         else
00554         {
00555             median = (column[colnum/2 - 1] + column[colnum/2])/2. ;
00556         }
00557 
00558         /* now select the pixels for the tilt calculation */
00559         
00560         colnum = 0 ;
00561         for ( j = 0; j < image -> ly ; j ++ )
00562         {
00563             if ( !qfits_isnan(image -> data[i+j*image->lx]) &&
00564                  fabs ( (image -> data[i+j*image->lx]) - median) <= noise )
00565             {
00566                 column[colnum] = image -> data[i+j*image->lx] ;
00567                 dat[colnum] = (float) j ;
00568                 sig[colnum] = 1. ;
00569                 colnum ++ ;
00570             }
00571         }
00572 
00573         if ( colnum == 0 )
00574         {
00575             /*for ( j = 0; j < image -> ly; j++ )
00576             {
00577                 im -> data[i+j*image->lx] -= median ;
00578             }
00579             cpl_free (column) ;
00580             cpl_free (sig)    ;
00581             cpl_free (dat)    ;
00582             continue ;*/
00583         a=0./0.;
00584         b=0./0.;
00585         }
00586     else
00587     {
00588         mwt = 0 ;
00589         myfit ( dat, column, colnum, sig, mwt, &a, &b, &siga, &sigb, &chi2, &q ) ;
00590     }
00591         if ( fabs(b) >= SLOPE || fabs(a) >= SATURATION  || qfits_isnan(b) || qfits_isnan(a))
00592         {
00593             cpl_msg_warning ("colTilt:","linear fit: slope is greater than limit: %f saturation level is reached: %f in column number %d ", b, a , i+1) ; 
00594         }
00595         
00596         /* subtract fit or median from data */
00597         for ( j = 0; j < image -> ly; j++ )
00598         {
00599             if ( !qfits_isnan(image -> data[i+j*image->lx]) &&
00600                  fabs(b) < SLOPE && fabs(a) < SATURATION )
00601             {
00602                 im -> data[i+j*image->lx] = image -> data[i+j*image->lx] - (a + b * (float)j) ;
00603             }
00604             else if ( qfits_isnan(image -> data[i+j*image->lx]) )
00605             {
00606                 im -> data[i+j*image->lx] = ZERO ;
00607             }
00608             else if ( (fabs(b) >= SLOPE || fabs(a) >= SATURATION || qfits_isnan(a) || qfits_isnan(b)) && 
00609                       !qfits_isnan(image -> data[i+j*image->lx]) )
00610             {
00611                 im -> data[i+j*image->lx] -= median ;
00612             }
00613             else
00614             {
00615                 cpl_msg_error ( "colTilt:"," case is not possible! %f %f", b,a) ;
00616                 /*cpl_free (column) ;
00617                 cpl_free (sig)    ;
00618                 cpl_free (dat)    ;
00619                 destroy_image(im) ;
00620                 return NullImage ;*/
00621             }
00622         }
00623         cpl_free (column) ;
00624         cpl_free (sig)    ;
00625         cpl_free (dat)    ;
00626     }
00627 
00628     return im     ;
00629 }
00630 
00631 
00632 /*----------------------------------------------------------------------------
00633    Function     :       medianImage()
00634    In           :       image, a median threshold parameter
00635    Out          :       resulting image
00636    Job          :       median filter, calculates the median for an image
00637                         by using the 8 closest pixels of every pixel.
00638                         The values in the output image are determined according
00639                         to the values of the input parameter. 
00640                         If fmedian = 0: always replace by median
00641                         if fmedian < 0: replace by median if |pixel - median| >
00642                                         -fmedian
00643                         if fmedian > 0: replace by median (fmedian as a factor of 
00644                                         the square root of the median itself)             
00645                                         if |pixel - median| >= fmedian * sqrt ( median ) 
00646                                         This can be used to consider photon noise.
00647                                         This considers a dependence of the differences on the 
00648                                         pixel values themselves.
00649    Notice       :       it is assumed that most of the 8 nearest neighbor pixels
00650                         are not bad pixels!
00651                         blank pixels are not replaced!
00652  ---------------------------------------------------------------------------*/
00653 
00654 OneImage * medianImage( OneImage * im, float fmedian )
00655 {
00656     OneImage *   image       ;
00657     pixelvalue * value       ;
00658     pixelvalue   median      ;
00659     int        * position    ;
00660     int          nposition   ;
00661     int          n, i, j ;
00662 
00663     if ( im == NullImage )
00664     {
00665         cpl_msg_error ("medianImage:","no image input") ;
00666         return NullImage ;
00667     }
00668     
00669     image = copy_image ( im ) ;
00670  
00671     /*---------------------------------------------------------------------- 
00672      * go through all pixels 
00673      */
00674 
00675     for ( i = 0 ; i < (int) im -> nbpix ; i++ )
00676     {
00677         /* blank pixels are not replaced */ 
00678         if ( qfits_isnan(im -> data[i]) )
00679         {
00680             continue ;
00681         }
00682 
00683         /* initialize the buffer variables for the 8 nearest neighbors */
00684         value = (pixelvalue * )cpl_calloc ( 8, sizeof ( pixelvalue * ) ) ;
00685         position = ( int * ) cpl_calloc ( 8, sizeof ( int * ) ) ;
00686 
00687         /*---------------------------------------------------------------------- 
00688          * determine the pixel position of the 8 nearest neighbors  
00689          */
00690 
00691         position[0] = i + im -> lx - 1 ; /* upper left  */
00692         position[1] = i + im -> lx     ; /* upper       */
00693         position[2] = i + im -> lx + 1 ; /* upper right */
00694         position[3] = i + 1            ; /* right       */
00695         position[4] = i - im -> lx + 1 ; /* lower right */
00696         position[5] = i - im -> lx     ; /* lower       */
00697         position[6] = i - im -> lx - 1 ; /* lower left  */
00698         position[7] = i - 1            ; /* left        */
00699        
00700         /*---------------------------------------------------------------------- 
00701          * determine the positions of the image margins, top positions are changed
00702          * to low positions and vice versa. Right positions are changed to left
00703          * positions and vice versa.
00704          */
00705 
00706         if ( i >= 0 && i < im -> lx )    /* bottom line */
00707         {
00708             position[4] += 2 * im -> lx ;  
00709             position[5] += 2 * im -> lx ;  
00710             position[6] += 2 * im -> lx ;  
00711         }
00712         else if ( i >= ((int) im -> nbpix - im -> lx ) && i < (int) im -> nbpix ) /* top line */
00713         {
00714             position[0] -= 2 * im -> lx ;  
00715             position[1] -= 2 * im -> lx ;  
00716             position[2] -= 2 * im -> lx ;  
00717         }
00718         else if ( i % im -> lx == 0 )    /* left side */
00719         {
00720             position[0] += 2 ;
00721             position[6] += 2 ;
00722             position[7] += 2 ;
00723         }
00724         else if ( i % im -> lx == im -> lx - 1 )    /* right side */
00725         {
00726             position[2] -= 2 ;
00727             position[3] -= 2 ;
00728             position[4] -= 2 ;
00729         }
00730 
00731         /* ----------------------------------------------------------------------------
00732          * read the pixel values of the neighboring pixels,
00733          * blanks are not considered 
00734          */
00735 
00736         nposition = 8 ;
00737         n = 0 ;
00738         for ( j = 0 ; j < nposition ; j ++ )
00739         {
00740             if ( !qfits_isnan(im -> data[position[j]]) )
00741             {
00742                 value[n] = im -> data[position[j]] ;    
00743                 n ++ ;
00744             }
00745         }
00746         nposition = n ;
00747 
00748         if ( nposition <= 1 )  /* almost all neighbors are blank */
00749         {
00750             image->data[i] = ZERO ;
00751             cpl_free(value) ;
00752             cpl_free(position) ;
00753             continue ;
00754         }
00755         
00756         /* sort the values and determine the median */
00757 
00758         pixel_qsort ( value, nposition ) ;
00759         if ( nposition % 2 == 1 )
00760         {
00761             median = value [ nposition/2 ] ;
00762         }
00763         else
00764         {
00765             median = ( value [nposition/2 - 1] + value [nposition/2] ) / 2. ;
00766         }
00767 
00768         /* -----------------------------------------------------------------
00769          * replace the pixel value by the median on conditions:
00770          * fmedian = 0: always replace with median.
00771          * fmedian < 0: interpret as absolute condition:
00772          *              if |pixel - median| > -fmedian
00773          *              replace with median.
00774          * fmedian > 0: replace by median (fmedian as a factor of 
00775          *              the square root of the median itself)                                  
00776          *              if |pixel - median| >= fmedian * sqrt ( median ) 
00777          *              considers a dependence on the pixel value.
00778          *              This can be used to consider photon noise.
00779          */
00780 
00781         if ( fmedian == 0 )
00782         {
00783             image -> data[i] = median ;
00784         }
00785         else if ( fmedian < 0 &&
00786                   fabs ( median - im -> data[i] ) >= -fmedian )
00787         {
00788             image -> data[i] = median ;
00789         }
00790         else if ( fmedian > 0 &&
00791                   fabs ( median - im -> data[i] ) >= fmedian * sqrt(fabs(median)) )
00792         {
00793             image -> data[i] = median ;
00794         }   
00795         else
00796         {
00797             cpl_free (value) ; 
00798             cpl_free (position) ; 
00799             continue ;
00800         }
00801 
00802         cpl_free (value) ;
00803         cpl_free (position) ;
00804     }       
00805     return image ;
00806 }
00807 
00808 
00809 /*----------------------------------------------------------------------------
00810    Function     :       compareImages()
00811    In           :       two images to be compared and the original image
00812    Out          :       resulting image
00813    Job          :       if a pixel value of one image (im1) equals 
00814                         the pixel value of the other image keep the 
00815                         pixel value of the original image otherwise replace 
00816                         it with ZEROs
00817  ---------------------------------------------------------------------------*/
00818 
00819 OneImage * compareImages( OneImage * im1, OneImage * im2, OneImage * origim )
00820 {
00821     OneImage * image ;
00822     int            i ;
00823 
00824     if ( im1 == NullImage || im2 == NullImage || origim == NullImage )
00825     {
00826         cpl_msg_error ("compareImages:","Null images as input" ) ;
00827         return NullImage ;
00828     }
00829 
00830     if ( im1 -> lx != im2 -> lx ||
00831          im1 -> ly != im2 -> ly )
00832     {
00833         cpl_msg_error ("compareImages:","incompatible image sizes" ) ;
00834         return NullImage ;
00835     }
00836     
00837     /* allocate memory */
00838     if ( NULL == (image = new_image ( im1 -> lx, im1 -> ly )) ) 
00839     {
00840         cpl_msg_error ("compareImages:","cannot allocate new image" ) ;
00841         return NullImage ;
00842     }
00843 
00844     for ( i = 0 ; i < (int) im1 -> nbpix ; i ++ )
00845     {
00846         if ( qfits_isnan(im1 -> data[i]) && qfits_isnan(im2 -> data[i]) )
00847         {
00848             image -> data[i] = ZERO ;
00849         }
00850         else
00851         {
00852             if ( im1 -> data[i] == im2 -> data[i] )
00853             {
00854                 image -> data[i] = origim -> data[i] ;
00855             }
00856             else
00857             {
00858                 image -> data[i] = ZERO ;
00859             }
00860         }
00861     }
00862     return image ;
00863 }
00864 
00865 
00866 /*----------------------------------------------------------------------------
00867    Function     :       threshImage()
00868    In           :       image, low cut pixel value, high cut pixel value
00869    Out          :       resulting image
00870    Job          :       simple search for static bad pixels for a flat field
00871                         or dark frame, values below and above the threshold
00872                         values are set to ZERO.
00873  ---------------------------------------------------------------------------*/
00874 
00875 OneImage * threshImage ( OneImage * im, float lo_cut, float hi_cut )
00876 {
00877     OneImage * image ;
00878     int            i ;
00879     
00880     if (im == NullImage)
00881     {
00882         cpl_msg_error ("threshImage:","null image given") ;
00883         return NullImage ;
00884     }
00885 
00886     image = copy_image(im) ;
00887     for ( i = 0 ; i < (int) im -> nbpix ; i ++ )
00888     {
00889         if ( im -> data[i] > (pixelvalue) hi_cut || 
00890              im -> data[i] < (pixelvalue) lo_cut )
00891         {
00892              image -> data[i] = ZERO ;
00893         }
00894     }
00895     return image ;
00896 }
00897 
00898     
00899 /*----------------------------------------------------------------------------
00900    Function     :       promoteImageToPixelmap()
00901    In           :       image
00902    Out          :       resulting pixelmap
00903    Job          :       changes an image with ZERO indicated bad pixels to
00904                         a bad pixel map.
00905  ---------------------------------------------------------------------------*/
00906 
00907 pixel_map * promoteImageToPixelmap ( OneImage * im )
00908 {
00909     pixel_map * newmap ;
00910     int              i ;
00911 
00912     if ( im == NullImage )
00913     {
00914         cpl_msg_error ( "promoteImageToPixelmap:","null image given" ) ;
00915         return NullMap ;
00916     }
00917 
00918     if ( NullMap == (newmap = new_pixelmap ( im->lx, im->ly )) )
00919     {
00920         cpl_msg_error ("promoteImageToPixelmap:","cannot allocate new pixelmap") ;
00921         return NullMap ;
00922     } 
00923     
00924     for ( i = 0 ; i < (int) im -> nbpix ; i ++ )
00925     {
00926         if ( qfits_isnan(im -> data[i]) )
00927         {
00928             newmap -> data[i] = ZERO ;
00929             newmap -> ngoodpix -- ;
00930         }
00931     }
00932 
00933     return newmap ;
00934 }
00935 
00936 /*----------------------------------------------------------------------------
00937    Function     :       promoteImageToMask()
00938    In           :       image with ZERO indicating bad pixels
00939    Out          :       resulting mask image that means 1 for good pixels
00940                         and 0 for bad pixel positions
00941                         n_badpixels: number of bad pixels
00942    Job          :       changes an image with ZERO indicated bad pixels to
00943                         a bad pixel mask image, that means the returned
00944                         image has values 1 at positions of good pixels and 
00945                         0 at positions of bad pixels.
00946  ---------------------------------------------------------------------------*/
00947 
00948 OneImage * promoteImageToMask ( OneImage * im, int * n_badpixels )
00949 {
00950     OneImage * reImage ;
00951     int        i ;
00952 
00953     if ( NullImage == im )
00954     {
00955         cpl_msg_error("promoteImageToMask:","no input image given!") ;
00956         return NullImage ;
00957     }
00958 
00959     /* allocate memory for the returned image */
00960     if ( NullImage == (reImage = new_image ( im->lx, im->ly )) )
00961     {
00962         cpl_msg_error ("promoteImageToMask:","cannot allocate new image!") ;
00963         return NullImage ;
00964     } 
00965     
00966     *n_badpixels = 0 ;
00967     for ( i = 0 ; i < (int) im -> nbpix ; i ++ )
00968     {
00969         if ( qfits_isnan(im -> data[i]) )
00970         {
00971             reImage -> data[i] = 0. ;
00972             (*n_badpixels)++ ;
00973         }
00974         else
00975         {
00976             reImage -> data[i] = 1. ;
00977         }
00978     }
00979     return reImage ;
00980 }
00981 
00982 
00983 /*----------------------------------------------------------------------------
00984    Function     :       multImageByMask()
00985    In           :       im:   input image
00986                         mask: mask image 
00987    Out          :       resulting image that means the input image with marked
00988                         static bad pixels (ZERO values)
00989    Job          :       changes an image to an image that has ZERO indicated 
00990                         static bad pixels 
00991  ---------------------------------------------------------------------------*/
00992 
00993 OneImage * multImageByMask ( OneImage * im, OneImage * mask )
00994 {
00995     OneImage * reImage ;
00996     int        i ;
00997 
00998     if ( NullImage == im )
00999     {
01000         cpl_msg_error("multImageByMask:","no input image given!") ;
01001         return NullImage ;
01002     }
01003     if ( NullImage == mask )
01004     {
01005         cpl_msg_error("multImageByMask:","no mask image given!") ;
01006         return NullImage ;
01007     }
01008     if ( im ->lx != mask->lx || im->ly != mask->ly)
01009     {
01010         cpl_msg_error("multImageByMask:","image sizes are not correspondent!") ;
01011         return NullImage ;
01012     }
01013 
01014     reImage = copy_image( im ) ;
01015 
01016     for ( i = 0 ; i < (int) im -> nbpix ; i ++ )
01017     {
01018         if ( mask -> data[i] == 0. )
01019         {
01020             reImage -> data[i] = ZERO ; 
01021         }
01022     }
01023 
01024     return reImage ;
01025 }
01026 
01027 /*----------------------------------------------------------------------------
01028    Function     :       interpolImage()
01029    In           :       im: raw image 
01030                         mask: bad pixel mask
01031                         max_radius: maximum x and y distance in pixels from the
01032                                     bad pixel within which valid pixels for 
01033                                     interpolation are searched.
01034                         n_pixels:   minimal number of pixels with which the bad
01035                                     pixel is interpolated if not enough 
01036                                     valid nearest neighbors are found.
01037    Out          :       resulting interpolated image without any ZEROS 
01038    Job          :       interpolates all bad pixels indicated by the bad pixel mask. 
01039                         Therefore, the mean of at least 2 valid values of 
01040                         the nearest 8 neighbors is taken. If too much 
01041                         neighbors are also bad pixels 
01042                         the neighbor radius is increased to a maximum of 
01043                         max_radius until n_pixels valid pixels are found.
01044                         The valid neighbors are searched by going through
01045                         the columns and rows around the central square that
01046                         was already searched.
01047                         The bad pixel is interpolated by the mean of these 
01048                         valid pixels (less than 9) or by the median of them
01049                         (more than 8).
01050  ---------------------------------------------------------------------------*/
01051 
01052 OneImage * interpolImage ( OneImage * im,
01053                            OneImage * mask,
01054                            int        max_radius,
01055                            int        n_pixels )
01056 {
01057     OneImage * returnImage ;
01058     float neighbors[4*max_radius*max_radius] ;
01059     float sum, mean ;
01060     int i, j, k ;
01061     int row, col ;
01062     int n_valid ;
01063     int agreed ;
01064 
01065     if ( NullImage == im )
01066     {
01067         cpl_msg_error("interpolImage:","sorry, no input image given!") ;
01068         return NullImage ;
01069     }
01070     if ( NullImage == mask )
01071     {
01072         cpl_msg_error("interpolImage:","sorry, no mask image given!") ;
01073         return NullImage ;
01074     }
01075     if ( mask->lx != im->lx || mask->ly != im->ly )
01076     {
01077         cpl_msg_error("interpolImage:","images not compatible !") ;
01078         return NullImage ;
01079     }
01080 
01081     if ( max_radius <= 0 )
01082     {
01083         cpl_msg_error("interpolImage:","wrong number of pixels for maximal search radius given!") ;
01084         return NullImage ;
01085     }
01086 
01087     if ( n_pixels <= 2 )
01088     {
01089         cpl_msg_error("interpolImage:","wrong number of pixels used for interpolation given!") ;
01090         return NullImage ;
01091     }
01092 
01093     returnImage = copy_image ( im ) ;
01094 
01095     /* go through the columns and rows of the input and mask image */
01096     for ( col = 0 ; col < im->lx ; col++ )
01097     {
01098         for ( row = 0 ; row < im->ly ; row++ )
01099         {
01100             /* look for the ZEROS that means the detected bad pixels */
01101             if ( qfits_isnan(mask->data[col+row*im->lx]) || mask->data[col+row*im->lx] == 0. )
01102             {
01103                 /* now the neighbors must be considered */
01104                 n_valid = 0 ;
01105                 agreed  = 0 ;
01106                 for ( j = 1 ; j <= max_radius ; j++ )
01107                 {
01108 
01109                     /* go through the left column */
01110                     for ( k = -j ; k < j ; k++ )
01111                     {
01112                         if ( col-j >= 0 && row+k < im->ly && row+k >= 0 )
01113                         {
01114                             if ( !qfits_isnan(mask->data[col-j+(row+k)*mask->lx]) || 
01115                                  mask->data[col-j+(row+k)*mask->lx] != 0 )
01116                             {
01117                                 neighbors[n_valid] = im->data[col-j+(row+k)*im->lx] ;
01118                                 n_valid++ ;
01119                             }
01120                         }
01121                     }
01122 
01123                     /* go through the upper row */
01124                     for ( k = -j ; k < j ; k++ )
01125                     {
01126                         if ( col+k < im->lx && col+k >= 0 && row+j < im->ly )
01127                         {
01128                             if ( !qfits_isnan(mask->data[col+k+(row+j)*mask->lx]) || 
01129                                  mask->data[col+k+(row+j)*mask->lx] != 0. )
01130                             {
01131                                 neighbors[n_valid] = im->data[col+k+(row+j)*im->lx] ;
01132                                 n_valid++ ;
01133                             }
01134                         }
01135                     }
01136 
01137                     /* go through the right column */
01138                     for ( k = -j ; k < j ; k++ )
01139                     {
01140                         if ( col+j < im->lx  && row-k >= 0 && row-k < im->ly )
01141                         {
01142                             if ( !qfits_isnan(mask->data[col+j+(row-k)*mask->lx]) || 
01143                                  mask->data[col+j+(row-k)*mask->lx] != 0. )
01144                             {
01145                                 neighbors[n_valid] = im->data[col+j+(row-k)*im->lx] ;
01146                                 n_valid++ ;
01147                             }
01148                         }
01149                     }
01150 
01151                     /* go through the lower row */
01152                     for ( k = -j ; k < j ; k++ )
01153                     {
01154                         if ( col-k >= 0 && col-k < im->lx && row-j < im->ly )
01155                         {
01156                             if ( !qfits_isnan(mask->data[col-k+(row-j)*mask->lx]) ||
01157                                  mask->data[col-k+(row-j)*mask->lx] != 0. )
01158                             {
01159                                 neighbors[n_valid] = im->data[col-k+(row-j)*im->lx] ;
01160                                 n_valid++ ;
01161                             }
01162                         }
01163                     }
01164 
01165                     /* control if the breaking criteria is fullfilled */
01166                     if ( n_valid >= n_pixels )
01167                     {
01168                         agreed = 1 ;
01169                         break ;
01170                     }
01171                     /* do a break if more than 2 nearest neighbors are found */
01172                     if ( j == 1 && n_valid >= 2 )
01173                     {
01174                         agreed = 1 ;
01175                         break ;
01176                     }
01177                 }
01178                 if ( n_valid < n_pixels && agreed == 0 )
01179                 {
01180                     cpl_msg_error("interpolImage:","not enough valid neighbors found to interpolate \
01181                                bad pixel in col: %d, row: %d\n", col, row ) ;
01182                     return NullImage ;
01183                 }
01184                 else
01185                 {
01186                     /* --------------------------------------------------------- 
01187                      * take the mean of the valid neighboring pixels if less 
01188                      * than 9 valid pixels are available else take the median.
01189                      */
01190                     if ( n_valid <= 8 )
01191                     {
01192                         sum = 0. ;
01193                          
01194                         for ( i = 0 ; i < n_valid ; i++ )
01195                         {
01196                             sum += neighbors[i] ;
01197                         }
01198                         mean = sum / n_valid ;
01199 
01200                         returnImage->data[col+row*im->lx] = mean ;
01201                     }
01202                     else
01203                     {
01204                         returnImage->data[col+row*im->lx] = median(neighbors, n_valid) ;
01205                     }
01206                 }
01207             }
01208         }
01209     }
01210 
01211     return returnImage ;
01212 }
01213 
01214 
01215 /*----------------------------------------------------------------------------
01216    Function     :       interpolSourceImage()
01217    In           :       im:         source raw image
01218                         mask:       bad pixel mask
01219                         max_rad:    maximum pixel distance from the bad pixel to 
01220                                     interpolate where valid pixel values are searched for. 
01221                         slit_edges: array of the edges of the slitlets. 
01222    Out          :       resulting interpolated image hopefully without any ZEROS 
01223    Job          :       interpolates all bad pixels indicated by the bad pixel mask. 
01224                         Therefore, the mean of the nearest 4 neighbors is taken, 
01225                         two in spectral direction and 2 in spatial direction.
01226                         The routine cares about the image and slitlet edges.
01227                         If there are no good pixel found within the nearest neighbors,
01228                         the next 4 nearest neighbors in spatial and spectral direction
01229                         are searched for valid pixels until a limit of max_rad.   
01230                         A maximum of 4 valid pixels are used for interpolation by their mean.
01231  ---------------------------------------------------------------------------*/
01232 
01233 OneImage * interpolSourceImage ( OneImage * im,
01234                                  OneImage * mask,
01235                                  int        max_rad,
01236                                  float   ** slit_edges )
01237 {
01238     OneImage * returnImage ;
01239     float validpixel[6] ;
01240     float sum ;
01241     int n, row, col ;
01242     int i, k ;
01243     int slitlet ;
01244     int n_slitlets ;
01245 
01246     if ( NullImage == im )
01247     {
01248         cpl_msg_error("interpolSourceImage:","sorry, no input image given!") ;
01249         return NullImage ;
01250     }
01251     if ( NullImage == mask )
01252     {
01253         cpl_msg_error("interpolSourceImage:","sorry, no bad pixel mask image given!") ;
01254         return NullImage ;
01255     }
01256     if ( mask->lx != im->lx || mask->ly != im->ly )
01257     {
01258         cpl_msg_error("interpolSourceImage:","images not compatible in size!") ;
01259         return NullImage ;
01260     }
01261 
01262     if ( max_rad < 1 )
01263     {
01264         cpl_msg_error("interpolSourceImage:","sorry, wrong maximum distance given!") ;
01265         return NullImage ;
01266     }
01267 
01268     if ( slit_edges == NULL )
01269     {
01270         cpl_msg_error("interpolSourceImage:","sorry, array slit_edges is empty!") ;
01271         return NullImage ;
01272     }
01273  
01274     /* determine the number of slitlets */
01275     n_slitlets = N_SLITLETS ;
01276  
01277     /* copy the original image in the image that will be returned */
01278     returnImage = copy_image( im ) ;
01279 
01280     /* go through the rows and columns of the image and search for the bad pixels */
01281     for ( row = 0 ; row < im -> ly ; row++ )
01282     {
01283         for ( col = 0 ; col < im -> lx ; col++ )
01284         {
01285             n = 0 ;
01286             if ( qfits_isnan(mask -> data[col + row*mask->lx]) || 
01287                  mask -> data[col + row*mask->lx] == 0. ||
01288                  qfits_isnan(im -> data[col + row*mask->lx]) )
01289             {
01290                 /* look for the slitlet where the bad pixel is found */
01291                 slitlet = -1000 ;
01292                 for ( k = 0 ; k < n_slitlets ; k++ )
01293                 {
01294                     if ( nint(slit_edges[k][0]) <= col && nint(slit_edges[k][1]) >= col )
01295                     {
01296                         slitlet = k ;
01297                     }
01298 /* The following else statement is wrong, because in the end slitlet will always be -1000
01299             else
01300                     {
01301                         slitlet = -1000 ;
01302                     } 
01303 */
01304                 }      
01305                 for ( i = 0 ; i < 6 ; i++ )
01306                 {
01307                     validpixel[i] = 0. ;
01308                 }
01309                 /* look for the valid nearest neighbors and collect them but only a maximum of 4 */
01310                 for ( i = 1 ; i <= max_rad ; i++ )
01311                 {
01312                     if ( row + i < im->ly)
01313                     {
01314                         if ( !qfits_isnan(mask -> data[col + (row+i) * mask->lx])
01315                              && mask -> data[col + (row+i) * mask->lx] != 0. &&
01316                                 !qfits_isnan(im -> data[col + (row+i) * im->lx]) )
01317                         {
01318                             validpixel[n] = im -> data[col + (row+i) * im->lx] ;
01319                             n++ ;
01320                         }
01321                     }
01322                     if ( row - i >= 0 )
01323                     {
01324                         if ( !qfits_isnan(mask -> data[col + (row-i) * mask->lx]) 
01325                              && mask -> data[col + (row-i) * mask->lx] != 0. &&
01326                                 !qfits_isnan(im -> data[col + (row-i) * im->lx]) )
01327                         {
01328                             validpixel[n] = im -> data[col + (row-i) * im->lx] ;
01329                             n++ ;
01330                         }
01331                     }
01332 
01333                     /* be aware of the slitlet edges in the spatial direction */
01334                     if ( col + i < im -> lx )  
01335                     {
01336                         if (  slitlet != -1000 )
01337                         {
01338                             if ( col + i <= nint(slit_edges[slitlet][1]) &&
01339                                  !qfits_isnan(mask -> data[col + i + row * mask->lx]) && 
01340                                  mask -> data[col + i + row * mask->lx] != 0. &&
01341                                  !qfits_isnan(im -> data[col + i + row * im->lx]) )
01342                             {
01343                                 validpixel[n] = im -> data[col + i + row * im->lx] ;
01344                                 n++ ;
01345                             }
01346                         }
01347                     }
01348                     if ( col - i >= 0 )  
01349                     {
01350                         if ( slitlet != -1000 )
01351                         {
01352                             if ( col - i >= nint(slit_edges[slitlet][0]) &&
01353                                  !qfits_isnan(mask -> data[col - i + row * mask->lx]) &&
01354                                  mask -> data[col - i + row * mask->lx] != 0. && 
01355                                  !qfits_isnan(im -> data[col - i + row * im->lx]) )  
01356                             {
01357                                 validpixel[n] = im -> data[col - i + row * im->lx] ;
01358                                 n++ ;
01359                             }
01360                         }
01361                     }
01362 
01363                     if ( i == 1 && n > 1 )
01364                     {
01365                         break ;
01366                     }
01367                     if ( n > 2 )
01368                     {
01369                         break ;
01370                     }
01371                 }
01372 
01373                 if ( n == 0 )
01374                 {
01375                     returnImage -> data[col + row*im->lx] = ZERO ; 
01376             /*          cpl_msg_warning("interpolSourceImage:","bad pixel in column: %d and row: %d could not be interpolated!", col, row) ; */
01377                 }
01378                 else
01379                 {
01380                     /* now compute the mean and replace the bad pixel value by the mean */ 
01381                     sum = 0. ;
01382                     for ( i = 0 ; i < n ; i++ )
01383                     {
01384                         sum += validpixel[i] ;
01385                     }
01386                     returnImage -> data[col + row*im->lx] = sum/n ; 
01387                 }
01388             }
01389         }
01390     }
01391 
01392     return returnImage ;
01393 }
01394 
01395 /*----------------------------------------------------------------------------
01396    Function     :       stackRowToImage()
01397    In           :       row: one image row as vector data structure
01398                         ly:  image length
01399    Out          :       resulting image
01400    Job          :       stack a given image row to build a whole image
01401  ---------------------------------------------------------------------------*/
01402 
01403 OneImage * stackRowToImage ( Vector * row, int ly )
01404 {
01405     OneImage * image ;
01406     int        col, ro ;
01407 
01408     if ( row == NullVector )
01409     {
01410         cpl_msg_error ("stackRowToImage:","Null vector as input" ) ;
01411         return NullImage ;
01412     }
01413     if ( ly <= 1 )
01414     {
01415         cpl_msg_error ("stackRowToImage:","wrong image length given" ) ;
01416         return NullImage ;
01417     }
01418 
01419     /* allocate memory */
01420     if ( NULL == (image = new_image ( row -> n_elements ,ly )) ) 
01421     {
01422         cpl_msg_error ("compareImages:","cannot allocate new image" ) ;
01423         return NullImage ;
01424     }
01425 
01426     for ( col = 0 ; col < row -> n_elements ; col++ )
01427     {
01428         for ( ro = 0 ; ro < ly ; ro++ )
01429         {
01430             image -> data[col + ro*ly] = row -> data[col] ;
01431         }
01432     }
01433     return image ;
01434 }
01435 
01436 /*----------------------------------------------------------------------------
01437    Function     :       imageStatsOnRectangle()
01438    In           :       im: flatfield image to search for bad pix
01439                         loReject, hiReject: percentage (0...100) of extrem values
01440                                             that should not be considered
01441                         llx, lly: lower left pixel position of rectangle
01442                         urx, ury: upper right pixel position of rectangle
01443    Out          :       data structure giving the mean and standard deviation
01444    Job          :       computes the mean and standard deviation of 
01445                         a given rectangle on an image by leaving the extreme
01446                         intensity values.
01447  ---------------------------------------------------------------------------*/
01448 
01449 Stats * imageStatsOnRectangle ( OneImage * im, 
01450                                 float      loReject,
01451                                 float      hiReject,
01452                                 int        llx, 
01453                                 int        lly, 
01454                                 int        urx, 
01455                                 int        ury )
01456 {
01457     Stats * retstats ;
01458     int i ;
01459     int row, col ;
01460     int n, npix ;
01461     int lo_n, hi_n ;
01462     double pix_sum, sqr_sum ;
01463     float * pix_array ;
01464 
01465     if ( NullImage == im )
01466     {
01467         cpl_msg_error("imageStatsOnRectangle:","sorry, no input image given!") ;
01468         return NULL ;
01469     }
01470     if ( loReject+hiReject >= 100. )
01471     {
01472         cpl_msg_error("imageStatsOnRectangle:","sorry, too much pixels rejected!") ;
01473         return NULL ;
01474     }
01475     if ( loReject < 0. || loReject >= 100. || hiReject < 0. || hiReject >=100. )
01476     {
01477         cpl_msg_error("imageStatsOnRectangle:","sorry, negative reject values!") ;
01478         return NULL ;
01479     }
01480     if ( llx < 0 || lly < 0 || urx < 0 || ury < 0 ||
01481          llx >= im->lx || lly >= im->ly || urx >= im->lx || 
01482          ury >= im->ly || ury <= lly || urx <= llx )
01483     {
01484         cpl_msg_error("imageStatsOnRectangle:","sorry, wrong pixel coordinates of rectangle!") ;
01485         return NULL ;
01486     }
01487 
01488      /* allocate memory */
01489     retstats = (Stats*) cpl_calloc(1, sizeof(Stats)) ;
01490     npix = (urx - llx + 1) * (ury - lly + 1) ;
01491     pix_array = (float*) cpl_calloc ( npix, sizeof(float) ) ;
01492 
01493     /*------------------------------------------------------------------------- 
01494      * go through the rectangle and copy the pixel values into an array.
01495      */
01496     n = 0 ;
01497     for ( row = lly ; row <= ury ; row++ )
01498     {
01499         for ( col = llx ; col <= urx ; col++ )
01500         {
01501             if ( !qfits_isnan(im -> data[col + row*im->lx]) )
01502             {
01503                 pix_array[n] = im -> data[col + row*im->lx] ;
01504                 n++ ;
01505             }
01506     }
01507     }
01508     
01509     npix = n;
01510     /*if (n != npix)
01511     {
01512         cpl_msg_error("imageStatsOnRectangle:","the computed number of pixel equals not the counted number, impossible!") ;
01513         cpl_free(retstats) ;
01514         cpl_free(pix_array) ;
01515         return NULL ;
01516     }*/
01517 
01518     /* determining the clean mean is already done in the recipes */
01519     if ( FLT_MAX == (retstats->cleanmean = clean_mean(pix_array, npix, loReject, hiReject)) )
01520     {    
01521         cpl_msg_error("imageStatsOnRectangle:","clean_mean() did not work!") ;
01522         cpl_free(retstats) ;
01523         cpl_free(pix_array) ;
01524         return NULL ;
01525     } 
01526 
01527     /* now the clean standard deviation must be calculated */
01528     /* initialize sums */
01529     lo_n = (int) (loReject / 100. * (float)npix) ;
01530     hi_n = (int) (hiReject / 100. * (float)npix) ;
01531     pix_sum = 0. ;
01532     sqr_sum = 0. ;
01533     n = 0 ;
01534     for ( i = lo_n ; i <= npix - hi_n ; i++ )
01535     {
01536         pix_sum += (double)pix_array[i] ;
01537         sqr_sum += ((double)pix_array[i] * (double)pix_array[i]) ;
01538         n++ ;
01539     }
01540 
01541     if ( n == 0 )
01542     {    
01543         cpl_msg_error("imageStatsOnRectangle:","number of clean pixels is zero!") ;
01544         cpl_free(retstats) ;
01545         cpl_free(pix_array) ;
01546         return NULL ;
01547     } 
01548     retstats -> npix = n ;
01549     pix_sum /= (double) n ;
01550     sqr_sum /= (double) n ;
01551     retstats -> cleanstdev = (float)sqrt(sqr_sum - pix_sum * pix_sum) ;
01552     cpl_free (pix_array) ;
01553     return retstats ;
01554 }
01555 
01556 /*----------------------------------------------------------------------------
01557    Function     :       normalize_to_central_pixel()
01558    In           :       image: image to normalize
01559    Out          :       resulting image
01560    Job          :       normalizes a raw flatfield image by dividing by the median of the central
01561                         spectral pixels to produce a master flatfield
01562  ---------------------------------------------------------------------------*/
01563 
01564 OneImage * normalize_to_central_pixel ( OneImage * image )
01565 {
01566     int col, row ;
01567     int i, n ;
01568     float array[2*image->lx] ;
01569     float divisor ;
01570     OneImage * retImage ;
01571 
01572     if ( image == NullImage )
01573     {
01574         cpl_msg_error("normalize_to_central_pixel:","no image given!") ;
01575         return NullImage ;
01576     }
01577     retImage = copy_image(image) ;
01578 
01579     n = 0 ;
01580     /* go through the central two image rows and store the values in an array */
01581     for ( row = image->ly/2 ; row < image->ly/2+1 ; row++ )
01582     {
01583         for ( col = 0 ; col < image->lx ; col++ )
01584         {
01585             if ( !qfits_isnan(image->data[col+image->lx*row]) )
01586             {
01587                 array[n] = image->data[col+image->lx*row] ;
01588                 n++ ;
01589             }
01590         }
01591     }
01592     /* compute the median of the central 2 spectral values of all spatial pixels*/
01593     if ( qfits_isnan(divisor = median(array, n) ) )
01594     {
01595         cpl_msg_error("normalize_to_central_pixel:","no median possible!") ;
01596         return NullImage ;
01597     }
01598     if ( 0 == divisor )
01599     {
01600         cpl_msg_error("normalize_to_central_pixel:","cannot divide by 0") ;
01601         return NullImage ;
01602     }
01603 
01604     for ( i = 0 ; i < (int) image->nbpix ; i++ )
01605     {
01606         if ( qfits_isnan(image->data[i]) )
01607         {
01608             retImage->data[i] = ZERO ;
01609         }
01610         else
01611         {
01612             retImage->data[i] = image->data[i]/divisor ;
01613         }
01614     }
01615     return retImage ;
01616 }
01617 
01618 /*-------------------------------------------------------------------------*/
01647 /*--------------------------------------------------------------------------*/
01648 
01649 OneImage *
01650 shiftImage(
01651     OneImage    *    image_in,
01652     double           shift_x,
01653     double           shift_y,
01654     double       *   interp_kernel)
01655 {
01656     OneImage    *       shifted ;
01657     pixelvalue  *       first_pass ;
01658     pixelvalue  *       second_pass ;
01659     int             samples = KERNEL_SAMPLES ;
01660     int          i, j ;
01661     double           fx, fy ;
01662     double           rx, ry ;
01663     int             px, py ;
01664     int             tabx, taby ;
01665     double           value ;
01666     size_t          pos ;
01667     register pixelvalue     *   pix ;
01668     register pixelvalue     *   pixint ;
01669     int             mid;
01670     double          norm ;
01671     double       *      ker ;
01672     int                         freeKernel = 1 ;
01673     /* error handling: test entries */
01674         if (image_in==NULL) return NULL ;
01675 
01676         /* Shifting by a zero offset returns a copy of the input image */
01677         if ((fabs(shift_x)<1e-2) && (fabs(shift_y)<1e-2))
01678                 return copy_image(image_in) ;
01679 
01680         /* See if a kernel needs to be generated */
01681     if (interp_kernel == NULL) {
01682         ker = generate_interpolation_kernel("default") ;
01683         if (ker == NULL) {
01684             cpl_msg_error("kernel generation failure:","aborting resampling") ;
01685             return NULL ;
01686         }
01687     } else {
01688         ker = interp_kernel ;
01689         freeKernel = 0 ;
01690     }
01691 
01692     mid = (int)samples/(int)2 ;
01693     first_pass = cpl_calloc(image_in->lx, image_in->ly*sizeof(pixelvalue)) ;
01694     shifted = new_image(image_in->lx, image_in->ly) ;
01695     second_pass = shifted->data ;
01696 
01697     pix = image_in->data ;
01698     if ( image_in->lx != 1 )
01699     { 
01700     for (j=0 ; j<image_in->ly ; j++) 
01701         {
01702          for (i=0 ; i<image_in->lx ; i++) {
01703            fx = (double)i-shift_x ;
01704            px = (int)fx ;
01705            rx = fx - (double)px ;
01706            pos = px + j * image_in->lx ;
01707 
01708            if ((px>1) && (px<(image_in->lx-2))) 
01709            {
01710            tabx = (int)(fabs((double)mid * rx)) ;
01711            /* exclude blank (ZERO) pixels from interpolation */
01712            if (qfits_isnan(pix[pos]))
01713                {
01714                value = ZERO ;
01715                }
01716            else
01717                {
01718                if (qfits_isnan(pix[pos-1]))
01719                {
01720                pix[pos-1] = 0. ;
01721                }
01722                if (qfits_isnan(pix[pos+1]))
01723                {
01724                pix[pos+1] = 0. ;
01725                }
01726                if (qfits_isnan(pix[pos+2]))
01727                {
01728                pix[pos+2] = 0. ;
01729                }
01730                
01731                /*
01732             * Sum up over 4 closest pixel values,
01733             * weighted by interpolation kernel values
01734             */
01735                value =     (double)pix[pos-1] * ker[mid+tabx] +
01736                (double)pix[pos] * ker[tabx] +
01737                (double)pix[pos+1] * ker[mid-tabx] +
01738                (double)pix[pos+2] * ker[samples-tabx-1] ;
01739                /*
01740             * Also sum up interpolation kernel coefficients
01741             * for further normalization
01742             */
01743                norm =      (double)ker[mid+tabx] +
01744                (double)ker[tabx] +
01745                (double)ker[mid-tabx] +
01746                (double)ker[samples-tabx-1] ;
01747                if (fabs(norm) > 1e-4) {
01748                value /= norm ;
01749                }
01750                }
01751            } else {
01752            value = ZERO ;
01753            }  
01754            /*
01755         * There may be a problem of rounding here if pixelvalue
01756         * has not enough bits to sustain the accuracy.
01757         */
01758            if ( qfits_isnan(value) )
01759            {
01760            first_pass[i+j*image_in->lx] = ZERO ;
01761            }
01762            else
01763            {
01764            first_pass[i+j*image_in->lx] = (pixelvalue)value ;
01765            }
01766          }
01767         }
01768     } 
01769     else
01770     {
01771     memcpy(first_pass,pix,image_in->ly*sizeof(pixelvalue));
01772     }
01773 
01774     pixint = first_pass ;
01775     for (i=0 ; i<image_in->lx ; i++) {
01776         for (j=0 ; j<image_in->ly ; j++) {
01777             fy = (double)j - shift_y ;
01778             py = (int)fy ;
01779             ry = fy - (double)py ;
01780             pos = i + py * image_in->lx ;
01781 
01782             taby = (int)(fabs((double)mid * ry)) ;
01783 
01784             if ((py>(int)1) && (py<(image_in->ly-2))) {
01785                 /* exclude blank (ZERO) pixels from interpolation */
01786                 if (qfits_isnan(pixint[pos]) && image_in->lx != 1 )
01787                 {
01788                     value = ZERO ;
01789                 }
01790                 else
01791                 {
01792                     if (qfits_isnan(pixint[pos-image_in->lx]))
01793                     {
01794                         pixint[pos-image_in->lx] = 0. ;
01795                     }
01796                     if (qfits_isnan(pixint[pos+image_in->lx]))
01797                     {
01798                         pixint[pos+image_in->lx] = 0. ;
01799                     }
01800                     if (qfits_isnan(pixint[pos+2*image_in->lx]))
01801                     {
01802                         pixint[pos+2*image_in->lx] = 0. ;
01803                     }
01804                     /*
01805                      * Sum up over 4 closest pixel values,
01806                      * weighted by interpolation kernel values
01807                      */
01808                     value = (double)pixint[pos-image_in->lx] * ker[mid+taby] +
01809                             (double)pixint[pos] * ker[taby] +
01810                             (double)pixint[pos+image_in->lx] * ker[mid-taby] +
01811                             (double)pixint[pos+2*image_in->lx]*ker[samples-taby-1];
01812                     /*
01813                      * Also sum up interpolation kernel coefficients
01814                      * for further normalization
01815                      */
01816                     norm =      (double)ker[mid+taby] +
01817                                 (double)ker[taby] +
01818                                 (double)ker[mid-taby] +
01819                                 (double)ker[samples-taby-1] ;
01820 
01821                     if (fabs(norm) > 1e-4) {
01822                         value /= norm ;
01823                     }
01824                 }
01825             } else {
01826                 value = ZERO ;
01827             }  
01828             if (qfits_isnan(value))
01829             {
01830                 second_pass[i+j*image_in->lx] = ZERO ;
01831             }
01832             else
01833             {
01834                 second_pass[i+j*image_in->lx] = (pixelvalue)value ;
01835             }
01836         }
01837     }
01838 
01839     cpl_free(first_pass) ;
01840     if (freeKernel)
01841         cpl_free(ker) ;
01842     return shifted ;
01843 }
01844 
01845 
01846 /*-------------------------------------------------------------------------*/
01877 /*--------------------------------------------------------------------------*/
01878 
01879 void
01880 shiftImageinCube(
01881     OneImage     *   image_in,
01882     double           shift_x,
01883     double           shift_y,
01884     double       *   interp_kernel,
01885     OneImage     *   shifted,
01886     pixelvalue   *   first_pass)
01887 {  
01888     pixelvalue  *       second_pass ;
01889     int             samples = KERNEL_SAMPLES ;
01890     int          i, j ;
01891     double           fx, fy ;
01892     double           rx, ry ;
01893     int             px, py ;
01894     int             tabx, taby ;
01895     double           value ;
01896     size_t          pos ;
01897     register pixelvalue     *   pix ;
01898     register pixelvalue     *   pixint ;
01899     int             mid;
01900     double          norm ;
01901     double       *      ker ;
01902     int                         freeKernel = 1 ;
01903     /* error handling: test entries */
01904         if (image_in==NULL) shifted = NULL ;
01905 
01906         /* Shifting by a zero offset returns a copy of the input image */
01907         if ((fabs(shift_x)<1e-2) && (fabs(shift_y)<1e-2))
01908                 memcpy(shifted->data, image_in->data, (size_t) shifted->nbpix * sizeof(pixelvalue)) ;
01909 
01910         /* See if a kernel needs to be generated */
01911     if (interp_kernel == NULL) {
01912         ker = generate_interpolation_kernel("default") ;
01913         if (ker == NULL) {
01914             cpl_msg_error("kernel generation failure:","aborting resampling") ;
01915             shifted = NULL ;
01916         }
01917     } else {
01918         ker = interp_kernel ;
01919         freeKernel = 0 ;
01920     }
01921 
01922     mid = (int)samples/(int)2 ;
01923     second_pass = shifted->data ;
01924 
01925     pix = image_in->data ;
01926     for (j=0 ; j<image_in->ly ; j++) {
01927         for (i=1 ; i<image_in->lx-2 ; i++) {
01928             fx = (double)i-shift_x ;
01929             px = (int)fx ;
01930             rx = fx - (double)px ;
01931 
01932             pos = px + j * image_in->lx ;
01933 
01934             if ((px>1) && (px<(image_in->lx-2))) {
01935                 tabx = (int)(fabs((double)mid * rx)) ;
01936                 /* exclude blank (ZERO) pixels from interpolation */
01937                 if (qfits_isnan(pix[pos]))
01938                 {
01939                     value = ZERO ;
01940                 }
01941                 else
01942                 {
01943                     if (qfits_isnan(pix[pos-1]))
01944                     {
01945                         pix[pos-1] = 0. ;
01946                     }
01947                     if (qfits_isnan(pix[pos+1]))
01948                     {
01949                         pix[pos+1] = 0. ;
01950                     }
01951                     if (qfits_isnan(pix[pos+2]))
01952                     {
01953                         pix[pos+2] = 0. ;
01954                     }
01955                    
01956                     /*
01957                      * Sum up over 4 closest pixel values,
01958                      * weighted by interpolation kernel values
01959                      */
01960                     value =     (double)pix[pos-1] * ker[mid+tabx] +
01961                                 (double)pix[pos] * ker[tabx] +
01962                                 (double)pix[pos+1] * ker[mid-tabx] +
01963                                 (double)pix[pos+2] * ker[samples-tabx-1] ;
01964                     /*
01965                      * Also sum up interpolation kernel coefficients
01966                      * for further normalization
01967                      */
01968                     norm =      (double)ker[mid+tabx] +
01969                                 (double)ker[tabx] +
01970                                 (double)ker[mid-tabx] +
01971                                 (double)ker[samples-tabx-1] ;
01972                     if (fabs(norm) > 1e-4) {
01973                         value /= norm ;
01974                     }
01975                 }
01976             } else {
01977                 value = 0.0 ;
01978             }  
01979             /*
01980              * There may be a problem of rounding here if pixelvalue
01981              * has not enough bits to sustain the accuracy.
01982              */
01983             if ( qfits_isnan(value) )
01984             {
01985                 first_pass[i+j*image_in->lx] = ZERO ;
01986             }
01987             else
01988             {
01989                 first_pass[i+j*image_in->lx] = (pixelvalue)value ;
01990             }
01991         }
01992     }
01993     pixint = first_pass ;
01994     for (i=0 ; i< image_in->lx ; i++) {
01995         for (j=1 ; j< image_in->ly-2 ; j++) {
01996             fy = (double)j - shift_y ;
01997             py = (int)fy ;
01998             ry = fy - (double)py ;
01999             pos = i + py * image_in->lx ;
02000 
02001             taby = (int)(fabs((double)mid * ry)) ;
02002 
02003             if ((py>(int)1) && (py<(image_in->ly-2))) {
02004                 /* exclude blank (ZERO) pixels from interpolation */
02005                 if (qfits_isnan(pixint[pos]))
02006                 {
02007                     value = ZERO ;
02008                 }
02009                 else
02010                 {
02011                     if (qfits_isnan(pixint[pos-image_in->lx]))
02012                     {
02013                         pixint[pos-image_in->lx] = 0. ;
02014                     }
02015                     if (qfits_isnan(pixint[pos+image_in->lx]))
02016                     {
02017                         pixint[pos+image_in->lx] = 0. ;
02018                     }
02019                     if (qfits_isnan(pixint[pos+2*image_in->lx]))
02020                     {
02021                         pixint[pos+2*image_in->lx] = 0. ;
02022                     }
02023                     /*
02024                      * Sum up over 4 closest pixel values,
02025                      * weighted by interpolation kernel values
02026                      */
02027                     value = (double)pixint[pos-image_in->lx] * ker[mid+taby] +
02028                             (double)pixint[pos] * ker[taby] +
02029                             (double)pixint[pos+image_in->lx] * ker[mid-taby] +
02030                             (double)pixint[pos+2*image_in->lx]*ker[samples-taby-1];
02031                     /*
02032                      * Also sum up interpolation kernel coefficients
02033                      * for further normalization
02034                      */
02035                     norm =      (double)ker[mid+taby] +
02036                                 (double)ker[taby] +
02037                                 (double)ker[mid-taby] +
02038                                 (double)ker[samples-taby-1] ;
02039 
02040                     if (fabs(norm) > 1e-4) {
02041                         value /= norm ;
02042                     }
02043                 }
02044             } else {
02045           /* value = 0.0 ; AMo: This affect slitlet #1 */
02046             }  
02047             if (qfits_isnan(value))
02048             {
02049                 second_pass[i+j*image_in->lx] = ZERO ;
02050             }
02051             else
02052             {
02053                 second_pass[i+j*image_in->lx] = (pixelvalue)value ;
02054             }
02055         }
02056     }
02057 
02058     if (freeKernel)
02059         cpl_free(ker) ;
02060 }
02061 
02062 /* function to delete the image statistics within python */
02063 void delStats( Stats * st)
02064 {
02065     cpl_free (st) ;
02066 }
02067 
02068 /*----------------------------------------------------------------------------
02069    Function     :       combineMasks()
02070    In           :       firstMask, secondMask: bad pixel masks to combine
02071    Out          :       resulting image
02072    Job          :       combines two bad pixel mask to one using an or relation
02073  ---------------------------------------------------------------------------*/
02074 
02075 OneImage * combineMasks ( OneImage * firstMask, OneImage * secondMask ) 
02076 {
02077     OneImage * retMask ;
02078     int n ;
02079 
02080     if ( firstMask == NullImage || secondMask == NullImage )
02081     {
02082         cpl_msg_error ("combineMasks:","no input mask image given!") ;
02083         return NullImage ;
02084     }
02085     retMask = copy_image (firstMask) ;
02086     
02087     for ( n = 0 ; n < (int) retMask->nbpix ; n++ )
02088     {
02089        if ( retMask -> data[n] == 0. || secondMask ->data[n] == 0. )
02090        {
02091            retMask ->data[n] = 0. ;
02092        }
02093        else
02094        {
02095            retMask ->data[n] = 1. ;
02096        }
02097     }
02098     return retMask ;
02099 }
02100 
02101 /*----------------------------------------------------------------------------
02102    Function     :       sliceCube()
02103    In           :       cube: input reduced data cube
02104                         x, y: x slice or y slice pixel value
02105    Out          :       resulting slice image
02106    Job          :       slices a data cube in x or y direction
02107  ---------------------------------------------------------------------------*/
02108 
02109 OneImage * sliceCube ( OneCube * cube, int x, int y ) 
02110 {
02111     OneImage * retImage ;
02112     int col, row, z ;
02113 
02114     if ( cube == NullCube )
02115     {
02116         cpl_msg_error("sliceCube:","no cube given!") ;
02117         return NullImage ;
02118     }
02119     if ( x > 31 || y > 31 )
02120     {
02121         cpl_msg_warning("sliceCube:","wrong x or y values!") ;
02122     }
02123 
02124     if ( x < 0 )
02125     {
02126         /* allocate memory */
02127         if ( NULL == (retImage = new_image(cube->lx, cube->np)) ) 
02128         {
02129             cpl_msg_error("sliceCube:","could not allocate memory!") ;
02130             return NullImage ;
02131         }
02132 
02133         for ( z = 0 ; z < cube->np ; z++ )
02134         {
02135             for ( col = 0 ; col < cube->lx ; col++ )  
02136             {
02137                 retImage -> data[col+z*cube->lx] = cube->plane[z]->data[col+y*cube->lx] ;
02138             }
02139         }
02140     }
02141     else if ( y < 0 )
02142     {
02143         /* allocate memory */
02144         if ( NULL == (retImage = new_image(cube->ly, cube->np)) ) 
02145         {
02146             cpl_msg_error("sliceCube:","could not allocate memory!") ;
02147             return NullImage ;
02148         }
02149 
02150         for ( z = 0 ; z < cube->np ; z++ )
02151         {
02152             for ( row = 0 ; row < cube->ly ; row++ )  
02153             {
02154                 retImage -> data[row+z*cube->ly] = cube->plane[z]->data[x+row*cube->ly] ;
02155             }
02156         }
02157     }
02158     else 
02159     {
02160         cpl_msg_error("sliceCube:","wrong input!") ;
02161         return NullImage ;
02162     }
02163     return retImage ;
02164 }
02165 
02166 /*----------------------------------------------------------------------------
02167    Function     :       divImagesRobust()
02168    In           :       im1, im2: input images im1/im2
02169    Out          :       resulting divided image
02170    Job          :       divides two images by considering blanks and
02171                         calculating first 1/im2 by
02172                         cutting the very high values and setting to 1,
02173                         then multiplying im1 * 1/im2.
02174  ---------------------------------------------------------------------------*/
02175 
02176 OneImage * divImagesRobust ( OneImage * im1, OneImage * im2 ) 
02177 {
02178     OneImage * retIm ;
02179     float help ;
02180     int i ;
02181   
02182     if ( im1 == NullImage || im2 == NullImage ) 
02183     {
02184         cpl_msg_error("divImagesRobust:","no input images given!") ;
02185         return NullImage ;
02186     }
02187     if ( im1 -> lx != im2 ->lx || im1->ly != im2->ly ) 
02188     {
02189         cpl_msg_error("divImagesRobust:","images not compatible!") ;
02190         return NullImage ;
02191     }
02192     if ( NullImage == (retIm = new_image(im1->lx, im1->ly)) ) 
02193     {
02194         cpl_msg_error("divImagesRobust:","could not allocate memory!") ;
02195         return NullImage ;
02196     }
02197 
02198     for ( i = 0 ; i < (int) im2 -> nbpix ; i++ )
02199     {
02200         if ( !qfits_isnan(im2 -> data[i]) )
02201         {
02202             help = 1./im2->data[i] ;
02203             if (fabs( help )> THRESH )
02204             {
02205                help = 1. ;
02206             }
02207         }
02208         else
02209         {
02210             help = ZERO ;
02211         }
02212         if ( qfits_isnan(help) || qfits_isnan(im1->data[i]) )
02213         {
02214             retIm->data[i] = ZERO ;
02215         }
02216         else
02217         {
02218             retIm->data[i] = im1->data[i] * help ;
02219         }
02220     }
02221     return retIm ;
02222 }
02223 
02224 OneImage * nullEdges ( OneImage * image) 
02225 {
02226     OneImage * new ;
02227     int i,j ;
02228 
02229     if ( image == NullImage )
02230     {
02231         cpl_msg_error ("nullEdges:","no input image given!\n") ;
02232         return NullImage ;
02233     }
02234     new = copy_image (image) ;
02235     
02236     for ( i = 0 ; i < new->lx ; i++ )
02237     {
02238         for ( j = 0 ; j < 4 ; j++)
02239     {
02240         new->data[i+j*new->lx]=0;
02241         new->data[i+(new->ly-j-1)*new->lx]=0;
02242     }
02243     }
02244     for ( i = 0 ; i < new->ly ; i++ )
02245     {
02246         for ( j = 0 ; j < 4 ; j++)
02247     {
02248         new->data[j+i*new->lx]=0;
02249         new->data[(new->lx-j-1)+i*new->lx]=0;
02250     }
02251     }
02252     return new ;
02253 }
02254 
02255 
02256 void usedcormap( OneImage *im, OneImage *map)
02257 {
02258     int i,j,index;
02259     float temp_array[2048];
02260 
02261     for( j=0; j<im->ly; j++)
02262     {
02263     for( i=0;i<im->lx;i++)
02264         {
02265           index = (int)map->data[i+j*im->lx];
02266           temp_array[i] = im->data[index+j*im->lx];
02267         }
02268     for( i=0;i<im->lx;i++)
02269         {
02270           im->data[i+j*im->lx]= temp_array[i];
02271         }
02272     }
02273 }
02274 /*--------------------------------------------------------------------------*/

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