/*
                                  coadd.c

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
*       Part of:        Swarp
*
*       Author:         E.BERTIN (IAP)
*
*       Contents:       Coaddition routines
*
*       Last modify:    15/03/2001
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/

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

#include "define.h"
#include "types.h"
#include "globals.h"
#include "fitscat.h"
#include "fitswcs.h"
#include "coadd.h"
#include "data.h"
#include "field.h"
#include "interpolate.h"
#include "prefs.h"
#include "weight.h"
#include "wcs/wcs.h"

char	padbuf[FBSIZE];

/****** resample_field *******************************************************
PROTO	void resample_field(fieldstruct **pinfield, fieldstruct **pinwfield,
			fieldstruct *outfield, fieldstruct *outwfield,
			interpenum *interptype)
PURPOSE	Resample an image.
INPUT	Input pointer to field structure pointer,
	Input pointer weight-map field structure pointer,
	Output total field structure pointer,
	Output total weight-map field structure pointer,
	Interpolation type.
	Output
OUTPUT	-.
NOTES	The structure pointers pointed by pinfield and and pinwfield are
	updated and point to the resampled fields on output.
AUTHOR	E. Bertin (IAP)
VERSION	07/11/2000
 ***/
void	resample_field(fieldstruct **pinfield, fieldstruct **pinwfield,
		fieldstruct *outfield, fieldstruct *outwfield,
		interpenum *interptype)
  {
   fieldstruct		*infield,*inwfield, *field, *wfield;
   wcsstruct		wcsloc,
			*wcs;
   char			filename[MAXCHAR],filename2[MAXCHAR],
			*pstr;
   PIXTYPE		*outbuf,*outwbuf, *out,*outw,
			pix,pixw;
   static double	rawmin[NAXIS], rawmax[NAXIS], rawpos[NAXIS],
			rawposover[NAXIS], stepover[NAXIS];
   double		*wcsbuf,*wcspos, *rawbuf,*rawbufc, 
			dynfac, val, dpix,dpixw;
   static int		nstepover[NAXIS];
   int			*oversamp,
			d, o, x,y, dynflag, naxis, width, height, size,
			noversamp, oversampflag, ninput;

  infield = *pinfield;
  inwfield = *pinwfield;

/* Create new file name */
  strcpy(filename2, infield->rfilename);
  *gstr = '\0';
  if (pstr=strrchr(filename2, '.'))
    {
    strcpy(gstr, pstr);
    *pstr = '\0';
    sprintf(filename, "%s/coadd.%s.%04d%s", prefs.swapdir_name, filename2,
				infield->frameno, gstr);
    }
  else
    sprintf(filename, "%s/coadd.%s.%04d", prefs.swapdir_name, filename2,
				infield->frameno);

/* We make a copy of the output field */
  field = inherit_field(filename, outfield, FIELD_WRITE);
  field->backmean = infield->backmean;
  field->backsig = infield->backsig;

/* Now modify some characteristics of the output file */
  if (open_cat(field->cat, WRITE_ONLY) != RETURN_OK)
      error(EXIT_FAILURE, "*Error*: cannot open for writing ", filename);
  if (prefs.removetmp_flag)
    add_cleanupfilename(filename);

/* We use a small dirty trick to frame the output */
  wcsloc = *(outfield->wcs);
  wcs = infield->wcs;
  naxis = wcsloc.naxis;
  for (d=0; d<naxis;d++)
    {
    wcsloc.outmin[d] = wcs->outmin[d] < 1.0? 1.0:wcs->outmin[d];
    wcsloc.outmax[d] = wcs->outmax[d] > wcsloc.naxisn[d]?
			(double)wcsloc.naxisn[d] : wcs->outmax[d];
    wcsloc.crpix[d] -= (wcsloc.outmin[d]-1.0);
    wcsloc.naxisn[d] = (wcsloc.outmax[d]-wcsloc.outmin[d]+1);
    }

  wcs = field->wcs = copy_wcs(&wcsloc);
  write_wcs(field->tab, wcs);

/* Update field characteristics */
  width = field->width = field->tab->naxisn[0];
  field->height = 1;
  for (d=1; d<naxis; d++)
    field->height *= field->tab->naxisn[d];
  height = field->height;
  field->npix = field->width*field->height;

/* Write image header */
  QFWRITE(field->tab->headbuf, field->tab->headnblock*FBSIZE,
	field->cat->file, filename);

/* Now go on with output weight-map */
  if (*gstr)
    sprintf(filename, "%s/coadd.%s.%04d.weight%s", prefs.swapdir_name,
					filename2,
					infield->frameno, gstr);
  else
    sprintf(filename, "%s/coadd.%s.%04d.weight", prefs.swapdir_name, filename2,
					infield->frameno);

  wfield = init_weight(filename, field);
  if (prefs.removetmp_flag)
    add_cleanupfilename(filename);

  if (inwfield)
    set_weightconv(inwfield);
/* Write image header */
  if (open_cat(wfield->cat, WRITE_ONLY) != RETURN_OK)
    error(EXIT_FAILURE, "*Error*: cannot open for writing ", filename);

  QFWRITE(wfield->tab->headbuf, wfield->tab->headnblock*FBSIZE,
	wfield->cat->file, filename);

/* Prepare oversampling stuff */
  noversamp = 1;
  oversamp = prefs.oversamp;
  for (d=0; d<naxis; d++)
    {
    if (!oversamp[d])
/*-- Automatic oversampling mode */
      oversamp[d] = wcs->wcsscale[d]>infield->wcs->wcsscale[d]?
	(int)(wcs->wcsscale[d]/infield->wcs->wcsscale[d]+0.5)
	: 1;
    noversamp *= oversamp[d];
    }
  oversampflag = (noversamp>1);

/* Allocate memory for the output buffers that contain the final data in */
/* internal format (PIXTYPE) */
  QMALLOC(outbuf, PIXTYPE, width);
  QMALLOC(outwbuf, PIXTYPE, width);
/* Provide memory space for the current astrometric line */
  QMALLOC(wcsbuf, double, naxis*width*noversamp);
  QMALLOC(rawbuf, double, naxis*width*noversamp);

/* Initialize the astrometric vector */
  for (d=0; d<naxis;d++)
    {
    rawmin[d] = rawpos[d] = 0.5 + 0.5*(stepover[d]=1.0/oversamp[d]);
    rawmax[d] = (double)wcs->naxisn[d];
    }

/* Initialize interpolation kernel */
  init_kernel(interptype, naxis);

/* Set some commonly used variables */
  dynflag = (infield->cflags & CONVERT_DYNCOMP);
  dynfac = infield->dyncomp_fac;

/* Go on with resampling!! */
/* Loop over output lines: this can be over more than 1 (Y) dimension */
  for (y=0; y<height; y++)
    {
    if (!y || !((y+1)%4))
      NPRINTF(OUTPUT, "\33[1M> Resampling line:%7d / %-7d\n\33[1A",
	y+1, height);
/*-- Compute all alpha's and delta's for the current line */
    wcspos = wcsbuf;
    rawpos[0] = rawmin[0];
    if (oversampflag)
      {
      for (x=width; x--; (*rawpos)+=1.0)
        {
        for (d=0; d<naxis; d++)
          rawposover[d] = rawpos[d];
        for (o=noversamp; o--; wcspos+=naxis)
          {
          raw_to_wcs(field->wcs, rawposover, wcspos);
/*----- Update coordinate vector */
          for (d=0; d<naxis; d++)
            {
            rawposover[d] += stepover[d];
            if (++nstepover[d]<oversamp[d])
              break;
            else
              {
              rawposover[d] = rawpos[d];
              nstepover[d] = 0;
              }
            }
          }
        }
      }
    else
      for (x=width; x--; wcspos+=naxis, (*rawpos)+=1.0)
        raw_to_wcs(field->wcs, rawpos, wcspos);
/*-- Interpolate input images over the next output line */
    wcspos = wcsbuf;
    rawbufc = rawbuf;
    for (x=width*noversamp; x--; wcspos+=naxis, rawbufc+=naxis)
      if (*wcspos == WCS_NOCOORD)
        *rawbufc = WCS_NOCOORD;
      else
        wcs_to_raw(infield->wcs, wcspos, rawbufc);
    rawbufc = rawbuf;
    out = outbuf;
    outw = outwbuf;
    if (oversampflag)
      for (x=width; x--;)
        {
        ninput = 0;
        dpix = dpixw = 0.0;
        for (o=noversamp; o--; rawbufc+=naxis)
          if (*rawbufc != WCS_NOCOORD)
            {
            interpolate_pix(infield, inwfield, rawbufc, &pix, &pixw);
            if (pixw<BIG)
            {
            dpix += (double)pix;
            dpixw += (double)pixw;
            ninput++;
            }
          }
        *(out++) = ninput? dpix/ninput: 0.0;
/*------ Expand dynamic if it has been compressed */
        if (dynflag)
          {
          val = *(--out)/(2.0*dynfac) + 1.0;
          *(out++) = dynfac * (val*val - 1.0);
          }
        *(outw++) = ninput? (ninput*ninput)/dpixw:0.0;
        }
    else
      for (x=width; x--; rawbufc+=naxis)
        if (*rawbufc == WCS_NOCOORD)
          *(out++) = *(outw++) = 0.0;
        else
          {
          interpolate_pix(infield, inwfield, rawbufc, out++, outw);
/*-------- Expand dynamic if it has been compressed */
          if (dynflag)
            {
            val = *(--out)/(2.0*dynfac) + 1.0;
            *(out++) = dynfac * (val*val - 1.0);
            }
/*------- Convert variance to weight */
          *outw = (*outw < BIG) ? 1.0/(*outw) : 0.0;
          outw++;
          }

    write_body(field->tab, outbuf, width);
    write_body(wfield->tab, outwbuf, width);

/*-- Update coordinate vector */
    for (d=1; d<naxis; d++)
      if ((rawpos[d]+=1.0)<= rawmax[d])
        break;
      else
        rawpos[d] = rawmin[d];
    }

/* FITS padding*/
  size = PADEXTRA(field->tab->tabsize);
  if (size)
    QFWRITE(padbuf, (size_t)size, field->cat->file, field->filename);
  size = PADEXTRA(wfield->tab->tabsize);
  if (size)
    QFWRITE(padbuf, (size_t)size, wfield->cat->file, wfield->filename);

/* Close files */
  close_cat(field->cat);
  close_cat(wfield->cat);

/* Free interpolation kernel */
  free_kernel();

/* Return back new field pointers */
  *pinfield = field;
  *pinwfield = wfield;

/* Free memory */
  free(outbuf);
  free(outwbuf);
  free(wcsbuf);
  free(rawbuf);
  end_field(infield);
  if (inwfield)
    end_field(inwfield);

  return;
  }


