public software.psfex

Compare Revisions

Rev 190 → Rev 191

/branches/rhl/man/psfex.1
1,4 → 1,4
.TH PSFEX "1" "June 2013" "PSFEx 3.15.2" "User Commands"
.TH PSFEX "1" "July 2013" "PSFEx 3.15.2" "User Commands"
.SH NAME
psfex \- generate a PSF super-tabulated model for SExtractor
.SH SYNOPSIS
/branches/rhl/src/field.h
80,6 → 80,7
extern void field_count(fieldstruct **fields, setstruct *set,
int counttype),
field_end(fieldstruct *field),
field_init_finalize(fieldstruct *field),
field_locate(fieldstruct *field),
field_psfsave(fieldstruct *field, char *filename),
field_stats(fieldstruct **fields, setstruct *set);
/branches/rhl/src/Makefile.am
34,8 → 34,8
endif
bin_PROGRAMS = psfex
psfex_SOURCES = check.c context.c $(CPLOTSOURCE) diagnostic.c fft.c \
field.c fitswcs.c homo.c main.c makeit.c misc.c \
pca.c prefs.c psf.c sample.c vignet.c xml.c \
field.c field_utils.c fitswcs.c homo.c main.c makeit.c makeit2.c misc.c \
pca.c prefs.c psf.c sample.c sample_utils.c vignet.c xml.c \
check.h context.h cplot.h define.h diagnostic.h \
fft.h field.h fitswcs.h globals.h homo.h key.h \
misc.h pca.h prefs.h preflist.h psf.h \
/branches/rhl/src/globals.h
25,10 → 25,20
* Last modified: 06/01/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
#include <time.h>
 
/*----------------------- miscellaneous variables ---------------------------*/
char gstr[MAXCHAR];
 
/*------------------------------- functions ---------------------------------*/
extern void makeit(void);
#include "fits/fitscat.h"
#include "psf.h"
#include "sample.h"
#include "field.h"
 
extern void error(int, char *, char *),
makeit(void),
makeit_body(fieldstruct **fields, contextstruct **context,
contextstruct **fullcontext, int free_sets);
psfstruct *make_psf(setstruct *set, float psfstep,
float *basis, int nbasis, contextstruct *context);
/branches/rhl/src/makeit2.c New file
0,0 → 1,532
/*
* makeit2.c
*
* Main loop, with the I/O parts removed
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
* This file part of: PSFEx
*
* Copyright: (C) 1997-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
*
* License: GNU General Public License
*
* PSFEx is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
* PSFEx is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
* You should have received a copy of the GNU General Public License
* along with PSFEx. If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 19/07/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
 
#ifdef USE_THREADS
#ifdef HAVE_MKL
#include MKL_H
#endif
#endif
 
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
 
#include "define.h"
#include "types.h"
#include "globals.h"
#include "fits/fitscat.h"
#include "check.h"
#include "context.h"
#include "cplot.h"
#include "diagnostic.h"
#include "field.h"
#include "homo.h"
#include "pca.h"
#include "prefs.h"
#include "psf.h"
#include "sample.h"
#include "xml.h"
 
psfstruct *make_psf(setstruct *set, float psfstep,
float *basis, int nbasis, contextstruct *context);
 
/********************************** makeit_body ******************************/
/*
*/
void
makeit_body(
fieldstruct **fields,
contextstruct **context_,
contextstruct **fullcontext_,
int free_sets /* if false, don't free any sets as we don't own them */
)
{
wcsstruct *wcs;
fieldstruct *field;
psfstruct **cpsf,
*psf;
setstruct *set, *set2;
contextstruct *context, *fullcontext;
char str[MAXCHAR];
char **incatnames;
float **psfbasiss,
*psfsteps, *psfbasis, *basis,
psfstep, step;
int c,i,p, ncat, ext, next, nmed, nbasis;
 
incatnames = prefs.incat_name;
ncat = prefs.ncat;
next = fields[0]->next;
 
psfstep = prefs.psf_step;
psfsteps = NULL;
nbasis = 0;
psfbasis = NULL;
psfbasiss = NULL;
 
/* Initialize context */
NFPRINTF(OUTPUT, "Initializing contexts...");
context = context_init(prefs.context_name, prefs.context_group,
prefs.ncontext_group, prefs.group_deg, prefs.ngroup_deg,
CONTEXT_REMOVEHIDDEN);
fullcontext = context->npc?
context_init(prefs.context_name, prefs.context_group,
prefs.ncontext_group, prefs.group_deg, prefs.ngroup_deg,
CONTEXT_KEEPHIDDEN)
: context;
*context_ = context;
*fullcontext_ = fullcontext;
 
if (context->npc && ncat<2)
warning("Hidden dependencies cannot be derived from",
" a single catalog");
else if (context->npc && prefs.stability_type == STABILITY_EXPOSURE)
warning("Hidden dependencies have no effect",
" in STABILITY_TYPE EXPOSURE mode");
 
/* Compute PSF steps */
if (!prefs.psf_step)
{
NFPRINTF(OUTPUT, "Computing optimum PSF sampling steps...");
if (prefs.newbasis_type==NEWBASIS_PCACOMMON
|| (prefs.stability_type == STABILITY_SEQUENCE
&& prefs.psf_mef_type == PSF_MEF_COMMON))
{
set = load_samples(incatnames, 0, ncat, ALL_EXTENSIONS, next, context);
psfstep = (float)((set->fwhm/2.35)*0.5);
if (free_sets) end_set(set);
}
/*-- Need to derive a common pixel step for each ext */
else if (prefs.newbasis_type == NEWBASIS_PCAINDEPENDENT
|| context->npc
|| (prefs.stability_type == STABILITY_SEQUENCE
&& prefs.psf_mef_type == PSF_MEF_INDEPENDENT))
{
/*-- Run through all samples to derive a different pixel step for each extension */
QMALLOC(psfsteps, float, next);
for (ext=0 ; ext<next; ext++)
{
set = load_samples(incatnames, 0, ncat, ext, next, context);
psfsteps[ext] = (float)(psfstep? psfstep : (set->fwhm/2.35)*0.5);
if (free_sets) end_set(set);
}
}
}
 
/* Derive a new common PCA basis for all extensions */
if (prefs.newbasis_type==NEWBASIS_PCACOMMON)
{
QMALLOC(cpsf, psfstruct *, ncat*next);
for (ext=0 ; ext<next; ext++)
for (c=0; c<ncat; c++)
{
sprintf(str, "Computing new PCA image basis from %s...",
fields[c]->rtcatname);
NFPRINTF(OUTPUT, str);
set = load_samples(incatnames, c, 1, ext, next, context);
step = psfstep;
cpsf[c+ext*ncat] = make_psf(set, psfstep, NULL, 0, context);
if (free_sets) end_set(set);
}
nbasis = prefs.newbasis_number;
psfbasis = pca_onsnaps(cpsf, ncat*next, nbasis);
for (i=0 ; i<ncat*next; i++)
psf_end(cpsf[i]);
free(cpsf);
}
/* Derive a new PCA basis for each extension */
else if (prefs.newbasis_type == NEWBASIS_PCAINDEPENDENT)
{
nbasis = prefs.newbasis_number;
QMALLOC(psfbasiss, float *, next);
for (ext=0; ext<next; ext++)
{
if (psfsteps)
step = psfsteps[ext];
else
step = psfstep;
QMALLOC(cpsf, psfstruct *, ncat);
for (c=0; c<ncat; c++)
{
sprintf(str, "Computing new PCA image basis from %s...",
fields[c]->rtcatname);
NFPRINTF(OUTPUT, str);
set = load_samples(incatnames, c, 1, ext, next, context);
cpsf[c] = make_psf(set, step, NULL, 0, context);
if (free_sets) end_set(set);
}
psfbasiss[ext] = pca_onsnaps(cpsf, ncat, nbasis);
for (c=0 ; c<ncat; c++)
psf_end(cpsf[c]);
free(cpsf);
}
}
 
if (context->npc && prefs.hidden_mef_type == HIDDEN_MEF_COMMON)
/*-- Derive principal components of PSF variation from the whole mosaic */
{
p = 0;
QMALLOC(cpsf, psfstruct *, ncat*next);
for (c=0; c<ncat; c++)
{
sprintf(str, "Computing hidden dependency parameter(s) from %s...",
fields[c]->rtcatname);
NFPRINTF(OUTPUT, str);
for (ext=0 ; ext<next; ext++)
{
set = load_samples(incatnames, c, 1, ext, next, context);
if (psfsteps)
step = psfsteps[ext];
else
step = psfstep;
basis = psfbasiss? psfbasiss[ext] : psfbasis;
cpsf[p++] = make_psf(set, step, basis, nbasis, context);
if (free_sets) end_set(set);
}
}
free(fullcontext->pc);
fullcontext->pc = pca_oncomps(cpsf, next, ncat, context->npc);
for (c=0 ; c<ncat*next; c++)
psf_end(cpsf[c]);
free(cpsf);
}
 
/* Compute "final" PSF models */
if (prefs.psf_mef_type == PSF_MEF_COMMON)
{
if (prefs.stability_type == STABILITY_SEQUENCE)
{
/*---- Load all the samples at once */
set = load_samples(incatnames, 0, ncat, ALL_EXTENSIONS, next, context);
step = psfstep;
basis = psfbasis;
field_count(fields, set, COUNT_LOADED);
psf = make_psf(set, step, basis, nbasis, context);
field_count(fields, set, COUNT_ACCEPTED);
if (free_sets) end_set(set);
NFPRINTF(OUTPUT, "Computing final PSF model...");
context_apply(context, psf, fields, ALL_EXTENSIONS, 0, ncat);
psf_end(psf);
}
else
for (c=0; c<ncat; c++)
{
/*------ Load the samples for current exposure */
sprintf(str, "Computing final PSF model from %s...",
fields[c]->rtcatname);
NFPRINTF(OUTPUT, str);
set = load_samples(incatnames, c, 1, ALL_EXTENSIONS, next, context);
if (psfstep)
step = psfstep;
else
step = (float)((set->fwhm/2.35)*0.5);
basis = psfbasis;
field_count(fields, set, COUNT_LOADED);
psf = make_psf(set, step, basis, nbasis, fullcontext);
field_count(fields, set, COUNT_ACCEPTED);
if (free_sets) end_set(set);
context_apply(fullcontext, psf, fields, ALL_EXTENSIONS, c, 1);
psf_end(psf);
}
}
else
for (ext=0 ; ext<next; ext++)
{
basis = psfbasiss? psfbasiss[ext] : psfbasis;
if (context->npc && prefs.hidden_mef_type == HIDDEN_MEF_INDEPENDENT)
/*------ Derive principal components of PSF components */
{
QMALLOC(cpsf, psfstruct *, ncat);
if (psfsteps)
step = psfsteps[ext];
else
step = psfstep;
for (c=0; c<ncat; c++)
{
if (next>1)
sprintf(str,
"Computing hidden dependency parameter(s) from %s[%d/%d]...",
fields[c]->rtcatname, ext+1, next);
else
sprintf(str, "Computing hidden dependency parameter(s) from %s...",
fields[c]->rtcatname);
NFPRINTF(OUTPUT, str);
set = load_samples(incatnames, c, 1, ext, next, context);
field_count(fields, set, COUNT_LOADED);
cpsf[c] = make_psf(set, step, basis, nbasis, context);
field_count(fields, set, COUNT_ACCEPTED);
if (free_sets) end_set(set);
}
free(fullcontext->pc);
fullcontext->pc = pca_oncomps(cpsf, 1, ncat, context->npc);
for (c=0 ; c<ncat; c++)
psf_end(cpsf[c]);
free(cpsf);
}
 
if (prefs.stability_type == STABILITY_SEQUENCE)
{
/*------ Load all the samples at once */
if (next>1)
{
sprintf(str, "Computing final PSF model for extension [%d/%d]...",
ext+1, next);
NFPRINTF(OUTPUT, str);
}
else
NFPRINTF(OUTPUT, "Computing final PSF model...");
set = load_samples(incatnames, 0, ncat, ext, next, fullcontext);
if (psfstep)
step = psfstep;
else if (psfsteps)
step = psfsteps[ext];
else
step = (float)((set->fwhm/2.35)*0.5);
field_count(fields, set, COUNT_LOADED);
psf = make_psf(set, step, basis, nbasis, fullcontext);
field_count(fields, set, COUNT_ACCEPTED);
if (free_sets) end_set(set);
context_apply(fullcontext, psf, fields, ext, 0, ncat);
psf_end(psf);
}
else
for (c=0; c<ncat; c++)
{
/*-------- Load the samples for current exposure */
if (next>1)
sprintf(str, "Reading data from %s[%d/%d]...",
fields[c]->rtcatname, ext+1, next);
else
sprintf(str, "Reading data from %s...",
fields[c]->rtcatname);
NFPRINTF(OUTPUT, str);
set = load_samples(incatnames, c, 1, ext, next, context);
if (psfstep)
step = psfstep;
else if (psfsteps)
step = psfsteps[ext];
else
step = (float)((set->fwhm/2.35)*0.5);
if (next>1)
sprintf(str, "Computing final PSF model for %s[%d/%d]...",
fields[c]->rtcatname, ext+1, next);
else
sprintf(str, "Computing final PSF model for %s...",
fields[c]->rtcatname);
NFPRINTF(OUTPUT, str);
field_count(fields, set, COUNT_LOADED);
psf = make_psf(set, step, basis, nbasis, context);
field_count(fields, set, COUNT_ACCEPTED);
if (free_sets) end_set(set);
context_apply(context, psf, fields, ext, c, 1);
psf_end(psf);
}
}
 
free(psfsteps);
if (psfbasiss)
{
for (ext=0; ext<next; ext++)
free(psfbasiss[ext]);
free(psfbasiss);
}
else if (psfbasis)
free(psfbasis);
 
/* Compute diagnostics and check-images */
#ifdef USE_THREADS
/* Force MKL using a single thread as diagnostic code is already multithreaded*/
#ifdef HAVE_MKL
if (prefs.context_nsnap>2)
mkl_set_num_threads(1);
#endif
#endif
QIPRINTF(OUTPUT,
" filename [ext] accepted/total samp. chi2/dof FWHM ellip."
" resi. asym.");
for (c=0; c<ncat; c++)
{
field = fields[c];
for (ext=0 ; ext<next; ext++)
{
psf = field->psf[ext];
wcs = field->wcs[ext];
if (next>1)
sprintf(str, "Computing diagnostics for %s[%d/%d]...",
field->rtcatname, ext+1, next);
else
sprintf(str, "Computing diagnostics for %s...",
field->rtcatname);
NFPRINTF(OUTPUT, str);
/*---- Check PSF with individual datasets */
set2 = load_samples(incatnames, c, 1, ext, next, context);
psf->samples_loaded = set2->nsample;
if (set2->nsample>1)
{
/*------ Remove bad PSF candidates */
psf_clean(psf, set2, prefs.prof_accuracy);
psf->chi2 = set2->nsample? psf_chi2(psf, set2) : 0.0;
}
psf->samples_accepted = set2->nsample;
/*---- Compute diagnostics and field statistics */
psf_diagnostic(psf);
psf_wcsdiagnostic(psf, wcs);
nmed = psf->nmed;
field_stats(fields, set2);
/*---- Display stats for current catalog/extension */
if (next>1)
sprintf(str, "[%d/%d]", ext+1, next);
else
str[0] = '\0';
QPRINTF(OUTPUT, "%-17.17s%-7.7s %5d/%-5d %6.2f %6.2f %6.2f %4.2f"
" %5.2f %5.2f\n",
ext==0? field->rtcatname : "",
str,
psf->samples_accepted, psf->samples_loaded,
psf->pixstep,
psf->chi2,
psf->moffat_fwhm,
psf->moffat_ellipticity,
psf->pfmoffat_residuals,
psf->sym_residuals);
/*---- Save "Check-images" */
for (i=0; i<prefs.ncheck_type; i++)
if (prefs.check_type[i])
{
sprintf(str, "Saving CHECK-image #%d...", i+1);
NFPRINTF(OUTPUT, str);
check_write(field, set2, prefs.check_name[i], prefs.check_type[i],
ext, next, prefs.check_cubeflag);
}
/*---- Free memory */
if (free_sets) end_set(set2);
}
}
 
#ifdef USE_THREADS
/* Back to multithreaded MKL */
#ifdef HAVE_MKL
mkl_set_num_threads(prefs.nthreads);
#endif
#endif
 
return;
}
 
