vircam_platesol.c

/* $Id: vircam_platesol.c,v 1.22 2008/07/10 13:05:53 jim Exp $
 *
 * This file is part of the VIRCAM Pipeline
 * Copyright (C) 2005 Cambridge Astronomy Survey Unit
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

/*
 * $Author: jim $
 * $Date: 2008/07/10 13:05:53 $
 * $Revision: 1.22 $
 * $Name:  $
 */

/* Includes */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <cpl.h>
#include <math.h>

#include "vircam_mods.h"
#include "vircam_utils.h"
#include "vircam_wcsutils.h"
#include "vircam_stats.h"
#include "vircam_pfits.h"

static double *work = NULL;
static unsigned char *iwork = NULL;

static int vircam_plate6(double *xpos, double *ypos, double *xi, double *eta,
                         unsigned char *flag, int nstds, double *a, double *b,
                         double *c, double *d, double *e, double *f);
static int vircam_plate4(double *xpos, double *ypos, double *xi, double *eta,
                         unsigned char *flag, int nstds, double *a, double *b,
                         double *c, double *d, double *e, double *f);
static void tidy(void);

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

extern int vircam_platesol(cpl_propertylist *plist, cpl_propertylist *tlist,
                           cpl_table *matchedstds, int nconst, int shiftan, 
                           int *status) {
    int nstds,nc2,i,niter,nrej,ngood,nreq=6,xcol,ycol;
    long n1,n2;
    const char *fctid = "vircam_platesol";
    float *tptr;
    double *xptr,*yptr,*xptr2,*yptr2,*ra,*dec,*xiptr,*etaptr,*wptr,averr;
    double r1,r2,d1,d2,newcrval1,newcrval2,a,b,c,d,e,f,xifit,etafit,dxi,deta;
    double crpix1,crpix2,xi,eta,rac_before,rac_after,decc_before,decc_after;
    double xcen,ycen,oldtheta,scale,oldcrpix1,oldcrpix2,*crdata;
    unsigned char *isbad,*wptr2;
    const char *reqcols[] = {"X_coordinate","Y_coordinate","xpredict",
                             "ypredict","RA","Dec"};
    char key[9];
    cpl_wcs *wcs;
    cpl_array *cr;
    cpl_matrix *crm;

    /* Inherited status */

    if (*status != VIR_OK)
        return(*status);

    /* Check the value of nconst */

    if (nconst != 4 && nconst != 6) {
        cpl_msg_error(fctid,"Value of nconst = %d is unsupported",nconst);
        FATAL_ERROR
    }

    /* How many standards are in the input matched standards table? */

    nstds = cpl_table_get_nrow(matchedstds);
    nc2 = nconst/2;
    if (nstds < nc2) {
        cpl_msg_error(fctid,
                      "Too few standards (%d) in table for %d coefficient fit",
                      nstds,nconst);
        FATAL_ERROR
    }

    /* Check that the matched standards table has all the required columns */

    for (i = 0; i < nreq; i++) {
        if (cpl_table_has_column(matchedstds,reqcols[i]) != 1) {
            cpl_msg_error(fctid,"Matched standards table missing column %s\n",
                          reqcols[i]);
            FATAL_ERROR
        }
    }
            
    /* Get some workspace now */

    work = cpl_malloc(10*nstds*sizeof(*work));
    iwork = cpl_calloc(3*nstds,sizeof(*isbad));
    xptr = work;
    yptr = work + nstds;
    xptr2 = work + 2*nstds;
    yptr2 = work + 3*nstds;
    ra = work + 4*nstds;
    dec = work + 5*nstds;
    xiptr = work + 6*nstds;
    etaptr = work + 7*nstds;
    wptr = work + 8*nstds;
    isbad = iwork;
    wptr2 = iwork + nstds;
    
    /* Get the data from the table and put it all into double precision
       arrays */

    tptr = cpl_table_get_data_float(matchedstds,"X_coordinate");
    for (i = 0; i < nstds; i++)
        xptr[i] = (double)tptr[i];
    tptr = cpl_table_get_data_float(matchedstds,"Y_coordinate");
    for (i = 0; i < nstds; i++)
        yptr[i] = (double)tptr[i];
    tptr = cpl_table_get_data_float(matchedstds,"xpredict");
    for (i = 0; i < nstds; i++)
        xptr2[i] = (double)tptr[i];
    tptr = cpl_table_get_data_float(matchedstds,"ypredict");
    for (i = 0; i < nstds; i++)
        yptr2[i] = (double)tptr[i];
    tptr = cpl_table_get_data_float(matchedstds,"RA");
    for (i = 0; i < nstds; i++)
        ra[i] = (double)tptr[i];
    tptr = cpl_table_get_data_float(matchedstds,"Dec");
    for (i = 0; i < nstds; i++)
        dec[i] = (double)tptr[i];

    /* If you want to shift the RA and Dec of the tangent point, then
       do that now */

    if (shiftan) {
        wcs = cpl_wcs_new_from_propertylist(plist);
        cr = cpl_wcs_get_crval(wcs);
        crdata = cpl_array_get_data_double(cr);
        for (i = 0; i < nstds; i++) {
            vircam_xytoradec(wcs,xptr[i],yptr[i],&r1,&d1);
            vircam_xytoradec(wcs,xptr2[i],yptr2[i],&r2,&d2);
            xiptr[i] = r2 - r1;
            etaptr[i] = d2 - d1;
        }
        r1 = vircam_dmed(xiptr,NULL,nstds);
        d1 = vircam_dmed(etaptr,NULL,nstds);
        newcrval1 = crdata[0] + r1;
        newcrval2 = crdata[1] + d1;
        cpl_propertylist_update_double(plist,"CRVAL1",newcrval1);
        cpl_propertylist_update_double(plist,"CRVAL2",newcrval2);
        cpl_wcs_delete(wcs);
    }

    /* Calculate the central RA and Dec */

    wcs = cpl_wcs_new_from_propertylist(plist);
    (void)vircam_pfits_get_naxis1(plist,&n1);
    (void)vircam_pfits_get_naxis2(plist,&n2);
    cr = cpl_wcs_get_crpix(wcs);
    crdata = cpl_array_get_data_double(cr);
    oldcrpix1 = crdata[0];
    oldcrpix2 = crdata[1];
    xcen = 0.5*(double)n1;
    ycen = 0.5*(double)n2;
    vircam_xytoradec(wcs,xcen,ycen,&rac_before,&decc_before);

    /* Right, calculate xi and eta for each of the standards */

    for (i = 0; i < nstds; i++) {
        vircam_radectoxieta(wcs,ra[i],dec[i],&xi,&eta);
        xiptr[i] = xi;
        etaptr[i] = eta;
    }

    /* Right, now loop for maximum number of iterations or until
       convergence */

    niter = 0;
    while (niter >= 0) {

        /* Do a plate solution */

        switch (nconst) {
        case 6:
            *status = vircam_plate6(xptr,yptr,xiptr,etaptr,isbad,nstds,&a,&b,
                                    &c,&e,&d,&f);
            break;
        case 4:
            *status = vircam_plate4(xptr,yptr,xiptr,etaptr,isbad,nstds,&a,&b,
                                    &c,&e,&d,&f);
            break;
        default:
            *status = vircam_plate6(xptr,yptr,xiptr,etaptr,isbad,nstds,&a,&b,
                                    &c,&e,&d,&f);
            break;
        }
        if (*status != VIR_OK) {
            cpl_msg_error(fctid,"Plate constant solution failed");
            tidy();
            FATAL_ERROR
        }

        /* Now look at the residuals and see if any should be rejected */

        for (i = 0; i < nstds; i++) {
            xifit = xptr[i]*a + yptr[i]*b + c;
            etafit = xptr[i]*d + yptr[i]*e + f;
            dxi = fabs(xifit - xiptr[i]);
            deta = fabs(etafit - etaptr[i]);
            wptr[i*2] = dxi;
            wptr[i*2+1] = deta;
            wptr2[i*2] = isbad[i];
            wptr2[i*2+1] = isbad[i];
        }
        averr = vircam_dmed(wptr,wptr2,2*nstds);
        averr *= 1.48;
        if (niter == 3)
            break;
        nrej = 0;
        ngood = 0;
        for (i = 0; i < nstds; i++) {
            if (!isbad[i] && (wptr[i*2] > 3.0*averr || wptr[i*2+1] > 3.0*averr))                nrej++;
            if (!isbad[i])
                ngood++;
        }
        ngood -= nrej;
        if (nrej == 0 || ngood < nconst)
            break;
        for (i = 0; i < nstds; i++) {
            if (!isbad[i] && (wptr[i*2] > 3.0*averr || wptr[i*2+1] > 3.0*averr))                isbad[i] = 1;
        }
        niter++;
    }

    /* Convert values to degrees now */

    crpix1 = (e*c - b*f)/(d*b - e*a);
    crpix2 = (a*f - d*c)/(d*b - e*a);
    a *= DEGRAD;
    b *= DEGRAD;
    d *= DEGRAD;
    e *= DEGRAD;

    /* Number of good points fit and average error in arcsec*/

    ngood = 0;
    for (i = 0; i < nstds; i++)
        if (! isbad[i])
            ngood++;
    averr *= DEGRAD*3600.0;

    /* Right, now update the header */

    cpl_propertylist_update_double(plist,"CRPIX1",crpix1);
    cpl_propertylist_update_double(plist,"CRPIX2",crpix2);
    cpl_propertylist_update_double(plist,"CD1_1",a);
    cpl_propertylist_update_double(plist,"CD1_2",b);
    cpl_propertylist_update_double(plist,"CD2_1",d);
    cpl_propertylist_update_double(plist,"CD2_2",e);
    cpl_propertylist_update_int(plist,"ESO DRS NUMBRMS",ngood);
    cpl_propertylist_set_comment(plist,"ESO DRS NUMBRMS",
                                 "Number of stars in WCS fit");
    cpl_propertylist_update_float(plist,"ESO DRS STDCRMS",(float)averr);
    cpl_propertylist_set_comment(plist,"ESO DRS STDCRMS",
                                 "[arcsec] Average error in WCS fit");

    /* Calculate the central RA and Dec again */

    crm = cpl_wcs_get_cd(wcs);
    crdata = cpl_matrix_get_data(crm);
    oldtheta = 0.5*(fabs(atan2(crdata[1],crdata[0])) + 
                    fabs(atan2(crdata[2],crdata[3])));
    cpl_wcs_delete(wcs);
    wcs = cpl_wcs_new_from_propertylist(plist);
    vircam_xytoradec(wcs,xcen,ycen,&rac_after,&decc_after);
    xcen = 3600.0*(rac_after - rac_before);
    ycen = 3600.0*(decc_after - decc_before);
    cpl_propertylist_update_float(plist,"ESO DRS WCSRAOFF",(float)xcen);
    cpl_propertylist_set_comment(plist,"ESO DRS WCSRAOFF",
                                 "[arcsec] cenpix RA_after - RA_before)");
    cpl_propertylist_update_float(plist,"ESO DRS WCSDECOFF",(float)ycen);
    cpl_propertylist_set_comment(plist,"ESO DRS WCSDECOFF",
                                 "[arcsec] cenpix Dec_after - Dec_before)");

    /* Update the table header */

    if (tlist != NULL) {
        xcol = vircam_findcol(tlist,"X");
        ycol = vircam_findcol(tlist,"Y");
        if (xcol != -1 && ycol != -1) {
            snprintf(key,9,"TCRPX%d",xcol);
            cpl_propertylist_update_double(tlist,key,crpix1);
            snprintf(key,9,"TCRPX%d",ycol);
            cpl_propertylist_update_double(tlist,key,crpix2);
            snprintf(key,9,"TC%d_%d",xcol,xcol);
            cpl_propertylist_update_double(tlist,key,a);
            snprintf(key,9,"TC%d_%d",xcol,ycol);
            cpl_propertylist_update_double(tlist,key,b);
            snprintf(key,9,"TC%d_%d",ycol,xcol);
            cpl_propertylist_update_double(tlist,key,d);
            snprintf(key,9,"TC%d_%d",ycol,ycol);
            cpl_propertylist_update_double(tlist,key,e);
            cpl_propertylist_update_int(tlist,"ESO DRS NUMBRMS",ngood);
            cpl_propertylist_set_comment(tlist,"ESO DRS NUMBRMS",
                                         "Number of stars in WCS fit");
            cpl_propertylist_update_float(tlist,"ESO DRS STDCRMS",(float)averr);
            cpl_propertylist_set_comment(tlist,"ESO DRS STDCRMS",
                                         "[arcsec] Average error in WCS fit");
            cpl_propertylist_update_float(tlist,"ESO DRS WCSRAOFF",(float)xcen);
            cpl_propertylist_set_comment(tlist,"ESO DRS WCSRAOFF",
                                         "[arcsec] cenpix RA_after - RA_before)");
            cpl_propertylist_update_float(tlist,"ESO DRS WCSDECOFF",(float)ycen);
            cpl_propertylist_set_comment(tlist,"ESO DRS WCSDECOFF",
                                         "[arcsec] cenpix Dec_after - Dec_before)");
        }
    }

    /* Back-calculate a crval for the old value of crpix. Compare to the 
       WCS crval and write to QC header */

    cr = cpl_wcs_get_crval(wcs);
    crdata = cpl_array_get_data_double(cr);
    vircam_xytoradec(wcs,oldcrpix1,oldcrpix2,&rac_after,&decc_after);
    rac_after -= crdata[0];
    decc_after -= crdata[1];
    cpl_propertylist_update_float(plist,"ESO QC WCS_DCRVAL1",(float)rac_after);
    cpl_propertylist_set_comment(plist,"ESO QC WCS_DCRVAL1",
                                 "[deg] change in crval1");
    cpl_propertylist_update_float(plist,"ESO QC WCS_DCRVAL2",(float)decc_after);
    cpl_propertylist_set_comment(plist,"ESO QC WCS_DCRVAL2",
                                 "[deg] change in crval2");

    /* Work out the change in the rotation */

    crm = cpl_wcs_get_cd(wcs);
    crdata = cpl_matrix_get_data(crm);
    oldtheta = 0.5*(fabs(atan2(crdata[1],crdata[0])) + 
                    fabs(atan2(crdata[2],crdata[3]))) - oldtheta;
    oldtheta *= DEGRAD;
    cpl_propertylist_update_float(plist,"ESO QC WCS_DTHETA",(float)oldtheta);
    cpl_propertylist_set_comment(plist,"ESO QC WCS_DTHETA",
                                 "[deg] change in rotation");
    
    /* Work out the mean plate scale */

    scale = 1800.0*(sqrt(pow(crdata[0],2.0) +  pow(crdata[1],2.0)) + 
                    sqrt(pow(crdata[2],2.0) +  pow(crdata[3],2.0)));
    cpl_propertylist_update_float(plist,"ESO QC WCS_SCALE",(float)scale);
    cpl_propertylist_set_comment(plist,"ESO QC WCS_SCALE",
                                 "[arcsec] mean plate scale");

    /* Work out the shear if this new WCS */

    oldtheta = fabs(atan2(crdata[1],crdata[0])) - 
        fabs(atan2(crdata[2],crdata[3]));
    cpl_propertylist_update_float(plist,"ESO QC WCS_SHEAR",(float)oldtheta);
    cpl_propertylist_set_comment(plist,"ESO QC WCS_SHEAR",
                                 "[deg] abs(xrot) - abs(yrot)");

    /* Now just add in the RMS */

    cpl_propertylist_update_float(plist,"ESO QC WCS_RMS",(float)averr);
    cpl_propertylist_set_comment(plist,"ESO QC WCS_RMS",
                                 "[arcsec] Average error in WCS fit");

    /* Right, get out of here now... */

    cpl_wcs_delete(wcs);
    tidy();
    GOOD_STATUS
}

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

