fors_image.c

00001 /* $Id: fors_image.c,v 1.54 2008/08/04 12:02:24 cizzo Exp $
00002  *
00003  * This file is part of the FORS Library
00004  * Copyright (C) 2002-2006 European Southern Observatory
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA
00019  */
00020 
00021 /*
00022  * $Author: cizzo $
00023  * $Date: 2008/08/04 12:02:24 $
00024  * $Revision: 1.54 $
00025  * $Name:  $
00026  */
00027 
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031 
00032 #include <fors_image.h>
00033 
00034 #include <fors_dfs.h>
00035 #include <fors_utils.h>
00036 #include <fors_pfits.h>
00037 #include <fors_double.h>
00038 
00039 #include <moses.h>
00040 
00041 #include <cpl.h>
00042 
00043 #include <math.h>
00044 #include <stdbool.h>
00045 #include <stdio.h>
00046 
00062 struct _fors_image
00063 {
00064     cpl_image *data;
00065     cpl_image *variance;
00066 
00067     /* Invariants:
00068        The CPL images are non-NULL. 
00069        The variance image is everywhere non-negative.
00070        The CPL image types are FORS_IMAGE_TYPE.
00071        The CPL image bad pixel masks are unused
00072     */
00073 };
00074 
00078 const cpl_type FORS_IMAGE_TYPE = CPL_TYPE_FLOAT;
00079 #define FORS_IMAGE_TYPE_MAX FLT_MAX  /* Use a #define rather than a variable here
00080                                         to avoid type casting */
00081 
00082 #undef cleanup
00083 
00084 /* 
00085  * The following static function passes a max filter of given box
00086  * size on the input data buffer. The output data buffer must be
00087  * pre-allocated. The box size must be a positive odd integer.
00088  * Returns 0 on success.
00089  */
00090  
00091 static int 
00092 max_filter(const float *ibuffer, float *obuffer, int length, int size)
00093 {
00094     float  max;
00095     int    start = size / 2;
00096     int    end   = length - size / 2;
00097     int    i, j;
00098 
00099 
00100     for (i = start; i < end; i++) {
00101         max = ibuffer[i-start];
00102         for (j = i - start + 1; j <= i + start; j++)
00103             if (max < ibuffer[j])
00104                 max = ibuffer[j];
00105         obuffer[i] = max;
00106     }
00107 
00108     for (i = 0; i < start; i++)
00109         obuffer[i] = obuffer[start];
00110  
00111     for (i = end; i < length; i++)
00112         obuffer[i] = obuffer[end-1];
00113 
00114     return 0;
00115 }
00116 
00117 #define cleanup
00118 
00127 fors_image *
00128 fors_image_new(cpl_image *data, cpl_image *variance)
00129 {
00130     fors_image *image = NULL;
00131 
00132     assure( data != NULL, return NULL, NULL );
00133     assure( variance != NULL, return NULL, NULL );
00134     
00135     assure( cpl_image_get_type(data) == FORS_IMAGE_TYPE, return NULL,
00136             "Provided data image type is %s, must be %s",
00137             fors_type_get_string(cpl_image_get_type(data)),
00138             fors_type_get_string(FORS_IMAGE_TYPE) );
00139 
00140     assure( cpl_image_get_type(variance) == FORS_IMAGE_TYPE, return NULL,
00141             "Provided weight image type is %s, must be %s",
00142             fors_type_get_string(cpl_image_get_type(variance)),
00143             fors_type_get_string(FORS_IMAGE_TYPE) );
00144 
00145     assure( cpl_image_get_size_x(data) == cpl_image_get_size_x(variance) &&
00146             cpl_image_get_size_y(data) == cpl_image_get_size_y(variance),
00147             return NULL,
00148             "Incompatible data and weight image sizes: %dx%d and %dx%d",
00149             cpl_image_get_size_x(data), cpl_image_get_size_y(data),
00150             cpl_image_get_size_x(variance), cpl_image_get_size_y(variance));
00151 
00152     assure( cpl_image_get_min(variance) >= 0, return NULL,
00153             "Variances must be non-negative, minimum is %f",
00154             cpl_image_get_min(variance));
00155 
00156     image = cpl_malloc(sizeof(*image));
00157 
00158     image->data = data;
00159     image->variance = variance;
00160 
00161     return image;    
00162 }
00163 
00164 #undef cleanup
00165 #define cleanup
00166 
00171 fors_image *
00172 fors_image_duplicate(const fors_image *image)
00173 {
00174     assure( image != NULL, return NULL, NULL );
00175 
00176     return fors_image_new(cpl_image_duplicate(image->data),
00177                           cpl_image_duplicate(image->variance));
00178 }
00179 
00184 void
00185 fors_image_delete(fors_image **image)
00186 {
00187     if (image && *image) {
00188         cpl_image_delete((*image)->data);
00189         cpl_image_delete((*image)->variance);
00190         cpl_free(*image); *image = NULL;
00191     }
00192     return;
00193 }
00194 
00199 void
00200 fors_image_delete_const(const fors_image **image)
00201 {
00202     fors_image_delete((fors_image **)image);
00203 
00204     return;
00205 }
00206 
00207 /* not used */
00208 #if 0
00209 
00213 static void
00214 fors_image_dump(const fors_image *image, FILE *file)
00215 {
00216     if (image == NULL) {
00217         fprintf(file, "Null image\n");
00218     }
00219     else {
00220         cpl_stats *stats;
00221 
00222         fprintf(file, "Data:\n");
00223         stats = cpl_stats_new_from_image(image->data, CPL_STATS_ALL);
00224         cpl_stats_dump(stats, CPL_STATS_ALL, file);
00225         cpl_stats_delete(stats);
00226             
00227         fprintf(file, "Variance:\n");
00228         stats = cpl_stats_new_from_image(image->variance, CPL_STATS_ALL);
00229         cpl_stats_dump(stats, CPL_STATS_ALL, file);
00230         cpl_stats_delete(stats);
00231     }
00232     
00233     return;
00234 }
00235 #endif
00236 
00237 #undef cleanup
00238 #define cleanup \
00239 do { \
00240     double_list_delete(&sat_percent, double_delete); \
00241 } while (0)
00242 
00256 fors_image_list *
00257 fors_image_load_list(const cpl_frameset *frames, const fors_image *bias,
00258                      const fors_setting *setting,
00259                      double *saturated)
00260 {
00261     fors_image_list *ilist = fors_image_list_new();
00262     double_list *sat_percent = double_list_new();
00263     
00264     assure( frames != NULL, return ilist, NULL );
00265     assure( !cpl_frameset_is_empty(frames), return ilist, "Empty frameset");
00266 
00267     {
00268         const cpl_frame *f;
00269         
00270         for (f = cpl_frameset_get_first_const(frames);
00271              f != NULL;
00272              f = cpl_frameset_get_next_const(frames)) {
00273 
00274             double saturated_one;
00275             
00276             fors_image *i = fors_image_load(f, bias, setting,
00277                                             saturated != NULL ? 
00278                                             &saturated_one : NULL);
00279             assure( !cpl_error_get_code(), return ilist, NULL );
00280             
00281             fors_image_list_insert(ilist, i);
00282             if (saturated != NULL) {
00283                 double_list_insert(sat_percent, 
00284                                    double_duplicate(&saturated_one));
00285             }
00286         }
00287     }
00288 
00289     /* fixme: make sure input is consistent */
00290 
00291     /* Compute overall percentage as mean of each percentage.
00292        Valid because input frames have same size. */
00293     if (saturated != NULL) {
00294         *saturated = double_list_mean(sat_percent, double_eval, NULL);
00295     }
00296 
00297     cleanup;
00298     return ilist;
00299 }
00300 
00312 const fors_image_list *
00313 fors_image_load_list_const(const cpl_frameset *frames, const fors_image *bias,
00314                            const fors_setting *setting,
00315                            double *saturated)
00316 {
00317     return (const fors_image_list *)
00318         fors_image_load_list(frames, bias, setting, saturated);
00319 }
00320 
00321 #undef cleanup
00322 #define cleanup \
00323 do { \
00324     cpl_propertylist_delete(header); \
00325     cpl_table_delete(overscans); \
00326     cpl_image_delete(temp); \
00327     cpl_mask_delete(non_saturated); \
00328     fors_setting_delete(&frame_setting); \
00329 } while (0)
00330 
00348 fors_image *
00349 fors_image_load(const cpl_frame *frame, const fors_image *bias,
00350                 const fors_setting *setting,
00351                 double *saturated)
00352 {
00353     fors_image *image           = NULL;
00354     fors_setting *frame_setting = NULL;
00355     cpl_image *data             = NULL;
00356     cpl_image *variance         = NULL;
00357     cpl_image *temp             = NULL;
00358     cpl_propertylist *header    = NULL;
00359     cpl_table *overscans        = NULL;
00360     cpl_mask *non_saturated     = NULL;
00361     const char *filename;
00362     int extension= 0;
00363     const int plane = 0;
00364 
00365     assure( frame != NULL, return image, NULL );
00366     /* bias may be NULL */
00367     assure( setting != NULL, return image, NULL );
00368     filename = cpl_frame_get_filename(frame);
00369     assure( filename != NULL, return image, 
00370             "NULL filename received");
00371     
00372     cpl_msg_info(cpl_func, "Loading %s: %s",
00373                  /* fors_frame_get_group_string(frame), */
00374          (cpl_frame_get_tag(frame) != NULL) ? 
00375          cpl_frame_get_tag(frame) : "NULL",
00376                  filename);
00377 
00378     /* Verify instrument setting */
00379     fors_setting_verify(setting, frame, &frame_setting);
00380     assure( !cpl_error_get_code(), return image, 
00381             "Could not verify %s setting", 
00382             filename);
00383     
00384 
00385     /* Get header */
00386     header = cpl_propertylist_load(filename, extension);
00387     assure( !cpl_error_get_code(), return image, 
00388             "Could not load %s extension %d header", 
00389             filename, extension);
00390 
00391     /* Get data */
00392     data = cpl_image_load(filename, 
00393                           FORS_IMAGE_TYPE, plane, extension);
00394     
00395     assure( !cpl_error_get_code(), return image, 
00396             "Could not load image from %s extension %d", 
00397             filename, extension);
00398 
00399     /* Read, remove overscan areas */
00400     {
00401         int xlow, ylow, xhig, yhig;
00402     bool check_consistency = false;
00403 
00404         overscans = mos_load_overscans_vimos(header, check_consistency);
00405     /* Uses always ESO.DET.OUT1.*  which is fine even
00406        for old FORS data */
00407         
00408         assure( !cpl_error_get_code(), return image, 
00409                 "Could not read overscan information from %s extension %d", 
00410                 filename, extension);
00411         
00412         xlow = cpl_table_get_int(overscans, "xlow", 0, NULL);
00413         ylow = cpl_table_get_int(overscans, "ylow", 0, NULL);
00414         xhig = cpl_table_get_int(overscans, "xhig", 0, NULL);
00415         yhig = cpl_table_get_int(overscans, "yhig", 0, NULL);
00416         
00417         cpl_table_delete(overscans); overscans = NULL;
00418         
00419         temp = cpl_image_duplicate(data);
00420         cpl_image_delete(data);
00421         data = cpl_image_extract(temp, xlow+1, ylow+1, xhig, yhig);
00422         
00423         cpl_image_delete(temp); temp = NULL;
00424     }
00425 
00426     cpl_propertylist_delete(header); header = NULL;       
00427 
00428     /* Define variance */
00429     if (cpl_frame_get_nextensions(frame) == 0 || bias != NULL) {
00430 
00431         /* No error bars provided, assume RON only.
00432            If frame is not a bias, add photon noise later.
00433         */
00434 
00435         /* Define read-out-noise */
00436         variance = cpl_image_new(
00437             cpl_image_get_size_x(data),
00438             cpl_image_get_size_y(data),
00439             FORS_IMAGE_TYPE);        
00440 
00441         if (cpl_frame_get_group(frame) == CPL_FRAME_GROUP_CALIB) {
00442             cpl_msg_warning(cpl_func, 
00443                             "No error bars provided for calibration frame %s, "
00444                             "assuming no errors. For complete error propagation, "
00445                             "you may recreate this frame with this pipeline", 
00446                             filename);
00447         } else {
00448             cpl_image_add_scalar(variance, frame_setting->ron*frame_setting->ron);
00449         }
00450     }
00451     else {
00452 
00453         extension = 1;
00454         
00455         /* Get error bars */
00456         variance = cpl_image_load(filename, 
00457                                   FORS_IMAGE_TYPE, plane, extension);
00458         
00459         assure( !cpl_error_get_code(), return image, 
00460                 "Could not load image from %s extension %d", 
00461                 filename, extension);
00462 
00463         cpl_image_power(variance, 2);
00464 
00465         assure( cpl_image_get_min(variance) >= 0, return image,
00466                 "Illegal minimum variance: %g",
00467                 cpl_image_get_min(variance));
00468 
00469         cpl_image_delete(temp); temp = NULL;
00470     }
00471     
00472     image = fors_image_new(data, variance);
00473 
00474     /* Count saturated pixels */
00475     if (saturated != NULL) {
00476         double lo_cut = 0.5;
00477         double hi_cut = 65534.5;
00478         non_saturated = cpl_mask_threshold_image_create(
00479             data, lo_cut, hi_cut);
00480         
00481         *saturated = (cpl_image_get_size_x(data) * 
00482                       cpl_image_get_size_y(data) - 
00483                       cpl_mask_count(non_saturated)) * 100.0 / 
00484             (cpl_image_get_size_x(data) * 
00485              cpl_image_get_size_y(data));
00486 
00487         cpl_mask_delete(non_saturated); non_saturated = NULL;
00488 
00489         cpl_msg_debug(cpl_func,
00490                      "%f %% saturated pixels", *saturated);
00491     }        
00492     
00493     if (bias == NULL) {
00494         /* Done. No bias to subtract */
00495     }
00496     else {
00497         assure( cpl_frame_get_group(frame) == CPL_FRAME_GROUP_RAW,
00498                 return image,
00499                 "Refusing to subtract bias from non-raw (%s) input frame: %s",
00500                 fors_frame_get_group_string(frame),
00501                 filename);
00502         
00503         cpl_msg_debug(cpl_func, "Subtracting bias from %s", filename);
00504         
00505         fors_image_subtract(image, bias);
00506         
00507         assure( !cpl_error_get_code(), return image, 
00508                 "Bias subtraction failed" );
00509         
00510         /* Variance is now (ron**2  +  bias_noise**2).
00511            Add photonic noise:
00512            
00513            variance := variance + |f|/conad
00514         */
00515 
00516         double conad = 1.0 / frame_setting->average_gain;
00517 
00518         temp = cpl_image_divide_scalar_create(image->data, conad);
00519         cpl_image_abs(temp);
00520 
00521         cpl_image_add(image->variance, temp);
00522     }
00523     
00524     cleanup;
00525     return image;
00526 }
00527 
00528 
00529 #undef cleanup
00530 #define cleanup \
00531 do { \
00532     cpl_image_delete(sigma); \
00533 } while(0)
00534 
00543 void
00544 fors_image_save(const fors_image *image, const cpl_propertylist *header,
00545                 const char *filename)
00546 {
00547     cpl_propertylist *extension_header = NULL;
00548     cpl_image *sigma = NULL;
00549 
00550     assure( image != NULL, return, NULL );
00551     /* header may be NULL */
00552     assure( filename != NULL, return, NULL );
00553     
00554     cpl_image_save(image->data, filename, CPL_BPP_IEEE_FLOAT, header,
00555                    CPL_IO_DEFAULT);
00556     assure( !cpl_error_get_code(), return, 
00557             "Cannot save product %s", filename);
00558     
00559     sigma = cpl_image_power_create(image->variance, 0.5);
00560     /* This would probably be faster if sqrt() is used rather than pow */
00561 
00562     cpl_image_save(sigma, filename, CPL_BPP_IEEE_FLOAT, extension_header,
00563                    CPL_IO_EXTEND);
00564     assure( !cpl_error_get_code(), return, 
00565             "Cannot save product %s", filename);
00566     
00567     cleanup;
00568     return;
00569 }
00570 
00571 
00572 #undef cleanup
00573 #define cleanup \
00574 do { \
00575     cpl_image_delete(var_bkg); \
00576     cpl_image_delete(sigma_bkg); \
00577 } while(0)
00578 
00589 void
00590 fors_image_save_sex(const fors_image *image, const cpl_propertylist *header,
00591                     const char *filename_dat,
00592                     const char *filename_var,
00593                     int radius)
00594 {
00595     cpl_propertylist *extension_header = NULL;
00596     cpl_image *sigma_bkg = NULL;
00597     cpl_image *var_bkg = NULL;
00598 
00599     assure( image != NULL, return, NULL );
00600     /* header may be NULL */
00601     assure( filename_dat != NULL, return, NULL );
00602     assure( filename_var != NULL, return, NULL );
00603 
00604     cpl_image_save(image->data, filename_dat, CPL_BPP_IEEE_FLOAT, header,
00605                    CPL_IO_DEFAULT);
00606     assure( !cpl_error_get_code(), return, 
00607             "Cannot save product %s", filename_dat);
00608     
00609     /* Sextractor wants as input the background error bars,
00610        i.e. excluding sources.
00611        Therefore filter away sources but keep the sharp edges
00612        between the illuminated / non-illuminated areas.
00613 
00614        I.e. use a median filter, average filter would not work.
00615     */
00616 
00617     cpl_msg_info(cpl_func, "Creating background error map");
00618 
00619     bool filter_data = false;  /* filter the variance image */
00620     int xstep = radius/2; /* 25 points sampling grid 
00621                              . . . . .
00622                              . . . . .
00623                              . . . . .
00624                              . . . . .
00625                              . . . . .
00626                            */
00627     int ystep = radius/2;
00628     int xstart = 1;
00629     int ystart = 1;
00630     int xend = fors_image_get_size_x(image);
00631     int yend = fors_image_get_size_y(image);
00632 
00633 
00634     var_bkg = fors_image_filter_median_create(image, 
00635                                               radius,
00636                                               radius,
00637                                               xstart, ystart,
00638                                               xend, yend,
00639                                               xstep, ystep,
00640                                               filter_data);
00641     assure( !cpl_error_get_code(), return, 
00642             "Median filtering failed");
00643 
00644     sigma_bkg = cpl_image_power_create(var_bkg, 0.5);
00645 
00646     cpl_image_save(sigma_bkg, filename_var,
00647                    CPL_BPP_IEEE_FLOAT, extension_header,
00648                    CPL_IO_DEFAULT);
00649     assure( !cpl_error_get_code(), return, 
00650             "Cannot save product %s", filename_var);
00651 
00652     cleanup;
00653     return;
00654 }
00655 
00656 #undef cleanup
00657 #define cleanup
00658 
00663 int fors_image_get_size_x(const fors_image *image)
00664 {
00665     assure( image != NULL, return -1, NULL );
00666     return cpl_image_get_size_x(image->data);
00667 }
00668 
00669 #undef cleanup
00670 #define cleanup
00671 
00676 int fors_image_get_size_y(const fors_image *image)
00677 {
00678     assure( image != NULL, return -1, NULL );
00679     return cpl_image_get_size_y(image->data);
00680 }
00681 
00682 #undef cleanup
00683 #define cleanup
00684 
00688 const float *fors_image_get_data_const(const fors_image *image)
00689 {
00690     assure( image != NULL, return NULL, NULL );
00691 
00692     assure( FORS_IMAGE_TYPE == CPL_TYPE_FLOAT, return NULL, NULL );
00693     /* This function (including API) would need to change
00694        if the pixel type changes */
00695 
00696     return cpl_image_get_data_float(image->data);
00697 }
00698 
00699 #undef cleanup
00700 #define cleanup
00701 
00708 void
00709 fors_image_abs(fors_image *image)
00710 {
00711     assure( image != NULL, return, NULL );
00712 
00713     cpl_image_abs(image->data);
00714 
00715     return;
00716 }
00717 
00718 #undef cleanup
00719 #define cleanup
00720 
00727 void
00728 fors_image_square(fors_image *image)
00729 {
00730     assure( image != NULL, return, NULL );
00731 
00732     cpl_image_multiply(image->data, image->data);
00733     /* It is an undocumented feature of CPL that you
00734        can pass the same image to cpl_image_multiply and get
00735        the right answer. Let us hope it does not change...
00736     */
00737     cpl_image_multiply_scalar(image->variance, 2);
00738 
00739     return;
00740 }
00741 
00742 
00743 #undef cleanup
00744 #define cleanup \
00745 do { \
00746     cpl_image_delete(temp); \
00747 } while(0)
00748 
00756 void
00757 fors_image_subtract(fors_image *left, const fors_image *right)
00758 {
00759     cpl_image *temp = NULL;
00760     assure( left != NULL, return, NULL );
00761     assure( right != NULL, return, NULL );
00762 
00763     cpl_image_subtract(left->data, right->data);
00764 
00765     /*  variance_left := variance_left + variance_right */
00766     cpl_image_add(left->variance, right->variance);
00767 
00768     cleanup;
00769     return;
00770 }
00771 
00772 #undef cleanup
00773 #define cleanup
00774 
00784 void
00785 fors_image_multiply_noerr(fors_image *left, const cpl_image *right)
00786 {
00787     assure( left != NULL, return, NULL );
00788     assure( right != NULL, return, NULL );
00789     assure( cpl_image_get_size_x(left->data) == cpl_image_get_size_x(right) &&
00790             cpl_image_get_size_y(left->data) == cpl_image_get_size_y(right),
00791             return, "Incompatible image sizes: %dx%d and %dx%d",
00792             cpl_image_get_size_x(left->data),
00793             cpl_image_get_size_y(left->data),
00794             cpl_image_get_size_x(right),
00795             cpl_image_get_size_y(right));
00796 
00797     cpl_image_multiply(left->data, right);
00798     cpl_image_multiply(left->variance, right);
00799     cpl_image_multiply(left->variance, right);
00800 
00801     return;
00802 }
00803 
00804 #undef cleanup
00805 #define cleanup
00806 
00822 void
00823 fors_image_divide_noerr(fors_image *left, cpl_image *right)
00824 {
00825     assure( left != NULL, return, NULL );
00826     assure( right != NULL, return, NULL );
00827     assure( cpl_image_get_size_x(left->data) == cpl_image_get_size_x(right) &&
00828             cpl_image_get_size_y(left->data) == cpl_image_get_size_y(right),
00829             return, "Incompatible image sizes: %dx%d and %dx%d",
00830             cpl_image_get_size_x(left->data),
00831             cpl_image_get_size_y(left->data),
00832             cpl_image_get_size_x(right),
00833             cpl_image_get_size_y(right));
00834 
00835     int x, y;
00836     int nx = cpl_image_get_size_x(right);
00837     int ny = cpl_image_get_size_y(right);
00838     float *datal = cpl_image_get_data_float(left->data);
00839     float *datav = cpl_image_get_data_float(left->variance);
00840     float *datar = cpl_image_get_data_float(right);
00841     for (y = 0; y < ny; y++) {
00842         for (x = 0; x < nx; x++) {
00843             if (datar[x + nx*y] == 0) {
00844                 datar[x + nx*y] = 1;
00845                 datal[x + nx*y] = 1;
00846 
00847                 datav[x + nx*y] = FORS_IMAGE_TYPE_MAX;
00848             }
00849         }
00850     }
00851 
00852     cpl_image_divide(left->data, right);
00853     cpl_image_divide(left->variance, right);
00854     cpl_image_divide(left->variance, right);
00855 
00856     return;
00857 }
00858 
00859 #undef cleanup
00860 #define cleanup \
00861 do { \
00862     fors_image_delete(&dupl); \
00863 } while(0)
00864 
00884 void
00885 fors_image_divide(fors_image *left, const fors_image *right)
00886 {
00887     fors_image *dupl = NULL;
00888 
00889     assure( left  != NULL, return, NULL );
00890     assure( right != NULL, return, NULL );
00891 
00892     dupl = fors_image_duplicate(right);
00893 
00894     cpl_image_divide(left->data, dupl->data); 
00895     /* This CPL function divides by zero by setting  x/0 = 1 for all x */
00896 
00897     cpl_image_multiply(dupl->variance, left->data);
00898     cpl_image_multiply(dupl->variance, left->data);
00899 
00900     /* Now  dupl->variance = sigma2^2 * data1^2 / data2^2 */
00901 
00902     cpl_image_add(left->variance, dupl->variance);
00903 
00904     /* Now  left->variance = sigma1^2 + sigma2^2 * data1^2 / data2^2 */
00905 
00906     cpl_image_divide(left->variance, dupl->data);
00907     cpl_image_divide(left->variance, dupl->data);
00908     /* QED */
00909 
00910     /* Handle division by zero */
00911     int x, y;
00912     int nx = cpl_image_get_size_x(left->data);
00913     int ny = cpl_image_get_size_y(left->data);
00914     float *datal = cpl_image_get_data_float(left->data);
00915     float *datav = cpl_image_get_data_float(left->variance);
00916     float *datar = cpl_image_get_data_float(right->data);
00917     for (y = 0; y < ny; y++) {
00918         for (x = 0; x < nx; x++) {
00919             if (datar[x + nx*y] == 0) {
00920                 datal[x + nx*y] = 1;
00921                 datav[x + nx*y] = FORS_IMAGE_TYPE_MAX;
00922             }
00923         }
00924     }
00925 
00926     cleanup;
00927     return;
00928 }
00929 
00930 #undef cleanup
00931 #define cleanup \
00932 do { \
00933     cpl_image_delete(s22d12); \
00934 } while(0)
00935 
00946 void
00947 fors_image_multiply(fors_image *left, const fors_image *right)
00948 {
00949     cpl_image *s22d12 = NULL;
00950 
00951     assure( left  != NULL, return, NULL );
00952     assure( right != NULL, return, NULL );
00953 
00954     s22d12 = cpl_image_duplicate(right->variance);
00955     cpl_image_multiply(s22d12, left->data);
00956     cpl_image_multiply(s22d12, left->data);
00957     
00958     cpl_image_multiply(left->variance, right->data);
00959     cpl_image_multiply(left->variance, right->data);
00960     cpl_image_add(left->variance, s22d12);
00961 
00962     cpl_image_multiply(left->data, right->data);
00963 
00964     cleanup;
00965     return;
00966 }
00967 
00968 
00969 #undef cleanup
00970 #define cleanup
00971 
00983 void fors_image_subtract_scalar(fors_image *image, double s, double ds)
00984 {
00985     assure( image != NULL, return, NULL );
00986     assure( ds <= 0, return, "Unsupported");
00987 
00988     cpl_image_subtract_scalar(image->data, s);
00989 
00990     return;
00991 }
00992 
00993 
00994 #undef cleanup
00995 #define cleanup
00996 
01008 void fors_image_divide_scalar(fors_image *image, double s, double ds)
01009 {
01010     assure( image != NULL, return, NULL );
01011     assure( s != 0, return, "Division by zero");
01012     assure( ds <= 0, return, "Unsupported");
01013 
01014     cpl_image_divide_scalar(image->data, s);
01015     cpl_image_divide_scalar(image->variance, s*s);
01016     
01017     return;
01018 }
01019 
01020 #undef cleanup
01021 #define cleanup
01022 
01034 void fors_image_multiply_scalar(fors_image *image, double s, double ds)
01035 {
01036     assure( image != NULL, return, NULL );
01037     assure( ds <= 0, return, "Unsupported");
01038 
01039     cpl_image_multiply_scalar(image->data, s);
01040     cpl_image_multiply_scalar(image->variance, s*s);
01041 
01042     return;
01043 }
01044 
01045 #undef cleanup
01046 #define cleanup \
01047 do { \
01048     cpl_image_delete(temp); \
01049 } while(0)
01050 
01063 void fors_image_exponential(fors_image *image, double b, double db)
01064 {
01065     cpl_image *temp = NULL;
01066 
01067     assure( image != NULL, return, NULL );
01068     assure( b >= 0, return, "Negative base: %f", b);
01069     assure( db <= 0, return, "Unsupported");
01070 
01071     cpl_image_exponential(image->data, b);
01072 
01073     double lnb = log(b);
01074     
01075     cpl_image_multiply_scalar(image->variance, lnb*lnb);
01076     cpl_image_multiply(image->variance, image->data);
01077     cpl_image_multiply(image->variance, image->data);
01078 
01079     return;
01080 }
01081 
01082 
01083 #undef cleanup
01084 #define cleanup
01085 
01090 double
01091 fors_image_get_min(const fors_image *image)
01092 {
01093     assure( image != NULL, return 0, NULL );
01094 
01095     return cpl_image_get_min(image->data);
01096 }
01097 
01098 #undef cleanup
01099 #define cleanup
01100 
01105 double
01106 fors_image_get_max(const fors_image *image)
01107 {
01108     assure( image != NULL, return 0, NULL );
01109 
01110     return cpl_image_get_max(image->data);
01111 }
01112 
01113 #undef cleanup
01114 #define cleanup
01115 
01121 double
01122 fors_image_get_mean(const fors_image *image, double *dmean)
01123 {
01124     assure( image != NULL, return 0, NULL );
01125     assure( dmean == NULL, return 0, "Unsupported");
01126 
01127     return cpl_image_get_mean(image->data);
01128 }
01129 
01130 #undef cleanup
01131 #define cleanup
01132 
01138 double
01139 fors_image_get_median(const fors_image *image, double *dmedian)
01140 {
01141     assure( image != NULL, return 0, NULL );
01142     assure( dmedian == NULL, return 0, "Unsupported");
01143 
01144     return cpl_image_get_median(image->data);
01145 }
01146 
01147 
01148 #undef cleanup
01149 #define cleanup
01150 
01164 void fors_image_crop(fors_image *image,
01165              int xlo, int ylo,
01166              int xhi, int yhi)
01167 {
01168     /* CPL is missing the function to locally extract an image,
01169        so this this inefficient CPL function */
01170     assure( image != NULL, return, NULL );
01171     assure( 1 <= xlo && xlo <= xhi && xhi <= fors_image_get_size_x(image) &&
01172             1 <= ylo && ylo <= yhi && yhi <= fors_image_get_size_y(image),
01173             return, "Cannot extraction region (%d, %d) - (%d, %d) of "
01174             "%dx%d image",
01175             xlo, ylo, xhi, yhi,
01176             fors_image_get_size_x(image),
01177             fors_image_get_size_y(image));
01178     
01179     cpl_image *new_data = cpl_image_extract(image->data,
01180                                             xlo, ylo,
01181                                             xhi, yhi);
01182     cpl_image_delete(image->data);
01183 
01184     cpl_image* new_variance = cpl_image_extract(image->variance,
01185                                                 xlo, ylo,
01186                                                 xhi, yhi);
01187     cpl_image_delete(image->variance);
01188 
01189     image->data = new_data;
01190     image->variance = new_variance;
01191 
01192     return;
01193 }
01194 
01220 cpl_image *
01221 fors_image_filter_median_create(const fors_image *image, 
01222                                 int xradius,
01223                                 int yradius,
01224                                 int xstart, 
01225                                 int ystart,
01226                                 int xend,
01227                                 int yend,
01228                                 int xstep,
01229                                 int ystep,
01230                                 bool use_data)
01231 {
01232     const cpl_image *input = NULL;
01233     cpl_image *smooth = NULL;
01234     int nx, ny;
01235     
01236     assure( image != NULL, return smooth, NULL );
01237     passure( image->data != NULL, return smooth );
01238     passure( image->variance != NULL, return smooth );
01239     
01240     input = (use_data) ? image->data : image->variance;
01241 
01242     nx = cpl_image_get_size_x(input);
01243     ny = cpl_image_get_size_y(input);
01244 
01245     if (xstep < 1) xstep = 1;
01246     if (ystep < 1) ystep = 1;
01247 
01248     assure( 1 <= xstart && xstart <= xend && xend <= nx &&
01249             1 <= ystart && ystart <= yend && yend <= ny, return smooth,
01250             "Illegal region (%d, %d) - (%d, %d) of %dx%d image",
01251             xstart, ystart,
01252             xend, yend,
01253             nx, ny);
01254     
01255     smooth = cpl_image_duplicate(input);
01256 
01257     /* For efficiency reasons, assume that the image type is float */
01258     assure( FORS_IMAGE_TYPE == CPL_TYPE_FLOAT, return smooth, NULL );
01259 
01260     const float *input_data  = cpl_image_get_data_float_const(input);
01261     float *smooth_data = cpl_image_get_data_float(smooth);
01262     float *data = cpl_malloc((2*yradius + 1)*(2*xradius + 1)*sizeof(*data));
01263     
01264     int y;
01265     for (y = ystart; y < yend; y++) {
01266         /*
01267           Sample kernel on grid which always contains the central pixel
01268           
01269           Trim window (note: this will cause fewer values to
01270           be used for the median near the region borders 
01271         */
01272         int ylo = y - (yradius/ystep) * ystep;
01273         int yhi = y + (yradius/ystep) * ystep;
01274         
01275         while (ylo < ystart) ylo += ystep;
01276         while (yhi > yend  ) yhi -= ystep;
01277         
01278         int x;
01279         for (x = xstart; x < xend; x++) {
01280             int xlo = x - (xradius/xstep) * xstep;
01281             int xhi = x + (xradius/xstep) * xstep;
01282             
01283             while (xlo < xstart) xlo += xstep;
01284             while (xhi > xend  ) xhi -= xstep;
01285             
01286             /* Collect data */
01287             int k = 0;
01288             int j, i;
01289             for (j = ylo; j <= yhi; j += ystep) {
01290                 for (i = xlo; i <= xhi; i += xstep) {
01291                     data[k++] = input_data[ (i-1) + (j-1)*nx ];
01292                 }
01293             }
01294         
01295             /* Get median */
01296             smooth_data[ (x-1) + (y-1)*nx ] = 
01297                                fors_tools_get_median_float(data, k);
01298         }
01299     }
01300 
01301     cpl_free(data);
01302     return smooth;
01303 }
01304 
01305 #undef cleanup
01306 #define cleanup \
01307 do { \
01308     cpl_image_delete(input); \
01309 } while(0)
01310 cpl_image *
01311 fors_image_flat_fit_create(fors_image *image, 
01312                            int step, 
01313                            int degree, 
01314                            float level)
01315 {
01316     cpl_image *temp = NULL;
01317     cpl_image *input = NULL;
01318     cpl_image *smooth = NULL;
01319     int nx, ny;
01320 
01321     assure( image != NULL, return smooth, NULL );
01322     passure( image->data != NULL, return smooth );
01323     assure( step > 0, return smooth, NULL );
01324     assure( degree >= 0, return smooth, NULL );
01325 
01326 
01327     temp = image->data;
01328 
01329     nx = cpl_image_get_size_x(temp);
01330     ny = cpl_image_get_size_y(temp);
01331 
01332     /* 
01333      * For efficiency reasons, assume that the image type is float 
01334      */
01335 
01336     assure( FORS_IMAGE_TYPE == CPL_TYPE_FLOAT, return smooth, NULL );
01337 
01338     /*
01339      * Apply light median filter, to eliminate big outliers from fit
01340      */
01341 
01342     cpl_matrix *kernel = cpl_matrix_new(3, 3);
01343     cpl_matrix_fill(kernel, 1.0);
01344     input = cpl_image_filter_median(image->data, kernel);
01345     cpl_matrix_delete(kernel);
01346     const float *input_data = cpl_image_get_data_float_const(input);
01347 
01348     /*
01349      * First of all, count how many points will have to be fitted
01350      */
01351 
01352     int x, y, pos;
01353     int count = 0;
01354     for (y = 0; y < ny; y += step) {
01355         pos = y*nx;
01356         for (x = 0; x < nx; x += step, pos += step) {
01357             if (input_data[pos] > level) {
01358                 count++;
01359             }
01360         }
01361     }
01362 
01363     if (count < (degree+1)*(degree+2)) {
01364         step = sqrt((nx*nx)/((degree+1)*(degree+2))) / 2;
01365         if (step == 0)
01366             step = 1;
01367         cpl_msg_error(cpl_func, "Flat field image too small (%dx%d). "
01368                       "Please provide a smaller resampling step (a good "
01369                       "value would be %d)", nx, ny, step);
01370         cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_INPUT);
01371         cleanup;
01372         return NULL;
01373     }
01374 
01375 
01376     /*
01377      * Fill position and flux vectors with appropriate values
01378      */
01379 
01380     cpl_bivector *positions = cpl_bivector_new(count);
01381     double *xpos = cpl_bivector_get_x_data(positions);
01382     double *ypos = cpl_bivector_get_y_data(positions);
01383     cpl_vector *fluxes = cpl_vector_new(count);
01384     double *flux = cpl_vector_get_data(fluxes);
01385 
01386     count = 0;
01387     for (y = 0; y < ny; y += step) {
01388         pos = y*nx;
01389         for (x = 0; x < nx; x += step, pos += step) {
01390             if (input_data[pos] > level) {
01391                 xpos[count] = x;
01392                 ypos[count] = y;
01393                 flux[count] = input_data[pos];
01394                 count++;
01395             }
01396         }
01397     }
01398 
01399     cpl_image_delete(input); input = NULL;
01400 
01401     /*
01402      * Do the fit, and fill the output image with the model
01403      * values in all pixels.
01404      */
01405 
01406     cpl_polynomial *model = cpl_polynomial_fit_2d_create(positions,
01407                                                          fluxes,
01408                                                          degree,
01409                                                          NULL);
01410 
01411     cpl_bivector_delete(positions);
01412     cpl_vector_delete(fluxes);
01413 
01414     smooth = cpl_image_new(nx, ny, FORS_IMAGE_TYPE);
01415     float *smooth_data = cpl_image_get_data_float(smooth);
01416 
01417     cpl_vector *point = cpl_vector_new(2);
01418     double *dpoint = cpl_vector_get_data(point);
01419 
01420     for (y = 0; y < ny; y++) {
01421         pos = y*nx;
01422         dpoint[1] = y;
01423         for (x = 0; x < nx; x++, pos++) {
01424             dpoint[0] = x;
01425             smooth_data[pos] = cpl_polynomial_eval(model, point);
01426         }
01427     }
01428 
01429     cpl_polynomial_delete(model);
01430     cpl_vector_delete(point);
01431 
01432     cleanup;
01433     return smooth;
01434 
01435 }
01436 
01437 #undef cleanup
01438 #define cleanup
01439 
01455 cpl_image *
01456 fors_image_filter_max_create(const fors_image *image, 
01457                              int xradius,
01458                              int yradius,
01459                              bool use_data)
01460 {
01461     const cpl_image *input = NULL;
01462     cpl_image *hmaxima = NULL;
01463     cpl_image *maxima = NULL;
01464     int nx, ny;
01465     
01466     assure( image != NULL, return maxima, NULL );
01467     passure( image->data != NULL, return maxima );
01468     passure( image->variance != NULL, return maxima );
01469     
01470     input = (use_data) ? image->data : image->variance;
01471 
01472     nx = cpl_image_get_size_x(input);
01473     ny = cpl_image_get_size_y(input);
01474 
01475     /*
01476      * Allocate space for horizontal max filter result.
01477      */
01478 
01479     hmaxima = cpl_image_duplicate(input);
01480 
01481     /* For efficiency reasons, assume that the image type is float */
01482     assure( FORS_IMAGE_TYPE == CPL_TYPE_FLOAT, return maxima, NULL );
01483 
01484     float *input_data  = (float *)cpl_image_get_data_float_const(input);
01485     float *maxima_data = cpl_image_get_data_float(hmaxima);
01486 
01487     int y;
01488     for (y = 0; y < ny; y++) {
01489         const float *irow = input_data + y * nx;
01490         float       *orow = maxima_data + y * nx;
01491         max_filter(irow, orow, nx, 2*xradius+1);
01492     }
01493 
01494     cpl_image_turn(hmaxima, 1);
01495 
01496     /*
01497      * Allocate space for vertical max filter result.
01498      */
01499 
01500     maxima = cpl_image_duplicate(hmaxima);
01501     input_data  = cpl_image_get_data_float(hmaxima);
01502     maxima_data = cpl_image_get_data_float(maxima);
01503 
01504     /*
01505      * Now nx is the y size of the rotated image...
01506      */
01507 
01508     int x;
01509     for (x = 0; x < nx; x++) {
01510         const float *irow = input_data + x * ny;
01511         float       *orow = maxima_data + x * ny;
01512         max_filter(irow, orow, ny, 2*yradius+1);
01513     }
01514 
01515     cpl_image_delete(hmaxima);
01516 
01517     cpl_image_turn(maxima, -1);
01518     
01519     return maxima;
01520 }
01521 
01522 #undef cleanup
01523 #define cleanup
01524 
01530 double
01531 fors_image_get_stdev(const fors_image *image, double *dstdev)
01532 {
01533     assure( image != NULL, return 0, NULL );
01534     assure( dstdev == NULL, return 0, "Unsupported");
01535 
01536     return cpl_image_get_stdev(image->data);
01537 }
01538 #undef cleanup
01539 #define cleanup \
01540 do { \
01541     cpl_mask_delete(rejected); \
01542     cpl_image_delete(im); \
01543 } while (0)
01544 
01552 double fors_image_get_stdev_robust(const fors_image *image, 
01553                    double cut,
01554                    double *dstdev)
01555 {
01556     cpl_mask *rejected = NULL;
01557     cpl_image *im = NULL;
01558 
01559     assure( image != NULL, return 0, NULL );
01560     assure( cut > 0, return 0, "Illegal cut: %f", cut );
01561     assure( dstdev == NULL, return 0, "Unsupported");
01562 
01563     double median = fors_image_get_median(image, NULL);
01564 
01565     im = cpl_image_duplicate(image->data);
01566     cpl_image_subtract_scalar(im, median); 
01567     cpl_image_power(im, 2);
01568     /* Now squared residuals wrt median */
01569     
01570     rejected = cpl_mask_threshold_image_create(image->data,
01571                                                median - cut,
01572                                                median + cut);
01573     cpl_mask_not(rejected);
01574     cpl_image_reject_from_mask(im, rejected);
01575 
01576     double robust_stdev = sqrt(cpl_image_get_mean(im));
01577 
01578     cleanup;
01579     return robust_stdev;
01580 }
01581 
01582 #undef cleanup
01583 #define cleanup
01584 
01594 double
01595 fors_image_get_error_mean(const fors_image *image, double *dmean)
01596 {
01597     double avg;
01598 
01599     assure( image != NULL, return 0, NULL );
01600     assure( dmean == NULL, return 0, "Unsupported");
01601 
01602     avg = cpl_image_get_mean(image->variance);
01603 
01604     /* This should never happen, but avoid sqrt of negative value in any case */
01605     assure( avg >= 0, return -1, "Average variance is %f", avg);
01606     
01607     return sqrt(avg);
01608 }
01609 
01610 
01611 #undef cleanup
01612 #define cleanup \
01613 do { \
01614     cpl_imagelist_delete(datlist); \
01615     cpl_imagelist_delete(varlist); \
01616 } while (0)
01617 
01626 fors_image *
01627 fors_image_collapse_create(const fors_image_list *images)
01628 {
01629     cpl_imagelist *datlist = NULL;
01630     cpl_imagelist *varlist = NULL;
01631     cpl_image *data = NULL;
01632     cpl_image *variance = NULL;
01633     const fors_image *i;
01634     int N = 0;
01635     
01636     assure( images != NULL, return NULL, NULL );
01637     assure( fors_image_list_size(images) > 0, return NULL, 
01638             "Cannot stack zero images");
01639 
01640     i = fors_image_list_first_const(images);
01641 
01642     datlist = cpl_imagelist_new();
01643     varlist = cpl_imagelist_new();
01644 
01645     while(i != NULL) {
01646 
01647         /* Append current image to image lists */
01648         cpl_imagelist_set(datlist, 
01649                           cpl_image_duplicate(i->data), 
01650                           cpl_imagelist_get_size(datlist));
01651         cpl_imagelist_set(varlist,
01652                           cpl_image_duplicate(i->variance),
01653                           cpl_imagelist_get_size(varlist));
01654         i = fors_image_list_next_const(images);
01655         N++;
01656     }
01657 
01658 #ifdef CPL_IS_NOT_CRAP
01659     data    = cpl_imagelist_collapse_create(datlist);
01660 
01661     variance = cpl_imagelist_collapse_create(varlist);
01662 #else
01663     data    = fors_imagelist_collapse_create(datlist);
01664 
01665     variance = fors_imagelist_collapse_create(varlist);
01666 #endif
01667 
01668     cpl_image_divide_scalar(variance, N);
01669 
01670     cleanup;
01671     return fors_image_new(data, variance);
01672 }
01673 
01674 
01675 #undef cleanup
01676 #define cleanup \
01677 do { \
01678     cpl_imagelist_delete(datlist); \
01679     cpl_imagelist_delete(varlist); \
01680 } while (0)
01681 
01692 fors_image *
01693 fors_image_collapse_minmax_create(const fors_image_list *images, 
01694                                   int low, int high)
01695 {
01696     cpl_imagelist *datlist = NULL;
01697     cpl_imagelist *varlist = NULL;
01698     cpl_image *data = NULL;
01699     cpl_image *variance = NULL;
01700     const fors_image *i;
01701     int N = 0;
01702     
01703     assure( images != NULL, return NULL, NULL );
01704     assure( fors_image_list_size(images) >  low + high, return NULL, 
01705             "Cannot reject more images than there are");
01706     assure( low*high >= 0 && low+high > 0, return NULL, 
01707             "Invalid minmax rejection criteria");
01708 
01709     i = fors_image_list_first_const(images);
01710 
01711     datlist = cpl_imagelist_new();
01712     varlist = cpl_imagelist_new();
01713 
01714     while(i != NULL) {
01715 
01716         /* Append current image to image lists */
01717         cpl_imagelist_set(datlist, 
01718                           cpl_image_duplicate(i->data), 
01719                           cpl_imagelist_get_size(datlist));
01720         cpl_imagelist_set(varlist,
01721                           cpl_image_duplicate(i->variance),
01722                           cpl_imagelist_get_size(varlist));
01723         i = fors_image_list_next_const(images);
01724         N++;
01725     }
01726 
01727     data     = cpl_imagelist_collapse_minmax_create(datlist, low, high);
01728     variance = cpl_imagelist_collapse_minmax_create(varlist, low, high);
01729 
01730     cpl_image_divide_scalar(variance, N);
01731 
01732     cleanup;
01733     return fors_image_new(data, variance);
01734 }
01735 
01747 fors_image *
01748 fors_image_collapse_median_create(const fors_image_list *images)
01749 {
01750     cpl_imagelist *datlist = NULL;
01751     cpl_imagelist *varlist = NULL;
01752     cpl_image *data = NULL;
01753     cpl_image *variance = NULL;
01754     const fors_image *i;
01755     int N = 0;
01756 
01757     assure( images != NULL, return NULL, NULL );
01758     assure( fors_image_list_size(images) > 0, return NULL, 
01759             "Cannot stack zero images");
01760 
01761     i = fors_image_list_first_const(images);
01762     
01763     datlist = cpl_imagelist_new();
01764     varlist = cpl_imagelist_new();
01765     while(i != NULL) {
01766         /* Append to image lists */
01767         cpl_imagelist_set(datlist, 
01768                           cpl_image_duplicate(i->data), 
01769                           cpl_imagelist_get_size(datlist));
01770         cpl_imagelist_set(varlist,
01771                           cpl_image_duplicate(i->variance),
01772                           cpl_imagelist_get_size(varlist));
01773 
01774         i = fors_image_list_next_const(images);
01775         N++;
01776     }
01777     
01778 #ifdef CPL_IS_NOT_CRAP
01779     data    = cpl_imagelist_collapse_median_create(datlist);
01780 
01781     variance = cpl_imagelist_collapse_create(varlist);
01782 #else
01783     data    = fors_imagelist_collapse_median_create(datlist);
01784 
01785     variance = fors_imagelist_collapse_create(varlist);
01786 #endif
01787 
01788     cpl_image_divide_scalar(variance, N);
01789 
01790     cpl_image_multiply_scalar(variance, 
01791                   fors_utils_median_corr(N) * 
01792                   fors_utils_median_corr(N));
01793     
01794     cleanup;
01795     return fors_image_new(data, variance);
01796 }
01797 
01798 #undef cleanup
01799 #define cleanup
01800 
01817 void fors_image_draw(fors_image *image, int type,
01818              double x, double y,
01819              int radius, double color)
01820 {
01821     assure( image != NULL, return, NULL );
01822 
01823     assure( type == 0 || type == 1 || type == 2,
01824             return , "Unsupported type %d", type);
01825 
01826     assure( radius > 0, return, NULL );
01827 
01828     if (type == 2) {
01829         int i;
01830         for (i = 0; i < 360; i++) {
01831             /* Step size of 1 degree is arbitrary */
01832 
01833             int px = x + radius*cos(i/(2*M_PI));
01834             int py = y + radius*sin(i/(2*M_PI));
01835             
01836             if (1 <= px && px <= cpl_image_get_size_x(image->data) &&
01837                 1 <= py && py <= cpl_image_get_size_y(image->data)) {
01838                 cpl_image_set(image->data, px, py, color);
01839                 cpl_image_set(image->variance, px, py, color > 0 ? color : 0);
01840             }
01841         }
01842     }
01843     else {
01844         int i;
01845 
01846         for (i = -radius; i <= radius; i++) {
01847 
01848             int px, py;
01849             
01850             if (type == 0) {
01851                 px = x + i;
01852                 py = y;
01853             }
01854             else {
01855                 px = x;
01856                 py = y + i;
01857             }
01858             
01859             if (1 <= px && px <= cpl_image_get_size_x(image->data) &&
01860                 1 <= py && py <= cpl_image_get_size_y(image->data)) {
01861                 cpl_image_set(image->data    , px, py, color);
01862                 cpl_image_set(image->variance, px, py, color > 0 ? color : 0);
01863             }
01864         }
01865     }
01866 
01867     return;
01868 }
01869 
01870 #define LIST_DEFINE
01871 #undef LIST_ELEM
01872 #define LIST_ELEM fors_image
01873 #include <list.h>
01874 

Generated on Wed Sep 10 07:31:51 2008 for FORS Pipeline Reference Manual by  doxygen 1.4.6