Skip to content
GitLab
Commit deaee39b authored by Emmanuel Bertin's avatar Emmanuel Bertin
Browse files

Fixed PCi_j WCS keyword parsing issue (thanks to O.Fors for reporting this issue). 

Fixed issue with magnitudes when used as PSF model parameters (thanks to A.Donnarumma for reporting this issue).
Added ID_PARENT catalog parameter to identify the common parent of deblended sources (thanks to R.McMahon for the suggestion).
Introduced regularization of Sersic model ellipticities based on AMALGAm prescriptions for the Great3 challenge (see Mandelbaum et al. 2015).
Fixed error in ecliptic coordinate computations (thanks to H.Bouy for reporting this issue).
Added support for differential geometry maps, with new DGEO_TYPE and DGEO_IMAGE configuration keywords and new DGEOX/Y_IMAGE and DGEOX/YWIN_IMAGE measurement parameters.
Added differential geometry map correction to model-fitting.
Added VIGNET_DGEOX and VIGNET_DGEOY differential geometry map cutouts in catalogs (for use in PSFEx).
Added support for PSF model dependencies in ASSOC vectors.
Fixed issue with single-line background meshes.
Fixed positional offsets of large galaxy models on MODELS and -MODELS check-images.
Pushed version number to 2.23.1.
parent eedbf805
Expand all Show whitespace changes
Inline Side-by-side
src/preflist.h
View file @ deaee39b
......@@ -7,7 +7,7 @@
*
* This file part of: SExtractor
*
* Copyright: (C) 1993-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 1993-2015 Emmanuel Bertin -- IAP/CNRS/UPMC
*
* License: GNU General Public License
*
......@@ -22,7 +22,7 @@
* You should have received a copy of the GNU General Public License
* along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 30/10/2012
* Last modified: 12/01/2015
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -100,6 +100,9 @@
{""}, 1, 2, &prefs.ndthresh},
{"DETECT_TYPE", P_KEY, &prefs.detect_type, 0,0, 0.0,0.0,
{"CCD","PHOTO",""}},
{"DGEO_IMAGE", P_STRING, prefs.dgeoimage_name},
{"DGEO_TYPE", P_KEY, &prefs.dgeo_type, 0,0, 0.0,0.0,
{"NONE","PIXEL", ""}},
{"FILTER", P_BOOL, &prefs.filter_flag},
{"FILTER_NAME", P_STRING, prefs.filter_name},
{"FILTER_THRESH", P_FLOATLIST, prefs.filter_thresh, 0,0,-BIG,BIG,
......@@ -225,6 +228,11 @@ char *default_prefs[] =
"*FLAG_TYPE OR # flag pixel combination: OR, AND, MIN, MAX",
"* # or MOST",
"*",
"*#----------------------- Differential Geometry Map ---------------------------",
"*",
"*DGEO_TYPE NONE # Differential geometry map type: NONE or PIXEL",
"*DGEO_IMAGE dgeo.fits # Filename for input differential geometry image",
"*",
"#------------------------------ Photometry -----------------------------------",
" ",
"PHOT_APERTURES 5 # MAG_APER aperture diameter(s) in pixels",
......
src/prefs.h
View file @ deaee39b
......@@ -7,7 +7,7 @@
*
* This file part of: SExtractor
*
* Copyright: (C) 1993-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 1993-2015 Emmanuel Bertin -- IAP/CNRS/UPMC
*
* License: GNU General Public License
*
......@@ -22,7 +22,7 @@
* You should have received a copy of the GNU General Public License
* along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 02/11/2011
* Last modified: 02/12/2015
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -98,6 +98,9 @@ typedef struct
int weightgain_flag; /* weight gain? */
int wscale_flag[2]; /* Weight rescaling */
int nwscale_flag; /* nb of params */
/*----- differential geometry map */
char dgeoimage_name[MAXCHAR]; /* diff.geo. filename */
enum {DGEO_NONE, DGEO_PIXELMAP} dgeo_type; /* diff.geo. scheme */
/*----- photometry */
enum {CAT_NONE, ASCII, ASCII_HEAD, ASCII_SKYCAT, ASCII_VO,
FITS_LDAC, FITS_TPX, FITS_10}
......@@ -178,6 +181,8 @@ typedef struct
char retina_name[MAXCHAR]; /* retina filename */
int vignetsize[2]; /* vignet size */
int vigshiftsize[2]; /* shift-vignet size */
int vignet_dgeoxsize[2]; /* dgeo-x vignet size */
int vignet_dgeoysize[2]; /* dgeo-y vignet size */
int cleanmargin; /* CLEANing margin */
char som_name[MAXCHAR]; /* SOM filename */
int somfit_flag; /* ASSOCiation flag */
......
src/profit.c
View file @ deaee39b
This diff is collapsed.
src/profit.h
View file @ deaee39b
......@@ -7,7 +7,7 @@
*
* This file part of: SExtractor
*
* Copyright: (C) 2006-2013 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 2006-2015 Emmanuel Bertin -- IAP/CNRS/UPMC
*
* License: GNU General Public License
*
......@@ -22,7 +22,7 @@
* You should have received a copy of the GNU General Public License
* along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 13/02/2013
* Last modified: 17/02/2015
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -73,9 +73,9 @@
#define PROFIT_MAXR2MAX 1e6 /* Maximum r2_max for truncating profiles */
#define PROFIT_DYNPARAM 10.0 /* Dynamic compression param. in sigma units */
#define PROFIT_SMOOTHR 4.0 /* Profile smoothing radius (pixels) */
#define PROFIT_MAXMODSIZE 512 /* Maximum size allowed for the model raster */
#define PROFIT_MAXMODSIZE 8192 /* Maximum size allowed for the model raster */
#define PROFIT_MAXSMODSIZE 64 /* Number of model planes */
#define PROFIT_MAXOBJSIZE 512 /* Maximum size allowed for the object raster */
#define PROFIT_MAXOBJSIZE 8192 /* Maximum size allowed for the object raster */
#define PROFIT_BARXFADE 0.1 /* Fract. of bar length crossfaded with arms */
#define PROFIT_MAXEXTRA 2 /* Max. nb of extra free params of profiles */
#define INTERP_MAXKERNELWIDTH 8 /* Max. range of kernel (pixels) */
......@@ -182,6 +182,8 @@ typedef struct
PIXTYPE *lmodpix2; /* 2nd low resolution pixmap of the model */
PIXTYPE *objpix; /* Copy of object pixmap */
PIXTYPE *objweight; /* Copy of object weight-map */
PIXTYPE *dgeopix[2]; /* Copy of object differential geometry maps */
int dgeoflag; /* Set if diff. geometry maps are provided */
int objnaxisn[2]; /* Dimensions along each axis */
int nobjpix; /* Total number of "final" pixels */
int ix, iy; /* Integer coordinates of object pixmap */
......@@ -227,7 +229,7 @@ float *profit_compresi(profitstruct *profit, float dynparam,
int profit_boundtounbound(profitstruct *profit,
float *param, double *dparam, int index),
profit_copyobjpix(profitstruct *profit, picstruct *field,
picstruct *wfield),
picstruct *wfield, picstruct *dgeofield),
profit_covarunboundtobound(profitstruct *profit,
double *dparam, float *param),
profit_minimize(profitstruct *profit, int niter),
......@@ -249,8 +251,8 @@ void profit_dfit(profitstruct *profit, profitstruct *dprofit,
prof_end(profstruct *prof),
profit_addparam(profitstruct *profit, paramenum paramindex,
float **param),
profit_fit(profitstruct *profit,
picstruct *field, picstruct *wfield,
profit_fit(profitstruct *profit, picstruct *field,
picstruct *wfield, picstruct *dgeofield,
objstruct *obj, obj2struct *obj2),
profit_convmoments(profitstruct *profit, obj2struct *obj2),
profit_convolve(profitstruct *profit, float *modpix),
......
src/psf.c
View file @ deaee39b
......@@ -7,7 +7,7 @@
*
* This file part of: SExtractor
*
* Copyright: (C) 1998-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 1998-2015 Emmanuel Bertin -- IAP/CNRS/UPMC
*
* License: GNU General Public License
*
......@@ -22,7 +22,7 @@
* You should have received a copy of the GNU General Public License
* along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 12/07/2012
* Last modified: 15/12/2015
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -37,6 +37,7 @@
#include "define.h"
#include "globals.h"
#include "key.h"
#include "prefs.h"
#include "fits/fitscat.h"
#include "check.h"
......@@ -89,6 +90,7 @@ void psf_end(psfstruct *psf, psfitstruct *psfit)
for (d=0; d<ndim; d++)
free(psf->contextname[d]);
free(psf->context);
free(psf->contextindex);
free(psf->contextname);
free(psf->contextoffset);
free(psf->contextscale);
......@@ -118,15 +120,17 @@ Read the PSF data from a FITS file.
*/
psfstruct *psf_load(char *filename, int ext)
{
extern tabstruct *objtab;
static objstruct saveobj;
static obj2struct saveobj2;
psfstruct *psf;
catstruct *cat;
tabstruct *tab;
keystruct *key;
char *head, *ci,*co;
char rkeyname[80],
*head, *ci,*co, *pstr, *pstrbuf;
int deg[POLY_MAXDIM], group[POLY_MAXDIM], ndim, ngroup,
e,i,k;
d,e,i,k,p,n, keyn;
/* Open the cat (well it is not a "cat", but simply a FITS file */
if (!(cat = read_cat(filename)))
......@@ -152,8 +156,7 @@ psfstruct *psf_load(char *filename, int ext)
head = tab->headbuf;
/*-- Dimension of the polynomial */
if (fitsread(head, "POLNAXIS", &ndim, H_INT,T_LONG) == RETURN_OK
&& ndim)
if (fitsread(head, "POLNAXIS", &ndim, H_INT,T_LONG) == RETURN_OK && ndim)
{
/*-- So we have a polynomial description of the PSF variations */
if (ndim > POLY_MAXDIM)
......@@ -165,6 +168,7 @@ psfstruct *psf_load(char *filename, int ext)
QMALLOC(psf->contextname, char *, ndim);
QMALLOC(psf->context, double *, ndim);
QMALLOC(psf->contextindex, int, ndim);
QMALLOC(psf->contexttyp, t_type, ndim);
QMALLOC(psf->contextoffset, double, ndim);
QMALLOC(psf->contextscale, double, ndim);
......@@ -194,18 +198,42 @@ psfstruct *psf_load(char *filename, int ext)
else
/*------ The context element is a dynamic object parameter */
{
if ((k = findkey(psf->contextname[i], (char *)objkey,
strncpy(rkeyname, psf->contextname[i],40);
strtok_r(rkeyname, "([{}])", &pstrbuf);
if ((k = findkey(rkeyname, (char *)objkey,
sizeof(keystruct)))==RETURN_ERROR)
{
sprintf(gstr, "*Error*: %s CONTEXT parameter in %s unknown",
psf->contextname[i], psf->name);
rkeyname, psf->name);
error(EXIT_FAILURE, gstr, "");
}
keyn = (pstr = strtok_r(NULL,"([{}])", &pstrbuf))? atoi(pstr) - 1 : 0;
key = objkey+k;
psf->context[i] = key->ptr;
if (key->naxis)
{
n = 1;
for (d=0; d<key->naxis; d++)
n *= key->naxisn[d];
if (keyn >= n)
{
/*---------- Increase vector size: we restrict ourselves to vectors (1D) */
key->nbytes = key->naxisn[0] = keyn+1;
key->naxis = 1;
changecatparamarrays(rkeyname, key->naxisn, key->naxis);
}
}
else
keyn = -1;
psf->context[i] = (double *)((char *)key->ptr);
psf->contextindex[i] = keyn;
psf->contexttyp[i] = key->ttype;
/*------ Declare the parameter "active" to trigger computation by SExtractor */
*((char *)key->ptr) = (char)'\1';
/*------ Add the key to the catalog */
if (add_key(key, objtab, 0) != RETURN_ERROR)
thecat.nparam++;
if (!cistrcmp("MAG", rkeyname, FIND_NOSTRICT))
psf->mag_flag = 1;
}
/*---- Scaling of the context parameter */
sprintf(gstr, "POLZERO%1d", i+1);
......@@ -1046,7 +1074,10 @@ void psf_build(psfstruct *psf)
ndim = psf->poly->ndim;
for (i=0; i<ndim; i++)
{
ttypeconv(psf->context[i], &pos[i], psf->contexttyp[i],T_DOUBLE);
ttypeconv(psf->contextindex[i]<0? (char *)psf->context[i]
: *((char **)psf->context[i])
+ psf->contextindex[i]*t_size[psf->contexttyp[i]],
&pos[i], psf->contexttyp[i],T_DOUBLE);
pos[i] = (pos[i] - psf->contextoffset[i]) / psf->contextscale[i];
}
poly_func(psf->poly, pos);
......
src/psf.h
View file @ deaee39b
......@@ -7,7 +7,7 @@
*
* This file part of: SExtractor
*
* Copyright: (C) 1998-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 1998-2014 Emmanuel Bertin -- IAP/CNRS/UPMC
*
* License: GNU General Public License
*
......@@ -22,7 +22,7 @@
* You should have received a copy of the GNU General Public License
* along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 13/06/2012
* Last modified: 15/12/2015
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -81,6 +81,7 @@ typedef struct psf
float *maskcomp; /* Complete pix. data (PSF components) */
float *maskloc; /* Local PSF */
double **context; /* Contexts */
int *contextindex; /* Context index (for arrays) */
t_type *contexttyp; /* Context types */
char **contextname; /* Array of context key-names */
double *contextoffset; /* Offset to apply to context data */
......@@ -90,6 +91,7 @@ typedef struct psf
double fwhm; /* Typical PSF FWHM */
float pixstep; /* PSF sampling step */
int build_flag; /* Set if the current PSF has been computed */
int mag_flag; /* Set if PSF contexts include magnitudes */
} psfstruct;
typedef struct
......
src/readimage.c
View file @ deaee39b
......@@ -7,7 +7,7 @@
*
* This file part of: SExtractor
*
* Copyright: (C) 1993-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 1993-2015 Emmanuel Bertin -- IAP/CNRS/UPMC
*
* License: GNU General Public License
*
......@@ -22,7 +22,7 @@
* You should have received a copy of the GNU General Public License
* along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 26/06/2012
* Last modified: 14/01/2015
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -58,11 +58,12 @@ void *loadstrip(picstruct *field, picstruct *wfield)
{
tabstruct *tab;
checkstruct *check;
int y, w, flags, interpflag;
int y, w, flags, interpflag, npixtot;
PIXTYPE *data, *wdata, *rmsdata;
tab = field->tab;
w = field->width;
npixtot = field->width*field->height;
flags = field->flags;
interpflag = (wfield && wfield->interp_flag);
wdata = NULL; /* To avoid gcc -Wall warnings */
......@@ -74,7 +75,7 @@ void *loadstrip(picstruct *field, picstruct *wfield)
nbpix = w*field->stripheight;
if (flags ^ FLAG_FIELD)
if (!(flags & (FLAG_FIELD|DGEO_FIELD)))
{
/*---- Allocate space for the frame-buffer */
if (!(field->strip=(PIXTYPE *)malloc(field->stripheight*field->width
......@@ -138,14 +139,36 @@ void *loadstrip(picstruct *field, picstruct *wfield)
}
}
}
else
else if (flags & FLAG_FIELD)
{
/*---- flag map */
if (!(field->fstrip=(FLAGTYPE *)malloc(field->stripheight*field->width
*sizeof(FLAGTYPE))))
error(EXIT_FAILURE,"Not enough memory for the flag buffer of ",
field->rfilename);
read_ibody(field->tab, field->fstrip, nbpix);
}
else
{
/*---- differential geometry map */
if (!(field->dgeostrip[0]=(PIXTYPE *)malloc(field->stripheight
* field->width * sizeof(PIXTYPE)))
|| !(field->dgeostrip[1]=(PIXTYPE *)malloc(field->stripheight
* field->width * sizeof(PIXTYPE))))
error(EXIT_FAILURE,
"Not enough memory for differential geometry buffers of ",
field->rfilename);
/*---- Read first data plane (x shift) */
read_body(tab, field->dgeostrip[0], nbpix);
/*---- Move to second data plane */
QFSEEK(tab->cat->file,tab->bodypos + npixtot*tab->bytepix, SEEK_SET,
field->rfilename);
/*---- Read second data plane (y shift) */
read_body(tab, field->dgeostrip[1], nbpix);
/*---- Move back to first data plane */
QFSEEK(tab->cat->file,tab->bodypos + nbpix*tab->bytepix, SEEK_SET,
field->rfilename);
}
field->ymax = field->stripheight;
if (field->ymax < field->height)
......@@ -154,7 +177,7 @@ void *loadstrip(picstruct *field, picstruct *wfield)
else
{
/*- other strips */
if (flags ^ FLAG_FIELD)
if (!(flags & (FLAG_FIELD|DGEO_FIELD)))
{
data = field->strip + field->stripylim*w;
/*---- We assume weight data have been read just before */
......@@ -216,17 +239,34 @@ void *loadstrip(picstruct *field, picstruct *wfield)
}
}
}
else
else if (flags & FLAG_FIELD)
read_ibody(tab, field->fstrip + field->stripylim*w, w);
else
{
/*---- differential geometry map */
/*---- Read first data plane (x shift) */
read_body(tab, field->dgeostrip[0] + field->stripylim*w, w);
/*---- Move to second data plane */
QFSEEK(tab->cat->file,tab->bodypos
+ (field->ymax * w + npixtot) * tab->bytepix, SEEK_SET,
field->rfilename);
/*---- Read second data plane (y shift) */
read_body(tab, field->dgeostrip[1] + field->stripylim*w, w);
/*---- Move back to first data plane */
QFSEEK(tab->cat->file, tab->bodypos
+ (field->ymax + 1) * w * tab->bytepix, SEEK_SET, field->rfilename);
}
field->stripylim = (++field->ymin)%field->stripheight;
if ((++field->ymax)<field->height)
field->stripysclim = (++field->stripysclim)%field->stripheight;
}
return (flags ^ FLAG_FIELD)?
return !(flags & (FLAG_FIELD|DGEO_FIELD))?
(void *)(field->strip + field->stripy*w)
: (void *)(field->fstrip + field->stripy*w);
: ((flags & FLAG_FIELD)?
(void *)(field->fstrip + field->stripy*w)
: (void *)(field->dgeostrip[0] + field->stripy*w));
}
......
src/refine.c
View file @ deaee39b
......@@ -7,7 +7,7 @@
*
* This file part of: SExtractor
*
* Copyright: (C) 1993-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 1993-2016 Emmanuel Bertin -- IAP/CNRS/UPMC
*
* License: GNU General Public License
*
......@@ -22,7 +22,7 @@
* You should have received a copy of the GNU General Public License
* along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 11/10/2010
* Last modified: 13/01/2016
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -166,9 +166,10 @@ int parcelout(objliststruct *objlistin, objliststruct *objlistout)
if (ok[k+1+xn*j] && obj[j].fdflux-obj[j].dthresh*obj[j].fdnpix
> value0)
{
objlist[k+1].obj[j].flag |= OBJ_MERGED /* Merge flag on */
obj[j].flag |= OBJ_MERGED /* Merge flag on */
| ((OBJ_ISO_PB|OBJ_APERT_PB|OBJ_OVERFLOW)
&debobjlist2.obj[0].flag);
obj[j].id_parent = debobjlist2.obj[0].id_parent;
if ((out = addobj(j, &objlist[k+1], &debobjlist2))
== RETURN_FATAL_ERROR)
goto exit_parcelout;
......
src/scan.c
View file @ deaee39b
......@@ -7,7 +7,7 @@
*
* This file part of: SExtractor
*
* Copyright: (C) 1993-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 1993-2016 Emmanuel Bertin -- IAP/CNRS/UPMC
*
* License: GNU General Public License
*
......@@ -22,7 +22,7 @@
* You should have received a copy of the GNU General Public License
* along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 21/12/2011
* Last modified: 13/01/2016
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -49,7 +49,7 @@
/****************************** scanimage ************************************
PROTO void scanimage(picstruct *field, picstruct *dfield, picstruct *ffield,
picstruct *wfield, picstruct *dwfield)
picstruct *wfield, picstruct *dwfield, picstruct *dgeofield)
PURPOSE Scan of the large pixmap(s). Main loop and heart of the program.
INPUT Measurement field pointer,
Detection field pointer,
......@@ -59,10 +59,11 @@ INPUT Measurement field pointer,
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 21/12/2011
VERSION 14/01/2015
***/
void scanimage(picstruct *field, picstruct *dfield, picstruct **pffield,
int nffield, picstruct *wfield, picstruct *dwfield)
int nffield, picstruct *wfield, picstruct *dwfield,
picstruct *dgeofield)
{
static infostruct curpixinfo, *info, *store,
......@@ -82,7 +83,7 @@ void scanimage(picstruct *field, picstruct *dfield, picstruct **pffield,
PIXTYPE thresh, relthresh, cdnewsymbol, cdwthresh,wthresh,
*scan,*dscan,*cdscan,*dwscan,*dwscanp,*dwscann,
*cdwscan,*cdwscanp,*cdwscann,*wscand,
*scant, *wscan,*wscann,*wscanp;
*scant, *wscan,*wscann,*wscanp, *dgeoscanx, *dgeoscany;
FLAGTYPE *pfscan[MAXFLAG];
status cs, ps, *psstack;
int *start, *end, ymax;
......@@ -90,7 +91,7 @@ void scanimage(picstruct *field, picstruct *dfield, picstruct **pffield,
/* Avoid gcc -Wall warnings */
scan = dscan = cdscan = cdwscan = cdwscann = cdwscanp
= dwscan = dwscann = dwscanp
= wscan = wscann = wscanp = NULL;
= wscan = wscann = wscanp = dgeoscanx = dgeoscany = NULL;
victim = NULL; /* Avoid gcc -Wall warnings */
blankh = 0; /* Avoid gcc -Wall warnings */
/*----- Beginning of the main loop: Initialisations */
......@@ -102,7 +103,7 @@ void scanimage(picstruct *field, picstruct *dfield, picstruct **pffield,
/* cdwfield is the detection weight-field if available */
cdwfield = dwfield? dwfield:(prefs.dweight_flag?wfield:NULL);
cdwthresh = cdwfield ? cdwfield->weight_thresh : 0.0;
if (cdwthresh>BIG*WTHRESH_CONVFAC);
if (cdwthresh>BIG*WTHRESH_CONVFAC)
cdwthresh = BIG*WTHRESH_CONVFAC;
wthresh = wfield? wfield->weight_thresh : 0.0;
......@@ -144,6 +145,13 @@ void scanimage(picstruct *field, picstruct *dfield, picstruct **pffield,
dwfield->stripysclim = 0;
}
if (dgeofield) {
dgeofield->y = dgeofield->stripy = 0;
dgeofield->ymin = dgeofield->stripylim = 0;
dgeofield->stripysclim = 0;
}
/*Allocate memory for buffers */
stacksize = w+1;
QMALLOC(info, infostruct, stacksize);
......@@ -301,6 +309,15 @@ void scanimage(picstruct *field, picstruct *dfield, picstruct **pffield,
else
dscan = scan;
if (dgeofield) {
dgeoscanx = (dgeofield->stripy==dgeofield->stripysclim) ?
(PIXTYPE *)loadstrip(dgeofield, (picstruct *)NULL)
: dgeofield->dgeostrip[0]
+ dgeofield->stripy * dgeofield->width;
dgeoscany = dgeofield->dgeostrip[1]
+ dgeofield->stripy * dgeofield->width;
}
if (prefs.filter_flag)
{
filter(cfield, cdscan, cfield->y);
......@@ -443,6 +460,10 @@ void scanimage(picstruct *field, picstruct *dfield, picstruct **pffield,
PLISTPIX(pixt, dthresh) = thresh;
if (PLISTEXIST(var))
PLISTPIX(pixt, var) = wscan[xl];
if (PLISTEXIST(dgeo)) {
PLISTPIX(pixt, dgeox) = dgeoscanx[xl];
PLISTPIX(pixt, dgeoy) = dgeoscany[xl];
}
if (cs != OBJECT)
/*------------------------------- Start Segment -----------------------------*/
......@@ -515,8 +536,8 @@ void scanimage(picstruct *field, picstruct *dfield, picstruct **pffield,
{
if ((int)info[co].pixnb >= prefs.ext_minarea)
{
sortit(field, dfield, wfield, cdwfield, &info[co], &objlist,
cdwscan, wscan);
sortit(field, dfield, wfield, cdwfield, dgeofield,
&info[co], &objlist, cdwscan, wscan);
}
/* ------------------------------------ free the chain-list */
......@@ -613,6 +634,8 @@ void scanimage(picstruct *field, picstruct *dfield, picstruct **pffield,
wfield->stripy = (wfield->y=yl)%wfield->stripheight;
if (dwfield)
dwfield->stripy = (dwfield->y=yl)%dwfield->stripheight;
if (dgeofield)
dgeofield->stripy = (dgeofield->y=yl)%dgeofield->stripheight;
/*-- Remove objects close to the ymin limit if ymin is ready to increase */
if (cfield->stripy==cfield->stripysclim)
......@@ -645,7 +668,8 @@ void scanimage(picstruct *field, picstruct *dfield, picstruct **pffield,
"Objects: %8d detected / %8d sextracted\n\33[1A",
yl>h? h:yl, thecat.ndetect, thecat.ntotal);
ontotal = thecat.ntotal;
endobject(field, dfield, wfield, cdwfield, i, cleanobjlist);
endobject(field, dfield, wfield, cdwfield, dgeofield,
i, cleanobjlist);
subcleanobj(i);
cleanobj = cleanobjlist->obj+i; /* realloc in subcleanobj() */
}
......@@ -691,7 +715,7 @@ void scanimage(picstruct *field, picstruct *dfield, picstruct **pffield,
"Objects: %8d detected / %8d sextracted\n\33[1A",
h, thecat.ndetect, thecat.ntotal);
ontotal = thecat.ntotal;
endobject(field, dfield, wfield, cdwfield, 0, cleanobjlist);
endobject(field, dfield, wfield, cdwfield, dgeofield, 0, cleanobjlist);
subcleanobj(0);
}
......@@ -745,7 +769,8 @@ void update(infostruct *infoptr1, infostruct *infoptr2, pliststruct *pixel)
build the object structure.
*/
void sortit(picstruct *field, picstruct *dfield, picstruct *wfield,
picstruct *dwfield, infostruct *info, objliststruct *objlist,
picstruct *dwfield, picstruct *dgeofield,
infostruct *info, objliststruct *objlist,
PIXTYPE *cdwscan, PIXTYPE *wscan)
{
......@@ -754,6 +779,7 @@ void sortit(picstruct *field, picstruct *dfield, picstruct *wfield,
static objstruct obj;
objstruct *cobj;
pliststruct *pixel;
static int id_parent;
int i,j,n;
cfield = dfield? dfield: field;
......@@ -775,6 +801,7 @@ void sortit(picstruct *field, picstruct *dfield, picstruct *wfield,
obj.flag = info->flag;
obj.dthresh = objlist->dthresh;
obj.thresh = objlist->thresh;
obj.id_parent = ++id_parent;
preanalyse(0, objlist, ANALYSE_FAST);
......@@ -852,7 +879,8 @@ void sortit(picstruct *field, picstruct *dfield, picstruct *wfield,
}
}
endobject(field, dfield, wfield, dwfield, victim, cleanobjlist);
endobject(field, dfield, wfield, dwfield, dgeofield,
victim, cleanobjlist);
subcleanobj(victim);
}
......@@ -877,16 +905,16 @@ INPUT objlist number,
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP & Leiden & ESO)
VERSION 28/11/2003
VERSION 21/01/2015
***/
void preanalyse(int no, objliststruct *objlist, int analyse_type)
{
objstruct *obj = &objlist->obj[no];
pliststruct *pixel = objlist->plist, *pixt;
PIXTYPE peak, cpeak, val, cval, minthresh, thresht;
PIXTYPE peak, cpeak, val, cval, minthresh, thresht, dx, dy;
double thresh,thresh2, t1t2,darea,
mx,my, mx2,my2,mxy, rv, tv,
mx,my, mx2,my2,mxy, rv, tv, mdx, mdy,
xm,ym, xm2,ym2,xym,
temp,temp2, theta,pmx2,pmy2;
int x, y, xmin,xmax, ymin,ymax,area2, fdnpix, dnpix;
......@@ -944,20 +972,26 @@ void preanalyse(int no, objliststruct *objlist, int analyse_type)
if (analyse_type & ANALYSE_FULL)
{
mx = my = tv = 0.0;
mx = my = tv = mdx = mdy = 0.0;
mx2 = my2 = mxy = 0.0;
thresh2 = (thresh + peak)/2.0;
area2 = 0;
for (pixt=pixel+obj->firstpix; pixt>=pixel; pixt=pixel+PLIST(pixt,nextpix))
{
x = PLIST(pixt,x)-xmin; /* avoid roundoff errors on big images */
y = PLIST(pixt,y)-ymin; /* avoid roundoff errors on big images */
cval = PLISTPIX(pixt, cdvalue);
tv += (val = PLISTPIX(pixt, dvalue));
if (val>thresh)
dnpix++;
if (val > thresh2)
area2++;
x = PLIST(pixt,x)-xmin; /* avoid roundoff errors on big images */
y = PLIST(pixt,y)-ymin; /* avoid roundoff errors on big images */
if PLISTEXIST(dgeo) {
x -= (dx = PLISTPIX(pixt, dgeox));
y -= (dy = PLISTPIX(pixt, dgeoy));
mdx += cval * dx;
mdy += cval * dy;
}
mx += cval * x;
my += cval * y;
mx2 += cval * x*x;
......@@ -968,6 +1002,8 @@ void preanalyse(int no, objliststruct *objlist, int analyse_type)
/*----- compute object's properties */
xm = mx / rv; /* mean x */
ym = my / rv; /* mean y */
mdx /= rv; /* mean Delta-x */
mdy /= rv; /* mean Delta-y */
/*-- In case of blending, use previous barycenters */
if ((analyse_type&ANALYSE_ROBUST) && (obj->flag&OBJ_MERGED))
......@@ -1014,6 +1050,8 @@ void preanalyse(int no, objliststruct *objlist, int analyse_type)
obj->dflux = tv;
obj->mx = xm+xmin; /* add back xmin */
obj->my = ym+ymin; /* add back ymin */
obj->dmx = (float)mdx;
obj->dmy = (float)mdy;
obj->mx2 = xm2;
obj->my2 = ym2;
obj->mxy = xym;
......
src/types.h
View file @ deaee39b
......@@ -7,7 +7,7 @@
*
* This file part of: SExtractor
*
* Copyright: (C) 1993-2013 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 1993-2015 Emmanuel Bertin -- IAP/CNRS/UPMC
*
* License: GNU General Public License
*
......@@ -22,7 +22,7 @@
* You should have received a copy of the GNU General Public License
* along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 05/07/2013
* Last modified: 02/12/2015
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -92,6 +92,7 @@ typedef struct
{
/* ---- basic parameters */
int number; /* ID */
int id_parent; /* Parent ID */
int fdnpix; /* nb of extracted pix */
int dnpix; /* nb of pix above thresh */
int npix; /* "" in measured frame */
......@@ -107,6 +108,7 @@ typedef struct
/* ---- astrometric data */
int peakx,peaky; /* pos of brightest pix */
double mx, my; /* barycenter */
float dmx, dmy; /* dgeo barycenter correction */
double poserr_mx2, poserr_my2,
poserr_mxy; /* Error ellips moments */
/* ---- morphological data */
......@@ -130,7 +132,6 @@ typedef struct
float mthresh; /* max. threshold (ADU) */
int iso[NISO]; /* isophotal areas */
float fwhm; /* IMAGE FWHM */
} objstruct;
/* II: "BLIND" parameters */
......@@ -175,6 +176,7 @@ typedef struct
double pixscale2; /* Local pixel area */
double mamaposx,mamaposy; /* "MAMA" pos. in pixels */
float sposx,sposy; /* single precision pos. */
float pos_dgeox, pos_dgeoy; /* d.geom. correction to pos. */
float poserr_a, poserr_b,
poserr_theta; /* Error ellips parameters */
float poserr_cxx, poserr_cyy,
......@@ -220,9 +222,12 @@ typedef struct
int assoc_number; /* nb of ASSOCiated objects */
float *vignet; /* Pixel data */
float *vigshift; /* (Shifted) pixel data */
float *vignet_dgeox; /* Differential x geom. data */
float *vignet_dgeoy; /* Differential y geom. data */
/* Windowed measurements */
double winpos_x,winpos_y; /* Windowed barycenter */
float winpos_dgeox, winpos_dgeoy; /* d.geom. corr. to win. pos */
double winposerr_mx2, winposerr_my2,
winposerr_mxy; /* Error ellips moments */
float winposerr_a, winposerr_b,
......@@ -574,6 +579,7 @@ typedef struct pic
int yblank; /* y blanking limit (highest+1) */
PIXTYPE *strip; /* pointer to the image buffer */
FLAGTYPE *fstrip; /* pointer to the FLAG buffer */
PIXTYPE *dgeostrip[2]; /* pointers to diff. geometry buffers */
int stripheight; /* height of a strip (in lines) */
int stripmargin; /* number of lines in margin */
int stripstep; /* number of lines at each read */
......
src/wcs/poly.c
View file @ deaee39b
/*
* poly.c
*
* Manage polynomials.
* Polynomial functions.
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
* This file part of: AstrOmatic WCS library
* This file part of: AstrOmatic software
*
* Copyright: (C) 1998-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 1998-2012 IAP/CNRS/UPMC
*
* License: GNU General Public License
*
......@@ -23,7 +23,7 @@
* along with AstrOmatic software.
* If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 20/12/2011
* Last modified: 20/11/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -36,16 +36,17 @@
#include <stdlib.h>
#include <string.h>
#include "poly.h"
#ifdef HAVE_ATLAS
#include ATLAS_LAPACK_H
#endif
#ifdef HAVE_LAPACKE
#include LAPACKE_H
//#define MATSTORAGE_PACKED 1
#endif
#include "poly.h"
#define QCALLOC(ptr, typ, nel) \
{if (!(ptr = (typ *)calloc((size_t)(nel),sizeof(typ)))) \
qerror("Not enough memory for ", \
......@@ -56,6 +57,13 @@
qerror("Not enough memory for ", \
#ptr " (" #nel " elements) !");;}
#define QMEMCPY(ptrin, ptrout, typ, nel) \
{if (ptrin) \
{if (!(ptrout = (typ *)malloc((size_t)(nel)*sizeof(typ)))) \
qerror("Not enough memory for ", \
#ptrout " (" #nel " elements) !"); \
memcpy(ptrout, ptrin, (size_t)(nel)*sizeof(typ));};;}
/********************************* qerror ************************************/
/*
I hope it will never be used!
......@@ -149,8 +157,8 @@ PURPOSE Free a polynom structure and everything it contains.
INPUT polystruct pointer.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
VERSION 09/04/2000
AUTHOR E. Bertin (IAP)
VERSION 04/11/2008
***/
void poly_end(polystruct *poly)
{
......@@ -158,16 +166,63 @@ void poly_end(polystruct *poly)
{
free(poly->coeff);
free(poly->basis);
free(poly->orthobasis);
free(poly->degree);
free(poly->group);
free(poly->orthomat);
free(poly->deorthomat);
free(poly);
}
return;
}
/****** poly_copy *************************************************************
PROTO polystruct *poly_copy(polystruct *poly)
PURPOSE Copy a polynom structure and everything it contains.
INPUT polystruct pointer.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 04/11/2008
***/
polystruct *poly_copy(polystruct *poly)
{
polystruct *newpoly;
if (poly)
{
QMALLOC(newpoly, polystruct, 1);
*newpoly = *poly;
if (poly->ncoeff)
{
QMEMCPY(poly->coeff, newpoly->coeff, double, poly->ncoeff);
QMEMCPY(poly->basis, newpoly->basis, double, poly->ncoeff);
}
if (poly->ndim)
QMEMCPY(poly->group, newpoly->group, int, poly->ndim);
if (poly->ngroup)
QMEMCPY(poly->degree, newpoly->degree, int, poly->ngroup);
if (poly->orthomat)
{
QMEMCPY(poly->orthomat, newpoly->orthomat, double,
poly->ncoeff*poly->ncoeff);
QMEMCPY(poly->deorthomat, newpoly->deorthomat, double,
poly->ncoeff*poly->ncoeff);
QMEMCPY(poly->orthobasis, newpoly->orthobasis, double, poly->ncoeff);
}
return newpoly;
}
else
return NULL;
}
/****** poly_func ************************************************************
PROTO double poly_func(polystruct *poly, double *pos)
PURPOSE Evaluate a multidimensional polynom.
PURPOSE Evaluate a multidimensional polynomial.
INPUT polystruct pointer,
pointer to the 1D array of input vector data.
OUTPUT Polynom value.
......@@ -243,9 +298,88 @@ double poly_func(polystruct *poly, double *pos)
}
/****** poly_cfunc ************************************************************
PROTO double poly_cfunc(polystruct *poly, double *pos)
PURPOSE Evaluate a multidimensional Chebyshev polynomial.
INPUT polystruct pointer,
pointer to the 1D array of input vector data.
OUTPUT Polynom value.
NOTES Values of the basis functions are updated in poly->basis.
AUTHOR E. Bertin (IAP)
VERSION 29/01/2013
***/
double poly_cfunc(polystruct *poly, double *pos)
{
double pol[POLY_MAXDIM*(POLY_MAXDEGREE+1)],
*polt, *post, *basis, *coeff, xval;
long double val;
int expo[POLY_MAXDIM+1], gexpo[POLY_MAXDIM+1];
int *expot, *degree,*degreet, *group,*groupt, *gexpot,
d,d2,g,t, ndim;
/* Prepare the vectors and counters */
ndim = poly->ndim;
basis = poly->basis;
coeff = poly->coeff;
group = poly->group;
degree = poly->degree;
if (ndim)
{
for (groupt=group, expot=expo, post=pos, d=0; d<ndim; d++)
{
*(expot++) = 0;
polt = pol + d*(POLY_MAXDEGREE+1);
*(polt++) = 1.0;
*(polt++) = xval = *(post++);
for (d2 = degree[*(groupt++)]; --d2 > 0; polt++)
*polt = 2.0*xval**(polt-1) - *(polt-2);
}
for (gexpot=gexpo, degreet=degree, g=poly->ngroup; g--;)
*(gexpot++) = *(degreet++);
if (gexpo[*group])
gexpo[*group]--;
}
/* The constant term is handled separately */
val = *(coeff++);
*(basis++) = 1.0;
*expo = 1;
/* Compute the rest of the polynom */
for (t=poly->ncoeff; --t; )
{
polt = pol;
expot = expo;
/*-- xval contains the current product of the polynomials */
xval = 1.0;
for (d=ndim; d--; polt += POLY_MAXDEGREE+1)
xval *= polt[*(expot++)];
val += (*(basis++)=xval)**(coeff++);
/*-- A complex recursion between terms of the polynom speeds up computations */
/*-- Not too good for roundoff errors (prefer Horner's), but much easier for */
/*-- multivariate polynomials: this is why we use a long double accumulator */
expot = expo;
groupt = group;
for (d=0; d<ndim; d++, groupt++)
if (gexpo[*groupt]--)
{
++*(expot++);
break;
}
else
{
gexpo[*groupt] = *expot;
*(expot++) = 0;
}
}
return (double)val;
}
/****** poly_fit *************************************************************
PROTO double poly_fit(polystruct *poly, double *x, double *y, double *w,
int ndata, double *extbasis)
PROTO int poly_fit(polystruct *poly, double *x, double *y, double *w,
int ndata, double *extbasis, double regul)
PURPOSE Least-Square fit of a multidimensional polynom to weighted data.
INPUT polystruct pointer,
pointer to the (pseudo)2D array of inputs to basis functions,
......@@ -253,23 +387,24 @@ INPUT polystruct pointer,
pointer to the 1D array of data weights,
number of data points,
pointer to a (pseudo)2D array of computed basis function values.
Tikhonov regularization parameter (0 = no regularization).
OUTPUT Chi2 of the fit.
NOTES If different from NULL, extbasis can be provided to store the
values of the basis functions. If x==NULL and extbasis!=NULL, the
precomputed basis functions stored in extbasis are used (which saves
CPU). If w is NULL, all points are given identical weight.
AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
VERSION 08/03/2005
AUTHOR E. Bertin (IAP)
VERSION 20/11/2012
***/
void poly_fit(polystruct *poly, double *x, double *y, double *w, int ndata,
double *extbasis)
int poly_fit(polystruct *poly, double *x, double *y, double *w, int ndata,
double *extbasis, double regul)
{
void qerror(char *msg1, char *msg2);
double /*offset[POLY_MAXDIM],*/x2[POLY_MAXDIM],
*alpha,*alphat, *beta,*betat, *basis,*basis1,*basis2, *coeff,
*extbasist,*xt,
val,wval,yval;
int ncoeff, ndim, matsize,
int ncoeff, ndim, matsize, info,
d,i,j,n;
if (!x && !extbasis)
......@@ -331,8 +466,13 @@ void poly_fit(polystruct *poly, double *x, double *y, double *w, int ndata,
}
}
if (regul>POLY_TINY)
/*-- Simple Tikhonov regularization */
for (i=0; i<ncoeff; i++)
alpha[i*(ncoeff+1)] += regul;
/* Solve the system */
poly_solve(alpha,beta,ncoeff);
info = poly_solve(alpha,beta,ncoeff);
free(alpha);
......@@ -346,7 +486,7 @@ void poly_fit(polystruct *poly, double *x, double *y, double *w, int ndata,
*/
free(beta);
return;
return info;
}
......@@ -433,33 +573,32 @@ void poly_addcste(polystruct *poly, double *cste)
return;
}
/****** poly_solve ************************************************************
PROTO void poly_solve(double *a, double *b, int n)
PROTO int poly_solve(double *a, double *b, int n)
PURPOSE Solve a system of linear equations, using Cholesky decomposition.
INPUT Pointer to the (pseudo 2D) matrix of coefficients,
pointer to the 1D column vector,
matrix size.
OUTPUT -.
OUTPUT 0 if solution OK, !=0 otherwise (e.g., singular matrix).
NOTES -.
AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
VERSION 20/12/2011
AUTHOR E. Bertin (IAP)
VERSION 20/11/2012
***/
void poly_solve(double *a, double *b, int n)
int poly_solve(double *a, double *b, int n)
{
#if defined(HAVE_LAPACKE)
LAPACKE_dposv(LAPACK_COL_MAJOR, 'L', n, 1, a, n, b, n);
return LAPACKE_dposv(LAPACK_COL_MAJOR, 'L', n, 1, a, n, b, n);
#elif defined(HAVE_ATLAS)
clapack_dposv(CblasRowMajor, CblasUpper, n, 1, a, n, b, n);
return clapack_dposv(CblasRowMajor, CblasUpper, n, 1, a, n, b, n);
#else
cholsolve(a,b,n);
return cholsolve(a,b,n);
#endif
return;
}
/****** cholsolve *************************************************************
PROTO int cholsolve(double *a, double *b, int n)
PROTO void cholsolve(double *a, double *b, int n)
PURPOSE Solve a system of linear equations, using Cholesky decomposition.
INPUT Pointer to the (pseudo 2D) matrix of coefficients,
pointer to the 1D column vector,
......@@ -582,3 +721,181 @@ int *poly_powers(polystruct *poly)
return powers;
}
/****** poly_initortho ********************************************************
PROTO void poly_initortho(polystruct *poly, double *data, int ndata)
PURPOSE Compute orthonormalization and de-orthonormalization matrices for a
polynomial basis on a data set.
INPUT polystruct pointer,
pointer to the 1D array of input vector data,
number of data vectors.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 10/07/2012
***/
void poly_initortho(polystruct *poly, double *data, int ndata)
{
double *basis, *coeff, *invec,*invect0,*invect,*invect02,*invect2,
*rdiag, *deortho,
scale,s, dval;
int c,i,j,m,n, ndmc, ndim,ncoeff;
/* Prepare the vectors and counters */
ndim = poly->ndim;
ncoeff = poly->ncoeff;
basis = poly->basis;
coeff = poly->coeff;
/* Allocate memory for orthonormalization matrix and vector */
QCALLOC(poly->deorthomat, double, ncoeff*ncoeff);
QMALLOC(poly->orthobasis, double, poly->ncoeff);
QMALLOC(rdiag, double, ncoeff);
/* Do a QR decomposition of input vector set */
/* Vectors are stored as rows to speed up the Householder transformation */
n = ncoeff;
m = ndata;
invec = data;
for (c=0; c<ncoeff; c++)
{
ndmc = ndata - c;
scale = 0.0;
invect = invect0 = data + c*(ndata+1);
for (i=ndmc; i--; invect++)
scale = sqrt(scale*scale + *invect**invect);
if (scale > POLY_TINY)
{
if (*invect0 < 0.0)
scale = -scale;
invect = invect0;
for (i=ndmc; i--;)
*(invect++) /= scale;
*invect0 += 1.0;
invect02 = invect0 + ndata;
for (j=ncoeff-c; --j; invect02+=ndata)
{
s = 0.0;
invect = invect0;
invect2 = invect02;
for (i=ndmc; i--;)
s += *(invect++)**(invect2++);
s /= -*invect0;
invect = invect0;
invect2 = invect02;
for (i=ndmc; i--;)
*(invect2++) += s**(invect++);
}
}
rdiag[c] = -scale;
}
/* Convert to deorthonormalization matrix */
deortho = poly->deorthomat;
for (j=0; j<ncoeff; j++)
for (i=0; i<ncoeff; i++)
deortho[j*ncoeff+i] = i<j? data[j*ndata+i] : (i==j?rdiag[i] : 0.0);
free(rdiag);
/* Compute the "unorthonormalization" matrix */
QMEMCPY(poly->deorthomat, poly->orthomat, double, ncoeff*ncoeff);
#if defined(HAVE_LAPACKE)
LAPACKE_dtrtri(LAPACK_ROW_MAJOR, 'L', 'N', ncoeff,poly->orthomat,ncoeff);
#elif defined(HAVE_ATLAS)
clapack_dtrtri(CblasRowMajor, CblasLower, CblasNonUnit, ncoeff,
poly->orthomat, ncoeff);
#else
qerror("*Internal Error*: no routine available", " for triangular inverse");
#endif
/* Transpose orthonormalization matrix to speed up later use */
deortho = poly->deorthomat;
for (j=0; j<ncoeff; j++)
for (i=j; i<ncoeff; i++)
{
dval = deortho[j*ncoeff+i];
deortho[j*ncoeff+i] = deortho[i*ncoeff+j];
deortho[i*ncoeff+j] = dval;
}
return;
}
/****** poly_ortho ************************************************************
PROTO double *poly_ortho(polystruct *poly, double *datain, double *dataout)
PURPOSE Apply orthonormalization to the poly basis vector ("ket>").
INPUT polystruct pointer,
pointer to the input vector,
pointer to the output vector.
OUTPUT Pointer to poly->orthobasis, or poly->basis if no ortho. matrix exists.
NOTES The poly->basis vector must have been updated with poly_func() first.
AUTHOR E. Bertin (IAP)
VERSION 04/11/2008
***/
double *poly_ortho(polystruct *poly, double *datain, double *dataout)
{
double *omat,*basis,*obasis,
dval;
int i,j, ncoeff;
if (!poly->orthomat)
return datain;
ncoeff = poly->ncoeff;
/* Compute matrix product */
omat = poly->orthomat;
obasis = dataout;
for (j=ncoeff; j--;)
{
basis = datain;
dval = 0.0;
for (i=ncoeff; i--;)
dval += *(omat++)**(basis++);
*(obasis++) = dval;
}
return dataout;
}
/****** poly_deortho **********************************************************
PROTO void poly_deortho(polystruct *poly, double *datain, double *dataout)
PURPOSE Apply deorthonormalization to the poly basis component vector("<bra|").
INPUT polystruct pointer,
pointer to the input vector,
pointer to the output vector.
OUTPUT Pointer to poly->basis.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 04/11/2008
***/
double *poly_deortho(polystruct *poly, double *datain, double *dataout)
{
double *omat,*basis,*obasis,
dval;
int i,j, ncoeff;
if (!poly->deorthomat)
return datain;
ncoeff = poly->ncoeff;
/* Compute matrix product */
omat = poly->deorthomat;
basis = dataout;
for (j=ncoeff; j--;)
{
obasis = datain;
dval = 0.0;
for (i=ncoeff; i--;)
dval += *(omat++)**(obasis++);
*(basis++) = dval;
}
return dataout;
}
src/wcs/poly.h
View file @ deaee39b
......@@ -5,9 +5,9 @@
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
* This file part of: AstrOmatic WCS library
* This file part of: AstrOmatic software
*
* Copyright: (C) 1998-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 1998-2011 IAP/CNRS/UPMC
*
* License: GNU General Public License
*
......@@ -23,7 +23,7 @@
* along with AstrOmatic software.
* If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 10/10/2010
* Last modified: 20/11/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -33,7 +33,8 @@
/*--------------------------------- constants -------------------------------*/
#define POLY_MAXDIM 4 /* Max dimensionality of polynom */
#define POLY_MAXDEGREE 10 /* Max degree of the polynom */
#define POLY_MAXDEGREE 40 /* Max degree of the polynom */
#define POLY_TINY 1e-30 /* A tiny number */
/*---------------------------------- macros ---------------------------------*/
......@@ -42,27 +43,40 @@
typedef struct poly
{
double *basis; /* Current values of the basis functions */
double *orthobasis; /* Curr orthonormalized basis function values */
double *coeff; /* Polynom coefficients */
int ncoeff; /* Number of coefficients */
int *group; /* Groups */
int ndim; /* dimensionality of the polynom */
int *degree; /* Degree in each group */
int ngroup; /* Number of different groups */
double *orthomat; /* Orthonormalization matrix */
double *deorthomat; /* "Deorthonormalization" matrix */
} polystruct;
/*---------------------------------- protos --------------------------------*/
extern polystruct *poly_init(int *group,int ndim,int *degree,int ngroup);
extern polystruct *poly_copy(polystruct *poly),
*poly_init(int *group,int ndim,int *degree,int ngroup);
extern double poly_func(polystruct *poly, double *pos);
extern int cholsolve(double *a, double *b, int n),
*poly_powers(polystruct *poly);
extern void poly_addcste(polystruct *poly, double *cste),
poly_end(polystruct *poly),
poly_fit(polystruct *poly, double *x, double *y,
double *w, int ndata, double *extbasis),
double *w, int ndata, double *extbasis,
double regul),
*poly_powers(polystruct *poly),
poly_solve(double *a, double *b, int n);
extern void poly_addcste(polystruct *poly, double *cste),
poly_end(polystruct *poly);
extern double *poly_deortho(polystruct *poly, double *datain,
double *dataout),
*poly_ortho(polystruct *poly, double *datain,
double *dataout);
extern void poly_initortho(polystruct *poly, double *data,
int ndata);
#endif
src/wcscelsys.h
View file @ deaee39b
......@@ -7,7 +7,7 @@
*
* This file part of: AstrOmatic software
*
* Copyright: (C) 1998-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 1998-2015 Emmanuel Bertin -- IAP/CNRS/UPMC
*
* License: GNU General Public License
*
......@@ -23,7 +23,7 @@
* along with AstrOmatic software.
* If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 10/10/2010
* Last modified: 23/11/2015
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -41,11 +41,11 @@ char celsysname[][2][8] = { {"RA--", "DEC-"},
double celsysorig[][2] = { {0.0, 0.0},
{266.40499625, -28.93617242},
{0.0, 0.0},
{42.29235, 59.52315}},
{42.308333, 59.528333}},
celsyspole[][2] = { {0.0, 90.0},
{192.85948123, 27.12825120},
{273.85261111, 66.99111111},
{283.7514, 15.70480}},
{270.00000000, 66.560709},
{283.754167, 15.708889}},
/* Note: the code to handle the rotation sign is not yet implemented!!! */
celsyssign[] = { 1.0,
1.0,
......
src/winpos.c
View file @ deaee39b
......@@ -7,7 +7,7 @@
*
* This file part of: SExtractor
*
* Copyright: (C) 2005-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 2005-2015 Emmanuel Bertin -- IAP/CNRS/UPMC
*
* License: GNU General Public License
*
......@@ -22,7 +22,7 @@
* You should have received a copy of the GNU General Public License
* along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 03/05/2011
* Last modified: 04/03/2015
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -50,23 +50,25 @@ INPUT Picture structure pointer,
OUTPUT -.
NOTES obj->posx and obj->posy are taken as initial centroid guesses.
AUTHOR E. Bertin (IAP)
VERSION 03/05/2011
VERSION 04/03/2015
***/
void compute_winpos(picstruct *field, picstruct *wfield, objstruct *obj)
void compute_winpos(picstruct *field, picstruct *wfield,
picstruct *dgeofield, objstruct *obj)
{
float r2,invtwosig2, raper,raper2, rintlim,rintlim2,rextlim2,
dx,dx1,dy,dy2, sig, invngamma, pdbkg,
dx,dx1,dy,dy2, ddx, ddy, sig, invngamma, pdbkg,
offsetx,offsety,scalex,scaley,scale2, locarea;
double tv, norm, pix, var, backnoise2, invgain, locpix,
dxpos,dypos, err,err2, emx2,emy2,emxy,
dxpos,dypos, ddxpos, ddypos, err,err2, emx2,emy2,emxy,
esum, temp,temp2, mx2, my2,mxy,pmx2, theta, mx,my,
mx2ph, my2ph;
dmx, dmy, mx2ph, my2ph, cx2, cy2, cxy;
int i,x,y, x2,y2, xmin,xmax,ymin,ymax, sx,sy, w,h,
fymin,fymax, pflag,corrflag, gainflag, errflag,
momentflag;
long pos;
PIXTYPE *strip,*stript, *wstrip,*wstript,
PIXTYPE *strip,*stript, *wstrip, *wstript,
*dgeoxstrip, *dgeoystrip, *dgeoxstript, *dgeoystript,
wthresh = 0.0;
if (wfield)
......@@ -113,6 +115,8 @@ void compute_winpos(picstruct *field, picstruct *wfield, objstruct *obj)
/* Use isophotal centroid as a first guess */
mx = obj2->posx - 1.0;
my = obj2->posy - 1.0;
dmx = obj2->pos_dgeox;
dmy = obj2->pos_dgeoy;
for (i=0; i<WINPOS_NITERMAX; i++)
{
......@@ -120,8 +124,8 @@ void compute_winpos(picstruct *field, picstruct *wfield, objstruct *obj)
xmax = (int)(mx+raper+1.499999);
ymin = (int)(my-raper+0.499999);
ymax = (int)(my+raper+1.499999);
mx2ph = mx*2.0 + 0.49999;
my2ph = my*2.0 + 0.49999;
mx2ph = (mx+dmx)*2.0 + 0.49999; // Pixel position, not the physical one
my2ph = (my+dmy)*2.0 + 0.49999;
if (xmin < 0)
{
......@@ -145,23 +149,35 @@ void compute_winpos(picstruct *field, picstruct *wfield, objstruct *obj)
}
tv = esum = emxy = emx2 = emy2 = mx2 = my2 = mxy = 0.0;
dxpos = dypos = 0.0;
dxpos = dypos = ddxpos = ddypos = 0.0;
strip = field->strip;
wstrip = wstript = NULL; /* To avoid gcc -Wall warnings */
if (wfield)
wstrip = wfield->strip;
if (dgeofield) { // Differential geometry map is present
dgeoxstrip = dgeofield->dgeostrip[0];
dgeoystrip = dgeofield->dgeostrip[1];
}
for (y=ymin; y<ymax; y++)
{
stript = strip + (pos = (y%h)*w + xmin);
if (wfield)
wstript = wstrip + pos;
for (x=xmin; x<xmax; x++, stript++, wstript++)
if (dgeofield) {
dgeoxstript = dgeoxstrip + pos;
dgeoystript = dgeoystrip + pos;
}
for (x=xmin; x<xmax; x++)
{
dx = x - mx;
dy = y - my;
if (dgeofield) {
dx -= (ddx = *(dgeoxstript++));
dy -= (ddy = *(dgeoystript++));
}
if ((r2=dx*dx+dy*dy)<rextlim2)
{
if (WINPOS_OVERSAMP>1 && r2> rintlim2)
if (WINPOS_OVERSAMP > 1 && r2 > rintlim2)
{
dx += offsetx;
dy += offsety;
......@@ -204,10 +220,16 @@ void compute_winpos(picstruct *field, picstruct *wfield, objstruct *obj)
}
if (pflag)
pix = expf(pix*invngamma);
dx = x - mx;
dy = y - my;
locpix = locarea*pix;
tv += locpix;
dx = x - mx;
dy = y - my;
if (dgeofield) {
dx -= ddx;
dy -= ddy;
ddxpos += locpix*ddx;
ddypos += locpix*ddy;
}
dxpos += locpix*dx;
dypos += locpix*dy;
if (errflag)
......@@ -228,20 +250,20 @@ void compute_winpos(picstruct *field, picstruct *wfield, objstruct *obj)
emy2 += err2*(dy*dy+0.0833); /* Finite pixel size */
emxy += err2*dx*dy;
}
if (momentflag)
{
mx2 += locpix*dx*dx;
my2 += locpix*dy*dy;
mxy += locpix*dx*dy;
}
}
stript++;
if (wfield)
wstript++;
}
}
if (tv>0.0)
{
mx += (dxpos /= tv)*WINPOS_FAC;
my += (dypos /= tv)*WINPOS_FAC;
mx += (dxpos /= tv) * WINPOS_FAC;
my += (dypos /= tv) * WINPOS_FAC;
}
else
break;
......@@ -250,11 +272,17 @@ void compute_winpos(picstruct *field, picstruct *wfield, objstruct *obj)
if (dxpos*dxpos+dypos*dypos < WINPOS_STEPMIN*WINPOS_STEPMIN)
break;
}
mx2 = mx2/tv - dxpos*dxpos;
my2 = my2/tv - dypos*dypos;
mxy = mxy/tv - dxpos*dypos;
obj2->winpos_x = mx + 1.0; /* The dreaded 1.0 FITS offset */
obj2->winpos_y = my + 1.0; /* The dreaded 1.0 FITS offset */
if (dgeofield) {
obj2->winpos_dgeox = (ddxpos / tv) * WINPOS_FAC;
obj2->winpos_dgeoy = (ddypos / tv) * WINPOS_FAC;
}
obj2->winpos_niter = i+1;
/* WINdowed flux */
......@@ -265,10 +293,24 @@ void compute_winpos(picstruct *field, picstruct *wfield, objstruct *obj)
obj2->snr_win = esum>(1.0/BIG)? obj2->flux_win / obj2->fluxerr_win: BIG;
}
temp2=mx2*my2-mxy*mxy;
obj2->win_flag = (tv <= 0.0)*4 + (mx2 < 0.0 || my2 < 0.0)*2 + (temp2<0.0);
cx2 = obj->cxx;
cy2 = obj->cyy;
cxy = obj->cxy;
dx = obj2->winpos_x - obj2->posx;
dy = obj2->winpos_y - obj2->posy;
obj2->win_flag = 8 * (dx*dx > 1.0 && dy*dy > 1.0
&& (cx2*dx*dx + cy2*dy*dy + cxy*dx*dy) > 1.0)
// Check that the difference vector with isophotal
// centroid does not lie outside of the 1-sigma
// footprint ellipse
| 4 * (tv <= 0.0) // Check negative flux
| 2 * (mx2 < 0.0 || my2 < 0.0) // Negative 2nd order moments
| (temp2<0.0);
if (obj2->win_flag)
{
/*--- Negative values: revert to isophotal estimates */
obj2->winpos_x = obj2->posx;
obj2->winpos_y = obj2->posy;
if (FLAG(obj2.winposerr_mx2))
{
obj2->winposerr_mx2 = obj->poserr_mx2;
......
src/winpos.h
View file @ deaee39b
......@@ -7,7 +7,7 @@
*
* This file part of: SExtractor
*
* Copyright: (C) 2005-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 2005-2015 Emmanuel Bertin -- IAP/CNRS/UPMC
*
* License: GNU General Public License
*
......@@ -22,7 +22,7 @@
* You should have received a copy of the GNU General Public License
* along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 11/10/2010
* Last modified: 26/01/2015
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -42,4 +42,4 @@ One must have:
/*------------------------------- functions ---------------------------------*/
extern void compute_winpos(picstruct *field, picstruct *wfield,
objstruct *obj);
picstruct *dgeofield, objstruct *obj);
tests/default.param
View file @ deaee39b
......@@ -3,13 +3,14 @@ X_IMAGE
Y_IMAGE
XWIN_IMAGE
YWIN_IMAGE
XMODEL_IMAGE
YMODEL_IMAGE
FLUX_AUTO
FLUX_MODEL
#XMODEL_IMAGE
#YMODEL_IMAGE
MAG_AUTO
MAG_MODEL
#MAG_MODEL
FLUX_RADIUS
FLAGS
NITER_MODEL
FLAGS_WIN
NITER_WIN
#SPREAD_MODEL
#SPREADERR_MODEL
ID_PARENT
tests/default.sex
View file @ deaee39b
......@@ -5,7 +5,7 @@
#-------------------------------- Catalog ------------------------------------
CATALOG_NAME test.cat # name of the output catalog
CATALOG_TYPE ASCII_HEAD # NONE,ASCII,ASCII_HEAD, ASCII_SKYCAT,
CATALOG_TYPE FITS_LDAC # NONE,ASCII,ASCII_HEAD, ASCII_SKYCAT,
# ASCII_VOTABLE, FITS_1.0 or FITS_LDAC
PARAMETERS_NAME default.param # name of the file containing catalog contents
......@@ -17,7 +17,7 @@ DETECT_THRESH 1.5 # <sigmas> or <threshold>,<ZP> in mag.arcsec-2
ANALYSIS_THRESH 2 # <sigmas> or <threshold>,<ZP> in mag.arcsec-2
FILTER Y # apply filter for detection (Y or N)?
FILTER_NAME default.conv # name of the file containing the filter
FILTER_NAME gauss_4.0_7x7.conv # name of the file containing the filter
DEBLEND_NTHRESH 32 # Number of deblending sub-thresholds
DEBLEND_MINCONT 0.005 # Minimum contrast parameter for deblending
......@@ -58,11 +58,11 @@ BACKPHOTO_TYPE GLOBAL # can be GLOBAL or LOCAL
#------------------------------ Check Image ----------------------------------
CHECKIMAGE_TYPE MODELS, -MODELS, -BACKGROUND # can be NONE, BACKGROUND, BACKGROUND_RMS,
#CHECKIMAGE_TYPE MODELS, -MODELS, -BACKGROUND # can be NONE, BACKGROUND, BACKGROUND_RMS,
# MINIBACKGROUND, MINIBACK_RMS, -BACKGROUND,
# FILTERED, OBJECTS, -OBJECTS, SEGMENTATION,
# or APERTURES
CHECKIMAGE_NAME prof.fits,subprof.fits,orig.fits # Filename for the check-image
#CHECKIMAGE_NAME prof.fits,subprof.fits,orig.fits # Filename for the check-image
#--------------------- Memory (change with caution!) -------------------------
......@@ -70,6 +70,17 @@ MEMORY_OBJSTACK 3000 # number of objects in stack
MEMORY_PIXSTACK 300000 # number of pixels in stack
MEMORY_BUFSIZE 1024 # number of lines in buffer
#------------------------------- ASSOCiation ---------------------------------
ASSOC_NAME sky.list # name of the ASCII file to ASSOCiate
ASSOC_DATA 2,3,4 # columns of the data to replicate (0=all)
ASSOC_PARAMS 2,3,4 # columns of xpos,ypos[,mag]
ASSOCCOORD_TYPE PIXEL # ASSOC coordinates: PIXEL or WORLD
ASSOC_RADIUS 2.0 # cross-matching radius (pixels)
ASSOC_TYPE NEAREST # ASSOCiation method: FIRST, NEAREST, MEAN,
# MAG_MEAN, SUM, MAG_SUM, MIN or MAX
ASSOCSELEC_TYPE MATCHED # ASSOC selection type: ALL, MATCHED or -MATCHED
#----------------------------- Miscellaneous ---------------------------------
VERBOSE_TYPE NORMAL # can be QUIET, NORMAL or FULL
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment