Commit 54bca520 authored by Emmanuel Bertin's avatar Emmanuel Bertin
Browse files

Fixed FLUX_HYBRID and MAG_HYBRID centering issues.

Added 2nd order moments measurements of the convolved model.
Fixed issues introduced in a recent updates with check-images.
Implemented new algorithms for SPREAD_MODEL and SPREADERR_MODEL (no optimization yet and the scaling is still variable however).
Added experimental SNR_WIN measurement parameter.
Fixed crash with empty ASSOC lists.
Pushed version number to 2.13.6.
parent 1b0fb5bd
#! /bin/sh #! /bin/sh
# Guess values for system-dependent variables and create Makefiles. # Guess values for system-dependent variables and create Makefiles.
# Generated by GNU Autoconf 2.63 for sextractor 2.13.2. # Generated by GNU Autoconf 2.63 for sextractor 2.13.6.
# #
# Report bugs to <bertin@iap.fr>. # Report bugs to <bertin@iap.fr>.
# #
...@@ -750,8 +750,8 @@ SHELL=${CONFIG_SHELL-/bin/sh} ...@@ -750,8 +750,8 @@ SHELL=${CONFIG_SHELL-/bin/sh}
# Identity of this package. # Identity of this package.
PACKAGE_NAME='sextractor' PACKAGE_NAME='sextractor'
PACKAGE_TARNAME='sextractor' PACKAGE_TARNAME='sextractor'
PACKAGE_VERSION='2.13.2' PACKAGE_VERSION='2.13.6'
PACKAGE_STRING='sextractor 2.13.2' PACKAGE_STRING='sextractor 2.13.6'
PACKAGE_BUGREPORT='bertin@iap.fr' PACKAGE_BUGREPORT='bertin@iap.fr'
   
