/*
* assoc.c
*
* Associate catalogues.
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
* This file part of: SExtractor
*
* Copyright: (C) 1997-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
*
* License: GNU General Public License
*
* SExtractor 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 3 of the License, or
* (at your option) any later version.
* SExtractor 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 SExtractor. If not, see .
*
* Last modified: 12/01/2011
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include
#include
#include
#include
#include "define.h"
#include "globals.h"
#include "prefs.h"
#include "assoc.h"
#include "fitswcs.h"
/********************************* comp_assoc ********************************/
/*
Comparison function for sort_assoc().
*/
int comp_assoc(const void *i1, const void *i2)
{
double *f1,*f2;
f1 = (double *)i1 + 1;
f2 = (double *)i2 + 1;
if (*f1<*f2)
return -1;
else return (*f1==*f2)?0:1;
}
/********************************* sort_assoc ********************************/
/*
Make the presentation histogram, order the assoc-list and build the hash-table.
*/
void sort_assoc(picstruct *field, assocstruct *assoc)
{
int comp_assoc(const void *i1, const void *i2);
double *list, rad;
int i,j, step,nobj, *hash;
step = assoc->ncol;
nobj = assoc->nobj;
list = assoc->list;
qsort(assoc->list, assoc->nobj, step*sizeof(double), comp_assoc);
/* Build the hash table that contains the first object in the sorted list */
/* which may be in reach from the current scanline */
QMALLOC(assoc->hash, int, field->height);
list = assoc->list+1; /* This is where the 1st y coordinate is stored */
hash = assoc->hash;
rad = assoc->radius;
for (i=0, j=0; iheight; i++)
{
/*-- For safety, we keep a 1-pixel margin */
for (;ji) ? -1 : j;
}
return;
}
/********************************* load_assoc ********************************/
/*
Read an assoc-list, and returns a pointer to the new assoc struct (or NULL if
no list was found).
*/
assocstruct *load_assoc(char *filename, wcsstruct *wcs)
{
assocstruct *assoc;
FILE *file;
double *list, val;
char str[MAXCHARL], str2[MAXCHARL], *sstr;
int *data,
i,ispoon,j,k,l, ncol, ndata, nlist, size,spoonsize,
xindex,yindex,mindex;
if (!(file = fopen(filename, "r")))
return NULL;
assoc = NULL; /* To avoid gcc -Wall warnings */
list = NULL; /* To avoid gcc -Wall warnings */
data = NULL; /* To avoid gcc -Wall warnings */
ispoon = ncol = ndata = nlist = size = spoonsize = xindex = yindex
= mindex = 0;
NFPRINTF(OUTPUT, "Reading ASSOC input-list...");
for (i=0; fgets(str, MAXCHARL, file);)
{
/*-- Examine current input line (discard empty and comment lines) */
if (!*str || strchr("#\t\n",*str))
continue;
if (!i)
{
strcpy(str2, str);
/*---- Let's count the number of columns in the first line */
for (ncol=0; strtok(ncol?NULL:str2, " \t\v\n\r\f"); ncol++);
if (!ncol)
error(EXIT_FAILURE, "*Error*: empty line in ", filename);
/*---- Build a look-up table containing the ordering of column data */
QCALLOC(data, int, ncol);
k = 1;
for (j=0; j= ncol)
error(EXIT_FAILURE, "*Error*: ASSOC_PARAMS #1 exceeds the number of ",
"fields in the ASSOC file");
if ((yindex = prefs.assoc_param[1]-1) >= ncol)
error(EXIT_FAILURE, "*Error*: ASSOC_PARAMS #2 exceeds the number of ",
"fields in the ASSOC file");
if (prefs.nassoc_param>2)
{
if ((mindex = prefs.assoc_param[2]-1) >= ncol)
error(EXIT_FAILURE,"*Error*: ASSOC_PARAMS #3 exceeds the number of ",
"fields in the ASSOC file");
}
else
{
mindex = -1;
prefs.assoc_type = ASSOC_FIRST;
}
nlist = ndata+3;
/*---- Allocate memory for the ASSOC struct and the filtered list */
QMALLOC(assoc, assocstruct, 1);
ispoon = ASSOC_BUFINC/(nlist*sizeof(double));
spoonsize = ispoon*nlist;
QMALLOC(assoc->list, double, size = spoonsize);
list = assoc->list;
}
else if (!(i%ispoon))
{
QREALLOC(assoc->list, double, size += spoonsize);
list = assoc->list + i*nlist;
}
if (!(++i%1000))
{
sprintf(str2, "Reading input list... (%d objects)", i);
NFPRINTF(OUTPUT, str2);
}
/*-- Read the data normally */
*(list+2) = 0.0;
for (sstr = str, j=0; jlist, double, i*nlist);
assoc->nobj = i;
assoc->radius = prefs.assoc_radius;
assoc->ndata = ndata;
assoc->ncol = nlist;
return assoc;
}
/********************************* init_assoc ********************************/
/*
Initialize the association procedure.
*/
void init_assoc(picstruct *field)
{
assocstruct *assoc;
/* Load the assoc-list */
if (!(assoc = field->assoc = load_assoc(prefs.assoc_name,
prefs.assoccoord_type==ASSOCCOORD_WORLD?
field->wcs : NULL)))
error(EXIT_FAILURE, "*Error*: Assoc-list file not found: ",
prefs.assoc_name);
/* Sort the assoc-list by y coordinates, and build the hash table */
sort_assoc(field, assoc);
/* Where data go for the current output pattern*/
assoc->data = outobj2.assoc;
return;
}
/********************************** end_assoc ********************************/
/*
Free memory related to the assoc operations.
*/
void end_assoc(picstruct *field)
{
/* Free the assoc-list */
if (field->assoc)
{
free((field->assoc)->list);
free((field->assoc)->hash);
free(field->assoc);
}
return;
}
/********************************** do_assoc *********************************/
/*
Perform the association task for a source and return the number of IDs.
*/
int do_assoc(picstruct *field, double x, double y)
{
assocstruct *assoc;
double aver, dx,dy, dist, rad, rad2, comp, wparam,
*list, *input, *data;
int h, step, i, flag, iy, nobj;
assoc = field->assoc;
/* Need to initialize the array */
memset(assoc->data, 0, prefs.assoc_size*sizeof(double));
aver = 0.0;
if (prefs.assoc_type == ASSOC_MIN || prefs.assoc_type == ASSOC_NEAREST)
comp = BIG;
else
comp = -BIG;
iy = (int)(y+0.499999);
if (iy<0 || iy>=field->height)
return 0;
/* If there is already no candidate in hash table, no need to go further */
if ((h=assoc->hash[iy])<0)
return 0;
/* Now loop over possible candidates */
nobj = assoc->nobj;
step = assoc->ncol;
list = assoc->list + step*h;
rad = assoc->radius;
rad2 = rad*rad;
for (flag = 0; (h++data, input, assoc->ndata*sizeof(double));
return 1;
}
wparam = *(list+2);
data = assoc->data;
switch(prefs.assoc_type)
{
case ASSOC_NEAREST:
if (distndata*sizeof(double));
comp = dist;
}
break;
case ASSOC_MEAN:
aver += wparam;
for (i=assoc->ndata; i--;)
*(data++) += *(input++)*wparam;
break;
case ASSOC_MAGMEAN:
wparam = fabs(wparam)<99.0?DEXP(-0.4*wparam): 0.0;
aver += wparam;
for (i=assoc->ndata; i--;)
*(data++) += *(input++)*wparam;
break;
case ASSOC_SUM:
for (i=assoc->ndata; i--;)
*(data++) += *(input++);
break;
case ASSOC_MAGSUM:
for (i=assoc->ndata; i--;)
*(data++) += fabs(wparam=*(input++))<99.0? DEXP(-0.4*wparam):0.0;
break;
case ASSOC_MIN:
if (wparamndata*sizeof(double));
comp = wparam;
}
break;
case ASSOC_MAX:
if (wparam>comp)
{
memcpy(data, input, assoc->ndata*sizeof(double));
comp = wparam;
}
break;
default:
error(EXIT_FAILURE, "*Internal Error*: unknown ASSOC type in ",
"pixlearn()");
}
}
}
/* No candidate found? exit! */
if (!flag)
return 0;
/* Terminate the computation of the mean */
if (prefs.assoc_type == ASSOC_MEAN || prefs.assoc_type == ASSOC_MAGMEAN)
{
if (aver<1e-30)
return 0;
data = assoc->data;
for (i=assoc->ndata; i--;)
*(data++) /= aver;
}
if (prefs.assoc_type == ASSOC_MAGSUM)
{
data = assoc->data;
for (i=assoc->ndata; i--; data++)
*data = *data>0.0? -2.5*log10(*data):99.0;
}
return flag;
}