public software.sextractor

/*

*                               check.c

*

* Manage "check-images".

*

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

*

*       This file part of:      SExtractor

*

*       Copyright:              (C) 1993,1998-2010 IAP/CNRS/UPMC

*                               (C) 1994,1997 ESO

*                               (C) 1995,1996 Sterrewacht Leiden

*

*       Author:                 Emmanuel Bertin (IAP)

*

*       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:          11/10/2010

*

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/

#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        "fits/fitscat.h"

#include        "fitswcs.h"

#include        "check.h"

/********************************* addcheck **********************************/

/*

Add a PSF to a CHECK-image (with a multiplicative factor).

Outside boundaries are taken into account.

*/

void    addcheck(checkstruct *check, float *psf,

                        int w,int h, int ix,int iy, float amplitude)

  {

   PIXTYPE      *pix;

   int          x,y, xmin,xmax,ymin,ymax,w2, dwpsf;

/* Set the image boundaries */

  w2 = w;

  ymin = iy-h/2;

  ymax = ymin + h;

  if (ymin<0)

    {

    psf -= ymin*w;

    ymin = 0;

    }

  if (ymax>check->height)

    ymax = check->height;

  xmin = ix-w/2;

  xmax = xmin + w;

  if (xmax>check->width)

    {

    w2 -= xmax-check->width;

    xmax = check->width;

    }

  if (xmin<0)

    {

    psf -= xmin;

    w2 += xmin;

    xmin = 0;

    }

  dwpsf = w-w2;

/* Subtract the right pixels to the destination */

  for (y=ymin; y<ymax; y++, psf += dwpsf)

    {

    pix = (float *)check->pix+y*check->width+xmin;

    for (x=w2; x--;)

      *(pix++) += amplitude**(psf++);

    }

  return;

  }

/********************************* blankcheck *******************************/

/*

Blank a part of the CHECK-image according to a mask.

*/

void    blankcheck(checkstruct *check, PIXTYPE *mask, int w,int h,

                int xmin,int ymin, PIXTYPE val)

  {

   PIXTYPE      *pixt;

   int          x,y, xmax,ymax,w2,wc;

/* Don't go further if out of frame!! */

  if (xmin+w<0 || xmin>=check->width

        || ymin+h<0 || ymin>=check->height)

    return;

/* Set the image boundaries */

  w2 = w;

  ymax = ymin + h;

  if (ymin<0)

    {

    mask -= ymin*w;

    ymin = 0;

    }

  if (ymax>check->height)

    ymax = check->height;

  xmax = xmin + w;

  if (xmax>check->width)

    {

    w2 -= xmax - check->width;

    xmax = check->width;

    }

  if (xmin<0)

    {

    mask += -xmin;

    w2 -= -xmin;

    xmin = 0;

    }

  w -= w2;

  wc = check->width;

  ymin = ymin*wc+xmin;

  ymax = ymax*wc+xmin;

/* Blank the right pixels in the image */

  for (y=ymin; y<ymax; y+=wc, mask += w)

    {

    pixt = (float *)check->pix + y;

    for (x=w2; x--; pixt++)

      if (*(mask++) > -BIG)

        *pixt = val;

    }

  return;

  }

/******************************** initcheck **********************************/

/*

initialize check-image.

*/

checkstruct     *initcheck(char *filename, checkenum check_type, int next)

  {

   catstruct    *fitscat;

   checkstruct  *check;

  QCALLOC(check, checkstruct, 1);

  strcpy(check->filename, filename);

  check->type = check_type;

  if (next>1)

/*-- Create a "pure" primary HDU */

    {

    fitscat = new_cat(1);

    init_cat(fitscat);

    strcpy(fitscat->filename, filename);

    fitsadd(fitscat->tab->headbuf, "NEXTEND ", "Number of extensions");

    fitswrite(fitscat->tab->headbuf, "NEXTEND ", &next, H_INT, T_LONG);

    if (open_cat(fitscat, WRITE_ONLY) != RETURN_OK)

      error(EXIT_FAILURE,"*Error*: cannot open for writing ", filename);

    save_tab(fitscat, fitscat->tab);

    check->file = fitscat->file;

    fitscat->file = NULL;

    free_cat(&fitscat, 1);

    }

  else

    if (!(check->file = fopen(check->filename, "wb")))

      error(EXIT_FAILURE, "*Error*: Cannot open for output ", check->filename);

  return check;

  }

/******************************** reinitcheck ********************************/

/*

initialize check-image (for subsequent writing).

*/