/****** make_psf *************************************************************
PROTO psfstruct *make_psf(setstruct *set, float psfstep,
float *basis, int nbasis, contextstruct *context)
PURPOSE Make PSFs from a set of FITS binary catalogs.
INPUT Pointer to a sample set,
PSF sampling step,
Pointer to basis image vectors,
Number of basis vectors,
Pointer to context structure.
OUTPUT Pointer to the PSF structure.
NOTES Diagnostics are computed only if diagflag != 0.
AUTHOR E. Bertin (IAP)
VERSION 03/09/2009
***/
psfstruct *make_psf(setstruct *set, float psfstep,
float *basis, int nbasis, contextstruct *context)
{
psfstruct *psf;
basistypenum basistype;
float pixsize[2];
 
pixsize[0] = (float)prefs.psf_pixsize[0];
pixsize[1] = (float)prefs.psf_pixsize[1];
// NFPRINTF(OUTPUT,"Initializing PSF modules...");
psf = psf_init(context, prefs.psf_size, psfstep, pixsize, set->nsample);
 
psf->samples_loaded = set->nsample;
psf->fwhm = set->fwhm;
/* Make the basic PSF-model (1st pass) */
// NFPRINTF(OUTPUT,"Modeling the PSF (1/3)...");
psf_make(psf, set, 0.2);
if (basis && nbasis)
{
QMEMCPY(basis, psf->basis, float, nbasis*psf->size[0]*psf->size[1]);
psf->nbasis = nbasis;
}
else
{
// NFPRINTF(OUTPUT,"Generating the PSF model...");
basistype = prefs.basis_type;
if (basistype==BASIS_PIXEL_AUTO)
basistype = (psf->fwhm < PSF_AUTO_FWHM)? BASIS_PIXEL : BASIS_NONE;
psf_makebasis(psf, set, basistype, prefs.basis_number);
}
psf_refine(psf, set);
 
/* Remove bad PSF candidates */
if (set->nsample>1)
{
psf_clean(psf, set, 0.2);
/*-- Make the basic PSF-model (2nd pass) */
// NFPRINTF(OUTPUT,"Modeling the PSF (2/3)...");
psf_make(psf, set, 0.1);
psf_refine(psf, set);
}
/* Remove bad PSF candidates */
if (set->nsample>1)
{
psf_clean(psf, set, 0.1);
/*-- Make the basic PSF-model (3rd pass) */
// NFPRINTF(OUTPUT,"Modeling the PSF (3/3)...");
psf_make(psf, set, 0.05);
psf_refine(psf, set);
}
/* Remove bad PSF candidates */
if (set->nsample>1)
psf_clean(psf, set, 0.05);
psf->samples_accepted = set->nsample;
/* Refine the PSF-model */
psf_make(psf, set, prefs.prof_accuracy);
psf_refine(psf, set);
/* Clip the PSF-model */
psf_clip(psf);
 
/*-- Just check the Chi2 */
psf->chi2 = set->nsample? psf_chi2(psf, set) : 0.0;
 
return psf;
}
/branches/rhl/src/Makefile.in
1,8 → 1,8
# Makefile.in generated by automake 1.12.2 from Makefile.am.
# Makefile.in generated by automake 1.10 from Makefile.am.
# @configure_input@
 
# Copyright (C) 1994-2012 Free Software Foundation, Inc.
 
# Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002,
# 2003, 2004, 2005, 2006 Free Software Foundation, Inc.
# This Makefile.in is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
# with or without modifications, as long as this notice is preserved.
45,27 → 45,9
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
VPATH = @srcdir@
am__make_dryrun = \
{ \
am__dry=no; \
case $$MAKEFLAGS in \
*\\[\ \ ]*) \
echo 'am--echo: ; @echo "AM" OK' | $(MAKE) -f - 2>/dev/null \
| grep '^AM OK$$' >/dev/null || am__dry=yes;; \
*) \
for am__flg in $$MAKEFLAGS; do \
case $$am__flg in \
*=*|--*) ;; \
*n*) am__dry=yes; break;; \
esac; \
done;; \
esac; \
test $$am__dry = yes; \
}
pkgdatadir = $(datadir)/@PACKAGE@
pkgincludedir = $(includedir)/@PACKAGE@
pkglibdir = $(libdir)/@PACKAGE@
pkglibexecdir = $(libexecdir)/@PACKAGE@
pkgincludedir = $(includedir)/@PACKAGE@
am__cd = CDPATH="$${ZSH_VERSION+.}$(PATH_SEPARATOR)" && cd
install_sh_DATA = $(install_sh) -c -m 644
install_sh_PROGRAM = $(install_sh) -c
82,8 → 64,7
host_triplet = @host@
bin_PROGRAMS = psfex$(EXEEXT)
subdir = src
DIST_COMMON = $(srcdir)/Makefile.am $(srcdir)/Makefile.in \
$(top_srcdir)/autoconf/depcomp
DIST_COMMON = $(srcdir)/Makefile.am $(srcdir)/Makefile.in
ACLOCAL_M4 = $(top_srcdir)/aclocal.m4
am__aclocal_m4_deps = $(top_srcdir)/acx_atlas.m4 \
$(top_srcdir)/acx_fftw.m4 $(top_srcdir)/acx_mkl.m4 \
96,28 → 77,30
mkinstalldirs = $(install_sh) -d
CONFIG_HEADER = $(top_builddir)/config.h
CONFIG_CLEAN_FILES =
CONFIG_CLEAN_VPATH_FILES =
am__installdirs = "$(DESTDIR)$(bindir)"
binPROGRAMS_INSTALL = $(INSTALL_PROGRAM)
PROGRAMS = $(bin_PROGRAMS)
am__psfex_SOURCES_DIST = check.c context.c cplot.c diagnostic.c fft.c \
field.c fitswcs.c homo.c main.c makeit.c misc.c pca.c prefs.c \
psf.c sample.c vignet.c xml.c check.h context.h cplot.h \
define.h diagnostic.h fft.h field.h fitswcs.h globals.h homo.h \
key.h misc.h pca.h prefs.h preflist.h psf.h sample.h threads.h \
types.h vignet.h wcscelsys.h xml.h
field.c field_utils.c fitswcs.c homo.c main.c makeit.c \
makeit2.c misc.c pca.c prefs.c psf.c sample.c sample_utils.c \
vignet.c xml.c check.h context.h cplot.h define.h diagnostic.h \
fft.h field.h fitswcs.h globals.h homo.h key.h misc.h pca.h \
prefs.h preflist.h psf.h sample.h threads.h types.h vignet.h \
wcscelsys.h xml.h
@USE_PLPLOT_TRUE@am__objects_1 = cplot.$(OBJEXT)
am_psfex_OBJECTS = check.$(OBJEXT) context.$(OBJEXT) $(am__objects_1) \
diagnostic.$(OBJEXT) fft.$(OBJEXT) field.$(OBJEXT) \
fitswcs.$(OBJEXT) homo.$(OBJEXT) main.$(OBJEXT) \
makeit.$(OBJEXT) misc.$(OBJEXT) pca.$(OBJEXT) prefs.$(OBJEXT) \
psf.$(OBJEXT) sample.$(OBJEXT) vignet.$(OBJEXT) xml.$(OBJEXT)
field_utils.$(OBJEXT) fitswcs.$(OBJEXT) homo.$(OBJEXT) \
main.$(OBJEXT) makeit.$(OBJEXT) makeit2.$(OBJEXT) \
misc.$(OBJEXT) pca.$(OBJEXT) prefs.$(OBJEXT) psf.$(OBJEXT) \
sample.$(OBJEXT) sample_utils.$(OBJEXT) vignet.$(OBJEXT) \
xml.$(OBJEXT)
psfex_OBJECTS = $(am_psfex_OBJECTS)
psfex_DEPENDENCIES = $(srcdir)/fits/libfits.a \
$(srcdir)/levmar/liblevmar.a $(srcdir)/wcs/libwcs_c.a
DEFAULT_INCLUDES = -I.@am__isrc@ -I$(top_builddir)
DEFAULT_INCLUDES = -I. -I$(top_builddir)@am__isrc@
depcomp = $(SHELL) $(top_srcdir)/autoconf/depcomp
am__depfiles_maybe = depfiles
am__mv = mv -f
COMPILE = $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) \
$(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS)
LTCOMPILE = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) \
136,45 → 119,12
install-pdf-recursive install-ps-recursive install-recursive \
installcheck-recursive installdirs-recursive pdf-recursive \
ps-recursive uninstall-recursive
am__can_run_installinfo = \
case $$AM_UPDATE_INFO_DIR in \
n|no|NO) false;; \
*) (install-info --version) >/dev/null 2>&1;; \
esac
RECURSIVE_CLEAN_TARGETS = mostlyclean-recursive clean-recursive \
distclean-recursive maintainer-clean-recursive
AM_RECURSIVE_TARGETS = $(RECURSIVE_TARGETS:-recursive=) \
$(RECURSIVE_CLEAN_TARGETS:-recursive=) tags TAGS ctags CTAGS \
distdir
ETAGS = etags
CTAGS = ctags
DIST_SUBDIRS = $(SUBDIRS)
DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST)
am__relativize = \
dir0=`pwd`; \
sed_first='s,^\([^/]*\)/.*$$,\1,'; \
sed_rest='s,^[^/]*/*,,'; \
sed_last='s,^.*/\([^/]*\)$$,\1,'; \
sed_butlast='s,/*[^/]*$$,,'; \
while test -n "$$dir1"; do \
first=`echo "$$dir1" | sed -e "$$sed_first"`; \
if test "$$first" != "."; then \
if test "$$first" = ".."; then \
dir2=`echo "$$dir0" | sed -e "$$sed_last"`/"$$dir2"; \
dir0=`echo "$$dir0" | sed -e "$$sed_butlast"`; \
else \
first2=`echo "$$dir2" | sed -e "$$sed_first"`; \
if test "$$first2" = "$$first"; then \
dir2=`echo "$$dir2" | sed -e "$$sed_rest"`; \
else \
dir2="../$$dir2"; \
fi; \
dir0="$$dir0"/"$$first"; \
fi; \
fi; \
dir1=`echo "$$dir1" | sed -e "$$sed_rest"`; \
done; \
reldir="$$dir2"
ACLOCAL = @ACLOCAL@
AMTAR = @AMTAR@
AM_CFLAGS = @AM_CFLAGS@
200,7 → 150,6
DATE3 = @DATE3@
DEFS = @DEFS@
DEPDIR = @DEPDIR@
DLLTOOL = @DLLTOOL@
DSYMUTIL = @DSYMUTIL@
DUMPBIN = @DUMPBIN@
ECHO_C = @ECHO_C@
227,7 → 176,6
LN_S = @LN_S@
LTLIBOBJS = @LTLIBOBJS@
MAKEINFO = @MAKEINFO@
MANIFEST_TOOL = @MANIFEST_TOOL@
MKDIR_P = @MKDIR_P@
MKL_CFLAGS = @MKL_CFLAGS@
MKL_LDFLAGS = @MKL_LDFLAGS@
263,7 → 211,6
abs_srcdir = @abs_srcdir@
abs_top_builddir = @abs_top_builddir@
abs_top_srcdir = @abs_top_srcdir@
ac_ct_AR = @ac_ct_AR@
ac_ct_CC = @ac_ct_CC@
ac_ct_DUMPBIN = @ac_ct_DUMPBIN@
am__include = @am__include@
314,8 → 261,8
SUBDIRS = fits levmar wcs
@USE_PLPLOT_TRUE@CPLOTSOURCE = cplot.c
psfex_SOURCES = check.c context.c $(CPLOTSOURCE) diagnostic.c fft.c \
field.c fitswcs.c homo.c main.c makeit.c misc.c \
pca.c prefs.c psf.c sample.c vignet.c xml.c \
field.c field_utils.c fitswcs.c homo.c main.c makeit.c makeit2.c misc.c \
pca.c prefs.c psf.c sample.c sample_utils.c vignet.c xml.c \
check.h context.h cplot.h define.h diagnostic.h \
fft.h field.h fitswcs.h globals.h homo.h key.h \
misc.h pca.h prefs.h preflist.h psf.h \
334,14 → 281,14
@for dep in $?; do \
case '$(am__configure_deps)' in \
*$$dep*) \
( cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh ) \
&& { if test -f $@; then exit 0; else break; fi; }; \
cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh \
&& exit 0; \
exit 1;; \
esac; \
done; \
echo ' cd $(top_srcdir) && $(AUTOMAKE) --foreign src/Makefile'; \
$(am__cd) $(top_srcdir) && \
$(AUTOMAKE) --foreign src/Makefile
echo ' cd $(top_srcdir) && $(AUTOMAKE) --gnu src/Makefile'; \
cd $(top_srcdir) && \
$(AUTOMAKE) --gnu src/Makefile
.PRECIOUS: Makefile
Makefile: $(srcdir)/Makefile.in $(top_builddir)/config.status
@case '$?' in \
359,54 → 306,35
cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh
$(ACLOCAL_M4): $(am__aclocal_m4_deps)
cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh
$(am__aclocal_m4_deps):
install-binPROGRAMS: $(bin_PROGRAMS)
@$(NORMAL_INSTALL)
@list='$(bin_PROGRAMS)'; test -n "$(bindir)" || list=; \
if test -n "$$list"; then \
echo " $(MKDIR_P) '$(DESTDIR)$(bindir)'"; \
$(MKDIR_P) "$(DESTDIR)$(bindir)" || exit 1; \
fi; \
for p in $$list; do echo "$$p $$p"; done | \
sed 's/$(EXEEXT)$$//' | \
while read p p1; do if test -f $$p || test -f $$p1; \
then echo "$$p"; echo "$$p"; else :; fi; \
done | \
sed -e 'p;s,.*/,,;n;h' -e 's|.*|.|' \
-e 'p;x;s,.*/,,;s/$(EXEEXT)$$//;$(transform);s/$$/$(EXEEXT)/' | \
sed 'N;N;N;s,\n, ,g' | \
$(AWK) 'BEGIN { files["."] = ""; dirs["."] = 1 } \
{ d=$$3; if (dirs[d] != 1) { print "d", d; dirs[d] = 1 } \
if ($$2 == $$4) files[d] = files[d] " " $$1; \
else { print "f", $$3 "/" $$4, $$1; } } \
END { for (d in files) print "f", d, files[d] }' | \
while read type dir files; do \
if test "$$dir" = .; then dir=; else dir=/$$dir; fi; \
test -z "$$files" || { \
echo " $(INSTALL_PROGRAM_ENV) $(LIBTOOL) $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=install $(INSTALL_PROGRAM) $$files '$(DESTDIR)$(bindir)$$dir'"; \
$(INSTALL_PROGRAM_ENV) $(LIBTOOL) $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=install $(INSTALL_PROGRAM) $$files "$(DESTDIR)$(bindir)$$dir" || exit $$?; \
} \
; done
test -z "$(bindir)" || $(MKDIR_P) "$(DESTDIR)$(bindir)"
@list='$(bin_PROGRAMS)'; for p in $$list; do \
p1=`echo $$p|sed 's/$(EXEEXT)$$//'`; \
if test -f $$p \
|| test -f $$p1 \
; then \
f=`echo "$$p1" | sed 's,^.*/,,;$(transform);s/$$/$(EXEEXT)/'`; \
echo " $(INSTALL_PROGRAM_ENV) $(LIBTOOL) --mode=install $(binPROGRAMS_INSTALL) '$$p' '$(DESTDIR)$(bindir)/$$f'"; \
$(INSTALL_PROGRAM_ENV) $(LIBTOOL) --mode=install $(binPROGRAMS_INSTALL) "$$p" "$(DESTDIR)$(bindir)/$$f" || exit 1; \
else :; fi; \
done
 
uninstall-binPROGRAMS:
@$(NORMAL_UNINSTALL)
@list='$(bin_PROGRAMS)'; test -n "$(bindir)" || list=; \
files=`for p in $$list; do echo "$$p"; done | \
sed -e 'h;s,^.*/,,;s/$(EXEEXT)$$//;$(transform)' \
-e 's/$$/$(EXEEXT)/' `; \
test -n "$$list" || exit 0; \
echo " ( cd '$(DESTDIR)$(bindir)' && rm -f" $$files ")"; \
cd "$(DESTDIR)$(bindir)" && rm -f $$files
@list='$(bin_PROGRAMS)'; for p in $$list; do \
f=`echo "$$p" | sed 's,^.*/,,;s/$(EXEEXT)$$//;$(transform);s/$$/$(EXEEXT)/'`; \
echo " rm -f '$(DESTDIR)$(bindir)/$$f'"; \
rm -f "$(DESTDIR)$(bindir)/$$f"; \
done
 
clean-binPROGRAMS:
@list='$(bin_PROGRAMS)'; test -n "$$list" || exit 0; \
echo " rm -f" $$list; \
rm -f $$list || exit $$?; \
test -n "$(EXEEXT)" || exit 0; \
list=`for p in $$list; do echo "$$p"; done | sed 's/$(EXEEXT)$$//'`; \
echo " rm -f" $$list; \
rm -f $$list
psfex$(EXEEXT): $(psfex_OBJECTS) $(psfex_DEPENDENCIES) $(EXTRA_psfex_DEPENDENCIES)
@list='$(bin_PROGRAMS)'; for p in $$list; do \
f=`echo $$p|sed 's/$(EXEEXT)$$//'`; \
echo " rm -f $$p $$f"; \
rm -f $$p $$f ; \
done
psfex$(EXEEXT): $(psfex_OBJECTS) $(psfex_DEPENDENCIES)
@rm -f psfex$(EXEEXT)
$(LINK) $(psfex_OBJECTS) $(psfex_LDADD) $(LIBS)
 
422,35 → 350,38
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/diagnostic.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/fft.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/field.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/field_utils.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/fitswcs.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/homo.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/main.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/makeit.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/makeit2.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/misc.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pca.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/prefs.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/psf.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/sample.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/sample_utils.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/vignet.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/xml.Po@am__quote@
 
.c.o:
@am__fastdepCC_TRUE@ $(COMPILE) -MT $@ -MD -MP -MF $(DEPDIR)/$*.Tpo -c -o $@ $<
@am__fastdepCC_TRUE@ $(am__mv) $(DEPDIR)/$*.Tpo $(DEPDIR)/$*.Po
@am__fastdepCC_TRUE@ mv -f $(DEPDIR)/$*.Tpo $(DEPDIR)/$*.Po
@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='$<' object='$@' libtool=no @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(COMPILE) -c $<
 
.c.obj:
@am__fastdepCC_TRUE@ $(COMPILE) -MT $@ -MD -MP -MF $(DEPDIR)/$*.Tpo -c -o $@ `$(CYGPATH_W) '$<'`
@am__fastdepCC_TRUE@ $(am__mv) $(DEPDIR)/$*.Tpo $(DEPDIR)/$*.Po
@am__fastdepCC_TRUE@ mv -f $(DEPDIR)/$*.Tpo $(DEPDIR)/$*.Po
@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='$<' object='$@' libtool=no @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(COMPILE) -c `$(CYGPATH_W) '$<'`
 
.c.lo:
@am__fastdepCC_TRUE@ $(LTCOMPILE) -MT $@ -MD -MP -MF $(DEPDIR)/$*.Tpo -c -o $@ $<
@am__fastdepCC_TRUE@ $(am__mv) $(DEPDIR)/$*.Tpo $(DEPDIR)/$*.Plo
@am__fastdepCC_TRUE@ mv -f $(DEPDIR)/$*.Tpo $(DEPDIR)/$*.Plo
@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='$<' object='$@' libtool=yes @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(LTCOMPILE) -c -o $@ $<
462,13 → 393,13
-rm -rf .libs _libs
 
# This directory's subdirectories are mostly independent; you can cd
# into them and run 'make' without going through this Makefile.
# To change the values of 'make' variables: instead of editing Makefiles,
# (1) if the variable is set in 'config.status', edit 'config.status'
# (which will cause the Makefiles to be regenerated when you run 'make');
# (2) otherwise, pass the desired values on the 'make' command line.
$(RECURSIVE_TARGETS) $(RECURSIVE_CLEAN_TARGETS):
@fail= failcom='exit 1'; \
# into them and run `make' without going through this Makefile.
# To change the values of `make' variables: instead of editing Makefiles,
# (1) if the variable is set in `config.status', edit `config.status'
# (which will cause the Makefiles to be regenerated when you run `make');
# (2) otherwise, pass the desired values on the `make' command line.
$(RECURSIVE_TARGETS):
@failcom='exit 1'; \
for f in x $$MAKEFLAGS; do \
case $$f in \
*=* | --[!k]*);; \
477,11 → 408,7
done; \
dot_seen=no; \
target=`echo $@ | sed s/-recursive//`; \
case "$@" in \
distclean-* | maintainer-clean-*) list='$(DIST_SUBDIRS)' ;; \
*) list='$(SUBDIRS)' ;; \
esac; \
for subdir in $$list; do \
list='$(SUBDIRS)'; for subdir in $$list; do \
echo "Making $$target in $$subdir"; \
if test "$$subdir" = "."; then \
dot_seen=yes; \
489,38 → 416,65
else \
local_target="$$target"; \
fi; \
($(am__cd) $$subdir && $(MAKE) $(AM_MAKEFLAGS) $$local_target) \
(cd $$subdir && $(MAKE) $(AM_MAKEFLAGS) $$local_target) \
|| eval $$failcom; \
done; \
if test "$$dot_seen" = "no"; then \
$(MAKE) $(AM_MAKEFLAGS) "$$target-am" || exit 1; \
fi; test -z "$$fail"
 
