profit.c

/*
*				profit.c
*
* Fit a range of galaxy models to an image.
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
*	This file part of:	SExtractor
*
*	Copyright:		(C) 2006-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
*
*	License:		GNU General Public License
*
*	SExtractor 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.
*	SExtractor 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 SExtractor. If not, see <http://www.gnu.org/licenses/>.
*
*	Last modified:		23/09/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/

#ifdef HAVE_CONFIG_H
#include        "config.h"
#endif

#ifndef HAVE_MATHIMF_H
#define _GNU_SOURCE
#endif

#include	<math.h>
#include	<stdio.h>
#include	<stdlib.h>
#include	<string.h>

#include	"define.h"
#include	"globals.h"
#include	"prefs.h"
#include	"fits/fitscat.h"
#include	"levmar/levmar.h"
#include	"fft.h"
#include	"fitswcs.h"
#include	"check.h"
#include	"image.h"
#include	"pattern.h"
#include	"psf.h"
#include	"profit.h"

static double	prof_gammainc(double x, double a),
		prof_gamma(double x);
static float	prof_interpolate(profstruct *prof, float *posin);
static float	interpolate_pix(float *posin, float *pix, int *naxisn,
		interpenum interptype);

static void	make_kernel(float pos, float *kernel, interpenum interptype);

/*------------------------------- variables ---------------------------------*/

const int	interp_kernwidth[5]={1,2,4,6,8};

const int	flux_flag[PARAM_NPARAM] = {0,
					1,0,0,
					1,0,0,0,0,
					1,0,0,0,
					1,0,0,0,0,0,0,0,
					1,0,0,
					1,0,0,
					1,0,0
					};

/* "Local" global variables for debugging purposes */
int theniter, the_gal;
static picstruct	*the_field, *the_wfield;
profitstruct		*theprofit,*thedprofit, *thepprofit, *theqprofit;

/****** profit_init ***********************************************************
PROTO	profitstruct profit_init(psfstruct *psf, unsigned int modeltype)
PURPOSE	Allocate and initialize a new profile-fitting structure.
INPUT	Pointer to PSF structure,
	Model type.
OUTPUT	A pointer to an allocated profit structure.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	22/04/2011
 ***/
profitstruct	*profit_init(psfstruct *psf, unsigned int modeltype)
  {
   profitstruct		*profit;
   int			t, nmodels;

  QCALLOC(profit, profitstruct, 1);
  profit->psf = psf;
  QMALLOC(profit->prof, profstruct *, MODEL_NMAX);
  nmodels = 0;
  for (t=1; t<(1<<MODEL_NMAX); t<<=1)
    if (modeltype&t)
      profit->prof[nmodels++] = prof_init(profit, t);
/* Allocate memory for the complete model */
  QMALLOC16(profit->modpix, float, PROFIT_MAXMODSIZE*PROFIT_MAXMODSIZE);
  QMALLOC16(profit->modpix2, float, PROFIT_MAXMODSIZE*PROFIT_MAXMODSIZE);
  QMALLOC16(profit->cmodpix, float, PROFIT_MAXMODSIZE*PROFIT_MAXMODSIZE);
  QMALLOC16(profit->psfpix, float, PROFIT_MAXMODSIZE*PROFIT_MAXMODSIZE);
  QMALLOC16(profit->objpix, PIXTYPE, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE);
  QMALLOC16(profit->objweight, PIXTYPE, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE);
  QMALLOC16(profit->lmodpix, PIXTYPE, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE);
  QMALLOC16(profit->lmodpix2, PIXTYPE, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE);
  QMALLOC16(profit->resi, float, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE);
  QMALLOC16(profit->covar, float, profit->nparam*profit->nparam);
  profit->nprof = nmodels;
  profit->fluxfac = 1.0;	/* Default */

  return profit;
  }  


/****** profit_end ************************************************************
PROTO	void prof_end(profstruct *prof)
PURPOSE	End (deallocate) a profile-fitting structure.
INPUT	Prof structure.
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	12/07/2012
 ***/
void	profit_end(profitstruct *profit)
  {
   int	p;

  for (p=0; p<profit->nprof; p++)
    prof_end(profit->prof[p]);
  free(profit->modpix);
  free(profit->modpix2);
  free(profit->cmodpix);
  free(profit->psfpix);
  free(profit->lmodpix);
  free(profit->lmodpix2);
  free(profit->objpix);
  free(profit->objweight);
  free(profit->resi);
  free(profit->prof);
  free(profit->covar);
  QFFTWF_FREE(profit->psfdft);
  free(profit);

  return;
  }


/****** profit_fit ************************************************************
PROTO	void profit_fit(profitstruct *profit, picstruct *field,
		picstruct *wfield, objstruct *obj, obj2struct *obj2)
PURPOSE	Fit profile(s) convolved with the PSF to a detected object.
INPUT	Array of profile structures,
	Number of profiles,
	Pointer to the profile-fitting structure,
	Pointer to the field,
	Pointer to the field weight,
	Pointer to the obj.
OUTPUT	Pointer to an allocated fit structure (containing details about the
	fit).
NOTES	It is a modified version of the lm_minimize() of lmfit.
AUTHOR	E. Bertin (IAP)
VERSION	12/07/2012
 ***/
