vircam_channel.c

/*

$Id: vircam_channel.c,v 1.12 2007/11/20 09:36:51 jim Exp $

*/

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <cpl.h>

#include "vircam_channel.h"
#include "vircam_utils.h"

/* Array of character values for all of the expected columns in a 
   channel table. Any other columns will be ignored */

#define NCHANTAB_COLS 15
static const char *chantab_columns[] = {"channum",
                                        "ixmin",
                                        "ixmax",
                                        "iymin",
                                        "iymax",
                                        "dcrpix1",
                                        "dcrpix2",
                                        "dcd1_1",
                                        "dcd1_2",
                                        "dcd2_1",
                                        "dcd2_2",
                                        "lin_10000",
                                        "lin_10000_err",
                                        "norder",
                                        "coeff_1"};

static const cpl_type chantab_types[] = {CPL_TYPE_INT,
                                         CPL_TYPE_INT,
                                         CPL_TYPE_INT,
                                         CPL_TYPE_INT,
                                         CPL_TYPE_INT,
                                         CPL_TYPE_INT,
                                         CPL_TYPE_INT,
                                         CPL_TYPE_INT,
                                         CPL_TYPE_INT,
                                         CPL_TYPE_INT,
                                         CPL_TYPE_INT,
                                         CPL_TYPE_DOUBLE,
                                         CPL_TYPE_DOUBLE,
                                         CPL_TYPE_INT,
                                         CPL_TYPE_DOUBLE};
#define NTESTPOSITIVE 7
static const char *testpositive[] = {"channum",
                                     "ixmin",
                                     "ixmax",
                                     "iymin",
                                     "iymax",
                                     "dcrpix1",
                                     "dcrpix2"};

#define NTEST1 4
static const char *test1[] = {"dcd1_1",
                              "dcd1_2",
                              "dcd2_1",
                              "dcd2_2"};

#define SZCOL 9

static int vircam_chan_init(int channum, int ixmin, int ixmax, int iymin,
                            int iymax, int dcrpix1, int dcrpix2, int dcd1_1,
                            int dcd1_2, int dcd2_1, int dcd2_2, parquet *p);