static int vircam_plate6(double *xpos, double *ypos, double *xi, double *eta,
                         unsigned char *flag, int nstds, double *a, double *b,
                         double *c, double *d, double *e, double *f) {
    double sx1sq,sy1sq,sx1y1,sx1x2,sy1x2;
    double sy1y2,sx1y2,xposmean,yposmean,ximean,etamean,xx1,yy1,xx2,yy2;
    int i,ngood,nbad;

    /* Is it worthwhile even being here? */

    (void)vircam_sumbpm(flag,nstds,&nbad);
    ngood = nstds - nbad;
    if (ngood < 2)
        return(VIR_FATAL);

    /* Initialise all the counters and summations */

    sx1sq = 0.0;
    sy1sq = 0.0;
    sx1y1 = 0.0;
    sx1x2 = 0.0;
    sy1x2 = 0.0;
    sy1y2 = 0.0;
    sx1y2 = 0.0;
    xposmean = 0.0;
    yposmean = 0.0;
    ximean = 0.0;
    etamean = 0.0;

    /* Find means in each coordinate system */

    xposmean = vircam_dmean(xpos,flag,nstds);
    yposmean = vircam_dmean(ypos,flag,nstds);
    ximean = vircam_dmean(xi,flag,nstds);
    etamean = vircam_dmean(eta,flag,nstds);

    /* Now accumulate the sums */

    for (i = 0; i < nstds; i++) {
        if (!flag[i]) {
            xx1 = xpos[i] - xposmean;
            yy1 = ypos[i] - yposmean;
            xx2 = xi[i] - ximean;
            yy2 = eta[i] - etamean;
            sx1sq += xx1*xx1;
            sy1sq += yy1*yy1;
            sx1y1 += xx1*yy1;
            sx1x2 += xx1*xx2;
            sy1x2 += yy1*xx2;
            sy1y2 += yy1*yy2;
            sx1y2 += xx1*yy2;
        }
    }

    /* Do solution for X */

    *a = (sx1y1*sy1x2 - sx1x2*sy1sq)/(sx1y1*sx1y1 - sx1sq*sy1sq);
    *b = (sx1x2*sx1y1 - sx1sq*sy1x2)/(sx1y1*sx1y1 - sx1sq*sy1sq);
    *c = -xposmean*(*a) - yposmean*(*b) + ximean;

    /* Now the solution for Y */

    *d = (sx1y1*sx1y2 - sy1y2*sx1sq)/(sx1y1*sx1y1 - sy1sq*sx1sq);
    *e = (sy1y2*sx1y1 - sy1sq*sx1y2)/(sx1y1*sx1y1 - sy1sq*sx1sq);
    *f = -xposmean*(*e) - yposmean*(*d) + etamean;

    /* Get outta here */

    return(VIR_OK);
}

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