void	profit_fit(profitstruct *profit,
		picstruct *field, picstruct *wfield,
		objstruct *obj, obj2struct *obj2)
  {
    profitstruct	*pprofit, *qprofit;
    patternstruct	*pattern;
    psfstruct		*psf;
    checkstruct		*check;
    double		emx2,emy2,emxy, a , cp,sp, cn, bn, n,
			sump,sumq, sumpw2,sumqw2,sumpqw, sump0,sumq0, err;
    PIXTYPE		valp,valq,sig2;
    float		param0[PARAM_NPARAM], param1[PARAM_NPARAM],
			param[PARAM_NPARAM],
			**list,
			*cov,
			psf_fwhm, dchi2, aspect, chi2;
    int			*index,
			i,j,p, nparam, nparam2, ncomp, nprof;

  nparam = profit->nparam;
  nparam2 = nparam*nparam;
  nprof = profit->nprof;

  if (profit->psfdft)
    QFFTWF_FREE(profit->psfdft);


  psf = profit->psf;
  profit->pixstep = psf->pixstep;
  obj2->prof_flag = 0;

/* 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;
  profit->objnaxisn[0] = (((int)((obj->xmax-obj->xmin+1) + psf_fwhm + 0.499)
		*1.2)/2)*2 + 1;
  profit->objnaxisn[1] = (((int)((obj->ymax-obj->ymin+1) + psf_fwhm + 0.499)
		*1.2)/2)*2 + 1;
  if (profit->objnaxisn[1]<profit->objnaxisn[0])
    profit->objnaxisn[1] = profit->objnaxisn[0];
  else
    profit->objnaxisn[0] = profit->objnaxisn[1];
  if (profit->objnaxisn[0]>PROFIT_MAXOBJSIZE)
    {
    profit->subsamp = ceil((float)profit->objnaxisn[0]/PROFIT_MAXOBJSIZE);
    profit->objnaxisn[1] = (profit->objnaxisn[0] /= (int)profit->subsamp);
    obj2->prof_flag |= PROFLAG_OBJSUB;
    }
  else
    profit->subsamp = 1.0;
  profit->nobjpix = profit->objnaxisn[0]*profit->objnaxisn[1];

/* Use (dirty) global variables to interface with lmfit */
  the_field = field;
  the_wfield = wfield;
  theprofit = profit;
  profit->obj = obj;
  profit->obj2 = obj2;

  profit->nresi = profit_copyobjpix(profit, field, wfield);
/* Check if the number of constraints exceeds the number of free parameters */
  if (profit->nresi < nparam)
    {
    if (FLAG(obj2.prof_vector))
      for (p=0; p<nparam; p++)
        obj2->prof_vector[p] = 0.0;
    if (FLAG(obj2.prof_errvector))
      for (p=0; p<nparam; p++)
        obj2->prof_errvector[p] = 0.0;
    if (FLAG(obj2.prof_errmatrix))
      for (p=0; p<nparam2; p++)
        obj2->prof_errmatrix[p] = 0.0;
    obj2->prof_niter = 0;
    obj2->prof_flag |= PROFLAG_NOTCONST;
    return;
    }

/* Create pixmap at PSF resolution */
  profit->modnaxisn[0] =
	((int)(profit->objnaxisn[0]*profit->subsamp/profit->pixstep
		+0.4999)/2+1)*2; 
  profit->modnaxisn[1] =
	((int)(profit->objnaxisn[1]*profit->subsamp/profit->pixstep
		+0.4999)/2+1)*2; 
  if (profit->modnaxisn[1] < profit->modnaxisn[0])
    profit->modnaxisn[1] = profit->modnaxisn[0];
  else
    profit->modnaxisn[0] = profit->modnaxisn[1];
  if (profit->modnaxisn[0]>PROFIT_MAXMODSIZE)
    {
    profit->pixstep = (double)profit->modnaxisn[0] / PROFIT_MAXMODSIZE;
    profit->modnaxisn[0] = profit->modnaxisn[1] = PROFIT_MAXMODSIZE;
    obj2->prof_flag |= PROFLAG_MODSUB;
    }
  profit->nmodpix = profit->modnaxisn[0]*profit->modnaxisn[1];

/* Compute the local PSF */
  profit_psf(profit);

/* Set initial guesses and boundaries */
  profit->guesssigbkg = profit->sigma = obj->sigbkg;
  profit->guessdx = obj->mx - (int)(obj->mx+0.49999);
  profit->guessdy = obj->my - (int)(obj->my+0.49999);
  if ((profit->guessflux = obj2->flux_auto) <= 0.0)
    profit->guessflux = 0.0;
  if ((profit->guessfluxmax = 10.0*obj2->fluxerr_auto) <= profit->guessflux)
    profit->guessfluxmax = profit->guessflux;
  if (profit->guessfluxmax <= 0.0)
    profit->guessfluxmax = 1.0;
  if ((profit->guessradius = 0.5*psf->fwhm) < obj2->hl_radius)
    profit->guessradius = obj2->hl_radius;
  profit->guessaspect = obj->b/obj->a;
  profit->guessposang = obj->theta;

  profit_resetparams(profit);

/* Actual minimisation */
  fft_reset();
the_gal++;

/*
char str[1024];
sprintf(str, "obj_%04d.fits", the_gal);
catstruct *bcat;
float *bpix, *opix,*lmodpix,*objpix;
bcat=read_cat("base.fits");
QMALLOC(bpix, float, field->npix);
QFSEEK(bcat->file, bcat->tab->bodypos, SEEK_SET, bcat->filename);
read_body(bcat->tab, bpix, field->npix); 
free_cat(&bcat,1);
bcat=read_cat(str);
QMALLOC(opix, float, profit->nobjpix);
QFSEEK(bcat->file, bcat->tab->bodypos, SEEK_SET, bcat->filename);
read_body(bcat->tab, opix, profit->nobjpix); 
free_cat(&bcat,1);
addfrombig(bpix, field->width, field->height,
		profit->objpix, profit->objnaxisn[0],profit->objnaxisn[1],
		profit->ix,profit->iy, -1.0);
objpix = profit->objpix;
lmodpix = opix;
for (i=profit->nobjpix; i--;)
*(objpix++) += *(lmodpix++);
free(bpix);
free(opix);
*/
  profit->niter = profit_minimize(profit, PROFIT_MAXITER);
/*
profit_residuals(profit,field,wfield, 0.0, profit->paraminit, NULL);
check=initcheck(str, CHECK_OTHER,1);
check->width = profit->objnaxisn[0];
check->height = profit->objnaxisn[1];
reinitcheck(field,check);
memcpy(check->pix, profit->lmodpix, profit->nobjpix*sizeof(float));
reendcheck(field,check);
endcheck(check);

  chi2 = profit->chi2;
  for (p=0; p<nparam; p++)
    param1[p] = profit->paraminit[p];
  profit_resetparams(profit);
  for (p=0; p<nparam; p++)
    profit->paraminit[p] = param1[p] + (profit->paraminit[p]<param1[p]?1.0:-1.0)
			* sqrt(profit->covar[p*(nparam+1)]);

  profit->niter = profit_minimize(profit, PROFIT_MAXITER);
  if (chi2<profit->chi2)
    for (p=0; p<nparam; p++)
      profit->paraminit[p] = param1[p];

list = profit->paramlist;
index = profit->paramindex;
for (i=0; i<PARAM_NPARAM; i++)
if (list[i] && i!= PARAM_SPHEROID_ASPECT && i!=PARAM_SPHEROID_POSANG)
profit->freeparam_flag[index[i]] = 0;
profit->niter = profit_minimize(profit, PROFIT_MAXITER);
*/

  if (profit->nlimmin)
    obj2->prof_flag |= PROFLAG_MINLIM;
  if (profit->nlimmax)
    obj2->prof_flag |= PROFLAG_MAXLIM;

  for (p=0; p<nparam; p++)
    profit->paramerr[p]= sqrt(profit->covar[p*(nparam+1)]);

/* CHECK-Images */
  if ((check = prefs.check[CHECK_PROFILES]))
    {
    profit_residuals(profit,field,wfield, 0.0, profit->paraminit, NULL);
    addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1],
		profit->ix,profit->iy, 1.0);
    }

  if ((check = prefs.check[CHECK_SUBPROFILES]))
    {
    profit_residuals(profit,field,wfield, 0.0, profit->paraminit, NULL);
    addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1],
		profit->ix,profit->iy, -1.0);
    }
  if ((check = prefs.check[CHECK_SPHEROIDS]))
    {
/*-- Set to 0 flux components that do not belong to spheroids */
    for (p=0; p<profit->nparam; p++)
      param[p] = profit->paraminit[p];
    list = profit->paramlist;
    index = profit->paramindex;
    for (i=0; i<PARAM_NPARAM; i++)
      if (list[i] && flux_flag[i] && i!= PARAM_SPHEROID_FLUX)
        param[index[i]] = 0.0;
    profit_residuals(profit,field,wfield, 0.0, param, NULL);
    addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1],
		profit->ix,profit->iy, 1.0);
    }
  if ((check = prefs.check[CHECK_SUBSPHEROIDS]))
    {
/*-- Set to 0 flux components that do not belong to spheroids */
    for (p=0; p<profit->nparam; p++)
      param[p] = profit->paraminit[p];
    list = profit->paramlist;
    index = profit->paramindex;
    for (i=0; i<PARAM_NPARAM; i++)
      if (list[i] && flux_flag[i] && i!= PARAM_SPHEROID_FLUX)
        param[index[i]] = 0.0;
    profit_residuals(profit,field,wfield, 0.0, param, NULL);
    addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1],
		profit->ix,profit->iy, -1.0);
    }
  if ((check = prefs.check[CHECK_DISKS]))
    {
/*-- Set to 0 flux components that do not belong to disks */
    for (p=0; p<profit->nparam; p++)
      param[p] = profit->paraminit[p];
    list = profit->paramlist;
    index = profit->paramindex;
    for (i=0; i<PARAM_NPARAM; i++)
      if (list[i] && flux_flag[i] && i!= PARAM_DISK_FLUX)
        param[index[i]] = 0.0;
    profit_residuals(profit,field,wfield, 0.0, param, NULL);
    addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1],
		profit->ix,profit->iy, 1.0);
    }
  if ((check = prefs.check[CHECK_SUBDISKS]))
    {
/*-- Set to 0 flux components that do not belong to disks */
    for (p=0; p<profit->nparam; p++)
      param[p] = profit->paraminit[p];
    list = profit->paramlist;
    index = profit->paramindex;
    for (i=0; i<PARAM_NPARAM; i++)
      if (list[i] && flux_flag[i] && i!= PARAM_DISK_FLUX)
        param[index[i]] = 0.0;
    profit_residuals(profit,field,wfield, 0.0, param, NULL);
    addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1],
		profit->ix,profit->iy, -1.0);
    }
 
