back.c

/*
*				back.c
*
* Functions dealing with the image background.
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
*	This file part of:	SExtractor
*
*	Copyright:		(C) 1993-2014 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:		17/06/2014
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/

#ifdef HAVE_CONFIG_H
#include        "config.h"
#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	"back.h"
#include	"field.h"
#include	"weight.h"

/******************************** makeback ***********************************/
/*
Background maps are established from the images themselves; thus we need to
make at least one first pass through the data.
*/
void	makeback(picstruct *field, picstruct *wfield, int wscale_flag)

  {
   backstruct	*backmesh,*wbackmesh, *bm,*wbm;
   tabstruct	*tab, *wtab;
   PIXTYPE	*buf,*wbuf, *buft,*wbuft;
   OFF_T2	fcurpos,wfcurpos, wfcurpos2,fcurpos2, bufshift, jumpsize;
   OFF_T2	currentElement, wcurrentElement, currentElement2, wcurrentElement2;
   size_t	bufsize, bufsize2,
		size,meshsize;
   int		i,j,k,m,n, step, nlines,
		w,bw, bh, nx,ny,nb,
		lflag, nr;
   float	*ratio,*ratiop, *weight, *sigma,
		sratio, sigfac;

/* If the weight-map is not an external one, no stats are needed for it */
  if (wfield && wfield->flags&(INTERP_FIELD|BACKRMS_FIELD))
    wfield= NULL;
  tab = field->tab;
  if (wfield)
    wtab = wfield->tab;
  else
    wtab = NULL;	/* to avoid gcc -Wall warnings */
  w = field->width;
  bw = field->backw;
  bh = field->backh;
  nx = field->nbackx;
  ny = field->nbacky;
  nb = field->nback;

  NFPRINTF(OUTPUT, "Setting up background maps");

/* Decide if it is worth displaying progress each 16 lines */

  lflag = (field->width*field->backh >= (size_t)65536);

/* Save current positions in files */

  wfcurpos = wfcurpos2 = 0;		/* to avoid gcc -Wall warnings */
  QFTELL(field->file, fcurpos, field->filename);
  currentElement = (tab->currentElement == 0) ? 1 : tab->currentElement; // CFITSIO

  if (wfield)
    {
    QFTELL(wfield->file, wfcurpos, wfield->filename);
    wcurrentElement = (wtab->currentElement == 0) ? 1 : wtab->currentElement; // CFITSIO
    }

/* Allocate a correct amount of memory to store pixels */

  bufsize = (OFF_T2)w*bh;
  meshsize = (size_t)bufsize;
  nlines = 0;
  if (bufsize > (size_t)BACK_BUFSIZE)
    {
    nlines = BACK_BUFSIZE/w;
    step = (field->backh-1)/nlines+1;
    bufsize = (size_t)(nlines = field->backh/step)*w;
    bufshift = (step/2)*(OFF_T2)w;
    jumpsize = (step-1)*(OFF_T2)w;
    }
  else
    bufshift = jumpsize = 0;		/* to avoid gcc -Wall warnings */

/* Allocate some memory */
  QMALLOC(backmesh, backstruct, nx);		/* background information */
  QMALLOC(buf, PIXTYPE, bufsize);		/* pixel buffer */
  free(field->back);
  QMALLOC(field->back, float, nb);		/* background map */
  free(field->backline);
  QMALLOC(field->backline, PIXTYPE, w);		/* current background line */
  free(field->sigma);
  QMALLOC(field->sigma, float, nb);		/* sigma map */
  if (wfield)
    {
    QMALLOC(wbackmesh, backstruct, nx);		/* background information */
    QMALLOC(wbuf, PIXTYPE, bufsize);		/* pixel buffer */
    free(wfield->back);
    QMALLOC(wfield->back, float, nb);		/* background map */
    free(wfield->backline);
    QMALLOC(wfield->backline, PIXTYPE, w);	/* current background line */
    free(wfield->sigma);
    QMALLOC(wfield->sigma, float, nb);		/* sigma map */
    wfield->sigfac = 1.0;
    }
  else
    {
    wbackmesh = NULL;
    wbuf = NULL;
    }

/* Loop over the data packets */

  for (j=0; j<ny; j++)
    {
    if (lflag && j)
      NPRINTF(OUTPUT, "\33[1M> Setting up background map at line:%5d\n\33[1A",
	      j*bh);
    if (!nlines)
      {
/*---- The image is small enough so that we can make exhaustive stats */
      if (j == ny-1 && field->npix%bufsize)
        bufsize = field->npix%bufsize;
      read_body(field->tab, buf, bufsize);
      if (wfield)
        {
        read_body(wfield->tab, wbuf, bufsize);
        weight_to_var(wfield, wbuf, bufsize);
        }
/*---- Build the histograms */
      backstat(backmesh, wbackmesh, buf, wbuf, bufsize,nx, w, bw,
	wfield?wfield->weight_thresh:0.0);
      bm = backmesh;
      for (m=nx; m--; bm++)
        if (bm->mean <= -BIG)
          bm->histo=NULL;
        else
          QCALLOC(bm->histo, LONG, bm->nlevels);
      if (wfield)
        {
        wbm = wbackmesh;
        for (m=nx; m--; wbm++)
          if (wbm->mean <= -BIG)
            wbm->histo=NULL;
          else
            QCALLOC(wbm->histo, LONG, wbm->nlevels);
        }
      backhisto(backmesh, wbackmesh, buf, wbuf, bufsize,nx, w, bw,
	wfield?wfield->weight_thresh:0.0);
      }
    else
      {
/*---- Image size too big, we have to skip a few data !*/
      QFTELL(field->file, fcurpos2, field->filename);
      currentElement2 = (tab->currentElement == 0) ? 1 : tab->currentElement; // CFITSIO
      if (wfield){
        QFTELL(wfield->file, wfcurpos2, wfield->filename);
        wcurrentElement2 = (wtab->currentElement == 0) ? 1 : wtab->currentElement; // CFITSIO
      }
      if (j == ny-1 && (n=field->height%field->backh))
        {
        meshsize = n*(size_t)w;
        nlines = BACK_BUFSIZE/w;
        step = (n-1)/nlines+1;
        bufsize = (nlines = n/step)*(size_t)w;
        bufshift = (step/2)*(OFF_T2)w;
        jumpsize = (step-1)*(OFF_T2)w;
        free(buf);
        QMALLOC(buf, PIXTYPE, bufsize);		/* pixel buffer */