vircam_imcombine.c

00001 /* $Id: vircam_imcombine.c,v 1.24 2007/05/14 15:26:03 jim Exp $
00002  *
00003  * This file is part of the VIRCAM Pipeline
00004  * Copyright (C) 2005 Cambridge Astronomy Survey Unit
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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00019  */
00020 
00021 /*
00022  * $Author: jim $
00023  * $Date: 2007/05/14 15:26:03 $
00024  * $Revision: 1.24 $
00025  * $Name:  $
00026  */
00027 
00028 /* Includes */
00029 
00030 #ifdef HAVE_CONFIG_H
00031 #include <config.h>
00032 #endif
00033 
00034 #include <cpl.h>
00035 
00036 #include "vircam_utils.h"
00037 #include "vircam_pfits.h"
00038 #include "vircam_stats.h"
00039 #include "vircam_fits.h"
00040 #include "vircam_mods.h"
00041 
00042 /* Macro definitions */
00043 
00044 #define FATAL_ERR(_f,_a) {cpl_msg_error(_f,_a); tidy(); *status = VIR_FATAL; return(*status);}
00045 #define WARN_ERR(_f,_a) {cpl_msg_error(_f,_a); tidy(); *status = VIR_WARN; return(VIR_WARN);}
00046 #define INFO_ERR(_f,_a) {cpl_msg_info(_f,_a);}
00047 #define MEDIANCALC 1
00048 #define MEANCALC 2
00049 #define SZBUF 1024
00050 
00051 /* Definition of litestruct which is used to hold useful information about
00052    each of the input frames and their associated images */
00053 
00054 typedef struct {
00055     vir_fits         *frame;
00056     float            exptime;
00057     float            expfudge;
00058     float            skylevel;
00059     float            skynoise;
00060     float            skyfudge;
00061     float            skyrenorm;
00062 } litestruct;
00063 
00064 /* Static Global variables */
00065 
00066 static litestruct *fileptrs = NULL;
00067 static float **datas = NULL;
00068 static float *odata;
00069 static long npts;
00070 static int nf;
00071 static float oskylevel;
00072 static long nx;
00073 static long ny;
00074 static unsigned char *rmask = NULL;
00075 static unsigned char *rplus = NULL;
00076 
00077 /* Static subroutine prototypes */
00078 
00079 static void skyest(float *data, float thresh, float *skymed, float *skynoise);
00080 static void medcalc(float, float, int);
00081 static void meancalc(float, float, int);
00082 static void xclip_med(float thresh, int scaletype);
00083 static void xclip_mean(float thresh, int scaletype);
00084 
00085 static void tidy(void);
00086 
00089 /*---------------------------------------------------------------------------*/
00149 /*---------------------------------------------------------------------------*/
00150 
00151 extern int vircam_imcombine(vir_fits **fset, int nfits, int combtype, 
00152                             int scaletype, int xrej, float thresh,
00153                             cpl_image **outimage, unsigned char **rejmask,
00154                             unsigned char **rejplus, cpl_propertylist **drs,
00155                             int *status) {
00156     int i,k,j;
00157     const char *ic_fctid = "vircam_imcombine";
00158     char msg[SZBUF];
00159     float sumsky,sumsig,exp1,exp2,expfudge,skylevel,skynoise,oskynoise;
00160     float *dat;
00161     litestruct *ff;
00162     cpl_propertylist *plist_p;
00163 
00164     /* Inherited status */
00165 
00166     *rejmask = NULL;
00167     *rejplus = NULL;
00168     *drs = NULL;
00169     *outimage = NULL;
00170     if (*status != VIR_OK) 
00171         return(*status);
00172 
00173     /* Check that there are any files in the first place...*/ 
00174 
00175     nf = nfits;
00176     if (nf == 0) 
00177         WARN_ERR(ic_fctid,"No files to combine")
00178 
00179     /* Get some file structures */
00180 
00181     fileptrs = cpl_calloc(nf,sizeof(litestruct));
00182 
00183     /* Get workspace for convenience array */
00184 
00185     datas = cpl_malloc(nf*sizeof(float*));
00186     npts = vircam_getnpts(vircam_fits_get_image(fset[0]));
00187     for (k = 0; k < nf; k++)
00188         datas[k] = cpl_malloc(npts*sizeof(float));
00189 
00190     /* Get pointers to the data arrays */
00191 
00192     for (k = 0; k < nf; k++) {
00193 
00194         dat = cpl_image_get_data_float(vircam_fits_get_image(fset[k]));
00195         if (dat == NULL) {
00196             snprintf(msg,SZBUF,"Failed to load data from extension %d in %s",
00197                      vircam_fits_get_nexten(fset[k]),
00198                      vircam_fits_get_filename(fset[k]));
00199             FATAL_ERR(ic_fctid,msg)
00200         }
00201         for (i = 0; i < npts; i++)
00202             datas[k][i] = dat[i];
00203     }
00204 
00205     /* Open each file in turn and fill in the necessary information. Start with
00206        the file name...*/
00207 
00208     sumsky = 0.0;
00209     sumsig = 0.0;
00210     for (i = 0; i < nf; i++) {
00211         ff = fileptrs + i;
00212         ff->frame = fset[i];
00213 
00214         /* If this is the first frame, then keep the size of the data
00215            array for future reference */
00216 
00217         if (i == 0) {
00218             nx = cpl_image_get_size_x(vircam_fits_get_image(fset[0]));
00219             ny = cpl_image_get_size_y(vircam_fits_get_image(fset[0]));
00220             npts = nx*ny;
00221         }
00222 
00223         /* Get the header and the exposure time */
00224 
00225         plist_p = vircam_fits_get_phu(ff->frame);
00226         if (plist_p == NULL) {
00227             snprintf(msg,SZBUF,"Failed to load primary property list %s",
00228                      vircam_fits_get_filename(ff->frame));
00229             INFO_ERR(ic_fctid,msg)
00230             exp2 = 1.0;
00231         } else {
00232             if (vircam_pfits_get_exptime(plist_p,&exp2) != VIR_OK) {
00233                 exp2 = 1.0;
00234                 snprintf(msg,SZBUF,"Failed to get exposure time from %s",
00235                          vircam_fits_get_filename(ff->frame));
00236                 INFO_ERR(ic_fctid,msg)
00237             }
00238         }
00239 
00240         /* Set a few properties */
00241 
00242         ff->exptime = exp2;
00243         exp1 = fileptrs->exptime;
00244         expfudge = exp1/exp2;
00245         ff->expfudge = expfudge;
00246 
00247         /* If scaling by relative exposure time, then do it now. NB: This
00248            isn't necessary for the first file as all the others are scaled
00249            relative to it */
00250 
00251         if (scaletype == 3 && i > 0) {
00252             for (j = 0; j < npts; j++) 
00253                 datas[i][j] *= ff->expfudge;
00254         }
00255 
00256         /* Get the background estimate and noise */
00257 
00258         skyest(datas[i],thresh,&skylevel,&skynoise);
00259         ff->skylevel = skylevel;
00260         ff->skynoise = skynoise;
00261         sumsky += skylevel;
00262         sumsig += skynoise;
00263     }
00264 
00265     /* Work out average background and noise. Then create background zeropoint
00266        or scale factor, depending upon which was requested in the call */
00267 
00268     sumsky /= (float)nf;
00269     sumsig /= (float)nf;
00270     switch (scaletype) {
00271     case 1:
00272         for (i = 0; i < nf; i++)
00273             (fileptrs+i)->skyfudge = sumsky - (fileptrs+i)->skylevel;
00274         break;
00275     case 2:
00276         for (i = 0; i < nf; i++)
00277             (fileptrs+i)->skyfudge = sumsky/(fileptrs+i)->skylevel;
00278         break;      
00279     case 3:
00280         for (i = 0; i < nf; i++)
00281             (fileptrs+i)->skyfudge = sumsky - (fileptrs+i)->skylevel;
00282         break;
00283     default:
00284         for (i = 0; i < nf; i++) 
00285             (fileptrs+i)->skyfudge = 0.0;
00286         break;
00287     }
00288 
00289     /* Open an output image based on the first frame in the frameset */
00290 
00291     *outimage = cpl_image_new(nx,ny,CPL_TYPE_FLOAT);
00292     odata = cpl_image_get_data_float(*outimage);
00293     if (*outimage == NULL || odata == NULL) 
00294         FATAL_ERR(ic_fctid,"Couldn't create output image")
00295     *rejmask = cpl_calloc(npts,sizeof(*rejmask));
00296     rmask = *rejmask;
00297     *rejplus = cpl_calloc(npts,sizeof(*rejplus));
00298     rplus = *rejplus;
00299 
00300     /* Now do the averaging/medianing */
00301 
00302     switch (combtype) {
00303     case MEDIANCALC:
00304         medcalc(thresh,sumsig,scaletype);
00305         break;
00306     case MEANCALC:
00307         meancalc(thresh,sumsig,scaletype);
00308         break;
00309     }
00310 
00311     /* Do the extra clipping here if you want it */
00312 
00313     if (xrej) {
00314         
00315         /* First get sky background and sigma from output data */
00316 
00317         skyest(odata,thresh,&oskylevel,&oskynoise);
00318 
00319         /* Now loop for all the files subtract off the mean frame (suitably
00320            scaled and zero pointed depending on what was done in the first 
00321            place) */
00322 
00323         for (i = 0; i < nf; i++) {
00324             ff = fileptrs + i;
00325             ff->skyrenorm = ff->skylevel/oskylevel;
00326             switch (scaletype) {
00327             case 1:
00328                 for (k = 0; k < npts; k++) 
00329                     datas[i][k] -= (odata[k] - oskylevel);
00330                 break;
00331             case 2:
00332                 for (k = 0; k < npts; k++)
00333                     datas[i][k] -= (odata[k] - oskylevel)*ff->skyrenorm;
00334                 break;
00335             case 3:
00336                 for (k = 0; k < npts; k++) 
00337                     datas[i][k] -= (odata[k] - oskylevel);
00338                 break;
00339             case 0:
00340                 for (k = 0; k < npts; k++) 
00341                     datas[i][k] -= (odata[k] - oskylevel);
00342                 break;
00343             }
00344 
00345             /* Re-estimate the noise for this image */
00346 
00347             skyest(datas[i],thresh,&skylevel,&skynoise);
00348             ff->skynoise = skynoise;
00349         }
00350 
00351         /* Now do the extra clip... */
00352 
00353         switch (combtype) {
00354         case MEDIANCALC:
00355             xclip_med(thresh,scaletype);
00356             break;
00357         case MEANCALC:
00358             xclip_mean(thresh,scaletype);
00359             break;
00360         }
00361     }
00362 
00363     /* Write provenance keywords */
00364 
00365     *drs = cpl_propertylist_new();
00366     vircam_prov(*drs,fset,nfits);
00367 
00368     /* Right, tidy and get out of here */
00369 
00370     tidy();
00371     return(VIR_OK);
00372 }
00373 
00374 /*---------------------------------------------------------------------------*/
00396 /*---------------------------------------------------------------------------*/
00397 
00398 static void xclip_med(float thresh, int scaletype) {
00399     int nf1,nf2,nfm,nrejmax,is_even,k,is_even2,nrej,nremain,nm,nmm,nplus;
00400     int nminus;
00401     float **work,*dork,value,cliplev;
00402     long i;
00403     int j;
00404     litestruct *ff;
00405 
00406     /* Set up a few useful variables */
00407 
00408     nf1 = nf/2 - 1;
00409     nf2 = nf1 + 1;
00410     nfm = (nf + 1)/2 - 1;
00411     nrejmax = nf/2;
00412     is_even = !(nf & 1);
00413 
00414     /* Get some workspace */
00415 
00416     work = cpl_malloc(3*sizeof(float *));
00417     for (i = 0; i < 3; i++)
00418         work[i] = cpl_malloc(nf*sizeof(float));
00419     dork = cpl_malloc(nf*sizeof(float));
00420 
00421     /* Loop for each input pixel now... */
00422 
00423     for (i = 0; i < npts; i++) {
00424 
00425         /* Scale or shift data */
00426 
00427         switch (scaletype) {
00428         case 0:
00429             for (k = 0; k < nf; k++) {
00430                 ff = fileptrs + k;
00431                 work[0][k] = datas[k][i];
00432                 work[1][k] = ff->skynoise;
00433                 work[2][k] = datas[k][i] + odata[i] - oskylevel;
00434             }
00435             break;
00436         case 1:
00437             for (k = 0; k < nf; k++) {
00438                 ff = fileptrs + k;
00439                 work[0][k] = datas[k][i] + ff->skyfudge;
00440                 work[1][k] = ff->skynoise;
00441                 work[2][k] = datas[k][i] + odata[i] - oskylevel + ff->skyfudge;
00442             }
00443             break;
00444         case 2:
00445             for (k = 0; k < nf; k++) {
00446                 ff = fileptrs + k;
00447                 work[0][k] = datas[k][i]*ff->skyfudge;
00448                 work[1][k] = ff->skynoise*ff->skyfudge;
00449                 work[2][k] = (datas[k][i] + odata[i]*ff->skyrenorm - 
00450                               ff->skylevel)*ff->skyfudge;
00451             }
00452             break;
00453         case 3:
00454             for (k = 0; k < nf; k++) {
00455                 ff = fileptrs + k;
00456                 work[0][k] = datas[k][i] + ff->skyfudge;
00457                 work[1][k] = ff->skynoise;
00458                 work[2][k] = datas[k][i] + odata[i] - oskylevel + ff->skyfudge;
00459             }
00460             break;
00461         }       
00462 
00463         /* Sort and get a first pass median */
00464 
00465         vircam_sort(work,nf,3);
00466         if (is_even)
00467             value = 0.5*(work[0][nf1] + work[0][nf2]);
00468         else 
00469             if (nf < 5) 
00470                 value = work[0][nfm];
00471             else 
00472                 value = 0.25*(work[0][nfm-1] + work[0][nfm+1]) + 0.5*work[0][nfm];
00473         
00474         /* Do clipping */
00475 
00476         nplus = 0;
00477         cliplev = value + thresh*work[1][nf-1];
00478         while (nplus < nrejmax && work[0][nf-nplus-1] > cliplev)
00479             nplus++;
00480         nminus = 0;
00481         cliplev = value - thresh*work[1][nf-1];
00482         while ((nplus+nminus) < nrejmax && work[0][nminus] < cliplev)
00483             nminus++;
00484         nrej = nplus + nminus;
00485 
00486         /* If there were any clipped out, the re-estimate the value */
00487 
00488         if (nrej > 0) {
00489             nremain = nf - nrej;
00490             if (nremain != 0) {
00491                 nm = nremain/2 - 1;
00492                 for (j = 0; j < nremain; j++) 
00493                     dork[j] = work[2][j+nminus];
00494                 nmm = (nremain + 1)/2 - 1;
00495                 is_even2 = !(nremain & 1);
00496                 vircam_sort(&dork,nm,1);
00497                 if (is_even2) 
00498                     value = 0.5*(dork[nm] + dork[nm+1]);
00499                 else 
00500                     if (nremain < 3) 
00501                         value = dork[nmm];
00502                     else 
00503                         value = 0.5*dork[nmm] + 0.25*(dork[nmm-1] + dork[nmm+1]);
00504             }
00505         
00506             /* Store the result away */
00507 
00508             odata[i] = value;
00509             rmask[i] = min(255,nrej);
00510             rplus[i] = min(255,nplus);
00511         } else {
00512             rmask[i] = 0;
00513             rplus[i] = 0;
00514         }
00515     }
00516 
00517     /* Ditch workspace and get out of here */
00518 
00519     for (i = 0; i < 3; i++)
00520         cpl_free(work[i]);
00521     cpl_free(work);
00522     cpl_free(dork);
00523 }
00524 
00525 /*---------------------------------------------------------------------------*/
00547 /*---------------------------------------------------------------------------*/
00548 
00549 static void xclip_mean(float thresh, int scaletype) {
00550     int k,nf2,nrej,nplus;
00551     float *work[3],value,cliplev1,cliplev2,value2;
00552     long i;
00553     litestruct *ff;
00554 
00555     /* Get some workspace */
00556 
00557     for (i = 0; i < 3; i++)
00558         work[i] = cpl_malloc(nf*sizeof(float));
00559 
00560     /* Loop for each input pixel now... */
00561 
00562     for (i = 0; i < npts; i++) {
00563 
00564         /* Scale or shift data */
00565 
00566         switch (scaletype) {
00567         case 0:
00568             for (k = 0; k < nf; k++) {
00569                 ff = fileptrs + k;
00570                 work[0][k] = datas[k][i];
00571                 work[1][k] = ff->skynoise;
00572                 work[2][k] = datas[k][i] + odata[i] - oskylevel;
00573             }
00574             break;
00575         case 1:
00576             for (k = 0; k < nf; k++) {
00577                 ff = fileptrs + k;
00578                 work[0][k] = datas[k][i] + ff->skyfudge;
00579                 work[1][k] = ff->skynoise;
00580                 work[2][k] = datas[k][i] + odata[i] - oskylevel + ff->skyfudge;
00581             }
00582             break;
00583         case 2:
00584             for (k = 0; k < nf; k++) {
00585                 ff = fileptrs + k;
00586                 work[0][k] = datas[k][i]*ff->skyfudge;
00587                 work[1][k] = ff->skynoise*ff->skyfudge;
00588                 work[2][k] = (datas[k][i] + odata[i]*ff->skyrenorm - 
00589                            ff->skylevel)*ff->skyfudge;
00590             }
00591             break;
00592         case 3:
00593             for (k = 0; k < nf; k++) {
00594                 ff = fileptrs + k;
00595                 work[0][k] = datas[k][i] + ff->skyfudge;
00596                 work[1][k] = ff->skynoise;
00597                 work[2][k] = datas[k][i] + odata[i] - oskylevel + ff->skyfudge;
00598             }
00599             break;
00600         }       
00601 
00602         /* Get a first pass mean */
00603 
00604         value = 0.0;
00605         for (k = 0; k < nf; k++) 
00606             value += work[0][k];
00607         value /= (float)nf;
00608         
00609         /* Do upper clipping */
00610 
00611         value2 = 0.0;
00612         nplus = 0;
00613         nf2 = 0;
00614         for (k = 0; k < nf; k++) {
00615             cliplev1 = value - thresh*work[1][k];
00616             cliplev2 = value + thresh*work[1][k];
00617             if (work[0][k] >= cliplev1 && work[0][k] <= cliplev2) {
00618                 value2 += work[2][k];
00619                 nf2++;
00620             } else if (work[0][k] > cliplev2) {
00621                 nplus++;
00622             }
00623         }
00624 
00625         /* If there were any clipped out, the re-estimate the value */
00626 
00627         if (nf2 < nf) {
00628             if (nf2 != 0)
00629                 value = value2/(float)nf2;
00630         
00631             /* Store the result away */
00632 
00633             odata[i] = value;
00634         }
00635         nrej = nf - nf2;
00636         rmask[i] = min(255,nrej);       
00637         rplus[i] = min(255,nplus);
00638     }
00639 
00640     /* Ditch workspace and get out of here */
00641 
00642     for (k = 0; k < 3; k++)
00643         cpl_free(work[k]);
00644 }
00645 
00646 
00647 /*---------------------------------------------------------------------------*/
00670 /*---------------------------------------------------------------------------*/
00671 
00672 static void medcalc(float thresh, float avskynoise, int scaletype) {
00673     int nf1,nf2,nfm,nrejmax,is_even,nrej,nremain,nm,nmm,is_even2,k,nminus;
00674     int nplus;
00675     long i;
00676     float value,cliplev,*work;
00677 
00678     /* Set up a few useful variables */
00679 
00680     nf1 = nf/2 - 1;
00681     nf2 = nf1 + 1;
00682     nfm = (nf + 1)/2 - 1;
00683     nrejmax = nf/2;
00684     is_even = !(nf & 1);
00685 
00686     /* Get a workspace */
00687 
00688     work = cpl_malloc(nf*sizeof(*work));
00689 
00690     /* Ok, loop for each pixel... */
00691 
00692     for (i = 0; i < npts; i++) {
00693 
00694         /* Scale or shift data */
00695 
00696         switch (scaletype) {
00697         case 0:
00698             for (k = 0; k < nf; k++) 
00699                 work[k] = datas[k][i];
00700             break;
00701         case 1:
00702             for (k = 0; k < nf; k++) 
00703                 work[k] = datas[k][i] + (fileptrs+k)->skyfudge;
00704             break;
00705         case 2:
00706             for (k = 0; k < nf; k++) 
00707                 work[k] = datas[k][i]*(fileptrs+k)->skyfudge;
00708             break;
00709         case 3:
00710             for (k = 0; k < nf; k++) 
00711                 work[k] = datas[k][i] + (fileptrs+k)->skyfudge;
00712             break;
00713         }
00714         
00715         /* Sort data and get the median */
00716 
00717         vircam_sort(&work,nf,1);
00718         if (is_even)
00719             value = 0.5*(work[nf1] + work[nf2]);
00720         else 
00721             if (nf < 5) 
00722                 value = work[nfm];
00723             else 
00724                 value = 0.25*(work[nfm-1] + work[nfm+1]) + 0.5*work[nfm];
00725 
00726         /* Enter a rejection loop */
00727 
00728         nplus = 0;
00729         cliplev = value + thresh*avskynoise;
00730         while (nplus < nrejmax && work[nf-nplus-1] > cliplev)
00731             nplus++;
00732         nminus = 0;
00733         cliplev = value - thresh*avskynoise;
00734         while ((nplus+nminus) < nrejmax && work[nminus] < cliplev)
00735             nminus++;
00736         nrej = nplus + nminus;
00737 
00738         /* If you've clipped any, then recalculate the median...*/
00739 
00740         if (nrej > 0) {
00741             nremain = nf - nrej;
00742             nm = nremain/2 - 1 + nminus;
00743             nmm = (nremain + 1)/2 - 1 + nminus;
00744             is_even2 = !(nremain & 1);
00745             if (is_even2) 
00746                 value = 0.5*(work[nm] + work[nm+1]);
00747             else 
00748                 if (nremain < 3) 
00749                     value = work[nmm];
00750                 else 
00751                     value = 0.5*work[nmm] + 0.25*(work[nmm-1] + work[nmm+1]);
00752         }
00753 
00754         /* Store the result away */
00755 
00756         odata[i] = value;
00757         rmask[i] = min(255,nrej);
00758         rplus[i] = min(255,nplus);
00759     }
00760 
00761     /* Get rid of workspace */
00762 
00763     cpl_free(work);
00764 }
00765 
00766 /*---------------------------------------------------------------------------*/
00789 /*---------------------------------------------------------------------------*/
00790 
00791 static void meancalc(float thresh, float avskynoise, int scaletype) {
00792     int nf2,k,nrej,nplus;
00793     long i;
00794     float *work,value,cliplev1,cliplev2,value2;
00795 
00796     /* Get a vector for workspace */
00797 
00798     work = cpl_malloc(nf*sizeof(*work));
00799 
00800     /* Ok, loop for each pixel... */
00801 
00802     for (i = 0; i < npts; i++) {
00803 
00804         /* Scale or shift data */
00805 
00806         switch (scaletype) {
00807         case 0:
00808             for (k = 0; k < nf; k++) 
00809                 work[k] = datas[k][i];
00810             break;
00811         case 1:
00812             for (k = 0; k < nf; k++) 
00813                 work[k] = datas[k][i] + (fileptrs+k)->skyfudge;
00814             break;
00815         case 2:
00816             for (k = 0; k < nf; k++) 
00817                 work[k] = datas[k][i]*(fileptrs+k)->skyfudge;
00818             break;
00819         case 3:
00820             for (k = 0; k < nf; k++) 
00821                 work[k] = datas[k][i] + (fileptrs+k)->skyfudge;
00822             break;
00823         }
00824         
00825         /* Get the mean */
00826         
00827         value = 0.0;
00828         for (k = 0; k < nf; k++) 
00829             value += work[k];
00830         value /= (float)nf;
00831 
00832         /* Enter a rejection loop  */
00833 
00834         value2 = 0.0;
00835         nf2 = 0;
00836         nplus = 0;
00837         cliplev1 = value - thresh*avskynoise;
00838         cliplev2 = value + thresh*avskynoise;
00839         for (k = 0; k < nf; k++) {
00840             if (work[k] >= cliplev1 && work[k] <= cliplev2) {
00841                 value2 += work[k];
00842                 nf2 += 1;
00843             } else if (work[k] > cliplev2) {
00844                 nplus++;
00845             }
00846         }
00847 
00848         /* If you've clipped any, then recalculate the mean...*/
00849 
00850         if (nf2 < nf && nf2 != 0)
00851             value = value2/(float)nf2;
00852         nrej = nf - nf2;
00853 
00854         /* Store the result away */
00855 
00856         odata[i] = value;
00857         rmask[i] = min(255,nrej);
00858         rplus[i] = min(255,nplus);
00859     }
00860 
00861     /* Get rid of workspace */
00862 
00863     cpl_free(work);
00864 }
00865             
00866 /*---------------------------------------------------------------------------*/
00889 /*---------------------------------------------------------------------------*/
00890 
00891 static void skyest(float *data, float thresh, float *skymed, float *skynoise) {
00892     unsigned char *bpm;
00893 
00894     /* Get a dummy bad pixel mask */
00895 
00896     bpm = cpl_calloc(npts,sizeof(*bpm));
00897 
00898     /* Get the stats */
00899 
00900     vircam_qmedsig(data,bpm,npts,thresh,3,-65535.0,65535.0,skymed,
00901                    skynoise);
00902 
00903     /* Clean up */
00904 
00905     cpl_free(bpm);
00906 }
00907 
00908 /*---------------------------------------------------------------------------*/
00912 /*---------------------------------------------------------------------------*/
00913 
00914 static void tidy(void) {
00915     int i;
00916 
00917     /* Free up work space associated with file structures */
00918 
00919     freespace(fileptrs);
00920     for (i = 0; i < nf; i++)
00921         freespace(datas[i]);
00922     freespace(datas);
00923 }
00924 
00928 /*
00929 
00930 $Log: vircam_imcombine.c,v $
00931 Revision 1.24  2007/05/14 15:26:03  jim
00932 Fixed bug in multiplicative scaling
00933 
00934 Revision 1.23  2007/03/29 12:19:39  jim
00935 Little changes to improve documentation
00936 
00937 Revision 1.22  2007/03/01 12:42:41  jim
00938 Modified slightly after code checking
00939 
00940 Revision 1.21  2006/10/02 13:47:33  jim
00941 Added missing .h file to include list
00942 
00943 Revision 1.20  2006/09/08 09:21:37  jim
00944 Fixed bug in median routine
00945 
00946 Revision 1.19  2006/08/21 09:08:10  jim
00947 Fixed a few subtle bugs in the median routines
00948 
00949 Revision 1.18  2006/07/04 09:19:05  jim
00950 replaced all sprintf statements with snprintf
00951 
00952 Revision 1.17  2006/05/08 10:56:14  jim
00953 Fixed bug where exposure time scaling was not happening to the data copy
00954 
00955 Revision 1.16  2006/04/20 11:26:45  jim
00956 Fixed so that original data isn't modified in the input vir_fits list.
00957 
00958 Revision 1.15  2006/03/23 21:18:48  jim
00959 Minor changes mainly to comment headers
00960 
00961 Revision 1.14  2006/03/22 13:58:31  jim
00962 Cosmetic fixes to keep lint happy
00963 
00964 Revision 1.13  2006/03/08 14:32:21  jim
00965 Lots of little modifications
00966 
00967 Revision 1.12  2006/03/03 14:29:46  jim
00968 Modified definition of vir_fits and channel table
00969 
00970 Revision 1.11  2006/03/01 10:31:28  jim
00971 Now uses new vir_fits objects
00972 
00973 Revision 1.10  2006/02/22 10:09:09  jim
00974 Added status variable to call
00975 
00976 Revision 1.9  2006/01/23 10:30:49  jim
00977 Mainly documentation mods
00978 
00979 Revision 1.8  2005/12/14 22:17:33  jim
00980 Updated docs
00981 
00982 Revision 1.7  2005/11/29 11:54:19  jim
00983 Modified call in skyest to allow for larger negative numbers as it is
00984 likely that the reset frames will be negative.
00985 
00986 Revision 1.6  2005/11/25 09:56:15  jim
00987 Tidied up some more documentation
00988 
00989 Revision 1.5  2005/11/08 12:47:44  jim
00990 Made garbage collection a little better
00991 
00992 Revision 1.4  2005/11/07 13:15:16  jim
00993 Fixed lots of bugs and added some error checking
00994 
00995 Revision 1.3  2005/11/03 13:28:49  jim
00996 All sorts of changes to tighten up error handling
00997 
00998 Revision 1.2  2005/10/14 13:17:54  jim
00999 Now spews out a mask of rejected pixels and rejected positive pixels
01000 (indicative of cosmic ray hits)
01001 
01002 Revision 1.1.1.1  2005/08/05 08:29:09  jim
01003 Initial import
01004 
01005 
01006 */
01007 

Generated on Wed Apr 10 04:01:55 2013 for VIRCAM Pipeline by  doxygen 1.5.1