/* Compute compressed residuals */
  profit_residuals(profit,field,wfield, 10.0, profit->paraminit,profit->resi);

/* Fill measurement parameters */
  if (FLAG(obj2.prof_vector))
    {
    for (p=0; p<nparam; p++)
      obj2->prof_vector[p]= profit->paraminit[p];
    }
  if (FLAG(obj2.prof_errvector))
    {
    for (p=0; p<nparam; p++)
      obj2->prof_errvector[p]= profit->paramerr[p];
    }
  if (FLAG(obj2.prof_errmatrix))
    {
    for (p=0; p<nparam2; p++)
      obj2->prof_errmatrix[p]= profit->covar[p];
    }

  obj2->prof_niter = profit->niter;
  obj2->flux_prof = profit->flux;
  if (FLAG(obj2.fluxerr_prof))
    {
    err = 0.0;
    cov = profit->covar;
    index = profit->paramindex;
    list = profit->paramlist;
    for (i=0; i<PARAM_NPARAM; i++)
      if (flux_flag[i] && list[i])
        {
        cov = profit->covar + nparam*index[i];
        for (j=0; j<PARAM_NPARAM; j++)
          if (flux_flag[j] && list[j])
            err += cov[index[j]];
        }
    obj2->fluxerr_prof = err>0.0? sqrt(err): 0.0;
    }

  obj2->prof_chi2 = (profit->nresi > profit->nparam)?
		profit->chi2 / (profit->nresi - profit->nparam) : 0.0;