static int vircam_plate4(double *xpos, double *ypos, double *xi, double *eta,
                         unsigned char *flag, int nstds, double *a, double *b,
                         double *c, double *d, double *e, double *f) {
    double sx1sq,sy1sq,sx1x2,sy1x2,sy1y2,sx1y2,xposmean,yposmean;
    double ximean,etamean,xx1,yy1,xx2,yy2,det,num,denom,theta,mag;
    double stheta,ctheta;
    int i,ngood,nbad;

    /* Is it worthwhile even being here? */

    (void)vircam_sumbpm(flag,nstds,&nbad);
    ngood = nstds - nbad;
    if (ngood < 2)
        return(VIR_FATAL);

    /* Initialise all the counters and summations */

    sx1sq = 0.0;
    sy1sq = 0.0;
    sx1x2 = 0.0;
    sy1x2 = 0.0;
    sy1y2 = 0.0;
    sx1y2 = 0.0;
    xposmean = 0.0;
    yposmean = 0.0;
    ximean = 0.0;
    etamean = 0.0;

    /* Find means in each coordinate system */

    xposmean = vircam_dmean(xpos,flag,nstds);
    yposmean = vircam_dmean(ypos,flag,nstds);
    ximean = vircam_dmean(xi,flag,nstds);
    etamean = vircam_dmean(eta,flag,nstds);

    /* Now accumulate the sums */

    for (i = 0; i < nstds; i++) {
        if (!flag[i]) {
            xx1 = xpos[i] - xposmean;
            yy1 = ypos[i] - yposmean;
            xx2 = xi[i] - ximean;
            yy2 = eta[i] - etamean;
            sx1sq += xx1*xx1;
            sy1sq += yy1*yy1;
            sx1x2 += xx1*xx2;
            sy1x2 += yy1*xx2;
            sy1y2 += yy1*yy2;
            sx1y2 += xx1*yy2;
        }
    }

    /* Compute the rotation angle */

    det = sx1x2*sy1y2 - sy1x2*sx1y2;
    if (det < 0.0) {
        num = sy1x2 + sx1y2;
        denom = -sx1x2 + sy1y2;
    } else {
        num = sy1x2 - sx1y2;
        denom = sx1x2 + sy1y2;
    }
    if (num == 0.0 && denom == 0.0) {
        theta = 0.0;
    } else {
        theta = atan2(num,denom);
        if (theta < 0.0)
            theta += M_TWOPI;
    }

    /* Compute magnification factor */

    ctheta = cos(theta);
    stheta = sin(theta);
    num = denom*ctheta  + num*stheta;
    denom = sx1sq + sy1sq;
    if (denom <= 0.0) {
        mag = 1.0;
    } else {
        mag = num/denom;
    }

    /* Compute coeffs */

    if (det < 0.0) {
        *a = -mag*ctheta;
        *e = mag*stheta;
    } else {
        *a = mag*ctheta;
        *e = -mag*stheta;
    }
    *b = mag*stheta;
    *d = mag*ctheta;
    *c = -xposmean*(*a) - yposmean*(*b) + ximean;
    *f = -xposmean*(*e) - yposmean*(*d) + etamean;

    /* Get outta here */

    return(VIR_OK);
}