/******* coadd_fields *********************************************************
PROTO	int coadd_fields(fieldstruct **infield, fieldstruct **inwfield,
			int ninput,
			fieldstruct *outfield, fieldstruct *outwfield,
			coaddenum coaddtype, PIXTYPE wthresh)
PURPOSE	Coadd images.
INPUT	Input field ptr array,
	Input weight field ptr array,
	number of input fields,
	Output field ptr,
	Output weight field ptr,
	Coaddition type.
OUTPUT	RETURN_OK if no error, or RETURN_ERROR in case of non-fatal error(s).
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 13/03/2001
 ***/
int coadd_fields(fieldstruct **infield, fieldstruct **inwfield,	int ninput,
			fieldstruct *outfield, fieldstruct *outwfield,
			coaddenum coaddtype, PIXTYPE wthresh)

  {
   fieldstruct		*field;
   tabstruct		*tab;
   wcsstruct		*wcs;
   PIXTYPE		*emptybuf,*emptywbuf, *inbuf,*inwbuf, *outbuf,*outwbuf,
			*pix,*wpix,
			*multibuf,*multiwbuf, *multi,*multiw;
   double		val;
   long			offbeg, offend;
   static int		rawmin[NAXIS], rawmax[NAXIS], rawpos[NAXIS],
			rawmaxmax[NAXIS];
   int			d,n, outwidth,multiwidth, width, height,
			naxis,nlines, x,y,yc,yp, size, offset;

  naxis = outfield->tab->naxis;
/* The output width is the useful length of the NAXIS1 axis */
  width = outfield->width;
/* The output ``height'' is the product of all other axis lengths */
  height = outfield->height;

/*-- Find global extrema of mapped image */
  for (d=0; d<naxis; d++)
    {
    rawmax[d] = -(rawmin[d] = 1<<30);
    for (n = 0; n<ninput; n++)
      {
      if (infield[n]->wcs->outmin[d] < rawmin[d])
        rawmin[d] = (int)infield[n]->wcs->outmin[d];
      if (infield[n]->wcs->outmax[d] > rawmax[d])
        rawmax[d] = (int)infield[n]->wcs->outmax[d];
      }
    rawpos[d] = 1;
    rawmaxmax[d] = outfield->wcs->naxisn[d];
    }

  if ((offbeg = (rawmin[0]-1)) < 0)
    offbeg = 0;
  if ((offend = (outfield->wcs->naxisn[0] - rawmax[0])) < 0)
    offend = 0;

  outwidth = outfield->wcs->naxisn[0] - offend - offbeg;

  nlines = 1;
  if (nlines<1)
    nlines = 1;
  multiwidth = outwidth*ninput;

/* Open output file and save header */
  if (open_cat(outfield->cat, WRITE_ONLY) != RETURN_OK)
      error(EXIT_FAILURE, "*Error*: cannot open for writing ",
		outfield->filename);
  QFWRITE(outfield->tab->headbuf, outfield->tab->headnblock*FBSIZE,
	outfield->cat->file, outfield->filename);

/* Open output weight file and save header */
  outwfield->sigfac = (float)ninput;
  if (open_cat(outwfield->cat, WRITE_ONLY) != RETURN_OK)
    error(EXIT_FAILURE, "*Error*: cannot open for writing ",
		outwfield->filename);
  QFWRITE(outwfield->tab->headbuf, outwfield->tab->headnblock*FBSIZE,
	outwfield->cat->file, outwfield->filename);
/* Only for WRITING the weights */
  set_weightconv(outwfield);

/* Re-open Input files */
  for (n = 0; n<ninput; n++)
    {
    if (open_cat(infield[n]->cat, READ_ONLY) != RETURN_OK)
      error(EXIT_FAILURE, "*Error*: cannot open for reading ",
		infield[n]->filename);
    tab = infield[n]->cat->tab;
    tab->bodypos = tab->headnblock*FBSIZE;
    QFSEEK(tab->cat->file, tab->bodypos, SEEK_SET, infield[n]->filename);
    if (open_cat(inwfield[n]->cat, READ_ONLY) != RETURN_OK)
      error(EXIT_FAILURE, "*Error*: cannot open for reading ",
		inwfield[n]->filename);
    tab = inwfield[n]->cat->tab;
    tab->bodypos = tab->headnblock*FBSIZE;
    QFSEEK(tab->cat->file, tab->bodypos, SEEK_SET, inwfield[n]->filename);
    }

/* Allocate memory for the output buffers that contain the input data in */
/* internal format (PIXTYPE) */
  QMALLOC(inbuf, PIXTYPE, nlines*outwidth);
  QMALLOC(inwbuf, PIXTYPE, nlines*outwidth);
/* Allocate memory for the output buffers that contain the final data in */
/* internal format (PIXTYPE) */
  QMALLOC(outbuf, PIXTYPE, nlines*outwidth);
  QMALLOC(outwbuf, PIXTYPE, nlines*outwidth);
/* Allocate memory for the output buffers that contain "empty data" */
  QCALLOC(emptybuf, PIXTYPE, width);
  QCALLOC(emptywbuf, PIXTYPE, width);
/* Allocate memory for the "multi-buffer" that stores resampled pixels from */
/* all images for the current line(s), prior to co-addition */
  QMALLOC(multibuf, PIXTYPE, nlines*multiwidth);
  QMALLOC(multiwbuf, PIXTYPE, nlines*multiwidth);
  multiw = NULL;

/* Loop over output ``lines'': this can be over more than 1 (Y) dimension */
  for (y = 0; y<height; y++)
    {
    if (!((y+1)%16))
      NPRINTF(OUTPUT, "\33[1M> Co-adding line:%7d / %-7d\n\33[1A",
	y+1, height);
/*---- Skip empty lines */
    for (d=naxis; --d;)
      if (rawpos[d]<rawmin[d] || rawpos[d]>rawmax[d])
        break;
    if (d>0)
      {
      write_body(outfield->tab, emptybuf, (long)width);
      write_body(outwfield->tab, emptywbuf, (long)width);
      }
    else
      {
/*---- Skip empty pixels */
      if (offbeg)
        {
        write_body(outfield->tab, emptybuf, offbeg);
        write_body(outwfield->tab, emptywbuf, offbeg);
        }

/*---- The current "physical" line is y modulo nlines */
      yc = y%nlines;

/*---- Initialize output data line */
      memset(multibuf+yc*multiwidth, 0, multiwidth*sizeof(PIXTYPE));
      multiw = multiwbuf+yc*multiwidth;
      for (x=multiwidth; x--;)
        *(multiw++) = BIG;

/*---- Loop over input images */
      for (n = 0; n<ninput; n++)
        {
        field = infield[n];
        wcs = field->wcs;
/*------ Check that the present line is contained in the field */
        for (d=naxis; --d;)
          if (rawpos[d]<wcs->outmin[d] || rawpos[d]>wcs->outmax[d])
            break;
        if (d>0)
          continue;
        read_body(field->tab, pix=inbuf, field->width);
        read_body(inwfield[n]->tab, wpix=inwbuf, field->width);
        offset = yc*multiwidth+(wcs->outmin[0]-1-offbeg)*ninput+n;
        multi = multibuf+offset;
        multiw = multiwbuf+offset;
        for (x=field->width; x--; wpix++)
          {
          *multi = *(pix++);
/*------- Convert weight to variance */
          *multiw = (*wpix> 0.0)? 1.0/(*wpix) : BIG;
          multi += ninput;
          multiw += ninput;
          }
        }

/*---- Now perform the coaddition itself */
      pix = outbuf;
      wpix = outwbuf;
      multi = multibuf+yc*multiwidth;
      multiw = multiwbuf+yc*multiwidth;
      for (x=outwidth; x--;)
        {
        coadd_pix(multi,multiw, ninput, pix++,wpix++,coaddtype,wthresh);
        multi += ninput;
        multiw += ninput;
        }

      write_body(outfield->tab, outbuf, outwidth);
      var_to_weight(outwbuf, outwidth);
      write_body(outwfield->tab, outwbuf, outwidth);

/*---- Skip empty pixels */
      if (offend)
        {
        write_body(outfield->tab, emptybuf, offend);
        write_body(outwfield->tab, emptywbuf, offend);
        }
      }

/*-- Update coordinate vector */
    for (d=1; d<naxis; d++)
      if ((++rawpos[d])<=rawmaxmax[d])
        break;
      else
        rawpos[d] = 1;

    }

/* FITS padding*/
  size = PADEXTRA(outfield->tab->tabsize);
  if (size)
    QFWRITE(padbuf, (size_t)size, outfield->cat->file, outfield->filename);
  size = PADEXTRA(outwfield->tab->tabsize);
  if (size)
    QFWRITE(padbuf, (size_t)size, outwfield->cat->file, outwfield->filename);

/* Close files */
  close_cat(outfield->cat);
  close_cat(outwfield->cat);
  for (n = 0; n<ninput; n++)
    {
    close_cat(infield[n]->cat);
    close_cat(inwfield[n]->cat);
    }

/* Free Buffers */
  free(emptybuf);
  free(emptywbuf);
  free(multibuf);
  free(multiwbuf);
  free(inbuf);
  free(inwbuf);
  free(outbuf);
  free(outwbuf);

  return RETURN_OK;
  }