/* Position */
  if (FLAG(obj2.x_prof))
    {
    i = profit->paramindex[PARAM_X];
    j = profit->paramindex[PARAM_Y];
/*-- Model coordinates follow the FITS convention (first pixel at 1,1) */
    if (profit->paramlist[PARAM_X])
      {
      obj2->x_prof = (double)profit->ix + *profit->paramlist[PARAM_X] + 1.0;
      obj2->poserrmx2_prof = emx2 = profit->covar[i*(nparam+1)];
      }
    else
      emx2 = 0.0;
    if (profit->paramlist[PARAM_Y])
      {
      obj2->y_prof = (double)profit->iy + *profit->paramlist[PARAM_Y] + 1.0;
      obj2->poserrmy2_prof = emy2 = profit->covar[j*(nparam+1)];
      }
    else
      emy2 = 0.0;
    if (profit->paramlist[PARAM_X] && profit->paramlist[PARAM_Y])
      obj2->poserrmxy_prof = emxy = profit->covar[i+j*nparam];
    else
      emxy = 0.0;

/*-- Error ellipse parameters */
    if (FLAG(obj2.poserra_prof))
      {
       double	pmx2,pmy2,temp,theta;

      if (fabs(temp=emx2-emy2) > 0.0)
        theta = atan2(2.0 * emxy,temp) / 2.0;
      else
        theta = PI/4.0;

      temp = sqrt(0.25*temp*temp+ emxy*emxy);
      pmy2 = pmx2 = 0.5*(emx2+emy2);
      pmx2+=temp;
      pmy2-=temp;

      obj2->poserra_prof = (float)sqrt(pmx2>0.0? pmx2 : 0.0);
      obj2->poserrb_prof = (float)sqrt(pmy2>0.0? pmy2 : 0.0);
      obj2->poserrtheta_prof = (float)(theta/DEG);
      }

    if (FLAG(obj2.poserrcxx_prof))
      {
       double	temp;

      obj2->poserrcxx_prof = (float)(emy2/(temp=emx2*emy2-emxy*emxy));
      obj2->poserrcyy_prof = (float)(emx2/temp);
      obj2->poserrcxy_prof = (float)(-2*emxy/temp);
      }
    }

/* Equivalent noise area */
  if (FLAG(obj2.prof_noisearea))
    obj2->prof_noisearea = profit_noisearea(profit);

/* Second order moments and ellipticities */
  if (FLAG(obj2.prof_mx2))
    profit_moments(profit, obj2);

/* Second order moments of the convolved model (used by other parameters) */
  if (FLAG(obj2.prof_convmx2))
    profit_convmoments(profit, obj2);

/* "Hybrid" magnitudes */
  if (FLAG(obj2.fluxcor_prof))
    {
    profit_residuals(profit,field,wfield, 0.0, profit->paraminit, NULL);
    profit_fluxcor(profit, obj, obj2);
    }

/* Do measurements on the rasterised model (surface brightnesses) */
  if (FLAG(obj2.fluxeff_prof))
    profit_surface(profit, obj2); 

/* Background offset */
  if (FLAG(obj2.prof_offset_flux))
    {
    obj2->prof_offset_flux = *profit->paramlist[PARAM_BACK];
    obj2->prof_offset_fluxerr=profit->paramerr[profit->paramindex[PARAM_BACK]];
    }

/* Point source */
  if (FLAG(obj2.prof_dirac_flux))
    {
    obj2->prof_dirac_flux = *profit->paramlist[PARAM_DIRAC_FLUX];
    obj2->prof_dirac_fluxerr =
		profit->paramerr[profit->paramindex[PARAM_DIRAC_FLUX]];
    }

