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

Fixed potential issues in ASSOCCOORD_TYPE WORLD mode. 

Fixed triggering of model surface brightness measurement.
Added new FLUX_DETMODEL,FLUXERR_DETMODEL, MAG_DETMODEL,MAGERR_DETMODEL, FLAGS_DE
TMODEL and CHI2_DETMODEL measurement parameters that apply on the measurement im
age a model fitted on the detection image (two PSF models required); these new p
arameters are meant for measuring faint galaxy colors.
Improved displayed info for memory allocation failures.
Increased default maximum number of threads set by configure to 1024.
Updated Copyright display line (2012).
Added support for "TPV" WCS projection.
Pushed version number to 2.16.0.
parent c3d4de93
Show whitespace changes
Inline Side-by-side
acx_urbi_resolve_dir.m4
View file @ 1755ca24
......@@ -7,7 +7,7 @@ dnl %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dnl
dnl This file part of: AstrOmatic software
dnl
dnl Copyrights: (C) 2007-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
dnl Copyrights: (C) 2007-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
dnl (C) 2007 Akim Demaille (original version)
dnl
dnl License: GPL
......@@ -24,7 +24,7 @@ dnl You should have received a copy of the GNU General Public License
dnl along with AstrOmatic software.
dnl If not, see <http://www.gnu.org/licenses/>.
dnl
dnl Last modified: 09/10/2010
dnl Last modified: 27/12/2011
dnl
dnl %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dnl
......@@ -55,7 +55,7 @@ urbi_resolve_dir ()
while true
do
eval ac_$0_res="$ac_$0_dir"
if test x"$ac_$0_dir" == x"$ac_$0_res"; then
if test x"$ac_$0_dir" = x"$ac_$0_res"; then
break
fi
ac_$0_dir=$ac_$0_res
......
configure
View file @ 1755ca24
#! /bin/sh
# Guess values for system-dependent variables and create Makefiles.
# Generated by GNU Autoconf 2.68 for sextractor 2.14.7.
# Generated by GNU Autoconf 2.68 for sextractor 2.16.0.
#
# Report bugs to <bertin@iap.fr>.
#
......@@ -714,8 +714,8 @@ MAKEFLAGS=
# Identity of this package.
PACKAGE_NAME='sextractor'
PACKAGE_TARNAME='sextractor'
PACKAGE_VERSION='2.14.7'
PACKAGE_STRING='sextractor 2.14.7'
PACKAGE_VERSION='2.16.0'
PACKAGE_STRING='sextractor 2.16.0'
PACKAGE_BUGREPORT='bertin@iap.fr'
PACKAGE_URL=''
 