static void tidy(void) {

    freespace(work);
    freespace(iwork);
}

/*

$Log: vircam_platesol.c,v $
Revision 1.22  2008/07/10 13:05:53  jim
Modified to use v4.2 version of cpl_wcs

Revision 1.21  2008/05/06 08:40:10  jim
Modified to use cpl_wcs interface

Revision 1.20  2008/01/22 19:46:10  jim
Fixed sign error in plate4

Revision 1.19  2007/10/25 17:34:01  jim
Modified to remove lint warnings

Revision 1.18  2007/05/03 11:15:33  jim
Fixed little problem with table wcs

Revision 1.17  2007/05/02 09:12:30  jim
Modified to update table header WCS keywords

Revision 1.16  2007/03/29 12:19:39  jim
Little changes to improve documentation

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

Revision 1.14  2007/01/17 23:54:01  jim
Plugged some memory leaks

Revision 1.13  2007/01/08 19:12:03  jim
Fixed shear comment so that it doesn't overflow

Revision 1.12  2006/10/02 13:43:50  jim
Added missing .h file

Revision 1.11  2006/08/11 12:45:41  jim
Modified to use wcslib

Revision 1.10  2006/06/14 14:33:31  jim
Fixed units of WCS_SCALE

Revision 1.9  2006/06/13 14:09:38  jim
Made some of the QC header values single precision

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

Revision 1.7  2006/05/26 15:05:46  jim
Fixed column names for input table

Revision 1.6  2006/03/23 21:18:53  jim
Minor changes mainly to comment headers

Revision 1.5  2006/03/22 14:05:37  jim
fixed a little bug

Revision 1.4  2006/03/22 13:58:32  jim
Cosmetic fixes to keep lint happy

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

Revision 1.2  2006/02/22 14:19:52  jim
Modified to check for the existence of the columns in the matched standards
table.

Revision 1.1  2006/02/18 11:52:34  jim
new file


*/

---