/* Spheroid */
  if (FLAG(obj2.prof_spheroid_flux))
    {
    if ((aspect = *profit->paramlist[PARAM_SPHEROID_ASPECT]) > 1.0)
      {
      *profit->paramlist[PARAM_SPHEROID_REFF] *= aspect;
      profit->paramerr[profit->paramindex[PARAM_SPHEROID_REFF]] *= aspect;
      profit->paramerr[profit->paramindex[PARAM_SPHEROID_ASPECT]]
			/= (aspect*aspect);
      *profit->paramlist[PARAM_SPHEROID_ASPECT] = 1.0 / aspect;
      *profit->paramlist[PARAM_SPHEROID_POSANG] += 90.0;
      }
    obj2->prof_spheroid_flux = *profit->paramlist[PARAM_SPHEROID_FLUX];
    obj2->prof_spheroid_fluxerr =
		profit->paramerr[profit->paramindex[PARAM_SPHEROID_FLUX]];
    obj2->prof_spheroid_reff = *profit->paramlist[PARAM_SPHEROID_REFF];
    obj2->prof_spheroid_refferr = 
		profit->paramerr[profit->paramindex[PARAM_SPHEROID_REFF]];
    obj2->prof_spheroid_aspect = *profit->paramlist[PARAM_SPHEROID_ASPECT];
    obj2->prof_spheroid_aspecterr = 
		profit->paramerr[profit->paramindex[PARAM_SPHEROID_ASPECT]];
    obj2->prof_spheroid_theta =
			fmod_m90_p90(*profit->paramlist[PARAM_SPHEROID_POSANG]);
    obj2->prof_spheroid_thetaerr = 
		profit->paramerr[profit->paramindex[PARAM_SPHEROID_POSANG]];
    if (FLAG(obj2.prof_spheroid_sersicn))
      {
      obj2->prof_spheroid_sersicn = *profit->paramlist[PARAM_SPHEROID_SERSICN];
      obj2->prof_spheroid_sersicnerr = 
		profit->paramerr[profit->paramindex[PARAM_SPHEROID_SERSICN]];
      }
    else
      obj2->prof_spheroid_sersicn = 4.0;
    if (FLAG(obj2.prof_spheroid_peak))
      {
      n = obj2->prof_spheroid_sersicn;
      bn = 2.0*n - 1.0/3.0 + 4.0/(405.0*n) + 46.0/(25515.0*n*n)
		+ 131.0/(1148175*n*n*n);	/* Ciotti & Bertin 1999 */
      cn = n * prof_gamma(2.0*n) * pow(bn, -2.0*n);
      obj2->prof_spheroid_peak = obj2->prof_spheroid_reff>0.0?
	obj2->prof_spheroid_flux
		/ (2.0 * PI * cn
		* obj2->prof_spheroid_reff*obj2->prof_spheroid_reff
		* obj2->prof_spheroid_aspect)
	: 0.0;
      if (FLAG(obj2.prof_spheroid_fluxeff))
        obj2->prof_spheroid_fluxeff = obj2->prof_spheroid_peak * exp(-bn);
      if (FLAG(obj2.prof_spheroid_fluxmean))
        obj2->prof_spheroid_fluxmean = obj2->prof_spheroid_peak * cn;
      }
    }