$(RECURSIVE_CLEAN_TARGETS):
@failcom='exit 1'; \
for f in x $$MAKEFLAGS; do \
case $$f in \
*=* | --[!k]*);; \
*k*) failcom='fail=yes';; \
esac; \
done; \
dot_seen=no; \
case "$@" in \
distclean-* | maintainer-clean-*) list='$(DIST_SUBDIRS)' ;; \
*) list='$(SUBDIRS)' ;; \
esac; \
rev=''; for subdir in $$list; do \
if test "$$subdir" = "."; then :; else \
rev="$$subdir $$rev"; \
fi; \
done; \
rev="$$rev ."; \
target=`echo $@ | sed s/-recursive//`; \
for subdir in $$rev; do \
echo "Making $$target in $$subdir"; \
if test "$$subdir" = "."; then \
local_target="$$target-am"; \
else \
local_target="$$target"; \
fi; \
(cd $$subdir && $(MAKE) $(AM_MAKEFLAGS) $$local_target) \
|| eval $$failcom; \
done && test -z "$$fail"
tags-recursive:
list='$(SUBDIRS)'; for subdir in $$list; do \
test "$$subdir" = . || ($(am__cd) $$subdir && $(MAKE) $(AM_MAKEFLAGS) tags); \
test "$$subdir" = . || (cd $$subdir && $(MAKE) $(AM_MAKEFLAGS) tags); \
done
ctags-recursive:
list='$(SUBDIRS)'; for subdir in $$list; do \
test "$$subdir" = . || ($(am__cd) $$subdir && $(MAKE) $(AM_MAKEFLAGS) ctags); \
test "$$subdir" = . || (cd $$subdir && $(MAKE) $(AM_MAKEFLAGS) ctags); \
done
cscopelist-recursive:
list='$(SUBDIRS)'; for subdir in $$list; do \
test "$$subdir" = . || ($(am__cd) $$subdir && $(MAKE) $(AM_MAKEFLAGS) cscopelist); \
done
 