void    reinitcheck(picstruct *field, checkstruct *check)

  {

   wcsstruct    *wcs;

   char         *buf;

   int          i, ival;

   size_t       padsize;

   double       dval;

   ULONG        *ptri;

   PIXTYPE      *ptrf;

/* Inherit the field FITS header */

  check->fitsheadsize = field->tab->headnblock*FBSIZE;

  QMEMCPY(field->tab->headbuf, check->fitshead, char, check->fitsheadsize);

  check->y = 0;

/* Neutralize possible scaling factors */

  dval = 1.0;fitswrite(check->fitshead, "BSCALE  ", &dval, H_FLOAT, T_DOUBLE);

  dval = 0.0;fitswrite(check->fitshead, "BZERO   ", &dval, H_FLOAT, T_DOUBLE);

  ival = 1;fitswrite(check->fitshead, "BITSGN  ", &ival, H_INT, T_LONG);

  if (field->tab->compress_type != COMPRESS_NONE)

    fitswrite(check->fitshead, "IMAGECOD", "NONE", H_STRING, T_STRING);

  fitswrite(check->fitshead, "ORIGIN  ", BANNER, H_STRING, T_STRING);

  switch(check->type)

    {

    case CHECK_IDENTICAL:

    case CHECK_BACKGROUND:

    case CHECK_FILTERED:

    case CHECK_SUBTRACTED:

      ival = -32;

      fitswrite(check->fitshead, "BITPIX  ", &ival, H_INT, T_LONG);

      check->width = field->width;

      check->height = field->height;

      check->npix = field->npix;

      QMALLOC(ptrf, PIXTYPE, check->width);

      check->pix = (void *)ptrf;

      QFWRITE(check->fitshead,check->fitsheadsize,check->file,check->filename);

      free(check->fitshead);

      break;

    case CHECK_BACKRMS:

    case CHECK_SUBOBJECTS:

      ival = -32;

      fitswrite(check->fitshead, "BITPIX  ", &ival, H_INT, T_LONG);

      check->width = field->width;

      check->height = field->height;

      check->npix = field->npix;

      QMALLOC(ptrf, PIXTYPE, check->width);

      check->pix = (void *)ptrf;

      QFWRITE(check->fitshead,check->fitsheadsize,check->file,check->filename);

      free(check->fitshead);

/*---- Allocate memory for replacing the blanked pixels by 0 */

      if (!check->line)

        QMALLOC(check->line, PIXTYPE, field->width);

      break;

    case CHECK_OBJECTS:

    case CHECK_APERTURES:

    case CHECK_PSFPROTOS:

    case CHECK_SUBPSFPROTOS:

    case CHECK_PCPROTOS:

    case CHECK_SUBPCPROTOS:

    case CHECK_PCOPROTOS:

    case CHECK_PROFILES:

    case CHECK_SUBPROFILES:

    case CHECK_SPHEROIDS:

    case CHECK_SUBSPHEROIDS:

    case CHECK_DISKS:

    case CHECK_SUBDISKS:

    case CHECK_PATTERNS:

      ival = -32;

      fitswrite(check->fitshead, "BITPIX  ", &ival, H_INT, T_LONG);

      check->width = field->width;

      check->height = field->height;

      check->npix = field->npix;

      check->overlay = 30*field->backsig;

      QCALLOC(ptrf, PIXTYPE, check->npix);

      check->pix = (void *)ptrf;

      QFWRITE(check->fitshead,check->fitsheadsize,check->file,check->filename);

      free(check->fitshead);

      break;

    case CHECK_SEGMENTATION:

      ival = 32;

      fitswrite(check->fitshead, "BITPIX  ", &ival, H_INT, T_LONG);

      check->width = field->width;

      check->height = field->height;

      check->npix = field->npix;

      QCALLOC(ptri, ULONG, check->npix);

      check->pix = (void *)ptri;

      QFWRITE(check->fitshead,check->fitsheadsize,check->file,check->filename);

      free(check->fitshead);

      break;

    case CHECK_ASSOC:

      ival = -32;

      fitswrite(check->fitshead, "BITPIX  ", &ival, H_INT, T_LONG);

      check->width = field->width;

      check->height = field->height;

      check->npix = field->npix;

      QMALLOC(ptrf, PIXTYPE, check->npix);

      check->pix = (void *)ptrf;

/*---- Initialize the pixmap to IEEE NaN */

      memset(ptrf, 0xFF, check->npix*sizeof(LONG));

      QFWRITE(check->fitshead,check->fitsheadsize,check->file,check->filename);

      free(check->fitshead);

      break;

    case CHECK_MINIBACKGROUND:

    case CHECK_MINIBACKRMS:

      ival = -32;

      fitswrite(check->fitshead, "BITPIX  ", &ival, H_INT, T_LONG);

      check->width = field->nbackx;

      fitswrite(check->fitshead, "NAXIS1  ", &check->width, H_INT, T_LONG);

      check->height = field->nbacky;

      fitswrite(check->fitshead, "NAXIS2  ", &check->height, H_INT, T_LONG);

/*---- Scale the WCS information if present */

      if ((wcs=field->wcs))

        {

        dval = wcs->cdelt[0]*field->backw;

        fitswrite(check->fitshead, "CDELT1  ", &dval, H_EXPO, T_DOUBLE);

        dval = wcs->cdelt[1]*field->backh;

        fitswrite(check->fitshead, "CDELT2  ", &dval, H_EXPO, T_DOUBLE);

        dval = (wcs->crpix[0]-0.5)/field->backw + 0.5;

        fitswrite(check->fitshead, "CRPIX1  ", &dval, H_EXPO, T_DOUBLE);

        dval = (wcs->crpix[1]-0.5)/field->backh + 0.5;

        fitswrite(check->fitshead, "CRPIX2  ", &dval, H_EXPO, T_DOUBLE);

        dval = wcs->cd[0]*field->backw;

        fitswrite(check->fitshead, "CD1_1   ", &dval, H_EXPO, T_DOUBLE);

        dval = wcs->cd[1]*field->backh;

        fitswrite(check->fitshead, "CD1_2  ", &dval, H_EXPO, T_DOUBLE);

        dval = wcs->cd[wcs->naxis]*field->backw;

        fitswrite(check->fitshead, "CD2_1   ", &dval, H_EXPO, T_DOUBLE);

        dval = wcs->cd[wcs->naxis+1]*field->backh;

        fitswrite(check->fitshead, "CD2_2  ", &dval, H_EXPO, T_DOUBLE);

        }

      check->npix = check->width*check->height;

      QMALLOC(ptrf, PIXTYPE, check->npix);

      check->pix = (void *)ptrf;

      if (check->type==CHECK_MINIBACKRMS)

        memcpy(check->pix, field->sigma, check->npix*sizeof(float));

      else

        memcpy(check->pix, field->back, check->npix*sizeof(float));

      QFWRITE(check->fitshead,check->fitsheadsize,check->file,check->filename);

      free(check->fitshead);

      if (bswapflag)

        swapbytes(check->pix, sizeof(float), (int)check->npix);

      QFWRITE(check->pix,check->npix*sizeof(float),check->file,

        check->filename);

/*---- Put the buffer back to its original state */

      if (bswapflag)

        swapbytes(check->pix, sizeof(float), (int)check->npix);

      free(check->pix);

      QCALLOC(buf, char, FBSIZE);

      padsize = (FBSIZE -((check->npix*sizeof(PIXTYPE))%FBSIZE))% FBSIZE;

      if (padsize)

        QFWRITE (buf, padsize, check->file, check->filename);

      free(buf);

      break;

    case CHECK_MAPSOM:

      ival = -32;

      fitswrite(check->fitshead, "BITPIX  ", &ival, H_INT, T_LONG);

      check->width = field->width;

      check->height = field->height;

      check->npix = field->npix;

      QMALLOC(ptrf, PIXTYPE, check->npix);

      check->pix = (void *)ptrf;

      for (i=check->npix; i--;)

        *(ptrf++) = -10.0;

      QFWRITE(check->fitshead,check->fitsheadsize,check->file,check->filename);

      free(check->fitshead);

      break;

    case CHECK_OTHER:

      ival = -32;

      fitswrite(check->fitshead, "BITPIX  ", &ival, H_INT, T_LONG);

      fitswrite(check->fitshead, "NAXIS1  ", &check->width, H_INT, T_LONG);

      fitswrite(check->fitshead, "NAXIS2  ", &check->height, H_INT, T_LONG);

      check->npix = check->width*check->height;

      QMALLOC(ptrf, PIXTYPE, check->npix);

      check->pix = (void *)ptrf;

      for (i=check->npix; i--;)

        *(ptrf++) = -10.0;

      QFWRITE(check->fitshead,check->fitsheadsize,check->file,check->filename);

      free(check->fitshead);

      break;

    default:

      error(EXIT_FAILURE, "*Internal Error* in ", "reinitcheck()!");

    }

  return;

  }

