Skip to content
GitLab
Commit 753719ab authored by Emmanuel Bertin's avatar Emmanuel Bertin
Browse files

Added automatic image rebinning for fitting very large objects (>512 pixels). 

Fixed several issues related to image and weight copy prior to fitting.
Fixed bug with symmetric models.
Fixed issue with new FFT planning scheme.
Pushed version number to 2.9.5.
parent d242beed
Show whitespace changes
Inline Side-by-side
config.h.in
View file @ 753719ab
...@@ -110,6 +110,9 @@ ...@@ -110,6 +110,9 @@
/* Define to 1 if you have the `munmap' function. */ /* Define to 1 if you have the `munmap' function. */
#undef HAVE_MUNMAP #undef HAVE_MUNMAP
/* Define to 1 if you have the `posix_memalign' function. */
#undef HAVE_POSIX_MEMALIGN
/* Define if you have POSIX threads libraries and header files. */ /* Define if you have POSIX threads libraries and header files. */
#undef HAVE_PTHREAD #undef HAVE_PTHREAD
......
configure
View file @ 753719ab
#! /bin/sh #! /bin/sh
# Guess values for system-dependent variables and create Makefiles. # Guess values for system-dependent variables and create Makefiles.
# Generated by GNU Autoconf 2.63 for sextractor 2.9.3. # Generated by GNU Autoconf 2.63 for sextractor 2.9.5.
# #
# Report bugs to <bertin@iap.fr>. # Report bugs to <bertin@iap.fr>.
# #
...@@ -750,8 +750,8 @@ SHELL=${CONFIG_SHELL-/bin/sh} ...@@ -750,8 +750,8 @@ SHELL=${CONFIG_SHELL-/bin/sh}
# Identity of this package. # Identity of this package.
PACKAGE_NAME='sextractor' PACKAGE_NAME='sextractor'
PACKAGE_TARNAME='sextractor' PACKAGE_TARNAME='sextractor'
PACKAGE_VERSION='2.9.3' PACKAGE_VERSION='2.9.5'
PACKAGE_STRING='sextractor 2.9.3' PACKAGE_STRING='sextractor 2.9.5'
PACKAGE_BUGREPORT='bertin@iap.fr' PACKAGE_BUGREPORT='bertin@iap.fr'
   
