Commit f8fe9609 authored by Emmanuel Bertin's avatar Emmanuel Bertin
Browse files

Added support for Gaussian priors on (transformed) model parameters (deactivated).

Fixed model sub-sampling issues with large objects, including check-image resizing.
Fixed initialization of DETMODEL size guess (thanks to Eli Rykoff).
Fixed image extension pb in double weighting/single image mode.
Improved debug info for memory allocations in FFT and FITS libraries.
Fixed crashes with DETMODEL magnitudes on null measurement weight maps.
Added Eli Rykoff's fix to discard sources with NaNs in the fitted positions.
Improved MKL autoconfiguration.
Improved vectorization in FITS data transformations and FFT-based convolutions.
Pushed SExtractor version number to 3.18.8.
parent b4bff35a
......@@ -7,8 +7,8 @@
#
# This file part of: AstrOmatic software
#
# Copyright: (C) 2007-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
## (C) 2004 Manolis Lourakis (original LevMar)
# Copyright: (C) 2007-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
## (C) 2004-2011 Manolis Lourakis (orig. LevMar)
#
# Licenses: GNU General Public License
#
......@@ -24,7 +24,7 @@
# along with AstrOmatic software.
# If not, see <http://www.gnu.org/licenses/>.
#
# Last modified: 10/10/2010
# Last modified: 10/07/2012
#
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
......@@ -35,5 +35,6 @@ EXTRA_liblevmar_a_SOURCES = Axb_core.c lmbc_core.c lm_core.c \
lmblec_core.c lmbleic_core.c lmlec_core.c \
misc_core.c \
LICENSE README README.txt \
Makefile.icc Makefile.vc expfit.c lmdemo.c lm.h matlab
Makefile.icc Makefile.vc expfit.c levmar.h.in \
lmdemo.c lm.h matlab
# Makefile.in generated by automake 1.11.1 from Makefile.am.
# Makefile.in generated by automake 1.12.2 from Makefile.am.
# @configure_input@
# Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002,
# 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation,
# Inc.
# Copyright (C) 1994-2012 Free Software Foundation, Inc.
# This Makefile.in is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
# with or without modifications, as long as this notice is preserved.
......@@ -24,7 +23,7 @@
#
# This file part of: AstrOmatic software
#
# Copyright: (C) 2007-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
# Copyright: (C) 2007-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
#
# Licenses: GNU General Public License
#
......@@ -40,11 +39,28 @@
# along with AstrOmatic software.
# If not, see <http://www.gnu.org/licenses/>.
#
# Last modified: 10/10/2010
# Last modified: 10/07/2012
#
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
VPATH = @srcdir@
am__make_dryrun = \
{ \
am__dry=no; \
case $$MAKEFLAGS in \
*\\[\ \ ]*) \
echo 'am--echo: ; @echo "AM" OK' | $(MAKE) -f - 2>/dev/null \
| grep '^AM OK$$' >/dev/null || am__dry=yes;; \
*) \
for am__flg in $$MAKEFLAGS; do \
case $$am__flg in \
*=*|--*) ;; \
*n*) am__dry=yes; break;; \
esac; \
done;; \
esac; \
test $$am__dry = yes; \
}
pkgdatadir = $(datadir)/@PACKAGE@
pkgincludedir = $(includedir)/@PACKAGE@
pkglibdir = $(libdir)/@PACKAGE@
......@@ -64,7 +80,9 @@ POST_UNINSTALL = :
build_triplet = @build@
host_triplet = @host@
subdir = src/levmar
DIST_COMMON = README $(srcdir)/Makefile.am $(srcdir)/Makefile.in
DIST_COMMON = README $(srcdir)/Makefile.am $(srcdir)/Makefile.in \
$(top_srcdir)/autoconf/depcomp \
$(top_srcdir)/autoconf/mkinstalldirs
ACLOCAL_M4 = $(top_srcdir)/aclocal.m4
am__aclocal_m4_deps = $(top_srcdir)/acx_atlas.m4 \
$(top_srcdir)/acx_fftw.m4 $(top_srcdir)/acx_mkl.m4 \
......@@ -101,6 +119,11 @@ LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) \
$(LDFLAGS) -o $@
SOURCES = $(liblevmar_a_SOURCES) $(EXTRA_liblevmar_a_SOURCES)
DIST_SOURCES = $(liblevmar_a_SOURCES) $(EXTRA_liblevmar_a_SOURCES)
am__can_run_installinfo = \
case $$AM_UPDATE_INFO_DIR in \
n|no|NO) false;; \
*) (install-info --version) >/dev/null 2>&1;; \
esac
ETAGS = etags
CTAGS = ctags
DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST)
......@@ -159,6 +182,7 @@ MAKEINFO = @MAKEINFO@
MANIFEST_TOOL = @MANIFEST_TOOL@
MKDIR_P = @MKDIR_P@
MKL_CFLAGS = @MKL_CFLAGS@
MKL_LDFLAGS = @MKL_LDFLAGS@
MKL_LIBS = @MKL_LIBS@
MKL_WARN = @MKL_WARN@
NM = @NM@
......@@ -245,7 +269,8 @@ EXTRA_liblevmar_a_SOURCES = Axb_core.c lmbc_core.c lm_core.c \
lmblec_core.c lmbleic_core.c lmlec_core.c \
misc_core.c \
LICENSE README README.txt \
Makefile.icc Makefile.vc expfit.c lmdemo.c lm.h matlab
Makefile.icc Makefile.vc expfit.c levmar.h.in \
lmdemo.c lm.h matlab
all: all-am
......@@ -284,7 +309,7 @@ $(am__aclocal_m4_deps):
clean-noinstLIBRARIES:
-test -z "$(noinst_LIBRARIES)" || rm -f $(noinst_LIBRARIES)
liblevmar.a: $(liblevmar_a_OBJECTS) $(liblevmar_a_DEPENDENCIES)
liblevmar.a: $(liblevmar_a_OBJECTS) $(liblevmar_a_DEPENDENCIES) $(EXTRA_liblevmar_a_DEPENDENCIES)
-rm -f liblevmar.a
$(liblevmar_a_AR) liblevmar.a $(liblevmar_a_OBJECTS) $(liblevmar_a_LIBADD)
$(RANLIB) liblevmar.a
......@@ -388,6 +413,20 @@ GTAGS:
&& $(am__cd) $(top_srcdir) \
&& gtags -i $(GTAGS_ARGS) "$$here"
cscopelist: $(HEADERS) $(SOURCES) $(LISP)
list='$(SOURCES) $(HEADERS) $(LISP)'; \
case "$(srcdir)" in \
[\\/]* | ?:[\\/]*) sdir="$(srcdir)" ;; \
*) sdir=$(subdir)/$(srcdir) ;; \
esac; \
for i in $$list; do \
if test -f "$$i"; then \
echo "$(subdir)/$$i"; \
else \
echo "$$sdir/$$i"; \
fi; \
done >> $(top_builddir)/cscope.files
distclean-tags:
-rm -f TAGS ID GTAGS GRTAGS GSYMS GPATH tags
......@@ -435,10 +474,15 @@ install-am: all-am
installcheck: installcheck-am
install-strip:
if test -z '$(STRIP)'; then \
$(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM="$(INSTALL_STRIP_PROGRAM)" \
install_sh_PROGRAM="$(INSTALL_STRIP_PROGRAM)" INSTALL_STRIP_FLAG=-s \
install; \
else \
$(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM="$(INSTALL_STRIP_PROGRAM)" \
install_sh_PROGRAM="$(INSTALL_STRIP_PROGRAM)" INSTALL_STRIP_FLAG=-s \
`test -z '$(STRIP)' || \
echo "INSTALL_PROGRAM_ENV=STRIPPROG='$(STRIP)'"` install
"INSTALL_PROGRAM_ENV=STRIPPROG='$(STRIP)'" install; \
fi
mostlyclean-generic:
clean-generic:
......@@ -524,7 +568,7 @@ uninstall-am:
.MAKE: install-am install-strip
.PHONY: CTAGS GTAGS all all-am check check-am clean clean-generic \
clean-libtool clean-noinstLIBRARIES ctags distclean \
clean-libtool clean-noinstLIBRARIES cscopelist ctags distclean \
distclean-compile distclean-generic distclean-libtool \
distclean-tags distdir dvi dvi-am html html-am info info-am \
install install-am install-data install-data-am install-dvi \
......
......@@ -7,7 +7,7 @@
*
* This file part of: SExtractor
*
* Copyright: (C) 1993-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 1993-2013 Emmanuel Bertin -- IAP/CNRS/UPMC
*
* License: GNU General Public License
*
......@@ -22,7 +22,7 @@
* You should have received a copy of the GNU General Public License
* along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 12/07/2012
* Last modified: 13/02/2013
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -337,7 +337,7 @@ void makeit()
field = newfield(prefs.image_name[0], DETECT_FIELD | MEASURE_FIELD,
nima0<0? ntab:nima0);
/*-- Prepare interpolation */
/*---- Prepare interpolation */
if ((prefs.dweight_flag || prefs.weight_flag)
&& prefs.interp_type[0] == INTERP_ALL)
init_interpolate(field, -1, -1); /* 0.0 or anything else */
......@@ -356,7 +356,8 @@ void makeit()
{
/*-------- First: the "measurement" weights */
wfield = newweight(prefs.wimage_name[1],field,prefs.weight_type[1],
nima1<0? ntab : (nweight1<0?1:nweight1));
(nima1<0 && prefs.image_name[1])?
ntab : (nweight1<0?1:nweight1));
wtype = prefs.weight_type[1];
interpthresh = prefs.weight_thresh[1];
/*-------- Convert the interpolation threshold to variance units */
......@@ -437,7 +438,7 @@ void makeit()
}
else if (dwfield && dwfield->flags^INTERP_FIELD)
{
makeback(field, dwfield, BACK_WSCALE);
makeback(field, dwfield, prefs.wscale_flag[0]);
QPRINTF(OUTPUT, "(D) "
"Background: %-10g RMS: %-10g / Threshold: %-10g \n",
field->backmean, field->backsig, field->dthresh);
......
......@@ -679,7 +679,6 @@
&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", ""},
......
......@@ -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/07/2012
* Last modified: 30/10/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -124,7 +124,7 @@
{"NONE","BLANK","CORRECT",""}},
{"MEMORY_BUFSIZE", P_INT, &prefs.mem_bufsize, 8, 65534},
{"MEMORY_OBJSTACK", P_INT, &prefs.clean_stacksize, 16,65536},
{"MEMORY_PIXSTACK", P_INT, &prefs.mem_pixstack, 1000, 10000000},
{"MEMORY_PIXSTACK", P_INT, &prefs.mem_pixstack, 1000, 100000000},
{"NTHREADS", P_INT, &prefs.nthreads, -THREADS_PREFMAX, THREADS_PREFMAX},
{"PARAMETERS_NAME", P_STRING, prefs.param_name},
{"PATTERN_TYPE", P_KEY, &prefs.pattern_type, 0,0, 0.0,0.0,
......
......@@ -7,7 +7,7 @@
*
* This file part of: SExtractor
*
* Copyright: (C) 2006-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 2006-2013 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: 23/09/2012
* Last modified: 16/02/2013
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -111,6 +111,7 @@ profitstruct *profit_init(psfstruct *psf, unsigned int modeltype)
QMALLOC16(profit->lmodpix, PIXTYPE, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE);
QMALLOC16(profit->lmodpix2, PIXTYPE, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE);
QMALLOC16(profit->resi, float, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE);
QMALLOC16(profit->presi, float, profit->nparam);
QMALLOC16(profit->covar, float, profit->nparam*profit->nparam);
profit->nprof = nmodels;
profit->fluxfac = 1.0; /* Default */
......@@ -143,6 +144,7 @@ void profit_end(profitstruct *profit)
free(profit->objpix);
free(profit->objweight);
free(profit->resi);
free(profit->presi);
free(profit->prof);
free(profit->covar);
QFFTWF_FREE(profit->psfdft);
......@@ -166,7 +168,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 12/07/2012
VERSION 26/02/2013
***/
void profit_fit(profitstruct *profit,
picstruct *field, picstruct *wfield,
......@@ -221,6 +223,25 @@ void profit_fit(profitstruct *profit,
profit->subsamp = 1.0;
profit->nobjpix = profit->objnaxisn[0]*profit->objnaxisn[1];
/* Create pixmap at model resolution */
profit->modnaxisn[0] =
((int)(profit->objnaxisn[0]*profit->subsamp/profit->pixstep
+0.4999)/2+1)*2;
profit->modnaxisn[1] =
((int)(profit->objnaxisn[1]*profit->subsamp/profit->pixstep
+0.4999)/2+1)*2;
if (profit->modnaxisn[1] < profit->modnaxisn[0])
profit->modnaxisn[1] = profit->modnaxisn[0];
else
profit->modnaxisn[0] = profit->modnaxisn[1];
if (profit->modnaxisn[0]>PROFIT_MAXMODSIZE)
{
profit->pixstep = (double)profit->modnaxisn[0] / PROFIT_MAXMODSIZE;
profit->modnaxisn[0] = profit->modnaxisn[1] = PROFIT_MAXMODSIZE;
obj2->prof_flag |= PROFLAG_MODSUB;
}
profit->nmodpix = profit->modnaxisn[0]*profit->modnaxisn[1];
/* Use (dirty) global variables to interface with lmfit */
the_field = field;
the_wfield = wfield;
......@@ -228,7 +249,11 @@ void profit_fit(profitstruct *profit,
profit->obj = obj;
profit->obj2 = obj2;
/* Compute the local PSF */
profit_psf(profit);
profit->nresi = profit_copyobjpix(profit, field, wfield);
profit->npresi = 0;
/* Check if the number of constraints exceeds the number of free parameters */
if (profit->nresi < nparam)
{
......@@ -246,28 +271,6 @@ void profit_fit(profitstruct *profit,
return;
}
/* Create pixmap at PSF resolution */
profit->modnaxisn[0] =
((int)(profit->objnaxisn[0]*profit->subsamp/profit->pixstep
+0.4999)/2+1)*2;
profit->modnaxisn[1] =
((int)(profit->objnaxisn[1]*profit->subsamp/profit->pixstep
+0.4999)/2+1)*2;
if (profit->modnaxisn[1] < profit->modnaxisn[0])
profit->modnaxisn[1] = profit->modnaxisn[0];
else
profit->modnaxisn[0] = profit->modnaxisn[1];
if (profit->modnaxisn[0]>PROFIT_MAXMODSIZE)
{
profit->pixstep = (double)profit->modnaxisn[0] / PROFIT_MAXMODSIZE;
profit->modnaxisn[0] = profit->modnaxisn[1] = PROFIT_MAXMODSIZE;
obj2->prof_flag |= PROFLAG_MODSUB;
}
profit->nmodpix = profit->modnaxisn[0]*profit->modnaxisn[1];
/* Compute the local PSF */
profit_psf(profit);
/* Set initial guesses and boundaries */
profit->guesssigbkg = profit->sigma = obj->sigbkg;
profit->guessdx = obj->mx - (int)(obj->mx+0.49999);
......@@ -358,14 +361,28 @@ profit->niter = profit_minimize(profit, PROFIT_MAXITER);
if ((check = prefs.check[CHECK_PROFILES]))
{
profit_residuals(profit,field,wfield, 0.0, profit->paraminit, NULL);
addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1],
if (profit->subsamp>1.0)
addcheck_resample(check, profit->lmodpix,
profit->objnaxisn[0],profit->objnaxisn[1],
profit->ix,profit->iy, 1.0/profit->subsamp,
1.0/(profit->subsamp*profit->subsamp));
else
addcheck(check, profit->lmodpix,
profit->objnaxisn[0],profit->objnaxisn[1],
profit->ix,profit->iy, 1.0);
}
if ((check = prefs.check[CHECK_SUBPROFILES]))
{
profit_residuals(profit,field,wfield, 0.0, profit->paraminit, NULL);
addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1],
if (profit->subsamp>1.0)
addcheck_resample(check, profit->lmodpix,
profit->objnaxisn[0],profit->objnaxisn[1],
profit->ix,profit->iy, 1.0/profit->subsamp,
-1.0/(profit->subsamp*profit->subsamp));
else
addcheck(check, profit->lmodpix,
profit->objnaxisn[0],profit->objnaxisn[1],
profit->ix,profit->iy, -1.0);
}
if ((check = prefs.check[CHECK_SPHEROIDS]))
......@@ -888,7 +905,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 12/07/2012
VERSION 16/02/2013
***/
void profit_dfit(profitstruct *profit, profitstruct *dprofit,
picstruct *field, picstruct *dfield,
......@@ -946,8 +963,9 @@ void profit_dfit(profitstruct *profit, profitstruct *dprofit,
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)
if (dprofit->nresi < nparam || profit->nresi < 1)
{
obj2->dprof_flag |= PROFLAG_NOTCONST;
return;
......@@ -986,8 +1004,8 @@ void profit_dfit(profitstruct *profit, profitstruct *dprofit,
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);
if ((dprofit->guessradius = sqrtf(obj->a*obj->b)*1.17)<0.5*dpsf->fwhm)
dprofit->guessradius = 0.5*dpsf->fwhm;
dprofit->guessaspect = obj->b/obj->a;
dprofit->guessposang = obj->theta;
......@@ -1024,7 +1042,7 @@ void profit_dfit(profitstruct *profit, profitstruct *dprofit,
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;
obj2->fluxerr_dprof = sumd!=0.0? dprofit->flux/sqrtf((float)sumd): 0.0f;
if (FLAG(obj2.dprof_chi2))
{
......@@ -1360,7 +1378,7 @@ INPUT Profile-fitting structure.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 07/07/2010
VERSION 13/02/2013
***/
void profit_psf(profitstruct *profit)
{
......@@ -1403,7 +1421,7 @@ void profit_psf(profitstruct *profit)
}
/* Normalize PSF flux (just in case...) */
flux *= profit->pixstep*profit->pixstep;
flux *= profit->pixstep*profit->pixstep / (profit->subsamp*profit->subsamp);
if (fabs(flux) <= 0.0)
error(EXIT_FAILURE, "*Error*: PSF model is empty or negative: ", psf->name);
......@@ -1424,7 +1442,7 @@ INPUT Pointer to the profit structure involved in the fit,
OUTPUT Number of iterations used.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 20/05/2011
VERSION 15/01/2013
***/
int profit_minimize(profitstruct *profit, int niter)
{
......@@ -1445,7 +1463,8 @@ int profit_minimize(profitstruct *profit, int niter)
nfree = profit_boundtounbound(profit, profit->paraminit, dparam,
PARAM_ALLPARAMS);
niter = dlevmar_dif(profit_evaluate, dparam, NULL, nfree, profit->nresi,
niter = dlevmar_dif(profit_evaluate, dparam, NULL, nfree,
profit->nresi + profit->npresi,
niter, lm_opts, info, NULL, dcovar, profit);
profit_unboundtobound(profit, dparam, profit->paraminit, PARAM_ALLPARAMS);
......@@ -1516,13 +1535,13 @@ INPUT Pointer to the vector of parameters,
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 20/05/2011
VERSION 16/01/2013
***/
void profit_evaluate(double *dpar, double *fvec, int m, int n, void *adata)
{
profitstruct *profit;
profstruct **prof;
double *dpar0, *dresi;
double *dpar0, *dresi, *fvect;
float *modpixt, *profpixt, *resi,
tparam, val;
PIXTYPE *lmodpixt,*lmodpix2t, *objpix,*weight,
......@@ -1682,8 +1701,14 @@ jflag = 0; /* Temporarily deactivated (until problems are fixed) */
profit_residuals(profit, the_field, the_wfield, PROFIT_DYNPARAM,
profit->param, profit->resi);
profit_presiduals(profit, dpar, profit->presi);
fvect = fvec;
for (p=0; p<profit->nresi; p++)
fvec[p] = profit->resi[p];
*(fvect++) = profit->resi[p];
for (p=0; p<profit->npresi; p++)
*(fvect++) = profit->presi[p];
}
// profit_printout(m, par, n, fvec, adata, 0, -1, 0 );
......@@ -1694,7 +1719,7 @@ jflag = 0; /* Temporarily deactivated (until problems are fixed) */
/****** profit_residuals ******************************************************
PROTO float *prof_residuals(profitstruct *profit, picstruct *field,
PROTO float *profit_residuals(profitstruct *profit, picstruct *field,
picstruct *wfield, float dynparam, float *param, float *resi)
PURPOSE Compute the vector of residuals between the data and the galaxy
profile model.
......@@ -1702,7 +1727,7 @@ INPUT Profile-fitting structure,
pointer to the field,
pointer to the field weight,
dynamic compression parameter (0=no compression),
pointer to the model parameters (output),
pointer to the model parameters
pointer to the computed residuals (output).
OUTPUT Vector of residuals.
NOTES -.
......@@ -1742,6 +1767,33 @@ float *profit_residuals(profitstruct *profit, picstruct *field,
}
/****** profit_presiduals *****************************************************
PROTO float *profit_presiduals(profitstruct *profit, float *param,
float *presi)
PURPOSE Compute the vector of prior "residuals" for the model parameters.
INPUT Profile-fitting structure,
pointer to the (unbound) model parameters,
pointer to the computed prior "residuals" (output).
OUTPUT Vector of residuals.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 15/01/2013
***/
float *profit_presiduals(profitstruct *profit, double *dparam, float *presi)
{
float *presit;
int p;
presit = presi;
for (p=0; p<profit->nparam; p++)
if (profit->dparampsig[p]>0.0)
*(presit++) = (float)((dparam[p] - profit->dparampcen[p])
/ profit->dparampsig[p]);
return presi;
}
/****** profit_compresi ******************************************************
PROTO float *profit_compresi(profitstruct *profit, float dynparam,
float *resi)
......@@ -1816,7 +1868,7 @@ INPUT Profile-fitting structure,
OUTPUT RETURN_ERROR if the rasters don't overlap, RETURN_OK otherwise.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 23/09/2012
VERSION 13/02/2013
***/
int profit_resample(profitstruct *profit, float *inpix, PIXTYPE *outpix,
float factor)
......@@ -1832,13 +1884,12 @@ int profit_resample(profitstruct *profit, float *inpix, PIXTYPE *outpix,
iysina, nyin, hmw,hmh, ix,iy, ixin,iyin;
invpixstep = profit->subsamp/profit->pixstep;
factor /= profit->subsamp*profit->subsamp;
xcin = (float)(profit->modnaxisn[0]/2);
xcout = ((int)(profit->subsamp*profit->objnaxisn[0])/2 + 0.5)
/ profit->subsamp - 0.5;
if ((dx=profit->paramlist[PARAM_X]))
xcout += *dx / profit->subsamp;
xcout += *dx/profit->subsamp;
xsin = xcin - xcout*invpixstep; /* Input start x-coord*/
......@@ -1865,7 +1916,7 @@ int profit_resample(profitstruct *profit, float *inpix, PIXTYPE *outpix,
ycout = ((int)(profit->subsamp*profit->objnaxisn[1])/2 + 0.5)
/ profit->subsamp - 0.5;
if ((dy=profit->paramlist[PARAM_Y]))
ycout += *dy / profit->subsamp;
ycout += *dy/profit->subsamp;
ysin = ycin - ycout*invpixstep; /* Input start y-coord*/
if ((int)ysin >= profit->modnaxisn[1])
......@@ -2986,19 +3037,19 @@ PROTO void profit_resetparam(profitstruct *profit, paramenum paramtype)
PURPOSE Set the initial, lower and upper boundary values of a profile parameter.
INPUT Pointer to the profit structure,
Parameter index.
OUTPUT -.
OUTPUT .
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 03/11/2011
VERSION 15/01/2013
***/
void profit_resetparam(profitstruct *profit, paramenum paramtype)
{
obj2struct *obj2;
float param, parammin,parammax, range;
float param, parammin,parammax, range, priorcen,priorsig;
parfitenum fittype;
obj2 = profit->obj2;
param = parammin = parammax = 0.0; /* Avoid gcc -Wall warnings*/
param = parammin = parammax = priorcen = priorsig = 0.0;
switch(paramtype)
{
......@@ -3012,8 +3063,8 @@ void profit_resetparam(profitstruct *profit, paramenum paramtype)
fittype = PARFIT_LINBOUND;
param = profit->guessdx;
range = profit->guessradius*4.0;
if (range>profit->objnaxisn[0]*2.0)
range = profit->objnaxisn[0]*2.0;
if (range>profit->objnaxisn[0]*profit->subsamp*2.0)
range = profit->objnaxisn[0]*profit->subsamp*2.0;
parammin = -range;
parammax = range;
break;
......@@ -3021,8 +3072,8 @@ void profit_resetparam(profitstruct *profit, paramenum paramtype)
fittype = PARFIT_LINBOUND;
param = profit->guessdy;
range = profit->guessradius*4.0;
if (range>profit->objnaxisn[1]*2)
range = profit->objnaxisn[1]*2;
if (range>profit->objnaxisn[1]*profit->subsamp*2)
range = profit->objnaxisn[1]*profit->subsamp*2;
parammin = -range;
parammax = range;
break;
......@@ -3050,6 +3101,8 @@ void profit_resetparam(profitstruct *profit, paramenum paramtype)
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;
priorcen = 0.3;
priorsig = 0.0 /*0.4*/;
break;
case PARAM_SPHEROID_POSANG:
fittype = PARFIT_UNBOUND;
......@@ -3062,6 +3115,8 @@ void profit_resetparam(profitstruct *profit, paramenum paramtype)
param = 4.0;
parammin = FLAG(obj2.prof_disk_flux)? 1.0 : 0.3;
parammax = 10.0;
priorcen = 1.0;
priorsig = 0.0 /*2.0*/;
break;
case PARAM_DISK_FLUX:
fittype = PARFIT_LOGBOUND;
......@@ -3207,7 +3262,8 @@ void profit_resetparam(profitstruct *profit, paramenum paramtype)
param = (parammin+parammax)/2.0;
else if (parammin==0.0 && parammax==0.0)
parammax = 1.0;
profit_setparam(profit, paramtype, param, parammin, parammax, fittype);
profit_setparam(profit, paramtype, param, parammin, parammax, fittype,
priorcen, priorsig);
return;
}
......@@ -3237,7 +3293,8 @@ void profit_resetparams(profitstruct *profit)
/****** profit_setparam ****************************************************
PROTO void profit_setparam(profitstruct *profit, paramenum paramtype,
float param, float parammin, float parammax,
parfitenum parfittype)
parfitenum parfittype,
float priorcen, float priorsig)
PURPOSE Set the actual, lower and upper boundary values of a profile parameter.
INPUT Pointer to the profit structure,
parameter index,
......@@ -3245,14 +3302,18 @@ INPUT Pointer to the profit structure,
lower boundary to the parameter,
upper boundary to the parameter,
parameter fitting type.
prior central value,
prior standard deviation.
OUTPUT RETURN_OK if the parameter is registered, RETURN_ERROR otherwise.
AUTHOR E. Bertin (IAP)
VERSION 20/05/2011
VERSION 16/01/2013
***/
int profit_setparam(profitstruct *profit, paramenum paramtype,
float param, float parammin,float parammax,
parfitenum parfittype)
parfitenum parfittype,
float priorcen, float priorsig)
{
double dtemp;
float *paramptr;
int index;
......@@ -3264,6 +3325,8 @@ int profit_setparam(profitstruct *profit, paramenum paramtype,
profit->parammin[index] = parammin;
profit->parammax[index] = parammax;
profit->parfittype[index] = parfittype;
profit_boundtounbound(profit, &priorcen, &profit->dparampcen[index], index);
profit->dparampsig[index] = priorsig;
return RETURN_OK;
}
else
......@@ -3758,7 +3821,7 @@ INPUT Profile-fitting structure,
OUTPUT Total (asymptotic) flux contribution.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 25/07/2011
VERSION 13/02/2013
***/
float prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag)
{
......
......@@ -7,7 +7,7 @@
*
* This file part of: SExtractor
*
* Copyright: (C) 2006-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
* Copyright: (C) 2006-2013 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: 08/12/2011
* Last modified: 13/02/2013
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
......@@ -157,6 +157,8 @@ typedef struct
float paraminit[PARAM_NPARAM];/* Parameter initial guesses */
float parammin[PARAM_NPARAM]; /* Parameter lower limits */
float parammax[PARAM_NPARAM]; /* Parameter upper limits */
double dparampcen[PARAM_NPARAM];/* Parameter prior center */
double dparampsig[PARAM_NPARAM];/* Parameter prior sigma */
int nlimmin; /* # of parameters that hit the min limit */
int nlimmax; /* # of parameters that hit the max limit */
float paramerr[PARAM_NPARAM]; /* Std deviations of parameters */
......@@ -185,6 +187,8 @@ typedef struct
int ix, iy; /* Integer coordinates of object pixmap */
float *resi; /* Vector of residuals */
int nresi; /* Number of residual elements */
float *presi; /* Vector of prior residuals */
int npresi; /* Number of prior residual elements */
float chi2; /* Std error per residual element */
float sigma; /* Standard deviation of the pixel values */
float flux; /* Total flux in final convolved model */
......@@ -209,6 +213,8 @@ profstruct *prof_init(profitstruct *profit, unsigned int modeltype);
float *profit_compresi(profitstruct *profit, float dynparam,
float *resi),
*profit_presiduals(profitstruct *profit, double *dparam,
float *presi),
*profit_residuals(profitstruct *profit, picstruct *field,
picstruct *wfield, float dynparam,
float *param, float *resi),
......@@ -231,7 +237,8 @@ int profit_boundtounbound(profitstruct *profit,
PIXTYPE *outpix, float factor),
profit_setparam(profitstruct *profit, paramenum paramtype,
float param, float parammin, float parammax,
parfitenum parfittype),
parfitenum parfittype,
float priorcen, float priorsig),
profit_unboundtobound(profitstruct *profit,
double *dparam, float *param, int index);
......
......@@ -127,6 +127,7 @@ psfstruct *psf_load(char *filename, int ext)
char *head, *ci,*co;
int deg[POLY_MAXDIM], group[POLY_MAXDIM], ndim, ngroup,
e,i,k;
/* Open the cat (well it is not a "cat", but simply a FITS file */
if (!(cat = read_cat(filename)))
error(EXIT_FAILURE, "*Error*: PSF file not found: ", filename);
......
......@@ -500,10 +500,10 @@ typedef struct
float prof_arms_scaleerrw; /* RMS error */
float prof_arms_posang; /* Arms true position angle */
float prof_arms_posangerr; /* RMS error */
// float prof_arms_thetaw; /* WORLD arms position angle */
// float prof_arms_thetas; /* Sky arms position angle */
// float prof_arms_theta2000; /* J2000 arms position angle */
// float prof_arms_theta1950; /* B1950 arms position angle */
float prof_arms_thetaw; /* WORLD arms position angle */
float prof_arms_thetas; /* Sky arms position angle */
float prof_arms_theta2000; /* J2000 arms position angle */
float prof_arms_theta1950; /* B1950 arms position angle */
float prof_arms_pitch; /* Arms pitch angle */
float prof_arms_pitcherr; /* RMS error */
float prof_arms_start; /* Arms starting radius */
......
# Makefile.in generated by automake 1.11.1 from Makefile.am.
# Makefile.in generated by automake 1.12.2 from Makefile.am.
# @configure_input@
# Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002,
# 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation,
# Inc.
# Copyright (C) 1994-2012 Free Software Foundation, Inc.
# This Makefile.in is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
# with or without modifications, as long as this notice is preserved.
......@@ -47,6 +46,23 @@
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
VPATH = @srcdir@
am__make_dryrun = \
{ \
am__dry=no; \
case $$MAKEFLAGS in \
*\\[\ \ ]*) \
echo 'am--echo: ; @echo "AM" OK' | $(MAKE) -f - 2>/dev/null \
| grep '^AM OK$$' >/dev/null || am__dry=yes;; \
*) \
for am__flg in $$MAKEFLAGS; do \
case $$am__flg in \
*=*|--*) ;; \
*n*) am__dry=yes; break;; \
esac; \
done;; \
esac; \
test $$am__dry = yes; \
}
pkgdatadir = $(datadir)/@PACKAGE@
pkgincludedir = $(includedir)/@PACKAGE@
pkglibdir = $(libdir)/@PACKAGE@
......@@ -66,7 +82,9 @@ POST_UNINSTALL = :
build_triplet = @build@
host_triplet = @host@
subdir = src/wcs
DIST_COMMON = README $(srcdir)/Makefile.am $(srcdir)/Makefile.in
DIST_COMMON = README $(srcdir)/Makefile.am $(srcdir)/Makefile.in \
$(top_srcdir)/autoconf/depcomp \
$(top_srcdir)/autoconf/mkinstalldirs
ACLOCAL_M4 = $(top_srcdir)/aclocal.m4
am__aclocal_m4_deps = $(top_srcdir)/acx_atlas.m4 \
$(top_srcdir)/acx_fftw.m4 $(top_srcdir)/acx_mkl.m4 \
......@@ -103,6 +121,11 @@ LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) \
$(LDFLAGS) -o $@
SOURCES = $(libwcs_c_a_SOURCES)
DIST_SOURCES = $(libwcs_c_a_SOURCES)
am__can_run_installinfo = \
case $$AM_UPDATE_INFO_DIR in \
n|no|NO) false;; \
*) (install-info --version) >/dev/null 2>&1;; \
esac
ETAGS = etags
CTAGS = ctags
DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST)
......@@ -161,6 +184,7 @@ MAKEINFO = @MAKEINFO@
MANIFEST_TOOL = @MANIFEST_TOOL@
MKDIR_P = @MKDIR_P@
MKL_CFLAGS = @MKL_CFLAGS@
MKL_LDFLAGS = @MKL_LDFLAGS@
MKL_LIBS = @MKL_LIBS@
MKL_WARN = @MKL_WARN@
NM = @NM@
......@@ -283,7 +307,7 @@ $(am__aclocal_m4_deps):
clean-noinstLIBRARIES:
-test -z "$(noinst_LIBRARIES)" || rm -f $(noinst_LIBRARIES)
libwcs_c.a: $(libwcs_c_a_OBJECTS) $(libwcs_c_a_DEPENDENCIES)
libwcs_c.a: $(libwcs_c_a_OBJECTS) $(libwcs_c_a_DEPENDENCIES) $(EXTRA_libwcs_c_a_DEPENDENCIES)
-rm -f libwcs_c.a
$(libwcs_c_a_AR) libwcs_c.a $(libwcs_c_a_OBJECTS) $(libwcs_c_a_LIBADD)
$(RANLIB) libwcs_c.a
......@@ -379,6 +403,20 @@ GTAGS:
&& $(am__cd) $(top_srcdir) \
&& gtags -i $(GTAGS_ARGS) "$$here"
cscopelist: $(HEADERS) $(SOURCES) $(LISP)
list='$(SOURCES) $(HEADERS) $(LISP)'; \
case "$(srcdir)" in \
[\\/]* | ?:[\\/]*) sdir="$(srcdir)" ;; \
*) sdir=$(subdir)/$(srcdir) ;; \
esac; \
for i in $$list; do \
if test -f "$$i"; then \
echo "$(subdir)/$$i"; \
else \
echo "$$sdir/$$i"; \
fi; \
done >> $(top_builddir)/cscope.files
distclean-tags:
-rm -f TAGS ID GTAGS GRTAGS GSYMS GPATH tags
......@@ -426,10 +464,15 @@ install-am: all-am
installcheck: installcheck-am
install-strip:
if test -z '$(STRIP)'; then \
$(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM="$(INSTALL_STRIP_PROGRAM)" \
install_sh_PROGRAM="$(INSTALL_STRIP_PROGRAM)" INSTALL_STRIP_FLAG=-s \
install; \
else \
$(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM="$(INSTALL_STRIP_PROGRAM)" \
install_sh_PROGRAM="$(INSTALL_STRIP_PROGRAM)" INSTALL_STRIP_FLAG=-s \
`test -z '$(STRIP)' || \
echo "INSTALL_PROGRAM_ENV=STRIPPROG='$(STRIP)'"` install
"INSTALL_PROGRAM_ENV=STRIPPROG='$(STRIP)'" install; \
fi
mostlyclean-generic:
clean-generic:
......@@ -515,7 +558,7 @@ uninstall-am:
.MAKE: install-am install-strip
.PHONY: CTAGS GTAGS all all-am check check-am clean clean-generic \
clean-libtool clean-noinstLIBRARIES ctags distclean \
clean-libtool clean-noinstLIBRARIES cscopelist ctags distclean \
distclean-compile distclean-generic distclean-libtool \
distclean-tags distdir dvi dvi-am html html-am info info-am \
install install-am install-data install-data-am install-dvi \
......
# Makefile.in generated by automake 1.11.1 from Makefile.am.
# Makefile.in generated by automake 1.12.2 from Makefile.am.
# @configure_input@
# Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002,
# 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation,
# Inc.
# Copyright (C) 1994-2012 Free Software Foundation, Inc.
# This Makefile.in is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
# with or without modifications, as long as this notice is preserved.
......@@ -15,6 +14,23 @@
@SET_MAKE@
VPATH = @srcdir@
am__make_dryrun = \
{ \
am__dry=no; \
case $$MAKEFLAGS in \
*\\[\ \ ]*) \
echo 'am--echo: ; @echo "AM" OK' | $(MAKE) -f - 2>/dev/null \
| grep '^AM OK$$' >/dev/null || am__dry=yes;; \
*) \
for am__flg in $$MAKEFLAGS; do \
case $$am__flg in \
*=*|--*) ;; \
*n*) am__dry=yes; break;; \
esac; \
done;; \
esac; \
test $$am__dry = yes; \
}
pkgdatadir = $(datadir)/@PACKAGE@
pkgincludedir = $(includedir)/@PACKAGE@
pkglibdir = $(libdir)/@PACKAGE@
......@@ -34,7 +50,8 @@ POST_UNINSTALL = :
build_triplet = @build@
host_triplet = @host@
subdir = tests
DIST_COMMON = $(srcdir)/Makefile.am $(srcdir)/Makefile.in
DIST_COMMON = $(srcdir)/Makefile.am $(srcdir)/Makefile.in \
$(top_srcdir)/autoconf/mkinstalldirs
ACLOCAL_M4 = $(top_srcdir)/aclocal.m4
am__aclocal_m4_deps = $(top_srcdir)/acx_atlas.m4 \
$(top_srcdir)/acx_fftw.m4 $(top_srcdir)/acx_mkl.m4 \
......@@ -50,8 +67,15 @@ CONFIG_CLEAN_FILES =
CONFIG_CLEAN_VPATH_FILES =
SOURCES =
DIST_SOURCES =
am__tty_colors = \
red=; grn=; lgn=; blu=; std=
am__can_run_installinfo = \
case $$AM_UPDATE_INFO_DIR in \
n|no|NO) false;; \
*) (install-info --version) >/dev/null 2>&1;; \
esac
am__tty_colors_dummy = \
mgn= red= grn= lgn= blu= brg= std=; \
am__color_tests=no
am__tty_colors = $(am__tty_colors_dummy)
DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST)
ACLOCAL = @ACLOCAL@
AMTAR = @AMTAR@
......@@ -108,6 +132,7 @@ MAKEINFO = @MAKEINFO@
MANIFEST_TOOL = @MANIFEST_TOOL@
MKDIR_P = @MKDIR_P@
MKL_CFLAGS = @MKL_CFLAGS@
MKL_LDFLAGS = @MKL_LDFLAGS@
MKL_LIBS = @MKL_LIBS@
MKL_WARN = @MKL_WARN@
NM = @NM@
......@@ -239,6 +264,8 @@ TAGS:
ctags: CTAGS
CTAGS:
cscope cscopelist:
check-TESTS: $(TESTS)
@failed=0; all=0; xfail=0; xpass=0; skip=0; \
......@@ -250,7 +277,7 @@ check-TESTS: $(TESTS)
if test -f ./$$tst; then dir=./; \
elif test -f $$tst; then dir=; \
else dir="$(srcdir)/"; fi; \
if $(TESTS_ENVIRONMENT) $${dir}$$tst; then \
if $(TESTS_ENVIRONMENT) $${dir}$$tst $(AM_TESTS_FD_REDIRECT); then \
all=`expr $$all + 1`; \
case " $(XFAIL_TESTS) " in \
*[\ \ ]$$tst[\ \ ]*) \
......@@ -321,14 +348,15 @@ check-TESTS: $(TESTS)
fi; \
dashes=`echo "$$dashes" | sed s/./=/g`; \
if test "$$failed" -eq 0; then \
echo "$$grn$$dashes"; \
col="$$grn"; \
else \
echo "$$red$$dashes"; \
col="$$red"; \
fi; \
echo "$$banner"; \
test -z "$$skipped" || echo "$$skipped"; \
test -z "$$report" || echo "$$report"; \
echo "$$dashes$$std"; \
echo "$${col}$$dashes$${std}"; \
echo "$${col}$$banner$${std}"; \
test -z "$$skipped" || echo "$${col}$$skipped$${std}"; \
test -z "$$report" || echo "$${col}$$report$${std}"; \
echo "$${col}$$dashes$${std}"; \
test "$$failed" -eq 0; \
else :; fi
......@@ -377,10 +405,15 @@ install-am: all-am
installcheck: installcheck-am
install-strip:
if test -z '$(STRIP)'; then \
$(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM="$(INSTALL_STRIP_PROGRAM)" \
install_sh_PROGRAM="$(INSTALL_STRIP_PROGRAM)" INSTALL_STRIP_FLAG=-s \
install; \
else \
$(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM="$(INSTALL_STRIP_PROGRAM)" \
install_sh_PROGRAM="$(INSTALL_STRIP_PROGRAM)" INSTALL_STRIP_FLAG=-s \
`test -z '$(STRIP)' || \
echo "INSTALL_PROGRAM_ENV=STRIPPROG='$(STRIP)'"` install
"INSTALL_PROGRAM_ENV=STRIPPROG='$(STRIP)'" install; \
fi
mostlyclean-generic:
clean-generic:
......
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