/* Disk */
  if (FLAG(obj2.prof_disk_flux))
    {
    if ((aspect = *profit->paramlist[PARAM_DISK_ASPECT]) > 1.0)
      {
      *profit->paramlist[PARAM_DISK_SCALE] *= aspect;
      profit->paramerr[profit->paramindex[PARAM_DISK_SCALE]] *= aspect;
      profit->paramerr[profit->paramindex[PARAM_DISK_ASPECT]]
			/= (aspect*aspect);
      *profit->paramlist[PARAM_DISK_ASPECT] = 1.0 / aspect;
      *profit->paramlist[PARAM_DISK_POSANG] += 90.0;
      }
    obj2->prof_disk_flux = *profit->paramlist[PARAM_DISK_FLUX];
    obj2->prof_disk_fluxerr =
		profit->paramerr[profit->paramindex[PARAM_DISK_FLUX]];
    obj2->prof_disk_scale = *profit->paramlist[PARAM_DISK_SCALE];
    obj2->prof_disk_scaleerr =
		profit->paramerr[profit->paramindex[PARAM_DISK_SCALE]];
    obj2->prof_disk_aspect = *profit->paramlist[PARAM_DISK_ASPECT];
    obj2->prof_disk_aspecterr =
		profit->paramerr[profit->paramindex[PARAM_DISK_ASPECT]];
    obj2->prof_disk_theta = fmod_m90_p90(*profit->paramlist[PARAM_DISK_POSANG]);
    obj2->prof_disk_thetaerr =
		profit->paramerr[profit->paramindex[PARAM_DISK_POSANG]];
    if (FLAG(obj2.prof_disk_inclination))
      {
      obj2->prof_disk_inclination = acos(obj2->prof_disk_aspect) / DEG;
      if (FLAG(obj2.prof_disk_inclinationerr))
        {
        a = sqrt(1.0-obj2->prof_disk_aspect*obj2->prof_disk_aspect);
        obj2->prof_disk_inclinationerr = obj2->prof_disk_aspecterr
					/(a>0.1? a : 0.1)/DEG;
        }
      }

    if (FLAG(obj2.prof_disk_peak))
      {
      obj2->prof_disk_peak = obj2->prof_disk_scale>0.0?
	obj2->prof_disk_flux
	/ (2.0 * PI * obj2->prof_disk_scale*obj2->prof_disk_scale
		* obj2->prof_disk_aspect)
	: 0.0;
      if (FLAG(obj2.prof_disk_fluxeff))
        obj2->prof_disk_fluxeff = obj2->prof_disk_peak * 0.186682; /* e^-(b_n)*/
      if (FLAG(obj2.prof_disk_fluxmean))
        obj2->prof_disk_fluxmean = obj2->prof_disk_peak * 0.355007;/* b_n^(-2)*/
      }

/* Disk pattern */
    if (prefs.pattern_flag)
      {
      profit_residuals(profit,field,wfield, PROFIT_DYNPARAM,
			profit->paraminit,profit->resi);
      pattern = pattern_init(profit, prefs.pattern_type,
		prefs.prof_disk_patternncomp);
      pattern_fit(pattern, profit);
      if (FLAG(obj2.prof_disk_patternspiral))
        obj2->prof_disk_patternspiral = pattern_spiral(pattern);
      if (FLAG(obj2.prof_disk_patternvector))
        {
        ncomp = pattern->size[2];
        for (p=0; p<ncomp; p++)
          obj2->prof_disk_patternvector[p] = (float)pattern->coeff[p];
        }
      if (FLAG(obj2.prof_disk_patternmodvector))
        {
        ncomp = pattern->ncomp*pattern->nfreq;
        for (p=0; p<ncomp; p++)
          obj2->prof_disk_patternmodvector[p] = (float)pattern->mcoeff[p];
        }
      if (FLAG(obj2.prof_disk_patternargvector))
        {
        ncomp = pattern->ncomp*pattern->nfreq;
        for (p=0; p<ncomp; p++)
          obj2->prof_disk_patternargvector[p] = (float)pattern->acoeff[p];
        }
      pattern_end(pattern);
      }

/* Bar */
    if (FLAG(obj2.prof_bar_flux))
      {
      obj2->prof_bar_flux = *profit->paramlist[PARAM_BAR_FLUX];
      obj2->prof_bar_fluxerr =
		profit->paramerr[profit->paramindex[PARAM_BAR_FLUX]];
      obj2->prof_bar_length = *profit->paramlist[PARAM_ARMS_START]
				**profit->paramlist[PARAM_DISK_SCALE];
      obj2->prof_bar_lengtherr = *profit->paramlist[PARAM_ARMS_START]
		  * profit->paramerr[profit->paramindex[PARAM_DISK_SCALE]]
		+ *profit->paramlist[PARAM_DISK_SCALE]
		  * profit->paramerr[profit->paramindex[PARAM_ARMS_START]];
      obj2->prof_bar_aspect = *profit->paramlist[PARAM_BAR_ASPECT];
      obj2->prof_bar_aspecterr =
		profit->paramerr[profit->paramindex[PARAM_BAR_ASPECT]];
      obj2->prof_bar_posang = 
			fmod_m90_p90(*profit->paramlist[PARAM_ARMS_POSANG]);
      obj2->prof_bar_posangerr =
		profit->paramerr[profit->paramindex[PARAM_ARMS_POSANG]];
      if (FLAG(obj2.prof_bar_theta))
        {
        cp = cos(obj2->prof_bar_posang*DEG);
        sp = sin(obj2->prof_bar_posang*DEG);
        a = obj2->prof_disk_aspect;
        obj2->prof_bar_theta = fmod_m90_p90(atan2(a*sp,cp)/DEG
				+ obj2->prof_disk_theta);
        obj2->prof_bar_thetaerr = obj2->prof_bar_posangerr*a/(cp*cp+a*a*sp*sp);
        }

/* Arms */
      if (FLAG(obj2.prof_arms_flux))
        {
        obj2->prof_arms_flux = *profit->paramlist[PARAM_ARMS_FLUX];
        obj2->prof_arms_fluxerr =
		profit->paramerr[profit->paramindex[PARAM_ARMS_FLUX]];
        obj2->prof_arms_pitch =
		fmod_m90_p90(*profit->paramlist[PARAM_ARMS_PITCH]);
        obj2->prof_arms_pitcherr =
		profit->paramerr[profit->paramindex[PARAM_ARMS_PITCH]];
        obj2->prof_arms_start = *profit->paramlist[PARAM_ARMS_START]
				**profit->paramlist[PARAM_DISK_SCALE];
        obj2->prof_arms_starterr = *profit->paramlist[PARAM_ARMS_START]
		  * profit->paramerr[profit->paramindex[PARAM_DISK_SCALE]]
		+ *profit->paramlist[PARAM_DISK_SCALE]
		  * profit->paramerr[profit->paramindex[PARAM_ARMS_START]];
        obj2->prof_arms_quadfrac = *profit->paramlist[PARAM_ARMS_QUADFRAC];
        obj2->prof_arms_quadfracerr =
		profit->paramerr[profit->paramindex[PARAM_ARMS_QUADFRAC]];
        obj2->prof_arms_posang =
			fmod_m90_p90(*profit->paramlist[PARAM_ARMS_POSANG]);
        obj2->prof_arms_posangerr =
		profit->paramerr[profit->paramindex[PARAM_ARMS_POSANG]];
        }
      }
    }

