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 /*--------------------------------------------------------------------------*/
1.2.13.1 written by Dimitri van Heesch,
© 1997-2001