......@@ -1465,7 +1465,7 @@ if test "$ac_init_help" = "long"; then
# Omit some internal or obsolete options to make the list less imposing.
# This message is too long to be a string in the A/UX 3.1 sh.
cat <<_ACEOF
\`configure' configures sextractor 2.14.7 to adapt to many kinds of systems.
\`configure' configures sextractor 2.16.0 to adapt to many kinds of systems.
 
Usage: $0 [OPTION]... [VAR=VALUE]...
 
......@@ -1535,7 +1535,7 @@ fi
 
if test -n "$ac_init_help"; then
case $ac_init_help in
short | recursive ) echo "Configuration of sextractor 2.14.7:";;
short | recursive ) echo "Configuration of sextractor 2.16.0:";;
esac
cat <<\_ACEOF
 
......@@ -1668,7 +1668,7 @@ fi
test -n "$ac_init_help" && exit $ac_status
if $ac_init_version; then
cat <<\_ACEOF
sextractor configure 2.14.7
sextractor configure 2.16.0
generated by GNU Autoconf 2.68
 
Copyright (C) 2010 Free Software Foundation, Inc.
......@@ -2296,7 +2296,7 @@ cat >config.log <<_ACEOF
This file contains any messages produced by compilers while
running configure, to aid debugging if configure makes a mistake.
 
It was created by sextractor $as_me 2.14.7, which was
It was created by sextractor $as_me 2.16.0, which was
generated by GNU Autoconf 2.68. Invocation command line was
 
$ $0 $@
......@@ -2972,7 +2972,7 @@ fi
 
# Define the identity of the package.
PACKAGE='sextractor'
VERSION='2.14.7'
VERSION='2.16.0'
 
 
cat >>confdefs.h <<_ACEOF
......@@ -20664,7 +20664,7 @@ urbi_resolve_dir ()
while true
do
eval ac_URBI_RESOLVE_DIR_PREPARE_res="$ac_URBI_RESOLVE_DIR_PREPARE_dir"
if test x"$ac_URBI_RESOLVE_DIR_PREPARE_dir" == x"$ac_URBI_RESOLVE_DIR_PREPARE_res"; then
if test x"$ac_URBI_RESOLVE_DIR_PREPARE_dir" = x"$ac_URBI_RESOLVE_DIR_PREPARE_res"; then
break
fi
ac_URBI_RESOLVE_DIR_PREPARE_dir=$ac_URBI_RESOLVE_DIR_PREPARE_res
......@@ -23325,7 +23325,7 @@ cat >>$CONFIG_STATUS <<\_ACEOF || ac_write_fail=1
# report actual input values of CONFIG_FILES etc. instead of their
# values after options handling.
ac_log="
This file was extended by sextractor $as_me 2.14.7, which was
This file was extended by sextractor $as_me 2.16.0, which was
generated by GNU Autoconf 2.68. Invocation command line was
 
CONFIG_FILES = $CONFIG_FILES
......@@ -23391,7 +23391,7 @@ _ACEOF
cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1
ac_cs_config="`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`"
ac_cs_version="\\
sextractor config.status 2.14.7
sextractor config.status 2.16.0
configured by $0, generated by GNU Autoconf 2.68,
with options \\"\$ac_cs_config\\"
 
......
configure.ac
View file @ 1755ca24
......@@ -7,7 +7,7 @@
#
# This file part of: SExtractor
#
# Copyright: (C) 2002-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
# Copyright: (C) 2002-2012 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: 06/09/2011
# Last modified: 12/04/2012
#
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
......@@ -31,7 +31,7 @@ define([AC_CACHE_LOAD],)
define([AC_CACHE_SAVE],)
# This is your standard Bertin source code...
AC_INIT(sextractor, 2.14.7, [bertin@iap.fr])
AC_INIT(sextractor, 2.16.0, [bertin@iap.fr])
AC_CONFIG_SRCDIR(src/makeit.c)
AC_CONFIG_AUX_DIR(autoconf)
AM_CONFIG_HEADER(config.h)
......
man/sex.1
View file @ 1755ca24
.TH SEXTRACTOR "1" "October 2011" "SExtractor 2.14.7" "User Commands"
.TH SEXTRACTOR "1" "April 2012" "SExtractor 2.16.0" "User Commands"
.SH NAME
sex \- extract a source catalogue from an astronomical FITS image
.SH SYNOPSIS
......
src/analyse.c
View file @ 1755ca24
......@@ -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: 19/05/2011
* Last modified: 02/11/2011
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -56,7 +56,7 @@
#include "winpos.h"
static obj2struct *obj2 = &outobj2;
extern profitstruct *theprofit;
extern profitstruct *theprofit,*thedprofit;
/********************************* analyse ***********************************/
void analyse(picstruct *field, picstruct *dfield, int objnb,
......@@ -73,6 +73,7 @@ void analyse(picstruct *field, picstruct *dfield, int objnb,
localback(field, obj);
else
obj->sigbkg = field->backsig;
obj->dsigbkg = dfield? dfield->backsig : field->backsig;
examineiso(field, dfield, obj, objlist->plist);
......@@ -416,7 +417,7 @@ void endobject(picstruct *field, picstruct *dfield, picstruct *wfield,
if (prefs.psf_flag)
thepsf->build_flag = 0; /* Reset PSF building flag */
if (prefs.dpsf_flag)
ppsf->build_flag = 0; /* Reset PSF building flag */
thedpsf->build_flag = 0; /* Reset PSF building flag */
if (FLAG(obj2.analtime))
analtime1 = counter_seconds();
......@@ -701,7 +702,7 @@ void endobject(picstruct *field, picstruct *dfield, picstruct *wfield,
if (prefs.psffit_flag)
{
if (prefs.dpsffit_flag)
double_psf_fit(ppsf, field, wfield, obj, thepsf, dfield, dwfield);
double_psf_fit(thepsf, field, wfield, obj, thedpsf, dfield, dwfield);
else
psf_fit(thepsf, field, wfield, obj);
obj2->npsf = thepsfit->npsf;
......@@ -724,6 +725,9 @@ void endobject(picstruct *field, picstruct *dfield, picstruct *wfield,
/*---- Express position error parameters in the FOCAL or WORLD frame */
if (FLAG(obj2.poserrmx2w_prof))
astrom_proferrparam(field, obj);
if (prefs.dprof_flag)
profit_dfit(theprofit, thedprofit, field, dfield, wfield, dwfield, obj,
obj2);
}
#endif
/*--- Express everything in magnitude units */
......
src/catout.c
View file @ 1755ca24
......@@ -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: 06/09/2011
* Last modified: 08/12/2011
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -414,6 +414,13 @@ void updateparamflags()
FLAG(obj2.prof_dirac_flux) |= FLAG(obj2.prof_dirac_fluxerr)
| FLAG(obj2.prof_dirac_mag);
FLAG(obj2.prof_offset_flux) |= FLAG(obj2.prof_offset_fluxerr);
FLAG(obj2.fluxerr_dprof) |= FLAG(obj2.magerr_dprof);
FLAG(obj2.flux_dprof) |= FLAG(obj2.mag_dprof)
| FLAG(obj2.fluxerr_dprof);
prefs.dprof_flag |= FLAG(obj2.dprof_chi2)
| FLAG(obj2.dprof_flag)
| FLAG(obj2.dprof_niter)
| FLAG(obj2.flux_dprof);
prefs.prof_flag |= FLAG(obj2.prof_chi2) | FLAG(obj2.prof_niter)
| FLAG(obj2.prof_vector) | FLAG(obj2.prof_errvector)
| FLAG(obj2.prof_errmatrix)
......@@ -425,7 +432,8 @@ void updateparamflags()
| FLAG(obj2.prof_disk_flux)
| FLAG(obj2.prof_spheroid_flux)
| FLAG(obj2.prof_dirac_flux)
| FLAG(obj2.prof_offset_flux);
| FLAG(obj2.prof_offset_flux)
| prefs.dprof_flag;
/* If only global parameters are requested, fit a Sersic model */
......@@ -535,7 +543,9 @@ void updateparamflags()
FLAG(obj2.mxf) |= FLAG(obj2.myf);
FLAG(obj2.mxw) |= FLAG(obj2.myw) | FLAG(obj2.mx2w) | FLAG(obj2.alphas)
| FLAG(obj2.poserr_mx2w);
| FLAG(obj2.poserr_mx2w)
| ((FLAG(obj2.assoc) | FLAG(obj2.assoc_number))
&& (prefs.assoccoord_type==ASSOCCOORD_WORLD));
FLAG(obj2.mamaposx) |= FLAG(obj2.mamaposy);
FLAG(obj2.fluxerr_win) |= FLAG(obj2.snr_win);
FLAG(obj2.flux_win) |= FLAG(obj2.mag_win)|FLAG(obj2.magerr_win)
......
src/define.h
View file @ 1755ca24
......@@ -7,7 +7,7 @@
*
* This file part of: SExtractor
*
* Copyright: (C) 1993-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 1993-2012 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/04/2011
* Last modified: 12/04/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -30,7 +30,7 @@
#ifndef HAVE_CONFIG_H
#define VERSION "2.x"
#define DATE "2009-03-31"
#define THREADS_NMAX 16 /* max. number of threads */
#define THREADS_NMAX 1024 /* max. number of threads */
#endif
/*------------------------ what, who, when and where ------------------------*/
......@@ -38,7 +38,7 @@
#define BANNER "SExtractor"
#define MYVERSION VERSION
#define EXECUTABLE "sex"
#define COPYRIGHT "2011 IAP/CNRS/UPMC"
#define COPYRIGHT "2012 IAP/CNRS/UPMC"
#define DISCLAIMER BANNER " comes with ABSOLUTELY NO WARRANTY\n" \
"You may redistribute copies of " BANNER "\n" \
"under the terms of the GNU General Public License."
......@@ -160,36 +160,62 @@
error(EXIT_FAILURE,"*Error*: file position unknown in ", \
fname)
#define QFREE(ptr) \
{free(ptr); \
ptr = NULL;}
#define QCALLOC(ptr, typ, nel) \
{if (!(ptr = (typ *)calloc((size_t)(nel),sizeof(typ)))) \
error(EXIT_FAILURE, "Not enough memory for ", \
#ptr " (" #nel " elements) !");;}
{ \
sprintf(gstr, #ptr " (" #nel "=%lld elements) " \
"at line %d in module " __FILE__ " !", \
(size_t)(nel)*sizeof(typ), __LINE__); \
error(EXIT_FAILURE, "Could not allocate memory for ", gstr);\
}; \
}
#define QMALLOC(ptr, typ, nel) \
{if (!(ptr = (typ *)malloc((size_t)(nel)*sizeof(typ)))) \
error(EXIT_FAILURE, "Not enough memory for ", \
#ptr " (" #nel " elements) !");;}
{ \
sprintf(gstr, #ptr " (" #nel "=%lld elements) " \
"at line %d in module " __FILE__ " !", \
(size_t)(nel)*sizeof(typ), __LINE__); \
error(EXIT_FAILURE, "Could not allocate memory for ", gstr);\
}; \
}
#define QMALLOC16(ptr, typ, nel) \
{if (posix_memalign((void **)&ptr, 16, (size_t)(nel)*sizeof(typ))) \
error(EXIT_FAILURE, "Not enough memory for ", \
#ptr " (" #nel " elements) !");;}
#define QFREE(ptr) \
{free(ptr); \
ptr = NULL;}
{ \
sprintf(gstr, #ptr " (" #nel "=%lld elements) " \
"at line %d in module " __FILE__ " !", \
(size_t)(nel)*sizeof(typ), __LINE__); \
error(EXIT_FAILURE, "Could not allocate memory for ", gstr);\
}; \
}
#define QREALLOC(ptr, typ, nel) \
{if (!(ptr = (typ *)realloc(ptr, (size_t)(nel)*sizeof(typ)))) \
error(EXIT_FAILURE, "Not enough memory for ", \
#ptr " (" #nel " elements) !");;}
{if (!(ptr = (typ *)realloc(ptr, (size_t)(nel)*sizeof(typ))))\
{ \
sprintf(gstr, #ptr " (" #nel "=%lld elements) " \
"at line %d in module " __FILE__ " !", \
(size_t)(nel)*sizeof(typ), __LINE__); \
error(EXIT_FAILURE, "Could not allocate memory for ", gstr);\
}; \
}
#define QMEMCPY(ptrin, ptrout, typ, nel) \
{if (ptrin) \
{if (!(ptrout = (typ *)malloc((size_t)(nel)*sizeof(typ)))) \
error(EXIT_FAILURE, "Not enough memory for ", \
#ptrout " (" #nel " elements) !"); \
memcpy(ptrout, ptrin, (size_t)(nel)*sizeof(typ));};}
{ \
sprintf(gstr, #ptrout " (" #nel "=%lld elements) " \
"at line %d in module " __FILE__ " !", \
(size_t)(nel)*sizeof(typ), __LINE__); \
error(EXIT_FAILURE,"Could not allocate memory for ",gstr);\
}; \
memcpy(ptrout, ptrin, (size_t)(nel)*sizeof(typ)); \
}; \
}
#define RINT(x) (int)(floor(x+0.5))
......
src/fitswcs.c
View file @ 1755ca24
......@@ -7,7 +7,7 @@
*
* This file part of: AstrOmatic software
*
* Copyright: (C) 1993-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 1993-2012 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: 22/07/2011
* Last modified: 12/04/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -328,7 +328,7 @@ INPUT tab structure.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 22/07/2011
VERSION 30/07/2010
***/
wcsstruct *read_wcs(tabstruct *tab)
......@@ -479,7 +479,6 @@ wcsstruct *read_wcs(tabstruct *tab)
FITSREADF(buf, "EPOCH", wcs->epoch, 2000.0);
FITSREADF(buf, "EQUINOX", wcs->equinox, wcs->epoch);
if (fitsread(buf, "RADESYS", str, H_STRING,T_STRING) != RETURN_OK)
FITSREADS(buf, "RADECSYS", str,
wcs->equinox >= 2000.0? "ICRS" : (wcs->equinox<1984.0? "FK4" : "FK5"));
if (!strcmp(str, "ICRS"))
......@@ -744,18 +743,18 @@ INPUT Proposed projection code name.
OUTPUT RETURN_OK if projection is supported, RETURN_ERROR otherwise.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 24/05/2000
VERSION 12/04/2012
***/
int wcs_supproj(char *name)
{
char projcode[26][5] =
{"AZP", "TAN", "SIN", "STG", "ARC", "ZPN", "ZEA", "AIR", "CYP", "CAR",
char projcode[28][5] =
{"TAN", "TPV", "AZP", "SIN", "STG", "ARC", "ZPN", "ZEA", "AIR", "CYP", "CAR",
"MER", "CEA", "COP", "COD", "COE", "COO", "BON", "PCO", "GLS", "PAR",
"AIT", "MOL", "CSC", "QSC", "TSC", "NONE"};
"AIT", "MOL", "CSC", "QSC", "TSC", "TNX", "NONE"};
int i;
for (i=0; i<26; i++)
for (i=0; i<28; i++)
if (!strcmp(name, projcode[i]))
return RETURN_OK;
......@@ -1050,23 +1049,23 @@ void range_wcs(wcsstruct *wcs)
/******* frame_wcs ***********************************************************
PROTO void frame_wcs(wcsstruct *wcsin, wcsstruct *wcsout)
PROTO int frame_wcs(wcsstruct *wcsin, wcsstruct *wcsout)
PURPOSE Find the x and y limits of an input frame in an output image.
INPUT WCS structure of the input frame,
WCS structure of the output frame.
OUTPUT -.
NOTES .
OUTPUT 1 if frames overlap, 0 otherwise.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 29/12/2004
VERSION 26/03/2012
***/
void frame_wcs(wcsstruct *wcsin, wcsstruct *wcsout)
int frame_wcs(wcsstruct *wcsin, wcsstruct *wcsout)
{
double rawin[NAXIS], rawout[NAXIS], world[NAXIS];
int linecount[NAXIS];
double worldc;
int *min, *max,
i,j, naxis, npoints, out, swapflag;
i,j, naxis, npoints, out, swapflag, overlapflag;
naxis = wcsin->naxis;
......@@ -1119,7 +1118,7 @@ void frame_wcs(wcsstruct *wcsin, wcsstruct *wcsout)
}
}
/* Just add a little margin, in case of... */
/* Add a little margin, just in case... */
for (i=0; i<naxis; i++)
{
if (min[i]>-2147483647)
......@@ -1128,7 +1127,16 @@ void frame_wcs(wcsstruct *wcsin, wcsstruct *wcsout)
max[i] += 2;
}
return;
/* Check overlap */
overlapflag = 1;
for (i=0; i<naxis; i++)
if (min[i]>wcsout->naxisn[i] || max[i]<0)
{
overlapflag = 0;
break;
}
return overlapflag;
}
......@@ -1570,7 +1578,7 @@ double wcs_scale(wcsstruct *wcs, double *pixpos)
/****** wcs jacobian *********************************************************
PROTO double wcs_jacobian(wcsstruct *wcs, double *pixpos, double *jacob)
PURPOSE Compute the local Jacobian matrice of the astrometric deprojection.
PURPOSE Compute the local Jacobian matrix of the astrometric deprojection.
INPUT WCS structure,
Pointer to the array of local raw coordinates,
Pointer to the jacobian array (output).
......@@ -1626,6 +1634,71 @@ double wcs_jacobian(wcsstruct *wcs, double *pixpos, double *jacob)
}
/****** wcs rawtoraw *********************************************************
PROTO double wcs_rawtoraw(wcsstruct *wcsin, wcsstruct *wcsout,
double *pixpos, double *jacob)
PURPOSE Convert raw coordinates from input wcs structure to raw coordinates
from output wcs structure, and the local Jacobian matrix of the
reprojection.
INPUT WCS input structure,
WCS output structure,
pointer to the array of local input raw coordinates,
pointer to the array of local output raw coordinates (output),
pointer to the jacobian array (output).
OUTPUT Determinant over spatial coordinates (ratio of pixel areas),
0.0 if jacob is NULL (Jacobian not computed in that case),
or -1.0 if mapping was unsuccesful.
NOTES Memory must have been allocated (naxis*naxis*sizeof(double)) for the
Jacobian array.
AUTHOR E. Bertin (IAP)
VERSION 31/03/2012
***/
double wcs_rawtoraw(wcsstruct *wcsin, wcsstruct *wcsout,
double *pixposin, double *pixposout, double *jacob)
{
double pixpos0[NAXIS], pixpos2[NAXIS], wcspos[NAXIS],
dpos;
int i,j, lng,lat,naxis;
naxis = wcsin->naxis;
for (i=0; i<naxis; i++)
pixpos0[i] = pixposin[i];
/* Coordinate transformation */
if (raw_to_wcs(wcsin, pixposin, wcspos) == RETURN_ERROR)
return -1.0;
if (wcs_to_raw(wcsout, wcspos, pixposout) == RETURN_ERROR)
return -1.0;
/* Jacobian */
if (!jacob)
return 0.0;
lng = wcsin->lng;
lat = wcsin->lat;
for (i=0; i<naxis; i++)
{
pixpos0[i] += 1.0;
if (raw_to_wcs(wcsin, pixpos0, wcspos) == RETURN_ERROR)
return -1.0;
if (wcs_to_raw(wcsout, wcspos, pixpos2) == RETURN_ERROR)
return -1.0;
pixpos0[i] -= 1.0;
for (j=0; j<naxis; j++)
jacob[j*naxis+i] = pixpos2[j] - pixposout[j];
}
if (lng==lat)
{
lng = 0;
lat = 1;
}
return fabs(jacob[lng+naxis*lng]*jacob[lat+naxis*lat]
- jacob[lat+naxis*lng]*jacob[lng+naxis*lat]);
}
/******* wcs_chirality *******************************************************
PROTO int wcs_chirality(wcsstruct *wcs)
PURPOSE Compute the chirality of a WCS projection.
......
src/fitswcs.h
View file @ 1755ca24
......@@ -7,7 +7,7 @@
*
* This file part of: AstrOmatic software
*
* Copyright: (C) 1993-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 1993-2012 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: 31/03/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -120,11 +120,15 @@ extern double fmod_0_p360(double angle),
double *wcspos1, double *wcspos2),
wcs_jacobian(wcsstruct *wcs, double *pixpos,
double *jacob),
wcs_rawtoraw(wcsstruct *wcsin, wcsstruct *wcsout,
double *pixposin, double *pixposout,
double *jacob),
wcs_scale(wcsstruct *wcs, double *pixpos);
extern int celsys_to_eq(wcsstruct *wcs, double *wcspos),
eq_to_celsys(wcsstruct *wcs, double *wcspos),
fcmp_0_p360(double anglep, double anglem),
frame_wcs(wcsstruct *wcsin, wcsstruct *wcsout),
raw_to_red(wcsstruct *wcs,
double *pixpos, double *redpos),
raw_to_wcs(wcsstruct *wcs,
......@@ -146,7 +150,6 @@ extern void b2j(double yearobs, double alphain, double deltain,
init_wcs(wcsstruct *wcs),
init_wcscelsys(wcsstruct *wcs),
invert_wcs(wcsstruct *wcs),
frame_wcs(wcsstruct *wcsin, wcsstruct *wcsout),
j2b(double yearobs, double alphain, double deltain,
double *alphaout, double *deltaout),
precess(double yearin, double alphain, double deltain,
......
src/makeit.c
View file @ 1755ca24
......@@ -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: 22/04/2011
* Last modified: 02/11/2011
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -57,7 +57,7 @@
static int selectext(char *filename);
time_t thetimet, thetimet2;
extern profitstruct *theprofit,*thepprofit,*theqprofit;
extern profitstruct *theprofit,*thedprofit, *thepprofit,*theqprofit;
extern char profname[][32];
double dtime;
......@@ -75,6 +75,7 @@ void makeit()
patternstruct *pattern;
static time_t thetime1, thetime2;
struct tm *tm;
unsigned int modeltype;
int nflag[MAXFLAG], nparam2[2],
i, nok, ntab, next, ntabmax, forcextflag,
nima0,nima1, nweight0,nweight1, npat;
......@@ -110,9 +111,13 @@ void makeit()
if (prefs.psf_flag)
{
NFPRINTF(OUTPUT, "Reading PSF information");
thepsf = psf_load(prefs.psf_name[0]);
if (prefs.dpsf_flag)
ppsf = psf_load(prefs.psf_name[1]);
{
thedpsf = psf_load(prefs.psf_name[0]);
thepsf = psf_load(prefs.psf_name[1]);
}
else
thepsf = psf_load(prefs.psf_name[0]);
/*-- Need to check things up because of PSF context parameters */
updateparamflags();
useprefs();
......@@ -124,19 +129,21 @@ void makeit()
fft_init(prefs.nthreads);
/* Create profiles at full resolution */
NFPRINTF(OUTPUT, "Preparing profile models");
theprofit = profit_init(thepsf,
(FLAG(obj2.prof_offset_flux)? MODEL_BACK : MODEL_NONE)
modeltype = (FLAG(obj2.prof_offset_flux)? MODEL_BACK : MODEL_NONE)
|(FLAG(obj2.prof_dirac_flux)? MODEL_DIRAC : MODEL_NONE)
|(FLAG(obj2.prof_spheroid_flux)?
(FLAG(obj2.prof_spheroid_sersicn)?
MODEL_SERSIC : MODEL_DEVAUCOULEURS) : MODEL_NONE)
|(FLAG(obj2.prof_disk_flux)? MODEL_EXPONENTIAL : MODEL_NONE)
|(FLAG(obj2.prof_bar_flux)? MODEL_BAR : MODEL_NONE)
|(FLAG(obj2.prof_arms_flux)? MODEL_ARMS : MODEL_NONE));
|(FLAG(obj2.prof_arms_flux)? MODEL_ARMS : MODEL_NONE);
theprofit = profit_init(thepsf, modeltype);
changecatparamarrays("VECTOR_MODEL", &theprofit->nparam, 1);
changecatparamarrays("VECTOR_MODELERR", &theprofit->nparam, 1);
nparam2[0] = nparam2[1] = theprofit->nparam;
changecatparamarrays("MATRIX_MODELERR", nparam2, 2);
if (prefs.dprof_flag)
thedprofit = profit_init(thedpsf, modeltype);
if (prefs.pattern_flag)
{
npat = prefs.prof_disk_patternvectorsize;
......@@ -445,11 +452,11 @@ void makeit()
if (prefs.psf_flag)
{
psf_readcontext(thepsf, field);
psf_init(thepsf);
psf_init();
if (prefs.dpsf_flag)
{
psf_readcontext(thepsf, dfield);
psf_init(thepsf); /*?*/
psf_init();
}
}
......@@ -546,6 +553,8 @@ void makeit()
if (prefs.prof_flag)
{
profit_end(theprofit);
if (prefs.dprof_flag)
profit_end(thedprofit);
if (FLAG(obj2.prof_concentration)|FLAG(obj2.prof_concentration))
{
profit_end(thepprofit);
......@@ -556,10 +565,10 @@ void makeit()
#endif
if (prefs.psf_flag)
psf_end(thepsf,thepsfit); /*?*/
psf_end(thepsf, thepsfit);
if (prefs.dpsf_flag)
psf_end(ppsf,ppsfit);
psf_end(thedpsf, thedpsfit);
if (FLAG(obj2.sprob))
neurclose();
......
src/paramprofit.h
View file @ 1755ca24
......@@ -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: 09/07/2011
* Last modified: 08/11/2011
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -679,3 +679,27 @@
&outobj2.prof_arms_quadfracerr, H_FLOAT, T_FLOAT, "%6.4f", "deg",
"stat.error;phot.count;arith.ratio;src.morph;stat.fit.param", "deg"},
*/
{"CHI2_DETMODEL", "Reduced Chi2 of the det. model fit to measurement image",
&outobj2.dprof_chi2, H_FLOAT, T_FLOAT, "%12.7g", "",
"stat.fit.chi2;src.morph", ""},
{"FLAGS_DETMODEL", "Detection model-fitting flags",
&outobj2.dprof_flag, H_INT, T_BYTE, "%3d", "",
"meta.code;stat.fit;src.morph", ""},
{"NITER_DETMODEL", "Number of iterations for detection model-fitting",
&outobj2.dprof_niter, H_INT, T_SHORT, "%3d", "",
"meta.number;stat.fit;src.morph", ""},
{"FLUX_DETMODEL", "Flux from detection model-fitting",
&outobj2.flux_dprof, H_FLOAT, T_FLOAT, "%12.7g", "count",
"phot.count;stat.fit.param", "ct"},
{"FLUXERR_DETMODEL", "RMS error on detection model-fitting flux",
&outobj2.fluxerr_dprof, H_FLOAT, T_FLOAT, "%12.7g", "count",
"stat.error;phot.count;stat.fit.param", "ct"},
{"MAG_DETMODEL", "Magnitude from detection model-fitting",
&outobj2.mag_dprof, H_FLOAT, T_FLOAT, "%8.4f", "mag",
"phot.mag;stat.fit.param", "mag"},
{"MAGERR_DETMODEL", "RMS error on detection model-fitting magnitude",
&outobj2.magerr_dprof, H_FLOAT, T_FLOAT, "%8.4f", "mag",
"stat.error;phot.mag;stat.fit.param", "mag"},
src/photom.c
View file @ 1755ca24
......@@ -7,7 +7,7 @@
*
* This file part of: SExtractor
*
* Copyright: (C) 1993-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 1993-2011 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: 19/10/2010
* Last modified: 08/11/2011
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -1025,6 +1025,15 @@ void computemags(picstruct *field, objstruct *obj)
/obj2->prof_arms_flux
:99.0;
if (FLAG(obj2.mag_dprof))
obj2->mag_dprof = obj2->flux_dprof>0.0?
-2.5*log10(obj2->flux_dprof) + prefs.mag_zeropoint
:99.0;
if (FLAG(obj2.magerr_dprof))
obj2->magerr_dprof = obj2->flux_dprof>0.0?
1.086*obj2->fluxerr_dprof/obj2->flux_dprof
:99.0;
/* Mag. WINdowed */
if (FLAG(obj2.mag_win))
obj2->mag_win = obj2->flux_win>0.0?
......
src/prefs.c
View file @ 1755ca24
......@@ -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: 19/05/2011
* Last modified: 02/11/2011
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -581,7 +581,10 @@ void useprefs()
if (FLAG(obj2.flux_psf) )
{
prefs.psffit_flag = 1;
prefs.dpsffit_flag = (prefs.npsf_name>1); /*?*/
/* We deactivate double-PSF fits for now */
/*
prefs.dpsffit_flag = (prefs.npsf_name>1);
*/
}
if (prefs.check_flag)
for (i=0; i<prefs.ncheck_type; i++)
......@@ -618,6 +621,8 @@ void useprefs()
prefs.prof_flag = 1;
if (prefs.prof_flag)
prefs.psf_flag = 1;
if (prefs.dprof_flag)
prefs.dpsf_flag = 1;
/*-------------------------- Tracking the PSF FWHM --------------------------*/
if (prefs.seeing_fwhm == 0 && FLAG(obj2.sprob) || FLAG(obj2.fwhm_psf))
......
src/prefs.h
View file @ 1755ca24
......@@ -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: 09/03/2011
* Last modified: 02/11/2011
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -228,6 +228,7 @@ typedef struct
int pc_flag; /* PC-fit needed */
int pc_vectorsize; /* nb of params */
int prof_flag; /* Profile-fitting */
int dprof_flag; /* Det. Prof.-fitting */
int pattern_flag; /* Pattern-fitting */
/*----- Profile-fitting */
int prof_vectorsize; /* nb of params */
......
src/profit.c
View file @ 1755ca24
......@@ -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: 08/10/2011
* Last modified: 08/12/2011
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -77,7 +77,7 @@ const int flux_flag[PARAM_NPARAM] = {0,
/* "Local" global variables for debugging purposes */
int theniter, the_gal;
static picstruct *the_field, *the_wfield;
profitstruct *theprofit, *thepprofit, *theqprofit;
profitstruct *theprofit,*thedprofit, *thepprofit, *theqprofit;
/****** profit_init ***********************************************************
PROTO profitstruct profit_init(psfstruct *psf, unsigned int modeltype)
......@@ -166,7 +166,7 @@ OUTPUT Pointer to an allocated fit structure (containing details about the
fit).
NOTES It is a modified version of the lm_minimize() of lmfit.
AUTHOR E. Bertin (IAP)
VERSION 06/09/2011
VERSION 03/11/2011
***/
void profit_fit(profitstruct *profit,
picstruct *field, picstruct *wfield,
......@@ -269,15 +269,19 @@ void profit_fit(profitstruct *profit,
profit_psf(profit);
/* Set initial guesses and boundaries */
profit->sigma = obj->sigbkg;
if ((profit->guessradius = 0.5*psf->fwhm) < obj2->hl_radius)
profit->guessradius = obj2->hl_radius;
profit->guesssigbkg = profit->sigma = obj->sigbkg;
profit->guessdx = obj->mx - (int)(obj->mx+0.49999);
profit->guessdy = obj->my - (int)(obj->my+0.49999);
if ((profit->guessflux = obj2->flux_auto) <= 0.0)
profit->guessflux = 0.0;
if ((profit->guessfluxmax = 10.0*obj2->fluxerr_auto) <= profit->guessflux)
profit->guessfluxmax = profit->guessflux;
if (profit->guessfluxmax <= 0.0)
profit->guessfluxmax = 1.0;
if ((profit->guessradius = 0.5*psf->fwhm) < obj2->hl_radius)
profit->guessradius = obj2->hl_radius;
profit->guessaspect = obj->b/obj->a;
profit->guessposang = obj->theta;
profit_resetparams(profit);
......@@ -341,6 +345,7 @@ if (list[i] && i!= PARAM_SPHEROID_ASPECT && i!=PARAM_SPHEROID_POSANG)
profit->freeparam_flag[index[i]] = 0;
profit->niter = profit_minimize(profit, PROFIT_MAXITER);
*/
if (profit->nlimmin)
obj2->prof_flag |= PROFLAG_MINLIM;
if (profit->nlimmax)
......@@ -755,9 +760,14 @@ profit->resi);
}
psf = pprofit->psf;
pprofit->pixstep = profit->pixstep;
pprofit->guessradius = profit->guessradius;
pprofit->guesssigbkg = profit->guesssigbkg;
pprofit->guessdx = profit->guessdx;
pprofit->guessdy = profit->guessdy;
pprofit->guessflux = profit->guessflux;
pprofit->guessfluxmax = profit->guessfluxmax;
pprofit->guessradius = profit->guessradius;
pprofit->guessaspect = profit->guessaspect;
pprofit->guessposang = profit->guessposang;
pprofit->ix = profit->ix;
pprofit->iy = profit->iy;
pprofit->objnaxisn[0] = profit->objnaxisn[0];
......@@ -790,9 +800,14 @@ profit->resi);
QFREE(qprofit->psfdft);
}
qprofit->pixstep = profit->pixstep;
qprofit->guessradius = profit->guessradius;
qprofit->guesssigbkg = profit->guesssigbkg;
qprofit->guessdx = profit->guessdx;
qprofit->guessdy = profit->guessdy;
qprofit->guessflux = profit->guessflux;
qprofit->guessfluxmax = profit->guessfluxmax;
qprofit->guessradius = profit->guessradius;
qprofit->guessaspect = profit->guessaspect;
qprofit->guessposang = profit->guessposang;
qprofit->ix = profit->ix;
qprofit->iy = profit->iy;
qprofit->objnaxisn[0] = profit->objnaxisn[0];
......@@ -859,6 +874,182 @@ profit->resi);
}
/****** profit_dfit ***********************************************************
PROTO void profit_dfit(profitstruct *profit, profitstruct *dprofit,
picstruct *field, picstruct *dfield,
picstruct *wfield, picstruct *dwfield,
objstruct *obj, obj2struct *obj2)
PURPOSE Fit profile(s) convolved with the PSF to a detected object on the
detection image, and use the measurement image to scale the flux.
INPUT Pointer to the measurement profile-fitting structure,
pointer to the detection profile-fitting structure,
pointer to the measurement field,
pointer to the detection field,
pointer to the measurement field weight,
pointer to the detection field weight,
pointer to the obj.
OUTPUT Pointer to an allocated fit structure (containing details about the
fit).
NOTES It is a modified version of the lm_minimize() of lmfit.
AUTHOR E. Bertin (IAP)
VERSION 08/12/2011
***/
void profit_dfit(profitstruct *profit, profitstruct *dprofit,
picstruct *field, picstruct *dfield,
picstruct *wfield, picstruct *dwfield,
objstruct *obj, obj2struct *obj2)
{
psfstruct *dpsf;
double emx2,emy2,emxy, a , cp,sp, cn, bn, n,
sumn,sumd;
PIXTYPE valn,vald,w2;
float param0[PARAM_NPARAM], param1[PARAM_NPARAM],
param[PARAM_NPARAM],
**list,
*cov, *pix,
psf_fwhm, dchi2, err, aspect, chi2, ffac;
int *index,
i,j,p, nparam, nparam2, ncomp, nprof;
nparam = dprofit->nparam;
nparam2 = nparam*nparam;
nprof = dprofit->nprof;
if (dprofit->psfdft)
{
QFREE(dprofit->psfdft);
}
dpsf = dprofit->psf;
dprofit->pixstep = dpsf->pixstep;
obj2->dprof_flag = 0;
/* Create pixmaps at image resolution */
dprofit->ix = (int)(obj->mx + 0.49999);/* internal convention: 1st pix = 0 */
dprofit->iy = (int)(obj->my + 0.49999);/* internal convention: 1st pix = 0 */
psf_fwhm = dpsf->masksize[0]*dpsf->pixstep;
dprofit->objnaxisn[0] = (((int)((obj->xmax-obj->xmin+1) + psf_fwhm + 0.499)
*1.2)/2)*2 + 1;
dprofit->objnaxisn[1] = (((int)((obj->ymax-obj->ymin+1) + psf_fwhm + 0.499)
*1.2)/2)*2 + 1;
if (dprofit->objnaxisn[1]<dprofit->objnaxisn[0])
dprofit->objnaxisn[1] = dprofit->objnaxisn[0];
else
dprofit->objnaxisn[0] = dprofit->objnaxisn[1];
if (dprofit->objnaxisn[0]>PROFIT_MAXOBJSIZE)
{
dprofit->subsamp = ceil((float)dprofit->objnaxisn[0]/PROFIT_MAXOBJSIZE);
dprofit->objnaxisn[1] = (dprofit->objnaxisn[0] /= (int)dprofit->subsamp);
obj2->dprof_flag |= PROFLAG_OBJSUB;
}
else
dprofit->subsamp = 1.0;
dprofit->nobjpix = dprofit->objnaxisn[0]*dprofit->objnaxisn[1];
/* Use (dirty) global variables to interface with lmfit */
the_field = dfield;
the_wfield = dwfield;
dprofit->obj = obj;
dprofit->obj2 = obj2;
dprofit->nresi = profit_copyobjpix(dprofit, dfield, dwfield);
/* Check if the number of constraints exceeds the number of free parameters */
if (dprofit->nresi < nparam)
{
obj2->dprof_flag |= PROFLAG_NOTCONST;
return;
}
/* Create pixmap at PSF resolution */
dprofit->modnaxisn[0] =
((int)(dprofit->objnaxisn[0]*dprofit->subsamp/dprofit->pixstep
+0.4999)/2+1)*2;
dprofit->modnaxisn[1] =
((int)(dprofit->objnaxisn[1]*dprofit->subsamp/dprofit->pixstep
+0.4999)/2+1)*2;
if (dprofit->modnaxisn[1] < dprofit->modnaxisn[0])
dprofit->modnaxisn[1] = dprofit->modnaxisn[0];
else
dprofit->modnaxisn[0] = dprofit->modnaxisn[1];
if (dprofit->modnaxisn[0]>PROFIT_MAXMODSIZE)
{
dprofit->pixstep = (double)dprofit->modnaxisn[0] / PROFIT_MAXMODSIZE;
dprofit->modnaxisn[0] = dprofit->modnaxisn[1] = PROFIT_MAXMODSIZE;
obj2->dprof_flag |= PROFLAG_MODSUB;
}
dprofit->nmodpix = dprofit->modnaxisn[0]*dprofit->modnaxisn[1];
/* Compute the local PSF */
profit_psf(dprofit);
/* Set initial guesses and boundaries */
dprofit->guesssigbkg = dprofit->sigma = obj->dsigbkg;
dprofit->guessdx = obj->mx - (int)(obj->mx+0.49999);
dprofit->guessdy = obj->my - (int)(obj->my+0.49999);
if ((dprofit->guessflux = obj->dflux) <= 0.0)
dprofit->guessflux = 0.0;
if ((dprofit->guessfluxmax = 10.0*obj->dnpix*obj->dsigbkg*obj->dsigbkg)
<= dprofit->guessflux)
dprofit->guessfluxmax = dprofit->guessflux;
if (dprofit->guessfluxmax <= 0.0)
dprofit->guessfluxmax = 1.0;
if ((dprofit->guessradius = 0.5*dpsf->fwhm) < sqrtf((float)obj->dnpix))
dprofit->guessradius = sqrtf((float)obj->dnpix);
dprofit->guessaspect = obj->b/obj->a;
dprofit->guessposang = obj->theta;
profit_resetparams(dprofit);
/* Actual minimisation */
fft_reset();
dprofit->niter = profit_minimize(dprofit, PROFIT_MAXITER);
if (dprofit->nlimmin)
obj2->dprof_flag |= PROFLAG_MINLIM;
if (dprofit->nlimmax)
obj2->dprof_flag |= PROFLAG_MAXLIM;
for (p=0; p<nparam; p++)
dprofit->paramerr[p]= sqrt(dprofit->covar[p*(nparam+1)]);
obj2->dprof_niter = dprofit->niter;
/* Now inject fitted parameters into the measurement model */
fft_reset();
profit_residuals(profit,field,wfield, 0.0, dprofit->paraminit, NULL);
/* Compute flux correction */
sumn = sumd = 0.0;
for (p=0; p<profit->nobjpix; p++)
if (profit->objweight[p]>0 && profit->objpix[p]>-BIG)
{
w2 = profit->objweight[p]*profit->objweight[p] * profit->lmodpix[p];
sumn += (double)(w2*profit->objpix[p]);
sumd += (double)(w2*profit->lmodpix[p]);
}
ffac = (float)(sumn/sumd);
obj2->flux_dprof = sumd!=0.0? dprofit->flux*ffac: 0.0f;
obj2->fluxerr_dprof = sumd!=0.0? 1.0f/sqrtf((float)sumd): 0.0f;
if (FLAG(obj2.dprof_chi2))
{
/*-- Compute reduced chi2 on measurement image */
pix = profit->lmodpix;
for (p=profit->nobjpix; p--;)
*(pix++) *= ffac;
profit_compresi(profit, 0.0, profit->resi);
obj2->dprof_chi2 = (profit->nresi > dprofit->nparam)?
profit->chi2 / (profit->nresi - dprofit->nparam) : 0.0;
}
/* clean up. */
fft_reset();
return;
}
/****** profit_noisearea ******************************************************
PROTO float profit_noisearea(profitstruct *profit)
PURPOSE Return the equivalent noise area (see King 1983) of a model.
......@@ -1558,7 +1749,7 @@ float *profit_residuals(profitstruct *profit, picstruct *field,
/****** profit_compresi ******************************************************
PROTO float *prof_compresi(profitstruct *profit, float dynparam,
PROTO float *profit_compresi(profitstruct *profit, float dynparam,
float *resi)
PURPOSE Compute the vector of residuals between the data and the galaxy
profile model.
......@@ -2800,16 +2991,14 @@ INPUT Pointer to the profit structure,
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 18/05/2011
VERSION 03/11/2011
***/
void profit_resetparam(profitstruct *profit, paramenum paramtype)
{
objstruct *obj;
obj2struct *obj2;
float param, parammin,parammax, range;
parfitenum fittype;
obj = profit->obj;
obj2 = profit->obj2;
param = parammin = parammax = 0.0; /* Avoid gcc -Wall warnings*/
......@@ -2818,12 +3007,12 @@ void profit_resetparam(profitstruct *profit, paramenum paramtype)
case PARAM_BACK:
fittype = PARFIT_LINBOUND;
param = 0.0;
parammin = -6.0*obj->sigbkg;
parammax = 6.0*obj->sigbkg;
parammin = -6.0*profit->guesssigbkg;
parammax = 6.0*profit->guesssigbkg;
break;
case PARAM_X:
fittype = PARFIT_LINBOUND;
param = obj->mx - (int)(obj->mx+0.49999);
param = profit->guessdx;
range = profit->guessradius*4.0;
if (range>profit->objnaxisn[0]*2.0)
range = profit->objnaxisn[0]*2.0;
......@@ -2832,7 +3021,7 @@ void profit_resetparam(profitstruct *profit, paramenum paramtype)
break;
case PARAM_Y:
fittype = PARFIT_LINBOUND;
param = obj->my - (int)(obj->my+0.49999);
param = profit->guessdy;
range = profit->guessradius*4.0;
if (range>profit->objnaxisn[1]*2)
range = profit->objnaxisn[1]*2;
......@@ -2854,19 +3043,19 @@ void profit_resetparam(profitstruct *profit, paramenum paramtype)
case PARAM_SPHEROID_REFF:
fittype = PARFIT_LOGBOUND;
param = FLAG(obj2.prof_disk_flux)? profit->guessradius
: profit->guessradius*sqrtf(obj->a/obj->b);
: profit->guessradius/sqrtf(profit->guessaspect);
parammin = 0.01;
parammax = param * 10.0;
break;
case PARAM_SPHEROID_ASPECT:
fittype = PARFIT_LOGBOUND;
param = FLAG(obj2.prof_disk_flux)? 1.0 : obj->b/obj->a;
param = FLAG(obj2.prof_disk_flux)? 1.0 : profit->guessaspect;
parammin = FLAG(obj2.prof_disk_flux)? 0.5 : 0.01;
parammax = FLAG(obj2.prof_disk_flux)? 2.0 : 100.0;
break;
case PARAM_SPHEROID_POSANG:
fittype = PARFIT_UNBOUND;
param = obj->theta;
param = profit->guessposang;
parammin = 90.0;
parammax = 90.0;
break;
......@@ -2884,19 +3073,19 @@ void profit_resetparam(profitstruct *profit, paramenum paramtype)
break;
case PARAM_DISK_SCALE: /* From scalelength to Re */
fittype = PARFIT_LOGBOUND;
param = profit->guessradius/1.67835*sqrtf(obj->a/obj->b);
param = profit->guessradius/(1.67835*sqrtf(profit->guessaspect));
parammin = 0.01/1.67835;
parammax = param * 10.0;
break;
case PARAM_DISK_ASPECT:
fittype = PARFIT_LOGBOUND;
param = obj->b/obj->a;
param = profit->guessaspect;
parammin = 0.01;
parammax = 100.0;
break;
case PARAM_DISK_POSANG:
fittype = PARFIT_UNBOUND;
param = obj->theta;
param = profit->guessposang;
parammin = 90.0;
parammax = 90.0;
break;
......@@ -3313,7 +3502,7 @@ INPUT Pointer to the profile-fitting structure,
OUTPUT A pointer to an allocated prof structure.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 08/10/2011
VERSION 08/12/2011
***/
profstruct *prof_init(profitstruct *profit, unsigned int modeltype)
{
......@@ -3341,6 +3530,7 @@ profstruct *prof_init(profitstruct *profit, unsigned int modeltype)
prof->naxisn[1] = PROFIT_MAXMODSIZE;
prof->naxisn[2] = 1;
prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2];
QMALLOC(prof->pix, float, prof->npix);
prof->typscale = 1.0;
profit_addparam(profit, PARAM_X, &prof->x[0]);
profit_addparam(profit, PARAM_Y, &prof->x[1]);
......@@ -3370,7 +3560,7 @@ profstruct *prof_init(profitstruct *profit, unsigned int modeltype)
prof->naxisn[1] = PROFIT_MAXMODSIZE;
prof->naxisn[2] = 1;
prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2];
QMALLOC(prof->pix, float, profit->nmodpix);
QMALLOC(prof->pix, float, prof->npix);
prof->typscale = 1.0;
profit_addparam(profit, PARAM_X, &prof->x[0]);
profit_addparam(profit, PARAM_Y, &prof->x[1]);
......@@ -3386,7 +3576,7 @@ profstruct *prof_init(profitstruct *profit, unsigned int modeltype)
prof->naxisn[1] = PROFIT_MAXMODSIZE;
prof->naxisn[2] = 1;
prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2];
QMALLOC(prof->pix, float, profit->nmodpix);
QMALLOC(prof->pix, float, prof->npix);
prof->typscale = 1.0;
profit_addparam(profit, PARAM_X, &prof->x[0]);
profit_addparam(profit, PARAM_Y, &prof->x[1]);
......@@ -3402,7 +3592,7 @@ profstruct *prof_init(profitstruct *profit, unsigned int modeltype)
prof->naxisn[1] = PROFIT_MAXMODSIZE;
prof->naxisn[2] = 1;
prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2];
QMALLOC(prof->pix, float, profit->nmodpix);
QMALLOC(prof->pix, float, prof->npix);
prof->typscale = 1.0;
profit_addparam(profit, PARAM_X, &prof->x[0]);
profit_addparam(profit, PARAM_Y, &prof->x[1]);
......@@ -3425,7 +3615,7 @@ profstruct *prof_init(profitstruct *profit, unsigned int modeltype)
prof->naxisn[1] = PROFIT_MAXMODSIZE;
prof->naxisn[2] = 1;
prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2];
QMALLOC(prof->pix, float, profit->nmodpix);
QMALLOC(prof->pix, float, prof->npix);
prof->typscale = 1.0;
profit_addparam(profit, PARAM_X, &prof->x[0]);
profit_addparam(profit, PARAM_Y, &prof->x[1]);
......@@ -3444,7 +3634,7 @@ profstruct *prof_init(profitstruct *profit, unsigned int modeltype)
prof->naxisn[1] = PROFIT_MAXMODSIZE;
prof->naxisn[2] = 1;
prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2];
QMALLOC(prof->pix, float, profit->nmodpix);
QMALLOC(prof->pix, float, prof->npix);
prof->typscale = 1.0;
profit_addparam(profit, PARAM_X, &prof->x[0]);
profit_addparam(profit, PARAM_Y, &prof->x[1]);
......@@ -3463,7 +3653,7 @@ profstruct *prof_init(profitstruct *profit, unsigned int modeltype)
prof->naxisn[1] = PROFIT_MAXMODSIZE;
prof->naxisn[2] = 1;
prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2];
QMALLOC(prof->pix, float, profit->nmodpix);
QMALLOC(prof->pix, float, prof->npix);
prof->typscale = 1.0;
profit_addparam(profit, PARAM_X, &prof->x[0]);
profit_addparam(profit, PARAM_Y, &prof->x[1]);
......@@ -3477,10 +3667,11 @@ profstruct *prof_init(profitstruct *profit, unsigned int modeltype)
case MODEL_TABULATED: /* An example of tabulated profile */
prof->name = "tabulated model";
prof->naxis = 3;
width = prof->naxisn[0] = PROFIT_PROFRES;
height = prof->naxisn[1] = PROFIT_PROFRES;
nsub = prof->naxisn[2] = PROFIT_PROFSRES;
QCALLOC(prof->pix, float, width*height*nsub);
width = prof->naxisn[0] = PROFIT_MAXMODSIZE;
height = prof->naxisn[1] = PROFIT_MAXMODSIZE;
nsub = prof->naxisn[2] = PROFIT_MAXSMODSIZE;
prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2];
QMALLOC(prof->pix, float, prof->npix);
ixc = width/2;
iyc = height/2;
rmax2 = (ixc - 1.0)*(ixc - 1.0);
......@@ -3488,7 +3679,7 @@ profstruct *prof_init(profitstruct *profit, unsigned int modeltype)
prof->typscale = re2;
re2 *= re2;
zero = prof->extrazero[0] = 0.2;
scale = prof->extrascale[0]= 8.0/PROFIT_PROFSRES;
scale = prof->extrascale[0]= 8.0/PROFIT_MAXSMODSIZE;
pix = prof->pix;
for (s=0; s<nsub; s++)
{
......@@ -3645,6 +3836,7 @@ float prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag)
switch(prof->code)
{
case MODEL_DIRAC:
memset(prof->pix, 0, profit->nmodpix);
prof->pix[profit->modnaxisn[0]/2
+ (profit->modnaxisn[1]/2)*profit->modnaxisn[0]] = 1.0;
prof->lostfluxfrac = 0.0;
......
src/profit.h
View file @ 1755ca24
......@@ -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: 20/05/2011
* Last modified: 08/12/2011
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -74,11 +74,10 @@
#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_MAXSMODSIZE 64 /* Number of model planes */
#define PROFIT_MAXOBJSIZE 512 /* 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 PROFIT_PROFRES 256 /* Pixmap size of model components */
#define PROFIT_PROFSRES 64 /* Number of model subcomponents */
#define INTERP_MAXKERNELWIDTH 8 /* Max. range of kernel (pixels) */
/* NOTES:
One must have: PROFIT_NITER > 0
......@@ -190,9 +189,13 @@ typedef struct
float sigma; /* Standard deviation of the pixel values */
float flux; /* Total flux in final convolved model */
float spirindex; /* Spiral index (>0 for CCW) */
float guessradius; /* Best guess for source half-light radius */
float guesssigbkg; /* Best guess for background noise sigma */
float guessdx,guessdy;/* Best guess for relative source coordinates */
float guessflux; /* Best guess for typical source flux (>0) */
float guessfluxmax; /* Best guess for source flux upper limit (>0)*/
float guessradius; /* Best guess for source half-light radius */
float guessaspect; /* Best guess for source aspect ratio */
float guessposang; /* Best guess for source position angle */
/* Buffers */
double dparam[PARAM_NPARAM];
} profitstruct;
......@@ -232,7 +235,11 @@ int profit_boundtounbound(profitstruct *profit,
profit_unboundtobound(profitstruct *profit,
double *dparam, float *param, int index);
void prof_end(profstruct *prof),
void profit_dfit(profitstruct *profit, profitstruct *dprofit,
picstruct *field, picstruct *dfield,
picstruct *wfield, picstruct *dwfield,
objstruct *obj, obj2struct *obj2),
prof_end(profstruct *prof),
profit_addparam(profitstruct *profit, paramenum paramindex,
float **param),
profit_fit(profitstruct *profit,
......
src/psf.c
View file @ 1755ca24
......@@ -7,7 +7,7 @@
*
* This file part of: SExtractor
*
* Copyright: (C) 1998-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 1998-2011 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/10/2011
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -55,19 +55,21 @@ extern objstruct outobj;
/*
Allocate memory and stuff for the PSF-fitting.
*/
void psf_init(psfstruct *psf)
void psf_init(void)
{
QMALLOC(thepsfit, psfitstruct, 1);
QMALLOC(thepsfit->x, double, prefs.psf_npsfmax);
QMALLOC(thepsfit->y, double, prefs.psf_npsfmax);
QMALLOC(thepsfit->flux, float, prefs.psf_npsfmax);
QMALLOC(thepsfit->fluxerr, float, prefs.psf_npsfmax);
QMALLOC(ppsfit, psfitstruct, 1); /*?*/
QMALLOC(ppsfit->x, double, prefs.psf_npsfmax);
QMALLOC(ppsfit->y, double, prefs.psf_npsfmax);
QMALLOC(ppsfit->flux, float, prefs.psf_npsfmax);
QMALLOC(ppsfit->fluxerr, float, prefs.psf_npsfmax);
if (prefs.dpsf_flag)
{
QMALLOC(thedpsfit, psfitstruct, 1);
QMALLOC(thedpsfit->x, double, prefs.psf_npsfmax);
QMALLOC(thedpsfit->y, double, prefs.psf_npsfmax);
QMALLOC(thedpsfit->flux, float, prefs.psf_npsfmax);
QMALLOC(thedpsfit->fluxerr, float, prefs.psf_npsfmax);
}
return;
}
......@@ -735,9 +737,9 @@ er */
/*******************************************************************************
****/
void double_psf_fit(psfstruct *ppsf, picstruct *pfield, picstruct *pwfield,
objstruct *obj, psfstruct *psf, picstruct *field,
picstruct *wfield)
void double_psf_fit(psfstruct *psf, picstruct *field, picstruct *wfield,
objstruct *obj, psfstruct *dpsf, picstruct *dfield,
picstruct *dwfield)
{
static double /* sum[PSF_NPSFMAX]*/ pdeltax[PSF_NPSFMAX],
pdeltay[PSF_NPSFMAX],psol[PSF_NPSFMAX], pcovmat[PSF_NPSFMAX*PSF_NPSFMAX],
......@@ -750,8 +752,8 @@ void double_psf_fit(psfstruct *ppsf, picstruct *pfield, picstruct *pwfield,
val, ppix,ppix2, /* dflux, */
gain, radmin2,radmax2,satlevel
,chi2,pwthresh,pbacknoise2, /* mr, */
r2=0, psf_fwhm,ppsf_fwhm ;
float **ppsfmasks, **ppsfmaskx,**ppsfmasky, *pps;
r2=0, dpsf_fwhm,psf_fwhm ;
float **psfmasks, **psfmaskx,**psfmasky, *pps;
float *pdata, *pdata2, *pdata3, *pweight, *pd, *pw,
*pdh, *pwh, pixstep,ppixstep;
PIXTYPE *pdatah, *pweighth;
......@@ -765,53 +767,53 @@ void double_psf_fit(psfstruct *ppsf, picstruct *pfield, picstruct *pwfield,
static obj2struct *obj2 = &outobj2;
pdx = pdy =dx = dy = 0.0;
ppixstep = 1.0/ppsf->pixstep;
pixstep = 1.0/psf->pixstep;
gain = (field->gain >0.0? field->gain: 1e30);
ppixstep = 1.0/psf->pixstep;
pixstep = 1.0/dpsf->pixstep;
gain = (dfield->gain >0.0? dfield->gain: 1e30);
npsfmax=prefs.psf_npsfmax;
pbacknoise2 = pfield->backsig*pfield->backsig;
satlevel = field->satur_level - obj->bkg;
pbacknoise2 = field->backsig*field->backsig;
satlevel = dfield->satur_level - obj->bkg;
gainflag = prefs.weightgain_flag;
dpsf_fwhm = dpsf->fwhm*dpsf->pixstep;
psf_fwhm = psf->fwhm*psf->pixstep;
ppsf_fwhm = ppsf->fwhm*ppsf->pixstep;
pwthresh = pwfield?pwfield->weight_thresh:BIG;
pwthresh = wfield?wfield->weight_thresh:BIG;
/* Initialize outputs */
ppsfit->niter = 0;
ppsfit->npsf = 0;
thepsfit->niter = 0;
thepsfit->npsf = 0;
for (j=0; j<npsfmax; j++)
{
ppsfit->x[j] = 999999.0;
ppsfit->y[j] = 999999.0;
ppsfit->flux[j] = 0.0;
ppsfit->fluxerr[j] = 0.0;
thepsfit->x[j] = 999999.0;
thepsfit->y[j] = 999999.0;
thepsfit->flux[j] = 0.0;
thepsfit->fluxerr[j] = 0.0;
pdeltax[j]= pdeltay[j]=psol[j]= pwmat[j]=pflux[j]=pfluxerr[j]=0.0;
}
ix = (obj->xmax+obj->xmin+1)/2;
iy = (obj->ymax+obj->ymin+1)/2;
width = obj->xmax-obj->xmin+1+psf_fwhm;
if (width < (ival=(int)(psf_fwhm*2)))
width = obj->xmax-obj->xmin+1+dpsf_fwhm;
if (width < (ival=(int)(dpsf_fwhm*2)))
width = ival;
height = obj->ymax-obj->ymin+1+psf_fwhm;
if (height < (ival=(int)(psf_fwhm*2)))
height = obj->ymax-obj->ymin+1+dpsf_fwhm;
if (height < (ival=(int)(dpsf_fwhm*2)))
height = ival;
npix = width*height;
radmin2 = PSF_MINSHIFT*PSF_MINSHIFT;
radmax2 = npix/2.0;
psf_fit(psf,field, wfield,obj);
npsf=thepsfit->npsf;
psf_fit(dpsf,dfield, dwfield,obj);
npsf=thedpsfit->npsf;
QMALLOC(ppsfmasks,float *,npsfmax);
QMALLOC(ppsfmaskx,float *,npsfmax);
QMALLOC(ppsfmasky,float *,npsfmax);
QMALLOC(psfmasks,float *,npsfmax);
QMALLOC(psfmaskx,float *,npsfmax);
QMALLOC(psfmasky,float *,npsfmax);
for (i=0; i<npsfmax; i++)
{
QMALLOC(ppsfmasks[i],float,npix);
QMALLOC(ppsfmaskx[i],float,npix);
QMALLOC(ppsfmasky[i],float,npix);
QMALLOC(psfmasks[i],float,npix);
QMALLOC(psfmaskx[i],float,npix);
QMALLOC(psfmasky[i],float,npix);
}
QMALLOC(pweighth, PIXTYPE, npix);
......@@ -824,19 +826,19 @@ void double_psf_fit(psfstruct *ppsf, picstruct *pfield, picstruct *pwfield,
for (j=0; j<npsf; j++)
{
pdeltax[j] =thepsfit->x[j]-ix-1 ;
pdeltay[j] =thepsfit->y[j]-iy-1 ;
ppsfit->flux[j] = 0;
ppsfit->fluxerr[j] = 0;
pdeltax[j] =thedpsfit->x[j]-ix-1 ;
pdeltay[j] =thedpsfit->y[j]-iy-1 ;
thepsfit->flux[j] = 0;
thepsfit->fluxerr[j] = 0;
}
/*------------------- Now the photometry fit ---------------------*/
copyimage(pfield, pdatah, width, height, ix, iy);
copyimage(field, pdatah, width, height, ix, iy);
/* Compute photometry weights */
wbad = 0;
if (pwfield)
if (wfield)
{
copyimage(pwfield, pweighth, width, height, ix, iy);
copyimage(wfield, pweighth, width, height, ix, iy);
for (pwh=pweighth, pw=pweight, pdh=pdatah,p=npix; p--;)
{
if ((ppix=*(pwh++)) < pwthresh && ppix>0
......@@ -877,7 +879,7 @@ void double_psf_fit(psfstruct *ppsf, picstruct *pfield, picstruct *pwfield,
/* Get the photmetry PSF */
psf_build(ppsf);
psf_build(psf);
for (j=1; j<=npsf; j++)
{
if (j>1)
......@@ -890,11 +892,11 @@ void double_psf_fit(psfstruct *ppsf, picstruct *pfield, picstruct *pwfield,
for (k=0; k<j-1; k++)
{
pd = pdata2;
pps = ppsfmasks[k];
pps = psfmasks[k];
for (p=npix; p--;)
*(pd++) -= pflux[k]**(pps++);
}
convolve_image(pfield, pdata2, pdata3, width,height);
convolve_image(field, pdata2, pdata3, width,height);
/*---- Ignore regions too close to stellar cores */
for (k=0; k<j-1; k++)
{
......@@ -914,12 +916,12 @@ void double_psf_fit(psfstruct *ppsf, picstruct *pfield, picstruct *pwfield,
for (k=0; k<j; k++)
{
/*------ Resample the PSFs here for the 1st iteration */
vignet_resample(ppsf->maskloc,
ppsf->masksize[0], ppsf->masksize[1],
ppsfmasks[k], width, height,
vignet_resample(psf->maskloc,
psf->masksize[0], psf->masksize[1],
psfmasks[k], width, height,
-pdeltax[k]*ppixstep, -pdeltay[k]*ppixstep,
ppixstep);
pm=compute_gradient_phot(pweight,width,height, ppsfmasks[k],pm);
pm=compute_gradient_phot(pweight,width,height, psfmasks[k],pm);
}
svdfit(pmat, pdata, npix, j, psol, pvmat, pwmat);
......@@ -942,7 +944,7 @@ void double_psf_fit(psfstruct *ppsf, picstruct *pfield, picstruct *pwfield,
for (pd=pdata,pw=pweight,p=0; p<npix; pw++,p++)
{
ppix = *(pd++);
ppix -= ppsfmasks[j][p]*pflux[j]**pw;
ppix -= psfmasks[j][p]*pflux[j]**pw;
chi2 += ppix*ppix;
if (chi2>1E29) chi2=1E28;
}
......@@ -962,7 +964,7 @@ void double_psf_fit(psfstruct *ppsf, picstruct *pfield, picstruct *pwfield,
for (pd=pdata,pw=pweight,p=0; p<npix; pw++,p++)
{
ppix = *(pd++)/pflux[j];
ppix -= ppsfmasks[j][p]**pw;
ppix -= psfmasks[j][p]**pw;
chi2 += ppix*ppix;
if (chi2>1E29) chi2=1E28;
}
......@@ -972,33 +974,33 @@ void double_psf_fit(psfstruct *ppsf, picstruct *pfield, picstruct *pwfield,
}
}
ppsfit->niter = thepsfit->niter;
ppsfit->npsf = npsf;
thepsfit->niter = thedpsfit->niter;
thepsfit->npsf = npsf;
for (j=0; j<npsf; j++)
{
thedpsfit->x[j] = ix+pdeltax[j]+1.0;
thedpsfit->y[j] = iy+pdeltay[j]+1.0;
thedpsfit->flux[j] = pflux[j];
thedpsfit->fluxerr[j] = pfluxerr[j];
thepsfit->x[j] = ix+pdeltax[j]+1.0;
thepsfit->y[j] = iy+pdeltay[j]+1.0;
thepsfit->flux[j] = pflux[j];
thepsfit->fluxerr[j] = pfluxerr[j];
ppsfit->x[j] = ix+pdeltax[j]+1.0;
ppsfit->y[j] = iy+pdeltay[j]+1.0;
ppsfit->flux[j] = pflux[j];
ppsfit->fluxerr[j] = pfluxerr[j];
}
for (i=0; i<npsfmax; i++)
{
QFREE(ppsfmasks[i]);
QFREE(ppsfmaskx[i]);
QFREE(ppsfmasky[i]);
QFREE(psfmasks[i]);
QFREE(psfmaskx[i]);
QFREE(psfmasky[i]);
}
QFREE(ppsfmasks);
QFREE(ppsfmaskx);
QFREE(ppsfmasky);
QFREE(psfmasks);
QFREE(psfmaskx);
QFREE(psfmasky);
QFREE(pdatah);
QFREE(pdata);
QFREE(pdata2);
......
src/psf.h
View file @ 1755ca24
......@@ -7,7 +7,7 @@
*
* This file part of: SExtractor
*
* Copyright: (C) 1998-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 1998-2011 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/10/2011
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -102,8 +102,8 @@ typedef struct
} psfitstruct;
/*----------------------------- Global variables ----------------------------*/
psfstruct *psf,*ppsf,*thepsf;
psfitstruct *thepsfit,*ppsfit,*psfit;
psfstruct *psf,*thedpsf,*thepsf;
psfitstruct *thepsfit,*thedpsfit;
PIXTYPE *checkmask;
/*-------------------------------- functions --------------------------------*/
......@@ -116,7 +116,7 @@ extern void compute_pos(int *pnpsf,int *pconvflag,int *pnpsfflag,
double *x2, double *y2,double *xy, int npsf),
psf_build(psfstruct *psf),
psf_end(psfstruct *psf, psfitstruct *psfit),
psf_init(psfstruct *psf),
psf_init(void),
svdfit(double *a, float *b, int m, int n, double *sol,
double *vmat, double *wmat),
svdvar(double *vmat, double *wmat, int n, double *covmat);
......
src/scan.c
View file @ 1755ca24
......@@ -7,7 +7,7 @@
*
* This file part of: SExtractor
*
* Copyright: (C) 1993-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 1993-2011 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: 21/12/2011
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -59,7 +59,7 @@ INPUT Measurement field pointer,
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 28/09/2010
VERSION 21/12/2011
***/
void scanimage(picstruct *field, picstruct *dfield, picstruct **pffield,
int nffield, picstruct *wfield, picstruct *dwfield)
......@@ -79,7 +79,7 @@ void scanimage(picstruct *field, picstruct *dfield, picstruct **pffield,
nposize, stacksize, w, h, blankh, maxpixnb,
varthreshflag, ontotal;
short trunflag;
PIXTYPE thresh, relthresh, cdnewsymbol, cdvar,cdwthresh,wthresh,
PIXTYPE thresh, relthresh, cdnewsymbol, cdwthresh,wthresh,
*scan,*dscan,*cdscan,*dwscan,*dwscanp,*dwscann,
*cdwscan,*cdwscanp,*cdwscann,*wscand,
*scant, *wscan,*wscann,*wscanp;
......@@ -350,7 +350,7 @@ void scanimage(picstruct *field, picstruct *dfield, picstruct **pffield,
curpixinfo.flag = trunflag;
if (varthreshflag)
thresh = relthresh*sqrt(cdvar = ((xl==w || yl==h)? 0.0:cdwscan[xl]));
thresh = relthresh*sqrt((xl==w || yl==h)? 0.0:cdwscan[xl]);
luflag = cdnewsymbol > thresh?1:0;
if (luflag)
......
src/types.h
View file @ 1755ca24
......@@ -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: 19/05/2011
* Last modified: 08/12/2011
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -124,7 +124,7 @@ typedef struct
float cxx,cyy,cxy; /* ellipse parameters */
int firstpix; /* ptr to first pixel */
int lastpix; /* ptr to last pixel */
float bkg, dbkg, sigbkg; /* Background stats (ADU) */
float bkg, dbkg, sigbkg, dsigbkg; /* Background stats (ADU) */
float thresh; /* measur. threshold (ADU) */
float dthresh; /* detect. threshold (ADU) */
float mthresh; /* max. threshold (ADU) */
......@@ -512,6 +512,13 @@ typedef struct
float prof_arms_starterrw; /* RMS error */
float prof_arms_quadfrac; /* Arms quadrature fraction */
float prof_arms_quadfracerr; /* RMS error */
float dprof_chi2; /* Det. model fit reduced chi2*/
BYTE dprof_flag; /* Detection model flags*/
short dprof_niter; /* # of detection model iter. */
float flux_dprof; /* Flux in detection model*/
float fluxerr_dprof; /* Error on detect model flux */
float mag_dprof; /* Mag from detection model */
float magerr_dprof; /* RMS mag from detect. model */
/* ---- MEF ----*/
short ext_number; /* FITS extension number */
/* ---- Time ---- */
......
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