static void vircam_chan_close(parquet *p);
static int vircam_chan_testpos(cpl_table *intab, const char *column);
static int vircam_chan_test1(cpl_table *intab, const char *column);
static int vircam_chan_fill_lin(cpl_table *intab, int row, parquet *p);

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern int vircam_chantab_verify(cpl_table *intab) {
    int i,missing[NCHANTAB_COLS],wrongtype[NCHANTAB_COLS],ierr,norder;
    int nrows,null;
    double *c1;
    const char *fctid = "vircam_chantab_verify";
    char errmsg[1024],colname[SZCOL],msg[1024];
    const char *col;

    /* First check that you've even supplied a channel table...*/

    errmsg[0] = '\0';
    if (intab == NULL) {
        cpl_msg_error(fctid,"Null input channel table");
        return(VIR_FATAL);
    }

    /* Loop through all the expected channel table columns and see if they
       are there... */

    ierr = 0;
    for (i = 0; i < NCHANTAB_COLS; i++) {
        missing[i] = 0;
        wrongtype[i] = 0;
        if (cpl_table_has_column(intab,chantab_columns[i]) != 1) {
            missing[i] = 1;
            ierr = 1;
        } else if (cpl_table_get_column_type(intab,chantab_columns[i]) !=
                   chantab_types[i]) {
            wrongtype[i] = 1;
            ierr = 1;
        }
    }
    if (ierr == 1) {
        (void)strcat(errmsg,"These columns are missing or have incorrect types");
        for (i = 0; i < NCHANTAB_COLS; i++) {
            if (missing[i] == 1) {
                (void)strcat(errmsg,"\n");
                (void)strcat(errmsg,chantab_columns[i]);
                (void)strcat(errmsg," Missing");
            } else if (wrongtype[i] == 1) {
                (void)strcat(errmsg,"\n");
                (void)strcat(errmsg,chantab_columns[i]);
                (void)strcat(errmsg," Wrong type");
            }
        }
        (void)strcat(errmsg,"\n");
    }

    /* Do a separate check for any coefficient columns to make sure they
       are all doubles. This might lead to coeff_1 with 2 error messages,
       but who cares...*/

    col = cpl_table_get_column_name(intab);
    while (col != NULL) {
        if (! strncmp(col,"coeff_",6)) {
            if (cpl_table_get_column_type(intab,col) != CPL_TYPE_DOUBLE) {
                (void)sprintf(msg,"Column %s must be double\n",col);
                (void)strcat(errmsg,msg);
            }
        }
        col = cpl_table_get_column_name(NULL);
    }

    /* Now do a few basic sanity tests. First check that some of the columns
       have only positive integers */

    ierr = 0;
    for (i = 0; i < NTESTPOSITIVE; i++) {
        missing[i] = 0;
        if (vircam_chan_testpos(intab,testpositive[i]) == VIR_OK)
            continue;
        ierr = 1;
        missing[i] = 1;
    }
    if (ierr == 1) {
        (void)strcat(errmsg,"These columns should be positive integers only:");
        for (i = 0; i < NTESTPOSITIVE; i++) {
            if (missing[i] == 1) {
                (void)strcat(errmsg,"\n");
                (void)strcat(errmsg,testpositive[i]);
            }
        }
        (void)strcat(errmsg,"\n");
    }

    /* Now check that some of the columns have only -1,0,1 as values */

    ierr = 0;
    for (i = 0; i < NTEST1; i++) {
        missing[i] = 0;
        if (vircam_chan_test1(intab,test1[i]) == VIR_OK)
            continue;
        ierr = 1;
        missing[i] = 1;
    }
    if (ierr == 1) {
        (void)strcat(errmsg,"These columns should be integers and (-1,0,1) only:");
        for (i = 0; i < NTEST1; i++) {
            if (missing[i] == 1) {
                (void)strcat(errmsg,"\n");
                (void)strcat(errmsg,test1[i]);
            }
        }
        (void)strcat(errmsg,"\n");
    }

    /* Finally check that the order for the fit is matched by the number
       of coefficients in the table. Also check that the first order 
       coefficient is always 1 */

    norder = cpl_table_get_int(intab,"norder",0,&null);
    nrows = cpl_table_get_nrow(intab);
    c1 = cpl_table_get_data_double(intab,"coeff_1");
    for (i = 0; i < nrows; i++) {
        if (c1[i] != 1.0) {
            (void)strcat(errmsg,"\ncoeff_1 must be all 1s!\n");
            break;
        }
    }
    for (i = 2; i <= norder; i++) {
        missing[i] = 0;
        (void)snprintf(colname,SZCOL,"coeff_%d",i);
        if (cpl_table_has_column(intab,colname) != 1) {
            missing[i] = 1;
            ierr = 1;
        }
    }
    if (ierr == 1) {
        (void)strcat(errmsg,
               "\nThese coefficient columns are missing from the channel table:");
        for (i = 2; i < norder; i++) {
            if (missing[i] == 1) {
                (void)snprintf(colname,SZCOL,"coeff_%d",i);
                (void)strcat(errmsg,"\n");
                (void)strcat(errmsg,colname);
            }
        }
        (void)strcat(errmsg,"\n");
    }

    /* Right, if there are errors to be reported, then do that now */

    if (errmsg[0] != '\0') {
        cpl_msg_error(fctid,"Errors in channel table\n%s",errmsg);
        return(VIR_FATAL);
    }

    /* Get out of here */

    return(VIR_OK);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern cpl_table *vircam_chantab_new(int nord, cpl_table *template) {
    cpl_table *tab;
    int i,oldnorder,null,j,nr;
    char colname[SZCOL];

    /* Create the table */

    tab = cpl_table_duplicate(template);

    /* Fix the table so that it's appropriate for the order in the 
       current fit */

    oldnorder = cpl_table_get_int(tab,"norder",0,&null);
    if (oldnorder > nord) {
        for (i = nord+1; i <= oldnorder; i++) {
            (void)snprintf(colname,SZCOL,"coeff_%d",i);
            cpl_table_erase_column(tab,colname);
        }
    } else if (oldnorder < nord) {
        for (i = oldnorder+1; i <= nord; i++) {
            (void)snprintf(colname,SZCOL,"coeff_%d",i);
            if (cpl_table_has_column(tab,colname)) 
                continue;
            cpl_table_new_column(tab,colname,CPL_TYPE_DOUBLE);
        }
    }
    
    /* Fill in some reasonable initial values for the linearity info */

    nr = cpl_table_get_nrow(tab);
    for (i = 0; i < nr; i++) {
        cpl_table_set_int(tab,"norder",i,nord);
        cpl_table_set_double(tab,"coeff_1",i,1.0);
        for (j = 2; j <= nord; j++) {
            (void)snprintf(colname,SZCOL,"coeff_%d",j);
            cpl_table_set_double(tab,colname,i,0.0);
        }
    }

    /* Get out of here */

    return(tab);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern int vircam_chan_fill(cpl_table *tab, parquet **p, long *np) {
    parquet *pp;
    int i,channum,ixmin,ixmax,iymin,iymax;
    int dcrpix1,dcrpix2,dcd1_1,dcd1_2,dcd2_1,dcd2_2,retval,null;

    /* Verify the table first */

    if (vircam_chantab_verify(tab) != VIR_OK)
        return(VIR_FATAL);

    /* Get the number of rows in the table. This will tell you how many 
       parquet structures to allocate */

    *np = cpl_table_get_nrow(tab);
    *p = cpl_malloc(*np*sizeof(parquet));

    /* Right, start reading them in... */

    for (i = 0; i < *np; i++) {
        pp = *p + i;
        channum = cpl_table_get_int(tab,"channum",i,&null);
        ixmin = cpl_table_get_int(tab,"ixmin",i,&null);
        ixmax = cpl_table_get_int(tab,"ixmax",i,&null);
        iymin = cpl_table_get_int(tab,"iymin",i,&null);
        iymax = cpl_table_get_int(tab,"iymax",i,&null);
        dcrpix1 = cpl_table_get_int(tab,"dcrpix1",i,&null);
        dcrpix2 = cpl_table_get_int(tab,"dcrpix2",i,&null);
        dcd1_1 = cpl_table_get_int(tab,"dcd1_1",i,&null);
        dcd1_2 = cpl_table_get_int(tab,"dcd1_2",i,&null);
        dcd2_1 = cpl_table_get_int(tab,"dcd2_1",i,&null);
        dcd2_2 = cpl_table_get_int(tab,"dcd2_2",i,&null);
            
        /* Right, now fill in the values */

        retval = vircam_chan_init(channum,ixmin,ixmax,iymin,iymax,dcrpix1,
                                  dcrpix2,dcd1_1,dcd1_2,dcd2_1,dcd2_2,pp);
        if (retval != VIR_OK) {
            freespace(*p);
            return(retval);
        }

        /* Fill in the linearity coefficient information */

        retval = vircam_chan_fill_lin(tab,i,pp);
        if (retval != VIR_OK) {
            freespace(*p);
            return(retval);
        }
    }

    /* Get out of here */

    return(VIR_OK);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern void vircam_chan_free(int np, parquet **p) {
    int i;

    /* Loop for and close each entry */

    for (i = 0; i < np; i++)
        vircam_chan_close(*p+i);

    /* Now free the array workspace */

    freespace(*p);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern long vircam_chan_d2r(parquet *p, long l) {
    int lx,ly,kx,ky;
    long k;

    /* Convert detector channel index to readout index. First cartesian 
       coords in detector coords */

    ly = (int)(l/(p->delta_i));
    lx = l - ly*(p->delta_i);

    /* Now cartesian in readout coords */

    kx = (lx - p->Lx)*(p->dcd1_1) + (ly - p->Ly)*(p->dcd1_2);
    ky = (lx - p->Lx)*(p->dcd2_1) + (ly - p->Ly)*(p->dcd2_2);

    /* Now the readout index */

    k = ky*(p->delta_x) + kx;
    return(k);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern long vircam_chan_r2d(parquet *p, long k) {
    int lx,ly,kx,ky;
    long l;

    /* Convert readout index to detector channel index. Do the calculation. 
       First cartesion readout coordinates */

    ky = (int)(k/(p->delta_x));
    kx = k - ky*(p->delta_x);

    /* Now cartesian detector coords */

    lx = (kx*(p->dcd2_2) - ky*(p->dcd1_2))/(p->det_cd) + p->Lx;
    ly = (ky*(p->dcd1_1) - kx*(p->dcd2_1))/(p->det_cd) + p->Ly;

    /* Now the detector index */
    
    l = ly*(p->delta_i) + lx;
    return(l);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern long vircam_chan_r2a(parquet *p, long naxis[2], long k) {
    int lx,ly,kx,ky;
    long l;

    /* Convert readout index to detector channel index. Do the calculation. 
       First cartesion readout coordinates */

    ky = (int)(k/(p->delta_x));
    kx = k - ky*(p->delta_x);

    /* Now cartesian detector coords */

    lx = (kx*(p->dcd2_2) - ky*(p->dcd1_2))/(p->det_cd) + p->Lx + p->ixmin - 1;
    ly = (ky*(p->dcd1_1) - kx*(p->dcd2_1))/(p->det_cd) + p->Ly + p->iymin - 1;

    /* Now the detector index */
    
    l = ly*naxis[0] + lx;
    return(l);
}


/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/


static int vircam_chan_init(int channum, int ixmin, int ixmax, int iymin, 
                            int iymax, int dcrpix1, int dcrpix2, int dcd1_1, 
                            int dcd1_2, int dcd2_1, int dcd2_2, parquet *p) {

    /* Fill in the relevant values */

    p->channum = channum;
    p->ixmin = ixmin;
    p->ixmax = ixmax;
    p->iymin = iymin;
    p->iymax = iymax;
    p->fpix[0] = ixmin;
    p->fpix[1] = iymin;
    p->lpix[0] = ixmax;
    p->lpix[1] = iymax;
    (void)snprintf(p->imsec,SECSIZ,"[%d:%d,%d:%d]",ixmin,ixmax,iymin,iymax);
    p->delta_i = ixmax - ixmin + 1;
    p->delta_j = iymax - iymin + 1;
    if (dcrpix1 < 0 || dcrpix1 > p->delta_i || 
        dcrpix2 < 0 || dcrpix2 > p->delta_j) {
        cpl_msg_error("vircam_chan_init",
                      "nonsense values of dcrpix");
        return(VIR_FATAL);
    }
    p->Lx = dcrpix1 - 1;
    p->Ly = dcrpix2 - 1;
    p->dcd1_1 = dcd1_1;
    p->dcd1_2 = dcd1_2;
    p->dcd2_1 = dcd2_1;
    p->dcd2_2 = dcd2_2;
    p->det_cd = dcd1_1*dcd2_2 - dcd1_2*dcd2_1;
    p->delta_x = abs(dcd1_1*(p->delta_i) + dcd1_2*(p->delta_j));
    p->delta_y = abs(dcd2_1*(p->delta_i) + dcd2_2*(p->delta_j));
    p->data = NULL;
    p->norder = 0;
    p->bb = NULL;
    return(VIR_OK);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

static void vircam_chan_close(parquet *p) {
    
    /* Check to see if the linearity array has been allocated. If so, then
       free the workspace. */

    freespace(p->data);
    freespace(p->bb);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

static int vircam_chan_testpos(cpl_table *intab, const char *column) {
    int nrows,i,*idata;

    /* First check the column type */

    if (cpl_table_get_column_type(intab,column) != CPL_TYPE_INT)
        return(VIR_FATAL);

    /* Now have a look at the actual values and see if they're all 
       positive */

    nrows = cpl_table_get_nrow(intab);
    idata = cpl_table_get_data_int(intab,column);
    for (i = 0; i < nrows; i++)
        if (idata[i] <= 0)
            return(VIR_FATAL);

    /* Otherwise we're OK */

    return(VIR_OK);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

static int vircam_chan_test1(cpl_table *intab, const char *column) {
    int nrows,i,*idata;

    /* First check the column type */

    if (cpl_table_get_column_type(intab,column) != CPL_TYPE_INT)
        return(VIR_FATAL);

    /* Now have a look at the actual values and see if they're all 
       0, 1 or -1 */

    nrows = cpl_table_get_nrow(intab);
    idata = cpl_table_get_data_int(intab,column);
    for (i = 0; i < nrows; i++)
        if (idata[i] != 0 && idata[i] != -1 && idata[i] != 1)
            return(VIR_FATAL);

    /* Otherwise we're OK */

    return(VIR_OK);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

static int vircam_chan_fill_lin(cpl_table *intab, int row, parquet *p) {
    int i,null;
    char colname[SZCOL];
    
    /* First get the order of the fit */

    p->norder = cpl_table_get_int(intab,"norder",row,&null);

    /* Get some workspace for the coefficients */

    p->bb = cpl_malloc(p->norder*sizeof(double));

    /* Now loop for the order numbers and read off the coefficents */

    for (i = 1; i <= p->norder; i++) {
        (void)snprintf(colname,SZCOL,"coeff_%d",i);
        p->bb[i-1] = cpl_table_get_double(intab,colname,row,&null);
    }
    return(VIR_OK);
}
        

/*

$Log: vircam_channel.c,v $
Revision 1.12  2007/11/20 09:36:51  jim
changed column 'qualfit' to 'lin_10000_err'

Revision 1.11  2007/10/25 17:34:00  jim
Modified to remove lint warnings

Revision 1.10  2007/10/19 09:25:09  jim
Fixed problems with missing includes

Revision 1.9  2007/03/01 12:42:41  jim
Modified slightly after code checking

Revision 1.8  2006/11/10 09:26:43  jim
Modified vircam_channel_new to copy the basic structure of a pre-existing
template.

Revision 1.7  2006/09/08 09:21:02  jim
Added routine vircam_chantab_new

Revision 1.6  2006/07/04 09:19:03  jim
replaced all sprintf statements with snprintf

Revision 1.5  2006/06/09 11:32:59  jim
A few more minor fixes for lint

Revision 1.4  2006/06/09 11:26:25  jim
Small changes to keep lint happy

Revision 1.3  2006/03/15 10:43:41  jim
Fixed a few things

Revision 1.2  2006/03/03 14:29:45  jim
Modified definition of vir_fits and channel table

Revision 1.1  2006/02/18 11:43:45  jim
new file


*/

---