ac_unique_file="src/makeit.c" ac_unique_file="src/makeit.c"
...@@ -1508,7 +1508,7 @@ if test "$ac_init_help" = "long"; then ...@@ -1508,7 +1508,7 @@ if test "$ac_init_help" = "long"; then
# Omit some internal or obsolete options to make the list less imposing. # 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. # This message is too long to be a string in the A/UX 3.1 sh.
cat <<_ACEOF cat <<_ACEOF
\`configure' configures sextractor 2.13.2 to adapt to many kinds of systems. \`configure' configures sextractor 2.13.6 to adapt to many kinds of systems.
   
Usage: $0 [OPTION]... [VAR=VALUE]... Usage: $0 [OPTION]... [VAR=VALUE]...
   
...@@ -1578,7 +1578,7 @@ fi ...@@ -1578,7 +1578,7 @@ fi
   
if test -n "$ac_init_help"; then if test -n "$ac_init_help"; then
case $ac_init_help in case $ac_init_help in
short | recursive ) echo "Configuration of sextractor 2.13.2:";; short | recursive ) echo "Configuration of sextractor 2.13.6:";;
esac esac
cat <<\_ACEOF cat <<\_ACEOF
   
...@@ -1711,7 +1711,7 @@ fi ...@@ -1711,7 +1711,7 @@ fi
test -n "$ac_init_help" && exit $ac_status test -n "$ac_init_help" && exit $ac_status
if $ac_init_version; then if $ac_init_version; then
cat <<\_ACEOF cat <<\_ACEOF
sextractor configure 2.13.2 sextractor configure 2.13.6
generated by GNU Autoconf 2.63 generated by GNU Autoconf 2.63
   
Copyright (C) 1992, 1993, 1994, 1995, 1996, 1998, 1999, 2000, 2001, Copyright (C) 1992, 1993, 1994, 1995, 1996, 1998, 1999, 2000, 2001,
...@@ -1725,7 +1725,7 @@ cat >config.log <<_ACEOF ...@@ -1725,7 +1725,7 @@ cat >config.log <<_ACEOF
This file contains any messages produced by compilers while This file contains any messages produced by compilers while
running configure, to aid debugging if configure makes a mistake. running configure, to aid debugging if configure makes a mistake.
   
It was created by sextractor $as_me 2.13.2, which was It was created by sextractor $as_me 2.13.6, which was
generated by GNU Autoconf 2.63. Invocation command line was generated by GNU Autoconf 2.63. Invocation command line was
   
$ $0 $@ $ $0 $@
...@@ -2428,7 +2428,7 @@ fi ...@@ -2428,7 +2428,7 @@ fi
   
# Define the identity of the package. # Define the identity of the package.
PACKAGE='sextractor' PACKAGE='sextractor'
VERSION='2.13.2' VERSION='2.13.6'
   
   
cat >>confdefs.h <<_ACEOF cat >>confdefs.h <<_ACEOF
...@@ -28593,7 +28593,7 @@ exec 6>&1 ...@@ -28593,7 +28593,7 @@ exec 6>&1
# report actual input values of CONFIG_FILES etc. instead of their # report actual input values of CONFIG_FILES etc. instead of their
# values after options handling. # values after options handling.
ac_log=" ac_log="
This file was extended by sextractor $as_me 2.13.2, which was This file was extended by sextractor $as_me 2.13.6, which was
generated by GNU Autoconf 2.63. Invocation command line was generated by GNU Autoconf 2.63. Invocation command line was
   
CONFIG_FILES = $CONFIG_FILES CONFIG_FILES = $CONFIG_FILES
...@@ -28656,7 +28656,7 @@ Report bugs to <bug-autoconf@gnu.org>." ...@@ -28656,7 +28656,7 @@ Report bugs to <bug-autoconf@gnu.org>."
_ACEOF _ACEOF
cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1 cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1
ac_cs_version="\\ ac_cs_version="\\
sextractor config.status 2.13.2 sextractor config.status 2.13.6
configured by $0, generated by GNU Autoconf 2.63, configured by $0, generated by GNU Autoconf 2.63,
with options \\"`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`\\" with options \\"`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`\\"
   
......
...@@ -22,7 +22,7 @@ ...@@ -22,7 +22,7 @@
# 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: 24/01/2011 # Last modified: 03/05/2011
# #
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
...@@ -31,7 +31,7 @@ define([AC_CACHE_LOAD],) ...@@ -31,7 +31,7 @@ define([AC_CACHE_LOAD],)
define([AC_CACHE_SAVE],) define([AC_CACHE_SAVE],)
# This is your standard Bertin source code... # This is your standard Bertin source code...
AC_INIT(sextractor, 2.13.2, [bertin@iap.fr]) AC_INIT(sextractor, 2.13.6, [bertin@iap.fr])
AC_CONFIG_SRCDIR(src/makeit.c) AC_CONFIG_SRCDIR(src/makeit.c)
AC_CONFIG_AUX_DIR(autoconf) AC_CONFIG_AUX_DIR(autoconf)
AM_CONFIG_HEADER(config.h) AM_CONFIG_HEADER(config.h)
......
.TH SEXTRACTOR "1" "March 2011" "SExtractor 2.13.2" "User Commands" .TH SEXTRACTOR "1" "May 2011" "SExtractor 2.13.6" "User Commands"
.SH NAME .SH NAME
sex \- extract a source catalogue from an astronomical FITS image sex \- extract a source catalogue from an astronomical FITS image
.SH SYNOPSIS .SH SYNOPSIS
......
...@@ -22,7 +22,7 @@ ...@@ -22,7 +22,7 @@
* 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: 12/01/2011 * Last modified: 04/05/2011
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
...@@ -109,7 +109,7 @@ assocstruct *load_assoc(char *filename, wcsstruct *wcs) ...@@ -109,7 +109,7 @@ assocstruct *load_assoc(char *filename, wcsstruct *wcs)
if (!(file = fopen(filename, "r"))) if (!(file = fopen(filename, "r")))
return NULL; return NULL;
assoc = NULL; /* To avoid gcc -Wall warnings */ QCALLOC(assoc, assocstruct, 1);
list = NULL; /* To avoid gcc -Wall warnings */ list = NULL; /* To avoid gcc -Wall warnings */
data = NULL; /* To avoid gcc -Wall warnings */ data = NULL; /* To avoid gcc -Wall warnings */
ispoon = ncol = ndata = nlist = size = spoonsize = xindex = yindex ispoon = ncol = ndata = nlist = size = spoonsize = xindex = yindex
...@@ -170,8 +170,7 @@ assocstruct *load_assoc(char *filename, wcsstruct *wcs) ...@@ -170,8 +170,7 @@ assocstruct *load_assoc(char *filename, wcsstruct *wcs)
nlist = ndata+3; nlist = ndata+3;
/*---- Allocate memory for the ASSOC struct and the filtered list */ /*---- Allocate memory for the filtered list */
QMALLOC(assoc, assocstruct, 1);
ispoon = ASSOC_BUFINC/(nlist*sizeof(double)); ispoon = ASSOC_BUFINC/(nlist*sizeof(double));
spoonsize = ispoon*nlist; spoonsize = ispoon*nlist;
QMALLOC(assoc->list, double, size = spoonsize); QMALLOC(assoc->list, double, size = spoonsize);
...@@ -210,11 +209,15 @@ assocstruct *load_assoc(char *filename, wcsstruct *wcs) ...@@ -210,11 +209,15 @@ assocstruct *load_assoc(char *filename, wcsstruct *wcs)
fclose(file); fclose(file);
free(data); free(data);
QREALLOC(assoc->list, double, i*nlist);
assoc->nobj = i; assoc->nobj = i;
assoc->radius = prefs.assoc_radius; if (i>0)
assoc->ndata = ndata; {
assoc->ncol = nlist; QREALLOC(assoc->list, double, i*nlist);
assoc->radius = prefs.assoc_radius;
assoc->ndata = ndata;
assoc->ncol = nlist;
}
return assoc; return assoc;
} }
...@@ -236,6 +239,9 @@ void init_assoc(picstruct *field) ...@@ -236,6 +239,9 @@ void init_assoc(picstruct *field)
error(EXIT_FAILURE, "*Error*: Assoc-list file not found: ", error(EXIT_FAILURE, "*Error*: Assoc-list file not found: ",
prefs.assoc_name); prefs.assoc_name);
if (assoc->nobj==0)
warning(prefs.assoc_name, " ASSOC input-list is empty");
/* Sort the assoc-list by y coordinates, and build the hash table */ /* Sort the assoc-list by y coordinates, and build the hash table */
sort_assoc(field, assoc); sort_assoc(field, assoc);
......
...@@ -7,7 +7,7 @@ ...@@ -7,7 +7,7 @@
* *
* This file part of: SExtractor * 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 * License: GNU General Public License
* *
...@@ -22,7 +22,7 @@ ...@@ -22,7 +22,7 @@
* 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: 11/10/2010 * Last modified: 23/03/2011
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
...@@ -60,7 +60,7 @@ void makeback(picstruct *field, picstruct *wfield, int wscale_flag) ...@@ -60,7 +60,7 @@ void makeback(picstruct *field, picstruct *wfield, int wscale_flag)
w,bw, bh, nx,ny,nb, w,bw, bh, nx,ny,nb,
lflag, nr; lflag, nr;
float *ratio,*ratiop, *weight, *sigma, float *ratio,*ratiop, *weight, *sigma,
sratio; sratio, sigfac;
/* If the weight-map is not an external one, no stats are needed for it */ /* If the weight-map is not an external one, no stats are needed for it */
if (wfield && wfield->flags&(INTERP_FIELD|BACKRMS_FIELD)) if (wfield && wfield->flags&(INTERP_FIELD|BACKRMS_FIELD))
...@@ -294,7 +294,7 @@ void makeback(picstruct *field, picstruct *wfield, int wscale_flag) ...@@ -294,7 +294,7 @@ void makeback(picstruct *field, picstruct *wfield, int wscale_flag)
filterback(wfield); filterback(wfield);
/* Compute normalization for variance- or weight-maps*/ /* Compute normalization for variance- or weight-maps*/
if (wfield && wscale_flag && wfield->flags&(VAR_FIELD|WEIGHT_FIELD)) if (wfield && wfield->flags&(VAR_FIELD|WEIGHT_FIELD))
{ {
nr = 0; nr = 0;
QMALLOC(ratio, float, wfield->nback); QMALLOC(ratio, float, wfield->nback);
...@@ -308,18 +308,27 @@ void makeback(picstruct *field, picstruct *wfield, int wscale_flag) ...@@ -308,18 +308,27 @@ void makeback(picstruct *field, picstruct *wfield, int wscale_flag)
*(ratiop++) = sratio; *(ratiop++) = sratio;
nr++; nr++;
} }
wfield->sigfac = fqmedian(ratio, nr); sigfac = fqmedian(ratio, nr);
for (i=0; i<nr && ratio[i]<=0.0; i++); for (i=0; i<nr && ratio[i]<=0.0; i++);
if (i<nr) if (i<nr)
wfield->sigfac = fqmedian(ratio+i, nr-i); sigfac = fqmedian(ratio+i, nr-i);
else else
{ {
warning("Null or negative global weighting factor:","defaulted to 1"); warning("Null or negative global weighting factor:","defaulted to 1");
wfield->sigfac = 1.0; sigfac = 1.0;
} }
free(ratio); free(ratio);
if (wscale_flag)
wfield->sigfac = sigfac;
else
{
wfield->sigfac = 1.0;
field->backsig /= sigfac;
}
} }
/* Compute 2nd derivatives along the y-direction */ /* Compute 2nd derivatives along the y-direction */
NFPRINTF(OUTPUT, "Computing background d-map"); NFPRINTF(OUTPUT, "Computing background d-map");
free(field->dback); free(field->dback);
......
...@@ -7,7 +7,7 @@ ...@@ -7,7 +7,7 @@
* *
* This file part of: SExtractor * 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 * License: GNU General Public License
* *
...@@ -22,7 +22,7 @@ ...@@ -22,7 +22,7 @@
* 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: 19/10/2010 * Last modified: 12/04/2011
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
...@@ -189,6 +189,17 @@ void updateparamflags() ...@@ -189,6 +189,17 @@ void updateparamflags()
FLAG(obj2.prof_disk_scale) |= prefs.pattern_flag; FLAG(obj2.prof_disk_scale) |= prefs.pattern_flag;
FLAG(obj2.prof_concentration) |= FLAG(obj2.prof_concentrationerr); FLAG(obj2.prof_concentration) |= FLAG(obj2.prof_concentrationerr);
FLAG(obj2.fluxcorerr_prof) |= FLAG(obj2.magcorerr_prof);
FLAG(obj2.fluxcor_prof) |= FLAG(obj2.magcor_prof)
| FLAG(obj2.fluxcorerr_prof);
FLAG(obj2.fluxerr_prof) |= FLAG(obj2.magerr_prof)
| FLAG(obj2.prof_concentrationerr)
| FLAG(obj2.fluxcorerr_prof);
FLAG(obj2.flux_prof) |= FLAG(obj2.mag_prof)
| FLAG(obj2.fluxerr_prof)
| FLAG(obj2.fluxcor_prof);
FLAG(obj2.poserraw_prof) |= FLAG(obj2.poserrbw_prof); FLAG(obj2.poserraw_prof) |= FLAG(obj2.poserrbw_prof);
FLAG(obj2.poserrcxxw_prof) |= FLAG(obj2.poserrcyyw_prof) FLAG(obj2.poserrcxxw_prof) |= FLAG(obj2.poserrcyyw_prof)
| FLAG(obj2.poserrcxyw_prof); | FLAG(obj2.poserrcxyw_prof);
...@@ -226,17 +237,9 @@ void updateparamflags() ...@@ -226,17 +237,9 @@ void updateparamflags()
| FLAG(obj2.poserrcxx_prof) | FLAG(obj2.poserrcxx_prof)
| FLAG(obj2.poserrmx2_prof) | FLAG(obj2.poserrmx2_prof)
| FLAG(obj2.prof_concentration) | FLAG(obj2.prof_concentration)
| FLAG(obj2.prof_class_star); | FLAG(obj2.prof_class_star)
FLAG(obj2.fluxcorerr_prof) |= FLAG(obj2.magcorerr_prof);
FLAG(obj2.fluxcor_prof) |= FLAG(obj2.magcor_prof)
| FLAG(obj2.fluxcorerr_prof);
FLAG(obj2.fluxerr_prof) |= FLAG(obj2.magerr_prof)
| FLAG(obj2.prof_concentrationerr)
| FLAG(obj2.fluxcorerr_prof);
FLAG(obj2.flux_prof) |= FLAG(obj2.mag_prof)
| FLAG(obj2.fluxerr_prof)
| FLAG(obj2.fluxcor_prof); | FLAG(obj2.fluxcor_prof);
FLAG(obj2.prof_dirac_mag) |= FLAG(obj2.prof_dirac_magerr); FLAG(obj2.prof_dirac_mag) |= FLAG(obj2.prof_dirac_magerr);
FLAG(obj2.prof_spheroid_mag) |= FLAG(obj2.prof_spheroid_magerr); FLAG(obj2.prof_spheroid_mag) |= FLAG(obj2.prof_spheroid_magerr);
FLAG(obj2.prof_spheroid_reff) |= FLAG(obj2.prof_spheroid_refferr); FLAG(obj2.prof_spheroid_reff) |= FLAG(obj2.prof_spheroid_refferr);
...@@ -365,6 +368,10 @@ void updateparamflags() ...@@ -365,6 +368,10 @@ void updateparamflags()
| FLAG(obj2.prof_e1) | FLAG(obj2.prof_pol1) | FLAG(obj2.prof_e1) | FLAG(obj2.prof_pol1)
| FLAG(obj2.prof_a) | FLAG(obj2.prof_cxx) | FLAG(obj2.prof_a) | FLAG(obj2.prof_cxx)
| FLAG(obj2.prof_mx2w); | FLAG(obj2.prof_mx2w);
FLAG(obj2.prof_conva) |= FLAG(obj2.prof_convb) | FLAG(obj2.prof_convtheta);
FLAG(obj2.prof_convcxx) |= FLAG(obj2.prof_convcyy)
/*| FLAG(obj2.fluxcor_prof) */;
FLAG(obj2.prof_convmx2) |= FLAG(obj2.prof_convcxx) | FLAG(obj2.prof_conva);
FLAG(obj2.fluxmean_prof) |= FLAG(obj2.mumean_prof); FLAG(obj2.fluxmean_prof) |= FLAG(obj2.mumean_prof);
FLAG(obj2.fluxeff_prof) |= FLAG(obj2.mueff_prof) FLAG(obj2.fluxeff_prof) |= FLAG(obj2.mueff_prof)
| FLAG(obj2.fluxmean_prof); | FLAG(obj2.fluxmean_prof);
...@@ -471,8 +478,7 @@ void updateparamflags() ...@@ -471,8 +478,7 @@ void updateparamflags()
FLAG(obj2.winposerr_mx2) |= FLAG(obj2.winposerr_my2) FLAG(obj2.winposerr_mx2) |= FLAG(obj2.winposerr_my2)
| FLAG(obj2.winposerr_mxy) | FLAG(obj2.winposerr_mxy)
| FLAG(obj2.winposerr_a) | FLAG(obj2.winposerr_cxx) | FLAG(obj2.winposerr_a) | FLAG(obj2.winposerr_cxx)
| FLAG(obj2.winposerr_mx2w) | FLAG(obj2.winposerr_mx2w);
| FLAG(obj2.fluxerr_win) | FLAG(obj2.magerr_win);
FLAG(obj2.winpos_alpha1950) |= FLAG(obj2.winpos_delta1950) FLAG(obj2.winpos_alpha1950) |= FLAG(obj2.winpos_delta1950)
| FLAG(obj2.win_theta1950) | FLAG(obj2.win_theta1950)
...@@ -531,8 +537,9 @@ void updateparamflags() ...@@ -531,8 +537,9 @@ void updateparamflags()
FLAG(obj2.mxw) |= FLAG(obj2.myw) | FLAG(obj2.mx2w) | FLAG(obj2.alphas) FLAG(obj2.mxw) |= FLAG(obj2.myw) | FLAG(obj2.mx2w) | FLAG(obj2.alphas)
| FLAG(obj2.poserr_mx2w); | FLAG(obj2.poserr_mx2w);
FLAG(obj2.mamaposx) |= FLAG(obj2.mamaposy); 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) FLAG(obj2.flux_win) |= FLAG(obj2.mag_win)|FLAG(obj2.magerr_win)
| FLAG(obj2.flux_win) | FLAG(obj2.fluxerr_win); | FLAG(obj2.flux_win) | FLAG(obj2.fluxerr_win);
FLAG(obj2.winpos_x) |= FLAG(obj2.winpos_y) FLAG(obj2.winpos_x) |= FLAG(obj2.winpos_y)
| FLAG(obj2.winposerr_mx2) | FLAG(obj2.win_mx2) | FLAG(obj2.winposerr_mx2) | FLAG(obj2.win_mx2)
| FLAG(obj2.winpos_xw) | FLAG(obj2.winpos_xf) | FLAG(obj2.winpos_xw) | FLAG(obj2.winpos_xf)
......
...@@ -22,7 +22,7 @@ ...@@ -22,7 +22,7 @@
* 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: 11/03/2011 * Last modified: 25/03/2011
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
...@@ -208,7 +208,6 @@ void reinitcheck(picstruct *field, checkstruct *check) ...@@ -208,7 +208,6 @@ void reinitcheck(picstruct *field, checkstruct *check)
check->y = 0; check->y = 0;
fitshead = tab->headbuf; fitshead = tab->headbuf;
/* Neutralize possible scaling factors */ /* Neutralize possible scaling factors */
tab->bytepix = 4;
tab->bscale = 1.0; tab->bscale = 1.0;
tab->bzero = 0.0; tab->bzero = 0.0;
fitswrite(fitshead, "BSCALE ", &tab->bscale, H_FLOAT, T_DOUBLE); fitswrite(fitshead, "BSCALE ", &tab->bscale, H_FLOAT, T_DOUBLE);
...@@ -227,7 +226,8 @@ void reinitcheck(picstruct *field, checkstruct *check) ...@@ -227,7 +226,8 @@ void reinitcheck(picstruct *field, checkstruct *check)
case CHECK_BACKGROUND: case CHECK_BACKGROUND:
case CHECK_FILTERED: case CHECK_FILTERED:
case CHECK_SUBTRACTED: case CHECK_SUBTRACTED:
tab->bitpix = -32; tab->bitpix = BP_FLOAT;
tab->bytepix = 4;
tab->bitsgn = 0; tab->bitsgn = 0;
tab->naxisn[0] = check->width = field->width; tab->naxisn[0] = check->width = field->width;
tab->naxisn[1] = check->height = field->height; tab->naxisn[1] = check->height = field->height;
...@@ -239,7 +239,8 @@ void reinitcheck(picstruct *field, checkstruct *check) ...@@ -239,7 +239,8 @@ void reinitcheck(picstruct *field, checkstruct *check)
case CHECK_BACKRMS: case CHECK_BACKRMS:
case CHECK_SUBOBJECTS: case CHECK_SUBOBJECTS:
tab->bitpix = -32; tab->bitpix = BP_FLOAT;
tab->bytepix = 4;
tab->bitsgn = 0; tab->bitsgn = 0;
tab->naxisn[0] = check->width = field->width; tab->naxisn[0] = check->width = field->width;
tab->naxisn[1] = check->height = field->height; tab->naxisn[1] = check->height = field->height;
...@@ -266,6 +267,7 @@ void reinitcheck(picstruct *field, checkstruct *check) ...@@ -266,6 +267,7 @@ void reinitcheck(picstruct *field, checkstruct *check)
case CHECK_SUBDISKS: case CHECK_SUBDISKS:
case CHECK_PATTERNS: case CHECK_PATTERNS:
tab->bitpix = -32; tab->bitpix = -32;
tab->bytepix = 4;
tab->bitsgn = 0; tab->bitsgn = 0;
tab->naxisn[0] = check->width = field->width; tab->naxisn[0] = check->width = field->width;
tab->naxisn[1] = check->height = field->height; tab->naxisn[1] = check->height = field->height;
...@@ -276,7 +278,8 @@ void reinitcheck(picstruct *field, checkstruct *check) ...@@ -276,7 +278,8 @@ void reinitcheck(picstruct *field, checkstruct *check)
break; break;
case CHECK_SEGMENTATION: case CHECK_SEGMENTATION:
tab->bitpix = 32; tab->bitpix = BP_LONG;
tab->bytepix = 4;
tab->bitsgn = 1; tab->bitsgn = 1;
tab->naxisn[0] = check->width = field->width; tab->naxisn[0] = check->width = field->width;
tab->naxisn[1] = check->height = field->height; tab->naxisn[1] = check->height = field->height;
...@@ -285,8 +288,23 @@ void reinitcheck(picstruct *field, checkstruct *check) ...@@ -285,8 +288,23 @@ void reinitcheck(picstruct *field, checkstruct *check)
save_head(cat, cat->tab); save_head(cat, cat->tab);
break; break;
case CHECK_MASK:
case CHECK_SUBMASK:
tab->bitpix = BP_BYTE;
tab->bytepix = 1;
tab->bitsgn = 1;
tab->naxisn[0] = check->width = field->width;
tab->naxisn[1] = check->height = field->height;
check->npix = field->npix;
save_head(cat, cat->tab);
/*---- Allocate memory */
if (!check->line)
QMALLOC(check->line, FLAGTYPE, check->width);
break;
case CHECK_ASSOC: case CHECK_ASSOC:
tab->bitpix = -32; tab->bitpix = BP_FLOAT;
tab->bytepix = 4;
tab->bitsgn = 0; tab->bitsgn = 0;
tab->naxisn[0] = check->width = field->width; tab->naxisn[0] = check->width = field->width;
tab->naxisn[1] = check->height = field->height; tab->naxisn[1] = check->height = field->height;
...@@ -299,7 +317,8 @@ void reinitcheck(picstruct *field, checkstruct *check) ...@@ -299,7 +317,8 @@ void reinitcheck(picstruct *field, checkstruct *check)
case CHECK_MINIBACKGROUND: case CHECK_MINIBACKGROUND:
case CHECK_MINIBACKRMS: case CHECK_MINIBACKRMS:
tab->bitpix = -32; tab->bitpix = BP_FLOAT;
tab->bytepix = 4;
tab->bitsgn = 0; tab->bitsgn = 0;
tab->naxisn[0] = check->width = field->nbackx; tab->naxisn[0] = check->width = field->nbackx;
tab->naxisn[1] = check->height = field->nbacky; tab->naxisn[1] = check->height = field->nbacky;
...@@ -336,7 +355,8 @@ void reinitcheck(picstruct *field, checkstruct *check) ...@@ -336,7 +355,8 @@ void reinitcheck(picstruct *field, checkstruct *check)
break; break;
case CHECK_MAPSOM: case CHECK_MAPSOM:
tab->bitpix = -32; tab->bitpix = BP_FLOAT;
tab->bytepix = 4;
tab->bitsgn = 0; tab->bitsgn = 0;
tab->naxisn[0] = check->width = field->width; tab->naxisn[0] = check->width = field->width;
tab->naxisn[1] = check->height = field->height; tab->naxisn[1] = check->height = field->height;
...@@ -349,10 +369,11 @@ void reinitcheck(picstruct *field, checkstruct *check) ...@@ -349,10 +369,11 @@ void reinitcheck(picstruct *field, checkstruct *check)
break; break;
case CHECK_OTHER: case CHECK_OTHER:
tab->bitpix = -32; tab->bitpix = BP_FLOAT;
tab->bytepix = 4;
tab->bitsgn = 0; tab->bitsgn = 0;
tab->naxisn[0] = check->width; tab->naxisn[0] = check->width = field->width;
tab->naxisn[1] = check->height; tab->naxisn[1] = check->height = field->height;
check->npix = check->width*check->height; check->npix = check->width*check->height;
QCALLOC(check->pix, PIXTYPE, check->npix); QCALLOC(check->pix, PIXTYPE, check->npix);
save_head(cat, cat->tab); save_head(cat, cat->tab);
...@@ -376,7 +397,7 @@ void writecheck(checkstruct *check, PIXTYPE *data, int w) ...@@ -376,7 +397,7 @@ void writecheck(checkstruct *check, PIXTYPE *data, int w)
if (check->type == CHECK_APERTURES || check->type == CHECK_SUBPSFPROTOS if (check->type == CHECK_APERTURES || check->type == CHECK_SUBPSFPROTOS
|| check->type == CHECK_SUBPCPROTOS || check->type == CHECK_PCOPROTOS || check->type == CHECK_SUBPCPROTOS || check->type == CHECK_PCOPROTOS
|| check->type == CHECK_SUBPROFILES || check->type == CHECK_SUBSPHEROIDS || check->type == CHECK_SUBPROFILES || check->type == CHECK_SUBSPHEROIDS
|| check->type == CHECK_SUBDISKS) || check->type == CHECK_SUBDISKS || check->type == CHECK_OTHER)
{ {
memcpy((PIXTYPE *)check->pix + w*(check->y++), data, w*sizeof(PIXTYPE)); memcpy((PIXTYPE *)check->pix + w*(check->y++), data, w*sizeof(PIXTYPE));
return; return;
...@@ -386,13 +407,33 @@ void writecheck(checkstruct *check, PIXTYPE *data, int w) ...@@ -386,13 +407,33 @@ void writecheck(checkstruct *check, PIXTYPE *data, int w)
int i; int i;
PIXTYPE *pixt; PIXTYPE *pixt;
pixt = check->line; pixt = (PIXTYPE *)check->line;
for (i=w; i--; data++) for (i=w; i--; data++)
*(pixt++) = (*data>-BIG)? *data:0.0; *(pixt++) = (*data>-BIG)? *data:0.0;
data = check->line; write_body(check->cat->tab, (PIXTYPE *)check->line, w);
}
else if (check->type == CHECK_MASK)
{
int i;
FLAGTYPE *pixt;
pixt = (FLAGTYPE *)check->line;
for (i=w; i--;)
*(pixt++) = (*(data++)>-BIG)?0:1;
write_ibody(check->cat->tab, (FLAGTYPE *)check->line, w);
} }
else if (check->type == CHECK_SUBMASK)
{
int i;
FLAGTYPE *pixt;
write_body(check->cat->tab, data, w); pixt = (FLAGTYPE *)check->line;
for (i=w; i--;)
*(pixt++) = (*(data++)>-BIG)?1:0;
write_ibody(check->cat->tab, (FLAGTYPE *)check->line, w);
}
else
write_body(check->cat->tab, data, w);
return; return;
} }
...@@ -452,6 +493,19 @@ void reendcheck(picstruct *field, checkstruct *check) ...@@ -452,6 +493,19 @@ void reendcheck(picstruct *field, checkstruct *check)
pad_tab(cat, check->npix*sizeof(FLAGTYPE)); pad_tab(cat, check->npix*sizeof(FLAGTYPE));
break; break;
case CHECK_MASK:
case CHECK_SUBMASK:
{
int y;
for (y=field->ymin; y<field->ymax; y++)
writecheck(check, &PIX(field, 0, y), field->width);
free(check->line);
check->line = NULL;
pad_tab(cat, check->npix*sizeof(unsigned char));
break;
}
case CHECK_SUBOBJECTS: case CHECK_SUBOBJECTS:
{ {
int y; int y;
......
...@@ -22,7 +22,7 @@ ...@@ -22,7 +22,7 @@
* 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: 02/12/2010 * Last modified: 25/03/2011
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
...@@ -41,7 +41,7 @@ typedef struct structcheck ...@@ -41,7 +41,7 @@ typedef struct structcheck
size_t npix; /* number of pixels in check-image */ size_t npix; /* number of pixels in check-image */
int y; /* current line in check-image */ int y; /* current line in check-image */
PIXTYPE overlay; /* intensity of the overlayed plots */ PIXTYPE overlay; /* intensity of the overlayed plots */
PIXTYPE *line; /* buffered image line */ void *line; /* buffered image line */
checkenum type; /* CHECKIMAGE_TYPE */ checkenum type; /* CHECKIMAGE_TYPE */
} checkstruct; } checkstruct;
......
...@@ -22,7 +22,7 @@ ...@@ -22,7 +22,7 @@
* 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: 09/01/2011 * Last modified: 12/04/2011
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
...@@ -65,10 +65,10 @@ ...@@ -65,10 +65,10 @@
#define OUTPUT stderr /* where all msgs are sent */ #define OUTPUT stderr /* where all msgs are sent */
#define PSF_NPSFMAX 9 /* Max number of fitted PSFs */ #define PSF_NPSFMAX 9 /* Max number of fitted PSFs */
#define DEG (PI/180.0) /* 1 deg in radians */
#ifndef PI #ifndef PI
#define PI 3.1415926535898 /* never met before? */ #define PI 3.1415926535898 /* never met before? */
#endif #endif
#define DEG (PI/180.0) /* 1 deg in radians */
/* NOTES: /* NOTES:
* *
......
...@@ -22,7 +22,7 @@ ...@@ -22,7 +22,7 @@
* 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: 09/03/2011 * Last modified: 22/04/2011
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
...@@ -57,7 +57,7 @@ ...@@ -57,7 +57,7 @@
static int selectext(char *filename); static int selectext(char *filename);
time_t thetimet, thetimet2; time_t thetimet, thetimet2;
extern profitstruct *theprofit; extern profitstruct *theprofit,*thepprofit,*theqprofit;
extern char profname[][32]; extern char profname[][32];
double dtime; double dtime;
...@@ -124,7 +124,15 @@ void makeit() ...@@ -124,7 +124,15 @@ void makeit()
fft_init(prefs.nthreads); fft_init(prefs.nthreads);
/* Create profiles at full resolution */ /* Create profiles at full resolution */
NFPRINTF(OUTPUT, "Preparing profile models"); NFPRINTF(OUTPUT, "Preparing profile models");
theprofit = profit_init(thepsf); theprofit = profit_init(thepsf,
(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));
changecatparamarrays("VECTOR_MODEL", &theprofit->nparam, 1); changecatparamarrays("VECTOR_MODEL", &theprofit->nparam, 1);
changecatparamarrays("VECTOR_MODELERR", &theprofit->nparam, 1); changecatparamarrays("VECTOR_MODELERR", &theprofit->nparam, 1);
nparam2[0] = nparam2[1] = theprofit->nparam; nparam2[0] = nparam2[1] = theprofit->nparam;
...@@ -161,9 +169,12 @@ void makeit() ...@@ -161,9 +169,12 @@ void makeit()
{ {
if (i) if (i)
QPRINTF(OUTPUT, "+"); QPRINTF(OUTPUT, "+");
QPRINTF(OUTPUT, "%s", profname[theprofit->prof[i]->code]); QPRINTF(OUTPUT, "%s", theprofit->prof[i]->name);
} }
QPRINTF(OUTPUT, "\n"); QPRINTF(OUTPUT, "\n");
if (FLAG(obj2.prof_concentration)|FLAG(obj2.prof_concentration))
thepprofit = profit_init(thepsf, MODEL_DIRAC);
theqprofit = profit_init(thepsf, MODEL_EXPONENTIAL);
#else #else
error(EXIT_FAILURE, error(EXIT_FAILURE,
"*Error*: model-fitting is not supported in this build.\n", "*Error*: model-fitting is not supported in this build.\n",
...@@ -535,6 +546,11 @@ void makeit() ...@@ -535,6 +546,11 @@ void makeit()
if (prefs.prof_flag) if (prefs.prof_flag)
{ {
profit_end(theprofit); profit_end(theprofit);
if (FLAG(obj2.prof_concentration)|FLAG(obj2.prof_concentration))
{
profit_end(thepprofit);
profit_end(theqprofit);
}
fft_end(); fft_end();
} }
#endif #endif
......
...@@ -22,7 +22,7 @@ ...@@ -22,7 +22,7 @@
* 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: 13/03/2011 * Last modified: 03/05/2011
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
...@@ -127,6 +127,9 @@ keystruct objkey[] = { ...@@ -127,6 +127,9 @@ keystruct objkey[] = {
{"MAGERR_WIN", "RMS error for MAG_WIN", {"MAGERR_WIN", "RMS error for MAG_WIN",
&outobj2.magerr_win, H_FLOAT, T_FLOAT, "%8.4f", "mag", &outobj2.magerr_win, H_FLOAT, T_FLOAT, "%8.4f", "mag",
"stat.stdev;phot.mag", "mag"}, "stat.stdev;phot.mag", "mag"},
{"SNR_WIN", "Gaussian-weighted SNR",
&outobj2.snr_win, H_FLOAT, T_FLOAT, "%10.4g", "",
"stat.snr", ""},
{"FLUX_SOMFIT", "Flux derived from SOM fit", {"FLUX_SOMFIT", "Flux derived from SOM fit",
&outobj2.flux_somfit, H_FLOAT, T_FLOAT, "%12.7g", "count", &outobj2.flux_somfit, H_FLOAT, T_FLOAT, "%12.7g", "count",
......
...@@ -22,7 +22,7 @@ ...@@ -22,7 +22,7 @@
* 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: 13/03/2011 * Last modified: 07/04/2011
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
...@@ -60,16 +60,16 @@ ...@@ -60,16 +60,16 @@
&outobj2.magerr_prof, H_FLOAT, T_FLOAT, "%8.4f", "mag", &outobj2.magerr_prof, H_FLOAT, T_FLOAT, "%8.4f", "mag",
"stat.error;phot.mag;stat.fit.param", "mag"}, "stat.error;phot.mag;stat.fit.param", "mag"},
{"FLUXHYBRID_MODEL", "Hybrid flux from model-fitting", {"FLUX_HYBRID", "Hybrid flux from model-fitting",
&outobj2.fluxcor_prof, H_FLOAT, T_FLOAT, "%12.7g", "count", &outobj2.fluxcor_prof, H_FLOAT, T_FLOAT, "%12.7g", "count",
"phot.count;stat.fit.param", "ct"}, "phot.count;stat.fit.param", "ct"},
{"FLUXHYBRIDERR_MODEL", "RMS error on hybrid flux", {"FLUXERR_HYBRID", "RMS error on hybrid flux",
&outobj2.fluxcorerr_prof, H_FLOAT, T_FLOAT, "%12.7g", "count", &outobj2.fluxcorerr_prof, H_FLOAT, T_FLOAT, "%12.7g", "count",
"stat.error;phot.count;stat.fit.param", "ct"}, "stat.error;phot.count;stat.fit.param", "ct"},
{"MAGHYBRID_MODEL", "Hybrid magnitude from model-fitting", {"MAG_HYBRID", "Hybrid magnitude from model-fitting",
&outobj2.magcor_prof, H_FLOAT, T_FLOAT, "%8.4f", "mag", &outobj2.magcor_prof, H_FLOAT, T_FLOAT, "%8.4f", "mag",
"phot.mag;stat.fit.param", "mag"}, "phot.mag;stat.fit.param", "mag"},
{"MAGHYBRIDERR_MODEL", "RMS error on hybrid magnitude", {"MAGERR_HYBRID", "RMS error on hybrid magnitude",
&outobj2.magcorerr_prof, H_FLOAT, T_FLOAT, "%8.4f", "mag", &outobj2.magcorerr_prof, H_FLOAT, T_FLOAT, "%8.4f", "mag",
"stat.error;phot.mag;stat.fit.param", "mag"}, "stat.error;phot.mag;stat.fit.param", "mag"},
...@@ -338,6 +338,11 @@ ...@@ -338,6 +338,11 @@
{"SPREADERR_MODEL", "Spread parameter error from model-fitting", {"SPREADERR_MODEL", "Spread parameter error from model-fitting",
&outobj2.prof_concentrationerr, H_FLOAT, T_FLOAT, "%8.5f", "", &outobj2.prof_concentrationerr, H_FLOAT, T_FLOAT, "%8.5f", "",
"src.morph.param", ""}, "src.morph.param", ""},
{"NOISEAREA_MODEL", "Equivalent noise area of the fitted model",
&outobj2.prof_noisearea, H_FLOAT, T_FLOAT, "%12.2f", "pixel**2",
"phys.area", "pix2"},
/* /*
{"CLASS_STAR_MODEL", "S/G classifier from model-fitting", {"CLASS_STAR_MODEL", "S/G classifier from model-fitting",
&outobj2.prof_class_star, H_FLOAT, T_FLOAT, "%7.4f", "", &outobj2.prof_class_star, H_FLOAT, T_FLOAT, "%7.4f", "",
......
...@@ -22,7 +22,7 @@ ...@@ -22,7 +22,7 @@
* 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: 09/03/2011 * Last modified: 25/03/2011
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
...@@ -81,7 +81,7 @@ ...@@ -81,7 +81,7 @@
{"NONE", "IDENTICAL", {"NONE", "IDENTICAL",
"BACKGROUND", "BACKGROUND_RMS", "MINIBACKGROUND", "BACKGROUND", "BACKGROUND_RMS", "MINIBACKGROUND",
"MINIBACK_RMS", "-BACKGROUND", "MINIBACK_RMS", "-BACKGROUND",
"FILTERED", "OBJECTS", "APERTURES", "SEGMENTATION", "ASSOC", "FILTERED", "OBJECTS", "APERTURES", "SEGMENTATION", "MASK", "-MASK", "ASSOC",
"-OBJECTS", "PSFS", "-PSFS", "-OBJECTS", "PSFS", "-PSFS",
"PC_CONVPROTOS", "-PC_CONVPROTOS", "PC_PROTOS", "MAP_SOM", "PC_CONVPROTOS", "-PC_CONVPROTOS", "PC_PROTOS", "MAP_SOM",
#ifdef USE_MODEL #ifdef USE_MODEL
......
...@@ -22,7 +22,7 @@ ...@@ -22,7 +22,7 @@
* 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: 11/03/2011 * Last modified: 29/04/2011
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
...@@ -62,12 +62,6 @@ static void make_kernel(float pos, float *kernel, interpenum interptype); ...@@ -62,12 +62,6 @@ static void make_kernel(float pos, float *kernel, interpenum interptype);
/*------------------------------- variables ---------------------------------*/ /*------------------------------- variables ---------------------------------*/
const char profname[][32]={"background offset", "point source",
"Sersic spheroid", "de Vaucouleurs spheroid",
"exponential disk", "spiral arms",
"bar", "inner ring", "outer ring", "tabulated model",
""};
const int interp_kernwidth[5]={1,2,4,6,8}; const int interp_kernwidth[5]={1,2,4,6,8};
const int flux_flag[PARAM_NPARAM] = {0, const int flux_flag[PARAM_NPARAM] = {0,
...@@ -83,71 +77,30 @@ const int flux_flag[PARAM_NPARAM] = {0, ...@@ -83,71 +77,30 @@ const int flux_flag[PARAM_NPARAM] = {0,
/* "Local" global variables for debugging purposes */ /* "Local" global variables for debugging purposes */
int theniter, the_gal; int theniter, the_gal;
static picstruct *the_field, *the_wfield; static picstruct *the_field, *the_wfield;
profitstruct *theprofit; profitstruct *theprofit, *thepprofit, *theqprofit;
/****** profit_init *********************************************************** /****** profit_init ***********************************************************
PROTO profitstruct profit_init(psfstruct *psf) PROTO profitstruct profit_init(psfstruct *psf, unsigned int modeltype)
PURPOSE Allocate and initialize a new profile-fitting structure. PURPOSE Allocate and initialize a new profile-fitting structure.
INPUT Pointer to PSF structure. INPUT Pointer to PSF structure,
Model type.
OUTPUT A pointer to an allocated profit structure. OUTPUT A pointer to an allocated profit structure.
NOTES -. NOTES -.
AUTHOR E. Bertin (IAP) AUTHOR E. Bertin (IAP)
VERSION 08/10/2010 VERSION 22/04/2011
***/ ***/
profitstruct *profit_init(psfstruct *psf) profitstruct *profit_init(psfstruct *psf, unsigned int modeltype)
{ {
profitstruct *profit; profitstruct *profit;
int p, nprof, int t, nmodels;
backflag, diracflag, spheroidflag, diskflag,
barflag, armsflag;
QCALLOC(profit, profitstruct, 1); QCALLOC(profit, profitstruct, 1);
profit->psf = psf; profit->psf = psf;
profit->psfdft = NULL; QMALLOC(profit->prof, profstruct *, MODEL_NMAX);
nmodels = 0;
profit->nparam = 0; for (t=1; t<(1<<MODEL_NMAX); t<<=1)
QMALLOC(profit->prof, profstruct *, PROF_NPROF); if (modeltype&t)
backflag = diracflag = spheroidflag = diskflag = barflag = armsflag = 0; profit->prof[nmodels++] = prof_init(profit, t);
nprof = 0;
for (p=0; p<PROF_NPROF; p++)
if (!backflag && FLAG(obj2.prof_offset_flux))
{
profit->prof[p] = prof_init(profit, PROF_BACK);
backflag = 1;
nprof++;
}
else if (!diracflag && FLAG(obj2.prof_dirac_flux))
{
profit->prof[p] = prof_init(profit, PROF_DIRAC);
diracflag = 1;
nprof++;
}
else if (!spheroidflag && FLAG(obj2.prof_spheroid_flux))
{
profit->prof[p] = prof_init(profit,
FLAG(obj2.prof_spheroid_sersicn)? PROF_SERSIC : PROF_DEVAUCOULEURS);
spheroidflag = 1;
nprof++;
}
else if (!diskflag && FLAG(obj2.prof_disk_flux))
{
profit->prof[p] = prof_init(profit, PROF_EXPONENTIAL);
diskflag = 1;
nprof++;
}
else if (diskflag && !barflag && FLAG(obj2.prof_bar_flux))
{
profit->prof[p] = prof_init(profit, PROF_BAR);
barflag = 1;
nprof++;
}
else if (barflag && !armsflag && FLAG(obj2.prof_arms_flux))
{
profit->prof[p] = prof_init(profit, PROF_ARMS);
armsflag = 1;
nprof++;
}
/* Allocate memory for the complete model */ /* Allocate memory for the complete model */
QMALLOC16(profit->modpix, float, PROFIT_MAXMODSIZE*PROFIT_MAXMODSIZE); QMALLOC16(profit->modpix, float, PROFIT_MAXMODSIZE*PROFIT_MAXMODSIZE);
QMALLOC16(profit->modpix2, float, PROFIT_MAXMODSIZE*PROFIT_MAXMODSIZE); QMALLOC16(profit->modpix2, float, PROFIT_MAXMODSIZE*PROFIT_MAXMODSIZE);
...@@ -159,7 +112,7 @@ profitstruct *profit_init(psfstruct *psf) ...@@ -159,7 +112,7 @@ profitstruct *profit_init(psfstruct *psf)
QMALLOC16(profit->lmodpix2, PIXTYPE, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE); QMALLOC16(profit->lmodpix2, PIXTYPE, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE);
QMALLOC16(profit->resi, float, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE); QMALLOC16(profit->resi, float, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE);
QMALLOC16(profit->covar, float, profit->nparam*profit->nparam); QMALLOC16(profit->covar, float, profit->nparam*profit->nparam);
profit->nprof = nprof; profit->nprof = nmodels;
profit->fluxfac = 1.0; /* Default */ profit->fluxfac = 1.0; /* Default */
return profit; return profit;
...@@ -213,17 +166,19 @@ OUTPUT Pointer to an allocated fit structure (containing details about the ...@@ -213,17 +166,19 @@ OUTPUT Pointer to an allocated fit structure (containing details about the
fit). fit).
NOTES It is a modified version of the lm_minimize() of lmfit. NOTES It is a modified version of the lm_minimize() of lmfit.
AUTHOR E. Bertin (IAP) AUTHOR E. Bertin (IAP)
VERSION 07/10/2010 VERSION 29/04/2011
***/ ***/
void profit_fit(profitstruct *profit, void profit_fit(profitstruct *profit,
picstruct *field, picstruct *wfield, picstruct *field, picstruct *wfield,
objstruct *obj, obj2struct *obj2) objstruct *obj, obj2struct *obj2)
{ {
profitstruct pprofit; profitstruct *pprofit, *qprofit;
patternstruct *pattern; patternstruct *pattern;
psfstruct *psf; psfstruct *psf;
checkstruct *check; checkstruct *check;
double emx2,emy2,emxy, a , cp,sp, cn, bn, n; double emx2,emy2,emxy, a , cp,sp, cn, bn, n,
sump,sumq, sumpw2,sumqw2,sumpqw, sump0,sumq0;
PIXTYPE valp,valq,sig2;
float param0[PARAM_NPARAM], param1[PARAM_NPARAM], float param0[PARAM_NPARAM], param1[PARAM_NPARAM],
param[PARAM_NPARAM], param[PARAM_NPARAM],
**list, **list,
...@@ -498,6 +453,7 @@ profit->niter = profit_minimize(profit, PROFIT_MAXITER); ...@@ -498,6 +453,7 @@ profit->niter = profit_minimize(profit, PROFIT_MAXITER);
obj2->prof_chi2 = (profit->nresi > profit->nparam)? obj2->prof_chi2 = (profit->nresi > profit->nparam)?
profit->chi2 / (profit->nresi - profit->nparam) : 0.0; profit->chi2 / (profit->nresi - profit->nparam) : 0.0;
/* Position */
if (FLAG(obj2.x_prof)) if (FLAG(obj2.x_prof))
{ {
i = profit->paramindex[PARAM_X]; i = profit->paramindex[PARAM_X];
...@@ -552,15 +508,25 @@ profit->niter = profit_minimize(profit, PROFIT_MAXITER); ...@@ -552,15 +508,25 @@ profit->niter = profit_minimize(profit, PROFIT_MAXITER);
} }
} }
/* Equivalent noise area */
if (FLAG(obj2.prof_noisearea))
obj2->prof_noisearea = profit_noisearea(profit);
/* Second order moments and ellipticities */
if (FLAG(obj2.prof_mx2))
profit_moments(profit, obj2);
/* Second order moments of the convolved model (used by other parameters) */
if (FLAG(obj2.prof_convmx2))
profit_convmoments(profit, obj2);
/* "Hybrid" magnitudes */
if (FLAG(obj2.fluxcor_prof)) if (FLAG(obj2.fluxcor_prof))
{ {
profit_residuals(profit,field,wfield, 0.0, profit->paraminit, NULL); profit_residuals(profit,field,wfield, 0.0, profit->paraminit, NULL);
profit_fluxcor(profit, obj, obj2); profit_fluxcor(profit, obj, obj2);
} }
if (FLAG(obj2.prof_mx2))
profit_moments(profit, obj2);
/* Do measurements on the rasterised model (surface brightnesses) */ /* Do measurements on the rasterised model (surface brightnesses) */
if (FLAG(obj2.peak_prof)) if (FLAG(obj2.peak_prof))
profit_surface(profit, obj2); profit_surface(profit, obj2);
...@@ -770,48 +736,108 @@ profit->niter = profit_minimize(profit, PROFIT_MAXITER); ...@@ -770,48 +736,108 @@ profit->niter = profit_minimize(profit, PROFIT_MAXITER);
/* Star/galaxy classification */ /* Star/galaxy classification */
if (FLAG(obj2.prof_class_star) || FLAG(obj2.prof_concentration)) if (FLAG(obj2.prof_class_star) || FLAG(obj2.prof_concentration))
{ {
pprofit = *profit; profit_residuals(profit,field,wfield, PROFIT_DYNPARAM, profit->paraminit,
memset(pprofit.paramindex, 0, PARAM_NPARAM*sizeof(int)); profit->resi);
memset(pprofit.paramlist, 0, PARAM_NPARAM*sizeof(float *)); pprofit = thepprofit;
pprofit.nparam = 0; nparam = pprofit->nparam;
QMALLOC(pprofit.prof, profstruct *, 1); if (pprofit->psfdft)
pprofit.prof[0] = prof_init(&pprofit, PROF_DIRAC); {
QMALLOC16(pprofit.covar, float, pprofit.nparam*pprofit.nparam); QFREE(pprofit->psfdft);
pprofit.nprof = 1; }
profit_resetparams(&pprofit); psf = pprofit->psf;
pprofit->pixstep = profit->pixstep;
pprofit->ix = profit->ix;
pprofit->iy = profit->iy;
pprofit->objnaxisn[0] = profit->objnaxisn[0];
pprofit->objnaxisn[1] = profit->objnaxisn[1];
pprofit->subsamp = profit->subsamp;
pprofit->nobjpix = profit->nobjpix;
pprofit->obj = obj;
pprofit->obj2 = obj2;
pprofit->nresi = profit_copyobjpix(pprofit, field, wfield);
pprofit->modnaxisn[0] = profit->modnaxisn[0];
pprofit->modnaxisn[1] = profit->modnaxisn[1];
pprofit->nmodpix = profit->nmodpix;
profit_psf(pprofit);
pprofit->sigma = obj->sigbkg;
for (p=0; p<pprofit->nparam; p++)
pprofit->freeparam_flag[p] = 1;
pprofit->nfreeparam = pprofit->nparam;
profit_resetparams(pprofit);
if (profit->paramlist[PARAM_X] && profit->paramlist[PARAM_Y]) if (profit->paramlist[PARAM_X] && profit->paramlist[PARAM_Y])
{ {
pprofit.paraminit[pprofit.paramindex[PARAM_X]] = *profit->paramlist[PARAM_X]; pprofit->paraminit[pprofit->paramindex[PARAM_X]] = *profit->paramlist[PARAM_X];
pprofit.paraminit[pprofit.paramindex[PARAM_Y]] = *profit->paramlist[PARAM_Y]; pprofit->paraminit[pprofit->paramindex[PARAM_Y]] = *profit->paramlist[PARAM_Y];
}
fft_reset();
pprofit->paraminit[pprofit->paramindex[PARAM_DIRAC_FLUX]] = profit->flux;
pprofit->niter = profit_minimize(pprofit, PROFIT_MAXITER);
profit_residuals(pprofit,field,wfield, PROFIT_DYNPARAM, pprofit->paraminit,
pprofit->resi);
qprofit = theqprofit;
nparam = qprofit->nparam;
if (qprofit->psfdft)
{
QFREE(qprofit->psfdft);
} }
pprofit.paraminit[pprofit.paramindex[PARAM_DIRAC_FLUX]] = profit->flux; qprofit->pixstep = profit->pixstep;
for (p=0; p<pprofit.nparam; p++) qprofit->ix = profit->ix;
pprofit.freeparam_flag[p] = 1; qprofit->iy = profit->iy;
pprofit.nfreeparam = pprofit.nparam; qprofit->objnaxisn[0] = profit->objnaxisn[0];
pprofit.niter = profit_minimize(&pprofit, PROFIT_MAXITER); qprofit->objnaxisn[1] = profit->objnaxisn[1];
profit_residuals(&pprofit,field,wfield, PROFIT_DYNPARAM, pprofit.paraminit, qprofit->subsamp = profit->subsamp;
pprofit.resi); qprofit->nobjpix = profit->nobjpix;
qprofit->obj = obj;
qprofit->obj2 = obj2;
qprofit->nresi = profit_copyobjpix(qprofit, field, wfield);
qprofit->modnaxisn[0] = profit->modnaxisn[0];
qprofit->modnaxisn[1] = profit->modnaxisn[1];
qprofit->nmodpix = profit->nmodpix;
profit_psf(qprofit);
qprofit->sigma = obj->sigbkg;
for (p=0; p<qprofit->nparam; p++)
qprofit->freeparam_flag[p] = 1;
qprofit->nfreeparam = qprofit->nparam;
profit_resetparams(qprofit);
fft_reset();
qprofit->paraminit[qprofit->paramindex[PARAM_X]] = pprofit->paraminit[pprofit->paramindex[PARAM_X]];
qprofit->paraminit[qprofit->paramindex[PARAM_Y]] = pprofit->paraminit[pprofit->paramindex[PARAM_Y]];
qprofit->paraminit[qprofit->paramindex[PARAM_DISK_FLUX]] = pprofit->paraminit[pprofit->paramindex[PARAM_DIRAC_FLUX]];
qprofit->paraminit[qprofit->paramindex[PARAM_DISK_SCALE]] = psf->fwhm/16.0;
qprofit->paraminit[qprofit->paramindex[PARAM_DISK_ASPECT]] = 1.0;
qprofit->paraminit[qprofit->paramindex[PARAM_DISK_POSANG]] = 0.0;
profit_residuals(qprofit,field,wfield, PROFIT_DYNPARAM, qprofit->paraminit,
qprofit->resi);
sump = sumq = sumpw2 = sumqw2 = sumpqw = sump0 = sumq0 = 0.0;
for (p=0; p<pprofit->nobjpix; p++)
if (pprofit->objweight[p]>0 && pprofit->objpix[p]>-BIG)
{
valp = pprofit->lmodpix[p];
sump += (double)(valp*pprofit->objpix[p]);
valq = qprofit->lmodpix[p];
sumq += (double)(valq*pprofit->objpix[p]);
sump0 += (double)(valp*valp);
sumq0 += (double)(valp*valq);
sig2 = 1.0f/(pprofit->objweight[p]*pprofit->objweight[p]);
sumpw2 += valp*valp*sig2;
sumqw2 += valq*valq*sig2;
sumpqw += valp*valq*sig2;
}
if (FLAG(obj2.prof_class_star)) if (FLAG(obj2.prof_class_star))
{ {
dchi2 = 0.5*(pprofit.chi2 - profit->chi2); dchi2 = 0.5*(pprofit->chi2 - profit->chi2);
obj2->prof_class_star = dchi2 < 50.0? obj2->prof_class_star = dchi2 < 50.0?
(dchi2 > -50.0? 2.0/(1.0+expf(dchi2)) : 2.0) : 0.0; (dchi2 > -50.0? 2.0/(1.0+expf(dchi2)) : 2.0) : 0.0;
} }
if (FLAG(obj2.prof_concentration)) if (FLAG(obj2.prof_concentration))
{ {
if (profit->flux > 0.0 && pprofit.flux > 0.0) obj2->prof_concentration = sump>0.0? (sumq/sump - sumq0/sump0) : 1.0;
obj2->prof_concentration = -2.5*log10(pprofit.flux / profit->flux);
else if (profit->flux > 0.0)
obj2->prof_concentration = 99.0;
else if (pprofit.flux > 0.0)
obj2->prof_concentration = -99.0;
if (FLAG(obj2.prof_concentrationerr)) if (FLAG(obj2.prof_concentrationerr))
obj2->prof_concentrationerr = (obj2->flux_prof > 0.0? obj2->prof_concentrationerr = sump>0.0?
1.086*(obj2->fluxerr_prof / obj2->flux_prof) : 99.0); sqrt(sumqw2*sump*sump+sumpw2*sumq*sumq-2.0*sumpqw*sump*sumq)
/ (sump*sump) : 0.0;
} }
prof_end(pprofit.prof[0]);
free(pprofit.prof);
free(pprofit.covar);
} }
/* clean up. */ /* clean up. */
...@@ -821,8 +847,36 @@ profit->niter = profit_minimize(profit, PROFIT_MAXITER); ...@@ -821,8 +847,36 @@ profit->niter = profit_minimize(profit, PROFIT_MAXITER);
} }
/****** profit_noisearea ******************************************************
PROTO float profit_noisearea(profitstruct *profit)
PURPOSE Return the equivalent noise area (see King 1983) of a model.
INPUT Profile-fitting structure,
OUTPUT Equivalent noise area, in pixels.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 19/10/2010
***/
float profit_noisearea(profitstruct *profit)
{
double dval, flux,flux2;
PIXTYPE *pix;
int p;
flux = flux2 = 0.0;
pix = profit->lmodpix;
for (p=profit->nobjpix; p--;)
{
dval = (double)*(pix++);
flux += dval;
flux2 += dval*dval;
}
return (float)(flux2>0.0? flux*flux / flux2 : 0.0);
}
/****** profit_fluxcor ****************************************************** /****** profit_fluxcor ******************************************************
PROTO void profit_fluxcor(profitstruct *profit objstruct *obj, PROTO void profit_fluxcor(profitstruct *profit, objstruct *obj,
obj2struct *obj2) obj2struct *obj2)
PURPOSE Integrate the flux within an ellipse and complete it with the wings of PURPOSE Integrate the flux within an ellipse and complete it with the wings of
the fitted model. the fitted model.
...@@ -832,15 +886,17 @@ INPUT Profile-fitting structure, ...@@ -832,15 +886,17 @@ INPUT Profile-fitting structure,
OUTPUT Model-corrected flux. OUTPUT Model-corrected flux.
NOTES -. NOTES -.
AUTHOR E. Bertin (IAP) AUTHOR E. Bertin (IAP)
VERSION 19/10/2010 VERSION 12/04/2011
***/ ***/
void profit_fluxcor(profitstruct *profit, objstruct *obj, obj2struct *obj2) void profit_fluxcor(profitstruct *profit, objstruct *obj, obj2struct *obj2)
{ {
double mx,my, dx,dy, cx2,cy2,cxy, klim2, tvobj,sigtvobj, checkstruct *check;
tvm,tvmin,tvmout; double mx,my, dx,dy, cx2,cy2,cxy, klim,klim2, tvobj,sigtvobj,
PIXTYPE *objpix,*objpixt,*objweight,*objweightt, *lmodpix, tvm,tvmin,tvmout, r1,v1;
PIXTYPE *objpix,*objpixt,*objweight,*objweightt, *lmodpix,
pix, weight,var; pix, weight,var;
int x,y, x2,y2, pos, w,h, area, corrflag; int x,y, x2,y2, pos, w,h, area, corrflag;
corrflag = (prefs.mask_type==MASK_CORRECT); corrflag = (prefs.mask_type==MASK_CORRECT);
w = profit->objnaxisn[0]; w = profit->objnaxisn[0];
...@@ -854,12 +910,13 @@ void profit_fluxcor(profitstruct *profit, objstruct *obj, obj2struct *obj2) ...@@ -854,12 +910,13 @@ void profit_fluxcor(profitstruct *profit, objstruct *obj, obj2struct *obj2)
if (profit->paramlist[PARAM_Y]) if (profit->paramlist[PARAM_Y])
my += *profit->paramlist[PARAM_Y]; my += *profit->paramlist[PARAM_Y];
} }
if (obj2->kronfactor>0.0) if (obj2->kronfactor>0.0)
{ {
cx2 = obj->cxx; cx2 = obj->cxx;
cy2 = obj->cyy; cy2 = obj->cyy;
cxy = obj->cxy; cxy = obj->cxy;
klim2 = 2.0; klim2 = 0.64*obj2->kronfactor*obj2->kronfactor;
} }
else else
/*-- ...if not, use the circular aperture provided by the user */ /*-- ...if not, use the circular aperture provided by the user */
...@@ -868,6 +925,33 @@ void profit_fluxcor(profitstruct *profit, objstruct *obj, obj2struct *obj2) ...@@ -868,6 +925,33 @@ void profit_fluxcor(profitstruct *profit, objstruct *obj, obj2struct *obj2)
cxy = 0.0; cxy = 0.0;
klim2 = (prefs.autoaper[1]/2.0)*(prefs.autoaper[1]/2.0); klim2 = (prefs.autoaper[1]/2.0)*(prefs.autoaper[1]/2.0);
} }
/*
cx2 = obj2->prof_convcxx;
cy2 = obj2->prof_convcyy;
cxy = obj2->prof_convcxy;
lmodpix = profit->lmodpix;
r1 = v1 = 0.0;
for (y=0; y<h; y++)
{
dy = y - my;
for (x=0; x<w; x++)
{
dx = x - mx;
pix = *(lmodpix++);
r1 += sqrt(cx2*dx*dx + cy2*dy*dy + cxy*dx*dy)*pix;
v1 += pix;
}
}
klim = r1/v1*2.0;
klim2 = klim*klim;
if ((check = prefs.check[CHECK_APERTURES]))
sexellips(check->pix, check->width, check->height,
obj2->x_prof-1.0, obj2->y_prof-1.0, klim*obj2->prof_conva,klim*obj2->prof_convb,
obj2->prof_convtheta, check->overlay, 0);
*/
area = 0; area = 0;
tvmin = tvmout = tvobj = sigtvobj = 0.0; tvmin = tvmout = tvobj = sigtvobj = 0.0;
...@@ -906,8 +990,8 @@ void profit_fluxcor(profitstruct *profit, objstruct *obj, obj2struct *obj2) ...@@ -906,8 +990,8 @@ void profit_fluxcor(profitstruct *profit, objstruct *obj, obj2struct *obj2)
} }
tvobj += pix; tvobj += pix;
sigtvobj += var; sigtvobj += var;
tvmin += *lmodpix; tvmin += *(lmodpix++);
*(lmodpix++) = pix; // *(lmodpix++) = pix;
} }
else else
tvmout += *(lmodpix++); tvmout += *(lmodpix++);
...@@ -929,6 +1013,10 @@ void profit_fluxcor(profitstruct *profit, objstruct *obj, obj2struct *obj2) ...@@ -929,6 +1013,10 @@ void profit_fluxcor(profitstruct *profit, objstruct *obj, obj2struct *obj2)
obj2->fluxcorerr_prof = sqrt(sigtvobj); obj2->fluxcorerr_prof = sqrt(sigtvobj);
} }
/*
if ((check = prefs.check[CHECK_OTHER]))
addcheck(check, profit->lmodpix, w, h, profit->ix,profit->iy, 1.0);
*/
return; return;
} }
...@@ -1041,17 +1129,17 @@ float profit_minradius(profitstruct *profit, float refffac) ...@@ -1041,17 +1129,17 @@ float profit_minradius(profitstruct *profit, float refffac)
{ {
switch (profit->prof[p]->code) switch (profit->prof[p]->code)
{ {
case PROF_BACK: case MODEL_BACK:
case PROF_DIRAC: case MODEL_DIRAC:
reff = 0.0; reff = 0.0;
break; break;
case PROF_SERSIC: case MODEL_SERSIC:
reff = *profit->paramlist[PARAM_SPHEROID_REFF]; reff = *profit->paramlist[PARAM_SPHEROID_REFF];
break; break;
case PROF_DEVAUCOULEURS: case MODEL_DEVAUCOULEURS:
reff = *profit->paramlist[PARAM_SPHEROID_REFF]; reff = *profit->paramlist[PARAM_SPHEROID_REFF];
break; break;
case PROF_EXPONENTIAL: case MODEL_EXPONENTIAL:
reff = *profit->paramlist[PARAM_DISK_SCALE]*1.67835; reff = *profit->paramlist[PARAM_DISK_SCALE]*1.67835;
break; break;
default: default:
...@@ -1266,7 +1354,7 @@ void profit_evaluate(double *dpar, double *fvec, int m, int n, void *adata) ...@@ -1266,7 +1354,7 @@ void profit_evaluate(double *dpar, double *fvec, int m, int n, void *adata)
jflag = 1; jflag = 1;
} }
jflag = 0; /* Temporarily deactivated (until problems are fixed) */ jflag = 0; /* Temporarily deactivated (until problems are fixed) */
if (jflag && !(profit->nprof==1 && profit->prof[0]->code == PROF_DIRAC)) if (jflag && !(profit->nprof==1 && profit->prof[0]->code == MODEL_DIRAC))
{ {
prof = profit->prof; prof = profit->prof;
nprof = profit->nprof; nprof = profit->nprof;
...@@ -1322,15 +1410,15 @@ jflag = 0; /* Temporarily deactivated (until problems are fixed) */ ...@@ -1322,15 +1410,15 @@ jflag = 0; /* Temporarily deactivated (until problems are fixed) */
case PARAM_SPHEROID_SERSICN: case PARAM_SPHEROID_SERSICN:
sflag = 0; /* We are in the same switch */ sflag = 0; /* We are in the same switch */
for (c=0; c<nprof; c++) for (c=0; c<nprof; c++)
if (prof[c]->code == PROF_SERSIC if (prof[c]->code == MODEL_SERSIC
|| prof[c]->code == PROF_DEVAUCOULEURS) || prof[c]->code == MODEL_DEVAUCOULEURS)
break; break;
case PARAM_DISK_SCALE: case PARAM_DISK_SCALE:
case PARAM_DISK_ASPECT: case PARAM_DISK_ASPECT:
case PARAM_DISK_POSANG: case PARAM_DISK_POSANG:
if (sflag) if (sflag)
for (c=0; c<nprof; c++) for (c=0; c<nprof; c++)
if (prof[c]->code == PROF_EXPONENTIAL) if (prof[c]->code == MODEL_EXPONENTIAL)
break; break;
sflag = 0; sflag = 0;
case PARAM_ARMS_QUADFRAC: case PARAM_ARMS_QUADFRAC:
...@@ -1342,14 +1430,14 @@ jflag = 0; /* Temporarily deactivated (until problems are fixed) */ ...@@ -1342,14 +1430,14 @@ jflag = 0; /* Temporarily deactivated (until problems are fixed) */
case PARAM_ARMS_WIDTH: case PARAM_ARMS_WIDTH:
if (sflag) if (sflag)
for (c=0; c<nprof; c++) for (c=0; c<nprof; c++)
if (prof[c]->code == PROF_ARMS) if (prof[c]->code == MODEL_ARMS)
break; break;
sflag = 0; sflag = 0;
case PARAM_BAR_ASPECT: case PARAM_BAR_ASPECT:
case PARAM_BAR_POSANG: case PARAM_BAR_POSANG:
if (sflag) if (sflag)
for (c=0; c<nprof; c++) for (c=0; c<nprof; c++)
if (prof[c]->code == PROF_ARMS) if (prof[c]->code == MODEL_ARMS)
break; break;
modpixt = profit->modpix; modpixt = profit->modpix;
profpixt = prof[c]->pix; profpixt = prof[c]->pix;
...@@ -1433,7 +1521,7 @@ float *profit_residuals(profitstruct *profit, picstruct *field, ...@@ -1433,7 +1521,7 @@ float *profit_residuals(profitstruct *profit, picstruct *field,
for (p=0; p<profit->nparam; p++) for (p=0; p<profit->nparam; p++)
profit->param[p] = param[p]; profit->param[p] = param[p];
/* Simple PSF shortcut */ /* Simple PSF shortcut */
if (profit->nprof == 1 && profit->prof[0]->code == PROF_DIRAC) if (profit->nprof == 1 && profit->prof[0]->code == MODEL_DIRAC)
{ {
profit_resample(profit, profit->psfpix, profit->lmodpix, profit_resample(profit, profit->psfpix, profit->lmodpix,
*profit->prof[0]->flux); *profit->prof[0]->flux);
...@@ -2232,7 +2320,7 @@ INPUT Profile-fitting structure, ...@@ -2232,7 +2320,7 @@ INPUT Profile-fitting structure,
OUTPUT -. OUTPUT -.
NOTES -. NOTES -.
AUTHOR E. Bertin (IAP) AUTHOR E. Bertin (IAP)
VERSION 20/08/2010 VERSION 22/04/2011
***/ ***/
void profit_moments(profitstruct *profit, obj2struct *obj2) void profit_moments(profitstruct *profit, obj2struct *obj2)
{ {
...@@ -2244,7 +2332,7 @@ void profit_moments(profitstruct *profit, obj2struct *obj2) ...@@ -2244,7 +2332,7 @@ void profit_moments(profitstruct *profit, obj2struct *obj2)
temp, temp2,invtemp2,invstemp2, temp, temp2,invtemp2,invstemp2,
pmx2,theta, flux, dval; pmx2,theta, flux, dval;
float *covart; float *covart;
int findex[PROF_NPROF], int findex[MODEL_NMAX],
i,j,p, nparam; i,j,p, nparam;
/* hw = (float)(profit->modnaxisn[0]/2);*/ /* hw = (float)(profit->modnaxisn[0]/2);*/
...@@ -2445,6 +2533,95 @@ void profit_moments(profitstruct *profit, obj2struct *obj2) ...@@ -2445,6 +2533,95 @@ void profit_moments(profitstruct *profit, obj2struct *obj2)
} }
/****** profit_convmoments ****************************************************
PROTO void profit_convmoments(profitstruct *profit, obj2struct *obj2)
PURPOSE Compute the 2nd order moments of the convolved object model.
INPUT Profile-fitting structure,
Pointer to obj2 structure.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 12/04/2011
***/
void profit_convmoments(profitstruct *profit, obj2struct *obj2)
{
double hw,hh, r2max, x,xstart,y, mx2,my2,mxy,mx,my,sum, dval,
temp,temp2,invtemp2, pmx2, theta;
PIXTYPE *pix;
int ix,iy, w,h;
w = profit->modnaxisn[0];
h = profit->modnaxisn[1];
hw = (double)(w/2);
hh = (double)(h/2);
r2max = hw<hh? hw*hw : hh*hh;
xstart = -hw;
y = -hh;
pix = profit->cmodpix;
mx2 = my2 = mxy = mx = my = sum = 0.0;
for (iy=h; iy--; y+=1.0)
{
x = xstart;
for (ix=w; ix--; x+=1.0)
if (y*y+x*x <= r2max)
{
dval = *(pix++);
sum += dval;
mx += dval*x;
my += dval*y;
mx2 += dval*x*x;
mxy += dval*x*y;
my2 += dval*y*y;
}
else
pix++;
}
if (sum <= 1.0/BIG)
sum = 1.0;
mx /= sum;
my /= sum;
obj2->prof_convmx2 = (mx2 = mx2/sum - mx*mx)*profit->pixstep*profit->pixstep;
obj2->prof_convmy2 = (my2 = my2/sum - my*my)*profit->pixstep*profit->pixstep;
obj2->prof_convmxy = (mxy = mxy/sum - mx*my)*profit->pixstep*profit->pixstep;
/* Handle fully correlated profiles (which cause a singularity...) */
if ((temp2=mx2*my2-mxy*mxy)<0.00694)
{
mx2 += 0.0833333;
my2 += 0.0833333;
temp2 = mx2*my2-mxy*mxy;
}
temp2 *= profit->pixstep*profit->pixstep;
if (FLAG(obj2.prof_convcxx))
{
invtemp2 = (temp2>=0.0) ? 1.0/temp2 : 0.0;
obj2->prof_convcxx = (float)(my2*invtemp2);
obj2->prof_convcyy = (float)(mx2*invtemp2);
obj2->prof_convcxy = (float)(-2*mxy*invtemp2);
}
if (1 /*FLAG(obj2.prof_conva)*/)
{
if ((fabs(temp=mx2-my2)) > 0.0)
theta = atan2(2.0 * mxy,temp) / 2.0;
else
theta = PI/4.0;
temp = sqrt(0.25*temp*temp+mxy*mxy);
pmx2 = 0.5*(mx2+my2);
obj2->prof_conva = (float)sqrt(pmx2 + temp)*profit->pixstep;
obj2->prof_convb = (float)sqrt(pmx2 - temp)*profit->pixstep;
obj2->prof_convtheta = theta/DEG;
}
return;
}
/****** profit_surface **************************************************** /****** profit_surface ****************************************************
PROTO void profit_surface(profitstruct *profit, obj2struct *obj2) PROTO void profit_surface(profitstruct *profit, obj2struct *obj2)
PURPOSE Compute surface brightnesses from the unconvolved object model. PURPOSE Compute surface brightnesses from the unconvolved object model.
...@@ -3065,16 +3242,16 @@ void profit_covarunboundtobound(profitstruct *profit, ...@@ -3065,16 +3242,16 @@ void profit_covarunboundtobound(profitstruct *profit,
/****** prof_init ************************************************************* /****** prof_init *************************************************************
PROTO profstruct prof_init(profitstruct *profit, proftypenum profcode) PROTO profstruct prof_init(profitstruct *profit, unsigned int modeltype)
PURPOSE Allocate and initialize a new profile structure. PURPOSE Allocate and initialize a new profile structure.
INPUT Pointer to the profile-fitting structure, INPUT Pointer to the profile-fitting structure,
profile type. model type.
OUTPUT A pointer to an allocated prof structure. OUTPUT A pointer to an allocated prof structure.
NOTES -. NOTES -.
AUTHOR E. Bertin (IAP) AUTHOR E. Bertin (IAP)
VERSION 08/10/2010 VERSION 22/04/2011
***/ ***/
profstruct *prof_init(profitstruct *profit, proftypenum profcode) profstruct *prof_init(profitstruct *profit, unsigned int modeltype)
{ {
profstruct *prof; profstruct *prof;
float *pix, float *pix,
...@@ -3083,16 +3260,18 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode) ...@@ -3083,16 +3260,18 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode)
d,s; d,s;
QCALLOC(prof, profstruct, 1); QCALLOC(prof, profstruct, 1);
prof->code = profcode; prof->code = modeltype;
switch(profcode) switch(modeltype)
{ {
case PROF_BACK: case MODEL_BACK:
prof->name = "background offset";
prof->naxis = 2; prof->naxis = 2;
prof->pix = NULL; prof->pix = NULL;
profit_addparam(profit, PARAM_BACK, &prof->flux); profit_addparam(profit, PARAM_BACK, &prof->flux);
prof->typscale = 1.0; prof->typscale = 1.0;
break; break;
case PROF_DIRAC: case MODEL_DIRAC:
prof->name = "point source";
prof->naxis = 2; prof->naxis = 2;
prof->naxisn[0] = PROFIT_MAXMODSIZE; prof->naxisn[0] = PROFIT_MAXMODSIZE;
prof->naxisn[1] = PROFIT_MAXMODSIZE; prof->naxisn[1] = PROFIT_MAXMODSIZE;
...@@ -3103,7 +3282,8 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode) ...@@ -3103,7 +3282,8 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode)
profit_addparam(profit, PARAM_Y, &prof->x[1]); profit_addparam(profit, PARAM_Y, &prof->x[1]);
profit_addparam(profit, PARAM_DIRAC_FLUX, &prof->flux); profit_addparam(profit, PARAM_DIRAC_FLUX, &prof->flux);
break; break;
case PROF_SERSIC: case MODEL_SERSIC:
prof->name = "Sersic spheroid";
prof->naxis = 2; prof->naxis = 2;
prof->naxisn[0] = PROFIT_MAXMODSIZE; prof->naxisn[0] = PROFIT_MAXMODSIZE;
prof->naxisn[1] = PROFIT_MAXMODSIZE; prof->naxisn[1] = PROFIT_MAXMODSIZE;
...@@ -3119,7 +3299,8 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode) ...@@ -3119,7 +3299,8 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode)
profit_addparam(profit, PARAM_SPHEROID_POSANG, &prof->posangle); profit_addparam(profit, PARAM_SPHEROID_POSANG, &prof->posangle);
profit_addparam(profit, PARAM_SPHEROID_SERSICN, &prof->extra[0]); profit_addparam(profit, PARAM_SPHEROID_SERSICN, &prof->extra[0]);
break; break;
case PROF_DEVAUCOULEURS: case MODEL_DEVAUCOULEURS:
prof->name = "de Vaucouleurs spheroid";
prof->naxis = 2; prof->naxis = 2;
prof->naxisn[0] = PROFIT_MAXMODSIZE; prof->naxisn[0] = PROFIT_MAXMODSIZE;
prof->naxisn[1] = PROFIT_MAXMODSIZE; prof->naxisn[1] = PROFIT_MAXMODSIZE;
...@@ -3134,7 +3315,8 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode) ...@@ -3134,7 +3315,8 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode)
profit_addparam(profit, PARAM_SPHEROID_ASPECT, &prof->aspect); profit_addparam(profit, PARAM_SPHEROID_ASPECT, &prof->aspect);
profit_addparam(profit, PARAM_SPHEROID_POSANG, &prof->posangle); profit_addparam(profit, PARAM_SPHEROID_POSANG, &prof->posangle);
break; break;
case PROF_EXPONENTIAL: case MODEL_EXPONENTIAL:
prof->name = "exponential disk";
prof->naxis = 2; prof->naxis = 2;
prof->naxisn[0] = PROFIT_MAXMODSIZE; prof->naxisn[0] = PROFIT_MAXMODSIZE;
prof->naxisn[1] = PROFIT_MAXMODSIZE; prof->naxisn[1] = PROFIT_MAXMODSIZE;
...@@ -3149,7 +3331,8 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode) ...@@ -3149,7 +3331,8 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode)
profit_addparam(profit, PARAM_DISK_ASPECT, &prof->aspect); profit_addparam(profit, PARAM_DISK_ASPECT, &prof->aspect);
profit_addparam(profit, PARAM_DISK_POSANG, &prof->posangle); profit_addparam(profit, PARAM_DISK_POSANG, &prof->posangle);
break; break;
case PROF_ARMS: case MODEL_ARMS:
prof->name = "spiral arms";
prof->naxis = 2; prof->naxis = 2;
prof->naxisn[0] = PROFIT_MAXMODSIZE; prof->naxisn[0] = PROFIT_MAXMODSIZE;
prof->naxisn[1] = PROFIT_MAXMODSIZE; prof->naxisn[1] = PROFIT_MAXMODSIZE;
...@@ -3171,7 +3354,8 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode) ...@@ -3171,7 +3354,8 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode)
profit_addparam(profit, PARAM_ARMS_POSANG, &prof->featposang); profit_addparam(profit, PARAM_ARMS_POSANG, &prof->featposang);
// profit_addparam(profit, PARAM_ARMS_WIDTH, &prof->featwidth); // profit_addparam(profit, PARAM_ARMS_WIDTH, &prof->featwidth);
break; break;
case PROF_BAR: case MODEL_BAR:
prof->name = "bar";
prof->naxis = 2; prof->naxis = 2;
prof->naxisn[0] = PROFIT_MAXMODSIZE; prof->naxisn[0] = PROFIT_MAXMODSIZE;
prof->naxisn[1] = PROFIT_MAXMODSIZE; prof->naxisn[1] = PROFIT_MAXMODSIZE;
...@@ -3189,7 +3373,8 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode) ...@@ -3189,7 +3373,8 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode)
profit_addparam(profit, PARAM_BAR_ASPECT, &prof->feataspect); profit_addparam(profit, PARAM_BAR_ASPECT, &prof->feataspect);
profit_addparam(profit, PARAM_ARMS_POSANG, &prof->featposang); profit_addparam(profit, PARAM_ARMS_POSANG, &prof->featposang);
break; break;
case PROF_INRING: case MODEL_INRING:
prof->name = "inner ring";
prof->naxis = 2; prof->naxis = 2;
prof->naxisn[0] = PROFIT_MAXMODSIZE; prof->naxisn[0] = PROFIT_MAXMODSIZE;
prof->naxisn[1] = PROFIT_MAXMODSIZE; prof->naxisn[1] = PROFIT_MAXMODSIZE;
...@@ -3207,7 +3392,8 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode) ...@@ -3207,7 +3392,8 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode)
profit_addparam(profit, PARAM_INRING_WIDTH, &prof->featwidth); profit_addparam(profit, PARAM_INRING_WIDTH, &prof->featwidth);
profit_addparam(profit, PARAM_INRING_ASPECT, &prof->feataspect); profit_addparam(profit, PARAM_INRING_ASPECT, &prof->feataspect);
break; break;
case PROF_OUTRING: case MODEL_OUTRING:
prof->name = "outer ring";
prof->naxis = 2; prof->naxis = 2;
prof->naxisn[0] = PROFIT_MAXMODSIZE; prof->naxisn[0] = PROFIT_MAXMODSIZE;
prof->naxisn[1] = PROFIT_MAXMODSIZE; prof->naxisn[1] = PROFIT_MAXMODSIZE;
...@@ -3224,7 +3410,8 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode) ...@@ -3224,7 +3410,8 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode)
profit_addparam(profit, PARAM_OUTRING_FLUX, &prof->flux); profit_addparam(profit, PARAM_OUTRING_FLUX, &prof->flux);
profit_addparam(profit, PARAM_OUTRING_WIDTH, &prof->featwidth); profit_addparam(profit, PARAM_OUTRING_WIDTH, &prof->featwidth);
break; break;
case PROF_SERSIC_TABEX: /* An example of tabulated profile */ case MODEL_TABULATED: /* An example of tabulated profile */
prof->name = "tabulated model";
prof->naxis = 3; prof->naxis = 3;
width = prof->naxisn[0] = PROFIT_PROFRES; width = prof->naxisn[0] = PROFIT_PROFRES;
height = prof->naxisn[1] = PROFIT_PROFRES; height = prof->naxisn[1] = PROFIT_PROFRES;
...@@ -3345,7 +3532,7 @@ float prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag) ...@@ -3345,7 +3532,7 @@ float prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag)
npix = profit->nmodpix; npix = profit->nmodpix;
if (prof->code==PROF_BACK) if (prof->code==MODEL_BACK)
{ {
amp = fabs(*prof->flux); amp = fabs(*prof->flux);
pixout = profit->modpix; pixout = profit->modpix;
...@@ -3357,7 +3544,7 @@ float prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag) ...@@ -3357,7 +3544,7 @@ float prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag)
scaling = profit->pixstep / prof->typscale; scaling = profit->pixstep / prof->typscale;
if (prof->code!=PROF_DIRAC) if (prof->code!=MODEL_DIRAC)
{ {
/*-- Compute Profile CD matrix */ /*-- Compute Profile CD matrix */
ctheta = cos(*prof->posangle*DEG); ctheta = cos(*prof->posangle*DEG);
...@@ -3395,24 +3582,24 @@ float prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag) ...@@ -3395,24 +3582,24 @@ float prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag)
switch(prof->code) switch(prof->code)
{ {
case PROF_DIRAC: case MODEL_DIRAC:
prof->pix[profit->modnaxisn[0]/2 prof->pix[profit->modnaxisn[0]/2
+ (profit->modnaxisn[1]/2)*profit->modnaxisn[0]] = 1.0; + (profit->modnaxisn[1]/2)*profit->modnaxisn[0]] = 1.0;
prof->lostfluxfrac = 0.0; prof->lostfluxfrac = 0.0;
threshflag = 0; threshflag = 0;
break; break;
case PROF_SERSIC: case MODEL_SERSIC:
case PROF_DEVAUCOULEURS: case MODEL_DEVAUCOULEURS:
case PROF_EXPONENTIAL: case MODEL_EXPONENTIAL:
/*---- Compute sharp/smooth transition radius */ /*---- Compute sharp/smooth transition radius */
rs = PROFIT_SMOOTHR*(xscale>yscale?xscale:yscale); rs = PROFIT_SMOOTHR*(xscale>yscale?xscale:yscale);
if (rs<=0) if (rs<=0)
rs = 1.0; rs = 1.0;
rs2 = rs*rs; rs2 = rs*rs;
/*---- The consequence of sampling on flux is compensated by PSF normalisation*/ /*---- The consequence of sampling on flux is compensated by PSF normalisation*/
if (prof->code==PROF_EXPONENTIAL) if (prof->code==MODEL_EXPONENTIAL)
bn = n = 1.0; bn = n = 1.0;
else if (prof->code==PROF_DEVAUCOULEURS) else if (prof->code==MODEL_DEVAUCOULEURS)
{ {
n = 4.0; n = 4.0;
bn = 7.66924944; bn = 7.66924944;
...@@ -3427,7 +3614,7 @@ float prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag) ...@@ -3427,7 +3614,7 @@ float prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag)
hinvn = 0.5/n; hinvn = 0.5/n;
k = -bn; k = -bn;
/*---- Compute central polynomial terms */ /*---- Compute central polynomial terms */
krspinvn = prof->code==PROF_EXPONENTIAL? -rs : k*expf(logf(rs)*invn); krspinvn = prof->code==MODEL_EXPONENTIAL? -rs : k*expf(logf(rs)*invn);
ekrspinvn = expf(krspinvn); ekrspinvn = expf(krspinvn);
p2 = krspinvn*invn*invn; p2 = krspinvn*invn*invn;
p1 = krspinvn*p2; p1 = krspinvn*p2;
...@@ -3438,7 +3625,7 @@ float prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag) ...@@ -3438,7 +3625,7 @@ float prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag)
x10 = -x1cout - dx1; x10 = -x1cout - dx1;
x2 = -x2cout - dx2; x2 = -x2cout - dx2;
pixin = prof->pix; pixin = prof->pix;
if (prof->code==PROF_EXPONENTIAL) if (prof->code==MODEL_EXPONENTIAL)
for (ix2=nx2; ix2--; x2+=1.0) for (ix2=nx2; ix2--; x2+=1.0)
{ {
x1 = x10; x1 = x10;
...@@ -3531,12 +3718,12 @@ float prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag) ...@@ -3531,12 +3718,12 @@ float prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag)
krpinvn *= dkrpinvn; krpinvn *= dkrpinvn;
} }
prof->lostfluxfrac = prof->code==PROF_EXPONENTIAL? prof->lostfluxfrac = prof->code==MODEL_EXPONENTIAL?
(1.0 + rmax)*exp(-rmax) (1.0 + rmax)*exp(-rmax)
:1.0 - prof_gammainc(2.0*n, bn*pow(r2max, hinvn)); :1.0 - prof_gammainc(2.0*n, bn*pow(r2max, hinvn));
threshflag = 0; threshflag = 0;
break; break;
case PROF_ARMS: case MODEL_ARMS:
r2min = *prof->featstart**prof->featstart; r2min = *prof->featstart**prof->featstart;
r2minxin = r2min * (1.0 - PROFIT_BARXFADE) * (1.0 - PROFIT_BARXFADE); r2minxin = r2min * (1.0 - PROFIT_BARXFADE) * (1.0 - PROFIT_BARXFADE);
r2minxout = r2min * (1.0 + PROFIT_BARXFADE) * (1.0 + PROFIT_BARXFADE); r2minxout = r2min * (1.0 + PROFIT_BARXFADE) * (1.0 + PROFIT_BARXFADE);
...@@ -3587,7 +3774,7 @@ width = 3.0; ...@@ -3587,7 +3774,7 @@ width = 3.0;
prof->lostfluxfrac = 0.0; prof->lostfluxfrac = 0.0;
threshflag = 1; threshflag = 1;
break; break;
case PROF_BAR: case MODEL_BAR:
r2min = *prof->featstart**prof->featstart; r2min = *prof->featstart**prof->featstart;
r2minxin = r2min * (1.0 - PROFIT_BARXFADE) * (1.0 - PROFIT_BARXFADE); r2minxin = r2min * (1.0 - PROFIT_BARXFADE) * (1.0 - PROFIT_BARXFADE);
r2minxout = r2min * (1.0 + PROFIT_BARXFADE) * (1.0 + PROFIT_BARXFADE); r2minxout = r2min * (1.0 + PROFIT_BARXFADE) * (1.0 + PROFIT_BARXFADE);
...@@ -3626,7 +3813,7 @@ width = 3.0; ...@@ -3626,7 +3813,7 @@ width = 3.0;
prof->lostfluxfrac = 0.0; prof->lostfluxfrac = 0.0;
threshflag = 1; threshflag = 1;
break; break;
case PROF_INRING: case MODEL_INRING:
rmin = *prof->featstart; rmin = *prof->featstart;
r2minxin = *prof->featstart-4.0**prof->featwidth; r2minxin = *prof->featstart-4.0**prof->featwidth;
if (r2minxin < 0.0) if (r2minxin < 0.0)
...@@ -3661,7 +3848,7 @@ width = 3.0; ...@@ -3661,7 +3848,7 @@ width = 3.0;
prof->lostfluxfrac = 0.0; prof->lostfluxfrac = 0.0;
threshflag = 1; threshflag = 1;
break; break;
case PROF_OUTRING: case MODEL_OUTRING:
rmin = *prof->featstart; rmin = *prof->featstart;
r2minxin = *prof->featstart-4.0**prof->featwidth; r2minxin = *prof->featstart-4.0**prof->featwidth;
if (r2minxin < 0.0) if (r2minxin < 0.0)
...@@ -3866,7 +4053,7 @@ int prof_moments(profitstruct *profit, profstruct *prof, double *jac) ...@@ -3866,7 +4053,7 @@ int prof_moments(profitstruct *profit, profstruct *prof, double *jac)
dc2 = ds2 = dcs = 0.0; /* To avoid gcc -Wall warnings */ dc2 = ds2 = dcs = 0.0; /* To avoid gcc -Wall warnings */
switch(prof->code) switch(prof->code)
{ {
case PROF_SERSIC: case MODEL_SERSIC:
n = fabs(*prof->extra[0]); n = fabs(*prof->extra[0]);
bn = 2.0*n - 1.0/3.0 + 4.0/(405.0*n) + 46.0/(25515.0*n*n) bn = 2.0*n - 1.0/3.0 + 4.0/(405.0*n) + 46.0/(25515.0*n*n)
+ 131.0/(1148175*n*n*n); /* Ciotti & Bertin 1999 */ + 131.0/(1148175*n*n*n); /* Ciotti & Bertin 1999 */
...@@ -3901,7 +4088,7 @@ int prof_moments(profitstruct *profit, profstruct *prof, double *jac) ...@@ -3901,7 +4088,7 @@ int prof_moments(profitstruct *profit, profstruct *prof, double *jac)
} }
index = profit->paramindex[PARAM_SPHEROID_FLUX]; index = profit->paramindex[PARAM_SPHEROID_FLUX];
break; break;
case PROF_DEVAUCOULEURS: case MODEL_DEVAUCOULEURS:
m20 = 10.83995 * *prof->scale**prof->scale; m20 = 10.83995 * *prof->scale**prof->scale;
if (jac) if (jac)
{ {
...@@ -3923,7 +4110,7 @@ int prof_moments(profitstruct *profit, profstruct *prof, double *jac) ...@@ -3923,7 +4110,7 @@ int prof_moments(profitstruct *profit, profstruct *prof, double *jac)
} }
index = profit->paramindex[PARAM_SPHEROID_FLUX]; index = profit->paramindex[PARAM_SPHEROID_FLUX];
break; break;
case PROF_EXPONENTIAL: case MODEL_EXPONENTIAL:
m20 = 3.0 * *prof->scale**prof->scale; m20 = 3.0 * *prof->scale**prof->scale;
if (jac) if (jac)
{ {
...@@ -3945,19 +4132,19 @@ int prof_moments(profitstruct *profit, profstruct *prof, double *jac) ...@@ -3945,19 +4132,19 @@ int prof_moments(profitstruct *profit, profstruct *prof, double *jac)
} }
index = profit->paramindex[PARAM_DISK_FLUX]; index = profit->paramindex[PARAM_DISK_FLUX];
break; break;
case PROF_ARMS: case MODEL_ARMS:
m20 = 1.0; m20 = 1.0;
index = profit->paramindex[PARAM_ARMS_FLUX]; index = profit->paramindex[PARAM_ARMS_FLUX];
break; break;
case PROF_BAR: case MODEL_BAR:
m20 = 1.0; m20 = 1.0;
index = profit->paramindex[PARAM_BAR_FLUX]; index = profit->paramindex[PARAM_BAR_FLUX];
break; break;
case PROF_INRING: case MODEL_INRING:
m20 = 1.0; m20 = 1.0;
index = profit->paramindex[PARAM_INRING_FLUX]; index = profit->paramindex[PARAM_INRING_FLUX];
break; break;
case PROF_OUTRING: case MODEL_OUTRING:
m20 = 1.0; m20 = 1.0;
index = profit->paramindex[PARAM_OUTRING_FLUX]; index = profit->paramindex[PARAM_OUTRING_FLUX];
break; break;
......
...@@ -7,7 +7,7 @@ ...@@ -7,7 +7,7 @@
* *
* This file part of: SExtractor * This file part of: SExtractor
* *
* Copyright: (C) 2006-2010 Emmanuel Bertin -- IAP/CNRS/UPMC * Copyright: (C) 2006-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
* *
* License: GNU General Public License * License: GNU General Public License
* *
...@@ -22,13 +22,28 @@ ...@@ -22,13 +22,28 @@
* 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: 24/01/2011 * Last modified: 22/04/2011
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
#ifndef _PROFIT_H_ #ifndef _PROFIT_H_
#define _PROFIT_H_ #define _PROFIT_H_
/*-------------------------------- models -----------------------------------*/
#define MODEL_NONE 0x0000
#define MODEL_BACK 0x0001
#define MODEL_DIRAC 0x0002
#define MODEL_SERSIC 0x0004
#define MODEL_DEVAUCOULEURS 0x0008
#define MODEL_EXPONENTIAL 0x0010
#define MODEL_ARMS 0x0020
#define MODEL_BAR 0x0040
#define MODEL_INRING 0x0080
#define MODEL_OUTRING 0x0100
#define MODEL_TABULATED 0x0200
#define MODEL_NMAX 11
/*-------------------------------- flags ------------------------------------*/ /*-------------------------------- flags ------------------------------------*/
#define PROFLAG_MODSUB 0x0001 #define PROFLAG_MODSUB 0x0001
...@@ -64,11 +79,6 @@ One must have: PROFIT_NITER > 0 ...@@ -64,11 +79,6 @@ One must have: PROFIT_NITER > 0
/*--------------------------------- typedefs --------------------------------*/ /*--------------------------------- typedefs --------------------------------*/
typedef enum {PROF_BACK, PROF_DIRAC, PROF_SERSIC, PROF_DEVAUCOULEURS,
PROF_EXPONENTIAL, PROF_ARMS, PROF_BAR, PROF_INRING,
PROF_OUTRING, PROF_SERSIC_TABEX, PROF_NPROF}
proftypenum; /* Profile code */
typedef enum {INTERP_NEARESTNEIGHBOUR, INTERP_BILINEAR, INTERP_LANCZOS2, typedef enum {INTERP_NEARESTNEIGHBOUR, INTERP_BILINEAR, INTERP_LANCZOS2,
INTERP_LANCZOS3, INTERP_LANCZOS4} interpenum; INTERP_LANCZOS3, INTERP_LANCZOS4} interpenum;
...@@ -90,7 +100,8 @@ typedef enum {PARAM_BACK, ...@@ -90,7 +100,8 @@ typedef enum {PARAM_BACK,
typedef struct typedef struct
{ {
proftypenum code; /* Model code */ unsigned int code; /* Model code */
char *name; /* Model name */
float *pix; /* Full pixmap of the model */ float *pix; /* Full pixmap of the model */
int naxis; /* Number of pixmap dimensions */ int naxis; /* Number of pixmap dimensions */
int naxisn[3]; /* Pixmap size for each axis */ int naxisn[3]; /* Pixmap size for each axis */
...@@ -174,9 +185,9 @@ typedef struct ...@@ -174,9 +185,9 @@ typedef struct
/*----------------------------- Global variables ----------------------------*/ /*----------------------------- Global variables ----------------------------*/
/*-------------------------------- functions --------------------------------*/ /*-------------------------------- functions --------------------------------*/
profitstruct *profit_init(struct psf *psf); profitstruct *profit_init(struct psf *psf, unsigned int modeltype);
profstruct *prof_init(profitstruct *profit, proftypenum profcode); profstruct *prof_init(profitstruct *profit, unsigned int modeltype);
float *profit_compresi(profitstruct *profit, float dynparam, float *profit_compresi(profitstruct *profit, float dynparam,
float *resi), float *resi),
...@@ -186,6 +197,7 @@ float *profit_compresi(profitstruct *profit, float dynparam, ...@@ -186,6 +197,7 @@ float *profit_compresi(profitstruct *profit, float dynparam,
prof_add(profitstruct *profit, profstruct *prof, prof_add(profitstruct *profit, profstruct *prof,
int extfluxfac_flag), int extfluxfac_flag),
profit_minradius(profitstruct *profit, float refffac), profit_minradius(profitstruct *profit, float refffac),
profit_noisearea(profitstruct *profit),
profit_spiralindex(profitstruct *profit); profit_spiralindex(profitstruct *profit);
int profit_copyobjpix(profitstruct *profit, picstruct *field, int profit_copyobjpix(profitstruct *profit, picstruct *field,
...@@ -206,6 +218,7 @@ void prof_end(profstruct *prof), ...@@ -206,6 +218,7 @@ void prof_end(profstruct *prof),
profit_fit(profitstruct *profit, profit_fit(profitstruct *profit,
picstruct *field, picstruct *wfield, picstruct *field, picstruct *wfield,
objstruct *obj, obj2struct *obj2), objstruct *obj, obj2struct *obj2),
profit_convmoments(profitstruct *profit, obj2struct *obj2),
profit_convolve(profitstruct *profit, float *modpix), profit_convolve(profitstruct *profit, float *modpix),
profit_covarunboundtobound(profitstruct *profit, profit_covarunboundtobound(profitstruct *profit,
double *dparam, float *param), double *dparam, float *param),
......
...@@ -7,7 +7,7 @@ ...@@ -7,7 +7,7 @@
* *
* This file part of: SExtractor * 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 * License: GNU General Public License
* *
...@@ -22,7 +22,7 @@ ...@@ -22,7 +22,7 @@
* 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: 11/10/2010 * Last modified: 25/03/2010
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
...@@ -165,6 +165,12 @@ void *loadstrip(picstruct *field, picstruct *wfield) ...@@ -165,6 +165,12 @@ void *loadstrip(picstruct *field, picstruct *wfield)
if ((flags & MEASURE_FIELD) && (check=prefs.check[CHECK_SUBOBJECTS])) if ((flags & MEASURE_FIELD) && (check=prefs.check[CHECK_SUBOBJECTS]))
writecheck(check, data, w); writecheck(check, data, w);
if ((flags & MEASURE_FIELD) && (check=prefs.check[CHECK_MASK]))
writecheck(check, data, w);
if ((flags & MEASURE_FIELD) && (check=prefs.check[CHECK_SUBMASK]))
writecheck(check, data, w);
if (flags & BACKRMS_FIELD) if (flags & BACKRMS_FIELD)
backrmsline(field, field->ymax, data); backrmsline(field, field->ymax, data);
else if (flags & INTERP_FIELD) else if (flags & INTERP_FIELD)
......
...@@ -7,7 +7,7 @@ ...@@ -7,7 +7,7 @@
* *
* This file part of: SExtractor * 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 * License: GNU General Public License
* *
...@@ -22,7 +22,7 @@ ...@@ -22,7 +22,7 @@
* 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: 12/01/2011 * Last modified: 03/05/2011
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
...@@ -71,8 +71,8 @@ typedef enum {BACK_RELATIVE, BACK_ABSOLUTE} ...@@ -71,8 +71,8 @@ typedef enum {BACK_RELATIVE, BACK_ABSOLUTE}
typedef enum {CHECK_NONE, CHECK_IDENTICAL, CHECK_BACKGROUND, typedef enum {CHECK_NONE, CHECK_IDENTICAL, CHECK_BACKGROUND,
CHECK_BACKRMS, CHECK_MINIBACKGROUND, CHECK_MINIBACKRMS, CHECK_BACKRMS, CHECK_MINIBACKGROUND, CHECK_MINIBACKRMS,
CHECK_SUBTRACTED, CHECK_FILTERED, CHECK_OBJECTS, CHECK_APERTURES, CHECK_SUBTRACTED, CHECK_FILTERED, CHECK_OBJECTS, CHECK_APERTURES,
CHECK_SEGMENTATION, CHECK_ASSOC, CHECK_SUBOBJECTS, CHECK_SEGMENTATION, CHECK_MASK, CHECK_SUBMASK, CHECK_ASSOC,
CHECK_PSFPROTOS, CHECK_SUBPSFPROTOS, CHECK_SUBOBJECTS, CHECK_PSFPROTOS, CHECK_SUBPSFPROTOS,
CHECK_PCPROTOS, CHECK_SUBPCPROTOS, CHECK_PCOPROTOS, CHECK_PCPROTOS, CHECK_SUBPCPROTOS, CHECK_PCOPROTOS,
CHECK_MAPSOM, CHECK_PROFILES, CHECK_SUBPROFILES, CHECK_MAPSOM, CHECK_PROFILES, CHECK_SUBPROFILES,
CHECK_SPHEROIDS, CHECK_SUBSPHEROIDS, CHECK_DISKS, CHECK_SUBDISKS, CHECK_SPHEROIDS, CHECK_SUBSPHEROIDS, CHECK_DISKS, CHECK_SUBDISKS,
...@@ -165,6 +165,7 @@ typedef struct ...@@ -165,6 +165,7 @@ typedef struct
float *magerr_aper; /* APER mag error vector */ float *magerr_aper; /* APER mag error vector */
float flux_win; /* WINdowed flux*/ float flux_win; /* WINdowed flux*/
float fluxerr_win; /* WINdowed flux error */ float fluxerr_win; /* WINdowed flux error */
float snr_win; /* WINdowed SNR */
float mag_win; /* WINdowed magnitude */ float mag_win; /* WINdowed magnitude */
float magerr_win; /* WINdowed magnitude error */ float magerr_win; /* WINdowed magnitude error */
/* ---- astrometric data */ /* ---- astrometric data */
...@@ -371,6 +372,12 @@ typedef struct ...@@ -371,6 +372,12 @@ typedef struct
prof_theta; /* Model-fitting ellip. params*/ prof_theta; /* Model-fitting ellip. params*/
float prof_cxx, prof_cyy, float prof_cxx, prof_cyy,
prof_cxy; /* Model-fitting ellip. params*/ prof_cxy; /* Model-fitting ellip. params*/
double prof_convmx2, prof_convmy2,
prof_convmxy; /* Convolved model moments */
float prof_convcxx, prof_convcyy,
prof_convcxy; /* Conv. model ellipse params */
float prof_conva, prof_convb,
prof_convtheta; /* Conv. model ellipse params */
float prof_pol1, prof_pol2; /* Model-fitting pol. vector*/ float prof_pol1, prof_pol2; /* Model-fitting pol. vector*/
float prof_pol1err, prof_pol2err, float prof_pol1err, prof_pol2err,
prof_pol12corr; /* Model-fitting pol. errors */ prof_pol12corr; /* Model-fitting pol. errors */
...@@ -395,6 +402,7 @@ typedef struct ...@@ -395,6 +402,7 @@ typedef struct
float prof_class_star; /* Mod.-fitting star/gal class*/ float prof_class_star; /* Mod.-fitting star/gal class*/
float prof_concentration; /* Model-fitting concentration*/ float prof_concentration; /* Model-fitting concentration*/
float prof_concentrationerr; /* RMS error */ float prof_concentrationerr; /* RMS error */
float prof_noisearea; /* Equivalent noise area */
float prof_offset_flux; /* Background offset */ float prof_offset_flux; /* Background offset */
float prof_offset_fluxerr; /* RMS error */ float prof_offset_fluxerr; /* RMS error */
float prof_dirac_flux; /* Point source total flux */ float prof_dirac_flux; /* Point source total flux */
......
...@@ -7,7 +7,7 @@ ...@@ -7,7 +7,7 @@
* *
* This file part of: SExtractor * This file part of: SExtractor
* *
* Copyright: (C) 2005-2010 Emmanuel Bertin -- IAP/CNRS/UPMC * Copyright: (C) 2005-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
* *
* License: GNU General Public License * License: GNU General Public License
* *
...@@ -22,7 +22,7 @@ ...@@ -22,7 +22,7 @@
* 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: 11/10/2010 * Last modified: 03/05/2011
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
...@@ -50,7 +50,7 @@ INPUT Picture structure pointer, ...@@ -50,7 +50,7 @@ INPUT Picture structure pointer,
OUTPUT -. OUTPUT -.
NOTES obj->posx and obj->posy are taken as initial centroid guesses. NOTES obj->posx and obj->posy are taken as initial centroid guesses.
AUTHOR E. Bertin (IAP) AUTHOR E. Bertin (IAP)
VERSION 16/07/2010 VERSION 03/05/2011
***/ ***/
void compute_winpos(picstruct *field, picstruct *wfield, objstruct *obj) void compute_winpos(picstruct *field, picstruct *wfield, objstruct *obj)
...@@ -79,7 +79,7 @@ void compute_winpos(picstruct *field, picstruct *wfield, objstruct *obj) ...@@ -79,7 +79,7 @@ void compute_winpos(picstruct *field, picstruct *wfield, objstruct *obj)
pflag = (prefs.detect_type==PHOTO)? 1:0; pflag = (prefs.detect_type==PHOTO)? 1:0;
corrflag = (prefs.mask_type==MASK_CORRECT); corrflag = (prefs.mask_type==MASK_CORRECT);
gainflag = wfield && prefs.weightgain_flag; gainflag = wfield && prefs.weightgain_flag;
errflag = FLAG(obj2.winposerr_mx2); errflag = FLAG(obj2.winposerr_mx2) | FLAG(obj2.fluxerr_win);
momentflag = FLAG(obj2.win_mx2) | FLAG(obj2.winposerr_mx2); momentflag = FLAG(obj2.win_mx2) | FLAG(obj2.winposerr_mx2);
var = backnoise2 = field->backsig*field->backsig; var = backnoise2 = field->backsig*field->backsig;
invgain = field->gain>0.0? 1.0/field->gain : 0.0; invgain = field->gain>0.0? 1.0/field->gain : 0.0;
...@@ -262,13 +262,14 @@ void compute_winpos(picstruct *field, picstruct *wfield, objstruct *obj) ...@@ -262,13 +262,14 @@ void compute_winpos(picstruct *field, picstruct *wfield, objstruct *obj)
{ {
obj2->flux_win = tv; obj2->flux_win = tv;
obj2->fluxerr_win = sqrt(esum); obj2->fluxerr_win = sqrt(esum);
obj2->snr_win = esum>(1.0/BIG)? obj2->flux_win / obj2->fluxerr_win: BIG;
} }
temp2=mx2*my2-mxy*mxy; temp2=mx2*my2-mxy*mxy;
obj2->win_flag = (tv <= 0.0)*4 + (mx2 < 0.0 || my2 < 0.0)*2 + (temp2<0.0); obj2->win_flag = (tv <= 0.0)*4 + (mx2 < 0.0 || my2 < 0.0)*2 + (temp2<0.0);
if (obj2->win_flag) if (obj2->win_flag)
{ {
/*--- Negative values: revert to isophotal estimates */ /*--- Negative values: revert to isophotal estimates */
if (errflag) if (FLAG(obj2.winposerr_mx2))
{ {
obj2->winposerr_mx2 = obj->poserr_mx2; obj2->winposerr_mx2 = obj->poserr_mx2;
obj2->winposerr_my2 = obj->poserr_my2; obj2->winposerr_my2 = obj->poserr_my2;
...@@ -308,7 +309,7 @@ void compute_winpos(picstruct *field, picstruct *wfield, objstruct *obj) ...@@ -308,7 +309,7 @@ void compute_winpos(picstruct *field, picstruct *wfield, objstruct *obj)
} }
else else
{ {
if (errflag) if (FLAG(obj2.winposerr_mx2))
{ {
norm = WINPOS_FAC*WINPOS_FAC/(tv*tv); norm = WINPOS_FAC*WINPOS_FAC/(tv*tv);
emx2 *= norm; emx2 *= norm;
......
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