vircam_match.c

/* $Id: vircam_match.c,v 1.13 2007/10/19 09:25:10 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: 2007/10/19 09:25:10 $
 * $Revision: 1.13 $
 * $Name:  $
 */

/* Includes */

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

#include <math.h>

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

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

#define NX 2048
#define NY 2048
#define NGRIDMAX 31

static cpl_table *vircam_mkmstd_table(cpl_table *objtab, cpl_table *stdstab);

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

extern int vircam_matchxy(cpl_table *progtab, cpl_table *template, float srad,
                          float *xoffset, float *yoffset, int *nm, 
                          int *status) {
    cpl_propertylist *p;
    float *xprog,*yprog,*xtemp,*ytemp,aveden,errlim,xoffbest,yoffbest,xoff;
    float yoff,x,y,*xoffs,*yoffs;
    const char *fctid = "vircam_matchxy";
    int nprog,ntemp,ngrid,ngrid2,ibest,ig,jg,nmatch,k,jm;

    /* Inherited status */

    *xoffset = 0.0;
    *yoffset = 0.0;
    *nm = 0;
    if (*status != VIR_OK)
        return(*status);

    /* Check the size of each of the tables */

    nprog = cpl_table_get_nrow(progtab);
    ntemp = cpl_table_get_nrow(template);
    if (nprog == 0) {
        cpl_msg_warning(fctid,"Program table has no rows");
        WARN_RETURN
    } else if (ntemp == 0) {
        cpl_msg_warning(fctid,"Template table has no rows");
        WARN_RETURN
    }

    /* First, sort the two tables by the Y coordinate */

    p = cpl_propertylist_new();
    cpl_propertylist_append_bool(p,"Y_coordinate",0);
    if (cpl_table_sort(progtab,p) != CPL_ERROR_NONE) {
        cpl_propertylist_delete(p);
        FATAL_ERROR
    }
    if (cpl_table_sort(template,p) != CPL_ERROR_NONE) {
        cpl_propertylist_delete(p);
        FATAL_ERROR
    }
    cpl_propertylist_delete(p);

    /* Get the x,y coordinates for each table */

    xprog = cpl_table_get_data_float(progtab,"X_coordinate");
    yprog = cpl_table_get_data_float(progtab,"Y_coordinate");
    xtemp = cpl_table_get_data_float(template,"X_coordinate");
    ytemp = cpl_table_get_data_float(template,"Y_coordinate");
    if (xprog == NULL || yprog == NULL || xtemp == NULL || ytemp == NULL)
        FATAL_ERROR

    /* Calculate the error limit and the number of grid points */

    aveden = (float)ntemp/(float)(NX*NY);
    errlim = 1.0/sqrt(4.0*M_PI*aveden);
    errlim = min(errlim,5.0);
    ngrid = (int)(srad/errlim);
    ngrid = (ngrid/2)*2 + 1;
    ngrid = max(5,min(NGRIDMAX,ngrid));
    ngrid2 = ngrid/2 + 1;

    /* Now search for the best solution */

    ibest = 0;
    xoffbest = 0.0;
    yoffbest = 0.0;
    for (ig = -ngrid2; ig <= ngrid2; ig++) {
        xoff = (float)ig*errlim*M_SQRT2;
        for (jg = -ngrid2; jg <= ngrid2; jg++) {
            yoff = (float)jg*errlim*M_SQRT2;
            nmatch = 0;
            for (k = 0; k < nprog; k++) {
                x = xprog[k] + xoff;
                y = yprog[k] + yoff;
                if (vircam_fndmatch(x,y,xtemp,ytemp,ntemp,errlim) > -1) 
                    nmatch++;
            }
            if (nmatch > ibest) {
                ibest = nmatch;
                xoffbest = xoff;
                yoffbest = yoff;
            }
        }
    }

    /* Now allocate some workspace so that you can calculate a good median
       coordinate difference */

    xoffs = cpl_malloc(nprog*sizeof(*xoffs));
    yoffs = cpl_malloc(nprog*sizeof(*yoffs));

    /* Now get the coordinate differences and find the medians */

    nmatch = 0;
    for (k = 0; k < nprog; k++) {
        x = xprog[k] + xoffbest;
        y = yprog[k] + yoffbest;
        jm = vircam_fndmatch(x,y,xtemp,ytemp,ntemp,errlim);
        if (jm > -1) {
            xoffs[nmatch] = xtemp[jm] - xprog[k];
            yoffs[nmatch] = ytemp[jm] - yprog[k];
            nmatch++;
        }
    }
    if (nmatch > 0) {
        *xoffset = vircam_med(xoffs,NULL,nmatch);
        *yoffset = vircam_med(yoffs,NULL,nmatch);
    } else {
        *xoffset = 0.0;
        *yoffset = 0.0;
    }
    *nm = nmatch;

    /* Tidy and exit */

    freespace(xoffs);
    freespace(yoffs);
    GOOD_STATUS
}

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

