From 5f63711aaf01a1ee403cbfc4fa60170c78869573 Mon Sep 17 00:00:00 2001 From: Emmanuel Bertin Date: Fri, 9 Jul 2010 14:41:51 +0000 Subject: [PATCH] Added VOTABLE "version" attribute to XML output. Reverted minimization algorithm back to double precision (while keeping single precision for data and parameters in "bounded" space). Added specific handling of "Jacobian calls" in model-fitting. Sped up the resampling process in model-fitting. Tuned up minimization algorithm in model-fitting. Added shortcuts to the computation of the local PSF model. Fixed triggering issues with some model-fitting parameters (thanks to S.Serrano). Removed CLASS_STAR_MODEL and DISK_PATTERN* measurement parameters. Added automatic tracking of the PSF FWHM from the local PSF model if SEEING_FWHM is 0 for CLASS_STAR. Fixed excessive WIN parameter computation time on large ring-shaped objects. Fixed model-fitting issue with null half-light radii (thanks to S.Desai and B.Armstrong). Fixed aperture display value error in FITS-LDAC catalogue header. Removed redundant update_tab() calls in FITS library. Fixed comment copy issue for slashes within strings in FITS library (thanks to F.Schuller). Added check of PSF values; SExtractor will now exit in error in the case where all PSF components sum up to 0. Replaced computation of model second moments and ellipticities with analytical estimates. Improved the accuracy of peak and average model surface brightness estimates. Fixed some compiler warnings. Pushed version number to 2.11.7. --- configure | 20 +- configure.ac | 4 +- man/sex.1 | 2 +- src/analyse.c | 8 +- src/catout.c | 10 +- src/check.c | 50 +- src/fits/fitsutil.c | 6 +- src/fits/fitswrite.c | 8 +- src/growth.c | 12 +- src/growth.h | 6 +- src/param.h | 4 +- src/paramprofit.h | 8 +- src/preflist.h | 11 +- src/prefs.c | 16 +- src/profit.c | 1280 +++++++++++++++++++++++++++++++----------- src/profit.h | 46 +- src/psf.c | 40 +- src/psf.h | 7 +- src/readimage.c | 10 +- src/sexhead1.h | 6 +- src/types.h | 10 +- src/xml.c | 6 +- 22 files changed, 1149 insertions(+), 421 deletions(-) diff --git a/configure b/configure index 0d81f85..254115e 100755 --- a/configure +++ b/configure @@ -1,6 +1,6 @@ #! /bin/sh # Guess values for system-dependent variables and create Makefiles. -# Generated by GNU Autoconf 2.63 for sextractor 2.11.0. +# Generated by GNU Autoconf 2.63 for sextractor 2.11.7. # # Report bugs to . # @@ -750,8 +750,8 @@ SHELL=${CONFIG_SHELL-/bin/sh} # Identity of this package. PACKAGE_NAME='sextractor' PACKAGE_TARNAME='sextractor' -PACKAGE_VERSION='2.11.0' -PACKAGE_STRING='sextractor 2.11.0' +PACKAGE_VERSION='2.11.7' +PACKAGE_STRING='sextractor 2.11.7' PACKAGE_BUGREPORT='bertin@iap.fr' ac_unique_file="src/makeit.c" @@ -1505,7 +1505,7 @@ if test "$ac_init_help" = "long"; then # Omit some internal or obsolete options to make the list less imposing. # This message is too long to be a string in the A/UX 3.1 sh. cat <<_ACEOF -\`configure' configures sextractor 2.11.0 to adapt to many kinds of systems. +\`configure' configures sextractor 2.11.7 to adapt to many kinds of systems. Usage: $0 [OPTION]... [VAR=VALUE]... @@ -1575,7 +1575,7 @@ fi if test -n "$ac_init_help"; then case $ac_init_help in - short | recursive ) echo "Configuration of sextractor 2.11.0:";; + short | recursive ) echo "Configuration of sextractor 2.11.7:";; esac cat <<\_ACEOF @@ -1706,7 +1706,7 @@ fi test -n "$ac_init_help" && exit $ac_status if $ac_init_version; then cat <<\_ACEOF -sextractor configure 2.11.0 +sextractor configure 2.11.7 generated by GNU Autoconf 2.63 Copyright (C) 1992, 1993, 1994, 1995, 1996, 1998, 1999, 2000, 2001, @@ -1720,7 +1720,7 @@ cat >config.log <<_ACEOF This file contains any messages produced by compilers while running configure, to aid debugging if configure makes a mistake. -It was created by sextractor $as_me 2.11.0, which was +It was created by sextractor $as_me 2.11.7, which was generated by GNU Autoconf 2.63. Invocation command line was $ $0 $@ @@ -2423,7 +2423,7 @@ fi # Define the identity of the package. PACKAGE='sextractor' - VERSION='2.11.0' + VERSION='2.11.7' cat >>confdefs.h <<_ACEOF @@ -28535,7 +28535,7 @@ exec 6>&1 # report actual input values of CONFIG_FILES etc. instead of their # values after options handling. ac_log=" -This file was extended by sextractor $as_me 2.11.0, which was +This file was extended by sextractor $as_me 2.11.7, which was generated by GNU Autoconf 2.63. Invocation command line was CONFIG_FILES = $CONFIG_FILES @@ -28598,7 +28598,7 @@ Report bugs to ." _ACEOF cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1 ac_cs_version="\\ -sextractor config.status 2.11.0 +sextractor config.status 2.11.7 configured by $0, generated by GNU Autoconf 2.63, with options \\"`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`\\" diff --git a/configure.ac b/configure.ac index f6df79b..d24c5c2 100644 --- a/configure.ac +++ b/configure.ac @@ -1,12 +1,12 @@ # configure.in for SExtractor -# (C) E.Bertin 2002-2009 +# (C) E.Bertin 2002-2010 # Process this file with autoconf to produce a configure script. # First, disable the annoying config.cache define([AC_CACHE_LOAD],) define([AC_CACHE_SAVE],) # This is your standard Bertin source code... -AC_INIT(sextractor, 2.11.0, [bertin@iap.fr]) +AC_INIT(sextractor, 2.11.7, [bertin@iap.fr]) AC_CONFIG_SRCDIR(src/makeit.c) AC_CONFIG_AUX_DIR(autoconf) AM_CONFIG_HEADER(config.h) diff --git a/man/sex.1 b/man/sex.1 index 01c11dd..d6e243b 100644 --- a/man/sex.1 +++ b/man/sex.1 @@ -1,4 +1,4 @@ -.TH SEXTRACTOR "1" "February 2010" "SWarp 2.11.0" "User Commands" +.TH SEXTRACTOR "1" "July 2010" "SWarp 2.11.7" "User Commands" .SH NAME sex \- extract a source catalog from an astronomical FITS image .SH SYNOPSIS diff --git a/src/analyse.c b/src/analyse.c index e2b6af2..7bc8765 100644 --- a/src/analyse.c +++ b/src/analyse.c @@ -395,6 +395,11 @@ void endobject(picstruct *field, picstruct *dfield, picstruct *wfield, analtime1; int i,j, ix,iy,selecflag, newnumber,nsub; + if (prefs.psf_flag || prefs.prof_flag) + thepsf->build_flag = 0; /* Reset PSF building flag */ + if (prefs.dpsf_flag) + ppsf->build_flag = 0; /* Reset PSF building flag */ + if (FLAG(obj2.analtime)) analtime1 = counter_seconds(); else @@ -589,7 +594,8 @@ void endobject(picstruct *field, picstruct *dfield, picstruct *wfield, { double fac2, input[10], output, fwhm; - fwhm = prefs.seeing_fwhm; + fwhm = (prefs.seeing_fwhm==0.0)? psf_fwhm(thepsf)*field->pixscale + : prefs.seeing_fwhm; fac2 = fwhm/field->pixscale; fac2 *= fac2; diff --git a/src/catout.c b/src/catout.c index 65eb970..4cc3af0 100644 --- a/src/catout.c +++ b/src/catout.c @@ -9,7 +9,7 @@ * * Contents: functions for output of catalog data. * -* Last modify: 07/10/2009 +* Last modify: 08/07/2010 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -193,6 +193,9 @@ void updateparamflags() | FLAG(obj2.poserraw_prof); FLAG(obj2.poserrcxx_prof) |= FLAG(obj2.poserrcyy_prof) | FLAG(obj2.poserrcxy_prof); + FLAG(obj2.poserrmx2_prof) |= FLAG(obj2.poserrmy2_prof) + | FLAG(obj2.poserrmxy_prof) + | FLAG(obj2.poserrmx2w_prof); FLAG(obj2.alpha1950_prof) |= FLAG(obj2.delta1950_prof) | FLAG(obj2.poserrtheta1950_prof); FLAG(obj2.alpha2000_prof) |= FLAG(obj2.delta2000_prof) @@ -208,6 +211,7 @@ void updateparamflags() | FLAG(obj2.xf_prof) | FLAG(obj2.poserra_prof) | FLAG(obj2.poserrcxx_prof) + | FLAG(obj2.poserrmx2_prof) | FLAG(obj2.prof_concentration) | FLAG(obj2.prof_class_star); @@ -336,12 +340,10 @@ void updateparamflags() | FLAG(obj2.prof_e1) | FLAG(obj2.prof_a) | FLAG(obj2.prof_cxx) | FLAG(obj2.prof_mx2w); - FLAG(obj2.mumax_prof) |= FLAG(obj2.mueff_prof) | FLAG(obj2.mumean_prof);/*!*/ FLAG(obj2.fluxmean_prof) |= FLAG(obj2.mumean_prof); FLAG(obj2.fluxeff_prof) |= FLAG(obj2.mueff_prof) | FLAG(obj2.fluxmean_prof); - FLAG(obj2.peak_prof) |= FLAG(obj2.mumax_prof) - | FLAG(obj2.fluxeff_prof); + FLAG(obj2.peak_prof) |= FLAG(obj2.mumax_prof); FLAG(obj2.prof_arms_flux) |= FLAG(obj2.prof_arms_fluxerr) | FLAG(obj2.prof_arms_mag) diff --git a/src/check.c b/src/check.c index 54e7f69..9f49e04 100644 --- a/src/check.c +++ b/src/check.c @@ -9,7 +9,7 @@ * * Contents: handling of "check-images". * -* Last modify: 17/09/2008 +* Last modify: 23/03/2010 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -236,13 +236,17 @@ void reinitcheck(picstruct *field, checkstruct *check) case CHECK_OBJECTS: case CHECK_APERTURES: - case CHECK_SUBPSFPROTOS: case CHECK_PSFPROTOS: - case CHECK_SUBPCPROTOS: + case CHECK_SUBPSFPROTOS: case CHECK_PCPROTOS: + case CHECK_SUBPCPROTOS: case CHECK_PCOPROTOS: - case CHECK_SUBPROFILES: case CHECK_PROFILES: + case CHECK_SUBPROFILES: + case CHECK_SPHEROIDS: + case CHECK_SUBSPHEROIDS: + case CHECK_DISKS: + case CHECK_SUBDISKS: case CHECK_PATTERNS: ival = -32; fitswrite(check->fitshead, "BITPIX ", &ival, H_INT, T_LONG); @@ -349,6 +353,20 @@ void reinitcheck(picstruct *field, checkstruct *check) free(check->fitshead); break; + case CHECK_OTHER: + ival = -32; + fitswrite(check->fitshead, "BITPIX ", &ival, H_INT, T_LONG); + fitswrite(check->fitshead, "NAXIS1 ", &check->width, H_INT, T_LONG); + fitswrite(check->fitshead, "NAXIS2 ", &check->height, H_INT, T_LONG); + check->npix = check->width*check->height; + QMALLOC(ptrf, PIXTYPE, check->npix); + check->pix = (void *)ptrf; + for (i=check->npix; i--;) + *(ptrf++) = -10.0; + QFWRITE(check->fitshead,check->fitsheadsize,check->file,check->filename); + free(check->fitshead); + break; + default: error(EXIT_FAILURE, "*Internal Error* in ", "reinitcheck()!"); } @@ -366,7 +384,8 @@ void writecheck(checkstruct *check, PIXTYPE *data, int w) { if (check->type == CHECK_APERTURES || check->type == CHECK_SUBPSFPROTOS || check->type == CHECK_SUBPCPROTOS || check->type == CHECK_PCOPROTOS - || check->type == CHECK_SUBPROFILES) + || check->type == CHECK_SUBPROFILES || check->type == CHECK_SUBSPHEROIDS + || check->type == CHECK_SUBDISKS) { memcpy((PIXTYPE *)check->pix + w*(check->y++), data, w*sizeof(PIXTYPE)); return; @@ -422,15 +441,21 @@ void reendcheck(picstruct *field, checkstruct *check) case CHECK_OBJECTS: case CHECK_APERTURES: - case CHECK_SUBPSFPROTOS: case CHECK_PSFPROTOS: - case CHECK_SUBPCPROTOS: + case CHECK_SUBPSFPROTOS: case CHECK_PCPROTOS: + case CHECK_SUBPCPROTOS: case CHECK_PCOPROTOS: - case CHECK_SUBPROFILES: case CHECK_PROFILES: + case CHECK_SUBPROFILES: + case CHECK_SPHEROIDS: + case CHECK_SUBSPHEROIDS: + case CHECK_DISKS: + case CHECK_SUBDISKS: case CHECK_ASSOC: case CHECK_PATTERNS: + case CHECK_MAPSOM: + case CHECK_OTHER: if (bswapflag) swapbytes(check->pix, sizeof(PIXTYPE), (int)check->npix); QFWRITE(check->pix,check->npix*sizeof(PIXTYPE), @@ -461,15 +486,6 @@ void reendcheck(picstruct *field, checkstruct *check) break; } - case CHECK_MAPSOM: - if (bswapflag) - swapbytes(check->pix, sizeof(PIXTYPE), (int)check->npix); - QFWRITE(check->pix,check->npix*sizeof(PIXTYPE), - check->file,check->filename); - free(check->pix); - padsize = (FBSIZE -((check->npix*sizeof(USHORT))%FBSIZE)) % FBSIZE; - break; - default: error(EXIT_FAILURE, "*Internal Error* in ", "endcheck()!"); } diff --git a/src/fits/fitsutil.c b/src/fits/fitsutil.c index cdecceb..0a32a33 100644 --- a/src/fits/fitsutil.c +++ b/src/fits/fitsutil.c @@ -9,7 +9,7 @@ * * Contents: functions for handling FITS keywords. * -* Last modify: 02/11/2009 +* Last modify: 02/07/2010 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -181,7 +181,7 @@ OUTPUT RETURN_OK if something was found, RETURN_ERROR otherwise. NOTES -. AUTHOR E. Bertin (IAP), E.R. Deul - Handling of NaN -VERSION 18/05/2009 +VERSION 02/07/2010 ***/ int fitspick(char *fitsline, char *keyword, void *ptr, h_type *htype, t_type *ttype, char *comment) @@ -274,7 +274,7 @@ int fitspick(char *fitsline, char *keyword, void *ptr, h_type *htype, for (fptr = fitsline + (i=j); i<80; i++) { if (*fptr == (char)'\'') - toggle^=toggle; + toggle ^= 1; if (*(fptr++) == (char)'/' && !toggle) { while (++i<80 && *fptr<=' ') diff --git a/src/fits/fitswrite.c b/src/fits/fitswrite.c index a00b249..3f15ca9 100644 --- a/src/fits/fitswrite.c +++ b/src/fits/fitswrite.c @@ -75,7 +75,7 @@ INPUT pointer to the catalog structure, OUTPUT -. NOTES -. AUTHOR E. Bertin (IAP & Leiden observatory) -VERSION 09/09/2003 +VERSION 02/11/2009 ***/ void save_tab(catstruct *cat, tabstruct *tab) @@ -90,8 +90,6 @@ void save_tab(catstruct *cat, tabstruct *tab) char *buf, *inbuf, *outbuf, *fptr,*ptr; int esize; -/* Make the table parameters reflect its content*/ - update_tab(tab); /* The header itself*/ tabflag = save_head(cat, tab)==RETURN_OK?1:0; /* Allocate memory for the output buffer */ @@ -265,7 +263,7 @@ INPUT catalog structure, OUTPUT -. NOTES -. AUTHOR E. Bertin (IAP & Leiden observatory) -VERSION 26/09/2004 +VERSION 02/11/2009 ***/ void init_writeobj(catstruct *cat, tabstruct *tab, char **pbuf) @@ -273,8 +271,6 @@ void init_writeobj(catstruct *cat, tabstruct *tab, char **pbuf) keystruct *key; int k; -/* Make the table parameters reflect its content*/ - update_tab(tab); /* The header itself */ if (save_head(cat, tab) != RETURN_OK) error(EXIT_FAILURE, "*Error*: Not a binary table: ", tab->extname); diff --git a/src/growth.c b/src/growth.c index cdd281e..8e02891 100644 --- a/src/growth.c +++ b/src/growth.c @@ -9,7 +9,7 @@ * * Contents: Make growth curves. * -* Last modify: 19/12/2007 +* Last modify: 02/07/2010 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -336,6 +336,8 @@ void makeavergrowth(picstruct *field, picstruct *wfield, objstruct *obj) i + (tv - *(growtht-1))/dg : i) : (*growth !=0.0 ?tv/(*growth) : 0.0)); + if (obj2->flux_radius[j] > rlim) + obj2->flux_radius[j] = rlim; } } @@ -346,10 +348,14 @@ void makeavergrowth(picstruct *field, picstruct *wfield, objstruct *obj) tv = 0.5*obj2->flux_auto; growtht = growth-1; for (i=0; ihl_radius = step*(i? ((dg=*growtht - *(growtht-1)) != 0.0 ? + obj2->hl_radius = fabs(step*(i? ((dg=*growtht - *(growtht-1)) != 0.0 ? i + (tv - *(growtht-1))/dg : i) - : (*growth !=0.0 ?tv/(*growth) : 0.0)); + : (*growth !=0.0 ?tv/(*growth) : 0.0))); + if (obj2->hl_radius > rlim) + obj2->hl_radius = rlim; + if (obj2->hl_radius < GROWTH_MINHLRAD) + obj2->hl_radius = GROWTH_MINHLRAD; } return; diff --git a/src/growth.h b/src/growth.h index fd79452..aadcfce 100644 --- a/src/growth.h +++ b/src/growth.h @@ -9,7 +9,7 @@ * * Contents: Include file for growth.c. * -* Last modify: 04/05/98 +* Last modify: 02/07/2010 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -18,8 +18,8 @@ #define GROWTH_NSTEP 64 /* number of growth curve samples */ #define GROWTH_OVERSAMP 5 /* pixel oversampling in each dimension */ -#define GROWTH_NSIG 3*MARGIN_SCALE /* MAG_AUTO analysis range (number */ - /* of sigma) */ +#define GROWTH_NSIG 3*MARGIN_SCALE /* MAG_AUTO analysis range (nsigmas) */ +#define GROWTH_MINHLRAD 0.5 /* Minimum internal half-light radius (pixels)*/ /* NOTES: One must have: GROWTH_SAMP >= 1 diff --git a/src/param.h b/src/param.h index 4490ee6..cbb2d52 100644 --- a/src/param.h +++ b/src/param.h @@ -9,7 +9,7 @@ * * Contents: parameter list for catalog data. * -* Last modify: 14/12/2009 +* Last modify: 18/05/2010 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -702,7 +702,7 @@ keystruct objkey[] = { &outobj2.win_polarw, H_FLOAT, T_FLOAT, "%7.5f", "", "src.ellipticity", ""}, {"CLASS_STAR", "S/G classifier output", - &outobj2.sprob, H_FLOAT, T_FLOAT, "%5.2f", "", + &outobj2.sprob, H_FLOAT, T_FLOAT, "%6.3f", "", "src.class.starGalaxy", ""}, {"VIGNET", "Pixel data around detection", &outobj2.vignet, H_FLOAT, T_FLOAT, "%12.7g", "count", diff --git a/src/paramprofit.h b/src/paramprofit.h index e16f8ba..9344c87 100644 --- a/src/paramprofit.h +++ b/src/paramprofit.h @@ -9,7 +9,7 @@ * * Contents: Model-fitting parameter list for catalog data. * -* Last modify: 01/12/2009 +* Last modify: 25/05/2010 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -272,10 +272,11 @@ {"SPREADERR_MODEL", "Spread parameter error from model-fitting", &outobj2.prof_concentrationerr, H_FLOAT, T_FLOAT, "%8.5f", "", "src.morph.param", ""}, +/* {"CLASS_STAR_MODEL", "S/G classifier from model-fitting", &outobj2.prof_class_star, H_FLOAT, T_FLOAT, "%7.4f", "", "src.class.starGalaxy", ""}, - +*/ {"FLUX_BACKOFFSET", "Background offset from fitting", &outobj2.prof_offset_flux, H_FLOAT, T_FLOAT, "%12.7g", "count", "instr.skyLevel;arith.diff;stat.fit.param", "ct"}, @@ -446,6 +447,7 @@ {"DISK_THETA_B1950", "Disk position angle (east of north, B1950)", &outobj2.prof_disk_theta1950, H_FLOAT, T_FLOAT, "%+7.3f", "deg", "pos.posAng;src.morph;stat.fit.param", "deg"}, +/* {"DISK_PATTERN_VECTOR", "Disk pattern fitted coefficients", &outobj2.prof_disk_patternvector, H_FLOAT, T_FLOAT, "%12.4g", "", "stat.fit.param;src.morph.param", "", 1, @@ -461,7 +463,7 @@ {"DISK_PATTERN_SPIRAL", "Disk pattern spiral index", &outobj2.prof_disk_patternspiral, H_FLOAT, T_FLOAT, "%12.4g", "", "stat.fit.param;src.morph.param", ""}, -/* + {"FLUX_BAR", "Bar total flux from fitting", &outobj2.prof_bar_flux, H_FLOAT, T_FLOAT, "%12.g", "count", "phot.count;stat.fit.param", "ct"}, diff --git a/src/preflist.h b/src/preflist.h index 0a44a5f..64bf318 100644 --- a/src/preflist.h +++ b/src/preflist.h @@ -9,7 +9,7 @@ * * Contents: Keywords for the configuration file. * -* Last modify: 05/02/2010 +* Last modify: 18/05/2010 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -69,9 +69,10 @@ "BACKGROUND", "BACKGROUND_RMS", "MINIBACKGROUND", "MINIBACK_RMS", "-BACKGROUND", "FILTERED", "OBJECTS", "APERTURES", "SEGMENTATION", "ASSOC", - "-OBJECTS", "-PSFS", "PSFS", - "-PC_CONVPROTOS", "PC_CONVPROTOS", "PC_PROTOS", - "MAP_SOM", "-MODELS", "MODELS", "PATTERNS", ""}, + "-OBJECTS", "PSFS", "-PSFS", + "PC_CONVPROTOS", "-PC_CONVPROTOS", "PC_PROTOS", + "MAP_SOM", "MODELS", "-MODELS", "SPHEROIDS", "-SPHEROIDS", + "DISKS", "-DISKS", "PATTERNS", "OTHER", ""}, 0, 17, &prefs.ncheck_type}, {"CLEAN", P_BOOL, &prefs.clean_flag}, {"CLEAN_PARAM", P_FLOAT, &prefs.clean_param, 0,0, 0.1,10.0}, @@ -129,7 +130,7 @@ {"PSF_NMAX", P_INT, &prefs.psf_npsfmax, 1, PSF_NPSFMAX}, {"SATUR_KEY", P_STRING, prefs.satur_key}, {"SATUR_LEVEL", P_FLOAT, &prefs.satur_level, 0,0, -1e+30, 1e+30}, - {"SEEING_FWHM", P_FLOAT, &prefs.seeing_fwhm, 0,0, 1e-10, 1e+10}, + {"SEEING_FWHM", P_FLOAT, &prefs.seeing_fwhm, 0,0, 0.0, 1e+10}, {"SOM_NAME", P_STRING, prefs.som_name}, {"STARNNW_NAME", P_STRING, prefs.nnw_name}, {"THRESH_TYPE", P_KEYLIST, prefs.thresh_type, 0,0, 0.0,0.0, diff --git a/src/prefs.c b/src/prefs.c index 229aa0b..fed2673 100644 --- a/src/prefs.c +++ b/src/prefs.c @@ -9,7 +9,7 @@ * * Contents: Functions to handle the configuration file. * -* Last modify: 11/09/2009 +* Last modify: 08/03/2010 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -588,13 +588,21 @@ void useprefs() prefs.pc_flag = 1; } -/*-------------------------- Profile-fitting -------------------------------*/ +/*----------------------------- Model-fitting -------------------------------*/ if (prefs.check_flag) for (i=0; imodpix, float, PROFIT_MAXMODSIZE*PROFIT_MAXMODSIZE); + QMALLOC16(profit->modpix2, float, PROFIT_MAXMODSIZE*PROFIT_MAXMODSIZE); + QMALLOC16(profit->cmodpix, float, PROFIT_MAXMODSIZE*PROFIT_MAXMODSIZE); + QMALLOC16(profit->psfpix, float, PROFIT_MAXMODSIZE*PROFIT_MAXMODSIZE); + QMALLOC16(profit->objpix, PIXTYPE, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE); + QMALLOC16(profit->objweight, PIXTYPE, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE); + 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->covar, float, profit->nparam*profit->nparam); profit->nprof = nprof; + profit->oversamp = PROFIT_OVERSAMP; profit->fluxfac = 1.0; /* Default */ return profit; @@ -142,7 +159,7 @@ INPUT Prof structure. OUTPUT -. NOTES -. AUTHOR E. Bertin (IAP) -VERSION 26/04/2008 +VERSION 06/04/2010 ***/ void profit_end(profitstruct *profit) { @@ -150,6 +167,15 @@ void profit_end(profitstruct *profit) for (p=0; pnprof; p++) prof_end(profit->prof[p]); + free(profit->modpix); + free(profit->modpix2); + free(profit->cmodpix); + free(profit->psfpix); + free(profit->lmodpix); + free(profit->lmodpix2); + free(profit->objpix); + free(profit->objweight); + free(profit->resi); free(profit->prof); free(profit->covar); free(profit->psfdft); @@ -173,7 +199,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 14/12/2009 +VERSION 08/03/2010 ***/ void profit_fit(profitstruct *profit, picstruct *field, picstruct *wfield, @@ -181,13 +207,15 @@ void profit_fit(profitstruct *profit, { profitstruct pprofit; profitstruct hdprofit; - patternstruct *pattern; + patternstruct *pattern; psfstruct *psf; checkstruct *check; double emx2,emy2,emxy, a , cp,sp, cn, bn, n; - float **list, + float param0[PARAM_NPARAM], param1[PARAM_NPARAM], + param[PARAM_NPARAM], + **list, *cov, - psf_fwhm, dchi2, err, aspect; + psf_fwhm, dchi2, err, aspect, chi2; int *index, i,j,p, nparam, nparam2, ncomp, nprof; @@ -223,6 +251,7 @@ void profit_fit(profitstruct *profit, } else profit->subsamp = 1.0; + profit->nobjpix = profit->objnaxisn[0]*profit->objnaxisn[1]; /* Use (dirty) global variables to interface with lmfit */ the_field = field; @@ -231,9 +260,6 @@ void profit_fit(profitstruct *profit, profit->obj = obj; profit->obj2 = obj2; - QMALLOC16(profit->objpix, PIXTYPE, profit->objnaxisn[0]*profit->objnaxisn[1]); - QMALLOC16(profit->objweight, PIXTYPE,profit->objnaxisn[0]*profit->objnaxisn[1]); - QMALLOC16(profit->lmodpix, PIXTYPE, profit->objnaxisn[0]*profit->objnaxisn[1]); profit->nresi = profit_copyobjpix(profit, field, wfield); /* Check if the number of constraints exceeds the number of free parameters */ if (profit->nresi < nparam) @@ -252,8 +278,6 @@ void profit_fit(profitstruct *profit, return; } - QMALLOC16(profit->resi, float, profit->nresi); - /* Create pixmap at PSF resolution */ profit->modnaxisn[0] = ((int)(profit->objnaxisn[0]*profit->subsamp/profit->pixstep @@ -271,13 +295,8 @@ void profit_fit(profitstruct *profit, profit->modnaxisn[0] = profit->modnaxisn[1] = PROFIT_MAXMODSIZE; obj2->prof_flag |= PROFLAG_MODSUB; } + profit->nmodpix = profit->modnaxisn[0]*profit->modnaxisn[1]; -/* Allocate memory for the complete model */ - QMALLOC16(profit->modpix, float, profit->modnaxisn[0]*profit->modnaxisn[1]); - memset(profit->modpix, 0, profit->modnaxisn[0]*profit->modnaxisn[1]*sizeof(float)); - QMALLOC16(profit->psfpix, float, profit->modnaxisn[0]*profit->modnaxisn[1]); -/* Allocate memory for the partial model */ - QMALLOC16(profit->pmodpix, float, profit->modnaxisn[0]*profit->modnaxisn[1]); /* Compute the local PSF */ profit_psf(profit); @@ -286,42 +305,117 @@ void profit_fit(profitstruct *profit, profit_resetparams(profit); -//the_gal++; - /* Actual minimisation */ fft_reset(); the_gal++; + for (p=0; pnparam; p++) + profit->freeparam_flag[p] = 1; + profit->nfreeparam = profit->nparam; + profit->niter = profit_minimize(profit, PROFIT_MAXITER); - profit_residuals(profit,field,wfield, 10.0, profit->param,profit->resi); +/* + chi2 = profit->chi2; + for (p=0; pparaminit[p]; + profit_resetparams(profit); + for (p=0; pparaminit[p] = param1[p] + (profit->paraminit[p]covar[p*(nparam+1)]); -/* Convert covariance matrix to bound space */ - profit_covarunboundtobound(profit); + profit->niter = profit_minimize(profit, PROFIT_MAXITER); + if (chi2chi2) + for (p=0; pparaminit[p] = param1[p]; + +list = profit->paramlist; +index = profit->paramindex; +for (i=0; ifreeparam_flag[index[i]] = 0; +profit->nfreeparam = 2; +profit->niter = profit_minimize(profit, PROFIT_MAXITER); +*/ for (p=0; pparamerr[p]= sqrt(profit->covar[p*(nparam+1)]); -/* Equate param and paraminit vectors to avoid confusion later on */ - for (p=0; pnparam; p++) - profit->param[p] = profit->paraminit[p]; - /* CHECK-Images */ + 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], + profit->ix,profit->iy, 1.0); + } if ((check = prefs.check[CHECK_SUBPROFILES])) { - profit_residuals(profit,field,wfield, 0.0, profit->param,profit->resi); + profit_residuals(profit,field,wfield, 0.0, profit->paraminit, NULL); addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1], profit->ix,profit->iy, -1.0); } - if ((check = prefs.check[CHECK_PROFILES])) + if ((check = prefs.check[CHECK_SPHEROIDS])) + { +/*-- Set to 0 flux components that do not belong to spheroids */ + for (p=0; pnparam; p++) + param[p] = profit->paraminit[p]; + list = profit->paramlist; + index = profit->paramindex; + for (i=0; ilmodpix, profit->objnaxisn[0],profit->objnaxisn[1], + profit->ix,profit->iy, 1.0); + } + if ((check = prefs.check[CHECK_SUBSPHEROIDS])) + { +/*-- Set to 0 flux components that do not belong to spheroids */ + for (p=0; pnparam; p++) + param[p] = profit->paraminit[p]; + list = profit->paramlist; + index = profit->paramindex; + for (i=0; ilmodpix, profit->objnaxisn[0],profit->objnaxisn[1], + profit->ix,profit->iy, -1.0); + } + if ((check = prefs.check[CHECK_DISKS])) { - profit_residuals(profit,field,wfield, 0.0, profit->param,profit->resi); +/*-- Set to 0 flux components that do not belong to disks */ + for (p=0; pnparam; p++) + param[p] = profit->paraminit[p]; + list = profit->paramlist; + index = profit->paramindex; + for (i=0; ilmodpix, profit->objnaxisn[0],profit->objnaxisn[1], profit->ix,profit->iy, 1.0); } + if ((check = prefs.check[CHECK_SUBDISKS])) + { +/*-- Set to 0 flux components that do not belong to disks */ + for (p=0; pnparam; p++) + param[p] = profit->paraminit[p]; + list = profit->paramlist; + index = profit->paramindex; + for (i=0; ilmodpix, profit->objnaxisn[0],profit->objnaxisn[1], + profit->ix,profit->iy, -1.0); + } + +/* Compute compressed residuals */ + profit_residuals(profit,field,wfield, 10.0, profit->paraminit,profit->resi); /* Fill measurement parameters */ if (FLAG(obj2.prof_vector)) { for (p=0; pprof_vector[p]= profit->param[p]; + obj2->prof_vector[p]= profit->paraminit[p]; } if (FLAG(obj2.prof_errvector)) { @@ -410,47 +504,12 @@ the_gal++; } } -/* Do measurements on the rasterised model (shear and surface brightnesses) */ - if (FLAG(obj2.prof_mx2) || FLAG(obj2.peak_prof)) - { - float scalefac, imsizefac, flux, lost, sum, lostfluxfrac; - -/*-- Allocate "high-definition" rasters only to make measurements */ - hdprofit.modnaxisn[0] = hdprofit.modnaxisn[1] = PROFIT_HIDEFRES; -/*-- Find best image size factor from fitting results */ - imsizefac = 2.0*profit_minradius(profit, PROFIT_REFFFAC)/profit->pixstep - / (float)profit->modnaxisn[0]; - if (imsizefac<0.01) - imsizefac = 0.01; - else if (imsizefac>100.0) - imsizefac = 100.0; - scalefac = (float)hdprofit.modnaxisn[0] / (float)profit->modnaxisn[0] - / imsizefac; - hdprofit.pixstep = profit->pixstep / scalefac; - hdprofit.fluxfac = scalefac*scalefac; - QCALLOC(hdprofit.modpix, float, - hdprofit.modnaxisn[0]*hdprofit.modnaxisn[1]*sizeof(float)); - QCALLOC(hdprofit.pmodpix, float, - hdprofit.modnaxisn[0]*hdprofit.modnaxisn[1]*sizeof(float)); - - lost = sum = 0.0; - for (p=0; pnprof; p++) - { - sum += (flux = prof_add(profit->prof[p], &hdprofit)); - lost += flux*profit->prof[p]->lostfluxfrac; - } - lostfluxfrac = sum > 0.0? lost / sum : 0.0; + if (FLAG(obj2.prof_mx2)) + profit_moments(profit, obj2); - if (FLAG(obj2.prof_mx2)) - profit_moments(&hdprofit, obj2); - - if (FLAG(obj2.peak_prof)) - profit_surface(&hdprofit, obj2, lostfluxfrac); - -/*-- Free rasters */ - free(hdprofit.modpix); - free(hdprofit.pmodpix); - } +/* Do measurements on the rasterised model (surface brightnesses) */ + if (FLAG(obj2.peak_prof)) + profit_surface(profit, obj2); /* Spheroid */ if (FLAG(obj2.prof_spheroid_flux)) @@ -556,7 +615,7 @@ the_gal++; if (prefs.pattern_flag) { profit_residuals(profit,field,wfield, PROFIT_DYNPARAM, - profit->param,profit->resi); + profit->paraminit,profit->resi); pattern = pattern_init(profit, prefs.pattern_type, prefs.prof_disk_patternncomp); pattern_fit(pattern, profit); @@ -657,8 +716,12 @@ the_gal++; pprofit.paraminit[pprofit.paramindex[PARAM_Y]] = *profit->paramlist[PARAM_Y]; } pprofit.paraminit[pprofit.paramindex[PARAM_DISK_FLUX]] = profit->flux; + for (p=0; pchi2); @@ -684,13 +747,6 @@ the_gal++; /* clean up. */ fft_reset(); - free(profit->modpix); - free(profit->psfpix); - free(profit->pmodpix); - free(profit->lmodpix); - free(profit->objpix); - free(profit->objweight); - free(profit->resi); return; } @@ -832,7 +888,7 @@ INPUT Profile-fitting structure. OUTPUT -. NOTES -. AUTHOR E. Bertin (IAP) -VERSION 05/10/2009 +VERSION 07/07/2010 ***/ void profit_psf(profitstruct *profit) { @@ -876,13 +932,13 @@ void profit_psf(profitstruct *profit) /* Normalize PSF flux (just in case...) */ flux *= profit->pixstep*profit->pixstep; - if (fabs(flux) > 0.0) - { - norm = 1.0/flux; - pixout = profit->psfpix; - for (i=profit->modnaxisn[0]*profit->modnaxisn[1]; i--;) - *(pixout++) *= norm; - } + if (fabs(flux) <= 0.0) + error(EXIT_FAILURE, "*Error*: PSF model is empty or negative: ", psf->name); + + norm = 1.0/flux; + pixout = profit->psfpix; + for (i=profit->modnaxisn[0]*profit->modnaxisn[1]; i--;) + *(pixout++) *= norm; return; } @@ -896,33 +952,33 @@ INPUT Pointer to the profit structure involved in the fit, OUTPUT Number of iterations used. NOTES -. AUTHOR E. Bertin (IAP) -VERSION 23/05/2008 +VERSION 02/04/2010 ***/ int profit_minimize(profitstruct *profit, int niter) { - float lm_opts[5], info[LM_INFO_SZ]; - int m,n; - -/* Allocate work space */ - n = profit->nparam; - m = profit->nresi; + double lm_opts[5], info[LM_INFO_SZ]; + double dcovar[PARAM_NPARAM*PARAM_NPARAM], + dparam[PARAM_NPARAM]; - memset(profit->resi, 0, profit->nresi*sizeof(float)); - memset(profit->covar, 0, profit->nparam*profit->nparam*sizeof(float)); - profit_boundtounbound(profit, profit->paraminit); + profit->iter = 0; + memset(dcovar, 0, profit->nparam*profit->nparam*sizeof(double)); /* Perform fit */ lm_opts[0] = 1.0e-3; - lm_opts[1] = 1.0e-17; - lm_opts[2] = 1.0e-17; - lm_opts[3] = 1.0e-17; - lm_opts[4] = 1.0e-6; + lm_opts[1] = 1.0e-12; + lm_opts[2] = 1.0e-12; + lm_opts[3] = 1.0e-12; + lm_opts[4] = 1.0e-3; - niter = slevmar_dif(profit_evaluate, profit->paraminit, profit->resi, - n, m, niter, lm_opts, info, NULL, profit->covar, profit); + profit_boundtounbound(profit, profit->paraminit, dparam, PARAM_ALLPARAMS); - profit_unboundtobound(profit, profit->paraminit); + niter = dlevmar_dif(profit_evaluate, dparam, NULL, profit->nfreeparam, + profit->nresi, niter, lm_opts, info, NULL, dcovar, profit); + profit_unboundtobound(profit, dparam, profit->paraminit, PARAM_ALLPARAMS); + +/* Convert covariance matrix to bounded space */ + profit_covarunboundtobound(profit, dcovar, profit->covar); return niter; } @@ -976,29 +1032,182 @@ void profit_printout(int n_par, float* par, int m_dat, float* fvec, /****** profit_evaluate ****************************************************** -PROTO void profit_evaluate(float *par, float *fvec, int m, int n, - void *adata) +PROTO void profit_evaluate(double *par, double *fvec, int m, int n, + void *adata) PURPOSE Provide a function returning residuals to levmar. INPUT Pointer to the vector of parameters, pointer to the vector of residuals (output), number of model parameters, number of data points, - pointer to a data structure (unused). + pointer to a data structure (we use it for the profit structure here). OUTPUT -. -NOTES Input arguments are there only for compatibility purposes (unused) +NOTES -. AUTHOR E. Bertin (IAP) -VERSION 18/09/2008 +VERSION 08/07/2010 ***/ -void profit_evaluate(float *par, float *fvec, int m, int n, - void *adata) +void profit_evaluate(double *dpar, double *fvec, int m, int n, void *adata) { - profitstruct *profit; + profitstruct *profit; + profstruct **prof; + double *dpar0, *dresi; + float *modpixt, *profpixt, *resi, + tparam, val; + PIXTYPE *lmodpixt,*lmodpix2t, *objpix,*weight, + wval; + int c,f,i,p,q, fd,pd, jflag,sflag, nprof; profit = (profitstruct *)adata; - profit_unboundtobound(profit, par); - profit_residuals(profit, the_field, the_wfield, PROFIT_DYNPARAM, par, fvec); - profit_boundtounbound(profit, par); - profit_printout(m, par, n, fvec, adata, 0, -1, 0 ); + +/* Detect "Jacobian-related" calls */ + jflag = pd = fd = 0; + dpar0 = profit->dparam; + if (profit->iter) + { + f = q = 0; + for (p=0; pnparam; p++) + { + if (dpar[f] - dpar0[f] != 0.0) + { + pd = p; + fd = f; + q++; + } + if (profit->freeparam_flag[p]) + f++; + } + if (f>0 && q==1) + jflag = 1; + } + + if (jflag) + { + prof = profit->prof; + nprof = profit->nprof; + +/*-- "Jacobian call" */ + tparam = profit->param[pd]; + profit_unboundtobound(profit, &dpar[fd], &profit->param[pd], pd); + sflag = 1; + switch(profit->paramrevindex[pd]) + { + case PARAM_X: + case PARAM_Y: + profit_resample(profit, profit->cmodpix, profit->lmodpix2, 1.0); + lmodpixt = profit->lmodpix; + lmodpix2t = profit->lmodpix2; + for (i=profit->nobjpix;i--;) + *(lmodpix2t++) -= *(lmodpixt++); + break; + case PARAM_SPHEROID_FLUX: + case PARAM_DISK_FLUX: + case PARAM_ARMS_FLUX: + case PARAM_BAR_FLUX: + if (nprof==1 && tparam != 0.0) + { + lmodpixt = profit->lmodpix; + lmodpix2t = profit->lmodpix2; + val = (profit->param[pd] - tparam) / tparam; + for (i=profit->nobjpix;i--;) + *(lmodpix2t++) = val**(lmodpixt++); + } + else + { + for (c=0; cflux == &profit->param[pd]) + break; + memcpy(profit->modpix, prof[c]->pix, profit->nmodpix*sizeof(float)); + profit_convolve(profit, profit->modpix); + profit_resample(profit, profit->modpix, profit->lmodpix2, + profit->param[pd] - tparam); + } + break; + case PARAM_SPHEROID_REFF: + case PARAM_SPHEROID_ASPECT: + case PARAM_SPHEROID_POSANG: + case PARAM_SPHEROID_SERSICN: + sflag = 0; /* We are in the same switch */ + for (c=0; ccode == PROF_SERSIC + || prof[c]->code == PROF_DEVAUCOULEURS) + break; + case PARAM_DISK_SCALE: + case PARAM_DISK_ASPECT: + case PARAM_DISK_POSANG: + if (sflag) + for (c=0; ccode == PROF_EXPONENTIAL) + break; + sflag = 0; + case PARAM_ARMS_QUADFRAC: + case PARAM_ARMS_SCALE: + case PARAM_ARMS_START: + case PARAM_ARMS_POSANG: + case PARAM_ARMS_PITCH: + case PARAM_ARMS_PITCHVAR: + case PARAM_ARMS_WIDTH: + if (sflag) + for (c=0; ccode == PROF_ARMS) + break; + sflag = 0; + case PARAM_BAR_ASPECT: + case PARAM_BAR_POSANG: + if (sflag) + for (c=0; ccode == PROF_ARMS) + break; + modpixt = profit->modpix; + profpixt = prof[c]->pix; + val = -*prof[c]->flux; + for (i=profit->nmodpix;i--;) + *(modpixt++) = val**(profpixt++); + memcpy(profit->modpix2, prof[c]->pix, profit->nmodpix*sizeof(float)); + prof_add(profit, prof[c], 0); + memcpy(prof[c]->pix, profit->modpix2, profit->nmodpix*sizeof(float)); + profit_convolve(profit, profit->modpix); + profit_resample(profit, profit->modpix, profit->lmodpix2, 1.0); + break; + default: + error(EXIT_FAILURE, "*Internal Error*: ", + "unknown parameter index in profit_jacobian()"); + break; + } + objpix = profit->objpix; + weight = profit->objweight; + lmodpixt = profit->lmodpix; + lmodpix2t = profit->lmodpix2; + resi = profit->resi; + dresi = fvec; + if (PROFIT_DYNPARAM > 0.0) + for (i=profit->nobjpix;i--; lmodpixt++, lmodpix2t++) + { + val = *(objpix++); + if ((wval=*(weight++))>0.0) + *(dresi++) = *(resi++) + *lmodpix2t + * wval/(1.0+wval*fabs(*lmodpixt - val)/PROFIT_DYNPARAM); + } + else + for (i=profit->nobjpix;i--; lmodpix2t++) + if ((wval=*(weight++))>0.0) + *(dresi++) = *(resi++) + *lmodpix2t * wval; + } + else + { +/*-- "Regular call" */ + for (p=0; pnparam; p++) + dpar0[p] = dpar[p]; + profit_unboundtobound(profit, dpar, profit->param, PARAM_ALLPARAMS); + + profit_residuals(profit, the_field, the_wfield, PROFIT_DYNPARAM, + profit->param, profit->resi); + + for (p=0; pnresi; p++) + fvec[p] = profit->resi[p]; + } + +// profit_printout(m, par, n, fvec, adata, 0, -1, 0 ); + profit->iter++; + return; } @@ -1017,15 +1226,15 @@ INPUT Profile-fitting structure, OUTPUT Vector of residuals. NOTES -. AUTHOR E. Bertin (IAP) -VERSION 24/03/2009 +VERSION 08/07/2010 ***/ float *profit_residuals(profitstruct *profit, picstruct *field, picstruct *wfield, float dynparam, float *param, float *resi) { - int p; + int p, nmodpix; - memset(profit->modpix, 0, - profit->modnaxisn[0]*profit->modnaxisn[1]*sizeof(float)); + nmodpix = profit->modnaxisn[0]*profit->modnaxisn[1]*sizeof(float); + memset(profit->modpix, 0, nmodpix); for (p=0; pnparam; p++) profit->param[p] = param[p]; /* Simple PSF shortcut */ @@ -1039,11 +1248,14 @@ float *profit_residuals(profitstruct *profit, picstruct *field, { profit->flux = 0.0; for (p=0; pnprof; p++) - profit->flux += prof_add(profit->prof[p], profit); - profit_convolve(profit, profit->modpix); - profit_resample(profit, profit->modpix, profit->lmodpix, 1.0); + profit->flux += prof_add(profit, profit->prof[p], 0); + memcpy(profit->cmodpix, profit->modpix, profit->nmodpix*sizeof(float)); + profit_convolve(profit, profit->cmodpix); + profit_resample(profit, profit->cmodpix, profit->lmodpix, 1.0); } - profit_compresi(profit, dynparam, resi); + + if (resi) + profit_compresi(profit, dynparam, resi); return resi; } @@ -1085,7 +1297,7 @@ float *profit_compresi(profitstruct *profit, float dynparam, float *resi) val = *(objpix++); if ((wval=*(objweight++))>0.0) { - val2 = (val - *lmodpix)*wval*invsig; + val2 = (*lmodpix - val)*wval*invsig; val2 = val2>0.0? logf(1.0+val2) : -logf(1.0-val2); *(resit++) = val2*dynparam; error += val2*val2; @@ -1100,7 +1312,7 @@ float *profit_compresi(profitstruct *profit, float dynparam, float *resi) val = *(objpix++); if ((wval=*(objweight++))>0.0) { - val2 = (val - *lmodpix)*wval; + val2 = (*lmodpix - val)*wval; *(resit++) = val2; error += val2*val2; } @@ -1113,82 +1325,206 @@ float *profit_compresi(profitstruct *profit, float dynparam, float *resi) /****** profit_resample ****************************************************** -PROTO void prof_resample(profitstruct *profit, float *inpix, - PIXTYPE *outpix) +PROTO int prof_resample(profitstruct *profit, float *inpix, + PIXTYPE *outpix, float factor) PURPOSE Resample the current full resolution model to image resolution. -INPUT Profile-fitting structure. -OUTPUT -. +INPUT Profile-fitting structure, + pointer to input raster, + pointer to output raster, + multiplicating factor. +OUTPUT RETURN_ERROR if the rasters don't overlap, RETURN_OK otherwise. NOTES -. AUTHOR E. Bertin (IAP) -VERSION 09/10/2009 +VERSION 01/05/2010 ***/ -void profit_resample(profitstruct *profit, float *inpix, PIXTYPE *outpix, +int profit_resample(profitstruct *profit, float *inpix, PIXTYPE *outpix, float factor) { - double flux; - float posin[2], posout[2], dnaxisn[2], - *dx,*dy, - xcout,ycout, xcin,ycin, invpixstep, pix, off,step; - int d,i, ix,iy,ns; + PIXTYPE *pixout,*pixout0; + float *pixin,*pixin0, *mask,*maskt, *pixinout, *dpixin,*dpixin0, + *dpixout,*dpixout0, *dx,*dy, + xcin,xcout,ycin,ycout, xsin,ysin, xin,yin, x,y, dxm,dym, val, + invpixstep; + int *start,*startt, *nmask,*nmaskt, + i,j,k,n,t, + ixsout,iysout, ixout,iyout, dixout,diyout, nxout,nyout, + iysina, nyin, hmw,hmh, ix,iy, ixin,iyin; - xcout = (float)(profit->objnaxisn[0]/2) + 1.0; /* FITS convention */ + invpixstep = profit->subsamp/profit->pixstep; + + xcin = (profit->modnaxisn[0]/2); + xcout = (float)(profit->objnaxisn[0]/2); if ((dx=(profit->paramlist[PARAM_X]))) xcout += *dx; - ycout = (float)(profit->objnaxisn[1]/2) + 1.0; /* FITS convention */ + + xsin = xcin - xcout*invpixstep; /* Input start x-coord*/ + if ((int)xsin >= profit->modnaxisn[0]) + return RETURN_ERROR; + ixsout = 0; /* Int. part of output start x-coord */ + if (xsin<0.0) + { + dixout = (int)(1.0-xsin/invpixstep); +/*-- Simply leave here if the images do not overlap in x */ + if (dixout >= profit->objnaxisn[0]) + return RETURN_ERROR; + ixsout += dixout; + xsin += dixout*invpixstep; + } + nxout = (int)((profit->modnaxisn[0]-xsin)/invpixstep);/* nb of interpolated + input pixels along x */ + if (nxout>(ixout=profit->objnaxisn[0]-ixsout)) + nxout = ixout; + if (!nxout) + return RETURN_ERROR; + + ycin = (profit->modnaxisn[1]/2); + ycout = (float)(profit->objnaxisn[1]/2); if ((dy=(profit->paramlist[PARAM_Y]))) ycout += *dy; - xcin = (profit->modnaxisn[0]/2) + 1.0; /* FITS convention */ - ycin = (profit->modnaxisn[1]/2) + 1.0; /* FITS convention */ - invpixstep = profit->subsamp/profit->pixstep; -/* Initialize multi-dimensional counters */ - for (d=0; d<2; d++) + ysin = ycin - ycout*invpixstep; /* Input start y-coord*/ + if ((int)ysin >= profit->modnaxisn[1]) + return RETURN_ERROR; + iysout = 0; /* Int. part of output start y-coord */ + if (ysin<0.0) + { + diyout = (int)(1.0-ysin/invpixstep); +/*-- Simply leave here if the images do not overlap in y */ + if (diyout >= profit->objnaxisn[1]) + return RETURN_ERROR; + iysout += diyout; + ysin += diyout*invpixstep; + } + nyout = (int)((profit->modnaxisn[1]-ysin)/invpixstep);/* nb of interpolated + input pixels along y */ + if (nyout>(iyout=profit->objnaxisn[1]-iysout)) + nyout = iyout; + if (!nyout) + return RETURN_ERROR; + +/* Set the yrange for the x-resampling with some margin for interpolation */ + iysina = (int)ysin; /* Int. part of Input start y-coord with margin */ + hmh = INTERPH/2 - 1; /* Interpolant start */ + if (iysina<0 || ((iysina -= hmh)< 0)) + iysina = 0; + nyin = (int)(ysin+nyout*invpixstep)+INTERPH-hmh;/* Interpolated Input y size*/ + if (nyin>profit->modnaxisn[1]) /* with margin */ + nyin = profit->modnaxisn[1]; +/* Express everything relative to the effective Input start (with margin) */ + nyin -= iysina; + ysin -= (float)iysina; + +/* Allocate interpolant stuff for the x direction */ + QMALLOC(mask, float, nxout*INTERPW); /* Interpolation masks */ + QMALLOC(nmask, int, nxout); /* Interpolation mask sizes */ + QMALLOC(start, int, nxout); /* Int. part of Input conv starts */ +/* Compute the local interpolant and data starting points in x */ + hmw = INTERPW/2 - 1; + xin = xsin; + maskt = mask; + nmaskt = nmask; + startt = start; + for (j=nxout; j--; xin+=invpixstep) { - posout[d] = 1.0; /* FITS convention */ - dnaxisn[d] = profit->objnaxisn[d]+0.5; + ix = (ixin=(int)xin) - hmw; + dxm = ixin - xin - hmw; /* starting point in the interpolation func */ + if (ix < 0) + { + n = INTERPW+ix; + dxm -= (float)ix; + ix = 0; + } + else + n = INTERPW; + if (n>(t=profit->modnaxisn[0]-ix)) + n=t; + *(startt++) = ix; + *(nmaskt++) = n; + for (x=dxm, i=n; i--; x+=1.0) + *(maskt++) = INTERPF(x); + } + + QCALLOC(pixinout, float, nxout*nyin); /* Intermediary frame-buffer */ + +/* Make the interpolation in x (this includes transposition) */ + pixin0 = inpix + iysina*profit->modnaxisn[0]; + dpixout0 = pixinout; + for (k=nyin; k--; pixin0+=profit->modnaxisn[0], dpixout0++) + { + maskt = mask; + nmaskt = nmask; + startt = start; + dpixout = dpixout0; + for (j=nxout; j--; dpixout+=nyin) + { + pixin = pixin0+*(startt++); + val = 0.0; + for (i=*(nmaskt++); i--;) + val += *(maskt++)**(pixin++); + *dpixout = val; + } } -/* Remap each pixel (and rebin if necessary) */ - flux = 0.0; - ns=(int)profit->subsamp; - if (ns>1) - { - step = 1.0/profit->subsamp; - off = 0.5*(step - 1.0); - xcin += off; - ycin += off; - for (i=profit->objnaxisn[0]*profit->objnaxisn[1]; i--;) +/* Reallocate interpolant stuff for the y direction */ + QREALLOC(mask, float, nyout*INTERPH); /* Interpolation masks */ + QREALLOC(nmask, int, nyout); /* Interpolation mask sizes */ + QREALLOC(start, int, nyout); /* Int. part of Input conv starts */ + +/* Compute the local interpolant and data starting points in y */ + hmh = INTERPH/2 - 1; + yin = ysin; + maskt = mask; + nmaskt = nmask; + startt = start; + for (j=nyout; j--; yin+=invpixstep) + { + iy = (iyin=(int)yin) - hmh; + dym = iyin - yin - hmh; /* starting point in the interpolation func */ + if (iy < 0) { - posin[0] = (posout[0] - xcout)*invpixstep + xcin; - posin[1] = (posout[1] - ycout)*invpixstep + ycin; - pix = 0.0; - for (iy=ns; iy--; posin[1] += step) - for (ix=ns; ix--; posin[0] += step) - pix += (PIXTYPE)(factor*interpolate_pix(posin, inpix, - profit->modnaxisn, INTERP_LANCZOS3)); - flux += (*(outpix++) = pix); - for (d=0; d<2; d++) - if ((posout[d]+=1.0) < dnaxisn[d]) - break; - else - posout[d] = 1.0; + n = INTERPH+iy; + dym -= (float)iy; + iy = 0; } - } - else - for (i=profit->objnaxisn[0]*profit->objnaxisn[1]; i--;) + else + n = INTERPH; + if (n>(t=nyin-iy)) + n=t; + *(startt++) = iy; + *(nmaskt++) = n; + for (y=dym, i=n; i--; y+=1.0) + *(maskt++) = INTERPF(y); + } + +/* Initialize destination buffer to zero */ + memset(outpix, 0, (size_t)profit->nobjpix*sizeof(PIXTYPE)); + +/* Make the interpolation in y and transpose once again */ + dpixin0 = pixinout; + pixout0 = outpix+ixsout+iysout*profit->objnaxisn[0]; + for (k=nxout; k--; dpixin0+=nyin, pixout0++) + { + maskt = mask; + nmaskt = nmask; + startt = start; + pixout = pixout0; + for (j=nyout; j--; pixout+=profit->objnaxisn[0]) { - posin[0] = (posout[0] - xcout)*invpixstep + xcin; - posin[1] = (posout[1] - ycout)*invpixstep + ycin; - flux += ((*(outpix++) = (PIXTYPE)(factor*interpolate_pix(posin, inpix, - profit->modnaxisn, INTERP_LANCZOS3)))); - for (d=0; d<2; d++) - if ((posout[d]+=1.0) < dnaxisn[d]) - break; - else - posout[d] = 1.0; + dpixin = dpixin0+*(startt++); + val = 0.0; + for (i=*(nmaskt++); i--;) + val += *(maskt++)**(dpixin++); + *pixout = (PIXTYPE)(factor*val); } + } - return; +/* Free memory */ + free(pixinout); + free(mask); + free(nmask); + free(start); + + return RETURN_OK; } @@ -1691,47 +2027,60 @@ INPUT Profile-fitting structure, OUTPUT -. NOTES -. AUTHOR E. Bertin (IAP) -VERSION 17/09/2009 +VERSION 07/07/2010 ***/ void profit_moments(profitstruct *profit, obj2struct *obj2) { - double mx,my, sum, mx2,my2,mxy, den, temp,temp2,pmx2,theta; - float *pix, - hw,hh, x,y, xstart, val, r2max; - int ix,iy; - - hw = (float)(profit->modnaxisn[0]/2); - hh = (float)(profit->modnaxisn[1]/2); - r2max = hwmodpix; - mx2 = my2 = mxy = mx = my = sum = 0.0; - for (iy=profit->modnaxisn[1]; iy--; y+=1.0) + profstruct *prof; + double m0, mx2,my2,mxy, den, temp,temp2,pmx2,theta; + int p; + +/* hw = (float)(profit->modnaxisn[0]/2);*/ +/* hh = (float)(profit->modnaxisn[1]/2);*/ +/* r2max = hwmodpix;*/ +/* mx2 = my2 = mxy = mx = my = sum = 0.0;*/ +/* for (iy=profit->modnaxisn[1]; iy--; y+=1.0)*/ +/* {*/ +/* x = xstart;*/ +/* for (ix=profit->modnaxisn[0]; ix--; x+=1.0)*/ +/* if (y*y+x*x <= r2max)*/ +/* {*/ +/* val = *(pix++);*/ +/* sum += val;*/ +/* mx += val*x;*/ +/* my += val*y;*/ +/* mx2 += val*x*x;*/ +/* mxy += val*x*y;*/ +/* my2 += val*y*y;*/ +/* }*/ +/* else*/ +/* pix++;*/ +/* }*/ + +/* if (sum <= 1.0/BIG)*/ +/* sum = 1.0;*/ +/* mx /= sum;*/ +/* my /= sum;*/ +/* obj2->prof_mx2 = mx2 = mx2/sum - mx*mx;*/ +/* obj2->prof_my2 = my2 = my2/sum - my*my;*/ +/* obj2->prof_mxy = mxy = mxy/sum - mx*my;*/ + + m0 = mx2 = my2 = mxy = 0.0; + for (p=0; pnprof; p++) { - x = xstart; - for (ix=profit->modnaxisn[0]; ix--; x+=1.0) - if (y*y+x*x <= r2max) - { - val = *(pix++); - sum += val; - mx += val*x; - my += val*y; - mx2 += val*x*x; - mxy += val*x*y; - my2 += val*y*y; - } - else - pix++; + prof = profit->prof[p]; + prof_moments(profit, prof); + m0 += *prof->flux; + mx2 += prof->mx2**prof->flux; + my2 += prof->my2**prof->flux; + mxy += prof->mxy**prof->flux; } - - if (sum <= 1.0/BIG) - sum = 1.0; - mx /= sum; - my /= sum; - obj2->prof_mx2 = mx2 = mx2/sum - mx*mx; - obj2->prof_my2 = my2 = my2/sum - my*my; - obj2->prof_mxy = mxy = mxy/sum - mx*my; + obj2->prof_mx2 = mx2 / m0; + obj2->prof_my2 = my2 / m0; + obj2->prof_mxy = mxy / m0; /* Handle fully correlated profiles (which cause a singularity...) */ if ((temp2=mx2*my2-mxy*mxy)<0.00694) @@ -1779,39 +2128,112 @@ void profit_moments(profitstruct *profit, obj2struct *obj2) obj2->prof_theta = theta*180.0/PI; } - return; } /****** profit_surface **************************************************** -PROTO void profit_surface(profitstruct *profit, obj2struct *obj2, - double lostfluxfrac) +PROTO void profit_surface(profitstruct *profit, obj2struct *obj2) PURPOSE Compute surface brightnesses from the unconvolved object model. INPUT Pointer to the profile-fitting structure, - Pointer to obj2 structure, - Fraction of total flux lost in model. + Pointer to obj2 structure. OUTPUT -. NOTES -. AUTHOR E. Bertin (IAP) -VERSION 21/09/2009 +VERSION 08/07/2010 ***/ -void profit_surface(profitstruct *profit, obj2struct *obj2, - double lostfluxfrac) +void profit_surface(profitstruct *profit, obj2struct *obj2) { + profitstruct hdprofit; double dsum,dhsum,dsumoff, dhval, frac, seff; - float *spix, *spixt; - int i, npix, neff; + float *spix, *spixt, + val,vmax, + scalefac, imsizefac, flux, lost, sum, lostfluxfrac; + int i,p, imax, npix, neff; + +/* Allocate "high-definition" raster only to make measurements */ + hdprofit.oversamp = PROFIT_OVERSAMP; + hdprofit.modnaxisn[0] = hdprofit.modnaxisn[1] = PROFIT_HIDEFRES; + npix = hdprofit.nmodpix = hdprofit.modnaxisn[0]*hdprofit.modnaxisn[1]; +/* Find best image size factor from fitting results */ + imsizefac = 2.0*profit_minradius(profit, PROFIT_REFFFAC)/profit->pixstep + / (float)profit->modnaxisn[0]; + if (imsizefac<0.01) + imsizefac = 0.01; + else if (imsizefac>100.0) + imsizefac = 100.0; + scalefac = (float)hdprofit.modnaxisn[0] / (float)profit->modnaxisn[0] + / imsizefac; + hdprofit.pixstep = profit->pixstep / scalefac; + hdprofit.fluxfac = scalefac*scalefac; + QCALLOC(hdprofit.modpix, float,npix*sizeof(float)); + + for (p=0; pnparam; p++) + profit->param[p] = profit->paraminit[p]; + lost = sum = 0.0; + + for (p=0; pnprof; p++) + { + sum += (flux = prof_add(&hdprofit, profit->prof[p],0)); + lost += flux*profit->prof[p]->lostfluxfrac; + } + lostfluxfrac = sum > 0.0? lost / sum : 0.0; - npix = profit->modnaxisn[0]*profit->modnaxisn[1]; -/* Sort model pixel values */ - spix = NULL; /* to avoid gcc -Wall warnings */ - QMEMCPY(profit->modpix, spix, float, npix); - hmedian(spix, npix); - obj2->peak_prof = spix[npix-1]; +/* +char filename[256]; +sprintf(filename, "raster_%02d.fits", the_gal); +check=initcheck(filename, CHECK_OTHER, 0); +check->width = hdprofit.modnaxisn[0]; +check->height = hdprofit.modnaxisn[1]; +reinitcheck(the_field, check); +memcpy(check->pix,hdprofit.modpix,check->npix*sizeof(float)); + + +int r,t; +double ratio,ratio0,ang,ang0, x,x0,y,y0; +list = profit->paramlist; +index = profit->paramindex; +for (p=0; pparaminit[p]; + +ratio0 = profit->paraminit[index[PARAM_SPHEROID_ASPECT]]; +ang0 = profit->paraminit[index[PARAM_SPHEROID_POSANG]]; +x0 = profit->paraminit[index[PARAM_X]]; +y0 = profit->paraminit[index[PARAM_Y]]; +for (r=0;rheight;r++) +for (t=0; twidth;t++) +{ +//x = (r-10.0)/100.0 + x0; +//y = (t-10.0)/100.0 + y0; +ratio = ratio0*exp((r-10.0)/400.0); +ang = ang0+(t-10.0)/3.0; + +for (i=0; iparaminit[index[i]]*sqrt(ratio0/ratio); +} +profit_residuals(profit,field,wfield, PROFIT_DYNPARAM, param,profit->resi); +*((float *)check->pix + t + r*check->width) = profit->chi2; +} +reendcheck(the_field, check); +endcheck(check); +*/ if (FLAG(obj2.fluxeff_prof)) { +/*-- Sort model pixel values */ + spix = NULL; /* to avoid gcc -Wall warnings */ + QMEMCPY(hdprofit.modpix, spix, float, npix); + hmedian(spix, npix); /*-- Build a cumulative distribution */ dsum = 0.0; spixt = spix; @@ -1848,9 +2270,34 @@ void profit_surface(profitstruct *profit, obj2struct *obj2, dsum += (double)*(spixt++); obj2->fluxmean_prof = seff > 0.0? dsum / seff : 0.0; } + free(spix); + } + +/* Compute model peak (overwrites oversampled model!!) */ + if (FLAG(obj2.peak_prof)) + { +/*-- Find position of maximum pixel in current hi-def raster */ + imax = 0; + vmax = -BIG; + spixt = hdprofit.modpix; + for (i=npix; i--;) + if ((val=*(spixt++))>vmax) + { + vmax = val; + imax = i; + } + imax = npix-1 - imax; +/*-- Recompute hi-def model raster without oversampling */ +/*-- and with the same flux correction factor */ + hdprofit.oversamp = 0; + memset(hdprofit.modpix,0, npix*sizeof(float)); + for (p=0; pnprof; p++) + prof_add(&hdprofit, profit->prof[p], 1); + obj2->peak_prof = hdprofit.modpix[imax]; } - free(spix); +/* Free hi-def model raster */ + free(hdprofit.modpix); return; } @@ -1866,7 +2313,7 @@ INPUT Pointer to the profit structure, OUTPUT -. NOTES -. AUTHOR E. Bertin (IAP) -VERSION 29/03/2007 +VERSION 29/03/2010 ***/ void profit_addparam(profitstruct *profit, paramenum paramindex, float **param) @@ -1879,7 +2326,8 @@ void profit_addparam(profitstruct *profit, paramenum paramindex, /*-- No */ { *param = profit->paramlist[paramindex] = &profit->param[profit->nparam]; - profit->paramindex[paramindex] = profit->nparam++; + profit->paramindex[paramindex] = profit->nparam; + profit->paramrevindex[profit->nparam++] = paramindex; } return; @@ -1894,7 +2342,7 @@ INPUT Pointer to the profit structure, OUTPUT -. NOTES -. AUTHOR E. Bertin (IAP) -VERSION 15/12/2009 +VERSION 02/07/2010 ***/ void profit_resetparam(profitstruct *profit, paramenum paramtype) { @@ -1915,7 +2363,7 @@ void profit_resetparam(profitstruct *profit, paramenum paramtype) break; case PARAM_X: param = obj->mx - (int)(obj->mx+0.49999); - range = fabs(obj2->hl_radius*4.0); + range = obj2->hl_radius*4.0; if (range>profit->objnaxisn[0]*2.0) range = profit->objnaxisn[0]*2.0; parammin = -range; @@ -1923,7 +2371,7 @@ void profit_resetparam(profitstruct *profit, paramenum paramtype) break; case PARAM_Y: param = obj->my - (int)(obj->my+0.49999); - range = fabs(obj2->hl_radius*4); + range = obj2->hl_radius*4.0; if (range>profit->objnaxisn[1]*2) range = profit->objnaxisn[1]*2; parammin = -range; @@ -1935,20 +2383,19 @@ void profit_resetparam(profitstruct *profit, paramenum paramtype) parammax = 2*obj2->flux_auto; break; case PARAM_SPHEROID_REFF: - param = fabs(obj2->hl_radius); + param = obj2->hl_radius; parammin = 0.1; parammax = param * 4.0; break; case PARAM_SPHEROID_ASPECT: param = FLAG(obj2.prof_disk_flux)? 1.0 : obj->b/obj->a; parammin = FLAG(obj2.prof_disk_flux)? 0.5 : 0.01; - parammax = 1.0; -parammax = FLAG(obj2.prof_disk_flux)? 2.0 : 100.0; + parammax = FLAG(obj2.prof_disk_flux)? 2.0 : 100.0; break; case PARAM_SPHEROID_POSANG: param = obj->theta; - parammin = 0.0; - parammax = 0.0; + parammin = 90.0; + parammax = 90.0; break; case PARAM_SPHEROID_SERSICN: param = 4.0; @@ -1961,7 +2408,7 @@ parammax = FLAG(obj2.prof_disk_flux)? 2.0 : 100.0; parammax = 2*obj2->flux_auto; break; case PARAM_DISK_SCALE: /* From scalelength to Re */ - param = fabs(obj2->hl_radius)/1.67835*sqrt(obj->a/obj->b); + param = obj2->hl_radius/1.67835*sqrt(obj->a/obj->b); parammin = param/4.0; parammax = param * 4.0; break; @@ -1972,8 +2419,8 @@ parammax = FLAG(obj2.prof_disk_flux)? 2.0 : 100.0; break; case PARAM_DISK_POSANG: param = obj->theta; - parammin = 0.0; - parammax = 0.0; + parammin = 90.0; + parammax = 90.0; break; case PARAM_ARMS_FLUX: param = obj2->flux_auto/2.0; @@ -2137,91 +2584,166 @@ int profit_setparam(profitstruct *profit, paramenum paramtype, /****** profit_boundtounbound ************************************************* -PROTO void profit_boundtounbound(profitstruct *profit, float *param) +PROTO void profit_boundtounbound(profitstruct *profit, + float *param, double *dparam, int index) PURPOSE Convert parameters from bounded to unbounded space. -INPUT Pointer to the profit structure. +INPUT Pointer to the profit structure, + input array of single precision parameters, + output (incomplete) array of double precision parameters, + parameter selection index (<0 = all parameters) OUTPUT -. NOTES -. AUTHOR E. Bertin (IAP) -VERSION 07/09/2009 +VERSION 29/03/2010 ***/ -void profit_boundtounbound(profitstruct *profit, float *param) +void profit_boundtounbound(profitstruct *profit, + float *param, double *dparam, int index) { double num,den; - int p; + float tparam; + int f,p, pstart,np; - for (p=0; pnparam; p++) - if (profit->parammin[p]!=profit->parammax[p]) + if (index<0) + { + pstart = 0; + np = profit->nparam; + } + else + { + pstart = index; + np = pstart+1; + } + + f = 0; + for (p=pstart ; pfreeparam_flag[p]) { - num = param[p] - profit->parammin[p]; - den = profit->parammax[p] - param[p]; - param[p] = num>1e-50? (den>1e-50? log(num/den): 50.0) : -50.0; + tparam = param[p-pstart]; + if (profit->parammin[p]!=profit->parammax[p]) + { + num = tparam - profit->parammin[p]; + den = profit->parammax[p] - tparam; + dparam[f] = num>1e-50? (den>1e-50? log(num/den): 50.0) : -50.0; + } + else if (profit->parammax[p] > 0.0 || profit->parammax[p] < 0.0) + dparam[f] = param[p-pstart] / profit->parammax[p]; + + f++; } return; - } /****** profit_unboundtobound ************************************************* -PROTO void profit_unboundtobound(profitstruct *profit, float *param) +PROTO void profit_unboundtobound(profitstruct *profit, + double *dparam, float *param, int index) PURPOSE Convert parameters from unbounded to bounded space. -INPUT Pointer to the profit structure. +INPUT Pointer to the profit structure, + input (incomplete) array of double precision parameters, + output array of single precision parameters. + parameter selection index (<0 = all parameters) OUTPUT -. NOTES -. AUTHOR E. Bertin (IAP) -VERSION 18/05/2009 +VERSION 08/07/2010 ***/ -void profit_unboundtobound(profitstruct *profit, float *param) +void profit_unboundtobound(profitstruct *profit, + double *dparam, float *param, int index) { - int p; + int f,p, pstart,np; - for (p=0; pnparam; p++) - if (profit->parammin[p]!=profit->parammax[p]) - param[p] = (profit->parammax[p] - profit->parammin[p]) - / (1.0 + exp(-(param[p]>50.0? 50.0 : param[p]))) - + profit->parammin[p]; + if (index<0) + { + pstart = 0; + np = profit->nparam; + } + else + { + pstart = index; + np = pstart+1; + } + + f = 0; + for (p=pstart; pfreeparam_flag[p]) + { + param[p-pstart] = (profit->parammin[p]!=profit->parammax[p])? + ((profit->parammax[p] - profit->parammin[p]) + / (1.0 + exp(-(dparam[f]>50.0? 50.0 + : (dparam[f]<-50.0? -50.0: dparam[f])))) + + profit->parammin[p]) + : dparam[f]*profit->parammax[p]; + f++; + } + else + param[p-pstart] = profit->paraminit[p]; + } return; } /****** profit_covarunboundtobound ******************************************** -PROTO void profit_covarunboundtobound(profitstruct *profit) +PROTO void profit_covarunboundtobound(profitstruct *profit, + double *dcovar, float *covar) PURPOSE Convert covariance matrix from unbounded to bounded space. -INPUT Pointer to the profit structure. +INPUT Pointer to the profit structure, + input (incomplete) matrix of double precision covariances, + output matrix of single precision covariances. OUTPUT -. NOTES -. AUTHOR E. Bertin (IAP) -VERSION 02/10/2009 +VERSION 08/07/2010 ***/ -void profit_covarunboundtobound(profitstruct *profit) +void profit_covarunboundtobound(profitstruct *profit, + double *dcovar, float *covar) { double *dxdy, dxmin,dxmax; - float *covar, *x,*xmin,*xmax; - int p,p1,p2, nparam; + float *x,*xmin,*xmax; + int *fflag, + f,f1,f2, nfree, p,p1,p2, nparam; nparam = profit->nparam; - QMALLOC16(dxdy, double, nparam); + fflag = profit->freeparam_flag; + nfree = profit->nfreeparam; + QMALLOC16(dxdy, double, nfree); x = profit->paraminit; xmin = profit->parammin; xmax = profit->parammax; + f = 0; for (p=0; pnparam; p++) - if (xmin[p]!=xmax[p]) + if (fflag[p]) { - dxmin = x[p] - xmin[p]; - dxmax= xmax[p] - x[p]; - dxdy[p] = (fabs(dxmin) < 1.0/BIG && fabs(dxmax) < 1.0/BIG) ? + if (xmin[p]!=xmax[p]) + { + dxmin = x[p] - xmin[p]; + dxmax= xmax[p] - x[p]; + dxdy[f++] = (fabs(dxmin) < 1.0/BIG && fabs(dxmax) < 1.0/BIG) ? 0.0 : dxmin*dxmax/(dxmin+dxmax); + } + else + dxdy[f++] = xmax[p]; } - else - dxdy[p] = 1.0; - covar = profit->covar; + memset(profit->covar, 0, nparam*nparam*sizeof(float)); + f2 = 0; for (p2=0; p2typscale = 1.0; break; case PROF_SERSIC: - prof->naxis = 3; - prof->pix = NULL; + prof->naxis = 2; + prof->naxisn[0] = PROFIT_MAXMODSIZE; + prof->naxisn[1] = PROFIT_MAXMODSIZE; + prof->naxisn[2] = 1; + prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2]; + QMALLOC(prof->pix, float, prof->npix); prof->typscale = 1.0; profit_addparam(profit, PARAM_X, &prof->x[0]); profit_addparam(profit, PARAM_Y, &prof->x[1]); @@ -2271,7 +2797,11 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode) break; case PROF_DEVAUCOULEURS: prof->naxis = 2; - prof->pix = NULL; + prof->naxisn[0] = PROFIT_MAXMODSIZE; + prof->naxisn[1] = PROFIT_MAXMODSIZE; + prof->naxisn[2] = 1; + prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2]; + QMALLOC(prof->pix, float, profit->nmodpix); prof->typscale = 1.0; profit_addparam(profit, PARAM_X, &prof->x[0]); profit_addparam(profit, PARAM_Y, &prof->x[1]); @@ -2282,7 +2812,11 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode) break; case PROF_EXPONENTIAL: prof->naxis = 2; - prof->pix = NULL; + prof->naxisn[0] = PROFIT_MAXMODSIZE; + prof->naxisn[1] = PROFIT_MAXMODSIZE; + prof->naxisn[2] = 1; + prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2]; + QMALLOC(prof->pix, float, profit->nmodpix); prof->typscale = 1.0; profit_addparam(profit, PARAM_X, &prof->x[0]); profit_addparam(profit, PARAM_Y, &prof->x[1]); @@ -2293,7 +2827,11 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode) break; case PROF_ARMS: prof->naxis = 2; - prof->pix = NULL; + prof->naxisn[0] = PROFIT_MAXMODSIZE; + prof->naxisn[1] = PROFIT_MAXMODSIZE; + prof->naxisn[2] = 1; + prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2]; + QMALLOC(prof->pix, float, profit->nmodpix); prof->typscale = 1.0; profit_addparam(profit, PARAM_X, &prof->x[0]); profit_addparam(profit, PARAM_Y, &prof->x[1]); @@ -2311,7 +2849,11 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode) break; case PROF_BAR: prof->naxis = 2; - prof->pix = NULL; + prof->naxisn[0] = PROFIT_MAXMODSIZE; + prof->naxisn[1] = PROFIT_MAXMODSIZE; + prof->naxisn[2] = 1; + prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2]; + QMALLOC(prof->pix, float, profit->nmodpix); prof->typscale = 1.0; profit_addparam(profit, PARAM_X, &prof->x[0]); profit_addparam(profit, PARAM_Y, &prof->x[1]); @@ -2325,7 +2867,11 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode) break; case PROF_INRING: prof->naxis = 2; - prof->pix = NULL; + prof->naxisn[0] = PROFIT_MAXMODSIZE; + prof->naxisn[1] = PROFIT_MAXMODSIZE; + prof->naxisn[2] = 1; + prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2]; + QMALLOC(prof->pix, float, profit->nmodpix); prof->typscale = 1.0; profit_addparam(profit, PARAM_X, &prof->x[0]); profit_addparam(profit, PARAM_Y, &prof->x[1]); @@ -2339,7 +2885,11 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode) break; case PROF_OUTRING: prof->naxis = 2; - prof->pix = NULL; + prof->naxisn[0] = PROFIT_MAXMODSIZE; + prof->naxisn[1] = PROFIT_MAXMODSIZE; + prof->naxisn[2] = 1; + prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2]; + QMALLOC(prof->pix, float, profit->nmodpix); prof->typscale = 1.0; profit_addparam(profit, PARAM_X, &prof->x[0]); profit_addparam(profit, PARAM_Y, &prof->x[1]); @@ -2352,7 +2902,10 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode) break; case PROF_DIRAC: prof->naxis = 2; - prof->pix = NULL; + prof->naxisn[0] = PROFIT_MAXMODSIZE; + prof->naxisn[1] = PROFIT_MAXMODSIZE; + prof->naxisn[2] = 1; + prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2]; prof->typscale = 1.0; profit_addparam(profit, PARAM_X, &prof->x[0]); profit_addparam(profit, PARAM_Y, &prof->x[1]); @@ -2403,7 +2956,9 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode) break; } - if (prof->pix) + QMALLOC(prof->pix, float, prof->npix); + + if (prof->naxis>2) { prof->kernelnlines = 1; for (d=0; dnaxis; d++) @@ -2443,16 +2998,18 @@ void prof_end(profstruct *prof) /****** prof_add ************************************************************** -PROTO float prof_add(profstruct *prof, profitstruct *profit) +PROTO float prof_add(profitstruct *profit, profstruct *prof, + int extfluxfac_flag) PURPOSE Add a model profile to an image. -INPUT Profile structure, - profile-fitting structure. -OUTPUT Corrected flux contribution. +INPUT Profile-fitting structure, + profile structure, + flag (0 if flux correction factor is to be computed internally) +OUTPUT Total (asymptotic) flux contribution. NOTES -. AUTHOR E. Bertin (IAP) -VERSION 14/12/2009 +VERSION 08/07/2010 ***/ -float prof_add(profstruct *prof, profitstruct *profit) +float prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag) { double xscale, yscale, saspect, ctheta,stheta, flux, scaling, bn, n; float posin[PROFIT_MAXEXTRA], posout[2], dnaxisn[2], @@ -2469,7 +3026,7 @@ float prof_add(profstruct *prof, profitstruct *profit) int npix, noversamp, threshflag, d,e,i, ix1,ix2, idx1,idx2, nx2, npix2; - npix = profit->modnaxisn[0]*profit->modnaxisn[1]; + npix = profit->nmodpix; if (prof->code==PROF_BACK) { @@ -2524,7 +3081,7 @@ float prof_add(profstruct *prof, profitstruct *profit) /*---- The consequence of sampling on flux is compensated by PSF normalisation*/ x10 = -x1cout - dx1; x2 = -x2cout - dx2; - pixin = profit->pmodpix; + pixin = prof->pix; for (ix2=nx2; ix2--; x2+=1.0) { x1 = x10; @@ -2539,7 +3096,7 @@ float prof_add(profstruct *prof, profitstruct *profit) continue; } val = expf(k*expf(logf(ra)*hinvn)); - noversamp = (int)(val*PROFIT_OVERSAMP+0.1); + noversamp = (int)(val*profit->oversamp+0.1); if (noversamp < 2) *(pixin++) = val; else @@ -2585,7 +3142,7 @@ float prof_add(profstruct *prof, profitstruct *profit) /*---- The consequence of sampling on flux is compensated by PSF normalisation*/ x10 = -x1cout - dx1; x2 = -x2cout - dx2; - pixin = profit->pmodpix; + pixin = prof->pix; for (ix2=nx2; ix2--; x2+=1.0) { x1 = x10; @@ -2600,7 +3157,7 @@ float prof_add(profstruct *prof, profitstruct *profit) continue; } val = expf(-7.66924944f*PROFIT_POWF(ra,0.125)); - noversamp = (int)(sqrt(val)*PROFIT_OVERSAMP+0.1); + noversamp = (int)(sqrt(val)*profit->oversamp+0.1); if (noversamp < 2) *(pixin++) = val; else @@ -2645,7 +3202,7 @@ float prof_add(profstruct *prof, profitstruct *profit) case PROF_EXPONENTIAL: x10 = -x1cout - dx1; x2 = -x2cout - dx2; - pixin = profit->pmodpix; + pixin = prof->pix; for (ix2=nx2; ix2--; x2+=1.0) { x1 = x10; @@ -2660,7 +3217,7 @@ float prof_add(profstruct *prof, profitstruct *profit) continue; } val = expf(-sqrtf(ra)); - noversamp = (int)(val*sqrt(PROFIT_OVERSAMP)+0.1); + noversamp = (int)(val*sqrt(profit->oversamp)+0.1); if (noversamp < 2) *(pixin++) = val; else @@ -2721,7 +3278,7 @@ float prof_add(profstruct *prof, profitstruct *profit) width = 3.0; x1 = -x1cout - dx1; x2 = -x2cout - dx2; - pixin = profit->pmodpix; + pixin = prof->pix; for (ix2=profit->modnaxisn[1]; ix2--; x2+=1.0) { x1t = cd12*x2 + cd11*x1; @@ -2768,7 +3325,7 @@ width = 3.0; sa = sinf(posang); x1 = -x1cout - dx1; x2 = -x2cout - dx2; - pixin = profit->pmodpix; + pixin = prof->pix; for (ix2=profit->modnaxisn[1]; ix2--; x2+=1.0) { x1t = cd12*x2 + cd11*x1; @@ -2806,7 +3363,7 @@ width = 3.0; cd21 /= *prof->feataspect; x1 = -x1cout - dx1; x2 = -x2cout - dx2; - pixin = profit->pmodpix; + pixin = prof->pix; for (ix2=profit->modnaxisn[1]; ix2--; x2+=1.0) { x1t = cd12*x2 + cd11*x1; @@ -2839,7 +3396,7 @@ width = 3.0; invwidth2 = 0.5 / (*prof->featwidth**prof->featwidth); x1 = -x1cout - dx1; x2 = -x2cout - dx2; - pixin = profit->pmodpix; + pixin = prof->pix; for (ix2=profit->modnaxisn[1]; ix2--; x2+=1.0) { x1t = cd12*x2 + cd11*x1; @@ -2862,9 +3419,8 @@ width = 3.0; threshflag = 1; break; case PROF_DIRAC: - memset(profit->pmodpix, 0, - profit->modnaxisn[0]*profit->modnaxisn[1]*sizeof(float)); - profit->pmodpix[profit->modnaxisn[0]/2 + memset(prof->pix, 0, npix*sizeof(float)); + prof->pix[profit->modnaxisn[0]/2 + (profit->modnaxisn[1]/2)*profit->modnaxisn[0]] = 1.0; prof->lostfluxfrac = 0.0; threshflag = 0; @@ -2901,7 +3457,7 @@ width = 3.0; } x1cin = (float)(prof->naxisn[0]/2); x2cin = (float)(prof->naxisn[1]/2); - pixin = profit->pmodpix; + pixin = prof->pix; for (i=npix; i--;) { x1 = posout[0] - x1cout - 1.0 - dx1; @@ -2937,7 +3493,7 @@ width = 3.0; /*-- Find best threshold (the max around the circle with radius rmax */ dx2 = -x2cout; - pixin = profit->pmodpix; + pixin = prof->pix; thresh = -BIG; for (ix2=profit->modnaxisn[1]; ix2--; dx2 += 1.0) { @@ -2949,7 +3505,7 @@ width = 3.0; /*-- Threshold and measure the flux */ flux = 0.0; - pixin = profit->pmodpix; + pixin = prof->pix; for (i=npix; i--; pixin++) if (*pixin >= thresh) flux += (double)*pixin; @@ -2959,27 +3515,95 @@ width = 3.0; else { flux = 0.0; - pixin = profit->pmodpix; + pixin = prof->pix; for (i=npix; i--;) flux += (double)*(pixin++); } - if (prof->lostfluxfrac < 1.0) - flux /= (1.0 - prof->lostfluxfrac); + if (extfluxfac_flag) + fluxfac = prof->fluxfac; + else + { + if (prof->lostfluxfrac < 1.0) + flux /= (1.0 - prof->lostfluxfrac); + + prof->fluxfac = fluxfac = fabs(flux)>0.0? profit->fluxfac/flux : 0.0; + } + + pixin = prof->pix; + for (i=npix; i--;) + *(pixin++) *= fluxfac; /* Correct final flux */ - fluxfac = profit->fluxfac ; - fluxfac *= fabs(flux)>0.0? *prof->flux / flux : 1.0; -// prof->fluxfac = fluxfac; - pixin = profit->pmodpix; + fluxfac = *prof->flux; + pixin = prof->pix; pixout = profit->modpix; for (i=npix; i--;) - *(pixout++) += fluxfac * *(pixin++); + *(pixout++) += fluxfac**(pixin++); return *prof->flux; } +/****** prof_moments ********************************************************** +PROTO void prof_moments(profitstruct *profit, profstruct *prof) +PURPOSE Computes (analytically or numerically) the 2nd moments of a profile + model. +INPUT Profile-fitting structure, + profile structure. +OUTPUT -. +NOTES -. +AUTHOR E. Bertin (IAP) +VERSION 08/07/2010 + ***/ +void prof_moments(profitstruct *profit, profstruct *prof) + { + double m20,m02, ct,st, bn, n; + + if (prof->posangle) + { + switch(prof->code) + { + case PROF_SERSIC: + n = fabs(*prof->extra[0]); + 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 */ + m20 = 0.5 * *prof->scale**prof->scale + * prof_gamma(4.0*n) / (prof_gamma(2.0*n) * pow(bn, 2.0*n)); + break; + case PROF_DEVAUCOULEURS: + m20 = 10.83995 * *prof->scale**prof->scale; + break; + case PROF_EXPONENTIAL: + m20 = 3.0 * *prof->scale**prof->scale; + break; + case PROF_ARMS: + case PROF_BAR: + case PROF_INRING: + case PROF_OUTRING: + m20 = 1.0 / 1.0; + break; + default: + m20 = 0.0; /* to avoid gcc -Wall warnings */ + error(EXIT_FAILURE, "*Internal Error*: Unknown oriented model in ", + "prof_moments()"); + break; + } + + m02 = m20**prof->aspect**prof->aspect; + ct = cos(*prof->posangle*DEG); + st = sin(*prof->posangle*DEG); + prof->mx2 = ct*ct*m20 + st*st*m02; + prof->my2 = st*st*m20 + ct*ct*m02; + prof->mxy = ct*st*(m20 - m02); + } + else + prof->mx2 = prof->my2 = prof->mxy = 0.0; + + return; + } + + /****** prof_interpolate ****************************************************** PROTO float prof_interpolate(profstruct *prof, float *posin) PURPOSE Interpolate a multidimensional model profile at a given position. diff --git a/src/profit.h b/src/profit.h index 5bf3945..4b41eff 100644 --- a/src/profit.h +++ b/src/profit.h @@ -9,7 +9,7 @@ * * Contents: Include file for profit.c. * -* Last modify: 09/10/2009 +* Last modify: 08/07/2010 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -30,11 +30,12 @@ /*----------------------------- Internal constants --------------------------*/ +#define PARAM_ALLPARAMS (-1) /* All parameters */ #define PROFIT_MAXITER 1000 /* Max. nb of iterations in profile fitting */ #define PROFIT_MAXPROF 8 /* Max. nb of profile components */ #define PROFIT_OVERSAMP 5 /* Max. profile oversamp. factor on each axis */ -#define PROFIT_HIDEFRES 201 /* Resolution of the high def. model raster */ -#define PROFIT_REFFFAC 6.0 /* Factor in r_eff for measurement radius*/ +#define PROFIT_HIDEFRES 201 /* Hi. def. model resol. (must be 0 for CCW) */ +/* Buffers */ + double dparam[PARAM_NPARAM]; } profitstruct; /*----------------------------- Global variables ----------------------------*/ @@ -156,39 +170,43 @@ float *profit_compresi(profitstruct *profit, float dynparam, *profit_residuals(profitstruct *profit, picstruct *field, picstruct *wfield, float dynparam, float *param, float *resi), - prof_add(profstruct *prof, profitstruct *profit), + prof_add(profitstruct *profit, profstruct *prof, + int extfluxfac_flag), profit_minradius(profitstruct *profit, float refffac), profit_spiralindex(profitstruct *profit); int profit_copyobjpix(profitstruct *profit, picstruct *field, picstruct *wfield), profit_minimize(profitstruct *profit, int niter), + profit_resample(profitstruct *profit, float *inpix, + PIXTYPE *outpix, float factor), profit_setparam(profitstruct *profit, paramenum paramtype, float param, float parammin, float parammax); void prof_end(profstruct *prof), + prof_moments(profitstruct *profit, profstruct *prof), profit_addparam(profitstruct *profit, paramenum paramindex, float **param), - profit_boundtounbound(profitstruct *profit, float *param), + profit_boundtounbound(profitstruct *profit, + float *param, double *dparam, int index), profit_fit(profitstruct *profit, picstruct *field, picstruct *wfield, objstruct *obj, obj2struct *obj2), profit_convolve(profitstruct *profit, float *modpix), - profit_covarunboundtobound(profitstruct *profit), + profit_covarunboundtobound(profitstruct *profit, + double *dparam, float *param), profit_end(profitstruct *profit), - profit_evaluate(float *par, float *fvec, int m, int n, + profit_evaluate(double *par, double *fvec, int m, int n, void *adata), profit_makedft(profitstruct *profit), profit_moments(profitstruct *profit, obj2struct *obj2), profit_printout(int n_par, float* par, int m_dat, float* fvec, void *data, int iflag, int iter, int nfev ), profit_psf(profitstruct *profit), - profit_resample(profitstruct *profit, float *inpix, - PIXTYPE *outpix, float factor), profit_resetparam(profitstruct *profit, paramenum paramtype), profit_resetparams(profitstruct *profit), - profit_surface(profitstruct *profit, obj2struct *obj2, - double lostfluxfrac), - profit_unboundtobound(profitstruct *profit, float *param); + profit_surface(profitstruct *profit, obj2struct *obj2), + profit_unboundtobound(profitstruct *profit, + double *dparam, float *param, int index); #endif diff --git a/src/psf.c b/src/psf.c index 12ac972..647a30b 100644 --- a/src/psf.c +++ b/src/psf.c @@ -10,7 +10,7 @@ * * Contents: Fit the PSF to a detection. * -* Last modify: 14/12/2009 +* Last modify: 07/07/2010 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -255,7 +255,8 @@ psfstruct *psf_load(char *filename) psf->fwhm = 3.0; /* PSF oversampling: defaulted to 1 */ - if (fitsread(head, "PSF_SAMP", &psf->pixstep,H_FLOAT,T_FLOAT) != RETURN_OK) + if (fitsread(head, "PSF_SAMP", &psf->pixstep,H_FLOAT,T_FLOAT) != RETURN_OK + || psf->pixstep <= 0.0) psf->pixstep = 1.0; /* Load the PSF mask data */ @@ -1016,6 +1017,9 @@ void psf_build(psfstruct *psf) float *ppc, *pl; int i,n,p, ndim, npix; + if (psf->build_flag) + return; + npix = psf->masksize[0]*psf->masksize[1]; /* Reset the Local PSF mask */ @@ -1042,10 +1046,42 @@ void psf_build(psfstruct *psf) *(pl++) += fac**(ppc++); } + psf->build_flag = 1; + return; } +/******************************** psf_fwhm **********************************/ +/* +Return the local PSF FWHM. +*/ +double psf_fwhm(psfstruct *psf) + { + float *pl, + val, max; + int n,p, npix; + + if (!psf->build_flag) + psf_build(psf); + + npix = psf->masksize[0]*psf->masksize[1]; + max = -BIG; + pl = psf->maskloc; + for (p=npix; p--;) + if ((val=*(pl++)) > max) + max = val; + pl = psf->maskloc; + max /= 2.0; + n = 0; + for (p=npix; p--;) + if (*(pl++) >= max) + n++; + + return 2.0*sqrt(n/PI)*psf->pixstep; + } + + /*****************************compute_gradient*********************************/ double *compute_gradient(float *weight,int width, int height, diff --git a/src/psf.h b/src/psf.h index 728b004..1141ad9 100644 --- a/src/psf.h +++ b/src/psf.h @@ -10,7 +10,7 @@ * * Contents: Include file for psffit.c. * -* Last modify: 14/12/2009 +* Last modify: 18/05/2010 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -78,6 +78,7 @@ typedef struct psf pcstruct *pc; /* PC components */ double fwhm; /* Typical PSF FWHM */ float pixstep; /* PSF sampling step */ + int build_flag; /* Set if the current PSF has been computed */ } psfstruct; typedef struct @@ -113,7 +114,8 @@ extern double *compute_gradient (float *weight,int width, int height, float *masks, float *maskx, float *masky, double *mat), *compute_gradient_phot(float *weight,int width, int height, - float *masks, double *pm); + float *masks, double *pm), + psf_fwhm(psfstruct *psf); extern psfstruct *psf_load(char *filename); @@ -129,3 +131,4 @@ extern void pc_end(pcstruct *pc), psf_readcontext(psfstruct *psf, picstruct *field); extern pcstruct *pc_load(catstruct *cat); + diff --git a/src/readimage.c b/src/readimage.c index bb85f9f..4f609b8 100644 --- a/src/readimage.c +++ b/src/readimage.c @@ -9,7 +9,7 @@ * * Contents: functions for input of image data. * -* Last modify: 29/09/2009 +* Last modify: 08/03/2010 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -113,6 +113,10 @@ void *loadstrip(picstruct *field, picstruct *wfield) writecheck(check, data, w); if ((check = prefs.check[CHECK_SUBPROFILES])) writecheck(check, data, w); + if ((check = prefs.check[CHECK_SUBSPHEROIDS])) + writecheck(check, data, w); + if ((check = prefs.check[CHECK_SUBDISKS])) + writecheck(check, data, w); } if ((flags&DETECT_FIELD) && (check=prefs.check[CHECK_BACKRMS])) { @@ -182,6 +186,10 @@ void *loadstrip(picstruct *field, picstruct *wfield) writecheck(check, data, w); if ((check = prefs.check[CHECK_SUBPROFILES])) writecheck(check, data, w); + if ((check = prefs.check[CHECK_SUBSPHEROIDS])) + writecheck(check, data, w); + if ((check = prefs.check[CHECK_SUBDISKS])) + writecheck(check, data, w); } if ((flags&DETECT_FIELD) && (check=prefs.check[CHECK_BACKRMS])) { diff --git a/src/sexhead1.h b/src/sexhead1.h index e909a1f..bbeb221 100644 --- a/src/sexhead1.h +++ b/src/sexhead1.h @@ -9,7 +9,7 @@ * * Contents: header (FITS format #1) and templates for catalog data. * -* Last modify: 19/12/2007 +* Last modify: 02/07/2010 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -95,8 +95,8 @@ keystruct headkey1[] = { &prefs.clean_param, H_FLOAT, T_DOUBLE, "%8f"}, {"SEXCLNST", "CLEANING OBJECT-STACK", &prefs.deblend_nthresh, H_INT, T_LONG, "%6d"}, - {"SEXAPERD", "APERTURE DIAMETER (PIXELS)", - &prefs.apert[0], H_INT, T_LONG, "%7.1f"}, + {"SEXAPERD", "1ST APERTURE DIAMETER (PIXELS)", + &prefs.apert[0], H_FLOAT, T_DOUBLE, "%8.2f"}, {"SEXAPEK1", "KRON PARAMETER", &prefs.autoparam[0], H_FLOAT, T_DOUBLE, "%4.1f"}, {"SEXAPEK2", "KRON ANALYSIS RADIUS", diff --git a/src/types.h b/src/types.h index 9475e5f..2e26872 100644 --- a/src/types.h +++ b/src/types.h @@ -9,7 +9,7 @@ * * Contents: global type definitions. * -* Last modify: 05/02/2010 +* Last modify: 23/03/2010 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -60,9 +60,11 @@ typedef enum {CHECK_NONE, CHECK_IDENTICAL, CHECK_BACKGROUND, CHECK_BACKRMS, CHECK_MINIBACKGROUND, CHECK_MINIBACKRMS, CHECK_SUBTRACTED, CHECK_FILTERED, CHECK_OBJECTS, CHECK_APERTURES, CHECK_SEGMENTATION, CHECK_ASSOC, CHECK_SUBOBJECTS, - CHECK_SUBPSFPROTOS, CHECK_PSFPROTOS, - CHECK_SUBPCPROTOS, CHECK_PCPROTOS, CHECK_PCOPROTOS, - CHECK_MAPSOM, CHECK_SUBPROFILES, CHECK_PROFILES, CHECK_PATTERNS, + CHECK_PSFPROTOS, CHECK_SUBPSFPROTOS, + CHECK_PCPROTOS, CHECK_SUBPCPROTOS, CHECK_PCOPROTOS, + CHECK_MAPSOM, CHECK_PROFILES, CHECK_SUBPROFILES, + CHECK_SPHEROIDS, CHECK_SUBSPHEROIDS, CHECK_DISKS, CHECK_SUBDISKS, + CHECK_PATTERNS,CHECK_OTHER, MAXCHECK} checkenum; /* CHECK_IMAGE type */ diff --git a/src/xml.c b/src/xml.c index be1fbf1..a5ffdc8 100644 --- a/src/xml.c +++ b/src/xml.c @@ -9,7 +9,7 @@ * * Contents: XML logging. * -* Last modify: 05/02/2010 +* Last modify: 25/05/2010 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -172,7 +172,7 @@ INPUT file or stream pointer. OUTPUT RETURN_OK if everything went fine, RETURN_ERROR otherwise. NOTES -. AUTHOR E. Bertin (IAP) -VERSION 14/07/2006 +VERSION 25/05/2010 ***/ int write_xml_header(FILE *file) { @@ -188,7 +188,7 @@ int write_xml_header(FILE *file) fprintf(file, "\n"); fprintf(file, "\n", prefs.xsl_name); - fprintf(file, "\n"); -- GitLab