/* Star/galaxy classification */
  if (FLAG(obj2.prof_class_star) || FLAG(obj2.prof_concentration))
    {
    profit_residuals(profit,field,wfield, PROFIT_DYNPARAM, profit->paraminit,
	FLAG(obj2.prof_class_star) ? profit->resi : NULL);
    pprofit = thepprofit;
    nparam = pprofit->nparam;
    if (pprofit->psfdft)
      QFFTWF_FREE(pprofit->psfdft);
    psf = pprofit->psf;
    pprofit->pixstep = profit->pixstep;
    pprofit->guesssigbkg = profit->guesssigbkg;
    pprofit->guessdx = profit->guessdx;
    pprofit->guessdy = profit->guessdy;
    pprofit->guessflux = profit->guessflux;
    pprofit->guessfluxmax = profit->guessfluxmax;
    pprofit->guessradius = profit->guessradius;
    pprofit->guessaspect = profit->guessaspect;
    pprofit->guessposang = profit->guessposang;
    pprofit->ix = profit->ix;
    pprofit->iy = profit->iy;
    pprofit->objnaxisn[0] = profit->objnaxisn[0];
    pprofit->objnaxisn[1] = profit->objnaxisn[1];
    pprofit->subsamp = profit->subsamp;
    pprofit->nobjpix = profit->nobjpix;
    pprofit->obj = obj;
    pprofit->obj2 = obj2;
    pprofit->nresi = profit_copyobjpix(pprofit, field, wfield);
    pprofit->modnaxisn[0] = profit->modnaxisn[0];
    pprofit->modnaxisn[1] = profit->modnaxisn[1];
    pprofit->nmodpix = profit->nmodpix;
    profit_psf(pprofit);
    pprofit->sigma = obj->sigbkg;
    profit_resetparams(pprofit);
    if (profit->paramlist[PARAM_X] && profit->paramlist[PARAM_Y])
      {
      pprofit->paraminit[pprofit->paramindex[PARAM_X]] = *profit->paramlist[PARAM_X];
      pprofit->paraminit[pprofit->paramindex[PARAM_Y]] = *profit->paramlist[PARAM_Y];
      }
    fft_reset();
    pprofit->paraminit[pprofit->paramindex[PARAM_DIRAC_FLUX]] = profit->flux;
    pprofit->niter = profit_minimize(pprofit, PROFIT_MAXITER);
    profit_residuals(pprofit,field,wfield, PROFIT_DYNPARAM, pprofit->paraminit,
			FLAG(obj2.prof_class_star)? pprofit->resi : NULL);
    qprofit = theqprofit;
    nparam = qprofit->nparam;
    if (qprofit->psfdft)
      QFFTWF_FREE(qprofit->psfdft);
    qprofit->pixstep = profit->pixstep;
    qprofit->guesssigbkg = profit->guesssigbkg;
    qprofit->guessdx = profit->guessdx;
    qprofit->guessdy = profit->guessdy;
    qprofit->guessflux = profit->guessflux;
    qprofit->guessfluxmax = profit->guessfluxmax;
    qprofit->guessradius = profit->guessradius;
    qprofit->guessaspect = profit->guessaspect;
    qprofit->guessposang = profit->guessposang;
    qprofit->ix = profit->ix;
    qprofit->iy = profit->iy;
    qprofit->objnaxisn[0] = profit->objnaxisn[0];
    qprofit->objnaxisn[1] = profit->objnaxisn[1];
    qprofit->subsamp = profit->subsamp;
    qprofit->nobjpix = profit->nobjpix;
    qprofit->obj = obj;
    qprofit->obj2 = obj2;
    qprofit->nresi = profit_copyobjpix(qprofit, field, wfield);
    qprofit->modnaxisn[0] = profit->modnaxisn[0];
    qprofit->modnaxisn[1] = profit->modnaxisn[1];
    qprofit->nmodpix = profit->nmodpix;
    profit_psf(qprofit);
    qprofit->sigma = obj->sigbkg;
    profit_resetparams(qprofit);
    fft_reset();
    qprofit->paraminit[qprofit->paramindex[PARAM_X]]
		= pprofit->paraminit[pprofit->paramindex[PARAM_X]];
    qprofit->paraminit[qprofit->paramindex[PARAM_Y]]
		= pprofit->paraminit[pprofit->paramindex[PARAM_Y]];
    qprofit->paraminit[qprofit->paramindex[PARAM_DISK_FLUX]]
		= pprofit->paraminit[pprofit->paramindex[PARAM_DIRAC_FLUX]];
    qprofit->paraminit[qprofit->paramindex[PARAM_DISK_SCALE]] = psf->fwhm/16.0;
    qprofit->paraminit[qprofit->paramindex[PARAM_DISK_ASPECT]] = 1.0;
    qprofit->paraminit[qprofit->paramindex[PARAM_DISK_POSANG]] = 0.0;
    profit_residuals(qprofit,field,wfield, PROFIT_DYNPARAM, qprofit->paraminit,
			FLAG(obj2.prof_class_star)? qprofit->resi : NULL);
    sump = sumq = sumpw2 = sumqw2 = sumpqw = sump0 = sumq0 = 0.0;
    for (p=0; p<pprofit->nobjpix; p++)
      if (pprofit->objweight[p]>0 && pprofit->objpix[p]>-BIG)
        {
        valp = pprofit->lmodpix[p];
        sump += (double)(valp*pprofit->objpix[p]);
	valq = qprofit->lmodpix[p];
        sumq += (double)(valq*pprofit->objpix[p]);
	sump0 += (double)(valp*valp);
	sumq0 += (double)(valp*valq);
        sig2 = 1.0f/(pprofit->objweight[p]*pprofit->objweight[p]);
        sumpw2 += valp*valp*sig2;
        sumqw2 += valq*valq*sig2;
        sumpqw += valp*valq*sig2;
        }

    if (FLAG(obj2.prof_class_star))
      {
      dchi2 = 0.5*(pprofit->chi2 - profit->chi2);
      obj2->prof_class_star = dchi2 < 50.0?
	(dchi2 > -50.0? 2.0/(1.0+expf(dchi2)) : 2.0) : 0.0;
      }
    if (FLAG(obj2.prof_concentration))
      {
      obj2->prof_concentration = sump>0.0? (sumq/sump - sumq0/sump0) : 1.0;
      if (FLAG(obj2.prof_concentrationerr))
        obj2->prof_concentrationerr = (sump>0.0 && (err = sumqw2*sump*sump
		+sumpw2*sumq*sumq-2.0*sumpqw*sump*sumq)>0.0)?
			sqrt(err) / (sump*sump) : 1.0;
      }
    }

/* clean up. */
  fft_reset();

  return;
  }


/****** profit_dfit ***********************************************************
PROTO	void profit_dfit(profitstruct *profit, profitstruct *dprofit,
		picstruct *field, picstruct *dfield,
		picstruct *wfield, picstruct *dwfield,
		objstruct *obj, obj2struct *obj2)
PURPOSE	Fit profile(s) convolved with the PSF to a detected object on the
	detection image, and use the measurement image to scale the flux.
INPUT	Pointer to the measurement profile-fitting structure,
	pointer to the detection profile-fitting structure,
	pointer to the measurement field,
	pointer to the detection field,
	pointer to the measurement field weight,
	pointer to the detection field weight,
	pointer to the obj.
OUTPUT	Pointer to an allocated fit structure (containing details about the
	fit).
NOTES	It is a modified version of the lm_minimize() of lmfit.
AUTHOR	E. Bertin (IAP)
VERSION	12/07/2012
 ***/