ac_unique_file="src/makeit.c" ac_unique_file="src/makeit.c"
...@@ -1505,7 +1505,7 @@ if test "$ac_init_help" = "long"; then ...@@ -1505,7 +1505,7 @@ if test "$ac_init_help" = "long"; then
# Omit some internal or obsolete options to make the list less imposing. # Omit some internal or obsolete options to make the list less imposing.
# This message is too long to be a string in the A/UX 3.1 sh. # This message is too long to be a string in the A/UX 3.1 sh.
cat <<_ACEOF cat <<_ACEOF
\`configure' configures sextractor 2.9.3 to adapt to many kinds of systems. \`configure' configures sextractor 2.9.5 to adapt to many kinds of systems.
   
Usage: $0 [OPTION]... [VAR=VALUE]... Usage: $0 [OPTION]... [VAR=VALUE]...
   
...@@ -1575,7 +1575,7 @@ fi ...@@ -1575,7 +1575,7 @@ fi
   
if test -n "$ac_init_help"; then if test -n "$ac_init_help"; then
case $ac_init_help in case $ac_init_help in
short | recursive ) echo "Configuration of sextractor 2.9.3:";; short | recursive ) echo "Configuration of sextractor 2.9.5:";;
esac esac
cat <<\_ACEOF cat <<\_ACEOF
   
...@@ -1706,7 +1706,7 @@ fi ...@@ -1706,7 +1706,7 @@ fi
test -n "$ac_init_help" && exit $ac_status test -n "$ac_init_help" && exit $ac_status
if $ac_init_version; then if $ac_init_version; then
cat <<\_ACEOF cat <<\_ACEOF
sextractor configure 2.9.3 sextractor configure 2.9.5
generated by GNU Autoconf 2.63 generated by GNU Autoconf 2.63
   
Copyright (C) 1992, 1993, 1994, 1995, 1996, 1998, 1999, 2000, 2001, Copyright (C) 1992, 1993, 1994, 1995, 1996, 1998, 1999, 2000, 2001,
...@@ -1720,7 +1720,7 @@ cat >config.log <<_ACEOF ...@@ -1720,7 +1720,7 @@ cat >config.log <<_ACEOF
This file contains any messages produced by compilers while This file contains any messages produced by compilers while
running configure, to aid debugging if configure makes a mistake. running configure, to aid debugging if configure makes a mistake.
   
It was created by sextractor $as_me 2.9.3, which was It was created by sextractor $as_me 2.9.5, which was
generated by GNU Autoconf 2.63. Invocation command line was generated by GNU Autoconf 2.63. Invocation command line was
   
$ $0 $@ $ $0 $@
...@@ -2423,7 +2423,7 @@ fi ...@@ -2423,7 +2423,7 @@ fi
   
# Define the identity of the package. # Define the identity of the package.
PACKAGE='sextractor' PACKAGE='sextractor'
VERSION='2.9.3' VERSION='2.9.5'
   
   
cat >>confdefs.h <<_ACEOF cat >>confdefs.h <<_ACEOF
...@@ -22428,8 +22428,9 @@ done ...@@ -22428,8 +22428,9 @@ done
   
   
   
for ac_func in atexit getenv memcpy memmove memset mkdir munmap strstr \ for ac_func in atexit getenv memcpy memmove memset mkdir munmap strstr \
sincosf logf gettimeofday sincosf logf gettimeofday posix_memalign
do do
as_ac_var=`$as_echo "ac_cv_func_$ac_func" | $as_tr_sh` as_ac_var=`$as_echo "ac_cv_func_$ac_func" | $as_tr_sh`
{ $as_echo "$as_me:$LINENO: checking for $ac_func" >&5 { $as_echo "$as_me:$LINENO: checking for $ac_func" >&5
...@@ -28290,7 +28291,7 @@ exec 6>&1 ...@@ -28290,7 +28291,7 @@ exec 6>&1
# report actual input values of CONFIG_FILES etc. instead of their # report actual input values of CONFIG_FILES etc. instead of their
# values after options handling. # values after options handling.
ac_log=" ac_log="
This file was extended by sextractor $as_me 2.9.3, which was This file was extended by sextractor $as_me 2.9.5, which was
generated by GNU Autoconf 2.63. Invocation command line was generated by GNU Autoconf 2.63. Invocation command line was
   
CONFIG_FILES = $CONFIG_FILES CONFIG_FILES = $CONFIG_FILES
...@@ -28353,7 +28354,7 @@ Report bugs to <bug-autoconf@gnu.org>." ...@@ -28353,7 +28354,7 @@ Report bugs to <bug-autoconf@gnu.org>."
_ACEOF _ACEOF
cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1 cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1
ac_cs_version="\\ ac_cs_version="\\
sextractor config.status 2.9.3 sextractor config.status 2.9.5
configured by $0, generated by GNU Autoconf 2.63, configured by $0, generated by GNU Autoconf 2.63,
with options \\"`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`\\" with options \\"`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`\\"
   
......
configure.ac
View file @ 753719ab
...@@ -6,7 +6,7 @@ define([AC_CACHE_LOAD],) ...@@ -6,7 +6,7 @@ define([AC_CACHE_LOAD],)
define([AC_CACHE_SAVE],) define([AC_CACHE_SAVE],)
# This is your standard Bertin source code... # This is your standard Bertin source code...
AC_INIT(sextractor, 2.9.3, [bertin@iap.fr]) AC_INIT(sextractor, 2.9.5, [bertin@iap.fr])
AC_CONFIG_SRCDIR(src/makeit.c) AC_CONFIG_SRCDIR(src/makeit.c)
AC_CONFIG_AUX_DIR(autoconf) AC_CONFIG_AUX_DIR(autoconf)
AM_CONFIG_HEADER(config.h) AM_CONFIG_HEADER(config.h)
...@@ -93,7 +93,7 @@ AC_TYPE_SIGNAL ...@@ -93,7 +93,7 @@ AC_TYPE_SIGNAL
AC_FUNC_STAT AC_FUNC_STAT
AC_FUNC_STRFTIME AC_FUNC_STRFTIME
AC_CHECK_FUNCS([atexit getenv memcpy memmove memset mkdir munmap strstr \ AC_CHECK_FUNCS([atexit getenv memcpy memmove memset mkdir munmap strstr \
sincosf logf gettimeofday]) sincosf logf gettimeofday posix_memalign])
# Check support for large files # Check support for large files
AC_SYS_LARGEFILE AC_SYS_LARGEFILE
......
man/sex.1
View file @ 753719ab
.TH SEXTRACTOR "1" "October 2009" "SWarp 2.9.3" "User Commands" .TH SEXTRACTOR "1" "October 2009" "SWarp 2.9.5" "User Commands"
.SH NAME .SH NAME
sex \- extract a source catalog from an astronomical FITS image sex \- extract a source catalog from an astronomical FITS image
.SH SYNOPSIS .SH SYNOPSIS
......
src/catout.c
View file @ 753719ab
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: functions for output of catalog data. * Contents: functions for output of catalog data.
* *
* Last modify: 01/10/2009 * Last modify: 07/10/2009
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -381,6 +381,8 @@ void updateparamflags() ...@@ -381,6 +381,8 @@ void updateparamflags()
| FLAG(obj2.prof_vector) | FLAG(obj2.prof_errvector) | FLAG(obj2.prof_vector) | FLAG(obj2.prof_errvector)
| FLAG(obj2.prof_errmatrix) | FLAG(obj2.prof_errmatrix)
| FLAG(obj2.x_prof) | FLAG(obj2.y_prof) | FLAG(obj2.x_prof) | FLAG(obj2.y_prof)
| FLAG(obj2.prof_flag)
| FLAG(obj2.flux_prof)
| FLAG(obj2.prof_mx2) | FLAG(obj2.prof_mx2)
| FLAG(obj2.peak_prof) | FLAG(obj2.peak_prof)
| FLAG(obj2.prof_disk_flux) | FLAG(obj2.prof_disk_flux)
......
src/define.h
View file @ 753719ab
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: global definitions. * Contents: global definitions.
* *
* Last modify: 31/03/2009 * Last modify: 02/10/2009
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -152,6 +152,11 @@ ...@@ -152,6 +152,11 @@
error(EXIT_FAILURE, "Not enough memory for ", \ error(EXIT_FAILURE, "Not enough memory for ", \
#ptr " (" #nel " elements) !");;} #ptr " (" #nel " elements) !");;}
#define QMALLOC16(ptr, typ, nel) \
{if (posix_memalign((void **)&ptr, 16, (size_t)(nel)*sizeof(typ))) \
error(EXIT_FAILURE, "Not enough memory for ", \
#ptr " (" #nel " elements) !");;}
#define QFREE(ptr) \ #define QFREE(ptr) \
{free(ptr); \ {free(ptr); \
ptr = NULL;} ptr = NULL;}
......
src/fft.c
View file @ 753719ab
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: Routines dealing with float precision FFT. * Contents: Routines dealing with float precision FFT.
* *
* Last modify: 26/06/2009 * Last modify: 08/10/2009
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -106,7 +106,7 @@ INPUT -. ...@@ -106,7 +106,7 @@ INPUT -.
OUTPUT -. OUTPUT -.
NOTES -. NOTES -.
AUTHOR E. Bertin (IAP) AUTHOR E. Bertin (IAP)
VERSION 26/06/2009 VERSION 08/10/2009
***/ ***/
void fft_reset(void) void fft_reset(void)
{ {
...@@ -133,11 +133,12 @@ OUTPUT -. ...@@ -133,11 +133,12 @@ OUTPUT -.
NOTES For data1 and fdata2, memory must be allocated for NOTES For data1 and fdata2, memory must be allocated for
size[0]* ... * 2*(size[naxis-1]/2+1) floats (padding required). size[0]* ... * 2*(size[naxis-1]/2+1) floats (padding required).
AUTHOR E. Bertin (IAP) AUTHOR E. Bertin (IAP)
VERSION 26/03/2007 VERSION 08/10/2009
***/ ***/
void fft_conv(float *data1, float *fdata2, int *size) void fft_conv(float *data1, float *fdata2, int *size)
{ {
float *fdata1p,*fdata2p, fftwf_complex *fdata0;
float *fdata1p,*fdata2p,*data0,
real,imag, fac; real,imag, fac;
int i, npix,npix2; int i, npix,npix2;
...@@ -153,7 +154,7 @@ void fft_conv(float *data1, float *fdata2, int *size) ...@@ -153,7 +154,7 @@ void fft_conv(float *data1, float *fdata2, int *size)
{ {
QFFTWMALLOC(fdata1, fftwf_complex, npix2); QFFTWMALLOC(fdata1, fftwf_complex, npix2);
fplan = fftwf_plan_dft_r2c_2d(size[1], size[0], data1, fplan = fftwf_plan_dft_r2c_2d(size[1], size[0], data1,
(fftwf_complex *)fdata1, FFTW_MEASURE|FFTW_DESTROY_INPUT); (fftwf_complex *)fdata1, FFTW_ESTIMATE);
} }
#ifdef USE_THREADS #ifdef USE_THREADS
QPTHREAD_MUTEX_UNLOCK(&fftmutex); QPTHREAD_MUTEX_UNLOCK(&fftmutex);
...@@ -188,7 +189,7 @@ void fft_conv(float *data1, float *fdata2, int *size) ...@@ -188,7 +189,7 @@ void fft_conv(float *data1, float *fdata2, int *size)
#endif #endif
if (!bplan) if (!bplan)
bplan = fftwf_plan_dft_c2r_2d(size[1], size[0], (fftwf_complex *)fdata1, bplan = fftwf_plan_dft_c2r_2d(size[1], size[0], (fftwf_complex *)fdata1,
data1, FFTW_MEASURE|FFTW_DESTROY_INPUT); data1, FFTW_ESTIMATE);
#ifdef USE_THREADS #ifdef USE_THREADS
QPTHREAD_MUTEX_UNLOCK(&fftmutex); QPTHREAD_MUTEX_UNLOCK(&fftmutex);
#endif #endif
...@@ -217,7 +218,7 @@ INPUT ptr to the image, ...@@ -217,7 +218,7 @@ INPUT ptr to the image,
OUTPUT Pointer to the compressed, memory-allocated Fourier transform. OUTPUT Pointer to the compressed, memory-allocated Fourier transform.
NOTES Input data may end up corrupted. NOTES Input data may end up corrupted.
AUTHOR E. Bertin (IAP) AUTHOR E. Bertin (IAP)
VERSION 26/03/2007 VERSION 08/10/2009
***/ ***/
float *fft_rtf(float *data, int *size) float *fft_rtf(float *data, int *size)
{ {
...@@ -234,7 +235,7 @@ float *fft_rtf(float *data, int *size) ...@@ -234,7 +235,7 @@ float *fft_rtf(float *data, int *size)
#endif #endif
QFFTWMALLOC(fdata, fftwf_complex, npix2); QFFTWMALLOC(fdata, fftwf_complex, npix2);
plan = fftwf_plan_dft_r2c_2d(size[1], size[0], data, plan = fftwf_plan_dft_r2c_2d(size[1], size[0], data,
fdata, FFTW_ESTIMATE|FFTW_DESTROY_INPUT); fdata, FFTW_ESTIMATE);
#ifdef USE_THREADS #ifdef USE_THREADS
QPTHREAD_MUTEX_UNLOCK(&fftmutex); QPTHREAD_MUTEX_UNLOCK(&fftmutex);
#endif #endif
......
src/fft.h
View file @ 753719ab
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: Include for fft.c. * Contents: Include for fft.c.
* *
* Last modify: 28/05/2009 * Last modify: 08/10/2009
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
......
src/profit.c
View file @ 753719ab
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: Fit an arbitrary profile combination to a detection. * Contents: Fit an arbitrary profile combination to a detection.
* *
* Last modify: 29/09/2009 * Last modify: 06/10/2009
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -127,7 +127,7 @@ profitstruct *profit_init(psfstruct *psf) ...@@ -127,7 +127,7 @@ profitstruct *profit_init(psfstruct *psf)
nprof++; nprof++;
} }
QMALLOC(profit->covar, float, profit->nparam*profit->nparam); QMALLOC16(profit->covar, float, profit->nparam*profit->nparam);
profit->nprof = nprof; profit->nprof = nprof;
profit->fluxfac = 1.0; /* Default */ profit->fluxfac = 1.0; /* Default */
...@@ -173,7 +173,7 @@ OUTPUT Pointer to an allocated fit structure (containing details about the ...@@ -173,7 +173,7 @@ OUTPUT Pointer to an allocated fit structure (containing details about the
fit). fit).
NOTES It is a modified version of the lm_minimize() of lmfit. NOTES It is a modified version of the lm_minimize() of lmfit.
AUTHOR E. Bertin (IAP) AUTHOR E. Bertin (IAP)
VERSION 24/09/2009 VERSION 06/10/2009
***/ ***/
void profit_fit(profitstruct *profit, void profit_fit(profitstruct *profit,
picstruct *field, picstruct *wfield, picstruct *field, picstruct *wfield,
...@@ -201,20 +201,29 @@ void profit_fit(profitstruct *profit, ...@@ -201,20 +201,29 @@ void profit_fit(profitstruct *profit,
psf = profit->psf; psf = profit->psf;
profit->pixstep = psf->pixstep; profit->pixstep = psf->pixstep;
obj2->prof_flag = 0;
/* Create pixmaps at image resolution */ /* Create pixmaps at image resolution */
profit->ix = (int)(obj->mx + 0.49999);/* internal convention: 1st pix = 0 */
profit->iy = (int)(obj->my + 0.49999);/* internal convention: 1st pix = 0 */
psf_fwhm = psf->masksize[0]*psf->pixstep; psf_fwhm = psf->masksize[0]*psf->pixstep;
profit->objnaxisn[0] = (((int)((obj->xmax-obj->xmin+1) + psf_fwhm + 0.499) profit->objnaxisn[0] = (((int)((obj->xmax-obj->xmin+1) + psf_fwhm + 0.499)
*1.2)/2)*2 + 1; *1.2)/2)*2 + 1;
profit->objnaxisn[1] = (((int)((obj->ymax-obj->ymin+1) + psf_fwhm + 0.499) profit->objnaxisn[1] = (((int)((obj->ymax-obj->ymin+1) + psf_fwhm + 0.499)
*1.2)/2)*2 + 1; *1.2)/2)*2 + 1;
profit->ix = (int)(obj->mx + 0.49999);/* internal convention: 1st pix = 0 */
profit->iy = (int)(obj->my + 0.49999);/* internal convention: 1st pix = 0 */
if (profit->objnaxisn[1]<profit->objnaxisn[0]) if (profit->objnaxisn[1]<profit->objnaxisn[0])
profit->objnaxisn[1] = profit->objnaxisn[0]; profit->objnaxisn[1] = profit->objnaxisn[0];
else else
profit->objnaxisn[0] = profit->objnaxisn[1]; profit->objnaxisn[0] = profit->objnaxisn[1];
if (profit->objnaxisn[0]>PROFIT_MAXOBJSIZE)
{
profit->nsubsamp = ceil((float)profit->objnaxisn[0]/PROFIT_MAXOBJSIZE);
profit->objnaxisn[1] = (profit->objnaxisn[0] /= (int)profit->nsubsamp);
profit->pixstep *= profit->nsubsamp;
obj2->prof_flag |= PROFLAG_OBJSUB;
}
else
profit->nsubsamp = 1.0;
/* Use (dirty) global variables to interface with lmfit */ /* Use (dirty) global variables to interface with lmfit */
the_field = field; the_field = field;
...@@ -223,10 +232,11 @@ void profit_fit(profitstruct *profit, ...@@ -223,10 +232,11 @@ void profit_fit(profitstruct *profit,
profit->obj = obj; profit->obj = obj;
profit->obj2 = obj2; profit->obj2 = obj2;
QMALLOC(profit->objpix, PIXTYPE, profit->objnaxisn[0]*profit->objnaxisn[1]); QMALLOC16(profit->objpix, PIXTYPE, profit->objnaxisn[0]*profit->objnaxisn[1]);
QMALLOC(profit->objweight, PIXTYPE,profit->objnaxisn[0]*profit->objnaxisn[1]); QMALLOC16(profit->objweight, PIXTYPE,profit->objnaxisn[0]*profit->objnaxisn[1]);
QMALLOC(profit->lmodpix, PIXTYPE, profit->objnaxisn[0]*profit->objnaxisn[1]); QMALLOC16(profit->lmodpix, PIXTYPE, profit->objnaxisn[0]*profit->objnaxisn[1]);
profit->nresi = profit_copyobjpix(profit, field, wfield); profit->nresi = profit_copyobjpix(profit, field, wfield);
/* Check if the number of constraints exceeds the number of free parameters */
if (profit->nresi < nparam) if (profit->nresi < nparam)
{ {
if (FLAG(obj2.prof_vector)) if (FLAG(obj2.prof_vector))
...@@ -239,10 +249,11 @@ void profit_fit(profitstruct *profit, ...@@ -239,10 +249,11 @@ void profit_fit(profitstruct *profit,
for (p=0; p<nparam2; p++) for (p=0; p<nparam2; p++)
obj2->prof_errmatrix[p] = 0.0; obj2->prof_errmatrix[p] = 0.0;
obj2->prof_niter = 0; obj2->prof_niter = 0;
obj2->prof_flag |= PROFLAG_NOTCONST;
return; return;
} }
QMALLOC(profit->resi, float, profit->nresi); QMALLOC16(profit->resi, float, profit->nresi);
/* Create pixmap at PSF resolution */ /* Create pixmap at PSF resolution */
profit->modnaxisn[0] = profit->modnaxisn[0] =
...@@ -253,19 +264,23 @@ void profit_fit(profitstruct *profit, ...@@ -253,19 +264,23 @@ void profit_fit(profitstruct *profit,
profit->modnaxisn[1] = profit->modnaxisn[0]; profit->modnaxisn[1] = profit->modnaxisn[0];
else else
profit->modnaxisn[0] = profit->modnaxisn[1]; 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;
}
/* Allocate memory for the complete model */ /* Allocate memory for the complete model */
QFFTWMALLOC(profit->modpix, float, profit->modnaxisn[0]*profit->modnaxisn[1]); QMALLOC16(profit->modpix, float, profit->modnaxisn[0]*profit->modnaxisn[1]);
memset(profit->modpix, 0, profit->modnaxisn[0]*profit->modnaxisn[1]*sizeof(float)); memset(profit->modpix, 0, profit->modnaxisn[0]*profit->modnaxisn[1]*sizeof(float));
QMALLOC(profit->psfpix, float, profit->modnaxisn[0]*profit->modnaxisn[1]); QMALLOC16(profit->psfpix, float, profit->modnaxisn[0]*profit->modnaxisn[1]);
/* Allocate memory for the partial model */ /* Allocate memory for the partial model */
QMALLOC(profit->pmodpix, float, profit->modnaxisn[0]*profit->modnaxisn[1]); QMALLOC16(profit->pmodpix, float, profit->modnaxisn[0]*profit->modnaxisn[1]);
/* Compute the local PSF */ /* Compute the local PSF */
profit_psf(profit); profit_psf(profit);
/* Set initial guesses and boundaries */ /* Set initial guesses and boundaries */
obj2->prof_flag = 0;
profit->sigma = obj->sigbkg; profit->sigma = obj->sigbkg;
profit_resetparams(profit); profit_resetparams(profit);
...@@ -274,6 +289,7 @@ void profit_fit(profitstruct *profit, ...@@ -274,6 +289,7 @@ void profit_fit(profitstruct *profit,
/* Actual minimisation */ /* Actual minimisation */
fft_reset(); fft_reset();
the_gal++;
profit->niter = profit_minimize(profit, PROFIT_MAXITER); profit->niter = profit_minimize(profit, PROFIT_MAXITER);
profit_residuals(profit,field,wfield, 10.0, profit->param,profit->resi); profit_residuals(profit,field,wfield, 10.0, profit->param,profit->resi);
...@@ -613,7 +629,7 @@ void profit_fit(profitstruct *profit, ...@@ -613,7 +629,7 @@ void profit_fit(profitstruct *profit,
pprofit.nparam = 0; pprofit.nparam = 0;
QMALLOC(pprofit.prof, profstruct *, 1); QMALLOC(pprofit.prof, profstruct *, 1);
pprofit.prof[0] = prof_init(&pprofit, PROF_DIRAC); pprofit.prof[0] = prof_init(&pprofit, PROF_DIRAC);
QMALLOC(pprofit.covar, float, pprofit.nparam*pprofit.nparam); QMALLOC16(pprofit.covar, float, pprofit.nparam*pprofit.nparam);
pprofit.nprof = 1; pprofit.nprof = 1;
profit_resetparams(&pprofit); profit_resetparams(&pprofit);
if (profit->paramlist[PARAM_X] && profit->paramlist[PARAM_Y]) if (profit->paramlist[PARAM_X] && profit->paramlist[PARAM_Y])
...@@ -648,6 +664,7 @@ void profit_fit(profitstruct *profit, ...@@ -648,6 +664,7 @@ void profit_fit(profitstruct *profit,
} }
/* clean up. */ /* clean up. */
fft_reset();
free(profit->modpix); free(profit->modpix);
free(profit->psfpix); free(profit->psfpix);
free(profit->pmodpix); free(profit->pmodpix);
...@@ -796,7 +813,7 @@ INPUT Profile-fitting structure. ...@@ -796,7 +813,7 @@ INPUT Profile-fitting structure.
OUTPUT -. OUTPUT -.
NOTES -. NOTES -.
AUTHOR E. Bertin (IAP) AUTHOR E. Bertin (IAP)
VERSION 07/09/2009 VERSION 05/10/2009
***/ ***/
void profit_psf(profitstruct *profit) void profit_psf(profitstruct *profit)
{ {
...@@ -839,7 +856,7 @@ void profit_psf(profitstruct *profit) ...@@ -839,7 +856,7 @@ void profit_psf(profitstruct *profit)
} }
/* Normalize PSF flux (just in case...) */ /* Normalize PSF flux (just in case...) */
flux *= psf->pixstep*psf->pixstep; flux *= profit->pixstep*profit->pixstep;
if (fabs(flux) > 0.0) if (fabs(flux) > 0.0)
{ {
norm = 1.0/flux; norm = 1.0/flux;
...@@ -1276,23 +1293,21 @@ INPUT Pointer to the profit structure, ...@@ -1276,23 +1293,21 @@ INPUT Pointer to the profit structure,
OUTPUT The number of valid pixels copied. OUTPUT The number of valid pixels copied.
NOTES Global preferences are used. NOTES Global preferences are used.
AUTHOR E. Bertin (IAP) AUTHOR E. Bertin (IAP)
VERSION 18/09/2008 VERSION 07/10/2009
***/ ***/
int profit_copyobjpix(profitstruct *profit, picstruct *field, int profit_copyobjpix(profitstruct *profit, picstruct *field,
picstruct *wfield) picstruct *wfield)
{ {
float dx2, dy2, dr2, rad2; float dx, dy2, dr2, rad2;
PIXTYPE *pixin,*pixout, *wpixin,*wpixout, PIXTYPE *pixin,*spixin, *wpixin,*swpixin, *pixout,*wpixout,
backnoise2, invgain, satlevel, wthresh, pix, wpix; backnoise2, invgain, satlevel, wthresh, pix,spix, wpix,swpix;
int i,x,y, xmin,xmax,ymin,ymax, w,h,w2,dw, npix, off, gainflag, int i,x,y, xmin,xmax,ymin,ymax, w,h,dw, npix, off, gainflag,
ix,iy; badflag, sflag, sx,sy,sn,sw, rem, ix,iy;
/* First put the image background to -BIG */ /* First put the image background to -BIG */
pixout = profit->objpix; pixout = profit->objpix;
wpixout = profit->objweight; wpixout = profit->objweight;
w = profit->objnaxisn[0]; for (i=profit->objnaxisn[0]*profit->objnaxisn[1]; i--;)
h = profit->objnaxisn[1];
for (i=w*h; i--;)
{ {
*(pixout++) = -BIG; *(pixout++) = -BIG;
*(wpixout++) = 0.0; *(wpixout++) = 0.0;
...@@ -1305,6 +1320,12 @@ int profit_copyobjpix(profitstruct *profit, picstruct *field, ...@@ -1305,6 +1320,12 @@ int profit_copyobjpix(profitstruct *profit, picstruct *field,
return 0; return 0;
backnoise2 = field->backsig*field->backsig; backnoise2 = field->backsig*field->backsig;
sn = (int)profit->nsubsamp;
sflag = (sn>1);
w = profit->objnaxisn[0]*sn;
h = profit->objnaxisn[1]*sn;
if (sflag)
backnoise2 *= (PIXTYPE)sn;
invgain = (field->gain > 0.0) ? 1.0/field->gain : 0.0; invgain = (field->gain > 0.0) ? 1.0/field->gain : 0.0;
satlevel = field->satur_level - profit->obj->bkg; satlevel = field->satur_level - profit->obj->bkg;
rad2 = h/2.0; rad2 = h/2.0;
...@@ -1312,7 +1333,6 @@ int profit_copyobjpix(profitstruct *profit, picstruct *field, ...@@ -1312,7 +1333,6 @@ int profit_copyobjpix(profitstruct *profit, picstruct *field,
rad2 = w/2.0; rad2 = w/2.0;
rad2 *= rad2; rad2 *= rad2;
/* Set the image boundaries */ /* Set the image boundaries */
pixout = profit->objpix; pixout = profit->objpix;
wpixout = profit->objweight; wpixout = profit->objweight;
...@@ -1320,40 +1340,102 @@ int profit_copyobjpix(profitstruct *profit, picstruct *field, ...@@ -1320,40 +1340,102 @@ int profit_copyobjpix(profitstruct *profit, picstruct *field,
ymax = ymin + h; ymax = ymin + h;
if (ymin<field->ymin) if (ymin<field->ymin)
{ {
off = (field->ymin-ymin)*w; off = (field->ymin-ymin-1)/sn + 1;
pixout += off; pixout += off*profit->objnaxisn[0];
wpixout += off; wpixout += off*profit->objnaxisn[0];
ymin = field->ymin; ymin += off*sn;
} }
if (ymax>field->ymax) if (ymax>field->ymax)
ymax = field->ymax; ymax -= ((ymax-field->ymax-1)/sn + 1)*sn;
xmin = ix-w/2; xmin = ix-w/2;
xmax = xmin + w; xmax = xmin + w;
w2 = w; dw = 0;
if (xmax>field->width) if (xmax>field->width)
{ {
w2 -= xmax-field->width; off = (xmax-field->width-1)/sn + 1;
xmax = field->width; dw += off;
xmax -= off*sn;
} }
if (xmin<0) if (xmin<0)
{ {
pixout -= xmin; off = (-xmin-1)/sn + 1;
wpixout -= xmin; pixout += off;
w2 += xmin; wpixout += off;
xmin = 0; dw += off;
xmin += off*sn;
}
/* Make sure the input frame size is a multiple of the subsampling step */
if (sflag)
{
/*
if (((rem=ymax-ymin)%sn))
{
ymin += rem/2;
ymax -= (rem-rem/2);
}
if (((rem=xmax-xmin)%sn))
{
xmin += rem/2;
pixout += rem/2;
wpixout += rem/2;
dw += rem;
xmax -= (rem-rem/2);
}
*/
sw = field->width;
} }
/* Copy the right pixels to the destination */ /* Copy the right pixels to the destination */
dw = w - w2;
npix = 0; npix = 0;
if (wfield) if (wfield)
{ {
wthresh = wfield->weight_thresh; wthresh = wfield->weight_thresh;
gainflag = prefs.weightgain_flag; gainflag = prefs.weightgain_flag;
/*-- Do the same for the weights */ if (sflag)
npix = 0; {
/*---- Sub-sampling case */
for (y=ymin; y<ymax; y+=sn, pixout+=dw,wpixout+=dw)
{
for (x=xmin; x<xmax; x+=sn)
{
pix = wpix = 0.0;
badflag = 0;
for (sy=0; sy<sn; sy++)
{
dy2 = (y+sy-iy);
dy2 *= dy2;
dx = (x-ix);
spixin = &PIX(field, x, y+sy);
swpixin = &PIX(wfield, x, y+sy);
for (sx=sn; sx--;)
{
dr2 = dy2 + dx*dx;
dx++;
spix = *(spixin++);
swpix = *(swpixin++);
if (dr2<rad2 && spix>-BIG && spix<satlevel && swpix<wthresh)
{
pix += spix;
wpix += swpix;
}
else
badflag=1;
}
}
*(pixout++) = pix;
if (!badflag) /* A single bad pixel ruins is all (saturation, etc.)*/
{
*(wpixout++) = 1.0 / sqrt(wpix+(pix>0.0?
(gainflag? pix*wpix/backnoise2:pix)*invgain : 0.0));
npix++;
}
else
*(wpixout++) = 0.0;
}
}
}
else
for (y=ymin; y<ymax; y++, pixout+=dw,wpixout+=dw) for (y=ymin; y<ymax; y++, pixout+=dw,wpixout+=dw)
{ {
dy2 = y-iy; dy2 = y-iy;
...@@ -1362,15 +1444,56 @@ int profit_copyobjpix(profitstruct *profit, picstruct *field, ...@@ -1362,15 +1444,56 @@ int profit_copyobjpix(profitstruct *profit, picstruct *field,
wpixin = &PIX(wfield, xmin, y); wpixin = &PIX(wfield, xmin, y);
for (x=xmin; x<xmax; x++) for (x=xmin; x<xmax; x++)
{ {
dx2 = (x-ix); dx = x-ix;
dr2 = dy2 + dx2*dx2; dr2 = dy2 + dx*dx;
pix = *(pixout++) = *(pixin++); pix = *(pixin++);
if (dr2<rad2 && pix>-BIG && pix<satlevel && (wpix=*(wpixin++))<wthresh) wpix = *(wpixin++);
if (dr2<rad2 && pix>-BIG && pix<satlevel && wpix<wthresh)
{ {
*(pixout++) = pix;
*(wpixout++) = 1.0 / sqrt(wpix+(pix>0.0? *(wpixout++) = 1.0 / sqrt(wpix+(pix>0.0?
(gainflag? pix*wpix/backnoise2:pix)*invgain : 0.0)); (gainflag? pix*wpix/backnoise2:pix)*invgain : 0.0));
npix++; npix++;
} }
else
*(pixout++) = *(wpixout++) = 0.0;
}
}
}
else
{
if (sflag)
{
/*---- Sub-sampling case */
for (y=ymin; y<ymax; y+=sn, pixout+=dw, wpixout+=dw)
{
for (x=xmin; x<xmax; x+=sn)
{
pix = 0.0;
badflag = 0;
for (sy=0; sy<sn; sy++)
{
dy2 = y+sy-iy;
dy2 *= dy2;
dx = x-ix;
spixin = &PIX(field, x, y+sy);
for (sx=sn; sx--;)
{
dr2 = dy2 + dx*dx;
dx++;
spix = *(spixin++);
if (dr2<rad2 && spix>-BIG && spix<satlevel)
pix += spix;
else
badflag=1;
}
}
*(pixout++) = pix;
if (!badflag) /* A single bad pixel ruins is all (saturation, etc.)*/
{
*(wpixout++) = 1.0 / sqrt(backnoise2 + (pix>0.0?pix*invgain:0.0));
npix++;
}
else else
*(wpixout++) = 0.0; *(wpixout++) = 0.0;
} }
...@@ -1384,16 +1507,18 @@ int profit_copyobjpix(profitstruct *profit, picstruct *field, ...@@ -1384,16 +1507,18 @@ int profit_copyobjpix(profitstruct *profit, picstruct *field,
pixin = &PIX(field, xmin, y); pixin = &PIX(field, xmin, y);
for (x=xmin; x<xmax; x++) for (x=xmin; x<xmax; x++)
{ {
dx2 = (x-ix); dx = x-ix;
dr2 = dy2 + dx2*dx2; dr2 = dy2 + dx*dx;
pix = *(pixout++) = *(pixin++); pix = *(pixin++);
if (dr2<rad2 && pix>-BIG) if (dr2<rad2 && pix>-BIG && pix<satlevel)
{ {
*(pixout++) = pix;
*(wpixout++) = 1.0 / sqrt(backnoise2 + (pix>0.0?pix*invgain : 0.0)); *(wpixout++) = 1.0 / sqrt(backnoise2 + (pix>0.0?pix*invgain : 0.0));
npix++; npix++;
} }
else else
*(wpixout++) = 0.0; *(pixout++) = *(wpixout++) = 0.0;
}
} }
} }
...@@ -1786,7 +1911,7 @@ void profit_resetparam(profitstruct *profit, paramenum paramtype) ...@@ -1786,7 +1911,7 @@ void profit_resetparam(profitstruct *profit, paramenum paramtype)
break; break;
case PARAM_DISK_SCALE: /* From scalelength to Re */ case PARAM_DISK_SCALE: /* From scalelength to Re */
param = obj2->hl_radius/1.67835*sqrt(obj->a/obj->b); param = obj2->hl_radius/1.67835*sqrt(obj->a/obj->b);
parammin = 0.1; parammin = param/4.0;
parammax = param * 4.0; parammax = param * 4.0;
break; break;
case PARAM_DISK_ASPECT: case PARAM_DISK_ASPECT:
...@@ -2017,7 +2142,7 @@ INPUT Pointer to the profit structure. ...@@ -2017,7 +2142,7 @@ INPUT Pointer to the profit structure.
OUTPUT -. OUTPUT -.
NOTES -. NOTES -.
AUTHOR E. Bertin (IAP) AUTHOR E. Bertin (IAP)
VERSION 07/09/2009 VERSION 02/10/2009
***/ ***/
void profit_covarunboundtobound(profitstruct *profit) void profit_covarunboundtobound(profitstruct *profit)
{ {
...@@ -2027,7 +2152,7 @@ void profit_covarunboundtobound(profitstruct *profit) ...@@ -2027,7 +2152,7 @@ void profit_covarunboundtobound(profitstruct *profit)
int p,p1,p2, nparam; int p,p1,p2, nparam;
nparam = profit->nparam; nparam = profit->nparam;
QMALLOC(dxdy, double, nparam); QMALLOC16(dxdy, double, nparam);
x = profit->paraminit; x = profit->paraminit;
xmin = profit->parammin; xmin = profit->parammin;
xmax = profit->parammax; xmax = profit->parammax;
...@@ -2237,7 +2362,7 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode) ...@@ -2237,7 +2362,7 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode)
(prof->kernelwidth[d] = interp_kernwidth[prof->interptype[d]]); (prof->kernelwidth[d] = interp_kernwidth[prof->interptype[d]]);
} }
prof->kernelnlines /= prof->kernelwidth[0]; prof->kernelnlines /= prof->kernelwidth[0];
QMALLOC(prof->kernelbuf, float, prof->kernelnlines); QMALLOC16(prof->kernelbuf, float, prof->kernelnlines);
} }
return prof; return prof;
...@@ -2274,7 +2399,7 @@ INPUT Profile structure, ...@@ -2274,7 +2399,7 @@ INPUT Profile structure,
OUTPUT Corrected flux contribution. OUTPUT Corrected flux contribution.
NOTES -. NOTES -.
AUTHOR E. Bertin (IAP) AUTHOR E. Bertin (IAP)
VERSION 29/09/2009 VERSION 08/10/2009
***/ ***/
float prof_add(profstruct *prof, profitstruct *profit) float prof_add(profstruct *prof, profitstruct *profit)
{ {
...@@ -2449,7 +2574,7 @@ float prof_add(profstruct *prof, profitstruct *profit) ...@@ -2449,7 +2574,7 @@ float prof_add(profstruct *prof, profitstruct *profit)
/*---- Copy the symmetric part */ /*---- Copy the symmetric part */
if ((npix2=(profit->modnaxisn[1]-nx2)*profit->modnaxisn[0]) > 0) if ((npix2=(profit->modnaxisn[1]-nx2)*profit->modnaxisn[0]) > 0)
{ {
pixin2 = pixin - profit->modnaxisn[0] - 1; pixin2 = pixin - profit->modnaxisn[0];
for (i=npix2; i--;) for (i=npix2; i--;)
*(pixin++) = *(pixin2--); *(pixin++) = *(pixin2--);
} }
...@@ -2504,7 +2629,7 @@ float prof_add(profstruct *prof, profitstruct *profit) ...@@ -2504,7 +2629,7 @@ float prof_add(profstruct *prof, profitstruct *profit)
/*---- Copy the symmetric part */ /*---- Copy the symmetric part */
if ((npix2=(profit->modnaxisn[1]-nx2)*profit->modnaxisn[0]) > 0) if ((npix2=(profit->modnaxisn[1]-nx2)*profit->modnaxisn[0]) > 0)
{ {
pixin2 = pixin - profit->modnaxisn[0] - 1; pixin2 = pixin - profit->modnaxisn[0];
for (i=npix2; i--;) for (i=npix2; i--;)
*(pixin++) = *(pixin2--); *(pixin++) = *(pixin2--);
} }
......
src/profit.h
View file @ 753719ab
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: Include file for profit.c. * Contents: Include file for profit.c.
* *
* Last modify: 24/09/2009 * Last modify: 07/10/2009
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -19,7 +19,9 @@ ...@@ -19,7 +19,9 @@
/*-------------------------------- flags ------------------------------------*/ /*-------------------------------- flags ------------------------------------*/
#define PROFIT_FLIPPED 0x0001 #define PROFLAG_MODSUB 0x0001
#define PROFLAG_OBJSUB 0x0002
#define PROFLAG_NOTCONST 0x0004
/*-------------------------------- macros -----------------------------------*/ /*-------------------------------- macros -----------------------------------*/
...@@ -34,6 +36,8 @@ ...@@ -34,6 +36,8 @@
#define PROFIT_HIDEFRES 201 /* Resolution of the high def. model raster */ #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_REFFFAC 6.0 /* Factor in r_eff for measurement radius*/
#define PROFIT_DYNPARAM 10.0 /* Dynamic compression param. in sigma units */ #define PROFIT_DYNPARAM 10.0 /* Dynamic compression param. in sigma units */
#define PROFIT_MAXMODSIZE 1024 /* Maximum size allowed for the model raster */
#define PROFIT_MAXOBJSIZE 512 /* Maximum size allowed for the object raster */
#define PROFIT_BARXFADE 0.1 /* Fract. of bar length crossfaded with arms */ #define PROFIT_BARXFADE 0.1 /* Fract. of bar length crossfaded with arms */
#define PROFIT_MAXEXTRA 2 /* Max. nb of extra free params of profiles */ #define PROFIT_MAXEXTRA 2 /* Max. nb of extra free params of profiles */
#define PROFIT_PROFRES 256 /* Pixmap size of model components */ #define PROFIT_PROFRES 256 /* Pixmap size of model components */
...@@ -121,6 +125,7 @@ typedef struct ...@@ -121,6 +125,7 @@ typedef struct
struct psf *psf; /* PSF */ struct psf *psf; /* PSF */
float pixstep; /* Model/PSF sampling step */ float pixstep; /* Model/PSF sampling step */
float fluxfac; /* Model flux scaling factor */ float fluxfac; /* Model flux scaling factor */
float nsubsamp; /* Subsampling factor */
float *psfdft; /* Compressed Fourier Transform of the PSF */ float *psfdft; /* Compressed Fourier Transform of the PSF */
float *psfpix; /* Full res. pixmap of the PSF */ float *psfpix; /* Full res. pixmap of the PSF */
float *modpix; /* Full res. pixmap of the complete model */ float *modpix; /* Full res. pixmap of the complete model */
......
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