/******************************** writecheck *********************************/

/*

Write ONE line of npix pixels of a check-image.

*/

void    writecheck(checkstruct *check, PIXTYPE *data, int w)

  {

  if (check->type == CHECK_APERTURES || check->type == CHECK_SUBPSFPROTOS

        || check->type == CHECK_SUBPCPROTOS || check->type == CHECK_PCOPROTOS

        || check->type == CHECK_SUBPROFILES || check->type == CHECK_SUBSPHEROIDS

        || check->type == CHECK_SUBDISKS)

    {

    memcpy((PIXTYPE *)check->pix + w*(check->y++), data, w*sizeof(PIXTYPE));

    return;

    }

  else if (check->type == CHECK_SUBOBJECTS)

    {

     int        i;

     PIXTYPE    *pixt;

    pixt = check->line;

    for (i=w; i--; data++)

      *(pixt++) = (*data>-BIG)? *data:0.0;

    data = check->line;

    }

  if (bswapflag)

    swapbytes(data, sizeof(PIXTYPE), w);

  QFWRITE(data, w*sizeof(PIXTYPE), check->file, check->filename);

  if (bswapflag)

/*-- Put the buffer back to its original state */

    swapbytes(data, sizeof(PIXTYPE), w);

  return;

  }

/********************************* reendcheck ********************************/

