/*
|
/*
|
* psf.c
|
* psf.c
|
*
|
*
|
* Fit a PSF model to an image.
|
* Fit a PSF model to an image.
|
*
|
*
|
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
*
|
*
|
* This file part of: SExtractor
|
* This file part of: SExtractor
|
*
|
*
|
* Copyright: (C) 1998-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
|
* Copyright: (C) 1998-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
|
*
|
*
|
* License: GNU General Public License
|
* License: GNU General Public License
|
*
|
*
|
* SExtractor is free software: you can redistribute it and/or modify
|
* SExtractor is free software: you can redistribute it and/or modify
|
* it under the terms of the GNU General Public License as published by
|
* it under the terms of the GNU General Public License as published by
|
* the Free Software Foundation, either version 3 of the License, or
|
* the Free Software Foundation, either version 3 of the License, or
|
* (at your option) any later version.
|
* (at your option) any later version.
|
* SExtractor is distributed in the hope that it will be useful,
|
* SExtractor is distributed in the hope that it will be useful,
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
* GNU General Public License for more details.
|
* GNU General Public License for more details.
|
* You should have received a copy of the GNU General Public License
|
* You should have received a copy of the GNU General Public License
|
* along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
|
* along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
|
*
|
*
|
* Last modified: 23/06/2011
|
* Last modified: 06/10/2011
|
*
|
*
|
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
|
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
|
|
|
#ifdef HAVE_CONFIG_H
|
#ifdef HAVE_CONFIG_H
|
#include "config.h"
|
#include "config.h"
|
#endif
|
#endif
|
|
|
#include <math.h>
|
#include <math.h>
|
#include <stdio.h>
|
#include <stdio.h>
|
#include <stdlib.h>
|
#include <stdlib.h>
|
#include <string.h>
|
#include <string.h>
|
|
|
#include "define.h"
|
#include "define.h"
|
#include "globals.h"
|
#include "globals.h"
|
#include "prefs.h"
|
#include "prefs.h"
|
#include "fits/fitscat.h"
|
#include "fits/fitscat.h"
|
#include "check.h"
|
#include "check.h"
|
#include "field.h"
|
#include "field.h"
|
#include "filter.h"
|
#include "filter.h"
|
#include "image.h"
|
#include "image.h"
|
#include "wcs/poly.h"
|
#include "wcs/poly.h"
|
#include "psf.h"
|
#include "psf.h"
|
|
|
/*------------------------------- variables ---------------------------------*/
|
/*------------------------------- variables ---------------------------------*/
|
|
|
|
|
extern keystruct objkey[];
|
extern keystruct objkey[];
|
extern objstruct outobj;
|
extern objstruct outobj;
|
|
|
/********************************* psf_init **********************************/
|
/********************************* psf_init **********************************/
|
/*
|
/*
|
Allocate memory and stuff for the PSF-fitting.
|
Allocate memory and stuff for the PSF-fitting.
|
*/
|
*/
|
void psf_init(psfstruct *psf)
|
void psf_init(psfstruct *psf)
|
{
|
{
|
QMALLOC(thepsfit, psfitstruct, 1);
|
QMALLOC(thepsfit, psfitstruct, 1);
|
QMALLOC(thepsfit->x, double, prefs.psf_npsfmax);
|
QMALLOC(thepsfit->x, double, prefs.psf_npsfmax);
|
QMALLOC(thepsfit->y, double, prefs.psf_npsfmax);
|
QMALLOC(thepsfit->y, double, prefs.psf_npsfmax);
|
QMALLOC(thepsfit->flux, float, prefs.psf_npsfmax);
|
QMALLOC(thepsfit->flux, float, prefs.psf_npsfmax);
|
QMALLOC(thepsfit->fluxerr, float, prefs.psf_npsfmax);
|
QMALLOC(thepsfit->fluxerr, float, prefs.psf_npsfmax);
|
QMALLOC(ppsfit, psfitstruct, 1); /*?*/
|
QMALLOC(ppsfit, psfitstruct, 1); /*?*/
|
QMALLOC(ppsfit->x, double, prefs.psf_npsfmax);
|
QMALLOC(ppsfit->x, double, prefs.psf_npsfmax);
|
QMALLOC(ppsfit->y, double, prefs.psf_npsfmax);
|
QMALLOC(ppsfit->y, double, prefs.psf_npsfmax);
|
QMALLOC(ppsfit->flux, float, prefs.psf_npsfmax);
|
QMALLOC(ppsfit->flux, float, prefs.psf_npsfmax);
|
QMALLOC(ppsfit->fluxerr, float, prefs.psf_npsfmax);
|
QMALLOC(ppsfit->fluxerr, float, prefs.psf_npsfmax);
|
|
|
return;
|
return;
|
}
|
}
|
|
|
|
|
/********************************* psf_end ***********************************/
|
/********************************* psf_end ***********************************/
|
/*
|
/*
|
Free memory occupied by the PSF-fitting stuff.
|
Free memory occupied by the PSF-fitting stuff.
|
*/
|
*/
|
void psf_end(psfstruct *psf, psfitstruct *psfit)
|
void psf_end(psfstruct *psf, psfitstruct *psfit)
|
{
|
{
|
int d, ndim;
|
int d, ndim;
|
|
|
if (psf->pc)
|
|
pc_end(psf->pc);
|
|
|
|
ndim = psf->poly->ndim;
|
ndim = psf->poly->ndim;
|
for (d=0; d<ndim; d++)
|
for (d=0; d<ndim; d++)
|
free(psf->contextname[d]);
|
free(psf->contextname[d]);
|
free(psf->context);
|
free(psf->context);
|
free(psf->contextname);
|
free(psf->contextname);
|
free(psf->contextoffset);
|
free(psf->contextoffset);
|
free(psf->contextscale);
|
free(psf->contextscale);
|
free(psf->contexttyp);
|
free(psf->contexttyp);
|
poly_end(psf->poly);
|
poly_end(psf->poly);
|
free(psf->maskcomp);
|
free(psf->maskcomp);
|
free(psf->maskloc);
|
free(psf->maskloc);
|
free(psf->masksize);
|
free(psf->masksize);
|
free(psf);
|
free(psf);
|
|
|
if (psfit)
|
if (psfit)
|
{
|
{
|
free(psfit->x);
|
free(psfit->x);
|
free(psfit->y);
|
free(psfit->y);
|
free(psfit->flux);
|
free(psfit->flux);
|
free(psfit->fluxerr);
|
free(psfit->fluxerr);
|
free(psfit);
|
free(psfit);
|
}
|
}
|
|
|
return;
|
return;
|
}
|
}
|
|
|
|
|
/********************************* psf_load *********************************/
|
/********************************* psf_load *********************************/
|
/*
|
/*
|
Read the PSF data from a FITS file.
|
Read the PSF data from a FITS file.
|
*/
|
*/
|
psfstruct *psf_load(char *filename)
|
psfstruct *psf_load(char *filename)
|
{
|
{
|
static objstruct saveobj;
|
|
static obj2struct saveobj2;
|
static obj2struct saveobj2;
|
psfstruct *psf;
|
psfstruct *psf;
|
catstruct *cat;
|
catstruct *cat;
|
tabstruct *tab;
|
tabstruct *tab;
|
keystruct *key;
|
keystruct *key;
|
char *head, *ci,*co;
|
char *head, *ci,*co;
|
int deg[POLY_MAXDIM], group[POLY_MAXDIM], ndim, ngroup,
|
int deg[POLY_MAXDIM], group[POLY_MAXDIM], ndim, ngroup,
|
i,k;
|
i,k;
|
|
|
/* Open the cat (well it is not a "cat", but simply a FITS file */
|
/* Open the cat (well it is not a "cat", but simply a FITS file */
|
if (!(cat = read_cat(filename)))
|
if (!(cat = read_cat(filename)))
|
error(EXIT_FAILURE, "*Error*: PSF file not found: ", filename);
|
error(EXIT_FAILURE, "*Error*: PSF file not found: ", filename);
|
|
|
/* OK, we now allocate memory for the PSF structure itself */
|
/* OK, we now allocate memory for the PSF structure itself */
|
QCALLOC(psf, psfstruct, 1);
|
QCALLOC(psf, psfstruct, 1);
|
|
|
/* Store a short copy of the PSF filename */
|
/* Store a short copy of the PSF filename */
|
if ((ci=strrchr(filename, '/')))
|
if ((ci=strrchr(filename, '/')))
|
strcpy(psf->name, ci+1);
|
strcpy(psf->name, ci+1);
|
else
|
else
|
strcpy(psf->name, filename);
|
strcpy(psf->name, filename);
|
|
|
if (!(tab = name_to_tab(cat, "PSF_DATA", 0)))
|
if (!(tab = name_to_tab(cat, "PSF_DATA", 0)))
|
error(EXIT_FAILURE, "*Error*: PSF_DATA table not found in catalog ",
|
error(EXIT_FAILURE, "*Error*: PSF_DATA table not found in catalog ",
|
filename);
|
filename);
|
|
|
head = tab->headbuf;
|
head = tab->headbuf;
|
|
|
/*-- Dimension of the polynomial */
|
/*-- Dimension of the polynomial */
|
if (fitsread(head, "POLNAXIS", &ndim, H_INT,T_LONG) == RETURN_OK
|
if (fitsread(head, "POLNAXIS", &ndim, H_INT,T_LONG) == RETURN_OK
|
&& ndim)
|
&& ndim)
|
{
|
{
|
/*-- So we have a polynomial description of the PSF variations */
|
/*-- So we have a polynomial description of the PSF variations */
|
if (ndim > POLY_MAXDIM)
|
if (ndim > POLY_MAXDIM)
|
{
|
{
|
sprintf(gstr, "*Error*: The POLNAXIS parameter in %s exceeds %d",
|
sprintf(gstr, "*Error*: The POLNAXIS parameter in %s exceeds %d",
|
psf->name, POLY_MAXDIM);
|
psf->name, POLY_MAXDIM);
|
error(EXIT_FAILURE, gstr, "");
|
error(EXIT_FAILURE, gstr, "");
|
}
|
}
|
|
|
QMALLOC(psf->contextname, char *, ndim);
|
QMALLOC(psf->contextname, char *, ndim);
|
QMALLOC(psf->context, double *, ndim);
|
QMALLOC(psf->context, double *, ndim);
|
QMALLOC(psf->contexttyp, t_type, ndim);
|
QMALLOC(psf->contexttyp, t_type, ndim);
|
QMALLOC(psf->contextoffset, double, ndim);
|
QMALLOC(psf->contextoffset, double, ndim);
|
QMALLOC(psf->contextscale, double, ndim);
|
QMALLOC(psf->contextscale, double, ndim);
|
|
|
/*-- We will have to use the outobj structs, so we first save their content */
|
/*-- We will have to use the flagobj2 struct, so we first save its content */
|
saveobj = outobj;
|
saveobj2 = flagobj2;
|
saveobj2 = outobj2;
|
/*-- flagobj2 is used as a FLAG array, so we initialize it to 0 */
|
/*-- outobj's are used as FLAG arrays, so we initialize them to 0 */
|
memset(&flagobj2, 0, sizeof(flagobj2));
|
memset(&outobj, 0, sizeof(outobj));
|
|
memset(&outobj2, 0, sizeof(outobj2));
|
|
for (i=0; i<ndim; i++)
|
for (i=0; i<ndim; i++)
|
{
|
{
|
/*---- Polynomial groups */
|
/*---- Polynomial groups */
|
sprintf(gstr, "POLGRP%1d", i+1);
|
sprintf(gstr, "POLGRP%1d", i+1);
|
if (fitsread(head, gstr, &group[i], H_INT,T_LONG) != RETURN_OK)
|
if (fitsread(head, gstr, &group[i], H_INT,T_LONG) != RETURN_OK)
|
goto headerror;
|
goto headerror;
|
|
|
/*---- Contexts */
|
/*---- Contexts */
|
QMALLOC(psf->contextname[i], char, 80);
|
QMALLOC(psf->contextname[i], char, 80);
|
sprintf(gstr, "POLNAME%1d", i+1);
|
sprintf(gstr, "POLNAME%1d", i+1);
|
if (fitsread(head,gstr,psf->contextname[i],H_STRING,T_STRING)!=RETURN_OK)
|
if (fitsread(head,gstr,psf->contextname[i],H_STRING,T_STRING)!=RETURN_OK)
|
goto headerror;
|
goto headerror;
|
if (*psf->contextname[i]==(char)':')
|
if (*psf->contextname[i]==(char)':')
|
/*------ It seems we're facing a FITS header parameter */
|
/*------ It seems we're facing a FITS header parameter */
|
psf->context[i] = NULL; /* This is to tell we'll have to load */
|
psf->context[i] = NULL; /* This is to tell we'll have to load */
|
/* a FITS header context later on */
|
/* a FITS header context later on */
|
else
|
else
|
/*------ The context element is a dynamic object parameter */
|
/*------ The context element is a dynamic object parameter */
|
{
|
{
|
if ((k = findkey(psf->contextname[i], (char *)objkey,
|
if ((k = findkey(psf->contextname[i], (char *)objkey,
|
sizeof(keystruct)))==RETURN_ERROR)
|
sizeof(keystruct)))==RETURN_ERROR)
|
{
|
{
|
sprintf(gstr, "*Error*: %s CONTEXT parameter in %s unknown",
|
sprintf(gstr, "*Error*: %s CONTEXT parameter in %s unknown",
|
psf->contextname[i], psf->name);
|
psf->contextname[i], psf->name);
|
error(EXIT_FAILURE, gstr, "");
|
error(EXIT_FAILURE, gstr, "");
|
}
|
}
|
key = objkey+k;
|
key = objkey+k;
|
psf->context[i] = key->ptr;
|
psf->context[i] = key->ptr;
|
psf->contexttyp[i] = key->ttype;
|
psf->contexttyp[i] = key->ttype;
|
/*------ Declare the parameter "active" to trigger computation by SExtractor */
|
/*------ Declare the parameter "active" to trigger computation by SExtractor */
|
*((char *)key->ptr) = (char)'\1';
|
*((char *)key->ptr) = (char)'\1';
|
}
|
}
|
/*---- Scaling of the context parameter */
|
/*---- Scaling of the context parameter */
|
sprintf(gstr, "POLZERO%1d", i+1);
|
sprintf(gstr, "POLZERO%1d", i+1);
|
if (fitsread(head, gstr, &psf->contextoffset[i], H_EXPO, T_DOUBLE)
|
if (fitsread(head, gstr, &psf->contextoffset[i], H_EXPO, T_DOUBLE)
|
!=RETURN_OK)
|
!=RETURN_OK)
|
goto headerror;
|
goto headerror;
|
sprintf(gstr, "POLSCAL%1d", i+1);
|
sprintf(gstr, "POLSCAL%1d", i+1);
|
if (fitsread(head, gstr, &psf->contextscale[i], H_EXPO, T_DOUBLE)
|
if (fitsread(head, gstr, &psf->contextscale[i], H_EXPO, T_DOUBLE)
|
!=RETURN_OK)
|
!=RETURN_OK)
|
goto headerror;
|
goto headerror;
|
}
|
}
|
|
|
/*-- Number of groups */
|
/*-- Number of groups */
|
if (fitsread(head, "POLNGRP ", &ngroup, H_INT, T_LONG) != RETURN_OK)
|
if (fitsread(head, "POLNGRP ", &ngroup, H_INT, T_LONG) != RETURN_OK)
|
goto headerror;
|
goto headerror;
|
|
|
for (i=0; i<ngroup; i++)
|
for (i=0; i<ngroup; i++)
|
{
|
{
|
/*---- Polynomial degree for each group */
|
/*---- Polynomial degree for each group */
|
sprintf(gstr, "POLDEG%1d", i+1);
|
sprintf(gstr, "POLDEG%1d", i+1);
|
if (fitsread(head, gstr, °[i], H_INT,T_LONG) != RETURN_OK)
|
if (fitsread(head, gstr, °[i], H_INT,T_LONG) != RETURN_OK)
|
goto headerror;
|
goto headerror;
|
}
|
}
|
|
|
psf->poly = poly_init(group, ndim, deg, ngroup);
|
psf->poly = poly_init(group, ndim, deg, ngroup);
|
|
|
/*-- Update the permanent FLAG arrays (that is, perform an "OR" on them) */
|
/*-- Restore previous flagobj2 content */
|
for (ci=(char *)&outobj,co=(char *)&flagobj,i=sizeof(objstruct); i--;)
|
flagobj2 = saveobj2;
|
*(co++) |= *(ci++);
|
|
for (ci=(char *)&outobj2,co=(char *)&flagobj2,i=sizeof(obj2struct); i--;)
|
|
*(co++) |= *(ci++);
|
|
|
|
/*-- Restore previous outobj contents */
|
|
outobj = saveobj;
|
|
outobj2 = saveobj2;
|
|
}
|
}
|
else
|
else
|
{
|
{
|
/*-- This is a simple, constant PSF */
|
/*-- This is a simple, constant PSF */
|
psf->poly = poly_init(group, 0, deg, 0);
|
psf->poly = poly_init(group, 0, deg, 0);
|
psf->context = NULL;
|
psf->context = NULL;
|
}
|
}
|
|
|
/* Dimensionality of the PSF mask */
|
/* Dimensionality of the PSF mask */
|
if (fitsread(head, "PSFNAXIS", &psf->maskdim, H_INT, T_LONG) != RETURN_OK)
|
if (fitsread(head, "PSFNAXIS", &psf->maskdim, H_INT, T_LONG) != RETURN_OK)
|
goto headerror;
|
goto headerror;
|
if (psf->maskdim<2 || psf->maskdim>3)
|
if (psf->maskdim<2 || psf->maskdim>3)
|
error(EXIT_FAILURE, "*Error*: wrong dimensionality for the PSF "
|
error(EXIT_FAILURE, "*Error*: wrong dimensionality for the PSF "
|
"mask in ", filename);
|
"mask in ", filename);
|
QMALLOC(psf->masksize, int, psf->maskdim);
|
QMALLOC(psf->masksize, int, psf->maskdim);
|
for (i=0; i<psf->maskdim; i++)
|
for (i=0; i<psf->maskdim; i++)
|
psf->masksize[i] = 1;
|
psf->masksize[i] = 1;
|
psf->masknpix = 1;
|
psf->masknpix = 1;
|
for (i=0; i<psf->maskdim; i++)
|
for (i=0; i<psf->maskdim; i++)
|
{
|
{
|
sprintf(gstr, "PSFAXIS%1d", i+1);
|
sprintf(gstr, "PSFAXIS%1d", i+1);
|
if (fitsread(head, gstr, &psf->masksize[i], H_INT,T_LONG) != RETURN_OK)
|
if (fitsread(head, gstr, &psf->masksize[i], H_INT,T_LONG) != RETURN_OK)
|
goto headerror;
|
goto headerror;
|
psf->masknpix *= psf->masksize[i];
|
psf->masknpix *= psf->masksize[i];
|
}
|
}
|
|
|
/* PSF FWHM: defaulted to 3 pixels */
|
/* PSF FWHM: defaulted to 3 pixels */
|
if (fitsread(head, "PSF_FWHM", &psf->fwhm, H_FLOAT,T_DOUBLE) != RETURN_OK)
|
if (fitsread(head, "PSF_FWHM", &psf->fwhm, H_FLOAT,T_DOUBLE) != RETURN_OK)
|
psf->fwhm = 3.0;
|
psf->fwhm = 3.0;
|
|
|
/* PSF oversampling: defaulted to 1 */
|
/* PSF oversampling: defaulted to 1 */
|
if (fitsread(head, "PSF_SAMP", &psf->pixstep,H_FLOAT,T_FLOAT) != RETURN_OK
|
if (fitsread(head, "PSF_SAMP", &psf->pixstep,H_FLOAT,T_FLOAT) != RETURN_OK
|
|| psf->pixstep <= 0.0)
|
|| psf->pixstep <= 0.0)
|
psf->pixstep = 1.0;
|
psf->pixstep = 1.0;
|
|
|
/* Load the PSF mask data */
|
/* Load the PSF mask data */
|
key = read_key(tab, "PSF_MASK");
|
key = read_key(tab, "PSF_MASK");
|
psf->maskcomp = key->ptr;
|
psf->maskcomp = key->ptr;
|
|
|
psf->pc = pc_load(cat);
|
|
|
|
QMALLOC(psf->maskloc, float, psf->masksize[0]*psf->masksize[1]);
|
QMALLOC(psf->maskloc, float, psf->masksize[0]*psf->masksize[1]);
|
|
|
/* But don't touch my arrays!! */
|
/* But don't touch my arrays!! */
|
blank_keys(tab);
|
blank_keys(tab);
|
|
|
free_cat(&cat, 1);
|
free_cat(&cat, 1);
|
|
|
return psf;
|
return psf;
|
|
|
headerror:
|
headerror:
|
error(EXIT_FAILURE, "*Error*: Incorrect or obsolete PSF data in ", filename);
|
error(EXIT_FAILURE, "*Error*: Incorrect or obsolete PSF data in ", filename);
|
return NULL;
|
return NULL;
|
}
|
}
|
|
|
|
|
/***************************** psf_readcontext *******************************/
|
/***************************** psf_readcontext *******************************/
|
/*
|
/*
|
Read the PSF context parameters in the FITS header.
|
Read the PSF context parameters in the FITS header.
|
*/
|
*/
|
void psf_readcontext(psfstruct *psf, picstruct *field)
|
void psf_readcontext(psfstruct *psf, picstruct *field)
|
{
|
{
|
static double contextval[POLY_MAXDIM];
|
static double contextval[POLY_MAXDIM];
|
int i, ndim;
|
int i, ndim;
|
|
|
ndim = psf->poly->ndim;
|
ndim = psf->poly->ndim;
|
for (i=0; i<ndim; i++)
|
for (i=0; i<ndim; i++)
|
if (!psf->context[i])
|
if (!psf->context[i])
|
{
|
{
|
psf->context[i] = &contextval[i];
|
psf->context[i] = &contextval[i];
|
psf->contexttyp[i] = T_DOUBLE;
|
psf->contexttyp[i] = T_DOUBLE;
|
if (fitsread(field->tab->headbuf, psf->contextname[i]+1, &contextval[i],
|
if (fitsread(field->tab->headbuf, psf->contextname[i]+1, &contextval[i],
|
H_FLOAT,T_DOUBLE) == RETURN_ERROR)
|
H_FLOAT,T_DOUBLE) == RETURN_ERROR)
|
{
|
{
|
sprintf(gstr, "*Error*: %s parameter not found in the header of ",
|
sprintf(gstr, "*Error*: %s parameter not found in the header of ",
|
psf->contextname[i]+1);
|
psf->contextname[i]+1);
|
error(EXIT_FAILURE, gstr, field->rfilename);
|
error(EXIT_FAILURE, gstr, field->rfilename);
|
}
|
}
|
}
|
}
|
|
|
return;
|
return;
|
}
|
}
|
|
|
|
|
/******************************** psf_fit ***********************************/
|
/******************************** psf_fit ***********************************/
|
/* standard PSF fit for one component */
|
/* standard PSF fit for one component */
|
/****************************************************************************/
|
/****************************************************************************/
|
|
|
void psf_fit(psfstruct *psf, picstruct *field, picstruct *wfield,
|
void psf_fit(psfstruct *psf, picstruct *field, picstruct *wfield,
|
objstruct *obj, obj2struct *obj2)
|
obj2struct *obj2)
|
{
|
{
|
checkstruct *check;
|
checkstruct *check;
|
static double x2[PSF_NPSFMAX],y2[PSF_NPSFMAX],xy[PSF_NPSFMAX],
|
static double x2[PSF_NPSFMAX],y2[PSF_NPSFMAX],xy[PSF_NPSFMAX],
|
deltax[PSF_NPSFMAX],
|
deltax[PSF_NPSFMAX],
|
deltay[PSF_NPSFMAX],
|
deltay[PSF_NPSFMAX],
|
flux[PSF_NPSFMAX],fluxerr[PSF_NPSFMAX],
|
flux[PSF_NPSFMAX],fluxerr[PSF_NPSFMAX],
|
deltaxb[PSF_NPSFMAX],deltayb[PSF_NPSFMAX],
|
deltaxb[PSF_NPSFMAX],deltayb[PSF_NPSFMAX],
|
fluxb[PSF_NPSFMAX],fluxerrb[PSF_NPSFMAX],
|
fluxb[PSF_NPSFMAX],fluxerrb[PSF_NPSFMAX],
|
sol[PSF_NTOT], covmat[PSF_NTOT*PSF_NTOT],
|
sol[PSF_NTOT], covmat[PSF_NTOT*PSF_NTOT],
|
vmat[PSF_NTOT*PSF_NTOT], wmat[PSF_NTOT];
|
vmat[PSF_NTOT*PSF_NTOT], wmat[PSF_NTOT];
|
float *data, *data2, *data3, *weight, *d, *w;
|
float *data, *data2, *data3, *weight, *d, *w;
|
double *mat,
|
double *mat,
|
*m, *var,
|
*m, *var,
|
dx,dy,
|
dx,dy,
|
pix,pix2, wthresh,val,
|
pix,pix2, wthresh,val,
|
backnoise2, gain, radmin2,radmax2,satlevel,chi2,
|
backnoise2, gain, radmin2,radmax2,satlevel,chi2,
|
r2, valmax, psf_fwhm;
|
r2, valmax, psf_fwhm;
|
float **psfmasks, **psfmaskx,**psfmasky,
|
float **psfmasks, **psfmaskx,**psfmasky,
|
*ps, *dh, *wh, pixstep;
|
*ps, *dh, *wh, pixstep;
|
PIXTYPE *datah, *weighth;
|
PIXTYPE *datah, *weighth;
|
int i,j,p, npsf,npsfmax, npix, nppix, ix,iy,niter,
|
int i,j,p, npsf,npsfmax, npix, nppix, ix,iy,niter,
|
width, height, pwidth,pheight, x,y,
|
width, height, pwidth,pheight, x,y,
|
xmax,ymax, wbad, gainflag, convflag, npsfflag,
|
xmax,ymax, wbad, gainflag, convflag, npsfflag,
|
ival,kill=0;
|
ival,kill=0;
|
|
|
dx = dy = 0.0;
|
dx = dy = 0.0;
|
niter = 0;
|
niter = 0;
|
npsfmax = prefs.psf_npsfmax;
|
npsfmax = prefs.psf_npsfmax;
|
pixstep = 1.0/psf->pixstep;
|
pixstep = 1.0/psf->pixstep;
|
gain = (field->gain >0.0? field->gain: 1e30);
|
gain = (field->gain >0.0? field->gain: 1e30);
|
backnoise2 = field->backsig*field->backsig;
|
backnoise2 = field->backsig*field->backsig;
|
satlevel = field->satur_level - obj->bkg;
|
satlevel = field->satur_level - obj2->bkg;
|
wthresh = wfield?wfield->weight_thresh:BIG;
|
wthresh = wfield?wfield->weight_thresh:BIG;
|
gainflag = prefs.weightgain_flag;
|
gainflag = prefs.weightgain_flag;
|
psf_fwhm = psf->fwhm*psf->pixstep;
|
psf_fwhm = psf->fwhm*psf->pixstep;
|
|
|
|
|
/* Initialize outputs */
|
/* Initialize outputs */
|
thepsfit->niter = 0;
|
thepsfit->niter = 0;
|
thepsfit->npsf = 0;
|
thepsfit->npsf = 0;
|
for (j=0; j<npsfmax; j++)
|
for (j=0; j<npsfmax; j++)
|
{
|
{
|
thepsfit->x[j] = obj2->posx;
|
thepsfit->x[j] = obj2->posx;
|
thepsfit->y[j] = obj2->posy;
|
thepsfit->y[j] = obj2->posy;
|
thepsfit->flux[j] = 0.0;
|
thepsfit->flux[j] = 0.0;
|
thepsfit->fluxerr[j] = 0.0;
|
thepsfit->fluxerr[j] = 0.0;
|
}
|
}
|
|
|
/* Scale data area with object "size" */
|
/* Scale data area with object "size" */
|
ix = (obj->xmax+obj->xmin+1)/2;
|
ix = (obj2->xmax+obj2->xmin+1)/2;
|
iy = (obj->ymax+obj->ymin+1)/2;
|
iy = (obj2->ymax+obj2->ymin+1)/2;
|
width = obj->xmax-obj->xmin+1+psf_fwhm;
|
width = obj2->xmax-obj2->xmin+1+psf_fwhm;
|
if (width < (ival=(int)(psf_fwhm*2)))
|
if (width < (ival=(int)(psf_fwhm*2)))
|
width = ival;
|
width = ival;
|
height = obj->ymax-obj->ymin+1+psf_fwhm;
|
height = obj2->ymax-obj2->ymin+1+psf_fwhm;
|
if (height < (ival=(int)(psf_fwhm*2)))
|
if (height < (ival=(int)(psf_fwhm*2)))
|
height = ival;
|
height = ival;
|
npix = width*height;
|
npix = width*height;
|
radmin2 = PSF_MINSHIFT*PSF_MINSHIFT;
|
radmin2 = PSF_MINSHIFT*PSF_MINSHIFT;
|
radmax2 = npix/2.0;
|
radmax2 = npix/2.0;
|
|
|
/* Scale total area with PSF FWHM */
|
/* Scale total area with PSF FWHM */
|
pwidth = (int)(psf->masksize[0]*psf->pixstep)+width;;
|
pwidth = (int)(psf->masksize[0]*psf->pixstep)+width;;
|
pheight = (int)(psf->masksize[1]*psf->pixstep)+height;
|
pheight = (int)(psf->masksize[1]*psf->pixstep)+height;
|
nppix = pwidth*pheight;
|
nppix = pwidth*pheight;
|
|
|
QMALLOC(weighth, PIXTYPE, npix);
|
QMALLOC(weighth, PIXTYPE, npix);
|
QMALLOC(weight, float, npix);
|
QMALLOC(weight, float, npix);
|
QMALLOC(datah, PIXTYPE, npix);
|
QMALLOC(datah, PIXTYPE, npix);
|
QMALLOC(data, float, npix);
|
QMALLOC(data, float, npix);
|
QMALLOC(data2, float, npix);
|
QMALLOC(data2, float, npix);
|
QMALLOC(data3, float, npix);
|
QMALLOC(data3, float, npix);
|
QMALLOC(mat, double, npix*PSF_NTOT);
|
QMALLOC(mat, double, npix*PSF_NTOT);
|
if (prefs.check[CHECK_SUBPSFPROTOS] || prefs.check[CHECK_PSFPROTOS]
|
if (prefs.check[CHECK_SUBPSFPROTOS] || prefs.check[CHECK_PSFPROTOS])
|
|| prefs.check[CHECK_SUBPCPROTOS] || prefs.check[CHECK_PCPROTOS]
|
|
|| prefs.check[CHECK_PCOPROTOS])
|
|
{
|
{
|
QMALLOC(checkmask, PIXTYPE, nppix);
|
QMALLOC(checkmask, PIXTYPE, nppix);
|
}
|
}
|
|
|
QMALLOC(psfmasks, float *, npsfmax);
|
QMALLOC(psfmasks, float *, npsfmax);
|
QMALLOC(psfmaskx, float *, npsfmax);
|
QMALLOC(psfmaskx, float *, npsfmax);
|
QMALLOC(psfmasky, float *, npsfmax);
|
QMALLOC(psfmasky, float *, npsfmax);
|
for (i=0; i<npsfmax; i++)
|
for (i=0; i<npsfmax; i++)
|
{
|
{
|
QMALLOC(psfmasks[i], float, npix);
|
QMALLOC(psfmasks[i], float, npix);
|
QMALLOC(psfmaskx[i], float, npix);
|
QMALLOC(psfmaskx[i], float, npix);
|
QMALLOC(psfmasky[i], float, npix);
|
QMALLOC(psfmasky[i], float, npix);
|
}
|
}
|
|
|
/* Compute weights */
|
/* Compute weights */
|
wbad = 0;
|
wbad = 0;
|
if (wfield)
|
if (wfield)
|
{
|
{
|
psf_copyobjpix(datah, weighth, width, height, ix,iy, obj2,
|
psf_copyobjpix(datah, weighth, width, height, ix,iy, obj2,
|
!(field->flags&MEASURE_FIELD));
|
!(field->flags&MEASURE_FIELD));
|
for (wh=weighth, w=weight, dh=datah,p=npix; p--;)
|
for (wh=weighth, w=weight, dh=datah,p=npix; p--;)
|
if ((pix=*(wh++)) < wthresh && pix>0
|
if ((pix=*(wh++)) < wthresh && pix>0
|
&& (pix2=*(dh++))>-BIG
|
&& (pix2=*(dh++))>-BIG
|
&& pix2<satlevel)
|
&& pix2<satlevel)
|
*(w++) = 1/sqrt(pix+(pix2>0.0?
|
*(w++) = 1/sqrt(pix+(pix2>0.0?
|
(gainflag? pix2*pix/backnoise2:pix2)/gain
|
(gainflag? pix2*pix/backnoise2:pix2)/gain
|
:0.0));
|
:0.0));
|
else
|
else
|
{
|
{
|
*(w++) = 0.0;
|
*(w++) = 0.0;
|
wbad++;
|
wbad++;
|
}
|
}
|
}
|
}
|
else
|
else
|
{
|
{
|
psf_copyobjpix(datah, NULL, width, height, ix,iy, obj2,
|
psf_copyobjpix(datah, NULL, width, height, ix,iy, obj2,
|
!(field->flags&MEASURE_FIELD));
|
!(field->flags&MEASURE_FIELD));
|
for (w=weight, dh=datah, p=npix; p--;)
|
for (w=weight, dh=datah, p=npix; p--;)
|
if ((pix=*(dh++))>-BIG && pix<satlevel)
|
if ((pix=*(dh++))>-BIG && pix<satlevel)
|
*(w++) = 1.0/sqrt(backnoise2+(pix>0.0?pix/gain:0.0));
|
*(w++) = 1.0/sqrt(backnoise2+(pix>0.0?pix/gain:0.0));
|
else
|
else
|
{
|
{
|
*(w++) = 0.0;
|
*(w++) = 0.0;
|
wbad++;
|
wbad++;
|
}
|
}
|
}
|
}
|
|
|
/* Special action if most of the weights are zero!! */
|
/* Special action if most of the weights are zero!! */
|
if (wbad>=npix-3)
|
if (wbad>=npix-3)
|
return;
|
return;
|
|
|
/* Weight the data */
|
/* Weight the data */
|
dh = datah;
|
dh = datah;
|
val = obj->dbkg; /* Take into account a local background change */
|
val = obj2->dbkg; /* Take into account a local background change */
|
d = data;
|
d = data;
|
w = weight;
|
w = weight;
|
for (p=npix; p--;)
|
for (p=npix; p--;)
|
*(d++) = (*(dh++)-val)**(w++);
|
*(d++) = (*(dh++)-val)**(w++);
|
|
|
/* Get the local PSF */
|
/* Get the local PSF */
|
psf_build(psf);
|
psf_build(psf, obj2);
|
|
|
npsfflag = 1;
|
npsfflag = 1;
|
r2 = psf_fwhm*psf_fwhm/2.0;
|
r2 = psf_fwhm*psf_fwhm/2.0;
|
fluxb[0] = fluxerrb[0] = deltaxb[0] = deltayb[0] = 0.0;
|
fluxb[0] = fluxerrb[0] = deltaxb[0] = deltayb[0] = 0.0;
|
|
|
for (npsf=1; npsf<=npsfmax && npsfflag; npsf++)
|
for (npsf=1; npsf<=npsfmax && npsfflag; npsf++)
|
{
|
{
|
kill=0;
|
kill=0;
|
/*-- First compute an optimum initial guess for the positions of components */
|
/*-- First compute an optimum initial guess for the positions of components */
|
if (npsf>1)
|
if (npsf>1)
|
{
|
{
|
/*---- Subtract previously fitted components */
|
/*---- Subtract previously fitted components */
|
d = data2;
|
d = data2;
|
dh = datah;
|
dh = datah;
|
for (p=npix; p--;)
|
for (p=npix; p--;)
|
*(d++) = (double)*(dh++);
|
*(d++) = (double)*(dh++);
|
for (j=0; j<npsf-1; j++)
|
for (j=0; j<npsf-1; j++)
|
{
|
{
|
d = data2;
|
d = data2;
|
ps = psfmasks[j];
|
ps = psfmasks[j];
|
for (p=npix; p--;)
|
for (p=npix; p--;)
|
*(d++) -= flux[j]**(ps++);
|
*(d++) -= flux[j]**(ps++);
|
}
|
}
|
convolve_image(field, data2, data3, width,height);
|
convolve_image(field, data2, data3, width,height);
|
/*---- Ignore regions too close to stellar cores */
|
/*---- Ignore regions too close to stellar cores */
|
for (j=0; j<npsf-1; j++)
|
for (j=0; j<npsf-1; j++)
|
{
|
{
|
d = data3;
|
d = data3;
|
dy = -((double)(height/2)+deltay[j]);
|
dy = -((double)(height/2)+deltay[j]);
|
for (y=height; y--; dy += 1.0)
|
for (y=height; y--; dy += 1.0)
|
{
|
{
|
dx = -((double)(width/2)+deltax[j]);
|
dx = -((double)(width/2)+deltax[j]);
|
for (x=width; x--; dx+= 1.0, d++)
|
for (x=width; x--; dx+= 1.0, d++)
|
if (dx*dx+dy*dy<r2)
|
if (dx*dx+dy*dy<r2)
|
*d = -BIG;
|
*d = -BIG;
|
}
|
}
|
}
|
}
|
/*---- Now find the brightest pixel (poor man's guess, to be refined later) */
|
/*---- Now find the brightest pixel (poor man's guess, to be refined later) */
|
d = data3;
|
d = data3;
|
valmax = -BIG;
|
valmax = -BIG;
|
xmax = width/2;
|
xmax = width/2;
|
ymax = height/2;
|
ymax = height/2;
|
for (y=0; y<height; y++)
|
for (y=0; y<height; y++)
|
for (x=0; x<width; x++)
|
for (x=0; x<width; x++)
|
{
|
{
|
if ((val = *(d++))>valmax)
|
if ((val = *(d++))>valmax)
|
{
|
{
|
valmax = val;
|
valmax = val;
|
xmax = x;
|
xmax = x;
|
ymax = y;
|
ymax = y;
|
}
|
}
|
}
|
}
|
deltax[npsf-1] = (double)(xmax - width/2);
|
deltax[npsf-1] = (double)(xmax - width/2);
|
deltay[npsf-1] = (double)(ymax - height/2);
|
deltay[npsf-1] = (double)(ymax - height/2);
|
}
|
}
|
else
|
else
|
{
|
{
|
/*---- Only one component to fit: simply use the barycenter as a guess */
|
/*---- Only one component to fit: simply use the barycenter as a guess */
|
deltax[npsf-1] = obj->mx - ix;
|
deltax[npsf-1] = obj2->mx - ix;
|
deltay[npsf-1] = obj->my - iy;
|
deltay[npsf-1] = obj2->my - iy;
|
}
|
}
|
|
|
niter = 0;
|
niter = 0;
|
convflag = 1;
|
convflag = 1;
|
for (i=0; i<PSF_NITER && convflag; i++)
|
for (i=0; i<PSF_NITER && convflag; i++)
|
{
|
{
|
convflag = 0,niter++,m=mat;
|
convflag = 0,niter++,m=mat;
|
for (j=0; j<npsf; j++)
|
for (j=0; j<npsf; j++)
|
{
|
{
|
/*------ Resample the PSFs here for the 1st iteration */
|
/*------ Resample the PSFs here for the 1st iteration */
|
vignet_resample(psf->maskloc, psf->masksize[0], psf->masksize[1],
|
vignet_resample(psf->maskloc, psf->masksize[0], psf->masksize[1],
|
psfmasks[j], width, height,
|
psfmasks[j], width, height,
|
-deltax[j]*pixstep, -deltay[j]*pixstep,
|
-deltax[j]*pixstep, -deltay[j]*pixstep,
|
pixstep);
|
pixstep);
|
m=compute_gradient(weight,width,height,
|
m=compute_gradient(weight,width,height,
|
psfmasks[j],psfmaskx[j],psfmasky[j],m);
|
psfmasks[j],psfmaskx[j],psfmasky[j],m);
|
}
|
}
|
|
|
|
|
svdfit(mat, data, npix, npsf*PSF_NA, sol, vmat, wmat);
|
svdfit(mat, data, npix, npsf*PSF_NA, sol, vmat, wmat);
|
|
|
compute_pos( &npsf, &convflag, &npsfflag,radmin2,radmax2,
|
compute_pos( &npsf, &convflag, &npsfflag,radmin2,radmax2,
|
r2, sol,flux, deltax, deltay,&dx,&dy);
|
r2, sol,flux, deltax, deltay,&dx,&dy);
|
}
|
}
|
|
|
/*-- Compute variances and covariances */
|
/*-- Compute variances and covariances */
|
svdvar(vmat, wmat, npsf*PSF_NA, covmat);
|
svdvar(vmat, wmat, npsf*PSF_NA, covmat);
|
var = covmat;
|
var = covmat;
|
for (j=0; j<npsf; j++, var += (npsf*PSF_NA+1)*PSF_NA)
|
for (j=0; j<npsf; j++, var += (npsf*PSF_NA+1)*PSF_NA)
|
{
|
{
|
/*---- First, the error on the flux estimate */
|
/*---- First, the error on the flux estimate */
|
fluxerr[j] = sqrt(*var)>0.0? sqrt(*var):999999.0;
|
fluxerr[j] = sqrt(*var)>0.0? sqrt(*var):999999.0;
|
//if (flux[j]<12*fluxerr && j>0)
|
//if (flux[j]<12*fluxerr && j>0)
|
// npsfmax--,flux[j]=0;
|
// npsfmax--,flux[j]=0;
|
if (flux[j]<12*fluxerr[j] && j>0)
|
if (flux[j]<12*fluxerr[j] && j>0)
|
{
|
{
|
flux[j]=0,kill++,npsfmax--;
|
flux[j]=0,kill++,npsfmax--;
|
//if(j==npsfmax-1)
|
//if(j==npsfmax-1)
|
// kill++;
|
// kill++;
|
}
|
}
|
}
|
}
|
if (npsfflag)
|
if (npsfflag)
|
{
|
{
|
/*--- If we reach this point we know the data are worth backuping */
|
/*--- If we reach this point we know the data are worth backuping */
|
for (j=0; j<npsf; j++)
|
for (j=0; j<npsf; j++)
|
{
|
{
|
deltaxb[j] = deltax[j];
|
deltaxb[j] = deltax[j];
|
deltayb[j] = deltay[j];
|
deltayb[j] = deltay[j];
|
fluxb[j] = flux[j];
|
fluxb[j] = flux[j];
|
fluxerrb[j]=fluxerr[j];
|
fluxerrb[j]=fluxerr[j];
|
}
|
}
|
}
|
}
|
}
|
}
|
npsf=npsf-1-kill;
|
npsf=npsf-1-kill;
|
|
|
/* Now keep only fitted stars that fall within the current detection area */
|
/* Now keep only fitted stars that fall within the current detection area */
|
i = 0;
|
i = 0;
|
for (j=0; j<npsf; j++)
|
for (j=0; j<npsf; j++)
|
{
|
{
|
x = (int)(deltaxb[j]+0.4999)+width/2;
|
x = (int)(deltaxb[j]+0.4999)+width/2;
|
y = (int)(deltayb[j]+0.4999)+height/2;
|
y = (int)(deltayb[j]+0.4999)+height/2;
|
if (x<0 || x>=width || y<0 || y>=height)
|
if (x<0 || x>=width || y<0 || y>=height)
|
continue;
|
continue;
|
if (weight[y*width+x] < 1/BIG)
|
if (weight[y*width+x] < 1/BIG)
|
continue;
|
continue;
|
if (10*fluxb[j]<fluxb[0] )
|
if (10*fluxb[j]<fluxb[0] )
|
continue;
|
continue;
|
if (fluxb[j]<=0 )
|
if (fluxb[j]<=0 )
|
continue;
|
continue;
|
|
|
if (FLAG(obj2.poserrmx2_psf))
|
if (FLAG(obj2.poserrmx2_psf))
|
compute_poserr(j,covmat,sol,obj2,x2,y2,xy, npsf);
|
compute_poserr(j,covmat,sol,obj2,x2,y2,xy, npsf);
|
|
|
deltax[i] = deltaxb[j];
|
deltax[i] = deltaxb[j];
|
deltay[i] = deltayb[j];
|
deltay[i] = deltayb[j];
|
flux[i] = fluxb[j];
|
flux[i] = fluxb[j];
|
fluxerr[i++] = fluxerrb[j];
|
fluxerr[i++] = fluxerrb[j];
|
}
|
}
|
|
|
npsf = i;
|
npsf = i;
|
|
|
/* Compute chi2 if asked to
|
/* Compute chi2 if asked to
|
if (FLAG(obj2.chi2_psf))
|
if (FLAG(obj2.chi2_psf))
|
{
|
{
|
for (j=0; j<npsf; j++)
|
for (j=0; j<npsf; j++)
|
{
|
{
|
chi2 = 0.0;
|
chi2 = 0.0;
|
for (d=data,w=weight,p=0; p<npix; w++,p++)
|
for (d=data,w=weight,p=0; p<npix; w++,p++)
|
{
|
{
|
pix = *(d++);
|
pix = *(d++);
|
pix -= psfmasks[j][p]*flux[j]**w;
|
pix -= psfmasks[j][p]*flux[j]**w;
|
chi2 += pix*pix;
|
chi2 += pix*pix;
|
if (chi2>1E29) chi2=1E28;
|
if (chi2>1E29) chi2=1E28;
|
}
|
}
|
obj2->chi2_psf = obj->sigbkg>0.?
|
obj2->chi2_psf = obj2->sigbkg>0.?
|
chi2/((npix - 3*npsf)*obj->sigbkg*obj->sigbkg):999999;
|
chi2/((npix - 3*npsf)*obj2->sigbkg*obj2->sigbkg):999999;
|
|
|
}
|
}
|
|
|
}*/
|
}*/
|
/* Compute relative chi2 if asked to */
|
/* Compute relative chi2 if asked to */
|
if (FLAG(obj2.chi2_psf))
|
if (FLAG(obj2.chi2_psf))
|
{
|
{
|
for (j=0; j<npsf; j++)
|
for (j=0; j<npsf; j++)
|
{
|
{
|
chi2 = 0.0;
|
chi2 = 0.0;
|
for (d=data,w=weight,p=0; p<npix; w++,p++)
|
for (d=data,w=weight,p=0; p<npix; w++,p++)
|
{
|
{
|
pix = *(d++)/flux[j];
|
pix = *(d++)/flux[j];
|
pix -= psfmasks[j][p]**w;
|
pix -= psfmasks[j][p]**w;
|
chi2 += pix*pix;
|
chi2 += pix*pix;
|
if (chi2>1E29) chi2=1E28;
|
if (chi2>1E29) chi2=1E28;
|
}
|
}
|
obj2->chi2_psf = flux[j]>0?
|
obj2->chi2_psf = flux[j]>0?
|
chi2/((npix - 3*npsf)*obj->sigbkg*obj->sigbkg):999999;
|
chi2/((npix - 3*npsf)*obj2->sigbkg*obj2->sigbkg):999999;
|
|
|
}
|
}
|
|
|
}
|
}
|
/* CHECK images */
|
/* CHECK images */
|
if (prefs.check[CHECK_SUBPSFPROTOS] || prefs.check[CHECK_PSFPROTOS])
|
if (prefs.check[CHECK_SUBPSFPROTOS] || prefs.check[CHECK_PSFPROTOS])
|
for (j=0; j<npsf; j++)
|
for (j=0; j<npsf; j++)
|
{
|
{
|
vignet_resample(psf->maskloc, psf->masksize[0], psf->masksize[1],
|
vignet_resample(psf->maskloc, psf->masksize[0], psf->masksize[1],
|
checkmask, pwidth, pheight,
|
checkmask, pwidth, pheight,
|
-deltax[j]*pixstep, -deltay[j]*pixstep, pixstep);
|
-deltax[j]*pixstep, -deltay[j]*pixstep, pixstep);
|
if ((check = prefs.check[CHECK_SUBPSFPROTOS]))
|
if ((check = prefs.check[CHECK_SUBPSFPROTOS]))
|
addcheck(check, checkmask, pwidth,pheight, ix,iy,-flux[j]);
|
addcheck(check, checkmask, pwidth,pheight, ix,iy,-flux[j]);
|
if ((check = prefs.check[CHECK_PSFPROTOS]))
|
if ((check = prefs.check[CHECK_PSFPROTOS]))
|
addcheck(check, checkmask, pwidth,pheight, ix,iy,flux[j]);
|
addcheck(check, checkmask, pwidth,pheight, ix,iy,flux[j]);
|
}
|
}
|
|
|
thepsfit->niter = niter;
|
thepsfit->niter = niter;
|
thepsfit->npsf = npsf;
|
thepsfit->npsf = npsf;
|
for (j=0; j<npsf; j++)
|
for (j=0; j<npsf; j++)
|
{
|
{
|
thepsfit->x[j] = ix+deltax[j]+1.0;
|
thepsfit->x[j] = ix+deltax[j]+1.0;
|
thepsfit->y[j] = iy+deltay[j]+1.0;
|
thepsfit->y[j] = iy+deltay[j]+1.0;
|
thepsfit->flux[j] = flux[j];
|
thepsfit->flux[j] = flux[j];
|
thepsfit->fluxerr[j] = fluxerr[j];
|
thepsfit->fluxerr[j] = fluxerr[j];
|
}
|
}
|
|
|
|
|
|
|
|
|
/* Now the morphology stuff */
|
|
if (prefs.pc_flag)
|
|
{
|
|
width = pwidth-1;
|
|
height = pheight-1;
|
|
npix = width*height;
|
|
|
|
/*-- Re-compute weights */
|
|
if (wfield)
|
|
{
|
|
psf_copyobjpix(datah, weighth, width, height, ix,iy, obj2,
|
|
!(field->flags&MEASURE_FIELD));
|
|
for (wh=weighth ,w=weight, p=npix; p--;)
|
|
*(w++) = (pix=*(wh++))<wthresh? sqrt(pix): 0.0;
|
|
}
|
|
else
|
|
{
|
|
psf_copyobjpix(datah, NULL, width, height, ix,iy, obj2,
|
|
!(field->flags&MEASURE_FIELD));
|
|
for (w=weight, dh=datah, p=npix; p--;)
|
|
*(w++) = ((pix = *(dh++))>-BIG && pix<satlevel)?
|
|
1.0/sqrt(backnoise2+(pix>0.0?pix/gain:0.0))
|
|
:0.0;
|
|
}
|
|
|
|
/*-- Weight the data */
|
|
dh = datah;
|
|
d = data;
|
|
w = weight;
|
|
for (p=npix; p--;)
|
|
*(d++) = *(dh++)*(*(w++));
|
|
|
|
pc_fit(psf, data, weight, width, height, ix,iy, dx,dy, npix,
|
|
field->backsig);
|
|
}
|
|
|
|
|
|
for (i=0; i<prefs.psf_npsfmax; i++)
|
for (i=0; i<prefs.psf_npsfmax; i++)
|
{
|
{
|
QFREE(psfmasks[i]);
|
QFREE(psfmasks[i]);
|
QFREE(psfmaskx[i]);
|
QFREE(psfmaskx[i]);
|
QFREE(psfmasky[i]);
|
QFREE(psfmasky[i]);
|
}
|
}
|
|
|
QFREE(psfmasks);
|
QFREE(psfmasks);
|
QFREE(psfmaskx);
|
QFREE(psfmaskx);
|
QFREE(psfmasky);
|
QFREE(psfmasky);
|
QFREE(datah);
|
QFREE(datah);
|
QFREE(data);
|
QFREE(data);
|
QFREE(data2);
|
QFREE(data2);
|
QFREE(data3);
|
QFREE(data3);
|
QFREE(weighth);
|
QFREE(weighth);
|
QFREE(weight);
|
QFREE(weight);
|
QFREE(data);
|
QFREE(data);
|
QFREE(mat);
|
QFREE(mat);
|
|
|
if (prefs.check[CHECK_SUBPSFPROTOS] || prefs.check[CHECK_PSFPROTOS]
|
if (prefs.check[CHECK_SUBPSFPROTOS] || prefs.check[CHECK_PSFPROTOS])
|
|| prefs.check[CHECK_SUBPCPROTOS] || prefs.check[CHECK_PCPROTOS]
|
|
|| prefs.check[CHECK_PCOPROTOS])
|
|
{
|
|
QFREE(checkmask);
|
QFREE(checkmask);
|
}
|
|
|
|
return;
|
return;
|
}
|
}
|
|
|
|
|
/****************************** double_psf_fit ******************************/
|
/****************************** double_psf_fit ******************************/
|
/* double fit to make psf detection on one image and photometry on another */
|
/* double fit to make psf detection on one image and photometry on another */
|
/****************************************************************************/
|
/****************************************************************************/
|
|
|
void double_psf_fit(psfstruct *ppsf, picstruct *pfield, picstruct *pwfield,
|
void double_psf_fit(psfstruct *ppsf, picstruct *pfield, picstruct *pwfield,
|
objstruct *obj, obj2struct *obj2,
|
obj2struct *obj2,
|
psfstruct *psf, picstruct *field, picstruct *wfield)
|
psfstruct *psf, picstruct *field, picstruct *wfield)
|
{
|
{
|
static double /* sum[PSF_NPSFMAX]*/ pdeltax[PSF_NPSFMAX],
|
static double /* sum[PSF_NPSFMAX]*/ pdeltax[PSF_NPSFMAX],
|
pdeltay[PSF_NPSFMAX],psol[PSF_NPSFMAX], pcovmat[PSF_NPSFMAX*PSF_NPSFMAX],
|
pdeltay[PSF_NPSFMAX],psol[PSF_NPSFMAX], pcovmat[PSF_NPSFMAX*PSF_NPSFMAX],
|
pvmat[PSF_NPSFMAX*PSF_NPSFMAX], pwmat[PSF_NPSFMAX],pflux[PSF_NPSFMAX],
|
pvmat[PSF_NPSFMAX*PSF_NPSFMAX], pwmat[PSF_NPSFMAX],pflux[PSF_NPSFMAX],
|
pfluxerr[PSF_NPSFMAX];
|
pfluxerr[PSF_NPSFMAX];
|
|
|
double *pmat,
|
double *pmat,
|
*pm, /* *pps, *px, *py,*/
|
*pm, /* *pps, *px, *py,*/
|
dx,dy,pdx,pdy, /* x1,y1, mx,my,mflux, */
|
dx,dy,pdx,pdy, /* x1,y1, mx,my,mflux, */
|
val, ppix,ppix2, /* dflux, */
|
val, ppix,ppix2, /* dflux, */
|
gain, radmin2,radmax2,satlevel
|
gain, radmin2,radmax2,satlevel
|
,chi2,pwthresh,pbacknoise2, /* mr, */
|
,chi2,pwthresh,pbacknoise2, /* mr, */
|
r2=0, psf_fwhm,ppsf_fwhm ;
|
r2=0, psf_fwhm,ppsf_fwhm ;
|
float **ppsfmasks, **ppsfmaskx,**ppsfmasky, *pps;
|
float **ppsfmasks, **ppsfmaskx,**ppsfmasky, *pps;
|
float *pdata, *pdata2, *pdata3, *pweight, *pd, *pw,
|
float *pdata, *pdata2, *pdata3, *pweight, *pd, *pw,
|
*pdh, *pwh, pixstep,ppixstep;
|
*pdh, *pwh, pixstep,ppixstep;
|
PIXTYPE *pdatah, *pweighth;
|
PIXTYPE *pdatah, *pweighth;
|
int i,j,k,p, npsf, npix,ix,iy,
|
int i,j,k,p, npsf, npix,ix,iy,
|
width, height, /* hw,hh, */
|
width, height, /* hw,hh, */
|
x,y, /* yb, */
|
x,y, /* yb, */
|
wbad, gainflag,
|
wbad, gainflag,
|
ival,npsfmax;
|
ival,npsfmax;
|
double *pvar;
|
double *pvar;
|
|
|
pdx = pdy =dx = dy = 0.0;
|
pdx = pdy =dx = dy = 0.0;
|
ppixstep = 1.0/ppsf->pixstep;
|
ppixstep = 1.0/ppsf->pixstep;
|
pixstep = 1.0/psf->pixstep;
|
pixstep = 1.0/psf->pixstep;
|
gain = (field->gain >0.0? field->gain: 1e30);
|
gain = (field->gain >0.0? field->gain: 1e30);
|
npsfmax=prefs.psf_npsfmax;
|
npsfmax=prefs.psf_npsfmax;
|
pbacknoise2 = pfield->backsig*pfield->backsig;
|
pbacknoise2 = pfield->backsig*pfield->backsig;
|
satlevel = field->satur_level - obj->bkg;
|
satlevel = field->satur_level - obj2->bkg;
|
gainflag = prefs.weightgain_flag;
|
gainflag = prefs.weightgain_flag;
|
psf_fwhm = psf->fwhm*psf->pixstep;
|
psf_fwhm = psf->fwhm*psf->pixstep;
|
ppsf_fwhm = ppsf->fwhm*ppsf->pixstep;
|
ppsf_fwhm = ppsf->fwhm*ppsf->pixstep;
|
pwthresh = pwfield?pwfield->weight_thresh:BIG;
|
pwthresh = pwfield?pwfield->weight_thresh:BIG;
|
|
|
/* Initialize outputs */
|
/* Initialize outputs */
|
ppsfit->niter = 0;
|
ppsfit->niter = 0;
|
ppsfit->npsf = 0;
|
ppsfit->npsf = 0;
|
for (j=0; j<npsfmax; j++)
|
for (j=0; j<npsfmax; j++)
|
{
|
{
|
ppsfit->x[j] = 999999.0;
|
ppsfit->x[j] = 999999.0;
|
ppsfit->y[j] = 999999.0;
|
ppsfit->y[j] = 999999.0;
|
ppsfit->flux[j] = 0.0;
|
ppsfit->flux[j] = 0.0;
|
ppsfit->fluxerr[j] = 0.0;
|
ppsfit->fluxerr[j] = 0.0;
|
pdeltax[j]= pdeltay[j]=psol[j]= pwmat[j]=pflux[j]=pfluxerr[j]=0.0;
|
pdeltax[j]= pdeltay[j]=psol[j]= pwmat[j]=pflux[j]=pfluxerr[j]=0.0;
|
|
|
}
|
}
|
|
|
ix = (obj->xmax+obj->xmin+1)/2;
|
ix = (obj2->xmax+obj2->xmin+1)/2;
|
iy = (obj->ymax+obj->ymin+1)/2;
|
iy = (obj2->ymax+obj2->ymin+1)/2;
|
width = obj->xmax-obj->xmin+1+psf_fwhm;
|
width = obj2->xmax-obj2->xmin+1+psf_fwhm;
|
if (width < (ival=(int)(psf_fwhm*2)))
|
if (width < (ival=(int)(psf_fwhm*2)))
|
width = ival;
|
width = ival;
|
height = obj->ymax-obj->ymin+1+psf_fwhm;
|
height = obj2->ymax-obj2->ymin+1+psf_fwhm;
|
if (height < (ival=(int)(psf_fwhm*2)))
|
if (height < (ival=(int)(psf_fwhm*2)))
|
height = ival;
|
height = ival;
|
npix = width*height;
|
npix = width*height;
|
radmin2 = PSF_MINSHIFT*PSF_MINSHIFT;
|
radmin2 = PSF_MINSHIFT*PSF_MINSHIFT;
|
radmax2 = npix/2.0;
|
radmax2 = npix/2.0;
|
psf_fit(psf, field, wfield, obj, obj2);
|
psf_fit(psf, field, wfield, obj2);
|
npsf=thepsfit->npsf;
|
npsf=thepsfit->npsf;
|
|
|
QMALLOC(ppsfmasks,float *,npsfmax);
|
QMALLOC(ppsfmasks,float *,npsfmax);
|
QMALLOC(ppsfmaskx,float *,npsfmax);
|
QMALLOC(ppsfmaskx,float *,npsfmax);
|
QMALLOC(ppsfmasky,float *,npsfmax);
|
QMALLOC(ppsfmasky,float *,npsfmax);
|
|
|
for (i=0; i<npsfmax; i++)
|
for (i=0; i<npsfmax; i++)
|
{
|
{
|
QMALLOC(ppsfmasks[i],float,npix);
|
QMALLOC(ppsfmasks[i],float,npix);
|
QMALLOC(ppsfmaskx[i],float,npix);
|
QMALLOC(ppsfmaskx[i],float,npix);
|
QMALLOC(ppsfmasky[i],float,npix);
|
QMALLOC(ppsfmasky[i],float,npix);
|
}
|
}
|
|
|
QMALLOC(pweighth, PIXTYPE, npix);
|
QMALLOC(pweighth, PIXTYPE, npix);
|
QMALLOC(pweight, float, npix);
|
QMALLOC(pweight, float, npix);
|
QMALLOC(pdatah, PIXTYPE, npix);
|
QMALLOC(pdatah, PIXTYPE, npix);
|
QMALLOC(pdata, float, npix);
|
QMALLOC(pdata, float, npix);
|
QMALLOC(pdata2, float, npix);
|
QMALLOC(pdata2, float, npix);
|
QMALLOC(pdata3, float, npix);
|
QMALLOC(pdata3, float, npix);
|
QMALLOC(pmat, double, npix*npsfmax);
|
QMALLOC(pmat, double, npix*npsfmax);
|
|
|
for (j=0; j<npsf; j++)
|
for (j=0; j<npsf; j++)
|
{
|
{
|
pdeltax[j] =thepsfit->x[j]-ix-1 ;
|
pdeltax[j] =thepsfit->x[j]-ix-1 ;
|
pdeltay[j] =thepsfit->y[j]-iy-1 ;
|
pdeltay[j] =thepsfit->y[j]-iy-1 ;
|
ppsfit->flux[j] = 0;
|
ppsfit->flux[j] = 0;
|
ppsfit->fluxerr[j] = 0;
|
ppsfit->fluxerr[j] = 0;
|
}
|
}
|
|
|
/*------------------- Now the photometry fit ---------------------*/
|
/*------------------- Now the photometry fit ---------------------*/
|
/* Compute photometry weights */
|
/* Compute photometry weights */
|
wbad = 0;
|
wbad = 0;
|
if (pwfield)
|
if (pwfield)
|
{
|
{
|
psf_copyobjpix(pdatah, pweighth, width, height, ix,iy, obj2, 0);
|
psf_copyobjpix(pdatah, pweighth, width, height, ix,iy, obj2, 0);
|
for (pwh=pweighth, pw=pweight, pdh=pdatah,p=npix; p--;)
|
for (pwh=pweighth, pw=pweight, pdh=pdatah,p=npix; p--;)
|
{
|
{
|
if ((ppix=*(pwh++)) < pwthresh && ppix>0
|
if ((ppix=*(pwh++)) < pwthresh && ppix>0
|
&& (ppix2=*(pdh++))>-BIG && ppix2<satlevel)
|
&& (ppix2=*(pdh++))>-BIG && ppix2<satlevel)
|
*(pw++) = 1/sqrt(ppix+(ppix2>0.0?
|
*(pw++) = 1/sqrt(ppix+(ppix2>0.0?
|
(gainflag? ppix2*ppix/pbacknoise2:ppix2)/gain : 0.0));
|
(gainflag? ppix2*ppix/pbacknoise2:ppix2)/gain : 0.0));
|
else
|
else
|
{
|
{
|
*(pw++) = 0.0;
|
*(pw++) = 0.0;
|
wbad++;
|
wbad++;
|
}
|
}
|
}
|
}
|
}
|
}
|
else
|
else
|
{
|
{
|
psf_copyobjpix(pdatah, NULL, width, height, ix,iy, obj2, 0);
|
psf_copyobjpix(pdatah, NULL, width, height, ix,iy, obj2, 0);
|
for (pw=pweight, pdh=pdatah, p=npix; p--;)
|
for (pw=pweight, pdh=pdatah, p=npix; p--;)
|
if ((ppix=*(pdh++))>-BIG && ppix<satlevel)
|
if ((ppix=*(pdh++))>-BIG && ppix<satlevel)
|
*(pw++) = 1.0/sqrt(pbacknoise2+(ppix>0.0?ppix/gain:0.0));
|
*(pw++) = 1.0/sqrt(pbacknoise2+(ppix>0.0?ppix/gain:0.0));
|
else
|
else
|
{
|
{
|
*(pw++) = 0.0;
|
*(pw++) = 0.0;
|
wbad++;
|
wbad++;
|
}
|
}
|
}
|
}
|
|
|
/* Special action if most of the weights are zero!! */
|
/* Special action if most of the weights are zero!! */
|
if (wbad>=npix-3)
|
if (wbad>=npix-3)
|
return;
|
return;
|
|
|
/* Weight the data */
|
/* Weight the data */
|
pdh = pdatah;
|
pdh = pdatah;
|
pd = pdata;
|
pd = pdata;
|
pw = pweight;
|
pw = pweight;
|
val = obj->dbkg;
|
val = obj2->dbkg;
|
for (p=npix; p--;)
|
for (p=npix; p--;)
|
*(pd++) = (*(pdh++)-val)**(pw++);
|
*(pd++) = (*(pdh++)-val)**(pw++);
|
|
|
|
|
/* Get the photmetry PSF */
|
/* Get the photmetry PSF */
|
psf_build(ppsf);
|
psf_build(ppsf, obj2);
|
for (j=1; j<=npsf; j++)
|
for (j=1; j<=npsf; j++)
|
{
|
{
|
if (j>1)
|
if (j>1)
|
{
|
{
|
/*---- Subtract //previously fitted components in photometry image */
|
/*---- Subtract //previously fitted components in photometry image */
|
pd = pdata2;
|
pd = pdata2;
|
pdh = pdatah;
|
pdh = pdatah;
|
for (p=npix; p--;)
|
for (p=npix; p--;)
|
*(pd++) = (double)*(pdh++);
|
*(pd++) = (double)*(pdh++);
|
for (k=0; k<j-1; k++)
|
for (k=0; k<j-1; k++)
|
{
|
{
|
pd = pdata2;
|
pd = pdata2;
|
pps = ppsfmasks[k];
|
pps = ppsfmasks[k];
|
for (p=npix; p--;)
|
for (p=npix; p--;)
|
*(pd++) -= pflux[k]**(pps++);
|
*(pd++) -= pflux[k]**(pps++);
|
}
|
}
|
convolve_image(pfield, pdata2, pdata3, width,height);
|
convolve_image(pfield, pdata2, pdata3, width,height);
|
/*---- Ignore regions too close to stellar cores */
|
/*---- Ignore regions too close to stellar cores */
|
for (k=0; k<j-1; k++)
|
for (k=0; k<j-1; k++)
|
{
|
{
|
pd = pdata3;
|
pd = pdata3;
|
dy = -((double)(height/2)+pdeltay[k]);
|
dy = -((double)(height/2)+pdeltay[k]);
|
for (y=height; y--; dy += 1.0)
|
for (y=height; y--; dy += 1.0)
|
{
|
{
|
dx = -((double)(width/2)+pdeltax[k]);
|
dx = -((double)(width/2)+pdeltax[k]);
|
for (x=width; x--; dx+= 1.0, pd++)
|
for (x=width; x--; dx+= 1.0, pd++)
|
if (dx*dx+dy*dy<r2) /*?*/
|
if (dx*dx+dy*dy<r2) /*?*/
|
*pd = -BIG;
|
*pd = -BIG;
|
}
|
}
|
}
|
}
|
}
|
}
|
|
|
pm=pmat;
|
pm=pmat;
|
for (k=0; k<j; k++)
|
for (k=0; k<j; k++)
|
{
|
{
|
/*------ Resample the PSFs here for the 1st iteration */
|
/*------ Resample the PSFs here for the 1st iteration */
|
vignet_resample(ppsf->maskloc,
|
vignet_resample(ppsf->maskloc,
|
ppsf->masksize[0], ppsf->masksize[1],
|
ppsf->masksize[0], ppsf->masksize[1],
|
ppsfmasks[k], width, height,
|
ppsfmasks[k], width, height,
|
-pdeltax[k]*ppixstep, -pdeltay[k]*ppixstep,
|
-pdeltax[k]*ppixstep, -pdeltay[k]*ppixstep,
|
ppixstep);
|
ppixstep);
|
pm=compute_gradient_phot(pweight,width,height, ppsfmasks[k],pm);
|
pm=compute_gradient_phot(pweight,width,height, ppsfmasks[k],pm);
|
}
|
}
|
|
|
svdfit(pmat, pdata, npix, j, psol, pvmat, pwmat);
|
svdfit(pmat, pdata, npix, j, psol, pvmat, pwmat);
|
compute_pos_phot( &j, psol,pflux);
|
compute_pos_phot( &j, psol,pflux);
|
|
|
for (k=0; k<j; k++)
|
for (k=0; k<j; k++)
|
{
|
{
|
svdvar(pvmat, pwmat, j, pcovmat);
|
svdvar(pvmat, pwmat, j, pcovmat);
|
pvar = pcovmat;
|
pvar = pcovmat;
|
pfluxerr[k]= sqrt(*pvar)>0.0 && sqrt(*pvar)<99?
|
pfluxerr[k]= sqrt(*pvar)>0.0 && sqrt(*pvar)<99?
|
sqrt(*pvar):99;
|
sqrt(*pvar):99;
|
}
|
}
|
}
|
}
|
/* Compute chi2 if asked to
|
/* Compute chi2 if asked to
|
if (FLAG(obj2.chi2_psf))
|
if (FLAG(obj2.chi2_psf))
|
{
|
{
|
for (j=0; j<npsf; j++)
|
for (j=0; j<npsf; j++)
|
{
|
{
|
chi2 = 0.0;
|
chi2 = 0.0;
|
for (pd=pdata,pw=pweight,p=0; p<npix; pw++,p++)
|
for (pd=pdata,pw=pweight,p=0; p<npix; pw++,p++)
|
{
|
{
|
ppix = *(pd++);
|
ppix = *(pd++);
|
ppix -= ppsfmasks[j][p]*pflux[j]**pw;
|
ppix -= ppsfmasks[j][p]*pflux[j]**pw;
|
chi2 += ppix*ppix;
|
chi2 += ppix*ppix;
|
if (chi2>1E29) chi2=1E28;
|
if (chi2>1E29) chi2=1E28;
|
}
|
}
|
obj2->chi2_psf = obj->sigbkg>0.?
|
obj2->chi2_psf = obj2->sigbkg>0.?
|
chi2/((npix - 3*npsf)*obj->sigbkg*obj->sigbkg):999999;
|
chi2/((npix - 3*npsf)*obj2->sigbkg*obj2->sigbkg):999999;
|
|
|
}
|
}
|
|
|
}
|
}
|
*/
|
*/
|
/* Compute relative error if asked to */
|
/* Compute relative error if asked to */
|
if (FLAG(obj2.chi2_psf))
|
if (FLAG(obj2.chi2_psf))
|
{
|
{
|
for (j=0; j<npsf; j++)
|
for (j=0; j<npsf; j++)
|
{
|
{
|
chi2 = 0.0;
|
chi2 = 0.0;
|
for (pd=pdata,pw=pweight,p=0; p<npix; pw++,p++)
|
for (pd=pdata,pw=pweight,p=0; p<npix; pw++,p++)
|
{
|
{
|
ppix = *(pd++)/pflux[j];
|
ppix = *(pd++)/pflux[j];
|
ppix -= ppsfmasks[j][p]**pw;
|
ppix -= ppsfmasks[j][p]**pw;
|
chi2 += ppix*ppix;
|
chi2 += ppix*ppix;
|
if (chi2>1E29) chi2=1E28;
|
if (chi2>1E29) chi2=1E28;
|
}
|
}
|
obj2->chi2_psf = pflux[j]>0?
|
obj2->chi2_psf = pflux[j]>0?
|
chi2/((npix - 3*npsf)*obj->sigbkg*obj->sigbkg):999999;
|
chi2/((npix - 3*npsf)*obj2->sigbkg*obj2->sigbkg):999999;
|
|
|
}
|
}
|
|
|
}
|
}
|
ppsfit->niter = thepsfit->niter;
|
ppsfit->niter = thepsfit->niter;
|
ppsfit->npsf = npsf;
|
ppsfit->npsf = npsf;
|
|
|
for (j=0; j<npsf; j++)
|
for (j=0; j<npsf; j++)
|
{
|
{
|
thepsfit->x[j] = ix+pdeltax[j]+1.0;
|
thepsfit->x[j] = ix+pdeltax[j]+1.0;
|
thepsfit->y[j] = iy+pdeltay[j]+1.0;
|
thepsfit->y[j] = iy+pdeltay[j]+1.0;
|
thepsfit->flux[j] = pflux[j];
|
thepsfit->flux[j] = pflux[j];
|
thepsfit->fluxerr[j] = pfluxerr[j];
|
thepsfit->fluxerr[j] = pfluxerr[j];
|
ppsfit->x[j] = ix+pdeltax[j]+1.0;
|
ppsfit->x[j] = ix+pdeltax[j]+1.0;
|
ppsfit->y[j] = iy+pdeltay[j]+1.0;
|
ppsfit->y[j] = iy+pdeltay[j]+1.0;
|
ppsfit->flux[j] = pflux[j];
|
ppsfit->flux[j] = pflux[j];
|
ppsfit->fluxerr[j] = pfluxerr[j];
|
ppsfit->fluxerr[j] = pfluxerr[j];
|
}
|
}
|
|
|
|
|
for (i=0; i<npsfmax; i++)
|
for (i=0; i<npsfmax; i++)
|
{
|
{
|
QFREE(ppsfmasks[i]);
|
QFREE(ppsfmasks[i]);
|
QFREE(ppsfmaskx[i]);
|
QFREE(ppsfmaskx[i]);
|
QFREE(ppsfmasky[i]);
|
QFREE(ppsfmasky[i]);
|
}
|
}
|
|
|
|
|
QFREE(ppsfmasks);
|
QFREE(ppsfmasks);
|
QFREE(ppsfmaskx);
|
QFREE(ppsfmaskx);
|
QFREE(ppsfmasky);
|
QFREE(ppsfmasky);
|
QFREE(pdatah);
|
QFREE(pdatah);
|
QFREE(pdata);
|
QFREE(pdata);
|
QFREE(pdata2);
|
QFREE(pdata2);
|
QFREE(pdata3);
|
QFREE(pdata3);
|
QFREE(pweighth);
|
QFREE(pweighth);
|
QFREE(pweight);
|
QFREE(pweight);
|
QFREE(pdata);
|
QFREE(pdata);
|
QFREE(pmat);
|
QFREE(pmat);
|
|
|
if (prefs.check[CHECK_SUBPSFPROTOS] || prefs.check[CHECK_PSFPROTOS]
|
if (prefs.check[CHECK_SUBPSFPROTOS] || prefs.check[CHECK_PSFPROTOS])
|
|| prefs.check[CHECK_SUBPCPROTOS] || prefs.check[CHECK_PCPROTOS]
|
|
|| prefs.check[CHECK_PCOPROTOS])
|
|
{
|
{
|
QFREE(checkmask);
|
QFREE(checkmask);
|
}
|
}
|
return;
|
return;
|
}
|
}
|
|
|
|
|
/****** psf_copyobjpix ******************************************************
|
/****** psf_copyobjpix ******************************************************
|
PROTO int psf_copyobjpix(PIXTYPE *data, PIXTYPE *weight,
|
PROTO int psf_copyobjpix(PIXTYPE *data, PIXTYPE *weight,
|
int wout, int hout, int ix, int iy,
|
int wout, int hout, int ix, int iy,
|
obj2struct *obj2, int detect_flag);
|
obj2struct *obj2, int detect_flag);
|
PURPOSE Copy a piece of the input object image/weights to local arrays.
|
PURPOSE Copy a piece of the input object image/weights to local arrays.
|
INPUT Pointer to the output data array,
|
INPUT Pointer to the output data array,
|
pointer to the output weight array,
|
pointer to the output weight array,
|
output frame width,
|
output frame width,
|
output frame height,
|
output frame height,
|
integer x coordinate,
|
integer x coordinate,
|
integer y coordinate,
|
integer y coordinate,
|
pointer to the obj2 structure,
|
pointer to the obj2 structure,
|
detection field flag (non 0 only for pure detection images).
|
detection field flag (non 0 only for pure detection images).
|
OUTPUT RETURN_ERROR if the coordinates are outside object image,
|
OUTPUT RETURN_ERROR if the coordinates are outside object image,
|
RETURN_OK otherwise.
|
RETURN_OK otherwise.
|
NOTES -.
|
NOTES -.
|
AUTHOR E. Bertin (IAP)
|
AUTHOR E. Bertin (IAP)
|
VERSION 23/06/2011
|
VERSION 23/06/2011
|
***/
|
***/
|
int psf_copyobjpix(PIXTYPE *data, PIXTYPE *weight,
|
int psf_copyobjpix(PIXTYPE *data, PIXTYPE *weight,
|
int wout, int hout, int ix, int iy,
|
int wout, int hout, int ix, int iy,
|
obj2struct *obj2, int detect_flag)
|
obj2struct *obj2, int detect_flag)
|
{
|
{
|
PIXTYPE *datat,*weightt, *imagein,*weightin;
|
PIXTYPE *datat,*weightt, *imagein,*weightin;
|
int i,y, win,hin, w2, xmin,xmax,ymin,ymax;
|
int i,y, win,hin, w2, xmin,xmax,ymin,ymax;
|
|
|
/* Set output to -BIG */
|
/* Set output to -BIG */
|
datat = data;
|
datat = data;
|
for (i=wout*hout; i--;)
|
for (i=wout*hout; i--;)
|
*(datat++) = -BIG;
|
*(datat++) = -BIG;
|
if (weight)
|
if (weight)
|
memset(weight, 0, wout*hout*sizeof(PIXTYPE));
|
memset(weight, 0, wout*hout*sizeof(PIXTYPE));
|
|
|
ix -= obj2->immin[0];
|
ix -= obj2->immin[0];
|
iy -= obj2->immin[1];
|
iy -= obj2->immin[1];
|
win = obj2->imsize[0];
|
win = obj2->imsize[0];
|
hin = obj2->imsize[1];
|
hin = obj2->imsize[1];
|
|
|
/* Don't go further if out of frame!! */
|
/* Don't go further if out of frame!! */
|
if (ix<0 || ix>=win || iy<0 || iy>=hin)
|
if (ix<0 || ix>=win || iy<0 || iy>=hin)
|
return RETURN_ERROR;
|
return RETURN_ERROR;
|
|
|
/* Set the image boundaries */
|
/* Set the image boundaries */
|
w2 = wout;
|
w2 = wout;
|
ymin = iy-hout/2;
|
ymin = iy-hout/2;
|
ymax = ymin + hout;
|
ymax = ymin + hout;
|
if (ymin<0)
|
if (ymin<0)
|
{
|
{
|
data -= ymin*wout;
|
data -= ymin*wout;
|
if (weight)
|
if (weight)
|
weight -= ymin*wout;
|
weight -= ymin*wout;
|
ymin = 0;
|
ymin = 0;
|
}
|
}
|
if (ymax>hin)
|
if (ymax>hin)
|
ymax = hin;
|
ymax = hin;
|
|
|
xmin = ix-wout/2;
|
xmin = ix-wout/2;
|
xmax = xmin + wout;
|
xmax = xmin + wout;
|
if (xmax>win)
|
if (xmax>win)
|
{
|
{
|
w2 -= xmax-win;
|
w2 -= xmax-win;
|
xmax = win;
|
xmax = win;
|
}
|
}
|
if (xmin<0)
|
if (xmin<0)
|
{
|
{
|
data -= xmin;
|
data -= xmin;
|
if (weight)
|
if (weight)
|
weight -= xmin;
|
weight -= xmin;
|
w2 += xmin;
|
w2 += xmin;
|
xmin = 0;
|
xmin = 0;
|
}
|
}
|
|
|
/* Copy the right pixels to the destination */
|
/* Copy the right pixels to the destination */
|
imagein = detect_flag? obj2->dimage:obj2->image;
|
imagein = detect_flag? obj2->dimage:obj2->image;
|
for (y=ymin; y<ymax; y++, data += wout)
|
for (y=ymin; y<ymax; y++, data += wout)
|
memcpy(data, imagein + xmin+y*win, w2*sizeof(PIXTYPE));
|
memcpy(data, imagein + xmin+y*win, w2*sizeof(PIXTYPE));
|
if (weight)
|
if (weight)
|
{
|
{
|
weightin = detect_flag? obj2->dweight:obj2->weight;
|
weightin = detect_flag? obj2->dweight:obj2->weight;
|
for (y=ymin; y<ymax; y++, data += wout)
|
for (y=ymin; y<ymax; y++, data += wout)
|
memcpy(weight, weightin + xmin+y*win, w2*sizeof(PIXTYPE));
|
memcpy(weight, weightin + xmin+y*win, w2*sizeof(PIXTYPE));
|
}
|
}
|
|
|
|
|
return RETURN_OK;
|
return RETURN_OK;
|
}
|
}
|
|
|
|
|
|
|
/******************************* psf_build **********************************/
|
/******************************* psf_build **********************************/
|
/*
|
/*
|
Build the local PSF (function of "context").
|
Build the local PSF (function of "context").
|
*/
|
*/
|
void psf_build(psfstruct *psf)
|
void psf_build(psfstruct *psf, obj2struct *obj2)
|
{
|
{
|
static double pos[POLY_MAXDIM];
|
static double pos[POLY_MAXDIM];
|
double *basis, fac;
|
double *basis, fac;
|
float *ppc, *pl;
|
float *ppc, *pl;
|
int i,n,p, ndim, npix;
|
int i,n,p, ndim, npix;
|
|
|
if (psf->build_flag)
|
if (psf->build_flag)
|
return;
|
return;
|
|
|
npix = psf->masksize[0]*psf->masksize[1];
|
npix = psf->masksize[0]*psf->masksize[1];
|
|
|
/* Reset the Local PSF mask */
|
/* Reset the Local PSF mask */
|
memset(psf->maskloc, 0, npix*sizeof(float));
|
memset(psf->maskloc, 0, npix*sizeof(float));
|
|
|
/* Grab the context vector */
|
/* Grab the context vector */
|
ndim = psf->poly->ndim;
|
ndim = psf->poly->ndim;
|
for (i=0; i<ndim; i++)
|
for (i=0; i<ndim; i++)
|
{
|
{
|
ttypeconv(psf->context[i], &pos[i], psf->contexttyp[i],T_DOUBLE);
|
ttypeconv((char *)obj2 + ((void *)psf->context[i]-(void *)&flagobj2),
|
|
&pos[i], psf->contexttyp[i],T_DOUBLE);
|
pos[i] = (pos[i] - psf->contextoffset[i]) / psf->contextscale[i];
|
pos[i] = (pos[i] - psf->contextoffset[i]) / psf->contextscale[i];
|
}
|
}
|
poly_func(psf->poly, pos);
|
poly_func(psf->poly, pos);
|
|
|
basis = psf->poly->basis;
|
basis = psf->poly->basis;
|
|
|
ppc = psf->maskcomp;
|
ppc = psf->maskcomp;
|
/* Sum each component */
|
/* Sum each component */
|
for (n = (psf->maskdim>2?psf->masksize[2]:1); n--;)
|
for (n = (psf->maskdim>2?psf->masksize[2]:1); n--;)
|
{
|
{
|
pl = psf->maskloc;
|
pl = psf->maskloc;
|
fac = *(basis++);
|
fac = *(basis++);
|
for (p=npix; p--;)
|
for (p=npix; p--;)
|
*(pl++) += fac**(ppc++);
|
*(pl++) += fac**(ppc++);
|
}
|
}
|
|
|
psf->build_flag = 1;
|
psf->build_flag = 1;
|
|
|
return;
|
return;
|
}
|
}
|
|
|
|
|
/******************************** psf_fwhm **********************************/
|
/******************************** psf_fwhm **********************************/
|
/*
|
/*
|
Return the local PSF FWHM.
|
Return the local PSF FWHM.
|
*/
|
*/
|
double psf_fwhm(psfstruct *psf)
|
double psf_fwhm(psfstruct *psf, obj2struct *obj2)
|
{
|
{
|
float *pl,
|
float *pl,
|
val, max;
|
val, max;
|
int n,p, npix;
|
int n,p, npix;
|
|
|
if (!psf->build_flag)
|
if (!psf->build_flag)
|
psf_build(psf);
|
psf_build(psf, obj2);
|
|
|
npix = psf->masksize[0]*psf->masksize[1];
|
npix = psf->masksize[0]*psf->masksize[1];
|
max = -BIG;
|
max = -BIG;
|
pl = psf->maskloc;
|
pl = psf->maskloc;
|
for (p=npix; p--;)
|
for (p=npix; p--;)
|
if ((val=*(pl++)) > max)
|
if ((val=*(pl++)) > max)
|
max = val;
|
max = val;
|
pl = psf->maskloc;
|
pl = psf->maskloc;
|
max /= 2.0;
|
max /= 2.0;
|
n = 0;
|
n = 0;
|
for (p=npix; p--;)
|
for (p=npix; p--;)
|
if (*(pl++) >= max)
|
if (*(pl++) >= max)
|
n++;
|
n++;
|
|
|
return 2.0*sqrt(n/PI)*psf->pixstep;
|
return 2.0*sqrt(n/PI)*psf->pixstep;
|
}
|
}
|
|
|
|
|
/*****************************compute_gradient*********************************/
|
/*****************************compute_gradient*********************************/
|
|
|
double *compute_gradient(float *weight,int width, int height,
|
double *compute_gradient(float *weight,int width, int height,
|
float *masks,float *maskx,float *masky
|
float *masks,float *maskx,float *masky
|
,double *m)
|
,double *m)
|
{
|
{
|
int x,y;
|
int x,y;
|
float *w, *ps,*px,*py;
|
float *w, *ps,*px,*py;
|
|
|
/*------ copy of the (weighted) PSF, with outer ring set to zero */
|
/*------ copy of the (weighted) PSF, with outer ring set to zero */
|
ps = masks;
|
ps = masks;
|
w = weight;
|
w = weight;
|
for (y=0; y<height; y++)
|
for (y=0; y<height; y++)
|
for (x=0; x<width; x++, ps++, w++)
|
for (x=0; x<width; x++, ps++, w++)
|
*(m++) = y?(y>=(height-1)?0:(x?(x>=(width-1)?0:*ps**w):0)):0;
|
*(m++) = y?(y>=(height-1)?0:(x?(x>=(width-1)?0:*ps**w):0)):0;
|
/*------ (weighted) PSF gradient in x (kernel for first moment in x) */
|
/*------ (weighted) PSF gradient in x (kernel for first moment in x) */
|
ps = masks;
|
ps = masks;
|
px = maskx;
|
px = maskx;
|
w = weight;
|
w = weight;
|
for (y=0; y<height; y++)
|
for (y=0; y<height; y++)
|
for (x=0; x<width; x++, ps++, w++)
|
for (x=0; x<width; x++, ps++, w++)
|
*(m++) = ((*px++) = (x?(x>=(width-1)?0:*(ps+1)-*(ps-1)):0))**w/2;
|
*(m++) = ((*px++) = (x?(x>=(width-1)?0:*(ps+1)-*(ps-1)):0))**w/2;
|
/*------ (weighted) PSF gradient in y (kernel for first moment in y) */
|
/*------ (weighted) PSF gradient in y (kernel for first moment in y) */
|
ps = masks;
|
ps = masks;
|
py = masky;
|
py = masky;
|
w = weight;
|
w = weight;
|
for (y=0; y<height; y++)
|
for (y=0; y<height; y++)
|
for (x=0; x<width; x++, ps++, w++)
|
for (x=0; x<width; x++, ps++, w++)
|
*(m++) = (*(py++)=(y?(y>=(height-1)?0:*(ps+width)-*(ps-width)):0))
|
*(m++) = (*(py++)=(y?(y>=(height-1)?0:*(ps+width)-*(ps-width)):0))
|
**w/2;
|
**w/2;
|
return m;
|
return m;
|
}
|
}
|
|
|
/*****************************compute_gradient_phot*****************************
|
/*****************************compute_gradient_phot*****************************
|
****/
|
****/
|
|
|
double *compute_gradient_phot(float *pweight,int width, int height,
|
double *compute_gradient_phot(float *pweight,int width, int height,
|
float *pmasks,double *pm)
|
float *pmasks,double *pm)
|
|
|
{
|
{
|
int x,y;
|
int x,y;
|
float *pw, *pps;
|
float *pw, *pps;
|
|
|
/*------ copy of the (weighted) PSF, with outer ring set to zero */
|
/*------ copy of the (weighted) PSF, with outer ring set to zero */
|
pps = pmasks;
|
pps = pmasks;
|
pw = pweight;
|
pw = pweight;
|
for (y=0; y<height; y++)
|
for (y=0; y<height; y++)
|
for (x=0; x<width; x++, pps++, pw++)
|
for (x=0; x<width; x++, pps++, pw++)
|
*(pm++) = y?(y>=(height-1)?0:(x?(x>=(width-1)?0:*pps**pw):0)):0;
|
*(pm++) = y?(y>=(height-1)?0:(x?(x>=(width-1)?0:*pps**pw):0)):0;
|
|
|
return pm;
|
return pm;
|
}
|
}
|
|
|
/**************************compute_pos********************************/
|
/**************************compute_pos********************************/
|
|
|
void compute_pos(int *pnpsf,int *pconvflag,int *pnpsfflag,double radmin2,
|
void compute_pos(int *pnpsf,int *pconvflag,int *pnpsfflag,double radmin2,
|
double radmax2,double r2,double *sol,double *flux
|
double radmax2,double r2,double *sol,double *flux
|
,double *deltax,double *deltay,double *pdx,double *pdy)
|
,double *deltax,double *deltay,double *pdx,double *pdy)
|
{
|
{
|
int j,k,convflag,npsfflag,npsf;
|
int j,k,convflag,npsfflag,npsf;
|
double dx,dy;
|
double dx,dy;
|
|
|
dx=*pdx;
|
dx=*pdx;
|
dy=*pdy;
|
dy=*pdy;
|
convflag=*pconvflag;
|
convflag=*pconvflag;
|
npsfflag=*pnpsfflag;
|
npsfflag=*pnpsfflag;
|
npsf=*pnpsf;
|
npsf=*pnpsf;
|
for (j=0; j<npsf; j++)
|
for (j=0; j<npsf; j++)
|
{
|
{
|
flux[j] = sol[j*PSF_NA];
|
flux[j] = sol[j*PSF_NA];
|
/*------ Update the PSF shifts */
|
/*------ Update the PSF shifts */
|
if (fabs(flux[j])>0.0)
|
if (fabs(flux[j])>0.0)
|
{
|
{
|
dx = -sol[j*PSF_NA+1]/((npsf>1?2:1)*flux[j]);
|
dx = -sol[j*PSF_NA+1]/((npsf>1?2:1)*flux[j]);
|
dy = -sol[j*PSF_NA+2]/((npsf>1?2:1)*flux[j]);
|
dy = -sol[j*PSF_NA+2]/((npsf>1?2:1)*flux[j]);
|
}
|
}
|
|
|
deltax[j] += dx;
|
deltax[j] += dx;
|
deltay[j] += dy;
|
deltay[j] += dy;
|
/*------ Continue until all PSFs have come to a complete stop */
|
/*------ Continue until all PSFs have come to a complete stop */
|
if ((dx*dx+dy*dy) > radmin2)
|
if ((dx*dx+dy*dy) > radmin2)
|
convflag = 1;
|
convflag = 1;
|
}
|
}
|
for (j=0; j<npsf; j++)
|
for (j=0; j<npsf; j++)
|
{
|
{
|
/*------ Exit if too much decentering or negative flux */
|
/*------ Exit if too much decentering or negative flux */
|
for (k=j+1; k<npsf; k++)
|
for (k=j+1; k<npsf; k++)
|
{
|
{
|
dx = deltax[j]-deltax[k];
|
dx = deltax[j]-deltax[k];
|
dy = deltay[j]-deltay[k];
|
dy = deltay[j]-deltay[k];
|
if (dx*dx+dy*dy<r2/4.0)
|
if (dx*dx+dy*dy<r2/4.0)
|
{
|
{
|
flux[j] = -BIG;
|
flux[j] = -BIG;
|
break;
|
break;
|
}
|
}
|
}
|
}
|
if (flux[j]<0.0
|
if (flux[j]<0.0
|
|| (deltax[j]*deltax[j] + deltay[j]*deltay[j]) > radmax2)
|
|| (deltax[j]*deltax[j] + deltay[j]*deltay[j]) > radmax2)
|
{
|
{
|
npsfflag = 0;
|
npsfflag = 0;
|
convflag = 0;
|
convflag = 0;
|
npsf--;
|
npsf--;
|
break;
|
break;
|
}
|
}
|
}
|
}
|
*pdx=dx;
|
*pdx=dx;
|
*pdy=dy;
|
*pdy=dy;
|
*pconvflag=convflag;
|
*pconvflag=convflag;
|
*pnpsfflag= npsfflag;
|
*pnpsfflag= npsfflag;
|
*pnpsf=npsf;
|
*pnpsf=npsf;
|
return;
|
return;
|
}
|
}
|
|
|
/**************************compute_pos_phot********************************/
|
/**************************compute_pos_phot********************************/
|
|
|
void compute_pos_phot(int *pnpsf,double *sol,double *flux)
|
void compute_pos_phot(int *pnpsf,double *sol,double *flux)
|
{
|
{
|
int j,npsf;
|
int j,npsf;
|
npsf=*pnpsf;
|
npsf=*pnpsf;
|
for (j=0; j<npsf; j++)
|
for (j=0; j<npsf; j++)
|
{
|
{
|
flux[j] = sol[j];
|
flux[j] = sol[j];
|
}
|
}
|
*pnpsf=npsf;
|
*pnpsf=npsf;
|
return;
|
return;
|
}
|
}
|
|
|
|
|
/************************************compute_poserr*****************************
|
/************************************compute_poserr*****************************
|
*********/
|
*********/
|
|
|
void compute_poserr( int j,double *var,double *sol,obj2struct *obj2,double *x2,
|
void compute_poserr( int j,double *var,double *sol,obj2struct *obj2,double *x2,
|
double *y2,double *xy, int npsf)
|
double *y2,double *xy, int npsf)
|
{
|
{
|
double vara,covab,varb, f2;
|
double vara,covab,varb, f2;
|
|
|
/*------ Variances and covariance along x and y */
|
/*------ Variances and covariance along x and y */
|
vara = *(var += (PSF_NA*npsf+1)*(j*PSF_NA+1));
|
vara = *(var += (PSF_NA*npsf+1)*(j*PSF_NA+1));
|
covab = *(++var);
|
covab = *(++var);
|
varb = *(var += PSF_NA*npsf);
|
varb = *(var += PSF_NA*npsf);
|
f2 = sol[PSF_NA*j];
|
f2 = sol[PSF_NA*j];
|
f2 *= f2;
|
f2 *= f2;
|
obj2->poserrmx2_psf = vara/f2;
|
obj2->poserrmx2_psf = vara/f2;
|
obj2->poserrmy2_psf = varb/f2;
|
obj2->poserrmy2_psf = varb/f2;
|
obj2->poserrmxy_psf = covab/f2;
|
obj2->poserrmxy_psf = covab/f2;
|
|
|
/*------ If requested, translate variances to major and minor error axes... */
|
/*------ If requested, translate variances to major and minor error axes... */
|
if (FLAG(obj2.poserra_psf))
|
if (FLAG(obj2.poserra_psf))
|
{
|
{
|
double pmx2,pmy2,temp,theta;
|
double pmx2,pmy2,temp,theta;
|
|
|
if (fabs(temp=obj2->poserrmx2_psf-obj2->poserrmy2_psf) > 0.0)
|
if (fabs(temp=obj2->poserrmx2_psf-obj2->poserrmy2_psf) > 0.0)
|
theta = atan2(2.0 * obj2->poserrmxy_psf,temp) / 2.0;
|
theta = atan2(2.0 * obj2->poserrmxy_psf,temp) / 2.0;
|
else
|
else
|
theta = PI/4.0;
|
theta = PI/4.0;
|
|
|
temp = sqrt(0.25*temp*temp+obj2->poserrmxy_psf*obj2->poserrmxy_psf);
|
temp = sqrt(0.25*temp*temp+obj2->poserrmxy_psf*obj2->poserrmxy_psf);
|
pmy2 = pmx2 = 0.5*(obj2->poserrmx2_psf+obj2->poserrmy2_psf);
|
pmy2 = pmx2 = 0.5*(obj2->poserrmx2_psf+obj2->poserrmy2_psf);
|
pmx2+=temp;
|
pmx2+=temp;
|
pmy2-=temp;
|
pmy2-=temp;
|
|
|
obj2->poserra_psf = (float)sqrt(pmx2);
|
obj2->poserra_psf = (float)sqrt(pmx2);
|
obj2->poserrb_psf = (float)sqrt(pmy2);
|
obj2->poserrb_psf = (float)sqrt(pmy2);
|
obj2->poserrtheta_psf = theta*180.0/PI;
|
obj2->poserrtheta_psf = theta*180.0/PI;
|
}
|
}
|
|
|
/*------ ...Or ellipse parameters */
|
/*------ ...Or ellipse parameters */
|
if (FLAG(obj2.poserr_cxx))
|
if (FLAG(obj2.poserr_cxx))
|
{
|
{
|
double xm2,ym2, xym, temp;
|
double xm2,ym2, xym, temp;
|
|
|
xm2 = obj2->poserrmx2_psf;
|
xm2 = obj2->poserrmx2_psf;
|
ym2 = obj2->poserrmy2_psf;
|
ym2 = obj2->poserrmy2_psf;
|
xym = obj2->poserrmxy_psf;
|
xym = obj2->poserrmxy_psf;
|
obj2->poserrcxx_psf = (float)(ym2/(temp=xm2*ym2-xym*xym));
|
obj2->poserrcxx_psf = (float)(ym2/(temp=xm2*ym2-xym*xym));
|
obj2->poserrcyy_psf = (float)(xm2/temp);
|
obj2->poserrcyy_psf = (float)(xm2/temp);
|
obj2->poserrcxy_psf = (float)(-2*xym/temp);
|
obj2->poserrcxy_psf = (float)(-2*xym/temp);
|
}
|
}
|
return;
|
return;
|
}
|
}
|
|
|
|
|
/******************************** svdfit ************************************/
|
/******************************** svdfit ************************************/
|
/*
|
/*
|
General least-square fit A.x = b, based on Singular Value Decomposition (SVD).
|
General least-square fit A.x = b, based on Singular Value Decomposition (SVD).
|
Loosely adapted from Numerical Recipes in C, 2nd Ed. (p. 671).
|
Loosely adapted from Numerical Recipes in C, 2nd Ed. (p. 671).
|
Note: the a and v matrices are transposed with respect to the N.R. convention.
|
Note: the a and v matrices are transposed with respect to the N.R. convention.
|
*/
|
*/
|
void svdfit(double *a, float *b, int m, int n, double *sol,
|
void svdfit(double *a, float *b, int m, int n, double *sol,
|
double *vmat, double *wmat)
|
double *vmat, double *wmat)
|
{
|
{
|
#define MAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
|
#define MAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
|
(maxarg1) : (maxarg2))
|
(maxarg1) : (maxarg2))
|
#define PYTHAG(a,b) ((at=fabs(a)) > (bt=fabs(b)) ? \
|
#define PYTHAG(a,b) ((at=fabs(a)) > (bt=fabs(b)) ? \
|
(ct=bt/at,at*sqrt(1.0+ct*ct)) \
|
(ct=bt/at,at*sqrt(1.0+ct*ct)) \
|
: (bt ? (ct=at/bt,bt*sqrt(1.0+ct*ct)): 0.0))
|
: (bt ? (ct=at/bt,bt*sqrt(1.0+ct*ct)): 0.0))
|
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
|
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
|
#define TOL 1.0e-11
|
#define TOL 1.0e-11
|
|
|
int flag,i,its,j,jj,k,l,nm,mmi,nml;
|
int flag,i,its,j,jj,k,l,nm,mmi,nml;
|
double c,f,h,s,x,y,z,
|
double c,f,h,s,x,y,z,
|
anorm, g, scale,
|
anorm, g, scale,
|
at,bt,ct,maxarg1,maxarg2,
|
at,bt,ct,maxarg1,maxarg2,
|
thresh, wmax,
|
thresh, wmax,
|
*w,*ap,*ap0,*ap1,*ap10,*rv1p,*vp,*vp0,*vp1,*vp10,
|
*w,*ap,*ap0,*ap1,*ap10,*rv1p,*vp,*vp0,*vp1,*vp10,
|
*tmpp, *rv1,*tmp;
|
*tmpp, *rv1,*tmp;
|
float *bp;
|
float *bp;
|
|
|
anorm = g = scale = 0.0;
|
anorm = g = scale = 0.0;
|
if (m < n)
|
if (m < n)
|
error(EXIT_FAILURE, "*Error*: Not enough rows for solving the system ",
|
error(EXIT_FAILURE, "*Error*: Not enough rows for solving the system ",
|
"in svdfit()");
|
"in svdfit()");
|
|
|
QMALLOC(rv1, double, n);
|
QMALLOC(rv1, double, n);
|
QMALLOC(tmp, double, n);
|
QMALLOC(tmp, double, n);
|
l = nm = nml = 0; /* To avoid gcc -Wall warnings */
|
l = nm = nml = 0; /* To avoid gcc -Wall warnings */
|
for (i=0;i<n;i++)
|
for (i=0;i<n;i++)
|
{
|
{
|
l = i+1;
|
l = i+1;
|
nml = n-l;
|
nml = n-l;
|
rv1[i] = scale*g;
|
rv1[i] = scale*g;
|
g = s = scale = 0.0;
|
g = s = scale = 0.0;
|
if ((mmi = m - i) > 0)
|
if ((mmi = m - i) > 0)
|
{
|
{
|
ap = ap0 = a+i*(m+1);
|
ap = ap0 = a+i*(m+1);
|
for (k=mmi;k--;)
|
for (k=mmi;k--;)
|
scale += fabs(*(ap++));
|
scale += fabs(*(ap++));
|
if (scale)
|
if (scale)
|
{
|
{
|
for (ap=ap0,k=mmi; k--; ap++)
|
for (ap=ap0,k=mmi; k--; ap++)
|
{
|
{
|
*ap /= scale;
|
*ap /= scale;
|
s += *ap**ap;
|
s += *ap**ap;
|
}
|
}
|
f = *ap0;
|
f = *ap0;
|
g = -SIGN(sqrt(s),f);
|
g = -SIGN(sqrt(s),f);
|
h = f*g-s;
|
h = f*g-s;
|
*ap0 = f-g;
|
*ap0 = f-g;
|
ap10 = a+l*m+i;
|
ap10 = a+l*m+i;
|
for (j=nml;j--; ap10+=m)
|
for (j=nml;j--; ap10+=m)
|
{
|
{
|
for (s=0.0,ap=ap0,ap1=ap10,k=mmi; k--;)
|
for (s=0.0,ap=ap0,ap1=ap10,k=mmi; k--;)
|
s += *(ap1++)**(ap++);
|
s += *(ap1++)**(ap++);
|
f = s/h;
|
f = s/h;
|
for (ap=ap0,ap1=ap10,k=mmi; k--;)
|
for (ap=ap0,ap1=ap10,k=mmi; k--;)
|
*(ap1++) += f**(ap++);
|
*(ap1++) += f**(ap++);
|
}
|
}
|
for (ap=ap0,k=mmi; k--;)
|
for (ap=ap0,k=mmi; k--;)
|
*(ap++) *= scale;
|
*(ap++) *= scale;
|
}
|
}
|
}
|
}
|
wmat[i] = scale*g;
|
wmat[i] = scale*g;
|
g = s = scale = 0.0;
|
g = s = scale = 0.0;
|
if (i < m && i+1 != n)
|
if (i < m && i+1 != n)
|
{
|
{
|
ap = ap0 = a+i+m*l;
|
ap = ap0 = a+i+m*l;
|
for (k=nml;k--; ap+=m)
|
for (k=nml;k--; ap+=m)
|
scale += fabs(*ap);
|
scale += fabs(*ap);
|
if (scale)
|
if (scale)
|
{
|
{
|
for (ap=ap0,k=nml;k--; ap+=m)
|
for (ap=ap0,k=nml;k--; ap+=m)
|
{
|
{
|
*ap /= scale;
|
*ap /= scale;
|
s += *ap**ap;
|
s += *ap**ap;
|
}
|
}
|
f=*ap0;
|
f=*ap0;
|
g = -SIGN(sqrt(s),f);
|
g = -SIGN(sqrt(s),f);
|
h=f*g-s;
|
h=f*g-s;
|
*ap0=f-g;
|
*ap0=f-g;
|
rv1p = rv1+l;
|
rv1p = rv1+l;
|
for (ap=ap0,k=nml;k--; ap+=m)
|
for (ap=ap0,k=nml;k--; ap+=m)
|
*(rv1p++) = *ap/h;
|
*(rv1p++) = *ap/h;
|
ap10 = a+l+m*l;
|
ap10 = a+l+m*l;
|
for (j=m-l; j--; ap10++)
|
for (j=m-l; j--; ap10++)
|
{
|
{
|
for (s=0.0,ap=ap0,ap1=ap10,k=nml; k--; ap+=m,ap1+=m)
|
for (s=0.0,ap=ap0,ap1=ap10,k=nml; k--; ap+=m,ap1+=m)
|
s += *ap1**ap;
|
s += *ap1**ap;
|
rv1p = rv1+l;
|
rv1p = rv1+l;
|
for (ap1=ap10,k=nml;k--; ap1+=m)
|
for (ap1=ap10,k=nml;k--; ap1+=m)
|
*ap1 += s**(rv1p++);
|
*ap1 += s**(rv1p++);
|
}
|
}
|
for (ap=ap0,k=nml;k--; ap+=m)
|
for (ap=ap0,k=nml;k--; ap+=m)
|
*ap *= scale;
|
*ap *= scale;
|
}
|
}
|
}
|
}
|
anorm=MAX(anorm,(fabs(wmat[i])+fabs(rv1[i])));
|
anorm=MAX(anorm,(fabs(wmat[i])+fabs(rv1[i])));
|
}
|
}
|
|
|
for (i=n-1;i>=0;i--)
|
for (i=n-1;i>=0;i--)
|
{
|
{
|
if (i < n-1)
|
if (i < n-1)
|
{
|
{
|
if (g)
|
if (g)
|
{
|
{
|
ap0 = a+l*m+i;
|
ap0 = a+l*m+i;
|
vp0 = vmat+i*n+l;
|
vp0 = vmat+i*n+l;
|
vp10 = vmat+l*n+l;
|
vp10 = vmat+l*n+l;
|
g *= *ap0;
|
g *= *ap0;
|
for (ap=ap0,vp=vp0,j=nml; j--; ap+=m)
|
for (ap=ap0,vp=vp0,j=nml; j--; ap+=m)
|
*(vp++) = *ap/g;
|
*(vp++) = *ap/g;
|
for (j=nml; j--; vp10+=n)
|
for (j=nml; j--; vp10+=n)
|
{
|
{
|
for (s=0.0,ap=ap0,vp1=vp10,k=nml; k--; ap+=m)
|
for (s=0.0,ap=ap0,vp1=vp10,k=nml; k--; ap+=m)
|
s += *ap**(vp1++);
|
s += *ap**(vp1++);
|
for (vp=vp0,vp1=vp10,k=nml; k--;)
|
for (vp=vp0,vp1=vp10,k=nml; k--;)
|
*(vp1++) += s**(vp++);
|
*(vp1++) += s**(vp++);
|
}
|
}
|
}
|
}
|
vp = vmat+l*n+i;
|
vp = vmat+l*n+i;
|
vp1 = vmat+i*n+l;
|
vp1 = vmat+i*n+l;
|
for (j=nml; j--; vp+=n)
|
for (j=nml; j--; vp+=n)
|
*vp = *(vp1++) = 0.0;
|
*vp = *(vp1++) = 0.0;
|
}
|
}
|
vmat[i*n+i]=1.0;
|
vmat[i*n+i]=1.0;
|
g=rv1[i];
|
g=rv1[i];
|
l=i;
|
l=i;
|
nml = n-l;
|
nml = n-l;
|
}
|
}
|
|
|
for (i=(m<n?m:n); --i>=0;)
|
for (i=(m<n?m:n); --i>=0;)
|
{
|
{
|
l=i+1;
|
l=i+1;
|
nml = n-l;
|
nml = n-l;
|
mmi=m-i;
|
mmi=m-i;
|
g=wmat[i];
|
g=wmat[i];
|
ap0 = a+i*m+i;
|
ap0 = a+i*m+i;
|
ap10 = ap0 + m;
|
ap10 = ap0 + m;
|
for (ap=ap10,j=nml;j--;ap+=m)
|
for (ap=ap10,j=nml;j--;ap+=m)
|
*ap=0.0;
|
*ap=0.0;
|
if (g)
|
if (g)
|
{
|
{
|
g=1.0/g;
|
g=1.0/g;
|
for (j=nml;j--; ap10+=m)
|
for (j=nml;j--; ap10+=m)
|
{
|
{
|
for (s=0.0,ap=ap0,ap1=ap10,k=mmi; --k;)
|
for (s=0.0,ap=ap0,ap1=ap10,k=mmi; --k;)
|
s += *(++ap)**(++ap1);
|
s += *(++ap)**(++ap1);
|
f = (s/(*ap0))*g;
|
f = (s/(*ap0))*g;
|
for (ap=ap0,ap1=ap10,k=mmi;k--;)
|
for (ap=ap0,ap1=ap10,k=mmi;k--;)
|
*(ap1++) += f**(ap++);
|
*(ap1++) += f**(ap++);
|
}
|
}
|
for (ap=ap0,j=mmi;j--;)
|
for (ap=ap0,j=mmi;j--;)
|
*(ap++) *= g;
|
*(ap++) *= g;
|
}
|
}
|
else
|
else
|
for (ap=ap0,j=mmi;j--;)
|
for (ap=ap0,j=mmi;j--;)
|
*(ap++)=0.0;
|
*(ap++)=0.0;
|
++(*ap0);
|
++(*ap0);
|
}
|
}
|
|
|
for (k=n; --k>=0;)
|
for (k=n; --k>=0;)
|
{
|
{
|
for (its=0;its<100;its++)
|
for (its=0;its<100;its++)
|
{
|
{
|
flag=1;
|
flag=1;
|
for (l=k;l>=0;l--)
|
for (l=k;l>=0;l--)
|
{
|
{
|
nm=l-1;
|
nm=l-1;
|
if (fabs(rv1[l])+anorm == anorm)
|
if (fabs(rv1[l])+anorm == anorm)
|
{
|
{
|
flag=0;
|
flag=0;
|
break;
|
break;
|
}
|
}
|
if (fabs(wmat[nm])+anorm == anorm)
|
if (fabs(wmat[nm])+anorm == anorm)
|
break;
|
break;
|
}
|
}
|
if (flag)
|
if (flag)
|
{
|
{
|
c=0.0;
|
c=0.0;
|
s=1.0;
|
s=1.0;
|
ap0 = a+nm*m;
|
ap0 = a+nm*m;
|
ap10 = a+l*m;
|
ap10 = a+l*m;
|
for (i=l; i<=k; i++,ap10+=m)
|
for (i=l; i<=k; i++,ap10+=m)
|
{
|
{
|
f=s*rv1[i];
|
f=s*rv1[i];
|
if (fabs(f)+anorm == anorm)
|
if (fabs(f)+anorm == anorm)
|
break;
|
break;
|
g=wmat[i];
|
g=wmat[i];
|
h=PYTHAG(f,g);
|
h=PYTHAG(f,g);
|
wmat[i]=h;
|
wmat[i]=h;
|
h=1.0/h;
|
h=1.0/h;
|
c=g*h;
|
c=g*h;
|
s=(-f*h);
|
s=(-f*h);
|
for (ap=ap0,ap1=ap10,j=m; j--;)
|
for (ap=ap0,ap1=ap10,j=m; j--;)
|
{
|
{
|
z = *ap1;
|
z = *ap1;
|
y = *ap;
|
y = *ap;
|
*(ap++) = y*c+z*s;
|
*(ap++) = y*c+z*s;
|
*(ap1++) = z*c-y*s;
|
*(ap1++) = z*c-y*s;
|
}
|
}
|
}
|
}
|
}
|
}
|
z=wmat[k];
|
z=wmat[k];
|
if (l == k)
|
if (l == k)
|
{
|
{
|
if (z < 0.0)
|
if (z < 0.0)
|
{
|
{
|
wmat[k] = -z;
|
wmat[k] = -z;
|
vp = vmat+k*n;
|
vp = vmat+k*n;
|
for (j=n; j--; vp++)
|
for (j=n; j--; vp++)
|
*vp = (-*vp);
|
*vp = (-*vp);
|
}
|
}
|
break;
|
break;
|
}
|
}
|
if (its == 99)
|
if (its == 99)
|
error(EXIT_FAILURE, "*Error*: No convergence in 100 SVD iterations ",
|
error(EXIT_FAILURE, "*Error*: No convergence in 100 SVD iterations ",
|
"in svdfit()");
|
"in svdfit()");
|
x=wmat[l];
|
x=wmat[l];
|
nm=k-1;
|
nm=k-1;
|
y=wmat[nm];
|
y=wmat[nm];
|
g=rv1[nm];
|
g=rv1[nm];
|
h=rv1[k];
|
h=rv1[k];
|
f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
|
f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
|
g=PYTHAG(f,1.0);
|
g=PYTHAG(f,1.0);
|
f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
|
f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
|
c=s=1.0;
|
c=s=1.0;
|
ap10 = a+l*m;
|
ap10 = a+l*m;
|
vp10 = vmat+l*n;
|
vp10 = vmat+l*n;
|
for (j=l;j<=nm;j++,ap10+=m,vp10+=n)
|
for (j=l;j<=nm;j++,ap10+=m,vp10+=n)
|
{
|
{
|
i=j+1;
|
i=j+1;
|
g=rv1[i];
|
g=rv1[i];
|
y=wmat[i];
|
y=wmat[i];
|
h=s*g;
|
h=s*g;
|
g=c*g;
|
g=c*g;
|
z=PYTHAG(f,h);
|
z=PYTHAG(f,h);
|
rv1[j]=z;
|
rv1[j]=z;
|
c=f/z;
|
c=f/z;
|
s=h/z;
|
s=h/z;
|
f=x*c+g*s;
|
f=x*c+g*s;
|
g=g*c-x*s;
|
g=g*c-x*s;
|
h=y*s;
|
h=y*s;
|
y=y*c;
|
y=y*c;
|
for (vp=(vp1=vp10)+n,jj=n; jj--;)
|
for (vp=(vp1=vp10)+n,jj=n; jj--;)
|
{
|
{
|
z = *vp;
|
z = *vp;
|
x = *vp1;
|
x = *vp1;
|
*(vp1++) = x*c+z*s;
|
*(vp1++) = x*c+z*s;
|
*(vp++) = z*c-x*s;
|
*(vp++) = z*c-x*s;
|
}
|
}
|
z=PYTHAG(f,h);
|
z=PYTHAG(f,h);
|
wmat[j]=z;
|
wmat[j]=z;
|
if (z)
|
if (z)
|
{
|
{
|
z=1.0/z;
|
z=1.0/z;
|
c=f*z;
|
c=f*z;
|
s=h*z;
|
s=h*z;
|
}
|
}
|
f=c*g+s*y;
|
f=c*g+s*y;
|
x=c*y-s*g;
|
x=c*y-s*g;
|
for (ap=(ap1=ap10)+m,jj=m; jj--;)
|
for (ap=(ap1=ap10)+m,jj=m; jj--;)
|
{
|
{
|
z = *ap;
|
z = *ap;
|
y = *ap1;
|
y = *ap1;
|
*(ap1++) = y*c+z*s;
|
*(ap1++) = y*c+z*s;
|
*(ap++) = z*c-y*s;
|
*(ap++) = z*c-y*s;
|
}
|
}
|
}
|
}
|
rv1[l]=0.0;
|
rv1[l]=0.0;
|
rv1[k]=f;
|
rv1[k]=f;
|
wmat[k]=x;
|
wmat[k]=x;
|
}
|
}
|
}
|
}
|
|
|
wmax=0.0;
|
wmax=0.0;
|
w = wmat;
|
w = wmat;
|
for (j=n;j--; w++)
|
for (j=n;j--; w++)
|
if (*w > wmax)
|
if (*w > wmax)
|
wmax=*w;
|
wmax=*w;
|
thresh=TOL*wmax;
|
thresh=TOL*wmax;
|
w = wmat;
|
w = wmat;
|
for (j=n;j--; w++)
|
for (j=n;j--; w++)
|
if (*w < thresh)
|
if (*w < thresh)
|
*w = 0.0;
|
*w = 0.0;
|
|
|
w = wmat;
|
w = wmat;
|
ap = a;
|
ap = a;
|
tmpp = tmp;
|
tmpp = tmp;
|
for (j=n; j--; w++)
|
for (j=n; j--; w++)
|
{
|
{
|
s=0.0;
|
s=0.0;
|
if (*w)
|
if (*w)
|
{
|
{
|
bp = b;
|
bp = b;
|
for (i=m; i--;)
|
for (i=m; i--;)
|
s += *(ap++)**(bp++);
|
s += *(ap++)**(bp++);
|
s /= *w;
|
s /= *w;
|
}
|
}
|
else
|
else
|
ap += m;
|
ap += m;
|
*(tmpp++) = s;
|
*(tmpp++) = s;
|
}
|
}
|
|
|
vp0 = vmat;
|
vp0 = vmat;
|
for (j=0; j<n; j++,vp0++)
|
for (j=0; j<n; j++,vp0++)
|
{
|
{
|
s=0.0;
|
s=0.0;
|
tmpp = tmp;
|
tmpp = tmp;
|
for (vp=vp0,jj=n; jj--; vp+=n)
|
for (vp=vp0,jj=n; jj--; vp+=n)
|
s += *vp**(tmpp++);
|
s += *vp**(tmpp++);
|
sol[j]=s;
|
sol[j]=s;
|
}
|
}
|
|
|
/* Free temporary arrays */
|
/* Free temporary arrays */
|
free(tmp);
|
free(tmp);
|
free(rv1);
|
free(rv1);
|
|
|
return;
|
return;
|
}
|
}
|
|
|
#undef SIGN
|
#undef SIGN
|
#undef MAX
|
#undef MAX
|
#undef PYTHAG
|
#undef PYTHAG
|
#undef TOL
|
#undef TOL
|
|
|
/******************************** svdvar ************************************/
|
/******************************** svdvar ************************************/
|
/*
|
/*
|
Computation of the covariance matrix from the SVD vmat and wmat matrices.A
|
Computation of the covariance matrix from the SVD vmat and wmat matrices.A
|
dapted from Numerical Recipes in C, 2nd Ed. (p. 679).
|
dapted from Numerical Recipes in C, 2nd Ed. (p. 679).
|
*/
|
*/
|
void svdvar(double *v, double *w, int n, double *cov)
|
void svdvar(double *v, double *w, int n, double *cov)
|
{
|
{
|
static double wti[PSF_NTOT];
|
static double wti[PSF_NTOT];
|
double sum;
|
double sum;
|
int i,j,k;
|
int i,j,k;
|
|
|
for (i=0; i<n; i++)
|
for (i=0; i<n; i++)
|
wti[i] = w[i]? 1.0/(w[i]*w[i]) : 0.0;
|
wti[i] = w[i]? 1.0/(w[i]*w[i]) : 0.0;
|
|
|
for (i=0; i<n; i++)
|
for (i=0; i<n; i++)
|
for (j=0; j<=i; j++)
|
for (j=0; j<=i; j++)
|
{
|
{
|
for (sum=0.0,k=0; k<n; k++)
|
for (sum=0.0,k=0; k<n; k++)
|
sum += v[k*n+i]*v[k*n+j]*wti[k];
|
sum += v[k*n+i]*v[k*n+j]*wti[k];
|
cov[j*n+i] = cov[i*n+j] = sum;
|
cov[j*n+i] = cov[i*n+j] = sum;
|
}
|
}
|
|
|
return;
|
return;
|
}
|
}
|
|
|
|
|