ID: $(HEADERS) $(SOURCES) $(LISP) $(TAGS_FILES)
list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \
unique=`for i in $$list; do \
if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \
done | \
$(AWK) '{ files[$$0] = 1; nonempty = 1; } \
END { if (nonempty) { for (i in files) print i; }; }'`; \
$(AWK) ' { files[$$0] = 1; } \
END { for (i in files) print i; }'`; \
mkid -fID $$unique
tags: TAGS
 
TAGS: tags-recursive $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \
$(TAGS_FILES) $(LISP)
set x; \
tags=; \
here=`pwd`; \
if ($(ETAGS) --etags-include --version) >/dev/null 2>&1; then \
include_option=--etags-include; \
532,58 → 486,40
list='$(SUBDIRS)'; for subdir in $$list; do \
if test "$$subdir" = .; then :; else \
test ! -f $$subdir/TAGS || \
set "$$@" "$$include_option=$$here/$$subdir/TAGS"; \
tags="$$tags $$include_option=$$here/$$subdir/TAGS"; \
fi; \
done; \
list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \
unique=`for i in $$list; do \
if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \
done | \
$(AWK) '{ files[$$0] = 1; nonempty = 1; } \
END { if (nonempty) { for (i in files) print i; }; }'`; \
shift; \
if test -z "$(ETAGS_ARGS)$$*$$unique"; then :; else \
$(AWK) ' { files[$$0] = 1; } \
END { for (i in files) print i; }'`; \
if test -z "$(ETAGS_ARGS)$$tags$$unique"; then :; else \
test -n "$$unique" || unique=$$empty_fix; \
if test $$# -gt 0; then \
$(ETAGS) $(ETAGSFLAGS) $(AM_ETAGSFLAGS) $(ETAGS_ARGS) \
"$$@" $$unique; \
else \
$(ETAGS) $(ETAGSFLAGS) $(AM_ETAGSFLAGS) $(ETAGS_ARGS) \
$$unique; \
fi; \
$(ETAGS) $(ETAGSFLAGS) $(AM_ETAGSFLAGS) $(ETAGS_ARGS) \
$$tags $$unique; \
fi
ctags: CTAGS
CTAGS: ctags-recursive $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \
$(TAGS_FILES) $(LISP)
tags=; \
here=`pwd`; \
list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \
unique=`for i in $$list; do \
if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \
done | \
$(AWK) '{ files[$$0] = 1; nonempty = 1; } \
END { if (nonempty) { for (i in files) print i; }; }'`; \
test -z "$(CTAGS_ARGS)$$unique" \
$(AWK) ' { files[$$0] = 1; } \
END { for (i in files) print i; }'`; \
test -z "$(CTAGS_ARGS)$$tags$$unique" \
|| $(CTAGS) $(CTAGSFLAGS) $(AM_CTAGSFLAGS) $(CTAGS_ARGS) \
$$unique
$$tags $$unique
 
GTAGS:
here=`$(am__cd) $(top_builddir) && pwd` \
&& $(am__cd) $(top_srcdir) \
&& gtags -i $(GTAGS_ARGS) "$$here"
&& cd $(top_srcdir) \
&& gtags -i $(GTAGS_ARGS) $$here
 
cscopelist: cscopelist-recursive $(HEADERS) $(SOURCES) $(LISP)
list='$(SOURCES) $(HEADERS) $(LISP)'; \
case "$(srcdir)" in \
[\\/]* | ?:[\\/]*) sdir="$(srcdir)" ;; \
*) sdir=$(subdir)/$(srcdir) ;; \
esac; \
for i in $$list; do \
if test -f "$$i"; then \
echo "$(subdir)/$$i"; \
else \
echo "$$sdir/$$i"; \
fi; \
done >> $(top_builddir)/cscope.files
 
distclean-tags:
-rm -f TAGS ID GTAGS GRTAGS GSYMS GPATH tags
 
603,41 → 539,29
if test -f $$file || test -d $$file; then d=.; else d=$(srcdir); fi; \
if test -d $$d/$$file; then \
dir=`echo "/$$file" | sed -e 's,/[^/]*$$,,'`; \
if test -d "$(distdir)/$$file"; then \
find "$(distdir)/$$file" -type d ! -perm -700 -exec chmod u+rwx {} \;; \
fi; \
if test -d $(srcdir)/$$file && test $$d != $(srcdir); then \
cp -fpR $(srcdir)/$$file "$(distdir)$$dir" || exit 1; \
find "$(distdir)/$$file" -type d ! -perm -700 -exec chmod u+rwx {} \;; \
cp -pR $(srcdir)/$$file $(distdir)$$dir || exit 1; \
fi; \
cp -fpR $$d/$$file "$(distdir)$$dir" || exit 1; \
cp -pR $$d/$$file $(distdir)$$dir || exit 1; \
else \
test -f "$(distdir)/$$file" \
|| cp -p $$d/$$file "$(distdir)/$$file" \
test -f $(distdir)/$$file \
|| cp -p $$d/$$file $(distdir)/$$file \
|| exit 1; \
fi; \
done
@list='$(DIST_SUBDIRS)'; for subdir in $$list; do \
list='$(DIST_SUBDIRS)'; for subdir in $$list; do \
if test "$$subdir" = .; then :; else \
$(am__make_dryrun) \
|| test -d "$(distdir)/$$subdir" \
|| $(MKDIR_P) "$(distdir)/$$subdir" \
|| exit 1; \
dir1=$$subdir; dir2="$(distdir)/$$subdir"; \
$(am__relativize); \
new_distdir=$$reldir; \
dir1=$$subdir; dir2="$(top_distdir)"; \
$(am__relativize); \
new_top_distdir=$$reldir; \
echo " (cd $$subdir && $(MAKE) $(AM_MAKEFLAGS) top_distdir="$$new_top_distdir" distdir="$$new_distdir" \\"; \
echo " am__remove_distdir=: am__skip_length_check=: am__skip_mode_fix=: distdir)"; \
($(am__cd) $$subdir && \
test -d "$(distdir)/$$subdir" \
|| $(MKDIR_P) "$(distdir)/$$subdir" \
|| exit 1; \
distdir=`$(am__cd) $(distdir) && pwd`; \
top_distdir=`$(am__cd) $(top_distdir) && pwd`; \
(cd $$subdir && \
$(MAKE) $(AM_MAKEFLAGS) \
top_distdir="$$new_top_distdir" \
distdir="$$new_distdir" \
top_distdir="$$top_distdir" \
distdir="$$distdir/$$subdir" \
am__remove_distdir=: \
am__skip_length_check=: \
am__skip_mode_fix=: \
distdir) \
|| exit 1; \
fi; \
660,22 → 584,16
 
installcheck: installcheck-recursive
install-strip:
if test -z '$(STRIP)'; then \
$(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM="$(INSTALL_STRIP_PROGRAM)" \
install_sh_PROGRAM="$(INSTALL_STRIP_PROGRAM)" INSTALL_STRIP_FLAG=-s \
install; \
else \
$(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM="$(INSTALL_STRIP_PROGRAM)" \
install_sh_PROGRAM="$(INSTALL_STRIP_PROGRAM)" INSTALL_STRIP_FLAG=-s \
"INSTALL_PROGRAM_ENV=STRIPPROG='$(STRIP)'" install; \
fi
$(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM="$(INSTALL_STRIP_PROGRAM)" \
install_sh_PROGRAM="$(INSTALL_STRIP_PROGRAM)" INSTALL_STRIP_FLAG=-s \
`test -z '$(STRIP)' || \
echo "INSTALL_PROGRAM_ENV=STRIPPROG='$(STRIP)'"` install
mostlyclean-generic:
 
clean-generic:
 
distclean-generic:
-test -z "$(CONFIG_CLEAN_FILES)" || rm -f $(CONFIG_CLEAN_FILES)
-test . = "$(srcdir)" || test -z "$(CONFIG_CLEAN_VPATH_FILES)" || rm -f $(CONFIG_CLEAN_VPATH_FILES)
 
maintainer-clean-generic:
@echo "This command is intended for maintainers to use"
696,8 → 614,6
 
html: html-recursive
 
html-am:
 
info: info-recursive
 
info-am:
706,28 → 622,18
 
install-dvi: install-dvi-recursive
 
install-dvi-am:
 
install-exec-am: install-binPROGRAMS
 
install-html: install-html-recursive
 
install-html-am:
 
install-info: install-info-recursive
 
install-info-am:
 
install-man:
 
install-pdf: install-pdf-recursive
 
install-pdf-am:
 
install-ps: install-ps-recursive
 
install-ps-am:
 
installcheck-am:
 
maintainer-clean: maintainer-clean-recursive
750,27 → 656,25
 
uninstall-am: uninstall-binPROGRAMS
 
.MAKE: $(RECURSIVE_CLEAN_TARGETS) $(RECURSIVE_TARGETS) \
cscopelist-recursive ctags-recursive install-am install-strip \
tags-recursive
.MAKE: $(RECURSIVE_CLEAN_TARGETS) $(RECURSIVE_TARGETS) install-am \
install-strip
 
.PHONY: $(RECURSIVE_CLEAN_TARGETS) $(RECURSIVE_TARGETS) CTAGS GTAGS \
all all-am check check-am clean clean-binPROGRAMS \
clean-generic clean-libtool cscopelist cscopelist-recursive \
ctags ctags-recursive distclean distclean-compile \
distclean-generic distclean-libtool distclean-tags distdir dvi \
dvi-am html html-am info info-am install install-am \
install-binPROGRAMS install-data install-data-am install-dvi \
install-dvi-am install-exec install-exec-am install-html \
install-html-am install-info install-info-am install-man \
install-pdf install-pdf-am install-ps install-ps-am \
install-strip installcheck installcheck-am installdirs \
installdirs-am maintainer-clean maintainer-clean-generic \
mostlyclean mostlyclean-compile mostlyclean-generic \
mostlyclean-libtool pdf pdf-am ps ps-am tags tags-recursive \
uninstall uninstall-am uninstall-binPROGRAMS
clean-generic clean-libtool ctags ctags-recursive distclean \
distclean-compile distclean-generic distclean-libtool \
distclean-tags distdir dvi dvi-am html html-am info info-am \
install install-am install-binPROGRAMS install-data \
install-data-am install-dvi install-dvi-am install-exec \
install-exec-am install-html install-html-am install-info \
install-info-am install-man install-pdf install-pdf-am \
install-ps install-ps-am install-strip installcheck \
installcheck-am installdirs installdirs-am maintainer-clean \
maintainer-clean-generic mostlyclean mostlyclean-compile \
mostlyclean-generic mostlyclean-libtool pdf pdf-am ps ps-am \
tags tags-recursive uninstall uninstall-am \
uninstall-binPROGRAMS
 
 
# Tell versions [3.59,3.63) of GNU make to not export all variables.
# Otherwise a system limit (for SysV at least) may be exceeded.
.NOEXPORT:
/branches/rhl/src/sample_utils.c New file
0,0 → 1,471
/*
* utils.c
*
* Miscellaneous functions.
*
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
* This file part of: PSFEx
*
* Copyright: (C) 1997-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
*
* License: GNU General Public License
*
* PSFEx is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
* PSFEx is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
* You should have received a copy of the GNU General Public License
* along with PSFEx. If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 14/07/2013, Bastille Day
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
 
#include <assert.h>
#include <ctype.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
 
#include "define.h"
#include "globals.h"
#include "prefs.h"
 
/*****************************************************************************/
/*
* Copied here from sample.c as the rest of that file's all about fits I/O
*/
#include "sample.h"
 
/****** malloc_samples *******************************************************
PROTO void malloc_samples(setstruct *set, int nsample)
PURPOSE Allocate memory for a set of samples.
INPUT set structure pointer,
desired number of samples.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
VERSION 26/03/2008
*/
void malloc_samples(setstruct *set, int nsample)
 
{
samplestruct *sample;
int n;
 
QMALLOC(set->sample, samplestruct *, nsample);
set->nsamplemax = nsample;
if (set != set->samples_owner)
return;
 
QMALLOC(set->s_sample, samplestruct, nsample);
for (n=0; n!=nsample; ++n)
{
sample = set->sample[n] = &set->s_sample[n];
QMALLOC(sample->vig, float, set->nvig);
QMALLOC(sample->vigresi, float, set->nvig);
QMALLOC(sample->vigweight, float, set->nvig);
QMALLOC(sample->vigchi, float, set->nvig);
if (set->ncontext)
QMALLOC(sample->context, double, set->ncontext);
}
 
return;
}
 
 
/****** realloc_samples ******************************************************
PROTO void realloc_samples(setstruct *set, int nsample)
PURPOSE Re-allocate memory for a set of samples.
INPUT set structure pointer,
desired number of samples.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
VERSION 26/03/2008
*/
void realloc_samples(setstruct *set, int nsample)
 
{
samplestruct *sample;
int n;
 
if (set != set->samples_owner)
{
assert(nsample <= set->nsamplemax);
return;
}
 
/* If we want to reallocate 0 samples, better free the whole thing! */
if (!nsample)
free_samples(set);
 
/* Two cases: either more samples are required, or the opposite! */
if (nsample>set->nsamplemax)
{
QREALLOC(set->s_sample, samplestruct, nsample);
free(set->sample);
QMALLOC(set->sample, samplestruct *, nsample);
for (n = set->nsamplemax; n<nsample; n++)
{
sample = &set->s_sample[n];
QMALLOC(sample->vig, float, set->nvig);
QMALLOC(sample->vigresi, float, set->nvig);
QMALLOC(sample->vigchi, float, set->nvig);
QMALLOC(sample->vigweight, float, set->nvig);
if (set->ncontext)
QMALLOC(sample->context, double, set->ncontext);
}
}
else if (nsample<set->nsamplemax)
{
for (n = nsample; n<set->nsamplemax; n++)
{
sample = &set->s_sample[n];
free(sample->vig);
free(sample->vigresi);
free(sample->vigchi);
free(sample->vigweight);
if (set->ncontext)
free(sample->context);
}
QREALLOC(set->s_sample, samplestruct, nsample);
free(set->sample);
QMALLOC(set->sample, samplestruct *, nsample);
}
 
set->nsamplemax = nsample;
 
for (n = 0; n<nsample; n++)
set->sample[n] = &set->s_sample[n];
 
return;
}
 
 
/****** free_samples *********************************************************
PROTO void free_samples(setstruct *set, int nsample)
PURPOSE free memory for a set of samples.
INPUT set structure pointer,
desired number of samples.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
VERSION 26/03/2008
*/
void free_samples(setstruct *set)
 
{
samplestruct *sample;
int n;
 
if (set->samples_owner == set)
{
for (n = 0; n<set->nsamplemax; ++n )
{
sample = set->sample[n];
free(sample->vig);
free(sample->vigresi);
free(sample->vigweight);
free(sample->vigchi);
if (set->ncontext)
free(sample->context);
}
free(set->s_sample);
}
 
free(set->sample);
set->sample = NULL;
set->nsample = set->nsamplemax = 0;
 
return;
}
 
 
/****** remove_sample ********************************************************
PROTO samplestruct *remove_sample(setstruct *set, int isample)
PURPOSE Remove an element from a set of samples.
INPUT set structure pointer,
sample number.
OUTPUT The new pointer for the element that replaced the removed one.
NOTES -.
AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
VERSION 01/03/99
*/
samplestruct *remove_sample(setstruct *set, int isample)
 
{
static samplestruct exsample;
samplestruct *sample;
int nsample;
 
/* If we want to reallocate 0 samples, better free the whole thing! */
nsample = set->nsample-1;
if (nsample>0)
{
sample = set->sample[isample];
exsample = *set->sample[nsample];
*set->sample[nsample] = *sample;
*sample = exsample;
}
else
nsample=0;
realloc_samples(set, nsample);
set->nsample = nsample;
 
return set->sample[isample];
}
 
 
/****** init_set ************************************************************
PROTO setstruct *init_set()
PURPOSE Allocate and initialize a set structure.
INPUT -.
OUTPUT -.
NOTES See prefs.h.
AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
VERSION 26/03/2008
*/
setstruct *init_set(contextstruct *context)
 
{
setstruct *set;
int i;
 
QCALLOC(set, setstruct, 1);
set->nsample = set->nsamplemax = 0;
set->samples_owner = set; /* we'll own the samples, so we'll need to free them */
set->vigdim = 2;
QMALLOC(set->vigsize, int, set->vigdim);
set->ncontext = context->ncontext;
if (set->ncontext)
{
QMALLOC(set->contextoffset, double, set->ncontext);
QMALLOC(set->contextscale, double, set->ncontext);
QMALLOC(set->contextname, char *, set->ncontext);
for (i=0; i<set->ncontext; i++)
QMALLOC(set->contextname[i], char, 80);
}
 
return set;
}
 
 
/****** end_set *************************************************************
PROTO void end_set(setstruct *set)
PURPOSE free memory allocated by a complete set structure.
INPUT set structure pointer,
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
VERSION 26/03/2008
*/
void end_set(setstruct *set)
 
{
int i;
 
free_samples(set);
free(set->vigsize);
if (set->ncontext)
{
for (i=0; i<set->ncontext; i++)
free(set->contextname[i]);
free(set->contextname);
free(set->contextoffset);
free(set->contextscale);
}
free(set->head);
free(set);
 
return;
}
 
 
/****** make_weights *********************************************************
PROTO void make_weights(setstruct *set, samplestruct *sample)
PURPOSE Produce a weight-map for each sample vignet.
INPUT set structure pointer,
sample structure pointer.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP,Leiden observatory & ESO)
VERSION 13/08/2007
*/
void make_weights(setstruct *set, samplestruct *sample)
 
{
float *vig, *vigweight,
backnoise2, gain, noise2, profaccu2, pix;
int i;
 
/* Produce a weight-map */
profaccu2 = prefs.prof_accuracy*prefs.prof_accuracy;
gain = sample->gain;
backnoise2 = sample->backnoise2;
for (vig=sample->vig, vigweight=sample->vigweight, i=set->nvig; i--;)
{
if (*vig <= -BIG)
*(vig++) = *(vigweight++) = 0.0;
else
{
pix = *(vig++);
noise2 = backnoise2 + profaccu2*pix*pix;
if (pix>0.0 && gain>0.0)
noise2 += pix/gain;
*(vigweight++) = 1.0/noise2;
}
}
 
return;
}
 
 
/****** recenter_sample ******************************************************
PROTO void recenter_samples(samplestruct sample,
setstruct *set, float fluxrad)
PURPOSE Recenter sample image using windowed barycenter.
INPUT sample structure pointer,
set structure pointer,
flux radius.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 16/11/2009
*/
void recenter_sample(samplestruct *sample, setstruct *set, float fluxrad)
 
{
double tv, dxpos, dypos;
float *ima,*imat,*weight,*weightt,
pix, var, locpix, locarea, sig,twosig2, raper,raper2,
offsetx,offsety, mx,my, mx2ph,my2ph,
rintlim,rintlim2,rextlim2, scalex,scaley,scale2,
dx,dy, dx1,dy2, r2;
int i, x,y,x2,y2, w,h, sx,sy, xmin,xmax,ymin,ymax, pos;
 
sig = fluxrad*2.0/2.35; /* From half-FWHM to sigma */
twosig2 = 2.0*sig*sig;
 
w = set->vigsize[0];
h = set->vigsize[1];
/* Integration radius */
raper = RECENTER_NSIG*sig;
raper2 = raper*raper;
/* Internal radius of the oversampled annulus (<r-sqrt(2)/2) */
rintlim = raper - 0.75;
rintlim2 = (rintlim>0.0)? rintlim*rintlim: 0.0;
/* External radius of the oversampled annulus (>r+sqrt(2)/2) */
rextlim2 = (raper + 0.75)*(raper + 0.75);
scaley = scalex = 1.0/RECENTER_OVERSAMP;
scale2 = scalex*scaley;
offsetx = 0.5*(scalex-1.0);
offsety = 0.5*(scaley-1.0);
/* Use isophotal centroid as a first guess */
mx = sample->dx + (float)(w/2);
my = sample->dy + (float)(h/2);
 
for (i=0; i<RECENTER_NITERMAX; i++)
{
xmin = (int)(mx-raper+0.499999);
xmax = (int)(mx+raper+1.499999);
ymin = (int)(my-raper+0.499999);
ymax = (int)(my+raper+1.499999);
mx2ph = mx*2.0 + 0.49999;
my2ph = my*2.0 + 0.49999;
 
if (xmin < 0)
xmin = 0;
if (xmax > w)
xmax = w;
if (ymin < 0)
ymin = 0;
if (ymax > h)
ymax = h;
tv = 0.0;
dxpos = dypos = 0.0;
ima = sample->vig;
weight = sample->vigweight;
for (y=ymin; y<ymax; y++)
{
imat= ima + (pos = y*w + xmin);
weightt = weight + pos;
dy = y - my;
for (x=xmin; x<xmax; x++, imat++, weightt++)
{
dx = x - mx;
if ((r2=dx*dx+dy*dy)<rextlim2)
{
if (RECENTER_OVERSAMP>1 && r2> rintlim2)
{
dx += offsetx;
dy += offsety;
locarea = 0.0;
for (sy=RECENTER_OVERSAMP; sy--; dy+=scaley)
{
dx1 = dx;
dy2 = dy*dy;
for (sx=RECENTER_OVERSAMP; sx--; dx1+=scalex)
if (dx1*dx1+dy2<raper2)
locarea += scale2;
}
}
else
locarea = 1.0;
locarea *= expf(-r2/twosig2);
/*-------- Here begin tests for pixel and/or weight overflows. Things are a */
/*-------- bit intricated to have it running as fast as possible in the most */
/*-------- common cases */
pix = *imat;
if ((var=*weightt)==0.0)
{
if ((x2=(int)(mx2ph-x))>=0 && x2<w
&& (y2=(int)(my2ph-y))>=0 && y2<h
&& (var=*(weight + (pos = y2*w + x2)))>0.0)
{
var = 1.0/var;
pix = *(ima+pos);
}
else
pix = var = 0.0;
}
dx = x - mx;
dy = y - my;
locpix = locarea*pix;
tv += locpix;
dxpos += locpix*dx;
dypos += locpix*dy;
}
}
}
 
if (tv>0.0)
{
mx += (dxpos /= tv)*RECENTER_GRADFAC;
my += (dypos /= tv)*RECENTER_GRADFAC;
}
else
break;
 
/*-- Stop here if position does not change */
if (dxpos*dxpos+dypos*dypos < RECENTER_STEPMIN*RECENTER_STEPMIN)
break;
}
 
sample->dx = mx - (float)(w/2);
sample->dy = my - (float)(h/2);
 
return;
}
/branches/rhl/src/sample.c
796,426 → 796,3
 
return set;
}
 
 
/****** recenter_sample ******************************************************
PROTO void recenter_samples(samplestruct sample,
setstruct *set, float fluxrad)
PURPOSE Recenter sample image using windowed barycenter.
INPUT sample structure pointer,
set structure pointer,
flux radius.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 16/11/2009
*/
void recenter_sample(samplestruct *sample, setstruct *set, float fluxrad)
 
{
double tv, dxpos, dypos;
float *ima,*imat,*weight,*weightt,
pix, var, locpix, locarea, sig,twosig2, raper,raper2,
offsetx,offsety, mx,my, mx2ph,my2ph,
rintlim,rintlim2,rextlim2, scalex,scaley,scale2,
dx,dy, dx1,dy2, r2;
int i, x,y,x2,y2, w,h, sx,sy, xmin,xmax,ymin,ymax, pos;
 
sig = fluxrad*2.0/2.35; /* From half-FWHM to sigma */
twosig2 = 2.0*sig*sig;
 
w = set->vigsize[0];
h = set->vigsize[1];
/* Integration radius */
raper = RECENTER_NSIG*sig;
raper2 = raper*raper;
/* Internal radius of the oversampled annulus (<r-sqrt(2)/2) */
rintlim = raper - 0.75;
rintlim2 = (rintlim>0.0)? rintlim*rintlim: 0.0;
/* External radius of the oversampled annulus (>r+sqrt(2)/2) */
rextlim2 = (raper + 0.75)*(raper + 0.75);
scaley = scalex = 1.0/RECENTER_OVERSAMP;
scale2 = scalex*scaley;
offsetx = 0.5*(scalex-1.0);
offsety = 0.5*(scaley-1.0);
/* Use isophotal centroid as a first guess */
mx = sample->dx + (float)(w/2);
my = sample->dy + (float)(h/2);
 
for (i=0; i<RECENTER_NITERMAX; i++)
{
xmin = (int)(mx-raper+0.499999);
xmax = (int)(mx+raper+1.499999);
ymin = (int)(my-raper+0.499999);
ymax = (int)(my+raper+1.499999);
mx2ph = mx*2.0 + 0.49999;
my2ph = my*2.0 + 0.49999;
 
if (xmin < 0)
xmin = 0;
if (xmax > w)
xmax = w;
if (ymin < 0)
ymin = 0;
if (ymax > h)
ymax = h;
tv = 0.0;
dxpos = dypos = 0.0;
ima = sample->vig;
weight = sample->vigweight;
for (y=ymin; y<ymax; y++)
{
imat= ima + (pos = y*w + xmin);
weightt = weight + pos;
dy = y - my;
for (x=xmin; x<xmax; x++, imat++, weightt++)
{
dx = x - mx;
if ((r2=dx*dx+dy*dy)<rextlim2)
{
if (RECENTER_OVERSAMP>1 && r2> rintlim2)
{
dx += offsetx;
dy += offsety;
locarea = 0.0;
for (sy=RECENTER_OVERSAMP; sy--; dy+=scaley)
{
dx1 = dx;
dy2 = dy*dy;
for (sx=RECENTER_OVERSAMP; sx--; dx1+=scalex)
if (dx1*dx1+dy2<raper2)
locarea += scale2;
}
}
else
locarea = 1.0;
locarea *= expf(-r2/twosig2);
/*-------- Here begin tests for pixel and/or weight overflows. Things are a */
/*-------- bit intricated to have it running as fast as possible in the most */
/*-------- common cases */
pix = *imat;
if ((var=*weightt)==0.0)
{
if ((x2=(int)(mx2ph-x))>=0 && x2<w
&& (y2=(int)(my2ph-y))>=0 && y2<h
&& (var=*(weight + (pos = y2*w + x2)))>0.0)
{
var = 1.0/var;
pix = *(ima+pos);
}
else
pix = var = 0.0;
}
dx = x - mx;
dy = y - my;
locpix = locarea*pix;
tv += locpix;
dxpos += locpix*dx;
dypos += locpix*dy;
}
}
}
 
if (tv>0.0)
{
mx += (dxpos /= tv)*RECENTER_GRADFAC;
my += (dypos /= tv)*RECENTER_GRADFAC;
}
else
break;
 
/*-- Stop here if position does not change */
if (dxpos*dxpos+dypos*dypos < RECENTER_STEPMIN*RECENTER_STEPMIN)
break;
}
 
sample->dx = mx - (float)(w/2);
sample->dy = my - (float)(h/2);
 
return;
}
 
 
/****** malloc_samples *******************************************************
PROTO void malloc_samples(setstruct *set, int nsample)
PURPOSE Allocate memory for a set of samples.
INPUT set structure pointer,
desired number of samples.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
VERSION 26/03/2008
*/
void malloc_samples(setstruct *set, int nsample)
 
{
samplestruct *sample;
int n;
 
QMALLOC(set->sample, samplestruct *, nsample);
set->nsamplemax = nsample;
if (set != set->samples_owner)
return;
 
QMALLOC(set->s_sample, samplestruct, nsample);
for (n=0; n!=nsample; ++n)
{
sample = set->sample[n] = &set->s_sample[n];
QMALLOC(sample->vig, float, set->nvig);
QMALLOC(sample->vigresi, float, set->nvig);
QMALLOC(sample->vigweight, float, set->nvig);
QMALLOC(sample->vigchi, float, set->nvig);
if (set->ncontext)
QMALLOC(sample->context, double, set->ncontext);
}
 
return;
}
 
 
/****** realloc_samples ******************************************************
PROTO void realloc_samples(setstruct *set, int nsample)
PURPOSE Re-allocate memory for a set of samples.
INPUT set structure pointer,
desired number of samples.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
VERSION 26/03/2008
*/
void realloc_samples(setstruct *set, int nsample)
 
{
samplestruct *sample;
int n;
 
if (set != set->samples_owner)
{
assert(nsample <= set->nsamplemax);
return;
}
 
/* If we want to reallocate 0 samples, better free the whole thing! */
if (!nsample)
free_samples(set);
 
/* Two cases: either more samples are required, or the opposite! */
if (nsample>set->nsamplemax)
{
QREALLOC(set->s_sample, samplestruct, nsample);
free(set->sample);
QMALLOC(set->sample, samplestruct *, nsample);
for (n = set->nsamplemax; n<nsample; n++)
{
sample = &set->s_sample[n];
QMALLOC(sample->vig, float, set->nvig);
QMALLOC(sample->vigresi, float, set->nvig);
QMALLOC(sample->vigchi, float, set->nvig);
QMALLOC(sample->vigweight, float, set->nvig);
if (set->ncontext)
QMALLOC(sample->context, double, set->ncontext);
}
}
else if (nsample<set->nsamplemax)
{
for (n = nsample; n<set->nsamplemax; n++)
{
sample = &set->s_sample[n];
free(sample->vig);
free(sample->vigresi);
free(sample->vigchi);
free(sample->vigweight);
if (set->ncontext)
free(sample->context);
}
QREALLOC(set->s_sample, samplestruct, nsample);
free(set->sample);
QMALLOC(set->sample, samplestruct *, nsample);
}
 
set->nsamplemax = nsample;
 
for (n = 0; n<nsample; n++)
set->sample[n] = &set->s_sample[n];
 
return;
}
 
 
/****** free_samples *********************************************************
PROTO void free_samples(setstruct *set, int nsample)
PURPOSE free memory for a set of samples.
INPUT set structure pointer,
desired number of samples.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
VERSION 26/03/2008
*/
void free_samples(setstruct *set)
 
{
samplestruct *sample;
int n;
 
if (set->samples_owner == set)
{
for (n = set->nsamplemax; n--; sample++)
{
sample = set->sample[n];
free(sample->vig);
free(sample->vigresi);
free(sample->vigweight);
free(sample->vigchi);
if (set->ncontext)
free(sample->context);
}
free(set->s_sample);
}
 
free(set->sample);
set->sample = NULL;
set->nsample = set->nsamplemax = 0;
 
return;
}
 
 
/****** remove_sample ********************************************************
PROTO samplestruct *remove_sample(setstruct *set, int isample)
PURPOSE Remove an element from a set of samples.
INPUT set structure pointer,
sample number.
OUTPUT The new pointer for the element that replaced the removed one.
NOTES -.
AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
VERSION 01/03/99
*/
samplestruct *remove_sample(setstruct *set, int isample)
 
{
static samplestruct exsample;
samplestruct *sample;
int nsample;
 
/* If we want to reallocate 0 samples, better free the whole thing! */
nsample = set->nsample-1;
if (nsample>0)
{
sample = set->sample[isample];
exsample = *set->sample[nsample];
*set->sample[nsample] = *sample;
*sample = exsample;
}
else
nsample=0;
realloc_samples(set, nsample);
set->nsample = nsample;
 
return set->sample[isample];
}
 
 
/****** init_set ************************************************************
PROTO setstruct *init_set()
PURPOSE Allocate and initialize a set structure.
INPUT -.
OUTPUT -.
NOTES See prefs.h.
AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
VERSION 26/03/2008
*/
setstruct *init_set(contextstruct *context)
 
{
setstruct *set;
int i;
 
QCALLOC(set, setstruct, 1);
set->nsample = set->nsamplemax = 0;
set->samples_owner = set; /* we'll own the samples, so we'll need to free them */
set->vigdim = 2;
QMALLOC(set->vigsize, int, set->vigdim);
set->ncontext = context->ncontext;
if (set->ncontext)
{
QMALLOC(set->contextoffset, double, set->ncontext);
QMALLOC(set->contextscale, double, set->ncontext);
QMALLOC(set->contextname, char *, set->ncontext);
for (i=0; i<set->ncontext; i++)
QMALLOC(set->contextname[i], char, 80);
}
 
return set;
}
 
 
/****** end_set *************************************************************
PROTO void end_set(setstruct *set)
PURPOSE free memory allocated by a complete set structure.
INPUT set structure pointer,
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
VERSION 26/03/2008
*/
void end_set(setstruct *set)
 
{
int i;
 
free_samples(set);
free(set->vigsize);
if (set->ncontext)
{
for (i=0; i<set->ncontext; i++)
free(set->contextname[i]);
free(set->contextname);
free(set->contextoffset);
free(set->contextscale);
}
free(set->head);
free(set);
 
return;
}
 
 
/****** make_weights *********************************************************
PROTO void make_weights(setstruct *set, samplestruct *sample)
PURPOSE Produce a weight-map for each sample vignet.
INPUT set structure pointer,
sample structure pointer.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP,Leiden observatory & ESO)
VERSION 13/08/2007
*/
void make_weights(setstruct *set, samplestruct *sample)
 
{
float *vig, *vigweight,
backnoise2, gain, noise2, profaccu2, pix;
int i;
 
/* Produce a weight-map */
profaccu2 = prefs.prof_accuracy*prefs.prof_accuracy;
gain = sample->gain;
backnoise2 = sample->backnoise2;
for (vig=sample->vig, vigweight=sample->vigweight, i=set->nvig; i--;)
{
if (*vig <= -BIG)
*(vig++) = *(vigweight++) = 0.0;
else
{
pix = *(vig++);
noise2 = backnoise2 + profaccu2*pix*pix;
if (pix>0.0 && gain>0.0)
noise2 += pix/gain;
*(vigweight++) = 1.0/noise2;
}
}
 
return;
}
 
/branches/rhl/src/field_utils.c New file
0,0 → 1,81
/*
* field.c
*
* Manage multiple PSFs.
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
* This file part of: PSFEx
*
* Copyright: (C) 2007-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
*
* License: GNU General Public License
*
* PSFEx is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
* PSFEx is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
* You should have received a copy of the GNU General Public License
* along with PSFEx. If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 25/06/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
 
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
 
#include "define.h"
#include "types.h"
#include "globals.h"
#include "fits/fitscat.h"
#include "check.h"
#include "fitswcs.h"
#include "misc.h"
#include "prefs.h"
#include "psf.h"
#include "field.h"
 
 
/****** field_init_finalize ************************************************************
PROTO void field_init(fieldstruct *)
PURPOSE Finish allocating and initializing a PSF MEF structure (groups of PSFs).
INPUT fieldstruct pointer
OUTPUT -.
NOTES .
AUTHOR E. Bertin (IAP)
VERSION 08/04/2010
***/
void
field_init_finalize(fieldstruct *field)
{
int countsize = prefs.context_nsnap*prefs.context_nsnap;
int e, next0 = field->next;
 
field_locate(field);
QCALLOC(field->ccat, catstruct *, MAXCHECK);
countsize = prefs.context_nsnap*prefs.context_nsnap;
QMALLOC(field->lcount, int *, next0);
QMALLOC(field->acount, int *, next0);
QMALLOC(field->count, int *, next0);
QMALLOC(field->modchi2, double *, next0);
QMALLOC(field->modresi, double *, next0);
for (e=0; e<next0; e++)
{
QCALLOC(field->lcount[e], int, countsize);
QCALLOC(field->acount[e], int, countsize);
QCALLOC(field->count[e], int, countsize);
QCALLOC(field->modchi2[e], double, countsize);
QCALLOC(field->modresi[e], double, countsize);
}
}
/branches/rhl/src/dummies.c New file
0,0 → 1,235
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include "define.h"
#include "sample.h"
 
struct cat; typedef struct cat catstruct;
struct field; typedef struct field fieldstruct;
struct key; typedef struct key keystruct;
struct tab; typedef struct tab tabstruct;
struct wcs; typedef struct wcs wcsstruct;
 
typedef enum {H_TYPE} h_type;
typedef enum {T_TYPE} t_type;
typedef enum {ACCESS_TYPE_T} access_type_t;
typedef float PIXTYPE;
 
const int t_size[] = {0};
time_t thetime, thetime2;
 
int
add_key(keystruct *key, tabstruct *tab, int pos)
{
abort();
return -1;
}
 
int
addkeywordto_head(tabstruct *tab, char *keyword, char *comment)
{
abort();
return -1;
}
 
int
blank_keys(tabstruct *tab)
{
abort();
return -1;
}
 
void
end_wcs(wcsstruct *wcs)
{
abort();
}
 
int
fitsread(char *fitsbuf, char *keyword, void *ptr, h_type htype,
t_type ttype)
{
abort();
return -1;
}
 
int
fitswrite(char *fitsbuf, char *keyword, void *ptr, h_type htype,
t_type ttype)
{
abort();
return -1;
}
 
int
init_cat(catstruct *cat)
{
abort();
return -1;
}
 
catstruct *
read_cat(char *filename)
{
abort();
return NULL;
}
 
keystruct *
read_key(tabstruct *tab, char *keyname)
{
abort();
return NULL;
}
 
keystruct *
new_key(char *keyname)
{
abort();
return NULL;
}
 
tabstruct *
new_tab(char *tabname)
{
abort();
return NULL;
}
 
wcsstruct *
read_wcs(tabstruct *tab)
{
abort();
return NULL;
}
 
double
wcs_dist(wcsstruct *wcs, double *wcspos1, double *wcspos2)
{
abort();
return -1.0;
}
 
double
wcs_scale(wcsstruct *wcs, double *pixpos)
{
abort();
return -1.0;
}
void
readbasic_head(tabstruct *tab)
{
abort();
}
 
catstruct *
new_cat(int ncat)
{
abort();
return NULL;
}
 
void
free_cat(catstruct **cat, int ncat)
{
abort();
}
 
void
free_tab(tabstruct *tab)
{
abort();
}
 
int
open_cat(catstruct *cat, access_type_t at)
{
abort();
return -1;
}
 
int
prim_head(tabstruct *tab)
{
abort();
return -1;
}
 
void
read_body(tabstruct *tab, PIXTYPE *ptr, size_t size)
{
abort();
}
 
void
save_tab(catstruct *cat, tabstruct *tab)
{
abort();
}
 
/*****************************************************************************/
/*
* These are not dummies. They replace psfex routines
*/
typedef enum {CHECKENUM} checkenum;
 
void
check_write(fieldstruct *field, setstruct *set,
char *checkname, checkenum checktype,
int ext, int next, int cubeflag)
{
;
}
 
setstruct *
load_samples(char **filenames, int catindex, int ncat, int ext,
int next, contextstruct *context)
{
/*
* The C version of this is called two ways:
* catindex == 0, ncat == ncat Read all catalogues
* catindex == c, ncat == 1 Read only catalogue c
*/
setstruct *completeSet = (setstruct *)(filenames[catindex + 0]);
/*
* Make a new set, which may be a subset of the completeSet
*/
setstruct *set = init_set(context);
set->fwhm = completeSet->fwhm;
for (int i = 0; i != completeSet->vigdim; ++i) {
set->vigsize[i] = completeSet->vigsize[i];
}
for (int i = 0; i != completeSet->ncontext; ++i) {
strcpy(set->contextname[i], completeSet->contextname[i]);
set->contextoffset[i] = completeSet->contextoffset[i];
set->contextscale[i] = completeSet->contextscale[i];
}
/*
* Count how many samples we'll be including
*/
int nsample_keep = 0;
for (int i = 0; i != ncat; ++i) {
setstruct *s = (setstruct *)(filenames[catindex + i]);
for (int j = 0; j != completeSet->nsample; ++j) {
samplestruct const *samp = s->sample[j];
if (ext == ALL_EXTENSIONS || ext == samp->extindex) {
++nsample_keep;
}
}
}
 
set->samples_owner = 0;
malloc_samples(set, nsample_keep);
for (int i = 0; i != ncat; ++i) {
setstruct *s = (setstruct *)(filenames[catindex + i]);
for (int j = 0; j != completeSet->nsample; ++j) {
samplestruct *samp = s->sample[j];
if (ext == ALL_EXTENSIONS || ext == samp->extindex) {
set->sample[set->nsample++] = samp;
}
}
}
 
return set;
}
/branches/rhl/src/makeit.c
40,7 → 40,6
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
 
#include "define.h"
#include "types.h"
58,10 → 57,8
#include "sample.h"
#include "xml.h"
 
psfstruct *make_psf(setstruct *set, float psfstep,
float *basis, int nbasis, contextstruct *context);
void write_error(char *msg1, char *msg2);
time_t thetime, thetime2;
void write_error(char *msg1, char *msg2);
 
/********************************** makeit ***********************************/
/*
147,361 → 144,20
}
QPRINTF(OUTPUT, "\n");
 
next = fields[0]->next;
 
if (prefs.xml_flag)
init_xml(ncat);
makeit_body(fields, &context, &fullcontext, 1);
next = fields[0]->next;
 
psfstep = prefs.psf_step;
psfsteps = NULL;
nbasis = 0;
psfbasis = NULL;
psfbasiss = NULL;
 
/* Initialize context */
NFPRINTF(OUTPUT, "Initializing contexts...");
context = context_init(prefs.context_name, prefs.context_group,
prefs.ncontext_group, prefs.group_deg, prefs.ngroup_deg,
CONTEXT_REMOVEHIDDEN);
fullcontext = context->npc?
context_init(prefs.context_name, prefs.context_group,
prefs.ncontext_group, prefs.group_deg, prefs.ngroup_deg,
CONTEXT_KEEPHIDDEN)
: context;
 
if (context->npc && ncat<2)
warning("Hidden dependencies cannot be derived from",
" a single catalog");
else if (context->npc && prefs.stability_type == STABILITY_EXPOSURE)
warning("Hidden dependencies have no effect",
" in STABILITY_TYPE EXPOSURE mode");
 
/* Compute PSF steps */
if (!prefs.psf_step)
/* Write XML */
if (prefs.xml_flag)
{
NFPRINTF(OUTPUT, "Computing optimum PSF sampling steps...");
if (prefs.newbasis_type==NEWBASIS_PCACOMMON
|| (prefs.stability_type == STABILITY_SEQUENCE
&& prefs.psf_mef_type == PSF_MEF_COMMON))
{
set = load_samples(incatnames, 0, ncat, ALL_EXTENSIONS, next, context);
psfstep = (float)((set->fwhm/2.35)*0.5);
end_set(set);
}
/*-- Need to derive a common pixel step for each ext */
else if (prefs.newbasis_type == NEWBASIS_PCAINDEPENDENT
|| context->npc
|| (prefs.stability_type == STABILITY_SEQUENCE
&& prefs.psf_mef_type == PSF_MEF_INDEPENDENT))
{
/*-- Run through all samples to derive a different pixel step for each extension */
QMALLOC(psfsteps, float, next);
for (ext=0 ; ext<next; ext++)
{
set = load_samples(incatnames, 0, ncat, ext, next, context);
psfsteps[ext] = (float)(psfstep? psfstep : (set->fwhm/2.35)*0.5);
end_set(set);
}
}
NFPRINTF(OUTPUT, "Writing XML file...");
write_xml(prefs.xml_name);
end_xml();
}
 
/* Derive a new common PCA basis for all extensions */
if (prefs.newbasis_type==NEWBASIS_PCACOMMON)
{
QMALLOC(cpsf, psfstruct *, ncat*next);
for (ext=0 ; ext<next; ext++)
for (c=0; c<ncat; c++)
{
sprintf(str, "Computing new PCA image basis from %s...",
fields[c]->rtcatname);
NFPRINTF(OUTPUT, str);
set = load_samples(incatnames, c, 1, ext, next, context);
step = psfstep;
cpsf[c+ext*ncat] = make_psf(set, psfstep, NULL, 0, context);
end_set(set);
}
nbasis = prefs.newbasis_number;
psfbasis = pca_onsnaps(cpsf, ncat*next, nbasis);
for (i=0 ; i<ncat*next; i++)
psf_end(cpsf[i]);
free(cpsf);
}
/* Derive a new PCA basis for each extension */
else if (prefs.newbasis_type == NEWBASIS_PCAINDEPENDENT)
{
nbasis = prefs.newbasis_number;
QMALLOC(psfbasiss, float *, next);
for (ext=0; ext<next; ext++)
{
if (psfsteps)
step = psfsteps[ext];
else
step = psfstep;
QMALLOC(cpsf, psfstruct *, ncat);
for (c=0; c<ncat; c++)
{
sprintf(str, "Computing new PCA image basis from %s...",
fields[c]->rtcatname);
NFPRINTF(OUTPUT, str);
set = load_samples(incatnames, c, 1, ext, next, context);
cpsf[c] = make_psf(set, step, NULL, 0, context);
end_set(set);
}
psfbasiss[ext] = pca_onsnaps(cpsf, ncat, nbasis);
for (c=0 ; c<ncat; c++)
psf_end(cpsf[c]);
free(cpsf);
}
}
 
if (context->npc && prefs.hidden_mef_type == HIDDEN_MEF_COMMON)
/*-- Derive principal components of PSF variation from the whole mosaic */
{
p = 0;
QMALLOC(cpsf, psfstruct *, ncat*next);
for (c=0; c<ncat; c++)
{
sprintf(str, "Computing hidden dependency parameter(s) from %s...",
fields[c]->rtcatname);
NFPRINTF(OUTPUT, str);
for (ext=0 ; ext<next; ext++)
{
set = load_samples(incatnames, c, 1, ext, next, context);
if (psfsteps)
step = psfsteps[ext];
else
step = psfstep;
basis = psfbasiss? psfbasiss[ext] : psfbasis;
cpsf[p++] = make_psf(set, step, basis, nbasis, context);
end_set(set);
}
}
free(fullcontext->pc);
fullcontext->pc = pca_oncomps(cpsf, next, ncat, context->npc);
for (c=0 ; c<ncat*next; c++)
psf_end(cpsf[c]);
free(cpsf);
}
 
/* Compute "final" PSF models */
if (prefs.psf_mef_type == PSF_MEF_COMMON)
{
if (prefs.stability_type == STABILITY_SEQUENCE)
{
/*---- Load all the samples at once */
set = load_samples(incatnames, 0, ncat, ALL_EXTENSIONS, next, context);
step = psfstep;
basis = psfbasis;
field_count(fields, set, COUNT_LOADED);
psf = make_psf(set, step, basis, nbasis, context);
field_count(fields, set, COUNT_ACCEPTED);
end_set(set);
NFPRINTF(OUTPUT, "Computing final PSF model...");
context_apply(context, psf, fields, ALL_EXTENSIONS, 0, ncat);
psf_end(psf);
}
else
for (c=0; c<ncat; c++)
{
/*------ Load the samples for current exposure */
sprintf(str, "Computing final PSF model from %s...",
fields[c]->rtcatname);
NFPRINTF(OUTPUT, str);
set = load_samples(incatnames, c, 1, ALL_EXTENSIONS, next, context);
if (psfstep)
step = psfstep;
else
step = (float)((set->fwhm/2.35)*0.5);
basis = psfbasis;
field_count(fields, set, COUNT_LOADED);
psf = make_psf(set, step, basis, nbasis, fullcontext);
field_count(fields, set, COUNT_ACCEPTED);
end_set(set);
context_apply(fullcontext, psf, fields, ALL_EXTENSIONS, c, 1);
psf_end(psf);
}
}
else
for (ext=0 ; ext<next; ext++)
{
basis = psfbasiss? psfbasiss[ext] : psfbasis;
if (context->npc && prefs.hidden_mef_type == HIDDEN_MEF_INDEPENDENT)
/*------ Derive principal components of PSF components */
{
QMALLOC(cpsf, psfstruct *, ncat);
if (psfsteps)
step = psfsteps[ext];
else
step = psfstep;
for (c=0; c<ncat; c++)
{
if (next>1)
sprintf(str,
"Computing hidden dependency parameter(s) from %s[%d/%d]...",
fields[c]->rtcatname, ext+1, next);
else
sprintf(str, "Computing hidden dependency parameter(s) from %s...",
fields[c]->rtcatname);
NFPRINTF(OUTPUT, str);
set = load_samples(incatnames, c, 1, ext, next, context);
field_count(fields, set, COUNT_LOADED);
cpsf[c] = make_psf(set, step, basis, nbasis, context);
field_count(fields, set, COUNT_ACCEPTED);
end_set(set);
}
free(fullcontext->pc);
fullcontext->pc = pca_oncomps(cpsf, 1, ncat, context->npc);
for (c=0 ; c<ncat; c++)
psf_end(cpsf[c]);
free(cpsf);
}
 
if (prefs.stability_type == STABILITY_SEQUENCE)
{
/*------ Load all the samples at once */
if (next>1)
{
sprintf(str, "Computing final PSF model for extension [%d/%d]...",
ext+1, next);
NFPRINTF(OUTPUT, str);
}
else
NFPRINTF(OUTPUT, "Computing final PSF model...");
set = load_samples(incatnames, 0, ncat, ext, next, fullcontext);
if (psfstep)
step = psfstep;
else if (psfsteps)
step = psfsteps[ext];
else
step = (float)((set->fwhm/2.35)*0.5);
field_count(fields, set, COUNT_LOADED);
psf = make_psf(set, step, basis, nbasis, fullcontext);
field_count(fields, set, COUNT_ACCEPTED);
end_set(set);
context_apply(fullcontext, psf, fields, ext, 0, ncat);
psf_end(psf);
}
else
for (c=0; c<ncat; c++)
{
/*-------- Load the samples for current exposure */
if (next>1)
sprintf(str, "Reading data from %s[%d/%d]...",
fields[c]->rtcatname, ext+1, next);
else
sprintf(str, "Reading data from %s...",
fields[c]->rtcatname);
NFPRINTF(OUTPUT, str);
set = load_samples(incatnames, c, 1, ext, next, context);
if (psfstep)
step = psfstep;
else if (psfsteps)
step = psfsteps[ext];
else
step = (float)((set->fwhm/2.35)*0.5);
if (next>1)
sprintf(str, "Computing final PSF model for %s[%d/%d]...",
fields[c]->rtcatname, ext+1, next);
else
sprintf(str, "Computing final PSF model for %s...",
fields[c]->rtcatname);
NFPRINTF(OUTPUT, str);
field_count(fields, set, COUNT_LOADED);
psf = make_psf(set, step, basis, nbasis, context);
field_count(fields, set, COUNT_ACCEPTED);
end_set(set);
context_apply(context, psf, fields, ext, c, 1);
psf_end(psf);
}
}
 
free(psfsteps);
if (psfbasiss)
{
for (ext=0; ext<next; ext++)
free(psfbasiss[ext]);
free(psfbasiss);
}
else if (psfbasis)
free(psfbasis);
 
/* Compute diagnostics and check-images */
#ifdef USE_THREADS
/* Force MKL using a single thread as diagnostic code is already multithreaded*/
#ifdef HAVE_MKL
if (prefs.context_nsnap>2)
mkl_set_num_threads(1);
#endif
#endif
QIPRINTF(OUTPUT,
" filename [ext] accepted/total samp. chi2/dof FWHM ellip."
" resi. asym.");
for (c=0; c<ncat; c++)
{
field = fields[c];
for (ext=0 ; ext<next; ext++)
{
psf = field->psf[ext];
wcs = field->wcs[ext];
if (next>1)
sprintf(str, "Computing diagnostics for %s[%d/%d]...",
field->rtcatname, ext+1, next);
else
sprintf(str, "Computing diagnostics for %s...",
field->rtcatname);
NFPRINTF(OUTPUT, str);
/*---- Check PSF with individual datasets */
set2 = load_samples(incatnames, c, 1, ext, next, context);
psf->samples_loaded = set2->nsample;
if (set2->nsample>1)
{
/*------ Remove bad PSF candidates */
psf_clean(psf, set2, prefs.prof_accuracy);
psf->chi2 = set2->nsample? psf_chi2(psf, set2) : 0.0;
}
psf->samples_accepted = set2->nsample;
/*---- Compute diagnostics and field statistics */
psf_diagnostic(psf);
psf_wcsdiagnostic(psf, wcs);
nmed = psf->nmed;
field_stats(fields, set2);
/*---- Display stats for current catalog/extension */
if (next>1)
sprintf(str, "[%d/%d]", ext+1, next);
else
str[0] = '\0';
QPRINTF(OUTPUT, "%-17.17s%-7.7s %5d/%-5d %6.2f %6.2f %6.2f %4.2f"
" %5.2f %5.2f\n",
ext==0? field->rtcatname : "",
str,
psf->samples_accepted, psf->samples_loaded,
psf->pixstep,
psf->chi2,
psf->moffat_fwhm,
psf->moffat_ellipticity,
psf->pfmoffat_residuals,
psf->sym_residuals);
/*---- Save "Check-images" */
for (i=0; i<prefs.ncheck_type; i++)
if (prefs.check_type[i])
{
sprintf(str, "Saving CHECK-image #%d...", i+1);
NFPRINTF(OUTPUT, str);
check_write(field, set2, prefs.check_name[i], prefs.check_type[i],
ext, next, prefs.check_cubeflag);
}
/*---- Free memory */
end_set(set2);
}
}
 
#ifdef USE_THREADS
/* Back to multithreaded MKL */
#ifdef HAVE_MKL
mkl_set_num_threads(prefs.nthreads);
#endif
#endif
 
/* Save result */
for (c=0; c<ncat; c++)
{
563,6 → 219,7
cplot_modchi2(fields[c]);
cplot_modresi(fields[c]);
#endif
 
/*-- Update XML */
if (prefs.xml_flag)
update_xml(fields[c]);
577,14 → 234,6
tm->tm_hour, tm->tm_min, tm->tm_sec);
prefs.time_diff = difftime(thetime2, thetime);
 
/* Write XML */
if (prefs.xml_flag)
{
NFPRINTF(OUTPUT, "Writing XML file...");
write_xml(prefs.xml_name);
end_xml();
}
 
/* Free memory */
for (c=0; c<ncat; c++)
field_end(fields[c]);
598,95 → 247,8
}
 
 
/****** make_psf *************************************************************
PROTO psfstruct *make_psf(setstruct *set, float psfstep,
float *basis, int nbasis, contextstruct *context)
PURPOSE Make PSFs from a set of FITS binary catalogs.
INPUT Pointer to a sample set,
PSF sampling step,
Pointer to basis image vectors,
Number of basis vectors,
Pointer to context structure.
OUTPUT Pointer to the PSF structure.
NOTES Diagnostics are computed only if diagflag != 0.
AUTHOR E. Bertin (IAP)
VERSION 03/09/2009
***/
psfstruct *make_psf(setstruct *set, float psfstep,
float *basis, int nbasis, contextstruct *context)
{
psfstruct *psf;
basistypenum basistype;
float pixsize[2];
 
pixsize[0] = (float)prefs.psf_pixsize[0];
pixsize[1] = (float)prefs.psf_pixsize[1];
// NFPRINTF(OUTPUT,"Initializing PSF modules...");
psf = psf_init(context, prefs.psf_size, psfstep, pixsize, set->nsample);
 
psf->samples_loaded = set->nsample;
psf->fwhm = set->fwhm;
/* Make the basic PSF-model (1st pass) */
// NFPRINTF(OUTPUT,"Modeling the PSF (1/3)...");
psf_make(psf, set, 0.2);
if (basis && nbasis)
{
QMEMCPY(basis, psf->basis, float, nbasis*psf->size[0]*psf->size[1]);
psf->nbasis = nbasis;
}
else
{
// NFPRINTF(OUTPUT,"Generating the PSF model...");
basistype = prefs.basis_type;
if (basistype==BASIS_PIXEL_AUTO)
basistype = (psf->fwhm < PSF_AUTO_FWHM)? BASIS_PIXEL : BASIS_NONE;
psf_makebasis(psf, set, basistype, prefs.basis_number);
}
psf_refine(psf, set);
 
/* Remove bad PSF candidates */
if (set->nsample>1)
{
psf_clean(psf, set, 0.2);
/*-- Make the basic PSF-model (2nd pass) */
// NFPRINTF(OUTPUT,"Modeling the PSF (2/3)...");
psf_make(psf, set, 0.1);
psf_refine(psf, set);
}
/* Remove bad PSF candidates */
if (set->nsample>1)
{
psf_clean(psf, set, 0.1);
/*-- Make the basic PSF-model (3rd pass) */
// NFPRINTF(OUTPUT,"Modeling the PSF (3/3)...");
psf_make(psf, set, 0.05);
psf_refine(psf, set);
}
/* Remove bad PSF candidates */
if (set->nsample>1)
psf_clean(psf, set, 0.05);
psf->samples_accepted = set->nsample;
/* Refine the PSF-model */
psf_make(psf, set, prefs.prof_accuracy);
psf_refine(psf, set);
/* Clip the PSF-model */
psf_clip(psf);
 
/*-- Just check the Chi2 */
psf->chi2 = set->nsample? psf_chi2(psf, set) : 0.0;
 
return psf;
}
 
 
/****** write_error ********************************************************
PROTO void write_error(char *msg1, char *msg2)
PURPOSE Manage files in case of a catched error
708,5 → 270,3
 
return;
}
 
 
/branches/rhl/src/field.c
62,7 → 62,7
tabstruct *tab, *imatab;
keystruct *key;
char *pstr;
int e, next, next0, ntab, countsize;
int next, next0, ntab;
 
QCALLOC(field, fieldstruct, 1);
/* Compute the number of valid input extensions */
130,27 → 130,11
 
free_cat(&cat, 1);
 
field_locate(field);
QCALLOC(field->ccat, catstruct *, MAXCHECK);
countsize = prefs.context_nsnap*prefs.context_nsnap;
QMALLOC(field->lcount, int *, next0);
QMALLOC(field->acount, int *, next0);
QMALLOC(field->count, int *, next0);
QMALLOC(field->modchi2, double *, next0);
QMALLOC(field->modresi, double *, next0);
for (e=0; e<next0; e++)
{
QCALLOC(field->lcount[e], int, countsize);
QCALLOC(field->acount[e], int, countsize);
QCALLOC(field->count[e], int, countsize);
QCALLOC(field->modchi2[e], double, countsize);
QCALLOC(field->modresi[e], double, countsize);
}
field_init_finalize(field);
 
return field;
}
 
 
/****** field_end *************************************************************
PROTO void field_end(fieldstruct *field)
PURPOSE Free a PSF MEF structure (groups of PSFs).
/branches/rhl/src/fits/fitscat.h
218,7 → 218,6
encode_checksum(unsigned int sum, char *str),
end_readobj(tabstruct *keytab, tabstruct *tab, char *buf),
end_writeobj(catstruct *cat, tabstruct *tab, char *buf),
error(int, char *, char *),
error_installfunc(void (*func)(char *msg1, char *msg2)),
fixexponent(char *s),
free_body(tabstruct *tab),