void	profit_dfit(profitstruct *profit, profitstruct *dprofit,
		picstruct *field, picstruct *dfield,
		picstruct *wfield, picstruct *dwfield,
		objstruct *obj, obj2struct *obj2)
  {
    psfstruct		*dpsf;
    double		emx2,emy2,emxy, a , cp,sp, cn, bn, n,
			sumn,sumd;
    PIXTYPE		valn,vald,w2;
    float		param0[PARAM_NPARAM], param1[PARAM_NPARAM],
			param[PARAM_NPARAM],
			**list,
			*cov, *pix,
			psf_fwhm, dchi2, err, aspect, chi2, ffac;
    int			*index,
			i,j,p, nparam, nparam2, ncomp, nprof;

  nparam = dprofit->nparam;
  nparam2 = nparam*nparam;
  nprof = dprofit->nprof;
  if (dprofit->psfdft)
    QFFTWF_FREE(dprofit->psfdft);

  dpsf = dprofit->psf;
  dprofit->pixstep = dpsf->pixstep;
  obj2->dprof_flag = 0;

/* Create pixmaps at image resolution */
  dprofit->ix = (int)(obj->mx + 0.49999);/* internal convention: 1st pix = 0 */
  dprofit->iy = (int)(obj->my + 0.49999);/* internal convention: 1st pix = 0 */
  psf_fwhm = dpsf->masksize[0]*dpsf->pixstep;
  dprofit->objnaxisn[0] = (((int)((obj->xmax-obj->xmin+1) + psf_fwhm + 0.499)
		*1.2)/2)*2 + 1;
  dprofit->objnaxisn[1] = (((int)((obj->ymax-obj->ymin+1) + psf_fwhm + 0.499)
		*1.2)/2)*2 + 1;
  if (dprofit->objnaxisn[1]<dprofit->objnaxisn[0])
    dprofit->objnaxisn[1] = dprofit->objnaxisn[0];
  else
    dprofit->objnaxisn[0] = dprofit->objnaxisn[1];
  if (dprofit->objnaxisn[0]>PROFIT_MAXOBJSIZE)
    {
    dprofit->subsamp = ceil((float)dprofit->objnaxisn[0]/PROFIT_MAXOBJSIZE);
    dprofit->objnaxisn[1] = (dprofit->objnaxisn[0] /= (int)dprofit->subsamp);
    obj2->dprof_flag |= PROFLAG_OBJSUB;
    }
  else
    dprofit->subsamp = 1.0;
  dprofit->nobjpix = dprofit->objnaxisn[0]*dprofit->objnaxisn[1];

/* Use (dirty) global variables to interface with lmfit */
  the_field = dfield;
  the_wfield = dwfield;
  dprofit->obj = obj;
  dprofit->obj2 = obj2;

  dprofit->nresi = profit_copyobjpix(dprofit, dfield, dwfield);
/* Check if the number of constraints exceeds the number of free parameters */
  if (dprofit->nresi < nparam)
    {
    obj2->dprof_flag |= PROFLAG_NOTCONST;
    return;
    }

/* Create pixmap at PSF resolution */
  dprofit->modnaxisn[0] =
	((int)(dprofit->objnaxisn[0]*dprofit->subsamp/dprofit->pixstep
		+0.4999)/2+1)*2; 
  dprofit->modnaxisn[1] =
	((int)(dprofit->objnaxisn[1]*dprofit->subsamp/dprofit->pixstep
		+0.4999)/2+1)*2; 
  if (dprofit->modnaxisn[1] < dprofit->modnaxisn[0])
    dprofit->modnaxisn[1] = dprofit->modnaxisn[0];
  else
    dprofit->modnaxisn[0] = dprofit->modnaxisn[1];
  if (dprofit->modnaxisn[0]>PROFIT_MAXMODSIZE)
    {
    dprofit->pixstep = (double)dprofit->modnaxisn[0] / PROFIT_MAXMODSIZE;
    dprofit->modnaxisn[0] = dprofit->modnaxisn[1] = PROFIT_MAXMODSIZE;
    obj2->dprof_flag |= PROFLAG_MODSUB;
    }
  dprofit->nmodpix = dprofit->modnaxisn[0]*dprofit->modnaxisn[1];

/* Compute the local PSF */
  profit_psf(dprofit);

/* Set initial guesses and boundaries */
  dprofit->guesssigbkg = dprofit->sigma = obj->dsigbkg;
  dprofit->guessdx = obj->mx - (int)(obj->mx+0.49999);
  dprofit->guessdy = obj->my - (int)(obj->my+0.49999);
  if ((dprofit->guessflux = obj->dflux) <= 0.0)
    dprofit->guessflux = 0.0;
  if ((dprofit->guessfluxmax = 10.0*obj->dnpix*obj->dsigbkg*obj->dsigbkg)
			 <= dprofit->guessflux)
    dprofit->guessfluxmax = dprofit->guessflux;
  if (dprofit->guessfluxmax <= 0.0)
    dprofit->guessfluxmax = 1.0;
  if ((dprofit->guessradius = 0.5*dpsf->fwhm) < sqrtf((float)obj->dnpix))
    dprofit->guessradius = sqrtf((float)obj->dnpix);
  dprofit->guessaspect = obj->b/obj->a;
  dprofit->guessposang = obj->theta;

  profit_resetparams(dprofit);

/* Actual minimisation */
  fft_reset();

  dprofit->niter = profit_minimize(dprofit, PROFIT_MAXITER);