vircam_utils.c

00001 /* $Id: vircam_utils.c,v 1.56 2008/07/10 13:04: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: 2008/07/10 13:04:03 $
00024  * $Revision: 1.56 $
00025  * $Name:  $
00026  */
00027 
00028 /* Includes */
00029 
00030 #ifdef HAVE_CONFIG_H
00031 #include <config.h>
00032 #endif
00033 
00034 #include <sys/time.h>
00035 #include <time.h>
00036 #include <libgen.h>
00037 #include <string.h>
00038 #include <unistd.h>
00039 
00040 #include <cpl.h>
00041 #include <math.h>
00042 
00043 #include "vircam_utils.h"
00044 #include "vircam_stats.h"
00045 #include "vircam_fits.h"
00046 #include "catalogue/imcore.h"
00047 
00048 #define SZKEY 32
00049 #define SZVAL 64
00050 
00051 /* Define some columns for illumination correction tables */
00052 
00053 #define NI_COLS 5
00054 static const char *illcor_cols[NI_COLS] = {"xmin","xmax","ymin","ymax",
00055                                            "illcor"};
00056 
00057 /* Static subroutine prototypes */
00058 
00059 static float madfunc(int npts, float *xt, float *yt, float b);
00060 
00074 /*---------------------------------------------------------------------------*/
00089 /*---------------------------------------------------------------------------*/
00090 
00091 extern const char *vircam_get_license(void) {
00092     const char  *vircam_license = 
00093         "This file is part of the VIRCAM Instrument Pipeline\n"
00094         "Copyright (C) 2006 Cambridge Astronomy Survey Unit\n"
00095         "\n"
00096         "This program is free software; you can redistribute it and/or modify\n"
00097         "it under the terms of the GNU General Public License as published by\n"
00098         "the Free Software Foundation; either version 2 of the License, or\n"
00099         "(at your option) any later version.\n"
00100         "\n"
00101         "This program is distributed in the hope that it will be useful,\n"
00102         "but WITHOUT ANY WARRANTY; without even the implied warranty of\n"
00103         "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\n"
00104         "GNU General Public License for more details.\n"
00105         "\n"
00106         "You should have received a copy of the GNU General Public License\n"
00107         "along with this program; if not, write to the Free Software\n"
00108         "Foundation, Inc., 59 Temple Place, Suite 330, Boston, \n"
00109         "MA  02111-1307  USA";
00110     return(vircam_license);
00111 }
00112 
00113 
00114 /*---------------------------------------------------------------------------*/
00138 /*---------------------------------------------------------------------------*/
00139 
00140 extern int vircam_compare_tags(const cpl_frame *frame1, 
00141                                const cpl_frame *frame2) {
00142     char *v1,*v2;
00143 
00144     /* Test entries */
00145 
00146     if (frame1 == NULL || frame2 == NULL) 
00147         return(-1);
00148 
00149     /* Get the tags */
00150 
00151     if ((v1 = (char *)cpl_frame_get_tag(frame1)) == NULL) 
00152         return(-1);
00153     if ((v2 = (char *)cpl_frame_get_tag(frame2)) == NULL) 
00154         return(-1);
00155 
00156     /* Compare the tags */
00157 
00158     if (strcmp(v1,v2)) 
00159         return(0);
00160     else 
00161         return(1);
00162 }
00163 
00164 /*---------------------------------------------------------------------------*/
00190 /*---------------------------------------------------------------------------*/
00191 
00192 extern cpl_frameset *vircam_frameset_subgroup(cpl_frameset *frameset, 
00193                                               int *labels, int nlab, 
00194                                               const char *tag) {
00195     int i;
00196     cpl_frameset *cur_set,*ret_set;
00197     cpl_frame *cur_frame;
00198     char *cur_tag;
00199 
00200     ret_set = NULL;
00201     for (i = 0; i < nlab; i++) {
00202         cur_set = cpl_frameset_extract(frameset,labels,i);
00203         if (cur_set == NULL)
00204             break;
00205         cur_frame = cpl_frameset_get_frame(cur_set,0);
00206         cur_tag = (char *)cpl_frame_get_tag(cur_frame);
00207         if (!strcmp(cur_tag,tag)) {
00208             ret_set = cur_set;
00209             break;
00210         }
00211         cpl_frameset_delete(cur_set);
00212     }
00213     return(ret_set);
00214 }
00215 
00216 /*---------------------------------------------------------------------------*/
00243 /*---------------------------------------------------------------------------*/
00244 
00245 extern cpl_frame *vircam_frameset_subgroup_1(cpl_frameset *frameset, 
00246                                              int *labels, int nlab, 
00247                                              const char *tag) {
00248     cpl_frameset *cur_set;
00249     cpl_frame *cur_frame,*new_frame;
00250 
00251     if ((cur_set = vircam_frameset_subgroup(frameset,labels,nlab,tag)) == NULL) {
00252         return(NULL);
00253     } else {
00254         cur_frame = cpl_frameset_get_frame(cur_set,0);
00255         new_frame = cpl_frame_duplicate(cur_frame);
00256         cpl_frameset_delete(cur_set);
00257         return(new_frame);
00258     }
00259 }
00260 
00261 
00262 /*---------------------------------------------------------------------------*/
00287 /*---------------------------------------------------------------------------*/
00288 
00289 extern void vircam_exten_range(int inexten, const cpl_frame *fr, int *out1, 
00290                                int *out2) {
00291     int nvircam = 16,n,nmax;
00292     char *fctid = "vircam_exten_range";
00293 
00294     /* Right, how many frames actually exist? */
00295 
00296     n = cpl_frame_get_nextensions(fr);
00297 
00298     /* If there are less than the cannonical number, then issue a warning
00299        message here */
00300 
00301     if (n < nvircam) {
00302         cpl_msg_warning(fctid,"Only %d extensions out of %d are present",
00303                         n,nvircam);
00304         nmax = n;
00305     } else {
00306         nmax = nvircam;
00307     }
00308 
00309     /* If the extension number requested is zero, then do all the extensions
00310        that are available */
00311     
00312     if (inexten == 0) {
00313         *out1 = 1;
00314         *out2 = nmax;
00315 
00316     /* Otherwise, check that the requested extension is actually present */
00317 
00318     } else {
00319 
00320         if (inexten > nmax) {
00321             cpl_msg_error(fctid,"Requested extension %d is not present",
00322                           inexten);
00323             *out1 = -1;
00324             *out2 = -1;
00325         } else {
00326             *out1 = inexten;
00327             *out2 = inexten;
00328         }
00329     }
00330     return;
00331 }
00332 
00333 /*---------------------------------------------------------------------------*/
00359 /*---------------------------------------------------------------------------*/
00360 
00361 extern void vircam_madfit(int npts, float *xdata, float *ydata, 
00362                           float *intercept, float *slope) {
00363     int j;
00364     float sx,sy,sxx,sxy,det,aa,bb,temp,chisq,sigb,b1,f1,b2,f2,f;
00365 
00366     /* Do a linear least squares for a first estimate */
00367 
00368     sx = 0.0;
00369     sy = 0.0;
00370     sxx = 0.0;
00371     sxy = 0.0;
00372     for (j = 0; j < npts; j++) {
00373         sx += xdata[j];
00374         sy += ydata[j];
00375         sxx += xdata[j]*xdata[j];
00376         sxy += xdata[j]*ydata[j];
00377     }
00378     det = (float)npts*sxx - sx*sx;
00379     if (det == 0.0) {
00380         *slope = 0.0;
00381         *intercept = 0.0;
00382         return;
00383     }
00384     aa = (sxx*sy - sx*sxy)/det;
00385     bb = ((float)npts*sxy - sx*sy)/det;
00386     chisq = 0.0;
00387     for (j = 0; j < npts; j++) {
00388         temp = ydata[j] - (aa + bb*xdata[j]);
00389         chisq += temp*temp;
00390     }
00391     sigb = sqrt(chisq/det);
00392     if (sigb == 0.0) {
00393         *slope = bb;
00394         *intercept = aa;
00395         return;
00396     }
00397 
00398     /* Now bracket the solution */
00399 
00400     b1 = bb;
00401     f1 = madfunc(npts,xdata,ydata,b1);
00402     b2 = bb + ((f1 > 0.0) ? fabs(3.0*sigb) : -fabs(3.0*sigb));
00403     f2 = madfunc(npts,xdata,ydata,b2);
00404     while (f1*f2 > 0.0) {
00405         bb = 2.0*b2 - b1;
00406         b1 = b2;
00407         f1 = f2;
00408         b2 = bb;
00409         f2 = madfunc(npts,xdata,ydata,b2);
00410     }
00411 
00412     /* Iterate to find the minimum value */
00413 
00414     sigb = 0.01*sigb;
00415     while (fabs(b2 - b1) > sigb) {
00416         bb = 0.5*(b1 + b2);
00417         if (bb == b1 || bb == b2) 
00418             break;
00419         f = madfunc(npts,xdata,ydata,bb);
00420         if (f*f1 >= 0.0) {
00421             f1 = f;
00422             b1 = bb;
00423         } else {
00424             f2 = f;
00425             b2 = bb;
00426         }
00427     }
00428     *intercept = aa;
00429     *slope = bb;
00430 }
00431 
00432 /*---------------------------------------------------------------------------*/
00455 /*---------------------------------------------------------------------------*/
00456 
00457 
00458 static float madfunc(int npts, float *xt, float *yt, float b) {
00459     float *arr,aa,d,sum;
00460     int j;
00461 
00462     arr = cpl_malloc(npts*sizeof(*arr));
00463     for (j = 0; j < npts; j++) 
00464         arr[j] = yt[j] - b*xt[j];
00465     aa = vircam_med(arr,NULL,(long)npts);
00466     sum = 0.0;
00467     for (j = 0; j < npts; j++) {
00468         d = yt[j] - (b*xt[j] + aa);
00469         sum += d > 0.0 ? xt[j] : -xt[j];
00470     }
00471     cpl_free(arr);
00472     return(sum);
00473 }
00474 
00475 /*---------------------------------------------------------------------------*/
00502 /*---------------------------------------------------------------------------*/
00503 
00504 extern void vircam_linfit(int npts, double *xdata, double *ydata, 
00505                           double *intercept, double *slope, double *sig) {
00506     int j;
00507     double sx,sy,sxx,sxy,det,aa,bb,temp,sum,sumsq;
00508 
00509     /* Do a linear least squares */
00510 
00511     sx = 0.0;
00512     sy = 0.0;
00513     sxx = 0.0;
00514     sxy = 0.0;
00515     for (j = 0; j < npts; j++) {
00516         sx += xdata[j];
00517         sy += ydata[j];
00518         sxx += xdata[j]*xdata[j];
00519         sxy += xdata[j]*ydata[j];
00520     }
00521     det = (double)npts*sxx - sx*sx;
00522     if (det == 0.0) {
00523         *slope = 0.0;
00524         *intercept = 0.0;
00525         *sig = 0.0;
00526         return;
00527     }
00528 
00529     /* Intercept and slope */
00530 
00531     aa = (sxx*sy - sx*sxy)/det;
00532     bb = ((double)npts*sxy - sx*sy)/det;
00533 
00534     /* Work out RMS of fit */
00535 
00536     sum = 0.0;
00537     sumsq = 0.0;
00538     for (j = 0; j < npts; j++) {
00539         temp = ydata[j] - (aa + bb*xdata[j]);
00540         sum += temp;
00541         sumsq += temp*temp;
00542     }
00543     sum /= (double)npts;
00544 
00545     /* Set return values */
00546 
00547     *sig = sqrt(sumsq/(double)npts - sum*sum);
00548     *slope = bb;
00549     *intercept = aa;
00550 }
00551 
00552 /*---------------------------------------------------------------------------*/
00576 /*---------------------------------------------------------------------------*/
00577 
00578 extern int vircam_solve_gauss(double **a, double *b, int m) {
00579     double temp,big,pivot,rmax;
00580     int i,iu,j,k,jl,ib,ir;
00581     int l = 0;
00582 
00583     iu = m - 1;
00584     for (i = 0; i < iu; i++) {
00585         big = 0.0;
00586 
00587         /* find largest remaining term in ith column for pivot */
00588 
00589         for (k = i; k < m; k++) {
00590             rmax = fabs(a[i][k]);
00591             if (rmax > big) {
00592                 big = rmax;
00593                 l = k;
00594             }
00595         }
00596 
00597         /* check for non-zero term */
00598 
00599         if (big == 0.0) {
00600             for (ib = 0; ib < m; ib++)
00601                 b[ib] = 0.0;
00602             cpl_msg_error("vircam_solve_gauss","Zero Determinant\n");
00603             return(VIR_FATAL);
00604         }
00605 
00606         if (i != l) {
00607 
00608             /* switch rows */
00609 
00610             for (j = 0; j < m; j++) {
00611                 temp = a[j][i];
00612                 a[j][i] = a[j][l];
00613                 a[j][l] = temp;
00614             }
00615             temp = b[i];
00616             b[i] = b[l];
00617             b[l] = temp;
00618         }
00619 
00620 
00621         /* pivotal reduction */
00622 
00623         pivot = a[i][i];
00624         jl = i+1;
00625 
00626         for (j = jl; j < m; j++) {
00627             temp = a[i][j]/pivot;
00628             b[j] -= temp*b[i];
00629             for (k = i; k < m; k++)
00630                 a[k][j] -= temp*a[k][i];
00631         }
00632     }
00633 
00634     /* back substitution for solution */
00635 
00636     for (i = 0; i < m; i++) {
00637         ir = m - 1 - i;
00638         if (a[ir][ir] != 0.0) {
00639             temp = b[ir];
00640             if (ir != m - 1) {
00641                 for (j = 1; j <= i; j++) {
00642                     k = m - j;
00643                     temp -= a[k][ir]*b[k];
00644                 }
00645             }
00646             b[ir] = temp/a[ir][ir];
00647         } else
00648             b[ir] = 0.0;
00649     }
00650     return(VIR_OK);
00651 }
00652 
00653 /*---------------------------------------------------------------------------*/
00691 /*---------------------------------------------------------------------------*/
00692 
00693 extern int vircam_polyfit(const cpl_array *xarray, const cpl_array *yarray, 
00694                           int ncoefs, int ilim, int niter, float lclip, 
00695                           float hclip, cpl_array **polycf, double *sigfit) {
00696     const char *fctid = "vircam_polyfit";
00697     int npts,iter,i,j,nnew,k,retval,n;
00698     double *xdata,*ydata,*pdata,*res,**a,*b,temp,sum,sumsq,val;
00699     double lcut,hcut;
00700     unsigned char *pm;
00701 
00702     /* Initialise a few things */
00703 
00704     *polycf = NULL;
00705     *sigfit = -1.0;
00706 
00707     /* How many data points do we have? */
00708 
00709     npts = cpl_array_get_size(xarray);
00710 
00711     /* Do we have enough points for the required order of the fit? */
00712 
00713     if (npts < ncoefs) {
00714         cpl_msg_warning(fctid,"Not data for fit, Npts = %d, Ncoefs = %d",
00715                         npts,ncoefs);
00716         return(VIR_FATAL);
00717     }
00718 
00719     /* Create some arrays */
00720 
00721     a = cpl_malloc(ncoefs*sizeof(double *));
00722     b = cpl_calloc(ncoefs,sizeof(double));
00723     for (i = 0; i < ncoefs; i++) 
00724         a[i] = cpl_calloc(ncoefs,sizeof(double));
00725     pm = cpl_calloc(npts,sizeof(unsigned char));
00726     res = cpl_malloc(npts*sizeof(double));
00727 
00728     /* Get pointers to the input arrays */
00729 
00730     xdata = (double *)cpl_array_get_data_double_const(xarray);
00731     ydata = (double *)cpl_array_get_data_double_const(yarray);
00732 
00733     /* Get an output array */
00734 
00735     *polycf = cpl_array_new(ncoefs,CPL_TYPE_DOUBLE);
00736     pdata = cpl_array_get_data_double(*polycf);
00737 
00738     /* Iteration loop */
00739 
00740     for (iter = 0; iter < niter; iter++) {
00741 
00742         /* Zero some accumulators */
00743 
00744         for (i = 0; i < ncoefs; i++)  {
00745             for (j = 0; j < ncoefs; j++) 
00746                 a[i][j] = 0.0;
00747             b[i] = 0.0;
00748         }
00749         nnew = 0;
00750 
00751         /* Cumulate sums */
00752  
00753         for (i = 0; i < npts; i++) {
00754             if (pm[i] == 1)
00755                 continue;
00756             for (k = 0; k < ncoefs; k++) {
00757                 temp = 1.0;
00758                 if (k + ilim != 0)
00759                     temp = pow(xdata[i],(double)(k+ilim));
00760                 b[k] += ydata[i]*temp;
00761                 for (j = 0; j <= k; j++) {
00762                     temp = 1.0;
00763                     if (k + j + 2*ilim != 0)
00764                         temp = pow(xdata[i],(double)(k+j+2*ilim));
00765                     a[j][k] += temp;
00766                 }
00767             }
00768         }
00769         for (k = 1; k < ncoefs; k++) 
00770             for (j = 0; j < k; j++) 
00771                 a[k][j] = a[j][k];
00772 
00773         /* Solve linear equations */
00774 
00775         retval = vircam_solve_gauss(a,b,ncoefs);
00776         if (retval != VIR_OK) {
00777             cpl_msg_warning(fctid,"Fit failed");
00778             freearray(*polycf);
00779             freespace2(a,ncoefs);
00780             freespace(b);
00781             freespace(pm);
00782             freespace(res);
00783             return(VIR_FATAL);
00784         }
00785 
00786         /* Ok, assuming this is OK, then fill the polynomial coefficients */
00787 
00788         for (i = 0; i < ncoefs; i++)
00789             pdata[i] = b[i];
00790 
00791         /* Calculate the fit quality */
00792 
00793         sum = 0.0;
00794         sumsq = 0.0;
00795         n = 0;
00796         for (i = 0; i < npts; i++) {
00797             if (pm[i] == 1)
00798                 continue;
00799             val = 0.0;
00800             for (j = ilim; j <= ncoefs; j++)
00801                 val += pdata[j-ilim]*pow(xdata[i],(double)j);
00802             res[i] = val - ydata[i];
00803             sum += res[i];
00804             sumsq += pow(res[i],2.0);
00805             n++;
00806         }
00807         sum /= (double)n;
00808         *sigfit = sqrt(sumsq/(double)n - sum*sum);
00809 
00810         /* If this is not the last iteration, then do some clipping */
00811 
00812         lcut = sum - lcut*(*sigfit);
00813         hcut = sum + hcut*(*sigfit);
00814         if (iter < niter - 1) {
00815             for (i = 0; i < npts; i++) {
00816                 if (pm[i] == 1) 
00817                     continue;
00818                 if (res[i] > hcut || res[i] < lcut) {
00819                     nnew++;
00820                     pm[i] = 1;
00821                 }
00822             }
00823         }
00824         
00825         /* If no new points have been clipped, then get out of here now... */
00826         
00827         if (nnew == 0) 
00828             break;
00829     }
00830 
00831     /* Tidy up and get out of here */
00832 
00833     freespace2(a,ncoefs);
00834     freespace(b);
00835     freespace(pm);
00836     freespace(res);
00837     return(VIR_OK);
00838 }
00839 
00840 /*---------------------------------------------------------------------------*/
00881 /*---------------------------------------------------------------------------*/
00882 
00883 extern void vircam_difference_image(cpl_image *master, cpl_image *prog, 
00884                                     unsigned char *bpm, cpl_table *chantab, 
00885                                     int ncells, int oper, float *global_diff, 
00886                                     float *global_rms, cpl_image **diffim, 
00887                                     cpl_table **diffimstats) {
00888     float *ddata,*work,mean,sig,med,mad;
00889     long nx,ny,npts;
00890     int nrows,i,nc1,nc2,nr,ixmin,ixmax,iymin,iymax,cnum,cx,cy,idx,idy;
00891     int icx,icy,indy1,indy2,indx1,indx2,jp,jcx,jj,jcy,ii,ncx,ncy;
00892     const char *fctid = "vircam_difference_image";
00893 
00894     /* Initialise the output */
00895 
00896     *diffim = NULL;
00897     *diffimstats = NULL;
00898     *global_diff = 0.0;
00899     *global_rms = 0.0;
00900 
00901     /* Is there a programme frame or a master? */
00902 
00903     if (prog == NULL || master == NULL)
00904         return;
00905 
00906     /* Start by subtracting the master image from the programme image */
00907 
00908     switch (oper) {
00909     case 1:
00910         *diffim = cpl_image_subtract_create(prog,master);
00911         break;
00912     case 2:
00913         *diffim = cpl_image_divide_create(prog,master);
00914         break;
00915     default:
00916         *diffim = NULL;
00917         cpl_msg_error(fctid,"Invalid operation requested %d",oper);
00918         break;
00919     }      
00920     if (*diffim == NULL)
00921         return;
00922 
00923     /* Work out a median difference */
00924 
00925     ddata = cpl_image_get_data_float(*diffim);
00926     nx = cpl_image_get_size_x(*diffim);
00927     ny = cpl_image_get_size_y(*diffim);
00928     npts = nx*ny;
00929     vircam_medmad(ddata,bpm,npts,global_diff,global_rms);
00930 
00931     /* Is there a channel table? */
00932 
00933     if (chantab == NULL)
00934         return;
00935 
00936     /* Before going any further, check that the channel table has all of 
00937        the columns we need */
00938 
00939     if (! cpl_table_has_column(chantab,"ixmin") ||
00940         ! cpl_table_has_column(chantab,"ixmax") ||
00941         ! cpl_table_has_column(chantab,"iymin") ||
00942         ! cpl_table_has_column(chantab,"iymax") ||
00943         ! cpl_table_has_column(chantab,"channum")) {
00944             cpl_msg_error(fctid,"Channel table is missing one of the required columns");
00945 
00946             return;
00947     }
00948 
00949     /* Work out how to divide the channels */
00950 
00951     switch (ncells) {
00952     case 1:
00953         nc1 = 1;
00954         nc2 = 1;
00955         break;
00956     case 2:
00957         nc1 = 2;
00958         nc2 = 1;
00959         break;
00960     case 4:
00961         nc1 = 4;
00962         nc2 = 1;
00963         break;
00964     case 8:
00965         nc1 = 8;
00966         nc2 = 1;
00967         break;
00968     case 16:
00969         nc1 = 16;
00970         nc2 = 1;
00971         break;
00972     case 32:
00973         nc1 = 16;
00974         nc2 = 2;
00975         break;
00976     case 64:
00977         nc1 = 32;
00978         nc2 = 2;
00979         break;
00980     default:
00981         nc1 = 32;
00982         nc2 = 2;
00983         break;
00984     }
00985 
00986     /* Create a difference image stats table */
00987 
00988     nrows = cpl_table_count_selected(chantab);
00989     *diffimstats = vircam_create_diffimg_stats(nrows*nc1*nc2);
00990 
00991     /* Loop for each data channel now */
00992 
00993     nr = 0;
00994     for (i = 0; i < nrows; i++) {
00995         ixmin = cpl_table_get_int(chantab,"ixmin",i,NULL);
00996         ixmax = cpl_table_get_int(chantab,"ixmax",i,NULL);
00997         iymin = cpl_table_get_int(chantab,"iymin",i,NULL);
00998         iymax = cpl_table_get_int(chantab,"iymax",i,NULL);
00999         cnum = cpl_table_get_int(chantab,"channum",i,NULL);
01000 
01001         /* Which is the long axis? If the channels are rectangular then
01002            divide the long axis by the greater number of cells. If the number
01003            of cells is less than 8, then just divide the long axis. If the
01004            channel is square, then divide the X axis more finely */
01005 
01006         cx = ixmax - ixmin + 1;
01007         cy = iymax - iymin + 1;
01008         if (cx > cy) {
01009             ncx = max(nc1,nc2);
01010             ncy = min(nc1,nc2);
01011         } else if (cx < cy) {
01012             ncy = max(nc1,nc2);
01013             ncx = min(nc1,nc2);
01014         } else {
01015             ncx = max(nc1,nc2);
01016             ncy = min(nc1,nc2);
01017         }
01018         
01019         /* How big is the cell? */
01020 
01021         idy = cy/ncy;
01022         idx = cx/ncx;
01023         work = cpl_malloc(idx*idy*sizeof(*work));
01024 
01025         /* Now loop for aach cell */
01026 
01027         for (icy = 0; icy < ncy; icy++) {
01028             indy1 = idy*icy;
01029             indy2 = min(iymax,indy1+idy-1);
01030             for (icx = 0; icx < ncx; icx++) {
01031                 indx1 = idx*icx;
01032                 indx2 = min(ixmax,indx1+idx-1);
01033                 jp = 0;
01034                 for (jcy = indy1; jcy < indy2; jcy++) {
01035                     jj = jcy*nx;
01036                     for (jcx = indx1; jcx < indx2; jcx++) {
01037                         ii = jj + jcx;
01038                         if (bpm != NULL && bpm[ii] == 0)
01039                             work[jp++] = ddata[ii];
01040                     }
01041                 }
01042                 (void)vircam_meansig(work,NULL,(long)jp,&mean,&sig);
01043                 (void)vircam_medmad(work,NULL,(long)jp,&med,&mad);
01044                 cpl_table_set_int(*diffimstats,"xmin",nr,indx1+1);
01045                 cpl_table_set_int(*diffimstats,"xmax",nr,indx2+1);
01046                 cpl_table_set_int(*diffimstats,"ymin",nr,indy1+1);
01047                 cpl_table_set_int(*diffimstats,"ymax",nr,indy2+1);
01048                 cpl_table_set_int(*diffimstats,"chan",nr,cnum);
01049                 cpl_table_set_float(*diffimstats,"mean",nr,mean);
01050                 cpl_table_set_float(*diffimstats,"median",nr,med);
01051                 cpl_table_set_float(*diffimstats,"variance",nr,(sig*sig));
01052                 cpl_table_set_float(*diffimstats,"mad",nr++,mad);
01053             }
01054         }
01055         cpl_free(work);
01056     }
01057 }               
01058 
01059 /*---------------------------------------------------------------------------*/
01076 /*---------------------------------------------------------------------------*/
01077 
01078 cpl_table *vircam_create_diffimg_stats(int nrows) {
01079     cpl_table *diffimstats;
01080 
01081     diffimstats = cpl_table_new(nrows);
01082     cpl_table_new_column(diffimstats,"xmin",CPL_TYPE_INT);
01083     cpl_table_set_column_unit(diffimstats,"xmin","pixels");
01084     cpl_table_new_column(diffimstats,"xmax",CPL_TYPE_INT);
01085     cpl_table_set_column_unit(diffimstats,"xmax","pixels");
01086     cpl_table_new_column(diffimstats,"ymin",CPL_TYPE_INT);
01087     cpl_table_set_column_unit(diffimstats,"ymin","pixels");
01088     cpl_table_new_column(diffimstats,"ymax",CPL_TYPE_INT);
01089     cpl_table_set_column_unit(diffimstats,"ymax","pixels");
01090     cpl_table_new_column(diffimstats,"chan",CPL_TYPE_INT);
01091     cpl_table_set_column_unit(diffimstats,"chan","pixels");
01092     cpl_table_new_column(diffimstats,"mean",CPL_TYPE_FLOAT);
01093     cpl_table_set_column_unit(diffimstats,"mean","ADU");
01094     cpl_table_new_column(diffimstats,"median",CPL_TYPE_FLOAT);
01095     cpl_table_set_column_unit(diffimstats,"median","ADU");
01096     cpl_table_new_column(diffimstats,"variance",CPL_TYPE_FLOAT);
01097     cpl_table_set_column_unit(diffimstats,"variance","ADU**2");
01098     cpl_table_new_column(diffimstats,"mad",CPL_TYPE_FLOAT);
01099     cpl_table_set_column_unit(diffimstats,"mad","ADU");
01100     return(diffimstats);
01101 }
01102 
01103 
01104 /*---------------------------------------------------------------------------*/
01128 /*---------------------------------------------------------------------------*/
01129 
01130 extern void vircam_sort(float **a, int n, int m) {
01131     int iii,ii,i,ifin,j,it;
01132     float *t;
01133 
01134     t = cpl_malloc(m*sizeof(*t));
01135 
01136     iii = 2;
01137     while (iii < n)
01138         iii *= 2;
01139     iii = min(n,(3*iii)/4 - 1);
01140 
01141     while (iii > 1) {
01142         iii /= 2;
01143         ifin = n - iii;
01144         for (ii = 0; ii < ifin; ii++) {
01145             i = ii;
01146             j = i + iii;
01147             if (a[0][i] > a[0][j]) {
01148                 for (it = 0; it < m; it++)
01149                    t[it] = a[it][j];
01150                 while (1) {
01151                     for (it = 0; it < m; it++) 
01152                         a[it][j] = a[it][i];
01153                     j = i;
01154                     i = i - iii;
01155                     if (i < 0 || a[0][i] <= t[0]) 
01156                         break;
01157                 }
01158                 for (it = 0; it < m; it++) 
01159                     a[it][j] = t[it];
01160             }
01161         }
01162     }
01163     cpl_free(t);
01164 }
01165 
01166 
01167 /*---------------------------------------------------------------------------*/
01184 /*---------------------------------------------------------------------------*/
01185 
01186 extern long vircam_getnpts(cpl_image *in) {
01187     int nx,ny;
01188     long npts;
01189     const char *fctid = "vircam_getnpts";
01190 
01191     if ((nx = cpl_image_get_size_x(in)) == -1) {
01192         cpl_msg_error(fctid,"NULL image input");
01193         return(0);
01194     }
01195     if ((ny = cpl_image_get_size_y(in)) == -1) {
01196         cpl_msg_error(fctid,"NULL image input");
01197         return(0);
01198     }
01199     npts = (long)nx*ny;
01200     return(npts);
01201 }
01202 
01203 /*---------------------------------------------------------------------------*/
01234 /*---------------------------------------------------------------------------*/
01235 
01236 extern int vircam_fndmatch(float x, float y, float *xlist, float *ylist, 
01237                            int nlist, float err) {
01238     int isp,ifp,index,i;
01239     float errsq,errmin,dx,dy,poserr;
01240 
01241     /* Find lower limit index */
01242 
01243     isp = 0;
01244     ifp = nlist - 1;
01245     errsq = err*err;
01246     index = (isp + ifp)/2;
01247     while (ifp-isp >= 2) {
01248         if (ylist[index] < y - err) {
01249             isp = index;
01250             index = (index+ifp)/2;
01251         } else if (ylist[index] > y - err) {
01252             ifp = index;
01253             index = (index+isp)/2;
01254         } else {
01255             isp = index;
01256             break;
01257         }
01258     }
01259 
01260     /* Now find nearest one within limit */
01261 
01262     index = -1;
01263     errmin = errsq;
01264     for (i = isp; i < nlist; i++) {
01265         if (ylist[i] > y+err)
01266             break;
01267         dx = x - xlist[i];
01268         dy = y - ylist[i];
01269         poserr = dx*dx + dy*dy;
01270         if (poserr < errsq) {
01271             if (poserr <= errmin) {
01272                 index = i;
01273                 errmin = poserr;
01274             }
01275         }
01276     }
01277     return(index);
01278 }
01279 
01280 /*---------------------------------------------------------------------------*/
01299 /*---------------------------------------------------------------------------*/
01300 
01301 extern int *vircam_dummy_confidence(long n) {
01302     int *cdata,i;
01303  
01304     cdata = cpl_malloc(n*sizeof(*cdata));
01305     for (i = 0; i < n; i++)
01306         cdata[i] = 100;
01307     return(cdata);
01308 }
01309 
01310 /*---------------------------------------------------------------------------*/
01333 /*---------------------------------------------------------------------------*/
01334 
01335 extern int vircam_compare_dims(cpl_image *im1, cpl_image *im2) {
01336     
01337     if ((cpl_image_get_size_x(im1) != cpl_image_get_size_x(im2)) ||
01338         cpl_image_get_size_y(im1) != cpl_image_get_size_y(im2))
01339         return(VIR_FATAL);
01340     else
01341         return(VIR_OK);
01342 }
01343    
01344 /*---------------------------------------------------------------------------*/
01365 /*---------------------------------------------------------------------------*/
01366 
01367 extern void vircam_prov(cpl_propertylist *p, vir_fits **inlist, int n) {
01368     int i;
01369     char keyword[SZKEY],value[SZVAL],*fn,*base;
01370 
01371     /* Delete all the provenance keywords that might already exist */
01372 
01373     cpl_propertylist_erase_regexp(p,"ESO DRS PROV*",0);
01374 
01375     /* Add the new provenance keywords */
01376 
01377     for (i = 0; i < n; i++) {
01378         (void)snprintf(keyword,SZKEY,"ESO DRS PROV%04d",i+1);
01379         fn = cpl_strdup(vircam_fits_get_fullname(inlist[i]));
01380         base = basename(fn);
01381         (void)snprintf(value,SZVAL,"%s",base);
01382         cpl_free(fn);
01383         cpl_propertylist_update_string(p,keyword,value);
01384         (void)snprintf(value,SZVAL,"Input file # %d",i+1);
01385         cpl_propertylist_set_comment(p,keyword,value);
01386     }
01387 }
01388 
01389 /*---------------------------------------------------------------------------*/
01407 /*---------------------------------------------------------------------------*/
01408 
01409 extern void vircam_merge_propertylists(cpl_propertylist *p1, 
01410                                        cpl_propertylist *p2) {
01411     int i;
01412     
01413     /* Return if either propertylist is NULL */
01414 
01415     if (p1 == NULL || p2 == NULL)
01416         return;
01417 
01418     /* Copy each property from the second list into the first one */
01419 
01420     for (i = 0; i < cpl_propertylist_get_size(p2); i++)
01421         cpl_propertylist_copy_property(p1,p2,cpl_property_get_name(cpl_propertylist_get(p2,i)));
01422 }
01423 
01424 
01425 /*---------------------------------------------------------------------------*/
01444 /*---------------------------------------------------------------------------*/
01445 
01446 extern void vircam_dummy_property(cpl_propertylist *p) {
01447 
01448     /* Check for silly input */
01449 
01450     if (p == NULL)
01451         return;
01452 
01453     /* Add the property now */
01454 
01455     cpl_propertylist_update_bool(p,"ESO DRS IMADUMMY",TRUE);
01456     cpl_propertylist_set_comment(p,"ESO DRS IMADUMMY",
01457                                  "This is a dummy product");
01458     return;
01459 }
01460 
01461 /*---------------------------------------------------------------------------*/
01477 /*---------------------------------------------------------------------------*/
01478 
01479 extern int vircam_is_dummy(cpl_propertylist *p) {
01480 
01481     /* Check for silly input */
01482 
01483     if (p == NULL)
01484         return(0);
01485 
01486     /* Check the propertylist and return the result */
01487 
01488     return(cpl_propertylist_has(p,"ESO DRS IMADUMMY"));
01489 }
01490 
01491 /*---------------------------------------------------------------------------*/
01521 /*---------------------------------------------------------------------------*/
01522 
01523 extern void vircam_overexp(vir_fits **fitslist, int *n, float lthr, 
01524                            float hthr, int ditch) {
01525     int i,m;
01526     cpl_image *im;
01527     double val;
01528 
01529     /* Loop for each of the fits items */
01530 
01531     m = 0;
01532     for (i = 0; i < *n; i++) {
01533         im = vircam_fits_get_image(fitslist[i]);
01534         val = cpl_image_get_mean_window(im,100,100,200,200);
01535         if (val > lthr && val < hthr)
01536             fitslist[m++] = fitslist[i];
01537         else
01538             if (ditch) 
01539                 vircam_fits_delete(fitslist[i]);
01540     }
01541     for (i = m; i < *n; i++)
01542         fitslist[i] = NULL;
01543     *n = m;
01544     
01545 }
01546 
01547 /*---------------------------------------------------------------------------*/
01564 /*---------------------------------------------------------------------------*/
01565 
01566 extern cpl_image *vircam_dummy_image(vir_fits *model) {
01567     cpl_image *im;
01568 
01569     /* Copy the input model image */
01570 
01571     im = cpl_image_duplicate(vircam_fits_get_image(model));
01572 
01573     /* Now change it all to zeros */
01574 
01575     cpl_image_multiply_scalar(im,0.0);
01576 
01577     /* Return the result */
01578 
01579     return(im);
01580 }
01581 
01582 /*---------------------------------------------------------------------------*/
01599 /*---------------------------------------------------------------------------*/
01600 
01601 extern cpl_table *vircam_dummy_catalogue(int type) {
01602 
01603     cattype = type;
01604     tabinit(NULL);
01605     return(tab);
01606 }
01607 
01608 
01609 /*---------------------------------------------------------------------------*/
01626 /*---------------------------------------------------------------------------*/
01627 
01628 extern cpl_table *vircam_illcor_newtab(int nrows) {
01629     cpl_table *illcor;
01630     int i;
01631 
01632     /* Create the table structure and fill in the columns */
01633 
01634     illcor = cpl_table_new(nrows);
01635     for (i = 0; i < NI_COLS; i++) 
01636         cpl_table_new_column(illcor,illcor_cols[i],CPL_TYPE_FLOAT);
01637     return(illcor);
01638 }
01639 
01640 /*---------------------------------------------------------------------------*/
01664 /*---------------------------------------------------------------------------*/
01665 
01666 extern void vircam_timestamp(char *out, int n) {
01667     struct timeval tv;
01668     struct tm *tm;
01669     float sec;
01670     
01671     /* Get the Greenwich Mean Time */
01672 
01673     (void)gettimeofday(&tv,NULL);
01674     tm = gmtime(&(tv.tv_sec));
01675     sec = (float)tm->tm_sec + 1.0e-6*(float)tv.tv_usec;
01676 
01677     /* Now format it */
01678 
01679     (void)snprintf(out,n,"%04d-%02d-%02dT%02d:%02d:%07.4f",1900+tm->tm_year,
01680                    tm->tm_mon+1,tm->tm_mday,tm->tm_hour,tm->tm_min,sec);
01681 }
01682 
01683 /*---------------------------------------------------------------------------*/
01712 /*---------------------------------------------------------------------------*/
01713 
01714 extern void vircam_backmap(vir_fits *in, vir_mask *mask, int nbsize, 
01715                            cpl_image **out, float *med) {
01716     int i,j,nx,ny,ifracx,ifracy,nbsizx,nbsizy,nbx,nby,*nps,jx,jy;
01717     int jy1,jy2,np,nout,jyp1,jxp1,jz1,jz2,jz3,jz4,nbsize2;
01718     float fracx,fracy,*bmap,**rowpoints,*data,*ptr,dely,delx,t1,t2;
01719     unsigned char *bpm;
01720 
01721     /* Check to see if nbsize is close to exact divisor */
01722 
01723     nx = cpl_image_get_size_x(vircam_fits_get_image(in));
01724     ny = cpl_image_get_size_y(vircam_fits_get_image(in));
01725     fracx = ((float)nx)/((float)nbsize);
01726     fracy = ((float)ny)/((float)nbsize);
01727     ifracx = (int)(fracx + 0.1);
01728     ifracy = (int)(fracy + 0.1);
01729     nbsizx = nx/ifracx;
01730     nbsizy = ny/ifracy;
01731     nbsize = max(vircam_nint(0.9*nbsize),min(nbsize,min(nbsizx,nbsizy)));
01732     nbsize = min(nx,min(ny,nbsize)); /* trap for small maps */
01733     nbsize2 = nbsize/2;
01734 
01735     /* Divide the map into partitions */
01736 
01737     nbx = nx/nbsize;
01738     nby = ny/nbsize;
01739 
01740     /* Get some space for to use for accumulating stats */
01741 
01742     bmap = cpl_malloc(nbx*nby*sizeof(float));
01743     rowpoints = cpl_malloc(nbx*sizeof(float *));
01744     nps = cpl_malloc(nbx*sizeof(int));
01745 
01746     /* Get the data from the input file and the mask */
01747 
01748     data = cpl_image_get_data_float(vircam_fits_get_image(in));
01749     bpm = vircam_mask_get_data(mask);
01750 
01751     /* Work out the global median */
01752 
01753     *med = vircam_med(data,bpm,(long)(nx*ny));
01754 
01755     /* For each cell allocate some nbsize*nbsize pixels in the row pointers.
01756        It's ok to do this here since this will be the maximum that each cell
01757        will have... */
01758 
01759     for (i = 0; i < nbx; i++) 
01760         rowpoints[i] = cpl_malloc(nbsize*nbsize*sizeof(float));
01761     
01762     /* Loop for each row of cells and work out which rows in the input
01763        map these relate to */
01764 
01765     for (jy = 0; jy < nby; jy++) {
01766         jy1 = jy*nbsize;
01767         jy2 = min((jy+1)*nbsize - 1,ny-1);
01768 
01769         /* Zero the counter for each cell */
01770 
01771         for (jx = 0; jx < nbx; jx++)
01772             nps[jx] = 0;
01773 
01774         /* Loop for the strip of the input map covered by this row of cells.
01775            Work out which cell the current pixel is in and then put it into
01776            the correct accumulator assuming the bad pixel mask says it's ok */
01777 
01778         for (j = jy1; j <= jy2; j++) {
01779             for (i = 0; i < nx; i++) {
01780                 jx = min(nbx-1,i/nbsize);
01781                 if (bpm[j*nx + i] == 0) {
01782                     ptr = rowpoints[jx];
01783                     np = nps[jx];
01784                     ptr[np++] = data[j*nx + i];
01785                     nps[jx] = np;
01786                 }
01787             }
01788         }
01789 
01790         /* Now do the medians for each of the cells in this row */
01791 
01792         for (jx = 0; jx < nbx; jx++) 
01793             bmap[jy*nbx+jx] = vircam_med(rowpoints[jx],NULL,(long)(nps[jx])) - *med;
01794     }
01795 
01796     /* Free up some workspace */
01797 
01798     for (jx = 0; jx < nbx; jx++) 
01799         freespace(rowpoints[jx]);
01800     freespace(rowpoints);
01801     freespace(nps);
01802 
01803     /* Create a new image for the output array */
01804 
01805     *out = cpl_image_new(nx,ny,CPL_TYPE_FLOAT);
01806     ptr = cpl_image_get_data_float(*out);
01807 
01808     /* Ok, we now have the block averaged background map. Do a two point
01809        interpolation to create the full map now */
01810 
01811     nout = 0;
01812     for (j = 1; j <= ny; j++) {
01813         jy = (j + nbsize2)/nbsize;
01814         jyp1 = jy + 1;
01815         jy = min(nby,max(1,jy));
01816         jyp1 = min(nby,jyp1);
01817         dely = (float)(j - nbsize*jy + nbsize2)/(float)nbsize;
01818         dely = max(0.0,min(1.0,dely));
01819         for (i = 1; i <= nx; i++) {
01820             jx = (i + nbsize2)/nbsize;
01821             jxp1 = jx + 1;
01822             jx = min(nbx,max(1,jx));
01823             jxp1 = min(nbx,jxp1);
01824             delx = (float)(i - nbsize*jx + nbsize2)/(float)nbsize;
01825             delx = max(0.0,min(1.0,delx));
01826             jz1 = (jy - 1)*nbx + jx;
01827             jz2 = (jyp1 - 1)*nbx + jx;
01828             if (jx == nbx) {
01829                 jz3 = jz1;
01830                 jz4 = jz2;
01831             } else {
01832                 jz3 = jz1 + 1;
01833                 jz4 = jz2 + 1;
01834             }
01835             t1 = (1.0 - delx)*bmap[jz1-1] + delx*bmap[jz3-1];
01836             t2 = (1.0 - delx)*bmap[jz2-1] + delx*bmap[jz4-1];
01837             ptr[nout++] = (1.0 - dely)*t1 + dely*t2;
01838         }
01839     }
01840 
01841     /* Tidy up before leaving... */
01842 
01843     freespace(bmap)
01844 }
01845 
01846 /*---------------------------------------------------------------------------*/
01867 /*---------------------------------------------------------------------------*/
01868 
01869 extern int vircam_findcol(cpl_propertylist *p, char *col) {
01870 
01871     if (!strcmp(col,"X")) {
01872         if (cpl_propertylist_has(p,"ESO DRS XCOL"))
01873             return(cpl_propertylist_get_int(p,"ESO DRS XCOL"));
01874         else 
01875             return(-1);
01876     } else if (!strcmp(col,"Y")) {
01877         if (cpl_propertylist_has(p,"ESO DRS YCOL"))
01878             return(cpl_propertylist_get_int(p,"ESO DRS YCOL"));
01879         else 
01880             return(-1);
01881     }
01882     return(-1);
01883 }
01884 
01885 /*---------------------------------------------------------------------------*/
01906 /*---------------------------------------------------------------------------*/
01907 
01908 extern void vircam_rename_property(cpl_propertylist *p, char *oldname,
01909                                    char *newname) {
01910     cpl_propertylist *temp;
01911     cpl_property *property;
01912 
01913     /* First get the old property. Note that we have to do this in a 
01914        particularly silly way since you cannot reference an individual property
01915        in a propertylist by its name. You can only do it by its index 
01916        number. Remeber to change this when CPL comes to it's senses... */
01917 
01918     if (! cpl_propertylist_has(p,oldname))
01919         return;
01920     temp = cpl_propertylist_new();
01921     cpl_propertylist_copy_property(temp,p,oldname);
01922     property = cpl_propertylist_get(temp,0);
01923 
01924     /* Now change its name */
01925 
01926     cpl_property_set_name(property,newname);
01927 
01928     /* Now insert this into the propertylist and delete the old one */
01929 
01930     cpl_propertylist_append(p,temp);
01931     cpl_propertylist_erase(p,oldname);
01932     cpl_propertylist_delete(temp);
01933 }
01934 
01935 /*---------------------------------------------------------------------------*/
01963 /*---------------------------------------------------------------------------*/
01964 
01965 extern int vircam_catpars(cpl_frame *index, char **catpath, char **catname) {
01966     cpl_propertylist *p;
01967     char *unk = "unknown",*fname;
01968     int status;
01969     const char *fctid = "vircam_catpars";
01970 
01971     /* Initialise a few things */
01972 
01973     *catpath = NULL;
01974     *catname = NULL;
01975 
01976     /* First get the full path to the index file and make sure it exists */
01977 
01978     fname = cpl_strdup(cpl_frame_get_filename(index));
01979     if (access((const char *)fname,R_OK) != 0) {
01980         cpl_msg_error(fctid,"Can't access index file %s",fname);
01981         cpl_free(fname);
01982         return(VIR_FATAL);
01983     }
01984     *catpath = cpl_strdup(dirname(fname));
01985 
01986     /* Load the propertylist if you can. If you can't then signal a fatal
01987        error since this probably means the whole file is messed up */
01988 
01989     if ((p = cpl_propertylist_load(cpl_frame_get_filename(index),0)) == NULL) {
01990         freespace(*catpath);
01991         cpl_msg_error(fctid,"Can't load index file header %s",fname);
01992         cpl_free(fname);
01993         return(VIR_FATAL);
01994     }
01995         
01996     /* If there is a catalogue name in the header then send it back. If there
01997        isn't then give a default name and send a warning */
01998 
01999     if (cpl_propertylist_has(p,"CATNAME")) {
02000         *catname = cpl_strdup(cpl_propertylist_get_string(p,"CATNAME"));
02001         status = VIR_OK;
02002     } else {
02003         *catname = cpl_strdup(unk);
02004         cpl_msg_warn(fctid,"Property CATNAME not in index file header %s",
02005                      fname);
02006         status = VIR_WARN;
02007     }
02008     cpl_free(fname);
02009     freepropertylist(p);
02010 
02011     /* Get out of here */
02012 
02013     return(status);
02014 }
02015 
02016 /*---------------------------------------------------------------------------*/
02045 /*---------------------------------------------------------------------------*/
02046 
02047 extern int vircam_gaincor_calc(cpl_frame *frame, int *n, float **cors,
02048                                int *status) {
02049     float sum,val;
02050     int i,ngood;
02051     unsigned char *iflag;
02052     cpl_propertylist *p;
02053 
02054     /* Inherited status */
02055 
02056     if (*status != VIR_OK)
02057         return(*status);
02058 
02059     /* Find the number of extensions in the file and allocate some workspace
02060        to hold the results. Allocate a small workspace to flag dummy 
02061        extensions */
02062 
02063     *n = cpl_frame_get_nextensions(frame);
02064     *cors = cpl_malloc(*n*sizeof(float));
02065     iflag = cpl_calloc(*n,sizeof(iflag));
02066 
02067     /* Ok, loop through the extensions and read the propertylists */
02068 
02069     sum = 0.0;
02070     ngood = 0;
02071     for (i = 0; i < *n; i++) {
02072         p = cpl_propertylist_load(cpl_frame_get_filename(frame),i+1);
02073         if (cpl_propertylist_has(p,"ESO DRS IMADUMMY")) {
02074             iflag[i] = 1;
02075         } else if (! cpl_propertylist_has(p,"ESO DRS MEDFLAT")) {
02076             iflag[i] = 1;
02077         } else {
02078             val = cpl_propertylist_get_double(p,"ESO DRS MEDFLAT");
02079             if (val == 0.0) {
02080                 iflag[i] = 1;
02081             } else {
02082                 sum += val;
02083                 (*cors)[i] = val;
02084                 ngood++;
02085             }
02086         }
02087         cpl_propertylist_delete(p);
02088     }
02089 
02090     /* If any of them are good, then work out what the average is and
02091        create the correction factors. If the image was a dummy or
02092        there was no MEDFLAT keyword in the header, then just give it
02093        a factor of 1 */
02094 
02095     if (ngood > 0)
02096         sum /= (float)ngood;
02097     for (i = 0; i < *n; i++) {
02098         if (iflag[i] == 0) {
02099             (*cors)[i] = sum/(*cors)[i];
02100         } else {
02101             (*cors)[i] = 1.0;
02102         }
02103     }
02104     cpl_free(iflag);
02105 
02106     /* Get out of here */
02107 
02108     GOOD_STATUS
02109 }
02110 
02113 /*
02114 
02115 $Log: vircam_utils.c,v $
02116 Revision 1.56  2008/07/10 13:04:03  jim
02117 Fixed some comments
02118 
02119 Revision 1.55  2008/05/06 12:15:02  jim
02120 Changed vircam_catpars to check that we can read the index file and to send
02121 out error messages if there are problems
02122 
02123 Revision 1.54  2008/05/06 08:38:29  jim
02124 Modified call to cpl_propertylist_get_ from float to double in
02125 vircam_gaincor_calc.
02126 
02127 Revision 1.53  2007/11/14 14:47:53  jim
02128 vircam_linfit now works only with doubles
02129 
02130 Revision 1.52  2007/11/14 10:46:07  jim
02131 Fixed bugs in vircam_polyfit
02132 
02133 Revision 1.51  2007/10/25 17:33:31  jim
02134 Modified to add vircam_polyfit and to remove lint warnings
02135 
02136 Revision 1.50  2007/10/19 09:25:10  jim
02137 Fixed problems with missing includes
02138 
02139 Revision 1.49  2007/10/19 06:55:06  jim
02140 Modifications made to use new method for directing the recipes to the
02141 standard catalogues using the sof
02142 
02143 Revision 1.48  2007/10/15 12:50:28  jim
02144 Modified for compatibility with cpl_4.0
02145 
02146 Revision 1.47  2007/07/09 13:21:06  jim
02147 Changed argument list to vircam_exten_range to include a cpl_frame
02148 
02149 Revision 1.46  2007/05/08 10:41:01  jim
02150 Added vircam_gaincor_calc
02151 
02152 Revision 1.45  2007/05/03 11:15:33  jim
02153 Fixed little problem with table wcs
02154 
02155 Revision 1.44  2007/05/02 09:14:48  jim
02156 Added vircam_findcol and vircam_rename_property
02157 
02158 Revision 1.43  2007/03/01 12:42:42  jim
02159 Modified slightly after code checking
02160 
02161 Revision 1.42  2007/02/14 12:53:34  jim
02162 Removed vircam_paf_print_header
02163 
02164 Revision 1.41  2007/02/07 10:12:15  jim
02165 Removed vircam_ndit_correct as it is now no longer necessary
02166 
02167 Revision 1.40  2007/01/15 12:36:27  jim
02168 Fixed bug in nditcor where factor was being divided by instead of multiplied
02169 
02170 Revision 1.39  2006/12/06 12:58:41  jim
02171 Fixed a small bug in vircam_backmap
02172 
02173 Revision 1.38  2006/12/01 14:08:33  jim
02174 added vircam_backmap
02175 
02176 Revision 1.37  2006/11/28 20:52:17  jim
02177 Added vircam_illcor_newtab
02178 
02179 Revision 1.36  2006/11/27 12:00:17  jim
02180 cosmetic changes
02181 
02182 Revision 1.35  2006/10/05 09:23:00  jim
02183 Small modifications to a couple of cpl calls to bring them into line with
02184 cpl v3.0
02185 
02186 Revision 1.34  2006/08/27 20:28:15  jim
02187 added vircam_is_dummy
02188 
02189 Revision 1.33  2006/08/21 09:09:29  jim
02190 Added routines vircam_create_diffimg_stats and vircam_dummy_property
02191 
02192 Revision 1.32  2006/07/11 14:53:58  jim
02193 Fixed vircam_ndit_correct so that images with bad status are not corrected
02194 
02195 Revision 1.31  2006/07/04 09:19:06  jim
02196 replaced all sprintf statements with snprintf
02197 
02198 Revision 1.30  2006/06/30 21:32:09  jim
02199 Fixed bug in vircam_prov so that sprintf is replaced by snprintf
02200 
02201 Revision 1.29  2006/06/20 18:59:51  jim
02202 Fixed small bug in vircam_ndit_correct
02203 
02204 Revision 1.28  2006/06/09 22:24:47  jim
02205 Added vircam_ndit_correct
02206 
02207 Revision 1.27  2006/06/09 11:26:26  jim
02208 Small changes to keep lint happy
02209 
02210 Revision 1.26  2006/06/06 13:07:54  jim
02211 Change the stats routine used in vircam_difference_image
02212 
02213 Revision 1.25  2006/05/30 13:31:47  jim
02214 Added vircam_timestamp
02215 
02216 Revision 1.24  2006/05/24 13:35:49  jim
02217 Added vircam_dummy_image and vircam_dummy_catalogue
02218 
02219 Revision 1.23  2006/05/04 15:16:33  jim
02220 Fixed memory bug in vircam_overexp
02221 
02222 Revision 1.22  2006/04/20 11:29:41  jim
02223 Added vircam_overexp
02224 
02225 Revision 1.21  2006/03/22 11:42:32  jim
02226 Moved some routines into vircam_mask
02227 
02228 Revision 1.20  2006/03/15 10:43:42  jim
02229 Fixed a few things
02230 
02231 Revision 1.19  2006/03/08 14:32:22  jim
02232 Lots of little modifications
02233 
02234 Revision 1.18  2006/03/01 10:31:29  jim
02235 Now uses new vir_fits objects
02236 
02237 Revision 1.17  2006/02/18 11:48:55  jim
02238 *** empty log message ***
02239 
02240 Revision 1.16  2006/01/23 10:30:50  jim
02241 Mainly documentation mods
02242 
02243 Revision 1.15  2005/12/14 22:17:33  jim
02244 Updated docs
02245 
02246 Revision 1.14  2005/12/01 16:26:03  jim
02247 Added vircam_dummy_confidence
02248 
02249 Revision 1.13  2005/11/25 15:33:22  jim
02250 Some code fixes to keep splint happy
02251 
02252 Revision 1.12  2005/11/25 09:56:15  jim
02253 Tidied up some more documentation
02254 
02255 Revision 1.11  2005/11/25 09:36:23  jim
02256 Added vircam_linfit
02257 
02258 Revision 1.10  2005/11/07 13:15:16  jim
02259 Fixed lots of bugs and added some error checking
02260 
02261 Revision 1.9  2005/11/03 15:16:28  jim
02262 Lots of changes mainly to strengthen error reporting
02263 
02264 Revision 1.8  2005/11/03 13:28:50  jim
02265 All sorts of changes to tighten up error handling
02266 
02267 Revision 1.7  2005/10/14 13:21:04  jim
02268 Added some new routines
02269 
02270 Revision 1.6  2005/09/20 15:07:46  jim
02271 Fixed a few bugs and added a few things
02272 
02273 Revision 1.5  2005/09/07 12:47:22  jim
02274 renamed a malloc and free to the cpl equivallent
02275 
02276 Revision 1.4  2005/08/10 09:55:05  jim
02277 Modified vircam_madfit so that if determinant is zero, then it sends a
02278 zero slope and intercept back
02279 
02280 Revision 1.3  2005/08/10 08:42:00  jim
02281 Modified vircam_madfit so that if the initial least squares fit gets a
02282 perfect result, then it just returns with those values rather than to try
02283 and do the median fit
02284 
02285 Revision 1.2  2005/08/09 15:30:00  jim
02286 vircam_exten_range had number of chips wrong!
02287 
02288 Revision 1.1.1.1  2005/08/05 08:29:09  jim
02289 Initial import
02290 
02291 
02292 */
02293 

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