fors_image.c

00001 /* $Id: fors_image.c,v 1.57 2009/06/22 08:18:05 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: 2009/06/22 08:18:05 $
00024  * $Revision: 1.57 $
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     double ocorr = 0.0;
00365     int has_overscans = 0;
00366     int i;
00367 
00368     assure( frame != NULL, return image, NULL );
00369     /* bias may be NULL */
00370     assure( setting != NULL, return image, NULL );
00371     filename = cpl_frame_get_filename(frame);
00372     assure( filename != NULL, return image, 
00373             "NULL filename received");
00374     
00375     cpl_msg_info(cpl_func, "Loading %s: %s",
00376                  /* fors_frame_get_group_string(frame), */
00377          (cpl_frame_get_tag(frame) != NULL) ? 
00378          cpl_frame_get_tag(frame) : "NULL",
00379                  filename);
00380 
00381     /* Verify instrument setting */
00382     fors_setting_verify(setting, frame, &frame_setting);
00383     assure( !cpl_error_get_code(), return image, 
00384             "Could not verify %s setting", 
00385             filename);
00386     
00387 
00388     /* Get header */
00389     header = cpl_propertylist_load(filename, extension);
00390     assure( !cpl_error_get_code(), return image, 
00391             "Could not load %s extension %d header", 
00392             filename, extension);
00393 
00394     /* Get data */
00395     data = cpl_image_load(filename, 
00396                           FORS_IMAGE_TYPE, plane, extension);
00397     
00398     assure( !cpl_error_get_code(), return image, 
00399             "Could not load image from %s extension %d", 
00400             filename, extension);
00401 
00402     /* Read, remove overscan areas */
00403     {
00404         int xlow, ylow, xhig, yhig;
00405     /* bool check_consistency = false; */
00406 
00407         overscans = mos_load_overscans_fors(header);
00408     /* Uses always ESO.DET.OUT1.*  which is fine even
00409        for old FORS data */
00410 
00411         assure( !cpl_error_get_code(), return image, 
00412                 "Could not read overscan information from %s extension %d", 
00413                 filename, extension);
00414 
00415         for (i = 1; i < cpl_table_get_nrow(overscans); i++) {
00416 
00417             /*
00418              * Overscan correction
00419              */
00420         
00421             xlow = cpl_table_get_int(overscans, "xlow", i, NULL);
00422             ylow = cpl_table_get_int(overscans, "ylow", i, NULL);
00423             xhig = cpl_table_get_int(overscans, "xhig", i, NULL);
00424             yhig = cpl_table_get_int(overscans, "yhig", i, NULL);
00425 
00426             temp = cpl_image_extract(data, xlow+1, ylow+1, xhig, yhig);
00427 
00428             ocorr += cpl_image_get_median(temp);
00429 
00430             cpl_image_delete(temp); temp = NULL;
00431         }
00432 
00433         if (cpl_table_get_nrow(overscans) > 1) {
00434             has_overscans = 1;
00435             ocorr /= cpl_table_get_nrow(overscans) - 1;
00436         }
00437 
00438         xlow = cpl_table_get_int(overscans, "xlow", 0, NULL);
00439         ylow = cpl_table_get_int(overscans, "ylow", 0, NULL);
00440         xhig = cpl_table_get_int(overscans, "xhig", 0, NULL);
00441         yhig = cpl_table_get_int(overscans, "yhig", 0, NULL);
00442         
00443         cpl_table_delete(overscans); overscans = NULL;
00444         
00445         temp = cpl_image_duplicate(data);
00446         cpl_image_delete(data);
00447         data = cpl_image_extract(temp, xlow+1, ylow+1, xhig, yhig);
00448         
00449         cpl_image_delete(temp); temp = NULL;
00450     }
00451 
00452     cpl_propertylist_delete(header); header = NULL;       
00453 
00454     /* Define variance */
00455     if (cpl_frame_get_nextensions(frame) == 0 || bias != NULL) {
00456 
00457         /* No error bars provided, assume RON only.
00458            If frame is not a bias, add photon noise later.
00459         */
00460 
00461         /* Define read-out-noise */
00462         variance = cpl_image_new(
00463             cpl_image_get_size_x(data),
00464             cpl_image_get_size_y(data),
00465             FORS_IMAGE_TYPE);        
00466 
00467         if (cpl_frame_get_group(frame) == CPL_FRAME_GROUP_CALIB) {
00468             cpl_msg_warning(cpl_func, 
00469                             "No error bars provided for calibration frame %s, "
00470                             "assuming no errors. For complete error propagation, "
00471                             "you may recreate this frame with this pipeline", 
00472                             filename);
00473         } else {
00474             cpl_image_add_scalar(variance, frame_setting->ron*frame_setting->ron);
00475         }
00476     }
00477     else {
00478 
00479         extension = 1;
00480         
00481         /* Get error bars */
00482         variance = cpl_image_load(filename, 
00483                                   FORS_IMAGE_TYPE, plane, extension);
00484         
00485         assure( !cpl_error_get_code(), return image, 
00486                 "Could not load image from %s extension %d", 
00487                 filename, extension);
00488 
00489         cpl_image_power(variance, 2);
00490 
00491         assure( cpl_image_get_min(variance) >= 0, return image,
00492                 "Illegal minimum variance: %g",
00493                 cpl_image_get_min(variance));
00494 
00495         cpl_image_delete(temp); temp = NULL;
00496     }
00497     
00498     image = fors_image_new(data, variance);
00499 
00500     /* Count saturated pixels */
00501     if (saturated != NULL) {
00502         double lo_cut = 0.5;
00503         double hi_cut = 65534.5;
00504         non_saturated = cpl_mask_threshold_image_create(
00505             data, lo_cut, hi_cut);
00506         
00507         *saturated = (cpl_image_get_size_x(data) * 
00508                       cpl_image_get_size_y(data) - 
00509                       cpl_mask_count(non_saturated)) * 100.0 / 
00510             (cpl_image_get_size_x(data) * 
00511              cpl_image_get_size_y(data));
00512 
00513         cpl_mask_delete(non_saturated); non_saturated = NULL;
00514 
00515         cpl_msg_debug(cpl_func,
00516                      "%f %% saturated pixels", *saturated);
00517     }        
00518     
00519     if (bias == NULL) {
00520         /* Done. No bias to subtract */
00521     }
00522     else {
00523         assure( cpl_frame_get_group(frame) == CPL_FRAME_GROUP_RAW,
00524                 return image,
00525                 "Refusing to subtract bias from non-raw (%s) input frame: %s",
00526                 fors_frame_get_group_string(frame),
00527                 filename);
00528         
00529         cpl_msg_debug(cpl_func, "Subtracting bias from %s", filename);
00530         
00531         fors_image_subtract(image, bias);
00532         
00533         assure( !cpl_error_get_code(), return image, 
00534                 "Bias subtraction failed" );
00535 
00536         if (has_overscans) {
00537 
00538             /*
00539              * Overscan correction
00540              */
00541 
00542             ocorr -= cpl_image_get_median(bias->data);
00543 
00544             cpl_msg_info(cpl_func, 
00545                          "Overscan correction applied: %.2f ADUs", ocorr);
00546 
00547             fors_image_subtract_scalar(image, ocorr, -1);
00548         }
00549         
00550         /* Variance is now (ron**2  +  bias_noise**2).
00551            Add photonic noise:
00552            
00553            variance := variance + |f|/conad
00554         */
00555 
00556         double conad = 1.0 / frame_setting->average_gain;
00557 
00558         temp = cpl_image_divide_scalar_create(image->data, conad);
00559         cpl_image_abs(temp);
00560 
00561         cpl_image_add(image->variance, temp);
00562     }
00563     
00564     cleanup;
00565     return image;
00566 }
00567 
00568 
00569 #undef cleanup
00570 #define cleanup \
00571 do { \
00572     cpl_image_delete(sigma); \
00573 } while(0)
00574 
00583 void
00584 fors_image_save(const fors_image *image, const cpl_propertylist *header,
00585                 const char *filename)
00586 {
00587     cpl_propertylist *extension_header = NULL;
00588     cpl_image *sigma = NULL;
00589 
00590     assure( image != NULL, return, NULL );
00591     /* header may be NULL */
00592     assure( filename != NULL, return, NULL );
00593     
00594     cpl_image_save(image->data, filename, CPL_BPP_IEEE_FLOAT, header,
00595                    CPL_IO_DEFAULT);
00596     assure( !cpl_error_get_code(), return, 
00597             "Cannot save product %s", filename);
00598     
00599     sigma = cpl_image_power_create(image->variance, 0.5);
00600     /* This would probably be faster if sqrt() is used rather than pow */
00601 
00602     cpl_image_save(sigma, filename, CPL_BPP_IEEE_FLOAT, extension_header,
00603                    CPL_IO_EXTEND);
00604     assure( !cpl_error_get_code(), return, 
00605             "Cannot save product %s", filename);
00606     
00607     cleanup;
00608     return;
00609 }
00610 
00611 
00612 #undef cleanup
00613 #define cleanup \
00614 do { \
00615     cpl_image_delete(var_bkg); \
00616     cpl_image_delete(sigma_bkg); \
00617 } while(0)
00618 
00629 void
00630 fors_image_save_sex(const fors_image *image, const cpl_propertylist *header,
00631                     const char *filename_dat,
00632                     const char *filename_var,
00633                     int radius)
00634 {
00635     cpl_propertylist *extension_header = NULL;
00636     cpl_image *sigma_bkg = NULL;
00637     cpl_image *var_bkg = NULL;
00638 
00639     assure( image != NULL, return, NULL );
00640     /* header may be NULL */
00641     assure( filename_dat != NULL, return, NULL );
00642     assure( filename_var != NULL, return, NULL );
00643 
00644     cpl_image_save(image->data, filename_dat, CPL_BPP_IEEE_FLOAT, header,
00645                    CPL_IO_DEFAULT);
00646     assure( !cpl_error_get_code(), return, 
00647             "Cannot save product %s", filename_dat);
00648     
00649     /* Sextractor wants as input the background error bars,
00650        i.e. excluding sources.
00651        Therefore filter away sources but keep the sharp edges
00652        between the illuminated / non-illuminated areas.
00653 
00654        I.e. use a median filter, average filter would not work.
00655     */
00656 
00657     cpl_msg_info(cpl_func, "Creating background error map");
00658 
00659     bool filter_data = false;  /* filter the variance image */
00660     int xstep = radius/2; /* 25 points sampling grid 
00661                              . . . . .
00662                              . . . . .
00663                              . . . . .
00664                              . . . . .
00665                              . . . . .
00666                            */
00667     int ystep = radius/2;
00668     int xstart = 1;
00669     int ystart = 1;
00670     int xend = fors_image_get_size_x(image);
00671     int yend = fors_image_get_size_y(image);
00672 
00673 
00674     var_bkg = fors_image_filter_median_create(image, 
00675                                               radius,
00676                                               radius,
00677                                               xstart, ystart,
00678                                               xend, yend,
00679                                               xstep, ystep,
00680                                               filter_data);
00681     assure( !cpl_error_get_code(), return, 
00682             "Median filtering failed");
00683 
00684     sigma_bkg = cpl_image_power_create(var_bkg, 0.5);
00685 
00686     cpl_image_save(sigma_bkg, filename_var,
00687                    CPL_BPP_IEEE_FLOAT, extension_header,
00688                    CPL_IO_DEFAULT);
00689     assure( !cpl_error_get_code(), return, 
00690             "Cannot save product %s", filename_var);
00691 
00692     cleanup;
00693     return;
00694 }
00695 
00696 #undef cleanup
00697 #define cleanup
00698 
00703 int fors_image_get_size_x(const fors_image *image)
00704 {
00705     assure( image != NULL, return -1, NULL );
00706     return cpl_image_get_size_x(image->data);
00707 }
00708 
00709 #undef cleanup
00710 #define cleanup
00711 
00716 int fors_image_get_size_y(const fors_image *image)
00717 {
00718     assure( image != NULL, return -1, NULL );
00719     return cpl_image_get_size_y(image->data);
00720 }
00721 
00722 #undef cleanup
00723 #define cleanup
00724 
00728 const float *fors_image_get_data_const(const fors_image *image)
00729 {
00730     assure( image != NULL, return NULL, NULL );
00731 
00732     assure( FORS_IMAGE_TYPE == CPL_TYPE_FLOAT, return NULL, NULL );
00733     /* This function (including API) would need to change
00734        if the pixel type changes */
00735 
00736     return cpl_image_get_data_float(image->data);
00737 }
00738 
00739 #undef cleanup
00740 #define cleanup
00741 
00748 void
00749 fors_image_abs(fors_image *image)
00750 {
00751     assure( image != NULL, return, NULL );
00752 
00753     cpl_image_abs(image->data);
00754 
00755     return;
00756 }
00757 
00758 #undef cleanup
00759 #define cleanup
00760 
00767 void
00768 fors_image_square(fors_image *image)
00769 {
00770     assure( image != NULL, return, NULL );
00771 
00772     cpl_image_multiply(image->data, image->data);
00773     /* It is an undocumented feature of CPL that you
00774        can pass the same image to cpl_image_multiply and get
00775        the right answer. Let us hope it does not change...
00776     */
00777     cpl_image_multiply_scalar(image->variance, 2);
00778 
00779     return;
00780 }
00781 
00782 
00783 #undef cleanup
00784 #define cleanup \
00785 do { \
00786     cpl_image_delete(temp); \
00787 } while(0)
00788 
00796 void
00797 fors_image_subtract(fors_image *left, const fors_image *right)
00798 {
00799     cpl_image *temp = NULL;
00800     assure( left != NULL, return, NULL );
00801     assure( right != NULL, return, NULL );
00802 
00803     cpl_image_subtract(left->data, right->data);
00804 
00805     /*  variance_left := variance_left + variance_right */
00806     cpl_image_add(left->variance, right->variance);
00807 
00808     cleanup;
00809     return;
00810 }
00811 
00812 #undef cleanup
00813 #define cleanup
00814 
00824 void
00825 fors_image_multiply_noerr(fors_image *left, const cpl_image *right)
00826 {
00827     assure( left != NULL, return, NULL );
00828     assure( right != NULL, return, NULL );
00829     assure( cpl_image_get_size_x(left->data) == cpl_image_get_size_x(right) &&
00830             cpl_image_get_size_y(left->data) == cpl_image_get_size_y(right),
00831             return, "Incompatible image sizes: %dx%d and %dx%d",
00832             cpl_image_get_size_x(left->data),
00833             cpl_image_get_size_y(left->data),
00834             cpl_image_get_size_x(right),
00835             cpl_image_get_size_y(right));
00836 
00837     cpl_image_multiply(left->data, right);
00838     cpl_image_multiply(left->variance, right);
00839     cpl_image_multiply(left->variance, right);
00840 
00841     return;
00842 }
00843 
00844 #undef cleanup
00845 #define cleanup
00846 
00862 void
00863 fors_image_divide_noerr(fors_image *left, cpl_image *right)
00864 {
00865     assure( left != NULL, return, NULL );
00866     assure( right != NULL, return, NULL );
00867     assure( cpl_image_get_size_x(left->data) == cpl_image_get_size_x(right) &&
00868             cpl_image_get_size_y(left->data) == cpl_image_get_size_y(right),
00869             return, "Incompatible image sizes: %dx%d and %dx%d",
00870             cpl_image_get_size_x(left->data),
00871             cpl_image_get_size_y(left->data),
00872             cpl_image_get_size_x(right),
00873             cpl_image_get_size_y(right));
00874 
00875     int x, y;
00876     int nx = cpl_image_get_size_x(right);
00877     int ny = cpl_image_get_size_y(right);
00878     float *datal = cpl_image_get_data_float(left->data);
00879     float *datav = cpl_image_get_data_float(left->variance);
00880     float *datar = cpl_image_get_data_float(right);
00881     for (y = 0; y < ny; y++) {
00882         for (x = 0; x < nx; x++) {
00883             if (datar[x + nx*y] == 0) {
00884                 datar[x + nx*y] = 1;
00885                 datal[x + nx*y] = 1;
00886 
00887                 datav[x + nx*y] = FORS_IMAGE_TYPE_MAX;
00888             }
00889         }
00890     }
00891 
00892     cpl_image_divide(left->data, right);
00893     cpl_image_divide(left->variance, right);
00894     cpl_image_divide(left->variance, right);
00895 
00896     return;
00897 }
00898 
00899 #undef cleanup
00900 #define cleanup \
00901 do { \
00902     fors_image_delete(&dupl); \
00903 } while(0)
00904 
00924 void
00925 fors_image_divide(fors_image *left, const fors_image *right)
00926 {
00927     fors_image *dupl = NULL;
00928 
00929     assure( left  != NULL, return, NULL );
00930     assure( right != NULL, return, NULL );
00931 
00932     dupl = fors_image_duplicate(right);
00933 
00934     cpl_image_divide(left->data, dupl->data); 
00935     /* This CPL function divides by zero by setting  x/0 = 1 for all x */
00936 
00937     cpl_image_multiply(dupl->variance, left->data);
00938     cpl_image_multiply(dupl->variance, left->data);
00939 
00940     /* Now  dupl->variance = sigma2^2 * data1^2 / data2^2 */
00941 
00942     cpl_image_add(left->variance, dupl->variance);
00943 
00944     /* Now  left->variance = sigma1^2 + sigma2^2 * data1^2 / data2^2 */
00945 
00946     cpl_image_divide(left->variance, dupl->data);
00947     cpl_image_divide(left->variance, dupl->data);
00948     /* QED */
00949 
00950     /* Handle division by zero */
00951     int x, y;
00952     int nx = cpl_image_get_size_x(left->data);
00953     int ny = cpl_image_get_size_y(left->data);
00954     float *datal = cpl_image_get_data_float(left->data);
00955     float *datav = cpl_image_get_data_float(left->variance);
00956     float *datar = cpl_image_get_data_float(right->data);
00957     for (y = 0; y < ny; y++) {
00958         for (x = 0; x < nx; x++) {
00959             if (datar[x + nx*y] == 0) {
00960                 datal[x + nx*y] = 1;
00961                 datav[x + nx*y] = FORS_IMAGE_TYPE_MAX;
00962             }
00963         }
00964     }
00965 
00966     cleanup;
00967     return;
00968 }
00969 
00970 #undef cleanup
00971 #define cleanup \
00972 do { \
00973     cpl_image_delete(s22d12); \
00974 } while(0)
00975 
00986 void
00987 fors_image_multiply(fors_image *left, const fors_image *right)
00988 {
00989     cpl_image *s22d12 = NULL;
00990 
00991     assure( left  != NULL, return, NULL );
00992     assure( right != NULL, return, NULL );
00993 
00994     s22d12 = cpl_image_duplicate(right->variance);
00995     cpl_image_multiply(s22d12, left->data);
00996     cpl_image_multiply(s22d12, left->data);
00997     
00998     cpl_image_multiply(left->variance, right->data);
00999     cpl_image_multiply(left->variance, right->data);
01000     cpl_image_add(left->variance, s22d12);
01001 
01002     cpl_image_multiply(left->data, right->data);
01003 
01004     cleanup;
01005     return;
01006 }
01007 
01008 
01009 #undef cleanup
01010 #define cleanup
01011 
01023 void fors_image_subtract_scalar(fors_image *image, double s, double ds)
01024 {
01025     assure( image != NULL, return, NULL );
01026     assure( ds <= 0, return, "Unsupported");
01027 
01028     cpl_image_subtract_scalar(image->data, s);
01029 
01030     return;
01031 }
01032 
01033 
01034 #undef cleanup
01035 #define cleanup
01036 
01048 void fors_image_divide_scalar(fors_image *image, double s, double ds)
01049 {
01050     assure( image != NULL, return, NULL );
01051     assure( s != 0, return, "Division by zero");
01052     assure( ds <= 0, return, "Unsupported");
01053 
01054     cpl_image_divide_scalar(image->data, s);
01055     cpl_image_divide_scalar(image->variance, s*s);
01056     
01057     return;
01058 }
01059 
01060 #undef cleanup
01061 #define cleanup
01062 
01074 void fors_image_multiply_scalar(fors_image *image, double s, double ds)
01075 {
01076     assure( image != NULL, return, NULL );
01077     assure( ds <= 0, return, "Unsupported");
01078 
01079     cpl_image_multiply_scalar(image->data, s);
01080     cpl_image_multiply_scalar(image->variance, s*s);
01081 
01082     return;
01083 }
01084 
01085 #undef cleanup
01086 #define cleanup \
01087 do { \
01088     cpl_image_delete(temp); \
01089 } while(0)
01090 
01103 void fors_image_exponential(fors_image *image, double b, double db)
01104 {
01105     cpl_image *temp = NULL;
01106 
01107     assure( image != NULL, return, NULL );
01108     assure( b >= 0, return, "Negative base: %f", b);
01109     assure( db <= 0, return, "Unsupported");
01110 
01111     cpl_image_exponential(image->data, b);
01112 
01113     double lnb = log(b);
01114     
01115     cpl_image_multiply_scalar(image->variance, lnb*lnb);
01116     cpl_image_multiply(image->variance, image->data);
01117     cpl_image_multiply(image->variance, image->data);
01118 
01119     return;
01120 }
01121 
01122 
01123 #undef cleanup
01124 #define cleanup
01125 
01130 double
01131 fors_image_get_min(const fors_image *image)
01132 {
01133     assure( image != NULL, return 0, NULL );
01134 
01135     return cpl_image_get_min(image->data);
01136 }
01137 
01138 #undef cleanup
01139 #define cleanup
01140 
01145 double
01146 fors_image_get_max(const fors_image *image)
01147 {
01148     assure( image != NULL, return 0, NULL );
01149 
01150     return cpl_image_get_max(image->data);
01151 }
01152 
01153 #undef cleanup
01154 #define cleanup
01155 
01161 double
01162 fors_image_get_mean(const fors_image *image, double *dmean)
01163 {
01164     assure( image != NULL, return 0, NULL );
01165     assure( dmean == NULL, return 0, "Unsupported");
01166 
01167     return cpl_image_get_mean(image->data);
01168 }
01169 
01170 #undef cleanup
01171 #define cleanup
01172 
01178 double
01179 fors_image_get_median(const fors_image *image, double *dmedian)
01180 {
01181     assure( image != NULL, return 0, NULL );
01182     assure( dmedian == NULL, return 0, "Unsupported");
01183 
01184     return cpl_image_get_median(image->data);
01185 }
01186 
01187 
01188 #undef cleanup
01189 #define cleanup
01190 
01204 void fors_image_crop(fors_image *image,
01205              int xlo, int ylo,
01206              int xhi, int yhi)
01207 {
01208     /* CPL is missing the function to locally extract an image,
01209        so this this inefficient CPL function */
01210     assure( image != NULL, return, NULL );
01211     assure( 1 <= xlo && xlo <= xhi && xhi <= fors_image_get_size_x(image) &&
01212             1 <= ylo && ylo <= yhi && yhi <= fors_image_get_size_y(image),
01213             return, "Cannot extraction region (%d, %d) - (%d, %d) of "
01214             "%dx%d image",
01215             xlo, ylo, xhi, yhi,
01216             fors_image_get_size_x(image),
01217             fors_image_get_size_y(image));
01218     
01219     cpl_image *new_data = cpl_image_extract(image->data,
01220                                             xlo, ylo,
01221                                             xhi, yhi);
01222     cpl_image_delete(image->data);
01223 
01224     cpl_image* new_variance = cpl_image_extract(image->variance,
01225                                                 xlo, ylo,
01226                                                 xhi, yhi);
01227     cpl_image_delete(image->variance);
01228 
01229     image->data = new_data;
01230     image->variance = new_variance;
01231 
01232     return;
01233 }
01234 
01260 cpl_image *
01261 fors_image_filter_median_create(const fors_image *image, 
01262                                 int xradius,
01263                                 int yradius,
01264                                 int xstart, 
01265                                 int ystart,
01266                                 int xend,
01267                                 int yend,
01268                                 int xstep,
01269                                 int ystep,
01270                                 bool use_data)
01271 {
01272     const cpl_image *input = NULL;
01273     cpl_image *smooth = NULL;
01274     int nx, ny;
01275     
01276     assure( image != NULL, return smooth, NULL );
01277     passure( image->data != NULL, return smooth );
01278     passure( image->variance != NULL, return smooth );
01279     
01280     input = (use_data) ? image->data : image->variance;
01281 
01282     nx = cpl_image_get_size_x(input);
01283     ny = cpl_image_get_size_y(input);
01284 
01285     if (xstep < 1) xstep = 1;
01286     if (ystep < 1) ystep = 1;
01287 
01288     assure( 1 <= xstart && xstart <= xend && xend <= nx &&
01289             1 <= ystart && ystart <= yend && yend <= ny, return smooth,
01290             "Illegal region (%d, %d) - (%d, %d) of %dx%d image",
01291             xstart, ystart,
01292             xend, yend,
01293             nx, ny);
01294     
01295     smooth = cpl_image_duplicate(input);
01296 
01297     /* For efficiency reasons, assume that the image type is float */
01298     assure( FORS_IMAGE_TYPE == CPL_TYPE_FLOAT, return smooth, NULL );
01299 
01300     const float *input_data  = cpl_image_get_data_float_const(input);
01301     float *smooth_data = cpl_image_get_data_float(smooth);
01302     float *data = cpl_malloc((2*yradius + 1)*(2*xradius + 1)*sizeof(*data));
01303     
01304     int y;
01305     for (y = ystart; y < yend; y++) {
01306         /*
01307           Sample kernel on grid which always contains the central pixel
01308           
01309           Trim window (note: this will cause fewer values to
01310           be used for the median near the region borders 
01311         */
01312         int ylo = y - (yradius/ystep) * ystep;
01313         int yhi = y + (yradius/ystep) * ystep;
01314         
01315         while (ylo < ystart) ylo += ystep;
01316         while (yhi > yend  ) yhi -= ystep;
01317         
01318         int x;
01319         for (x = xstart; x < xend; x++) {
01320             int xlo = x - (xradius/xstep) * xstep;
01321             int xhi = x + (xradius/xstep) * xstep;
01322             
01323             while (xlo < xstart) xlo += xstep;
01324             while (xhi > xend  ) xhi -= xstep;
01325             
01326             /* Collect data */
01327             int k = 0;
01328             int j, i;
01329             for (j = ylo; j <= yhi; j += ystep) {
01330                 for (i = xlo; i <= xhi; i += xstep) {
01331                     data[k++] = input_data[ (i-1) + (j-1)*nx ];
01332                 }
01333             }
01334         
01335             /* Get median */
01336             smooth_data[ (x-1) + (y-1)*nx ] = 
01337                                fors_tools_get_median_float(data, k);
01338         }
01339     }
01340 
01341     cpl_free(data);
01342     return smooth;
01343 }
01344 
01345 #undef cleanup
01346 #define cleanup \
01347 do { \
01348     cpl_image_delete(input); \
01349 } while(0)
01350 cpl_image *
01351 fors_image_flat_fit_create(fors_image *image, 
01352                            int step, 
01353                            int degree, 
01354                            float level)
01355 {
01356     cpl_image *temp = NULL;
01357     cpl_image *input = NULL;
01358     cpl_image *smooth = NULL;
01359     int nx, ny;
01360 
01361     assure( image != NULL, return smooth, NULL );
01362     passure( image->data != NULL, return smooth );
01363     assure( step > 0, return smooth, NULL );
01364     assure( degree >= 0, return smooth, NULL );
01365 
01366 
01367     temp = image->data;
01368 
01369     nx = cpl_image_get_size_x(temp);
01370     ny = cpl_image_get_size_y(temp);
01371 
01372     /* 
01373      * For efficiency reasons, assume that the image type is float 
01374      */
01375 
01376     assure( FORS_IMAGE_TYPE == CPL_TYPE_FLOAT, return smooth, NULL );
01377 
01378     /*
01379      * Apply light median filter, to eliminate big outliers from fit
01380      */
01381 
01382     cpl_matrix *kernel = cpl_matrix_new(3, 3);
01383     cpl_matrix_fill(kernel, 1.0);
01384     input = cpl_image_filter_median(image->data, kernel);
01385     cpl_matrix_delete(kernel);
01386     const float *input_data = cpl_image_get_data_float_const(input);
01387 
01388     /*
01389      * First of all, count how many points will have to be fitted
01390      */
01391 
01392     int x, y, pos;
01393     int count = 0;
01394     for (y = 0; y < ny; y += step) {
01395         pos = y*nx;
01396         for (x = 0; x < nx; x += step, pos += step) {
01397             if (input_data[pos] > level) {
01398                 count++;
01399             }
01400         }
01401     }
01402 
01403     if (count < (degree+1)*(degree+2)) {
01404         step = sqrt((nx*nx)/((degree+1)*(degree+2))) / 2;
01405         if (step == 0)
01406             step = 1;
01407         cpl_msg_error(cpl_func, "Flat field image too small (%dx%d). "
01408                       "Please provide a smaller resampling step (a good "
01409                       "value would be %d)", nx, ny, step);
01410         cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_INPUT);
01411         cleanup;
01412         return NULL;
01413     }
01414 
01415 
01416     /*
01417      * Fill position and flux vectors with appropriate values
01418      */
01419 
01420     cpl_bivector *positions = cpl_bivector_new(count);
01421     double *xpos = cpl_bivector_get_x_data(positions);
01422     double *ypos = cpl_bivector_get_y_data(positions);
01423     cpl_vector *fluxes = cpl_vector_new(count);
01424     double *flux = cpl_vector_get_data(fluxes);
01425 
01426     count = 0;
01427     for (y = 0; y < ny; y += step) {
01428         pos = y*nx;
01429         for (x = 0; x < nx; x += step, pos += step) {
01430             if (input_data[pos] > level) {
01431                 xpos[count] = x;
01432                 ypos[count] = y;
01433                 flux[count] = input_data[pos];
01434                 count++;
01435             }
01436         }
01437     }
01438 
01439     cpl_image_delete(input); input = NULL;
01440 
01441     /*
01442      * Do the fit, and fill the output image with the model
01443      * values in all pixels.
01444      */
01445 
01446     cpl_polynomial *model = cpl_polynomial_fit_2d_create(positions,
01447                                                          fluxes,
01448                                                          degree,
01449                                                          NULL);
01450 
01451     cpl_bivector_delete(positions);
01452     cpl_vector_delete(fluxes);
01453 
01454     smooth = cpl_image_new(nx, ny, FORS_IMAGE_TYPE);
01455     float *smooth_data = cpl_image_get_data_float(smooth);
01456 
01457     cpl_vector *point = cpl_vector_new(2);
01458     double *dpoint = cpl_vector_get_data(point);
01459 
01460     for (y = 0; y < ny; y++) {
01461         pos = y*nx;
01462         dpoint[1] = y;
01463         for (x = 0; x < nx; x++, pos++) {
01464             dpoint[0] = x;
01465             smooth_data[pos] = cpl_polynomial_eval(model, point);
01466         }
01467     }
01468 
01469     cpl_polynomial_delete(model);
01470     cpl_vector_delete(point);
01471 
01472     cleanup;
01473     return smooth;
01474 
01475 }
01476 
01477 #undef cleanup
01478 #define cleanup
01479 
01495 cpl_image *
01496 fors_image_filter_max_create(const fors_image *image, 
01497                              int xradius,
01498                              int yradius,
01499                              bool use_data)
01500 {
01501     const cpl_image *input = NULL;
01502     cpl_image *hmaxima = NULL;
01503     cpl_image *maxima = NULL;
01504     int nx, ny;
01505     
01506     assure( image != NULL, return maxima, NULL );
01507     passure( image->data != NULL, return maxima );
01508     passure( image->variance != NULL, return maxima );
01509     
01510     input = (use_data) ? image->data : image->variance;
01511 
01512     nx = cpl_image_get_size_x(input);
01513     ny = cpl_image_get_size_y(input);
01514 
01515     /*
01516      * Allocate space for horizontal max filter result.
01517      */
01518 
01519     hmaxima = cpl_image_duplicate(input);
01520 
01521     /* For efficiency reasons, assume that the image type is float */
01522     assure( FORS_IMAGE_TYPE == CPL_TYPE_FLOAT, return maxima, NULL );
01523 
01524     float *input_data  = (float *)cpl_image_get_data_float_const(input);
01525     float *maxima_data = cpl_image_get_data_float(hmaxima);
01526 
01527     int y;
01528     for (y = 0; y < ny; y++) {
01529         const float *irow = input_data + y * nx;
01530         float       *orow = maxima_data + y * nx;
01531         max_filter(irow, orow, nx, 2*xradius+1);
01532     }
01533 
01534     cpl_image_turn(hmaxima, 1);
01535 
01536     /*
01537      * Allocate space for vertical max filter result.
01538      */
01539 
01540     maxima = cpl_image_duplicate(hmaxima);
01541     input_data  = cpl_image_get_data_float(hmaxima);
01542     maxima_data = cpl_image_get_data_float(maxima);
01543 
01544     /*
01545      * Now nx is the y size of the rotated image...
01546      */
01547 
01548     int x;
01549     for (x = 0; x < nx; x++) {
01550         const float *irow = input_data + x * ny;
01551         float       *orow = maxima_data + x * ny;
01552         max_filter(irow, orow, ny, 2*yradius+1);
01553     }
01554 
01555     cpl_image_delete(hmaxima);
01556 
01557     cpl_image_turn(maxima, -1);
01558     
01559     return maxima;
01560 }
01561 
01562 #undef cleanup
01563 #define cleanup
01564 
01570 double
01571 fors_image_get_stdev(const fors_image *image, double *dstdev)
01572 {
01573     assure( image != NULL, return 0, NULL );
01574     assure( dstdev == NULL, return 0, "Unsupported");
01575 
01576     return cpl_image_get_stdev(image->data);
01577 }
01578 #undef cleanup
01579 #define cleanup \
01580 do { \
01581     cpl_mask_delete(rejected); \
01582     cpl_image_delete(im); \
01583 } while (0)
01584 
01592 double fors_image_get_stdev_robust(const fors_image *image, 
01593                    double cut,
01594                    double *dstdev)
01595 {
01596     cpl_mask *rejected = NULL;
01597     cpl_image *im = NULL;
01598 
01599     assure( image != NULL, return 0, NULL );
01600     assure( cut > 0, return 0, "Illegal cut: %f", cut );
01601     assure( dstdev == NULL, return 0, "Unsupported");
01602 
01603     double median = fors_image_get_median(image, NULL);
01604 
01605     im = cpl_image_duplicate(image->data);
01606     cpl_image_subtract_scalar(im, median); 
01607     cpl_image_power(im, 2);
01608     /* Now squared residuals wrt median */
01609     
01610     rejected = cpl_mask_threshold_image_create(image->data,
01611                                                median - cut,
01612                                                median + cut);
01613     cpl_mask_not(rejected);
01614     cpl_image_reject_from_mask(im, rejected);
01615 
01616     double robust_stdev = sqrt(cpl_image_get_mean(im));
01617 
01618     cleanup;
01619     return robust_stdev;
01620 }
01621 
01622 #undef cleanup
01623 #define cleanup
01624 
01634 double
01635 fors_image_get_error_mean(const fors_image *image, double *dmean)
01636 {
01637     double avg;
01638 
01639     assure( image != NULL, return 0, NULL );
01640     assure( dmean == NULL, return 0, "Unsupported");
01641 
01642     avg = cpl_image_get_mean(image->variance);
01643 
01644     /* This should never happen, but avoid sqrt of negative value in any case */
01645     assure( avg >= 0, return -1, "Average variance is %f", avg);
01646     
01647     return sqrt(avg);
01648 }
01649 
01650 
01651 #undef cleanup
01652 #define cleanup \
01653 do { \
01654     cpl_imagelist_delete(datlist); \
01655     cpl_imagelist_delete(varlist); \
01656 } while (0)
01657 
01666 fors_image *
01667 fors_image_collapse_create(const fors_image_list *images)
01668 {
01669     cpl_imagelist *datlist = NULL;
01670     cpl_imagelist *varlist = NULL;
01671     cpl_image *data = NULL;
01672     cpl_image *variance = NULL;
01673     const fors_image *i;
01674     int N = 0;
01675     
01676     assure( images != NULL, return NULL, NULL );
01677     assure( fors_image_list_size(images) > 0, return NULL, 
01678             "Cannot stack zero images");
01679 
01680     i = fors_image_list_first_const(images);
01681 
01682     datlist = cpl_imagelist_new();
01683     varlist = cpl_imagelist_new();
01684 
01685     while(i != NULL) {
01686 
01687         /* Append current image to image lists */
01688         cpl_imagelist_set(datlist, 
01689                           cpl_image_duplicate(i->data), 
01690                           cpl_imagelist_get_size(datlist));
01691         cpl_imagelist_set(varlist,
01692                           cpl_image_duplicate(i->variance),
01693                           cpl_imagelist_get_size(varlist));
01694         i = fors_image_list_next_const(images);
01695         N++;
01696     }
01697 
01698 #ifdef CPL_IS_NOT_CRAP
01699     data    = cpl_imagelist_collapse_create(datlist);
01700 
01701     variance = cpl_imagelist_collapse_create(varlist);
01702 #else
01703     data    = fors_imagelist_collapse_create(datlist);
01704 
01705     variance = fors_imagelist_collapse_create(varlist);
01706 #endif
01707 
01708     cpl_image_divide_scalar(variance, N);
01709 
01710     cleanup;
01711     return fors_image_new(data, variance);
01712 }
01713 
01714 
01715 #undef cleanup
01716 #define cleanup \
01717 do { \
01718     cpl_imagelist_delete(datlist); \
01719     cpl_imagelist_delete(varlist); \
01720 } while (0)
01721 
01732 fors_image *
01733 fors_image_collapse_minmax_create(const fors_image_list *images, 
01734                                   int low, int high)
01735 {
01736     cpl_imagelist *datlist = NULL;
01737     cpl_imagelist *varlist = NULL;
01738     cpl_image *data = NULL;
01739     cpl_image *variance = NULL;
01740     const fors_image *i;
01741     int N = 0;
01742     
01743     assure( images != NULL, return NULL, NULL );
01744     assure( fors_image_list_size(images) >  low + high, return NULL, 
01745             "Cannot reject more images than there are");
01746     assure( low*high >= 0 && low+high > 0, return NULL, 
01747             "Invalid minmax rejection criteria");
01748 
01749     i = fors_image_list_first_const(images);
01750 
01751     datlist = cpl_imagelist_new();
01752     varlist = cpl_imagelist_new();
01753 
01754     while(i != NULL) {
01755 
01756         /* Append current image to image lists */
01757         cpl_imagelist_set(datlist, 
01758                           cpl_image_duplicate(i->data), 
01759                           cpl_imagelist_get_size(datlist));
01760         cpl_imagelist_set(varlist,
01761                           cpl_image_duplicate(i->variance),
01762                           cpl_imagelist_get_size(varlist));
01763         i = fors_image_list_next_const(images);
01764         N++;
01765     }
01766 
01767     data     = cpl_imagelist_collapse_minmax_create(datlist, low, high);
01768     variance = cpl_imagelist_collapse_minmax_create(varlist, low, high);
01769 
01770     cpl_image_divide_scalar(variance, N);
01771 
01772     cleanup;
01773     return fors_image_new(data, variance);
01774 }
01775 
01787 fors_image *
01788 fors_image_collapse_median_create(const fors_image_list *images)
01789 {
01790     cpl_imagelist *datlist = NULL;
01791     cpl_imagelist *varlist = NULL;
01792     cpl_image *data = NULL;
01793     cpl_image *variance = NULL;
01794     const fors_image *i;
01795     int N = 0;
01796 
01797     assure( images != NULL, return NULL, NULL );
01798     assure( fors_image_list_size(images) > 0, return NULL, 
01799             "Cannot stack zero images");
01800 
01801     i = fors_image_list_first_const(images);
01802     
01803     datlist = cpl_imagelist_new();
01804     varlist = cpl_imagelist_new();
01805     while(i != NULL) {
01806         /* Append to image lists */
01807         cpl_imagelist_set(datlist, 
01808                           cpl_image_duplicate(i->data), 
01809                           cpl_imagelist_get_size(datlist));
01810         cpl_imagelist_set(varlist,
01811                           cpl_image_duplicate(i->variance),
01812                           cpl_imagelist_get_size(varlist));
01813 
01814         i = fors_image_list_next_const(images);
01815         N++;
01816     }
01817     
01818 #ifdef CPL_IS_NOT_CRAP
01819     data    = cpl_imagelist_collapse_median_create(datlist);
01820 
01821     variance = cpl_imagelist_collapse_create(varlist);
01822 #else
01823     data    = fors_imagelist_collapse_median_create(datlist);
01824 
01825     variance = fors_imagelist_collapse_create(varlist);
01826 #endif
01827 
01828     cpl_image_divide_scalar(variance, N);
01829 
01830     cpl_image_multiply_scalar(variance, 
01831                   fors_utils_median_corr(N) * 
01832                   fors_utils_median_corr(N));
01833     
01834     cleanup;
01835     return fors_image_new(data, variance);
01836 }
01837 
01838 #undef cleanup
01839 #define cleanup
01840 
01857 void fors_image_draw(fors_image *image, int type,
01858              double x, double y,
01859              int radius, double color)
01860 {
01861     assure( image != NULL, return, NULL );
01862 
01863     assure( type == 0 || type == 1 || type == 2,
01864             return , "Unsupported type %d", type);
01865 
01866     assure( radius > 0, return, NULL );
01867 
01868     if (type == 2) {
01869         int i;
01870         for (i = 0; i < 360; i++) {
01871             /* Step size of 1 degree is arbitrary */
01872 
01873             int px = x + radius*cos(i/(2*M_PI));
01874             int py = y + radius*sin(i/(2*M_PI));
01875             
01876             if (1 <= px && px <= cpl_image_get_size_x(image->data) &&
01877                 1 <= py && py <= cpl_image_get_size_y(image->data)) {
01878                 cpl_image_set(image->data, px, py, color);
01879                 cpl_image_set(image->variance, px, py, color > 0 ? color : 0);
01880             }
01881         }
01882     }
01883     else {
01884         int i;
01885 
01886         for (i = -radius; i <= radius; i++) {
01887 
01888             int px, py;
01889             
01890             if (type == 0) {
01891                 px = x + i;
01892                 py = y;
01893             }
01894             else {
01895                 px = x;
01896                 py = y + i;
01897             }
01898             
01899             if (1 <= px && px <= cpl_image_get_size_x(image->data) &&
01900                 1 <= py && py <= cpl_image_get_size_y(image->data)) {
01901                 cpl_image_set(image->data    , px, py, color);
01902                 cpl_image_set(image->variance, px, py, color > 0 ? color : 0);
01903             }
01904         }
01905     }
01906 
01907     return;
01908 }
01909 
01910 #define LIST_DEFINE
01911 #undef LIST_ELEM
01912 #define LIST_ELEM fors_image
01913 #include <list.h>
01914 

Generated on Wed Jun 24 14:11:17 2009 for FORS Pipeline Reference Manual by  doxygen 1.4.7