/*

Finish current check-image.

*/

void    reendcheck(picstruct *field, checkstruct *check)

  {

   char         *buf;

   size_t       padsize;

  padsize = 0;                           /* To avoid gcc -Wall warnings */

  switch(check->type)

    {

    case CHECK_MINIBACKGROUND:

    case CHECK_MINIBACKRMS:

      return;

    case CHECK_IDENTICAL:

    case CHECK_BACKGROUND:

    case CHECK_BACKRMS:

    case CHECK_FILTERED:

    case CHECK_SUBTRACTED:

      free(check->pix);

      free(check->line);

      check->line = NULL;

      padsize = (FBSIZE -((check->npix*sizeof(PIXTYPE))%FBSIZE)) % FBSIZE;

      break;

    case CHECK_OBJECTS:

    case CHECK_APERTURES:

    case CHECK_PSFPROTOS:

    case CHECK_SUBPSFPROTOS:

    case CHECK_PCPROTOS:

    case CHECK_SUBPCPROTOS:

    case CHECK_PCOPROTOS:

    case CHECK_PROFILES:

    case CHECK_SUBPROFILES:

    case CHECK_SPHEROIDS:

    case CHECK_SUBSPHEROIDS:

    case CHECK_DISKS:

    case CHECK_SUBDISKS:

    case CHECK_ASSOC:

    case CHECK_PATTERNS:

    case CHECK_MAPSOM:

    case CHECK_OTHER:

      if (bswapflag)

        swapbytes(check->pix, sizeof(PIXTYPE), (int)check->npix);

      QFWRITE(check->pix,check->npix*sizeof(PIXTYPE),

                check->file,check->filename);

      free(check->pix);

      padsize = (FBSIZE-((check->npix*sizeof(PIXTYPE))%FBSIZE)) % FBSIZE;

      break;

    case CHECK_SEGMENTATION:

      if (bswapflag)

        swapbytes(check->pix, sizeof(ULONG), (int)check->npix);

      QFWRITE(check->pix,check->npix*sizeof(ULONG),

                check->file,check->filename);

      free(check->pix);

      padsize = (FBSIZE -((check->npix*sizeof(ULONG))%FBSIZE)) % FBSIZE;

      break;

    case CHECK_SUBOBJECTS:

      {

       int      y;

      for (y=field->ymin; y<field->ymax; y++)

        writecheck(check, &PIX(field, 0, y), field->width);

      free(check->pix);

      free(check->line);

      check->line = NULL;

      padsize = (FBSIZE -((check->npix*sizeof(PIXTYPE))%FBSIZE)) % FBSIZE;

      break;

      }

    default:

      error(EXIT_FAILURE, "*Internal Error* in ", "endcheck()!");

    }

  QCALLOC(buf, char, FBSIZE);

  if (padsize)

    QFWRITE (buf, padsize, check->file, check->filename);

  free(buf);

  return;

  }

/********************************* endcheck **********************************/

/*

close check-image.

*/

void    endcheck(checkstruct *check)

  {

  fclose(check->file);

  free(check);

  return;

  }