/******* coadd_pix ************************************************************
PROTO	int coadd_pix(PIXTYPE *inpix, PIXTYPE *inwpix, int ninput,
		PIXTYPE *outpix , PIXTYPE *outwpix, coaddenum coaddtype,
		PIXTYPE wthresh)
PURPOSE	Coadd images.
INPUT	Input pixel array ptr,
	Input weight ptr array,
	number of input fields,
	Output pixel ptr,
	Output weight ptr,
	Coaddition type.
OUTPUT	RETURN_OK if no error, or RETURN_ERROR in case of non-fatal error(s).
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 15/03/2001
 ***/
int coadd_pix(PIXTYPE *inpix, PIXTYPE *inwpix, int ninput,
		PIXTYPE *outpix , PIXTYPE *outwpix, coaddenum coaddtype,
		PIXTYPE wthresh)

  {
   static PIXTYPE	pixstack[MAXINFIELD];
   PIXTYPE		*pixt;
   double		val, val2, wval,wval2,wval3;
   int			i, ninput2;

  switch(coaddtype)
    {
    case COADD_MEDIAN:
      ninput2 = 0;
      wval = wval3 = 0.0;
      pixt = pixstack;
      for (i=ninput; i--;)
        {
        val2 = *(inpix++);
        wval2 = *(inwpix++);        
        if (wval2 < wthresh)
          {
          *(pixt++) = val2;
          wval += 1.0/sqrt(wval2);
          wval3 += wval2;
          ninput2++;
          }      
        }
      if (ninput2)
        {
        *outpix = hmedian(pixstack, ninput2);
/*------ We assume Gaussian input noise */
        *outwpix = ninput2>2?
		(PI*ninput2*ninput2/(2*wval*wval*(ninput2+((ninput2&1)?(PI/2-1)
			:(PI-2)))))
	  : wval3/4.0;
        }
      else
        {
        *outpix = 0.0;
        *outwpix = BIG;
        }
      break;
    case COADD_AVERAGE:
      ninput2 = 0;
      wval = val = 0.0;
      for (i=ninput; i--;)
        {
        val2 = *(inpix++);
        wval2 = *(inwpix++);
        if (wval2<wthresh)
          {
          ninput2++;
          val += val2;
          wval += wval2;
          }
        }
      if (ninput2)
        {
        *outpix = val/ninput2;
        *outwpix = wval/(ninput2*ninput2);
        }
      else
        {
        *outpix = 0.0;
        *outwpix = BIG;
        }
      break;
    case COADD_MIN:
      ninput2 = 0;
      val = BIG;
      wval = 0.0;
      for (i=ninput; i--;)
        {
        val2 = *(inpix++);
        wval2 = *(inwpix++);
        if (wval2<wthresh && val2<val)
          {
          ninput2++;
          val = val2;
          }
        }
      if (ninput2)
        {
        *outpix = val;
        *outwpix = 1.0;
        }
      else
        {
        *outpix = 0.0;
        *outwpix = BIG;
        }
      break;
    case COADD_MAX:
      ninput2 = 0;
      val = -BIG;
      wval = 0.0;
      for (i=ninput; i--;)
        {
        val2 = *(inpix++);
        wval2 = *(inwpix++);
        if (wval2<wthresh && val2>val)
          {
          ninput2++;
          val = val2;
          }
        }
      if (ninput2)
        {
        *outpix = val;
        *outwpix = 1.0;
        }
      else
        {
        *outpix = 0.0;
        *outwpix = BIG;
        }
      break;
    case COADD_WEIGHTED:
      ninput2 = 0;
      wval = val = 0.0;
      for (i=ninput; i--;)
        {
        val2 = *(inpix++);
        wval2 = *(inwpix++);
        if (wval2<wthresh)
          {
          ninput2++;
          val += val2/wval2;
          wval += 1.0/wval2;
          }
        }
      if (ninput2)
        {
        *outpix = val/wval;
        *outwpix = 1.0/wval;
        }
      else
        {
        *outpix = 0.0;
        *outwpix = BIG;
        }
      break;
    case COADD_CHI2:
      ninput2 = 0;
      wval = val = 0.0;
      for (i=ninput; i--;)
        {
        val2 = *(inpix++);
        wval2 = *(inwpix++);
        if (wval2<wthresh)
          {
          ninput2++;
          val += val2*val2/wval2;
          }
        }
      if (ninput2)
        {
        *outpix = sqrt(val/ninput2);
        *outwpix = 1.0;
        }
      else
        {
        *outpix = 0.0;
        *outwpix = BIG;
        }
      break;
    default:
      error(EXIT_FAILURE, "*Internal Error*: Unknown Combine option in ",
			"coadd_pix()");
    }

  return RETURN_OK;
  }