extern int vircam_matchstds(cpl_table *objtab, cpl_table *stdstab, float srad,
                            cpl_table **outtab, int *status) {
    const char *fctid = "vircam_matchstds";
    char *colname;
    int nobj,nstd,ngrid,ngrid2,ibest,ig,jg,nmatch,k,*matches,jm,l,dont,null;
    float *xstd,*ystd,*xobj,*yobj,aveden,errlim,xoffbest,yoffbest,*xoffs;
    float *yoffs,x,y,x2,y2,r2,x1,y1,r1,xoffmed,sigx,yoffmed,sigy,xoff,yoff;
    cpl_propertylist *p;
    cpl_table *mstds;

    /* Inherited status */

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

    /* Check the size of each of the tables */

    nobj = cpl_table_get_nrow(objtab);
    nstd = cpl_table_get_nrow(stdstab);
    if (nobj == 0) {
        cpl_msg_warning(fctid,"Object table has no rows");
        mstds = vircam_mkmstd_table(objtab,stdstab);
        *outtab = cpl_table_extract_selected(mstds);
        cpl_table_delete(mstds);
        WARN_RETURN
    } else if (nstd == 0) {
        cpl_msg_warning(fctid,"Standards RA/DEC table has no rows");
        mstds = vircam_mkmstd_table(objtab,stdstab);
        *outtab = cpl_table_extract_selected(mstds);
        cpl_table_delete(mstds);
        WARN_RETURN
    }

    /* First, sort the two tables by the Y coordinate */

    p = cpl_propertylist_new();
    cpl_propertylist_append_bool(p,"Y_coordinate",0);
    if (cpl_table_sort(objtab,p) != CPL_ERROR_NONE) {
        cpl_propertylist_delete(p);
        FATAL_ERROR
    }
    cpl_propertylist_erase(p,"Y_coordinate");
    cpl_propertylist_append_bool(p,"ypredict",0);
    if (cpl_table_sort(stdstab,p) != CPL_ERROR_NONE) {
        cpl_propertylist_delete(p);
        FATAL_ERROR
    }
    cpl_propertylist_delete(p);

    /* Get the x,y coordinates for each table */

    xobj = cpl_table_get_data_float(objtab,"X_coordinate");
    yobj = cpl_table_get_data_float(objtab,"Y_coordinate");
    xstd = cpl_table_get_data_float(stdstab,"xpredict");
    ystd = cpl_table_get_data_float(stdstab,"ypredict");
    if (xstd == NULL || ystd == NULL || xobj == NULL || yobj == NULL)
        FATAL_ERROR

    /* Calculate the error limit and the number of grid points */

    aveden = (float)nstd/(float)(NX*NY);
    errlim = 1.0/sqrt(4.0*M_PI*aveden);
    errlim = min(errlim,5.0);
    ngrid = (int)(srad/errlim);
    ngrid = (ngrid/2)*2 + 1;
    ngrid = max(5,min(NGRIDMAX,ngrid));
    ngrid2 = ngrid/2 + 1;

    /* Now search for the best solution */

    ibest = 0;
    xoffbest = 0.0;
    yoffbest = 0.0;
    for (ig = -ngrid2; ig <= ngrid2; ig++) {
        xoff = (float)ig*errlim*M_SQRT2;
        for (jg = -ngrid2; jg <= ngrid2; jg++) {
            yoff = (float)jg*errlim*M_SQRT2;
            nmatch = 0;
            for (k = 0; k < nobj; k++) {
                x = xobj[k] + xoff;
                y = yobj[k] + yoff;
                if (vircam_fndmatch(x,y,xstd,ystd,nstd,errlim) > -1) 
                    nmatch++;
            }
            if (nmatch > ibest) {
                ibest = nmatch;
                xoffbest = xoff;
                yoffbest = yoff;
            }
        }
    }

    /* Now allocate some workspace so that you can calculate a good median
       coordinate difference */

    xoffs = cpl_malloc(nstd*sizeof(*xoffs));
    yoffs = cpl_malloc(nstd*sizeof(*yoffs));
    matches = cpl_malloc(nstd*sizeof(*matches));
    for (k = 0; k < nstd; k++)
        matches[k] = -1;

    /* Now get the best matches */

    nmatch = 0;
    for (k = 0; k < nstd; k++) {
        x = xstd[k] - xoffbest;
        y = ystd[k] - yoffbest;
        jm = vircam_fndmatch(x,y,xobj,yobj,nobj,errlim);
        if (jm > -1) {
            dont = 0;
            x2 = xobj[jm] - x;
            y2 = yobj[jm] - y;
            r2 = sqrt(x2*x2 + y2*y2);
            for (l = 0; l < nstd; l++) {
                if (matches[l] == jm) {
                    x1 = xobj[jm] - (xstd[l] - xoffbest);
                    y1 = yobj[jm] - (ystd[l] - yoffbest);
                    r1 = sqrt(x1*x1 + y1*y1);
                    if (r2 < r1) 
                        matches[l] = -1;
                    else
                        dont = 1;
                    break;
                }
            }
            if (dont == 0)
                matches[k] = jm;
        }
    }

    /* Now get the coordinate difference for the best matches */

    for (k = 0; k < nstd; k++) {
        jm = matches[k];
        if (jm != -1) {
            xoffs[nmatch] = xobj[jm] - xstd[k];
            yoffs[nmatch] = yobj[jm] - ystd[k];
            nmatch++;
        }
    }
    if (nmatch == 0) {
        xoffmed = 0.0;
        sigx = 1.0;
        yoffmed = 0.0;
        sigy = 1.0;
    } else {
        vircam_medmad(xoffs,NULL,nmatch,&xoffmed,&sigx);
        sigx *= 1.48;
        vircam_medmad(yoffs,NULL,nmatch,&yoffmed,&sigy);
        sigy *= 1.48;
    }

    /* Now go through one final time with a reduced error box and get
       the final matches */

    errlim = 3.0*max(sigx,sigy);
    for (k = 0; k < nstd; k++)
        matches[k] = -1;
    for (k = 0; k < nstd; k++) {
        x = xstd[k] + xoffmed;
        y = ystd[k] + yoffmed;
        jm = vircam_fndmatch(x,y,xobj,yobj,nobj,errlim);
        if (jm > -1) {
            dont = 0;
            x2 = xobj[jm] - x;
            y2 = yobj[jm] - y;
            r2 = sqrt(x2*x2 + y2*y2);
            for (l = 0; l < nstd; l++) {
                if (matches[l] == jm) {
                    x1 = xobj[jm] - (xstd[l] + xoffmed);
                    y1 = yobj[jm] - (ystd[l] + yoffmed);
                    r1 = sqrt(x1*x1 + y1*y1);
                    if (r2 < r1) 
                        matches[l] = -1;
                    else
                        dont = 1;
/*                  break; */
                }
            }
            if (dont == 0)
                matches[k] = jm;
        }
    }
    jm = matches[1];

    /* Make a copy of the standards table and add all the columns from the
       object catalogue to it. Ingore the RA and DEC columns in the catalogue
       as the standards table will already have these.*/

    mstds = cpl_table_duplicate(stdstab);
    colname = (char *)cpl_table_get_column_name(objtab);
    while (colname != NULL) {
        if (strcmp(colname,"RA") && strcmp(colname,"DEC"))
            cpl_table_new_column(mstds,colname,
                                 cpl_table_get_column_type(objtab,colname));
        colname = (char *)cpl_table_get_column_name(NULL);
    }
    cpl_table_unselect_all(mstds);

    /* Now go through and find the matches */

    for (k = 0; k < nstd; k++) {
        jm = matches[k];
        if (jm != -1) {
            colname = (char *)cpl_table_get_column_name(objtab);
            while (colname != NULL) {
                if (!strcmp(colname,"RA") || !strcmp(colname,"DEC")) {
                    colname = (char *)cpl_table_get_column_name(NULL);
                    continue;
                }
                null = 0;
                switch (cpl_table_get_column_type(objtab,colname)) {
                case CPL_TYPE_INT:
                    cpl_table_set_int(mstds,colname,k,
                                      cpl_table_get_int(objtab,colname,jm,
                                                        &null));
                    break;
                case CPL_TYPE_FLOAT:
                    cpl_table_set_float(mstds,colname,k,
                                        cpl_table_get_float(objtab,colname,jm,
                                                            &null));
                    break;
                case CPL_TYPE_DOUBLE:
                    cpl_table_set_double(mstds,colname,k,
                                         cpl_table_get_double(objtab,colname,
                                                              jm,&null));
                    break;
                default:
                    cpl_table_set_float(mstds,colname,k,
                                        cpl_table_get_float(objtab,colname,jm,
                                                            &null));
                    break;
                }
                colname = (char *)cpl_table_get_column_name(NULL);
            }
            cpl_table_select_row(mstds,k);
        }
    }

    /* Now extract the selected rows into the output table */

    *outtab = cpl_table_extract_selected(mstds);
    cpl_table_delete(mstds);

    /* Tidy up */

    freespace(matches);
    freespace(xoffs);
    freespace(yoffs);
    GOOD_STATUS
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
    
static cpl_table *vircam_mkmstd_table(cpl_table *objtab, cpl_table *stdstab) {
    cpl_table *mstds;
    char *colname;

    /* Copy the input standards table */

    mstds = cpl_table_duplicate(stdstab);

    /* Loop throught the object table columns and copy all of them over, 
       except for the RA and DEC tables */

    colname = (char *)cpl_table_get_column_name(objtab);
    while (colname != NULL) {
        if (strcmp(colname,"RA") && strcmp(colname,"DEC"))
            cpl_table_new_column(mstds,colname,
                                 cpl_table_get_column_type(objtab,colname));
        colname = (char *)cpl_table_get_column_name(NULL);
    }
    cpl_table_unselect_all(mstds);

    /* Get out of here */

    return(mstds);
}
        
/*

$Log: vircam_match.c,v $
Revision 1.13  2007/10/19 09:25:10  jim
Fixed problems with missing includes

Revision 1.12  2007/04/23 13:02:02  jim
Added static routine to make the matched standards table

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

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

Revision 1.9  2007/02/25 06:34:20  jim
Plugged memory leak

Revision 1.8  2006/07/03 09:33:18  jim
Fixed a few things to keep the compiler happy

Revision 1.7  2006/05/26 15:04:07  jim
Fixed sign error

Revision 1.6  2006/05/18 12:33:38  jim
Fixed bugs in the definition of the matched standard catalogue

Revision 1.5  2006/05/17 12:06:56  jim
Modified to take new definition of matched standards catalogue into account

Revision 1.4  2006/03/23 21:18:51  jim
Minor changes mainly to comment headers

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

Revision 1.2  2006/02/22 14:10:21  jim
Fixed omission in docs

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


*/

---