/******* fast_median **********************************************************
PROTO	PIXTYPE fast_median(PIXTYPE *arr, int n)
PURPOSE	Fast median from an input array, optimized version based on the
	select() routine (Numerical Recipes, 2nd ed. Section 8.5 and
	http://www.eso.org/~ndevilla/median/). If n is even, then the result
	is the average of the 2 "central" values.
INPUT	Input pixel array ptr,
	number of input pixels,
OUTPUT	Value of the median.
NOTES	n must be >0. Warning: changes the order of data (but does not sort
	them)!
AUTHOR	E. Bertin (IAP)
VERSION	19/08/2000
 ***/
#define MEDIAN_SWAP(a,b) { register PIXTYPE t=(a);(a)=(b);(b)=t; }

PIXTYPE fast_median(PIXTYPE *arr, int n) 
  {
   PIXTYPE	*alow, *ahigh, *amedian, *amiddle, *all, *ahh,
		val, valmax, valmax2;
   int		i;

  alow = arr;
  ahigh = arr + n - 1;
  amedian = arr + n/2;
  while (ahigh > (all=alow + 1))
    {
/*-- Find median of low, middle and high items; swap into position low */
    amiddle = alow + (ahigh-alow)/2;
    if (*amiddle > *ahigh)
      MEDIAN_SWAP(*amiddle, *ahigh);
    if (*alow > *ahigh)
      MEDIAN_SWAP(*alow, *ahigh);
    if (*amiddle > *alow)
      MEDIAN_SWAP(*amiddle, *alow);

/*-- Swap low item (now in position middle) into position (low+1) */
    MEDIAN_SWAP(*amiddle, *all);

/*-- Nibble from each end towards middle, swapping items when stuck */
    ahh = ahigh;
    for (;;)
      {
      while (*alow > *(++all));
      while (*(--ahh) > *alow);

      if (ahh < all)
        break;

      MEDIAN_SWAP(*all, *ahh);
      }

/*-- Swap middle item (in position low) back into correct position */
    MEDIAN_SWAP(*alow, *ahh) ;

/*-- Re-set active partition */
    if (ahh <= amedian)
      alow = all;
    if (ahh >= amedian)
      ahigh = ahh - 1;
    }

/* One or two elements left */
  if (*alow > *ahigh)
    MEDIAN_SWAP(*alow, *ahigh);

  if (n&1)
/*-- Odd case */
    return *amedian;
  else
    {
/*-- Even case */
    valmax2 = *amedian;
    valmax = -BIG;
    alow = arr;
    for (i=n;i--;)
      if ((val=*(alow++)) > valmax && val<valmax2)
        valmax = val;
    return (*amedian+valmax)/2.0;
    }

  }

#undef MEDIAN_SWAP


/******************************** hmedian ***********************************/
/*
Median using Heapsort algorithm (for float arrays) (based on Num.Rec algo.).
Warning: changes the order of data!
*/
PIXTYPE   hmedian(PIXTYPE *ra, int n)

  {
   int          l, j, ir, i;
   float        rra;


  if (n<2)
    return *ra;
  ra--;
  for (l = ((ir=n)>>1)+1;;)
    {
    if (l>1)
      rra = ra[--l];
    else
      {
      rra = ra[ir];
      ra[ir] = ra[1];
      if (--ir == 1)
        {
        ra[1] = rra;
        return n&1? ra[n/2+1] : (ra[n/2]+ra[n/2+1])/2.0;
        }
      }
    for (j = (i=l)<<1; j <= ir;)
      {
      if (j < ir && ra[j] < ra[j+1])
        ++j;
      if (rra < ra[j])
        {
        ra[i] = ra[j];
        j += (i=j);
        }
      else
        j = ir + 1;
      }
    ra[i] = rra;
    }

/* (the 'return' is inside the loop!!) */
  }